# A comparison of lifetime prediction methods for a thermal fatigue experiment

Sébastien Amiable, Stéphane Chapuliot, Andreï Constantinescu, Antoine

Fissolo

► To cite this version:

Sébastien Amiable, Stéphane Chapuliot, Andreï Constantinescu, Antoine Fissolo. A comparison of lifetime prediction methods for a thermal fatigue experiment. International Journal of Fatigue, 2006, 28, pp.692-706. 10.1016/j.ijfatigue.2005.09.002. hal-00111459

HAL Id: hal-00111459

https://hal.science/hal-00111459v1

Submitted on 23 Jul 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A comparison of lifetime prediction methods for a thermal fatigue experiment

Sébastien Amiable<sup>a</sup>, Stéphane Chapuliot<sup>b</sup>, Andreï Constantinescu<sup>a,*</sup>, Antoine Fissolo<sup>b</sup>

**Abstract:** This paper is dedicated to the comparison of several numerical models for estimating the lifetime in a fatigue experiment. The models simulate the SPLASH experiment, which produces thermal fatigue by locally quenching stainless steel specimens. All models predict first a stabilized mechanical state (plastic shakedown) and then a lifetime prediction using several fatigue crack initiation criteria.

The numerical methods are either completely nonlinear or combine approximate elastic solutions obtained from minimizing a potential energy or closed form solutions with a Neuber or Zarka technique to estimate directly the elastoplastic state.

The fatigue criteria used are Manson, dissipated energy and dissipated energy combined with a hydrostatic pressure term. The latter had provided a best prediction over a series of anisothermal and isothermal LCF experiments in a classical fatigue analysis.

The analysis shows that for fatigue criteria taking into account the triaxiality of the mechanical response we obtain a systematic and conservative error. As a consequence of this work, we show that simplified models can be used for lifetime prediction. Moreover the paper provides a general technique to assess from the point of view of the design engineer the combination between a numerical method and a fatigue criterion.

**Keywords:** Thermo-mechanical fatigue; Quenching; Finite elements method; Numerical methods

## 1. Introduction

The lifetime of a large class of industrial parts is estimated by numerical computations during their design cycles in order to assess their response under normal and exceptional loading. The existing design methods can be composed from closed-form, finite element solutions, boundary element solutions,... and fatigue criteria. The choice is rather large after decades of developments in solid mechanics as well as in fatigue experiments and criteria. It is important to notice that the choice of the design method is not only based on the global accuracy of the technique, but also by engineering environment: time constraints, available numerical, lifetime assessment codes...

The aim of this paper is to provide a comparison of a series of numerical techniques ranging from full nonlinear

elastoplastic finite element solutions lasting one day on existing computers, to closed-form or approximate solution obtained almost instantaneously. The comparison will assess the precision of the lifetime predictions through different methods and prove that simple models can still be used provided the error is systematic, conservative and a priori known.

At the basis of the comparison we shall consider the SPLASH experiment [13], in which a specimen heated by Joule effect is cyclically cooled down by a water spray on its surface. The cooling, similar to quenching, creates large temperature gradients over time and space, and a confined plastic zone. After cycling, an evolving crack network can be observed under certain loading conditions and observations have been reported in [26,27].

Initially results have been interpreted according to the classical industrial lifting codes ASME, R5, RCC-M, by a simple closed-form solution combined with a Manson type criterion, with an equivalent plastic strain amplitude as the damage parameter. However it leads to unsatisfactory lifetime predictions. Recently complete 3D finite element computations [4] provided the complete mechanical fields during the cycle. They showed an important evolution of the stress triaxiality

<footnote><sup>a</sup> Solid Mechanics Laboratory, CNRS UMR 7649, Department of Mechanics, Ecole Polytechnique, 91128 Palaiseau cedex, France</footnote>

<footnote><sup>b</sup> DM2S/SEMT, CEA Saclay, 91191 Gif sur Yvette cedex, France</footnote>

<footnote>* Corresponding author. Tel.: +33 1 69 33 33 30; fax: +33 1 69 33 30 26.
E-mail address: andrei constantinescu@lms.polytechnique.fr
(A. Constantinescu).</footnote>
during the experiment. Therefore a fatigue criteria combining
the hydrostatic stress and the dissipated energy per cycle has
successfully been proposed.

Fatigue analysis for structures with thermal quenching
have already been reported in the literature. Srivastava et al.
[41] analyzes the lifetime of a die-casting tool using a Manson-
Coffin fatigue criteria. A similar study is provided by
Mackin et al. for the fatigue of automotive brakes in [25].
Kerezsi et al. [21,22] reviewed a series of quenching
experiments and analyzed the results, i.e. the evolution of the
crack network in terms of stress intensity factors of the
involved cracks.

However, to our knowledge, no complete comparison of
different lifing methods involving elastoplastic finite elements,
closed-form solutions, estimation of plastic strains and
different fatigue criteria has been reported for such problems.
The aim of the paper is to propose such a comparison for a
thermal shock fatigue experiment.

The analysis of the SPLASH experiment proposed in the
sequel, more precisely the construction of the numerical
techniques, is based on a series of assumptions: uncoupling of
the heating and the mechanical power and uncoupling of the
constitutive law from damage evolution. Moreover, we shall
suppose that the cyclic loading leads to a plastic shakedown
state on the structure. As a consequence the lifetime prediction,
understood as a crack initiation prediction, is done in three
steps:

* computation of the evolution of the temperature field;

* computation of the stabilized cycle, corresponding to the plastic shakedown state (with the preceding temperature field as the loading parameter);

* lifetime prediction through a *fatigue crack initiation* criteria, based on the characteristics of the stabilized cycle.

For the computation of the stabilized cycle, we shall
compare the mechanical states obtained with five different
methods. The first two involve elastoplastic models:

* M0: a finite element solution with an elastoplastic constitutive law with non-linear kinematic hardening.

* M1: a finite element solution with an elastoplastic constitutive law with linear kinematic hardening.

The next three methods, combine an elastic computation
with an estimation of the plastic strain obtained through
different techniques:

* M2: a finite element solution with elastic constitutive law and Zarka's technique [43] for the computation of the bounds of the stabilized cycle

* M3: an elastic closed-form solution for a half-space submitted to a cylindrical temperature gap at its center and a plastic strain estimated from a Neuber [24] type formula (K<sub>v</sub> technique [32])

* M4: an approximate solution obtained by minimizing the complementary potential energy on a polynomial basis and

plastic strain estimated from a Neuber [24] type formula
(K<sub>v</sub> technique [32]).

For the lifetime prediction, three criteria have been
employed: Manson–Coffin, dissipated energy per cycle and
dissipated energy corrected with a hydrostatic pressure term.
The comparison is done for the stress and strain values of the
stabilized cycle and for the complete lifetime prediction.

The paper is organized as follows: Section 2 describes the material and the test rig of SPLASH experiment, the third section presents general assumptions of the analysis and the fourth section is then dedicated to the description of the different methods. The comparison of the stabilized mechanical states and the lifetime predictions are done in the last two sections.

## 2. Material and experimental procedure

The material studied here is a 304 stainless steel with a classical chemical composition given in Table 1 and thermo-mechanical properties in the studied temperature range (20–300 °C) given in Table 2.

The SPLASH test is a thermo-mechanical fatigue experiment designed by the CEA (similar to the test rig proposed by Marsh [29]), for the understanding of the apparition of crack networks under periodic thermal shocks. Another aim is to reproduce crack networks similar to those observed on real structures [14, 15]. Experimental observations and analysis of these crack networks have already been studied and some important results are reported in [27].

The test facility is composed of a specimen (240×30×20 mm³) continuously heated by an electrical DC current and periodically cooled by water sprayed on two opposites faces of the specimen as displayed in Fig. 1.

The cycles start with a homogeneous temperature T<sub>max</sub> in the hole specimen, followed by a local short cooling period (0.25 s) to T<sub>min</sub> and the heating up period (7.5 s) to T<sub>max</sub>. The temperature difference ΔT = T<sub>max</sub> - T<sub>min</sub>, measured at the center of the quenched zone is considered as the loading parameter of the experiment.

Thermal down-shocks produced by water-spraying induce large temporal and spatial gradients: the cooling-rate is about 600 °C/s and the gradient along depth of specimen is approximatively 100 °C/mm. The thermal loading zone is confined to some square millimeters on the surface and of a few millimeters in depth, therefore, the experiment conducts to a small localized plastic zone which develops into a crack network with growing number of cycles.

3. Modeling assumptions

The modeling assumptions of the SPLASH experiment can
be splitted in two classes: the uncoupling assumptions between
heating, mechanical behavior and damage evolution; and the
assumptions on the cyclic material behavior which will permit
to compute directly the stabilized cycle corresponding to the
plastic shakedown state of the structure.
<table><caption>Table 1
Chemical compositions (wt%), 304L stainless steel</caption><thead><tr><th></th><th>C</th><th>Mn</th><th>Si</th><th>Cr</th><th>Ni</th><th>Mo</th><th>S</th><th>P</th><th>Cu</th><th>Al</th><th>B</th><th>N</th><th>Fe</th><th>Bal.</th></tr></thead><tbody><tr><td>304L</td><td>0.031</td><td>1.48</td><td>0.55</td><td>19.4</td><td>8.6</td><td>0.23</td><td>0.003</td><td>0.028</td><td>0.17</td><td>0.025</td><td>0.0015</td><td>0.058</td><td></td><td></td></tr></tbody></table>

<table><caption>Table 2
Material properties, 304L stainless steel</caption><thead><tr><th></th><th>E<br>(GPa)</th><th>α (K<sup>-1</sup>)</th><th>R<sub>ρ0.2</sub><br>(Mpa)</th><th>ν</th><th>ρ (Kg/m<sup>3</sup>)</th><th>k (W/m/K)</th><th>CP (J/<br>K<sub>β</sub>/K)</th></tr></thead><tbody><tr><td>20 °C</td><td>196</td><td>15.9×<br>10<sup>-6</sup></td><td>175</td><td>0.3</td><td>7800</td><td>18</td><td>550</td></tr><tr><td>350 °C</td><td>172</td><td>17.8×<br>10<sup>-6</sup></td><td>105</td><td>0.3</td><td>–</td><td>–</td><td>–</td></tr></tbody></table>

### 3.1. Uncoupling assumptions

The first assumption is the uncoupling of the thermal analysis from the mechanical analysis. We have considered that the thermal analysis is independent and the heat source and the flux boundary conditions have been identified in order to match the temperature measurements. Therefore, the heat source provided by the mechanical power $\sigma : \dot{\epsilon}$ is included directly in the identified terms.

The second assumption states the uncoupling of the constitutive law and the damage evolution, i.e. the uncoupling of the mechanical analysis from the fatigue analysis.

This uncoupling assumption on damage has already successfully been used in lifetime predictions. In the high-cycle fatigue domain, it is justified by the elastic shakedown of the structures [9] and is the base concept to criteria like Dang-Van [10], Papadopoulos [35] or Morel [30].

In the low-cycle thermal fatigue domain, this hypothesis has been successfully used on 3D automotive structures [7,6,42,8] and therefore justified by the a posteriori results.

Moreover, biaxial experiments on cruciform specimens reported in [38] showed that lifetime predictions obtained with damage uncoupled from mechanical behavior were as precise as predictions obtained with coupled models.

<figure><img src="image_2.png" /><figcaption>Fig. 1. SPLASH specimen.</figcaption></figure>
## 3.2. *The cyclic behavior and the stabilized cycle*

The cyclic behavior of the 304 stainless steel [31] shows a complex evolution which can be described using a series of superposed hardening laws as proposed for example in [24, 23]. However, taking into account that the structure under consideration is submitted to $\approx 10^4–10^6$ cycles, and that the elastoplastic computation of a load increment takes $\approx 3$ min this would conduct to heavy numerical computations.

In order to simplify the analysis, we shall also suppose that the material behavior is stabilized during the complete lifetime of the structure, i.e. that a plastic shakedown state is reached and that it characterizes the lifetime. As a consequence the different hardening and softening periods are therefore neglected in the present approach (see Fig. 2). This leads in an approximate error of about 5% in the computed stress level, which is of the order of magnitude of other uncertainties in the model. We would also like to recall that this hypothesis is coherent with the interpretations given by Skelton [8,39] of the cyclic behavior of the material and its fatigue damage under anisothermal conditions. The steps defined in these papers are delimited by the *stabilization*, the *tangent* and the *final* points (Fig. 2). Therefore, a reasonable choice is to admit that the mechanical response at $N_f/2$ lies within the *stabilized period* and will correspond to the material behavior identified for computational needs. It is important to remark that our choice is compatible with the fatigue analysis where characteristics used in the interpretation are equally taken after $N_f/2$ cycles.

The experiments used for material parameter identification were uniaxial strain controlled tension experiments on cylindrical specimens reported in [31]. The variations of the stress level obtained with parameters identified from different loading conditions (which will also imply different $N_f/2$) were in the expected error of about 5%.

# 4. Numerical methods for the computation of the stabilized cycle

We describe in this section five numerical methods with various complexity in order to compute the stabilized mechanical response of the SPLASH specimen.

The thermal analysis of the SPLASH test has been modeled by finite element. Numerical results matched experimental measurements from thermocouples placed in the quenched zone at surface and in depth of specimen as displayed in Fig. 3. A complete analysis of the thermal computations is provided in [3].

As one of the motivations of this work is the comparison of computational burden of numerical methods, we mention that all computations presented in the sequel have been performed on a 3GHz bi-processor with 4Go RAM running under Linux.

Finite element models have been realized using the object-oriented finite element program Cast3M [1]. The same mesh has been used with 17,000 nodes and combined linear cubical and prismatic elements for all FEM computations. The finest zone of the mesh corresponds to the quenched zone of the specimen and the elements in this part of the mesh are 100 µm long. The spatial gradients are well represented over at least a dozen of element. In order to estimate the possible error due to large spatial gradients finer meshes have been tested and we can report that refining mesh above the present limit causes an important increase of computing time without a noticeable improvement of the results.

## 4.1. M0: finite element model with a non-linear kinematic hardening

The reference analysis will consider that the material obeys to an elastoplastic constitutive law. To respect to the modeling assumptions described in the previous section, we shall choose a constitutive law without isotropic hardening. This is coherent

<figure><img src="image_3.png" /><figcaption>Fig. 2. Evolution of the global load in isothermal uniaxial experiments under fixed strain range.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Comparison of experiments with the computed temperature field: evolution of temperature with depth (left panel) and evolution of temperature with time in the center of the quenched zone (right panel).</figcaption></figure>

with the assumption of a stabilized cyclic behavior during all lifetime of the specimen. Moreover, we can report that in the temperature range of the SPLASH test 150–350 °C, we can neglect viscous effects.

The chosen elastoplastic constitutive law is the Armstrong–Fredericks’ law [5] which is based on a classical additive decomposition of the strain into an elastic and a plastic part: $\epsilon^t = \epsilon^e + \epsilon^p$, a von Mises yield criterion and a non-linear kinematic hardening rule expressed as

$$ \sigma = \mathbf{A} : \boldsymbol{\epsilon}^e = \mathbf{A} : (\boldsymbol{\epsilon}^t - \boldsymbol{\epsilon}^p) \quad (1) $$

$$ \dot{\boldsymbol{\epsilon}}^p = \lambda \frac{\partial f}{\partial \boldsymbol{\sigma}} \quad (2) $$

$$ f = \sqrt{3J_2(\boldsymbol{\sigma} - \mathbf{X})} - \sigma_Y \quad (3) $$

$$ \dot{\mathbf{X}} = c \left( \frac{2}{3} a \dot{\boldsymbol{\epsilon}}^p - \mathbf{X} \dot{\boldsymbol{p}} \right) \quad (4) $$

where a and c are material parameters and p denotes the cumulated plastic strain defined as:

$$ p = \int_{0}^{t} \sqrt{\frac{2}{3} \dot{\boldsymbol{\epsilon}}^p : \dot{\boldsymbol{\epsilon}}^p} dt \quad (5) $$

The four material parameters, Young’s modulus E, elastic limit $\sigma_Y$ and two hardening coefficients a and c, are identified on isothermal uniaxial strain-controlled experiments on cylindrical specimens reported in [31]. The parameters have been obtained by minimization of a least squares distance between computed and measured stresses. The procedure is based on a Levenberg–Marquart minimization algorithm (from Mathematica [2]) and a gradient computed using finite differences based on direct computations with Cast3M [1]. Without discussing further details of the identification procedure here, let us just observe the match between the computed and experimental material response as displayed in Fig. 4. The identified parameters (Table 3) were in good

agreement with the orders of magnitudes reported in the literature [24].

Let us mention that the identified parameters correspond to 320 °C and are in negligible distance with those at 165 °C. The relative errors for E, a, c, $\sigma_Y$ are less than % between the two temperatures. As SPLASH test temperature range lies within this interval we shall further consider that the material parameters are temperature independents.

The computations done in the case of the SPLASH experiment correspond to three different temperature drops at the center of the quenched zone. For each case up to 30 cycles have been computed in order to assure that the stabilized cycles

<figure><img src="image_5.png" /><figcaption>Fig. 4. Comparison of an isothermal (320 °C) uniaxial LCF test with computed results using the identified constitutive law with non-linear kinematic hardening.</figcaption></figure>

<table><caption>Table 3
Materials parameters identified for the non-linear kinematic hardening rule at 320 °C</caption><thead><tr><th>E (Gpa)</th><th>ν</th><th>σ<sub>Y</sub> (Mpa)</th><th>A (MPa)</th><th>C</th></tr></thead><tbody><tr><td>174</td><td>0.3</td><td>102</td><td>114</td><td>532</td></tr></tbody></table>
has correctly been computed. One can report that the mechanical behavior of the structure is stabilized after 5 cycles. A load increment for this computations takes about 190 s which leads to 20 h to compute the 5 cycles needed for the mechanical stabilization (shakedown state).

## 4.2. M1: finite element model with a linear kinematic hardening

In order to decrease the computational burden of the reference model described in the previous subsection, we perform a finite element analysis with a simpler constitutive law then presented before in Eqs (1)–(3), with a linear kinematic hardening:

$$X = H\epsilon^p \quad (6)$$

where $H$ is the linear kinematic hardening parameter.

We still use the isothermal uniaxial strain-controlled experiments reported in [31] to identify the model. In these experiments, the observed behavior of the 304 stainless steel is non-linear and can, therefore, hardly be represented by a linear hardening rule, however experiments are matched by the computation in the strain range of 0.2–0.4% which is the domain of interest for the SPLASH experiment (see Fig. 5). The identification of linear hardening is done by trial and error and material parameters identified are given in Table 4.

As for the reference model, the mechanical response of the specimen is stabilized after 5 cycles. The main advantage of linear kinematic hardening comes from the smaller computational burden, as the computation of a load increment takes ‘only’ 120 s. This leads to 13 h for the computation of 5 cycles.

## 4.3. M2: finite element model based on Zarka's method

We shall use an elastic FEM solution and the Zarka's method [43] to estimate bounds of stress and plastic strain ranges on the stabilized cycle. This technique is essentially based on shakedown theory [28,16] for standard generalized

<table><caption>Table 4
Materials parameters identified for the linear kinematic hardening rule at 320 °C</caption><thead><tr><td>E (Gpa)</td><td>ν</td><td>σ<sub>Y</sub> (Mpa)</td><td>H (GPa)</td></tr></thead><tbody><tr><td>174</td><td>0.3</td><td>167</td><td>11</td></tr></tbody></table>

materials formalism [17] and is shortly presented in the Appendix A.

Two important hypothesis should be verified to apply successfully the Zarka method for the SPLASH experiment: (i) an elastoplastic material behavior with linear kinematic hardening and (ii) a radial loading path.

The first hypothesis (i) has already been discussed in the previous subsection and we shall consider here the behavior already defined and displayed in Table 4. With respect to the radial loading path, we displayed in Fig. 6 the stress path at the center of the quenched zone in the $\sigma_{yy} - \sigma_{zz}$ space. This shows that we can reasonably consider that it is radial as we recall that $\sigma_{xx} = 0$ at this point.

## 4.4. M3: closed-form solution of semi-infinite wall coupled with the Kv's method

We consider an elastic half-space (see Fig. 7). The temperature field is homogeneous with the exception of a cold inclusion, with temperature difference $\Delta T$, in the form of a semi-infinite cylinder, with radius $R$ defined as the set: $\{(x, y, z)|x \le 0, y^2 + z^2 < R\}$

The deformation of the material is characterized by the strain tensor:

$$\epsilon = \alpha \Delta T I \quad (7)$$

where $\alpha$ is the dilatation coefficient and $I$ the identity matrix.

As at the free surface $\sigma \cdot n|_{z=0} = 0$, assuming that $\sigma_{xx} = 0$ it follows that the elastic solution is:

<figure><img src="image_6.png" /><figcaption>Fig. 5. Comparison of an isothermal (320 °C) uniaxial LCF test with computed results using the identified constitutive law with linear kinematic hardening.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Evolution from $\sigma_{zz}$ with $\sigma_{yy}$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 7. Schema of the semi-infinite wall model (Ox axis oriented towards depth of the sample).</figcaption></figure>

$$ \epsilon = \begin{pmatrix} \frac{-2\nu}{1-\nu}\alpha\Delta T & 0 & 0 \\ 0 & \alpha\Delta T & 0 \\ 0 & 0 & \alpha\Delta T \end{pmatrix} \qquad (8) $$

$$ \sigma = \begin{pmatrix} 0 & 0 & 0 \\ 0 & \frac{E\alpha\Delta T}{1-\nu} & 0 \\ 0 & 0 & \frac{E\alpha\Delta T}{1-\nu} \end{pmatrix} $$

In order to avoid any confusion, we recall that equivalent stresses and strains are defined as

$$ \sigma_{\text{eq}} = \sqrt{\frac{3}{2}} \mathbf{s} : \mathbf{s} \qquad \varepsilon_{\text{eq}} = \sqrt{\frac{2}{3}} \mathbf{e} : \mathbf{e} \qquad (9) $$

where s and e are deviatoric part of stress and, respectively, strain tensors:

$$ \text{dev } \mathbf{t} = \mathbf{t} - \frac{1}{3} \text{Tr}(\mathbf{t}) \cdot \mathbf{I} $$

The equivalent stresses and strain for this solution are then obtained as:

$$ \sigma_{\text{eq}} = \frac{E\alpha\Delta T}{1-\nu} \qquad \varepsilon_{\text{eq}} = \frac{2}{3} \left( \frac{1+\nu}{1-\nu} \right) \alpha\Delta T \qquad (10) $$

This closed-form solution, denoted as the ‘semi-infinite wall’, provides a rough but instantaneous evaluation of the stress level and of the elastic strain range in the center of the quenched zone from the SPLASH specimen. In order to compute the plastic strains from the elastic solution, we shall use an extension of Neuber’s rule [24] known as the $K_v$ method [32]. The total strain range is calculated from the elastic strain range by multiplying it by a coefficient denoted $K_v$. This coefficient depends on the elastic strain range and on a

modified biaxiality ratio. The plastic strain range is obtained by using a classical additive decomposition of total strain:

$$ \Delta\varepsilon^p = \Delta\varepsilon^t - \Delta\varepsilon^e = \Delta\varepsilon^e (K_v - 1) \qquad (11) $$

### 4.5. M4: elastic solution by minimizing the complementary energy coupled with the $K_v$'s method

The elastic approximate solution used in this section stems from a direct minimization of the elastic complementary energy potential over a particular set of statically admissible fields. The method provides two advantages: (i) one can use a temperature field closer to reality than in the previous section, (ii) the statically admissible stress fields can depend on a small number of parameters providing an instantaneous solution.

The plastic strains in the SPLASH experiment are localized in a tiny elliptical zone around the quenched part of specimen (some square-millimeters at the surface and a few millimeters in depth). The rest of specimen remains at the same thermo mechanical state during a cycle. Therefore we will only consider the solution in a cylinder of radius $h$ and height $p$, using axisymmetrical coordinates (see Fig. 8).

We admit that the temperature distribution at maximum cooling point is expressed as a second order polynomial in $r$ and $z$.

$$ T(r, z) = T_0 - \Delta T + \Phi z - \frac{\Phi^2}{4\Delta T} z^2 + \frac{\Delta T}{h^2} r^2 - \frac{\Phi}{h^2} z r^2 + \frac{\Phi^2}{4\Delta T h^2} r^2 z^2 \qquad (12) $$

where $T_0$ is the initial temperature of the body and $\Delta T$ and $\Phi$ are, respectively, the temperature difference and the maximal flux at the surface of specimen.

The height $p$ of the choose cylinder is estimated from boundary conditions and is related to $\Delta T$ and $\Phi$ by:

$$ p = \frac{2\Delta T}{\Phi} \qquad (13) $$

The temperature field in the cylinder is then completely described with only four parameters which are $h$, $\Delta T$, $\Phi$ and $T_0$ and can be obtained from measurements or numerical estimations using for example a least squares fit. It is therefore not an exact solution of the heat equation but satisfies however the boundary conditions.

We recall that, the elastic stress field solution of the problem minimizes the elastic complementary energy potential [37]

$$ u(\boldsymbol{\sigma}) = \frac{1}{2} \int_{\Omega} ((\boldsymbol{\sigma} - \boldsymbol{\sigma}^{\text{th}}) : \mathbf{\Lambda}) : (\boldsymbol{\sigma} - \boldsymbol{\sigma}^{\text{th}})) dV - \int_{\partial \Omega^d} u^d \cdot ((\boldsymbol{\sigma} - \boldsymbol{\sigma}^{\text{th}}) \cdot \mathbf{n}) dS \qquad (14) $$

over statically admissible stress fields $\boldsymbol{\sigma}$, i.e. satisfying the balance equation and the boundary conditions.

The first term in Eq. (14) denotes the strain energy with $\mathbf{\Lambda}$ the fourth order tensor of elastic compliances and the second
<figure><img src="image_9.png" /><figcaption>Fig. 8. Quarter of real structure (left panel) and axisymmetrical model (right panel). The Ox axis (left panel) and the Oz axis (right panel) are oriented towards the depth of the sample.</figcaption></figure>

term is the boundary mechanical work of the traction vector on the imposed displacement fields.

The proposed admissible stress field has the following form:

$$ \sigma_{rr}(r, z) = A + Bz + Cz^2 + Dr^2 + Ezr^2 + Fr^2z^2 \quad (15) $$

$$ \sigma_{\theta\theta} = r \frac{\partial \sigma_{rr}}{\partial r} + \sigma_{rr} \quad (16) $$

$$ \sigma_{zz} = \sigma_{r\theta} = \sigma_{z\theta} = 0 \quad (17) $$

where A, B, C, D, E, F are the parameters determined by minimization of the complementary energy potential. The form of the $\sigma_{\theta\theta}$ component is imposed by the balance equation

$$ \mathrm{div} \ \boldsymbol{\sigma} = 0 $$

The minimization of the elastic complementary energy provides the stress field in the cylinder. Elastic strains are then obtained with the Hooke's law and plastic strains are estimated using the $K_v$ method as in Section 4.4.

## 5. Comparison of the plastic shakedown state

In this section we shall provide a comparison of the mechanical fields during the shakedown cycle. For the methods M0, M1, M2 we shall compare the complete solution whether for M3, M4 we shall only discuss the elastic estimations. Finally all methods will be compared by their predictions of stress and plastic strain ranges in the center of the quenched zone.

### 5.1. Results from FEM computations (M0, M1, M2)

Let us first remark the multiaxiality of the mechanical response as displayed on the left panel of Fig. 9. We remark a negligible value of $\sigma_{xx}$ compared to $\sigma_{yy}$ or $\sigma_{zz}$, opposite to the observations for strain where one can remark a huge value of $\epsilon_{xx}$ compared to $\epsilon_{yy}$ or $\epsilon_{zz}$.

On the right panel of Fig. 9 we plot the stress-strain curves obtained on the stabilized cycle at the center of the quenched

zone with the linear and non-linear hardening laws. A series of characteristics of the stabilized cycle such as stress-strain amplitude, cumulative plastic strain per cycle, dissipated energy per cycle, have similar values with both models.

From the complete elastoplastic solution obtained by applying Zarka's method (M2) we present in the right panel of Fig. 10 the equivalent plastic strain range field. As expected, one can remark a good agreement with the shape of the same field displayed in the left panel and obtained with the reference model M0.

### 5.2. Discussion on elastic estimations (M3, M4)

On Fig. 11 we compare evolutions of elastic equivalent stresses estimated with M3 and M4 with the solution obtained by an elastic finite element model. We remark that the closed-form solution of the semi-infinite wall predicts a constant stress in all directions of the specimen whether the elastic finite element solution shows a decay of stress with depth or growing distance from the quenched zone, which is equally the solution provided by the M0, M1 methods.

The elastic energy minimization provides an excellent match for the evolution within depth of specimen. The evolution for growing distance of the quenched zone (right panel of Fig. 11) is less precise. This is probably a consequence of the assumption that the state is axial symmetric and could be corrected in a future version by adding further parameters for the admissible fields.

### 5.3. Values at center of the quenched zone

In Table 5 we compare equivalent stress and plastic strain at the center of the quenched zone obtained with the five different methods.

First, concerning the complete elastoplastic finite element computations M0 and M1, we remark that for small values of strain the difference of stress results obtained by the two models is important (close to 33%). This can be explained on the
<figure><img src="image_10.png" /><figcaption>Fig. 9. Stress-strain curves obtained with M0 (left panel) and comparison of $\sigma_{yy} - \epsilon_{yy}$ curve computed with M0 and M1 (right panel).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Plastic strain range fields obtained with M0 (left panel) and M1 (right panel).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Comparison of the evolution of elastic equivalent stress with depth (left panel) and height (right panel) of specimen.</figcaption></figure>

uniaxial curves of these models (see Fig. 4 and 5) by estimating the gap between linear and non-linear curves. However the difference closes to 3% for larger values of strains.

Second, concerning finite element models, it is important to notice that equivalent stress and plastic strain ranges computed

with M1 (linear kinematic hardening) enter exactly into the bounds predicted with M2 which proves the accuracy of Zarka's method for prediction for this type of loading.

From elastic estimations associated with the *K<sub>v</sub>* method (M3 and M4), it is important to recall that stresses estimated are
<table><caption>Table 5
Results on stress range and plastic strain range from the five numerical methods at center of the quenched zone</caption><thead><tr><th></th><th></th><th>ΔT=125 °C</th><th>ΔT=150 °C</th><th>ΔT=200 °C</th></tr></thead><tbody><tr><td rowspan="2">M0: FEM<br>non-linear</td><td>Δσ (MPa)</td><td>245</td><td>270</td><td>321</td></tr><tr><td>Δε<sup>p</sup></td><td>0.07%</td><td>0.11%</td><td>0.21%</td></tr><tr><td rowspan="2">M1: FEM<br>linear</td><td>Δσ (MPa)</td><td>335</td><td>341</td><td>356</td></tr><tr><td>Δε<sup>p</sup></td><td>0.01%</td><td>0.07%</td><td>0.20%</td></tr><tr><td rowspan="2">M2: FEM<br>Zarka</td><td>Δσ (MPa)</td><td>331–338</td><td>339–349</td><td>365–380</td></tr><tr><td>Δε<sup>p</sup></td><td>0.02% - 0.<br>0.05%</td><td>0.05% - 0.<br>10%</td><td>0.19% - 0.<br>29%</td></tr><tr><td rowspan="2">M3: Semi-<br>infinite wall<br>+ K<sub>v</sub></td><td>Δσ (MPa)</td><td>550</td><td>660</td><td>880</td></tr><tr><td>Δε<sup>p</sup></td><td>0.07%</td><td>0.11%</td><td>0.17%</td></tr><tr><td rowspan="2">M4: Elastic<br>minimiz-<br>ation + K<sub>v</sub></td><td>Δσ (MPa)</td><td>371</td><td>345</td><td>594</td></tr><tr><td>Δε<sup>p</sup></td><td>0.02%</td><td>0.04%</td><td>0.08%</td></tr></tbody></table>

elastic stresses. This conducts to an overestimation of stress range when compared to elastoplastic finite element computations. But if we compare it with an elastic finite element model (see Fig. 11), we remark that the elastic energy minimization (M4) provides an excellent estimation of elastic stress while the closed-form solution (M3) gives a largely overestimate elastic stress at the surface of specimen. However, plastic strain ranges estimated from M3 are closed from the reference model M0 while those obtained from M4 are underestimated. This shows the roughness of the $K_v$ method and as a consequence, results obtained with such a method should be interpreted carefully.

## 6. Comparison of lifetime predictions

In the previous sections we have discussed a series of models to compute the stabilized mechanical cycle of the SPLASH experiment. In this section we shall use the results described previously in order to estimate the fatigue lifetime using different criteria.

The fatigue crack initiation criterion is classically defined as a local relation, i.e. in each spatial point of the structure ($x \in \Omega$)) between the values of the mechanical fields $\epsilon, \epsilon^p, \sigma, ...$ computed for the stabilized cycle and the number of cycles to failure $N_f$ of the structure:

$$ \max_{x \in \Omega} [f(\epsilon, \epsilon^p, \sigma, ...)] N_f^\beta = c \quad (18) $$

where $\beta$ and $c$ are material parameters.

The analysis presented next is based on three main assumptions. First, we suppose that the complete finite element solution with nonlinear hardening law provides the 'best' values for the mechanical fields and will be further considered as the reference. We have been obliged to proceed in this way, because no direct experimental measurement of

the mechanical fields in the quenched zone was possible for SPLASH. Second, we identify the parameters $\beta$ and $c$ for each studied fatigue criterion using the preceding elastoplastic reference computation M0. Finally, the number of cycles to failure $N_f$ has been defined as the initiation of a surface crack of length 100 µm. This value should be understood as an order of magnitude and corresponds to the mean value of the observed crack length 50–150 µm reported in [26,27].

It is important to remark that if both the elastoplastic model M0 and the fatigue criteria are precise we would have:

$$ \beta_{M0} = \beta_{test} \quad c_{M0} = c_{test} $$

where the subscripts M0 and test define if the parameters $\beta$ and $c$ have been identified using the computed results from M0 or, respectively, the uniaxial experiments. The error between M0 and test will be estimated in the end of the paper for the best fatigue criterion.

Starting from this reference data, we shall estimate the error in lifetime estimations for each chosen fatigue criteria when using the different methods to estimate the mechanical solution. The discussed criteria are a classical Manson's, a dissipated energy per cycle and a modified energy with a maximal hydrostatic pressure damage parameter for the function $f$.

We recall that the Manson's fatigue criterion is an uniaxial law based on plastic strain range defined as

$$ \Delta\epsilon^p N_f^\beta = c \quad (19) $$

where the amplitude of the plastic strain $\Delta\epsilon^p$ is directly measured from experiments. As the SPLASH experiment is multiaxial, we have to extend the amplitude to this case. The presented extension measures the diameter of the plastic strain path and is defined by

$$ \Delta\epsilon^p = \max_{t_1} \max_{t_2} \sqrt{\frac{2}{3}(\epsilon^p(t_1) - \epsilon^p(t_2)) : (\epsilon^p(t_1) - \epsilon^p(t_2))} \quad (20) $$

where $t_1, t_2$ represent different time instants of the stabilized cycle.

The corresponding values for the different models have been displayed in Table 5.

The comparison between the observed and the predicted number of cycles to failure for the Manson criterion is given in Fig. 12. We remark that all points lie between in the zone of half and twice the lifetime. An exception are the points obtained for the smallest $\Delta T$, computed with the linear kinematic hardening using the complete FEM solution or Zarka's method, were the prediction is not conservative. In the case of linear hardening, the identification was done such that the simulated and computed strain-stress curves overlap for strains of approximately 0.2–0.4%. The difference at small strains will conduct to differences in the computed cycles at small $\Delta T$ as commented in a previous section. We can therefore conclude that errors in plastic strain are not amplified by the Manson criterion. Moreover, if using this criterion we can use
<figure><img src="image_13.png" /><figcaption>Fig. 12. Comparison of experimental vs. predicted lifetime for various numerical methods using Manson's fatigue parameter.</figcaption></figure>

The simplified Neuber calculation for plastic strains, as errors are compensated.

The criterion based on the dissipated energy per cycle as a damage parameter is defined as

$$W_p N_f^\beta = c \qquad (21)$$

where the dissipated energy per cycle is integrated in each point of the structure over the complete cycle and is approximated as the product of the amplitude of equivalent stress and plastic strain rate for simplified methods:

$$W_p = \int_{cycle} \sigma : \varepsilon^p dt \approx \Delta\sigma \cdot \Delta\varepsilon^p \qquad (22)$$

The results using the dissipated energy as a damage parameter are very scattered (see Fig. 13). Let us first observe, that the stress computed of the simplified wall solution is largely overestimated, as already commented and this conducts to the observed underestimation of the lifetime in this case. As for the Manson criterion we obtain an overestimation of the lifetime in the case of the complete FEM solution with linear kinematic hardening and a small ΔT, which is again due to the underestimation of the plastic strains in this case. If other

<figure><img src="image_14.png" /><figcaption>Fig. 13. Comparison of experimental vs. predicted lifetime for various numerical methods using the dissipated energy per cycle as a fatigue criterion.</figcaption></figure>

computed points lie more or less in the zone of half and twice the lifetime, it is important to understand that sometimes this is only due to a balancing of errors. For example the simplified elastic model combined with a Neuber rule provides a overestimated stress and a underestimated plastic strain which produce a 'correct' result.

We can also observe that the errors are not systematic, with the exception of the two simplified models combined with the Neuber rule.

The modified dissipated energy with a maximal hydrostatic pressure is defined as:

$$(W_p + \alpha P_{max}) N_f^\beta = c \qquad (23)$$

with α an additional material parameter and $P_{max}$ the maximal hydrostatic pressure attended during the stabilized cycle. Again, $W_p + \alpha P_{max}$ is computed for each point of the structure.

As shown in Fig. 14, this criterion provides the best results when we compare the SPLASH test with other multiaxial experiments from the literature [20,11,40]. The fact that all points lie within half and twice lifetime of the best fit line shows the important influence of the hydrostatic pressure as well the pertinence of the modified energy criterion.

If we compare the numerical methods of this paper using the modified energy as a fatigue criterion we obtain the results displayed in Fig. 15. If the results are globally still scattered we can remark that the error is systematic. However, we remark on one hand that the lifetime is always underestimated meaning that the estimation is always conservative; and on the other hand that each method has a deviation which is similar for all computed loads.

The relative error between $\beta_{M0}$ and $\beta_{test}$ is % and for $c_{M0}$ and $c_{test}$ of %. The plot in Fig. 14 has been obtained using $\beta_{test}$ and $c_{test}$, and in Fig. 15 using $\beta_{M0}$ and $c_{M0}$. Changing between the M0 and test values would only shift the 'diagonal', i.e. the continuous line, with a small negligible distance which can be

<figure><img src="image_15.png" /><figcaption>Fig. 14. Comparison of experimental vs. predicted lifetime for various experiments and the complete elastoplastic analysis for SPLASH using the modified energy fatigue parameter. The experiments are: uniaxial isothermal tests from CEA [31], different multiaxial and uniaxial experimentations from literature [20,11].</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 15. Comparison of experimental vs. predicted lifetime for various numerical methods using the modified energy as a fatigue parameter.</figcaption></figure>

easily estimated by least square fit through the star-shaped points in Fig. 14.

The final conclusion from a fatigue point of view, drawn from Fig. 14 and 15, is that within the tested methods and criteria, the only accurate combination is the complete nonlinear finite element computations M0 combined with a modified dissipated energy with a maximal hydrostatic pressure as a damage parameter.

## 7. Conclusion

This paper presented different methods to perform a mechanical analysis of the SPLASH thermal shock test rig. From closed-form solution coupled with simple Neuber's rule to complete FEM elastoplastic with non-linear kinematic hardening, we tested five models and analyzed the obtained mechanical results on a stabilized cycle. Afterwards, we performed a fatigue lifetime estimation with three different damage parameters as fatigue criteria for each method.

From the mechanical point of view, as expected we showed that only complete FEM computations permitted the knowledge of all mechanical fields on the stabilized cycle. Closed-form solutions and simplified models provided only the ranges of the mechanical fields.

From the fatigue point of view, we showed that the most precise estimation was obtained by the complete nonlinear finite element computations combined with a modified dissipated energy with a maximal hydrostatic pressure as a damage parameter.

However, simplified methods can still be used in combination with the modified dissipated energy criteria, if one takes into account the systematic error. It is now clear that the final choice of the designer is a balance between the demanded accuracy and the computational burden one is disposed to spend on the design cycle.

It is also important to underline, that the simplified methods presented here, can still be improved without increasing their complexity. For example one can take into account the exact biaxial state of stresses in the center of the quenched zone when computing the elastic solution, or use Zarka's method instead

of Neuber's rule to provide better estimation of the mechanical cycle.

## Acknowledgements

The authors would like to thank M. Payen (IRSN, Institut de Radioprotection et de Sureté Nucléaire) for his particular support of this project.

## Appendix A. A description of Zarka's method

This Appendix is dedicated to Zarka's method [43,19] to estimate plastic strain and the stabilized cyclic. It is based on Halphen's shakedown theory [28,16,33,34] and Halphen and Son's Standard Generalized Materials formalism [17]. The presentation done next correspond to the programming of Zarka's method we have implemented in the Cast3M finite element program with help of an iterative procedure suggested in [18,19].

Let us consider a structure submitted to a quasi-static loading, occupying in the reference configuration the domain $\Omega$ with boundary $\partial\Omega$. The loading, depending on time, is characterized by body forces $\mathbf{f}$, surface tractions $\mathbf{t}^d$ on $\partial_t\Omega$ and imposed displacements $\mathbf{u}^d$ on $\partial_u\Omega$, where $\partial_t\Omega$ and $\partial_u\Omega$ are two complementary parts of the boundary ($\Omega$, i.e. $\partial_t\Omega \cap \partial_u\Omega = \emptyset$ and $\partial_t\Omega \cup \partial_u\Omega = \partial\Omega$). The problems considered here are subject to the hypothesis of small strains, therefore:

$$ \varepsilon[\mathbf{u}] = \frac{1}{2}(\nabla + \nabla^T)\mathbf{u} \qquad (A.1) $$

### A.1. Elastic and elastoplastic problem settings

In order to describe the behavior of the structure, we shall introduce the sets of statically admissible stress fields and respectively kinematically admissible strain fields as:

$$ S(\mathbf{f}, \mathbf{t}^d) = \{\boldsymbol{\sigma}|\mathrm{div} \ \boldsymbol{\sigma} + \mathbf{f} = 0 \text{ in } \Omega \text{ and} \\ \boldsymbol{\sigma} \cdot \mathbf{n} = \mathbf{t}^d \text{ on } \partial_t \Omega \} \qquad (A.2) $$

$$ C(\mathbf{u}^d) = \{\mathbf{u}|\mathbf{u} = \mathbf{u}^d \text{ on } \partial_u \Omega\} \qquad (A.3) $$

The elastic problem can be, therefore, summarized as follows:

$$ \begin{gathered} \boldsymbol{\sigma}^{el} \in S(\mathbf{f}, \mathbf{t}^d), \quad \mathbf{u}^{el} \in C(\mathbf{u}^d), \quad \varepsilon^{el} = \mathbf{\Lambda} : \boldsymbol{\sigma}^{el} + \varepsilon_0, \\ \varepsilon^{el} = \varepsilon[\mathbf{u}^{el}] \end{gathered} \qquad (A.4) $$

In order to describe the elastoplastic behavior of the structure, we assume an additive decomposition of strains into an elastic, a plastic part and eventually an initial strain:

$$ \varepsilon = \varepsilon^e + \varepsilon^p + \varepsilon_0 \qquad (A.5) $$

The set of plastically admissible stress fields is described as:

$$ \wp(\mathbf{X}, \sigma_Y) = \{\boldsymbol{\sigma}|f(\boldsymbol{\sigma}, \mathbf{X}, \sigma_Y) \le 0\} \qquad (A.6) $$

where $f$ is a von Mises yield criterion
$$f(\boldsymbol{\sigma}, \mathbf{X}, \sigma_Y) = \frac{3}{2}(\operatorname{dev} \boldsymbol{\sigma} - \mathbf{X}) : (\operatorname{dev} \boldsymbol{\sigma} - \mathbf{X}) - \sigma_Y^2 \quad (A.7)$$

with $\sigma_Y$ the yield limit in tension. The evolution of the plastic strain is determined from the normal plastic flow rule:

$$\varepsilon^p = \lambda \frac{\partial f}{\partial \boldsymbol{\sigma}} \quad (A.8)$$

where $\lambda$ denotes the plastic multiplicator.

The elastoplastic problem is then summarized as follows:

$$\boldsymbol{\sigma} \in S(\mathbf{f}, \mathbf{t}^d), \quad \mathbf{u} \in C(\mathbf{u}^d), \quad \boldsymbol{\sigma} \in \wp(\mathbf{X}, \sigma_Y), \quad (A.9)$$

$$\boldsymbol{\varepsilon} = \mathbf{\Lambda} : \boldsymbol{\sigma} + \boldsymbol{\varepsilon}^p + \boldsymbol{\varepsilon}_0, \quad \boldsymbol{\varepsilon} = \boldsymbol{\varepsilon}[\mathbf{u}]$$

In the sequel, we shall consider the difference between the elastoplastic and the purely elastic solution of the problem, i.e. the difference between (A.9) and (A.4). This problem denoted as the inelastic problem, is described by the stress field: $\mathbf{R} = \boldsymbol{\sigma} - \boldsymbol{\sigma}^e$, a strain field: $\boldsymbol{\varepsilon}^{\text{ine}} = \boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}^e$ and a displacement field: $\mathbf{u}^{\text{ine}} = \mathbf{u} - \mathbf{u}^{\text{el}}$. We note that:

$$\boldsymbol{\varepsilon}^{\text{el}} = \mathbf{\Lambda} : \boldsymbol{\sigma}^{\text{el}} + \boldsymbol{\varepsilon}_0 \quad \boldsymbol{\varepsilon}^{\text{ine}} = \mathbf{\Lambda} : \mathbf{R} + \boldsymbol{\varepsilon}^p \quad (A.10)$$

The difference problem is therefore described by:

$$\mathbf{R} \in S(0, 0^d), \quad \mathbf{u}^{\text{ine}} \in C(0^d), \quad \boldsymbol{\varepsilon}^{\text{ine}} = \mathbf{\Lambda} : \mathbf{R} + \boldsymbol{\varepsilon}^p, \quad (A.11)$$

$$\boldsymbol{\varepsilon}^{\text{ine}} = \boldsymbol{\varepsilon}[\mathbf{u}^{\text{ine}}]$$

## A.2. Transformed parameters of the structure

In the deviatoric stresses space, the Von Mises yield criterion (Eq. A.7) is a sphere centered in $\mathbf{X} = H\boldsymbol{\varepsilon}^p$, it is, therefore, convenient to consider the transformed parameter space represented by

$$\mathbf{Y} = \mathbf{X} - \operatorname{dev} \mathbf{R} = H\boldsymbol{\varepsilon}^p - \operatorname{dev} \mathbf{R} \quad \mathbf{Y} \in \wp(\boldsymbol{\sigma}^{\text{el}}, \sigma_Y) \quad (A.12)$$

where the yield criterion becomes a sphere centered in $\operatorname{dev} \boldsymbol{\sigma}^{\text{el}}$:

$$f(\boldsymbol{\sigma}^{\text{el}}, \mathbf{Y}, \sigma_Y) = \frac{3}{2}(\operatorname{dev} \boldsymbol{\sigma}^{\text{el}} - \mathbf{Y}) : (\operatorname{dev} \boldsymbol{\sigma}^{\text{el}} - \mathbf{Y}) - \sigma_Y^2 \quad (A.13)$$

In the transformed parameters space, the difference problem described by (A.11), can be written under the following form

$$\mathbf{R} \in S(0, 0^d), \quad \mathbf{u}^{\text{ine}} \in C(0^d), \quad \boldsymbol{\varepsilon}^{\text{ine}} = \mathbf{\Lambda}^* : \mathbf{R} + \frac{1}{H}\mathbf{Y}, \quad \boldsymbol{\varepsilon}^{\text{ine}} = \boldsymbol{\varepsilon}[\mathbf{u}^{\text{ine}}] \quad (A.14)$$

where $\mathbf{\Lambda}^*$, denoted as the modified elasticity compliance tensor is a 4th order tensor defined by:

$$\mathbf{\Lambda}_{ijkl}^* = \left( \frac{1 + \nu}{E} + \frac{1}{H} \right) (\delta_{ik} \delta_{jl} + \delta_{ij} \delta_{kl}) - \left( \frac{\nu}{E} + \frac{1}{3H} \right) \delta_{ij} \delta_{kl} \quad (A.15)$$

## A.3. Asymptotic response of a structure under cyclic loading

In order to assess the elastic and plastic shakedown zone of the structure, one has to refer to the following elastic shakedown theorem proved by Mandel, Zarka and Halphen in [28].

For a Standard Generalized Material with linear kinematic hardening, in the transformed parameters space, if on the periodic loading path, the intersection of all local yield surfaces centered on dev $\boldsymbol{\sigma}^{\text{el}}$ is non-void, we have an elastic shakedown and the transformed parameters of the structure at the ultimate state are inside the intersection of all yield surfaces.

For a general presentation of shakedown results see [28,16, 17,36].

As a consequence of the preceding result, in each point $x \in \Omega$ the variation of the elastic solution defined as:

$$k^* = \min_{\mathbf{Y}} \left\{ \max_{t \in [0, T]} f(\boldsymbol{\sigma}^{\text{el}}, \mathbf{Y}, \sigma_Y) + \sigma_Y^2 \right\} \quad (A.16)$$

assures an elastic shakedown if $k^* \le \sigma_Y^2$, and a plastic shakedown if $k^* > \sigma_Y^2$.

In the particular case of a radial path, the length of the loading path in the deviatoric stress space

$$d^{\text{el}} = \sqrt{\frac{3}{2} (\operatorname{dev} \boldsymbol{\sigma}_{\text{max}}^{\text{el}} - \operatorname{dev} \boldsymbol{\sigma}_{\text{min}}^{\text{el}}) : (\operatorname{dev} \boldsymbol{\sigma}_{\text{max}}^{\text{el}} - \operatorname{dev} \boldsymbol{\sigma}_{\text{min}}^{\text{el}})} \quad (A.17)$$

assures an elastic shakedown if $d^{\text{el}} \le 2\sigma_Y$ and otherwise a plastic shakedown.

## A.4. Construction of the shakedown state

The construction of the shakedown state is performed iteratively from elastic computations and the solutions of the difference problem (37) using the following algorithm:

(1) compute $\boldsymbol{\sigma}_{\text{min}}^{\text{el}}$ and $\boldsymbol{\sigma}_{\text{max}}^{\text{el}}$ as elastic solution of (A.4)

(2) construct the transformed parameters in each point according to the nature of the estimated limit state

(3) compute the solution of the problem (A.14)

(4) estimate from the preceding solution the plastic strains.

**Elastic shakedown** - First we compute a step by step analysis of the first half cycle with the real elastoplastic behavior. We obtain the transformed parameters $\mathbf{Y}_1$ at the end of this first half cycle.

Second since the final transformed parameters $\mathbf{Y}_1$ will be such as $\mathbf{Y}_1 \in \wp(\boldsymbol{\sigma}_{\text{min}}^{\text{el}}, \sigma_Y)$ and $\mathbf{Y}_1 \in \wp(\boldsymbol{\sigma}_{\text{max}}^{\text{el}}, \sigma_Y)$, $\mathbf{Y}_1$ are locally projected orthogonally on the local intersection of the plastic yield surfaces, giving $\mathbf{Y}_1$. Then the steps 3 and 4 of the general procedure are performed.

**Plastic shakedown** - In this case, only mean the value and a bounding range of the transformed parameters can be estimated, as displayed in (Fig. A1). For the mean value and respectively the lower bound we obtain
<figure><img src="image_17.png" /><figcaption>Fig. A1. Plastic range bounding $\Delta Y_{\min}$ (—), $\Delta Y_{\min}$ (---)</figcaption></figure>

$$ \mathbf{Y}_m = \frac{1}{2} (\mathrm{dev} \ \boldsymbol{\sigma}_{\max}^{el} + \mathrm{dev} \ \boldsymbol{\sigma}_{\min}^{el}) \qquad (\mathrm{A.18}) $$

$$ \begin{aligned} \Delta \mathbf{Y}_{\min} = & (\mathrm{dev} \ \boldsymbol{\sigma}_{\max}^{el} - \mathrm{dev} \ \boldsymbol{\sigma}_{\min}^{el}) \\ & - 2\sigma_Y \frac{\mathrm{dev} \ \boldsymbol{\sigma}_{\max}^{el} - \mathrm{dev} \ \boldsymbol{\sigma}_{\min}^{el}}{\|\mathrm{dev} \ \boldsymbol{\sigma}_{\max}^{el} - \mathrm{dev} \ \boldsymbol{\sigma}_{\min}^{el}\|} \end{aligned} \qquad (\mathrm{A.19}) $$

where the norm $||\cdot||$ is defined as in Eq. (A.17). The upper bound is obtained in the following way. Suppose the two extremal yield surfaces, centered in $(\mathrm{dev} \ \boldsymbol{\sigma}_{\min}^{el})$ and, respectively, in $(\mathrm{dev} \ \boldsymbol{\sigma}_{\max}^{el})$ in the transformed parameters space, move away from each indefinitely, then the amplitude of the transformed parameters grows indefinitely and its rate $\dot{\mathbf{Y}}$ tends toward the elastic loading rate $(\mathrm{dev} \ \dot{\boldsymbol{\sigma}})$ in the direction $(\mathrm{dev} \ \boldsymbol{\sigma}_{\max}^{el} - \mathrm{dev} \ \boldsymbol{\sigma}_{\min}^{el})$. Therefore, we shall assume that:

$$ \dot{\mathbf{Y}} = (\mathrm{dev} \ \boldsymbol{\sigma}_{\max}^{el} - \mathrm{dev} \ \boldsymbol{\sigma}_{\min}^{el}) \qquad (\mathrm{A.20}) $$

The ultimate plastic flow rate $\epsilon_{\infty}^p$ by performing the steps 3 and 4 of the general procedure with $\dot{\mathbf{Y}}$.

From the plastic flow rule we deduce that the transformed parameter $\mathbf{Y}^u$ at the ultimate state is on the yield surface, and the plastic strain rate is an internal normal to it. As a consequence, the value of the ultimate transformed parameter is:

$$ \mathbf{Y}^u = \mathrm{dev} \ \boldsymbol{\sigma}_{\max}^{el} - \sigma_Y \frac{\dot{\boldsymbol{\epsilon}}_{p}^{\infty}}{\|\dot{\boldsymbol{\epsilon}}_{p}^{\infty}\|} \qquad (\mathrm{A.21}) $$

The upper bound of the transformed parameter is deduced from $\mathbf{Y}^u$ and $\mathbf{Y}_m$:

$$ \begin{aligned} \Delta \mathbf{Y}_{\max} = & 2(\mathbf{Y}^u - \mathbf{Y}_m) \\ = & (\mathrm{dev} \ \boldsymbol{\sigma}_{\max}^{el} - \mathrm{dev} \ \boldsymbol{\sigma}_{\min}^{el}) - 2\sigma_Y \frac{\dot{\boldsymbol{\epsilon}}_{p}^{\infty}}{\|\dot{\boldsymbol{\epsilon}}_{p}^{\infty}\|} \end{aligned} \qquad (\mathrm{A.22}) $$

### A.5. Iterative procedure

The preceding procedure has generally to be applied iteratively in order to take into account the fact that some regions of the structure remain elastic during the entire evolution of the structure, whether other are submitted to elastic shakedown.

The procedure is initiated with $\{\epsilon_0^p, \mathbf{X}_0, \mathbf{Y}_0\}$. We the structure is partitioned in regions where plasticity mechanisms are active and inactive. This is decided by the relative position of $\mathbf{Y}_0$ with respect to the intersection of all local yield surfaces centered in $\mathrm{dev} \ \boldsymbol{\sigma}^{el}$ at the computed time instants. If the $\mathbf{Y}_0$ the plasticity

mechanisms are inactive and active otherwise. If no there is plastic evolution, $\epsilon_l^p = \epsilon_0^p$ and the elastic compliances tensor $\Delta$ is used in order to solve the inelastic problem (A.11).

If the plasticity mechanisms are active $\mathbf{Y}_l$ are chosen equal to the orthogonal projection of $\mathbf{Y}_0$ on the intersection of yield surfaces. Moreover, the modified elastic compliances tensor $\Delta^*$ is used to solve the modified inelastic (A.14).

A global elastic computation provides the unknown fields: $\mathbf{Y}_l$ for the inactive zone and $\epsilon_l^p$ for the active zone.

Finally the position of $\mathbf{Y}_l$ with respect to the intersection of yield surfaces is checked with respect with the initial assumption that the point is in an active or inactive zone and correction are done.

Iterations are then performed up to the stabilization of the active and inactive zone of the structure.

### References

[1] Documentation Cast3M. http://www-cast3m.cea.fr.

[2] Documentation mathematica. http://www.wolfram.com.

[3] Amiable S, Chapuliot S, Constantinescu A, Fissolo A. A computational lifetime prediction of a thermal shock experiment; 2005.

[4] Amiable S, Chapuliot S, Constantinescu A, Fissolo A. Splash test numerical modelisations. In: 11th International conference on fracture, Turin, Italy; 2005.

[5] Armstrong PJ, Fredericks CO. A mathematical representation of the multiaxial baushinger effect. Technical report RD/B/N731. Janvier: Berkeley Nuclear Laboratory; 1966.

[6] Charkaluk E, Bignonnet A, Constantinescu A, Dang Van K. Fatigue design of structures under thermomechanical loadings. Fatigue Fract Eng Mater Struct 2002;25:1199-206.

[7] Charkaluk E, Constantinescu A. An energetic approach in thermomechanical fatigue for silicon molybdenum cast iron. Mater High Temp 2000; 17(3):373-80.

[8] Constantinescu A, Charkaluk E, Lederer G, Verger L. A global computational fatigue design method for structures under thermomechanical loading: application to cast iron exhaust manifolds. Int J Fatigue 2004;26:805-18.

[9] Constantinescu A, Dang Van K, Maitournam H. A unified approach for high and low cycle fatigue based on shakedown concepts. Fatigue Fract Eng Mater Struct 2003;26:561-8.

[10] Dang Van K, Papadopoulos I. High cycle metal fatigue. New York: Springer; 1999.

[11] Fatemi A, Socie D. A critical plane approach to multiaxial fatigue damage including out-of-phase loading. Fatigue Fract Eng Mater Struct 1988; 11(3):149-65.

[13] Fissolo A, Marini B, Nais G, Wident P. Thermal fatigue behaviour for a 316l type steel. J Nucl Mater 1996;233-237:156-61.

[14] Haddar N, Fissolo A. 2D simulation of the initiation and propagation of crack array under thermal fatigue. Nucl Eng Design 2005;235: 945-64.

[15] Haddar N, Fissolo A, Maillot V. Thermal fatigue crack networks: an computational study. Int J Solids Struct 2005;42:771-88.

[16] Halphen B. L'accomodation des structures elastoplastiques à écrouissage cinématique. Compte Rendu à l'Académie des Sciences de Paris 1976; 283:799-802.

[17] Halphen B, Nguyen QS. Plastic and visco-plastic materials with generalized potential. Mech Res Commun 1974;1(1):43-7.

[18] Inglebert G. Analyse simplifiée des structures élasto-visco-plastiques sous chargements cycliques. PhD Thesis, Université Pierre et Maris Curie, Paris VI; 1984.

[19] Inglebert G, Frelat J. Quick analysis of inelastic structures using a simplified method. Nucl Eng Design 1989;116:281-91.
[20] Itoh T, Sakane M, Ohnami M, Socie D. Nonproportional low cycle fatigue criterion for type 304 stainless steel. J Eng Mater Technol 1995;117: 285-92.

[21] Kereszi BB, Price JWH, Ibrahim RN. Using s-n curves to analyse cracking due to repeated thermal shock. J Mater Process Technol 2004; 145(1):118-25.

[22] Kereszi BB, Price JWH, Kotousov AG. Features of fatigue crack growth due to repeated thermal shock. J Fatigue Fract Eng Mater Struct 2002; 25(2):215-22.

[23] Krempl E, Khan F. Rate (time)-dependant deformation behavior: an overview of some properties of metal and solid polymers. Int J Plasticity 2003;19:1069-95.

[24] Lemaitre J, Chaboche JL. Mécanique des matériaux solides. Paris, France: Dunod; 1988.

[25] Mackin TJ, et al. Thermal cracking in disc brakes. Eng Fail Anal 2002;9: 63-76.

[26] Maillot V. Amorçage et propagation de réseaux de fissures de fatigue thermique dans un acier inoxydable austenitique de type X2 CrNi18-09 (AISI 304L). PhD Thesis. Ecole Centrale de Lille; 2003.

[27] Maillot V, Fissolo A, Degallaix G, Degallaix S. Thermal fatigue crack networks parameters and stability: an experimental study. Int J Solids Struct 2005;42:759-69.

[28] Mandel J, Halphen B, Zarka J. Adaptation d'une structure elastoplastique à écrouissage cinématique. Mech Res Commun 1977; 4(5):309-14.

[29] Marsh DJ. A thermal-shock fatigue study of type-304 and type-316 stainless-steels. Fatigue Eng Mater Struct 1981;4(2):179-95.

[30] Morel F. A critical plane approach for life prediction of high cycle fatigue under multiaxial variable amplitude loading. Int J Fatigue 2000;22: 101-19.

[31] Mottot M. Etude du comportement en fatigue oligocyclique à 165 et 320°C du 304I (tole t112) pour de faibles niveaux de deformation (δε<sub>r</sub>1%) TN SRMA 2001-2403. Janvier: CEA; 2001.

[32] Moulin D, Roche RL. Correction of the poisson effect in the elastic analysis of low-cycle fatigue. Int J Pressure Vessel Piping 1985;19:213-33.

[33] Nguyen QS. On shakedown theorems in hardening plasticity. J Mech Phys Solids 2004;51:101-25.

[34] Nguyen QS, Pham DC. On shakedown theorems in hardening plasticity. Compte Rendu à l'Académie des Sciences de Paris 2001; 329:307-14.

[35] Papadopoulos IV. Exploring the high-cycle fatigue behaviour of metals from the mesoscopic scale. J Mech Behav Mater 1996;6(2):93-118.

[36] Sakae C, Keer LM. Application of direct method for a nonlinear-kinematic-hardening material under rolling/sliding line contact: constant ratchetting rate. J Mech Phys Solids 1997;45(9):1577-94.

[37] Salençon J. Mécanique des milieux continus. Paris: Ecole Polytechnique; 2004.

[38] Sermage JP, Lemaitre J, Desmorat R. Multiaxial creep-fatigue under anisothermal conditions. Fatigue Fract Eng Mater Struct 2000;23(3): 241-53.

[39] Skelton RP. Energy criterion for high temperature low cycle fatigue failure. Mater Sci Technol 1991;7:427-39.

[40] Socie D. Multiaxial fatigue damage models. J Eng Mater Technol 1987; 109:293-8.

[41] Srivastava A, Joshi V, Shivpuri R. Computer modeling and prediction of thermal fatigue cracking in die-casting tooling. Wear 2004;256:38-43.

[42] Thomas JJ, Verger L, Bignonnet A, Charkaluk E. Thermomechanical design in the automotive industry. Fatigue Fract Eng Mater Struct 2004; 27(10):887-95.

[43] Zarka J, Frelat J, Inglebert G, Kasmai-Navidi P. A new approach in inelastic analysis of structures. Paris: Ecole Polytechnique; 1990.