Dynamic crack propagation in plates weakened by inclined cracks: an investigation based on peridynamics

Abstract

Peridynamics is a theory in solid mechanics that uses integral equations instead of partial differential equations as governing equations. It can be applied to fracture problems in contrast to the approach of fracture mechanics. In this paper by using peridynamics, the crack path for inclined crack under dynamic loading were investigated. The peridynamics solution for this problem represents the main features of dynamic crack propagation such as crack bifurcation. The problem is solved for various angles and different stress values. In addition, the influence of geometry on inclined crack growth is studied. The results are compared with molecular dynamic solutions that seem to show reasonable agreement in branching position and time.

Keywords

peridynamics / inclined crack / dynamic fracture / crack branching

Introduction

The classical solid mechanics is based on partial differential equations. These equations are not valid on discontinuities. On the other hand, many problems of solid mechanics involve discontinuities in geometry or displacement field. So the equations of classical solid mechanics can not be applied on the whole of such problems. Fracture Mechanics involving criteria such as strain energy density criterion [¹,²] have been used to treat problems with initial defects. These criteria must be used in solution procedure. The situation in dynamic fracture is more difficult. The time dependence of a dynamic fracture problem results in equations that are more complicated than the same static problems [³–⁷]. The challenges in simulating dynamic fracture problems such as crack propagation speed and the time of crack bifurcation are many [⁸–¹³].

Recently, Silling [¹⁴] as a nonlocal expansion of classical continuum mechanics proposes a new theory called peridynamics (PD). PD is a theory in solid mechanics that uses integral equations instead of partial differential equations as governing equations [¹⁵]. Because the same equations are used on either continuous or discontinuous points, there is no need to study the points which are on sides of the crack surface, separately. In this theory particles of a continuum interact with each other across a finite distance, similar to molecular dynamics (MD).It has been shown that when horizon size in peridynamic formulation approach zero, molecular dynamics can be recovered from peridynamics as a special case. Therefore, peridynamics has the potential for a multiscale method [¹⁶–¹⁸]. Because both of discretized peridynamics and molecular dynamics have the same computational structures, more recently the peridynamics is implemented in LAMMPS, a molecular dynamic code [¹⁹,²⁰].

In PD theory, internal forces are formulated in terms of interaction between the material particles. As interactions between the particles stop, cracks propagate. Hence, a failure has been implicated at the molecular scale. Therefore, PD does not require any special criterion for crack growth [²¹,²²] and integral equations in PD theory are directly useable in cracked surfaces. Kilic and Madenci [²³] used this theory to predict crack growth patterns in quenched glass plates containing single and multiple pre-existing cracks. Using peridynamic model, Ha and Bobaru [²⁴,²⁵] were able to capture main characteristics of dynamic brittle fracture for a straight crack such as crack path shapes and crack propagation speeds.

In this study by using a bond-based peridynamic model, the crack path for inclined crack is investigated under dynamic loading. The PD solution for this problem represents the main features of dynamic crack propagation like crack propagation speed and the time of crack branching. The problem is solved for various angles and different stress values. In addition, the influence of geometry on inclined crack growth is studied. The results are compared with molecular dynamic solutions and show good agreement in branching position and time.

Basic equations of peridynamics

The equation of motion in PD is [¹⁴]:

(1)\rho \ddot{u}(x,t)={\displaystyle \underset{H_x}{\overset{}{\int }}f(u(x^{\prime },t)-u(x,t),x^{\prime }-\mathbf{x})d{V}_{x^{\prime }}+b(\mathbf{x},t)},

wherefis defined as a pairwise force vector that the particle x^{\prime } exerts on the particle x, and H_x is a neighborhood of the particle x. This neighborhood is a spherical (in 2D a circular) region around xwith radius d. The d is called horizon ofx. In classical Peridynamics horizon required to be constant. For no uniform refinement and adaptive modelling the dual horizon peridynamics approach can resolve the issue of ghost forces and spurious waves present in the classical peridynamics formulation [²⁶,²⁷].

The interaction between particles x and x^{\prime } is called bond. There is no interaction between particles that are not in horizon of x. Therefore, with the notation\xi =x^{\prime }-x as the relative position vector in the reference configuration and \eta =u(x^{\prime })-u(x) as the current relative displacement, we have:

(2)\left|\xi \right|>\delta \stackrel{}{\Rightarrow }f(\eta ,\xi )=\mathrm{0}\forall \xi ,\eta ,

where f(h,x) determines the behavior of material. For a linear isotropic material, the pairwise force function has the following form:

(3)f(\eta ,\xi )=cs{\displaystyle \frac{\partial \left|\xi +\eta \right|}{\partial \eta }},

wherec is called the micromodulus parameter and s is the bond relative elongation:

(4)s={\displaystyle \frac{\left|\xi +\eta \right|-\left|\xi \right|}{\left|\xi \right|}}.

Bond relative elongation is the other definition of one-dimensional strain. We can see from Eq. (4) that a bond at its equilibrium length has no stretch and s=\mathrm{0}, and when a bond stretches twice its equilibrium length has s=\mathrm{1}. For a 2D problem with plane stress conditions, c is proposed as [²⁴]:

(5)c_{\mathrm{2}\mathrm{D}}={\displaystyle \frac{\mathrm{6}E}{\pi {\delta }^{\mathrm{3}}(\mathrm{1}-\nu )}},

with E being the Young’s modulus and v the Poisson’s ratio.

Damage

When bonds stretch beyond a given limit, they break. When a bond breaks, it is failed permanently. In addition, during the course of a simulation new bonds are never created. Now in this part we discuss one criterion for bond break which is called the critical stretch criterion. With the entrance of failure into the peridynamics, the pairwise force will be [²⁸]:

(6)f(\eta ,\xi )={\displaystyle \frac{\xi +\eta }{\left|\xi +\eta \right|}}cs\mu (t,\left|\xi \right|),

where\mu (t,\left|\xi \right|) is a history-dependent scalar function:

(7)

The s_{\mathrm{0}} is the critical relative elongation and in 2D is derived as [²⁴]:

(8)s_{\mathrm{0}}=\sqrt{{\displaystyle \frac{\mathrm{4}\pi G_{\mathrm{0}}}{\mathrm{9}E\delta }}},

whereG_{\mathrm{0}} is the fracture energy. The Eq. (6) states, if the relative elongation of a bond exceeds s_{\mathrm{0}}, the bond breaks forever. With this notion of bond failure in peridynamics, the local damage at a material point is defined as:

(9)\phi (x,t)=\mathrm{1}-{\displaystyle \frac{{{\displaystyle \int }}_{H_x}^{}\mu (x,t,\xi )\mathrm{d}{V}_{\xi }}{{{\displaystyle \int }}_{H_x}^{}\mathrm{d}{V}_{\xi }}}.

This is the ratio of the broken bonds to the total initial bonds.

Numerical scheme

For the purpose of numerical solutions, a meshfree method is proposed in [²⁸]. In this method, the body is discretized into nodes. Using the mid-point integration scheme, the equation of motion (Eq. (1)) is written as:

(10)\rho {\ddot{u}}_i^n={\displaystyle \sum _{j}f(u_j^n-u_i^n,x_j-x_i)V_j+b_i^n},

whereu_i^n is the displacement of node i at nth timestep, j is a particle within horizon of i, and V_j=\Delta x^{\mathrm{3}} is the volume of node j. In 2D the volume of a node is its area. The acceleration term on the left-hand side of this equation is replaced with an explicit central-difference formula:

(11){\ddot{u}}_i^n={\displaystyle \frac{u_i^{n+\mathrm{1}}-\mathrm{2}{u}_i^n+u_i^{n-\mathrm{1}}}{\Delta t^{\mathrm{2}}}}.

In this paper, instead of central-difference formula the velocity-Verlet algorithm [²⁹] is used which is faster and numerically more stable. This algorithm is implemented in three steps:

(12){\dot{u}}_{n+\mathrm{1}/\mathrm{2}}={\dot{u}}_n+{\displaystyle \frac{\Delta t}{\mathrm{2}}}\ddot{u}{}_n,

(13)u_{n+\mathrm{1}}=u_n+\Delta t{\dot{u}}_{n+\mathrm{1}/\mathrm{2}},

(14){\dot{u}}_{n+\mathrm{1}}={\dot{u}}_{n+\mathrm{1}/\mathrm{2}}+{\displaystyle \frac{\Delta t}{\mathrm{2}}}\ddot{u}{}_{n+\mathrm{1}},

where\Delta t is the timestep size. The stable timestep size is:

(15)

The crack propagation speed at t_i is computed by

(16)v_i={\displaystyle \frac{\Vert x_i-x_{i-\mathrm{1}}\Vert }{t_i-t_{i-\mathrm{1}}}},

where x_i and x_{i-\mathrm{1}} denote the crack tip positions at the times t_i andt_{i-\mathrm{1}}, respectively.

Numerical examples

To validate the current PD code accuracy, at first we consider a rectangular plate with a long pre-notch (Fig. 1). The material properties are density r=2235 kg/m³, Young’s modulus E=65 GPa, Poisson’s ratio n=0.2 and energy release rate G₀=204 J/m². The top and bottom boundaries are subjected to a uniform tensile stress s=12 MPa. The stress is applied as a step function. For numerical solution, we discretize the plate into nodes, each node has a grid spacing Dx=0.0005 m. The horizon size is d=0.002 m and the stable time step size is Dt=50 ns (as given by Eq. (15)). To introduce the pre-crack, all the bonds that cross the pre-crack line are broken.

The crack path at 46 ms for this research, Ref [²⁴], and peridynamic in LAMMPS results is shown Fig. 2. In LAMMPS, at the start of simulation, no bond is considered between the upper and lower sides of the pre-crack. So the damage of particles that are near the initial crack, is computed zero. As can be seen clearly from Fig. 2, the results are close to each other. They predict the branching event and a considerable growth of crack before branching. For more comparison, Table 1 represents the branching initiation time and position predicted by these three approaches.

Inclined crack growth

Square plate

Consider a square plate with 0.05 m length. A pre-crack of length 0.01 m and angle q with respect to x axis is at the center of the plate (Fig. 3). The plate is subjected to the tensile stresses s_x and s_y in x and y directions, respectively. The material properties are the same as used in the former problem.

Apply the stresses σ_x=σ_y=17 MPa and detect the crack path up to 22 μs. Figure 4 shows the crack path at this time for angles θ=0°, 15°, 30°, 45° which are compared with the LAMMPS results. The branching events do not occur in these cases. In addition, the crack paths are not exactly straight except for θ=0°.

The problem is considered again by applying the stresses σ_x=σ_y=23 MPa. The crack paths at 22 μs are shown in Fig. 5. As observed from Fig. 5, the pre-crack at first is bifurcated into two cracks and then each crack propagates. There is no prediction of crack growth before bifurcation.

Rectangular plate

In this section, the influence of plate geometry on the crack path is studied. Consider a rectangular plate with 0.1 m length and 0.04 m width. An inclined pre-crack of length 0.01 m is at the center of the plate in the direction of the rectangle diagonal (Fig. 6). The plate is subjected to the tensile stresses s_x and s_y in x and y directions, respectively. The material properties are the same as used in the former problems.

Apply the stresses s_x=s_y=15 MPa to the plate. Figure 7 shows the crack path at 30 ms. As observed from Fig. 7, after a growth parallel the horizontally axis, the crack path deviates from its direction and tends to the top and bottom boundaries.

The problem is solved again by applying the stresses s_x=s_y=18 MPa. The crack path at 30 ms is shown in Fig. 8. As Fig. 8 shown, after the deviation of the crack path, the branching is occurred. Table 2 represents the time and position of the crack deviation and branching. In addition, Figure 9 shows the crack propagation speed for the lower right-hand side path. The crack propagation speed is computed by determining the crack-tip position at any time. The crack-tip is the most lower-right node, which its damage is larger than 0.3. With regard to Table 2 and Fig. 9, at crack deviation and crack branching points, the crack propagation speed decreases.

Conclusions

In this paper by using a bond-based peridynamic model, some characteristics of dynamic fracture in inclined cracks are shown. Furthermore, the effect of different factors such as geometry of crack and the boundary condition on the crack propagation were investigated.

In the case of a square plate, it is noteworthy that simulation predicted there is no crack growth before possible bifurcation. However, in the case of a rectangular plate, there is a crack growth and then possibly crack bifurcates.

The most advantage of the peridynamics, is that it does not use any additional criterion in the simulation and solution procedure of a problem. Also, because the same equations are used on either continuous or discontinuous points, there is no need to study the points which are on sides of crack surface, separately. The numerical solution of peridynamics is done by the mesh free method. This method do not need to generate mesh, therefore this model has no difficulty to simulate the complex geometries.

References

[1] Rabczuk T. Computational Methods for Fracture in Brittle and Quasi-Brittle Solids: State-of-the-Art Review and Future Perspectives. ISRN Applied Mathematics, 2013: 1–38

[2] Zehnder A. Fracture Mechanics. Springer Netherlands, 2012

[3] Rabczuk T, Bordas S, Zi G. On three-dimensional modelling of crack growth using partition of unity methods. Computers & Structures, 2010, 88(23–24): 1391–1411

[4] Areias P, Rabczuk T, Camanho P P. Initially rigid cohesive laws and fracture based on edge rotations. Computational Mechanics, 2013, 52(4): 931–947

[5] Amiri F, Anitescu C, Arroyo M, Bordas S, Rabczuk T. XLME interpolants, a seamless bridge between XFEM and enriched meshless methods. Computational Mechanics, 2014, 53(1): 45–57

[6] Rabczuk T, Gracie R, Song J H, Belytschko T. Immersed particle method for fluid-structure interaction. International Journal for Numerical Methods in Engineering, 2010, (81): 48–71

[7] Ravi-Chandar. Dyanamic Fracture. Elsevier, 2004

[8] Areias P M A, Rabczuk T, Camanho P P. Finite strain fracture of 2D problems with injected anisotropic softening elements. Theoretical and Applied Fracture Mechanics, 2014, 72: 50–63

[9] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A simple and robust three-dimensional cracking-particle method without enrichment. Computer Methods in Applied Mechanics and Engineering, 2010, 199(37–40): 2437–2455

[10] Song J, Wang H, Belytschko T. A comparative study on finite element method for dynamic fracture. Computational Mechanics, 2008, 42(2): 239–250

[11] Rabczuk T, Belytschko T. Cracking particles: a simplified meshfree method for arbitrary evolving cracks. International Journal for Numerical Methods in Engineering, 2004, 61(13): 2316–2343

[12] Rabczuk T, Bordas S, Zi G. A three-dimensional meshfree method for continuous multiple-crack initiation, propagation and junction in statics and dynamics. Computational Mechanics, 2007, 40(3): 273–495

[13] Areias P, Rabczuk T. Finite strain fracture of plates and shells with configurational forces and edge rotation. International Journal for Numerical Methods in Engineering, 2013, 94(12): 1099–1122

[14] Silling S. Reformulation of elasticity theory for discontinuities and long-rang forces. Journal of the Mechanics and Physics of Solids, 2000, 48(1): 175–209

[15] Silling S, Lehoucq R. Peridynamic theory of solid mechanics. Advances in Applied Mechanics, 2010, 44(10): 73–168

[16] Talebi H, Silani M, Bordas S P A, Kerfriden P, Rabczuk T. Molecular dynamics/XFEM coupling by a three-dimensional extended bridging domain with applications to dynamic brittle fracture. International Journal for Multiscale Computational Engineering, 2013, 11(6): 527–541

[17] Budarapu P, Gracie R, Bordas S, Rabczuk T. An adaptive multiscale method for quasi-static crack growth. Computational Mechanics, 2014, 53(6): 1129–1148

[18] Rahman R, Foster J T, Haque A. A multiscale modeling scheme based on peridynamic theory. International Journal of Multiscale Computational Engineering, 2014, 12(3): 223–248

[19] Parks M, Lehoucq R, Plimpton S, Silling S. Implementing peridynamics within a molecular dynamics code. Computer Physics Communications, 2008, 179: 777–783

[20] Parks M, Seleson P, Plimpton S, Lehoucq R, Silling S. Peridynamics with LAMMPS: A User Guide v0.2 Beta, Sandia Report, 2010

[21] Silling S, Weckner O, Askari E, Bobaru F. Crack nucleation in a peridynamic solid. International Journal of Fracture, 2010, 162(1–2): 219–227

[22] Madenci E, Oterkus E. Peridynamics theory and its applications. New York: Springer-Verlag, 2014

[23] Kilic B, Madenci E. Peridiction of crack paths in a quenched glass plate by using peridynamic theory. International Journal of Fracture, 2009, 156(2): 165–177

[24] Ha Y D, Bobaru F. Studies of dynamic crack propagation and crack branching with Peridynamics. International Journal of Fracture, 2010, 162(1–2): 229–244

[25] Ha Y D, Bobaru F. Characteristics of dynamic brittle fracture captured with Peridynamics. Engng Fract Mech, 2011 (78): 1156–1168

[26] Ren H, Zhuang X, Rabczuk T. Dual-horizon peridynamics: a stable solution to varying horizons. Computer Methods in Applied Mechanics and Engineering, 2017, 318: 762–782

[27] Ren H, Zhuang X, Cai Y, Rabczuk T. Dual-Horizon Peridynamics. International Journal for Numerical Methods in Engineering, 2016, 108(12): 1451–1476

[28] Silling S, Askari E. A mesh free method based on the peridynamic model of solid Mechanics. Comput Struct, 2005, 83(17): 1526–1535

[29] Sticker B, Schachinger E. Basic Concepts in Computational Physics. Springer, 2014

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature