A phase field model for fracture based on the strain gradient elasticity theory with hybrid formulation

https://doi.org/10.1016/j.engfracmech.2021.107975Get rights and content

Highlights

  • A PF fracture model is developed based on the strain gradient elasticity theory.

  • An energy decomposition method is proposed for the new PF fracture model.

  • The PF model is numerically implemented through Abaqus subroutine UEL.

  • The proposed PF model can characterize the toughening effect of gradient elasticity.

  • The proposed PF model can well predict the curvilinear crack propagation path.

Abstract

In this paper, a novel phase field (PF) model for fracture is developed in the framework of strain gradient elasticity. The strain energy decomposition methods initially proposed for linear elastic fracture problems are extended to the gradient elasticity situation. The PF model is numerically implemented through Abaqus subroutine UEL (user defined element) with nine-node quatrilateral C1-continous element. The finite element implementation is validated by studying the stress field near the static crack tip. When the length-scale parameter of the gradient elasticity is zero, the stress field is consistent with the analytical solution predicted by linear elastic fracture mechanics (LEFM). For non-zero length-scale parameters, the Cauchy stress at the crack tip is found to be non-singular and the numerical results exhibit less mesh sensitivity in the crack tip region compared with the linear elastic case. After that, four different strain energy decomposition methods (Method I-IV) are compared and discussed by studying Mode I fracture behavior under tensile load. It is found that Method I can properly characterize the toughening effect of gradient elasticity. Based on Method I, the influences of the fracture length-scale parameter lc and the strain gradient length-scale parameter l on the characteristics of the smeared crack model and the load–displacement curves are systematically discussed. The specimen size effect of Mode I crack propagation is also investigated. It is found that the length-scale parameter l has a larger influence on the load–displacement curve as the specimen size decreases. Finally, crack propagation behavior in pure shear test, three point bending test and tensile test of a notched plate with hole is investigated. It is demonstrated that the proposed PF model can well predict the curvilinear crack propagation path and properly characterize the effect of gradient elasticity. Thus, the PF model developed in this paper can be applied to study complex fracture behaviors in the framework of gradient elasticity, where the effect of internal structure on the mechanical responses become significant.

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 Ω=Γ×L where L=[-,+]. The cross section Γ at the position x=0 represents the crack surface. As shown in Fig. 1(a), an auxiliary phase field (PF) variable d[0,1] is used to describe the sharp crack withd=1(ifx=0)0(ifx=1)where 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δL=Ω[(1-d)2+k](-δεTDε-l2δηTDηη)+2(1-d)Hδd-Gc(lcd·δd+1lcdδd)dΩ+Ωb*δudΩ+ΩTt¯*δudS=0

We use the standard vector–matrix notation and denote the nodal values of the displacement and phase field as ui and d, respectively. The discretization is given byu=Nuû,d=Ndd̂where 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.ux,uyandd). 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)

  • J. Fang et al.

    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

    (2019)
  • T.T. Nguyen et al.

    A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure

    Eng Fract Mech

    (2015)
  • J.-Y. Ye et al.

    Large strained fracture of nearly incompressible hyperelastic materials: enhanced assumed strain methods and energy decomposition

    J Mech Phys Solids

    (2020)
  • T.K. Mandal et al.

    A length scale insensitive anisotropic phase field fracture model for hyperelastic composites

    Int J Mech Sci

    (2020)
  • Y. Sun et al.

    Ferroelastic toughening of single crystalline yttria-stabilized t’zirconia: A phase field study

    Eng Fract Mech

    (2020)
  • T. Zhao et al.

    Study of crack propagation behavior in single crystalline tetragonal zirconia with the phase field method

    Eng Fract Mech

    (2016)
  • J. Zhu et al.

    Study of transformation induced intergranular microcracking in tetragonal zirconia polycrystals with the phase field method

    Mater Sci Eng A

    (2017)
  • R. Makvandi et al.

    A phase-field fracture model based on strain gradient elasticity

    Eng Fract Mech

    (2019)
  • P. Isaksson et al.

    Approximation of mode I crack-tip displacement fields by a gradient enhanced elasticity theory

    Eng Fract Mech

    (2014)
  • R.D. Mindlin

    Second gradient of strain and surface-tension in linear elasticity

    Int J Solids Struct

    (1965)
  • R.D. Mindlin et al.

    On first strain-gradient theories in linear elasticity

    Int J Solids Struct

    (1968)
  • E.C. Aifantis

    On the role of gradients in the localization of deformation and fracture

    Int J Eng Sci

    (1992)
  • H. Askes et al.

    Gradient elasticity in statics and dynamics: an overview of formulations, length scale identification procedures, finite element implementations and new results

    Int J Solids Struct

    (2011)
  • B.S. Altan et al.

    Longitudinal vibrations of a beam: a gradient elasticity approach

    Mech Res Commun

    (1996)
  • K.A. Lazopoulos

    On the gradient strain elasticity theory of plates

    Euro J Mech A/Solids

    (2004)
  • S. Khakalo et al.

    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

    (2018)
  • S. Khakalo et al.

    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

    (2018)
  • A. Carpinteri et al.

    Asymptotic analysis in Linear Elasticity: From the pioneering studies by Wieghardt and Irwin until today

    Eng Fract Mech

    (2009)
  • G.C. Sih et al.

    Scaling of volume energy density function reflecting damage by singularities at macro-, meso-and microscopic level

    Theoret Appl Fract Mech

    (2005)
  • E.C. Aifantis

    On non-singular GRADELA crack fields

    Theor Appl Mech Lett

    (2014)
  • G. Exadaktylos

    Gradient elasticity with surface energy: Mode-I crack problem

    Int J Solids Struct

    (1998)
  • S.B. Altan et al.

    On the structure of the mode III crack-tip in gradient elasticity

    Scripta Metall Mater

    (1992)
  • Cited by (8)

    • A thermodynamically consistent phase-field regularized cohesive fracture model with strain gradient elasticity and surface stresses

      2022, Engineering Fracture Mechanics
      Citation 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 Mechanics
      Citation 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.

    View all citing articles on Scopus
    View full text