Efficient analysis of dynamic fracture mechanics in various media by a novel meshfree approach

Accepted Manuscript Efficient Analysis of Dynamic Fracture Mechanics in Various Media by a Novel Meshfree Approach A. Aghahosseini, A. Khosravifard, Tinh Quoc Bui PII: DOI: Reference:

S0167-8442(18)30440-3 https://doi.org/10.1016/j.tafmec.2018.12.002 TAFMEC 2140

To appear in:

Theoretical and Applied Fracture Mechanics

Received Date: Revised Date: Accepted Date:

13 September 2018 6 November 2018 6 December 2018

Please cite this article as: A. Aghahosseini, A. Khosravifard, T. Quoc Bui, Efficient Analysis of Dynamic Fracture Mechanics in Various Media by a Novel Meshfree Approach, Theoretical and Applied Fracture Mechanics (2018), doi: https://doi.org/10.1016/j.tafmec.2018.12.002

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Efficient Analysis of Dynamic Fracture Mechanics in Various Media by a Novel Meshfree Approach A. Aghahosseinia, A. Khosravifarda,1, Tinh Quoc Buib a b

Department of Mechanical Engineering, Shiraz University, Shiraz 71936, Iran

Department of Civil and Environmental Engineering, Tokyo Institute of Technology, 2-12-1W8-22, Ookayama, Meguro-ku, Tokyo 152-8552, Japan

Abstract A novel approach for meshfree analysis of various fracture mechanics problems is introduced. This approach is based on a weak-form meshfree method, which utilizes a refined nodal configuration in the vicinity of the crack tip and also a special meshfree integration technique. In this method, not only the computational cost is reduced, but also superior accuracy and convergence is achieved. In the proposed method, instead of enriching the basis functions for modelling the singular behavior of the stress fields near the crack tip, a few nodes are added to this region, and then by applying an optimal domain integration method a solution with high speed and accuracy is obtained. In the present work, stress intensity factors (SIFs) are used as a criterion to analyze crack’s behavior under dynamic loading conditions. The domain form of the interaction integral is used for calculation of the SIFs. Several benchmark problems are analyzed by the proposed method and the results are compared with those obtained from other methods. The example problems cover a large range of material behavior, from homogenous to nonhomogeneous, and from isotropic to generally anisotropic. Validity of the results obtained for the dynamic stress intensities of cracked bodies asserts the usefulness of the proposed method for the efficient and accurate numerical analysis of fracture mechanics problems in homogenous, nonhomogeneous, isotropic and anisotropic media in dynamic conditions. Keywords: Fracture, Dynamic stress intensity factor, Meshfree, Radial point interpolation method, Background decomposition method

1. Introduction The presence of flaws leads to considerable reduction in the strength and integrity of a structural component under various loadings. Cracks are among the most common flaws that may initiate during the manufacturing stage, or later under working conditions of the structure.

1

Corresponding author at: Department of Mechanical Engineering, Shiraz University, Shiraz 71936, Iran. E-mail address: [email protected]

1

Stress and displacement fields that arise in a structure have remarkable differences under static and dynamic loadings. Therefore, elastodynamic analysis of cracks has crucial importance particularly in bodies that experience instant loadings. Generally, the problems of the dynamic fracture mechanics are divided into two major categories: the first one analyzes the behavior of stationary cracks under dynamic loading and the other one considers conditions that the crack reaches the point of instability and starts to grow. The present study focuses on the former type of dynamic fracture. The rapid development and practical importance of newfound materials in several industries such as aerospace, transportation, and military attracts the interests toward the investigation of mechanical behavior of such media in the presence of cracks under static and dynamic loadings. Early scientific investigations on fracture phenomenon are certified to Inglis [1] and Griffith [2]. During his researches, Irwin [3] realized that stress fields near the crack tip in various conditions have the same order of singularity; consequently, he developed a new concept for characterizing crack’s behavior which is called stress intensity factor and since then became very useful in fracture analysis. Although several kinds of research have been done on crack behavior by analytical methods [49], this approach is limited to cases with simple loading and basic geometries. Due to difficulties of obtaining analytical solutions in general conditions, particularly in intricate engineering problems, numerical methods are among the most widely used approaches. Finite elements is a powerful tool for modeling complicated solid mechanics problems, however, its dependency on the quality of the elements, induces challenges in the process of analysis of fracture mechanics problems [10]. To overcome this obstacle, novel numerical methods based on finite elements were developed; ES-FEM [11-14] and XFEM [15-18] are among the most popular methods that are used for analysis of cracked bodies. The Extended isogeometric analysis (XIGA) is also a powerful tool for analysis of the problems that have flaws in their domain. Bhardwaj et al. [19, 20] studied the fatigue life of an interfacially cracked plate in the presence of flaws by homogenized XIGA. Other methods such as the boundary element method (BEM) [21-24] and meshfree methods [25-28] are also developed to eliminate the need of domain meshing. Early studies on fracture mechanics by meshfree methods is accredited to Belytschko et al. [29]. They employed MLS interpolation to create shape functions and analyzed static and propagating cracks in isotropic materials. Using improved shape functions with the aid of enriched basis functions is a convenient method to capture singular behavior of the stress fields at the crack tip. Several studies have been performed on basis functions to find a compromise between accurate results and computational efficiency. Wen and Aliabadi [28] used enriched basis functions in a meshfree Galerkin method to evaluate SIFs in a cracked body under static and dynamic loadings. Peng et al. [30] analyzed cracks in stiffened plates by using an improved meshless method and investigated the stress and displacement fields in a plate in presence and absence of the stiffener. Pathak et al. [31] simulated the three-dimensional interfacial crack growth by coupled finite element and element free Galerkin (FE-EFG) approach under mechanical and thermoelastic loading. 2

Cai et al. [32] employed Shepard function and partition of unity (MSPU) to introduce an improved meshless method for evaluation of SIFs and prediction of propagation path of cracks. Pant et al. [33] implement the EFG method for the stress analysis of structures having cracks at the interface of two dissimilar materials. The material discontinuity at the interface was modeled using a jump function, which enriches the approximation by the addition of special shape functions that contain discontinuities in their derivative. Pathak et al. [34] studied quasi-static fatigue crack growth of homogeneous and bi-material interfacial cracks by the EFG method under mechanical and thermo-elastic loads. In their study, discontinuities in the temperature and displacement fields were captured by extrinsic partition of unity enrichment technique. Applying the MLS interpolation for calculation of shape functions in the EFG method, due to its lack of Kronecker delta property, makes it impossible to enforce essential boundary conditions directly. Various strategies have been considered to overcome this shortcoming. For instance, Rao and Rahman [35-36] introduced a new weight function for accurate implementation of essential boundary conditions in the problem of evaluating the SIFs in cracked bodies with various types of materials. Penalty or Lagrange multiplier techniques are also used to enforce essential boundary conditions, however, obtaining proper results by these methods require further computational effort [26]. In addition, alternative interpolation techniques such as the radial point interpolation method (RPIM) and the moving Kriging method have been developed, which possess the Kronecker delta function property. RPIM is an interpolation technique that is widely used in various engineering fields. This method was introduced by Liu and Gu to analyze solid mechanics problems [37-38] and was later extended to study fracture mechanics problems [39]. Farahani et al. [40] used the RPIM and the FEM to obtain stress distribution and to calculate SIF for a compact tension specimen (CT) during fatigue crack growth test with thermoelastic stress analysis. A new meshless sub-region radial point interpolation method (MS-RPIM) was proposed by Zhuang et al. [41] to obtain the SIFs directly by adding a few nodes around the crack tip. Nguyen et al. [27] proposed a new approach based on the local partition of unity for modeling of crack growth in elastic solids. In their method, RPIM was utilized for producing the shape functions. The natural neighbor radial point interpolation method (NNRPIM) was utilized by Azevedo et al. [42] to predict crack propagation path in elastic solids. Nguyen et al. [26] studied dynamic fracture mechanics problems for isotropic and orthotropic materials using extended radial point interpolation method (X-RPIM). In their study, uniform distribution of nodes on the domain and the enriched basis functions for nodes around the crack tip were used. Using improved shape functions in a meshfree method, Bui et al. [43] calculated fracture parameters in the functionally graded materials under dynamic loading. In all of the mentioned studies, a somewhat complicated process is proposed to model the singular stress fields at the crack tip. Khosravifard et al. [10] used the meshfree RPIM to investigate stationary and propagating cracks in 2D bodies. In their study, instead of employing enriched basis functions to capture stress singularities, a few nodes were added around the crack tip, and with the assistance of an accurate domain integration method, called the background 3

decomposition method (BDM), accurate results were attained. In this method, a non-uniform nodal distribution in the entire domain is used, therefore, adopting RPIM is a suitable choice due to its low sensitivity to nonuniform arrangements and density of nodes in the domain. The BDM is a domain integration technique, proposed by Hematiyan et al. [44] which is most useful for problems with variable nodal density in their domains. The meshfree method proposed by Khosravifard et al. [10] was utilized for analysis of quasi-static fracture mechanics problems in isotropic media . In the present work, this technique is extended for determination of dynamic SIFs in a 2D stationary crack under dynamic loadings in isotropic, anisotropic and FG media. In the proposed technique, only a few nodes are added near the crack tip, and the domain integrals are efficiently computed by the BDM. In this way, not only the computational cost is minimized, but also accurate results are obtained. Solution of several benchmark problems by the proposed meshfree approach demonstrates that the total number of nodes required to reach a specific amount of accuracy in this method is less than other novel numerical methods. The rest of this paper is organized as follows: In Section 2, a brief review of fracture mechanics in various media is presented. The formulation of the RPIM for dynamic analysis of fracture mechanics is given in section 3. To assess validity and adequacy of the proposed method, several example problems are analyzed in Section 4. Eventually, some conclusions drawn from this study are stated in Section 5.

2. A brief review of dynamic fracture mechanics in anisotropic and nonhomogeneous media 2.1

Governing equations

Fig. 1 depicts a typical cracked body and the corresponding boundary conditions in 2D dynamic fracture mechanics. The equations of motion, boundary conditions, and initial conditions for homogeneous materials are written as:

(1)

4

Fig. 1: A 2D body with a traction-free crack and the associated coordinate systems. where the physical domain is enclosed by the boundary . is the traction boundary and is the displacement boundary. In Eq. (1), is the stress tensor, is the body force per unit volume, and is the displacement vector. and are pre-defined displacement and traction vectors, respectively, and is the outward unit normal to the boundary. and are density and acceleration, respectively, and and are the initial displacement and velocity. In order to use these equations for a non-homogenous medium, spatial variations of the material properties should be considered. For 2D elastic solids, the stress-strain relations are determined by the generalized Hook's law as: (2) where are components of the fourth-order elasticity tensor which after simplification for 2D isotropic elastic bodies under plane stress conditions can be written in matrix form as: (3)

and in Eq. (3) are the Young’s modulus and Poisson’s ratio respectively. For isotropic solids in the plane strain condition, the matrix of material constants D can be written by replacing and , respectively, with and . The elasticity tensor for non-homogenous and orthotropic materials under plane stress state can be expressed as follows, respectively: (3)

5

(4)

Eqs. (3) and (4) present the elasticity tensor for non-homogenous and orthotropic materials, respectively. In Eq. (4), the subscripts 1 and 2 denote the principal material directions. 2.2

Evaluation of the dynamic stress intensity factors

The stress intensity factor can be considered as one of the most important parameters in the analysis of fracture mechanics, by means of which not only the stress and displacement fields near the crack tip can be characterized, but also it can be used as a criterion for crack growth and prediction of the crack propagation path. There are several methods for numerical computation of the stress intensity factor, such as the displacement correlation technique (DCT) [45], the virtual crack closure technique (VCCT) [46], conventional and extended path independent J integral [27] and interaction integral [47]. In the present work, the domain form of the interaction integral, which is derived from the J integrals is employed. The concept of the path independent J-integral was developed by Rice in 1968 [48] for a non-linear homogeneous elastic material as follows: (5) In Eq. (5), is an arbitrary closed contour around the crack tip, is the strain energy density and for linear elastic materials can be written as , is the Kronecker delta, refers to the components of the outward unit normal vector. All the quantities in Eq. (5) are defined in the local crack tip coordinates (see Fig. 2) with the -axis parallel to the crack face.

Fig. 2: The arbitrary closed contour around the crack tip for evaluation of the J-integral. Although the J-integral is a very successful and efficient method in many problems, its conventional form has some drawbacks, such as the necessity to compute singular stresses at the crack tip, simultaneous dependency on both SIFs in Mode-I and II, and being path independent for functionally graded materials. To overcome these drawbacks, the interaction integral [49] and the extended J-integral, called J* were developed [50]. 6

The two-state interaction integral was proposed by Yau et al. [51] to evaluate SIFs in homogenous structures. The two states correspond to actual and auxiliary states. Crack tip asymptotic stress and displacement fields in the auxiliary problem associated with each mode are available in the literature. The interaction integral for homogeneous materials under static loading is written as: (6) In Eq. (6), the superscripts (1) and (2) denote the actual and auxiliary states, respectively. The auxiliary state can be the either of the modes I or II. With the aid of a smoothing function , which takes a value of unity on the inner contour and vanishes on the outer contour and the use of divergence theorem, the domain form of the interaction integral for homogeneous and isotropic media is obtained as: (7) For the case of dynamic fracture mechanics, inertia effects play an important role in the crack behavior and the interaction integral is written as: (8) In the case of fracture mechanics in heterogeneous materials, using developed auxiliary fields for homogeneous materials causes contravention of one of the basic mechanics laws, i.e., equilibrium, compatibility or constitutive laws, that leads to three independent formulations: “non-equilibrium,” “incompatibility,” and “constant-constitutive-tensor”. Kim and Paulino [52] employed these basic laws to obtain three stable forms for the interaction integral in functionally graded materials. Their study demonstrates that among these three formulations those, which are based on the non-equilibrium and incompatibility relations, are more accurate than the one based on the constant-constitutive-tensor. In the present work, for the sake of simplicity the nonequilibrium formulation of the interaction integral is taken for dynamic cases.

(9)

For isotropic media the following equation is used to calculate the SIFs using the interaction integral value. (10) Where and are the stress intensity factors corresponding to mod-I and II, respectively, and for plane stress and strain is equal to ,and , respectively. The SIFs in Mode-I 7

and II respectively can be obtained by considering and then using the following relation:

,

and

,

, (11)

The relation between interaction integral and the SIFs for orthotropic materials can be expressed as: (12) where

(13)

In Eq. (13), are components of the compliance matrix and following equation [23]:

are non-conjugated roots of (14)

Finally, the SIFs are obtained from the following relation: (15)

3. Formulation of the meshfree RPIM for dynamic fracture mechanics analysis Herein, the meshfree RPIM is used as the numerical method to obtain the solution of fracture mechanics problems. In this method, radial basis functions along with the polynomial basis functions are used to approximate the field variables and produce the shape functions. By using this interpolation method, the displacement approximation for a typical point can be written as: (16) In this equation, and are radial and polynomial basis functions, respectively. is the number of nodes in the support domain of point , and is the number of monomial terms; finally, and are unknown constants.

8

Several radial basis functions can be found in the literature. In the present study, the thin plate spline (TPS) function is employed: (17) where is the distance between the point of interest and a node in the support domain and the proposed values for are 3.001, 4.001, and 5.001 [53]. After performing a series of mathematical operations on Eq. (16), the following relation is obtained for the displacement field in terms of the nodal values of the function and the shape functions: (18) where is the vector of shape functions in the meshfree RPIM and is obtained from the following relation: (19) The moment matrix

is expressed as:

where

(20)

Using the Hamilton’s principle, the weak form of Eq. (1) can be written as:

Γ

Γ

(21)

Where is variation of the displacement vector u and is the applied traction vector. Subsequently, by substitution of the approximate displacement field, Eq. (18), into the Galerkin weak form, Eq. (21), the discrete form of the governing equations for the meshfree RPIM [26] is obtained as: (22) where u is the nodal displacement vector and the other matrices and vectors are expressed in the following expressions:

9

(23) (24) (25) where

(26)

Herein, for solution of the obtained system of ordinary differential equations, Eq. (22), the Newmark time integration scheme is adopted. In this method, the solution procedure starts by writing Eq. (22) for each time step, i.e.: (27) Then, by applying the Newmark algorithm to Eq. (27), the following algebraic system of equations at time is obtained: (28) in this relation, is the time step size. is calculated from Eq. (28) and subsequently, the acceleration and velocity vectors are calculated at the same time step by the following relations: (29) (30) In every new time-step, the displacements, velocities, and accelerations are calculated based on their previous time-step values. Considering and guarantees the unconditional stability of the Newmark method. 3.1

The Background decomposition method for evaluation of the domain integrals

As previously stated, in the present work, a few nodes are added to the region around the crack tip in order to capture the singular fields at the crack tip. Therefore, the density of the nodes varies in different parts of the domain. Obviously, it is necessary to take a special strategy for efficient and accurate evaluation of the domain integrals, without using too much integration points. Herein, the BDM is used to perform the domain integrals of the weak-form. In this method, the distribution of integration points conforms to that of the nodes in the domain. In other words, where the nodes are closer to each other, automatically a denser distribution of the 10

integration points will be created. This technique is very useful for analysis of problems such as fracture mechanics, stress concentration, and problems including concentrated loads. Utilizing BDM for evaluation of domain integrals was proposed by Hematiyan et al. [44] for weak-form meshfree methods with non-uniform nodal distribution in the problem domain. Calculation of the integrals using the BDM in a domain that has significant changes in nodal density is done in four steps. In the first step, scattered domain nodes are categorized into several groups based on their local density. Subsequently, the domain is divided into several sections such that the nodal distribution in each section is uniform. In the third step, based on the local density of the nodes, the integration points and their weights are determined and finally, the global vectors for the calculation of the domain integral are obtained in the fourth step [44].

4. Numerical results and discussions In this section, to examine the accuracy and efficiency of the proposed method, seven numerical example problems are presented. The first three cases deal with the evaluation of the DSIFs in isotropic and homogeneous media. The next two examples determine DSIFs in isotropic and non-homogenous materials (i.e., in FGMs), and the last two ones deal with evaluation of the DSIFs in orthotropic media and investigate the effect of inclination angle between the principal material axes and the crack-line in such media. It should be mentioned that for analysis of the example problems, the codes developed by the authors in MATLAB are used. For the numerical simulations, step and ramp loadings, as depicted in Fig. 3, are adopted. In fact, in these analyses, the purpose of the study is to investigate the time variations of the DSIFs in cracked bodies during the loading time. All of the results are validated in comparison with those obtained from the literature.

Fig. 3: (a) the Step loading and (b) the ramp loading, used in the analyses. 4.1

Example 1: A benchmark example problem

The purpose of this example is to investigate the accuracy of the results obtained by the present method. To this end, an example problem for which analytical solutions can be found in the literature is adopted. The example problem deals with a semi-infinite crack in an infinite isotropic and homogeneous domain that is subjected to tensile stresses at the top and bottom 11

edges. In the numerical procedure, the infinite domain is modeled by a 10×4 rectangle. The geometry, loading and boundary conditions are shown in Fig. 4. As shown in Fig. 4(b), due to the symmetry, only half of the problem domain is modeled in the numerical analysis. The analytical solution for mode-I DSIF is valid only for , i.e., up to the time when the reflected stress wave reaches the crack tip. The parameter is the time at which the first stress wave strikes the crack tip and is the dilatational wave speed. In the case of step loading, the DSIF can be obtained by the following equation as the stress wave reaches the crack tip at [55]: (31) In this problem, the state of plane strain prevails and the material properties are as follows: (32)

Fig. 4: Geometry and loading condition of an edge-cracked specimen, Example 1, (a) the problem domain, (b) the computational domain. Firstly, a convergence study is performed in order to investigate the sensitivity of the results with respect to the distribution of the nodal points. The domain is modeled with three sets of nodal arrangements as shown in Fig. 5. The number of nodes in these arrangements is 157, 304, and 515. It is seen that in all of the nodal arrangements, the density of the nodal points vary drastically in different parts of the domain.

12

Fig. 5: Distribution of the nodes in the domain of example 1, (a) arrangement 1, with 157 nodes, (b) arrangement 2, with 304 nodes, (c) arrangement 3, with 515 nodes. The normalized mode-I DSIF ( ) versus the normalized time ( ) is calculated by the proposed method, and the results are compared with the analytical solution. Fig. 6 illustrates the convergence of the computed DSIFs to the analytical solution, as the number of nodes increases. It is interesting to see that even with a small number of nodes, acceptable results are obtained. It is also seen that the calculated DSIFs converge quickly to the analytical solution.

Fig. 6: The mode-I DSIF versus time for different nodal arrangements in Example 1. 13

In Table 1, the overall errors of the computed DSIFs are reported according to the following expression: (33) Where and are mode-I DSIFs evaluated by the analytical solution and the proposed method, respectively, and M is the number of time steps. Table 1: The overall errors of mode-I DSIF for different nodal arrangements. Arrangements

Number of nodes

Error

Arrangement 1

157

6.19×10–4

Arrangement 2

304

6.68×10–5

Arrangement 3

515

4.64×10–5

Fig. 6 and Table 1 clearly show the excellent agreement between the results obtained by the proposed method and the analytical solution. Although the results obtained by all nodal arrangements indicate an acceptable accuracy, they become more accurate when more nodes are utilized in the meshfree method. Fig. 6 implies that by using the second nodal arrangement, efficiently accurate results can be obtained. As the result, this nodal arrangement will be used for further analyses of this example problem. Many other researchers have studied this problem using other numerical methods. In Fig. 7 results obtained by the present method are compared with those obtained by the XFEM [56], and the meshfree XRPIM based on linear ramp function [26].

Fig. 7: The mode-I DSIF versus time for different numerical methods in Example 1. 14

From Fig. 7 it can be observed that the results obtained by the present method match the analytical solution most closely in comparison with other numerical techniques. For instance, compared with the analytical solutions, the present method gives a smaller error near than the other approaches. Remarkably, the results of the present work are obtained by only 304 nodes, while 2400 nodes were used in the XRPIM, and 4000 nodes in the XFEM. Finally, the effect of loading types and load intensities on mode-I DSIF is investigated. In the ramp loading, the applied traction to the boundary varies according to the following relation: (34) In Eq. (34), is the rise time of the ramp loading. The DSIF can be obtained by the following equation as the stress wave strikes the crack tip at [57]:

(35)

In Eq. (35), . The same problem with ramp loading is analyzed by the present method and the results are depicted in Fig. 8. An excellent agreement with the analytical solution is observed. In the analyzed problem, the rise time of the ramp loading is .

Fig. 8: The mode-I DSIF versus time for ramp loading in Example 1. Lastly, the problem is solved with various load intensities for both loading types. The selected load intensities are 25 MPa, 50 MPa, and 100 MPa. The DSIFs for the step and ramp loading are plotted in Fig. 9(a) and 9(b), respectively. As expected, as the intensity of the applied loading increases, so does the value of the DSIF. Fig. 10 compares the DSIFs obtained from step and

15

ramp loading with different rise times, and with the same load intensity. It is observed that the step loading produces larger DSIFs in comparison with the ramp loading with any rise time.

Fig. 9: The mode-I DSIF versus time for (a) step and (b) ramp loading with different load intensity in Example 1.

Fig. 10: Comparison of the DSIF values corresponding to step and ramp loading, Example 1. 4.2

Example 2: Mode-I DSIF in an isotropic medium

The second numerical example for dynamic fracture analysis of isotropic materials deals with a rectangular body with a central crack, which is subjected to a tensile step loading , where is the loading amplitude and is the Heaviside function, at the top and bottom edges. Fig. 11 depicts the geometry, boundary conditions, and the applied loading of this example. 16

This example was previously studied by Chen [58] and since then, it has been regarded as a benchmark problem. Chen evaluated the mode-I SIFs for a center cracked tension specimen by using a time-dependent finite difference method (FDM). The geometry and loading condition are shown in Fig. 11. In this example, plane strain condition is assumed and the material and geometrical properties are as follows: (36)

Fig. 11: Geometry and loading of a center-cracked specimen, Example 1. The problem domain is modeled by a total of 360 nodes with sufficient nodal refinement at the crack tips. Fig. 12(a) depicts the overall distribution of the nodes within the problem domain. A closer view of the nodal distribution near the crack tips is shown in Fig. 12(c). The BDM is used to calculate the domain integrals. Fig. 12(b) shows distribution of the BDM integration points for the nodal arrangement shown in Fig. 12(a). It can be seen that the density of the integration points is consistent with the density of the nodes in the domain. Fig. 12(d) shows a closer view of the distribution of the integration points near the crack tips. It is clear that the density of the integration points near the crack tips is larger than other parts of the domain.

17

Fig. 12: Distribution of the nodes and integration points in example 2; (a) overall nodal distribution in the domain, (b) distribution of the BDM integration points, (c) nodal refinement at the crack tips and (d) BDM integration points at the crack tip. In the process of computing mode-I DSIF, the Newmark method with time steps of

was

used to obtain the field variables. The obtained DSIF is normalized by and plotted over normalized time in Fig. 13, where is the longitudinal wave velocity and is written as: (37) The obtained results are compared with those obtained with the FDM [58] and also the BEM [23]. The oscillations of the SIF through the loading time can be attributed to the changes in the stress field at the crack tip due to the interaction between elastic waves and geometrical boundaries or crack faces. In order to discretize the domain with the BEM, 34 elements have been used for the internal domain, 24 elements for the external boundary and 10 elements for the crack, while a gird of 50 zones/cm was used in FDM.

18

Fig. 13: the mode-I dynamic SIF versus time in Example 2. According to Fig. 13, both of the normalized DSIFs are zero at the beginning of the dynamic loading until at which the elastic waves reach the crack tip, suddenly, the mode-I DSIF increases, while mode-II DSIF remains close to zero. This implies that in the present example the opening mode dominates. This occurs due to the normal angle between the crack and direction of the applied loads. Because of the geometrical and loading symmetry, and because of material homogeneity, the obtained DSIFs at both crack tips are almost identical. As illustrated in Fig. 13, the DSIF calculated by the present method has a good agreement with those obtained with other methods. 4.3

Example 3: Mixed-mode DSIFs in a slanted edge-cracked rectangular domain

This example of dynamic fracture analysis in an isotropic material is dedicated to a mixed-mode fracture problem. For this purpose, a rectangular domain that includes a slanted edge crack at an angle of with the x-axis is considered. As shown in Fig. 14, three edges of the domain are supported with rollers, and a uniform tensile load of , see Fig. 3, is applied to the other edge. The material and geometrical properties of this example are: (38)

19

Fig. 14: Geometry, boundary conditions and loading of slanted edge cracked specimen, Example 3. As shown in Fig. 15, a set of 250 nodes are used to model the problem domain and operation time is set to . The proposed meshfree RPIM is employed to solve the mixed-mode fracture problem and the obtained numerical results of DSIFs for both mode-I and mode-II are depicted in Fig. 16. In this figure, results of the same problem that have been obtained by the EFG method [28], and the dual BEM [22] are also shown.

Fig. 15: Distribution of nodes on problem’s domain. As shown in Fig. 16, there is a close agreement between results of the proposed method and those of the DBEM [22], while the results of the EFG method [28] have a small deviation. It is noteworthy to mention that the results of the present method have been obtained with only 250 20

nodes and these results are more accurate than those obtained by 600 EFG nodes of [28]. Also Nguyen et al. [26], utilized 2790 meshfree nodes to obtain the same level of accuracy. Therefore, it is concluded that by using the proposed method, the computational cost of analysis of fracture mechanics problems can be reduced considerably.

Fig. 16: Time variations of both Mode-I and Mode-II DSIFs in example 3. To study the effect of the angle between the crack line and the force direction on the variation of the stress intensity factor, we consider an oblique center-cracked body under the action of an axial tensile loading (see Fig. 11). Results for mode-I DSIFs with different crack angles are represented in Fig. 17. This figure shows a significant effect of the crack angle on the DSIFs. In particular, the results indicate that by increasing the angle between the crack and force direction from 0 to 90 degrees, the effect of the elastic wave propagation on the fluctuations of the DSIF decreases in mode-I, which might be due to the difference in the reflection angle of these waves from the crack surface.

21

Fig. 17: Effect of crack inclination angle on the time variations of the DSIF, example 3.

4.4

Example 4: Mode-I DSIF in non-homogenous media

In this example, dynamic behavior of a cracked isotropic and non-homogeneous medium is studied. As shown in Fig. 11, a rectangular domain with a center crack under the application of uniform tensile loads is considered. Similar to the previous examples, the applied load is , see Fig. 3. The medium is made of an FGM, for which the Young's modulus varies in both directions, while the Poisson’s ratio is a constant, i.e. . The material and geometrical properties of this example are:

(39) where and respectively, are the Young’s modulus and mass density of Example 2 (see Eq. (36)). is the non-homogeneity parameter, i.e., represents a homogeneous material. Plane strain condition is also assumed for this case. The numerical model and the nodal arrangement are the same as those of Example 2. This problem is examined to demonstrate the applicability of the present method in problems that deal with non-homogenous media, especially FGMs. The computed DSIFs are compared with those available in the literature. Because of the non-homogeneity of the medium, the values of the DSIFs at the two crack tips are not the same. The time variations of the DSIFs at the right crack tip are depicted in Fig. 18. In this figure, results of the FEM [49], obtained by a mesh of 816 Q8 and 142 T6 2D plane strain elements and X-RPIM [43] with 1800 nodes are also shown. A good agreement between the results is seen.

22

Fig. 18: Time variations of the DSIFs at the right crack tip in a FG medium, example 4. Fig. 19 compares the time variations of the DSIFs at the right and left crack tips. It is seen that the non-homogeneity of the medium has led to different stress intensities at the two crack tips. Also, by comparing the results of this example with those of example 2, it is concluded that unlike the case of the homogeneous medium, in the FG medium, the mode-II DSIF has a nonzero value. An important result of this phenomenon is that the direction of crack growth in homogeneous and nonhomogeneous media under the same geometry and loading conditions, is not the same.

Fig. 19: Comparison of the time variation of the DSIFs at the right and left crack tips.

23

4.5

Example 5: Mixed-mode DSIFs in non-homogenous materials

This example of dynamic fracture analysis in non-homogenous media deals with a mixed-mode fracture problem. A rectangular FG domain with an inclined center crack at an angle of with the -axis is considered. As shown in Fig. 11, this specimen is subjected to a uniform tensile stress at the top and bottom edges. Young's modulus and mass density have an exponential variation along the -axis with constant ratio. The material and geometrical properties of this example are:

(40)

As shown in Fig. 20, a set of 350 nodes are utilized for modeling of the problem domain. Also, the time step size of the Newmark method is . The normalized DSIFs at the right crack tip in comparison with the reported results using the FEM by Song and Paulino [49] are shown in Fig. 21. It is worth mentioning that the FEM results were obtained by a mesh of 206 Q8 and 198 T6 elements. Therefore, the total number of FEM nodes is much more than that of the present meshfree technique.

Fig. 20: Distribution of the nodes in problem’s domain, example 5.

24

Fig. 21: Time variation of mode-I and II normalized DSIFs at the right crack tip of example 5. Here again, the non-homogeneity of the medium has led to different DSIFs at the two crack tips. Fig. 22 depicts the time variations of the normalized DSIFs for the right and left crack tips.

Fig. 22: Dynamic stress intensity factors for both crack tips, example 5.

4.6

Example 5: Mode-I DSIF in an orthotropic medium

In this example, a dynamic fracture mechanics problem in an orthotropic medium is studied. As shown in Fig. 11 a rectangular domain with a central crack subjected to a uniform tensile stress at the top and bottom edges is considered in this example. The principal material axis and the crack’s orientation are the same. Geometrical properties are the same as Example 2, and material properties are given as: 25

(41) where, and are the Young’s modulus in the principal material axes, is the shear modulus, and is the Poisson’s ratio. This problem is analyzed using the same nodal arrangement as the one used for example 2, see Fig. 12, and the time step size of the problem solution is . Fig. 23 shows the normalized DSIF versus normalized time velocity along the material axis 2.

, where

is the longitudinal wave

Fig. 23: Time variation of mode-I DSIF in the orthotropic medium of Example 6. In Fig. 23, the mode-I DSIF calculated using other numerical methods such as the BEM by Sánchez et al. [23], and also Wünsche et al. [59] and the DBEM by Albuquerque et al. [60] are also shown. Here again, a very close agreement between the results of different numerical approaches can be seen. Considering the homogeneity of the orthotropic medium in this example, the computed DSIFs at the both crack tips are the same.

4.7

Example 7: Mixed-mode DSIFs in an orthotropic medium

In the final example, an orthotropic rectangular domain with various crack inclinations is considered (see Fig. 24). In this case, material and spatial coordinates (body coordinate) are not aligned, and therefore, the material behaves as an anisotropic medium. The plane stress condition is assumed for this case. The geometrical and material properties of this example are the same as the previous example.

26

Fig. 24: Geometry and boundary conditions of center cracked specimen in an anisotropic medium, example 7. Fig. 25 depicts the normalized DSIFs of this example for the case in which the angle between the principal material direction and the crack line is . In this figure, the results of a BEM analysis performed by Sánchez et al. [23] are also shown. Here again, a very good agreement between the results of the two methods is observed.

Fig. 25: Variation of mode-I and II normalized DSIFs of the anisotropic cracked medium, Example 7. 27

In order to investigate the effect of crack inclination with respect to principal material direction on the dynamic stress intensity factors, various angles from to are considered. The mode-I and II DSIFs are plotted in Fig. 26 and Fig. 27 respectively. These two figures show that material anisotropy has a great influence on the time variations and magnitude of the DSIFs. For instance, when the crack line and the principal material direction are aligned, the mode-II DSIF vanished and the problem is in pure mode I. However, for other inclination angles, the problem is in a mixed mode condition.

Fig. 26: Mode-I DSIFs for various crack inclination angles with respect to principal material direction.

Fig. 27: Mode-II DSIFs for various crack inclination angles with respect to principal material direction.

28

5. Conclusions and outlook This paper presents an improved meshfree method for numerical analysis of dynamic fracture mechanics problems in various media, including isotropic, homogeneous, nonhomogeneous, and anisotropic media. In this method, a few nodes are added to the area near the crack tip in order to capture the singular behavior of the stress field. The domain integrals of the weak-form are computed with the efficient BDM. In this integration technique, the density of the integration points conforms to that of the nodal points. The meshfree RPIM is used to discretize the governing equations. This method possesses the kronecker-delta function property and therefore, imposition of the essential boundary conditions is straightforward. As shown in various examples, the prominent features of this method are its independency from the nodal distribution in problem’s domain and the decrease in the number of nodes required for accurate prediction of the field variables in comparison with other numerical methods such as the finite element, extended finite element, and other meshfree methods. Using fewer nodes will reduce the computational cost of obtaining appropriate results, and that is considered as an advantage for this approach in dealing with problems that have complexity in their geometry or loading conditions. By using the BDM for computation of the domain integrals in the weak-form of the governing equations, an accurate solution for determining fracture parameters is achieved. Proper distribution of the integration points along with the appropriate variation of the density of nodes in the problem’s domain and accurate identification of the crack boundaries are among the advantages of the present approach. Furthermore, the following remarks can be stated:  In an isotropic medium with a slanted crack, by increasing the angle between the crack and force direction, the effect of the elastic wave propagation on the fluctuations of the DSIF in mode-I decreases.  The inhomogeneity of the FGMs has a great influence on the time variation of the DSIFs, similar crack tips a nonhomogeneous medium can have completely different SIFs.  Even in geometrically symmetric problems with symmetric loading, and when the load is perpendicular to the crack line, unlike the homogeneous materials, the mode-II DSIF does not vanish in FGMs. Therefore, material inhomogeneity not only changes the values of the DSIF, but also changes the crack growth path.  The angle between the principal material direction and the crack line has a strong effect on the DSIFs in an orthotropic medium as it plays an important role in the time variation of the DSIFs and determination of fracture modes that occur in anisotropic specimens.

29

References 1. Inglis, C. E. (1913). Stresses in a plate due to the presence of cracks and sharp corners. Transactions of the institute of naval architects, 55(219-241), 193-198. 2. Griffith, A. A. (1921). The phenomena of rupture and flow in solids. Philosophical transactions of the royal society of london. Series A, containing papers of a mathematical or physical character, 221, 163-198. 3. Irwin, G. R. (1957). Analysis of stresses and strains near the end of a crack traversing a plate. Journal of applied mechanics, 24(3), 361-364. 4. Sih, G. C., Paris, P., & Irwin, G. R. (1965). On cracks in rectilinearly anisotropic bodies. International Journal of Fracture Mechanics, 1(3), 189-203. 5. Rice, J., & Rosengren, G. F. (1968). Plane strain deformation near a crack tip in a powerlaw hardening material. Journal of the Mechanics and Physics of Solids, 16(1), 1-12. 6. Freund, L. (1976). Dynamic crack propagation(in elastic solids). Mechanics of fracture, 105-134. 7. Delale, F., & Erdogan, F. (1983). The crack problem for a nonhomogeneous plane. Journal of Applied Mechanics, 50(3), 609-614. 8. Viola, E., Piva, A., & Radi, E. (1989). Crack propagation in an orthotropic medium under general loading. Engineering Fracture Mechanics, 34(5-6), 1155-1174. 9. Konda, N., & Erdogan, F. (1994). The mixed mode crack problem in a nonhomogeneous elastic medium. Engineering Fracture Mechanics, 47(4), 533-545. 10. Khosravifard, A., Hematiyan, M., Bui, T., & Do, T. (2017). Accurate and efficient analysis of stationary and propagating crack problems by meshless methods. Theoretical and Applied Fracture Mechanics, 87, 21-34. 11. Jiang, Y., Liu, G., Zhang, Y., Chen, L., & Tay, T. (2011). A singular ES-FEM for plastic fracture mechanics. Computer Methods in Applied Mechanics and Engineering, 200(45), 2943-2955. 12. Liu, G., Nourbakhshnia, N., & Zhang, Y. (2011). A novel singular ES-FEM method for simulating singular stress fields near the crack tips for linear fracture problems. Engineering Fracture Mechanics, 78(6), 863-876. 13. Liu, P., Bui, T., Zhang, C., Yu, T., Liu, G., & Golub, M. (2012). The singular edge-based smoothed finite element method for stationary dynamic crack problems in 2D elastic solids. Computer Methods in Applied Mechanics and Engineering, 233, 68-80. 14. Nguyen-Xuan, H., Liu, G., Nourbakhshnia, N., & Chen, L. (2012). A novel singular ESFEM for crack growth simulation. Engineering Fracture Mechanics, 84, 41-66. 15. Bhardwaj, G., Singh, S., Singh, I., Mishra, B., & Rabczuk, T. (2016). Fatigue crack growth analysis of an interfacial crack in heterogeneous materials using homogenized XIGA. Theoretical and Applied Fracture Mechanics, 85, 294-319. 16. Menouillard, T., & Belytschko, T. (2010). Dynamic fracture with meshfree enriched XFEM. Acta mechanica, 213(1), 53-69. 17. Motamedi, D., & Mohammadi, S. (2010). Dynamic analysis of fixed cracks in composites by the extended finite element method. Engineering Fracture Mechanics, 77(17), 3373-3393. 18. Tanaka, S., Sannomaru, S., Imachi, M., Hagihara, S., Okazawa, S., & Okada, H. (2015). Analysis of dynamic stress concentration problems employing spline-based wavelet Galerkin method. Engineering Analysis with Boundary Elements, 58, 129-139. 30

19. Bhardwaj, G., Singh, I., & Mishra, B. (2015). Fatigue crack growth in functionally graded material using homogenized XIGA. Composite Structures, 134, 269-284. 20. Bhardwaj, G., Singh, I., & Mishra, B. (2015). Stochastic fatigue crack growth simulation of interfacial crack in bi-layered FGMs using XIGA. Computer Methods in Applied Mechanics and Engineering, 284, 186-229. 21. Fedeliński, P. (2004). Boundary element method in dynamic analysis of structures with cracks. Engineering Analysis with Boundary Elements, 28(9), 1135-1147. 22. Fedelinski, P., Aliabadi, M., & Rooke, D. (1996). The Laplace transform DBEM for mixed-mode dynamic crack analysis. Computers & structures, 59(6), 1021-1031. 23. García-Sánchez, F., Zhang, C., & Sáez, A. (2008). A two-dimensional time-domain boundary element method for dynamic crack problems in anisotropic solids. Engineering Fracture Mechanics, 75(6), 1412-1430. 24. Santana, E., & Portela, A. (2016). Dual boundary element analysis of fatigue crack growth, interaction and linkup. Engineering Analysis with Boundary Elements, 64, 176195. 25. Ghorashi, S. S., Mohammadi, S., & Sabbagh-Yazdi, S.-R. (2011). Orthotropic enriched element free Galerkin method for fracture analysis of composites. Engineering Fracture Mechanics, 78(9), 1906-1927. 26. Nguyen, N. T., Bui, T. Q., & Truong, T. T. (2017). Transient dynamic fracture analysis by an extended meshfree method with different crack-tip enrichments. Meccanica, 52(10), 2363-2390. 27. Nguyen, N. T., Bui, T. Q., Zhang, C., & Truong, T. T. (2014). Crack growth modeling in elastic solids by the extended meshfree Galerkin radial point interpolation method. Engineering analysis with boundary elements, 44, 87-97. 28. Wen, P., & Aliabadi, M. (2009). Evaluation of mixed-mode stress intensity factors by the mesh-free Galerkin method: static and dynamic. The Journal of Strain Analysis for Engineering Design, 44(4), 273-286. 29. Belytschko, T., Gu, L., & Lu, Y. (1994). Fracture and crack growth by element free Galerkin methods. Modelling and Simulation in Materials Science and Engineering, 2(3A), 519. 30. Peng, L., Tao, Y., Liang, N., Li, L., Qin, X., Zeng, Z., & Teng, X. (2017). Simulation of a crack in stiffened plates via a meshless formulation and FSDT. International Journal of Mechanical Sciences, 131, 880-893. 31. Pathak, H., Singh, A., & Singh, I. V. (2016). Three-dimensional quasi-static interfacial crack growth simulations in thermo-mechanical environment by coupled FE-EFG approach. Theoretical and Applied Fracture Mechanics, 86, 267-283. 32. Cai, Y., Han, L., Tian, L., & Zhang, L. (2016). Meshless method based on Shepard function and partition of unity for two-dimensional crack problems. Engineering Analysis with Boundary Elements, 65, 126-135. 33. Pant, M., Singh, I., & Mishra, B. (2011). Evaluation of mixed mode stress intensity factors for interface cracks using EFGM. Applied Mathematical Modelling, 35(7), 34433459. 34. Pathak, H., Singh, A., & Singh, I. V. (2014). Fatigue crack growth simulations of homogeneous and bi-material interfacial cracks using element free Galerkin method. Applied Mathematical Modelling, 38(13), 3093-3123.

31

35. Rao, B., & Rahman, S. (2000). An efficient meshless method for fracture analysis of cracks. Computational mechanics, 26(4), 398-408. 36. Rao, B., & Rahman, S. (2003). Mesh-free analysis of cracks in isotropic functionally graded materials. Engineering fracture mechanics, 70(1), 1-27. 37. Liu, G., & Gu, Y. (2001). A local radial point interpolation method (LRPIM) for free vibration analyses of 2-D solids. Journal of Sound and vibration, 246(1), 29-46. 38. Liu, G., & Gu, Y. (2003). A meshfree method: meshfree weak–strong (MWS) form method, for 2-D solids. Computational Mechanics, 33(1), 2-14. 39. Liu, G., Jiang, Y., Chen, L., Zhang, G., & Zhang, Y. (2011). A singular cell-based smoothed radial point interpolation method for fracture problems. Computers & Structures, 89(13), 1378-1396. 40. Farahani, B. V., Tavares, P. J., Moreira, P., & Belinha, J. (2017). Stress intensity factor calculation through thermoelastic stress analysis, finite element and RPIM meshless method. Engineering Fracture Mechanics. 41. Zhuang, X., Cai, Y., & Augarde, C. (2014). A meshless sub-region radial point interpolation method for accurate calculation of crack tip fields. Theoretical and Applied Fracture Mechanics, 69, 118-125. 42. Azevedo, J., Belinha, J., Dinis, L., & Jorge, R. N. (2015). Crack path prediction using the natural neighbour radial point interpolation method. Engineering Analysis with Boundary Elements, 59, 144-158. 43. Bui, T. Q., Nguyen, N. T., Van Lich, L., Nguyen, M. N., & Truong, T. T. (2017). Analysis of transient dynamic fracture parameters of cracked functionally graded composites by improved meshfree methods. Theoretical and Applied Fracture Mechanics. 44. Hematiyan, M., Khosravifard, A., & Liu, G. (2014). A background decomposition method for domain integration in weak-form meshfree methods. Computers & Structures, 142, 64-78. 45. Khoei, A., Yasbolaghi, R., & Biabanaki, S. (2015). A polygonal finite element method for modeling crack propagation with minimum remeshing. International Journal of Fracture, 194(2), 123-148. 46. Rybicki, E. F., & Kanninen, M. F. (1977). A finite element calculation of stress intensity factors by a modified crack closure integral. Engineering fracture mechanics, 9(4), 931938. 47. Zhang, H., & Ma, G. (2014). Fracture modeling of isotropic functionally graded materials by the numerical manifold method. Engineering Analysis with Boundary Elements, 38, 61-71. 48. Rice, J. R. (1968). A path independent integral and the approximate analysis of strain concentration by notches and cracks. 49. Song, S. H., & Paulino, G. H. (2006). Dynamic stress intensity factors for homogeneous and smoothly heterogeneous materials using the interaction integral method. International Journal of Solids and Structures, 43(16), 4830-4866. 50. Chen, J., Wu, L., & Du, S. (2000). A modified J integral for functionally graded materials. Mechanics research communications, 27(3), 301-306. 51. Yau, J., Wang, S., & Corten, H. (1980). A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity. Journal of Applied Mechanics, 47(2), 335-341.

32

52. Kim, J. H., & Paulino, G. H. (2002). Finite element evaluation of mixed mode stress intensity factors in functionally graded materials. International Journal for Numerical Methods in Engineering, 53(8), 1903-1935. 53. Liu, G.-R., & Gu, Y.-T. (2005). An introduction to meshfree methods and their programming: Springer Science & Business Media. 54. Khosravifard, A., Hematiyan, M., & Marin, L. (2011). Nonlinear transient heat conduction analysis of functionally graded materials in the presence of heat sources using an improved meshless radial point interpolation method. Applied Mathematical Modelling, 35(9), 4157-4174. 55. Freund, L. B. (1998). Dynamic fracture mechanics: Cambridge university press. 56. Kumar, S., Singh, I., Mishra, B., & Singh, A. (2016). New enrichments in XFEM to model dynamic crack response of 2-D elastic solids. International Journal of Impact Engineering, 87, 198-211. 57. Ravi-Chandar, K. Dynamic Fracture–Amsterdam, Netherlands: Elsevier.–2004.–254 p. 58. Chen, Y. (1975). Numerical computation of dynamic stress intensity factors by a Lagrangian finite-difference method (the HEMP code). Engineering Fracture Mechanics, 7(4), 653-660. 59. Wünsche, M., Zhang, C., Kuna, M., Hirose, S., Sladek, J., & Sladek, V. (2009). A hypersingular time‐ domain BEM for 2D dynamic crack analysis in anisotropic solids. International journal for numerical methods in engineering, 78(2), 127-150. 60. Albuquerque, E., Sollero, P., & Fedelinski, P. (2003). Dual reciprocity boundary element method in Laplace domain applied to anisotropic dynamic crack problems. Computers & structures, 81(17), 1703-1713.

33

Highlights  A novel meshfree approach for analysis of dynamic fracture in various media types is presented. 

Using the BDM integration technique, the computational burden is reduced significantly.



The number of nodes required for obtaining accurate results is less than similar numerical methods.



Effect of crack inclination angle with respect to principal material direction in anisotropic media is investigated.



Effect of gradation scheme and crack orientation in FG media is studied.

34