Brittle Anisotropy based on a Tensor Damage Phase Field Model

Abstract

A consistent model that properly captures the behavior of materials when submitted to loads and at the same time considers the influence of damage growth in the material properties is a challenging task. Specially in brittle materials, due to the mechanism of cleavage, as the damage grows, the material mechanical response is differently degraded according to the load direction considered, characterizing an anisotropic damage. A fourth-order degradation tensor is used with a phase field framework in order to model the induced damage anisotropy without the need of defining the damage principal direction a priori. The transient non-linear coupled system of equations is solved using appropriate time integration procedures. Due to the mild non-linearities, standard linearization methods are not considered. The numerical examples illustrate the model’s features and conclude that the adopted approach is able to simulate anisotropic damage totally driven by the strain state of the material.

Keywords
Phase field; Anisotropic damage; Damage-induced anisotropy; Finite element method

Graphical Abstract

1 INTRODUCTION

As damage increases, the mechanical properties of materials are affected and this influence can manifest differently depending on the direction analyzed. This phenomenon is called damage induced anisotropy (Jaric et al., 2013; Brunig, 2004; Chow and Wang, 1987).

Since damage is a localized phenomenon, only the properties in the nearby zone are affected. This means that initially isotropic materials become locally anisotropic in the presence of damage and that non-isotropic materials have their original anisotropy changed. Due to the geometry of microdefects, this phenomenon is more pronounced for flat defects as cleavage in brittle materials. Nevertheless, the mechanical behavior of the damaged body may depend on the direction and distribution of microcracks (Skrzypek and Ganczarski, 2015). Such anisotropy of current damage state can also influence its development when non-proportional loading is considered (Chaboche, 1981).

Mathematical models have been developed aiming to describe realistically the mechanical, thermal and chemical behavior of many types of materials. In the context of damage and fracture models, the phase field approach is a competitive alternative to the conventional discrete methods (Haveroth et al., 2018; Petrini et al., 2020; Vale et al., 2023). However, anisotropic damage has been mainly addressed by models based on the Continuum Damage Mechanics (CDM). The description of damage by an internal variable is the base of both theories (Clayton and McDowell, 2003).

The difference between phase field and classical CDM models is that the onset and continuation of damage is governed by the evolution equation of the damage phase field variable considering an energetically motivated fracture criterion, which depends on the strain state generated by the external loading instead of phenomenological kinetic equations specifying the rate of damage accumulation, as pointed out by Clayton and Knap (2014). The model considered in this work uses the ideas of CDM into a phase field framework in order to obtain a general anisotropic damage model.

Some models, mainly of fourth-order, have focused on reproducing some specific anisotropic behavior of damaged materials by defining its final symmetry (Jaric et al., 2013). In most of these models, the principal directions are defined from the beginning of the simulation based on the load direction, and all points of the material have the same damage orientation (Mozaffari and Voyiadjis, 2015), rather than different ones according to the local strain state. Some models adopt the principal strains as the principal damage directions (Pituba, 2003). However, in such cases, each component of the degrading tensor are scalar functions similarly to multi-scalar damage approaches. To account for different directions, a rotation matrix is also included in some works (Chow and Wang, 1987). But, if the loading changes over time, it is still not clear how the scalar damage variables are evolved to represent the damage associated with the new direction.

The present model of this work uses a fourth-order degradation tensor, defined in the global coordinate system, to represent the damage effect on the material elasticity. The fourth-order degradation tensor is thermodynamically consistent according to the damage phase field framework proposed by Boldrini et al. (2016). It is an internal variable with its own evolution equation with an energetically motivated fracture criterion. The governing equations were obtained based on the principle of virtual power (PVP) and on the balance of energy and the Clausius–Duhem inequality for the entropy. The degradation tensor affects the elasticity tensor using the CDM hypothesis of complementary elastic energy equivalence, Petrini (2023).

The transient system of equation is given by the equation of motion coupled with the equation for the evolution of the degradation tensor. They are solved by appropriate implicit time integration procedures. Due to the mild non-linearities, explicit linearization of both equations is not considered. The finite element method (FEM) is used for spatial approximation of equations. Plane stress state is considered for the examples.

The effects of tension-compression asymmetry and damage irreversibility were considered in the numerical examples. However, the details of these features are not here presented and more details can be found in Petrini (2023).

This paper is organized as follows. The next two sections present an overview of the phase field model based on the degradation tensor and the time and spatial discretizations by the FEM of the motion and degradation equations. Various cases are considered illustrating the main features of the proposed model. Finally, conclusions are drawn from the obtained results.

2 TENSOR DAMAGE MODEL

To include the damage induced anisotropy, the stiffness loss due to the damaging process is given by the symmetric degradation tensor ($\mathbb{G}$), introduced in a phase field framework combined with the ideas from the CDM.

The traditional phase field method allows the modeling of discontinuous interfaces of cracks as a continuum by introducing the concept of smeared cracks represented by an additional scalar variable $\mathrm{d}$. The phase field $\mathrm{d}$ assumes value 0 for virgin material; 1 for fractured material; and 0 < $\mathrm{d}$ < 1 is associated to a damaged material between these two extreme states meaning the volumetric fraction of damaged material. The degradation of the material mechanical response is originally isotropic and given by the degradation function, which depends of $\mathrm{d}$ and ranges between 1 and 0 for scalar models and multiplies the material stiffness, degrading it in an isotropic manner.

Based on a thermodynamically consistent non-isothermal damage phase field framework presented by Boldrini et al. (2016), we substituted the scalar phase field by a fourth-order variable, the degradation tensor $\mathbb{G}$. In this case, the material degradation is no more given by the scalar degradation function and the elastic energy of the damaged material is obtained by the strain energy equivalence principle of the CDM.

Due to the nature of brittle damage, $\mathbb{G}$ is considered an internal-type variable, not allowing sources of damage from the external environment to be included. The governing equations were obtained based on the PVP, the balance of energy and the Clausius–Duhem inequality for the entropy. This methodology has the advantage of not requiring the definition of degradation functions. Moreover, general anisotropy can be achieved.

The governing equations obtained are written in terms of the specific density of the free-energy $\mathrm{ø}$ which can be defined such that it is possible to obtain a thermodynamically consistent fracture models for many types of materials.

To include the possibility of damage, we take the expression of the free-energy density, $\mathrm{ø}$, as the sum of the elastic energy density of the damaged body, $\mathcal{E}(\mathbb{G},$E), and the energy density required to create the crack, $\mathcal{I}(\mathbb{G},\nabla \mathbb{G})$. Therefore,

$\rho \mathrm{}\psi (\mathbb{G},\nabla \mathbb{G},\mathrm{}E)\mathcal{}=\mathcal{}\mathcal{E}(\mathbb{G},E)\mathcal{}+\mathcal{I}(\mathbb{G},\nabla \mathbb{G})$ (1)

where $\mathrm{ñ}$ is the volumetric density.

The degraded elastic energy $\mathcal{E}(\mathbb{G},\mathrm{E})$ is determined based on the CDM. The main idea is that the damaged material can be represented by a fictitious undamaged state whose mechanical behavior is described by the effective stress. The effective stress and strain are those on the bulk material of the damaged state taking into account the area reduction and the crack opening (Pituba, 2003). To obtain a simple but rather general model, we use these ideas and take the relation between the actual damaged state and the effective state given by the hypothesis of complementary strain energy equivalence. These CDM ideas lead us to write the elastic energy ($\mathcal{E}$) for the actual damage state, in terms of the Hooke’s laws for small strains, as

$\mathcal{E}(\mathbb{G},\mathrm{E})\mathrm{}=\mathrm{}\frac{1}{2}\mathrm{E}\mathrm{}:\mathbb{G}\mathrm{}:\mathrm{}{\mathbb{C}}^0\mathrm{}:\mathrm{}\mathbb{G}\mathrm{}:\mathrm{}\mathrm{E}\mathrm{}.$ (2)

Equation (2) is written in coordinates associated to an orthonormal basis; ${\mathrm{G}}_{\mathrm{i}\mathrm{j}\mathrm{k}\mathrm{l}}$ denotes the components of $\mathbb{G}$, and to guarantee the usual standard symmetries of the effective stress tensor, it is required that they satisfy ${\mathrm{G}}_{\mathrm{j}\mathrm{i}\mathrm{k}\mathrm{l}}{=\mathrm{G}}_{\mathrm{i}\mathrm{j}\mathrm{k}\mathrm{l}}$, ${\mathrm{G}}_{\mathrm{i}\mathrm{j}\mathrm{l}\mathrm{k}}{=\mathrm{G}}_{\mathrm{i}\mathrm{j}\mathrm{k}\mathrm{l}}$ and ${\mathrm{G}}_{\mathrm{k}\mathrm{l}\mathrm{i}\mathrm{j}}{=\mathrm{G}}_{\mathrm{i}\mathrm{j}\mathrm{k}\mathrm{l}}$. Moreover, ${\mathbb{C}}^0$ denotes the elasticity tensor of the virgin material and $\mathrm{E}=\frac{1}{2}(\nabla \mathrm{u}\mathrm{}+{\nabla }^{\mathrm{T}}\mathrm{u}\mathrm{})$ is the infinitesimal strain tensor field obtained from the displacement vector field $\mathrm{u}$. It is worth mentioning that, in a state without damage, $\mathbb{G}$ is taken as the fourth-order symmetric identity tensor ${\mathbb{I}}_{\mathrm{s}}$ whose components are given by ${\left({\mathbb{I}}_{\mathrm{s}}\right)}_{\mathrm{i}\mathrm{j}\mathrm{k}\mathrm{l}}\mathrm{}=\frac{1}{2}({\mathrm{ä}}_{\mathrm{i}\mathrm{k}}{\mathrm{ä}}_{\mathrm{j}\mathrm{k}}+{\mathrm{ä}}_{\mathrm{i}\mathrm{l}}{\mathrm{ä}}_{\mathrm{j}\mathrm{k}})$.

Since the effect of damage growth on the material response is represented by the degradation tensor $\mathbb{G}$, its contribution to the free-energy density is introduced by the term $\mathcal{I}(\mathbb{G},\nabla \mathbb{G})$ defined as

$\mathcal{I}(\mathbb{G},\nabla \mathbb{G})=\frac{{\mathrm{g}}_{\mathrm{c}}}{2\mathrm{ã}}\mathrm{}({\mathbb{I}}_{\mathrm{s}}-\mathbb{G})::({\mathbb{I}}_{\mathrm{s}}-\mathbb{G})\mathrm{}+\frac{{\mathrm{g}}_{\mathrm{c}}\mathrm{ã}}{2}\mathrm{}(\nabla \mathbb{G}):\cdot :(\nabla \mathbb{G}).$ (3)

Here, $\mathrm{ã}$ is a positive constant related to the width of the fracture layers and ${\mathrm{g}}_{\mathrm{c}}$ is the critical Griffith fracture energy assumed to be a positive constant. $\mathcal{I}$ was defined inspired on phase field models such as (Miehe et al. 2010), that define a surface crack density function for isotropic scalar damage based on a variational approach.

By considering the free energy defined in (1), with the expressions (2) and (3), in the case of an isothermal state and a nearly-incompressible and homogeneous material, the governing equations are given by (Petrini, 2023)

$\begin{array}{c}\ddot{\mathrm{u}}=\frac{1}{{\mathrm{ñ}}_0}divT+f,\\ T=G:{\mathbb{C}}^0:G:E-sym({\mathrm{g}}_{\mathrm{c}}ã\nabla G::\nabla G)+\widehat{\mathrm{b}}D,\\ \dot{\mathbb{G}}=-\widehat{\mathrm{F}}{\left[\mathrm{E}⨂\left({\mathbb{C}}^0\mathrm{}:\mathrm{}\mathbb{G}\mathrm{}:\mathrm{}\mathrm{E}\right)\right]}^{\tilde{\mathrm{s}}}-\widehat{\mathrm{F}}\frac{{\mathrm{g}}_{\mathrm{c}}}{\mathrm{ã}}(G-{\mathbb{I}}_{\mathrm{s}})+\widehat{\mathrm{F}}{\mathrm{g}}_{\mathrm{c}}ãdiv(\nabla G),\end{array}$ (4)

where f denotes the body force vector field, D is the symmetric part of the velocity gradient tensor field, $\mathrm{D}=\frac{1}{2}(\nabla \mathrm{v}\mathrm{}+{\nabla }^{\mathrm{T}}\mathrm{v}\mathrm{})$, T is the Cauchy stress tensor field, and $\widehat{\mathrm{b}}$ and $\widehat{\mathrm{F}}$ are positive constants. The parameter $\widehat{\mathrm{b}}$ is related to the elastic damping and $\widehat{\mathrm{F}}$ controls how much the strain state influences in the degradation rate.

Another important aspect to be considered, as the model does not have an explicit damage variable (due to the use of $\mathbb{G}$), is to define a suitable damage criterion. From $\mathbb{G}$, such a criterion must give the level of damage at any point in the material body and be used as an irreversibility damage criterion.

According to Neeman et al. (2008), there is a relationship between the direction of the eigenvector associated to the smallest eigenvalue and the direction of fracture for brittle materials. Based on that, we established that the degradation tensor would locally affect the elasticity tensor by lowering one of its eigenvalues and that failure would happen when one of them reached a value close to zero. Despite needing only one eigenvalue, there is no efficient way to track which eigenvalue reaches zero consistently. Not only that, but there might be changes in the other eigenvalues that also represent the damage level. Therefore, we chose to use the geometric mean of all eigenvalues to obtain an overall representation of the damage state. This approach also preserves the ability to keep track whenever one of the eigenvalues is close to zero.

Based on these considerations, we define $\mathbf{C}$ as the matrix representation of the reduced degraded elasticity $\mathbb{C}$ for plane stress and similarly $\mathbf{\epsilon }$ is the vector representation of the reduced strain tensor. Recalling the min-max theorem, the minimum eigenvalue ${\mathrm{ë}}_1$for the degraded elasticity can be obtained by:

${\mathrm{ë}}_1\mathrm{}=\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{}\left\{\frac{\mathbf{\epsilon }\cdot \mathbf{C}\mathbf{\epsilon }}{{\left|\right|\mathbf{\epsilon }\left|\right|}^2}\mathrm{}:\mathrm{}\mathrm{å}\ne 0\right\}.$ (5)

This eigenvalue represents the elastic energy for any strain state. This can be understood as the strain that minimizes the elastic energy and therefore the eigenvalues (and their geometric mean) are seen as a measure of energy density, as also described in Mehrabadi and Cowin (1991) which assumed a unitary strain state.

In the present work, the damage irreversibility criterion is constructed as follows. We compute the eigenvalues of the original elastic tensor ${\mathbb{C}}^0$ denoted by ${{\mathrm{ë}}_1}^{\left(0\right)}\le {{\mathrm{ë}}_2}^{\left(0\right)}\le {{\mathrm{ë}}_3}^{\left(0\right)}$. Next, consider the degraded elastic tensor $\mathbb{C}=\mathbb{G}\mathrm{}:\mathrm{}{\mathbb{C}}^0\mathrm{}:\mathrm{}\mathbb{G}$ and its eigenvalues, say ${\mathrm{ë}}_1\le {\mathrm{ë}}_2\le {\mathrm{ë}}_3$. Then, we compute the quotient of the geometric averages of the eigenvalues of the degraded and the original elastic tensors to obtain:

$\mathrm{\ae }\mathrm{}=\mathrm{}\frac{{\left({\mathrm{ë}}_1{\mathrm{ë}}_2{\mathrm{ë}}_3\right)}^{1/3}}{{\left({{\mathrm{ë}}_1}^{\left(0\right)}{{\mathrm{ë}}_2}^{\left(0\right)}{{\mathrm{ë}}_3}^{\left(0\right)}\right)}^{1/3}}$.(6)

Due to the degradation effect of $\mathbb{G}$, we have that $0\le \mathrm{\ae }\le 1$. This gives a scalar field with the required damage level as in standard models with one scalar phase field; in particular, when $\mathrm{\ae }=0$ the material is in broken state. To impose damage irreversibility, it is enough to require that $\dot{\mathrm{\zeta }}\le 0$; that is, in the approximated version of the model, $\mathbb{G}$ is updated only when the discretized version of $\dot{\mathrm{\ae }}\le 0$ is satisfied.

By considering a plane stress state, we also obtain a plane damage state. In this case, there are six components of $\mathbb{G}$ to be evolved independently. Defining ${\mathrm{G}}_{11}:={\mathrm{G}}_{1111}$,${\mathrm{G}}_{12}:={\mathrm{G}}_{1122}$, ${\mathrm{G}}_{13}:={\mathrm{G}}_{1112}$, ${\mathrm{G}}_{22}:={\mathrm{G}}_{2222}$, ${\mathrm{G}}_{23}:={\mathrm{G}}_{1222}$, ${\mathrm{G}}_{33}:={\mathrm{G}}_{1212}$, we can define the matrix form of G for plane stress as:

$\left[{\mathbb{G}}^{\mathrm{p}\mathrm{s}}\right]=\left[\begin{array}{ccc}{\mathrm{G}}_{11}& {\mathrm{G}}_{12}& {\mathrm{G}}_{13}\\ {\mathrm{G}}_{12}& {\mathrm{G}}_{22}& {\mathrm{G}}_{23}\\ {\mathrm{G}}_{13}& \mathrm{}{\mathrm{G}}_{23}& {\mathrm{G}}_{33}\end{array}\right].$ (7)

As $\mathbb{G}$ has both major and minor symmetries, while solving the equation of the degradation evolution, it is more interesting to write it in a vector form as

$\left\{{\mathbb{G}}^{\mathrm{p}\mathrm{s}}\right\}={\left\{{\mathrm{G}}_{11},\mathrm{}{\mathrm{G}}_{22},\mathrm{}{\mathrm{G}}_{33},\mathrm{}{\mathrm{G}}_{12},\mathrm{}{\mathrm{G}}_{13},\mathrm{}{\mathrm{G}}_{23}\right\}}^{\mathrm{T}}.$ (8)

The final governing equations for plane-stress are the following:

$\left\{\begin{array}{c}\left\{\ddot{\mathrm{u}}\right\}=\frac{1}{{\mathrm{\rho }}_0}div\left\{\mathrm{T}\right\}+\left\{f\right\},\\ \left\{T\right\}=\left[{\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}\mathrm{}:\mathrm{}{\mathbb{C}}^0\mathrm{}:\mathrm{}{\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}\right]:\left\{\mathrm{E}\right\}-\left\{\mathrm{s}\mathrm{y}\mathrm{m}({\mathrm{g}}_{\mathrm{c}}\mathrm{}\mathrm{\gamma }\nabla {\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}::\nabla {\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}})\right\}+\widehat{\mathrm{b}}\left\{\mathrm{D}\right\},\\ \left\{{\dot{\mathbb{G}}}^{\mathrm{p}\mathrm{s}}\right\}=-\widehat{\mathrm{F}}\left\{{\left[\mathrm{E}⨂\left({\mathbb{C}}^0\mathrm{}:\mathrm{}{\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}\mathrm{}:\mathrm{}\mathrm{E}\right)\right]}^{\tilde{\mathrm{s}}}\right\}-\widehat{\mathrm{F}}\frac{{\mathrm{g}}_{\mathrm{c}}}{\mathrm{\gamma }}(\left\{{\mathbb{G}}^{\mathrm{p}\mathrm{s}}\right\}-d)+\widehat{\mathrm{F}}{\mathrm{g}}_{\mathrm{c}}\gamma \left\{\mathrm{d}\mathrm{i}\mathrm{v}(\nabla {\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}})\mathrm{}\right\},\end{array}\right.$ (9)

where $\mathrm{d}\mathrm{}={\left\{1\mathrm{}1\mathrm{}0.5\mathrm{}0\mathrm{}0\mathrm{}0\right\}}^{\mathrm{T}}$. Note that curly braces { } were used for vectors and brackets [ ] for matrices in Voigt notation in order to distinguish it from the original tensors.

As a consequence of the procedure used to deduce the governing equations (Eq. (4)) for the plane stress state, the strain component ${\mathrm{E}}_{33}$ can be calculated as a function of ${\mathrm{E}}_{11}$, ${\mathrm{E}}_{22}$ and the components of [$\mathbb{G}$]. As for isotropic materials, this relation is taken into account while constructing the stiffness matrix of the plane stress state.

Additionally, we adopt an energy decomposition criterion where ${\dot{\mathbb{G}}}^{\mathrm{p}\mathrm{s}}=0$ when the material is under compression. To determine if an element is under tension or compression, the sign of the strain tensor is analyzed.

In the next section, the numerical implementation of the system of equations given in Eq. (9) is described.

3 NUMERICAL APROXIMATION

The time interval $[0,\mathrm{T}]$ is split into $\mathrm{N}+1$ equally spaced time instants ${\mathrm{t}}_{\mathrm{i}}\mathrm{}(\mathrm{i}=0,...,\mathrm{N})$. All variables of the problem are assumed to be known at time ${\mathrm{t}}_{\mathrm{n}}$ and we want to determine them at time ${\mathrm{t}}_{\mathrm{n}+1}$, including the displacement and degradation fields, $\mathrm{u}(\mathrm{x},{\mathrm{t}}_{\mathrm{n}+1})$ and $\mathbb{G}(\mathrm{x},{\mathrm{t}}_{\mathrm{n}+1})$. By using an appropriated time integration method, we solve each equation separately. First, we solve the equation of motion to obtain the updated deformation field and then we use this result as input in the degradation equation. The two-dimensional spatial discretization for plane stress state is done using the FEM.

By assuming infinitesimal deformation, $\mathrm{\Omega }\simeq {\mathrm{\Omega }}_{\mathrm{n}}\simeq {\mathrm{\Omega }}_{\mathrm{n}+1}$, and a test function $\mathrm{w}$, the weak formulation of the equation of motion, Eq. (9a), after the substitution of $\mathrm{T}$ by Eq. (9b) and integration by parts is

${\int }_{\mathrm{\Omega }}\ddot{\mathrm{u}}\cdot \mathrm{w}\mathrm{}=\mathrm{}-\frac{1}{\mathrm{\rho }}{\int }_{\mathrm{\Omega }}{\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}\mathrm{}:\mathrm{}{\mathbb{C}}^0\mathrm{}:\mathrm{}{\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}:\mathrm{}\mathrm{E}\mathrm{}:\nabla \mathrm{w}+\frac{1}{\mathrm{\rho }}{\mathrm{g}}_{\mathrm{c}}\mathrm{\gamma }{\int }_{\mathrm{\Omega }}\left(\nabla {\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}::\nabla {\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}\right):\nabla \mathrm{w}-\frac{1}{\mathrm{\rho }}\widehat{\mathrm{b}}{\int }_{\mathrm{\Omega }}\mathrm{D}\mathrm{}:\nabla \mathrm{w}+{\int }_{\mathrm{\Omega }}\mathrm{f}\mathrm{}\cdot \mathrm{w}+\mathrm{}\frac{1}{\mathrm{\rho }}{\int }_{{\mathrm{\Gamma }}_{\mathrm{t}}}\mathrm{t}\mathrm{}\cdot \mathrm{w},$ (10)

where $\mathrm{t}$ is the surface traction on the boundary ${\mathrm{}\mathrm{\Gamma }}_{\mathrm{t}}$.

The conforming global approximations for the displacement, velocity, acceleration, and test function in a mesh of K elements are indicated as

$\mathrm{u}(\mathrm{x},\mathrm{t})\approx {\bigcup }_{\mathrm{e}=1}^{\mathrm{K}}{\mathrm{u}}_{\mathrm{h}}^{\mathrm{e}}(\mathrm{x},\mathrm{t}),\mathrm{}\dot{\mathrm{u}}(\mathrm{x},\mathrm{t})\approx {\bigcup }_{\mathrm{e}=1}^{\mathrm{K}}{\dot{\mathrm{u}}}_{\mathrm{h}}^{\mathrm{e}}(\mathrm{x},\mathrm{t}),$ $\ddot{\mathrm{u}}(\mathrm{x},\mathrm{t})\approx {\bigcup }_{\mathrm{e}=1}^{\mathrm{K}}{\ddot{\mathrm{u}}}_{\mathrm{h}}^{\mathrm{e}}(\mathrm{x},\mathrm{t}),\mathrm{}\mathrm{w}(\mathrm{x},\mathrm{t})\approx {\bigcup }_{\mathrm{e}=1}^{\mathrm{K}}{\mathrm{w}}_{\mathrm{h}}^{\mathrm{e}}(\mathrm{x},\mathrm{t}),$ (11)

The approximations in each element is given by the linear combinations of the local nodal basis functions as follows:

${\mathrm{u}}_{\mathrm{h}}^{\mathrm{e}}\mathrm{}={\mathrm{N}}_{\mathrm{u}}{\widehat{\mathrm{u}}}^{\mathrm{e}},\mathrm{}{\dot{\mathrm{u}}}_{\mathrm{h}}^{\mathrm{e}}\mathrm{}={\mathrm{N}}_{\mathrm{u}}{\mathrm{}\widehat{\dot{\mathrm{u}}}}^{\mathrm{e}},\mathrm{}{\ddot{\mathrm{u}}}_{\mathrm{h}}^{\mathrm{e}}\mathrm{}={\mathrm{N}}_{\mathrm{u}}{\mathrm{}\widehat{\ddot{\mathrm{u}}}}^{\mathrm{e}},\mathrm{}{\mathrm{w}}_{\mathrm{h}}^{\mathrm{e}}\mathrm{}={\mathrm{N}}_{\mathrm{u}}{\mathrm{}\widehat{\mathrm{w}}}^{\mathrm{e}}\mathrm{},\mathrm{}{\left\{\mathrm{E}\right\}}^{\mathrm{e}}\mathrm{}={\mathrm{B}}_{\mathrm{u}}{\mathrm{}\widehat{\mathrm{u}}}^{\mathrm{e}},\mathrm{}\nabla {\mathrm{w}}^{\mathrm{e}}\mathrm{}={\mathrm{B}}_{\mathrm{u}}{\mathrm{}\widehat{\mathrm{w}}}^{\mathrm{e}},\mathrm{}{\left\{\mathrm{D}\right\}}^{\mathrm{e}}\mathrm{}={\mathrm{B}}_{\mathrm{u}}{\mathrm{}\widehat{\dot{\mathrm{u}}}}^{\mathrm{e}},$ (12)

where ${\mathrm{N}}_{\mathrm{u}}$ is the matrix with the local shape functions and ${\mathrm{B}}_{\mathrm{u}}$ is the matrix with the global derivatives of the shape functions.

Substituting the previous approximations in Eq. (10), we can obtain the following semi-discrete system of equation at element level:

${\mathbb{M}}_{\mathrm{u}}^{\mathrm{e}}{\mathrm{}\widehat{\ddot{\mathrm{u}}}}^{\mathrm{e}}\mathrm{}={\mathbb{K}}_{\mathrm{u}}^{\mathrm{e}}{\mathrm{}\widehat{\mathrm{u}}}^{\mathrm{e}}+\mathrm{}{\mathbb{K}}_{\mathrm{v}}^{\mathrm{e}}{\mathrm{}\widehat{\dot{\mathrm{u}}}}^{\mathrm{e}}+{\mathrm{w}}_{\mathrm{a}}^{\mathrm{e}}.$ (13)

where

${\mathbb{M}}_{\mathrm{u}}^{\mathrm{e}}={\int }_{\mathrm{\Omega }}{{\mathrm{N}}_{\mathrm{u}}}^{\mathrm{T}}{\mathrm{N}}_{\mathrm{u}},\mathrm{}\mathrm{}$ (14)

${\mathbb{K}}_{\mathrm{u}}^{\mathrm{e}}=-\frac{1}{{\mathrm{\rho }}_0}{\int }_{\mathrm{\Omega }}{{\mathrm{B}}_{\mathrm{u}}}^{\mathrm{T}}\left[{\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}\mathrm{}:\mathrm{}{\mathbb{C}}^0\mathrm{}:\mathrm{}{\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}\right]{\mathrm{B}}_{\mathrm{u}},$ (15)

${\mathbb{K}}_{\mathrm{v}}^{\mathrm{e}}=-\frac{1}{{\mathrm{\rho }}_0}\widehat{\mathrm{b}}{\int }_{\mathrm{\Omega }}{{\mathrm{B}}_{\mathrm{u}}}^{\mathrm{T}}{\mathrm{B}}_{\mathrm{u}},\mathrm{}$ (16)

${\mathrm{w}}_{\mathrm{a}}^{\mathrm{e}}=\frac{1}{{\mathrm{\rho }}_0}{\mathrm{g}}_{\mathrm{c}}\mathrm{\gamma }{\int }_{\mathrm{\Omega }}{{\mathrm{B}}_{\mathrm{u}}}^{\mathrm{T}}\left\{{\nabla {\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}}^{\mathrm{T}}::\nabla {\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}\right\},$ (17)

The global system of equations is obtained using the standard assembling procedure of the element operators and given by

${\mathbb{M}}_{\mathrm{u}}\mathrm{}\widehat{\ddot{\mathrm{u}}}={\mathbb{K}}_{\mathrm{u}}\widehat{\mathrm{u}}+{\mathbb{K}}_{\mathrm{v}}\widehat{\dot{\mathrm{u}}}+{\mathrm{w}}_{\mathrm{a}}.$ (18)

In order to approximate the solution of the displacement ${\mathrm{u}}_{\mathrm{n}+1}$, we adopt the implicit standard Newmark procedure. In this method, the acceleration and velocity are evaluated using the updated values of the displacement as

${\widehat{\ddot{\mathrm{u}}}}_{\mathrm{n}+1}={\mathrm{\alpha }}_1({\widehat{\mathrm{u}}}_{\mathrm{n}+1}-{\widehat{\mathrm{u}}}_{\mathrm{n}})-{\mathrm{\alpha }}_2{\mathrm{}\widehat{\dot{\mathrm{u}}}}_{\mathrm{n}}-{\mathrm{\alpha }}_3{\widehat{\ddot{\mathrm{u}}}}_{\mathrm{n}},$ (19)

${\widehat{\dot{\mathrm{u}}}}_{\mathrm{n}+1}={\mathrm{á}}_4({\widehat{\mathrm{u}}}_{\mathrm{n}+1}-{\widehat{\mathrm{u}}}_{\mathrm{n}})+{\mathrm{á}}_5{\mathrm{}\widehat{\dot{\mathrm{u}}}}_{\mathrm{n}}+{\mathrm{á}}_6{\widehat{\ddot{\mathrm{u}}}}_{\mathrm{n}}$,(20)

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

Substituting them into Eq. (18) and considering the body forces equal to zero, we obtain the system of equations to be solved at each time step as

$\left[{\mathrm{\alpha }}_1{\mathbb{M}}_{\mathrm{u}}-{\mathbb{K}}_{\mathrm{u}}-{\mathrm{\alpha }}_4{\mathbb{K}}_{\mathrm{v}}\right]{\widehat{\mathrm{u}}}_{\mathrm{n}+1}={\mathbb{M}}_{\mathrm{u}}\left[{\mathrm{\alpha }}_1{\widehat{\mathrm{u}}}_{\mathrm{n}}+{\mathrm{\alpha }}_2{\mathrm{}\widehat{\dot{\mathrm{u}}}}_{\mathrm{n}}+{\mathrm{\alpha }}_3{\widehat{\ddot{\mathrm{u}}}}_{\mathrm{n}}\right]+{\mathbb{K}}_{\mathrm{v}}\left[-{\mathrm{\alpha }}_4{\widehat{\mathrm{u}}}_{\mathrm{n}}+{\mathrm{\alpha }}_5{\mathrm{}\widehat{\dot{\mathrm{u}}}}_{\mathrm{n}}+{\mathrm{\alpha }}_6{\widehat{\ddot{\mathrm{u}}}}_{\mathrm{n}}\right]+{\mathrm{w}}_{\mathrm{a}}.$ (21)

In order to solve the degradation equation, we need to rewrite Eq. (9c) as follows:

$\left\{{\dot{\mathbb{G}}}^{\mathrm{p}\mathrm{s}}\right\}\mathrm{}=\mathrm{}-\widehat{\mathrm{F}}\left[\mathrm{M}\right]\left\{{\mathbb{G}}^{\mathrm{p}\mathrm{s}}\right\}-\widehat{\mathrm{F}}\frac{{\mathrm{g}}_{\mathrm{c}}}{\mathrm{ã}}(\left\{{\mathbb{G}}^{\mathrm{p}\mathrm{s}}\right\}-\mathrm{d})\mathrm{}+\mathrm{}\widehat{\mathrm{F}}{\mathrm{g}}_{\mathrm{c}}\mathrm{ã}\left\{\mathrm{d}\mathrm{i}\mathrm{v}(\nabla {\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}})\right\}\mathrm{}.$ (22)

This is done to factorize the components of the damage tensor. The components of matrix [$\mathrm{M}$] are such that $\left\{{\left[\mathrm{E}⨂\left({\mathbb{C}}^0\mathrm{}:\mathrm{}{\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}\mathrm{}:\mathrm{}\mathrm{E}\right)\right]}^{\tilde{\mathrm{s}}}\right\}=\left[\mathrm{M}\right]\left\{{\mathbb{G}}^{\mathrm{p}\mathrm{s}}\right\}$.

Multiplying the previous equation by the test function ${\mathrm{w}}_2$ and integrating by parts on the domain $\mathrm{\Omega }$, we derive the weak form of the degradation equation and by applying the backward Euler procedure, we obtain

${\int }_{\mathrm{\Omega }}{\left\{{\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}\right\}}^{\mathrm{n}+1}\cdot {\mathrm{w}}_2+\mathrm{\Delta }\mathrm{t}\mathrm{}\widehat{\mathrm{F}}{\int }_{\mathrm{\Omega }}{\left[\mathrm{M}\right]\left\{{\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}\right\}}^{\mathrm{n}+1}\cdot {\mathrm{w}}_2+\mathrm{\Delta }\mathrm{t}\widehat{\mathrm{F}}\frac{{\mathrm{g}}_{\mathrm{c}}}{\mathrm{\gamma }}{\int }_{\mathrm{\Omega }}{\left\{{\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}\right\}}^{\mathrm{n}+1}\cdot {\mathrm{w}}_2+\mathrm{\Delta }\mathrm{t}\widehat{\mathrm{F}}{\mathrm{g}}_{\mathrm{c}}\mathrm{\gamma }{\int }_{\mathrm{\Omega }}{\nabla \left\{{\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}\right\}}^{\mathrm{n}+1}:{\nabla \mathrm{w}}_2=$

$\mathrm{\Delta }\mathrm{t}\widehat{\mathrm{F}}\frac{{\mathrm{g}}_{\mathrm{c}}}{\mathrm{\gamma }}{\int }_{\mathrm{\Omega }}\mathrm{d}\cdot {\mathrm{w}}_2+{\int }_{\mathrm{\Omega }}{\left\{{\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}\right\}}^{\mathrm{n}}\cdot {\mathrm{w}}_2.$ (23)

Considering the approximations in each element

${\left\{{\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}\right\}}_{\mathrm{h}}^{\mathrm{e}}={\mathrm{N}}_{\mathrm{g}}{\left\{\widehat{{\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}}\right\}}^{\mathrm{e}},\mathrm{}\nabla {\left\{{\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}\right\}}_{\mathrm{h}}^{\mathrm{e}}={\mathrm{B}}_{\mathrm{g}}{\left\{\widehat{{\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}}\right\}}^{\mathrm{e}},\mathrm{}{\mathrm{w}}_2^{\mathrm{e}}={\mathrm{N}}_{\mathrm{g}}{\widehat{\mathrm{w}}}_2^{\mathrm{e}},\mathrm{}{\nabla \mathrm{w}}_2^{\mathrm{e}}={\mathrm{B}}_{\mathrm{g}}{\widehat{\mathrm{w}}}_2^{\mathrm{e}}\mathrm{},$ (24)

we can define the operators at the element level as

${\mathbb{M}}_{\mathrm{g}}^{\mathrm{e}}={\int }_{\mathrm{\Omega }}{{\mathrm{N}}_{\mathrm{g}}}^{\mathrm{T}}{\mathrm{N}}_{\mathrm{g}},\mathrm{}{\mathbb{Q}}_{\mathrm{g}}^{\mathrm{e}}={\int }_{\mathrm{\Omega }}{{\mathrm{N}}_{\mathrm{g}}}^{\mathrm{T}}\mathrm{M}{\mathrm{N}}_{\mathrm{g}},{\mathbb{K}}_{\mathrm{g}}^{\mathrm{e}}={\int }_{\mathrm{\Omega }}{{\mathrm{B}}_{\mathrm{g}}}^{\mathrm{T}}{\mathrm{B}}_{\mathrm{g}},{\mathrm{w}}_{\mathrm{b}}^{\mathrm{e}}=\widehat{\mathrm{F}}\frac{{\mathrm{g}}_{\mathrm{c}}}{\mathrm{\gamma }}{\int }_{\mathrm{\Omega }}{{\mathrm{N}}_{\mathrm{g}}}^{\mathrm{T}}\mathrm{d}.$ (25)

Using the standard assembling procedure, the equation of degradation to be implemented and solved numerically is

$\left[{\mathbb{M}}_{\mathrm{g}}+\mathrm{\Delta }\mathrm{t}\widehat{\mathrm{F}}\left({\mathbb{Q}}_{\mathrm{g}}+\frac{{\mathrm{g}}_{\mathrm{c}}}{\mathrm{\gamma }}{\mathbb{M}}_{\mathrm{g}}+{\mathrm{g}}_{\mathrm{c}}\mathrm{\gamma }{\mathbb{K}}_{\mathrm{g}}\right)\right]{\left\{\widehat{{\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}}\right\}}^{\mathrm{n}+1}={\mathbb{M}}_{\mathrm{g}}{\left\{\widehat{{\mathbb{G}}^{\mathrm{p}\mathrm{s}\mathrm{}}}\right\}}^{\mathrm{n}}+\mathrm{\Delta }\mathrm{t}{\mathrm{w}}_{\mathrm{b}}.$ (26)

4 RESULTS

In this section, we present some numerical examples aiming to evaluate the evolution of the material degradation and its relation to the material constants. Furthermore, the model ability to qualitatively represent the material response under different load conditions is also evaluated. In the following examples, we use $\widehat{b}=0$ for all analyses. This term is related to the elastic damping due to the contribution of the velocity gradient and its effect is not relevant for the stability of the implicit time integration methods considered. The other phase field parameters have been adjusted to qualitatively match experimental and numerical results from the literature.

4.1 I-shaped tensile test with unloading

The first numerical example is inspired on the experimental procedure used to indirectly measure the material damage. After loading the specimen until a certain damage level, it is unloaded and the inclination of the stress-strain diagram is measured to obtain the Young’s modulus ($\mathrm{E}$) of the damaged material. The Poisson’s ratio can also be obtained by its definition

$\mathrm{\nu }=-\frac{{\mathrm{E}}_{11}}{{\mathrm{E}}_{22}}$ .(27)

The Young’s modulus $\mathrm{E}$ and the Poisson’s ratio $\nu$ can then be used to compute the damage variable depending on the theoretical model considered (Lemaitre and Dufailly, 1987). For the present simulation, we used a rectangular tension test specimen from ASTM Standard E8/E8M-13a (ASTM Standard, 2013), shown in Fig. 1. The elasticity modulus will be evaluated in only one direction.

Figure 1
– Dimension of the specimen (in millimeters), boundary conditions and finite element mesh of the I-shaped tensile test.

First, in the damaging phase, a displacement was prescribed at the right edge of the specimen and linearly increased by the rate of 10⁻³ m/s until the degradation matrix component ${\mathrm{G}}_{11}$ reached the value of 0.7 as show in Figure 2. After that, the specimen was unloaded at the same displacement rate. It is worth mentioning that ${\mathrm{G}}_{11}$ was used as a parameter for measuring the damage state of the material because for a tensile test of a specimen oriented along the x-direction, we can associate the degradation process with the component ${\mathrm{G}}_{11}$ as discussed in Petrini (2023). For a more complex strain state or internal load not aligned to one of the global axes, different components of the degradation tensor can have a more significative influence in the fracture behavior.

Figure 2
– Contour plot of the component ${\mathrm{G}}_{11}$ after the damaging phase.

The domain was discretized using 4066 triangular linear elements and the mesh is shown in Fig. 1. The initial material properties used are Young’s modulus $\mathrm{E}$ = 180 GPa, Poisson’s ratio $\nu$ = 0.3, density $\rho$ = 7300 kg/m³ and Griffith fracture energy ${\mathrm{g}}_{\mathrm{c}}$ = 2700 N/m. The model parameters for this simulation are $\gamma$ = 2.5×10⁻⁴ m, $\widehat{\mathrm{F}}$ = 1. Triangular linear elements were used for this analysis to validate the model with an element different from the commonly used quadrilateral linear element. A time-step of 10^-2 s was chosen over 64 steps.

The stress-strain diagram is shown in Fig. 3. It was obtained by plotting the stress and strain components in x-direction (${\mathrm{T}}_{11}$ and ${\mathrm{E}}_{11}$) of a point in the geometric center of the specimen. The change in the slope between the load and unloading paths evidences the reduction in the stiffness of the material after damage. Note that the unloading path is straight as no damage occurs. As expected, there is no residual strain since we are not considering the effects of plasticity.

Figure 3
– Stress-strain diagram of the I-shaped tensile test with unloading.

By dividing the stress component ${\mathrm{T}}_{11}$ by the strain component ${\mathrm{E}}_{11}$ for each time-step we obtain the actual Young’s modulus as depicted in Fig. 4(a). We can see that $\mathrm{E}$ changes non-linearly from its original value (180 GPa) until the value of 859 MPa at the end of the loading phase. As expected, in the unloading phase the Young’s modulus $\mathrm{E}$ remains constant. The Poisson’s ratio has a similar behavior as shown in Fig. 4(b) and the final value is 0.289. However, these values are not sufficient to calculate the elasticity matrix of the material. Due to the tensor nature of the damage model, the material becomes locally anisotropic and therefore we need the material parameters related to the y-direction to determine the material matrix for plane stress state.

Figure 4
– Time evolution of (a) Young-modulus $\mathrm{E}$ and (b) Poisson’s ratio.

From the results presented above, we may observe the effect of the degradation on the material properties analyzed in the x-direction for a tension test.

4.2 L-shaped panel static test

The L-shaped specimen test is a method used to determine the fracture toughness of concrete. The test involves creating an L-shaped specimen from a concrete mixture and submitting it to a concentrated load at its tip as illustrated in Fig. 5(a). The crack propagation is then monitored, and the fracture toughness is calculated based on the load at which the crack propagates a certain distance. This test has become a popular benchmark test for the validating damage models for the numerical simulation of plain concrete cracking as pointed out by Arruda (2011).

Figure 5
– (a) Geometry and boundary conditions of the L-shaped panel example (in millimeters) (b) Finite element mesh.

The geometry and boundary conditions used in the simulation are depicted in Fig. 5(a). An incremental displacement rate of 0.01mm/s was considered. The domain was discretized with 11663 quadrilateral linear elements and the mesh is shown in Fig. 5(b). The material properties used for concrete are Young’s modulus $\mathrm{E}$ = 25.985 GPa, Poisson’s ratio $\nu$= 0.18, density, $\mathrm{\rho }$ = 2400 kg/m³ and Griffith fracture energy ${\mathrm{g}}_{\mathrm{c}}$ = 100 N/m. The model parameters for this simulation are $\mathrm{\gamma }$= 1.0×10⁻⁶ m, $\widehat{\mathrm{F}}$ = 1. A time-step of 10^-2 s was chosen for this example.

Figure 6 presents the normalized geometric mean of the elasticity matrix eigenvalues, which gives the crack path over the L-shaped specimen considering the tensor damage model at time t = 9.00 s.

Figure 6
– Numerical result obtained with the tensor damage model.

The crack pattern obtained for the tensor damage model is in agreement with the experimental result presented by Arruda (2011) in terms of the location and extent of the damage in the panel.

4.3 Plate with elongated hole

The next numerical example is the tensile test of a plate with an elongated hole. The geometry and boundary conditions, illustrated in Fig. 7(a), are the same as considered by Fassin et al. (2019) where a gradient-extended continuum mechanical framework with a second order damage tensor is used to model anisotropic brittle damage. Only one fourth of the specimen was simulated due to symmetry. Therefore, the displacement in the x-direction was constrained on the right edge and in the y-direction on bottom edge of the specimen. An incremental displacement rate of 10⁻⁴ mm/s was applied on the top edge.

Figure 7
– (a) Geometry and boundary conditions of the plate with elongated hole example (in millimeters). Adapted from: Fassin et al. (2019) (b) Finite element mesh.

The domain was discretized with 26266 quadrilateral linear elements as shown in Fig. 7(b). The material properties used were Young’s modulus $\mathrm{E}$ = 200 GPa, Poisson’s ratio $\mathrm{\nu }$ = 0.24, density $\mathrm{\rho }$ = 7800 kg/m³ and Griffith fracture energy ${\mathrm{g}}_{\mathrm{c}}$ = 2700 N/m. The model parameters for this simulation are $\mathrm{\gamma }$ = 10⁻⁶ m, $\widehat{\mathrm{F}}$ = 1 and time-step of ${10}^{-3}$s. Figure 8 illustrates the crack evolution represented by the normalized geometric mean of the elasticity matrix eigenvalues (ith the tensor damage model. The crack pattern obtained for the tensor damage model is in agreement with the experimental result presented by Arruda (2011) in terms of the location and extent of the damage in the panel. 4.3 Plate with elongated hole The next numerical example is the tensile test of a plate with an elongated hole. The geometry and boundary conditions, illustrated in Fig. 7(a), are the same as considered by Fassin et al. (2019) where a gradient-extended continuum mechanical framework with a second order damage tensor is used to model anisotropic brittle damage. Only one fourth of the specimen was simulated due to symmetry. Therefore, the displacement in the x-direction was constrained on the right edge and in the y-direction on bottom edge of the specimen. An incremental displacement rate of 10−4 mm/s was applied on the top edge. Figure 7 – (a) Geometry and boundary conditions of the plate with elongated hole example (in millimeters). Adapted from: Fassin et al. (2019) (b) Finite element mesh. The domain was discretized with 26266 quadrilateral linear elements as shown in Fig. 7(b). The material properties used were Young’s modulus E = 200 GPa, Poisson’s ratio ν = 0.24, density ρ = 7800 kg/m3 and Griffith fracture energy g_c = 2700 N/m. The model parameters for this simulation are γ = 10−6 m, F ̂ = 1 and time-step of 10^(-3)s. Figure 8 illustrates the crack evolution represented by the normalized geometric mean of the elasticity matrix eigenvalues ( ).

Figure 8
– Contour plot of the normalized geometric mean of elasticity matrix eigenvalues at three different time steps.

The results show that the crack nucleates on the edge of the elongated hole near the transition between the straight and the circular sections. The crack then grows in the x-axis towards the external side of the specimen. The crack width is smaller close to the nucleation and becomes thicker until it approaches the edges of the specimen, where we can observe a crack bifurcation. Diffuse damage was present in the lateral of the specimen.

The location of crack nucleation and the growth direction are in agreement with the results reported by Fassin et al. (2019). However, they obtained a crack with approximately constant width and observed different crack propagation angles for different damage thresholds.

4.4 Cross specimen

The anisotropy induced by the damaging process is particularly important in situations where sequential loads are applied. As showed by Fassin et al. (2019), isotropic and anisotropic damage models can lead to qualitatively different results. Two loading steps were applied in a cross specimen composed by two different materials (see Figure 9). First, a pre-damaging step was applied through a tensile load on the upper edge of the specimen. After that, the specimen is unloaded and a second load is applied in the perpendicular direction until the crack propagates completely, see Fassin et al. (2019).

Figure 9
– (a) Cross specimen dimensions (in millimeters). (b) Finite element mesh.

Both isotropic and anisotropic gradient-extended brittle damage models were considered by Fassin et al. (2019), with the latter incorporating a second-order damage tensor to model anisotropy. In the present work, the same specimen geometry and material constants ($\mathrm{E}$ and $\mathrm{í}$) were used to simulate the same example with damage process given by the tensor damage model.

In this example, the crack is assumed to propagate only in the inner circle of the specimen. The outer material has a larger Young’s modulus (${\mathrm{E}}_1$) and is not subjected to damage. Therefore, its critical Griffith energy is chosen to be very high. The inner material is assumed to be softer (${\mathrm{E}}_2=0.1{\mathrm{E}}_1$) and to define its critical Griffith energy ${\mathrm{g}}_{\mathrm{c}}$ and the model constants, we compared the obtained damage state at ${\mathrm{t}}_0=0.5$ s with that presented by Fassin et al. (2019).

The boundary conditions are shown in Figure 10. Only one fourth of the specimen is simulated. In the pre-damaging step ($\mathrm{t}\mathrm{}\le 0.5$), a displacement rate of 1.28×10⁻³ m/s is applied on the top edge of the specimen. After that, the specimen is totally unloaded in the same rate. Then, the displacement is prescribed on the left edge at the same rate of 1.28×10⁻³ m/s as shown in Figure 10. A mesh of 17410 quadrilateral elements was used as illustrated in Figure 9(b). The material properties and the model parameters for this simulation are presented in Table 1. A time-step of ${10}^{-4}$ s was used.

Figure 10
– Boundary conditions for the cross specimen simulation.

Table 1
– Material parameters and model constants.

Material constants
	Outer Material	Inner material
$E$	100000 MPa	10000 MPa
$\nu$	0.3	0.3
$g_c$	1000 N/m	100 N/m
$\rho$	2400 Kg/m3	2400 Kg/m3
Model constants
$\widehat{F}$	1000
$\mathrm{\gamma }$	10−6 m

Figures 11 to 13 illustrate the contour plots of the normalized geometric mean of elasticity matrix eigenvalues $\mathrm{\ae }$ and the components ${\mathrm{G}}_{11}$ and ${\mathrm{G}}_{22}$ at different time steps, where the first image represents the contour plot at the end of the pre-damaging step.

Figure 11
– Evolution of the normalized geometric mean of elasticity matrix eigenvalues $\mathrm{\ae }$ in the final crack propagation step.

Figure 13
– Evolution of degradation matrix ${\mathrm{G}}_{22}$ in the final crack propagation step.

From Figure 11, we can observe that the crack nucleates near the interface of the materials and starts to propagate from both sides of the crack until they meet at the horizontal line of symmetry. Comparing Figures 12 and 13 with Figure 11 we observe that, for the crack path obtained, the largest contribution for the crack propagation comes from the component ${\mathrm{G}}_{11}$. Moreover, during the application of the load in x-direction, the component ${\mathrm{G}}_{22}$ remains almost constant.

Figure 12
– Evolution of degradation matrix component ${\mathrm{G}}_{11}$ in the final crack propagation step.

The crack path obtained is similar to the one reported by Fassin et al. (2019) for the anisotropic damage model except by the fact that they obtained a crack in the upper and bottom parts of the inner circle, which did not happen for the final time step considered in the simulation of the present work. This difference may be a consequence of the chosen energy split for tension-compression asymmetry. The missing pattern occurs in a region of the domain that is under compression which, according to our model, does not propagate cracks. Further adjustments could be needed in the model to match with the original results.

5 CONCLUSION

The present work aimed to contribute to the development of a general anisotropic damage model that eliminates the need of knowing the principal damage directions a priori. The formulation is based on a thermodynamically consistent non-isothermal damage phase field framework. An internal fourth-order variable, called degradation tensor, is introduced to represent the damage effect on the material properties.

The results have shown that the new local material symmetry resulting from the degradation process was totally driven by the strain state obtained from the applied load. The definition of the degradation tensor in the global axes was successful in simulating anisotropic damage in different numerical examples.

The tensile test results show the gradual degradation of the material properties both indirectly through the stress-strain relationship and directly by each individual property over time. This is an interesting result as the degradation tensor transforms the elasticity tensor rather than having an explicit degradation function for them.

The model was also capable of tracking straight and curved cracks for different sets of boundary conditions and symmetries. For the elongated hole, the rate of the applied load may have had influence on the curvature of the crack as we consider inertia effects.

Lastly, the irreversibility criterion herein adopted was shown to properly prevent the healing of the degraded material. Although the results were different in compressive regions for the cross example, there were very good similarities in the tractive regions among the degradation components and the damage variables of the original paper.

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 140214/2022-4.

• https://doi.org/10.1590/1679-78257589

References

• ASTM Standard E8/E8M-13a, ”Standard Test Methods for Tension Testing of Metallic Materials”. ASTM International, i:1–27, 2013.
• J. L. Boldrini, E. A. B. D. Moraes, L. R. Chiarelli, F. G. Fumes, M. L. Bittencourt, E. A. Barros de Moraes, L. R. Chiarelli, F. G. Fumes, and M. L. Bittencourt. A non-isothermal thermodynamically consistent phase field framework for structural damage and fatigue. Computer Methods in Applied Mechanics and Engineering, 312:395–427, 2016.
• M. Brunig. An anisotropic continuum damage model: Theory and numerical analyses. Latin American Journal of Solids and Structures, 1(2):185–218, 2004.
• J. L. Chaboche. Continuous damage mechanics - A tool to describe phenomena before crack initiation. Nuclear Engineering and Design, 64(2):233–247, 1981.
• C. L. Chow and J. Wang. An anisotropic theory of continuum damage mechanics for ductile fracture. Engineering Fracture Mechanics, 27(5):547–558, 1987.
• J. D. Clayton and J. Knap. A geometrically nonlinear phase field theory of brittle fracture. International Journal of Fracture, 189(2):139–148, 2014.
• J. D. Clayton and D. L. McDowell. Finite polycrystalline elastoplasticity and damage: Multiscale kinematics. International Journal of Solids and Structures, 40(21):5669–5688, 2003.
• M. Fassin, R. Eggersmann, S. Wulfinghoff, and S. Reese. 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, 2019.
• J. Jaric, D. Kuzmanovi ́c, and D. Sumarac. On anisotropic elasticity damage mechanics. International Journal of Damage Mechanics, 22(7):1023–1038, 2013. ISSN 10567895.
• J. Lemaitre and J. Dufailly. Damage measurements. Engineering Fracture Mechanics, 28(5-6):643–661, 1987.
• N. Mozaffari and G. Z. Voyiadjis. Phase field based nonlocal anisotropic damage mechanics model. Physica D: Nonlinear Phenomena, 308:11–25, 2015.
• J. J. d. C. Pituba. Sobre a Formulação De Um Modelo De Dano para o Concreto. PhD thesis, Universidade de São Paulo, 2003.
• J. J. Skrzypek and Artur W. Ganczarski. Mechanics of Anisotropic Materials. Springer, Ganczarski, 2015. ISBN 978-3-319-17159-3.
• A. L. E. R. Petrini, A thermodynamically consistent phase field framework for anisotropic damage, Ph.D. thesis, Universidade Estadual de Campinas, 2023.
• M. R. T. Arruda. Static and dynamic analysis of concrete structures using damage mechanics. Tese, PhD thesis, Universidade de Lisboa, 2011.
• C. Miehe, F. Welschinger, M. Hofacker. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations, International Journal for Numerical Methods in Engineering, 2010.
• G. A. Haveroth, E. A. B. de Moraes, J. L. Boldrini and M. L. Bittencourt. Comparison of semi and fully-implicit time integration schemes applied to a damage and fatigue phase field mode. Latin American Journal of Solids and Structures, 2018.
• A. L. E. R. Petrini, J. L. Boldrini and M. L. Bittencourt. A Thermodynamically Consistent Phase Field Framework for Anisotropic Damage Propagation. Latin American Journal of Solids and Structures, 2020.
• M. G. Vale, J. A. Avila, J. L. Boldrini and M. L. Bittencourt. Efficient crack length measurement using A* shortest path methodology for a phase-field fracture framework. Latin American Journal of Solids and Structures, 2023.
• A. G. Neeman, R. Brannon, B. Jeremic, A. V. Gelder and A. Pang. Decomposition and visulatization of fourth-order elastic-plastic tensors. IEEE/EG Symposium on Volume and Point-Based Graphics, 2008.
• M. M. Mehrabadi, S. C. Cowin. Eigentensors of linear anisotropic elastic materials. The quaterly Journal of Mechanics and Applied Mathematics, 1991.

Edited by

Editor: Marco L. Bittencourt and Josué Labaki

Publication Dates

• Publication in this collection
  17 Nov 2023
• Date of issue
  2023

History

• Received
  21 Mar 2023
• Reviewed
  02 July 2023
• Accepted
  01 Aug 2023
• Published
  21 Aug 2023