A Thermodynamically Consistent Phase Field Framework for Anisotropic Damage Propagation

Corresponding author https://doi.org/10.1590/1679-78255959 Abstract In the present work, a thermodynamically consistent damage phase field formulation is adapted to include the effect of preferential directions in the damage evolution. A scalar damage variable is associated with each predefined preferential direction of crack propagation. Any other direction is penalized by a parameter ( 1) β  that represents the degree of anisotropy of the fracture. When 0 β= , the isotropic case is recovered. When there is more than one preferential direction, the material is considered totally fractured when any of the damage variables reaches value one. Simulations using the developed model show that it can reproduce the expected crack propagation pattern for materials with one and two preferential directions. In particular, the model was successful in simulating a zigzag crack pattern commonly obtained in double cantilever beam of spinel 2 4 MgAl O crystals. The model is fully dynamic in the sense that describes the actual time evolution of the unknown variables, in particular of damage growth. Moreover, anisotropic mechanical response can be easily included in the model by


INTRODUCTION
Anisotropy is used to indicate variations in physical properties or certain material behavior along different directions. In particular, the term anisotropic damage has been used in literature to refer to different phenomena. Damage models with tension-compression split to prevent cracking under compressive loading, found in (Wu et al., 2020;Ambati et al., 2014;Levitas et al., 2018), are examples of anisotropic damage.
Another use is in models for damage-induced material anisotropy, i.e., when the damage introduces new symmetry planes in the initially isotropic material, making it anisotropic or changing its original anisotropy as in (Jarić et al., 2013;Brünig, 2004;Chow et al., 1987). In phase field models, this type of anisotropic damage has been modeled either with second order tensors or with multiple scalar damage variables which degrades the elastic part given in terms of the constitutive tensors as in (Fassin et al., 2019;Mozaffari and Voyiadjis, 2015). This idea follows the hypothesis that the damage comes from planar microvoids which are dependent on the stress or strain directions. Therefore, cracks may create symmetry planes in the material making it anisotropic. It is also known that crack evolution depends, in the microstructure level, on the material symmetries and grain orientation. Orientation of the crystalline planes can define preferential cleavage planes as observed in crystalline materials or inside the grains in polycrystals. This phenomenon is also referred in the literature by anisotropic crack propagation or directional fracture (Clayton and Knap, 2015). This directional effect on damage growth can also be motivated purely by the macrostructure, as in honeycomb materials considered in Réthoré et al. (2017). This kind of material has isotropic bulk mechanical behavior but has two preferential admissible crack angles. Its structure has different geometric strength and does not allow the crack to grow in any direction. Other materials in which the macrostructure has influence on the directional damage propagation are fiber reinforced unidirectional composites. Their distinct interfaces are weak planes with favorable local stress conditions (Talreja, 2014), which can lead to damage dependency on the lamina orientation; in such case, they also have anisotropic mechanical response.
Within the context of phase field models, damage anisotropy caused by preferential cleavage directions has been considered by penalization parameters (anisotropy parameters) in the crack surface energy associated to predefined set of preferential directions for the crack growth (Nguyen et al., 2017;Clayton and Knap, 2015). This approach has been applied for example to composites (Denli et al., 2020;Dal et al., 2017), polycrystals (Clayton and Knap, 2015), organic tissues Teichtmeister and Miehe, 2015) and thin films (Li et al., 2015). Beyond the ability to capture the crack path orientation dependency, another advantage of this method is that, by the use of multiple phase-field variables, damage associated with different directions or with distinct damage mechanisms can be included. For composites materials, it is possible to associate different fracture energies to distinct phase-field variables to model fiber and matrix fracture and a third damage mechanism specifically for shear, as done by Bleyer and Alessi (2018). The works derive the isothermal governing equations using variational methods, but the way the damage variables couple with the equilibrium equation does not make clear the thermodynamically consistency.
The aim of the present work is to obtain a rather general thermodynamically consistent isothermal phase field model with the directional effect associated to preferential orientations given in the driving fracture energy. Tensioncompression asymmetry is included in order to prevent damage growth and crack healing for compression loads.
In addition, the model is fully dynamic and describes the actual time evolution of the unknown variables, in particular the damage growth. This differs from several other models presented in the literature that assume quasi-static damage growth (sometimes a pseudo-time is introduced to compute crack growth). This dynamic feature of predicting damage growth at earlier stages may help devise actions to prevent the final failure of structures.
Our model is developed by extending the framework presented in Boldrini et al. (2016) to include the directional fracture effect. In Boldrini et al. (2016), a general thermodynamically consistent non-isothermal continuum thermomechanic phase field framework for the evolution of damage in materials under the hypothesis of small deformation was developed. This framework is based on the use of the principle of virtual power (PVP), the balance of energy and the Clausius-Duhem inequality for the entropy. It is rather general and allows the use of different free-energy potentials to describe different material behaviors.
For simplicity of exposition, we focus here in the case of materials with isotropic and isothermal response with preferential damage directions, such as for honeycomb materials. Anisotropic material response is not considered here, but can be easily implemented just by modifying the elasticity tensor. Similarly, non-isothermal situations could be considered in our model and for details see (Boldrini et al., 2016).
The paper is organized as follows. In Section 2, the governing equations and the finite element discretization are presented followed by the algorithm used to solve the problem. In Section 3, the numerical results are considered. The conclusions are addressed in Section 4.

Governing equations
The phase field methodology was originally developed to solve separation of fluids (Cahn, 1959). Later, it was applied to damage processes in materials using a smeared representation of the crack geometry by a scalar variable ϕ depending on the position x of the body Ω and meaning the volumetric fraction of damaged material. For virgin material, 0 ϕ = ; for fractured material 1 ϕ = , and 0 1 ϕ < < is associated to damaged material between these two extreme states.
Although our approach is rather general, the model presented here focuses on the anisotropic damage propagation due to geometrical effects from the material structure, such as preferential cleavage planes and fibers in composites. It can also be applied to honeycomb materials, in which failure occurs and propagates along the corners of cells, following an angle that depends on the material orientation .
The crack orientation dependence is introduced in the diffusion term of the damage surface energy. Preferential directions i m ( ) 1 , 2,…, N i = for the crack growth are defined and a scalar damage variable i ϕ is associated to each direction. Consequently, instead of having only one equation for the damage evolution, there will be N equations for scalar damage phase fields, one for each preferential direction. As in the isotropic case developed in Boldrini et al. (2016), these equations are obtained by applying the principle of virtual power and the entropy inequality based on suitable specific free energy potential Ψ . By making use of the Coleman-Noll arguments, the relevant constitutive relations for the model are obtained and the thermodynamically consistency is guaranteed.
The first physical law to be considered is the conservation of mass. For a closed system defined by an arbitrary region t D , with boundary surface t D ∂ , moving with the body that occupies a domain 3 Ω R ⊂ of mass m continuously distributed, the conservation of mass is characterized by the following statement: where ρ is the density scalar field and v represents the velocity vector field.
The virtual power of the interior forces int P gives the internal response of any region t D at time t . It is given here by the sum of the classical stress power and the contribution of generalized loads associated to microscopic motions related to damage and expressed in the spatial configuration as follows: where T is the Cauchy stress tensor field, D is the gradient of the virtual velocity vector field, b is the volume density of energy exchanged by variation of a unit of ϕ in a unit of time (represented by ĉ ), and h is the flux of energy associated to the spatial variation of a unit of ϕ in a unit of time (represented by ĉ ∇ ), see (Frémond, 2013). Integrating by parts the first and third integrals of the previous expression, we can rewrite the virtual power of the interior loads as where n is the normal vector field on the domain boundary.
The virtual power of external actions for any region t D at time t is given by The first two integrals in the previous equation comes from the classical formulation, where t is the macroscopic stress vector field and f is the body force vector field per unit of mass. The last two integrals are associated to the virtual powers of actions at the distance and the contact force associated to the damage. Variable h t is the superficial density of energy supplied to the material by the flux h and a is the specific density of energy supplied to the material from the exterior acting on the microscopic level and affecting the damage level of the material (for instance, energy supplied by irradiation, electrical or chemical actions which modifies the microscopic bonds).
The PVP states that the power of external loads must be equal to the power of internal loads plus the change of the kinetic energy K for any virtual field v and ˆi c . Therefore, Substituting Equations 3, 4 and 6 in 5 and rearranging the terms, we obtain As the virtual actions � and ̂ and the domain are arbitrary, the terms in parentheses of the previous equations must be simultaneously zero resulting in the following equilibrium equations and boundary conditions: div 0 From now on, we will only use the letter " k " to indicate a preferential direction. The principles of conservation of mass and linear momentum balance must be supplemented by the first and second principles of thermodynamics. The first law of thermodynamics establishes the balance of mechanical and thermal energy as The thermal power ( ) Q t is given in terms of the specific heat source density scalar field ( r ) and the heat flux vector and the internal energy E is given by where e is the specific internal energy density.
Substituting Equations 3, 12 and 13 in Equation 11, we achieve the local form of the first law of thermodynamics The second principle, or entropy inequality principle, postulates that the total entropy production per unit of time for all processes can never be negative. It is expressed in terms of the Clausius-Duhem differential form, based on the arguments presented by Fabrizio, Giorgi and Morro (2006) and Boldrini et al. (2016), as div r where η is the specific entropy density and θ + q k and r θ are, respectively, the entropy flux and the specific entropy production term and θ is the absolute temperature.
In order to obtain thermodynamically consistent constitutive relations, we assume that the previous relation is written in terms of the specific Helmholtz free energy density functional Ψ . By expressing this functional in terms of the internal energy ( e ) and the entropy (η ) as -e θη Ψ = and using the balance of energy given in Equation 14, the previous expression may be rewritten as We assume now that the constitutive properties of the considered material are expressed in terms of a free-energy and φ and ∇φ represent the dependency on all damage scalar variables and their gradients and E the infinitesimal strain tensor given by the symmetric part s ∇ u of the displacement gradient.
As in Frémond (2013), we assume that the terms T and b can be decomposed in their non-dissipative (reversible) and dissipative (irreversible) parts indicated, respectively, by the superscripts ( ) r and ( ) ir and that q is totally irreversible ( From now on we make use of the Coleman-Noll arguments to obtain the constitutive relations. The non-dissipative terms must be such that for any admissible process they produce no entropy increase. Therefore, Since that the dependence of the previous equation on θ  , ∇v , are linear and that such quantities can assume arbitrary values different of zero depending on the adopted forcing terms f and r , their respective coefficients must be zero. Therefore, Moreover, we take the reversible parts of b , h and T respectively as Consequently, the irreversible part of 21 is The previous inequality is satisfied assuming Latin American Journal of Solids and Structures, 2021, 18(1), e314 7/17 where k λ  , b  and c  are non-negative coefficients.
By substituting Equations (28) and (32) in (17), we obtain the stress tensor together with the following restriction: Similarly, from Equations 18, 26 and 31, we have By collecting all the previous results, we finally obtain the governing equation of the model as div 0 ρ ρ Considering a homogeneous material under small strain with material density  By choosing suitable free-energy functionals Ψ in the previous system, it is possible to obtain thermodynamically consistent damage models for different types of materials. To include the possibility of damage, we take the expression of the free-energy density as the sum of the elastic energy density of the cracked body, ( ) , φ E  , and the energy density required to create the crack, ( ) , , ∇ φ φ  ω , according to the Griffith criterion, as The elastic energy density of the cracked body is given by : where denotes the degraded elastic tensor of the isotropic material  is the elasticity tensor of the virgin material; 0 k is the bulk modulus for the undamaged material, and ( ) The second term in Equation 47 was included to assure that there is no degradation under compression but allowing damage evolution under traction and shear. I is the identity tensor and ( ) g φ is the degradation function, which is commonly taken as a quadratic function (Miehe et al., 2010;Nguyen et al., 2017). Therefore, where 1 κ  is a small parameter used to maintain the well-posedness of the system in case any of the damage variables reaches value one. The second term of the free-energy density (Equation 45) is based on the concept of crack density function per unit of volume γ and the Griffith fracture energy c g as follows: To explain this last expression, let us firstly recall that, for the isotropic case, the crack density function γ is defined as where l is a positive parameter related to the width of the damage phase field layers. As in Nguyen et al. (2017), this formulation was extended to the anisotropic case, including the influence of the direction, according to the following expressions: where ω is the second order structural tensor defined by the unit vector m normal to the preferential cleavage plane with respect to the material coordinates (Clayton and Knap, 2015). It is important to notice that the degraded elastic tensor does not become anisotropic, but makes the energy release rate orientation dependent. In fact, β is the parameter used to control how strong is this dependency on one specific direction and, consequently, the tendency of the crack to follow that direction. High values of ( 1) β >> avoid the crack to propagate on planes not normal to m . To recover the isotropic case, we take only one damage variable and 0 β = .
By substituting the defined free-energy potential in Equation 44 and proceeding as in (Miehe et al., 2010;Nguyen et al., 2017) to ensure the non-reversibility condition of the damage, the strong formulation becomes 1 div : sym div div 1 : for where k  is the strain history functional defined as (Miehe et al., 2010;Nguyen et al., 2017) (

Finite Element Approximation
The first step of the finite element approximation of the given model is to obtain the weak forms of Equations 53a and 53b. For this purpose, they are multiplied by suitable vector 1 w and scalar 2 w test functions and integrated over the domain Ω . After integrating by parts and applying the boundary conditions, we obtain the following weak forms: Two-dimensional spatial finite element discretization is considered for each of the previous equations. The conforming global approximations for the displacement, velocity, acceleration, damage and test functions for a mesh of K elements are obtained from the standard assembling procedure of the element approximations, respectively, by The global derivatives of the element approximations for the displacement and velocity, in terms of infinitesimal strain tensor and the rate of deformation tensor, and the damage gradient are given, respectively, by Analogously for the test functions, where u B and φ B are matrices with the global derivatives of the shape functions given for 2D problems and each node i, respectively, by The expressions for the shape functions of the triangular and quadrangular elements used in the present work can be found in (Bittencourt, 2014).
A semi-implicit time discretization scheme is adopted where each equation is separately solved by an implicit time integration method. Firstly, the damage evolution equation for each preferential direction k is solved using the backward Euler method. After that, the obtained damage solutions . For more details see (Chiarelli et al., 2017). The time-marching rule for the damage equation related to a preferential direction k at time step 1 n + is given by The global operators of the previous equation are obtained by the standard assembling procedure of the following element operators: The approximate solution of the displacement at time step 1 n + , The previous global operators are also obtained by the standard assembling procedure of the following element operators: After solving the evolution equation for the displacement (Equation 65), the acceleration and velocity vectors are updated by using the following expressions: The overall procedure is given in Algorithm

Example1: I-shaped-specimen tensile test
The simulation of an I-shaped tensile test specimen was performed to verify the capability of the model to obtain anisotropic damage. The bottom edge of the specimen is fixed and a vertical prescribed displacement of  that the crack growth behaved anisotropically as the damage increased predominantly at o 45 .

Example2: Notched tensile test
The second example is the single notched tensile test of a 1 m 1 m × plate with a horizontal notch in the middle of the left edge with half width length. The plate is constrained at its bottom edge while a prescribed displacement of  For stronger damage anisotropy, more time steps were required for the complete crack propagation. To obtain a smoother crack path, it is necessary to increase the number of elements in the region of possible crack propagation.

Example 3: Zig-zag crack propagation
The influence of orientation of preferred cleavage planes in the crack propagation can be observed in a double cantilever beam with a guiding groove. In this case, a distinct zig-zag crack propagation pattern appears. This same pattern happens in some single crystals, such as 2 4 MgAl O , as shown in Wu et al. (1995).
The present example aims to reproduce this crack behavior by simulating the boundary conditions of the experiment. A square domain with 2 mm side was submitted to the boundary conditions illustrated in Figure 4 with incrementally applied displacement of . We recall that in case of more than one preferential direction, the material is considered totally fractured when one of the damage variables reaches one. The results are shown in Figure 5.  ϕ reaches the bottom of the central band, where the crack kinks, following another energetically preferred direction associated with the second phase field 2 ϕ ( Figure 5 (c-f)). Another crack kinking occurs near the end of the plate before the second damage variable propagates until the domain boundary, producing a crack pattern similar to the observed in experimental tests in Wu et al. (1995). The results obtained by using the time steps

CONCLUSION
This paper presented a thermodynamically consistent phase field model that takes into account the directionality effect resulting from the existence of cleavage planes or geometric resistance in the material macroscale. This effect was considered by introducing anisotropic parameters in the crack density function. The dynamic equations for the displacement and damage are solved by a spatial finite element approximation and a semi-implicit time scheme.
The numerical examples showed that the model was able to recover the isotropic damage evolution by using the anisotropic parameter equal to zero. It was also able to reproduce the expected crack propagation pattern for a material with initial damage at o 0 and preferential cleavage plane oriented at o -45 . The influence of parameter β , that regulates how strong is the fracture anisotropy, was studied for a notched tensile test. For larger β , the crack angle tends to the given preferential direction. Simulations with different meshes also revealed that, except for the obvious question of accuracy, crack directions were basically independent of mesh refinement.
Finally, by choosing appropriate boundary conditions and time step ( -3 1 10 s ≤ × ), the model was also able to reproduce the zig-zag crack pattern with two preferential directions. For time step -2 1 10 s × , the expected zig-zag pattern was smeared, probably due to excessive numerical diffusion.
From the simulations results, it is possible to conclude that the presented approach was successful in modeling directional crack propagation. In order to evaluate separately the effect of the anisotropy due to the directional effect, we only presented the case of materials with isotropic mechanical response but with preferential damage directions as, for instance, honeycomb materials. However other types of anisotropies, such as anisotropic mechanical response could be easily included in the model just by modifying the elasticity tensor. For simplicity of exposition, we just presented isothermal numerical simulations; again, non-isothermal situations could also be easily implemented in our framework. In fact, the thermodynamic consistency of our approach encourages future analyses considering effects of temperature and correlation with experiments. Finally, the important fully dynamic feature of our model makes possible to describe the actual time evolution of the unknown variables and in particular the damage growth. By predicting damage growth at earlier stages, devise actions may be taken to prevent the final failure of structures.