Original article
Comparison of semi and fullyimplicit time integration schemes applied to a damage and fatigue phase field model
Geovane A. Haveroth^{a
}
Eduardo A. Barros de Moraes^{b
}
José L. Boldrini^{a
}
Marco L. Bittencourt^{b
}^{*
}
^{a} Instituto de Matemática, Estatística e Computação Científica, Universidade Estadual de Campinas, Campinas, SP, Brasil. Email: geovaneah@gmail.com, boldrini@ime.unicamp.br
^{b} Departamento de Sistemas Integrados, Faculdade de Engenharia Mecânica, Universidade Estadual de Campinas, Campinas, SP, Brasil. Email: eduardo.moraes15@yahoo.com.br, mlb@fem.unicamp.br
Abstract
In this work, we apply semi and fullyimplicit time integration schemes to the damage and fatigue phase field presented in Boldrini et al. (2016). The damage phase field is considered a continuous dynamic variable whose evolution equation is obtained by the principle of virtual power. The fatigue phase field is a continuous internal variable whose evolution equation is considered as a constitutive relation to be determined in a thermodynamically consistent way. In the semiimplicit scheme, each equation is solved separately by suited implicit method. The Newton’s method is used to linearize the equations in the fullyimplicit scheme. The time integration methods are compared and the results of damage and fracture evolution under the influence of fatigue effects are presented. The computational cost associated to the semiimplicit scheme showed be lower than the fully counterpart.
Keywords Damage; fatigue; phase field; finite element method
1 INTRODUCTION
The effects of the fatigue are of great concern when designing mechanical components, as they account for the majority of mechanical failures. Modeling and simulation of fatigue accumulation can lower development time and reduce prototype and testing costs ( ^{Lemaitre and Desmorat, 2005} ). However, some aspects such as crack initiation are rather difficult to handle in the classic models. The correct and robust assessment of such issues is challenging, but the introduction of phase field models for the evolution of damage and fatigue can overcome the main difficulties.
Originally, the phase field model was developed by John W. Cahn and John E. Hilliard to describe the phase separation in fluids through the fourthorder CahnHilliard equation ( ^{Cahn and Hilliard, 1958} ). Later on, Sam Allen and John Cahn developed the classical AllenCahn equation ( ^{Allen and Cahn, 1972} ), a secondorder nonlinear partial differential equation modeling phase separation in iron alloys. Both equations are based on the GinzburgLandau free energy functional. Recently, researchers have given more attention to phase field models, using them in a wide range of applications such as biology and medicine fields. Some of the different uses of phase field models include the study of tumor growth ( ^{Silva, 2009} ), vesiclefluid interation ( ^{Du et al., 2007} ) and fluidstructure interaction problems ( ^{Sun et al., 2014} ).
Damage, crack propagation and fatigue have also used phase field models in recent years. The works on brittle fracture showed accurate crack propagation patterns for different benchmarks ( ^{Miehe et al., 2010a} ; ^{Miehe et al., 2010b} ; ^{Borden et al., 2014} ), see ^{Ambati et al. (2015a)} for an enlightening review on the phase field models of brittle fracture in elastic solids. ^{Duda et al. (2015)}
extends the phase field damage model to describe brittle fracture in elastoplastic materials. Effects of ductile fracture are captured by the models proposed in ^{Ambati et al. (2015b)} and ^{Borden et al. (2016)} . Impressive numerical results broadening the previous phase field models are presented. Damage and fatigue phase field models were developed by ^{Fabrizio et al. (2006)} , ^{Amendola and Fabrizio (2014)} and, ^{Caputo and Fabrizio (2015)} . In those models, an expression for fatigue is given from the beginning making generalization for more complex situations rather difficult. More recently, ^{Boldrini et al. (2016)} proposed a damage and fatigue model similar to ^{Amendola and Fabrizio (2014)} . However, in this case, fatigue is treated as an internal variable with a constitutive law found in a thermodynamically consistent way. Disregarding the fatigue variable in ^{Boldrini et al. (2016)} , a model which is closely related to, but more general than ^{Miehe et al. (2010a)} , is obtained.
^{Boldrini et al. (2016)} present a thermodynamically consistent model for nonisothermal evolution of damage and fatigue. The difference between the damage phase field
φ
and the fatigue internal variables
F
is in their scales: the first variable is related to meso/macrocracks and the second one associated to microcracks. The results are for 1D examples, solved using highorder finite elements ( ^{Bittencourt, 2014} ), and the explicit RungeKutta time integration scheme. The results had expected convergence rates for the
p
refinement by using manufactured solution. The simulations also showed that fatigue effects play a major role to drive the crack initiation under cyclic loading with small stresses (we remark, however, that fatigue may also evolve even under monotonic loading as discussed in Section 5.1). Temperature rises were also observed in simulations of the nonisothermal equations, in agreement with experimental observations.
Although the 1D results showed expected behavior and the explicit RungeKutta time integration method performed well, the solution of 2D and 3D problems is more complex due to the high computational costs. The use of an explicit time integration scheme requires smaller time increments to achieve stability, compared to the implicit methods. In this way, high dimensional problems demand finer meshes and smallest time increments for better approximations, making the use of explicit time integration methods in most cases unfeasible. On the other hand, a fullyimplicit method is generally costly due to the problem and the involved nonlinearities.
In the present work, we propose a semiimplicit scheme to solve the isothermal model for the evolution of damage and fatigue and which provides the desired stability, even for large time steps. In this scheme, each equation is treated separately by the suited implicit method, keeping the others variables frozen. More specifically, for each time step the evolution of one variable is solved and the results are used as input for equations of the other variables. In order to show the efficiency and accuracy of this scheme, the results are compared with the fullyimplicit method obtained using the Newton’s procedure.
Section 2 presents the phase field model and the governing equations for the isothermal case and linear elastic brittle material. Section 3 shows the finite element spatial discretization for the proposed model. Section 4 describes the semi and fullyimplicit schemes. For the semiimplicit case, different time integration methods are presented. Among the methods covered are Newmark, backward Euler, trapezoidal method and an implicit 2ndorder RungeKutta scheme called TRBDF2. In the fullyimplicit scheme, the residue and the Jacobian matrix of the Newton’s procedure are described based on the Newmark and trapezoidal methods. Section 5 describe the results for an Ishaped specimen under monotonic and cyclic loadings. From these numerical results, a convergence rate analysis is considered for different combinations of the proposed methods using the semiimplicit scheme. A comparison of the accuracy and computational time cost between the semi and fullyimplicit closes this section. Finally, Section 6 outlines the conclusions of the numerical example and the aspects of the proposed schemes.
2 GOVERNING EQUATIONS
The work of ^{Boldrini et al. (2016)} presents a framework for structural damage and fatigue using thermodynamically consistent and nonisothermal phase field models. The coupled dynamic equations for the evolution of displacement
u
, velocity
u˙
and damage phase field
φ
are obtained by applying the principle of virtual power and the entropy inequality based on the specific free energy potential. The damage phase field
φ
represents the volumetric fraction of damaged material with
φ=0
for a virgin material,
φ=1
for a fractured material and
0<φ<1
for damaged material between virgin and fully broken states. The evolution equation for the internal fatigue variable
F
is derived by finding a differential constitutive relation for
F
that satisfies the entropy inequality for all possible admissible processes.
From suitable choices of freeenergy, initial and boundary conditions, and disregarding the temperature evolution, the isotropic isothermal case described in Eq. (62) of ^{Boldrini et al. (2016)} is
{u¨=div((1−φ)2Cρ0E)+bρ0div(D)−γgcρ0div(∇φ⊗∇φ)+fφ˙=γgcλΔφ+1λ(1−φ)ETCE−1λγ[gcH'(φ)+FHf'(φ)]F˙=−F^γHf(φ) ,
(1)
where
E=∇su
,
D=∇su˙
and
C
are, respectively, the infinitesimal strain tensor, the rate of strain tensor and the elasticity tensor (using Voigt notation) given in terms of the Young’s modulus
E
and Poisson’s coefficient
ν
. Coefficient
b
refers to the viscous damping of the material,
gc
is the Griffith fracture energy,
ρ0
corresponds to the material’s density and
γ>0
is a parameter related to the width of the damage phase field layers. The parameter
λ
is associated to the rate of damage change and should increase proportional to the damage as described in ^{Lemaitre and Desmorat (2005)} . It is given as
1λ=c(1+δ−φ)ζ ,
(2)
where the parameters
c
and
ζ
are positive and dependent on the material and
δ
is a small positive constant used to avoid singularity. In the present work, it is assumed
ζ=1
for simplicity.
Additionally, suitable potentials
H
and
Hf
are needed. The role of these potentials is to drive the damage from zero to one as the fatigue changes from zero to
gc
. Different potentials change the way that damage varies by intermediate values of fatigue. A suitable choice for the potentials is such that damage changes in a monotonic and continuous way. Here we adopt
H(φ)={0.5φ2for0≤φ≤10.5+δ(φ−1)forφ>1−δφforφ<0andHf(φ)={−φfor0≤φ≤1−1forφ>10forφ<0 .
(3)
Remark 1: The complete definitions of the potentials
H(φ)
and
Hf(φ)
, that is, their expressions for all
φ
, are necessary so that the hypotheses of Theorem 1 presented in ^{Boldrini et al. (2016)}on page 405 can be verified. That theorem states that necessarily
0≤φ≤1
for all
t>t0
for sufficiently regular solutions for Eq. (1), with initial condition values in
[0,1]
at time
t0
, homogeneous Neumann boundary condition (
∂φ/∂n=0
) and under some additional hypotheses on the potentials (which the previously given
H(φ)
and
Hf(φ)
fulfill).
This means that for solutions of the continuous model (1), the damage variable
φ
assumes values only in the interval
[0,1]
, and thus only the parts of the potentials
H(φ)
and
Hf(φ)
associated to values of
φ
in the interval
[0,1]
are in fact used during the evolution.
However, another reason for given the complete definition of the potentials is that, when numerically solvingEq. (1)there is the possibility that some computed values of the approximation of
φ
fall slightly outside
[0,1]
. This may happen due to discretization error, especially for rather large time steps, and thus, to guarantee that the computations may proceed, the complete definitions of the potentials are required.
Finally, the parameter
F^
for the fatigue evolution equation is based on the microcracks that are formed by cyclic loadings, even small ones. The microcrack formation depends on how large is the stress associated to the virgin material stress (see Remark 1). Choosing a linear dependence in terms of a parameter
a
, we have
F^=a(1−φ)‖CE‖ .
(4)
From the choice of the potential
Hf
in Eq. (3b) and
F^
as above, the evolution equation of the fatigue in Eq. (1c) assume
F˙≥0
in any point of the material. In particular if
φ=1
, the fatigue
F
variable stops increasing.
There are many possibilities for the boundary conditions in Eq. (1). For Eq. (1a), there are the cases where either the displacement or the stresses are given at certain parts of the boundary. In the case of the damage variable
φ
, the standard boundary condition is a homogeneous Neumann condition (null flux for
φ
at the boundary). There is no boundary condition associated to the fatigue
F
since Eq. (1c) is locally an ordinary differential equation.
Remark 2: Eq. (4)differs slightly of the similar expression presented originally in^{Boldrini et al. (2016)}, which considers the absolute value of the power associated to the virgin material stress instead of the stress norm. This new choice of
F^
is due to the following physical observation: the fatigue of a material subjected to alternating loads increases proportionally to the number of cycles, rather than the velocity at which the cycles were performed. In fact, there are many possibilities to choose
F^
. For example, instead of the stress norm, the von Mises stress or a quadratic dependence (instead a linear as inEq. (4)) may be used. These choices should be made in order to better describe the problem of interest.
Also remark that the main advantage of using a model including the fatigue variable as in^{Boldrini et al. (2016)}is the possibility of reaching a rather general result for the possible effects of fatigue, still preserving thermodynamical consistency. In fact, in the derivation of the model presented in^{Boldrini et al. (2016)}, from the beginning it was just assumed that fatigue had to satisfy a differential constitutive relation, but whose specific expression and how it would interact with the other physical variables had to be determined at the final steps of the entropy condition argument. This process leaded to a general expression depending on the freeenergy potential and pseudopotential of dissipation, and thus can be applied to several classes of materials. If we had restricted the form of the constitutive expression or specified how the fatigue should interacted with the other constitutive relations, the results would be also much more restricted.
Remark 3: We observe thatEq. (1)does not guarantee
φ˙≥0
. Then, the presented model allows the possibility of healing of meso and macro cracks. It is possible to adapt the model in order to prevent this healing and then have irreversible damage. For this, a history of elastic energy similar to^{Miehe et al. (2010b)}or penalty arguments can be applied. In this work, the simpler strategy proposed by the first one is adopted.
3 SPATIAL DISCRETIZATION
This section presents the finite element approximation of Eq. (1) . For this purpose, each equation is multiplied by suitable test functions
ω1
,
ω2
and
ω3
, and integrated over the domain
Ω
of the body. By using the Gauss theorem and after the application of the boundary conditions, we obtain the following weak form for Eq. (1):
{∫Ωu¨ω1dΩ=−∫Ω1ρ0(1−φ)2CE:∇ω1dΩ−∫Ωbρ0D:∇ω1dΩ++∫Ωγgcρ0(∇φ⊗∇φ):∇ω1dΩ+∫Ωf.ω1dΩ+B.T.∫Ωφ˙ω2dΩ=−∫Ωγgcλ∇φ.∇ω2dΩ+∫Ω1λ(1−φ)ETCEω2dΩ+−∫Ω1λγ[gcH'(φ)+FHf'(φ)]ω2dΩ∫ΩF˙ω3dΩ=−∫ΩF^γHf(φ)ω3dΩ .
(5)
The boundary terms that appear in the above weak form are denoted by
B.T.
. They depend on the boundary conditions being considered. Due the assumed homogeneous boundary condition for the damage
φ
, there are no extra boundary terms in Eq. (5b).
Following the similar 2D spatial discretization presented in ^{Chiarelli et al. (2017)} , we use the finite element method to approximate the above equations. For this purpose, we consider a finite element mesh such that
Ω=Ae=1nelemΩeandΓ=Ae=1nelemΓe ,
(6)
where
Ωe
and
Γe
are the domain and the boundary of each element
e
,
A
represents the assembly procedure and
nelem
is the number of elements. In this way, the approximation in each element
e
is written as the linear combination of the
η
local nodal basis functions
Nj
as
ue=Neu^e,u˙e=Neu˙^e,fe=Nef^e,φe=Nφeφ^eandFe=NFeF^e ,
(7)
where
Ne
,
Nφe
,
NFe
,
u^e
,
u^˙e
,
φ^e
and
F^e
are defined by
Ne=[N10N20...Nη00N10N2...0Nη],
(8)
Nφe=NFe=[N1N2...Nη],
(9)
u^e=[u1xeu1yeu2xeu2ye...uηxeuηye]T,
(10)
u^˙e=[u˙1xeu˙1yeu˙2xeu˙2ye...u˙ηxeu˙ηye]T,
(11)
φ^e=[φ1eφ2e...φηe]T,
(12)
F^e=[F1eF2e...Fηe]T.
(13)
The finite element interpolations for the gradients of displacement, velocity and damage for plane state are given in terms of linear combinations of the shape function global derivatives by
Ee=Bueu^e,De=Bveu^˙eand∇φe=Bφeφ^e ,
(14)
where the derivative matrices are defined as
Bue=[N1,x0N2,x0...Nη,x00N1,y0N2,y...0Nη,yN1,yN1,xN2,yN2,x...Nη,yNη,x] ,
(15)
Bve=diag(1,1,2/2)Bue ,
(16)
Bφe=[N1,xN2,x...Nη,xN1,yN2,y...Nη,y] .
(17)
Replacing the above approximations in Eq. (5) , we obtain the semidiscrete form for the
e
th element as
{Meu¨^e=Kueu^e+Kveu˙^e+wae+Mef^eMφeφ˙^e=(Pφe+Kce)φ^e+wbe+wceMFeF˙^e=wde
(18)
with the operators defined by
Me=∫Ωe(Ne)T(Ne)dΩe ,
(19)
Mφe=∫Ωe(Nφe)T(Nφe)dΩe ,
(20)
MFe=∫Ωe(NFe)T(NFe)dΩe ,
(21)
Kue=−∫Ωe1ρ0(1−Nφeφ^e)2(Bue)TC(Bue)dΩe ,
(22)
Kve=−∫Ωebρ0(Bve)T(Bve)dΩe
(23)
Pφe=−∫Ωeγgcλ(Bφe)T(Bφe)dΩe−∫Ωegcγλ(Nφe)T(Nφe)dΩe ,
(24)
Kce=−∫Ωe1λ(Bueu^e)TC(Bueu^e)(Nφe)T(Nφe)dΩe ,
(25)
wae=∫Ωeγgcρ0(Bφeφ^e⊗Bφeφ^e):BuedΩe ,
(26)
wbe=∫Ωe1λ(Bueu^e)TC(Bueu^e)(Nφe)TdΩe ,
(27)
wce=∫Ωe1λγ(NFeF^e)(Nφe)TdΩe ,
(28)
wde=∫Ωeaγ(1−Nφeφ^e)(Nφeφ^e)‖CBueu^e‖(NFe)TdΩe .
(29)
Note that the integrand in Eq. (26) results in a secondorder tensor which is expressed as a vector using Voigt notation. Also, the potentials
H(φ)
,
Hf(φ)
and
Hf'(φ)
that appear in the damage and fatigue equations are given in Eq. (3) for
0≤φ≤1
. The global form of the above matrices and vectors are obtained by applying the standard assembly procedure
A
and indicated without the superscript
e
. Appropriate choice of the elasticity tensor
C
define particularization for plane stress or strain problems. Both could be considered; however, just the last one is used in this work.
4 TIME INTEGRATION SCHEMES
This section presents different time integration methods used by the semi and the fullyimplicit procedures to solve Eq. (18) . In the semiimplicit scheme, each equation is solved separately by an implicit method. Firstly, the damage evolution equation is solved using one of three methods presented below, namely the backward Euler, trapezoidal and the 2stage secondorder accurate diagonally implicit RungeKutta method (TRBDF2). The damage solution is replaced in the equation of motion which is solved by the standard Newmark method. Finally, the current displacement, velocity and damage fields are used to evaluate the fatigue by the backward Euler or trapezoidal methods. This procedure is repeated for each time step. The fullyimplicit procedure is based on the iterative technique of the Newton’s method, by applying the Newmark for Eqs. (18a) and trapezoidal method for Eqs. (18b) and (18c).
Consider the split of the solution time interval
[0,T]
in discrete time steps
tn
with time increments given by
Δt=tn+1−tn>0
for
n=0,1,...
. The global approximations for the variables are denoted at time step
tn+1
as
un+1=u^(tn+1),u˙n+1=u˙^(tn+1),φn+1=φ^(tn+1)andFn+1=F^(tn+1) .
(30)
The integration procedures used in Eq. (18) are presented in the following sections. The semi and fullyimplicit algorithms are summarized at the end of their respective Sections.
4.1 Semiimplicit Scheme
This section presents the time integration procedures for Eq. (18) used in the semiimplicit scheme.
4.1.1 Damage Evolution Equation
The simplest method for time integration of Eq. (18b) in global form is the forward Euler, whose explicit scheme of evolution is given by
φn+1=[I+Δt(Mφ)−1(Pφ,n+Kc,n)]φn+Δt(Mφ)−1(wb,n+wc,n) ,
(31)
where
I
is an appropriate identity matrix. In a very similar way, now considering the righthand side in the current time step, we have the implicit backward Euler. In this case, the timemarching rule is given by the following system to be solved at each time step
[Mφ−Δt(Pφ,n+1+Kc,n+1)]φn+1=Mφφn+Δt(wb,n+1+wc,n+1) .
(32)
Performing an average between the above two Euler methods, we obtain the implicit trapezoidal method. Unlike the Euler methods that are only first order accurate for linear problems, this method is second order accurate due to its symmetric approximation. Applying the trapezoidal method in Eq. (18b) results in the following time evolution rule:
[Mφ−Δt2(Pφ,n+1+Kc,n+1)]φn+1=[Mφ+Δt2(Pφ,n+Kc,n)]φn+−Δt2(wb,n+1+wb,n)−Δt2(wc,n+1+wc,n).
(33)
This work also proposes the use of the 2stage secondorder accurate diagonally implicit method, or simply TRBDF2 method. As the backward Euler, the TRBDF2 is a Lstable method, indicated for the time stiff problems ( ^{LeVeque, 2007} ). The first stage consists of the trapezoidal method applied over
Δt
/2 obtaining the solution
φ*
. In the second stage, the solution is updated by solving the following system:
[Mφ−Δt3(Pφ,n+1+Kc,n+1)]φn+1=Mφ[43φ*−13φn]+Δt3(wb,n+1+wc,n+1) .
(34)
Unlike the others, this method requires the solution of two linear systems at each time step.
The above operators depend on the current variables, an exception is the forward Euler. However, we use the values at the previous time steps in the semiimplicit scheme. For example, in Eq. (32) is used
φn
,
un
and
Fn
instead of
φn+1
,
un+1
and
Fn+1
, respectively. This makes possible to perform the assembly of the operators, and thus the time integration. In this way, the backward Euler method applied to the damage equation becomes
[Mφ−Δt(Pφ,n+Kc,n)]φn+1=Mφφn+Δt(wb,n+wc,n) .
(35)
The same idea is used for all the above methods.
4.1.2 Equation of Motion
Consider the system of ordinary differential equations Eq. (18a) given in the global form by
Mu¨^=Kuu^+Kvu˙^+wa+Mf^ .
(36)
In order to approximate the solution for the displacement
un+1
, the implicit standard Newmark procedure is adopted ( ^{Newmark, 1952} ). In this method, the acceleration and velocity are evaluated using the updated values of the displacement by
u¨n+1=α1(un+1−un)−α2u˙n−α3u¨n ,
(37)
u˙n+1=α4(un+1−un)+α5u˙n+α6u¨n ,
(38)
where the coefficients
αj
(j=1,…,6)
are given in terms of the Newmark coefficients
γ˜
and
β˜
(we assume the values
γ˜=0.5
and
β˜=0.25
) by
α1=1β˜Δt2,α2=1β˜Δt,α3=1−2β˜2β˜,α4=γ˜β˜Δt,α5=1−γ˜β˜,andα6=[1−γ˜2β˜]Δt .
(39)
Substituting the above approximations in Eq. (36) , we obtain the following timemarching system:
[α1M−Ku−α4Kv]un+1=M[α3u¨n+α2u˙n+α1un]+Kv[α6u¨n+α5u˙n−α4un]+wa+Mfn+1 .
(40)
After computing the displacement, the acceleration and velocity fields are updated using Eqs. (37) and (38) . The imposition of prescribed displacement
u¯(tn+1)
produces the
B.T.
terms in Eq. (5a). This results in an additional step of update, where we should impose the velocity and acceleration at the prescribed boundary as
u˙¯n+1=ddtu¯(tn+1)andu¨¯n+1=d2dt2u¯(tn+1) ,
(41)
where the bar symbol represents the prescribed degrees of freedom.
4.1.3 Fatigue Equation
For the fatigue equation Eq. (18c), we use the backward Euler and trapezoidal methods whose time evolution expressions are given respectively by
MFFn+1=MFFn+Δtwd,n+1 ,
(42)
and
MFFn+1=MFFn+Δt2(wd,n+1+wd,n) .
(43)
Note that both methods have basically the same computational cost, since
wd,n
can be treated as a stored vector.
The complete algorithm for the solution of the system is summarized in Algorithm 1 .
Algorithm 1 Algorithm of the semiimplicit time integration scheme.
Semiimplicit time integration scheme 
1: for
t
from
0
to
tf
do

2: Given
un
,
u˙n
and
λn
use one of the Eqs. (32 – 34) to obtain
φn+1
; 
3: Given the updated damage
φn+1
, solve Eq. (40) to obtain
un+1
; 
4: From the current displacement
un+1
, update the acceleration
u¨n+1
and velocity
u˙n+1
using 
Eqs. (37) and (38), respectively; 5: Update the prescribed velocity and acceleration using Eq. (41); 
6: Given
un+1
and
φn+1
, update the fatigue
Fn+1
using Eq. (42) or (43); 
7: Update the time step. 
8: end

4.2 Fullyimplicit Scheme
This section presents the fullyimplicit scheme to solve Eq. (18) obtained from the Newton’s method. Based on the idea of linear approximation, such method is generally employed to solve nonlinear equations by an iterative sequence of linear systems ( ^{LeVeque, 2007} ).
Consider a general nonlinear system given by
r(v)=0 .
(44)
This system is solved employing some iterative procedure instead of a direct method. Supposing that
v(i)
is the approximation to
v
at iteration
i
, Newton’s method estimates the solution at iteration
i
via the Taylor series expansion
r(v(i+1))=r(v(i))+∂r(v(i))∂v(i)(v(i+1)−v(i))+... .
(45)
Dropping the higher order terms and setting
r(v(i+1))=0
results in
r(v(i))+∂r(v(i))∂v(i)(v(i+1)−v(i))=0 ,
(46)
giving the following update step
v(i+1)=v(i)+δv(i) ,
(47)
where
δv(i)
solves the linear system
J(v(i))δv(i)=−r(v(i)) .
(48)
The Jacobian matrix
J
and its coefficients are respectively defined as
J(v)=∂r(v(i))∂v(i)andJpq(v(i))=∂rp(v(i))∂vq(i) .
(49)
Newton’s method is summarized in Algorithm 2 .
Algorithm 2 Algorithm of the Newton’s method.
Fullyimplicit time integration scheme 
1: Specify an initial guess
v(0)=v0
and a prescribed tolerance tol; 2: while tol < error do

3: Evaluate the Jacobian
J(v(i))
and the residue
r(v(i))
; 
4: Solve the system
J(v(i))δv(i)=−r(v(i))
for
δv(i)
; 
5: Update the variable
v(i+1)←v(i)+δv(i)
; 6: Update the iteration
i←i+1
; 7: Evaluate the error as
error=max[‖r(v(i))‖,‖δv(i)‖/‖v(i+1)‖]
; 8: end

9:
v←v(i+1)

In order to solve Eq. (18) using the Newton’s method, the nonlinear system (residue) is defined by applying the Newmark in Eq. (18a) and trapezoidal method in Eqs. (18b) and (18c). This choice of time integration methods is due their high orders of accuracy presented. From Eqs. (33) , (40) , (43) and (44) the discretized nonlinear system is
r(un+1,φn+1,Fn+1):=[r1r2r3]T=0 ,
(50)
whose components are
r1:=[α1M−Ku,n+1−α4Kv,n+1]un+1−M[α3u¨n+α2u˙n+α1un]+−Kv,n+1[α6u¨n+α5u˙n−α4un]−wa,n+1−Mfn+1,
(51)
r2:=[Mφ−Δt2(Pφ,n+1+Kc,n+1)]φn+1−[Mφ+Δt2(Pφ,n+Kc,n)]φn+−Δt2(wb,n+1+wb,n)−Δt2(wc,n+1+wc,n),
(52)
r3:=MF(Fn+1−Fn)−Δt2(wd,n+1+wd,n) .
(53)
The Jacobian matrix can be written in the following block matrix form:
J=[J11J12J13J21J22J23J31J32J33] ,
(54)
with coefficients defined by
Jp1=∂rp∂un+1,Jp2=∂rp∂φn+1andJp3=∂rp∂Fn+1,p=1,2,3 .
(55)
For the
e
th element, the coefficients of the Jacobian matrix are
J11e=α1Me−Kue−α4Kve ,
(56)
J12e=−2∫Ωe1ρ0(1−Nφeφ^e)(Bue)TCBueu^eNφedΩe−∫Ωeγgcρ0(Bue)TB˜eBφedΩe ,
(57)
J13e=0 ,
(58)
J21e=−Δt∫Ωe1λ(1−Nφeφ^e)(Ne)T(Bueu^e)TC(Bue)dΩe,
(59)
J22e=Mφe−Δt2Pφe+Δt2∫Ωeγgcλ¯(Bφe)T(Bφeφ^e)(Nφe)dΩe++Δt2∫Ωe1γλ¯(gcNφeφ^e−NFeF^e)(Nφe)T(Nφe)dΩe,
(60)
J23e=−Δt2∫Ωe1γλ(Nφe)T(Nφe)dΩe ,
(61)
J31e=−Δt2∫Ωeaγ‖CBueu^e‖(1−Nφeφ^e)(Nφeφ^e)(NFe)T(CBueu^e)TCBuedΩe(u^e≠0) ,
(62)
J32e=Δt2∫Ωeaγ(2Nφeφ^e−1)‖CBueu^e‖(NFe)T(NFe)dΩe ,
(63)
J33e=MFe ,
(64)
where
B˜e
and 1/
λ¯
in Eqs. (57) and (60) are given respectively by
B˜e=2∑i=1η[Ni,xeφ^ie00Ni,yeφ^ieNi,yeφ^ieNi,xeφ^ie]and1λ¯=c(1+δ−φ)2 .
(65)
The Jacobian matrix is nonsymmetric, since
J31≠J13T
for any
u^≠0
.
5 NUMERICAL RESULTS
This section presents the numerical results obtained using the semi and fullyimplicit time integration schemes. The problem consists of the Ishaped specimen with dimensions, boundary conditions and cartesian system of coordinates given in Figure 1 . A similar problem disregarding the effects of fatigue was proposed in ^{Chiarelli et al. (2017)} . We use a mesh of quadratic quadrilaterals (serendipity element of
η=8
nodes) totalizing 2700 nodes. In the present simulation, the following material and phase field parameters are used:
E=
180GPa,
ν=
0.3,
ρ=
7300Kg/m^{3},
γ=
0.5mm,
b=5×107
Ns/m^{2},
c=5×10−4
m/Ns,
gc=
6000N/m and
a=8×10−9
m^{2}. For the iterative procedure in the fullyimplicit scheme, the tolerance is
tol=10−8
.
In order to verify the evolution of damage and fatigue, two prescribed displacement cases are considered. In the first one, a displacement of 0.1mm/s is incrementally applied through the time steps on the bottom edge, denominated case A. In the second, a sinusoidal displacement (given in millimeters) of
u¯(t)={−0.167sin(10πt)for0≤t<0.05−0.075sin(20πt−1.6)−0.092fort≥0.05 ,
(66)
is applied on the same edge, called case B. A convergence study of the time integration schemes is considered in the first case. For future use, the effective and softened stress are defined as
σeff=CE ,
(67)
σsof=(1−φ)2CE+bD−γgcsym(∇φ⊗∇φ) ,
(68)
where the last term in Eq. (68) is expressed as a vector using Voigt notation. The last two terms in Eq. (68) have a small contribution to the values of the softened stress. The results are presented in the following sections.
5.1 Case A
The monotonic displacement of 0.1mm/s is incrementally applied along the time range
[0,1.6]s
on the bottom edge aims to verify the growth of the crack path in the specimen without a remarkable influence of fatigue. For illustrative purposes, Figure 2 exhibits the crack path evolution for particular time steps. Figure 3 illustrates the damage, fatigue and stress evolution for the considered time interval in the central node of the specimen. Note the fast growth of the damage levels for a small increment of time close to
t=1.4s
, characterizing a time stiff behavior. Also, observe that the damage variable increases quickly over the time without any significant fatigue values, while the effective stress rise as time progresses. In this case, failure occurs due to the high stress levels, as expected from a tensile test. If the fatigue equation is disregard, the damage response is qualitatively the same, their values, however, are slightly smaller.
This problem is solved using different time integration schemes. It was verified that the use of any RungeKutta explicit method results in numerical instability even for small increments of time as
Δt=10−8
. In order to avoid the use of small time increments, an alternative is adopting implicit time methods. Therefore, we use the Newmark procedure to solve the motion equation and the procedures previously presented in Section 4.1 for the damage and fatigue equations. In particular for the fullyimplicit scheme, these last two equations are solved by using the trapezoidal method.
The use of different combinations of methods and schemes resulted in the convergence rates and CPU times given in Table 1 . The convergence rates are obtained using the relative error evaluated at
t=1.5s
, given for a generic vector
s
and the seminorm by
ε(s)=‖s−sref‖2‖sref‖2withs=∫Ω∇s⋅∇sdΩ ,
(69)
where the reference values
sref
are computed using the fullyimplicit scheme with time increment of
Δt=10−4
.
Table 1 Convergence rates for different combinations of time integration methods and the related CPU time obtained in
t=1.5s
for time increment of
Δt=10−3
.
Scheme 
Equation 
Convergence ratio 
CPU Time 
Fatigue 
Damage 
Displacement 
Velocity 
Damage 
Fatigue 
Semi 
BE 
BE 
1.62 
1.30 
1.64 
0.78 
2580 
TR 
1.98 
1.71 
1.97 
0.21 
2558 
TRBDF2 
1.58 
1.38 
1.76 
0.99 
2601 
TR 
BE 
1.63 
1.33 
1.64 
0.53 
2585 
TR 
1.92 
1.45 
1.15 
0.28 
2587 
TRBDF2 
1.58 
1.40 
1.36 
0.97 
2626 
Fully 
TR 
TR 
2.07 
2.06 
2.06 
2.06 
22512 
Due to the decoupled time integration in the semiimplicit scheme, the convergence rates are below the expected theoretical values, even with the use of secondorder methods. Arising from the stiff behavior, the use of the Lstable TRBDF2 method presented the highest average rates of convergence when compared to the others methods for the damage equation. No significant difference is observed in the rates of the trapezoidal method in the fatigue equation. In general, all combinations of methods presented the same computational cost. The error values when using trapezoidal and TRBDF2 for fatigue and damage equations are illustrated in Figure 4 a. In order to obtain error below
5%
for the velocity, a time increment of
Δt=10−4
is required. For the other variables as displacement, damage and fatigue, the increment of
Δt=5×10−4
is sufficient.
The full scheme shows quadratic rate, as expected due the use of secondorder accurate methods. The error values illustrated in Figure 4 b are below
0.5%
for all variables. However, the computational cost is approximately
850%
higher than the semiimplicit scheme. In this analysis, the Jacobian matrix conditioning number had values of order
K(J)=1017
.
5.2 Case B
In this case, the sinusoidal displacement given in Eq. (66) is applied along the time range
[0,2.5]s
. The objective of this new boundary condition is to verify the growth of the crack path in the specimen with influence of the fatigue. For illustrative purposes, Figure 5 shows the damage evolution in the Ishaped specimen for particular time steps.
Conversely to the monotonic displacement, the growth of damage is largely due to fatigue. This can be seen in Figure 6 a, which illustrates the evolution of damage and fatigue along the time interval for the central node. Note that the fatigue grows at an approximately constant rate over the time, reaching large values close to failure. This results in an abruptly growth of the damage close to 2.35s, characterizing again the stiff behavior of the damage equation. In absence of fatigue, the damage grows only due the contribution of the elastic energy and their value is limited by
≈0.1
. Figure 6 b shows the effective and softened stress levels evaluated at the central node of the specimen. The effective stress does not exhibit remarkable growth in magnitude over time when compared to the case A, except near failure.
6 CONCLUSIONS
This work proposed the use of semi and fullyimplicit schemes to perform the time integration of the phase field equations of damage and fatigue presented in ^{Boldrini et al. (2016)} . The explicit RungeKutta methods showed to be unfeasible for 2D examples, presenting instability with time increments as small as
Δt=10−8
. An alternative to avoid the use of small time steps maintaining stability is to use implicit methods. The fullyimplicit method presented quadratic convergence and lower errors when compared to the semiimplicit procedures. However, the fullyimplicit scheme has some drawbacks:

The Jacobian of Newton’s method is hard to obtain and becomes particular to the problem;
The Jacobian matrix is nonsymmetric;
There is a high computational cost and memory demand, since it is necessary to solve a large system in each Newton’s iteration;
The semiimplicit procedure overcomes these difficulties. In this scheme, the equations are decoupled and solved separately by suited implicit time integration procedures. Different alternatives were used for the damage equation and the TRBDF2 method showed better accuracy with no significant additional cost. This method is Lstable and therefore suitable to solve a stiff problem as the damage equation. For the fatigue evolution equation, the trapezoidal and backward Euler methods performed similarly. For the equation of motion, only the standard Newmark procedure was used. The convergence rates for the semiimplicit scheme were below ranging close to the linear rate in most of the equations. Accurate responses (error below
5%
) for displacement, damage and fatigue are possible to obtain by using time increments of
Δt=5×10−4
. Although the velocity has high errors, this fact did not affect the displacement solution and consequently the damage and fatigue variables.
Although not applied in this work, an iteration of correction in the semiimplicit procedure would reduce errors. In particular, the prediction/correction AdamsBashforth/Moulton methods can be used to improve accuracy. Further study on this subject will be carried out in the future.
In general, the numerical results for the Ishaped specimen for tensile and fatigue tests show the expected behavior qualitatively. The model proposed in ^{Boldrini et al. (2016)} makes possible to represent the evolution of the damage under cyclic loads due to the fatigue. This feature provides the identification of regions of initiation and propagation of cracks efficiently with a reasonable computational cost.
As the next step, we will adjust this model quantitatively with experimental tests.
7 ACKNOWLEDGEMENTS
Authors would like to thank the São Paulo Research Foundation (FAPESP) under grants 2013/502383, 2015/103102 and 2015/201880, the Coordination for the Improvement of Higher Education Personnel (CAPES) under grant 33003017, and also the National Council for Research and Development (CNPq) under grant 306182/20149.
References
Allen, S.M., Cahn, J.W., (1972). Ground state structures in ordered binary alloys with second neighbor interactions. Acta Metallurgica, 20(3):423–433.
[ Links ]
Ambati, M., Gerasimov, T., de Lorenzis, L., (2015a). A review on phasefield models of brittle fracture and a new fast hybrid formulation. Journal of Computational Mechanics, 55(2):383–405.
[ Links ]
Ambati, M., Gerasimov, T., de Lorenzis, L., (2015b). Phasefield modeling of ductile fracture. Journal of Computational Mechanics, 55(5):1017–1040. doi:10.1007/s0046601511514.
[ Links ]
Amendola, G., Fabrizio, M., (2014). Thermomechanics of damage and fatigue by a phase field model. arXiv preprint arXiv:1410.7042.
[ Links ]
Bittencourt, M.L., (2014). Computational solid mechanics: Variational formulation and high order approximation. CRC Press.
[ Links ]
Boldrini, J.L., de Moraes, E.A.B., Chiarelli, L.R., Fumes, F.G., Bittencourt, M.L., (2016). A nonisothermal thermodynamically consistent phase field framework for structural damage and fatigue. Computer Methods in Applied Mechanics and Engineering, 312(395–427).
[ Links ]
Borden, M.J., Hughes, T.J.R., Landis, C.M., Verhoosel, C.V., (2014). A higherorder phasefield model for brittle fracture: Formulation and analysis within the isogeometric analysis framework. Computer Methods in Applied Mechanics and Engineering, 273:100–118.
[ Links ]
Borden, M.J., Hughes, T.J.R., Landis, C.M., Anvari, A., Lee, I.J., (2016). A phasefield formulation for fracture in ductile materials: Finite deformation, balance law derivation, plastic degradation, and stress triaxiality effects. Computer Methods in Applied Mechanics and Engineering, 312:130–166.
[ Links ]
Cahn, J.W., Hilliard, J.E., (1958). Free energy of a nonuniform system. I. Interfacial free energy. The Journal of chemical physics, 28(2):258–267.
[ Links ]
Caputo, M., Fabrizio, M., (2015). Damage and fatigue described by a fractional derivative model. Journal of Computational Physics, 293:400–408.
[ Links ]
Chiarelli, L.R., Fumes, F.G., de Moraes, E.A.B., Haveroth, G.A., Boldrini, J.L., Bittencourt, M.L., (2017). Comparison of high order finite element and discontinuous Galerkin methods applied to phase field equations: Application to structural damage. Computers & Mathematics with Applications, http://dx.doi.org/10.1016/j.camwa.2017.05.003 .
[ Links ]
Du, Q., Li, M., Liu, C., (2007). Analysis of a phase field NavierStokes vesiclefluid interaction model. Discrete and Continuous Dynamical Systems Series B, 8(3):539.
[ Links ]
Duda, F. P., Ciarbonetti, A. Sánchez, P. J. Huespe, A. E., (2015). A phasefield/gradient damage model for brittle fracture in elasticplastic solids. International Journal of Plasticity, 65:269296
[ Links ]
Fabrizio, M., Giorgi, C., Morro, A., (2006). A thermodynamic approach to nonisothermal phasefield evolution in continuum physics. Physica D: Nonlinear Phenomena, 214(2):144–156.
[ Links ]
Lemaitre, J., Desmorat, R., (2005). Engineering damage mechanics: ductile, creep, fatigue and brittle failures. Springer Science & Business Media.
[ Links ]
LeVeque, R.J., (2007). Finite difference methods for ordinary and partial differential equations: steadystate and timedependent problems, Society for Industrial and Applied Mathematics.
[ Links ]
Miehe, C., Hofacker, M., Welschinger, F., (2010a). A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering, 199(45):2765–2778.
[ Links ]
Miehe, C., Welschinger, F., Hofacker, M., (2010b). Thermodynamically consistent phasefield models of fracture: Variational principles and multifield FE implementations. International Journal for Numerical Methods in Engineering, 83(10):1273–1311.
[ Links ]
Newmark, N.M., (1952). Computation of dynamic structural response in the range approaching failure. Department of Civil Engineering, United States.
[ Links ]
Silva, S.P.L.D.d., (2009). Phasefield models for tumor growth in the avascular phase. Ms. Thesis. University of Coimbra, Portugal.
[ Links ]
Sun, P., Xu, J., Zhang, L., (2014). Full Eulerian finite element method of a phase field model for fluid–structure interaction problem.” Computers & Fluids, 90:1–8.
[ Links ]