Modeling nonlinear elastic solid with correlated lattice bond cell for dynamic fracture simulation

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347 www.elsevier.com/locate/cma

Modeling nonlinear elastic solid with correlated lattice bond cell for dynamic fracture simulation Zhennan Zhang ∗ , Yaxiong Chen School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai, 200240, China Received 9 February 2014; received in revised form 12 June 2014; accepted 25 June 2014 Available online 3 July 2014

Abstract The correlated lattice bond cell is a kind of lattice model developed for modeling nonlinear elastic solid. Unlike the usual lattice model, the present one is composed of unit cells, which can adopt any geometry with any number of bonds. The unit cells are compiled together via the particles along their boundaries to form a discrete system. The energy of unit cell is characterized by the modified Stillinger–Weber (SW) potential, in which the initial value of each bond angle in the reference configuration, not the ‘ideal’ tetrahedral angle, is taken as the reference value of this bond angle in the current configuration. This makes the modified SW potential applicable to other materials than silicon. Because the modified SW potential can simultaneously account for the bond stretch and rotation effect, the present lattice model can represent variable Poisson ratios. The bond parameters are calibrated as the function of material constants (Young’s modulus, Poisson ratio), the volume and the total bond number of a unit cell based on the ideal unit cell conception. Enriched with a more elaborate fracture mechanism than the two-body potential through the modified SW potential, the present lattice model can capture most characters of dynamic fracture, providing an efficient nonlinear elasticity modeling method for dynamic fracture simulation. For the present method can directly simulate the dynamic fracture initiation, arrest, propagation and branching in the finite deformation cases without any separate criterion, its perspective should be inspiring. c 2014 Elsevier B.V. All rights reserved. ⃝

Keywords: Nonlinear elastic solid; Dynamic fracture simulation; Lattice model; Modified Stillinger–Weber potential; Unit cell

1. Introduction The lattice model, pioneered by Hrennikoff [1], represents the continuum with a discrete lattice system. Through this discrete system, the 3D fracture of continuum is interpreted into the 1D linkage rupture. The rupture mechanism of the 1D linkage is much easier to capture than the 3D continuum. Hence, the lattice method can avoid many difficulties that the continuum mechanics method encounters in fracture simulation and presents many advantages to continuum mechanics method. For the lattice method, there are two key problems that have not been well addressed. One is the Poisson ratio problem and another is the parameter calibration. The linkage between particles in a lattice model could be a spring (bond) or a beam. When the linkage is a spring, which has no bending stiffness, the lattice model is essentially identical to the particle methods [2–4]. The Poisson ratio ∗ Corresponding author. Tel.: +86 21 34204346; fax: +86 21 34204346.

E-mail addresses: [email protected], [email protected] (Z. Zhang). http://dx.doi.org/10.1016/j.cma.2014.06.036 c 2014 Elsevier B.V. All rights reserved. 0045-7825/⃝

326

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

they present is fixed, namely 1/3 for 2D case and 1/4 for 3D case. To solve the fixed Poisson ratio problem, an effective approach is to replace the spring with a beam, e.g. [5,6]. Because the beam can account for the shear deformation effect between particles due to its bending stiffness, the variable Poisson ratios can be presented. However, the additional degree of freedom to the particle displacement, i.e. the particle rotation, is introduced into the lattice beam model. Thus, the beam conception could not be extended to the interatomic interaction because the degree of freedom of ‘atom’ is only characterized by the displacement, with no rotation. In order to account for the bond rotation effect in the context of the particle displacement, many approaches have been proposed. In the 3D distinct lattice spring model [7], the shear and the normal deformation between particles are considered separately. However, the shear deformation directly calculated by the particle displacements fails to exclude the rigid body rotation. To remedy this drawback, a local strain-based method has to be adopted in this method. Wang et al. [8,9] creatively introduced an angular spring that restricts the variation of bond angle in the lattice model. By this measure, the Lennard-Jones type of potential can be used to describe the axial and angular spring. For detailed review on this method, refer to [10]. But, in the axial–angular spring model, how to define the energy potential of the two kinds of springs is still an important and open issue. In the parameter calibration, the usual lattice model relies on the specific topology of the lattice structure. Through the equivalence of strain energy or stiffness coefficients, the bond parameters are related to the material constants. However, such parameter calibration method is case-dependent. When the lattice structure is changed, the corresponding bond parameters will have to be re-calibrated. For the 3D cases, such calibration method becomes highly complicated. The discretized virtual internal bond method (discretized-VIB) [11] is a newly-developed lattice method based on the VIB theory [12]. Different from the usual lattice bond methods, the discretized-VIB is composed of unit cells. Each cell has finite number of bonds and can take arbitrary geometry. The bond parameter calibration is based on an ideal unit cell rather than a specific one. So, a general micro–macro parameter relationship is derived. This general micro–macro parameter relationship can make the discretized-VIB reach a very high precision when applied to the specific unit cell, e.g. the triangular unit cell, as reported in [11]. However, the critical limitation of discretized-VIB is that the Poisson ratio it presents is fixed since it employs the two-body potential to describe the bond. From the more physical standpoint, the fixed Poisson ratio suggests that the two-body potential is inadequate to characterize the micro energy mechanism of material. It must ignore some other basic micro energy mechanisms. The Stillinger–Weber (SW) potential [13,14] suggests that the bond energy is not only related to the bond stretch, but also related to the bond angle subtended by this given bond and other bonds. It presents more elaborate micro energy mechanism than the two-body potential. So, the SW potential should be a feasible choice to characterize the unit cell energy of the discretized-VIB. However, the original SW potential was specially developed for the silicon materials. It takes the ‘ideal’ tetrahedral angle as the reference value of each bond angle in the current configuration. Hence, the original SW potential could not be directly applied to describe the unit cell energy. Recently, Zhang et al. [15] have modified the SW potential and embedded the modified SW potential into the constitutive relation of continuum. In this modified SW potential, the initial value of each bond angle in the reference configuration, not the ‘ideal’ tetrahedral angle, is taken as the reference value of this bond angle in the current configuration. By this modification, not only is the micro energy mechanism presented in the original SW potential preserved, the SW potential also can be applied to other materials than the silicon. Hence, in the present paper, we will use the modified SW potential to characterize the unit cell energy and then use it to simulate the dynamic fracture propagation. Because the bonds in a unit cell are correlated together through the modified SW potential, the present method is here termed as the ‘correlated lattice bond cell’ method. Different micro energy mechanism reveals different micro fracture mechanism. Through the modified SW potential, more micro fracture mechanisms are incorporated into the lattice model and the fixed Poisson ratio problem of lattice model is addressed. Another merit of using the modified SW potential to describe the micro bond is that the response of dynamic fracture is more likely to be captured since the hyperelasticity governs the dynamic fracture according to [16,17]. 2. Idea of correlated lattice bond method 2.1. Modeling method Material is usually composed of grains on the meso scale, e.g. rock, shown in Fig. 1a. Inspired by such mesostructure, the discretized-VIB [11] models the continuum with a discrete bond system (Fig. 1b). The discrete system is

327

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

a

b

c

Micro bond

Fig. 1. Correlated lattice bond cell modeling method for materials (a) mesostructure of material (SEM image of marble); (b) discrete bond system; (c) the correlated lattice unit cell.

Fig. 2. The conjugate bond vectors subtended at a given particle i in a discrete system with N particles. (The total particle number is N . The bond number of ri j is N (N − 1) and its conjugate bond rik number is N − 2 with ri j ̸= r ji .)

composed of lattice bond cells (Fig. 1c). Each unit cell has finite number of particles, thus has finite number of bonds. It can take any geometry. Such structural character makes the unit cell quite suitable for modeling the grains of materials on the meso-scale. Many unit cells are compiled together by the particles along their boundaries to form a discrete system. 2.2. Characterization of unit bond cell energy In the discretized-VIB, each bond is described by a two-body potential. The two-body potential ignores the contribution of bond rotation, leading to the fixed Poisson ratio on the macro scale. Compared with the two-body potential, the SW potential can simultaneously account for the bond stretch and bond rotation effect. The original SW potential has the form       Φ3 ri j , rik , θ jik (1) Φ= Φ2 ri j + i< j

i̸= j j
in which Φ2 , Φ3 stands for the two- and three-body interaction, respectively; ri j the bond length of bond vector ri j ; θ jik the bond angle subtended by the bond vector ri j , rik at vertex i. According to Eq. (1), the energy of the bond ri j is not only related to its length, but also related to the bond angles subtended by itself and its conjugate bonds rik . For a system with N particles, the total number of rik is N − 2, shown in Fig. 2. The specific expressions of Φ2 and Φ3 are respectively    p  q      ℓ0 ℓ0 1 − exp Φ2 ri j = A B ri j ri j ri j /ℓ0 − r0∗ (2)  2       1 γ γ Φ3 ri j , rik , θ jik = λ cos θ jik + exp exp 3 ri j /ℓ0 − r0∗ rik /ℓ0 − r0∗

328

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

˜ 3 /∂θIJ · Fig. 3. The normalized resistance force versus the bond angle deviation from its reference state. (The normalized force is defined as F = ∂ Φ ˜ 3 = λ (cos θIJ − cos θIJ0 )2 . The ITA denotes the ideal tetrahedral angle, ITA = arccos (−1/3).) (2λ)−1 = − sin θIJ (cos θIJ − cos θIJ0 ) with Φ

in which r0∗ is the cutoff radii; the 1/3 being the minus cosine of the ‘ideal’ tetrahedral angle. The parameters A, λ, here, implicitly contain the energy scaling parameter ε of the original two- and three-body interactions of [13]. From Eq. (2), it is noted that the original SW potential takes the ‘ideal’ tetrahedral angle, whose cosine is −1/3, as the reference value of each bond angle in the current configuration. This makes the original SW potential specially applicable to the silicon, not to the other materials. Moreover, governed by Φ3 of Eq. (2), the randomized particle configuration is not an equilibrium system and will be forced into a structured configuration. Hence, the original SW potential cannot be directly used to describe the unit cell. Recently, Zhang et al. [15] have modified the original SW potential and incorporated the modified SW potential into the constitutive relation of continuum. The topology information in the original SW potential has been wiped off in the modified SW potential, which makes the modified SW potential applicable to other materials than silicon, e.g. the quasi-brittle materials. The two- and three-body interaction of the modified SW potential read    −1  −q −p Φ2 (˜r I ) = A B r˜I − r˜I exp r˜I − r0∗   (3) −1  −1  Φ3 (˜r I , r˜ J , θIJ ) = λ (cos θIJ − cos θIJ0 )2 exp γ r˜I − r0∗ + γ r˜ J − r0∗ in which r˜I , r˜ J are the normalized bond length, r˜I = ri j /ℓ0 , r˜ J = rik /ℓ0 ; θIJ0 stands for the bond angle subtended by the bond ri j and rik in the reference configuration whereas θIJ its value in the current configuration. From Eq. (3) it is noted that the initial value of the bond angle in the reference configuration, not the ‘ideal’ tetrahedral angle, is taken as the reference value of this bond angle in the current configuration. This indicates that the reference configuration is the most stable one. The resistance to the variation of bond angle by Φ3 of Eq. (3) is shown in Fig. 3. From Fig. 3 it is seen that the resistance force generally increases with increasing the deviation of bond angle from its initial state. But the magnitude of resistance force is relevant to the initial value of bond angle. Besides the nonlinear modified SW potential, Zhang et al. [15] also suggested the linear modified SW potential, whose two- and three-body interaction read 1 2 Aℓ (˜r I − 1)2 2 0 1 Φ3 = λ (θIJ − θIJ0 )2 . 2

Φ2 =

(4)

Zhang et al. [15] have verified that both the Hookean matrixes derived from the linear and nonlinear SW potential are consistent with that of the linear elastic continuum. This provides the physical foundation for applying the modified SW potential to characterize the lattice bond cell energy.

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

329

3. Force–displacement constitutive relation of lattice bond cell Let u denote the particle displacement vector of a unit cell, f the corresponding force vector and K the stiffness matrix. By Eq. (1), the total energy of a unit cell is identified as W =

N −1) N −1) N −2  1 (N 1 (N Φ2 (˜r I ) + Φ3 (˜r I , r˜ J , θIJ ) 2 2 I I J

(5)

in which N denotes the total particle number in a unit cell with ri j ̸= r ji . By Eq. (5), the particle force vector f is derived as fi =

 N −1) N −1) N −2   ∂W ∂Φ2 ∂ r˜I ∂Φ3 ∂ r˜I 1 (N 1 (N ∂Φ3 ∂ r˜ J ∂Φ3 ∂θIJ = · + · + · + · ∂u i 2 ∂ r˜I ∂u i 2 ∂ r˜I ∂u i ∂ r˜ J ∂u i ∂θIJ ∂u i I I J

(6)

and the stiffness matrix of a unit cell is   N −1) ∂2W ∂ 2 Φ2 ∂ r˜I ∂ r˜I 1 (N ∂Φ2 ∂ 2r˜I Ki j = · = · + · ∂u i ∂u j 2 ∂ r˜I ∂u i ∂u j ∂ r˜I2 ∂u j ∂u i I  N −1) N −2  ∂ 2 Φ3 ∂ r˜I ∂ r˜I 1 (N ∂Φ3 ∂ 2r˜I · + · + · 2 ∂ r˜I ∂u i ∂u j ∂ r˜I2 ∂u j ∂u i I J  ∂ 2 Φ3 ∂ r˜ J ∂ r˜ J ∂Φ3 ∂ 2r˜ J ∂ 2 Φ3 ∂θIJ ∂θIJ ∂Φ3 ∂ 2 θIJ · + · · + · + · + · . (7) ∂ r˜ J ∂u i ∂u j ∂θIJ ∂u i ∂u j ∂ r˜ J2 ∂u j ∂u i ∂θIJ2 ∂u j ∂u i   According to Eq. (3), there are totally five parameters, i.e. A, λ, r0∗ , γ , B , in the force–displacement constitutive relation. The key point is that what values these parameters should take so that the discrete system is equivalent to a continuum in mechanics. 4. Determination of bond parameters   To determine the bond parameters A, λ, r0∗ , γ , B , we firstly define an ideal unit cell [11] in which the bond number is large enough and the bonds are uniformly distributed in space. A uniform strain is imposed on this ideal unit cell and the bonds have the identical initial length. The strain energy density of the ideal unit cell is   N −1) N −1) N −2  1 (N 1 (N 1 · Φ2 (˜r I ) + Φ3 (˜r I , r˜ J , θIJ ) (8) W = V 2 2 I I J in which V is the volume of the ideal unit cell. According to the character of an ideal unit cell, we can use a continuous function to describe the bond distribution. Let D denote the bond distribution density function, for an isotropic unit cell, D = D (α, ϕ) = constant with α, ϕ being the coordinates in a spherical coordinate system. The relationship between the bond distribution density and the total bond number can be derived as  2π  π nob (9) D · sin α · dαdϕ = nob ⇒ D = 4π 0 0 in which nob means the total bond number. Let ξ , η respectively denote the orientation vector of a pair of conjugate bonds r I , r J connected at one particle, ξ = (sin α cos ϕ, sin α sin ϕ, cos α)T and η = (sin α¯ cos ϕ, ¯ sin α¯ sin ϕ, ¯ cos α) ¯ T in the spherical coordinates. With the bond distribution density D, the unit cell energy due to Φ2 is N (N −1)  I

Φ2 (˜r I ) =

 0

2π

 0

π

Φ2 (˜r I ) · D I sin α · dαdϕ

(10)

330

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

and that due to Φ3 is N (N −1) N −2   I

Φ3 (˜r I , r˜ J , θIJ ) =

2π

0

I

J



N (N −1)  2π 

π





2π

π



0

0

¯ ϕ¯ Φ3 (˜r I , r˜ J , θIJ ) · D J sin α¯ · d αd

0

 Φ3 (˜r I , r˜ J , θIJ ) · D J sin α¯ · d αd ¯ ϕ¯ · D I sin α · dαdϕ

= 0

π



(11)

0

in which D I is the bond r˜I distribution density, D I = N (N − 1) /4π and D J is the bond r˜ J distribution density, D J = (N − 2) /4π, according to Eq. (9). By Eqs. (8), (10) and (11), the energy density W of the ideal unit cell is written as 1 1 ⟨Φ2 · D I ⟩2 + ⟨Φ3 · D I D J ⟩3 2V 2V  2π  π in which the operator ⟨· · · ⟩2 = 0 0 (· · · ) · sin α · dαdϕ and

(12)

W =



2π

π



2π



π



⟨· · · ⟩3 = 0

0

0

¯ ϕdαdϕ. ¯ (· · · ) · sin α¯ · sin α · d αd

0

In the small deformation cases, according to Cauchy–Born rule [18], the normalized bond length and the bond angle can be calculated as r˜I = ξ T εξ + 1 r˜ J = η T εη + 1  θIJ εi j 



 = arccos  

(13)

 r˜ TI r˜ J  r˜ TI r˜ I r˜ TJ r˜ J

  ≈ θIJ (0) + Θi j εi j

in which  ∂θIJ (ε)  = − Θi j = ∂εi j ε=0

  1 T T · 2ξ η − ξ η · ξ ξ − ξ η · η η i j i j i j  2 1 − ξTη

(14)

with ε being the strain tensor of the ideal unit cell. By Eq. (13), we have ∂ 2r˜I /∂εi j ∂εkl = 0, ∂ 2r˜ J /∂εi j ∂εkl = 0, ∂ 2 θIJ /∂εi j ∂εkl = 0. So, by Eq. (7), the tangent modulus of the unit cell is derived as   2W 2Φ ∂ D ∂ I 2 K idis = · ξi ξ j ξk ξl jkl = ∂εi j ∂εkl 2V ∂ r˜I2 2   2 2 D I D J ∂ Φ3 ∂ Φ3 ∂ 2 Φ3 + · ξi ξ j ξk ξl + · ηi η j ηk ηl + · Θi j Θkl . (15) 2V ∂ r˜I2 ∂ r˜ J2 ∂θIJ2 3 Under the undeformed state (i.e. ε = 0), Eq. (15) becomes       D I ∂ 2 Φ2  D I D J ∂ 2 Φ3  dis  · ξi ξ j ξk ξl + K i jkl  =   ε=0 2V 2V ∂ r˜I2  ∂θIJ2  ε=0

2

 · Θi j Θkl ε=0

since ∂ 2 Φ3 /∂ r˜I2 = ∂ 2 Φ3 /∂ r˜ J2 = 0 according to Eqs. (3) and (4) when ε = 0.

(16) 3

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

331

4.1. Parameter calibration for nonlinear SW potential Substituting Eq. (3) into Eq. (16), we have        D I Ac1 exp h −1 D I D J λ exp 2γ h −1   0 0  = K idis ξi ξ j ξk ξl 2 + sin2 θIJ0 · Θi j Θkl jkl  ε=0 3 2V V in which the coefficient c1 is     −3 −2 −3 2 c1 = B p 2 + p + ph −2 + 2h − q + q + qh + 2h 0 0 0 0     −2 B = h −2 sinceΦ2′ (1) = 0. 0 + q / h0 + p

(17)

(18)

   ⌢ ⌢  The sub-Hookean matrixes obtained by integrating the term ξi ξ j ξk ξl 2 and sin2 θIJ0 · Θ i j Θ kl 3 are respectively 3 Ω dis L

4π = 15

1  1  0  0 0

1 3 1 0 0 0

1 1 3 0 0 0

0 0 0 1 0 0

0 0 0 0 1 0

0 0  0 , 0  0 1

4 −2  128π 2  −2 = 225  0 0 0 

Ω dis R

−2 4 −2 0 0 0

−2 0 −2 0 4 0 0 3 0 0 0 0

So, the full Hookean matrix by Eq. (17) is calculated as     −1 D I Ac1 exp h −1 D D λ exp 2γ h I J 0 0 Ω dis = · Ω dis · Ω dis L + R . 2V V

0 0 0 0 3 0

 0 0  0 . 0  0 3

(19)

(20)

The nonzero components of the symmetrical Hookean matrix of the continuum are c c c Ω11 = Ω22 = Ω33 =

E (1 − ν) (1 + ν) (1 − 2ν)

c c c Ω12 = Ω13 = Ω23 =

Eν (1 + ν) (1 − 2ν)

c c c Ω44 = Ω55 = Ω66 =

E . 2 (1 + ν)

(21)

By equaling Ω dis to Ω c , the bond parameters are calibrated as A=

2V N (N − 1) c1 exp



h −1 0

·

3E , (1 − 2ν)

45E (1 − 4ν)  · λ= . −1 16 (1 + ν) (1 − 2ν) N (N − 1) (N − 2) exp 2γ h 0 V

(22)

4.2. Parameter calibration for linear SW potential Substituting Eq. (4) into Eq. (16), we have   K idis jkl 

ε=0

=

  D I ℓ20 A  DI D J λ  ξi ξ j ξk ξl 2 + Θi j Θkl 3 . 2V 2V

(23)

332

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

Fig. 4. Displacement of a pair of conjugate bonds. 2 Here, it is interesting to note that the integrand corresponding to the three-body interaction   in Eq. (17) is sin θIJ0 · Θi j Θkl while is Θi j Θkl in Eq. (23). The Hookean matrixes by integrating the term Θi j Θkl 3 is

Ω dis R

4 −2 32π 2  −2 =  45  0 0 0

−2 4 −2 0 0 0

−2 −2 4 0 0 0

0 0 0 3 0 0

0 0 0 0 3 0

0 0  0 . 0 0 3

(24)

According to Eq. (24) and Ω dis L of Eq. (19), we have the full Hookean matrix Ω dis =

D I ℓ20 A dis D I D J λ dis ΩL + ΩR . 2V 2V

(25)

By equaling Eq. (25) to Eq. (21), we can obtain the micro–macro parameter relationship A=

3E 2V · , N (N − 1) ℓ20 (1 − 2ν)

λ=

9E (1 − 4ν) V · . N (N − 1) (N − 2) 2 (1 + ν) (1 − 2ν)

(26)

5. Implementation Generally, when using the present method, the first step is to ‘mesh’ the interesting zone with the specific unit cells, e.g. the tetrahedron cells, and then compute the particle force vector and the stiffness matrix of the unit cell according to Eqs. (6) and (7). Compile the force vector and stiffness matrix of unit cells together into the global ones. Then this system can be solved as the usual finite element method. 5.1. Computation of bond length and its derivatives ¯ = Take a pair of conjugate vectors as example (Fig. 4). The bond vectors in the reference configuration are X ¯ = X3 − X1 with X1 , X2 , X3 being the coordinate vectors of the particles 1, 2, 3. Their corresponding X2 − X1 and Y bond vectors in the current configuration are ¯ + u¯ I x¯ = X ¯ + u¯ J y¯ = Y in which u¯ I = u2 − u1 and u¯ J = u3 − u1 with u1 , u2 , u3 being the displacement vectors of the particles 1, 2, 3.

(27)

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

333

For the bond x¯ with local particle index 1, 2, its particle displacement vector can be written as u¯ = [u¯ 1 , u¯ 2 , u¯ 3 , u¯ 4 , u¯ 5 , u¯ 6 ]T . The normalized bond length and its derivatives with respect to the node displacement vector are 1 x¯m x¯m ℓ0 ∂ r˜I 1 = · ti ∂ u¯ i ℓ0 ℓ 1 ℓ0 ∂ r˜I ∂ r˜I ∂ 2r˜I = pi j − · · ∂ u¯ i ∂ u¯ j ℓ0 ℓ ℓ ∂ u¯ i ∂ u¯ j √ √ ¯ ℓ = x¯ T x¯ ; the matrix p and the vector t are respectively ¯ T X, in which ℓ0 = X       pi j = (δ4i − δ1i ) δ4 j − δ1 j + (δ5i − δ2i ) δ5 j − δ2 j + (δ6i − δ3i ) δ6 j − δ3 j r˜I =

ti = x¯1 (δ4i − δ1i ) + x¯2 (δ5i − δ2i ) + x¯3 (δ6i − δ3i )

(28)

(29)

with δi j as the Kronecker delta. 5.2. Computation of bond angle and its derivatives ¯ Y ¯ in the reference and current configuration can be respectively calculated The cosine of the bond angle between X, as    ¯TX Y ¯TY ¯ ¯ T Y/ ¯ X cos θIJ0 = X    cos θIJ = x¯ T y¯ / x¯ T x¯ y¯ T y¯ .

(30)

To compute the derivative of the bond √ angle with respect to node displacement vector, let the scalar Q 11 = x¯m x¯m , Q 22 = y¯m y¯m , Q 12 = x¯m y¯m , T = Q 12 / Q 11 Q 22 , and then the bond angle is calculated as θIJ = arccos (T ). The first derivative of bond angle with respect to displacement vector is identified as   ∂θIJ (u) 1 ∂ Q 11 ∂ Q 22 −1 ∂ Q 12 T T θIJ′ (u) = · · =√ · − − (31) √ ∂u i 2Q 11 ∂u i 2Q 22 ∂u i Q 11 Q 22 ∂u i 1 − T2 in which  T ∂ Q 12 /∂u i = −x¯1 − y¯1 , −x¯2 − y¯2 , −x¯3 − y¯3 , y¯1 , y¯2 , y¯3 , x¯1 , x¯2 , x¯3  T (32) ∂ Q 11 /∂u i = 2 −x¯1 , −x¯2 , −x¯3 , x¯1 , x¯2 , x¯3 , 0, 0, 0  T ∂ Q 22 /∂u i = 2 − y¯1 , − y¯2 , − y¯3 , 0, 0, 0, y¯1 , y¯2 , y¯3 . According to Eq. (7), the stiffness matrix of unit cell involves with the second derivatives of the bond angle with respect to node displacement. However, Eq. (31) indicates that second derivative of the bond angle would be extremely complicated. To simplify this problem, we expand the bond angle into Taylor series and take the linear term as its approximate value. That is θIJ (u) ≈ θIJ (0) + Θi u i

(33)

in which Θ = θIJ′ (0). Thus, the variation of bond angle can be directly calculated as θIJ (u) − θIJ (0) ≈ Θi u i in the small deformation cases. In the stiffness matrix, the usage of the linearized bond angle does not affect the convergence nature of computation according to the principle of the Newton–Raphson method. 5.3. Constitutive relation with determined parameters By Eqs. (3), (6), (7) and (22), we have the force vector of the unit cell f i=

 N −1) N −1) N −2   ∂ r˜I 1 (N ∂ r˜I ∂ r˜ J ∂θIJ 1 (N s1 · + s2 · + s3 · + s4 · 2 ∂u i 2 ∂u i ∂u i ∂u i I I J

(34)

334

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

and the stiffness matrix  N −1)  N −1) N −2   ∂ r˜I ∂ r˜I ∂ r˜I ∂ r˜I 1 (N ∂ 2r˜I 1 (N g1 · g2 · Ki j = · + s1 · + · 2 ∂u j ∂u i ∂u i ∂u j 2 ∂u j ∂u i I I J  ∂ 2r˜I ∂ r˜ J ∂ r˜ J ∂ 2r˜ J ∂θIJ ∂θIJ + s2 · + g3 · · + s3 · + g4 · · . ∂u i ∂u j ∂u j ∂u i ∂u i ∂u j ∂u j ∂u i The coefficients s1 , s2 , s3 , s4 , g1 , g2 , g3 , g4 appearing in Eqs. (34) and (35) are respectively       −q−1 −q − p−1 −p −1 s1 = ∂Φ2 /∂ r˜I = a · exp h −1 + q r˜I − h −2 −Bpr˜I B r˜I − r˜I I − h0 I     −1 −1 −2 2 θ − cos θ s2 = ∂Φ3 /∂ r˜I = b · exp γ h −1 + γ h − 2γ h −γ h (cos ) IJ IJ0 I J I 0     −1 −1 −1 2 s3 = ∂Φ3 /∂ r˜ J = b · exp γ h I + γ h J − 2γ h 0 (cos θIJ − cos θIJ0 ) −γ h −2 J   −1 −1 s4 = ∂Φ3 /∂θIJ = 2b · exp γ h −1 + γ h − 2γ h (cos θIJ − cos θIJ0 ) (− sin θIJ ) I J 0   − p−2 −q−2 −1 g1 = ∂ 2 Φ2 /∂ r˜I2 = a · exp h −1 Bp ( p + 1) r˜I − q (q + 1) r˜I I − h0      − p−1 −q−1 −p −q −4 −3 + 2h −2 Bp r ˜ − q r ˜ + h + 2h B r ˜ − r ˜ I I I I I I I     −1 −1 −1 −3 2 2 g2 = ∂ Φ3 /∂ r˜I = b · exp γ h I + γ h J − 2γ h 0 (cos θIJ − cos θIJ0 )2 γ 2 h −4 I + 2γ h I     −1 −1 −3 2 2 −4 g3 = ∂ 2 Φ3 /∂ r˜ J2 = b · exp γ h −1 + γ h − 2γ h θ − cos θ γ h + 2γ h (cos ) IJ IJ0 I J J J 0    −1 −1 −1 2 2 2 1 − 2 cos θIJ + cos θIJ0 cos θIJ g4 = ∂ Φ3 /∂θIJ = 2b · exp γ h I + γ h J − 2γ h 0 h I = r˜I − r0∗ ;

h J = r˜ J − r0∗ ;

(35)

(36)

h 0 = 1 − r0∗

in which the quantities     −1 3E 2V −3 −2 −3 2 · · B p 2 + p + ph −2 + 2h − q + q + qh + 2h 0 0 0 0 (1 − 2ν) N (N − 1) 45E (1 − 4ν) V b= · 16 (1 + ν) (1 − 2ν) N (N − 1) (N − 2)

a=

B=

h −2 0 +q h −2 0 + p

(37)

.

6. Precision check for continuum representation In the present method, the parameter calibration is based on the ideal unit cell. However, in using the present method, the specific unit cell deviates from the ideal one both in the bond number and bond distribution. To check how precisely the discrete system can represent the continuum, we ‘mesh’ a cubic body (Fig. 5) with the unstructured tetrahedral unit cell. Before simulation, we firstly prescribe a Young’s modulus and Poisson ratio. Then we will use Eqs. (34), (35) to simulate the initial apparent Young’s modulus and Poisson ratio of this discrete system. The aim is to check whether the simulated apparent Young’s modulus and Poisson ratio are equal to the prescribed values. We take the shape parameters p = 12, q = 6, r0∗ = 1.05, γ = 0.8. To compute the apparent Young’s modulus and Poisson ratio, a uniaxial tensile displacement is applied on the cubic body. To guarantee the initial state, a very small displacement step, say 1.0 × 10−6 m, is imposed on the top of the cubic body. For the volume of the cubic body is √ 3 3 1 m , the quantity S = 1/N E can be taken as a measurement of unit cell size in which NE means the total number of unit cells. Table 1 lists the simulated Young’s modulus E ∗ , Poisson ratio νx∗ , ν y∗ and their prescribed values. The simulated Poisson ratios νx∗ , ν y∗ are respectively computed by the x − z and y − z deformation ratio. ν ∗ is the mean   value of νx∗ , ν y∗ , ν ∗ = 0.5 νx∗ + ν y∗ . In all cases, the prescribed Young’s modulus is E = 10.0 GPa. From Table 1,

335

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

Table 1 Comparison between the simulated Poisson ratios νx∗ , ν y∗ , Young’s modulus E ∗ (unit: GPa) and their prescribed values with different unit cell numbers. (The prescribed Young’s modulus E = 10.0 GPa in all cases. NUC stands for the total number of unit cells.) Meshing case NUC ν = 0.05

ν = 0.10

ν = 0.15

ν = 0.20

ν = 0.25

ν = 0.30

Case1

Case2

Case3

Case4

Case5

Case6

Case7

Case8

256

2046

12 288

24 373

52 565

115 262

231 846

433 047

E∗

9.37

9.33

9.55

9.99

9.66

9.71

9.73

9.72

νx∗ ν y∗

0.0378 0.0414

0.0421 0.0468

0.0469 0.0479

0.0454 0.0441

0.0441 0.0455

0.0454 0.0448

0.0447 0.0447

0.0445 0.0446

E∗

9.33

9.36

9.62

νx∗ ν y∗

0.0821 0.0833

0.0879 0.0887

0.0914 0.0922

E∗

9.41

9.37

9.68

νx∗ ν y∗

0.1268 0.125

0.1344 0.1303

0.1366 0.137

10.03 0.0936 0.0922 10.05 0.1427 0.1412

9.65

9.76

9.78

9.77

0.0923 0.093

0.0913 0.0906

0.0906 0.0906

0.0906 0.0907

9.61

9.80

9.82

9.80

0.142 0.1414

0.1379 0.137

0.1371 0.1371

0.1373 0.1376

E∗

9.45

9.36

9.73

9.48

9.80

9.83

9.80

νx∗ ν y∗

0.172 0.1665

0.1823 0.1719

0.1829 0.183

0.1932 0.1915

0.1938 0.1917

0.1858 0.1846

0.1849 0.1848

0.1853 0.1856

10.03

E∗

9.45

9.27

9.71

9.91

9.16

9.72

9.75

9.70

νx∗ ν y∗

0.2182 0.2081

0.2324 0.2137

0.2316 0.231

0.2461 0.2442

0.2497 0.2449

0.236 0.2344

0.2352 0.2349

0.2355 0.2361

E∗

9.27

8.93

9.43

9.46

8.27

9.30

9.33

9.26

νx∗ ν y∗

0.2665 0.2499

0.2863 0.2569

0.2851 0.2838

0.3037 0.3019

0.3144 0.3046

0.2912 0.2892

0.2907 0.2902

0.2908 0.2917

Fig. 5. A cubic body discretized into discrete system with the tetrahedral unit cell. (The cubic body dimensions are 1 m × 1 m × 1 m.)

it is seen that both the simulated Young’s modulus and Poisson ratio converge to fixed values with increasing the unit cell number, namely decreasing the unit cell size. Fig. 6 shows the relative errors of the simulated Young’s modulus and Poisson ratio. From Fig. 6a, it is seen that the relative errors of Young’s modulus approach to fixed values with decreasing the unit cell size. For the cases of Poisson ratio smaller than 0.25, the final relative error is about −2% while for the cases of Poisson ratio 0.3, the final relative error is about −7%. From Fig. 6b, it is seen that the relative errors of the simulated Poisson ratio also approach to fixed values with decreasing unit cell size. The bigger the prescribed Poisson ratio is, the smaller the relative error of the simulated Poisson ratio. But for different prescribed Poisson ratios, the final relative errors of the simulated Poisson ratio are scattered within the range of −3% ∼ −10%. Fig. 6 suggests that the discrete system with specific lattice bond cells can reach a very high precision with decreasing unit cell size.

336

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

Fig. 6. Errors of the lattice bond cell system in representing continuum (a) the relative error of the simulated Young’s modulus versus the unit cell size; (b) the relative error of the simulated Poisson ratio versus the unit cell size.

Fig. 7. Notched beam subjected to impact loading (a) beam dimensions (unit: mm) and boundary conditions (the thickness B = 25.4 mm); (b) illustration of meshing scheme.

7. Simulation examples of dynamic fracture 7.1. Notched beam subjected to impact loading To examine the performance of the present method in simulating dynamic fracture growth, we use it to simulate a notched beam subjected to impact loading, which was reported in [19,20]. The beam specimen is shown in Fig. 7a and the discretizing scheme with tetrahedral unit cell is shown in Fig. 7b. The total node number is 37 462 and the

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

337

Fig. 8. The reproduced uniaxial stress–strain curves for the three cases.

Fig. 9. Crack facet created in an over-stretched tetrahedral unit cell (a) mode I and (b) mode II.

total unit cell number is 210 653. The beam is subjected to a velocity [19]  v t/t for t ≤ t1 v0 (t) = 1 1 v1 for t > t1

(38)

at the upper edge shown as Fig. 7a, where t1 = 1.96 × 10−4 s and v1 = 0.065 m/s. To check the influence of material parameters on the dynamic fracture propagation, the following three cases with different input parameters are simulated. Case I: E = 50.0 GPa, r0∗ = 1.0114, the reproduced σb = 2.99 MPa, εt = 0.15 × 10−3 ; Case II: E = 50.0 GPa, r0∗ = 1.0144, the reproduced σb = 4.78 MPa, εt = 0.24 × 10−3 ; Case III: E = 31.4 GPa, r0∗ = 1.0144, the reproduced σb = 2.99 MPa, εt = 0.24 × 10−3 in which σb denotes the uniaxial tensile strength, εt the uniaxial strain value corresponding to σb . In all cases, the material density ρ = 2400 kg/m3 , Poisson ratio ν = 0.2, p = 12, q = 6, γ = 0.8. The reproduced √uniaxial stress–strain curves of the three cases are shown in Fig. 8. According to the shear wave velocity cs = E/ (2 (1 + ν) ρ), the Rayleigh wave velocity can be estimated as c R ≈ 0.9cs . So, the Rayleigh wave velocity of the Case I, Case II and Case III are c R = 2652 m/s, c R = 2652 m/s and c R = 2100 m/s, respectively. The explicit integration scheme is adopted with the time interval ∆t = 0.05 µs. Since the present model is a ‘strain’ softening one and no separate fracture criterion is adopted here, we will use a criterion to identify the over-stretched unit cells. The array of these over-stretched unit cells could be approximately considered as the crack trajectory. For

338

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

Fig. 10. Fracture propagation process in the Case III at the time (a) 363 µs; (b) 388.5 µs; (c) 402 µs; (d) 417 µs; (e) 430.5 µs; (f) 49.75 µs. The left column is the crack at different times whereas the right column is the corresponding node configuration. (The node displacements are magnified 1000 times.)

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

339

Fig. 11. The created rough crack faces in (a) Case I; (b) Case II; (c) Case III.

the tetrahedral unit cell, the unit cell in which at least three bonds exceed the following critical deformation will be identified as an over-stretched one. r˜ = 1.0 + 10.0εt .

(39)

In an over-stretched unit cell, a crack facet is created, shown in Fig. 9. Fig. 10 shows the fracture propagation process in Case III. From Fig. 10 it is seen that the dynamic crack begins to grow approximately at the 363 µs after the impact load was imposed. The dynamic crack fast propagates towards the loading point. Fig. 11 shows the final created crack surfaces in the three cases. From Fig. 11 it is noted that the crack surfaces are not smooth, but rough, which agrees with the topography of fracture surface in the usual experiment of rock or concrete. Fig. 12 shows the detailed kinetic of crack frontier during propagation in Case III. In Fig. 12 an interesting phenomenon is found that the crack frontier does not move forwards continuously, but sometimes discontinuously. At the beginning of crack growth, the discontinuous growth behavior is very obvious, shown as Fig. 12a–e. However, after a certain time of propagation, the frontier tends to move continuously, shown as Fig. 12f. The same phenomenon also can be found in other cases. Referring to Fig. 14, it is found that the discontinuous propagation happens in the first stage in which the crack speed is being accelerated and the continuous propagation happens in the later stage in which the crack speed is being decelerated. The fact that amount of energy stored in the specimen tends to quickly release at the beginning may be responsible for the discontinuous propagation phenomenon.

340

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

Fig. 12. Crack frontier at the time (a) 364 µs; (b) 378 µs; (c) 383.25 µs; (d) 393.75 µs; (e) 399.75 µs; (f) 484.5 µs. (Case III.)

Fig. 13. Crack trajectory viewed from the lateral side of beam in (a) Case I; (b) Case II; (c) Case III; (d) crack paths reported in [19,20].

Fig. 13a–c shows the fracture trajectory viewed from the lateral side of specimen. The simulated fracture trajectories generally agree with the observation in the experiment (Fig. 13d) in [20]. Because the fracture in 3D cases is a surface rather than a line, it is improper to characterize the crack growth with the line speed as does in 2D cases. However, in the present simulation, we could use an apparent line speed to describe the crack growth. To obtain the line speed of crack growth, we firstly compute the average crack face creation speed over the time interval 15∆t = 0.75 µs and then divide it with the thickness of the beam. The simulated line speed of crack growth is shown in Fig. 14. From Fig. 14a it is seen that the bigger the critical strain εt is, the later the crack begins to grow. In contrast to εt , the bigger the Young’s modulus is, the earlier the crack begins to grow. Once the crack begins to grow, the general trend of the crack growth speed with the crack length are almost the same, shown in Fig. 14b. But in all cases, the crack begins to grow earlier than that simulated in [21]. This is because the fracture energy is not preserved in the present simulation. The present model is a ‘strain’ softening one, which suffers from the ‘element’ size sensitivity. If no regularization technique is adopted, the fracture energy will approach to zero with the unit cell size approaching to zero. Consequently, the fracture will propagate with few strain energy release rate.

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

341

Fig. 14. Crack growth speed in three cases (a) crack speed versus time; (b) crack speed versus fracture length.

According to [22], the product h S is approximately equal to the intrinsic J -integral of material, in which h is the ‘element’ height at fracture tip and S is the area under the uniaxial tensile stress–strain curve. According to Fig. 8, S ≈ 0.73 kPa for Case I, S ≈ 10.4 kPa for Case II and S ≈ 1.2 kPa for Case III. The mean height of unit cell in the trajectory of the extended fracture is about 0.001 m. So, the fracture energy preserved in the present simulation is only about G f = 0.73–10.4 J/m2 , which is much lower than the usual concrete. Hence, the simulated crack begins to grow much earlier than that reported in [21]. According to [23], a dynamic crack cannot branch when its maximum velocity is below a threshold value 0.35c R . For the Case I, Case II and Case III, the threshold value is 928 m/s, 928 m/s and 735 m/s, respectively. Fig. 14 shows that the crack velocity in different cases is below their corresponding threshold value. Hence, the fracture should not branch under such crack speed. Just as predicted, we have not observed the branching phenomenon in this simulation, shown as Fig. 11. 7.2. Effect of Poisson ratio on the dynamic fracture Compared with the usual lattice method, the correlated lattice bond cell model can present variable Poisson ratios. To examine the necessity of presenting the variable Poisson ratios, the beam of Case I in Section 7.1 is re-simulated respectively with the Poisson ratios ν = 0.1 and ν = 0.25. The simulated results demonstrate that a lower Poisson ratio can lead to a higher crack speed, shown in Fig. 15a. The maximum crack speed with ν = 0.1 is about 1000 m/s whereas that with ν = 0.25 is about 800 m/s. The relative difference between them is about 1/5. Moreover, the lower the Poisson ratio is, the earlier the crack begins to grow, shown in Fig. 15b. Besides the crack speed, the Poisson ratio also influences the roughness of crack face. The higher the Poisson ratio is, the rougher the created crack face, shown in Fig. 16. So, the Poisson ratio effect on the dynamic fracture could not be neglected. The representation of variable Poisson ratios is quite necessary. 7.3. Crack branching To further examine the performance of the present method in dynamic fracture simulation, we use it to simulate the classical crack branching example. The simulating specimen is a notched thin slab subjected to a tensile stress on its two sides, shown in Fig. 17a. The tensile stress is applied on the slab as  σ t/t for t ≤ t1 σt (t) = 0 1 (40) σ0 for t > t1 in which t1 = 50.0 µs. In different cases, σ0 takes different values. The ‘meshing’ scheme of the thin slab is illustrated as Fig. 17b, in which the tetrahedron unit cell is adopted. The total node number is 32964. Take the Case I of Section 7.1 as the material parameter input. The explicit time integration scheme is adopted and the time interval is taken as ∆t = 0.05 µs. The simulated crack branching patterns with different notch depths H and tensile stresses σ0 are shown in Fig. 18. From Fig. 18 it is seen that the fracture

342

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

Fig. 15. Effect of Poisson ratio on crack speed (a) crack speed versus crack length; (b) crack speed versus time.

a

b

Fig. 16. Effects of Poisson ratio on the created crack face (a) ν = 0.10; (b) ν = 0.25. (The left column is the detailed crack face while the right column is the crack trajectory viewed from the lateral side of specimen.)

firstly propagates in a stable manner and then branches. The higher the imposed traction is, the more seriously the fracture branches. During the whole propagation process, the dynamic fracture experiences a series of branching points. For example, in the case of (H = 50 mm, σ0 = 1.5 MPa), the fracture successively experiences the first and the second branching point, shown in Fig. 19. If we pick up the crack speed from the propagation initiation to the first branching point, we can find that the crack speeds in the four cases are almost the same, shown as Fig. 20a. This indicates that when the crack speed is below a certain value, the crack does not branch. Fig. 20a also suggests that the higher the tensile stress is, the shorter the crack extended in the stable manner. The length of the stable crack is almost irrelevant to the notch depth, but relevant to the applied traction. Fig. 20b shows the crack speed versus the propagation time, from which it is found that the time to fracture is only related to the imposed tensile stress, almost unrelated to the notch depth. The bigger the tensile stress is, the earlier the crack begins to grow. According to [24–26], the crack speed range for instability onset is approximately 0.4c R –0.57c R , namely 1061 m/s–1512 m/s for the material of Case I. The crack speeds shown in Fig. 20 just fall into this range. From

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

343

Fig. 17. The simulation specimen for crack branching (a) specimen dimensions (H means the notch depth); (b) ‘meshing’ scheme.

Fig. 18. Simulated crack branching patterns in the specimen (a) H = 50 mm, σ0 = 1.5 MPa; (b) H = 50 mm, σ0 = 2.9 MPa; (c) H = 20 mm, σ0 = 1.5 MPa; (d) H = 20 mm, σ0 = 2.9 MPa. (The time to propagate in (a)–(d) is 29.0 µs, 21.5 µs, 29.5 µs, 21.5 µs, respectively. The crack patterns shown in (a)–(d) are at the time t = 70.0 µs, t = 60.5 µs, t = 95.5 µs, t = 66.5 µs, respectively.)

Fig. 19. Detailed crack branching pattern in the case of (H = 50 mm, σ0 = 1.5 MPa).

344

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

Fig. 20. Simulated crack speed before the first branching point (a) crack speed versus extended fractured length; (b) crack speed versus time.

Fig. 20 it is also interesting to note that the crack speed suddenly drops down just before crack branching. This is consistent with the prediction of [23] that the crack branching occurs not at the point of maximum velocity, but some time later, after the crack velocity has dropped below the maximum. In the crack branching examples, the discontinuous propagation phenomena are also observed. Fig. 21 shows the detailed propagation process of the case of (H = 20 mm, σ0 = 2.9 MPa). From Fig. 21 it is seen that the crack does not always propagate continuously, sometimes discontinuously (Fig. 21b, c). In the dynamic fracture experiments carried in [27–32], the ‘mirror–mist–hackle’ phenomenon is observed in that the dynamic crack growth successively experiences the smooth, rough and zigzag stage. In the molecular dynamic simulation conducted by Abraham et al. [24,25], such ‘mirror–mist–hackle’ phenomenon is also reproduced. To examine whether the present method can capture this typical feature of dynamic fracture, we will simulate an edge-notched rectangular plate under a constant strain rate ε˙ yy = 40.0 s−1 , shown in Fig. 22a. The lateral sides of specimen are stress free. The GMESH [33] is adopted to discretize the specimen with the tetrahedral unit cells, shown in Fig. 22b. The total node number is 35 536 and the total unit cell number is 161 935. Take the Case III of Section 7.1 as the material parameter input and the time interval is taken as ∆t = 0.0001 µs. The simulated crack propagation process is shown in Fig. 23, from which it is seen that the ‘mirror–mist–hackle’ phenomenon is basically represented. Through the simulation examples, it is found that the present method can capture the characters of the dynamic fracture propagation. However, we have to point out that the modified SW potential is a hyperelastic potential.

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

345

Fig. 21. Discontinuous propagation and branching process of dynamic crack at the time (a) 32 µs; (b) 42 µs; (c) 44 µs; (d) 53 µs; (e) 61 µs. (H = 20 mm, σ0 = 2.9 MPa.)

Fig. 22. An edge-notched rectangular specimen (a) dimensions (unit: mm); (b) meshing scheme.

Governed by such hyperelastic potential, the bond deformation is reversible. So, the present method cannot account for the plastic deformation effect. But most solids exhibit some plasticity near the crack tips even the solids are nominally brittle. To make the SW potential capable of accounting for the plastic deformation effect at fracture tip, the SW potential will have to be further modified to consider the irreversible deformation. 8. Conclusion remarks The nonlinear elastic solid is modeled with the correlated lattice bond cell method. Governed by the modified SW potential, the discrete system composed of the correlated lattice bond cells can represent the variable Poisson ratios. The bond parameters in a unit cell are calibrated as the function of macro material constants, the volume of unit cell and the total bond number in this unit cell. Because more elaborate fracture mechanisms are incorporated into the force–displacement constitutive relation of a unit cell through the modified SW potential, the present method can directly simulate the dynamic fracture initiation, arrest, propagation and branching behaviors in the finite deformation case without any separate fracture criterion. The simulation examples demonstrate that the present method can capture most characters of the dynamic fracture observed in the physical experiments and the molecular dynamics simulations. It is suggested that the present method is an efficient approach to dynamic fracture simulation.

346

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347

Fig. 23. Simulated result for ‘mirror–mist–hackle’ examination (a) fracture process; (b) overview at the time t = 3.50 µs.

Acknowledgments The present work is supported by the National Natural Science Foundation of China (No. 11172172) and the National Basic Research Program of China (Grant No. 2011CB013505), which are gratefully acknowledged. Appendix A. Supplementary data Supplementary material related to this article can be found online at http://dx.doi.org/10.1016/j.cma.2014.06.036. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]

A. Hrennikoff, Solution of problems of elasticity by the framework method, J. Appl. Mech. (1941) A169–A175. D. Greenspan, Computer-Oriented Mathematical Physics, Pergamon Press, Oxford, 1981. D. Greenspan, Supercomputer simulation of cracks and fractures by quasi-molecular dynamics, J. Phys. Chem. Solids 50 (1989) 1245–1249. G. Wang, M. Ostoja-Starzewski, Particle modeling of dynamic fragmentation—I: theoretical considerations, Comput. Mater. Sci. 33 (2005) 429–442. D.V. Griffiths, G.G.W. Mustoe, Modeling of elastic continua using a grillage of structural elements based on discrete element concepts, Internat. J. Numer. Methods Engrg. 50 (2001) 1759–1775. G. Lilliu, J.G.M. Van Mier, 3D lattice type of fracture model for concrete, Eng. Fract. Mech. 70 (2003) 927–941. G.F. Zhao, J. Fang, J. Zhao, A 3D distinct lattice spring model for elasticity and dynamic failure, Int. J. Numer. Anal. Methods Geomech. 35 (2011) 859–885. G. Wang, A. Al-Ostaz, A.H.D. Cheng, P.R. Mantena, Hybrid lattice particle modeling: theoretical considerations for a 2D elastic spring network for dynamic fracture simulations, Comput. Mater. Sci. 44 (2009) 1126–1134. G. Wang, A. Al-Ostaz, A.H.D. Cheng, P.R. Mantena, Hybrid lattice particle modeling of wave propagation induced fracture of solids, Comput. Methods Appl. Mech. Engrg. 199 (2009) 197–209. M. Ostoja-Starzewski, Lattice models in micromechanics, Appl. Mech. Rev. 55 (2002) 35–60. Z.N. Zhang, Discretized virtual internal bond model for nonlinear elasticity, Int. J. Solids Struct. 50 (2013) 3618–3625. H.J. Gao, P. Klein, Numerical simulation of crack growth in an isotropic solid with randomized internal cohesive bond, J. Mech. Phys. Solids 46 (1998) 187–218. F.H. Stillinger, T.A. Weber, Computer simulation of local order in condensed phases of silicon, Phys. Rev. B 31 (1985) 5262–5271. L. Pizzagalli, J. Godet, J. Guenole, S. Brochard, E. Holmstrom, K. Nordlund, T. Albaret, A new parametrization of the Stillinger–Weber potential for an improved description of defects and plasticity of silicon, J. Phys. Condens. Matter. 25 (2013) 055801. Z.N. Zhang, Y.X. Chen, H. Zheng, A modified Stillinger–Weber potential-based hyperelastic constitutive model for nonlinear elasticity, Int. J. Solids Struct. 51 (2014) 1542–1554.

Z. Zhang, Y. Chen / Comput. Methods Appl. Mech. Engrg. 279 (2014) 325–347 [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33]

347

M.J. Buehler, F.F. Abraham, H. Gao, Hyperelasticity governs dynamic fracture at a critical length scale, Nature 426 (2003) 141–146. M.J. Buehler, H. Gao, Dynamical fracture instabilities due to local hyperelasticity at crack tips, Nature 439 (2006) 307–310. E.B. Tadmor, M. Ortiz, R. Phillips, Quasicontinuum analysis of defects in solids, Phil. Mag. A 73 (1996) 1529–1563. T. Belytschko, M. Tabbara, Dynamic fracture using element-free Galerkin methods, Internat. J. Numer. Methods Engrg. 39 (1996) 923–938. R. John, S.P. Shah, Mixed-mode fracture of concrete subjected to impact loading, J. Struct. Eng. 116 (1990) 585–602. J.H. Song, P.M.A. Areias, T. Belytschko, A method for dynamic crack and shear band propagation with phantom nodes, Internat. J. Numer. Methods Engrg. 67 (2006) 868–893. Z.N. Zhang, H.J. Gao, Simulating fracture propagation in rock and concrete by an augmented virtual internal bond method, Int. J. Numer. Anal. Methods Geomech. 36 (2012) 459–482. S.J. Zhou, P.S. Lomdahl, R. Thomson, B.L. Holian, Dynamic crack processes via molecular dynamics, Phys. Rev. Lett. 76 (1996) 2318–2321. F.F. Abraham, D. Brodbeck, R.A. Rafey, W.E. Rudge, Instability dynamics of fracture: a computer simulation investigation, Phys. Rev. Lett. 73 (1994) 272–275. F.F. Abraham, D. Brodbeck, W.E. Rudge, X.P. Xu, A molecular dynamics investigation of rapid fracture mechanics, J. Mech. Phys. Solids 45 (1997) 1595–1619. E. Sharon, J. Fineberg, Confirming the continuum theory of dynamic brittle fracture for fast cracks, Nature 397 (1999) 333–335. J. Fineberg, S.P. Gross, M. Marder, H.L. Swinney, Instability in dynamic fracture, Phys. Rev. Lett. 67 (1991) 457–460. T. Cramer, A. Wanner, P. Gumbsch, Energy dissipation and path instabilities in dynamic fracture of silicon single crystals, Phys. Rev. Lett. 85 (2000) 788–791. K. Ravi-Chandar, W.G. Knauss, An experimental investigation into dynamic fracture: I. Crack initiation and arrest, Int. J. Fract. 25 (1984) 247–262. K. Ravi-Chandar, W.G. Knauss, An experimental investigation into dynamic fracture: II. Microstructural aspects, Int. J. Fract. 26 (1984) 65–80. K. Ravi-Chandar, W.G. Knauss, An experimental investigation into dynamic fracture: III. On steady-state crack propagation and crack branching, Int. J. Fract. 26 (1984) 141–154. K. Ravi-Chandar, W.G. Knauss, An experimental investigation into dynamic fracture: IV. On the interaction of stress waves with propagating cracks, Int. J. Fract. 26 (1984) 189–200. C. Geuzaine, J.-F. Remacle, Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities, Internat. J. Numer. Methods Engrg. 79 (2009) 1309–1331.