# A computational approach to thermomechanical fatigue

Andrei Constantinescu, Eric Charkaluk, Guy Lederer, Laetitia Verger

► To cite this version:

Andrei Constantinescu, Eric Charkaluk, Guy Lederer, Laetitia Verger. A computational approach to thermomechanical fatigue. International Journal of Fatigue, 2004, 26, pp.805-818. 10.1016/j.ijfatigue.2004.01.006. hal-00022603

HAL Id: hal-00022603

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

Submitted on 1 Jun 2021

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 computational approach to thermomechanical fatigue

Andrei Constantinescu$^{a,*}$, Eric Charkaluk$^{b}$, Guy Lederer$^{c}$, Laetitia Verger$^{c}$

$^{a}$ Laboratoire de Mécanique des Solides, CNRS UMR 7649, École Polytechnique, 91128 Palaiseau cedex, France

$^{b}$ Laboratoire de Mécanique des Lille, CNRS UMR 8107, boulevard Paul Langevin, 59655 Villeneuve d'Ascq cedex, France

$^{c}$ PSA Peugeot-Citroën, Direction des Plateformes Techniques et Achats, 18 Rue des Fauvelles, 92256 La Garenne Colombes cedex, France

This paper presents a computational approach for the lifetime assessment of structures under thermomechanical loading. One of the main features of the work is the search for simplicity and robustness in all steps of the modeling, in order to match the proposed method with industrial constraints. Among the peculiarities of the engineering environment are the imposed commercial computer programs and the short time delays for manufacturing. The proposed method is composed of a fluid flow, a thermal and a mechanical finite element computation, as well as a final fatigue analysis. The efficiency of the approach is proved by comparison with experimental results on prototype manifolds. The numerical results match both the localization of the damage zone and the lifetime.

**Keywords:** Cast iron; Viscoplasticity; Low cycle fatigue; Thermomechanical fatigue; Dissipated energy

# 1. Introduction

A series of industrial products like internal-combustion engines, turbines, nuclear reactors, etc. are subject to thermomechanical loading cycles. The actual design process is constrained by a series of specifications such as energy consumption, lifetime, manufacturing processes, materials, production rates, etc. All are summing up to the final cost of the product. The choice left to the design engineer is therefore very scarce.

In this context, materials and production processes are a priori defined and the steps of the design cycle with respect to fatigue lifetime predictions can therefore be presumed as [3]:

(1) definition of the representative load and the desired lifetime (see for example [27])

(2) modeling and testing in order to obtain the desired lifetime, generally only geometrical optimization and adjustment of boundary conditions and stiffnesses.

Moreover, the chosen design method should be compatible with the engineering environment (computing tools and testing techniques) and the time requirements of the industry, typically a couple of hours.

In this paper, we shall show that one can define a computational design method for assessing the lifetime of parts under thermomechanical loading cycles compatible with industrial specifications, and as such we shall only discuss the computational part of the second step of the design cycle presented before. This work is focused on cast iron parts for internal-combustion engines with application to the exhaust manifold. Even so, the underlying method is a general approach to this type of problem, such as the aluminum alloy cylinder head design (see [29]).

The originality of the present study lies essentially in adapting and interpreting classical notions, models and techniques in order to obtain a reliable fatigue prediction for complete 3D structures.

The damaging cycle of an exhaust manifold corresponds to the start-full load-stop cycle of the engine. A schematic view of its functioning is presented in Fig. 1. In the presented conditions, the manifold is subjected to huge temperature changes, implying that

<footnote>* Corresponding author. Tel.: +33-1-69-33-33-30; fax: +33-1-69-33-30-26.</footnote>

<footnote>E-mail address: andrei constantinescu@lms.polytechnique.fr (A. Constantinescu).</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic view of the position and the functioning of an exhaust manifold.</figcaption></figure>

at low temperatures (below 500 °C), the structure has elastoplastic behavior and that at high temperatures (above 500 °C), the structure has viscoplastic behavior.

The classical fatigue criteria, for example, Manson-Coffin [9,16], Smith-Watson-Topper [23] or strain range partitioning [10], have initially been defined in an isothermal framework for uniaxial specimens. Therefore, their ideas are not easily generalized in a multi-axial anisothermal framework where plasticity and viscosity occur simultaneously.

The fatigue problems in the nuclear industry have been generally reduced to pure creep problems at high temperatures (around 600 °C) and to crack propagation [18]. The direct application of this approach is not possible in our case as temperature changes in automotive engines imply the simultaneous existence of plasticity and creep and as no crack monitoring can be imagined for these parts.

A series of precise coupled constitutive-damage models have been developed and successfully applied for turbine blades in aerospace engines (see [5,12]). The success of these models is based, on the one hand, on fine control of the material manufacturing process, and, on the other hand, on a large number of degrees of freedom representing the material microstructure, leading to large computing time. Both arguments are not compatible with the constraints of the automotive industry. There are a large series of products with less control of the material quality than in the aerospace industry and very sharp constraints on the design period.

The method proposed in this work is based on a series of simple assumptions. The assumptions permit one to simplify the existing models and to eliminate a series of drawbacks in order to respect the industrial constraints.

The main assumption is that the cyclic material behavior is stabilized for most of the lifetime and that the evolution of the damage parameters can be neglected in the constitutive law and therefore in the

structural computation. The simple elastoviscoplastic constitutive law has been chosen to represent the material behavior at both low and high temperatures and identified from the stabilized cycle. The same behavior will determine the fatigue criterion as a function of the dissipated energy per cycle.

The paper presents, in the first section, the general material properties of the studied cast iron. The second and third sections describe the thermal and mechanical finite element analysis and the underlying assumptions. The last section discusses the fatigue analysis and presents the global lifetime prediction for structures.

## 2. The material

The studied material is a silicon molybdenum (SiMo) spherical graphite cast iron (SiMo cast iron, with the composition presented in Table 1). In the case of exhaust manifolds, the material is subject to temperatures ranging between -20 and 800 °C with heating rates attaining 10-20 °C/s.

The microstructure at room temperature is presented in Fig. 2. It exhibits a ferritic matrix (gray background) with the perlitic zones (white spots) around grain boundaries and a distribution of graphite noduli (black spots). At high temperatures, growth of the perlitic zones lead to a general decrease of the mechanical and fatigue properties of the material.

In isothermal tension-compression experiments, the material presents a general decrease of the standard material properties. The precise evolution of Young's modulus *E*, the yield limit *R*<sub>p0.2</sub> at 0.2% and the ultimate strength UTS are plotted in Fig. 3. Experiments also present an increase of the influence of viscosity, understood as strain rate dependence, with increasing temperature. Fig. 4 shows the evolution of the viscous part of the stress with temperature compared to the total stress. The viscous part of the stress corresponds to the relaxed stress in a relaxation test. As can be observed, this part represents more than 80% at high temperature for the SiMo cast iron.

One of the key points of the present work is the interpretation of the cyclic behavior of the material presented in the schematic plot in Fig. 11. The cyclic behavior of the SiMo cast iron in uniaxial isothermal experiments under constant strain range presents either a short period of hardening, if the temperature is below 500 °C, or one of softening if the temperature is above 600 °C. After this first softening or hardening period,

<table><caption>Table 1
Composition of the SiMo cast iron: elements as a percentage of weight</caption><thead><tr><th>C</th><th>Si</th><th>Mo</th><th>Mg</th><th>Ni</th></tr></thead><tbody><tr><td>3.47%</td><td>3.86%</td><td>0.63%</td><td>0.22%</td><td>0.56%</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>Fig. 2. Microstructure of the SiMo cast iron at room temperature.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Evolution of the conventional mechanical properties: yield limit $R_{p0.2}$ ultimate tensile stress limit UTS and Young's modulus E.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Evolution of the viscous part of the stress with temperature compared to the total stress. The viscous part of the stress corresponds to the relaxed stress in a relaxation test.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Evolution of the load during a thermal fatigue test corresponding to the test rig of Fig. 8. The observed saturation generalizes the stabilization obtained during isothermal cyclic test to a non-isothermal loading.</figcaption></figure>

extreme stresses describe a long slightly softening plateau, where mechanical behavior can be considered as stabilized. This stabilized period ends with a short period describing the final failure. As those observations are extended to anisothermal cycles (see Fig. 5), we can conclude that the major part of the lifetime can be described by damage independent constitutive law if the short starting and final transient periods are neglected. This justifies the choice of a stabilized limit cycle for the identification of the elastoviscoplastic constitutive law and will be discussed from the point of view of the fatigue criterion in the following sections.

### 3. The thermal and the mechanical analysis

The thermomechanical loading cycles considered for an exhaust manifold correspond to an idle-full load cycle of the engine. More precisely, this implies an inflow of hot gases through the manifold heating it up, a dwell time at maximum temperature and a cooling period [1,14]. Therefore, three computations will determine the global behavior of the structure:

* a *fluid flow* simulating the flow exhaust gases and supplying the heat exchange coefficients and the thermal fluxes;

* a *transient thermal diffusion computation* simulating heating and cooling down of the part;

* a *mechanical computation* with the temperature field as the main load parameter.

The main assumption is the uncoupling of the three computations.

First, the uncoupling between the thermal fluid computation and the thermal diffusion problem finds its justification in the two following elements.
On the one hand, the time scale of the heating of the part is of the order of $10^2$ s, whereas the engine cycles representing the explosion are of the order of $10^{-2}$ s.

On the other hand, the energetic balance of the gas during its flow through the exhaust manifold shows that the energy loss due to the heating of the part is negligible with respect to the transported energy from the inlet to the outlet. Therefore, the fluid flow is supposed to be independent of the temperature distribution of the part and an estimated average temperature will be applied as a boundary condition for the fluid flow computation.

The uncoupling of the thermal problem from the mechanical one is justified by the negligible contribution of the inelastic dissipation to heating. The thermal energy involved in the cycle can be estimated from the temperature evolution presented in Fig. 7.

Next, we shall extensively discuss only the transient thermal computation and the mechanical computation supposing that the fluid flow is given.

Let $\theta$, $u$, $\varepsilon$, $\sigma$ denote, respectively, the temperature, the displacement, the strain and the stress fields.

### 3.1. The thermal analysis

The aim of the thermal computations is the prediction of the temperature field distribution on the manifold from engine data (i.e. exhaust gas temperature and mass flow rate) which can already be estimated in the early phases of the project. The experimental temperature cycle studied here is transient and as such determined by the classical equation

$$ \rho c \frac{\partial \theta}{\partial t} - \nabla \cdot k \nabla \theta = 0 \qquad (1) $$

given an initial temperature field and a series of boundary conditions which will be discussed next. Here, $\rho$ denotes the volume mass, $c$ the specific heat of the mass and $k$ the thermal conductivity.

If the thermal conductivity $k$ is given as a function of temperature, the steady-state diffusion problem can easily be solved provided that the boundary conditions are precisely evaluated. Three families of boundary conditions can be observed: internal surface, external surface and contact conduction on contact surfaces.

The major heat flux contribution to the heating of the manifold comes from the forced convection with the exhaust gas on the inner surface of the manifold. The difficulties of the problem arise from the estimation of local convection coefficients, which demand a fair knowledge of both the fluid flow near the manifold wall and the gas temperature. The resolution of this fully turbulent anisothermal fluid problem is performed using the commercial FIRE SOFTWARE CODE [2]. Its description is out of the scope of this paper. How-

<figure><img src="image_7.png" /><figcaption>Fig. 6. Steady-state temperature field at full load. The maximal temperature, around 970 K, is obtained near the exhaust pipe, close to the turbo-charger.</figcaption></figure>

ever, let us summarize the assumptions at the basis of these computations.

* The mesh of the fluid part is generated automatically and the cell thickness in the wall neighborhood is adapted during the computations in order to get a correct evaluation of the boundary layer.

* The boundary conditions for the fluid computation are directly derived from engine and turbo-charger data, i.e. imposed gas temperature and fluid flow at the intake of the manifold and imposed pressure at the exhaust given by operating parameters of the turbo-charger. The presented temperature distributions (see Figs. 6 and 7) are obtained under full load conditions. The anisothermal $k_{\varepsilon}$ analysis of the fluid flow does not present important differences in gas temperature between intake and exhaust for the manifolds discussed here.

* The local heat transfer coefficient between the gas and the wall at full load is obtained from the Diwa-kar law [26] based on a logarithmic description of the wall model. The resulting heat transfer coefficients at full load are then projected on the internal surface of the FEM mesh. The values of these

<figure><img src="image_8.png" /><figcaption>Fig. 7. Temperature evolution at two measuring points.</figcaption></figure>
coefficients can vary from approximately 200 W m⁻² K⁻¹ to more than 1500 W m⁻² K⁻¹ near the turbo-charger and depend on the local flow rate and the temperature of the gas.

The heat exchange on the external surface has been modeled using both natural convection with air at room temperature and radiation:

$$ \phi = \phi_{\text{conv}} + \phi_{\text{rad}} \qquad (2) $$

The natural convection coefficient, h, is related to the heat exchange rate per unit surface by

$$ \phi_{\text{conv}} = h(\theta_{\text{wall}} - \theta_{\text{ext}}) \qquad (3) $$

and is derived from the classical empirical correlation [15]

$$ h = 0.53 \frac{Ra(\theta_f)^{1/4} \lambda_f}{D_{\text{ext}}} \qquad (4) $$

where λf is the air conductivity at an estimated “film” temperature and

$$ \theta_f = \frac{1}{2}(\theta_{\text{wall}} + \theta_{\text{ext}}) \qquad (5) $$

$\theta_{\text{wall}}$ and $\theta_{\text{ext}}$ denote, respectively, the wall temperature of the manifold and the surrounding air temperature; $D_{\text{ext}}$ is a characteristic external diameter and $Ra(\theta_f)$ is the Rayleigh number evaluated at $\theta_f$. Eq. (4) is valid for a laminar flow ($103 < Ra(\theta_f) < 109$) and leads in this case to values of h typically around 20 W m⁻² K⁻¹. The radiation part, $\phi_{\text{rad}}$, of the heat flow is computed by the classical radiation law

$$ \phi_{\text{rad}} = \omega(\theta)\beta(\theta_{\text{wall}}^4 - \theta_{\text{ext}}^4) \qquad (6) $$

where $\omega(\theta)$ is the emissivity of the cast iron as function of temperature and $\beta$ is Stefan's constant. The radiation part of the heat flux is quite significant in hot zones of the manifold as $\omega(\theta)$ reaches 0.9 for a wall temperature, $\theta_{\text{wall}}$, around 700 °C.

The thermal conduction between the manifold and the cylinder head or the turbo-charger is taken into account by means of thermal contact conduction, $\phi_{\text{cont}}$. The temperature field of the cylinder head in the neighborhood of the manifold, $\theta_{\text{cyl}}$, is known from earlier experimental analyses. The boundary condition, $\phi_{\text{cont}}$, on the contact surface between the cylinder head and the manifold is

$$ \phi_{\text{cont}} = \frac{1}{R_c}(\theta_{\text{wall}} - \theta_{\text{cyl}}) \qquad (7) $$

where $\theta_{\text{wall}}$ is as defined above and $R_c$ denotes the contact resistance. $R_c$ globally integrates the thermal characteristics of the gasket and of the interface and is obtained from measurements performed under representative clamping conditions.

The variation of the convection coefficient on the internal surface between low and full gas load is neglected. This can be accepted, as the revolution speed of the engine remains almost constant during the test, and as a first approximation, the convection coefficient also follows the revolution speed of the engine. Similarly, it is assumed that the conduction coefficient on the external surface as well as the contact conduction are constant during the cycle.

Supposing that all boundary conditions are known at full gas load, the transient boundary conditions are obtained under the following assumptions. At the time scale of the experiment, which is of an order of magnitude of minutes, the exchange temperature on the internal surface is linearly dependent on the evolution of load. The variation of the convection coefficient on the internal surface between small and full gas load is neglected. This can be deduced from the constant engine speed and from the linear dependence, in a first approximation, between the convection coefficient and the gas flow rate. The conduction coefficient on the external surface as well as the contact conduction are considered constant.

The computations are performed using the transient analysis in ABAQUS/standard [11]. The CPU time on an HP-V class computer is about 800 s with an FE mesh of the manifolds containing about 15 000 linear 3D brick elements, corresponding to about 66 000 degrees of freedom.

The resulting steady-state temperature field corresponding to a full load condition is plotted in Fig. 6. It should be noted that the temperature distribution exceeds 970 K near the turbo-charger and features a spatial gradient of more than 400 K between the input and the output of the manifold in less than 0.2 m.

A comparison between the computed and measured temperature evolution is presented in Fig. 7. A fair agreement between measured and computed temperature is observed, notwithstanding that the temperature rate reaches up to 15 K s⁻¹.

As a consequence, we can conclude that these results also assess ‘a posteriori’ the basic assumptions of our thermal computations.

## 3.2. The mechanical analysis

The mechanical fields are determined by the equilibrium equation

$$ \mathrm{div} \, \sigma = 0 \qquad (8) $$

and an elastoviscoplastic constitutive law under the hypothesis of small strains.

The constitutive law has to characterize the material from 20 to 700 °C. At 20 °C, it presents an elastoplastic behavior and at 700 °C, an almost purely viscous one, as described in Fig. 4. The fact that one has to
deal with a simple constitutive law over the whole temperature domain is dictated by the FE analysis and the fact that the temperature distribution is not uniform in the part. This leads to the choice of an elastoviscoplastic material model. A similar conclusion has also been drawn by Sermage et al. in [19].

The constitutive law used in this work is based on an additive decomposition of the stress into a plastic and a viscous part (for a complete description, see [13]):

$$ \sigma = \sigma^p + \sigma^v \quad (9) $$

The mechanical fields of the constitutive law, representing the two parts, i = p, v, plastic and viscous, respectively, are related through the following relation:

$$ \sigma^i = C_i(\varepsilon - \varepsilon^i - \varepsilon^{th}) \quad (10) $$

$$ \dot{\varepsilon}_i = \gamma_i \nabla f(A_i) \quad (11) $$

where $\varepsilon^{th}$ is the thermal strain tensor field given by the heating of the material. $C_p$ and $C_v$ are the fourth-order tensors of elasticity considered as isotropic and represented as such by the Young moduli $E_p$ and $E_v$ and the Poisson ratios $\nu_p$ and $\nu_v$.

$\varepsilon^p$ and $\varepsilon^v$ are two inelastic strain tensors corresponding to the plastic and viscous part of the model, where $f(A_i)$ represents the yield criterion and $A_i$ is a function of stresses and strains. The exact expression of $A_i$ will be detailed later. The preceding equations are complemented with the Kuhn–Tucker complementarity and the consistency conditions.

The plastic flow is given by a von Mises equivalent stress yield criterion, $J_2$, corresponding to the second invariant of the deviatoric part of the stress tensor, and a linear kinematic hardening rule, expressed as

$$ A_p = \operatorname{dev} \sigma_p - \frac{3}{2} H \varepsilon_p \quad (12) $$

$$ \gamma_p = \sqrt{\frac{3}{2} \langle J_2(A_p) - \sigma_Y \rangle} \quad (13) $$

where the expression of consistency parameter $\gamma_p$ corresponds to a plastic loading. $H$ and $\sigma_Y$ denote the kinematic hardening and the yield limit, respectively.

The viscous flow is represented by a three-dimensional Norton–Hoff flow rule, written as

$$ A_v = \operatorname{dev} \sigma_v \quad (14) $$

$$ \gamma_v = \sqrt{\frac{3}{2} \left( \frac{J_2(A_v)}{\eta} \right)^m} \quad (15) $$

where the expression of consistency parameter $\gamma_v$ corresponds to a viscoelastic loading. $\eta$ denotes the viscosity.

The model has finally six temperature-depending parameters, $E_p$, $E_v$, $H$, $\sigma_Y$, $\eta$, $m$.

The Poisson ratios $\nu_p$ and $\nu_v$ have been considered as constant and equal to 0.28.

The other parameters of the model have been identified from cyclic tension–relaxation tests at several temperatures and the parameters have been considered to vary linearly between these temperatures. An overview of the variation of the parameters with temperature can be gained from Fig. 3, which displays the variation of the standard yield limit $R_c$ at 0.2%. A discussion of the identification procedure and its implication on the final computational result is presented by Verger et al. in [28].

Regarding the identification procedure, the Young modulus $E_p$ and the hardening $H$ have been obtained from the slope at the origin and the work hardening of the experimental curve. The yield limit $\sigma_Y$ and the parameters $\eta$, $E_v$ and $m$ describing the viscous stress $\sigma_v$ are obtained from the relaxation curve and the sensitivity of the tension curve with respect to the strain rate. The Young modulus of the viscous part, $E_v$, is very important as it determines the characteristic time of the viscous relaxation.

One can remark that the yield limit lowers tremendously with increasing temperature and this will demand special care in the numerical integration algorithm. A comparison between three different integration algorithms for this constitutive law—forward Euler, semi-implicit and backward Euler with radial return— has already been discussed in a previous communication [8].

The computations presented next have been done using the Euler backward algorithm with a radial return, which proved to be unconditionally stable (for a complete reference on algorithms, see [20]). This solution showed to be essentially faster for temperature-dependent yield limit and viscosity.

The mechanical computations have been performed using the ABAQUS software [11] as for the thermal computations. From a practical point of view, the integration algorithm has been implemented by means of a User Material subroutine (UMAT) and some details of the radial return algorithm are presented in Appendix A.

The mechanical computations have been performed on two types of structures:

* a thermomechanical specimen submitted to heating by Joule effect, described in Fig. 8;

* Two different types of exhaust manifolds prototypes described as A and B, with two successive design versions for type B.

The thermomechanical tests (Fig. 8) were performed on axial symmetric clamped specimens heated by the Joule effect. The heating procedure produced a thermal gradient along the specimen axis, which also induces a variation of the mechanical fields in the same direction. This experimental machine has also been used to provide the thermomechanical fatigue test and the corre-
<figure><img src="image_9.png" /><figcaption>Fig. 8. The anisothermal fatigue machine.</figcaption></figure>

sponding results in the lifetime prediction as described in the next section.

Maximum temperatures have been varied between 40 and 700 °C with a heating rate of 20 °C s⁻¹. The maximum temperature is obtained in a region of approximately 10 mm in the center of the specimen. The maximum temperature gradient was 30–40 °C/mm. The parameters of the test are the clamp value due to the variable stiffness (proportional to the flexural moment of the beams: 183 000 and 227 000 N/mm, see Fig. 8) and the dwell time at 700 °C (60 and 900 s). This stiffness corresponds to a linear relation between the measured displacement (in mm) with the strain gauge and an applied load (in N) on the top of the test rig.

The thermal load of this experiment was numerically estimated by an electric and a thermal FEM computation. The results of the temperature distribution were compared with measurements obtained by infrared thermography and thermocouples. The temperature distribution at the maximal temperature of the cycle is represented in Fig. 9.

The mechanical computation was used to validate the accuracy of the constitutive behavior of this simple structure. In Fig. 10, we represent the measured and computed nominal axial stress versus temperature during some cycles.

The thermal loading cycle for the exhaust manifolds has been described in the previous section and mechanical computation has been performed for the same mesh. The cylinder head has been modeled as a thick plate and the coupling between the cylinder head and the manifold has been represented by a number of screws and a standard small sliding friction contact model of ABAQUS.

The mechanical analysis has been performed in two steps: one simulating the manufacturing of the engine

<figure><img src="image_10.png" /><figcaption>Fig. 9. Computed and measured temperature on the axes of the specimen in the anisothermal fatigue machine at maximal heating temperature.</figcaption></figure>

and one representing a functioning cycle of the engine. In the first step, the manifold is bolted on the cylinder head and in the second, it is submitted to the idle–full load thermal cycles described in the previous section.

In both cases, i.e. thermal fatigue specimens and exhaust manifolds, the response of the structure is multiaxial; more precisely, the computed $\sigma_{xx}$ and $\sigma_{yy}$ stresses exhibits a nonproportional affine loading as shown by Charkaluk et al. [6].

The thermomechanical analyses for both manifold versions were completed within approximately 6 h of CPU time on an HP-V class computer. This CPU time is very sensitive to the amount of inelastic deformation experienced by the structure as well as to the nonlinearities induced by frictional contact.

<figure><img src="image_11.png" /><figcaption>Fig. 10. Computed and measured nominal axial stress versus temperature for the specimen in the anisothermal fatigue machine.</figcaption></figure>
# 4. Fatigue analysis

## 4.1. Fatigue criterion

Once the mechanical fields are obtained, the major challenge is to determine a low cycle fatigue (LCF) criterion in the anisothermal multiaxial framework. We shall define the fatigue criterion as a relation between a criterion function $\Phi$ and the number of cycles to failure of the structure $N_f$.

A series of classical fatigue criteria can be expressed under a general form [24,25] as:

$$ \Phi(\varepsilon, \varepsilon^P, \sigma, \dots) N_f^\beta = c \quad (16) $$

where $N_f$ denotes the number of cycles to failure and $\beta$, $c$ are material dependent constants. $\varepsilon$, $\varepsilon^P$, $\sigma$ are the values of the fields during a complete cycle. One can remark that an underlying hypothesis in the framework is the existence of a stabilized cycle. This assumption will be verified in our case through direct computations.

The following paragraphs review some classical fatigue criteria and some engineering applications. For the Manson-Coffin [9,16], Smith-Watson-Topper [23] and Ostergren [17] criteria, the function $\Phi$ is represented in Table 2. Here $\Delta\varepsilon^P$, $\Delta\varepsilon$, $\sigma_{max}$ denote, respectively, the plastic strain range, the total strain range and the maximum tensile stress over a uniaxial cycle. $E$ is the Young modulus and $v$ the frequency of the test cycle. $\beta$, $k$, $c$ are material parameters characterizing the fatigue behavior.

Another classical criterion is the strain range partitioning (SRP) proposed by Halford and Manson [10], which gives an a priori separation of the stress-strain cycle between creep and plastic parts. The major drawbacks of the partitioning method are that the method is empirical and that the proposed separation is not unique as pointed out by Chaboche et al. [4]. Moreover, multiaxial anisothermal computations with elastoviscoplastic constitutive laws show that creep, i.e. viscous effects, and plasticity are coupled and simultaneous phenomena.

In the present work, the first step concerns the definition of $\Phi$, and to that end, it is essential to discuss the previous classical definitions. In a multiaxial case, one has to replace the strain range with a well chosen amplitude of the strain cycle computed using a “norm” such as the von Mises equivalent stress $J_2$. Similarly,

<table><caption>Table 2
Classical low cycle fatigue criteria</caption><thead><tr><th>Criterion</th><th>Φ</th></tr></thead><tbody><tr><td>Manson-Coffin</td><td>Δε<sup>P</sup></td></tr><tr><td>Smith-Watson-Topper</td><td>√Eσ<sub>max</sub>Δε</td></tr><tr><td>Ostergren</td><td>(σ<sub>max</sub>Δε)v<sup>β(k-1)</sup></td></tr></tbody></table>

the maximum tensile stress over a cycle can be replaced by the maximum hydrostatic pressure, as the mean stress has an influence on the lifetime.

If these changes seem feasible for the passage to multiaxial isothermal problems, the passage to the multiaxial anisothermal problems can hardly be justified, as over a large range of temperature, the material parameters $\beta$, $k$, $c$ are temperature dependent (see [7]).

The notion of cycle frequency, usually denoted by $v$ [17], is easily defined for an experiment on specimens. In this case, it correctly features the effect of the loading rate. During a multiaxial anisothermal cycle, the strain rate depends not only on the spatial location of the considered point, but also on the time instant when it has been computed. Therefore, $v$ is only a global mean value and cannot be directly associated with the creep phenomenon.

As a consequence, one has to introduce another intrinsic definition for the criterion function $\Phi$.

The second step concerns the definition of the number of cycles to failure $N_f$ in the case of a structure. In order to better understand the underlying phenomena, let us revisit the interpretation of the cyclic behavior of the material proposed by Skelton [21,22]. The evolution of the apparent load during a LCF test under imposed strain range generally displays three important phases (see Fig. 11). These steps are defined by the stabilization, the tangent and the final points, denoted by subscripts: sta, tan, and fin.

The stabilization phase corresponds to a cyclic hardening or softening behavior until $N_{sta}$ cycles. Typically for the studied cast iron, we have cyclic hardening up to 500 °C and cyclic softening above 600 °C. During this period, the microstructure is subject to a plastic shakedown cycle under the macroscopic imposed

<figure><img src="image_12.png" /><figcaption>Fig. 11. Evolution of the global load in isothermal uniaxial experiments under fixed strain range.</figcaption></figure>
strain. Infinitesimal cracks at the length scale of µm can initiate and propagate [22].

After $N_{\text{sta}}$ cycles, the material behavior is quite often stabilized and small cracks grow until $N_{\text{tan}}$ cycles. At $N_{\text{tan}}$ cycles, existing cracks reach a macroscopic length scale and become mechanically important for the structure. This means that the macroscopic fields are influenced by the presence of the crack and failure is to be expected after a small number of cycles, at $N_{\text{fin}}$.

Skelton [21,22] showed that the accumulated dissipated energy until stabilization $N_{\text{sta}}$ can be considered constant for a given material. These results also suggested that the accumulated dissipated energy till the hardening or the softening stabilization point is reached can be used as a crack initiation criterion in LCF. It is important to remember that the crack initiation at $N_{\text{sta}}$ correspond to a microscopic crack.

The previous considerations are being extended for the definition of the criterion function $\Phi$ and the number of cycles to failure $N_f$ in the case of structures. In the design of structures, we are interested in detecting the initiation of a macroscopic crack which may endanger the integrity of the component and which is definitely not the final failure of the structure. From the point of view of the cyclic material behavior, we can then state that the initiation of the macroscopic crack corresponds to the point at $N_{\text{tan}}$ cycles as explained before. However, as $N_{\text{tan}}$ is not easily detected, we prefer to define the failure point $N_f$ in our criterion by a 10% load drop with respect to the value of the stabilized cycles at $N_{\text{sta}}$. As a consequence, $N_f$ lies between $N_{\text{tan}}$ and $N_{\text{fin}}$.

The next step in the definition of the fatigue criteria is to specify the function $\Phi$. We note that some of the classical criteria (see Table 2) are related to the dissipated energy per cycle as they combine products of strains and stresses. Therefore, as explained in detail in [6,7], we propose

$$ \Phi(x) = \int_{\text{cycle}} \sigma(x) : \dot{\epsilon}(x)dt \quad (17) $$

where the stress and strain fields should be considered at the stabilized cycle $N_{\text{sta}}$ and $x$ is the point where the fields have been computed.

As the fields are not homogeneous on the structure, we shall rewrite the criterion as

$$ \left[ \max_x \Phi(x) \right] N_f^\beta = c \quad (18) $$

The parameters $\beta$, $c$ of the fatigue criterion have been established using classical isothermal LCF. The LCF tests have been carried out on classical specimens at different constant temperatures between 350 and 700 °C under strain control. The total strain rate was $10^{-3}$ s$^{-1}$, the strain ranges were 0, 25%, 0, 5% and 1% and the proportion between extreme strains was:

$\varepsilon_{\text{max}}/\varepsilon_{\text{min}} = R = -1$. The dissipated energy has been evaluated on the stabilized cycle by a numerical integration procedure from the stress-strain experimental curves (for details, see [7]). The analysis of those classical LCF tests show that temperature does not have an important influence on the value of the energy [6]. Therefore, $\beta$ and $c$ are identified using all the realized tests at all temperatures and correspond to the mean value in a linear regression sense. These tests have been completed with similar experiments at 700 °C which included a dwell period at maximal and minimal strain.

A series of complementary thermomechanical fatigue tests were conducted in order to assess the fatigue criterion for the machine described in Fig. 8 and in the previous section. For these tests, the dissipated energy has been evaluated on the stabilized strain-stress cycles computed by the FEM.

The experimental versus predicted lifetime for the tests is presented in a log-log plot in Fig. 12 and shows the efficiency of the chosen criterion.

## 4.2. Fatigue predictions for the exhaust manifolds

For the exhaust manifolds, the mechanical fields have been obtained by the FE simulation described in the previous section and the fatigue criterion function $\Phi$ has been computed as post processing. Two prototypes, A and B, of exhaust manifolds have been tested.

The first remark concerns the number of cycles needed to achieve a stabilized limit cycle for the mechanical field, leading to a stabilization of the cyclic dissipated energy. It appeared that three cycles were enough to get a representative value of the dissipated energy. Therefore, the lifetime results are computed

<figure><img src="image_13.png" /><figcaption>Fig. 12. Estimated number of cycles to failure versus experimental number of cycles to failure on a log-log scale. The two dashed lines represent one and three standard deviation, respectively, of twice the lifetime.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 13. Experimental crack observation on the type A manifold prototype.</figcaption></figure>

using the third loading cycle. It is important to note that this is a pure observation and that to our knowledge, no mathematical results are available to predict the existence of the limit cycle and the number of cycles needed to attain it in such a nonisothermal context.

The second remark concerns the experimental number of cycles to failure $N_f$, which corresponds to the initiation of a macroscopic crack on the exhaust manifold. As the inspection of the manifold for cracks is only done at the exterior surface and not at every cycle, the number should be understood only as an upper bound with a certain error.

For both types of manifolds, we can see from Fig. 12 that the lifetime predictions are in good agreement with the experimental observations.

For the type A manifold, we have represented in Fig. 13 the position of the detected crack and in Fig. 14 the distribution of the criterion function $\Phi$.

<figure><img src="image_15.png" /><figcaption>Fig. 14. Distribution of the computed fatigue criterion on the type A manifold prototype.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15. Experimental crack observation on the type B manifold prototype (version 1).</figcaption></figure>

Similar results are displayed for type B manifold in Figs. 15 and 16. Figs. 17 and 18 show two damaged zones of the type B (version 1) manifold. The criterion

<figure><img src="image_17.png" /><figcaption>Fig. 16. Distribution of the computed fatigue criterion on the type B manifold prototype (version 1).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. Experimental crack observation on the type A manifold prototype (version 1).</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. 18. Distribution of the computed fatigue criterion on the type A manifold prototype (version 1).</figcaption></figure>

function in Fig. 16 corresponds to the inner surface of the manifold where it presents a maximal value. The fact that the maximum is on the interior surface suggests that the crack will initiate at the interior and will propagate toward the exterior, which was confirmed by the direct material observations.

Both zones observed on the type B manifold prototype present approximately the same level of cracking danger, and the lifetime estimation falls within the standard deviation of the criterion as shown in Fig. 12.

Fig. 19 presents an example of von Mises equivalent stress distribution at full load on the type B manifold. It appears clearly that the highest von Mises equivalent stresses are located in the coldest parts of the manifold (see Fig. 2). This is coherent with the evolution of yield limit with temperature sketched in Fig. 3 but it gives no indication about the location of the critical zone. It can therefore not be used as criterion function Φ.

Computations performed on the type B (version 2) manifold exhibited an important increase in lifetime: the rear zone is not critical anymore, and the criterion shows a significant lifetime increase for the first zone, as detailed in Fig. 20. This predicted lifetime increase

<figure><img src="image_20.png" /><figcaption>Fig. 19. Distribution of the computed von Mises equivalent stresses on the type B manifold prototype (left: type B version 1, right: type B version 2).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 20. Distribution of the fatigue criterion with the critical zone on two consecutive design versions of the type B manifold prototype.</figcaption></figure>

of 150% between versions 1 and 2 can be considered satisfactory compared with the experimental lifetime increase of 110% from the point of view of industrial design.

These results show that the proposed criterion function Φ provides a good spatial localization of the crack initiation zones and that the criterion predicts the lifetime of prototypes within reasonable limits. Moreover, we have observed that the fine changes in geometry between the two versions of the prototype are correctly interpreted in terms of lifetime increase by the criterion.

## 5. Conclusion

A complete lifetime estimation approach has been derived for structures undergoing thermomechanical fatigue and applied to the design of cast iron exhaust manifolds. The structural approach integrates a simple elastoviscoplastic description of the high temperature behavior of the material. This material law allows sufficient precision without major CPU time increase, implying that this approach can be included in a classical development schedule of such parts. The results of these computations, associated with a reasonable simple anisothermal LCF criterion, show good agreement between predicted and observed critical zones for two versions of a given manifold. Numerical and experimental lifetime increase between the two versions is also satisfactorily correlated. This approach is not limited to exhaust manifolds or cast iron: further applications could concern steel or aluminum alloys and other parts in the automotive industry or elsewhere.

Some aspects of the actual approach could be elaborated, enhancing the results, but should keep the simplicity and robustness of the present version.

For example, a better fatigue criterion should take into account the influence of the mean stress. The mean stress was automatically eliminated in the present
approach by the deviatoric character of the constitutive law and implicit character of the dissipated energy.

Another example is the constitutive law, where non-linear hardening could be included in the formulation. In this respect, one should be aware that this will probably improve the representation of the limit cycle on simple experiments, but will influence in an a priori unpredictable way the cycles of the 3D structure.

These two aspects will be considered in future research on the subject.

## Acknowledgements

The authors would like to thank Mr K. Dang Van, LMS-CNRS-École Polytechnique, and Mr A. Bignonnet, PSA Peugeot-Citroën, for fruitful discussions and continuous support during the fulfillment of this project.

## Appendix A. The radial return algorithm for the constitutive law

This section presents a series of results regarding the radial return algorithm and the consistent tangent moduli for the plastic and viscous branches of the chosen constitutive law. The computations follow the general presentation given in Simo and Hughes [20].

### Notations

$A$ $A = s - H \cdot \varepsilon_{vp}$

$A_{n+1}^*$ trial stress

$C$ elastic moduli tensor (4th rank)

$J$ deviatoric identity tensor (4th rank)

$\Delta\varepsilon$ total strain increment

$\Delta\varepsilon^v$ viscous strain increment

$\Delta\varepsilon^p$ plastic strain increment

$\gamma$ inelastic strain rate norm

$H$ hardening modulus

$\eta$ viscosity norm

$J_2$ second invariant of the deviatoric stress tensor

$\lambda, \mu$ Lame's coefficients

$m$ viscosity exponent

$n$ actual increment

$n+1$ estimated increment

$s$ deviatoric stress tensor

$\sigma_y$ yield limit

$\Delta t$ time increment

### A.1. Radial return in elastoplasticity with linear hardening

The trial stress $A_{n+1}^*$ is defined using the stress $A$:

$$A_n = \operatorname{dev}(C_n: \varepsilon_n) - (2\mu_n + H_n)\varepsilon_n^p$$

$$A_{n+1} = \operatorname{dev}(C_{n+1}: \varepsilon_{n+1}) - (2\mu_{n+1} + H_{n+1})\varepsilon_{n+1}^p$$

$$A_{n+1} = \operatorname{dev}(C_{n+1}: \varepsilon_n) + \operatorname{dev}(C_{n+1}: \Delta\varepsilon) - (2\mu_{n+1} + H_{n+1})\varepsilon_n^p - (2\mu_{n+1} + H_{n+1})\Delta\varepsilon^p$$

$$A_{n+1} = A_n + \operatorname{dev}(\Delta C: \varepsilon_n) - 2(\mu_{n+1} - \mu_n)\varepsilon_n^p + \operatorname{dev}(C_{n+1}: \Delta\varepsilon) - (2\mu_{n+1} + H_{n+1})\Delta\varepsilon^p - \Delta H \varepsilon_n^p$$

The trial stress is defined, supposing that on plastic flow occurs during the increment, i.e. $\Delta\varepsilon^p = 0$:

$$A_{n+1}^* = A_n + \operatorname{dev}(\Delta C: \varepsilon_n) - 2(\mu_{n+1} - \mu_n)\varepsilon_n^p + \operatorname{dev}(C_{n+1}: \Delta\varepsilon) - \Delta H \varepsilon_n^p$$

and

$$A_{n+1} = A_{n+1}^* - (2\mu_{n+1} + H_{n+1})\Delta t \dot{\varepsilon}_{n+1}^p$$

The normality rule at instant $n+1$ is given by

$$\dot{\varepsilon}_{n+1}^p = \gamma_{n+1} \frac{A_{n+1}}{\|A_{n+1}\|} = \gamma_{n+1} \frac{A_{n+1}^*}{\|A_{n+1}^*\|} = \frac{\Delta\varepsilon^p}{\Delta t} \quad (A.1)$$

This leads to:

$$A_{n+1} = A_{n+1}^* - (2\mu_{n+1} + H_{n+1})\Delta t \gamma_{n+1} \frac{A_{n+1}^*}{\|A_{n+1}^*\|}$$

The norm of the preceding expression is then written as

$$\|A_{n+1}\| = \|A_{n+1}^*\| - (2\mu_{n+1} + H_{n+1})\Delta t \gamma_{n+1} \quad (A.2)$$

The consistency condition $f_{n+1} = 0$ can be written as

$$J_2(A_{n+1}) = \sigma_y$$

which implies that

$$\sqrt{\frac{3}{2}} (\|A_{n+1}^*\| - (2\mu_{n+1} + H_{n+1})\Delta t \gamma_{n+1}) = \sigma_y$$

and translates into

$$\gamma_{n+1} = \frac{\sqrt{\frac{3}{2}} \|A_{n+1}^*\| - \sigma_y}{\sqrt{\frac{3}{2}}(2\mu_{n+1} + H_{n+1})\Delta t} \quad (A.3)$$

Finally, the form of the plastic strain increment is given by

$$\Delta\varepsilon^p = \gamma_{n+1} \Delta t \frac{A_{n+1}^*}{\|A_{n+1}^*\|} \quad (A.4)$$

### A.2. Radial return in nonlinear viscoelasticity

As in the elastoplastic case, we define a trial stress $A_{n+1}^*$, supposing an elastic loading, i.e. without a viscous increment: $\Delta\varepsilon^v = 0$:

$$A_{n+1}^* = A_n + \operatorname{dev}(\Delta C: \varepsilon_n) - 2\mu_{n+1}\varepsilon_n^v + \operatorname{dev}(C_{n+1}: \Delta\varepsilon)$$
or

$$A_{n+1} = A_{n+1}^{*} - 2\mu_{n+1}\Delta t \varepsilon_{n+1}^{v}$$

The normality rule is expressed at the n+1th time instant as

$$\dot{\varepsilon}_{n+1}^{v} = \gamma_{n+1} \frac{A_{n+1}}{\|A_{n+1}\|} = \gamma_{n+1} \frac{A_{n+1}^{*}}{\|A_{n+1}^{*}\|} = \frac{\Delta \varepsilon^{v}}{\Delta t} \quad (A.5)$$

Using Eq. (A.5), we obtain

$$A_{n+1} = A_{n+1}^{*} - 2\mu_{n+1}\Delta t \gamma_{n+1} \frac{A_{n+1}^{*}}{\|A_{n+1}^{*}\|}$$

The computation of the norm leads to

$$\|A_{n+1}\| = \|A_{n+1}^{*}\| - 2\mu_{n+1}\Delta t \gamma_{n+1} \quad (A.6)$$

and the flow rule can be written under an incremental form:

$$\frac{\Delta \varepsilon^{v}}{\Delta t} = \frac{3}{2} \left( \frac{J_2(A_{n+1})}{\eta} \right)^m \frac{A_{n+1}}{J_2(A_{n+1})} = \gamma_{n+1} \frac{A_{n+1}^{*}}{\|A_{n+1}^{*}\|}$$

where

$$J_2(A_{n+1}) = \sqrt{\frac{3}{2}} \|A_{n+1}\|$$

The last expression can also be expressed in the following form:

$$J_2(A_{n+1}) = J_2(A_{n+1}^*, \gamma_{n+1}) = \sqrt{\frac{3}{2}} (\|A_{n+1}^*\| - 2\mu_{n+1} \Delta t \gamma_{n+1})$$

The consistency parameter $\gamma_{n+1}$ is then computed as the solution of the following implicit equation:

$$\gamma_{n+1} = \sqrt{\frac{3}{2}} \left( \frac{J_2(A_{n+1}^*, \gamma_{n+1})}{\eta} \right)^m \quad (A.7)$$

As a practical method of solution, the classical Newton-Raphson algorithm has been used.

### A.3. Computation of the tangent modulus

In order to compute the tangent modulus, we start from the following relation:

$$\sigma = C: \varepsilon^e$$

expressed in time discretized form as

$$\Delta \sigma = C_{n+1}: \Delta \varepsilon^e + \Delta C: \varepsilon_n^e$$

The determination of the tangent modulus is reduced to the computation of

$$\frac{\partial \Delta \sigma}{\partial \Delta \varepsilon}$$

Using the previous notations and formulas, we obtain, after lengthy computations a tangent modulus for the

plastic and, respectively, viscous branch of the following forms:

$$\frac{\partial \Delta \sigma^p}{\partial \Delta \varepsilon} = C - 2\mu \Delta t \times \left[ \frac{C:J}{(2\mu + H)\Delta t} \frac{A_{n+1}^* \otimes A_{n+1}^*}{\|A_{n+1}^*\|^2} + \gamma_{n+1} \times \frac{1}{\|A_{n+1}^*\|} \frac{C:J - C:J}{C:J} \frac{A_{n+1}^* \otimes A_{n+1}^*}{\|A_{n+1}^*\|^3} \right]$$

$$\frac{\partial \Delta \sigma^v}{\partial \Delta \varepsilon} = C - 2\mu \Delta t \times \left[ T_1 C:J \frac{A_{n+1}^* \otimes A_{n+1}^*}{\|A_{n+1}^*\|^2} + \gamma_{n+1} \times \left( \frac{1}{\|A_{n+1}^*\|} \frac{C:J - C:J}{C:J} \frac{A_{n+1}^* \otimes A_{n+1}^*}{\|A_{n+1}^*\|^3} \right) \right]$$

## References

[1] Anderson DH, Bisaro DR, Hann DM, Olree M. A thermoviscoplastic FE model for the strain prediction in high temperature, thermal cycling applications for silicon molybdenum nodular. SAE paper N 980697, 1998.

[2] AVL List Gmbh. FIRE menu reference, Graz, 1996.

[3] Bignonnet A. Fatigue design in the automotive industry. In: Dang Van K, Papadopoulos IV, editors. High-cycle fatigue. New York: Springer, Wien; 1999 (CISM lecture no. 392).

[4] Chaboche JL, Policella H, Kaczmarek H. Applicability of the SRP method and creep-fatigue approach to the LCFHT life prediction of IN-100 alloy. Technical Report, Réunion AGARD, 1978.

[5] Chaboche JL, Gallerneau F. An overview of the damage approach of durability modelling at elevated temperature. Fatigue Fract Eng Mater Struct 2002;24(6):405-18.

[6] Charkaluk E, Bignonnet A, Constantinescu A, Dang Van K. Fatigue design of structures under thermomechanical loading. Fatigue Fract Eng Mater Struct 2002;25(12):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] Charkaluk E, Verger L, Constantinescu A, Lederer G, Stolz C. Lois de comportement viscoplastiques anisothermes pour calculs cycliques sur structures. In: Guédra-Degeorges D, Ladevèze P, Raous M, editors. 4eme Colloque National en Calcul des Structures. Teknea; 1999. p. 575-80.

[9] Coffin LF. A study of the effects of cyclic thermal stresses on a ductile metal. ASME paper 1953;53-A76:931-50.

[10] Halford GR, Manson S. Life prediction of thermal-mechanical fatigue using strain range partitioning. Therm Fatigue Mater Components, ASTM STP 1976;612:239-54.

[11] Hibbitt, Karlsson & Sorensen Inc. ABAQUS/Standard user manual, 1998.

[12] Kanoute P, Pacou D, Poirier S, Gallernau F, Cardonna JM. Thin wall thermal gradient: experimental study and fatigue life analysis on multiple perforated components. In: Rémy, Petit, editors. Temperature-fatigue interaction. Elsevier; 2001. p. 341-50.
[13] Kichenin J. Comportement mécanique du polyéthylène. PhD Thesis, École Polytechnique, Palaiseau, France, 1992.

[14] Lederer G, Charkaluk E, Verger L, Constantinescu A. Numerical life assessment of engine parts submitted to thermo-mechanical fatigue, application to exhaust manifolds. SAE Technical paper series 2000-01-0789, 2001.

[15] McAdams WH. Heat transmission, 3rd ed. Mc Graw Hill; 1954.

[16] Manson S. Behavior of materials under conditions of thermal stress. Technical Report NACA TN 2933, 1953.

[17] Ostergren WJ. A damage function and associated failure equations for predicting hold time and frequency effects in elevated temperature, low cycle fatigue. J Test Eval 1976;4:327-39.

[18] RCC—MR. Technical appendix A3-1z, Technical Report, AFCEN, Paris, 1985.

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

[20] Simo JP, Hughes TJR. Computational inelasticity. Springer Verlag; 1998.

[21] Skelton RP. Energy criteria for high temperature low cycle fatigue. Mater Sci Tech 1991;7(1991):427-39.

[22] Skelton RP, Vilhemsen T, Webster GA. Energy criteria and cumulative damage during fatigue crack growth. Int J Fatigue 1998;20(9):641-9.

[23] Smith KN, Watson P, Topper TH. A stress-strain function for the fatigue of metals. J Mater 1970;5(4):767-78.

[24] Socie DF, Marquis GB. Multiaxial fatigue. Society of Automotive Engineers; 2000.

[25] Suresh S. Fatigue of materials. Cambridge University Press; 1998.

[26] Taine J, Petit JP. Transferts thermiques et mécanique des fluides anisothermes. Paris: Dunod; 1989.

[27] Thomas JJ, Bignonnet A. L'équivalence fatigue—etat de l'art. 18ème journées de printemps, Société Française de Metallurgie et des Matériaux, Paris, 1999.

[28] Verger L, Constantinescu A, Charkaluk E. On the simulation of large viscoplastic structures under anisothermal cyclic loading. In: Murakami, Ohno, editors. IUTAM creep in structures. Kluwer; 2000.

[29] Verger L, Constantinescu A, Charkaluk E. Thermomechanical fatigue design of aluminium components. In: Rémy, Petit, editors. Fatigue-temperature interaction. Elsevier, 2001. p. 309-18.