Comparison of semi and fully-implicit time integration schemes applied to a damage and fatigue phase field model

In this work, we apply semi and fully-implicit 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 semi-implicit scheme, each equation is solved separately by suited implicit method. The Newton’s method is used to linearize the equations in the fully-implicit 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 semi-implicit scheme showed be lower than the fully counterpart.


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 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 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  is derived by finding a differential constitutive relation for  that satisfies the entropy inequality for all possible admissible processes.
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  [0,1] are in fact used during the evolution.
However, another reason for given the complete definition of the potentials is that, when numerically solving Eq. (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 micro-cracks that are formed by cyclic loadings, even small ones. The micro-crack 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 From the choice of the potential f  in Eq. (3b) and F as above, the evolution equation of the fatigue in Eq.
(1c) assume 0    in any point of the material. In particular if 1   , the fatigue  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  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 in Eq. (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 free-energy 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 that Eq. (1) does not guarantee 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.

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): (1 ) 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 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 j N as where e N , e  N , e N  , ˆe u , ˆe u  , ˆe  and ˆe  are defined by 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 ˆˆ, e e e e e e e e e u v  where the derivative matrices are defined as 1, 2 Replacing the above approximations in Eq. (5), we obtain the semi-discrete form for the e -th element as 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 ( ) that appear in the damage and fatigue equations are given in 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.

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 [0, ] T in discrete time steps n t with time increments given by 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.

Semi-implicit Scheme
This section presents the time integration procedures for Eq. (18) used in the semi-implicit scheme.

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 I is an appropriate identity matrix. In a very similar way, now considering the right-hand side in the current time step, we have the implicit backward Euler. In this case, the time-marching rule is given by the following system to be solved at each time step 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 t  /2 obtaining the solution *  .
In the second stage, the solution is updated by solving the following system: 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) The same idea is used for all the above methods.

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 1 n  u , 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 where the coefficients j  ( 1, , 6) j   are given in terms of the Newmark coefficients   and   (we assume the Substituting the above approximations in Eq. (36), we obtain the following time-marching system: where the bar symbol represents the prescribed degrees of freedom.

Fatigue Equation
For the fatigue equation Eq. (18c), we use the backward Euler and trapezoidal methods whose time evolution expressions are given respectively by Note that both methods have basically the same computational cost, since ,

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 ( ) i v is the approximation to v at iteration i , Newton's method estimates the solution at iteration i via the Taylor series expansion Dropping the higher order terms and setting giving the following update step where ( ) i v solves the linear system The Jacobian matrix J and its coefficients are respectively defined as Newton's method is summarized in Algorithm 2.
Algorithm 2: Algorithm of the Newton's method. 4: Solve the system 7: Evaluate the error as 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 ] , u n v n n n n n v n n n n a n n The Jacobian matrix can be written in the following block matrix form: .
For the e -th element, the coefficients of the Jacobian matrix are where e B  and 1/ in Eqs. (57) and (60)

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 0.167 sin (10 ) for 0 0.05 ( ) 0.075 sin(20 1.6) 0.092 for 0.05 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.

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 1.4 s t  , 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 Runge-Kutta explicit method results in numerical instability even for small increments of time as . 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 fully-implicit 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 1.5 t s  , given for a generic vector s and the semi-norm by  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 4a. In order to obtain error below 5% for the velocity, a time increment of The full scheme shows quadratic rate, as expected due the use of second-order accurate methods. The error values illustrated in Figure 4b are below 0.5% for all variables. However, the computational cost is approximately 850% higher than the semi-implicit scheme. In this analysis, the Jacobian matrix conditioning number had values of order 17 ( ) 10  J  .

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 I-shaped 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 6a, 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 6b 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.

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 . An alternative to avoid the use of small time steps maintaining stability is to use implicit methods. The fully-implicit method presented quadratic convergence and lower errors when compared to the semi-implicit procedures. However, the fully-implicit scheme has some drawbacks: • The Jacobian of Newton's method is hard to obtain and becomes particular to the problem; This model deals with phenomena from micro to macro-scale, then the resulting iterative system of the Newton's method is poorly conditioned, even with the use of preconditioners; • 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 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 semi-implicit 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 4 5 10 t     . 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 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.