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 fourth-order Cahn-Hilliard equation ( Cahn and Hilliard, 1958 ). Later on, Sam Allen and John Cahn developed the classical Allen-Cahn equation ( Allen and Cahn, 1972 ), a second-order nonlinear partial differential equation modeling phase separation in iron alloys. Both equations are based on the Ginzburg-Landau 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 ), vesicle-fluid interation ( Du et al., 2007 ) and fluid-structure 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 non-isothermal evolution of damage and fatigue. The difference between the damage phase field
Although the 1D results showed expected behavior and the explicit Runge-Kutta 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 fully-implicit method is generally costly due to the problem and the involved non-linearities.
In the present work, we propose a semi-implicit 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 fully-implicit 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 fully-implicit schemes. For the semi-implicit case, different time integration methods are presented. Among the methods covered are Newmark, backward Euler, trapezoidal method and an implicit 2nd-order Runge-Kutta scheme called TR-BDF2. In the fully-implicit 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 I-shaped 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 semi-implicit scheme. A comparison of the accuracy and computational time cost between the semi and fully-implicit 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 non-isothermal phase field models. The coupled dynamic equations for the evolution of displacement
From suitable choices of free-energy, initial and boundary conditions, and disregarding the temperature evolution, the isotropic isothermal case described in Eq. (62) of Boldrini et al. (2016) is
where
where the parameters
Additionally, suitable potentials
Remark 1: The complete definitions of the potentials
This means that for solutions of the continuous model (1), the damage variable
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
Finally, the parameter
From the choice of the potential
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
Remark 2: Eq. (4)differs slightly of the similar expression presented originally inBoldrini 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
Also remark that the main advantage of using a model including the fatigue variable as inBoldrini 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 inBoldrini 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 free-energy potential and pseudo-potential 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
3 SPATIAL DISCRETIZATION
This section presents the finite element approximation of Eq. (1) . For this purpose, each equation is multiplied by suitable test functions
The boundary terms that appear in the above weak form are denoted by
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
where
where
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
where the derivative matrices are defined as
Replacing the above approximations in Eq. (5) , we obtain the semi-discrete form for the
with the operators defined by
Note that the integrand in Eq. (26) results in a second-order tensor which is expressed as a vector using Voigt notation. Also, the potentials
4 TIME INTEGRATION SCHEMES
This section presents different time integration methods used by the semi and the fully-implicit procedures to solve Eq. (18) . In the semi-implicit 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 2-stage second-order accurate diagonally implicit Runge-Kutta method (TR-BDF2). 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 fully-implicit 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
The integration procedures used in Eq. (18) are presented in the following sections. The semi and fully-implicit algorithms are summarized at the end of their respective Sections.
4.1 Semi-implicit Scheme
This section presents the time integration procedures for Eq. (18) used in the semi-implicit 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
where
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:
This work also proposes the use of the 2-stage second-order accurate diagonally implicit method, or simply TR-BDF2 method. As the backward Euler, the TR-BDF2 is a L-stable method, indicated for the time stiff problems ( LeVeque, 2007 ). The first stage consists of the trapezoidal method applied over
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 semi-implicit scheme. For example, in Eq. (32) is used
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
In order to approximate the solution for the displacement
where the coefficients
Substituting the above approximations in Eq. (36) , we obtain the following time-marching system:
After computing the displacement, the acceleration and velocity fields are updated using Eqs. (37) and (38) . The imposition of prescribed displacement
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
and
Note that both methods have basically the same computational cost, since
The complete algorithm for the solution of the system is summarized in Algorithm 1 .
Algorithm 1 Algorithm of the semi-implicit time integration scheme.
Semi-implicit time integration scheme |
---|
1: for
|
2: Given
|
3: Given the updated damage
|
4: From the current displacement
|
Eqs. (37) and (38), respectively; 5: Update the prescribed velocity and acceleration using Eq. (41); |
6: Given
|
7: Update the time step. |
8: end |
4.2 Fully-implicit Scheme
This section presents the fully-implicit 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
This system is solved employing some iterative procedure instead of a direct method. Supposing that
Dropping the higher order terms and setting
giving the following update step
where
The Jacobian matrix
Newton’s method is summarized in Algorithm 2 .
Algorithm 2 Algorithm of the Newton’s method.
Fully-implicit time integration scheme |
---|
1: Specify an initial guess
2: while tol < error do |
3: Evaluate the Jacobian
|
4: Solve the system
|
5: Update the variable
6: Update the iteration 7: Evaluate the error as 8: end |
9:
|
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
whose components are
The Jacobian matrix can be written in the following block matrix form:
with coefficients defined by
For the
where
The Jacobian matrix is non-symmetric, since
5 NUMERICAL RESULTS
This section presents the numerical results obtained using the semi and fully-implicit time integration schemes. The problem consists of the I-shaped 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

Figure 1 I-shaped specimen geometry and boundary conditions. Attention to detail A, where there is a reduction of 2.5% in the original width.
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
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
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

Figure 2 Damage evolution in the specimen under monotonic displacement. Results along the time: (a)

Figure 3 (a) Damage and fatigue. (b) Vertical stress. Results along the time at the central node of the specimen.
This problem is solved using different time integration schemes. It was verified that the use of any Runge-Kutta explicit method results in numerical instability even for small increments of time as
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
where the reference values
Table 1 Convergence rates for different combinations of time integration methods and the related CPU time obtained in
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 | ||
TR-BDF2 | 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 | ||
TR-BDF2 | 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 semi-implicit scheme, the convergence rates are below the expected theoretical values, even with the use of second-order methods. Arising from the stiff behavior, the use of the L-stable TR-BDF2 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

Figure 4 Error by using: (a) Trapezoidal and TR-BDF2 in the semi-implicit. (b) Fully-implicit. Results along the time increment in log scale.
The full scheme shows quadratic rate, as expected due the use of second-order accurate methods. The error values illustrated in Figure 4 b are below
5.2 Case B
In this case, the sinusoidal displacement given in Eq. (66) is applied along the time range

Figure 5 Damage evolution in the specimen under sinusoidal displacement. Results along the time: (a)
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
6 CONCLUSIONS
This work proposed the use of semi and fully-implicit schemes to perform the time integration of the phase field equations of damage and fatigue presented in Boldrini et al. (2016) . The explicit Runge-Kutta methods showed to be unfeasible for 2D examples, presenting instability with time increments as small as
-
The Jacobian of Newton’s method is hard to obtain and becomes particular to the problem;
The Jacobian matrix is non-symmetric;
There is a high computational cost and memory demand, since it is necessary to solve a large system in each Newton’s iteration;
The semi-implicit 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 TR-BDF2 method showed better accuracy with no significant additional cost. This method is L-stable 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 semi-implicit scheme were below ranging close to the linear rate in most of the equations. Accurate responses (error below
Although not applied in this work, an iteration of correction in the semi-implicit procedure would reduce errors. In particular, the prediction/correction Adams-Bashforth/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 I-shaped 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.