Crack propagation simulation in brittle elastic materials by a phase field method

Theoretical & Applied Mechanics Letters 9 (2019) 339-352  

Contents lists available at ScienceDirect

Theoretical & Applied Mechanics Letters journal homepage: www.elsevier.com/locate/taml    

Article

Crack propagation simulation in brittle elastic materials by a phase field method Xingxue Lua,b, Cheng Lia,*, Ying Tiea, Yuliang Houa, Chuanzeng Zhangb,*  

a b

School of Mechanical Engineering, Zhengzhou University, Zhengzhou 450001, China Department of Civil Engineering, University of Siegen, D-57076 Siegen, Germany

H  I  G  H  L  I  G  H  T  S

• Implementation and application of a phase field method for crack propagation simulation in brittle elastic materials. • Investigation of the effects of the loading type (tension and shear), boundary conditions, as well as size, location and orientation of the initial crack on the crack propagation path and force-displacement relation in homogenous and linear elastic plates. • Investigation of the effects of the loading type (tension and shear), boundary conditions, bi-material interface, and initial crack location on the crack propagation path and force-displacement relation in linear elastic bi-material plates. • Enhancement of the implementation and application convenience of the phase field method (ParaView as postprocessor, use of the colon function in MATLAB, etc.). A  R  T  I  C  L  E      I  N  F  O

A  B  S  T  R  A  C  T

 

 

Article history: Received 19 June 2019 Received in revised form 20 August 2019 Accepted 22 August 2019 Available online 24 August 2019 *This article belongs to the Solid Mechanics.

   

Keywords: Brittle fracture Phase field method Crack propagation Finite element method

To overcome the difficulties of re-meshing and tracking the crack-tip in other computational methods for crack propagation simulations, the phase field method based on the minimum energy principle is introduced by defining a continuous phase field variable Á (x)∈[0,1] to characterize discontinuous cracks in brittle materials. This method can well describe the crack initiation and propagation without assuming the shape, size and orientation of the initial crack in advance. In this paper, a phase field method based on Miehe's approach [Miehe et al., Comp. Meth. App. Mech. Eng. (2010)] is applied to simulate different crack propagation problems in twodimensional (2D), isotropic and linear elastic materials. The numerical implementation of the phase field method is realized within the framework of the finite element method (FEM). The validity, accuracy and efficiency of the present method are verified by comparing the numerical results with other reference results in literature. Several numerical examples are presented to show the effects of the loading type (tension and shear), boundary conditions, and initial crack location and orientation on the crack propagation path and force-displacement curve. Furthermore, for a single edge-cracked bi-material specimen, the influences of the loading type and the crack location on the crack propagation trajectory and force-displacement curve are also investigated and discussed. It is demonstrated that the phase field method is an efficient tool for the numerical simulation of the crack propagation problems in brittle elastic materials, and the corresponding results may have an important relevance for predicting and preventing possible crack propagations in engineering applications. ©2019 The Authors. Published by Elsevier Ltd on behalf of The Chinese Society of Theoretical and Applied Mechanics. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).  

 

    

 

* Corresponding author. E-mail address: [email protected] (C. Li), [email protected] (Ch. Zhang).

  http://dx.doi.org/10.1016/j.taml.2019.06.001 2095-0349/© 2019 The Authors. Published by Elsevier Ltd on behalf of The Chinese Society of Theoretical and Applied Mechanics. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

340

X.X. Lu et al. / Theoretical & Applied Mechanics Letters 9 (2019) 339-352

1     Introduction Fracture is a major failure mode in engineering materials and structures. However, fracture is usually a rather complex phenomenon and its adequate description is quite difficult in practice. Therefore, it is an important task to simulate and predict the initiation, growth and branching of possible cracks by computational models for practical engineering applications. The classical theory of brittle fracture in elastic solids, that a crack propagates if the energy release rate reaches a critical value, was first proposed by Griffith [1] and Irwin [2]. Even though it provides a criterion for crack propagation, it fails to account for the crack nucleation, curvilinear crack paths, crack kinking or branching. In last decades, many numerical simulation methods were developed to model different crack propagation problems. However, each has its own advantages and drawbacks. The well-established classical finite element method (FEM) [3, 4 ] usually needs to adjust the mesh or re-meshing in the crack propagation simulation. On the other hand, the node splitting [5], interface (cohesive) element formulation [6–8], and discrete element model [9] are strongly mesh dependent. Moës et al. [10] introduced the extended FEM (X-FEM) to avoid the abovementioned problems by using a local enrichment in the shape functions of the FE formulation [11, 12]. A computational scheme based on the configurational forces was adopted for crack propagation simulation by Gürses and Miehe [13]. Moreover, Francfort and Marigo [14] proposed a variational model of quasi-static crack evolution based on the energy minimization which can overcome the weaknesses of the original approach of Griffith. Some other related research works have been also done in the past [15–19]. Bourdin et al. [17, 20] studied the regularized setting of their framework in Γ − convergence based on the work on image segmentation of Mumford and Shah [21]. A regularized approximation of the brittle fracture by a scalar auxiliary variable was described in the phase field approaches [22–24], which assumes that a sharp crack is continuous in the solid. The variable connecting the unbroken and the broken states of the material is considered as a phase field. In recent years, the research related to the phase field formulation is receiving more and more attention. The main reason lies in the fact that in the phase field method, no additional ad-hoc criteria are required to simulate the crack initiation, propagation, kinking, branching and so on. Among many recent research works on this topic, we just mention a few of them, such as the quasistatic fracture [25–29], dynamic crack propagation [30–33], finite deformation [34, 35], multi-physics problems [36, 37], and so on [38–42]. Moreover, the phase field method has been also extended to simulate the crack growth problems in the anisotropic materials [43–45]. Although many research works on the phase field method have been reported in literature so far, few works have been yet devoted to a detailed study on the effects of the key factors such as the size, location and orientation of the initial crack, and the loading type (tensions and shear) on the crack propagation characteristics. The main purpose of this paper is to present an efficient and accurate phase field method for crack propagation simulation in two-dimensional (2D), isotropic and linear elastic materials, and to reveal the aforementioned essential effects. The corresponding results can be used to predict and avoid possible crack propagations in practical engineering applications,

which is an important and indispensable issue for reliable design and optimization of advanced engineering materials and structures. As previously mentioned, the continuous phase field variable ϕ is introduced to characterize the discontinuous cracks. So, we write the phase field, displacement field and mesh data at each loading step in the visualization toolkit (VTK) file format for plots to be viewed by using ParaView, which is an open-source, multi-platform data analysis and visualization software. We can optionally check and analyze the data and conveniently follow the process of the crack initiation and growth. It avoids the output of the graphical results which cannot be further analyzed. Besides, the present phase field code is implemented in MATLAB. We use the colon function to create a vector of indices to select rows, columns, or elements of arrays, especially select the number of the elements, and then multiply the arrays. In this way, a higher computational efficiency can also be obtained on the computers which cannot perform a parallel computing. The remainder of this paper is structured as follows. In the next section, the key idea and the theoretical background of the phase field method based on Miehe's approach [22] are briefly reviewed. Section 3 describes the numerical implementation. The validation of the implemented phase field method is first presented in Sect. 4. Then, some representative numerical examples are given and discussed to reveal the influences of the loading type (tension and shear), and the size, location and orientation of an initial crack on the crack propagation path and the force-displacement curve, respectively. Moreover, the effects of the loading type, and the location of an initial crack in a bi-material plate on the crack propagation characteristics are also studied. Finally, Sect. 5 gives some conclusions. 2     Description of the phase field method 2.1    Phase field approximation of the crack topology Let us consider a homogeneous, isotropic and linear elastic domain Ω ⊂ R 2 with a crack Γ and the boundary ∂Ω as depicted in Fig. 1. In a regularized framework, the crack geometry is a "smeared" representation defined by an auxiliary field variable ϕ(x) ∈ [0, 1], x ∈ Ω, which is denoted as the crack phase field as illustrated in Fig. 1(b). The un-cracked and fully cracked states of the considered domain are characterized by ϕ = 0 and ϕ = 1, respectively. According to Griffith's theory [1], an existing crack will propagate when the energy release rate G associated with the crack extension exceeds a critical value equal to the material fracture toughness G c . However, the Griffith's theory has a shortcoming which cannot accommodate the crack nucleation or predict the branching of fracture. The variational approach to fracture mechanics proposed by Francfort and Marigo [14] introduces the following energy functional for the cracked body ∫ E (u, Γ) = E d (u, Γ) +E s (Γ) =

ψ (ε (u)) dΩ +G c H1 (Γ),

(1)

Ω

where E d (u, Γ) is the elastic energy stored in the cracked body, E s (Γ) represents the energy required to create the new crack according to the Griffith theory, ψ is the elastic energy density function, ε is the strain tensor, u is the displacement vector, G c is

X.X. Lu et al. / Theoretical & Applied Mechanics Letters 9 (2019) 339-352 a

 

b

t

t

∂Ωt

1 ∂Ωt 0.5

Φ(x)

Γ ∂Ωu

∂Ωu

Ω

u

341

Ω

u

0

 

Fig. 1.   2D sharp and diffusive crack topology: a sharp crack model, b regularized representation of a crack by phase field approximation.

the critical energy release rate, and H1 is the Hausdorff surface measure giving the crack length. With the assumption of small deformation,( the strain tensor is related to the displacement ) vector by ε = ∇u + ∇uT /2 . To consider the crack phase field, the energy functional Eq. (1) can be expressed as the following functional ( ) E u, ϕ ≈

∫

( ) ψ ε (u) , ϕ dΩ +G c

Ω

∫

=

∫

Ω

∫ (

Ω

) 1 2 l ϕ + ∇ϕ · ∇ϕ dΩ, 2l 2

ε = ε+ +ε− ,

(6)

with the superscripts "±" describing the tension and compression modes. Based on the spectral decomposition of the

γ(ϕ)dΩ Ω

( ) ψ ε (u) , ϕ dΩ +G c

positive parameter used to maintain the well-posedness of the system for the fully cracked state ϕ = 1. The strain field is decomposed into positive and negative parts in order to account for the stress degradation only in tension, i.e.,

strain tensor ε = (2) ε± :=

where γ(ϕ, ∇ϕ) denotes the crack density function per unit volume, l is a regularization parameter describing the width of the smeared crack and recovers for l → 0 the sharp crack topology. It should be noted here that the regularization parameter l in the conventional phase field model may significantly affect the crack propagation path and the corresponding peak-load if it is not appropriately chosen. Frequently, two ways are suggested to determine l , either by considering it as a pure numerical parameter of the regularized model of the brittle fracture or as a real material parameter of a gradient damage model [19]. More detailed descriptions on how to choose an adequate value of l can be found in Refs. [19, 26, 31, 42, 46 ]. Moreover, the element size h may also have significant influences on the numerical convergence of the mechanical responses. Miehe et al. [23] suggested that the condition h ≤ l /2 has to be fulfilled. In this aspect, the length-scale insensitive phase field model developed by Wu [28], Wu and Nguyen [42] should be particularly mentioned.

2 ∑

i =1

εi ni ⊗ ni we obtain

2 ⟨ ⟩ ∑ ε i ± ni ⊗ ni ,

(7)

i =1

where εi are the eigenvalues and ni are the corresponding eigenvectors of ε. The bracket operators in Eqs. (5) and (7) are defined by 〈x〉± = (x ± |x|)/2. Two projection tensors are defined for the derivative of Eq. (6) with respect to the total strains as [ ] P+ := ∂ε ε+ (ε) , P− := I − P+ ,

(8)

where I is the identity matrix. The fourth-order tensors in Eq. (8) project the total strains onto their positive and negative parts, i.e. ε+ = P+ : ε and ε− = P− : ε. For more details one can refer to the Refs. [22, 24, 47, 48]. So the total potential energy function Eq. (2) can be rewritten as ( ) E u, ϕ =

∫ Ω

[ ] (1 − ϕ)2 + k ψ+ (ε) dΩ + ∫

+ Ω

∫ ψ− (ε) dΩ Ω

[ ] Gc 1 l ∇ϕ · ∇ϕ + ϕ2 dΩ. 2 l

(9)

2.2    Strain energy degradation in the cracked solid For a homogeneous, isotropic and linear elastic solid without damage, the elastic energy density function can be written as [ ] ψ0 (ε) = λtr2 [ε]/2 + µtr ε2 ,

(3)

where λ and µ are the Lamé constants. Since the elastic energy is released only in tension but not in compression, the degradation of the elastic energy due to fracturing can be described by ( ) ( ) ψ ε, ϕ = g (ϕ)+k ψ+0 (ε) + ψ−0 (ε), [( ) ] 2 ψ±0 (ε) = λ 〈tr [ε]〉2± /2 + µtr ε± ,

(4)

2.3    Variational formulation and phase field evolution equation As shown in Fig. 1, the exterior boundary of the considered domain can be decomposed into two parts ∂Ωu and ∂Ωt , ∂Ω = ∂Ωu ∪ ∂Ωt . The Dirichlet- and Neumann-type boundary conditions are applied. The first variation of the external work, δWext , can be written as ∫ δWext = Ω

(5)

with g (ϕ) = (1 − ϕ)2 as given in Refs. [22–24], k ≪ 1 is a small

∫ bδudΩ +

tδu∂Ω,

(10)

∂Ω

where b and t are the body force vector and boundary traction vector on the surfaces ∂Ω , respectively. For the first variation of the internal work or elastic strain energy, δWint, we have

342

X.X. Lu et al. / Theoretical & Applied Mechanics Letters 9 (2019) 339-352

∫ δWint = Ω

[ ] (1 − ϕ)2 + k σ+ δε+ dΩ + ∫

+

σ− δε− dΩ Ω

( ) −2 1 − ϕ δϕψ+0 (ε)dΩ +

Ω

 ∫

∫ ∫ Ω

( ) 1 G c l ∇ϕ · ∇δϕ + ϕδϕ dΩ, l

(11)

in which σ± = λ〈tr [ε]〉± 1 + 2µε± and [1] ={1, 1, 0}T . According to the principle of virtual work, we have δWi nt − δWext = 0, from which as well as Eqs. (10) and (11) we can obtain the following equations ) ( ) Gc ( ϕ − l 2 ∆ϕ − 2 1 − ϕ ψ+0 (ε) =0, l ∫

[ ] (1 − ϕ)2 + k σ+ δε+ dΩ +

Ω

∫

( ) −2 1 − ϕ δϕψ+0 (ε)dΩ +

Ω

x ∈ Ω, ∫

σ− δε− dΩ −

∫

Ω

(12) ∫ bδudΩ −

Ω

tδudΩ = 0,

∂Ω

( ) 1 G c l ∇ϕ · ∇δϕ + ϕδϕ dΩ=0. l

(13)

Equation (12) is the phase field evolution equation which implies that the phase field variable ϕ is a function of ψ+0 (ε). 3     Numerical implementation The coupled system of equations given in Eq. (13) can be solved numerically in a classical finite element framework. The displacement field and the phase field in 2D cases are approximated by u = Nu uu , ϕ = Nϕ ϕϕ ,

(14)

where uu and ϕϕ are the vectors consisting of the nodal displacement and phase field values, respectively, Nu and Nϕ are the shape functions associated with the nodes. The strain field and the gradient of the phase field can be obtained as ε = Bu u, ∇ϕ = Bϕ ϕϕ ,

(15)

in which Bu and Bϕ contain the derivatives of the shape functions. Then, Eq. (13) can be written as ∫ [ ∫ ∫ ∫ ] ( )2 1 − ϕ + k BTu σ+ dΩ+ BTu σ− dΩ− NTu bdΩ− NTu t∂Ω=0, Ω

Ω

Ω

∂Ω

( ) ] ∫ [ Gc G c l BTϕ ∇ϕ + + 2ψ+0 (ε) NTϕ ϕ − 2NTϕ ψ+0 (ε) dΩ=0. l

(16)

Ω

The staggered algorithm [22, 24] is used to solve Eq. (16). In order to guarantee the irreversible evolution of the phase field which corresponds to the crack healing, a local history strain energy function H is defined here as described in the work of Miehe et al. [22]. In each step, we first solve the phase field equation   ( ) ] ∫ ∫ [  Gc T T G c l Bϕ Bϕ + + 2Hn Nϕ Nϕ dΩ · ϕn+1 = 2NTϕ Hn dΩ,   l Ω

and then the momentum equation

Ω

(17)

[(

]( )2 ) 1 − ϕn+1 + k λR+n [1]T [1] + 2µP+n Bu dΩ

Ω

  ( ) + BTu λR−n [1]T [1] + 2µP−n Bu dΩ · un+1  Ω ∫ ∫ = NTu bdΩ + NTu tdΩ, ∫

Ω

(18)

∂Ω

the load step, R±n = 〈tr[εn ]〉± and Hn (x)= + , which is evaluated at the Gaussian points. max ψ0 (x, N )

where n {

N ∈[0,n]

∫ Ω



BTu

is

}

Before the first step, we assume H0 = 0. Then the staggered update of the phase field ϕ and the displacement field u can be calculated step by step. The phase field model based on Miehe's approach has been implemented in MATLAB. For the calculations of the strains, stresses and strain energy, we did not adopt the parallel calculation because of the available computer configuration. We used the colon function in MATLAB to create a vector of indices to select rows, columns, or elements of arrays, especially select the number of elements, and then multiplied the arrays. For instance, the following MATLAB commands were used to calculate the strains in each element at the Gaussian points: for i =1:ncompon for j =1:ndof strain(:, igaus, i) = strain(:, igaus, i) +B_u(:, i, j).*displa(:, j); end end Here, "ncompon" is the number of the strain components, "ndof" is the number of the degree of freedom in each element, the colon ":" denotes the number of the element, "igaus" is the number of the Gaussian point, "B_u" is the matrix of the derivatives of the shape functions in Eq. (15), and "displa" is the nodal displacement vector in each element. In this manner, a high computational efficiency can also be achieved on the computers which cannot perform a parallel computing. For the data post-processing, the phase field, displacement field and mesh data at each loading step are written in the VTK file format according to the format requirements, which can be viewed by using ParaView. The open-source, multi-platform data analysis and visualization software ParaView is easy to learn and master. Using the "Hover points on" icon, we can optionally check the data at each node. The animation view, color maps for pseudo-coloring datasets and crack pattern can be achieved by using the "Animation View", "Choose Preset" and "Save Screenshot Options" icons, respectively. In this way, a convenient and efficient data post-processing can be realized. The values of the phase field variable are obtained based on the element nodes, thus the mesh is refined in the area where the crack is expected to extend. The linear triangular and quadrilateral elements are used in this work. However, the disadvantage of the low computational accuracy of the triangular elements can be improved by the local mesh refinement. Moreover, the triangular elements are more suitable to simulate the crack propagation of inclined cracks and require lower computational cost than the quadrilateral elements. As the staggered algorithm is extremely robust, the present phase field method is a convenient and efficient numerical tool for the simulation of the crack propagation problems.

X.X. Lu et al. / Theoretical & Applied Mechanics Letters 9 (2019) 339-352

4     Numerical examples

method can provide a higher resolution. The force-displacement curves computed by both methods show a slight difference when the force reaches the peak value, which is presumably caused by using different meshing in the calculations. Then, we use the same specimen to simulate the crack propagation problems with different shear loading angles α between the positive x -axis and the loading direction as shown in Fig. 6(a). The angle β between the tangent line of the crackpath at the crack-tip and the positive horizontal direction is defined as the propagation direction of the crack. The results are illustrated in Fig. 6, which show good agreement with the results of Molnár et al. [27] and Bourdin et al. [17].

4.1    Model verification In the first numerical example, the main purpose is to validate the phase field modeling as described in the previous sections. We consider a square plate with a horizontal rectangular notch as a crack. In order to capture the crack trajectory properly, the mesh is refined in the area where the crack is expected to extend. Following Ref. [23], the maximum element size in this refined zone is chosen as one half of the length scale parameter l . The minimum element size is h ≈ 0.001 mm in the critical zone. The width of the notch is also taken as one half of the length scale parameter l . If not explicitly stated otherwise, the following numerical examples are assumed to be in plane strain, homogeneous, isotropic and linear elastic with the material parameters: λ=121.15 kN/mm2, µ=80.77 kN/mm2, G c =2.7 × 10−3 kN/mm , and the length scale parameter is chosen as l = 0.0075 mm [22-24]. The geometry and the boundary conditions of the considered problem are depicted in Fig. 2. For the tension load as shown in Fig. 2(a), the displacement on the lower boundary (y=0) of the specimen along the y -direction is fixed, while a uniform displacement load on the upper boundary is increased continuously. The displacements along the x -direction on both upper and lower boundaries are free. However, the node (x=0, y =0) is fixed to prevent the rigid body motion. The boundary conditions for the shear deformation as depicted Fig. 2(b) are given as follows: the bottom side of the plate is fixed and the top side is subjected to a uniform displacement load ∆u along the x-direction, while the displacement along the y -direction is fixed. Besides, the vertical displacements of the left side (x=0) and right side (x=1) are set to zero. In the tension loading case, a loading increment ∆u = 10−5 mm is applied for 500 steps, then the uniform loading increment is changed to ∆u = 10−6 mm for the rest steps of the simulation. In the shear loading case, the displacement loading increment ∆u = 10−5 mm is applied for 1500 steps. As shown in Figs. 3–5, the crack patterns and force-displacement curves predicted by the present phase field method are in excellent agreement with the results in Ref. [22]. By a closer look at the crack pattern it can be seen that the present phase field  

a

y

4.2    Single edge-cracked specimen with different lengths of the initial crack In this subsection, the influence of the initial crack-length on the crack propagation is studied. The geometry of the specimen with different lengths (L=0.3 mm, 0.4 mm, 0.5 mm, 0.6 mm, and 0.7 mm, respectively) of the initial crack parallel to the x-axis is shown in Fig. 7. The vertical coordinate of the initial crack is constant as y =0.5 mm. The same material parameters, length scale parameter and boundary conditions as described in the Sect. 4.1 are used. 4.2.1    Pure tension loading In the tension loading case, the loading increment ∆u = 10−5 mm is applied for 500 steps, then it is changed to ∆u = 10−6 mm for the rest steps of the simulation. Figure 8(a) and

8(b) show the crack patterns and force-displacement curves. In all considered cases, the crack propagates nearly along the horizontal direction. The peak force and displacement decrease when the length of the initial crack increases. However, the respective peak forces have almost the same displacement value. 4.2.2    Pure shear loading In

this

case,

the

b

y u

0.5

0.5

0.5

0.5

0.5

loading

increment

plied displacement loading are depicted in Fig. 7(b). Figure 9(a) and 9(b) illustrate the crack patterns and force-displacement curves, respectively. Here, we can observe the curved crack paths initiating with almost the same propagation direction β.

u

0.5

displacement

∆u =10−5 mm is applied for 2000 steps. The geometry and the ap-

x

 

343

x 0.5

0.5

Fig. 2.   Schematic of a single edge-cracked specimen: a pure tension loading, b pure shear loading.

344

X.X. Lu et al. / Theoretical & Applied Mechanics Letters 9 (2019) 339-352

 

a

1.00

b

0.75

0.50

0.25

0

 

Fig. 3.   Pure tension loading: a crack pattern of the phase field simulation in Ref. [22], b crack pattern of the phase field simulation in this paper.  

a

b

1.00

0.75

0.50

0.25

0

 

Fig. 4.   Pure shear loading: a crack pattern of the phase field simulation in Ref. [22], b crack pattern of the phase field simulation in this paper.  

a

b 0.6

0.8 0.7

0.5

0.5

Force F (kN)

Force F (kN)

0.6

0.4 0.3 0.2

Miehe et al. [22] This work

0.1 0

 

0

0.002 0.004 0.006 Displacement u (mm)

0.4 0.3 0.2

Miehe et al. [22] This work

0.1

0.008

0

0

0.015 0.005 0.010 Displacement u (mm)

0.020

Fig. 5.   Force-displacement curves compared with the results in Ref. [22]. a Pure tension loading, b pure shear loading.

The force decreases first and then increases after the first peak value. For each curve in Fig. 9(b), the first circle-mark indicates that the crack has been initiated and started to propagate, and the second circle-mark designates that the crack has been extended to the bottom or right side of the specimen. Moreover, as the length of the initial crack increases, the peak force gradually changes from large to small.

4.3    Single edge-cracked specimen with different locations of the initial crack In this subsection, the influence of the initial crack location on the crack propagation is studied. The geometry of the edgecracked specimen with different y -coordinates (y=0.2 mm, 0.3 mm, 0.4 mm, 0.5 mm, 0.6 mm, 0.7 mm, and 0.8 mm, respect-

X.X. Lu et al. / Theoretical & Applied Mechanics Letters 9 (2019) 339-352  

a

α

1.0

u b 70

0.9

60

y (mm)

0.7

Direction of crack β (°)

0.8 β

0.6 0.5 0.4 0.3 0.2

 

50 40 30 Molnar et al. [27] Bourdin et al. [17] This work Linear fit

20 10

0.1 0

345

0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 x (mm)

0

10

20 30 40 50 60 70 Angle of shear loading α (°)

80

90

Fig. 6.   a Crack patterns with different angles of shear loading α, b propagation direction of the crack β versus the angle of shear loading α. a

 

y

b

u

y u

L 1.0

L

1.0 y

y x

x

1.0

 

1.0

Fig. 7.   Schematic of a single edge-cracked specimen with different lengths of the initial crack. a Pure tension loading, b pure shear loading. a

1.0 0.9 0.8

y (mm)

0.7

b 1.0 L = 0.3

L = 0.4

0.5 0.4 L = 0.5

L = 0.6

L = 0.7

0.7 0.6 0.5 0.4 0.3

0.2

0.2

0.1

0.1

0

L = 0.3 L = 0.4 L = 0.5 L = 0.6 L = 0.7

0.9 0.8

0.6

0.3

 

L = 0.3 L = 0.4 L = 0.5 L = 0.6 L = 0.7

Force F (kN)

 

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 x (mm)

0

0

1

2 3 4 4 Displacement u (mm)

6

×10−3 7

Fig. 8.   Single edge-cracked specimen with different lengths of the initial crack under pure tension loading. a Crack patterns of the phase field simulation, b force-displacement curves.

ively) of the initial crack parallel to the x-axis is shown in Fig. 7. The length of the initial crack is taken as L = 0.2 mm.

4.3.1    Pure tension loading Figure 10(a) and 10(b) show the crack patterns and force-dis-

346  

X.X. Lu et al. / Theoretical & Applied Mechanics Letters 9 (2019) 339-352 a

1.0

L = 0.3 L = 0.4 L = 0.5 L = 0.6 L = 0.7

0.9 0.8

L = 0.4

L = 0.6

L = 0.7

b 0.9

0.6 0.5 0.4 0.3

0.7

L = 0.5

0.2

0.6 0.5 0.4 0.3 0.2

0.1

0.1

0

0

 

L = 0.3 L = 0.4 L = 0.5 L = 0.6 L = 0.7

0.8

Force F (kN)

y (mm)

0.7

L = 0.3

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 x (mm)

0

×10−3 0.020

0.005 0.010 0.015 Displacement u (mm)

Fig. 9.   Single edge-cracked specimen with different lengths of the initial crack under pure shear loading. a Crack patterns of the phase field simulation, b force-displacement curves.  

a

1.0

y = 0.2 y = 0.3 y = 0.4

0.9 0.8

y = 0.5 y = 0.6 y = 0.7

y = 0.8

y

0.6 0.5 0.4 0.3 0.2 0.1 0

 

y = 0.2 y = 0.3 y = 0.4 y = 0.5 y = 0.6 y = 0.7 y = 0.8

Force F (kN)

0.7

b 1.2 1.1 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 x

0

1

2 3 4 5 Displacement u (mm)

6

7

×10−3

Fig. 10.   Single edge-cracked specimen with different y-coordinates of the initial crack under pure tension loading. a Crack patterns of the phase field simulation, b force-displacement curves.

placement curves. We can see a symmetrical propagation tendency of the cracks towards the axis of symmetry (y=0.5) when the initial cracks are symmetric with respect to the symmetry axis, and they have almost the same force-displacement curve. The farther the initial crack away from the symmetry axis, the larger the peak force is. 4.3.2    Pure shear loading For the pure shear loading case, the displacement loading increment ∆u = 10−5 mm is also applied for 2000 steps. The crack patterns and force-displacement curves are illustrated in Fig. 11(a) and 11(b), respectively. The propagation tendencies of the cracks are very similar when a pure shear loading is applied. The force-displacement curves in Fig. 11(b) have nearly the same first peak value and the same corresponding displacement as in the case when the initial cracks are symmetric with respect to the axis of symmetry (y=0.5). An initial crack farther away from the symmetry axis shows a larger first peak force, which is similar to the tension loading case. The force decreases first and then increases after the first peak value. As the distance of the initial crack to the bottom side (y=0) increases, the force becomes la-

ger gradually. The reason is that the initial crack closer to the bottom side has a shorter crack trajectory. When the crack is closer to the bottom side, the constraint of the left bottom side below the crack decreases gradually, but it keeps the same constraint at the right bottom side. These conditions will make the specimen can bear a larger force. 4.4    Single edge-cracked specimen with different inclination angles of the initial crack This subsection investigates the propagation of an inclined crack in the same square plate as considered in the previous examples. The geometry of the plate and the inclination angles of the initial crack (θ=15°, 30°, 45°, 60°, and 75°, respectively) are shown in Fig. 12. The crack-length is taken as L=0.2 mm. 4.4.1    Pure tension loading For the pure tension loading as depicted in Fig. 12(a), the crack patterns and the force-displacement curves are given in Fig. 13(a) and 13(b) , respectively. The smaller the inclination angle of the initial crack, the smaller the corresponding peak

X.X. Lu et al. / Theoretical & Applied Mechanics Letters 9 (2019) 339-352

force. But in all considered cases, the cracks propagate nearly toward the horizontal direction parallel to the top and bottom boundaries of the specimen. a

1.0

y = 0.2 y = 0.3 y = 0.4 y = 0.5 y = 0.6 y = 0.7 y = 0.8

0.9 0.8

y (mm)

0.7 0.6 0.5

y = 0.2

y = 0.4

0.3

b 1.4

y = 0.5

y = 0.7

1.0 0.8 0.6

0.2

0.1

 

y = 0.2 y = 0.3 y = 0.4 y = 0.5 y = 0.6 y = 0.7 y = 0.8

0.4

y = 0.8

0.2 0

We also consider an equivalent horizontal crack, which is the projection of the inclined crack on the horizontal axis as shown in Fig. 12(a). The processes of the crack propagation for differ-

1.2

y = 0.6

0.4

y = 0.3

Force F (kN)

 

347

0

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 x (mm)

0

0.005 0.010 0.015 Displacement u (mm)

0.02

Fig. 11.   Single edge-cracked specimen with different y-coordinates of the initial crack under pure shear loading. a Crack patterns of the phase field simulation, b force-displacement curves. a

 

b

y

y

u

0.5

u

0.5

θ

θ

Equivalent notch L 0.5

0.5

L x

x 1.0

1.0

 

Fig. 12.   Schematic of the single edge-cracked specimen with different inclination angles of the initial crack. a Pure tension loading, b pure shear loading. a

1.0 0.9 0.8

y (mm)

0.7

θ = 15° θ = 30 ° θ = 45° θ = 60° θ = 75°

θ = 15°

0.5 0.4

0

 

1.4 1.2 1.0 0.8 0.6

0.3 0.1

θ = 15° θ = 30 ° θ = 45° θ = 60° θ = 75°

1.8 1.6

0.6

0.2

b 2.0

Force F (kN)

 

θ = 30 °

θ = 45°

θ = 60°

θ = 75°

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 x (mm)

0.4 0.2 0

0

1

2

3 4 5 6 7 Displacement u (mm)

8

9

×10−3 10

Fig. 13.   Single edge-cracked specimen with different inclination angles of the initial crack under pure tension loading. a Crack patterns of the phase field simulation, b force-displacement curves.

348

X.X. Lu et al. / Theoretical & Applied Mechanics Letters 9 (2019) 339-352

ent inclination angles (θ = 15◦, 30°, 45°, 60° and 75°, respectively) of the initial crack and the corresponding equivalent horizontal crack are simulated, respectively. Figure 14(a) and 14(b) show that the inclined crack has a very similar crack propagation pattern and force-displacement curve as the corresponding equivalent horizontal crack. But as the inclination angle of the initial crack increases, the inclined crack shows a larger peak force. The crack propagation patterns of the equivalent cracks are always slightly above that of the corresponding inclined cracks. In fact, the tip of the equivalent horizontal cracks (notches) is always slightly above the tip of the initial inclined cracks (notches) as illustrated in Fig. 12(a).

when the angle is larger than 45°. Besides, all the inclined initial cracks show the same force-displacement curves at the beginning, as shown in Fig. 15(b). With the increase of the displacement loading ∆u , the force gradually increases. It is noted here that when the displacement loading ∆u reaches a certain level, the force suddenly drops. When the inclination angle of the initial crack is 45°, the maximum peak force appears. 4.5    Single edge-cracked bi-material specimen with different y-coordinates of the initial crack Next, single edge-cracked bi-material specimens made of two different materials with different y-coordinates (y=0.10 mm, 0.20 mm, 0.25 mm, 0.30 mm, 0.40 mm, and 0.48 mm, respectively) of the initial crack in the material B are studied. The length of the initial crack is given by L=0.2 mm. We assume that the interface is perfectly bonded. The geometry and the applied displacement loading are depicted in Fig. 16. The upper material has better material properties (E a = 210 GPa, νa = 0.3 , G ca = 2.7 × 10−3 kN/mm ), while the lower material has poor ma-

4.4.2    Pure shear loading After the tension loading case we consider the shear loading case as illustrated in Fig. 12(b). In Fig. 15(a), the simulated crack patterns are depicted. The crack propagates from the right-upper tip of the initial crack when the angle is less than 45°. The starting point of the crack propagation is in the lower-middle part of the inclined crack for 45°. However, the starting point of the crack growth changes to the left mouth of the inclined crack a

1.0

θ = 15° θ = 30° θ = 45° θ = 60° θ = 75°

0.9 0.8

y (mm)

0.7

b 2.0

Equivalent crack for θ = 15° Equivalent crack for θ = 30° Equivalent crack for θ = 45° Equivalent crack for θ = 60° Equivalent crack for θ = 75°

1.6

0.6 0.5 0.4 0.3

θ = 15°E

θ = 30°E

θ = 45°E

θ = 60°E

θ = 75°E

0.1

 

1.4 1.2 1.0 0.8 0.6 0.4

0.2 0

θ = 15° θ = 15° equivalent crack θ = 30° θ = 30° equivalent crack θ = 45° θ = 45° equivalent crack θ = 60° θ = 60° equivalent crack θ = 75° θ = 75° equivalent crack

1.8

Force F (kN)

 

θ = 15°

θ = 30°

θ = 45°

θ = 60°

θ = 75°

0.2 0

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 x (mm)

0

1

2

3 4 5 6 7 Displacement u (mm)

8

9

×10−3 10

Fig. 14.   Single edge-cracked specimen with different inclination angles of the initial crack and the equivalent horizontal crack under pure tension loading. a Crack patterns of the phase field simulation, b force-displacement curves. a

1.0 0.9 0.8

y (mm)

0.7 0.6

b 2.0 θ = 15°

θ = 30°

0.4

θ = 45°

θ = 60°

θ = 75°

1.4 1.2 1.0 0.8 0.6

0.2

0.4

0.1

0.2

0

θ = 15° θ = 30° θ = 45° θ = 60° θ = 75°

1.8 1.6

0.5 0.3

 

θ = 15° θ = 30° θ = 45° θ = 60° θ = 75°

Force F (kN)

 

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 x (mm)

0

0

0.005

0.010 0.015 0.020 Displacement u (mm)

0.025

0.030

Fig. 15.   Single edge-cracked specimen with different inclination angles of the initial crack under pure shear loading. a Crack patterns of the phase field simulation, b force-displacement curves.

X.X. Lu et al. / Theoretical & Applied Mechanics Letters 9 (2019) 339-352

terial properties (E b = 21 GPa, νb = 0.3, G cb = 2.7 × 10−4 kN/mm ), and the length-scale parameter is also taken as l = 0.075 mm .

∆u = 10−5 mm is applied for 1500 steps. Figure 18(a) shows the

simulated crack patterns for different y-coordinates of the initial crack. They have almost the same propagation direction β as shown in Fig. 11(a) for the crack propagation in a single material specimen. The force-displacement curves are depicted in Fig. 18(b). The slope of each curve decreases first and then increases, but the slope change for the initial crack closer to the bottom side (y=0) is not obvious. The global trends of the force-displacement curves are very similar to that of the crack propagation in a single material specimen as shown in Fig. 11(b): the closer is the initial crack to the bottom side, the larger the ultimate force will be.

4.5.1    Pure tension loading The

displacement

loading

increment

is

taken

as

∆u = 10−5 mm for the first 350 steps, and then ∆u = 10−6 mm for

the following simulation. Figure 17(a) and 17(b) illustrate the simulated crack patterns and force-displacement curves for different y -coordinates of the initial crack. As we can see in Fig. 17(a), the crack far away from the interface propagates towards the interface, while the crack close to the interface extends away from the interface. However, for some special crack locations the crack propagates horizontally. The peak force gradually increases when the initial crack is farther away from the interface, as shown in Fig. 17(b). Compared to the crack propagation in a single material specimen as illustrated in Fig. 10, different phenomena regarding the crack growth path and the force-displacement curve can be observed in Fig. 17. These differences are caused by the interface and different material properties.

5     Conclusions In this paper, a phase field method based on Miehe's approach is applied for the crack propagation simulation in 2D brittle elastic materials. The method can be used to simulate the crack initiation, propagation, kinking, branching and so on. The simulation results are in good agreement with the reference results [22]. Several numerical examples are presented and dis-

4.5.2    Pure shear loading In

this

case,

the a

 

displacement

increment b

y

y

u

u

Material A

0.5

0.5

loading

Material A

0.5

Material B

L

0.5

y

L

Material B y

x

x

1.0

 

349

1.0

Fig. 16.   Schematic of the single edge-cracked bi-material specimen with different y-coordinates of the initial crack. a Pure tension loading, b pure shear loading. a

1.0 0.9 0.8

y (mm)

0.7 0.6

y = 0.10 y = 0.20 y = 0.25 y = 0.30 y = 0.40 y = 0.48

y = 0.10 y = 0.20 y = 0.25

y = 0.30 y = 0.40 y = 0.48

0.5 0.4 0.3

θ = 75°

0.2

 

y = 0.10 y = 0.20 y = 0.25 y = 0.30 y = 0.40 y = 0.48

0.14 0.12 0.10 0.08 0.06 0.04 0.02

0.1 0

b 0.16

Force F (kN)

 

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 x (mm)

0

0

×10−3 1.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 Displacement u (mm)

Fig. 17.   Single edge-cracked bi-material specimen with different y-coordinates of the initial crack under pure tension loading. a Crack patterns of the phase field simulation, b force-displacement curves.

 

X.X. Lu et al. / Theoretical & Applied Mechanics Letters 9 (2019) 339-352 a

1.0 0.9 0.8

y (mm)

0.7 0.6

y = 0.10 y = 0.20 y = 0.25 y = 0.30 y = 0.40 y = 0.48

y = 0.10

y = 0.20

y = 0.25

y = 0.30

y = 0.40

y = 0.48

0.3

0.12 0.10 0.08 0.06

0.2

0.04

0.1

0.02

0

y = 0.10 y = 0.20 y = 0.25 y = 0.30 y = 0.40 y = 0.48

0.16 0.14

0.5 0.4

b 0.18

Force F (kN)

350

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 x (mm)

0

0

0.003

 

0.006 0.009 0.012 Displacement u (mm)

0.015

Fig. 18.   Single edge-cracked bi-material specimen with different y-coordinates of the initial crack under pure shear loading. a Crack patterns of the phase field simulation, b force-displacement curves.

cussed to show the effects of different geometrical and material parameters on the crack propagation paths and force-displacement curves. The following conclusions can be drawn from the present investigation. (1) In the pure tension loading case, initial horizontal cracks of different lengths located on the horizontal symmetry-axis (y=0.5) of the specimen propagate straightforward in the horizontal direction, and the peak force and maximum displacement decrease when the length of the initial crack increases. The horizontal cracks located symmetrically with respect to the symmetry-axis show a propagation tendency toward the symmetryaxis, and they exhibit almost the same force-displacement curves for the same initial crack-length. Moreover, the farther is the initial crack away from the symmetry-axis, the larger the peak force is. For inclined cracks of the same length, as the crack inclination angle increases, the peak force increases. They have almost the same crack trajectory and force-displacement curves as their equivalent horizontal cracks. But when the angle is over 45°, the peak force for the equivalent crack is smaller than that for the original inclined crack. (2) When the pure shear loading is applied, initial horizontal cracks of different lengths located on the horizontal symmetryaxis (y=0.5) of the specimen have the same propagation trend. The horizontal cracks with different y -coordinates also have a similar propagation trend, and the peak forces are the same if the cracks are located symmetrically with respect to the symmetry-axis. The ultimate force is larger for the crack closer to the bottom boundary. On the other hand, the crack propagation tendencies as well as the force-displacement curves are different for the inclined cracks with different inclination angles. The peak force is maximum and minimum at 45° and 0°, respectively. The crack propagates from the tip of the initial crack when the inclination angle is less than 45°. However, when the angle is larger than 45°, the crack propagates from the mouth of the initial crack. While the starting point of the crack propagation is in the middle part of the initial crack for 45°. (3) For the bi-material plate containing a horizontal crack with different y -coordinates under a tension loading, an initial crack far away from the interface propagates towards the interface, while a crack close to the interface extends away from the

interface. For some special crack locations, the crack propagates horizontally. However, an axis of symmetry for the crack paths does not exist. For the shear loading case, the results are quite similar to that for an initial crack in a single material specimen. The results presented in this paper can be used to estimate or prevent possible crack propagations in brittle elastic materials, which are significant in practical engineering applications. As the phase field model based on Miehe's approach is somehow sensitive to the length-scale (regularization parameter l and element-size h), the length-scale insensitive phase field model developed by Wu [28, 42 ] should be implemented in our future work. Acknowledgments This work was supported by the National Natural Science Foundation of China (Grant U1833116). The authors would like to express their sincere thanks to Professor Jian-Ying Wu (South China University of Technology, Guangzhou, China), Dr. Tianxue Ma and Mr. Zhen Yan (University of Siegen, Siegen, Germany) for the helpful discussions in the course of this study. The first author is also grateful to the financial support by the China Scholarship Council (CSC) for the joint PhD scholarship at the Chair of Structural Mechanics, Department of Civil Engineering, University of Siegen, Germany.

References [1] A.A. Griffith, VI. The phenomena of rupture and flow in solids, Philosophical Transactions of the Royal Society London 221 (1921) 163–198. [2] G.R. Irwin, Fracture, In: Handbuch der physik, Band VI, Elastizität und Plastizität, Springer-Verlag, Berlin (1958). (in German) [3] A.R. Shahani, M.R. Amini Fasakhodi, Finite element analysis of dynamic crack propagation using remeshing technique, Materials and Design 30 (2009) 1032–1041. [4] R. Rangarajan, M.M. Chiaramonte, M.J. Hunsweck, et al., Simulating curvilinear crack propagation in two dimensions with universal meshes, International Journal for Numerical Meth-

X.X. Lu et al. / Theoretical & Applied Mechanics Letters 9 (2019) 339-352

ods in Engineering 102 (2015) 632–670.

[5] G.L. Peng, Y.H. Wang, A node split method for crack growth

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15] [16]

[17]

[18]

[19]

[20] [21]

[22]

problem, Applied Mechanics and Materials 182-183 (2012) 1524–1528. X.P. Xu, A. Needleman, Numerical simulations of fast crack growth in brittle solids, Journal of the Mechanics and Physics of Solids 42 (1994) 1397–1434. F. Zhou, J.F. Molinari, Dynamic crack propagation with cohesive elements: A methodology to address mesh dependency, International Journal for Numerical Methods in Engineering 59 (2004) 1–24. Z.J. Yang, A.J. Deeks, Modelling cohesive crack growth using a two-step finite element-scaled boundary finite element coupled method, International Journal of Fracture 143 (2007) 333–354. B.D. Le, G. Koval, C. Chazallon, Discrete element model for crack propagation in brittle materials, International Journal for Numerical and Analytical Methods in Geomechanics 40 (2016) 583–595. N. Moës, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, International Journal for Numerical Methods in Engineering 46 (1999) 131–150. C. Daux, N. Moes, J. Dolbow, et al., Arbitrary branched and intersecting cracks with the extended finite element method, International Journal for Numerical Methods in Engineering 48 (2000) 1741–1760. N. Moës, C. Stolz, P.E. Bernard, et al., A level set based model for damage growth: The thick level set approach, International Journal for Numerical Methods in Engineering 86 (2011) 358–380. E. Gürses, C. Miehe, A computational framework of three-dimensional configurational-force-driven brittle crack propagation, Computer Methods in Applied Mechanics and Engineering 198 (2009) 1413–1428. G.A. Francfort, J.J. Marigo, Revisiting brittle fracture as an energy minimization problem, Journal of the Mechanics and Physics of Solids 46 (1998) 1319–1342. M. Buliga, Energy minimizing brittle crack propagation, Journal of Elasticity 52 (1999) 201–238. G. Dal Maso, R. Toader, A Model for the quasi-static growth of brittle fractures: Existence and approximation results, Archive for Rational Mechanics and Analysis 162 (2002) 101–135. B. Bourdin, G.A. Francfort, J.J. Marigo, Numerical experiments in revisited brittle fracture, Journal of the Mechanics and Physics of Solids 48 (2000) 797–826. B. Bourdin, C.J. Larsen, C.L. Richardson, A time-discrete model for dynamic fracture based on crack regularization, International Journal of Fracture 168 (2011) 133–143. H. Amor, J.-J. Marigo, C. Maurini, Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments, Journal of the Mechanics and Physics of Solids 57 (2009) 1209–1229. B. Bourdin, G.A. Francfort, J.-J. Marigo, The variational approach to fracture, Journal of Elasticity 91 (2008) 5–148. D. Mumford, J. Shah, Optimal approximations by piecewise smooth functions and associated variational problems, Communications on Pure and Applied Mathematics 42 (1989) 577–685. C. Miehe, M. Hofacker, F. Welschinger, A phase field model for rate-independent crack propagation: Robust algorithmic im-

[23]

[24]

[25]

[26]

[27]

[28]

[29]

[30]

[31]

[32]

[33]

[34]

[35]

[36]

[37]

351

plementation based on operator splits, Computer Methods in Applied Mechanics and Engineering 199 (2010) 2765–2778. C. Miehe, F. Welschinger, M. Hofacker, Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations, International Journal for Numerical Methods in Engineering 83 (2010) 1273–1311. T.T. Nguyen, J. Yvonnet, Q.Z. Zhu, et al., A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure, Engineering Fracture Mechanics 139 (2015) 18–39. N. Singh, C.V. Verhoosel, R. De Borst, et al., A fracture-controlled path-following technique for phase-field modeling of brittle fracture, Finite Elements in Analysis and Design 113 (2016) 14–29. T.T. Nguyen, J. Yvonnet, M. Bornert, et al, On the choice of parameters in the phase field method for simulating crack initiation with experimental validation, International Journal of Fracture 197 (2016) 213–226. G. Molnár, A. Gravouil, 2D and 3D Abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture, Finite Elements in Analysis and Design 130 (2017) 27–38. J.-Y. Wu, A unified phase-field theory for the mechanics of damage and quasi-brittle failure, Journal of the Mechanics and Physics of Solids 103 (2017) 72–99. P.F. Li, Q.Z. Zhu, S.T. Gu, et al., A phase field method to simulate crack nucleation and crack propagation in rock-like materials, Engineering Mechanics 35 (2018) 41–48. M. Hofacker, C. Miehe, Continuum phase field modeling of dynamic fracture: Variational principles and staggered FE implementation, International Journal of Fracture 178 (2012) 113–129. M.J. Borden, C.V. Verhoosel, M.A. Scott, et al., A phase-field description of dynamic brittle fracture, Computer Methods in Applied Mechanics and Engineering 217-220 (2012) 77–95. T. Li, J.-J. Marigo, D. Guilbaud, et al., Gradient damage modeling of brittle fracture in an explicit dynamics context, International Journal for Numerical Methods in Engineering 108 (2016) 1381–1405. V.P. Nguyen, J.-Y. Wu, Modeling dynamic fracture of solids with a phase-field regularized cohesive zone model, Computer Methods in Applied Mechanics and Engineering 340 (2018) 1000–1022. C. Hesch, K. Weinberg, Thermodynamically consistent algorithms for a finite-deformation phase-field approach to fracture, International Journal for Numerical Methods in Engineering 99 (2014) 906–924. C. Hesch, A.J. Gil, R. Ortigosa, et al., A framework for polyconvex large strain phase-field methods to fracture, Computer Methods in Applied Mechanics and Engineering 317 (2017) 649–683. C. Miehe, L.-M. Schänzel, H. Ulmer, Phase field modeling of fracture in multi-physics problems. Part I. Balance of crack surface and failure criteria for brittle crack propagation in thermoelastic solids, Computer Methods in Applied Mechanics and Engineering 294 (2015) 449–485. C. Miehe, M. Hofacker, L.-M. Schänzel, et al., Phase field modeling of fracture in multi-physics problems. Part II. Coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic-plastic solids, Computer Methods in Applied Mechanics and Engineering 294 (2015) 486–522.

352

X.X. Lu et al. / Theoretical & Applied Mechanics Letters 9 (2019) 339-352

[38] C. Miehe, L.-M. Schänzel, Phase field modeling of fracture in

[43] S. Teichtmeister, D. Kienle, F. Aldakheel, et al., Phase field

rubbery polymers. Part I: Finite elasticity coupled with brittle failure, Journal of the Mechanics and Physics of Solids 65 (2014) 93–113. T. Heister, M.F. Wheeler, T. Wick, A primal-dual active set method and predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field approach, Computer Methods in Applied Mechanics and Engineering 290 (2015) 466–495. C. Miehe, F. Aldakheel, A. Raina, Phase field modeling of ductile fracture at finite strains: A variational gradient-extended plasticity-damage theory, International Journal of Plasticity 84 (2016) 1–32. L. Xia, J. Yvonnet, S. Ghabezloo, Phase field modeling of hydraulic fracturing with interfacial damage in highly heterogeneous fluid-saturated porous media, Engineering Fracture Mechanics 186 (2017) 158–180. J.-Y. Wu, V.P. Nguyen, A length scale insensitive phase-field damage model for brittle fracture, Journal of the Mechanics and Physics of Solids 119 (2018) 20–42.

modeling of fracture in anisotropic brittle solids, International Journal of Non-Linear Mechanics 97 (2017) 1–21. T.T. Nguyen, J. Réthoré, M.-C. Baietto, Phase field modelling of anisotropic crack propagation, European Journal of Mechanics - A/Solids 65 (2017) 279–288. Z.-J. Yang, B.-B. Li , J.-Y. Wu, X-ray computed tomography images based phase-field modeling of mesoscopic failure in concrete, Engineering Fracture Mechanics 208 (2019) 151–170. K. Pham, H. Amor, J.-J. Marigo, et al., Gradient damage models and their use to approximate brittle fracture, International Journal of Damage Mechanics 20 (2011) 618–652. C. Miehe, M. Lambrecht, Algorithms for computation of stresses and elasticity moduli in terms of Seth-Hill's family of generalized strain tensors, Communications in Numerical Methods in Engineering 17 (2001) 337–353. C. Miehe, Comparison of two algorithms for the computation of fourth-order isotropic tensor functions, Computers & Structures 66 (1998) 37–43.

[39]

[40]

[41]

[42]

[44]

[45]

[46]

[47]

[48]