A phase field model for fracture based on the strain gradient elasticity theory with hybrid formulation
Introduction
Fracture is one of the most common failure modes of engineering materials and structural components. In the process of engineering structural designs, the prevention of crack-induced fracture is a major concern. It is an important subject to predict crack nucleation and propagation with complex geometric configurations and loading modes which can help the designers to gain insights into the fracture mechanisms of the structural components and improve the structural design.
In recent years, the phase field (PF) model for fracture has emerged as a robust numerical approach to simulate complex fracture behaviors in engineering materials and structures [1], [2], [3], [4], [5]. The PF model for fracture is based on a regularization of the variational formulation of brittle fracture proposed by Francfort and Marigo [6], where the sharp crack surface is approximated as a diffusive crack zone [7], [8] by introducing an auxiliary scalar variable (or phase field variable). The fracture energy is characterized by a volume integral in the diffusive crack zone in the sense of Gamma convergence [7]. Crack initiation and propagation are then regarded as an energy minimization process with respect to the PF field variable and displacement field subject to certain boundary conditions [6]. Compared with previous numerical methods based on the discrete crack model like the cohesive zone model (CZM) and the extended finite element method (XFEM), the PF approach needs no extra ad-hoc criteria to simulate crack nucleation, propagation, deflection and branching as well as the coalescence of multiple cracks [9], [10], [11], [12], [13], [14], which is especially appealing to model fracture behavior with complex crack patterns.
The regularized model proposed by Bourdin et al. [15], [16] predicts symmetric fracture behavior in traction and in compression and allows for material interpenetration in cracks under compressive actions. In order to avoid unrealistic fracture under compression, various strain energy decomposition methods have been proposed to differentiate the deformation state [1], [3], [5], [9], [17]. For example, Amor et al. [17] decomposed the strain energy density into the volumetric and deviatoric contributions. Only the strain energy associated with expansion and shear contributes to the evolution of the crack phase. Miehe et al. [3], [9] adopted the spectral decomposition of the strain tensor to split the strain energy density, where the positive tensile part degrades with the increasing phase-field. In addition, a hybrid formulation was proposed by Ambati et al. [5] to dealing with the asymmetric tensile and compressive behavior for quasi-static brittle fracture. In the hybrid formulation, the constitutive relation between stress and strain is still derived from the isotropic elastic energy density but the evolution of the crack phase field is driven merely by the tensile elastic energy to avoid cracking in the compressed regions. Though the hybrid phase field formulation is variationally inconsistent, it does not violate the second law of thermodynamics [4]. Within a staggered implementation, the hybrid formulation is much more computationally efficient compared with the anisotropic PF models [5]. Recently, Wu introduced a modified reference energy in terms of an equivalent effective stress, where the ratio between the uniaxial compressive strength and the tensile strength is characterized by a material parameter. Zhou et al. [1] proposed a novel PF model to simulate compressive-shear fractures in rock-like materials, where a new driving force was constructed by considering the influence of cohesion and internal friction angle. With the above strain energy decomposition methods, the asymmetric fracture behavior under tension and compression can be well characterized. In the past decades, the PF model for fracture has received fast development and has been widely adopted to study complex fracture behaviors in brittle materials [3], [5], ductile materials [13], [18], heterogeneous materials [19], hyperelastic materials [20], [21], polycrystals and phase transformation materials [22], [23], [24], [25], etc.
Although the PF method has been shown to be successful to predict complex fracture behavior in various materials, till now, most of the PF models are based on the classical Cauchy continuum theory, which neglects the influence of the underlying microstructure. The Cauchy continuum theory is developed to describe the macroscopic material responses on condition that the specimen size is much larger than the internal micro-structure [26]. The size effect in micron- and nano-dimensions cannot be captured with the classical continuum theory due to the lack of a length scale parameter. In addition, the presence of stress singularity at the crack tip is also known as one of the drawbacks of the classical continuum theory [27], [28], [29], [30]. To remedy these drawbacks, higher order continuum theories were proposed to characterize the effect of microstructure. Mindlin [31], [32], [33] developed a general gradient elasticity theory where eighteen new constants are introduced to reconcile the dimensions of strains and higher order gradients of strains and to correlate the micro-strains with macro-strains. This gradient theory was later simplified into three different forms where only two Lamé constants and other five material parameters are involved [33], [34]. In Form-I, the strain energy density is a quadratic function of the classical strains and the second gradient of displacement; in Form-II, the second displacement gradient is substituted by the strain gradient and in Form-III the strain energy function is expressed as a function of the strain, the symmetric part of the strain gradient and the gradient of rotation [34]. Although the final governing partial differential equations, written in terms of the displacement field, are identical, Form-II leads to a symmetric total stress tensor [34]. Aifantis [35] and Ru and Aifantis [36] further provided a simplified version of Mindlin’s theory in Form-II with only one length-scale parameter besides the Lamé constants.
The simplified strain gradient elasticity theory is widely adopted to characterize the size effects of the mechanical responses of miniaturized structural components [37], [38], [39]. It has been also shown that the strain gradient theory can be applied to homogenize the micro-structure of metamaterials [40], [41], [42]. Besides that, the gradient elasticity theory provides a regularization of the governing equations of linear elasticity [43], [44]. In the framework of gradient elasticity, the stress singularity around micro-defects or cracks predicted by linear elasticity can be mitigated or eliminated [28], [29], [45]. The first gradient theory has been applied by many researchers to characterize the stress field of Mode I, II, III cracks [28], [34], [46], [47], [48], [49], [50], [51]. Although there are still some debates on the exact form of the analytical solution of the stress field near the crack tip, recent theoretical and numerical studies have shown that the gradient elasticity theory predicts non-singular stress field near 2D crack tips [28], [45]. Giannakopoulos and Stamoulis [51] derived the energy release rate of a gradient elastic cracked bar under uniaxial tension by adopting Form II strain gradient theory. They found that the normalized strain energy release rate decreases monotonically with the length scale parameter, which implies a toughening effect of the gradient elasticity.
Recently, Placidi and Barchiesi [52], Makvandi et al. [27] attempted to develop PF fracture models based on the strain gradient elasticity theory. It has been shown that the stress singularity at the crack tip can be eliminated by adopting the strain gradient elasticity theory. In the meanwhile, the numerical solution of the stress field around the crack tip exhibits less mesh sensitivity compared with the linear elastic case. This character is very appealing from the aspect of numerical simulation as the mesh size requirement is less strict in the crack tip region. However, in these studies, the total elastic strain energy, which is a function of the strain and strain gradient, degrades with the PF variable without discriminating the stress state. As a result, the PF models proposed by Placidi and Barchiesi [52] and Makvandi et al. [27] are only applicable for modelling fracture behavior under tensile stresses. Actually, only Mode-I fracture behavior was discussed in their studies. Proper characterization of the asymmetry fracture behavior under tension and compression is a prerequisite for modelling crack propagation in gradient medium with complex geometries. Currently, this topic has not been tackled yet.
In this paper, a novel PF fracture model is developed in the framework of strain gradient elasticity. Inspired by the strain energy decomposition methods proposed for PF models of linear elastic fracture, four different strain energy decomposition methods are proposed for the gradient medium. These four methods are compared and discussed by studying Mode I fracture. Eventually, the energy decomposition method based on the spectral decomposition of the strain and the strain gradient is chosen to characterize the asymmetric fracture behavior under tension and compression. The proposed PF model is based on the hybrid formulation, which is computationally efficient. In the hybrid formulation, the stress tensor, higher order stress tensor, strain and strain gradient tensor are still derived from the isotropic strain energy. But the phase field evolution has considered the energy decomposition. Although the hybrid PF formulation is variationally inconsistent, it obeys the second law of thermodynamics [4]. To the best knowledge of the authors, such a hybrid PF model is the first time to be presented in the framework of gradient elasticity and applied to predict the crack propagation trajectory.
For PF models developed based on the strain gradient elasticity, higher-order derivatives of displacement are involved in the governing equations. Conventional finite elements satisfying C0 continuity across the element boundary are not appropriate to solve equations with higher order derivatives. As the first gradient theory is adopted in this study, C1-continous element is needed to solve the governing partial differential equations. Another possible approach to solve partial differential equations with higher order derivatives is to use the so-called mixed finite element method [53], [54]. However, the computational cost of this approach is usually very high. In this paper, the higher order partial differential equations are solved through Abaqus subroutine UEL (user defined element) with nine-node quadrilateral C1-continous element.
The content of this paper is structured as follows. The hybrid PF model based on the strain gradient elasticity for brittle fractures is described in Section 2. Subsequently, in Section 3, the finite element discretization of the governing equations is presented. Section 4 describes the numerical implementation of the proposed PF model based on strain gradient elasticity in Abaqus through user defined element (UEL). The numerical implementation of the PF model as well as the comparison of different energy decomposition methods are presented in Section 5. In Section 6, the rationality of the proposed PF model is demonstrated by studying representative numerical examples. Finally, conclusions are drawn in Section 7.
Section snippets
Phase field approximation of the smeared crack topology
The PF fracture model is based on the smeared crack model. For the purpose of demonstration, a simple one-dimension model is considered. In the 1-D setting, an infinite bar with cross section occupies the domain where . The cross section at the position represents the crack surface. As shown in Fig. 1(a), an auxiliary phase field (PF) variable is used to describe the sharp crack withwhere the d = 0 and d = 1 represent the unbroken state and the
Finite element formulation
In the present work, the finite element method (FEM) is used to solve the governing equations. The weak form of Eq. (30) can be written as
We use the standard vector–matrix notation and denote the nodal values of the displacement and phase field as and , respectively. The discretization is given bywhere and are the vectors consisting of node values and ; and are shape function matrices. The
Numerical implementation through Abaqus user subroutine UEL
To take advantage of the powerful nonlinear solver in Abaqus[57], [58], the PF fracture model proposed above is implemented through Abaqus subroutine UEL (user defined element) with nine-node quadrilateral C1-continous element. In addition, each node has three degrees of freedom (i.e.,and). By using the staggered solution method, the displacement field and phase field minimization problems are solved separately. The associated nonlinear equations are solved iteratively with the
Crack tip stress field
In this section, the stress field near a static crack tip is examined to validate the numerical implementation of the finite element formulation. The influence of the mesh size and length scale parameter of gradient elasticity on the crack tip stress field is also studied.
As shown in Fig. 4(a), a squared plate of side 1 mm containing an edge crack of length 0.5 mm is considered. The geometric setup and the boundary conditions are depicted in Fig. 4(a). In order to bypass the influence of the PF
Numerical examples based on Method I
In this section, some benchmark numerical examples are presented to demonstrate the suitability of Method I to characterize the fracture behavior of the gradient medium. The influence of the length scale parameter of the gradient elasticity on the fracture behavior is also discussed. These numerical examples include the tensile and shear fracture of single edge-cracked specimen, the symmetric three point bending test and the fracture of a notched plate with hole under tensile load.
Conclusions
In this work, a PF model based on the strain gradient elasticity with a hybrid formulation is proposed to study the fracture behavior in the gradient medium. A strain energy decomposition method is proposed to properly characterize the asymmetric fracture behavior under tension and compression. The numerical accuracy and validity of the proposed PF model is demonstrated by studying the stress field near the crack tip and several benchmark problems in the open literature. The major conclusions
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgments
This work is supported by the National Natural Science Foundation of China (Grant Nos. 11872023 and 11372114).
References (64)
- et al.
Phase field modeling of brittle compressive-shear fractures in rock-like materials: A new driving force and a hybrid formulation
Comput Methods Appl Mech Eng
(2019) - et al.
A fracture-controlled path-following technique for phase-field modeling of brittle fracture
Finite Elem Anal Des
(2016) - et al.
A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits
Comput Methods Appl Mech Eng
(2010) - et al.
A length scale insensitive phase-field damage model for brittle fracture
J Mech Phys Solids
(2018) - et al.
Revisiting brittle fracture as an energy minimization problem
J Mech Phys Solids
(1998) - et al.
A two-set order parameters phase-field modeling of crack deflection/penetration in a heterogeneous microstructure
Comput Methods Appl Mech Eng
(2019) - et al.
An open-source Abaqus implementation of the phase-field method to study the effect of plasticity on the instantaneous fracture toughness in dynamic crack propagation
Comput Method Appl Mech Eng
(2020) - et al.
On the crack onset and growth in martensitic micro-structures; a phase-field approach
Int J Mech Sci
(2021) - et al.
Numerical experiments in revisited brittle fracture
J Mech Phys Solids
(2000) - et al.
Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments
J Mech Phys Solids
(2009)
Phase field fracture in elasto-plastic solids: Variational formulation for multi-surface plasticity and effects of plastic yield surfaces and hardening
Int J Mech Sci
A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure
Eng Fract Mech
Large strained fracture of nearly incompressible hyperelastic materials: enhanced assumed strain methods and energy decomposition
J Mech Phys Solids
A length scale insensitive anisotropic phase field fracture model for hyperelastic composites
Int J Mech Sci
Ferroelastic toughening of single crystalline yttria-stabilized t’zirconia: A phase field study
Eng Fract Mech
Study of crack propagation behavior in single crystalline tetragonal zirconia with the phase field method
Eng Fract Mech
Study of transformation induced intergranular microcracking in tetragonal zirconia polycrystals with the phase field method
Mater Sci Eng A
A phase-field fracture model based on strain gradient elasticity
Eng Fract Mech
Approximation of mode I crack-tip displacement fields by a gradient enhanced elasticity theory
Eng Fract Mech
Second gradient of strain and surface-tension in linear elasticity
Int J Solids Struct
On first strain-gradient theories in linear elasticity
Int J Solids Struct
On the role of gradients in the localization of deformation and fracture
Int J Eng Sci
Gradient elasticity in statics and dynamics: an overview of formulations, length scale identification procedures, finite element implementations and new results
Int J Solids Struct
Longitudinal vibrations of a beam: a gradient elasticity approach
Mech Res Commun
On the gradient strain elasticity theory of plates
Euro J Mech A/Solids
Modelling size-dependent bending, buckling and vibrations of 2D triangular lattices by strain gradient elasticity models: applications to sandwich beams and auxetics
Int J Eng Sci
Form II of Mindlin’s second strain gradient theory of elasticity with a simplification: For materials and structures from nano-to macro-scales
Euro J Mech A/Solids
Asymptotic analysis in Linear Elasticity: From the pioneering studies by Wieghardt and Irwin until today
Eng Fract Mech
Scaling of volume energy density function reflecting damage by singularities at macro-, meso-and microscopic level
Theoret Appl Fract Mech
On non-singular GRADELA crack fields
Theor Appl Mech Lett
Gradient elasticity with surface energy: Mode-I crack problem
Int J Solids Struct
On the structure of the mode III crack-tip in gradient elasticity
Scripta Metall Mater
Cited by (8)
The implicit stabilized dual-horizon peridynamics-based strain gradient damage model
2024, Applied Mathematical ModellingPhase field study of crack growth in t′ yttria stabilized zirconia with initial domain structures
2023, Materials Today CommunicationsPhase field study of the thermo-electro-mechanical fracture behavior of flexoelectric solids
2023, Theoretical and Applied Fracture MechanicsA double-phase field method for mixed mode crack modelling in 3D elasto-plastic solids with crack-direction-based strain energy decomposition
2023, Computer Methods in Applied Mechanics and EngineeringA thermodynamically consistent phase-field regularized cohesive fracture model with strain gradient elasticity and surface stresses
2022, Engineering Fracture MechanicsCitation Excerpt :In addition, Rao et al. [44] rigorously derived a new macroscopic brittle fracture model for materials with micro-cracks inside using a two-scale asymptotic analysis procedure, and strain gradient elasticity is naturally involved in the model in the absence of any phenomenological assumptions. Furthermore, in materials with many microstructures, strain gradient elasticity is included in fracture criterion or damage evolution equations [22,45–52]. Mandal et al. [16] proposed a PFM for brittle fracture with strain gradient elasticity for Mode I fracture.
A phase field model for electromechanical fracture in flexoelectric solids
2022, Engineering Fracture MechanicsCitation Excerpt :Thus the results presented in Fig. 8 might be universal, i.e., the flexoelectricity around the crack tip can have significant influences on the macro fracture behavior of dielectric materials. In our previous study [67], it has been shown that the stress concentration near the crack tip is more intensive when the length scale parameter l takes a small value, for example, 10 nm. In this case, the numerical solutions of the crack tip fields are more sensitive to the mesh size.