- Email: [email protected]
A new phase field fracture model for brittle materials that accounts for elastic anisotropy
Available online at www.sciencedirect.com
ScienceDirect Comput. Methods Appl. Mech. Engrg. 358 (2020) 112643 www.elsevier.com/locate/cma
A new phase field fracture model for brittle materials that accounts for elastic anisotropy Shuaifang Zhanga,b , Wen Jiangc , Michael R. Tonksa,b ,∗ b
a University of Florida, Gainesville, FL, United States of America The Pennsylvania State University, State College, PA, United States of America c Idaho National Laboratory, Idaho Falls, ID, United States of America
Received 22 June 2018; received in revised form 10 September 2019; accepted 11 September 2019 Available online xxxx
Highlights • • • •
The new phase field fracture model includes anisotropic elasticity. Compressive stress does not contribute to crack propagation. Single crystal and bicrystal simulations show the proper relationship with crystal orientation. Polycrystal simulations show the impact of texture on the stress–strain curves.
Abstract In this manuscript, a new phase field model of brittle fracture is proposed that accounts for any elastic anisotropy using spectral decomposition of the stress. To verify the model, both Mode I and Mode II fracture simulations were performed to ensure that the correct crack paths were predicted. Next, the fracture of face-centered cubic (FCC) and hexagonal close packed (HCP) single crystals was modeled, comparing the impact of crystal orientation on the stress–strain curve. Third, fracture in an HCP bicrystal was simulated and the impact of crystal orientation on the crack paths and stress–strain curves was studied. Finally, a polycrystal structure was simulated to study the impact of crystallographic texture on transgranular fracture. c 2019 Elsevier B.V. All rights reserved. ⃝ Keywords: Phase field; Fracture mechanics; Anisotropic materials; Brittle fracture; MOOSE
1. Introduction Various methods have been developed to model discrete crack initiation and growth, including cohesive zone models [1] and the extended finite element method [2], and these methods have been developed for isotropic and anisotropic materials [3]. Recently, the phase field method has emerged as a powerful means of modeling discrete crack nucleation and growth. A large part of its value comes from its simplicity; it can predict crack propagation solely by solving a system of partial differential equations (PDEs). This is because the system of PDEs simultaneously predicts crack formation and creates new cracked free surface. ∗ Corresponding author at: University of Florida, Gainesville, FL, United States of America.
E-mail addresses: [email protected] (S. Zhang), [email protected] (W. Jiang), [email protected] (M.R. Tonks). https://doi.org/10.1016/j.cma.2019.112643 c 2019 Elsevier B.V. All rights reserved. 0045-7825/⃝
2
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
The foundations for phase field methods to model crack propagation in brittle materials was established in the late 1990s, with the variational formulation of brittle fracture by Francfort [4] and Bourdin [5]. In these original models, the quasi-static brittle fracture process is governed by the minimization of the total energy, which is a function of strain energy and crack dissipation. Aranson et al. [6] presented a phase field fracture model to solve mode I fracture based on Ginzburg–Landau phase transition. Karma and collaborators [7,8] modified the strain energy density function from Aranson et al. [6] to simulate mode III fracture. Henry et al. [9] modified the energy function of Aranson et al. [6] to work for mode I and II fracture. Kuhn and Muller [10] were the first to propose the thermodynamics-based phase field model. In this model, a continuous variable field was used to represent cracked material. The system strain energy density was determined by solving the stress divergence equation and the evolution of the crack variable field was defined by a partial differential equation that minimized the energy. Amor [11] separated the strain energy density into two parts, one representing the tensile part of the strain energy, and one representing the compressive part, which did not contribute to the crack propagation. Miehe and collaborators [12,13] proposed a phase field fracture model based on thermodynamic consistency and continuum mechanics in which the strain energy was also decomposed into the tensile and compressive parts. Borden [14,15] developed a higher order equation for the evolution of the phase field parameter that required advanced finite element formulations to solve. Such phase field fracture models have been applied to model the impact of microstructural features on fracture behavior [16,17]. Initial phase field fracture models were developed assuming isotropic material behavior, however newer models have focused on anisotropic materials. With respect to brittle fracture, materials are anisotropic in two ways: anisotropic elastic response and anisotropic surface energy. Hakim and Karma (2005) [18] derived a force balance condition that was used to account for anisotropic surface energy, in which the path is determined by the directional fracture energy, but not elastic anisotropy. Abdollahi and Arias [19] also developed a model for anisotropic surface energy but not anisotropic elasticity, but in ferroelectrics under indentation with electro-mechanical coupling. Bin Li et al. [20] proposed a phase field fracture model that introduced a strongly anisotropic surface energy and anisotropic elasticity, though the contribution of compressive stress to crack propagation is not avoided. Clayton and Knap developed models with anisotropic surface energy using an isotropic elasticity tensor [21] and then an anisotropic elasticity tensor [22]; however, they also did not avoid the contribution of compressive stress to crack propagation. Teichtmeister and Miehe [23] employed a fourth-order structural tensor that was specific for cubic anisotropy. Shanthraj [24] developed a finite-strain viscoplastic phase field method that considered anisotropic surface energy, but not elastic anisotropy. Nguyen et al. [25] also considered surface energy anisotropy but not elastic anisotropy. None of these anisotropic phase field fracture models provide both a means of avoiding the impact of compressive stress on crack propagation and a means of considering arbitrary elastic anisotropy. In this paper, a new phase field fracture model is proposed that employs the full fourth-order elastic tensor, such that any form of anisotropy can be included. The model decomposes the strain energy into compressive and tensile parts and only the tensile stress contributes to crack propagation. The necessary equations are derived and then verified by modeling mode I and II fracture in an isotropic material. Then, the model is demonstrated on anisotropic materials. The model was developed using the finite element method with the open source Multiphysics Object-Oriented Simulation Environment (MOOSE) [26]. 2. Phase field fracture model in anisotropic brittle materials 2.1. Phase field approximation of crack In contrast to traditional finite element methods that are used to simulate discrete cracks in which the cracks are sharp, phase field fracture models employ a diffuse-interface description of cracks. The continuous variable d is used to characterize unbroken material (d = 0) and fully cracked material (d = 1), and the material is considered as damaged when 0 < d < 1. The shape of d across the crack width is described by the differential equation d(x) − l 2
∂ 2d =0 ∂x2
(1)
with solution d(x) = e−|x|/l ,
(2)
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
3
Fig. 1. Phase field parameter as a function of spatial coordinate: d(x) = e−|x|/l .
where l is the length parameter used to represent the crack width, as shown in Fig. 1. The crack approximation converges to a sharp crack as l approaches zero. In multiple dimensions, the differential equation that governs the shape of d is d − l 2 ∆d = 0
(3)
2.2. Crack propagation with a phase field model in anisotropic brittle materials The total potential energy of a cracked brittle material is defined as a sum of elastic energy and crack surface energy [14]: ∫ ∫ gc dΓ , (4) ψe dΩ + Ψtot = Ω
Γ
where ψe is the elastic energy density, gc is the crack surface energy, Ω represents the entire material, and Γ is the crack surface area. As the crack in the phase field model now has a finite thickness, the crack energy can be expressed as a volume integral using the crack surface density function per unit volume from Miehe [13], which is an extension of Bourdin et al. [27]: 1 l γ = d 2 + |∇d|2 , (5) 2l 2 ∫ ∫ such that Γ gc dΓ = Ω gc ( 2l1 d 2 + 2l |∇d|2 )dΩ . Thus, Eq. (4) becomes ) ∫ ( gc 2 gc l 2 ψe + d + Ψtot = |∇d| dΩ . (6) 2l 2 Ω As we are considering only brittle materials, it is reasonable to assume small deformation, so that the strain tensor is defined as ) 1( ∇u + ∇uT . (7) ϵ= 2 We also assume a linear elastic material response, such that the stress tensor of the unbroken material σ 0 = Cϵ,
(8)
where C is the fourth-order elasticity tensor, with as much symmetry as desired. In our model, we divide the elastic energy density into its tensile ψ + and compressive ψ − parts, such that ψe = ((1 − d)2 (1 − k) + k)ψ + + ψ − ,
(9)
where k is a small value that is used to ensure positive definiteness after the material is fully broken (d = 1). The energy is split for two reasons: first, to ensure that only tensile stress contributes to crack propagation; second, to ensure cracks only support compressive stress. We apply spectral decomposition of the unbroken stress tensor to determine values for ψ + and ψ − σ 0 = QΛ Q T ,
(10)
4
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
where Λ = diag(λ1 , λ2 , λ3 , ) is a diagonal matrix containing the three eigenvalues of the stress tensor (or the principal stress tensor) and Q is the eigenvector matrix of the stress. We then define the tensile Λ+ and compressive Λ− principal stress tensors as Λ+ = I+ ΛI+ Λ− = I− ΛI−
(11) (12)
with I+ and I− being similar to the identify matrix, except that Iii+ = 0 when λi < 0 and Iii− = 0 when λi > 0. From the principal stress tensors, we can determine the compressive and tensile stress tensors σ + = QΛ+ QT σ − = QΛ− QT ,
(13) (14)
and thus 1 + σ :ϵ (15) 2 1 ψ − = σ − : ϵ. (16) 2 The primary differences between this model and that from Miehe [13] are that this model includes the full rank four elasticity tensor, while Miehe’s model only includes the two isotropic elastic constants, and that this model directly employs a spectral decomposition of the elastic stress tensor, while Miehe’s model decomposes the elastic stress tensor using a function of the eigenvalues and eigenvectors of the strain tensor and the two isotropic elastic constants. We also define projection tensors that separate the stress into its compressive and tensile parts ψ+ =
σ + = P+ σ 0 σ − = P− σ 0 ,
(17) (18)
where ∂σ + ∂σ 0 − P = I − P+ . P+ =
(19) (20)
These projection tensors are calculated from the spectral decomposition of the unbroken stress tensor using the algorithms defined by Miehe [28] and Miehe and Lambrecht [29]. This relationship involving the projection tensor has been widely used in damage mechanics [30,31] and was also used in phase field fracture by Miehe [12]. 3. Governing equations The fracture model is defined by two equations, the stress equilibrium equation and the evolution equation. In this section, we will derive these two equations. 3.1. Equilibrium equation Eq. (9) defines the strain energy density as a function of the phase field parameter d, which represents the damaged region of the material. The stress tensor in the damaged anisotropic material is defined as the derivative of the strain energy density with respect to the strain tensor, ∂ψe ∂ψ + ∂ψ − = ((1 − d)2 (1 − k) + k) + . ∂ϵ ∂ϵ ∂ϵ The partial derivative of Eq. (15) with respect to the strain is σ =
∂ψ + 1 dσ + = ϵ+ ∂ϵ 2 dϵ + 1 dσ = ϵ+ 2 dϵ
1 + dϵ σ 2 dϵ 1 + σ 2
(21)
(22) (23)
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
5
We take the partial derivative of σ + with respect to the strain and substitute in Eq. (8) to obtain dσ + dσ + dσ 0 dσ + = = C dϵ dσ 0 dϵ dσ 0 and the partial derivative of Eq. (17) with respect to the unbroken stress to get
(24)
dσ + = P+ . (25) dσ 0 The expressions involving the compressive stress are obtained in a similar manner. If we combine these expressions with Eq. (19) we obtain ∂ψ + = σ+ ∂ϵ ∂ψ − = σ −. ∂ϵ Then, we substitute these into Eq. (21) to obtain σ = ((1 − d)2 (1 − k) + k)σ + + σ − .
(26) (27)
(28)
This expression for the stress is then used in the stress equilibrium equation to solve the mechanics problem. 3.2. Phase field evolution equation The evolution of the damage parameter d minimizes the total free energy of the system using an Allen–Cahn equation: ) ( ∂ψe ∂d δΨtot gc = −L = −L + d − ∇ · gc l∇d . (29) ∂t δd ∂d l Combining this expression with Eq. (9) gives ( ) ∂d gc = −L (−2(1 − d)(1 − k) + k)ψ + + d − ∇ · gc l∇d . (30) ∂t l This equation does not guarantee that cracks do not heal. Therefore, we alter it by introducing the history variable H , which is the maximum positive tensile energy experienced during the simulation, H = max ψ + for t = [0, tc ],
(31)
where tc is the current time. This is the same approach taken in Miehe et al. [12]. Thus, the evolution equation becomes ) 1( gc ∂d = − (−2(1 − d)(1 − k) + k)H + d − ∇ · gc l∇d . (32) ∂t η l The expression L = 1/η, with the rate parameter η (called viscosity by Miehe et al. [12]), has been used to obtain the form from Miehe et al. [13]. Note that this equation is identical to the Allen–Cahn equation used in traditional phase field simulations, ( ) ∂d ∂ floc = −L − ∇ · κ∇d , (33) ∂t ∂d with floc = ((1 − d)2 (1 − k) + k)H + ψ − + 1 η κ = gc l.
L =
gc 2 d 2l
(34) (35) (36)
This model has been implemented with the finite element method using MOOSE [32], as was done for isotropic materials by Chakroborty et al. [17]. It is solved using implicit time integration.
6
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
Fig. 2. Geometry and boundary conditions used in the mode I fracture simulations shown in Figs. 3, 7 and 9. Table 1 Material properties used in the Mode I and II fracture simulations of an isotropic material. E (GPa) 120
ν
gc (GPa mm)
η (s/mm)
k
0.3
1.0 × 10−3
1.0 × 10−7
1.0 × 10−6
4. Numerical results To validate that our fracture model adequately represents the fracture of brittle materials, we perform a series of 2D plane strain simulations. First we model mode I and II fracture of an isotropic material (with two unique elastic constants, the elastic modulus E and Poisson’s ratio ν). Then, we model fracture in hexagonal and cubic single crystals, as well as hexagonal bicrystals. Finally, we model fracture in hexagonal polycrystals. 4.1. Mode I fracture in isotropic materials To model mode I fracture in an isotropic material, we used a single edge notch tension specimen shown in Fig. 2. The dimensions and boundary conditions of the specimen are also shown in the figure. The material properties used for the phase field fracture simulations do not represent a specific material and are shown in Table 1. The crack length parameter is l = 0.005 mm. Chakraborty et al. [16] studied the mesh size requirement for phase field fracture simulations and found that the interface was adequately resolved when the smallest element size was equal to half the length parameter (l/2). In this work, we employed mesh adaptivity to reduce computational expense but ensure that the elements across the crack are always l/2 or smaller. Different length parameters were used in the different simulations in order to obtain a good balance between computational cost (which increases as l decreases) and similarity to a sharp crack model (which also generally increases as l decreases). In all the simulations, the time step size dt = 1 × 10−4 s. The final crack path and stress–strain curve for mode I fracture are shown in Figs. 3a and 3b, respectively. The stress field at different displacements of the top boundary are shown in Figs. 3b and 3c. From the stress evolution and crack propagation plots we can see that the crack propagates parallel to the initial crack, consistent with the mode I fracture behavior shown in Miehe et al. [12]. The material behavior under a compressive load was also investigated, where v = −t. The boundary conditions were slightly modified, such that the x-displacements are set to zero on the right side and the y-displacements are set to zero on the bottom. A larger compressive displacement load was applied than was applied in the tensile simulations and yet no fracture occurred, as shown in Fig. 4. This verifies that cracks do not propagate due to compressive stress. 4.2. Mode II fracture test in isotropic materials Shear force induced cracks are also very common in engineering structures, so we conducted simulations with an applied shear load (mode II fracture). The geometry and boundary conditions of the mode II crack simulation
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
7
Fig. 3. Crack propagation in the Mode I fracture simulation of an isotropic material with l = 0.005 mm. The final crack path is shown in (a); the σ yy stress distribution in GPa at displacements of v = 3.65 µm and v = 3.68 µm are shown in (b) and (c), respectively, with regions where d = 1 colored white. The resultant stress–strain curve is shown in (d).
Fig. 4. Material response under a compressive load of the isotropic material with l = 0.005 mm at v = −10 µm, where the σ yy stress distribution in GPa is shown in (a) and the phase field variable field is shown in (b). Note that no crack propagation occurs due to the uniform compressive stress distribution.
are shown in Fig. 5, in which the single notched specimen is used and the bottom boundary is fixed while the top side has an applied x-displacement. The material properties from Table 1 are used; the length parameter l = 0.01 mm. Figs. 6b and 6c show the crack evolution and stress evolution at different x-displacements of the top surface. As shown in Fig. 6a, the crack path of mode II fracture is consistent with classical mode II fracture results and with the results from other phase field fracture models [5,12], proving that our model works well with mode II fracture.
8
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
Fig. 5. Boundary conditions and geometry parameters for the Mode II fracture simulation shown in Fig. 6.
Fig. 6. Crack propagation in the Mode II fracture simulation of the isotropic material, where the final crack path is shown in (a) and the σx x stress distribution in GPa at displacements u = 6.43 µm and u = 8.92 µm are shown in (b) and (c), respectively, with regions where d = 1 colored white. Table 2 Picked material properties of HCP structure. c11 (GPa) 115.8
c12 (GPa) 39.8
c13 (GPa) 40.6
c33 (GPa) 51.4
c44 (GPa)
gc (GPa mm)
η (s/mm)
20.4
1.0 × 10−3
1.0 × 10−6
4.3. Brittle fracture in anisotropic materials 4.3.1. Fracture simulation in single crystals with anisotropic elasticity In this section, we demonstrate the influence of crystal orientation on the fracture behaviors in anisotropic single crystals. For all of these simulations, we employed the geometry and boundary conditions shown in Fig. 2. Firstly, a transversely isotropic material was studied, where hexagonal close-packed (hcp) crystals are the most common example. Second, we modeled a cubic material. There are five unique elastic constants in hcp crystals, and here we used cadmium (Cd) as the example material, with the material properties shown in Table 2. Note that this material is symmetric along the x- and y-directions, so the 2D simulations are done in the x-z plane. To verify that the crystal orientation influences the fracture behavior in hcp Cd, simulations were run for single crystals with different orientations, defined by a rotation angle θ around the out-of-plane axis. Examples of the σ yy stress field and the crack field after a given amount of deformation for two representative crystal orientations are shown Figs. 7a and 7b. The magnitude of the stress at the crack tip increases with increasing θ , and the amount of deformation required for the crack to propagate through the entire sample decreases with increasing θ . However, the crack path does not change. This is because the direction of the maximum principal stress at the crack tip is not changing, just the magnitude, such that the crack path that minimizes the energy is
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
9
Fig. 7. Impact of crystal orientation on fracture behavior in hcp Cd single crystals, where (a) shows the σ yy stress distribution in GPa for θ = 0◦ and v = 8.79 µm with regions where d = 1 colored white, (b) shows the same stress distribution for θ = 90◦ and v = 6.31 µm, and (c) shows the average stress–strain curves for the various crystal orientations. Table 3 Material parameters of fcc copper used in simulations. c11 (GPa)
c44 (GPa)
c12 (GPa)
gc (GPa mm)
η (s/mm)
127.0
73.55
70.8
1.0 × 10−3
1.0 × 10−6
unchanged. In order for the crack path to change with orientation, gc would need to be anisotropic, such that certain crack angles would create larger surface energies than others. The stress–strain curves for different orientations are shown in Fig. 7c, with θ = 0◦ resulting in the lowest stress and the largest strain at fracture and θ = 90◦ resulting in the highest stress and lowest strain. The area under the stress–strain curves represents the total external work, including the impact of the strain energy and the fracture energy. The fracture energy should be independent of orientation, since gc is isotropic in this model, but the strain energy should increase with orientation. This is consistent with our results, with the area under the stress–strain curve increasing with orientation such that θ = 90◦ has an area that is 14% larger than θ = 0◦ . To verify that our approach properly accounts for cubic symmetry, fracture in face-centered cubic (fcc) copper was also modeled using the parameters from Table 3. The σ yy stress field and crack path for various values of the orientation angle θ are shown in Figs. 8a–8e. The stress field is identical for θ = 0◦ and 90◦ . The stress field for θ = 30◦ and 60◦ are similar, but antisymmetric across the crack plane. The crack paths are again independent of crystal orientation, for the same reason as for the hcp example. The average stress–strain curves are shown in Fig. 8f, and they are equivalent for θ = 0◦ and 90◦ and also for θ = 30◦ and 60◦ . θ = 45◦ exhibits the highest stress and lowest strain at fracture while θ = 0◦ and 90◦ exhibit the lowest stress and largest strain. The area under the stress–strain curves for 45◦ is 10% larger than for 0◦ and 90◦ . These results indicate that our model properly accounts for the impact of crystal orientation in a cubic material on mode I fracture. Mode II fracture simulations of fcc copper were also conducted for comparison, where the boundary conditions and geometry sizes are shown in Fig. 5. Crystal orientations ranging from 0◦ to 90◦ were performed for comparison and the resultant crack paths and stress distributions are shown in Figs. 9a–9e. For mode II fracture, the stress field for θ = 0◦ and 90◦ is identical and the crack propagates at the largest deformation. The crack propagates at the lowest deformation for θ = 60◦ . Unlike mode I, in mode II the orientation impacts the crack path, as shown in Fig. 9f. Again, cubic symmetry is evident in the response of the material. 4.3.2. Fracture simulation in bicrystals Having established that our model accurately accounts for elastic anisotropy in single crystals, we now investigate the behavior of hcp Cd bicrystals with an initial crack. The properties from Table 2 were used in these simulations. The geometry and boundary conditions for these simulations are shown in Fig. 10. In our first bicrystal simulations, we considered two crystal orientations, θ = 0◦ and 90◦ . Two bicrystals were modeled, one with the θ = 90◦ crystal on the top (case I) and one with it on the bottom (case II). Note that in these simulations, the grain boundary
10
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
Fig. 8. Impact of crystal orientation on mode I fracture behavior in fcc single crystals, where (a)–(e) show the σ yy stress distribution in GPa for various orientations and applied displacements, with regions where d = 1 colored white, and (f) shows the stress–strain curves for the various crystal orientations. Note that the θ = 0◦ line is obscured by the θ = 90◦ because their behaviors are identical.
Fig. 9. Impact of crystal orientation on mode II fracture behavior in fcc copper single crystals, where (a)–(e) show the σx x stress distribution in GPa for various orientations and applied displacements, with regions where d = 1 colored white, and (f) overlays the different crack paths, showing that they change with orientation.
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
11
Fig. 10. Boundary conditions and geometry parameters for the bicrystal fracture simulations.
Fig. 11. Crack propagation in Cd bicrystals. (a) and (b) show the normalized elastic energy (|ψe | = ψe / max(ψe )) prior to crack propagation (left) and the σ yy stress field in GPa and post-processed crack path at a displacement of v = 7.9 µm (right) for case I and case II, respectively, with regions where d = 1 colored white. Note that the elastic energy field is not symmetric across the grain boundary, causing the crack to propagate through the θ = 0◦ grain. The stress–strain curves for these two cases are shown in (c). The single crystal stress–strain curves for θ = 0◦ and θ = 90◦ are also shown for comparison. They are plotted using a dashed line and labeled Single in the legend.
behaves in the same manner as the rest of the material. It is only distinguished by a transition in the elasticity tensor, due to the different crystal orientations. The elastic energy density field and the σ yy stress field and crack path from both cases are shown in Figs. 11a–11b. The results show different crack paths occurring in the two bicrystals, with the crack always propagating through the θ = 0◦ grain. The crack path changes because the elastic energy density at the crack tip is not symmetric around the grain boundary. At face value it is surprising that the crack propagates through the θ = 0◦ grain, as the single crystal simulations showed that the stress is higher and fracture occurred at lower strains when θ = 90◦ , as
12
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
Fig. 12. Results from Cd bicrystal fracture simulations with different crystal orientations of the top grain. (a)–(d) show the σ yy stress field in GPa and the post-processed crack path for various orientations of the top grain and approximately the same crack length (showing the corresponding displacement), with regions where d = 1 colored white. The stress–strain curves for the bicrystals with different angles are shown in (e).
shown in Fig. 7c. However, when the tensile load is applied to the material, the softer grain (θ = 0◦ ) will deform more than the θ = 90◦ grain. For this reason, it has a higher deformation energy and so the energy decreases faster through fracture in the θ = 0◦ grain. The stress–strain curves for both bicrystals are shown in Fig. 11c, and the two curves are nearly identical. This is because in both bicrystals, the crack is propagating through the same crystal orientation. The curves do not match those for the single crystal fracture with θ = 90◦ , nor do they match the single crystal curve with θ = 0◦ , even though the crack propagates through the θ = 0◦ grain. This is because the deformation energy throughout the bicrystal is impacted by both grains. To study the influence of crystal orientation on the crack path in bicrystals, we study four polycrystals in which the orientation of the top grain varies from θ = 30◦ to 90◦ , and the bottom grain has no rotation. The σ yy stress fields and crack paths are shown in Figs. 12a–12d. In each case, the crack propagates through the bottom, unrotated grain. In addition, as the orientation angle increases, the crack curves further from the grain boundary. Thus, as the difference in the elasticity tensor between the two grains increases, the energy drives the crack further away from the grain boundary. The stress–strain curves, shown in Fig. 12e change with crystal orientation, but not nearly to the magnitude observed in the single crystal simulations (Fig. 7c). 4.3.3. Fracture simulation in polycrystal structures As our final example, we investigate the impact of elastic anisotropy on the fracture behavior in 2D polycrystals. The polycrystal is composed of Cd, using the properties from Table 2, and has 36 grains generated using a Voronoi tessellation. The orientations of the grains are randomly generated, ranging from 0◦ to 90◦ , and three different sets of random orientations are generated to compare the impact on the crack behavior. The three polycrystals have an initial crack of length 0.3 mm and each have the same grain structure but different grain orientations as shown in Figs. 13a–13c. As with the bicrystal simulations, the grain boundaries behave in the same manner as the rest of the material. The polycrystals have zero x-displacements on the left side of the structure, zero y-displacements on the bottom surface, and they are loaded on the top surface with a positive y-displacement load with a rate of
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
13
Fig. 13. Polycrystal simulation initial conditions, where (a)–(c) show the crystal orientations and initial crack for the three different cases, and (d)–(f) shows the corresponding C1111 elasticity tensor component in GPa.
1 × 10−3 mm/s. To show the influence of crystal orientation on the elastic constants, the C1111 component of the heterogeneous elasticity tensor is shown in Figs. 13d–13f. The σ yy stress field and crack path for the three polycrystals are shown in Figs. 14a–14c. The impact of the grain orientations is clearly evident on the stress distribution, resulting in some deviation in the crack path near grain boundaries. In addition, the crack propagation rate is directly impacted by the crystal orientations, as seen in the release portion of the stress–strain curves (shown in Fig. 14d). Large changes in the decreasing slope result from the crack moving between grains with very different orientations. For example, this occurs very early in the crack propagation in Case I, it occurs late in the propagation in Case II, and it occurs twice in Case III. 5. Conclusions In this paper, we present a new phase field model of brittle fracture that includes the impact of anisotropic elastic behavior. The stress is decomposed into its tensile and compressive parts using a spectral decomposition, such that only the tensile energy contributes to fracture. To verify that it properly captures the behavior of isotropic materials, this new model was applied to mode I and mode II fracture in isotropic materials, and was shown to predict the expected crack paths. To investigate its ability to model fracture in anisotropic materials, fracture in single crystal, bicrystal, and polycrystal domains was simulated. Single crystal fracture with different crystal orientations was investigated in both hcp and fcc materials. In each case, the crack path was unaffected by the crystal orientation but the stress–strain curves were significantly impacted. The fcc material demonstrated cubic symmetry in its behavior. The bicrystal simulations showed that the fracture behavior was impacted by the orientations of both grains and that the bicrystal stress field diverted the crack path. The polycrystal simulations were conducted with three different
14
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
Fig. 14. Results from polycrystal fracture simulations with different configurations of crystal orientations , where (a)–(c) show the postprocessed crack path and σ yy stress field in GPa for the different orientation distributions, with regions where d = 1 colored white. The resultant stress–strain curves are shown in (d).
sets of crystal orientations, and the final crack paths and stress–strain curves showed the impact of the different orientations. Acknowledgments Zhang and Tonks gratefully acknowledge financial support from the United States Department of Energy, Office of Nuclear Energy under Nuclear Energy University Program (NEUP) Grant 15-8243. Zhang and Jiang gratefully acknowledge financial support from the United States Department of Energy, Office of Nuclear Energy under Nuclear Energy University Program (NEUP) Grant 15-8229. References [1] M. Elices, G. Guinea, J. Gomez, J. Planas, The cohesive zone model: advantages, limitations and challenges, Eng. Fract. Mech. 69 (2) (2002) 137–163. [2] N. Moës, T. Belytschko, Extended finite element method for cohesive crack growth, Eng. Fract. Mech. 69 (7) (2002) 813–833. [3] A. Menk, S.P. Bordas, Numerically determined enrichment functions for the extended finite element method and applications to bi-material anisotropic fracture and polycrystals, Internat. J. Numer. Methods Engrg. 83 (7) (2010) 805–828. [4] G.A. Francfort, J.-J. Marigo, Revisiting brittle fracture as an energy minimization problem, J. Mech. Phys. Solids 46 (8) (1998) 1319–1342. [5] B. Bourdin, G.A. Francfort, J.-J. Marigo, Numerical experiments in revisited brittle fracture, J. Mech. Phys. Solids 48 (4) (2000) 797–826.
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643 [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32]
15
I. Aranson, V. Kalatsky, V. Vinokur, Continuum field description of crack propagation, Phys. Rev. Lett. 85 (1) (2000) 118. A. Karma, D.A. Kessler, H. Levine, Phase-field model of mode III dynamic fracture, Phys. Rev. Lett. 87 (4) (2001) 045501. V. Hakim, A. Karma, Laws of crack motion and phase-field models of fracture, J. Mech. Phys. Solids 57 (2) (2009) 342–368. H. Henry, H. Levine, Dynamic instabilities of fracture under biaxial strain using a phase field model, Phys. Rev. Lett. 93 (10) (2004) 105504. C. Kuhn, R. Müller, A phase field model for fracture, PAMM 8 (1) (2008) 10223–10224. H. Amor, J.-J. Marigo, C. Maurini, Regularized formulation of the variational brittle fracture with unilateral contact: numerical experiments, J. Mech. Phys. Solids 57 (8) (2009) 1209–1229. C. Miehe, M. Hofacker, F. Welschinger, A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits, Comput. Methods Appl. Mech. Engrg. 199 (45) (2010) 2765–2778. C. Miehe, F. Welschinger, M. Hofacker, Thermodynamically consistent phase-field models of fracture: variational principles and multi-field FE implementations, Internat. J. Numer. Methods Engrg. 83 (10) (2010) 1273–1311. M.J. Borden, C.V. Verhoosel, M.A. Scott, T.J. Hughes, C.M. Landis, A phase-field description of dynamic brittle fracture, Comput. Methods Appl. Mech. Engrg. 217 (2012) 77–95. M.J. Borden, T.J. Hughes, C.M. Landis, C.V. Verhoosel, A higher-order phase-field model for brittle fracture: Formulation and analysis within the isogeometric analysis framework, Comput. Methods Appl. Mech. Engrg. 273 (2014) 100–118. P. Chakraborty, Y. Zhang, M.R. Tonks, Multi-scale modeling of microstructure dependent intergranular brittle fracture using a quantitative phase-field based method, Comput. Mater. Sci. 113 (2016) 38–52. P. Chakraborty, P. Sabharwall, M.C. Carroll, A phase-field approach to model multi-axial and microstructure dependent fracture in nuclear grade graphite, J. Nucl. Mater. 475 (2016) 200–208. V. Hakim, A. Karma, Crack path prediction in anisotropic brittle materials, Phys. Rev. Lett. 95 (23) (2005) 235501. A. Abdollahi, I. Arias, Phase-field modeling of fracture in ferroelectric materials, Arch. Comput. Methods Eng. 22 (2) (2015) 153–181. B. Li, C. Peco, D. Millán, I. Arias, M. Arroyo, Phase-field modeling and simulation of fracture in brittle materials with strongly anisotropic surface energy, Internat. J. Numer. Methods Engrg. 102 (3–4) (2015) 711–727. J. Clayton, J. Knap, A geometrically nonlinear phase field theory of brittle fracture, Int. J. Fract. 189 (2) (2014) 139–148. J. Clayton, J. Knap, Phase field modeling of directional fracture in anisotropic polycrystals, Comput. Mater. Sci. 98 (2015) 158–169. S. Teichtmeister, C. Miehe, Phase-field modeling of fracture in anisotropic media, PAMM 15 (1) (2015) 159–160. P. Shanthraj, B. Svendsen, L. Sharma, F. Roters, D. Raabe, Elasto-viscoplastic phase field modelling of anisotropic cleavage fracture, J. Mech. Phys. Solids 99 (2017) 19–34. T.T. Nguyen, J. Réthoré, M.-C. Baietto, Phase field modelling of anisotropic crack propagation, Eur. J. Mech. A Solids 65 (2017) 279–288. D. Gaston, C. Newman, G. Hansen, D. Lebrun-Grandie, Moose: A parallel computational framework for coupled systems of nonlinear equations, Nucl. Eng. Des. 239 (10) (2009) 1768–1778. B. Bourdin, G.A. Francfort, J.-J. Marigo, The variational approach to fracture, J. Elasticity 91 (1–3) (2008) 5–148. C. Miehe, Comparison of two algorithms for the computation of fourth-order isotropic tensor functions, Comput. Struct. 66 (1) (1998) 37–43. C. Miehe, M. Lambrecht, Algorithms for computation of stresses and elasticity moduli in terms of seth–hill’s family of generalized strain tensors, Int. J. Numer. Methods Biomed. Eng. 17 (5) (2001) 337–353. V. Lubarda, D. Krajcinovic, S. Mastilovic, Damage model for brittle elastic solids with unequal tensile and compressive strengths, Eng. Fract. Mech. 49 (5) (1994) 681–697. S. Murakami, Continuum Damage Mechanics: a Continuum Mechanics Approach to the Analysis of Damage and Fracture, Vol. 185, Springer Science & Business Media, 2012. M.R. Tonks, D. Gaston, P.C. Millett, D. Andrs, P. Talbot, An object-oriented finite element framework for multiphysics phase field simulations, Comput. Mater. Sci. 51 (1) (2012) 20–29.
ScienceDirect Comput. Methods Appl. Mech. Engrg. 358 (2020) 112643 www.elsevier.com/locate/cma
A new phase field fracture model for brittle materials that accounts for elastic anisotropy Shuaifang Zhanga,b , Wen Jiangc , Michael R. Tonksa,b ,∗ b
a University of Florida, Gainesville, FL, United States of America The Pennsylvania State University, State College, PA, United States of America c Idaho National Laboratory, Idaho Falls, ID, United States of America
Received 22 June 2018; received in revised form 10 September 2019; accepted 11 September 2019 Available online xxxx
Highlights • • • •
The new phase field fracture model includes anisotropic elasticity. Compressive stress does not contribute to crack propagation. Single crystal and bicrystal simulations show the proper relationship with crystal orientation. Polycrystal simulations show the impact of texture on the stress–strain curves.
Abstract In this manuscript, a new phase field model of brittle fracture is proposed that accounts for any elastic anisotropy using spectral decomposition of the stress. To verify the model, both Mode I and Mode II fracture simulations were performed to ensure that the correct crack paths were predicted. Next, the fracture of face-centered cubic (FCC) and hexagonal close packed (HCP) single crystals was modeled, comparing the impact of crystal orientation on the stress–strain curve. Third, fracture in an HCP bicrystal was simulated and the impact of crystal orientation on the crack paths and stress–strain curves was studied. Finally, a polycrystal structure was simulated to study the impact of crystallographic texture on transgranular fracture. c 2019 Elsevier B.V. All rights reserved. ⃝ Keywords: Phase field; Fracture mechanics; Anisotropic materials; Brittle fracture; MOOSE
1. Introduction Various methods have been developed to model discrete crack initiation and growth, including cohesive zone models [1] and the extended finite element method [2], and these methods have been developed for isotropic and anisotropic materials [3]. Recently, the phase field method has emerged as a powerful means of modeling discrete crack nucleation and growth. A large part of its value comes from its simplicity; it can predict crack propagation solely by solving a system of partial differential equations (PDEs). This is because the system of PDEs simultaneously predicts crack formation and creates new cracked free surface. ∗ Corresponding author at: University of Florida, Gainesville, FL, United States of America.
E-mail addresses: [email protected] (S. Zhang), [email protected] (W. Jiang), [email protected] (M.R. Tonks). https://doi.org/10.1016/j.cma.2019.112643 c 2019 Elsevier B.V. All rights reserved. 0045-7825/⃝
2
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
The foundations for phase field methods to model crack propagation in brittle materials was established in the late 1990s, with the variational formulation of brittle fracture by Francfort [4] and Bourdin [5]. In these original models, the quasi-static brittle fracture process is governed by the minimization of the total energy, which is a function of strain energy and crack dissipation. Aranson et al. [6] presented a phase field fracture model to solve mode I fracture based on Ginzburg–Landau phase transition. Karma and collaborators [7,8] modified the strain energy density function from Aranson et al. [6] to simulate mode III fracture. Henry et al. [9] modified the energy function of Aranson et al. [6] to work for mode I and II fracture. Kuhn and Muller [10] were the first to propose the thermodynamics-based phase field model. In this model, a continuous variable field was used to represent cracked material. The system strain energy density was determined by solving the stress divergence equation and the evolution of the crack variable field was defined by a partial differential equation that minimized the energy. Amor [11] separated the strain energy density into two parts, one representing the tensile part of the strain energy, and one representing the compressive part, which did not contribute to the crack propagation. Miehe and collaborators [12,13] proposed a phase field fracture model based on thermodynamic consistency and continuum mechanics in which the strain energy was also decomposed into the tensile and compressive parts. Borden [14,15] developed a higher order equation for the evolution of the phase field parameter that required advanced finite element formulations to solve. Such phase field fracture models have been applied to model the impact of microstructural features on fracture behavior [16,17]. Initial phase field fracture models were developed assuming isotropic material behavior, however newer models have focused on anisotropic materials. With respect to brittle fracture, materials are anisotropic in two ways: anisotropic elastic response and anisotropic surface energy. Hakim and Karma (2005) [18] derived a force balance condition that was used to account for anisotropic surface energy, in which the path is determined by the directional fracture energy, but not elastic anisotropy. Abdollahi and Arias [19] also developed a model for anisotropic surface energy but not anisotropic elasticity, but in ferroelectrics under indentation with electro-mechanical coupling. Bin Li et al. [20] proposed a phase field fracture model that introduced a strongly anisotropic surface energy and anisotropic elasticity, though the contribution of compressive stress to crack propagation is not avoided. Clayton and Knap developed models with anisotropic surface energy using an isotropic elasticity tensor [21] and then an anisotropic elasticity tensor [22]; however, they also did not avoid the contribution of compressive stress to crack propagation. Teichtmeister and Miehe [23] employed a fourth-order structural tensor that was specific for cubic anisotropy. Shanthraj [24] developed a finite-strain viscoplastic phase field method that considered anisotropic surface energy, but not elastic anisotropy. Nguyen et al. [25] also considered surface energy anisotropy but not elastic anisotropy. None of these anisotropic phase field fracture models provide both a means of avoiding the impact of compressive stress on crack propagation and a means of considering arbitrary elastic anisotropy. In this paper, a new phase field fracture model is proposed that employs the full fourth-order elastic tensor, such that any form of anisotropy can be included. The model decomposes the strain energy into compressive and tensile parts and only the tensile stress contributes to crack propagation. The necessary equations are derived and then verified by modeling mode I and II fracture in an isotropic material. Then, the model is demonstrated on anisotropic materials. The model was developed using the finite element method with the open source Multiphysics Object-Oriented Simulation Environment (MOOSE) [26]. 2. Phase field fracture model in anisotropic brittle materials 2.1. Phase field approximation of crack In contrast to traditional finite element methods that are used to simulate discrete cracks in which the cracks are sharp, phase field fracture models employ a diffuse-interface description of cracks. The continuous variable d is used to characterize unbroken material (d = 0) and fully cracked material (d = 1), and the material is considered as damaged when 0 < d < 1. The shape of d across the crack width is described by the differential equation d(x) − l 2
∂ 2d =0 ∂x2
(1)
with solution d(x) = e−|x|/l ,
(2)
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
3
Fig. 1. Phase field parameter as a function of spatial coordinate: d(x) = e−|x|/l .
where l is the length parameter used to represent the crack width, as shown in Fig. 1. The crack approximation converges to a sharp crack as l approaches zero. In multiple dimensions, the differential equation that governs the shape of d is d − l 2 ∆d = 0
(3)
2.2. Crack propagation with a phase field model in anisotropic brittle materials The total potential energy of a cracked brittle material is defined as a sum of elastic energy and crack surface energy [14]: ∫ ∫ gc dΓ , (4) ψe dΩ + Ψtot = Ω
Γ
where ψe is the elastic energy density, gc is the crack surface energy, Ω represents the entire material, and Γ is the crack surface area. As the crack in the phase field model now has a finite thickness, the crack energy can be expressed as a volume integral using the crack surface density function per unit volume from Miehe [13], which is an extension of Bourdin et al. [27]: 1 l γ = d 2 + |∇d|2 , (5) 2l 2 ∫ ∫ such that Γ gc dΓ = Ω gc ( 2l1 d 2 + 2l |∇d|2 )dΩ . Thus, Eq. (4) becomes ) ∫ ( gc 2 gc l 2 ψe + d + Ψtot = |∇d| dΩ . (6) 2l 2 Ω As we are considering only brittle materials, it is reasonable to assume small deformation, so that the strain tensor is defined as ) 1( ∇u + ∇uT . (7) ϵ= 2 We also assume a linear elastic material response, such that the stress tensor of the unbroken material σ 0 = Cϵ,
(8)
where C is the fourth-order elasticity tensor, with as much symmetry as desired. In our model, we divide the elastic energy density into its tensile ψ + and compressive ψ − parts, such that ψe = ((1 − d)2 (1 − k) + k)ψ + + ψ − ,
(9)
where k is a small value that is used to ensure positive definiteness after the material is fully broken (d = 1). The energy is split for two reasons: first, to ensure that only tensile stress contributes to crack propagation; second, to ensure cracks only support compressive stress. We apply spectral decomposition of the unbroken stress tensor to determine values for ψ + and ψ − σ 0 = QΛ Q T ,
(10)
4
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
where Λ = diag(λ1 , λ2 , λ3 , ) is a diagonal matrix containing the three eigenvalues of the stress tensor (or the principal stress tensor) and Q is the eigenvector matrix of the stress. We then define the tensile Λ+ and compressive Λ− principal stress tensors as Λ+ = I+ ΛI+ Λ− = I− ΛI−
(11) (12)
with I+ and I− being similar to the identify matrix, except that Iii+ = 0 when λi < 0 and Iii− = 0 when λi > 0. From the principal stress tensors, we can determine the compressive and tensile stress tensors σ + = QΛ+ QT σ − = QΛ− QT ,
(13) (14)
and thus 1 + σ :ϵ (15) 2 1 ψ − = σ − : ϵ. (16) 2 The primary differences between this model and that from Miehe [13] are that this model includes the full rank four elasticity tensor, while Miehe’s model only includes the two isotropic elastic constants, and that this model directly employs a spectral decomposition of the elastic stress tensor, while Miehe’s model decomposes the elastic stress tensor using a function of the eigenvalues and eigenvectors of the strain tensor and the two isotropic elastic constants. We also define projection tensors that separate the stress into its compressive and tensile parts ψ+ =
σ + = P+ σ 0 σ − = P− σ 0 ,
(17) (18)
where ∂σ + ∂σ 0 − P = I − P+ . P+ =
(19) (20)
These projection tensors are calculated from the spectral decomposition of the unbroken stress tensor using the algorithms defined by Miehe [28] and Miehe and Lambrecht [29]. This relationship involving the projection tensor has been widely used in damage mechanics [30,31] and was also used in phase field fracture by Miehe [12]. 3. Governing equations The fracture model is defined by two equations, the stress equilibrium equation and the evolution equation. In this section, we will derive these two equations. 3.1. Equilibrium equation Eq. (9) defines the strain energy density as a function of the phase field parameter d, which represents the damaged region of the material. The stress tensor in the damaged anisotropic material is defined as the derivative of the strain energy density with respect to the strain tensor, ∂ψe ∂ψ + ∂ψ − = ((1 − d)2 (1 − k) + k) + . ∂ϵ ∂ϵ ∂ϵ The partial derivative of Eq. (15) with respect to the strain is σ =
∂ψ + 1 dσ + = ϵ+ ∂ϵ 2 dϵ + 1 dσ = ϵ+ 2 dϵ
1 + dϵ σ 2 dϵ 1 + σ 2
(21)
(22) (23)
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
5
We take the partial derivative of σ + with respect to the strain and substitute in Eq. (8) to obtain dσ + dσ + dσ 0 dσ + = = C dϵ dσ 0 dϵ dσ 0 and the partial derivative of Eq. (17) with respect to the unbroken stress to get
(24)
dσ + = P+ . (25) dσ 0 The expressions involving the compressive stress are obtained in a similar manner. If we combine these expressions with Eq. (19) we obtain ∂ψ + = σ+ ∂ϵ ∂ψ − = σ −. ∂ϵ Then, we substitute these into Eq. (21) to obtain σ = ((1 − d)2 (1 − k) + k)σ + + σ − .
(26) (27)
(28)
This expression for the stress is then used in the stress equilibrium equation to solve the mechanics problem. 3.2. Phase field evolution equation The evolution of the damage parameter d minimizes the total free energy of the system using an Allen–Cahn equation: ) ( ∂ψe ∂d δΨtot gc = −L = −L + d − ∇ · gc l∇d . (29) ∂t δd ∂d l Combining this expression with Eq. (9) gives ( ) ∂d gc = −L (−2(1 − d)(1 − k) + k)ψ + + d − ∇ · gc l∇d . (30) ∂t l This equation does not guarantee that cracks do not heal. Therefore, we alter it by introducing the history variable H , which is the maximum positive tensile energy experienced during the simulation, H = max ψ + for t = [0, tc ],
(31)
where tc is the current time. This is the same approach taken in Miehe et al. [12]. Thus, the evolution equation becomes ) 1( gc ∂d = − (−2(1 − d)(1 − k) + k)H + d − ∇ · gc l∇d . (32) ∂t η l The expression L = 1/η, with the rate parameter η (called viscosity by Miehe et al. [12]), has been used to obtain the form from Miehe et al. [13]. Note that this equation is identical to the Allen–Cahn equation used in traditional phase field simulations, ( ) ∂d ∂ floc = −L − ∇ · κ∇d , (33) ∂t ∂d with floc = ((1 − d)2 (1 − k) + k)H + ψ − + 1 η κ = gc l.
L =
gc 2 d 2l
(34) (35) (36)
This model has been implemented with the finite element method using MOOSE [32], as was done for isotropic materials by Chakroborty et al. [17]. It is solved using implicit time integration.
6
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
Fig. 2. Geometry and boundary conditions used in the mode I fracture simulations shown in Figs. 3, 7 and 9. Table 1 Material properties used in the Mode I and II fracture simulations of an isotropic material. E (GPa) 120
ν
gc (GPa mm)
η (s/mm)
k
0.3
1.0 × 10−3
1.0 × 10−7
1.0 × 10−6
4. Numerical results To validate that our fracture model adequately represents the fracture of brittle materials, we perform a series of 2D plane strain simulations. First we model mode I and II fracture of an isotropic material (with two unique elastic constants, the elastic modulus E and Poisson’s ratio ν). Then, we model fracture in hexagonal and cubic single crystals, as well as hexagonal bicrystals. Finally, we model fracture in hexagonal polycrystals. 4.1. Mode I fracture in isotropic materials To model mode I fracture in an isotropic material, we used a single edge notch tension specimen shown in Fig. 2. The dimensions and boundary conditions of the specimen are also shown in the figure. The material properties used for the phase field fracture simulations do not represent a specific material and are shown in Table 1. The crack length parameter is l = 0.005 mm. Chakraborty et al. [16] studied the mesh size requirement for phase field fracture simulations and found that the interface was adequately resolved when the smallest element size was equal to half the length parameter (l/2). In this work, we employed mesh adaptivity to reduce computational expense but ensure that the elements across the crack are always l/2 or smaller. Different length parameters were used in the different simulations in order to obtain a good balance between computational cost (which increases as l decreases) and similarity to a sharp crack model (which also generally increases as l decreases). In all the simulations, the time step size dt = 1 × 10−4 s. The final crack path and stress–strain curve for mode I fracture are shown in Figs. 3a and 3b, respectively. The stress field at different displacements of the top boundary are shown in Figs. 3b and 3c. From the stress evolution and crack propagation plots we can see that the crack propagates parallel to the initial crack, consistent with the mode I fracture behavior shown in Miehe et al. [12]. The material behavior under a compressive load was also investigated, where v = −t. The boundary conditions were slightly modified, such that the x-displacements are set to zero on the right side and the y-displacements are set to zero on the bottom. A larger compressive displacement load was applied than was applied in the tensile simulations and yet no fracture occurred, as shown in Fig. 4. This verifies that cracks do not propagate due to compressive stress. 4.2. Mode II fracture test in isotropic materials Shear force induced cracks are also very common in engineering structures, so we conducted simulations with an applied shear load (mode II fracture). The geometry and boundary conditions of the mode II crack simulation
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
7
Fig. 3. Crack propagation in the Mode I fracture simulation of an isotropic material with l = 0.005 mm. The final crack path is shown in (a); the σ yy stress distribution in GPa at displacements of v = 3.65 µm and v = 3.68 µm are shown in (b) and (c), respectively, with regions where d = 1 colored white. The resultant stress–strain curve is shown in (d).
Fig. 4. Material response under a compressive load of the isotropic material with l = 0.005 mm at v = −10 µm, where the σ yy stress distribution in GPa is shown in (a) and the phase field variable field is shown in (b). Note that no crack propagation occurs due to the uniform compressive stress distribution.
are shown in Fig. 5, in which the single notched specimen is used and the bottom boundary is fixed while the top side has an applied x-displacement. The material properties from Table 1 are used; the length parameter l = 0.01 mm. Figs. 6b and 6c show the crack evolution and stress evolution at different x-displacements of the top surface. As shown in Fig. 6a, the crack path of mode II fracture is consistent with classical mode II fracture results and with the results from other phase field fracture models [5,12], proving that our model works well with mode II fracture.
8
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
Fig. 5. Boundary conditions and geometry parameters for the Mode II fracture simulation shown in Fig. 6.
Fig. 6. Crack propagation in the Mode II fracture simulation of the isotropic material, where the final crack path is shown in (a) and the σx x stress distribution in GPa at displacements u = 6.43 µm and u = 8.92 µm are shown in (b) and (c), respectively, with regions where d = 1 colored white. Table 2 Picked material properties of HCP structure. c11 (GPa) 115.8
c12 (GPa) 39.8
c13 (GPa) 40.6
c33 (GPa) 51.4
c44 (GPa)
gc (GPa mm)
η (s/mm)
20.4
1.0 × 10−3
1.0 × 10−6
4.3. Brittle fracture in anisotropic materials 4.3.1. Fracture simulation in single crystals with anisotropic elasticity In this section, we demonstrate the influence of crystal orientation on the fracture behaviors in anisotropic single crystals. For all of these simulations, we employed the geometry and boundary conditions shown in Fig. 2. Firstly, a transversely isotropic material was studied, where hexagonal close-packed (hcp) crystals are the most common example. Second, we modeled a cubic material. There are five unique elastic constants in hcp crystals, and here we used cadmium (Cd) as the example material, with the material properties shown in Table 2. Note that this material is symmetric along the x- and y-directions, so the 2D simulations are done in the x-z plane. To verify that the crystal orientation influences the fracture behavior in hcp Cd, simulations were run for single crystals with different orientations, defined by a rotation angle θ around the out-of-plane axis. Examples of the σ yy stress field and the crack field after a given amount of deformation for two representative crystal orientations are shown Figs. 7a and 7b. The magnitude of the stress at the crack tip increases with increasing θ , and the amount of deformation required for the crack to propagate through the entire sample decreases with increasing θ . However, the crack path does not change. This is because the direction of the maximum principal stress at the crack tip is not changing, just the magnitude, such that the crack path that minimizes the energy is
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
9
Fig. 7. Impact of crystal orientation on fracture behavior in hcp Cd single crystals, where (a) shows the σ yy stress distribution in GPa for θ = 0◦ and v = 8.79 µm with regions where d = 1 colored white, (b) shows the same stress distribution for θ = 90◦ and v = 6.31 µm, and (c) shows the average stress–strain curves for the various crystal orientations. Table 3 Material parameters of fcc copper used in simulations. c11 (GPa)
c44 (GPa)
c12 (GPa)
gc (GPa mm)
η (s/mm)
127.0
73.55
70.8
1.0 × 10−3
1.0 × 10−6
unchanged. In order for the crack path to change with orientation, gc would need to be anisotropic, such that certain crack angles would create larger surface energies than others. The stress–strain curves for different orientations are shown in Fig. 7c, with θ = 0◦ resulting in the lowest stress and the largest strain at fracture and θ = 90◦ resulting in the highest stress and lowest strain. The area under the stress–strain curves represents the total external work, including the impact of the strain energy and the fracture energy. The fracture energy should be independent of orientation, since gc is isotropic in this model, but the strain energy should increase with orientation. This is consistent with our results, with the area under the stress–strain curve increasing with orientation such that θ = 90◦ has an area that is 14% larger than θ = 0◦ . To verify that our approach properly accounts for cubic symmetry, fracture in face-centered cubic (fcc) copper was also modeled using the parameters from Table 3. The σ yy stress field and crack path for various values of the orientation angle θ are shown in Figs. 8a–8e. The stress field is identical for θ = 0◦ and 90◦ . The stress field for θ = 30◦ and 60◦ are similar, but antisymmetric across the crack plane. The crack paths are again independent of crystal orientation, for the same reason as for the hcp example. The average stress–strain curves are shown in Fig. 8f, and they are equivalent for θ = 0◦ and 90◦ and also for θ = 30◦ and 60◦ . θ = 45◦ exhibits the highest stress and lowest strain at fracture while θ = 0◦ and 90◦ exhibit the lowest stress and largest strain. The area under the stress–strain curves for 45◦ is 10% larger than for 0◦ and 90◦ . These results indicate that our model properly accounts for the impact of crystal orientation in a cubic material on mode I fracture. Mode II fracture simulations of fcc copper were also conducted for comparison, where the boundary conditions and geometry sizes are shown in Fig. 5. Crystal orientations ranging from 0◦ to 90◦ were performed for comparison and the resultant crack paths and stress distributions are shown in Figs. 9a–9e. For mode II fracture, the stress field for θ = 0◦ and 90◦ is identical and the crack propagates at the largest deformation. The crack propagates at the lowest deformation for θ = 60◦ . Unlike mode I, in mode II the orientation impacts the crack path, as shown in Fig. 9f. Again, cubic symmetry is evident in the response of the material. 4.3.2. Fracture simulation in bicrystals Having established that our model accurately accounts for elastic anisotropy in single crystals, we now investigate the behavior of hcp Cd bicrystals with an initial crack. The properties from Table 2 were used in these simulations. The geometry and boundary conditions for these simulations are shown in Fig. 10. In our first bicrystal simulations, we considered two crystal orientations, θ = 0◦ and 90◦ . Two bicrystals were modeled, one with the θ = 90◦ crystal on the top (case I) and one with it on the bottom (case II). Note that in these simulations, the grain boundary
10
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
Fig. 8. Impact of crystal orientation on mode I fracture behavior in fcc single crystals, where (a)–(e) show the σ yy stress distribution in GPa for various orientations and applied displacements, with regions where d = 1 colored white, and (f) shows the stress–strain curves for the various crystal orientations. Note that the θ = 0◦ line is obscured by the θ = 90◦ because their behaviors are identical.
Fig. 9. Impact of crystal orientation on mode II fracture behavior in fcc copper single crystals, where (a)–(e) show the σx x stress distribution in GPa for various orientations and applied displacements, with regions where d = 1 colored white, and (f) overlays the different crack paths, showing that they change with orientation.
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
11
Fig. 10. Boundary conditions and geometry parameters for the bicrystal fracture simulations.
Fig. 11. Crack propagation in Cd bicrystals. (a) and (b) show the normalized elastic energy (|ψe | = ψe / max(ψe )) prior to crack propagation (left) and the σ yy stress field in GPa and post-processed crack path at a displacement of v = 7.9 µm (right) for case I and case II, respectively, with regions where d = 1 colored white. Note that the elastic energy field is not symmetric across the grain boundary, causing the crack to propagate through the θ = 0◦ grain. The stress–strain curves for these two cases are shown in (c). The single crystal stress–strain curves for θ = 0◦ and θ = 90◦ are also shown for comparison. They are plotted using a dashed line and labeled Single in the legend.
behaves in the same manner as the rest of the material. It is only distinguished by a transition in the elasticity tensor, due to the different crystal orientations. The elastic energy density field and the σ yy stress field and crack path from both cases are shown in Figs. 11a–11b. The results show different crack paths occurring in the two bicrystals, with the crack always propagating through the θ = 0◦ grain. The crack path changes because the elastic energy density at the crack tip is not symmetric around the grain boundary. At face value it is surprising that the crack propagates through the θ = 0◦ grain, as the single crystal simulations showed that the stress is higher and fracture occurred at lower strains when θ = 90◦ , as
12
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
Fig. 12. Results from Cd bicrystal fracture simulations with different crystal orientations of the top grain. (a)–(d) show the σ yy stress field in GPa and the post-processed crack path for various orientations of the top grain and approximately the same crack length (showing the corresponding displacement), with regions where d = 1 colored white. The stress–strain curves for the bicrystals with different angles are shown in (e).
shown in Fig. 7c. However, when the tensile load is applied to the material, the softer grain (θ = 0◦ ) will deform more than the θ = 90◦ grain. For this reason, it has a higher deformation energy and so the energy decreases faster through fracture in the θ = 0◦ grain. The stress–strain curves for both bicrystals are shown in Fig. 11c, and the two curves are nearly identical. This is because in both bicrystals, the crack is propagating through the same crystal orientation. The curves do not match those for the single crystal fracture with θ = 90◦ , nor do they match the single crystal curve with θ = 0◦ , even though the crack propagates through the θ = 0◦ grain. This is because the deformation energy throughout the bicrystal is impacted by both grains. To study the influence of crystal orientation on the crack path in bicrystals, we study four polycrystals in which the orientation of the top grain varies from θ = 30◦ to 90◦ , and the bottom grain has no rotation. The σ yy stress fields and crack paths are shown in Figs. 12a–12d. In each case, the crack propagates through the bottom, unrotated grain. In addition, as the orientation angle increases, the crack curves further from the grain boundary. Thus, as the difference in the elasticity tensor between the two grains increases, the energy drives the crack further away from the grain boundary. The stress–strain curves, shown in Fig. 12e change with crystal orientation, but not nearly to the magnitude observed in the single crystal simulations (Fig. 7c). 4.3.3. Fracture simulation in polycrystal structures As our final example, we investigate the impact of elastic anisotropy on the fracture behavior in 2D polycrystals. The polycrystal is composed of Cd, using the properties from Table 2, and has 36 grains generated using a Voronoi tessellation. The orientations of the grains are randomly generated, ranging from 0◦ to 90◦ , and three different sets of random orientations are generated to compare the impact on the crack behavior. The three polycrystals have an initial crack of length 0.3 mm and each have the same grain structure but different grain orientations as shown in Figs. 13a–13c. As with the bicrystal simulations, the grain boundaries behave in the same manner as the rest of the material. The polycrystals have zero x-displacements on the left side of the structure, zero y-displacements on the bottom surface, and they are loaded on the top surface with a positive y-displacement load with a rate of
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
13
Fig. 13. Polycrystal simulation initial conditions, where (a)–(c) show the crystal orientations and initial crack for the three different cases, and (d)–(f) shows the corresponding C1111 elasticity tensor component in GPa.
1 × 10−3 mm/s. To show the influence of crystal orientation on the elastic constants, the C1111 component of the heterogeneous elasticity tensor is shown in Figs. 13d–13f. The σ yy stress field and crack path for the three polycrystals are shown in Figs. 14a–14c. The impact of the grain orientations is clearly evident on the stress distribution, resulting in some deviation in the crack path near grain boundaries. In addition, the crack propagation rate is directly impacted by the crystal orientations, as seen in the release portion of the stress–strain curves (shown in Fig. 14d). Large changes in the decreasing slope result from the crack moving between grains with very different orientations. For example, this occurs very early in the crack propagation in Case I, it occurs late in the propagation in Case II, and it occurs twice in Case III. 5. Conclusions In this paper, we present a new phase field model of brittle fracture that includes the impact of anisotropic elastic behavior. The stress is decomposed into its tensile and compressive parts using a spectral decomposition, such that only the tensile energy contributes to fracture. To verify that it properly captures the behavior of isotropic materials, this new model was applied to mode I and mode II fracture in isotropic materials, and was shown to predict the expected crack paths. To investigate its ability to model fracture in anisotropic materials, fracture in single crystal, bicrystal, and polycrystal domains was simulated. Single crystal fracture with different crystal orientations was investigated in both hcp and fcc materials. In each case, the crack path was unaffected by the crystal orientation but the stress–strain curves were significantly impacted. The fcc material demonstrated cubic symmetry in its behavior. The bicrystal simulations showed that the fracture behavior was impacted by the orientations of both grains and that the bicrystal stress field diverted the crack path. The polycrystal simulations were conducted with three different
14
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643
Fig. 14. Results from polycrystal fracture simulations with different configurations of crystal orientations , where (a)–(c) show the postprocessed crack path and σ yy stress field in GPa for the different orientation distributions, with regions where d = 1 colored white. The resultant stress–strain curves are shown in (d).
sets of crystal orientations, and the final crack paths and stress–strain curves showed the impact of the different orientations. Acknowledgments Zhang and Tonks gratefully acknowledge financial support from the United States Department of Energy, Office of Nuclear Energy under Nuclear Energy University Program (NEUP) Grant 15-8243. Zhang and Jiang gratefully acknowledge financial support from the United States Department of Energy, Office of Nuclear Energy under Nuclear Energy University Program (NEUP) Grant 15-8229. References [1] M. Elices, G. Guinea, J. Gomez, J. Planas, The cohesive zone model: advantages, limitations and challenges, Eng. Fract. Mech. 69 (2) (2002) 137–163. [2] N. Moës, T. Belytschko, Extended finite element method for cohesive crack growth, Eng. Fract. Mech. 69 (7) (2002) 813–833. [3] A. Menk, S.P. Bordas, Numerically determined enrichment functions for the extended finite element method and applications to bi-material anisotropic fracture and polycrystals, Internat. J. Numer. Methods Engrg. 83 (7) (2010) 805–828. [4] G.A. Francfort, J.-J. Marigo, Revisiting brittle fracture as an energy minimization problem, J. Mech. Phys. Solids 46 (8) (1998) 1319–1342. [5] B. Bourdin, G.A. Francfort, J.-J. Marigo, Numerical experiments in revisited brittle fracture, J. Mech. Phys. Solids 48 (4) (2000) 797–826.
S. Zhang, W. Jiang and M.R. Tonks / Computer Methods in Applied Mechanics and Engineering 358 (2020) 112643 [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32]
15
I. Aranson, V. Kalatsky, V. Vinokur, Continuum field description of crack propagation, Phys. Rev. Lett. 85 (1) (2000) 118. A. Karma, D.A. Kessler, H. Levine, Phase-field model of mode III dynamic fracture, Phys. Rev. Lett. 87 (4) (2001) 045501. V. Hakim, A. Karma, Laws of crack motion and phase-field models of fracture, J. Mech. Phys. Solids 57 (2) (2009) 342–368. H. Henry, H. Levine, Dynamic instabilities of fracture under biaxial strain using a phase field model, Phys. Rev. Lett. 93 (10) (2004) 105504. C. Kuhn, R. Müller, A phase field model for fracture, PAMM 8 (1) (2008) 10223–10224. H. Amor, J.-J. Marigo, C. Maurini, Regularized formulation of the variational brittle fracture with unilateral contact: numerical experiments, J. Mech. Phys. Solids 57 (8) (2009) 1209–1229. C. Miehe, M. Hofacker, F. Welschinger, A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits, Comput. Methods Appl. Mech. Engrg. 199 (45) (2010) 2765–2778. C. Miehe, F. Welschinger, M. Hofacker, Thermodynamically consistent phase-field models of fracture: variational principles and multi-field FE implementations, Internat. J. Numer. Methods Engrg. 83 (10) (2010) 1273–1311. M.J. Borden, C.V. Verhoosel, M.A. Scott, T.J. Hughes, C.M. Landis, A phase-field description of dynamic brittle fracture, Comput. Methods Appl. Mech. Engrg. 217 (2012) 77–95. M.J. Borden, T.J. Hughes, C.M. Landis, C.V. Verhoosel, A higher-order phase-field model for brittle fracture: Formulation and analysis within the isogeometric analysis framework, Comput. Methods Appl. Mech. Engrg. 273 (2014) 100–118. P. Chakraborty, Y. Zhang, M.R. Tonks, Multi-scale modeling of microstructure dependent intergranular brittle fracture using a quantitative phase-field based method, Comput. Mater. Sci. 113 (2016) 38–52. P. Chakraborty, P. Sabharwall, M.C. Carroll, A phase-field approach to model multi-axial and microstructure dependent fracture in nuclear grade graphite, J. Nucl. Mater. 475 (2016) 200–208. V. Hakim, A. Karma, Crack path prediction in anisotropic brittle materials, Phys. Rev. Lett. 95 (23) (2005) 235501. A. Abdollahi, I. Arias, Phase-field modeling of fracture in ferroelectric materials, Arch. Comput. Methods Eng. 22 (2) (2015) 153–181. B. Li, C. Peco, D. Millán, I. Arias, M. Arroyo, Phase-field modeling and simulation of fracture in brittle materials with strongly anisotropic surface energy, Internat. J. Numer. Methods Engrg. 102 (3–4) (2015) 711–727. J. Clayton, J. Knap, A geometrically nonlinear phase field theory of brittle fracture, Int. J. Fract. 189 (2) (2014) 139–148. J. Clayton, J. Knap, Phase field modeling of directional fracture in anisotropic polycrystals, Comput. Mater. Sci. 98 (2015) 158–169. S. Teichtmeister, C. Miehe, Phase-field modeling of fracture in anisotropic media, PAMM 15 (1) (2015) 159–160. P. Shanthraj, B. Svendsen, L. Sharma, F. Roters, D. Raabe, Elasto-viscoplastic phase field modelling of anisotropic cleavage fracture, J. Mech. Phys. Solids 99 (2017) 19–34. T.T. Nguyen, J. Réthoré, M.-C. Baietto, Phase field modelling of anisotropic crack propagation, Eur. J. Mech. A Solids 65 (2017) 279–288. D. Gaston, C. Newman, G. Hansen, D. Lebrun-Grandie, Moose: A parallel computational framework for coupled systems of nonlinear equations, Nucl. Eng. Des. 239 (10) (2009) 1768–1778. B. Bourdin, G.A. Francfort, J.-J. Marigo, The variational approach to fracture, J. Elasticity 91 (1–3) (2008) 5–148. C. Miehe, Comparison of two algorithms for the computation of fourth-order isotropic tensor functions, Comput. Struct. 66 (1) (1998) 37–43. C. Miehe, M. Lambrecht, Algorithms for computation of stresses and elasticity moduli in terms of seth–hill’s family of generalized strain tensors, Int. J. Numer. Methods Biomed. Eng. 17 (5) (2001) 337–353. V. Lubarda, D. Krajcinovic, S. Mastilovic, Damage model for brittle elastic solids with unequal tensile and compressive strengths, Eng. Fract. Mech. 49 (5) (1994) 681–697. S. Murakami, Continuum Damage Mechanics: a Continuum Mechanics Approach to the Analysis of Damage and Fracture, Vol. 185, Springer Science & Business Media, 2012. M.R. Tonks, D. Gaston, P.C. Millett, D. Andrs, P. Talbot, An object-oriented finite element framework for multiphysics phase field simulations, Comput. Mater. Sci. 51 (1) (2012) 20–29.