A Phase-Field Approach to Continuum Damage Mechanics

This paper develops a phase-field approach to describe the damage within continuum mechanics. The body is associated with the standard stresses and body forces of macroscopic character. As is the case in many contexts, the phase field is a scalar variable whose time rate is governed by a constitutive equation. The generality of the approach allows the modeling of non-stationary heat conduction, mechanical hysteretic effects, and the macroscopic damage. The thermodynamic consistency is investigated through the constraint of the Clausius–Duhem inequality following the standard procedure of Rational Thermodynamics. The entropy production is considered as a constitutive function; this view was proved to be essential in the elaboration of hysteretic models. Here, the scheme is improved by viewing the entropy production as a sum of two terms, one induced by the other constitutive equations and one given by a constitutive equation of the entropy production per se.


Introduction
This paper is devoted to a thermodynamically consistent description of damage evolution in continuum mechanics. Damage in continuum mechanics is of interest in many respects. It is a way of characterizing the aging of materials, and it is also a measure used to describe the nucleation of new cracks in the absence of macroscopic initial cracks.
Perhaps motivated by the search of models more akin to the microscopic setting of damage, some approaches have been developed where micro stresses and micro forces are framed along with the standard stresses and forces of macroscopic character (see, e.g., [1][2][3][4][5]). Though this idea leads to a more detailed scheme, it involves additional unknown fields to be determined in an operative way. To make the approach simpler, micro kinetic terms are neglected and hence, micro forces are subject to equilibrium conditions. The purpose of this paper is to set up a phase-field approach to describe the damage within continuum mechanics. No microscopic fields are considered; the body is associated with the standard stresses and body forces of macroscopic character. As is the case in many contexts, the phase field (or order parameter or internal variable) is a scalar variable whose time rate is governed by a constitutive equation [6]. The generality of the present approach allows the modeling of non-stationary heat conduction, mechanical hysteretic effects, and macroscopic damage.
The thermodynamic consistency is investigated through the constraint of the Clausius-Duhem inequality following the standard procedure of Rational Thermodynamics [7]. Lately, a more general scheme has been applied in that the entropy production is considered as a constitutive function; this view has proved essential in the elaboration of hysteretic models [8][9][10]. Here, the scheme is improved by viewing the entropy production as a sum of two terms, one induced by the other constitutive equations and one given by a constitutive equation of the entropy production per se (extended entropy production). We will see that the occurrence of the extended entropy production is essential in connection with the damage variable with causes such as, e.g., excesses of temperature, strain, and stress.
Recent developments in material modeling show a clear distinction between different causes of damage. For instance, the displacement field is decomposed into translation, rotation, plastic stretches, elastic stretches, and volumetric and shear stretches [11,12]. The present approach is consistent with such more refined descriptions provided only that the entropy production, as any constitutive function, allows for the pertinent dependencies. Notation 1. We consider a solid occupying a time-dependent region Ω ⊂ E 3 . Throughout, v is the velocity, ∇ denotes the gradient operator, ∂ t is the partial time derivative at a point x ∈ Ω, while a superposed dot stands for the total time derivative,ḟ = ∂ t f + v · ∇ f . Cartesian coordinates are used: L is the velocity gradient, L ij = ∂ x j v i , and D = symL is the stretching tensor. We let R be a reference configuration; the motion χ is a function that maps each point vector X ∈ R into a point x ∈ Ω. The deformation gradient F is defined by F iK = ∂ X K χ i and J = det F > 0. The other mathematical characters are defined at the first stage of usage.

Balance Equations and Admissible Processes
Let ϕ be the damage variable. According to the literature, there are several measures associated with the scalar ϕ. For instance, ϕ may be the fraction of damaged area [13] or the volumetric fraction of damaged material [5]. Anyway ϕ is a scalar with values in [0, 1]; for definiteness, ϕ = 0 represents an undamaged material, ϕ = 1 a fully damaged material. There are cases in which the damage is anisotropic as it happens in elastostatics [12]. The scalar character of the damage variable might be maintained by considering a set of degradation functions, as seen in [14].
The body is modeled as a material with an internal variable (or phase field), the scalar ϕ. Hence, the balance equations are those of a simple body [15]. In the local form, the continuity equation, the balance of linear momentum, and the balance of energy are written asρ where ρ is the mass density, v the velocity, T the symmetric stress tensor, b the body force (per unit mass), ε the internal energy density, q the heat flux, and r the external energy supply (per unit mass). Let η be the entropy density (per unit mass). The balance of entropy is assumed in the form where j is the entropy flux; we let j = q/θ + k with k the extra-entropy flux to be determined. In general, the entropy flux j need not equal q/θ, and is given by a constitutive equation. This view traces back to I. Müller [16] (see [17] for a detailed exposition). We define the entropy production σ as and, by the balance of entropy, σ ≥ 0. Furthermore, we let σ be given by a constitutive equation, per se, as is the case for η and j. We let the set of fields ρ, v, T, b, ε, r, q and η, j (or k), and σ, subject to (1), be an admissible process. We take as the second law of thermodynamics the following statement. Second law of thermodynamics. For every admissible process, the inequality (2) must hold for all times t ∈ R and points x ∈ Ω.
Using the Helmholtz free energy ψ = ε − θη, we obtain the inequality in the form Additionally, for a connection with the literature, we mention other approaches to the modeling of damage.

Other Balance Equations in the Literature
Motivated by a microscopic picture of damage, additional power terms are associated with ϕ viaφ and ∇φ. By appealing to the principle of virtual power [18] the damage variable ϕ is taken to occur in the internal and external virtual powers P i , P e . With reference to [19], we take Let ϕ → ϕ +φ, v → v +ṽ. Hence, the requirement The microscopic fields α, β, h, γ are then subject to the local balance Equations (4). This in turn means that appropriate constitutive equations are needed.
The power βφ + h · ∇φ would also occur in the entropy inequality [5]. We now go back to the present approach, which is free from fields of microscopic character.

Constitutive Equations and Thermodynamic Restrictions
The interest in constitutive dependencies on ϕ,φ, and ∇φ indicates that the Lagrangian description is more convenient than the Eulerian one. Hence, we consider the referential quantities [15] where T RR is the second Piola (or Piola-Kirchhoff) stress tensor, while for any function f (x). The multiplication of (3) by J = det F then results in the Lagrangian version of the second-law inequality, where E = 1 2 (F T F − 1) is the Green-St. Venant strain tensor. The use of the referential quantities T RR , E, and q R is mathematically advantageous whenever we describe the material behavior through rate equations, in thatṪ RR ,Ė, andq R are objective, whereasṪ anḋ q are not. Of course (3) and (5) are equivalent. Though we use E, rather than C = F T , F, to describe the strain, the present approach applies to finite deformation in that no linear approximation is considered. Let be the set of independent variables. Computeψ and substitute in (5) to obtain The linearity and arbitrariness ofË,T RR ,q R ,θ,φ, ∇ Rθ , andθ imply that Hence, the free energy is a function of a reduced number of variables, namely subject to the standard relation (6) for the entropy. For definiteness, we now examine further thermodynamic restrictions by specifying the type of continuum we have in mind. Divide the remaining inequality by θ to get In light of the identity we can write inequality (7) in the form is the variational derivative of ψ with respect to ϕ. This suggests that we let whereφ is yet to be determined. For definiteness, we take ∇ R θ andq R as independent ofĖ,Ṫ RR ,φ and let where σ| · denotes the appropriate restriction of the function σ. Hence, we can split the inequality to get

Further Thermodynamic Restrictions
The damaged material has properties affected by the level of damage, formally characterized by ϕ. We then look for a modeling with the free energy in the form We begin with Equation (9), which governs heat conduction.

Heat Conduction
The function ψ q is a scalar function of q R . We then let ψ q depend on q R through ξ = |q R | 2 . Hence, Equation (9) can be written in the form where σ q is independent ofĖ,Ṫ RR ,φ. The requirement (10) implies that q R ,q R , and ∇ R θ cannot be independent. Among the relations consistent with (10), we select so that in stationary conditions, we have Fourier-like laws We then obtain, respectively, where κ is allowed to depend on θ and ϕ. The non-negative character of σ q implies that κ > 0.
To obtain the corresponding Eulerian version, we observe that where q is the Truesdell derivative [20,21] q=q − Wq − Dq + (∇ · v)q.
As to (11) 2 we have Left multiplication of (13) by FC −1 results in In both cases, we can view τ = 2θρ R κ∂ ξ ψ q as the relaxation time. Both Equations (11) can be viewed as nonlinear generalizations of the Maxwell-Cattaneo equation [22,23]. The occurrence of the left Cauchy-Green tensor B in (12) and (14) shows possible effects of deformation on the heat conduction. In both cases, in stationary conditions (q= 0), we have the classical Fourier law, q = −κ∇θ, or the modification due to deformation, q = −κB∇θ, whereκ = κ/J. There are infinitely many free energies ψ q consistent with this scheme. If, for simplicity, we let We now go back to Equation (8) and assumeφ is independent ofĖ andṪ RR so that we can write the independent equations

Hysteretic Mechanical Effects
Some classes of materials models described by (16) are now investigated. First, we consider the particular case σ visco = 0 and ∂ T RR ψ el = 0. The resulting equation is The arbitrariness ofĖ implies that This relation models an elastic solid parameterized by the temperature θ and the damage variable ϕ.
If instead σ visco = 0, but ∂ T RR ψ el = 0, then we have Consequently,Ṫ RR can be given a linear representation inĖ. Indeed, we use a representation formula of tensors [10]; for any tensor Z with a known value of the inner product Z · N, and N · n = 1, we have where I is the fourth-order identity and G is an arbitrary second-order tensor. Here, we let Z =Ṫ RR and N = ρ R ∂ T RR ψ el /|ρ R ∂ T RR ψ el | to obtaiṅ Observe that by letting G = G RR (θ, ϕ, E, T RR )Ė, we can writė where The representation (18) describes the constitutive properties of hypoelastic solids. Elastic-plastic models are characterized by an entropy production-here, σ viscowhich depends linearly on |Ė| or |Ṫ RR |. Back to (16); let Jθσ visco = γ E |Ė| to obtain Assume ∂ T RR ψ el = 0. Hence, we can writė and applying the representation formula (17), we obtaiṅ The analogue holds if we let σ visco = γ T |Ṫ RR |.
As is shown in refs [8][9][10], in one-dimensional settings, the simultaneous occurrence ofĖ and |Ė| allows the modeling of hysteretic processes. Here, we have proved that the structure (19) follows directly from the entropy inequality by simply letting the positive quantity σ equal γ E |Ė| or γ T |Ṫ RR |. Moreover, this value of σ is not identically equal to the left-hand side of the entropy inequality, as instead it happens for the model of heat conduction, where σ q is determined by the left-hand side of (11). This conceptual aspect will be more refined, in the next section, in connection with the modeling of damage.
We expect that the thermoelastic properties decay with the increase in the damage variable ϕ, and henceĠ This qualitative conclusion is consistent with the observation that, by (20), sgn{G φ} = sgn{−σ dam }, not necessarily but consistently. With this in mind, we can view G as a known function subject to the monotone character G ≤ 0. Our attention is then restricted to ψ dam , a function of the values ϕ and ∇ R ϕ; the dependence on ∇ R ϕ represents the possible effects of spatial inhomogeneities.
We now consider Equation (15) in the form The left-hand side is non-negative iḟ Equation (21) allows for a further contribution toφ so thaṫ This is so that which shows that the rate Equation (23) yields an entropy production Look at the two effects separately. First, letσ dam = 0. Additionally, with reference to the literature (see e.g., [4,5] and refs therein) we may take λ as a constant, possibly related to parameters of the material (here, θ, E, T RR , q R ), and Hence, the evolution Equation (22) readṡ If, in particular, ψ 2 (|∇ R ϕ| 2 ) = 1 2 c|∇ R ϕ| 2 and ρ, θ are space independent, then we havė Incidentally, in these cases, and hence the requirement σ dam ≥ 0 holds for every function ψ dam . We now look for specific forms ofσ dam . Suppose that an increase of ϕ is due to high temperatures, θ > θ h , freeze-thaw transitions at θ = θ 0 , and large strains |E| > κ. Large cycles are also of interest. Letting E = Ee ⊗ e, we have damage effects if E > E + > 0 or E < E − < 0 for suitable values E − , E + . Effects of large values are modeled through terms proportional to H being the Heaviside step function. The whole effect onσ dam can then be written via appropriate coefficients in the form

History Effects on Damage
Damage may be the effect of phase transformations, as, e.g., in the solid-liquid case, or cyclic processes, as e.g., in periodic or hysteretic evolutions. This indicates that the present value of ϕ depends on the history of appropriate fields. The idea is not new (see, e.g., [5]). The subject deserves a detailed treatment both for conceptual questions associated with the thermodynamic consistency and for the development of definite models.
We follow a Lagrangian description and suppose the constitutive properties are expressed by the set of variables where θ t , E t denote the histories of θ, E up to time t, We then let ψ, η, k, σ, and the rateφ be functionals of Γ. Upon computingψ we replace in the Clausius-Duhem inequality (5) to obtain where dψ(θ t |θ t ) and dψ(E t |Ė t ) denote the Fréchet derivatives of ψ at θ t and E t in the directionθ t andĖ t . The reduced inequality is Observe that allows us to describe hysteretic effects, models heat conduction, and −ρ R ∂ ϕ ψφ indicates the evolution of damage. A sufficient condition to satisfy (25) is to let σ hys , σ cond , σ dam being non-negative in that they are particular cases of σ. Hence, in addition to the equations for σ hys , σ cond , we have For definiteness, we look for specific forms of σ dam . As observed above, large values effects are modeled through terms proportional to A model for the effect of a freeze-thaw transition is made formal by letting δ > 0 and taking a contribution, at time t, of the form which works for freeze-thaw and for thaw-freeze. The whole change in ϕ as t ∈ [t 1 , t 2 ] is then expressed by A model for the effect of a freeze-thaw transition is made formal by letting δ > 0 and taking a contribution, at time t, of the form which works for freeze-thaw and for thaw-freeze.
For the sake of simplicity, sometimes the damage variable is evaluated by a cumulative assessment (Miner's rule [24,25]). So, if the damage is produced by a number of stress cycles, at a given stress level, and N is the number of cycles producing failure, then ϕ associated with n ≤ N cycles is determined by Incidentally, for each freeze-thaw transition at timet, we have the number of freeze-thaw transitions as t ∈ [t 1 , t 2 ]. Once we fix the derivative ∂ ϕ ψ and the coefficients c h , c E , c f t , c + , c − , we obtain the variation of the damage variable ϕ in [t 1 , t 2 ].

Remarks on the Phase-Field Models
The evolution Equation (24) generalizes known model equations appeared in the literature. The rateφ consists of two parts, −λδ ϕ Ψ, − θ ρδ ϕ Ψσ dam .
The first term involves the present value of ϕ, possibly through the Laplacian ∆ R ϕ. With reference to the review paper [26], as particular cases we mention the models by Karma et al., Henry and Levine, Ambati et al., and Miehe et al., whereφ is affected by ϕ, ∆ϕ, and is governed by the strain.
The second term allows specific causes of damage. Equation (24), and the analogue (26) for modeling through histories, indicates how the model equation ofφ may account for high temperatures, freeze-thaw transitions, large strains, and cycles. The possible dependence of the coefficients c h , c E , c + , c − on ϕ itself and the temperature θ allows an effective modeling of damage effects. Experimental evidence might give insights into the constitutive dependencies of the coefficients.

Relation to Other Approaches
A number of papers involve the use of microforces in the balance equations. Nevertheless, some similarities appear. In [3], the free energy is assumed in the form which shows the correspondencẽ G and as parameters. The function describes the double well potential. This potential is widely used in the modeling of phase transitions. To our mind, in the modeling of damage, a monotone dependence would be more appropriate; the monotone character is advocated, e.g., in [27]. As in [5], the free energy is sometimes taken to depend on the history of E t and, moreover, the free energy potential is considered in terms of a fractional derivative. The thermodynamic consistency requires that the derivative dψ(E t |Ė t ) is among the contributions to the non-negative entropy production. Now, while a free energy potential can be written for the stress tensor, problems arise when we investigate the thermodynamic consistency, mainly because the kernel of fractional derivatives is singular at the origin [28].

Conclusions
A model for the characterization of damage in continuum mechanics is developed for a dissipative and heat-conducting solid. The damage is described by a scalar variable ϕ (phase field) whose evolution is governed by a rate equation consistent with thermodynamics. The consistency is investigated within a differential equation of the form (23), where Ψ is the free energy density. The term λδ ϕ Ψ is similar to models in the literature. The second term allows an account of effects such as those arising from large temperatures, large strains, and freeze-thaw transitions. Conceptually, the two terms have a different origin. From δ ϕ Ψφ ≤ 0, we conclude thatφ may be given by −λδ ϕ Ψ, λ > 0. This in turn implies that the corresponding entropy production is ρ/θ times λ(δ ϕ Ψ) 2 . This value of entropy production is induced by the constitutive function Ψ. The second term containsσ, which has a constitutive equation per se, subject to the positive valuedness and the dependence of variables chosen for all of the constitutive equations. This scheme is likely to allow some improvements of the modeling of damage. In this sense, appropriate generalizations can account for more involved effects through the entropy productionσ dam . Indeed, the use ofσ dam looks more flexible than the recourse to the dissipation potential. For instance, in [29] the damage rate, here,φ, is related to a dissipation potential F byφ = −∂ Y F, where Y is an appropriate density release rate. Hence, terms such as those in (24) and (26) account for the rate as (nine) of [29] for the accumulated plastic strain.
Future developments might investigate the use of degradation functions induced by the damage variable so as to describe, e.g., possible anisotropic effects. Data Availability Statement: The study did not report any data.