A simple and unified implementation of phase field and gradient damage models

In this work, the analogous treatment between coupled temperature–displacement problems and material failure models is explored within the context of a commercial software (Abaqus®). The implicit gradient Lemaitre damage and phase field models are implemented utilizing the software underlying capabilities for coupled temperature–displacement problems. The heat conduction equation is made compatible with the diffusive regularization of such material models and calculations are carried out at the material point level. This bypasses the need to implement explicitly the weak form resultant from the coupling between the momentum conservation and the evolution of the diffusive field. Throughout benchmarking examples, the proposed methodology is assessed and validated by investigating typical issues risen from the considered local inelastic-based deformation models, such as mesh dependency and the difficulties to predict cracked regions.


Introduction
The emergence of the so-called 'regularized' solutions for damage and failure in engineering materials has evolved considerably in recent years. They constitute appealing simplifications of the microstructural complexity, usually resorting to non-local and gradient theories. To ensure well-posedness of the set of partial differential equations to be solved, which is affected by the presence of softening at the constitutive level, these solutions abandon the principle of local action by coupling a diffusion equation of non-local variables with the momentum balance equation. The resultant procedure acts as a spatial regularization of local constitutive relations for failure, usually with the recourse of a characteristic length associated with the size of the non-local support region definition. Herein, two types of regularization are explored, namely: the non-local implicit gradient of the Lemaitre damage model [1] and the phase field model [2] that approximate crack evolution under inelastic deformations.
Prediction of initiation and propagation of cracks defines a highly crucial task in behavior assessment of ductile metallic structures, which requires the development of accurate and robust constitutive models that trigger the formation of meso-cracks initiated by the material internal progressive degradation due to nucleation, growth, and coalescence of microvoids. Within the context of coupled models, which take into account that degra-et al. [35], the phase field elastoplastic models by Duda et al. [36] and Ambati et al. [37], the crystal plasticity phase field coupled model of Padilla et al. [38], the phase field based anisotropic damage model of Mozaffari et al. [39] and the variational based, large deformation plasticity-phase field model by Miehe et al. [40,41] and rate-dependent plasticity phase field model of Badnava et al. [42] mostly concentrated on approximating the postpeak material behavior via a phase field methodology. Also in [43], de Borst et al. studied several similar and different aspects of gradient enhanced damage-based approaches and the phase field framework and a discussion has made over the solution of the broadening of the damaged zone in the wake of the crack tip with these methods.
The present study addresses a methodology for solving the diffusion equation of the phase field model and gradient-enhanced non-local damage model, in the context of an existing commercial software. It takes the advantage of built-in heat equation solver in an existing software code, Abaqus/Standard in the present case, bypassing the need to establish explicitly the weak form of the governing equations containing the momentum balance and the diffusion equation to solve separately for displacement and non-local or phase field variables. In fact, this procedure could be used with other commercial softwares with similar capabilities. This work is organized as follows: "Formulation" section builds the relation between the heat conduction problem and diffusion equation of gradient problems, and explains the gradient non-local model based on Lemaitre phenomenological ductile damage model and the phase field model. Numerical implementation and results of some illustrative examples are presented in "Results and discussion" section, and finally conclusions are drawn in "Concluding remarks" section and a shortened version of routine to implement the non-local gradient damage model is provided.

Formulation
In both gradient damage and phase field models, the mechanical equilibrium is performed for the coupling between displacements u and a field variable ψ. In a body , the strong form of the boundary value problem can be written, in general, as where σ is the stress tensor and b is the body force per unit of volume, with the following boundary conditions: being ∂ = ∂ u ∪ ∂ t the body surface boundary, n outward unit normal vector to the boundary surface at ∂ t , t n the traction force, u d prescribed displacements at the boundary surface ∂ u and a a set of possible material and numerical parameters. Second part of Eq. (1) depends on a "diffusion" coefficient Γ (a), on the "flux", or gradient, term ∇ψ, and on a "source" term h s (u, ψ, a), and is similar to steady-state heat conduction equation in solids. The identification of the different parameters of the two models with those in the heat conduction equation can be found in Table 1 and will be clarified in the following sections. The fact that almost always this equation is implemented in common Table 1 Diffusion equation forms for coupled temperature displacement, gradientdamage, and phase-field models   Name Field variable "Diffusion" coefficient Γ "Flux" term ∇ψ "Source" term h s Thermal conductivity, h s heat source, l n characteristic length for gradient damage model, G c fracture energy, H elastic energy density, l d characteristic length for phase field-model commercial software packages in the computational mechanics field hints that it may be used in establishing gradient and phase field models within these frameworks, taking advantage of the software built-in implementations. In Table 1, the parameters l n and l d are respectively length parameters associated to the gradient damage and phase field approaches, commonly denoted as the characteristic length scale of the material and are mostly considered as numerical parameters. Since its meaning with the microstructure is inherently complex, its precise determination usually relies on inverse methods [44].

Lemaitre damage model coupled with plasticity
The constitutive modeling of the fully-coupled damage-plasticity model in [1] is recovered in this section, considering its extensive application in describing ductile damage in metallic materials. This model is often categorized within the Continuous Damage Mechanics material behavior description as the damage is phenomenologically introduced at the macroscopic constitutive level by internal variables whose evolution mimic the nucleation, growth and coalescence of internal micro-voids. Mostly it is based on the hypothesis of strain equivalence and the associated concept of effective stress. Herein, the simplified version of the Lemaitre thermodynamically-consistent damage model proposed by De Souza Neto et al. [45], in the context of small strains, is adopted for simplicity. The decomposition of total strain tensor into elastic (ε e ) and plastic (ε p ) contributions yields: Accordingly, the free energy can be stated as the sum of the elastic damage ψ ed and plastic (ψ p ) potentials: where α is an isotropic hardening (internal) variable and D is the damage (internal) variable defined as a scalar, 0 ≤ D ≤ 1, evolving continuously from zero (virgin material) to one (fully damaged material). Assuming the elastic damage contribution as: the constitutive law reads: where σ is the Cauchy stress tensor and D e is the isotropic fourth order elasticity tensor, which may be defined as: where G and K are respectively shear and bulk moduli, I is the second order identity tensor, I d = I s − 1 3 I ⊗ I is the fourth order deviatoric projection tensor and I s is the symmetric identity tensor, given by: The yield condition to define the plastically admissible stress states is defined as: whereσ (s) = 3J 2 (s) is the von Mises equivalent stress with J 2 (s) as the second invariant of the deviatoric stress tensor, s. The yield stress σ y (α) is a function of the internal hardening variable, α, which is considered herein as the equivalent plastic strain,ε p , written as:ε whereε p is the rate of the plastic strain. The associative plastic flow rule is assumed and then:ε whereγ is the plastic multiplier and N is the flow vector, defined as: These relations are supplemented by the definition of damage-related energy derived from the total potential as: where Y is the damage energy release rate which can be alternatively written as: being p = 1 3 tr (σ) the hydrostatic pressure. The evolution of the internal variables reads: where r 1 and r 2 are material parameters which can be calibrated by experimental observations. Finally, the usual loading/unloading conditions of rate-independent plasticity follow as:

Gradient regularization of the local damage model
The non-local implicit gradient formulation is here recovered. Based on this approach, a non-local fieldf is related to the corresponding local field f , as in Eq. (1), through the equation: This approach bypasses the calculation of second order gradients of the local field of explicit nonlocal models, and establishes a more robust solution. Here, we focus on the gradient regularization of Lemaitre based local damage field, resorting to an alternative numerical treatment in the context of the treatment of multi-field problem within a commercial finite element code Abaqus/Standard. As a point of departure, the following set of governing equations defines the implicit gradient elastoplastic framework coupled with the internal damage: which includes the stress equilibrium and the additional partial differential equation corresponds to the implicit gradient regularization of the local damage field, D. In Eq. (18) 2 , D = ∇ 2 D denotes the Laplacian of the nonlocal damage field. With this strong form at hand, the following homogenous Neumann boundary conditions are prescribed: Unlike the explicit gradient formulation, this boundary condition is homogenous and does not require the gradient of the local field as an additional nodal degree of freedom. Following the constitutive relations of the previous section, the yield function of Eq. (9) can be rewritten based on the nonlocal field: Accordingly, the flow vector may be modified as: This is followed by the rate forms of the internal hardening variable and damage parameter: together with the damage driving force written as:

Integration algorithm of the non-local resolution
The integration procedure of the explained non-local gradient model is based on the conventional elastic predictor-plastic corrector, which is implemented via an operatorsplit scheme. The solution of the multi-field problem of Eq. (18) consists of the solution of mechanical problem based on the updated stresses and the local damage field, and the calculation of the nonlocal field from the gradient equation. The plastic correction in the mechanical problem is performed via the known iterative Newton Raphson algorithm to solve the return mapping equation, as is proved to have quadratic convergence rates. Given ε and nonlocal damage variable from the previous increment (D n ) and considering the pseudo time interval t ∈ (t n , t n+1 ), elastic predictor stage gives the following elastic trial states: where ε e trial v n+1 and ε e trial d n+1 are, respectively, the trial values of elastic volumetric strain and elastic deviatoric strain tensor: The plastic admissibility is evaluated by the following incremental yield function: If this condition is satisfied, the elastic stage proceeds by taking the trial predictor states as the updated values: Otherwise, the single equation return mapping algorithm is activated and triggers the solution of the following equation to obtain the plastic multiplier, γ : considering the so-called material integrity in the previous and current increments as: The damage energy release rate is given by: The Newton Raphson solution strategy is utilized to solve the Eq. (29) iteratively, taking the following initial value for the plastic multiplier rate: By achieving the prescribed tolerance, the state variables are updated based on the new material integrity: Using Abaqus/Standard as a support software, the above integration algorithm is coded into a UMAT subroutine to reach the updated stress and local damage fields. This is processed based on the nonlocal damage fields obtained from the solution of diffusion Eq.
with "thermal" constants (specific heat c p and density ρ), "pseudo time" and time increment so that the steady-state conditions are closely attained in each load increment. Given the local damage field, the nonlocal solution is achieved by the following assumptions (see Table 1): The layout of the multifield algorithm is represented in Fig. 1. Accordingly, the mechanical problem consisting of calculation of stresses and Jacobian is performed with the UMAT code and the solution dependent variables (SDVs) are stored in array STATEV (including temperature). In the subsequent increment, this data is used in the HETVAL routine, whereas the flux is defined within the FLUX array based on the values of local damage and nonlocal damage (temperature) to compute the updated nonlocal field. In "Appendix", a brief version of UMAT and HETVAL routines to implement the above gradient damage model are provided with concentration on the definition of nonlocality in the algorithm definition.

Phase field approximation in fracture
Recently, in the context of fracture analysis of structures, the phase field approach appeared as an alternative to classical discrete crack approaches due to their intrinsic difficulty in dealing with more complex crack topologies such as kinking or branching patterns. This methodology utilizes an order parameter called phase field, which distinguishes between the intact and fully-damaged regions of the deformed body, using physical phase transformation concepts and conceptually it may be more closely related to continuous evolution of internal degradation in damage mechanics. In its implementation, this method avoids the necessity to deal with displacement jumps due to a sharp description of cracks in discrete-based approaches, as all the calculations are performed on the original mesh without any ad hoc criteria. Adding to these the straightforward extensibility of the problem to higher dimensions in the same manner of lower dimensions, this approach has been widely used in the fracture analyses and more recently in prediction of evolution of inelastic material imperfections. As for the numerical implementation, the standard finite element shape functions can be utilized to interpolate the displacement and the diffusive field which, comparatively, bypasses the complexities of employment of enriched shape functions in methods such as XFEM. Also the treatment of the diffusion equation can be closely related to non-local gradient approach as both models use spatial gradients of the regularizing field.

Phase field model of brittle fracture
Let ⊂ R n be the reference configuration of an arbitrary body with n as the space dimensions, ∂ = ∂ u ∪ ∂ t ⊂ R n−1 be its boundary and Γ ⊂ R n−1 as the sharp crack discontinuity as is schematically depicted in Fig. 2. According to the classical Griffith theory of brittle fracture and assuming small and quasi-static deformations, the total potential may be expressed as the sum of the elastic strain energy, external work and fracture energy given by where ψ e is the elastic energy density function of the second order infinitesimal strain tensor, ε (u) = 1 2 ∇u + ∇ T u , and G c is the threshold value of energy releasing from the formation of meso-cracks. The elastic energy density can be written as: with the elastic stiffness tensor D e and without loss of generality of Eq. (7), and λ and μ as the Lamé constants, which are defined in terms of the elastic (E) and shear (μ) moduli: The phase field methodology may be regarded as a specific case of non-local gradient model, in which the regularization is performed on sharp crack interfaces with a pure geometrical representation (Fig. 2b), as the sharp discontinuity of Fig. 2a) is smoothened by an auxiliary phase field variable d ∈ [0, 1], discriminates between the intact (d = 0) and fully-broken (d = 1) phases. The width of the regularized region is driven by a length scale l d , which may be determined based on the micromechanical observations such as particle shape, size or strength. The sharp discontinuity pattern is recovered for the vanishing length scale. This crack treatment was emerged based on the variational formulation of Francfort and Marigo [21], subsequently complemented by the regularization of Bourdin et al. [22], whereas the classical surface integral of Eq. (36) is replaced by a volume integral, reads the following potential: One commonly used form of the stress degrading function g (d) is the following monotonically decreasing function: where η 1 is a dimensionless residual stiffness at the total failure to avoid numerical difficulties. This function should satisfy the following properties [25]: Based on the work of Miehe et al. [26], the strain tensor can be additively decomposed to the positive (tensile) and negative (compressive) parts in order to avoid cracking under compression: with ε a and n a as the eigenvalues and eigenvectors of ε in n dimensions, and ε a ± = 1 2 (ε a ± |ε a |). Accordingly, the degraded isotropic elastic energy function can be decomposed to positive and negative parts: with ψ ± e (ε) = λ 2 tr ε ± 2 + μtr ε ± 2 , where the negative part remains undegraded. This leads to the following constitutive law: where σ ± 0 = λ tr (ε) ± + 2μ ε a ± n a ⊗ n a and ∂ x (.) is the derivative of function (.) with respect to x. Applying Eqs. (40) and (43) to the potential in Eq. (39), the modified functional may be written as: According to the history-field-based formulation in [26], the crack irreversibility condition (ḋ ≤ 0) is applied to the model at time t on each point x: the strong form of the initial boundary value problem can be obtained based on the minimization principles by taking the variation of the functional of Eq. (45), with the following boundary conditions on the displacement field and the phase field:

Phase field model in a ductile failure context
Ductile failure occurs in conjunction with extensive plastic deformation. The above phase field approach is applied herein to the plasticity model based on the von Mises hardening criterion. To achieve a thermodynamically-consistent model, the inclusion of the inelastic deformations into the total potential should be considered. This necessitates adding the plastic potential, ψ p , to the functional in Eq. (45), which yields the following relation: This introduces a modified elastoplastic history field: where β e and β p are the constants that are regulating the contribution rate of the elastic and plastic works respectively and the Macaulay bracket of an arbitrary variable x is defined as: The second term in Eq. (50) is the part of total energy dissipated by the plastic mechanisms, represented by the following rate form [42,46]: whereas the constitutive law expressed in Eq. (44) is adopted. The strong form of Eq. (47) 2 may be replaced by the following relation:

Integration algorithm
The integration procedure of this phase field modelling is constructed upon the elastic predictor-plastic corrector algorithm in an analogous way to the presented nonlocal damage model. The calculation of the stress field is carried out based on the von Mises rateindependent plasticity model that may be closely related to the isotropically-hardening plasticity model in [45]. The coupled plasticity model is implemented with UMAT and HETVAL subroutines in Abaqus/Standard (in the same manner that is explained earlier and in "Appendix" for the gradient damage model). In each increment, the phase field value is frozen and retrieved in routine as a temperature field. Given the flux value based on the updated history field from the previous increment in the UMAT code, the updated phase field value is computed through the HETVAL subroutine. This staggered-type solution procedure bypasses the need to calculate the coupled phase-field-displacement tangent terms as is presented in the UEL monotonic scheme in [47].
Having ε and the phase field variable from the previous increment (d n ), the predicted elastic trial states are given by: The yield condition to check the plastic admissiblity is written as: If the yield condition is not satisfied, the following return mapping equation is solved iteratively for γ : By having the solution at hand, the update procedure followed as: The implementation of this algorithm is performed by the same strategy in dealing with the nonlocal gradient damage model discussed before, using steady state heat problem analogy referring to Eq. (34) and assuming the following associations:

Results and discussion
The capability of the presented diffusive strategy is assessed herein, through the analysis of two common benchmarks and by observation of force-deflection diagrams, damage evolution, and crack propagation prediction.

Notched specimen tensile test
First benchmark is a notched specimen under tension, with the geometry and boundary conditions depicted in Fig. 3. This example aims to observe the propagation of the damage and assess the nonlocal gradient damage framework robustness and applicability of the proposed methodology in the post-peak regime to circumvent mesh sensitivity. Due to the symmetry, only one quarter of the specimen is modelled and discretized, assuming distinct mesh discretizations (identically adopted from [48]) and using 4-noded plane strain quadrilateral elements. The material properties used for analysis are represented in Table 2, considering that the same length scale is used and treated as a material parameter for all mesh topologies. Assuming the imposed displacement of u 1 = 0.2 mm, the propagation of the local and non-local damage fields are illustrated in Figs. 4 and 5 in a same time increment and for three mesh sizes. Expectedly the localization of the damage is at the center of the specimen and, as it proceeds, a 45 • localization band is revealed. The local damage contours show different values by varying the mesh density (denoted by h), whereas the non-local model shows less discrepancy between different meshes at the same chosen time increment. The more sparse damage propagation in non-local case is attributed to the lower values obtained with the non-local model, and it is in a close agreement with the non-local integral-type Lemaitre damage model results presented in [48] (therein stated by L-D).
This can be more plainly observed in the force-deflection diagrams, by comparing the curve discrepancies in Fig. 6a, b for the local and non-local damage frameworks respectively, where it is visible that mesh sensitivity is attenuated in the case of the non-local model. The results are compared with the results obtained in [48] for local and non-local integral type Lemaitre damage models, showing some discrepancy in the non-local case in which the gradient model predicts less softening than the integral type.
Attenuation of mesh dependency also is observable in Fig. 7a, b showing the evolution of the local and non-local damage fields at the center of the specimen and the effect of the non-local model in restoring mesh objectivity. It is also worth noting the smoother evolution of the non-local damage field compared to the local one for the more refined mesh (Fig. 8). More importantly, these diagrams prove the applicability of parallel treatment of  the non-local field with the heat equation solver and therefore the proposed methodology may be utilized with reliability, for this purpose, with this commercial software in the conditions described.

Single edged notched shear test
The proposed methodology for the phase field model is analyzed via a single edged notched test (SENT) specimen as can be found in [40,43], among other brittle fracture based articles. The geometric setup of the 2D specimen is adopted from [37] with L = 5 mm and a horizontal displacement-controlled loading applied on the top edge of the specimen, as is depicted in Fig. 9a. Notice that the notch is considered as a geometrical entity with its tip located at the center of the specimen where, based on experimental evidence, the crack begins to propagate. The aim here is to observe the ability of the proposed algorithm that  from the undamaged zones. The material properties, mainly, taken from reference [37], are summarized in Table 3 while the same contribution rates of elastic and plastic works are considered herein, β e = β p = 1.0.
The contour plots of equivalent plastic strain and phase field at various displacement loading states are depicted in Fig. 10. The crack initiation, as expected, occurs at the initial crack tip in the center of the specimen, and proceeds horizontally to the right edge of the specimen. The pattern of the "regularized crack" in Fig. 10b is obtained resorting to the embedded solution of the heat conduction equation in Abaqus solver package for thermomechanical problems. Performing a denser mesh discretization with an effective element size of 0.04 mm, the mesh sensitivity of the model is indicated via the force deflection diagram in Fig. 11 and it is in a close agreement with the result reported in [37]. Also notable that no delayed softening is visible due to the chosen null value of plastic work threshold.

Concluding remarks
In this study, the evolution of damage and failure, via non-local gradient damage and phase field models, respectively, was simulated resorting to an analogy to heat diffusion in solids, using the built-in thermo-mechanically coupled finite-element solution procedure in Abaqus. The developed procedure was justified and assessed in terms of qualitative verification of crack propagation patterns and its capability of avoiding mesh dependence pathologies.
The utilization of this approach may be viewed as a simple alternative to include the regularization procedures associated with both gradient and phase field models in commercial codes with no need of cumbersome implementations of explicit weak form derivations within required lengthy user files.