Thermomechanical fatigue life prediction of metallic materials by a gradient‐enhanced viscoplastic damage approach

Fatigue failure plays an important role in engineering applications, especially when structural components experience significant cyclic thermal loading and complex force loading simultaneously. During the last decades, several post‐processing techniques have been developed based on empirical investigations of experimental evidence to predict the fatigue life of materials. The work at hand postulates a conventional continuum damage theory for thermomechanical fatigue failure modeling. In particular, an implicit gradient‐enhanced approach is employed to address the ill‐posedness of the partial differential equation system when the damage onsets. An internal fatigue variable is phenomenologically defined based on the accumulation of viscoplasticity. In the sequel, a regularized fatigue variable is obtained to further yield the damage softening function, which straightforwardly applies to the stress, material tangent, and viscoplastic dissipation. A multi‐field problem, consisting of the strain field, the temperature, and the non‐local variable, is taken into consideration, leading to a fully coupled system. This numerical methodology is consistently derived and implemented into the context of the finite element method. Several representative and demonstrative examples are performed, which yield good numerical stability and agreement with experimental data. Conclusive findings and further perspectives close this article.


INTRODUCTION
crack initiation, stable crack growth is observed according to experimental investigations. Fracture evolution starts from the initiated macroscopic crack and propagates steadily during each loading cycle. During this phase, the crack surface gradually increases until the rest of the force-transferring area is not capable of bearing the loading. Eventually, rupture is observed, which reflects an unstable crack propagation. A number of pioneers have contributed their investigations of fatigue failure to mechanical engineering applications in different fields, referring to the work of Wöhler, 1 Orowan, 2 and Paris 3 for the early literature.
Another categorization classifies fatigue failure into low cycle fatigue and high cycle fatigue. High cycle fatigue is driven only by elastic deformations, whereas low cycle fatigue is governed by elastoplastic properties. With respect to low cycle fatigue, several inelastic constitutive models have been developed. In particular, a viscoplastic model utilizing a nonlinear kinematic hardening law and a nonlinear isotropic hardening law have been proposed by Chaboche. 4,5 For the consideration of time-dependent plasticity in the Chaboche model, see the recent work of Ahmed et al.. 6 In addition to inelastic deformations, metals are often subjected to high-temperature cycles, 7 for instance, in the cases of turbine structures, engine components, and braking discs. This factor leads to more complex considerations with respect to a reliable fatigue life prediction. Hence, thermomechanical fatigue analysis [8][9][10][11][12] has been used to address thermal-cycle-induced fatigue damage.
Several sophisticated numerical models have been developed to study fatigue life criteria, which are mainly based on a parametric lifetime evaluation. Classical parametric formulations 13,14 consist of plastic or plastic-creep strain models, stress-based models, and energy-based models. The damage approximation is characterized as an evolution independent of the increase of the load cycles, which decouples the material response from the damage status. With the aim of a straightforward numerical implementation, the damage evolution is usually characterized by a post-processing technique. Another classification of the fatigue modeling is based on a continuous damage evolution coupled with the deformation response during each cycle. This approach incorporates an internal variable, which is governed by an empirical damage law, to approximate the accumulative damage and to identify the instantaneous nominal stress state throughout the fatigue loading history. Unfortunately, a localized description of the damage formulation has limited application capabilities, especially within a classical finite element method framework. This approach generally suffers from unstable or diverged numerical solutions due to an ill-posed nature of the partial differential equation (PDE) system when the local damage evolution is triggered. 15 In the meantime, it also leads to a pathological mesh sensitivity with respect to the finite element results. Accounting for rate-dependent characteristics, for example, viscoplasticity, the ill-posedness of PDEs can be reduced to some extent in case of strain softening material models. Nevertheless, the viscoplastic regularization cannot always guarantee mesh independent solutions, in particular for fully evolved damage. An earlier study by Needleman 16 has investigated this feature by studying a one-dimensional problem. In the sequel, as explained by the work of Niazi et al., 17 the effective regularization length gradually decreases with damage growth. It eventually gets smaller than the local element size, and the numerical results naturally become mesh-dependent again. Therefore, to address this issue, an additional regularization technique becomes indispensable. The regularization is commonly based on a non-local integral 18 or a gradient expansion 19 approach. The application of the latter one, that is, the gradient-enhanced approach, for modeling high cycle fatigue is presented by Peelings et al.. 20,21 In addition to an implicit gradient approximation of the internal damage variable, a stress-gradient-based model by Tovo and Livieri 22 has been applied to predict the fatigue life of complex structures. The recent variational phase-field method has also been considered for fatigue fracture analysis. The reader is referred to the recently published work [23][24][25][26][27][28] for representative literature.
In the work at hand, a conventional continuum damage formulation is expanded by incorporating an implicit gradient approach into a thermomechanical fatigue model. Considering the significant influence of high-temperature loading, the thermal governing law is derived based on a classical Fourier description. For the sake of simplicity, only the material dissipation due to the established viscoplastic characteristics is assumed to drive the temperature evolution, while the damage dissipation is omitted for the thermal evolution. The behavior of metallic materials is constituted by a rate-dependent plastic model, 5 which utilizes nonlinear kinematic hardening and nonlinear isotropic hardening simultaneously. Furthermore, with respect to the fatigue evaluation, an empirical definition of an internal quantity depending upon the effective deviatoric stress and the plastic multiplier is considered, which is inspired by Sommitsch et al.. 29 This empirical definition is also one of the most important reasons for the omission of damage dissipation for the thermal evolution in this work, since it brings unexpected difficulties to formulate the damage potential with respect to the empirical algorithm and to the further derivations of the fatigue damage evolution. Instead of directly using the internal fatigue variable to evaluate the damage state, this work employs a gradient enhanced regularization based on an implicit mathematical description. Thus, the damage status is eventually identified by an exponential function of the regularized variable or the non-local variable. In this regard, not only the numerical stability is largely improved, but also the issue of mesh-dependency in terms of the finite element discretization is appropriately addressed. In contrast to the classical post-processing evaluation, the present model fully couples three independent fields, namely, the deformation, the temperature, and the non-local variable, and yields a monolithic solution for all of them.
The structure of this work is outlined as follows. In Section 2, the fundamental constitutive description of the thermomechanical fatigue damage law is presented. In particular, the governing law of the temperature evolution, the constitutive Chaboche viscoplasticity, and the gradient-enhanced approach are introduced. Section 3 describes the numerical realization of the return mapping, the stress, and the consistent material tangent. Furthermore, the finite element implementations of the Galerkin-type discretization with respect to the element residual and its linearized tangent are presented. In the sequel, a couple of representative examples are numerically studied in Section 4, and Section 5 closes the article with findings and potential perspectives.

CONSTITUTIVE LAW OF THERMOMECHANICAL FATIGUE DAMAGE
In order to understand and derive the governing laws of the multi-field coupled problem, one can start from the definitions of the global field variables. Three primary fields, namely, the displacement u, the temperature Θ, and the non-local quantity are considered within the continuous solid domain ℬ ∈ R 3 along with the time ∈ R 1 + , which are characterized by respectively. The small strain tensor is defined by the gradient of the displacement field with a symmetric operation, reading For a general decomposition of the total strain, an additive relation of the volumetric and the deviatoric contributions is considered as The fourth order tensors V, P, and I represent the volumetric, deviatoric, and identity projections, respectively. In addition to the volumetric deviatoric decomposition, the total strain tensor can be decomposed in a different manner by the consideration of the viscoplastic tensor and the thermomechanical strain tensor Taking a classical J 2 plasticity model into consideration, the viscoplastic strain tensor is characterized by a deviatoric tensor with the relations Meanwhile, the assumption of homogeneous thermal expansion defines the thermal strain as Based on the definition in Equation (7), one can easily conclude that the thermal strain tensor is purely volumetric, that is, The coefficient th describes the expansion effects and Θ 0 is given as a reference temperature. As a result, the effective stress is obtained as The fourth order tensor C e is the linear elastic tangent, which is governed by the bulk modulus and the shear modulus , such that C e = 3 V + 2 P.
Furthermore, regarding the temperature field, a classical Fourier law is taken into account in this work, representing the heat flow oriented from hot to cold. The spatial heat flux vector is defined as The second order tensor  describes the thermal conductivity, which could be isotropic or anisotropic depending on material characteristics. This work only considers a simple isotropic feature governed by the parameter  0 , namely,  =  0 1.

Temperature evolution
In order to derive the temperature evolution law, one recalls the first law of thermodynamicṡ where the scalar quantity  (x, t) represents the specific heat source. The mass-specific internal energy density  and its time derivative are defined through the partial Legendre transformation of the Gibbs equation as respectively. The quantity is the Helmholtz free energy density function, and is the entropy density. Considering a general material constitutive model involving inelasticity , the time derivative of the free energy density yieldṡ The non-negative local dissipation is presented by the Clausius-Planck inequality, reading which, according to the relation in Equation (11), is alternatively expressed as With further algorithmic manipulation of Equations (12) 2 , (13), and (15), a relation is obtained, which yields the relations see the work of Holzapfel 30 for a detailed explanation. Due to the characterization of and the relations in Equations (13) and (17) 2 , one obtains the rate of entropy expanded aṡ where the heat capacity is defined as As a result, according to Equations (14), (17) 3 , and (19), the spatial description of the heat conduction evolution law readsΘ According to the investigation of Aldakheel and Miehe, 31 the quantity  involving second order mixed derivatives characterizes the latent heating (thermoelastic) effect, while  loc straightforwardly represents the local dissipation due to inelastic constitutive properties. For most metallic materials in the context of thermoplasticity applications, the latent heating effect is considerably smaller than the inelastic dissipation, that is,  ≪  loc . The interested reader is referred to the work of Miehe 32 for a demonstration in entropic thermoelasticity. The work at hand neglects the latent heating effect and only considers the local inelastic dissipation to drive the temperature evolution. Hence, it largely reduces the complexity of the thermomechanical coupled problem. The governing equation for the heat conduction is eventually written asΘ

Classical Chaboche viscoplastic model
This work considers a classical Chaboche model, which decomposes the hardening into kinematic and isotropic components governed by their respective constitutive laws. A conceptual von Mises yield surface function for the consideration of incompressible plasticity is defined by where the hardening variable is denoted by r and the coefficient Q r is a material constant. The deviatoric stress is dev = P ∶ , and two back stress tensors X 1 and X 2 are governed by the nonlinear rule aṡ The model parameters are denoted by C i and i . In particular, a linear kinematic hardening law can be recovered by setting i = 0 to exclude the plastic multiplier and the back stress from the description. For applications of this model, the reader is referred to the works of Chaboche 4 and Armstrong and Frederick, 33 to name a few. This work takes a general superposition algorithm of X i with i = 1, 2 to achieve the modeling flexibility, which follows the work of Chaboche. 5 Regarding the isotropic hardening, it depicts the relatively slow plastic evolution and is governed by the accumulated plastic strain. A simple and classical function of the isotropic hardening variable reads which naturally satisfies the initial condition r = 0 before triggering plasticity (p = 0). As a consequence, this definition yields the isotropic hardening rate aṡr For a further extension, the definition of the hardening rate from Equation (26) is expressed in this work aṡ where b r , g r , and a r are the model parameters for the nonlinear isotropic hardening. The governing functions of the viscoplastic strain rate and the viscoplastic multiplier are defined bẏ respectively. The parameters K and n in Equation (28) are material constants to describe the viscoplastic evolution.
According to the definition of the yield function in Equation (23) and the viscoplastic strain evolution in Equation (28), the constitutive equation of the viscoplastic strain is simplified aṡ characterizes the flow direction of the viscoplastic strain. The plastic flow direction in Equation (29) at an arbitrary loading state is a unit vector normal to the yield surface of conventional J 2 plasticity. The Kuhn-Tucker conditionṡ are the governing criteria for the viscoplastic evolution. In particular, it is notable that this work accounts for small strain kinematics for the constitutive description of metallic materials. Nevertheless, when metallic materials are subjected to excessive external loading, they often undergo large deformations, which requires a finite strain kinematics to constitute the numerical model due to inevitable geometrical nonlinearity. Therefore, in order to obtain reasonable results, the numerical applications within this work are restricted to small strain model problems without excessive deformation.
Furthermore, the viscoplastic yield function for the majority of metallic materials is typically temperature-dependent, since the material parameters are characterized as temperature-dependent. Nevertheless, this work does not consider the temperature dependency of material parameters based on two main reasons. First, the parametric study of temperature dependency requires substantial experimental support, which is not the scope of the current work. Then, it definitely brings considerable complexities to the constitutive derivation for the coupling terms in the consistently linearized tangent for the FE implementation, if one considers a temperature-dependent viscoplastic yield function in a fully coupled thermomechanical model. Therefore, this work does not attempt to take finite kinematics and temperature-dependent parameters into consideration, but keeps these two issues for the next priorities in the future research.

Implicit gradient-enhanced damage approach
Singularities induced by local softening in conventional continuum theories may yield mathematically ill-posed differential equations. 15 As a result, numerical solutions may not converge, which is not reliable from a mathematic point of view. Furthermore, finite element simulations suffer from the issue of pathological mesh sensitivity after the onset of localized damage. To overcome this issue, pioneers developed regularization approaches to avoid singularity issues, namely the non-local integral algorithm 18 and the gradient-enhanced approach. 19 The work at hand incorporates the latter one, a classical gradient-enhanced damage model, to simulate thermomechanical fatigue failure. Considering an internal length effect parameter c, the governing relation between the local internal variable and its non-local representation inside the continuum domain Ω reads and the Neumann boundary condition of the non-local variable at the domain surface Ω is given as where is a normal vector outward at the domain surface. The thermodynamically consistent derivation of the non-local governing equation in Equations (31) and (32) is already an established knowledge, which is not presented in this work. The reader is referred to the works of, for example, Forest, 34 Dimitrijevic and Hackl, 35 and Poh and Sun, 36 for further interests and detailed insights. Following the empirical rule, 29 the internal quantity is constituted as It is necessary to point out that Sommitsch et al. 29 use a pure post-processing technique to evaluate the fatigue life, and the variable is directly used to evaluate the damage status of the materials. In contrast, this work extends it to a continuum damage model to simulate thermomechanical fatigue failure. The damage status is evaluated depending on another regularization function. The normalization constant of the inelastic multiplier is denoted bẏ0. The stress dependency of the lifetime behavior is governed by 0 and ℳ 0 . The power index 0 represents the time dependence of the lifetime. According to the statement by Sommitsch et al., 29 0 = 1 yields rate independent behavior, and 0 = 0 leads to a fully time dependent lifetime behavior. In general, 0 is identified to be a positive but rather small value ( 0 ≪ 1) for high temperature loading. Another simplification in this work is that these material parameters are temperature-independent, namely 0 , ℳ 0 , 0 , anḋ0 are always constant for the subsequent formulations and simulations. This significantly simplifies the consistent derivation and parameter identification. Furthermore, the equivalent stress eq in Equation (33) is a scalar defined as As a result, without considering direct temperature influence, the internal quantity is driven only by the equivalent stress eq and the viscoplastic multiplierṗ.
As mentioned previously, instead of using in Equation (33) to directly evaluate the damage status, a softening function is defined to depend on the non-local variable . In this work, a simple exponential function is considered through for the continuum damage evolution. This damage function depicts the softening that starts when the non-local variable reaches the threshold cr . The coefficient  describes the rate at which the damage grows. The operator ⟨ * ⟩ + = ( * + | * |) ∕2 is known as Macauley bracket. The model presented within this manuscript aims at simulating low-cycle thermomechanical fatigue. It follows the established approaches in this field, namely, the use of a viscoplastic model accounting for nonlinear kinematic and isotropic hardening, and a fatigue law accounting for the plastic variables and the stresses. In several existing viscoplastic-damage approaches with mixed hardening laws, the viscoplastic model is used separately to calculate the plastic strains for certain cycles, and, then, the post-processing step for fatigue calculation is subsequently performed with the assumption of linear fatigue accumulation. This depicts an oversimplified interdependence between fatigue damage and viscoplasticity evolution. So far, very limited publications have considered or implemented such an approach that combines these two formulations into one comprehensive constitutive model. The main scope of this work attempts to introduce a comprehensive model with thermomechanical coupling, complex viscoplasticity, and gradient regularization, which allows the evaluation of fatigue damage under random temperature and stress variations. Meanwhile, it is capable of describing the evolution of fatigue damage in a straightforward manner, without presupposing a linear accumulation in further post-processing calculations.

Return mapping algorithm
This section presents the detailed algorithm to obtain the material stress and tangent using a return mapping scheme. First of all, the yield surface needs to be evaluated for the sake of determining the onset of viscoplasticity. The condition  ≤ 0 yields a pure elastic response and  > 0 triggers viscoplasticity. With respect to the inelastic response, the goal is to find the solution for the evolved viscoplastic strain and the accumulated internal hardening. Due to the complex nonlinearity of the yield function, obtaining an analytical solution does not seem to be always possible. Therefore, based on a numerical approximation, the classical radial return mapping is a powerful tool to solve the internal variables. 37 An initial guess defines the trial isotropic hardening, the trial back stress, and the trial deviatoric stress at the current step t n+1 using the values at previous loading state t n , namely, Therefore, the trial yield surface reads If the condition  tr t n+1 ≤ 0 is fulfilled, the trial assumption is understood as the real state of the material. Hence, no additional effort is required and the variables are updated using the trial state in Equation (36) as the history variables for the next loading state. On the other hand, for the case of  tr t n+1 > 0, a common approach to obtain the solution is applying a local Newton procedure. Within a pseudo-time increment [t n , t n+1 ], the internal variables are implicitly discretized as The back stress tensor is expressed as and the isotropic hardening yields Then, Equation (39) leads to the relation where The Kuhn-Tucker conditions in Equation (30) are rewritten as The deviatoric stress at the current step is dev With further algorithmic manipulations, the explicit J 2 quantity at the current loading step is expanded by trial quantities, reading In the meantime, the normalized second order tensor for the viscoplastic flow direction, which is defined in Equation (29), can be equivalently expressed as The proof steps of the relations in Equations (45) and (46) are presented in the Appendix. As a result, incorporating the aforementioned expression in Equation (45) into Equation (23), the yield function is rewritten as To obtain the numerical solution for the unknowns Δp and Δr, one can consider Equations (38) and (40) as two coupled equations . In the sequel, the local residuals are defined by Their consistent tangent is derived by a straightforward linearization, reading The variable a 2 is defined by the underbrace in Equation (49) and the derivative term  ∕ Δp according to Equation (47) is expanded as In the sequel, an internal Newton-Raphson scheme, shown in Table 1, can be conducted to obtain the unknown Δp and Δr at the current loading step t n+1 . The viscoplastic strain, the back stress terms, and the isotropic hardening variable are updated accordingly.

Material stress and tangent
The deviatoric stress is eventually obtained by Equation (44) as long as the plastic multiplier is resolved. Then, the consistent deviatoric material tangent reads

TA B L E 1 Local Newton-Raphson iteration scheme
Initiation According to the relation in Equation (46), a straightforward partial derivative of with respect to yields Herein, the only unknown quantity in both Equations (51) and (52) remains the term Δp∕ , which has been consistently derived in Appendix B. Thereby, the consistent material tangent for the deviatoric part is obtained by substituting Δp∕ into Equation (52), and further in Equation (51). Consequently, the total stress and the consistent material tangent read respectively.

Finite element residual and linearized tangent
This subsection depicts the numerical implementation within the conventional finite element context. Naturally, instead of using the stress obtained in Equation (53), the actual stress tensor̃needs to be considered for the strong forms of the strain field and temperature field, where the actual stress̃and the viscoplastic dissipation rate are rewritten as respectively. Regarding the term loc considered in Equation (54) 2 , only the established viscoplasticity dissipation is considered for the sake of simplicity, even though other internal variables, for example, d, X 1 , X 2 , r, , may yield a portion of inelastic dissipation. The reason behind it is that the explicit dissipation rate terms associated with these internal quantities cannot be straightforwardly derived due to the complexity of the phenomenological and empirical definition within the present constitutive model. This simplification is also inspired by the work of Dittmann et al., 38 who define the total dissipation by artificially scaling the plasticity, fracture, among other inelastic dissipations. As a result, the coupled multi-field problem is governed by the strong forms of for the strain field,Θ for the temperature field and for the evolution of the non-local variable. It is noteworthy that the surface traction, the body force contribution, and the Neumann boundary conditions are not presented for the sake of simplicity. In the sequel, according to a straightforward manipulation, the weak forms of these three fields are obtained by multiplying with their virtual quantities, and are further expanded as for the strain field, for the temperature evolution and for the non-local variable evolution. Neglecting the surface integration terms in Equations (58)-(60), which depict the Neumann boundary condition for each field quantity, the volumetric integrations can be straightforwardly implemented by considering the isoparametric element approach. The values u, u, Θ, Θ, and , are interpolated using shape functions N I u , N I Θ , N I and the nodal values u I , u I , Θ I , Θ I and I , I , that is, where the index n represents the number of nodes in one element. It is noteworthy that the interpolations of the three different field variables are independent of each other. One can use identical or nonidentical shape functions based on the necessity of continuity. This work only considers a simple linear interpolation function, namely a 4-node isoparametric quadratic element and an 8-node isoparametric brick element for two-dimensional problems and three-dimensional problems, respectively. As a result, only 2 Gauss integration points along each spatial direction are necessary to be taken into account for the volume integral in an element. It is also possible to apply selective integration algorithms, which is nevertheless not the main scope of this work. Meanwhile, the interpolation of the three aforementioned fields is based on the same linear shape functions, that is, N I u = N I Θ = N I = N I with I = 1, … , n. The residuals are eventually obtained as Their consistent tangents are obtained by a straightforward linearization, reading It is noteworthy that all the rate terms are numerically approximated by an implicit backward Euler algorithm, for example, loc = 1 Δt̃∶ Δ p and = t n + Δt The derivatives shown in Equation (63) are given in Appendix C. Having the implementation of the sub-matrices for each node, the element integration of the total tangent and residual reads respectively. By means of assembling these sub-matrices over all elements within the global domain, the global tangent and residual are obtained. The formulation is implemented by Fortran coding into an in-house platform and a user interface (UI) of Ansys APDL. Furthermore, it is apparent that the derivation of the coupling components in Equation (63) shows non-symmetry. Therefore, an unsymmetric solver must be employed in order to guarantee accuracy of the simulation outcomes.

A cube with uniaxial loading
Before applying the proposed thermomechanical fatigue damage model to structural simulations, a simple material test is first carried out to verify the constitutive response. In general, an individual cubic element with the size of 100 mm subjected to uniaxial strain loading is taken into account. Basically, two classes of strain applications are considered in order to investigate the different responses. A published work 39 has also studied such a model problem by showing corresponding experimental data. In this work, for the sake of simplicity in regard to the material parameter identification, nonlinear kinematic hardening is excluded, and only nonlinear isotropic hardening evolution is taken into account. The material parameters are given in Table 2.
The first class of strain application, see Figure 1A, studies the material relaxation. Here, a tensile strain rapidly increases to 2% and stays constant for 3600 s. This simulation neglects damage ( = 0), which intends to verify the correct constitutive material law and identify its parameters. Analyzing the simulation results, the stress-time relationship shows good agreement with the experimental results of Schicker et al., 39 see Figure 1B.
The second application is cyclic loading, which comprises 25 load cycles in total. Each cycle consists of rapid tensile and compressive loading, followed by material relaxation at the compressed state for 180 s. As representatives, the first three cycles are depicted in Figure 2A. The resultant stress-strain curve can be seen in Figure 2B, which shows a relatively fast stabilization after a couple of load cycles. This modeling does not consider the damage evolution yet. In particular, the simulation provides a satisfactory result of the stress-strain relationship compared to the experimental data with respect to the first and the third cycle, see Figure 2C,D, respectively.
In the sequel, damage-like softening induced by is triggered when the condition  > 0 is fulfilled, yielding a gradually decreasing material stress along with increasing loading cycles. Considering a zero threshold for the non-local variable, that is, cr = 0, the stress-strain relationships for the cases with  = 0.01 and  = 0.02 are shown in Figure 3A,B, respectively. By comparing the results in Figure 3A,B, one can deduce that for a larger factor  damage grows faster and the fully damaged state is achieved earlier. By setting a positive threshold value cr > 0, the initiation of the damage evolution is postponed. A larger value of cr triggers the damage softening comparatively later. As a result, less damage has evolved after 25 loading cycles. Meanwhile, a larger portion of the effective stress remains after the final (25th) cycle, see Figure 3C,D. A published work of Da Costa Mattos 40 also models such cyclic stress-stain relationship coupled with damage evolution, without taking the relaxation into account. In general, a similar profile of the actual stress with respect to the strain is obtained in Figure 3 compared to the result of Da Costa Mattos. 40 According to the examination of the two important parameters  and cr of the softening function in Equation (35), expected simulation results are obtained. Hence, the present model is capable of capturing an expected fatigue damage behavior for material failure.
Furthermore, a numerical investigation on the solution convergence based on the Newton-Raphson method is presented. In this work, the main numerical convergence criterion is to check if the coupled residual norm is less than a predefined tolerance. In particular, taking the simulation of Figure 2C as one example, the first 55 loading steps consist of a monotonic compressive load and a compressive relaxation, which require a total of 180 Newton-Raphson iterations. The logarithmic value of the global residual norm, that is, log 10 ||R|| = log 10 , is shown  Figure 4A. Evidently, for both the compressive loading phase and the compressive relaxation phase, quadratic convergence behavior is obtained. In addition, an internal iteration is also employed for the return mapping algorithm based on the Newton-Raphson approach. To investigate the internal convergence behavior, four arbitrary points (AP) from the global iterations are selected as representative examples. The logarithmic value of the local residual norm, that is, log 10 ||ℛ|| = log 10

TA B L E 2 Material parameters for the cubic element test
, versus the local iteration step at one certain integration point is shown in Figure 4B.
The quadratic convergence behavior is again observed. According to the global and local convergence study, one can deduce that the numerical derivation consistency and the implementation robustness of the fully coupled model are successfully verified. It is necessary to clarify that the residual norm shown in the aforementioned figures is based on its absolute value, even if the convergence criterion of the Newton solution is based on the approach of checking the relative error of the residual norm. The purpose is to distinguish between the compressive loading phase and the relaxation phase, since their absolute residual values are distinctive and can be easily identified. Another issue is the normalization of the global residual based on three components, namely, temperature, displacement, and non-local fields. They are characterized to be of different orders as well as distinctive units. However, regarding a monolithic solution of this fully coupled system, the mechanical problem is eventually converted to resolve the non-linear mathematical equations. For the sake of simplicity, a simple global residual from the mathematical description is examined in this example, since a converged solution for the global nonlinear equations definitely yields simultaneously converged solutions for the three different evolution responses. This mechanism also applies to the local Newton solution, for example, ℛ p and ℛ r .

Fatigue failure of a dumbbell structure
As already known, compared to conventional local damage modeling, the present gradient-enhanced damage formulation yields a damage profile with spatial regularization, which is characterized as mesh insensitive. Furthermore, the numerical stability to resolve the coupled PDEs, even in the material softening phase, is sufficiently enhanced. In this example, the advantages of the continuum damage model with a gradient-enhanced algorithm versus a conventional one is elaborated by a cyclic tension test on a dumbbell-shaped specimen. The model problem in this example does not follow any conducted experiment, but for similar geometry and loading setup, the reader is referred to existing literature, for example, Broggiato et al.. 41 Herein, considering a three-dimensional boundary value problem (BVP), the bottom surface of the finite element model is fully fixed. The top surface is subjected to a uniaxial displacement loading u consisting of five cycles, see Figure 5. The amplitude of the load cycles increases stepwise for the first three cycles, while the rest of the cycles have the same amplitude as the third one. During this process, the loading or unloading rate stays constant. In order to effectively trigger a localized damage initiation, the material with a length of 10 mm in the middle part is artificially weakened by assigning a value of 0 that is ten times that of the unweakened material. For the purpose of a comparative study, the geometry is discretized by four different meshes, namely, 4032 elements (Mesh 1), 5712 elements (Mesh 2), 9072 elements (Mesh 3), and 15,792 elements (Mesh 4). The mesh refinement is only considered in the cylinder part of the specimen. Meanwhile, several length effect parameters are considered for the damage regularization investigation. The material parameters are given in Table 3. As one illustrative simulation for Mesh 3 with c = 20 mm 2 , the full damage initiates at the center of specimen, and subsequently evolves outwards. Based upon a post-processing technique, one generates an isosurface of the damage variable d = 0.99 within the specimen, shown in Figure 6A. For each cycle, the damage variable d along the central axis of the structure is depicted in Figure 6B, showing a stepwise increment in the weakened region. First of all, the local damage model, that is, c = 0 mm 2 , is considered based on the aforementioned four FE discretizations. Not surprisingly, the numerical stability is lost at an early stage during the onset of the localized damage evolution. Regarding the simulations of all four different meshes, the damage profiles before the fully diverged solution are shown in Figure 7A-D. The pre-weakened elements are 2 element layers for Mesh 1, 4 element layers for Mesh 2, 8 element layers for Mesh 3, and 16 element layers for Mesh 4. After the simulation, the fully damaged elements (d max ≥ 0.99) are 2 layers for Mesh 1, 2 layers for Mesh 2, 6 layers for Mesh 3, and 14 layers for Mesh 4. The fact of mesh-dependent characteristics of a continuum local damage model is evidently observed, for example, the damage localization is strongly affected by element refinement. This can be observed from another perspective, namely, evaluating the damage distribution along the central axis. Figure 7E shows that a comparatively localized damage zone is Fatigue damage distribution based on a gradient-enhanced formulation with c = 10 mm 2 concentrated on the pre-weakened region, and Figure 7F depicts the magnified damage distribution in the pre-weakened region. As a consequence, one can deduce that the local damage does not exhibit a strict spatial convergence, since a convincing spatial convergence should yield a gradual approaching tendency from the coarsest Mesh 1 to finest Mesh 4. The magnified plot in Figure 7F, unfortunately, does not show such characteristics for local damage. Furthermore, some unexpected local oscillations of the damage variable from the fully broken zone towards its neighboring region are observed as well. Meanwhile, the potential reasons, why the localized damage plateau does not become narrower as a denser mesh is applied, are the influence of the pre-defined weakening zone band and the different viscosity regularization effects that occur before evolving full damage. Regarding a gradient-enhanced damage model, for example, c = 10 mm 2 , the issue of mesh sensitivity for the damage evolution is addressed to a great extent. Based on the same simulation setup, the results of non-local damage modeling yield similar damage profiles for different meshes, see Figure 8. In fact, the magnified plot in Figure 8F shows that the damage profiles for Mesh 2, 3, and 4 are almost identical to each other. Furthermore, the numerical stability is largely enhanced, and all the solutions converge quadratically even at fully damaged states.
In addition, the influence of the length effect parameter on the damage response is also investigated. Using the same FE discretization, for example, Mesh 2, four different length effect parameters are adopted to compare the results, that is, c = [0, 10, 20, 40] mm 2 corresponding to Figure 9A-D, respectively. In general, a larger c leads to a broader and a more diffusive damaged zone, also see Figure 9E. Compared to the local damage result (c = 0 mm 2 ), the simulations with regularization (c ≠ 0 mm 2 ) are capable of enduring more deformation with a smooth necking effect. Another interesting finding is that the damage zone still becomes broadened with further loadings, even if the pre-weakened zone is fully damaged. Due to the mathematical characteristics of the present damage model, a heuristic (1 − d) type continuum damage theory does not yield a vanishing derivative of the softening function when the material integrity is completely lost, for example, d ≈ 1, which naturally results in a broadening and spreading damage zone. This phenomenon has been explained in details in the work of De Borst and Verhoosel. 42 Another particular interest is the relationship between the resultant force and the displacement with respect to the local and gradient-enhanced damage models. Taking three discretizations (Mesh 2, 3, and 4) as illustrations, the cyclic loading-displacement curves are depicted in Figure 10A,B for local (c = 0 mm 2 ) and non-local modeling (c = 40 mm 2 ), respectively. For the local case, the first and second cycles do not yield an obvious deviation for all meshes, since the viscoplastic regularization plays an important role before full damage occurring (d ≤ 0.9). Starting from the third one onward, the further accumulated provides a sufficient driving force to yield full damage. Thereafter, the viscoplastic regularization effect vanishes and the growth of local damage results in an obvious mesh dependency in Figure 10A. Additionally, all simulations consequently fail to converge stably from the numerical point of view. The solid points in Figure 10A identify the loading states where the solutions fail. For detailed insights, the reader is referred to Figure 11 for the converged or diverged behavior of the global algebraic solution. With respect to the gradient damage case, the load-displacement curves, shown in Figure 10B, do not show apparent deviations for each cycle. Even for the fourth and the fifth cycles, where the damage has been fully developed, the differences are still negligible. Moreover, the numerical stability for resolving the global algebraic equations does not show any issues or difficulties, yielding a satisfactory convergence behavior.
Furthermore, the fully coupled thermomechanical formulation also leads to temperature evolution due to inelastic dissipation. Inspired by Aldakheel, 43 who formulates a gradient plasticity model coupled to the thermomechanical problem, this example also investigates the temperature evolution along with the cyclic loading history. As expected, the temperature, shown in Figure 12A-D, starts to rise in the weakened material, where the accumulative plasticity concentrates. The temperature distribution along the axial direction, illustratively depicted in Figure 12E, also outlines this finding. In regard to the simulation results in Figure 12, from the phenomenological point of view, physical characteristics of the temperature increase are captured appropriately, which is in a general agreement with the simulations of Aldakheel. 43

A specimen with a hole under random loading
The third example is inspired by Rosa et al., 44 who concentrate on a strain-life approach for the thermomechanical fatigue lifetime prediction using a conceptual post-processing technique. Herein, a structure, which is similar to the modified ASTM International Specimen, 44 is analyzed by adopting the fully coupled thermomechanical fatigue damage model. In this regard, the geometry and boundary setup are depicted in Figure 13A. The three-dimensional BVP is numerically discretized by 17,520 hexahedral elements, where a relatively fine mesh is applied in the central hole region. The material parameters in Table 4 are considered for the subsequent analyses. Three classes of load are taken into account: pure cyclic displacement loading, pure cyclic temperature loading, and random cyclic mixed displacement-temperature loading.
The first analysis investigates cyclic loading governed by displacement application with constant amplitude, which includes both tensile and compressive loads. The period of each cycle is 40 s, and the simulation consists of 250 cycles in total. The loading amplitude for each cycle is schematically shown in Figure 13B. As an interesting information with respect to the structural response, the reaction force monitored at the top surface is plotted versus the loading time in Figure 14. One can clearly observe in Figure 14A that the envelope of the cyclic reaction force smoothly decreases along with the increase of the number of cycles. The underlying mechanism is that the accumulated fatigue-damage evolution leads to a gradually decreasing nominal stress for the whole domain. For a detailed investigation, three distinct sections are presented to show a magnified load-time relationship, see Figure 14B-D. Damage evolution of the specimen is shown in Figure 15 with respect to the 50th, 125th, 200th, and 250th cycle. The damage is nucleated from the vicinity of the central  Figure 16 is considered as a representative location for damage nucleation, for which the damage evolution as well as the temperature evolution over the entire loading period is shown in Figure 17. As a result, both damage and temperature increase monotonically at Point A. The second analysis studies cyclic loading governed by temperature application with constant amplitude. It is noteworthy that not the entire solid domain is directly affected by temperature loading. In contrast, only the central region highlighted in Figure 13A is constrained by temperature Dirichlet boundary condition, see Figure 13C. Meanwhile, the upper and lower displacement boundaries are fully fixed. As a consequence, the thermo-effect induced expansion and shrinking deformation triggers plastic accumulation, and eventually leads to fatigue damage. As one can see in Figure 18, the potential damage may take place in the central area, from two external edges, but not from the internal hole. The reason behind this is that the lateral expansion is larger at the two edges in comparison to the deformation of the hole-region. Furthermore, the deformations of the upper and lower boundary regions are also indirectly influenced by thermal conduction from the heated zone. Due to the displacement constraint, a portion of plastic accumulation also concentrates in these regions, which naturally yields slight damage. All in all, the cyclic temperature loading with an amplitude ΔΘ = 300 K does not involve considerable damage at all for the whole loading period, that is, d max ≤ 0.024. The reaction forces versus the loading time do not show an apparent decrease of the envelope, which is not presented here.
For the third analysis, a mixed loading pattern, which consists of both displacement and temperature loadings, is considered by means of simultaneous random excitations, see Figure 13D. Similar to the case with pure displacement cyclic loading, one observes in Figure 19 the gradual decrease of the reaction force amplitude with respect to time, and in Figure 20 the damage profiles during the fatigue loading.

Thermomechanical fatigue analysis of a turbo charger
This numerical example studies a representative engineering component, a turbo charger, whose fatigue life estimation is strongly affected by cyclic temperature loading. A complete turbo charger device comprises several components, and is assembled by various complex mechanical joints. A delicate simulation definitely requires incredible computational efforts considering different aspects. For the purpose of an academic insight, a simplified model problem is taken into account, which does not correspond to any commercial product in terms of geometrical dimensions and loading conditions. The model geometry is generated on a CAD platform, which is subsequently discretized by ANSYS Workbench. A further pre-processing manipulation imports the discretized model to ANSYS APDL. Thereafter, the simulation is performed using the UI via ANSYS APDL platform, which consists of the present thermomechanical fatigue damage modeling.
The FE discretization of the three-dimensional BVP is shown in Figure 21A, which consists of 15,453 hexahedral or prismatic elements. Apparently, this structural component does not have any symmetric characteristics. The inner surface of the central axial hollow and the four boltholes are fully fixed, see Figure 21B. The turbo charger is assumed to be subjected to a homogeneous thermal loading. In other words, the same temperature amplitude ΔΘ = 500 K is applied to all integration points of the entire FE model of the turbo charger. As an alternative and equivalent setup, all nodes are fixed for the temperature degree of freedom, to which the same temperature amplitude is applied. In this regard, the thermal evolution is not taken into account. In total, 200 cycles are performed with a loading/unloading rate ofΘ = 500 K/min. The model parameters are given in Table 5. The simulation results in Figure 22 depict the temperature induced fatigue damage evolution for several representative cycles. Although the magnitude of the damage variable is considerably small (d max = 3.43 e −4 versus full failure state d = 1) within only 200 loading cycles, the overall distribution of damage plays a meaningful role in the predictive estimation where the fatigue failure potentially initiates. As expected, some of the potential damage zones are around the fully fixed boundaries, namely, the boltholes on the upper plate and the central axial hollow region. Analyzing Figure 22C, damage starts to gradually concentrate in the fully fixed boundary locations from the 80th cycle, and in the sequel, damage keeps increasing in these regions.

F I G U R E 22
Damage profile for 6 representative cycles As long as the loading cycles are sufficient, fatigue damage will localize in the aforementioned regions, which eventually yields structural failure. Many experimental evidences have shown that the joint regions endure more failure risks. The failure mechanisms are complex, which include the overloaded stress concentration, fatigue induced accumulative defects, and frequent friction due to relative motion of the joint pairs and so forth. Ductile phase-field fracture modeling by Borden et al. 45 simulates a bolted steel plate subjected to a blast load. According to their work, although the main fracture region is in the center part, where the external loading is applied in a concentrated manner, partial damage also exists in the bolt regions. Another potential damage zone is the joint region between the circular pipe and the straight one, see the Region B in Figure 22F. This connective joint leads to a sharp angle of the material geometry, which naturally yields stress concentration and plastic strain accumulation. For a similar investigation, the reader is referred to a practical simulation of a turbo charger component, 39 which, nevertheless, is based on a conventional post-processing algorithm. Therein, the geometry and the cyclic thermal loading of the structure are realistically defined. Supported by experimental validation, the result successfully captures the potential location of fatigue damage nucleation induced by cyclic thermal loading, which is also the transferring zone from the straight pipe to the circular one.

SUMMARY
The work at hand constitutes a thermomechanical fatigue damage model incorporating an implicit gradient-enhanced approach. In particular, three distinct fields are fully coupled, resulting in a globally monolithic equation with respect to temperature, deformation, and a non-local fatigue variable. The temperature evolution is derived out of the first and second law of thermodynamics. The governing equation is simplified by excluding the thermoelastic effect for inelastic materials. For the sake of further simplicity, only plastic dissipation is assumed to be taken into account as the temperature driving term. The constitutive behavior follows a classical Chaboche model, which considers nonlinear kinematic and isotropic hardening in the viscoplastic response. A conventional return mapping scheme is adopted for the solution for the internal plastic variable. Furthermore, an established gradient-enhanced approach is employed to regularize the internal fatigue variable, which is empirically defined based on the accumulative viscoplastic deformation. Thereafter, an exponential softening function based on the regularized field variable degrades the stress and material tangent, to model phenomenologically damage-like behavior. Comprehensively, the fully coupled model is derived in a consistent manner, and is numerically implemented into an in-house platform and ANSYS APDL (via a UI) using standard FE methodology. Evaluating the simulation results of different model problems, the present work provides a satisfactory prediction capability of fatigue damage problems. On the one hand, several existing publications decouple the material deformation and damage into independent evolution laws for the purpose of fatigue life estimation based on an assumption of linear fatigue accumulation using post-processing techniques. This depicts an oversimplified interdependence among the temperature, fatigue damage, and viscoplasticity responses. A scientific increment of this manuscript is that a continuum damage modeling, which simultaneously resolves the temperature response, material deformation and the fatigue damage evolution, is developed. The present model is capable of obtaining an actual and instantaneous response for three distinctive fields in a straightforward manner. On the other hand, compared to a conventional localized damage theory, the gradient-enhanced approach successfully overcomes the issue of damage-softening-induced ill-posed equation. The issues of numerical instability and mesh sensitivity are improved to a large extent, and a couple of illustrative examples are studied in detail for insights into the numerical convergence behavior using an iterative Newton-Raphson solver for nonlinear algebra. Both the global multi-field coupled problem and the local return mapping show a quadratic convergence behavior, which has verified the consistent derivation and implementation. Furthermore, a Fourier-type temperature evolution law is formulated within this work, and inelastic dissipation drives thermal evolution. The thermal coupling allows the evaluation of fatigue damage under random temperature variations in addition to a stress loading. In general, the present methodology plays a meaningful role in solving practical engineering problems.
In general, the present model is capable of capturing the fundamental characteristics of thermomechanical fatigue damage failure. However, many open and interesting questions certainly still exist. One of the next priorities for future research is an extension towards finite strain kinematics to constitute the material model, which allows to predict the fatigue failure at a relatively large deformation. Then, more studies on parameters, especially those accounting for temperature-dependent characteristics, are of importance for practical and realistic engineering applications. Furthermore, in order to obtain a relatively sharp and thin damage profile within the gradient damage framework, several approaches, for example, the conceptual transient length effect model, 46,47 the localizing gradient damage approach 36,[48][49][50] and the gradient damage model based on a micromorphic extension, 51,52 can be considered for further developments.

ACKNOWLEDGMENT
The authors would like to acknowledge the financial and the technical support of ANSYS Inc., Canonsburg, USA, the financial support of the China Scholarship Council (No. 202008080328), as well as of the Deutsche Forschungsgemeinschaft (DFG) under Grant KA 1163/42.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.

APPENDIX A
According to Equations (39) and (44), the deviatoric stress excluding the back stresses at current time step t n+1 is expanded as , see Equation (39) = ( dev,tr t n+1 For an arbitrary non-zero second order tensor A (‖A‖ ≠ 0), an established tensorial algebraic relation

APPENDIX B
To obtain Δp∕ , the internal balances ℛ p = 0 and ℛ r = 0 are again taken into account. In detail, when the numerical solutions Δp and Δr fulfill the conditions ℛ p = 0 and ℛ r = 0 simultaneously, the complete variation of ℛ r with respect to Δp should vanish, reading 1 + b r Δp + a r g r ( r t n + Δr ) a r −1 ) Δr As a consequence, Equation (B1) leads to 1 + b r Δp + a r g r ( r t n + Δr ) a r −1 .
Moreover, the complete variation of ℛ p with respect to the strain is also 0, which is expanded as Hence, the term Δp is derived as where the partial derivative terms  Δp and Δr Δp in Equation (B4) are known in Equations (50) and (B2), respectively.

APPENDIX C
The derivatives shown in Equation (63), which have not been stated, are derived as