Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh

Engineering Fracture Mechanics xxx (2017) xxx–xxx

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh Tao Ni, Qi-zhi Zhu ⇑, Lun-Yang Zhao, Peng-Fei Li Key Laboratory of Ministry of Education for Geomechanics and Embankment Engineering, Hohai University, Nanjing 210098 China Jiangsu Research Center for Geotechnical Engineering Technology, Hohai University, Nanjing 210098 China

a r t i c l e

i n f o

Article history: Received 21 February 2017 Received in revised form 13 August 2017 Accepted 17 August 2017 Available online xxxx Keywords: Peridynamic method Irregular grids Local encryption technology Crack propagation Brittle solids

a b s t r a c t When subjected to increasing external loads, quasi brittle solids usually experience crack nucleation and propagation, leading to final failure of material. Numerical modeling and simulation of cracking process is an essentially important issue in mechanics and engineering science. This paper aims at applying the peridynamic method to simulate crack propagation using irregularly distributed material grids generated from finite element mesh. For post-treatment of physical fields such as strain and stress, a relevant visualization technique is developed. Various numerical tests are conducted for validation and application. Firstly, a square plate with a centrally located circular hole is considered for uniaxial tension tests performed with different non-uniform distributions of material grids. Peridynamic predictions are compared to finite element results and the influence of grid density, distribution and orientation on fracturing process is evaluated and the nonlocality of numerical results is demonstrated. Secondly, a pre-cracked plate with an off-centre circular hole is simulated under tension condition to assess the performance of locally refined grid. Finally, a square plate containing two parallel pre-existing cracks under biaxial tension is tested under different combinations of geometric parameters. The present work is in favour of developing a combined peridynamics/finite element numerical technique for efficient modelling and simulation of complex cracking problems. Ó 2017 Published by Elsevier Ltd.

1. Introduction Numerical modelling of damage and failure induced by fracturing is one of major concerns in solid mechanics and engineering science. In recent two decades, due to the growing demand of simulating crack initiation and propagation process in quasi brittle materials, various efforts have been made to develop suitable and robust numerical technique. According to their difference in basic assumptions, these approaches can be divided grossly into two categories. The first one is based on the classical theory of continuum mechanics, such as the extended finite element method (XFEM) [1,2] and the virtual multidimensional internal bonds (VMIB) method [3,4]. For the XFEM, additional degrees of freedom have to be introduced to formulate correctly displacement field near discontinuities like cracks and imperfect interfaces [6,7,5]. Also, crack prediction requires path tracing technique [8], which makes the XFEM less efficient in dealing with complex cracking problems like interactions between multiple cracks, crack branching and coalescence. The VMIB incorporates fracture criterion directly into the constitutive formulations of material and thus requires no external fracture criterion, making it very effective in ⇑ Corresponding author. E-mail address: [email protected] (Q.-z. Zhu). http://dx.doi.org/10.1016/j.engfracmech.2017.08.028 0013-7944/Ó 2017 Published by Elsevier Ltd.

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

2

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Nomenclature b c d t dc D E f G0 M s s0 S; V t Dt u € u x; x0 d dr

e g m q . r

n

body force density local stiffness modulus of a bond accounting for distance effect thickness of structure models under plane stress or plane strain assumption damping coefficient damage parameter describing bond breaking state Young’s modulus of solids pairwise force between two material particles connected by a bond macroscopic fracture energy fictitious diagonal density coefficient matrix for time integration local normal deformation (stretch) of bond critical stretch of bonds beyond which bond is broken surface and volume occupied by a material particle, respectively time instant time increment displacement vector acceleration vector of material particle location vector of material particles material horizon serving as a representative length reference horizon related to the biggest element of irregular mesh macroscopic strain tensor relative displacement vector of two material particles Poisson’s ratio of solids mass density of material particle characteristic function describing the connection status of bonds macroscopic stress tensor relative position vector of two material particles

simulating fracture growth [9]. However, the VMIB is still faced with some theoretical challenges in describing and simulating propagation, bifurcation and coalescence of cracks in brittle and quasi-brittle materials. The second category of numerical methods for cracking problems is based on the assumption of discrete theories, such as the particle discrete element method. This type of discrete approaches deal with failure process from diffused cracks to macro fracture without using complex constitutive equations, while they present an obvious shortcoming that there is no satisfactory way to calibrate the model’s parameters in terms of the inherent material properties [10–12]. Alternatively, a nonlocal theory of continuum called peridynamics was initiated by Silling [13] and later generalized by Silling et al. [14] and some other investigators. The peridynamic method (PDM) formulates mechanical problems in terms of integral equations rather than partial differential equations, making it more suitable for dealing with both continuous and discontinuous problems than the traditional finite element method (FEM). It should be emphasized that in the bond-based peridynamic model [13], even in the linear isotropic case, the effective Poisson’s ratio has been theoretically proved to be limited to be 1=3 in plain stress conditions and 1=4 in 3D or plain strain context [14,15]. In order to generalize the peridynamic theory by eliminating the limitation in Poisson’s ratio, the state-based peridynamics was proposed [14], and the bondbased peridynamics is proved to be a special case of the former. By involving some local failure criteria, the peridynamic model is able to efficiently predict the process of damage and fracture at different scales [16,17]. On the other hand, by making specific extensions, peridynamics has been applied successfully to model plastic deformation [18,19], dynamic brittle fracturing [17,20–22], quasi-static brittle fracturing problems [26,23–25], fatigue crack growth problems [27–29] and some other nonlinear problems [30–33]. By introducing brittle constitutive equations, it can be applied to simulate damage and failure processes [35,36,34,15,37]. It is worthy noticing that during peridynamic implementation material body is commonly discretized with a uniform distribution of material particles and using a constant horizon [17,38]. In fracture simulation with uniform material grids, the increase in accuracy would require a denser grid and a sufficiently large horizon, which makes computation less efficient for exact description of crack geometry. For improvement, several routines have been proposed in literatures to refine local particles in the area(s) of particular concern. For example, the adaptive grid refinement method [39] has been proved to be an effective way of predicting cracking path correctly. However, its computational cost becomes higher in view that particles’ redistribution should be done repeatedly and thus is not in favour of potential coupling with the finite element method. In this work, the peridynamic method was applied to simulate fracturing process in quasi brittle materials like concrete and rocks where direct use has been made of irregular finite element mesh to generate non uniform distribution of material particles. In that process, a relevant visualization technique is developed for post-treatment of physical fields such as strain and stress. The influence of the particle density and distribution orientation on the local and global behavior of the tested

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

3

structures is investigated to assess the non-locality of peridynamic predictions. In Section 4, a plate with a circular hole at its centre is used to validate the numerical implementation with irregular meshes in purely elastic and elastic brittle fracture problems. In order to further present the advantages of non uniform distribution of material particles by local refinement, an example of a pre-cracked plate with an off-centre circular hole under unilateral tensile loading is tested in Section 5.1, and the results are compared to some existing works. Finally, as an application example, fracturing process of a square plate containing two parallel pre-existing cracks is simulated under biaxial tension condition, in which numerical oscillation during the loading process is eliminated. Finally, some conclusions are addressed in Section 6. 2. Theoretical basis of the peridynamic method As a nonlocal numerical approach, the peridynamic theory initiated by Silling [13] describes mechanical behavior of continuum media by solving a space integral equation rather than partial differential equations in the finite element method. It assumes that two material points in a solid body interact once their distance is less than a given range d (usually called horizon), even if they are not in direct contact. In the classical peridynamics, interaction between material particles is assumed to be transmitted by means of a bond (see Fig. 1). At instant t, the pairwise constitutive force between material particles x and x0 belonging to the region R can be characterized by the following general formula [13]:

f ¼ f ðx; x0 ; uðx; tÞ; uðx0 ; t ÞÞ

ð1Þ

0

0

where uðx; t Þ and uðx ; tÞ denote the displacements of material particles x and x , respectively. According to the Newton’s second law, the motion equation of the material particle x at instant t is given as [17]:

qu€ ¼

Z

Hx

f ðx; x0 ; uðx; t Þ; uðx0 ; t ÞÞdV x0 þ bðx; t Þ

ð2Þ

€ is the acceleration vector, b is the body force density, that is, the external force applied to per where q is the mass density, u unit reference volume. Hx is the horizon associated with material particle x, which is usually a circle in 2D models and a sphere in 3D models. Mathematically Hx can be defined by [15]

Hx ¼ Hðx; dÞ ¼ fkx0  xk < d : x0 2 Rg

ð3Þ

The following relative position vector is usually introduced to simplify the motion equation and the constitutive force function

n ¼ x0  x

ð4Þ

Accordingly, a relative displacement vector at instant t is defined

g ¼ uðx0 ; tÞ  uðx; tÞ

ð5Þ

The constitutive force function is then rewritten as f ðg; nÞ. In order to describe the normal deformation (stretch) of a bond, the following scalar variable can be defined

Fig. 1. Illustration of the bond-based peridynamic model.

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

4

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

sðg; n; t Þ ¼

kg þ nk  knk knk

ð6Þ

where knk and kg þ nk are the norms of n and ðg þ nÞ, which represent the lengths of the bond in its initial and deformed state, respectively. In constitutive formulations, s and f are respectively analogous to the strain and the stress in the traditional continuum mechanics. By taking the bond-broken effect into consideration, the relationship between the constitutive force f and the stretch of the bond s is given as

f ðg; n; t Þ ¼ csðg; n; tÞ.ðg; n; tÞ

gþn kg þ nk

ð7Þ

where . is a scalar function recording the connection status of all bonds adhered to the particle x. The parameter c denotes the micro stiffness that describes the local resistance to the bond stretch. This parameter can be derived by relating the macroscopic strain energy obtained from the elasticity theory to that from the peridynamic theory. For simplicity, the parameter c is assumed to be uniform for all bonds inside the horizon. According to Huang et al. [15,37], the micromodulus c should gradually get weaker as the length of bond n ¼ knk increases, and its value should be zero when the length reaches the horizon size. In this view, the following relation in plain stress condition is used in this paper

cðn; dÞ ¼

8 < :

h

315E 8pdd3

1

 n  2 i2 d

0;

; when n 6 d

ð8Þ

when n > d

in which d is the thickness of structure model and E is the Young’s modulus of the solid under investigation, while the effective Poisson’s ratio is theoretically limited to be 1=3 in the plane stress condition [15]. In order to simulate brittle fracturing process, a local characteristic function is applied to identify the connection state of each bond

.ðg; n; tÞ ¼



s 6 s0 s > s0

1; 0;

ð9Þ

where s0 is the critical stretch of bonds. When s reaches s0 , the bond is broken and the two particles initially connected via this bond will no longer interact each other. s0 can be determined in terms of the macroscopic fracture energy G0 . By following the procedure developed by Ha and Bobaru (2010) [40] under plane stress condition, the fracture energy G0 can be formulated as:

Z G0 ¼ 2

d

0

Z

d

z

Z

cos1 ðz=nÞ

0

cs20 n ndhdndz 2

ð10Þ

By substituting Eq. (8) into Eq. (10) and after some simplifications, G0 finally reads

G0 

8Eds20 5p d

ð11Þ

It follows the expression of s0

s0 ¼

rffiffiffiffiffiffiffiffiffiffiffiffiffiffi 5pdG0 8Ed

ð12Þ

3. Numerical implementation using finite element mesh 3.1. Discretization process In order to solve the peridynamic equations of motion, the solid body is often discretized by material particles. More precisely, the spatial integration in peridynamics is transformed into a summation over a finite number of material grids. Structure models are generally geometrically irregular. Therefore, it is usually unrealistic to use regularly distributed material grids to match outer boundaries and inner discontinuities (pores, cracks, imperfect interfaces, etc.) satisfactorily. We propose here to use finite element mesh to generate irregular distribution of material particles for peridynamic analyses. As a potential, this process is in favour of peridynamics (PD)/finite element (FE) hybrid modelling and simulation. 3.1.1. Transformation from FE meshes to PD material points Take the body in Fig. 2 as an example. First, the material region is discretized by quadrilateral elements and then we take the barycentre ðxb ; yb Þ of each element as the position of a material particle. In 2D peridynamic context, given the area of an element Sb and its thickness d, the volume occupied by this particle is V b ¼ Sb d. Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

5

Fig. 2. A sketch of peridynamic material particles generated from finite element mesh.

3.1.2. Choice of the horizon dr If the material region is discretized by uniform grids spacing Dx; d for each material particle can be related to Dx, and d ¼ 3Dx is a common choice [17]. However, for the discretization based on finite element mesh, the spacing between adjacent particles is generally not a constant. Peridynamic computation becomes stable when the number of material particles belonging to every horizon is big enough. It is then reasonable to focus on the element that has the largest area among all the pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi elements. Denoting this area by Smax , we define a reference spacing Dxr ¼ Smax =p, to which the horizon dr is related. 3.2. Numerical implementation For numerical purpose, the integral Eq. (2) is cast into a summation

qu€ ðx; tÞ ¼

N Hx X f ðg; n; t ÞV cj þ bðx; t Þ

ð13Þ

j¼1

€ ðx; tÞ is the acceleration of point x at instant t; N Hx is the number of material particles that are located inside the horiwhere u zon of point x; bðx; tÞ is the density of the body force applied to point x at instant t. 3.2.1. Modification to boundary volume and description of damage state As illustrated in Fig. 3, the particles near the boundary actually belong partly to the horizon. In this sense, it is necessary to correct the volume fraction of these particles when performing the weighted summation. Inspired by Madenci and Oterkus [41], the following correction scheme for particle xj is used

Fig. 3. An example of volume correction.

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

6

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

V cj

8 > > <

kn þ gk P d   kn þ gk 6 d  Rj ¼   >   > : dþRj jjnþgjj V j ; when d  Rj 6 kn þ gk 6 d 2R 0; V j;

when when

ð14Þ

j

where Rj and V j respectively represent the reference radius and the volume of the particle xj . Given the area Sj issued from finite element mesh data, the following relations are defined:

Rj ¼

qffiffiffiffiffiffiffiffiffiffi Sj =p;

V j ¼ Sj d

ð15Þ

Moreover, when the local brittle constitutive law (9) is adopted, the damage parameter of particle x can be defined as

PN H x Dðx; t Þ ¼ 1 

.ðg; n; tÞV cj PNHx c j¼1 V j

j¼1

ð16Þ

3.2.2. Time integration Kilic and Madenci [42] introduced the adaptive dynamic relaxation method to simulate the quasi-static loading condition in the peridynamic framework, and this method has been detailed by Madenci and Oterkus [41] in their book. This technique t

is also suitable for the simulation of quasi-static mechanical behaviors of solids. By adding a damping coefficient dc and a € the velocity vector U _ fictitious diagonal density coefficient matrix M, the relationship between the acceleration vector U, and the force density vector F is set up

€ þ d MU _ ¼F MU c t

ð17Þ t

According to Madenci and Oterkus [41], dc and kii (the entries on the main diagonal of M) in 2D context take the following expressions t dc

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u t  T t t u U KU ¼ 2t  T Ut Ut

kii ¼

ð18Þ

pd2r dC Dt2

ð19Þ

4 Dxr

Recall that d is the thickness of the 2D structure model, and C takes the maximum value of formula (8). In Eq. (18), Kt is the ’local’ diagonal stiffness matrix, the entries on its main diagonal are given as

K tii ¼ 

F ti  F t1 i kii Dtu_ it1=2

ð20Þ

When the initial condition and the boundary condition of the problem are imposed, the velocity of particle x at instant t can be calculated with the following central-difference explicit integration formula [42]:



_ U

tþ1=2

¼

 t 2  dc Dt U_ t1=2 þ 2M1 Ft Dt t

2 þ dc Dt

ð21Þ

The displacement is then given as

Utþ1 ¼ Ut þ U_ tþ1=2 Dt

ð22Þ

At t ¼ 0, the time integration begins with: 1 0 _ 1=2 ¼ M F Dt U 2

ð23Þ

where Dt represents the time increment, which does not need to refer to the actual time. According to the previous researches [42], Dt ¼ 1 is a simple acceptable choice. 4. Elastic response and elastic brittle fracturing in a plate with a circular hole In this section, numerical tests about linearly elastic and brittle fracture problems are performed successively to validate the numerical algorithm and to assess the FE mesh-based discretization method. The geometric and mechanical parameters of the structure are provided in Fig. 4. The particles used to discretize the material region of the plate are generated from finite element meshes. Fig. 5 presents four different discretization schemes for specific purposes, in which different mesh Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

7

Fig. 4. Geometry and mechanical parameters of a plate under uniaxial tension.

Fig. 5. Different discretization schemes using finite element meshes.

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

8

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

generation modes are applied, and all the elements are of convex quadrilateral form. The peridynamic structure model under quasi-static tension is established with different micro constitutive laws. In each case, the peridynamic predictions are compared with the finite element results and/or analytic solutions. 4.1. Elastic problem Assume that the bond connecting two particles behaves linearly elastic. The external load r ¼ 30 MPa is applied to the particles at the outermost boundary layer, as shown in Fig. 5a. The reference horizon in this elastic example is chosen as dr ¼ 3Dxr . Denoting Lx as the length of the outer boundary layer, the force per unit volume applied to the material particle x is given by

Bforce ðxÞ ¼

rLx d Vx

ð24Þ

Fig. 6. Comparison of the displacement predictions between FEM and PDM with the micro-elastic model.

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

9

The irregular FE meshes presented in Fig. 5 are used to generate material particles for peridynamic simulations. Finite element analyses are carried out using the commercial software AbaqusÒ. In order to compare the PD predictions to the FE results, the displacements on the top and left edges of the plate are depicted, which are the values at the boundary nodes in finite element method (FEM) and outermost particles in peridynamics, respectively. It is seen from Fig. 6a and b that the results calculated by the FEM and the PDM are in a good agreement. To further investigate the influence of the parameter dr , a convergence study was conducted with the locally refined discretization scheme. A series of numerical results are shown in Fig. 7. It is seen that with the increase of dr , the PD solutions converge gradually to the FEM solution. It is known by comparing the solutions in different cases that the numerical result with dr ¼ 3Dxr is relatively close to those with dr ¼ 4Dxr and dr ¼ 5Dxr . In view of computational efficiency, the choice of dr ¼ 3Dxr is satisfactorily acceptable. A direct visualization of computational results can make the presentation more vivid and intuitive. During peridynamic simulations, material region is discretized with particles, and the grid information only contains the coordinate of these particles, bringing out some difficulty in making direct use of finite element software to visualize physical fields such as stress

Fig. 7. d-related convergence test using the locally refined discretization scheme (Fig. 5c).

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

10

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 8. Illustration of the winding node method.

and strain. By contrast, grid data obtained from finite element MeshTool to discretize material region can generate the information at nodes by using averaging method and applying linear extrapolation of peridynamic simulation data at material particles, that is to say, the barycentre of each finite element. As shown in Fig. 8, the calculation formula in 2D context is given in the form of surface-weighted averaging

PNn Si Di Dn ¼ Pi¼1 Nn i¼1 Si

ð25Þ

where Dn is the value at a FE node, Di are the values at material particles which are predicted by the peridynamic model. N n is the number of particles surrounding the FE node under consideration. In this way, the field data at FE nodes is achieved and on this basis, data visualization can be realized via commercial software. The displacement fields obtained by the PDM and the FEM are presented respectively in Figs. 9 and 10, which exhibit a very high similarity. It is seen that the proposed visualization technique can meet the requirement in accuracy. 4.2. Brittle fracturing tests We are now concerned with the fracturing process of the structure described in Fig. 4 by implementing the micro brittle constitutive model. The simulation procedure is the same as the test using the micro-elastic model. Two reference horizons dr ¼ 2Dxr and dr ¼ 3Dxr are chosen for comparison, as detailed in Table 1. The external load r ¼ 30 MPa is applied to the particles located at the boundary layer, and the body force allocated to each particle is calculated according to Eq. (24). In

Fig. 9. Contour of the displacement from peridynamic predictions.

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

11

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 10. Contour of the displacement from finite element predictions.

Table 1 Calculation schemes of brittle fracture examples. Mesh type

Reference horizon

Number of elements

Irregular sparse mesh

2Dxr 3Dxr 2Dxr 3Dxr 2Dxr 3Dxr 3Dxr

11488

Irregular dense mesh Locally refined mesh Oriented mesh

18412 14077 10371

addition, the parameter describing fracturing characteristics of the brittle materials G0 is chosen as 110 J=m2 for a kind of rock-like materials [7]. By inserting the value of G0 into the formula (12), the critical stretch of the bond s0 is obtained. As described above, when the stretch of a bond reaches its critical value, the bond is broken and a crack initiates and propagates between the particles connected by this broken bond. The deviation degree of crack path is estimated using a weighted indicator introduced in this paper. 4.2.1. Simulation with irregular and oriented grids The dependence of cracking path on finite element meshes is a perpetual numerical issue for fracture simulation. For the peridynamic method, this numerical issue has been discussed by Dipasquale et al. [43]. As a non-local theory of continuum, by nature, the dependence of cracking paths on grid distribution should be negligible. For validation, discretization schemes with irregular and locally aligned grids generated from specifically designed finite element meshes (see Fig. 5a–d) are used to predict fracturing process in the plate under tensile loading. When external loads are applied, some bonds connecting particles are deforming, approaching the critical value s0 and finally breaking up. The coalescence of broken bonds leads to crack propagation. The simulation results of cracking path using the irregular grids of Fig. 5a and b are shown in Fig. 11, and the peridynamic prediction of cracking path using the locally aligned grids of Fig. 5d is presented in Fig. 12. For dr ¼ 2Dxr , Fig. 11a and c present an asymmetry of initiation and propagation of two cracks; while for dr ¼ 3Dxr , the cracking paths are nearly symmetric with respect to the central horizontal line of the plate, as shown in Fig. 11b and d. In Fig. 12, the initiation and propagation of two cracks is perfectly symmetric, which conforms to the analytical solution by fracture mechanics. 4.2.2. Simulation with locally refined grids A locally refined FE mesh is applied to generate material particles for peridynamic analyses. The discretization scheme is shown in Fig. 5c. The potential region where cracks may appear is refined, while the rest of the structure is discretized with coarse mesh. The results predicted using the micro-elastic constitutive law with locally refined grid have been compared Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

12

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 11. Peridynamic prediction of crack initiation and propagation under different combinations of material grid and horizon.

Fig. 12. Crack initiation and propagation with oriented grids and dr ¼ 3Dxr .

with the finite element results in Section 4.1. When the micro-brittle constitutive law is concerned, two rectilinear cracks initiate and propagate along the central axis of the plate. Two levels of dr ; 2Dxr and 3Dxr , have been chosen in order to investigate the influence of dr on crack initiation and propagation. The simulation results are presented in Fig. 13.

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

13

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 13. Prediction of crack initiation and propagation using material grids generated from the locally refined mesh.

It is seen that the cracks initiate and propagate symmetrically in both the cases of dr ¼ 2Dxr and dr ¼ 3Dxr . Crack paths in both Fig. 13a and b are approaching to that predicted by fracture mechanics. It is intuitive to observe that the cracks in Fig. 13a are more localized than those in Fig. 13b, implying that the value of dr is a deterministic factor for the thickness of cracking zones. 4.2.3. Quantitative comparison between predicted crack paths We now focus on the deviation of cracking paths from the symmetry axis of the structure. To this end, a weighted indicator is introduced and detailed as follows. Theoretically, the fracturing path in the structure shown in Fig. 4 under quasistatic boundary tension should be rectilinear and located at y ¼ 0. As have been described above, when the stretch of the bond connecting the particles x and x0 reaches its critical value, the bond will be broken up with the nucleation of a new microcrack. Accordingly, with the initiation and propagation of cracks associated with particle x, the damage parameter Dðx; tÞ will increase. Let us collect all particles at which damage values are greater than zero. Suppose that the y-direction coordinate of the particle x is C yx , and the expected location of the crack is C ax ¼ 0, the position deviation of particle

10-4

1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0

1

2

3

4

5

6

Fig. 14. The values of error indicator for different combinations of material grids and the horizon size.

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

14

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

    x-related cracks is measured as C yx  C ax . The averaged position deviation of the particles in the set is used as an error indicator to describe the global degree of deviation from the expected value. We thus define the error indicator

Dc ¼

Nc   1 X  x x C y  C a  N c i¼1

ð26Þ

where N c is the number of the particles involved in the calculation. The values of the error indicator for the discretization schemes in Table 1 are calculated and depicted in Fig. 14. Based on the results shown in Figs. 11–14, the following remarks can be made: (1) The density of material particles affects the initiation and propagation of cracks. The relative error can be reduced by refining the particles in potential fracturing zones, as illustrated in Figs. 11 and 13. (2) The value of the reference horizon dr affects greatly the initiation and propagation of cracks. As the value of dr increases, simulation results becomes more accurate and thicker crack bands appear. A sufficiently big value of dr can ensure that the material horizon contains a sufficient number of particles. It allows to eliminate the deviation of initiation positions and propagation path. This has been confirmed by comparing Fig. 11a and b, or by comparing Fig. 11c and d. (3) Analogously to the above observations, the dependence of cracking path on the orientation of grid distribution encountered by Dipasquale et al. [43] can be avoided only if dr is big enough. (4) The local refinement on finite element mesh is an effective way to improve the prediction of crack propagation. This method can reduce computational cost and maintain numerical accuracy, as shown in Figs. 13 and 14. This point will be further discussed in the following section.

5. Application to crack predictions In this section, the peridynamic method together with the proposed discretization process are applied to simulate two representative cracking problems: a pre-cracked plate with an off-centre circular hole under unilateral tensile loading as

Fig. 15. Geometry and finite element discretization of a pre-cracked plate with an off-centre circular hole.

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

15

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx Table 2 Model schemes of pre-cracked plate with an off-centre circular hole. Models

h (mm)

Horizon dr

Grid number

A B C

5 10 15

3Dxr 3Dxr 3Dxr

4607 6845 5566

Fig. 16. Peridynamic cracking prediction for different pre-cracking positions.

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

16

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 17. The cracking paths reported by Tabiei and Wu [45] and Dipasquales and Zaccariotto [39].

well as a square plate with two parallel pre-fabricated fissures under biaxial tensile loading. During the tests, different simulation schemes are taken into account for specific purpose. 5.1. Unilateral tension test upon a pre-cracked plate with an off-center circular hole This test has been simulated by various investigators. We mention, for example, Rashid [44] using the arbitrary local mesh replacement method on the basis of finite element method, Tabiei and Wu [45] using the Dyna3D method and Dipasquale et al. [39] using 2D peridynamics with adaptive grid refinement. The geometry parameters and discretization details are presented in Fig. 15. The structure is subjected to a uniform normal traction at the top and bottom edges. The purpose of this test is to assess the influence of the off-centre circular hole on the cracking path. The location of the initial crack is described by the parameter h in Fig. 15a. Three tests are performed respectively with h ¼ 0:005 m; h ¼ 0:01 m and h ¼ 0:015 m. The other parameters used in the simulation are collected in Table 2. The mechanical parameters of the material are identical to Dipasquale et al. (2014) [39]: E ¼ 71:4 GPa, q ¼ 2700 kg=m3 ; m ¼ 1=3 and G0 ¼ 1000 J=m2 . The traction applied to the top edge is such that r ¼ 22 MPa in model A, r ¼ 26 MPa in model B and r ¼ 28 MPa in model C. The material particles at the bottom edge are constrained, as shown in Fig. 15b. The simulation results in the three cases are presented in Fig. 16. For each distance h, there are three damage diagrams corresponding to different loading steps, which indicate the process of crack propagation in the plate. It is clearly seen that smaller the value of h is, more greatly the crack deviates from the hole.

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

17

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 18. Geometry and mesh of a rock specimen containing two pre-existing cracks.

Table 3 Model schemes of square plate containing two parallel pre-existing cracks. Model scheme

a (mm)

b (mm)

h (mm)

Crack angle ( )

Grid number

A B C A1 A2 A3 A4 A5 A6

10 10 10 10 10 10 10 10 10

10 20 10 10 10 10 10 10 10

0 0 5 0 0 0 0 0 0

45 45 45 0 15 30 60 75 90

4607 6845 5566 3418 4956 5050 5110 5287 4879

Fig. 19. Prediction of the crack path for the models A, B and C.

For comparison purpose, the numerical results reported by Tabiei and Wu [45] and Dipasquale et al. [39] are presented in Fig. 17a and b, respectively. The present peridynamic predictions in Fig. 16 shows more similarities with Fig. 17a than Fig. 17b. The satisfactory prediction is attributed to the discretization scheme using irregular finite element mesh. The profile of the hole’s geometry is described more accurately in the present model than Dipasquale et al. [39].

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

18

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

8 7 6 5 4 3 2 1 0

0

0.5

1

1.5

2

2.5

3 10-4

8 7 6 5 4 3 2 1 0

0

0.5

1

1.5

2

2.5

3 10

-4

Fig. 20. Stress–strain curves predicted using the models A, B and C.

5.2. Fracture in rock containing two parallel cracks under bi-directional tensile loading In this application, the fracturing path in a square rock plate containing two parallel pre-existing cracks is simulated under biaxial tensile loading condition. The material region is discretized with a locally refined mesh from which are generated material particles for peridynamic analyses. The geometric parameters and discretization scheme are presented in Fig. 18. In order to study the effect of the location of the pre-fabricated fissures, different combinations of the geometric parameters are considered and collected in Table 3.

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

19

Fig. 21. Prediction of the crack path using the models A A1 ; A2 ; A3 ; A4 ; A5 and A6 .

5.2.1. Mechanical parameters and loading condition The mechanical parameters of the material are such that E ¼ 30 GPa, m ¼ 1=3, q ¼ 2700 kg=m3 and G0 ¼ 30 J=m2 for the limestone used in Ha et al. [46]. In order to obtain complete stress–strain curves, the numerical tests are displacement-controlled, and the distribution of the boundary particles is shown in Fig. 18b. Two loading rates, 106 mm=s and 107 mm=s, are considered.

5.2.2. Simulation results According to Lehoucq and Silling [47], denoting the displacement and the resultant of the force density of the particles at boundary layer as Dtx and F tx , respectively, the axial strain ea and the axial stress ra on macro scale can be given by:

ea ¼

2Dtx L

ð27Þ

b 1 X FtV i Wd i¼1 i

N

ra ¼

ð28Þ

where L and W is the length and width of the specimen in the axial direction, N b is the index number of particles located at the boundary layer, d is the thickness of the structure. Note that Eq. (28) is consistent with the definition of the Piola–Kirchhoff stress. The numerical tests have been conducted with the combinations of parameters listed in Table 3. The results for the schemes A; B and C performed at the loading rate of 107 mm per step are shown in Fig. 19. The curves of ea versus ra are plotted in Fig. 20. It is seen from comparison between Fig. 19a–c that the influence of the pre-existing cracks’ locations on crack propagation is significant: the parameter h affects the type of rock bridge failure, while b determines whether the rock bridge breaks through earlier than lateral rock.

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

20

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

16 14 12 10 8 6 4 2 0

0

1

2

3

4

5 10-4

16 14 12 10 8 6 4 2 0

0

1

2

3

4

5 10-4

Fig. 22. Stress–strain curves obtained using the models A1 ; A2 ; A3 ; A4 ; A5 and A6 .

The simulation results for the schemes A; A1 ; A2 ; A3 ; A4 ; A5 and A6 are shown in Fig. 21, and the curves of ea and ra are plotted and compared in Fig. 22. As the inclined angle increases, the pre-existing cracks propagate similarly except for a ¼ 90 . It can be seen from Fig. 19 and 21 that all the pre-cracking patterns nearly present the same failure configuration: a failure region appears between the two initial cracks while the other two cracks develop towards the left and right sides of the plate. The predicted crack propagation paths show a high degree of consistency with the results reported in Tian [48] and Zhang and Chen [3]. It is concluded that the peridynamic method equipped with irregular grids from finite element mesh allows predicting correctly crack propagation process in quasi brittle materials. 5.2.3. A possible way of eliminating numerical oscillation The curves shown in Figs. 20 and 22 present a slight difference between subfigures ðaÞ and ðbÞ that are obtained at different loading rates. The curves with lower loading rate are smoother than those with higher loading rate. Also, the peak

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

21

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

15 14 13 12 11 10 9 8 7 6 5

0

10

20

30

40

50

60

70

80

90

80

90

Fig. 23. Peak stress values for the models A; A1 ; A2 ; A3 ; A4 ; A5 and A6 .

31.5 31 30.5 30 29.5 29 28.5 28 27.5 27

0

10

20

30

40

50

60

70

Fig. 24. Slope values for the models A; A1 ; A2 ; A3 ; A4 , A5 and A6 .

values of the latter are greater than the former. These phenomena are caused by the numerical oscillation in the peridynamic method. To be more explicit, the slope and peak values of the curves in Figs. 20 and 22 are plotted and compared in Figs. 23 and 24, respectively. The values of slope approximately equals to the elastic modulus of the material (E ¼ 30 GPa), which demonstrates the validity of numerical simulation and the post-treatment on the global stress and strain. The uniaxial tensile strength of the specimen increases progressively within the range from 6 MPa to 15 MPa as the inclined angle increases, as shown in Figs. 22 and 23. It is worthy noticing that even with the adaptive dynamic relaxation technique, numerical oscillation still occurs in the peridynamic method. In fact, we showed that when the loading rate is low enough, the numerical problem of oscillation can be limited to an acceptable extent.

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

22

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

Fig. 25. Contour of stress components: comparisons between PDM and FEM.

6. Conclusions and discussions In this paper, the peridynamic method was applied to simulate cracking problems in brittle solids. Material grids generated directly from finite element mesh have been applied for peridynamic modelling of quasi-static crack propagation. An efficient visualization process of numerical results was also developed. With 2D structure examples, the proposed implementation technique was evaluated and a first stage of applications was conducted. We confirmed that the dependence

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

23

of numerical results on the density and orientation of material grids can be removed when material horizon is sufficiently large. Such results conform to the nonlocal nature of the peridynamic theory. In the bond-based peridynamics, although the expression of the Piola–Kirchhoff stress can be given directly, the expression of Cauchy’s stress can not be derived theoretically, mainly due to the discrete nature of peridynamic theory. However, the displacement field can be produced just like in finite element framework. As mentioned in Section 4.1, the values of displacement on each node of finite element mesh can be calculated by applying the averaging method of winding nodes. In turn, the displacement field u helps to generate the field of Cauchy stress r by making use of finite element mesh. As an example, Fig. 25 presents the comparison of stress contours over the plate between the peridynamic method and the finite element method. The result issued from peridynamics simulation agrees quite well with finite element analyses. In addition, there are also no significant differences between the two methods for both the stress components rx and ry (see Fig. 25a–d), while sxy shows a slight difference. In summary, the method we developed to generate strain and stress field is effective, which provides a good foundation for further hybrid PD/FEM analyses. Acknowledgements This work has been jointly supported by the National Natural Science Foundation of China (Grant No. 51679068), the Fundamental Research Funds for the Central Universities (Grant No. 2016B20214) and the 111 Project (Grant No. B13024). Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.engfracmech.2017.08.028. References [1] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Meth Eng 1999;45:601–20. [2] Moes N, Dolbow J. A finite element method for crack growth without remeshing. Int J Numer Meth Eng 1999;46:131–50. B.T.. [3] Zhang ZN, Chen YQ. Numerical simulation for fracture propagation of multi-cracked rock materials using virtual multidimensional internal bonds. Chinese J Geotech Eng 2008;30(10):1490–5. [4] Zhang ZN, Ge XR. Micromechanical consideration of tensile crack behavior based on virtual internal bond in contrast to cohesive stress. Theor Appl Fract Mec 2005;43(3):342–59. [5] Zhu QZ, Gu ST, Yvonnet J, Shao JF, He QC. Three-dimensional numerical modelling by xfem of spring-layer imperfect curved interfaces with applications to linearly elastic composite materials. Int J Numer Meth Eng 2011;88(4):307–28. [6] Yvonnet J, Quang HL, He QC. An xfem/level set approach to modelling surface/interface effects and to computing the size-dependent effective properties of nanocomposites. Comput Mech 2008;42(1):119–31. [7] Cox JV. An extended finite element method with analytical enrichment for cohesive crack modeling. Int J Numer Meth Eng 2009;78(1):48–83. [8] Dumstorff P, Meschke G. Crack propagation criteria in the framework of x-fem-based structural analyses. Int J Numer Anal Met 2007;31(2):239–59. [9] Zhang ZN, Ge XR. A new quasi-continuum constitutive model for crack growth in an isotropic solid. Eur J Mech A – Solid 2005;24(2):243–52. [10] Wu SC, Xu XL. A study of three intrinsic problems of the classic discrete element method using flat-joint model. Rock Mech Rock Eng 2016;49 (5):1813–30. [11] Wang YN, Tonon F. Modeling lac du bonnet granite using a discrete element model. Int J Rock Mech Min 2009;46(7):1124–35. [12] Hazzard JF, Young RP. Simulating acoustic emissions in bonded-particle models of rock. Int J Rock Mech Min 2000;37(5):867–72. [13] Silling SA. Reformulation of elasticity theory for discontinuities and long-range forces. J Mech Phys Solids 2000;48(1):175–209. [14] Silling SA, Epton M, Weckner O, Xu J, Askari E. Peridynamic states and constitutive modeling. J Elasticity 2007;88(2):151–84. [15] Huang D, Lu GD, Qiao PZ. An improved peridynamic approach for quasi-static elastic deformation and brittle fracture analysis. Int J Mech Sci 2015;9495:111–22. [16] Askari E, Bobaru F, Lehoucq RB, Parks ML, Silling SA, Weckner O. Peridynamics for multiscale materials modeling. J Phys Conf Ser 2008;125(1):12–78. [17] Silling SA, Askari E. A meshfree method based on the peridynamic model of solid mechanics. Comput Struct 2005;83(17-18):1526–35. [18] Foster JT, Silling SA, Chen WW. Viscoplasticity using peridynamics. Int J Numer Meth Eng 2010;81(10):1242–58. [19] Madenci E, Oterkus S. Ordinary state-based peridynamics for plastic deformation according to von mises yield criteria with isotropic hardening. J Mech Phys Solids 2016;86:192–219. [20] Gu X, Zhang Q, Huang D, Yv Y. Wave dispersion analysis and simulation method for concrete shpb test in peridynamics. Eng Fract Mech 2016;160 (160):124–37. [21] Lee J, Liu W, Hong JW. Impact fracture analysis enhanced by contact of peridynamic and finite element formulations. Int J Impact Eng 2016;87 (10):108–19. [22] Shojaei A, Mudric T, Zaccariotto M, Galvanetto U. A coupled meshless finite point/peridynamic method for 2d dynamic fracture analysis. Int J Mech Sci 2016;119:419–31. [23] Zhou X, Wang Y, Xu X. Numerical simulation of initiation, propagation and coalescence of cracks using the non-ordinary state-based peridynamics. Int J Fract 2016;203(1-2):1–22. [24] Wang Y, Zhou X, Xu X. Numerical simulation of propagation and coalescence of flaws in rock materials under compressive loads using the extended non-ordinary state-based peridynamics. Eng Fract Mech 2016;163:248–73. [25] Lee J, Ha YD, Hong JW. Crack coalescence morphology in rock-like material under compression. Int J Fract 2017:1–26. [26] Zhou XP, Bi J, Qian QH. Numerical simulation of crack growth and coalescence in rock-like materials containing multiple pre-existing flaws. Rock Mech Rock Eng 2015;48(3):1097–114. [27] Oterkus E, Guven I, Madenci E. Fatigue failure model with peridynamic theory. In: Ther Thermomech Phenomena Electr Syst. p. 1–6. [28] Zaccariotto M, Galvanetto U. Modeling of fatigue crack propagation with a peridynamics approach. In: Esmc 2012 Euromech Solid Mechanics Conference. [29] Zaccariotto M, Luongo F, Sarego G, Dipasquale D, Galvanetto U. Fatigue crack propagation with peridynamics: a sensitivity study of paris law parameters. In: 4:th CEAS 2013 Conference. [30] Liu W, Hong J. Discretized peridynamics for brittle and ductile solids. Int J Numer Meth Eng 2012;89(8):1028–46. [31] Silling SA, Weckner O, Askari E, Bobaru F. Crack nucleation in a peridynamic solid. Int J Fract 2010;162(1-2):219–27.

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028

24

T. Ni et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx

[32] Breitenfeld MS, Geubelle PH, Weckner O, Silling SA. Non-ordinary state-based peridynamic analysis of stationary crack problems. Comput Method Appl M 2014;272(11):233–50. [33] Mossaiby F, Shojaei A, Zaccariotto M, Galvanetto U. Opencl implementation of a high performance 3d peridynamic model on graphics accelerators. Comput Math Appl 2017. doi: http://dx.doi.org/10.1016/j.camwa.2017.06.045. in press. [34] Zhou XP, Gu XB, Wang YT. Numerical simulations of propagation, bifurcation and coalescence of cracks in rocks. Int J Rock Mech Min 2015;80:241–54. [35] Gerstle WH, Sau N, Sakhavand N. On peridynamic computational simulation of concrete structures. Aci Special Publication 2009:245–64. [36] Kilic B, Madenci E. Structural stability and failure analysis using peridynamic theory. Int J Nonlin Mech 2009;44(8):845–54. [37] Huang D, Lu GD, Wang CW, Qiao PZ. An extended peridynamic approach for deformation and fracture analysis. Eng Fract Mech 2015;141:196–211. [38] Silling SA, Lehoucq RB. Peridynamic theory of solid mechanics. Adv Appl Mech 2010;44(10):73–168. [39] Dipasquale D, Zaccariotto M, Galvanetto U. Crack propagation with adaptive grid refinement in 2d peridynamics. Int J Fract 2014;190(1-2):1–22. [40] Ha YD, Bobaru F. Studies of dynamic crack propagation and crack branching with peridynamics. Int J Fract 2010;162(1-2):229–44. [41] Madenci E, Oterkus E. Peridynamic theory and its applications. New York: Springer; 2014. [42] Kilic B, Madenci E. An adaptive dynamic relaxation method for quasi-static simulations using the peridynamic theory. Theor Appl Fract Mec 2010;53 (3):194–204. [43] Dipasquale D, Sarego G, Zaccariotto M, Galvanetto U. Dependence of crack paths on the orientation of regular 2d peridynamic grids. Eng Fract Mech 2016;160:248–63. [44] Rashid MM. The arbitrary local mesh replacement method: an alternative to remeshing for crack propagation analysis. Comput Method Appl M 1998;154(1):133–50. [45] Ala T, Jin W. Development of the dyna3d simulation code with automated fracture procedure for brick elements. Int J Numer Meth Eng 2003;57 (14):1979–2006. [46] Ha YD, Lee J, Hong JW. Fracturing patterns of rock-like materials in compression captured with peridynamics. Eng Fract Mech 2015;144:176–93. [47] Lehoucq RB, Silling SA. Force flux and the peridynamic stress tensor. J Mech Phys Solids 2008;56(4):1566–77. [48] Tian R. Finite-cover-based element-free method for continuous and discontinuous deformation analysis with applications in geotechnical engineering (in chinese). Ph.D. thesis. Dalian University of Technology; 2001.

Please cite this article in press as: Ni T et al. Peridynamic simulation of fracture in quasi brittle solids using irregular finite element mesh. Engng Fract Mech (2017), http://dx.doi.org/10.1016/j.engfracmech.2017.08.028