A Thermodynamically Consistent Phase Field Framework for Anisotropic Damage Propagation

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 $(\beta \gg 1)$ that represents the degree of anisotropy of the fracture. When $\beta =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 $MgA{l}_2{O}_4$ 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 modifying the elasticity tensor.

Keywords:
Phase field models; Damage; Anisotropy; Finite element method

Graphical Abstract

1 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 (Gültekin et al., 2017; 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. Tension-compression 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 thermo-mechanic 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.

2 PHASE FIELD MODEL FOR ANISOTROPIC CRACK PROPAGATION

2.1 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 $\phi$ depending on the position $x$ of the body $\Omega$ and meaning the volumetric fraction of damaged material. For virgin material, $\phi =0$; for fractured material $\phi =1$, and $0<\phi <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 (Réthoré et al., 2017).

The crack orientation dependence is introduced in the diffusion term of the damage surface energy. Preferential directions $m_i$ $\left(i=\mathrm{1,}\mathrm{2,}\dots ,N\right)$ for the crack growth are defined and a scalar damage variable ${\phi }_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 $\Psi$. 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 $D_t$, with boundary surface $\partial D_t$, moving with the body that occupies a domain $\Omega \subset R^3$ of mass $m$ continuously distributed, the conservation of mass is characterized by the following statement:

$\dot{\rho }+\rho di{v}_{x}v=\mathrm{0,}$ (1)

where $\rho$ is the density scalar field and $v$ represents the velocity vector field.

The virtual power of the interior forces $P_{\mathrm{int}}$ gives the internal response of any region $D_t$ 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:

$P_{\mathrm{int}}=\underset{D_t}{\overset{}{{\displaystyle \int }}}T:\widehat{D}d{D}_t+{\displaystyle \sum }_{i=1}^N\left(\underset{D_t}{\overset{}{{\displaystyle \int }}}\left(b_i{\widehat{c}}_i+h_i\cdot \nabla {\widehat{c}}_i\right)d{D}_t\right),$ (2)

where $T$ is the Cauchy stress tensor field, $\widehat{D}$ is the gradient of the virtual velocity vector field, $b$ is the volume density of energy exchanged by variation of a unit of $\phi$ in a unit of time (represented by $\widehat{c}$), and $h$is the flux of energy associated to the spatial variation of a unit of $\phi$ in a unit of time (represented by $\nabla \widehat{c}$), 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

$P_{\mathrm{int}}=\underset{\partial D_t}{\overset{}{{\displaystyle \int }}}Tn\cdot \widehat{v}d\partial D_t-\underset{D_t}{\overset{}{{\displaystyle \int }}}\left(divT\right)\cdot \widehat{v}d{D}_t+{\displaystyle \sum }_{i=1}^N\left(\underset{D_t}{\overset{}{{\displaystyle \int }}}{b}_i{\widehat{c}}_{i}d{D}_t+\underset{\partial D_t}{\overset{}{{\displaystyle \int }}}\left(h_i\cdot n\right){\widehat{c}}_{i}d\partial D_t-\underset{D_t}{\overset{}{{\displaystyle \int }}}\left(div{h}_i\right){\widehat{c}}_{i}d{D}_t\right),$ (3)

where $n$ is the normal vector field on the domain boundary.

The virtual power of external actions for any region$D_t$ at time $t$ is given by

$P_{ext}=\underset{\partial D_t}{\overset{}{{\displaystyle \int }}}t\cdot \widehat{v}d\partial D_t+\underset{D_t}{\overset{}{{\displaystyle \int }}}\rho f\cdot \widehat{v}d{D}_t+{\displaystyle \sum }_{i=1}^N\left(\underset{\partial D_t}{\overset{}{{\displaystyle \int }}}{t}_{h_i}{\widehat{c}}_{i}d\partial D_t+\underset{D_t}{\overset{}{{\displaystyle \int }}}\rho a{\widehat{c}}_{i}d{D}_t\right).$ (4)

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 $t_h$ 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 $\widehat{v}$ and ${\widehat{c}}_i$. Therefore,

$P_{ext}=\frac{DK\left(t\right)}{Dt}+P_{\mathrm{int}},$ (5)

where

$\frac{DK\left(t\right)}{Dt}=\underset{D_t}{\overset{}{{\displaystyle \int }}}\rho \dot{v}\cdot \widehat{v}d{D}_t.$ (6)

Substituting Equations 3, 4 and 6 in 5 and rearranging the terms, we obtain

$\begin{array}{l}\underset{\partial D_t}{\overset{}{{\displaystyle \int }}}\left(t-Tn\right)\cdot \widehat{v}d\partial D_t+\underset{D_t}{\overset{}{{\displaystyle \int }}}\left(\rho f+divT-\rho \dot{v}\right)\cdot \widehat{v}d{D}_t\\ +{\displaystyle \sum }_{i=1}^N\left(\underset{\partial D_t}{\overset{}{{\displaystyle \int }}}(t_{h_i}-h_i\cdot n){\widehat{c}}_{i}d\partial D_t+\underset{D_t}{\overset{}{{\displaystyle \int }}}\left(\rho a-b_i+div{h}_i\right){\widehat{c}}_{i}d{D}_t\right)=0.\end{array}$

As the virtual actions $\widehat{v}$ and $\widehat{c}$ and the domain $D_t$are arbitrary, the terms in parentheses of the previous equations must be simultaneously zero resulting in the following equilibrium equations and boundary conditions:

$\rho \dot{v}=div\left(T\right)+\rho f,in\Omega ,$ (7)

$div{h}_k+\rho a_k-b_k=0,withk=\mathrm{1,2,}\dots ,N,$ (8)

$Tn=t,in\partial \Omega ,$ (9)

$h_k\cdot n=t_{h_k},in\partial \Omega ,$ (10)

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

$P_{\mathrm{int}}+Q\left(t\right)=\frac{DE\left(t\right)}{Dt}.$ (11)

The thermal power $Q\left(t\right)$ is given in terms of the specific heat source density scalar field ($r$) and the heat flux vector field $(q)$ as

$Q=\underset{D_t}{\overset{}{{\displaystyle \int }}}\rho rd{D}_t-\underset{\partial D_t}{\overset{}{{\displaystyle \int }}}q\cdot \widehat{n}d\partial D_t$ (12)

and the internal energy $E$ is given by

$E=\underset{D_t}{\overset{}{{\displaystyle \int }}}ed{D}_t,$ (13)

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

$\rho \dot{e}=T:D+{\displaystyle \sum }_{k=1}^N\left(h_k\cdot \nabla {\dot{\phi }}_k+b_k{\dot{\phi }}_k\right)+\rho r-divq.$ (14)

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

$\rho \dot{\eta }-\rho \frac{r}{\theta }\ge -div\left(\frac{q}{\theta }+k\right),$ (15)

where $\eta$ is the specific entropy density and $\frac{q}{\theta }+k$ and $\frac{r}{\theta }$ are, respectively, the entropy flux and the specific entropy production term and $\theta$ 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 $\Psi$. By expressing this functional in terms of the internal energy ($e$) and the entropy ($\eta$) as $\Psi =e–\theta \eta$ and using the balance of energy given in Equation 14, the previous expression may be rewritten as

$-\rho \left(\dot{\Psi }+\dot{\theta }\eta \right)+T:D+{\displaystyle \sum }_{k=1}^N\left(h_k\cdot \nabla {\dot{\phi }}_k+b_k{\dot{\phi }}_k\right)+\theta divk-q\cdot \frac{1}{\theta }\nabla \theta \ge 0.$ (16)

We assume now that the constitutive properties of the considered material are expressed in terms of a free-energy $\Psi \left(\Gamma \right)$, where $\Gamma =\left(\rho ,\theta ,\phi ,\nabla \rho ,\nabla \theta ,\nabla \phi ,E\right)$ and $\phi$ and $\nabla \phi$ represent the dependency on all damage scalar variables and their gradients and $E$ the infinitesimal strain tensor given by the symmetric part ${\nabla }^{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 $\left(r\right)$ and $\left(ir\right)$ and that $q$ is totally irreversible ($q^{\left(r\right)}=0$) and $h_k^{\left(ir\right)}=0$, we have

$T=T^{\left(r\right)}+T^{\left(ir\right)},$ (17)

$b_k^{}=b_k^{\left(r\right)}+b_k^{\left(ir\right)},$ (18)

$h_k^{}=h_k^{\left(r\right)},$ (19)

$q=q^{\left(ir\right)}.$ (20)

By using the chain rule and representing the derivatives of $\Psi$ with respect to the variables given in $\Gamma$, as for example $\frac{\partial \Psi }{\partial \theta }={\Psi }_{\theta }$, and considering small strain, we obtain

$\begin{array}{r}-\rho \left({\Psi }_{\theta }+\eta \right)\dot{\theta }+{\displaystyle \sum }_{k=1}^N\left[\left(-\rho {\Psi }_{{\phi }_k}+b_k{}^{\left(r\right)}+b_k{}^{\left(ir\right)}\right){\dot{\phi }}_k\right]+{\rho }^2{\Psi }_{\nabla \rho }\nabla \left(div\left(v\right)\right)+\rho {\Psi }_{\nabla \rho }\left(div\left(v\right)I+{\left(\nabla v\right)}^T\right)\nabla \rho \\ -\rho {\Psi }_{\nabla \theta }\stackrel{·}{\overline{\nabla \theta }}+{\displaystyle \sum _{k=1}^N\left[\left(-\rho {\Psi }_{\nabla {\phi }_k}+h_k^{\left(r\right)}\right)\cdot \stackrel{·}{\overline{\nabla {\phi }_k}}\right]-q^{\left(ir\right)}\cdot \frac{1}{\theta }}\nabla \theta +\theta divk\\ +\left(T^{\left(r\right)}+T^{\left(ir\right)}-\rho {\Psi }_E+{\rho }^2{\Psi }_{\rho }I+{\displaystyle \sum _{k=1}^N\nabla {\phi }_k}\otimes h_k^{\left(r\right)}\right):\nabla v+\rho {\Psi }_E{\left(\nabla v\right)}^T{\nabla }^{s}u\ge 0.\end{array}$ (21)

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,

$\begin{array}{r}-\rho \left({\Psi }_{\theta }+\eta \right)\dot{\theta }+{\displaystyle \sum }_{k=1}^N\left[\left(-\rho {\Psi }_{{\phi }_k}+b_k{}^{\left(r\right)}\right){\dot{\phi }}_k\right]+{\rho }^2{\Psi }_{\nabla \rho }\nabla \left(div\left(v\right)\right)-\rho {\Psi }_{\nabla \theta }\stackrel{·}{\overline{\nabla \theta }}+{\displaystyle \sum _{k=1}^N\left[\left(-\rho {\Psi }_{\nabla {\phi }_k}+h_k^{\left(r\right)}\right)\cdot \stackrel{·}{\overline{\nabla {\phi }_k}}\right]}\\ +\left(T^{\left(r\right)}+T^{\left(ir\right)}-\rho {\Psi }_E+{\rho }^2{\Psi }_{\rho }I+{\displaystyle \sum _{k=1}^N\nabla {\phi }_k}\otimes h_k^{\left(r\right)}\right):\nabla v=0.\end{array}$ (22)

Since that the dependence of the previous equation on $\dot{\theta }$,$\nabla v$, $\nabla \left(div\left(v\right)\right)$ and $\stackrel{·}{\overline{\nabla \theta }}$ 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,

${\Psi }_{\nabla \rho }=\mathrm{0,}$ (23)

${\Psi }_{\nabla \theta }=\mathrm{0,}$ (24)

$\eta =-{\Psi }_{\theta }.$ (25)

Moreover, we take the reversible parts of $b$, $h$ and $T$ respectively as

$b_k^{\left(r\right)}=\rho {\Psi }_{{\phi }_k},$ (26)

$h_k^{\left(r\right)}=\rho {\Psi }_{\nabla {\phi }_k},$ (27)

$T^{\left(r\right)}=\rho {\Psi }_E-{\rho }^2{\Psi }_{\rho }I-{\displaystyle \sum }_{k=1}^N\nabla {\phi }_k\otimes h_k^r.$ (28)

Consequently, the irreversible part of 21 is

${\displaystyle \sum }_{k=1}^N{b}_k^{\left(ir\right)}{\dot{\phi }}_k+T^{\left(ir\right)}:\nabla v-q^{\left(ir\right)}\cdot \frac{1}{\theta }\nabla \theta +\theta divk\ge 0.$ (29)

For $k=0$,

${\displaystyle \sum }_{k=1}^N{b}_k^{\left(ir\right)}{\dot{\phi }}_k+T^{\left(ir\right)}:\nabla v-q^{\left(ir\right)}\cdot \frac{1}{\theta }\nabla \theta \ge 0.$ (30)

The previous inequality is satisfied assuming

$b_k^{\left(ir\right)}={\tilde{\lambda }}_k{\dot{\phi }}_k,$ (31)

$T^{\left(ir\right)}=\tilde{b}D,$ (32)

$q^{\left(ir\right)}=-\tilde{c}\frac{1}{\theta }\nabla \theta ,$ (33)

where ${\tilde{\lambda }}_k$, $\tilde{b}$ and $\tilde{c}$ are non-negative coefficients.

By substituting Equations (28) and (32) in (17), we obtain the stress tensor

$T=T^{\left(r\right)}+T^{\left(ir\right)}=\rho {\Psi }_E-{\rho }^2{\Psi }_{\rho }I-\rho sym\left({\displaystyle \sum }_{k=1}^N\nabla {\phi }_k\otimes {\Psi }_{\nabla {\phi }_k}\right)+\tilde{b}D,$ (34)

together with the following restriction:

$skw\left({\displaystyle \sum }_{k=1}^N\nabla {\phi }_k\otimes h_k^r\right)=0.$ (35)

Similarly, from Equations 18, 26 and 31, we have

$b_k=b_k^{\left(r\right)}+b_k^{\left(ir\right)}=\rho {\Psi }_{{\phi }_k}+{\tilde{\lambda }}_k{\dot{\phi }}_k.$ (36)

By collecting all the previous results, we finally obtain the governing equation of the model as

$\dot{\rho }+\rho di{v}_x\dot{u}=\mathrm{0,}$ (37)

$\dot{u}=v,$ (38)

$\rho \dot{v}=divT+\rho f,in\Omega ,$ (39)

${\tilde{\lambda }}_k{\dot{\phi }}_k=div\left({\rho }_0{\Psi }_{\nabla {\phi }_k}\right)-\rho {\Psi }_{{\phi }_k}+{\rho }_0{a}_k,fork=\mathrm{1,}\dots ,N,$ (40)

$\rho \dot{e}-\rho r=T:D+{\displaystyle \sum }_{k=1}^N\left[\rho {\Psi }_{\nabla {\phi }_k}\cdot \nabla {\dot{\phi }}_k+{\rho }_0{\Psi }_{{\phi }_k}{\dot{\phi }}_k+{\tilde{\lambda }}_k{\dot{\phi }}_k{\dot{\phi }}_k\right]-div\left(\tilde{c}\frac{1}{\theta }\nabla \theta \right),$ (41)

$T=\rho {\Psi }_E-{\rho }^2{\Psi }_{\rho }I-\rho sym\left({\displaystyle \sum }_{k=1}^N\nabla {\phi }_k\otimes {\Psi }_{\nabla {\phi }_k}\right)+\tilde{b}D,$ (42)

$e=\Psi -\theta {\Psi }_{\theta }.$ (43)

Considering a homogeneous material under small strain with material density $\rho ={\rho }_0$assumed constant, we have the following system of governing equations for the isothermal case:

$\{\begin{array}{l}{\dot{\rho }}_0=0in\Omega \\ {\rho }_0\ddot{u}=divT+{\rho }_{0}f\\ {\tilde{\lambda }}_k{\dot{\phi }}_k=div\left({\rho }_0\frac{\partial \Psi }{\partial \nabla {\phi }_k}\right)-{\rho }_0\frac{\partial \Psi }{\partial {\phi }_k}fork=\mathrm{1,}\mathrm{2,}\dots ,N\\ T={\rho }_0{\Psi }_E-{\rho }_{0}sym\left({\displaystyle \sum }_{k=1}^N\nabla {\phi }_k\otimes {\Psi }_{\nabla {\phi }_k}\right)+\tilde{b}D\end{array}$ (44)

where $E={\nabla }^{s}u=\frac{1}{2}\left(\nabla u+{\left(\nabla u\right)}^T\right)$ and $D={\nabla }^{s}v=\frac{1}{2}\left(\nabla v+{\left(\nabla v\right)}^T\right)$. The coefficient $\tilde{b}$is associated to the viscous damping of the material, while ${\tilde{\lambda }}_k$ is associated to the rate of damage change.

By choosing suitable free-energy functionals $\Psi$ 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, ${\rm E}\left(\phi ,E\right)$, and the energy density required to create the crack, $\mathcal{J}\left(\phi ,\nabla \phi ,\omega \right)$, according to the Griffith criterion, as

${\rho }_0\Psi \left(E,\phi ,\nabla \phi ,\omega \right)={\rm E}\left(\phi ,E\right)+\mathcal{J}\left(\phi ,\nabla \phi ,\omega \right).$ (45)

The elastic energy density of the cracked body is given by

${\rm E}\left(\phi ,E\right)=\frac{1}{2}\left[E:ℂ\left(\phi \right):E\right],$ (46)

where $ℂ$denotes the degraded elastic tensor of the isotropic material

$ℂ\left(\phi \right)=g\left(\phi \right){ℂ}^0+k_{0}I\otimes I\left[1-g\left(\phi \right)\right]sig{n}^{-}tr\left(E\right).$ (47)

${ℂ}^0$is the elasticity tensor of the virgin material; $k_0$is the bulk modulus for the undamaged material, and $sig{n}^{-}\left(x\right)=1$ if $x<0$ and $sig{n}^{-}\left(x\right)=0$ if $x\ge 0$.

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\left(\phi \right)$is the degradation function, which is commonly taken as a quadratic function (Miehe et al., 2010; Nguyen et al., 2017). Therefore,

$g\left(\phi \right)=\left(1-\kappa \right){\displaystyle \prod }_k^N{\left(1-{\phi }_k\right)}^2+\kappa$ (48)

where $\kappa \ll 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 $\gamma$ and the Griffith fracture energy $g_c$ as follows:

$\mathcal{J}\left(\phi ,\nabla \phi ,\omega \right)=g_c{\displaystyle \sum }_{k=1}^N{\gamma }_k\left({\phi }_k,\nabla {\phi }_k\right).$ (49)

To explain this last expression, let us firstly recall that, for the isotropic case, the crack density function $\gamma$ is defined as

$\gamma \left(\phi ,\nabla \phi \right)=\frac{1}{2l}{\phi }^2+\frac{l}{2}\nabla \phi \cdot \nabla \phi$ (50)

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:

${\gamma }_k\left({\phi }_k,\nabla {\phi }_k\right)=\frac{1}{2l}{\phi }_k^2+\frac{l}{2}{\omega }_k:\left(\nabla {\phi }_k\otimes \nabla {\phi }_k\right),$ (51)

${\omega }_k=I+\beta \left(I-m_k\otimes m_k\right),$ (52)

where $\omega$ 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, $\beta$ 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 $(\beta >>1)$ avoid the crack to propagate on planes not normal to $m$. To recover the isotropic case, we take only one damage variable and $\beta =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

$\{\begin{array}{l}\ddot{u}=\frac{1}{{\rho }_0}div\left(ℂ:E\right)-\frac{l{g}_c}{{\rho }_0}sym\left({\displaystyle \sum }_{k=1}^{N}div\left(\left({\omega }_k\nabla {\phi }_k\right)\otimes \nabla {\phi }_k\right)\right)+\frac{\tilde{b}}{{\rho }_0}divD+f\\ {\dot{\phi }}_k=\frac{g_{c}l}{{\tilde{\lambda }}_k}{\omega }_k:{\nabla }^2{\phi }_k+\frac{1}{{\tilde{\lambda }}_k}{\mathcal{H}}_k-\frac{g_c}{{\tilde{\lambda }}_{k}l}{\phi }_{k}fork=\mathrm{1,2,}\dots ,N\end{array}$ (53)

where ${\mathcal{H}}_k$ is the strain history functional defined as (Miehe et al., 2010; Nguyen et al., 2017)

${\mathcal{H}}_k\left(x,t\right)=\underset{\tau \in \left[\mathrm{0,}t\right]}{\max }\left\{\left(1-\kappa \right){\displaystyle \prod }_{i\ne k}{\left(1-{\phi }_i\right)}^{2}E\left(x,\tau \right):{ℂ}_h\left(x,\tau \right):E\left(x,\tau \right)\right\},$, (54)

with ${ℂ}_h\left(x,\tau \right)={ℂ}^0-k_{0}I\otimes Isig{n}^{-}tr\left(E\left(x,\tau \right)\right).$

2.2 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 $w_1$ and scalar $w_2$ test functions and integrated over the domain $\Omega$. After integrating by parts and applying the boundary conditions, we obtain the following weak forms:

$\begin{array}{r}{\displaystyle \underset{D_t}{\int }\ddot{u}}{w}_{1}d{D}_t=-{\displaystyle \underset{D_t}{\int }\frac{1}{{\rho }_0}ℂ:E:}\nabla w_{1}d{D}_t+\frac{l{g}_c}{{\rho }_0}{\displaystyle \underset{D_t}{\int }sym\left({\displaystyle \sum }_{k=1}^N\left(\left({\omega }_k\nabla {\phi }_k\right)\otimes \nabla {\phi }_k\right)\right)}:\nabla w_{1}d{D}_t-{\displaystyle \underset{D_t}{\int }\frac{\tilde{b}}{{\rho }_0}divD}:\nabla w_{1}d{D}_t\\ +{\displaystyle \underset{D_t}{\int }f\cdot w_{1}d{D}_t}+\frac{1}{{\rho }_0}{\displaystyle \underset{\partial D_t}{\int }t\cdot }{w}_{1}d\partial D_t\end{array}$ (55)

${\displaystyle \underset{D_t}{\int }{\dot{\phi }}_k{w}_{2}d{D}_t}=-{\displaystyle \underset{D_t}{\int }\frac{g_{c}l}{{\tilde{\lambda }}_k}({\omega }_k\nabla {\phi }_k)\cdot \nabla w_{2}d{D}_t}+{\displaystyle \underset{D_t}{\int }\frac{1}{{\tilde{\lambda }}_k}{\mathcal{H}}_k{w}_{2}d{D}_t}-{\displaystyle \underset{D_t}{\int }\frac{g_c}{{\tilde{\lambda }}_{k}l}{\phi }_k{w}_{2}d{D}_t.}$ (56)

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

$\begin{array}{l}u\left(x,t\right)\approx {\cup }_{e=1}^K{u}_h^e\left(x,t\right),\dot{u}\left(x,t\right)\approx {\cup }_{e=1}^K{\dot{u}}_h^e\left(x,t\right),\ddot{u}\left(x,t\right)\approx {\cup }_{e=1}^K{\ddot{u}}_h^e\left(x,t\right),\phi \left(x,t\right)\approx {\cup }_{e=1}^K{\phi }_h^e\left(x,t\right)\\ w_1\left(x,t\right)\approx {\cup }_{e=1}^K{w}_1{}_h^e\left(x,t\right)and{w}_2\left(x,t\right)\approx {\cup }_{e=1}^K{w}_2{}_h^e\left(x,t\right),\end{array}$

where $x=\left(x,y\right)$ are the global coordinates which are mapped to the standard reference system $\xi \times \eta$ by the mapping $\left(x,y\right)=\left(x\left(\xi ,\eta \right),y\left(\xi ,\eta \right)\right).$

The element approximation for each variable is given by the linear combinations of the local nodal basis functions, $N_1\left(x\right),N_2\left(x\right),\dots ,N_n\left(x\right)$, and the nodal time dependent values. Therefore, the element approximations for the displacement, velocity, acceleration, damage and test functions are given, respectively, by

$\begin{array}{l}{u}_h^e\left(x,t\right)=N_u\left(x\right){\widehat{u}}^e\left(t\right),{\dot{u}}_h^e\left(x,t\right)=N_u\left(x\right){\widehat{\dot{u}}}^e\left(t\right),{\ddot{u}}_h^e\left(x,t\right)=N_u\left(x\right){\widehat{\ddot{u}}}^e\left(t\right),{\phi }_h^e\left(x,t\right)=N_{\phi }\left(x\right){\widehat{\phi }}^e\left(t\right),\\ w_1{}_h^e\left(x,t\right)=N_u\left(x\right){\widehat{w}}_1^e\left(t\right)and{w}_2{}_h^e\left(x,t\right)=N_{\phi }\left(x\right){\widehat{w}}_2{}^e\left(t\right)\end{array}$

where $N_u$ and $N_{\phi }$ are the matrices of shape functions in terms of Lagrange polynomials. For an element with n nodes, we have

$N_u=\left[\begin{array}{ccc}{N}_1\left(x\right)& 0& \cdots \\ 0& N_1\left(x\right)& \cdots \end{array}\begin{array}{cc}{N}_n\left(x\right)& 0\\ 0& N_n\left(x\right)\end{array}\right]and{N}_{\phi }=\left[\begin{array}{ccc}{N}_1\left(x\right)& \cdots & N_n\left(x\right)\end{array}\right].$

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

$E_h^e\left(x,t\right)=B_u\left(x\right){\widehat{u}}^e\left(t\right),D_h^e\left(x,t\right)=B_u\left(x\right){\widehat{\dot{u}}}^e\left(t\right)and\nabla {\phi }_h^e\left(x,t\right)=B_{\phi }\left(x\right){\widehat{\phi }}^e\left(t\right).$

Analogously for the test functions,

$\nabla w_1^e\left(x,t\right)=B_u\left(x\right){\widehat{w}}_1^e\left(t\right)and\nabla w_2{}_h^e\left(x,t\right)=B_{\phi }\left(x\right){\widehat{w}}_2\left(t\right),$

where $B_u$ and $B_{\phi }$ are matrices with the global derivatives of the shape functions given for 2D problems and each node i, respectively, by

$B_{u_i}=\left[\begin{array}{c}{N}_{i,x}\\ 0\\ N_{i,y}\end{array}\begin{array}{c}0\\ N_{i,y}\\ N_{i,x}\end{array}\right]and{B}_{{\phi }_i}=\left[\begin{array}{c}{N}_{i,x}\\ N_{i,y}\end{array}\right].$

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 ${\phi }_k^{n+1}$ is replaced in the equation of motion and then solved by the standard Newmark method. This procedure is repeated for each time step $t_n$ in the time interval $\left[\mathrm{0,}T\right]$ with time increments$\Delta t=t_{n+1}-t_n>0$ for $n=\mathrm{0,1,}\dots ,T$. 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 $n+1$ is given by

$\left[{\mathbb{M}}_{\phi }^{}+\Delta t\left({\mathbb{P}}_{\phi }^{n+1}+{\mathbb{K}}_c^{n+1}\right)\right]{\phi }_k^{n+1}={\mathbb{M}}_{\phi }{\phi }_k^n+\Delta t{w}_b^{n+1}.$ (57)

The global operators of the previous equation are obtained by the standard assembling procedure of the following element operators:

${\mathbb{M}}_{\phi }^e=\underset{{\Omega }_e}{\overset{}{{\displaystyle \int }}}{\tilde{\lambda }}_k{N}_{\phi }^T{N}_{\phi }d\Omega ,$ (58)

${\mathbb{P}}_{\phi }^e=\underset{{\Omega }_e}{\overset{}{{\displaystyle \int }}}{g}_{c}l{B}_{\phi }^T\left({\omega }_k{B}_{\phi }\right)d\Omega +\underset{{\Omega }_e}{\overset{}{{\displaystyle \int }}}\frac{g_c}{l}{N}_{\phi }^T{N}_{\phi }d\Omega ,$ (59)

${\mathbb{K}}_c^e=\underset{{\Omega }_e}{\overset{}{{\displaystyle \int }}}{\tilde{\mathcal{H}}}_k\left({\phi }_i^n,{\widehat{u}}_n^e\right)N_{\phi }^T{N}_{\phi }d\Omega ,$ (60)

$w_b^e=\underset{{\Omega }_e}{\overset{}{{\displaystyle \int }}}{\tilde{\mathcal{H}}}_k\left({\phi }_i^n,{\widehat{u}}_n^e\right)N_{\phi }^{T}d\Omega .$ (61)

The approximate solution of the displacement at time step $n+1$, $u^{n+1}$, is computed by solving

$\left[{\alpha }_1\mathbb{M}+{\mathbb{K}}_u+{\alpha }_4{\mathbb{K}}_v\right]u^{n+1}=\mathbb{M}\left[{\alpha }_1{u}^n+{\alpha }_2{\dot{u}}^n+{\alpha }_3{\ddot{u}}^n\right]+w_a+\left[{\alpha }_4{u}^n+{\alpha }_5{\dot{u}}^n+{\alpha }_6{\ddot{u}}^n\right]+\mathbb{M}{f}^{n+1}.$ (62)

The previous global operators are also obtained by the standard assembling procedure of the following element operators:

${\mathbb{M}}_{}^e=\underset{{\Omega }_e}{\overset{}{{\displaystyle \int }}}{\rho }_0{N}_u^T{N}_{u}d\Omega ,$ (63)

${\mathbb{K}}_u^e=\underset{{\Omega }_e}{\overset{}{{\displaystyle \int }}}{B}_u^T:ℂ\left(N_{\phi }{\widehat{\phi }}_k^e\right):B_u^{}d\Omega ,$ (64)

$w_a^e=\underset{{\Omega }_e}{\overset{}{{\displaystyle \int }}}{g}_{c}l{B}_u^{T}sym\left({\displaystyle \sum }_{k=1}^N\left({\omega }_k{B}_{\phi }{\widehat{\phi }}_k^e{}^n\right)\otimes B_{\phi }{\widehat{\phi }}_k^e{}^n\right)d\Omega ,$ (65)

${\mathbb{K}}_v^e=\underset{{\Omega }_e}{\overset{}{{\displaystyle \int }}}\tilde{b}{B}_u^T{B}_{u}d\Omega .$ (66)

After solving the evolution equation for the displacement (Equation 65), the acceleration and velocity vectors are updated by using the following expressions:

${\ddot{u}}^{n+1}={\alpha }_1(u^{n+1}-u^n)-{\alpha }_2{\dot{u}}^n-{\alpha }_3{\ddot{u}}^n,$ (67)

${\dot{u}}^{n+1}={\alpha }_4(u^{n+1}-u^n)-{\alpha }_5{\dot{u}}^n-{\alpha }_6{\ddot{u}}^n,$ (68)

where the coefficients ${\alpha }_j\left(j=\mathrm{1,...,6}\right)$ are given, in terms of the Newmark coefficients $\tilde{\gamma }=0.5$ and $\tilde{\beta }=0.25$, by ${\alpha }_1=\frac{1}{\tilde{\beta }\Delta t^2}$, ${\alpha }_2=\frac{1}{\tilde{\beta }\Delta t}$, ${\alpha }_3=\frac{1-2\tilde{\beta }}{2\tilde{\beta }}$, ${\alpha }_4=\frac{\tilde{\gamma }}{\tilde{\beta }\Delta t}$, ${\alpha }_5=1-\frac{\tilde{\gamma }}{\tilde{\beta }}$ and ${\alpha }_6=\left[1-\frac{\tilde{\gamma }}{2\tilde{\beta }}\right]\Delta t$.

The overall procedure is given in Algorithm 1.

Algorithm 1 Semi-implicit time integration procedure for the system of equations:

1. for$t=0\to T$do

2. for $k=1:N$do

3. Given $u^n$, ${\dot{u}}^n$, ${\phi }_k^{n}and$ ${\tilde{\lambda }}_n$, solve Equation 57 for ${\phi }_k^{n+1}$

4. end

5. Given the updated damage ${\phi }_{k=\mathrm{1,..,}N}^{n+1}$, solve Equation 62 for $u^{n+1}$;

6. Update the acceleration ${\ddot{u}}^{n+1}$and velocity ${\dot{u}}^{n+1}$using Equations 67 and 71;

7. Increment the time step by $\Delta t$.

8. end

3 NUMERICAL RESULTS

3.1 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 $1\times {10}^{-4}\Delta tm/s$ is applied to the upper edge incrementally along the time steps ($\Delta t$). The adopted time increment was $1\times {10}^{-3}s$. The material parameters for this case are the following: Young modulus $E=160GPa$, Poisson coefficient $\nu =0.3$, density ${\rho }_0=7300kg/m^3$, Griffith fracture energy $g_c=2.7\times {10}^{3}N/m$ and $l=1.0\times {10}^{-3}m$. A mesh of linear triangles with $4961$ nodes was used. An initial damage of $\phi =0.999$ was set at some nodes on middle of the left vertical edge of the specimen as illustrated in Figure 1(b).

Figure 1
(a) Mesh used for the I-shaped tensile test. (b) Initial horizontal damage of 0.99 imposed on nodes. (c) Final damage distribution and crack path for the isotropic case. (d) Damage distribution and crack path for the anisotropic propagation. (e) Legend.

First an isotropic crack behavior ($\beta =0$) was considered and illustrated in Figure 1(c). As expected, the material has an isotropic damage behavior and the crack remains at $0^o$. The value $\beta =20$ was used to simulate an anisotropic damage propagation with preferential direction oriented at ${45}^o$ with respect to the horizontal axis. Figure 1(d) shows that the crack growth behaved anisotropically as the damage increased predominantly at ${45}^o.$

3.2 Example2: Notched tensile test

The second example is the single notched tensile test of a $1m\times 1m$ 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 $1\times {10}^{-3}\Delta tm/s$ is applied on the upper edge. The following parameters were used: $E=200GPa$, $\nu =0.3$, ${\rho }_0=7800kg/m^3$, $l=2.5\times {10}^{-3}m$, $b=1$ and $g_C=2.7\times {10}^{3}N/m$. The parameter $\beta$ was first taken zero to recover the isotropic case. For the anisotropic crack propagation with preferential direction oriented at $-{45}^o$, three values were considered,$\beta =10$, $\beta =50$ and $\beta =150$, in order to verify its influence on the crack angle. For this simulation, a time step $\Delta t=1\times {10}^{-3}s$ and a mesh of $12053$linear triangular elements and $6155$nodes were used (Figure 2).

Figure 2
Mesh used for the single notched tensile test.

Figure 3
(a) Final damage distribution for the isotropic case. (b-d) Damage distribution for the anisotropic case.

As shown in Figure 3, larger $\beta$ makes the crack propagation angle closer to the material preferential direction. For $\beta =10$, the crack propagation angle is $-{\mathrm{22,2}}^o$, for $\beta =50$, the angle increases to $-{\mathrm{33,6}}^o$ and for $\beta =150$, it is $-{\mathrm{38,9}}^o.$ 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.

3.3 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 $MgA{l}_2{O}_4$, 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 $2mm$ side was submitted to the boundary conditions illustrated in Figure 4 with incrementally applied displacement of $1\times {10}^{-4}tm/s$. The idea is to restrict the crack path to the central band of $0.8mm$ by imposing the displacement in the upper and lower edges. The domain was discretized with $9086$ linear square elements. The material properties were $E=10GPa$, $\nu =0.2$ and ${\rho }_0=7300kg/m^3$; the model parameters used were $l=1\times {10}^{-5}m$, $b=1$ and $c=1\times {10}^{-4}$. Since the same damage mechanism is in both preferential directions (${45}^o$ and $-{45}^o$), the same parameter $\beta =20$ and Griffith fracture energy $g_c=250N/m$ were adopted for both directions. In order to evaluate the dependency of predicted crack path with the time step, three values were considered: $1\times {10}^{-2}s$, $1\times {10}^{-3}s$ and $5\times {10}^{-4}s$. 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.

Figure 4
Boundary conditions considered for the zigzag crack propagation.

Figure 5
(a-b) Damage for time step $1\times {10}^{-2}s$, (c-d) Damage for time step $1\times {10}^{-3}s$, (e-f) Damage for time step $5\times {10}^{-4}s$.

For all time steps the crack initiates with an angle of$-{45}^o$, which is the preferential direction of the first damage variable. As the crack propagates, the angle becomes smaller for time step $1\times {10}^{-2}s$ and before the crack reaches the boundary condition barrier, it changes the propagation direction. In this case the phase field ${\phi }_1$ is dominant and can represent alone the crack. As can be seen in the Figure 5(a), the expected sharp zig-zag crack pattern could not be obtained for this time step. For the smaller time steps, $1\times {10}^{-3}s$ and $5\times {10}^{-4}s$, the crack represented by the phase-field variable ${\phi }_1$ reaches the bottom of the central band, where the crack kinks, following another energetically preferred direction associated with the second phase field ${\phi }_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 $1\times {10}^{-3}s$ and $5\times {10}^{-4}s$ are qualitatively very similar and the total crack occurs with prescribed displacement $u_{\max }=0.0093mm$ and $u_{\max }=0.0086mm$, respectively. The maximum displacement with time step $1\times {10}^{-2}s$ was $u_{\max }=0.0138mm$.

4 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 $0^o$ and preferential cleavage plane oriented at$-{45}^o$. The influence of parameter $\beta$, that regulates how strong is the fracture anisotropy, was studied for a notched tensile test. For larger $\beta$, 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 ($\le 1\times {10}^{-3}s$), the model was also able to reproduce the zig-zag crack pattern with two preferential directions. For time step $1\times {10}^{-2}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.

Acknowledgements

This work was financed by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brazil (CAPES) - Finance Code 001 and Conselho Nacional de Pesquisa Científica e Tecnológica under grant number 310351/2019-7.

References

• Ambati, M., Gerasimov, T., & De Lorenzis, L. (2014). A review on phase-field models of brittle fracture and a new fast hybrid formulation. Computational Mechanics, 55(2), 383-405.
• Bleyer, J. and Alessi, R. (2018), Phase-field modeling of anisotropic brittle fracture including several damage mechanisms. Computer Methods in Applied Mechanics and Engineering 336, 213-236.
• Brünig, M. (2004). An anisotropic continuum damage model: Theory and numerical analyses. Latin American Journal of Solids and Structures, 1(2), 185-218.
• Bittencourt, M.L. (2014), Computational Solid Mechanics, CRC Press.
• Boldrini, J. L., Moraes, E. A. B. D., Chiarelli, L. R., Fumes, F. G., Bittencourt, M. L. (2016). A non-isothermal thermodynamically consistent phase field framework for structural damage and fatigue. Computer Methods in Applied Mechanics and Engineering 312, 395-427.
• Cahn, J. W. (1959). Free energy of a nonuniform system. II. Thermodynamic basis. The Journal of Chemical Physics 30(5), 1121-1124.
• Chiarelli, L. R., Fumes, F. G., Moraes, E. A. d. B. D., Haveroth, G. A., Boldrini, J. L. and Bittencourt, M. L. (2017). Comparison of high order finite element and discontinuous Galerkin methods for phase field equations: Application to structural damage. Computers and Mathematics with Applications 74(7), 1542-1564.
• Chow, C. L. and Wang, J. (1987). An anisotropic theory of continuum damage mechanics for ductile fracture. Engineering Fracture Mechanics 27(5), 547-558.
• Clayton, J. D. and Knap, J. (2015). Phase field modeling of directional fracture in anisotropic polycrystals. Computational Materials Science 98, 158-169.
• Dal, H., Gültekin, O., Denli, F. A. and Holzapfel, G. A. (2017). Phase-Field Models for the Failure of Anisotropic Continua. PAMM Proc. Appl. Math. Mech. 17, 91 - 94.
• Denli, F. A., Gültekin, O., Holzapfel, G. A., & Dal, H. (2020). A phase-field model for fracture of unidirectional fiber-reinforced polymer matrix composites. Computational Mechanics, 65(4), 1149-1166.
• Fabrizio, M., Giorgi, C., & Morro, A. (2006). A thermodynamic approach to non-isothermal phase-field evolution in continuum physics. Physica D: Nonlinear Phenomena, 214(2), 144-156.
• Fassin, M., Eggersmann, R., Wulfinghoff, S., & Reese, S. (2019). Gradient-extended anisotropic brittle damage modeling using a second order damage tensor - Theory, implementation and numerical examples. International Journal of Solids and Structures, 167, 93-126.
• Frémond, M. (2013). Non-Smooth Thermomechanics. Springer Science & Business Media.
• Gültekin, O., Dal, H., Holzapfel, G. A. (2017). Numerical aspects of anisotropic failure in soft biological tissues favor energy-based criteria: A rate-dependent anisotropic crack phase-field model. Computer Methods in Applied Mechanics and Engineering, 23-52.
• Jarić, J., Kuzmanović, D., & Šumarac, D. (2013). On anisotropic elasticity damage mechanics. International Journal of Damage Mechanics, 22(7), 1023-1038.
• Levitas, V. I., Jafarzadeh, H., Farrahi, G. H., & Javanbakht, M. (2018). Thermodynamically consistent and scale-dependent phase field approach for crack propagation allowing for surface stresses. In International Journal of Plasticity (Vol. 111). Elsevier Ltd.
• Li, B., Peco, C., Millán, D., Arias, I. and Arroyo, M. (2015). Phase-field modeling and simulation of fracture in brittle materials with strongly anisotropic surface energy. International Journal for Numerical Methods in Engineering 102, 711-727.
• Miehe, C., Hofacker, M. and Welschinger, F. (2010). A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering 199(45-48), 2765-2778.
• Mozaffari, N. and Voyiadjis, G. Z. (2015). Phase field based nonlocal anisotropic damage mechanics model. Physica D: Nonlinear Phenomena 308, 11-25.
• Nguyen, T. T., Réthoré, J. and Baietto, M.-C. (2017). Phase field modelling of anisotropic crack propagation. European Journal of Mechanics - A/Solids 65, 279-288.
• Réthoré, J., Dang, T. B. T., & Kaltenbrunner, C. (2017). Anisotropic failure and size effects in periodic honeycomb materials: A gradient-elasticity approach. Journal of the Mechanics and Physics of Solids, 99, 35-49.
• Talreja, R. (2014). Assessment of the fundamentals of failure theories for composite materials. Composites Science and Technology 105, 190-201.
• Teichtmeister, S. and Miehe, C. (2015). Phase-Field Modeling of Fracture in Anisotropic Media. PAMM 15(1), 159-160.
• Wu, C. C., McKinney, K. R. and Rice, R. W. (1995). Zig-zag crack propagation in MgAl2O4 crystals. Journal of Materials Science Letters 14(7), 474-477.
• Wu, J. Y., Nguyen, V. P., Zhou, H., & Huang, Y. (2020). A variationally consistent phase-field anisotropic damage model for fracture. Computer Methods in Applied Mechanics and Engineering, 358, 112629.

Edited by

Guest Editors:

Volnei Tita and Nicholas Fantuzzi.

Publication Dates

• Publication in this collection
  12 Feb 2021
• Date of issue
  2021

History

• Received
  30 Jan 2020
• Reviewed
  16 July 2020
• Accepted
  28 Aug 2020
• Published
  02 Sept 2020