Lev Truskinovsky

Anna Vainchtein

QUASICONTINUUM MODELS OF DYNAMIC PHASE TRANSITIONS

Keywords: Quasicontinuum approximations, martensitic phase transitions, kinetic relations, lattice models, nonlocal interactions, gradient models, microkinetics, configurational forces

We propose a series of quasicontinuum approximations for the simplest lattice model of dynamic phase transition in one dimension. The approximations are dispersive and include various non-classical corrections to both kinetic and potential energies. We show that the well-posed quasicontinuum theory can be chosen in such a way that the associated closed-form kinetic relation is in excellent agreement with the predictions of the discrete theory.

Introduction.

A recent trend in replacing continuum mechanical phenomenology with microscopic simulations raised an interest in developing lattice models for propagating defects (see [START_REF] Abeyaratne | A lattice-based model of the kinetics of twin boundary motion[END_REF][START_REF] Kresse | Mobility of lattice defects: discrete and continuum approaches[END_REF][START_REF] Marder | Origin of crack tip instabilities[END_REF][START_REF] Slepyan | Models and phenomena in Fracture Mechanics[END_REF][START_REF]Kinetics of martensitic phase transitions: Lattice model[END_REF] and the literature cited therein). The large nonlinear systems of ordinary differential equations used in such studies are rather complex and are typically not amenable to a detailed parametric analysis. Continuum theories, which assume that the lattice parameter is equal to zero, are much more transparent than their discrete prototypes. However, in order to handle defects they still require external phenomenological relations governing the dynamics of singularities.

Since outside singularities continuum theories usually deliver reasonably accurate reproduction of the behavior of the underlying discrete models, efforts have been made to augment continuum models by including various internal length scales bringing in terms that describe behavior at short wave lengths. The resulting models can be called quasicontinuum since unlike the truly continuum theories they explicitly account for the cutoff distance. Among the most well known examples of the quasicontinuum models are strain-gradient elasticity, incorporating into the energy function higher derivatives of displacements (e.g. [START_REF] Mindlin | Second gradient of strain and surface tension in linear elasticity[END_REF][START_REF] Toupin | Theories of elasticity with couple-stresses[END_REF]), and fully nonlocal theory of elasticity leading to formulation in terms of integral equations (e.g. [START_REF] Kunin | Elastic Media with Microstructure I: One-Dimensional Models[END_REF][START_REF] Rogula | Introduction to nonlocal theory of material media[END_REF]). While in linear regimes the relation between discrete, continuum and quasicontinuum theories can be considered rather well understood (e.g. [START_REF] Gustaffsson | Time-dependent problems and difference methods[END_REF]), in nonlinear cases the situation is far from being clear as indicated, for example, by the status of discrete approximations for general hyperbolic equations (e.g. [START_REF] Lax | The generation and propagation of oscillations in dispersive initial value problems and their limiting behavior[END_REF]). The problem is even more difficult when the essential nonlinearity manifests itself through the loss of convexity of the energy as in the situations involving formation and propagation of defects in solids.

In the discrete problems with nonconvex energy long-wave continuum approximations develop singularities whose dynamics remains largely arbitrary. This leads to nonuniqueness and requires an additional input from the discrete theory. To overcome this problem, one can try to develop a quasicontinuum theory that combines the advantages of a homogenized description with an ability to capture all important physics at small scales. To test the effectiveness of such quasicontinuum theories it is natural to begin with the study of a motion of an isolated defect, which may be a dislocation, a crack tip or a phase boundary.

Since the dynamics of all these defects exhibit considerable similarities (e.g. [START_REF]Lattice friction for crystalline defects: from dislocations to cracks[END_REF])

we focus in this paper on a martensitic phase boundary as the simplest example of a singularity that can be adequately described already in one dimension. Similar to shocks, martensitic phase boundaries are modeled at the continuum level as propagating strain discontinuities. However, unlike shocks, they require for their unique determination an explicit specification of the rate of dissipation. This is usually accomplished by complementing the system of Rankine-Hugoniot jump conditions by a kinetic relation between the driving force acting on a phase boundary and its velocity [START_REF] Abeyaratne | A continuum model of a thermoelastic solid capable of undergoing phase transitions[END_REF][START_REF] Truskinovsky | Dynamics of nonequilibrium phase boundaries in a heat conducting elastic medium[END_REF]. Lattice models generates such relations automatically due to the phenomenon of radiative damping (e.g. [START_REF] Kevrekidis | Dynamics of lattice kinks[END_REF]). In the present paper we choose for determinacy a prototypical lattice model first introduced in [START_REF]Kinetics of martensitic phase transitions: Lattice model[END_REF] and try to reproduce the corresponding discrete kinetic relation quantitatively in the framework of various quasicontinuum approximations. The approximate kinetic relations can then provide examples of particular admissibility conditions for the long-wave continuum models.

More specifically, we develop a set of quasicontinuum approximations for a rather general class of one-dimensional fully inertial lattice models exhibiting essential nonlinearity (bi-stability) of the interatomic potential combined with long-range interactions. An important feature of our quasicontinuum theories is that the associated energy functionals contain higher-order spatial derivatives of both strain and velocity fields, introducing dispersive corrections to both kinetic and potential energies. In the context of lattice theories first models of this type were introduced by Rosenau [START_REF] Rosenau | Dynamics of nonlinear mass-spring chains near the continuum limit[END_REF][START_REF]Dynamics of dense discrete systems[END_REF][START_REF]Hamiltonian dynamics of dense chains and lattices: or how to correct the continuum[END_REF] (see [START_REF] Kevrekidis | Continuum approach to discreteness[END_REF] for a recent review). In the phenomenological theories of solids velocity gradients were previously incorporated into kinetic energy as "microkinetic" terms [START_REF] Theil | A study of a hamiltonian model for martensitic phase transformations including microkinetic energy[END_REF], following their earlier microscopic derivation in the theories of rods [START_REF] Rayleigh | The theory of sound, v. I[END_REF] and models of bubbly fluids [START_REF] Caflich | Wave propagation in bubbly liquid[END_REF][START_REF] Gavrilyuk | Mathematical and numerical modeling of two-phase compressible flows with micro-inertia[END_REF]. One of our main conclusions is that the quasicontinuum approximations of this type can be chosen to provide closed-form kinetic relations which are in excellent agreement with the predictions of the discrete theory for sufficiently fast moving phase boundaries. We show how this goal can be accomplished without sacrificing the well-posedness of an initial value problem, which has always been an an unpleasant signature of the straightforward quasicontinuum approximations based on truncated Taylor series.

The structure of the paper is as follows. In Section 2 we briefly review the phenomenological origin of kinetic relations in the continuum theory of dynamic phase transitions. We then introduce in Section 3 a prototypical discrete model which leads to a canonical advance-delay differential equation governing the dynamics of an isolated phase boundary. In Section 4 we present a nonstandard reformulation of the discrete problem as a continuous problem with strong spatial nonlocality of both kinetic and potential energies. Various quasicontinuum approximations based on either Taylor or Padé approximations of the corresponding kernels are presented in Section 5. In Section 6 we use Fourier transform to obtain a series of analytical solutions for the traveling waves in the discrete and quasicontinuum problems. In Section 7 we show how to extract from a fully conservative solution of the discrete/quasicontinuum problem a kinetic relation governing the rate of dissipation at the continuum scale. Examples of quasicontinuum approximations generating explicit kinetic relations are discussed in Section 8, where we also provide a detailed comparison of the resulting kinetics with the exact solution for the discrete theory. Finally, in Section 9 we discuss the issue of well-posedness of the truncated quasicontinuum approximations and offer a rather universal way of dealing with unboundedness of the approximate operators at short wave lengths. Our conclusions are presented in Section 10.

Continuum model.

To present our ideas in the simplest setting, consider longitudinal motions of a homogeneous elastic bar with the total energy

E = ρ u2 2 + φ(u x ) dx, (2.1)
where u(x, t) is the displacement field, u ≡ ∂u/∂t is the velocity, u x ≡ ∂u/∂x is the strain, ρ is the constant mass density, and φ(u x ) is the elastic energy density. When motion is smooth, the function u(x, t) satisfies the nonlinear wave equation

ρü = (σ(u x )) x , (2.2)
where σ(u x ) = φ (u x ) is the stress-strain relation. The parameters on a discontinuity moving with velocity V must satisfy both the classical Rankine-Hugoniot jump conditions

[[ u]] + V [[u x ]] = 0, ρV[[ u]] + [[σ(u x )]] = 0,
and the entropy inequality R = GV ≥ 0, where R is the rate of dissipation and

G = [[φ]] -{σ}[[u x ]] (2.3)
is the configurational (driving) force. Here [[f ]] ≡ f +f -denotes the jump and {f } ≡ (f + + f -)/2 -the average of f across the discontinuity. Unlike conventional shocks, subsonic discontinuities violate the Lax condition. The resulting nonuniqueness can be remedied by supplementing the above jump conditions by a kinetic relation specifying the dependence of the configurational force on the velocity of the phase boundary G = G(V ) [START_REF] Abeyaratne | A continuum model of a thermoelastic solid capable of undergoing phase transitions[END_REF][START_REF] Truskinovsky | Dynamics of nonequilibrium phase boundaries in a heat conducting elastic medium[END_REF]. Note that the continuum model provides no information about this function. One way to recover the kinetic relation is to consider a lattice model representing physically justified "dehomogenization" of the continuum model (2.2).

Discrete model.

As a micro-model for a bar undergoing phase transition, consider an infinite chain of particles, each interacting with q nearest neighbors on both sides. The total energy of the system may be written as

E = ε ∞ n=-∞ ρ u2 n 2 + q p=1 pφ p u n+p -u n pε , (3.1)
Here u n (t) is the displacement of the nth particle, un is the particle velocity, ε is the reference interparticle distance, ρ is the mass density and φ p (w) is the energy density of the interaction between pth nearest neighbors.

To represent martensitic phase transitions in the simplest way we may assume that the nearest-neighbor (NN) interactions are governed by a bi-parabolic potential representing an example of the simplest essential nonlinearity:

φ 1 (w) = ⎧ ⎨ ⎩ 1 2 Ψ(1)w 2 , w ≤ w c 1 2 Ψ(1)(w -a) 2 + aΨ(1) w c - a 2 , w ≥ w c , (3.2)
The potentials governing long-range interactions were assumed to be quadratic:

φ p (w) = 1 2 pΨ(p)w 2 , p = 2, ..., q. (3.3)
In terms of dimensionless variables t = t(Ψ(1)/ρ) 1/2 /ε, ūn = u n /(aε), Ψ(p) = Ψ(p)/Ψ(1), p = 1, . . . , q, but with bars dropped, the total energy can be rewritten as

E = ∞ n=-∞ u2 n 2 - 1 2 |k-n|≤q u n Ψ(k -n)u k -(w n -w c )θ(w n -w c ) . (3.4)
Here w n = u nu n-1 is the strain in NN springs, θ(w) is a unit step function, and we define

Ψ(0) = -2 q p=1 Ψ(p), Ψ(-p) = Ψ(p) (3.5)
The dynamic equations can be formulated for the strain variables w n :

ẅn

- |k-n|≤q Ψ(k -n)w k = 2θ(w n -w c ) -θ(w n+1 -w c ) -θ(w n-1 -w c ). (3.6)
Note that in the classical continuum limit the discrete model converges to the rescaled equations (2.1), (2.2) with ρ = 1 and the bilinear macroscopic stress given by

σ(w) = c 2 w -θ(w -w c ), (3.7)
where

c = q p=1 p 2 Ψ(p) 1/2 (3.8)
is the dimensionless macroscopic sound speed.

To model the motion of an isolated phase boundary, we now consider solutions of (3.6) in the form of the traveling waves w n (t) = w(ξ), ξ = n -V t, with all springs are in phase I (w < w c ) for ξ > 0 and in phase II (w > w c ) otherwise. For these solutions the system of equations (3.6) reduces to a single advance-delay differential equation for w(x):

V 2 w - |p|≤q Ψ(p)w(ξ + p) = 2θ(-ξ) -θ(-ξ -1) -θ(1 -ξ).
(3.9)

4. Reformulation of the problem. Observe that equation (3.9) may be written in the form

Λ D (w -w + ) = θ(-ξ), (4.1)
where the operator Λ D is defined by

Λ D ≡ -Δ -1 d V 2 Δ c - |p|≤q Ψ(p) exp(pD) (4.2)
Here D is the operator of differentiation, Dw = w ξ . We further defined the continuum Laplacian operator Δ c = D 2 and the shift operator exp(pD) acting according to exp(pD)w = w(ξ + p). Finally Δ d = exp(D) -2I + exp(-D) is the discrete Laplacian, defined by its action on w as Δ d w = w(ξ + 1) -2w(ξ) + w(ξ -1). The Fourier image of (4.2) is given by

Λ D (k, V ) = ω 2 (k) -V 2 k 2 4 sin 2 (k/2) , (4.3)
where

ω 2 (k) = 4 q p=1 Ψ(p) sin 2 pk 2 . (4.4)
is the dispersion relation of the discrete model. Following some insights from [START_REF]Hamiltonian dynamics of dense chains and lattices: or how to correct the continuum[END_REF], we introduce a sufficiently smooth continuum displacement field u(x, t) such that

Du(n, t) = (I -exp(-D))u n (t). (4.5)
This definition can be rewritten as u n (t) = Mu(n, t), where the operator M is defined by M = (Iexp(-D)) -1 D. In most derivations of quasicontinuum approximations it is assumed that M = I. To see the effect of this alternative (or "indirect", in the language of [START_REF]Hamiltonian dynamics of dense chains and lattices: or how to correct the continuum[END_REF]) definition of the continuum variable it is instructive to compute the kinetic energy. Using (4.5) and replacing the sum by an integral, we obtain

1 2 ∞ n=-∞ un (t) • un (t) ⇒ 1 2 ∞ -∞ Mu t (x, t) • Mu t (x, t)dx = 1 2 ∞ -∞ u t (x, t)Ku t (x, t)dx (4.6)
Here the operator K is the continuum analog of the operator Δ -1 d Δ c , in the sense that their Fourier images coincide:

K(k) = F[Δ -1 d Δ c ] = k 2 4 sin 2 (k/2)
.

One can see that contrary to the more traditional approach (e.g. [START_REF] Kunin | Elastic Media with Microstructure I: One-Dimensional Models[END_REF]), in our continuum variables the kinetic energy is nonlocal. By rewriting in the new variables the whole discrete energy (3.4) we obtain

E QC = ∞ -∞ 1 2 u t Ku t + 1 2 u x Lu x -(u x -w c )θ(u x -w c ) dx, (4.7)
where the operator L is the continuum analog of Δ -1

d |p|≤q Ψ(p) exp(pD), in the sense that L(k) = F[Δ -1 d |p|≤q Ψ(p) exp(pD)] = q p=1 Ψ(p) sin 2 (pk/2) sin 2 (k/2) .
The equation of motion generated by (4.7) is

∂ ∂t [Ku t ] = ∂ ∂x [Lu x -θ(u x -w c )].
Applying the traveling wave ansatz u(x, t) = ũ(ξ), ξ = x -V t and integrating the resulting equation, we obtain

(L -V 2 K)w -θ(w -w c ) = const,
where w(ξ) = ũ (ξ) is the strain. Now, assuming as before that θ(ww c ) = θ(-ξ) and applying the boundary condition w → w + as ξ → ∞, we arrive at

(L -V 2 K)(w -w + ) = θ(-ξ), (4.8)
Note that (4.8) is the exact analog of (4.1) due to the identity

Λ D (k, V ) = L(k) -V 2 K(k) (4.9)
The continuum model with the energy (4.7) is equivalent to the discrete model with the energy (3.4). Notice that in contrast to the quasicontinuum models proposed in [START_REF] Kunin | Elastic Media with Microstructure I: One-Dimensional Models[END_REF], which preserve the structure of the kinetic energy while modifying only the elastic energy, we propose to modify expressions for both kinetic and elastic energies. This reflects the physical understanding that the microscopic lattice-scale vibrations cannot be treated at the mesoscale exclusively as thermodynamic free energy. Here we follow a tradition which dates back to Rayleigh who in his studies of vibrating beams understood the necessity of incorporating the concept of transversal inertia [START_REF] Rayleigh | The theory of sound, v. I[END_REF]. Similar ideas appeared in many areas, for instance in the homogenized theory of bubbly fluids where the extra terms in the kinetic energy describe inertia of individual bubbles [START_REF] Caflich | Wave propagation in bubbly liquid[END_REF][START_REF] Gavrilyuk | Mathematical and numerical modeling of two-phase compressible flows with micro-inertia[END_REF].

Quasicontinuum approximations.

The quadratic part of the nonlocal functional (4.7) can be approximated by local functionals in several ways. This gives rise to significantly different quasicontinuum theories.

Strain-gradient approximation.

The most common way of approximating discrete energy functional is based on retention of only the non-dispersive continuum expression for the kinetic energy. This means that K = Δ -1 d Δ c is replaced by I, which is equivalent to replacing the discrete Laplacian Δ d by Δ c . For consistency the operator L should then be replaced by Δ -1 c |p|≤q Ψ(p) exp(pD). Expanding the Fourier image of the resulting approximate operator in Taylor series around k = 0 and retaining the first few terms, we obtain

K(k) ≈ 1, L(k) ≈ c 2 + α 1 k 2 + α 2 k 4 , (5.1)
where

α 1 = - 1 12 q p=1 p 4 Ψ(p), α 2 = 1 360 q p=1 p 6 Ψ(p). (5.2)
Observe that in the physical space the associated quadratic part of potential energy reduces to

1 2 ∞ -∞ u x Lu x dx ≈ 1 2 ∞ -∞ [c 2 u 2 x + α 1 u 2 xx + α 2 u 2 xxx ]dx.
This approximation is often referred to as the strain-gradient (SG) model (e.g. [START_REF] Mindlin | Second gradient of strain and surface tension in linear elasticity[END_REF][START_REF] Triantafyllidis | On higher order gradient continuum theories in 1-D nonlinear elasticity. Derivation from and comparison to the corresponding discrete models[END_REF]). Its advantage is relative simplicity. However, the Fourier image of the approximate operator agrees with Λ D (k, V ) only up to O(k 2 ). This order of approximation is shared even by strongly nonlocal quasicontinuum models that are focused exclusively on the approximation of the elastic energy (e.g. [START_REF] Kunin | Elastic Media with Microstructure I: One-Dimensional Models[END_REF][START_REF] Rogula | Introduction to nonlocal theory of material media[END_REF]).

Mixed gradient approximation.

To obtain a higher order of approximation one must include the dispersive corrections in both kinetic and elastic energies. For instance, by expanding the Fourier images of the two operators K and L to the same order in Tailor series we obtain

K(k) ≈ 1 + k 2 12 + k 4 240 , L(k) ≈ c 2 + a 1 k 2 + a 2 k 4 , (5.3)
where

a 1 = 1 12 q p=2 (p 2 -p 4 )Ψ(p), a 2 = 1 720 q p=2 p 2 (3 -5p 2 + 2p 4 )Ψ(p).
(5.4)

In this approximation the kinetic energy and the quadratic part of the potential energy take the form

1 2 ∞ -∞ u t Ku t dx ≈ 1 2 ∞ -∞ u 2 t + 1 12 u 2 xt + 1 240 u 2 xxt dx, 1 2 ∞ -∞ u x Lu x dx ≈ 1 2 ∞ -∞ c 2 u 2 x + a 1 u 2 xx + a 2 u 2 xxx dx.
(5.5)

We will refer to this class of models as mixed gradient (MG) models. First approximation of this type with only quadratic term in k preserved in the kinetic energy and only zero-order term retained in the elastic energy was suggested by Rosenau [START_REF] Rosenau | Dynamics of nonlinear mass-spring chains near the continuum limit[END_REF]; more recent attempts along the same lines can be found in [START_REF] Kevrekidis | Continuum approach to discreteness[END_REF][START_REF]Hamiltonian dynamics of dense chains and lattices: or how to correct the continuum[END_REF][START_REF] Sharma | The kinetic relation of a Peierls dislocation in a higher-gradient dispersive continuum[END_REF][START_REF] Wattis | Approximations to solitary waves on lattices, III: the monatomic lattice with second-neighbour interactions[END_REF].

The advantage of the Rosenau model is that it is uniformly well-posed in the sense that it avoids unphysical short-wave instabilities that are characteristic for generic truncated polynomial approximations. This issue will be addressed in more detail in Section 9, while here we simply mention that by well-posedness we mean the discrete analog of hyperbolicity of the linear operator of the model in the sense that ω 2 (k) = k 2 L(k)/K(k) is positive for all nonzero k (see [START_REF] Kreiss | Time-dependent partial differential equations and their numerical solution[END_REF] for more detail).

To construct an improved Rosenau approximation with higher resolution we assume that L(k) and K(k) are still approximated by polynomials. Then the MG model leads to a well-posed problem if (i) the approximation of L(k) has no real zeroes and (ii) k 2 L(k) and K(k) are approximated by polynomials of the same order. This leads to the next simplest Rosenau-type well-posed model

K(k) ≈ 1 + k 2 12 + k 4 240 , L(k) ≈ c 2 + a 1 k 2 (5.6)
with a 1 as in (5.4). Both of the above requirements are satisfied provided that a 1 > 0, so that c 2 + a 1 k 2 = 0 has no real roots. Note that unlike the general MG model (5.3), where both kinetic and elastic operators are expanded up to O(k 6 ), this model approximates the discrete operator only up to O(k 4 ) because of lower order accuracy in approximating L(k). In what follows we refer to the model (5.6) as MG1/2 model because of its intermediate character.

Mixed gradient-Padé approximation.

Another possibility to preserve the well-posedness of the linear problem is to approximate L(k) by a (2, 2) Padé (rational) approximation accurate up to O(k 6 ), while keeping for K(k) a quartic polynomial expansion:

K(k) ≈ 1 + k 2 12 , L(k) ≈ c 2 + b 1 k 2 1 + b 2 k 2 , (5.7)
with

b 1 = 1 12 q p=2 (p 2 -p 4 )Ψ(p) - c 2 q p=2 p 2 (3 -5p 2 + 2p 4 )Ψ(p) 60 q p=2 (p 2 -p 4 )Ψ(p) , b 2 = - q p=2 p 2 (3 -5p 2 + 2p 4 )Ψ(p) 60 q p=2 (p 2 -p 4 )Ψ(p)
.

(5.8)

We will refer to this approximation as MGP model. In this case the dispersion relation (see (9.1)) remains to be the ratio of two polynomials of the same order, and if b 1 > 0, the MGP model leads to a well-posed problem. Note that like the well-posed version of the MG model (5.6), this approximation is accurate only up to O(k 4 ). In this case it is due to the reduced accuracy in the expansion of K(k).

6. Traveling waves. The traveling wave equation for both discrete and quasicontinuum models can be written as

Λ(w -w + ) = θ(-ξ), (6.1)
where

Λ = Λ D (k, V ) (6.2)
for the discrete model,

Λ = (c 2 -V 2 )I -α 1 D 2 + α 2 D 4 (6.3)
for the SG model,

Λ = (c 2 -V 2 )I -a 1 - V 2 12 D 2 + a 2 - V 2 240 D 4 (6.4)
for the MG model and

Λ = (I -b 2 D 2 ) -1 (c 2 I -b 1 D 2 ) -V 2 1 - 1 12 D 2 (6.5)
for the MGP model. For consistency with the assumptions that led to (6.1), we must also require that solution of equation (6.1) satisfies the following constraints:

w(0) = w c (6.6)
and w(ξ) < w c for ξ > 0, w(ξ) > w c for ξ < 0 (6.7)

The average states at ξ = ±∞ must correspond to stable constant-strain equilibria w(ξ) ≡ lim We solve (6.1) by writing w(ξ) = h(ξ) + w + and applying the complex Fourier transform ĥ

(k) = ∞ -∞ h(ξ)e i(k-iα)ξ dξ, h(ξ) = 1 2π ∞-iα -∞-iα ĥ(k)e -ikξ dk,
where α > 0 is a small parameter which guarantees convergence of the integrals.

After inverting the Fourier transform and letting α → 0, we obtain the integral representation

w(ξ) = w + - 1 2πi Γ e ikξ kΛ(k, V ) dk, (6.9)
where the contour Γ coincides with the real axis everywhere except near the singular points. Specifically, it goes above k = 0 and, in order to satisfy the radiation condition, passes below the nonzero real roots of Λ(k, V ) = 0 if the group velocity exceeds V and above otherwise (in this case the radiation condition reduces to requiring that the contour goes below a nonzero real root if kΛ k (k, V ) > 0 and above otherwise).

To compute the integral (6.9) explicitly, we use the residue method closing the contour in the upper half-plane when ξ > 0 and in the lower half-plane when ξ < 0.

w(ξ) = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ w -+ k∈M -(V ) e ikξ kΛ k (k, V ) for ξ < 0 w + - k∈M + (V ) e ikξ kΛ k (k, V ) for ξ > 0, (6.10)
The summations in (6.10) are over the sets of roots of Λ(k, V ) = 0 defined as

M ± (V ) = {k : Λ(k, V ) = 0, Imk ≷ 0} N ± (V ),
where

N ± (V ) = {k : Λ(k, V ) = 0, Imk = 0, kΛ k (k, V ) ≷ 0}
are the sets of real roots responsible for radiative modes whose amplitude does not decrease at infinity. These waves are placed either behind or in front of the moving phase boundary according to the radiation condition.

Observe that by construction all plane waves with real wave numbers k involved in the traveling wave solution (6.1) have real frequency ω(k) = V k, so they lie inside the region in Fourier space where the models are well-posed.

7. Kinetic relation. The parameters entering (6.10) are not all independent. The average strains at ±∞ must satisfy a macroscopic Rankine-Hugoniot condition

w + = w -- 1 c 2 -V 2 , (7.1)
where we recall that c is the dimensionless macroscopic sound speed defined in (3.8).

Note that since we assumed that w -< w + , the solutions are necessarily subsonic (V < c). The continuity of strain at ξ = 0 and (6.6) imply that

w ± = w c ∓ 1 2(c 2 -V 2 ) + k∈N pos (V ) 1 |kΛ k (k, V )| , (7.2)
where N pos (V ) = {k ∈ N + (V ) N -(V ) : k > 0} is the set of all positive real roots of Λ(k, V ) = 0. Computing microscopic energy fluxes carried by radiative modes to ±∞ (see [START_REF] Kresse | Mobility of lattice defects: discrete and continuum approaches[END_REF][START_REF]Kinetics of martensitic phase transitions: Lattice model[END_REF] for more detail), we obtain the following expression for the driving force:

G(V ) = k∈N pos (V ) 1 |kΛ k (k, V )| (7.3)
Since the sets N pos depend on V , (7.3) furnishes, although not in an explicit form, the desired kinetic relation between the driving force and the velocity of the moving boundary. Note that formula (7.3) can also be obtained from the macroscopic assessment of dissipation (2.3) yielding in the bilinear case

G = 1 2 (w -+ w + ) -w c . (7.4)
Indeed, substituting (7.2) in (7.4), we recover (7.3).

Examples.

To illustrate the general analysis presented in the previous sections, below we focus on a specific case of a lattice with q = 2 as the simplest nontrivial example of a lattice model accounting for bilinear local (nearest neighbor, or NN) and linear nonlocal (next to nearest neighbor, or NNN) interactions. In the nondimensional variables the model is fully characterized by a single dimensionless parameter β = 4Ψ(2)/Ψ(1) measuring the relative strength of NNN interactions. The total energy of the system can be written as

E = ∞ n=-∞ u2 n 2 + 1 + β 2 w 2 n -θ(w n -w c )(w n -w c ) - β 4 (w n+1 -w n ) 2 . (8.1)
The macroscopic sound speed (3.8) reduces in this case to

c = (1 + β) 1/2 . (8.2)
For this discrete model one can obtain explicit traveling wave solutions describing the motion of an isolated phase boundary and present the kinetic relation in the form of a numerical procedure generating the function G = G(V, β) [START_REF]Kinetics of martensitic phase transitions: Lattice model[END_REF].

Our task in what follows is to choose a quasicontinuum approximation which approximates this function best. Following [START_REF]Kinetics of martensitic phase transitions: Lattice model[END_REF], we limit our analysis to the interval

β 1 < β ≤ 0, (8.3)
where β 1 ≈ -0.9539. The lower bound ensures the existence of a fast subsonic traveling wave solution with monotone leading edge and a single radiative mode behind the phase boundary, which will be the prime target of our modeling, while the upper bound is suggested by the linearization of the potentials of the Lennard-Jones type [START_REF] Truskinovsky | Peierls-Nabarro landscape for martensitic phase transitions[END_REF][START_REF]The origin of nucleation peak in transformational plasticity[END_REF]. It is known that the kinetic relations generated by fully inertial lattice models are quite complex in the small velocity range due to a large number of lattice resonances [START_REF] Kresse | Mobility of lattice defects: discrete and continuum approaches[END_REF][START_REF]Kinetics of martensitic phase transitions: Lattice model[END_REF]. Since the stability of a steadily moving isolated front in this regime is questionable (e.g. [START_REF]Kinetics of martensitic phase transitions: Lattice model[END_REF]), we intentionally leave it outside our quasicontinuum modeling.

Quadratic expansions.

We begin by considering the simplest quadratic expansions in both SG and MG models. This corresponds to setting α 2 = 0 in (5.1) and a 2 = 0 in (5.3); the first model will be called SG1 and the second one MG1. In both cases equation Λ(k, V ) = 0 has two symmetric roots which are either real, k = ±r, or purely imaginary, k = ±is. When the roots are real the solution takes the form

w(ξ) = ⎧ ⎨ ⎩ w c + 1 -cos(rξ) c 2 -V 2 for ξ < 0 w c for ξ > 0. (8.4)
and the kinetic relation then reduces to

G = 1 2(c 2 -V 2 ) . (8.5)
In the case of purely imaginary roots the (subsonic) solution takes the form

w(ξ) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ w c + 1 -e sξ 2(c 2 -V 2 )
for ξ < 0

w c - 1 -e -sξ 2(c 2 -V 2 )
for ξ > 0.

(8.6)

The corresponding kinetic relation simply states that the dissipation is absent:

G = 0. (8.7)
To be more specific, consider first the quadratic strain-gradient approximation (SG1 model). In this case we have

Λ(k, V ) = c 2 -V 2 - 1 + 4β 12 k 2 . (8.8)
It is easy to see that at V < c the roots of Λ(k, V ) = 0 are real at -1/4 < β ≤ 0 and imaginary otherwise. This implies that SG1 model predicts zero dissipation for β < -1/4 and is therefore unable to capture the phenomenon of radiative damping predicted in this parameter range by the lattice model. We observe that it is exactly the region where SG1 model is usually used (e.g. [START_REF] Triantafyllidis | On higher order gradient continuum theories in 1-D nonlinear elasticity. Derivation from and comparison to the corresponding discrete models[END_REF]) since it is uniformly well-posed there.

In the rest of the domain -1/4 < β ≤ 0 the SG1 model predicts a nontrivial kinetic relation, although it can only be used if the real wavenumbers are in the interval |k| ≤ k * = c 12/(1 + 4β). In this case there are two real roots of the dispersion relation and the formulas (8.4), (8.5) provide the solution of the problem. In the MG1 approximation the dispersion relation takes the form

Λ(k, V ) = c 2 -V 2 - β 4 + V 2 12 k 2 . (8.9)
In the interval -1/4 < β ≤ 0 the MG1 model gives the same kinetics as SG1 model for (-3β) 1/2 < V ≤ c and predicts zero dissipation for V < (-3β) 1/2 . Like SG1 model, it is unable to capture radiative damping for β < -1/4 and predicts zero dissipation.

8.2.

Comparison with discrete theory. The results of the numerical evaluation of the solutions of the quasicontinuum models are summarized in Figure 8.1, where SG1 (dotted curves) and MG1 (dashed) approximations are compared to the discrete model in the interval -1/4 < β ≤ 0. One can see that both SG1 and MG1 model overestimate the frequency and the amplitude of the radiation behind the front in the near-sonic domain (Figure 8.1c). In addition, the approximations require the strain in front (ξ > 0) to be constant and equal to the critical value w c , while in the discrete model the strain exponentially decreases to a finite value w + < w c . Although at velocities sufficiently close to sonic limit the models capture the qualitative character of the discrete kinetic relation (see Figure 8.1d), they overestimate the driving force. Similar comparison in the interval β < -1/4 does not make much sense because both quasicontinuum models predict zero dissipation while the discrete theory generate a nontrivial kinetic relation (see Fig 8.3). Overall, one can conclude that none of the quadratic models capture the discrete kinetics adequately in the whole range of parameters. To achieve a better approximation one needs to include higher order terms in the long-wave expansion, which is done in the next subsection.

Quartic approximations.

If we keep fourth-order terms in our Taylor expansions, we obtain what we call SG2 model (5.1) and MG2 model (5.3). Assume as before (8.3) and subsonic V . With the exception of the interval -1/16 < β ≤ 0, where SG2 model leads to four real roots, in both models the characteristic equation, Λ(k, V ) = 0 has two real roots, k 1,2 = ±r, and two imaginary roots, k 3,4 = ±ip. The traveling wave solution (6.1) at subsonic velocities then reduces to

w(ξ) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ w c + 2p 2 (1 -cos(rξ)) + r 2 (1 -e pξ ) 2(c 2 -V 2 )(p 2 + r 2 ) for ξ < 0 w c + r 2 (e -pξ -1) 2(c 2 -V 2 )(p 2 + r 2 )
for ξ > 0 (8.10)

The kinetic relation (7.3) takes the form

G = p 2 2(c 2 -V 2 )(p 2 + r 2 ) . (8.11)
To specialize these formulas to SG2 model, we observe that in our example

Λ(k, V ) = 1 + β -V 2 - 1 + 4β 12 k 2 + 1 + 16β 360 k 4 = 0, (8.12)
If -1/16 < β ≤ 0 the characteristic equation has four real roots k = ±r 1 and k = ±r 2 , r 1 < r 2 . The roots k = ±r 1 approximate the discrete radiation, whereas the roots k = ±r 2 are an artifact of the polynomial expansion with a positive higher-order coefficient and do not have a discrete analog. Thus, SG2 model cannot be used in this interval. At β = -1/16 the model reduces to SG1 approximation. Assume now that β 1 < β < -1/16. By solving the dispersion relation we obtain two imaginary roots k = ±ip and two real roots k = ±r with the magnitudes

p = - 180 1 + 16β 1 + 4β 12 + 1 + 4β 12 2 - 1 + 16β 90 (1 + β -V 2 ) , r = - 180 1 + 16β 1 + 4β 12 2 - 1 + 16β 90 (1 + β -V 2 ) - 1 + 4β 12 .
Observe that this approximation is valid (V is real) for real k satisfying |k| ≤ k * , where

k * = 15 + 60β -3 √ 15 -1 -32β -16β 2 1 + 16β . (8.13)
The explicit kinetic relation in the SG2 approximation takes the form

G(V ) = 1 4(c 2 -V 2 ) 1 -1 - 8(1 + 16β) 5(1 + 4β) 2 (c 2 -V 2 ) -1/2 (8.14)
To obtain the kinetic relation in the mixed gradient approximation we first adopt it to our example and get

Λ = c 2 -V 2 - β 4 + V 2 12 k 2 + β 48 - V 2 240 k 4 (8.15)
In this case β ≤ 0 ensures that the characteristic equation always has two imaginary and two real roots. Their magnitudes are given by the formulas

p = √ 2 -(15β + 5V 2 ) --75β(4 + β) + 30(2 + 17β)V 2 -35V 4 5β -V 2 , r = √ 2 15β + 5V 2 --75β(4 + β) + 30(2 + 17β)V 2 -35V 4 5β -V 2
In this model the cutoff wavenumber is

k * = 6 + 2 √ 3(4 + β) -β(4 + β)
.

Note that it is greater than k * for SG2 model given in (8.13) for all β where the latter is real, so that the domain of applicability of the mixed gradient approximation is wider than that of the strain gradient approximation. The kinetic relation (7.3) is given in the mixed gradient approximation by the following explicit formula are compared to each other and the discrete kinetics in Figure 8.2d. As expected, the kinetic curve for MG2 model (dashed curve) is much closer to the discrete kinetic relation (solid) than the SG2 approximation (dotted curve).

G(V ) = 1 4(c 2 -V 2 ) 1 + √ 5(3β + V 2 ) -15β(4 + β) + 6(2 + 17β)V 2 -7V 4 . (8.16)
Next consider the interval β 1 < β < -1/4. In this parameter range the positive real root corresponding to radiation behind the front is not close to k = 0, so that both SG2 and MG2 approximations cannot be expected to capture the discrete kinetics quantitatively. Nevertheless, as Figure 8.3 illustrates, MG2 model (dashed curves) provides a substantially better approximation of the discrete model (solid curves) than the SG2 model (dotted curves) in terms of both the kinetic relations G(V ) and the near-sonic traveling wave profiles. Note in particular that the MG2 kinetic curve stays close to the discrete fast branch over the whole velocity range, while the SG2 curve deviates away from it significantly as V approaches c. Although MG2 approximation naturally becomes worse for higher |β| (the roots corresponding to radiation move further away from zero), it remains significantly closer to the discrete curve than the SG2 approximation.

The better performance of MG2 model is due to the fact that it directly approximates Λ D (k, V ) near k = 0 together with its first two derivatives, while SG2 model effectively replaces 4 sin 2 (k/2) in (4.3) by k 2 before expanding it in Taylor series. Even though the Taylor expansion is no longer valid at larger k, the optimal "tuning" of the function and its derivatives near k = 0 apparently results in a better agreement with the discrete theory in a broad range of wave lengths.

Next consider the well-posed version of the MG model which we referred to as MG1/2 and MGP approximations. The corresponding roots of the dispersion relations 2 and β = -0.5, respectively, which also compares them with the roots for the discrete problem. One can see that at smaller |β| the MG1/2 model gives a slightly better prediction of the real roots than the MGP model, while at larger |β| the opposite is true. The reason for this becomes clear when one observes that the error of MG1/2 approximation is -βk 4 /48 + O(k 6 ) and therefore the quality of approximation is getting better at β gets closer to zero. Meanwhile, the error for MGP model is independent of β (to the lowest order) but depends instead on the velocity V 2 k 4 /240 + O(k 6 ), so this approximation improves at smaller velocities for |k| close to zero. Figures 8.4c and 8.5c compare the strain profiles in both models to the discrete case. One can see again that MG1/2 model works better at smaller β, while MGP model gives a better prediction of the discrete strain profile for larger β. Both models overestimate the frequency of radiation.

The kinetic relations generated by both models are compared with the discrete kinetics from [START_REF]Kinetics of martensitic phase transitions: Lattice model[END_REF] in Figures 8.4d and 8.5d. As expected, kinetic relation for MG1/2 model is closer to the discrete one at smaller β (notice, however, that unlike the discrete model [START_REF]Kinetics of martensitic phase transitions: Lattice model[END_REF] it predicts G(0) = 0), and the MGP model gives a better prediction at higher |β|. Quite expectedly, both models gives a worse approximation than MG2 model which is of higher order, but their advantage is that they are linearly stable at all wavelengths. In the next section we show how the best so far MG2 approximation can be adjusted to enjoy the same level of well-posedness. 9. Well-posed approximations. In constructing our quasicontinuum approximations we were guided so far exclusively by the desire to provide an optimal fit for the discrete dispersion relation at long waves. In order to be able to use the model outside the narrow framework of the traveling wave solutions, one needs to impose an additional condition that the model produces well-posed initial value problem. At least necessary conditions can be formulated quite easily. Indeed, we recall that the governing equations for all the approximate models considered so far are linear in each phase. Therefore to test the well-posedness of an initial value problem in the pure phases it is sufficient to consider the plane waves of the form exp[i(ωt + kx)]. The dispersion relation for a generic quasicontinuum approximation can be written as

ω 2 (k) = k 2 L A (k) K A (k) , (9.1)
where L A (k) and K A (k) are the corresponding approximations of L(k) and K(k).

The wavenumbers k at which ω 2 (k) < 0 correspond to modes of instability for the homogeneous states away from the front. In view of (9.1) and the fact that K A (k) > 0 in all approximations considered above, these are the wavenumbers at which L A (k) < 0. For instance, if the coefficients in front of the highest-order term in L A (k) are negative, homogeneous deformation is unstable with respect to perturbations with sufficiently short wave length, |k| > k * , where L A (k * ) = 0. As one can readily see, this is indeed what happens with our "best" polynomial approximation (MG2 model) in at least one of the relevant parameter ranges. Although this formally implies that this model is ill-posed, one has to keep in mind that the polynomial expansions are expected to be valid only for sufficiently long waves. The apparent short-wave instability is the result of the rigid, low parametric nature of the approximation and is unphysical given that the homogeneous deformation is stable in the original discrete problem. To eliminate this unphysical instability, the approximate models must be restricted to functions whose Fourier image is confined in the interval |k| < k * and the corresponding projection of the operator leads to a well-posed problem. In statics this restriction can be achieved by setting ω(k) = 0 for real wave numbers that satisfy |k| > π and thus correspond to length scales less than the lattice size [START_REF] Kunin | Elastic Media with Microstructure I: One-Dimensional Models[END_REF][START_REF]Quasicontinuum modeling of short-wave instabilities in crystal lattices[END_REF].

In dynamics we can similarly replace L(k) by

L ε (k) = L A (k)η ε (k/k * ), (9.2)
where ε > 0 is a small parameter, k * satisfies L A (k) = 0 and

η ε (y) = 1 2 1 -tanh y 2 -1 ε (9.3)
is a smooth and analytic cut-off function, which is positive for all real y. Notice that for sufficiently small ε this function is close to one for |y| < 1 and decays exponentially fast to zero outside this interval. Thus for |k| > k * the function ω 2 ε (k) = k 2 L ε (k)/K(k) is negative but approaches zero exponentially fast at large k. It tends to zero together with ε, so that the growth rate for short-wave instabilities also disappears in the limit. For real k one can show that η ε (k/k * ) → θ(1k 2 /k 2 * ) as ε tends to zero, effectively cutting off the unwanted wave-numbers. For instance, applying this procedure to the MG2 model we obtain in the limit ε → 0 a nonlocal energy operator:

L(k) ≈ (c 2 + a 1 k 2 + a 2 k 4 )θ(1 -k 2 /k 2 * ) (9.4)
In this case the quadratic part of the elastic energy takes the form is the oscillatory kernel that decays as 1/x 2 at infinity. We remark that at sufficiently long waves the nonlocal approximation (9.4) of the potential energy operator is equivalent to the local gradient approximation (5.3) 2 . The two models are expected to be close in at least some nonlinear regimes as well; in particular, they generate identical kinetic relations for all nonzero V .

1 2 ∞ -∞ u x Lu x dx ≈ 1 2 ∞ -∞ ∞ -∞ Φ(x -

Conclusions.

In this paper we proposed a series of quasicontinuum approximations for the simplest one-dimensional lattice model describing dynamic phase transitions. The approximations are dispersive and include various non-classical corrections to both kinetic and potential energies. We showed that the quasicontinuum theory can be chosen in such a way that the associated closed-form kinetic relation is in excellent agreement with the predictions of the discrete theory for a broad range of phase boundary velocities. The essential ingredient making our approximations adequate in a wide range of parameters is the dispersive modification of the kinetic energy which introduces in the constitutive model the effects of micro-inertia. The presence of somewhat unusual kinetic energy terms at the mesoscale, where quasicontinuum models actually operate, reflects the oscillatory nature of the typical solutions in the discrete phase transition problem. To deal with the ill-posedness of the straightforward high-gradient quasicontinuum expansions we proposed to use a cutoff in the spectrum of the linearized operator. This operation, detailed in the paper, furnishes a canonical replacement of an ill-posed high-gradient model with a well-posed strongly nonlocal model with an explicit kernel. An unusual feature of this kernel is that the strong nonlocality of the model remains invisible in the case of long wave initial data until the nonlinearity is able to generate sufficiently short waves. The fact that this does not happen in the solutions describing isolated phase boundaries gives a hope that in most cases of interest one should be able to safely replace integral equations by partial differential equations.

  dζ → w ± , as ξ → ±∞.(6.8)

Fig. 8 . 1 .

 81 Fig. 8.1. Comparison of discrete (solid curves), MG1 (dashed curves) and SG1 (dotted curves) models at -1/4 < β ≤ 0: the structure of (a) real and (b) imaginary roots; (c) strain profiles at V /c = 0.92, w c = 1; (d) the mobility curves. In all pictures β = -0.2.

Fig. 8 . 2 .

 82 Fig. 8.2. Comparison of discrete (solid curves), MG2 (dashed curves) and SG2 (dotted curves) models at -1/4 < β < -1/16: the structure of (a) real and (b) imaginary roots; (c) strain profiles at V /c = 0.85, w c = 1; (d) the mobility curves. In all pictures β = -0.2.

Fig. 8 . 3 .

 83 Fig. 8.3. Comparison of discrete (solid curves), MG2 (dashed curves) and SG2 (dotted curves) models at β 1 < β < -1/4: the structure of (a) real and (b) imaginary roots; (c) strain profiles at V /c = 0.85, w c = 1; (d) the mobility curves. In all pictures β = -0.5.

Fig. 8 . 4 .

 84 Fig. 8.4. Comparison of discrete (solid curves), MG1/2 (dashed curves) and MGP (dotted curves) models at -1/4 < β ≤ 0: the structure of (a) real and (b) imaginary roots; (c) strain profiles at V /c = 0.85, w c = 1; (d) the mobility curves. In all pictures β = -0.2.

Fig. 8 . 5 .

 85 Fig. 8.5. Comparison of discrete (solid curves), MG1/2 (dashed curves) and MGP (dotted curves) models at β 1 < β < -1/4: the structure of (a) real and (b) imaginary roots; (c) strain profiles at V /c = 0.85, w c = 1; (d) the mobility curves. In all pictures β = -0.5.

  y)u x (x, t)u x (y, t)dxdy, x((a 1 + 2a 2 k 2 * )x 2 -12a 2 ) cos(k * x) + (24a 2 -2(a 1 + 6a 2 k 2 * )x 2 ) sin(k * x)

	where
	Φ(x) =	1 πx 5 2k

*

Acknowledgements. This work was supported by the National Science Foundation grants DMS-0102841 (L.T.) and DMS-0137634 (A.V.).

Next consider the well-posed MG1/2 and MGP models defined in (5.6) and (5.7), respectively. For the MG1/2 model the characteristic equation is given by

In MGP approximation we obtain

The characteristic equation in both models has two real and two imaginary roots with

for MG1/2 model and

for MGP model. The kinetic equations are again explicit taking the form

for MG1/2 model and

for MGP model.

Comparison with discrete theory.

We can now compute the solutions in quartic models explicitly and compare the approximate kinetic relations with the exact one for the discrete model (see [START_REF]Kinetics of martensitic phase transitions: Lattice model[END_REF]). The first series of results is summarized in Figure 8.2a,b which compares in the interval -1/4 < β < -1/16 the following three models: discrete (solid curves), MG2 (dashed) and SG2 (dotted). One can see that MG2 model gives a much better approximation of the both real and imaginary roots. Near-sonic strain profiles for all three models are juxtaposed in Figure 8.2c. The profiles for discrete (solid curve) and MG2 (dashed) models nearly coincide, while the SG2 approximation (dotted) gives a noticeably larger amplitude and smaller frequency of radiation behind the front. The kinetic relations for both approximations


François Comte
email: comte@lms.polytechnique.fr

Habibou Maitournam
email: habibou@lms.polytechnique.fr

Pierre Burry
email: pierre.burry@mpsa.com

T Mac

Lan Nguyen
email: maclan.nguyen@mpsa.com

A direct method for the solution of evolution problems

Keywords: Computational solid mechanics, Solids and structures, Large increment, Wavelets Mécanique des solides numérique, Solides et structures, Grand incrément, Ondelettes

An alternative numerical method to solve structural evolution problems is presented. It belongs to the family of Large Time Increment Methods. The problem is solved on the whole time interval separating a global stage where a SA and KA solution is searched and a local stage where the constitutive law is satisfied. To solve the global stage, consuming the most CPU time, the mechanical fields are decomposed in a wavelet basis and the equilibrium is solved only for the largest coefficients.

Résumé

La méthode directe. Cette Note présente une méthode numérique alternative de résolution de problèmes transitoires. Elle s'inscrit dans le cadre des méthodes à Grand Incrément de Temps ; la solution est recherchée sur la totalité du trajet de chargement en résolvant alternativement un problème dit global où l'on recherche une solution SA et CA et un problème local où on satisfait la loi de comportement. L'étape globale étant la plus coûteuse en temps de calcul, sa résolution se fait en décomposant le chargement sur une base d'ondelettes et en résolvant l'équilibre uniquement pour les coefficients d'ondelettes les plus importants.

Version française abrégée

Le travail présenté porte sur un nouvel algorithme de résolution numérique développé dans le cadre de la simulation d'un procédé de traitement thermique pour les culasses en aluminium. La prise en compte du process, et par conséquent la détermination de l'état mécanique résiduel, est importante pour prédire la tenue en service de la culasse.

Introduction

Nowadays, numerical modeling is mandatory in process development and validation of automotive components. Because of the large number of degrees of freedom of the F.E. models (about one million for a cylinder head), the complex loading and the nonlinear constitutive laws, calculation time is often prohibitive. In that case, computations cannot been used in an optimisation process.

The purpose of this Note is to propose a new numerical method to reduce calculation time. This method is based on the Large Time Increment Method and the wavelet decomposition for the equilibrium solution.

Problem position

Let us consider the quasi-static evolution, on a time interval [0, T ], of a structure made from a non-linear material. From its initial state, the structure is subjected at each time t to volumic body force f (x, t) on its volume Ω and to surface forces T d on ∂Ω T , part of the structure surface ∂Ω. On surface ∂Ω u , complementary to ∂Ω T , displacements u(x, t) are imposed. A temperature field T (x, t), result of a previous thermal calculation, is also given. The evolution satisfies the hypothesis of infinitesimal transformations. For simplicity, the material is assumed to be thermoelastoplastic.

The solution of the problem, defined on [0, T ] × Ω, consists of the displacements u, strains ε, plastic strains ε p and stresses σ , whose evolutions are obtained on [0, T ] × Ω. σ is statically admissible (SA), u kinematically admissible (KA) and ε, ε p and σ are related through the constitutive equations (PA).

Principle of the solution algorithm

From a general point of view, the approach is based on two main methods: the Large Time Increment Method developed by Ladevèze [START_REF] Ladevèze | Nonlinear Computational Structural Methods: New Approaches and Non Incremental Methods of Calculation[END_REF] and the Direct Cyclic Method by Maitournam et al. [START_REF] Maitournam | Détermination de la réponse asymptotique d'une structure anélastique sous chargement thermomécanique cyclique[END_REF]. From the first one, we retain the fact that the problem is directly treated on the whole time interval by separating a linear global stage where a statically and kinematically admissible (SKA) response is looked for and a local stage where a plastically admissible (PA) solution is obtained. This solution scheme accelerates the solution convergence. The second suggests that in the global stage, consuming the most CPU time, the loading is decomposed in a specific basis. In this method, suitable for cyclic loadings, since a periodic limit response of the structure is sought, Fourier decomposition is used. For the general evolution problem, where the solution is not periodic, the Fourier expansion is not appropriate and a better projection has to be found.

The iterative process can then be broken down into:

• a global stage: on the whole interval [0, T ], a SKA solution is searched. The problem, written in its discretized residual form K u(t) = R eq (t), is solved after a projection on a suitable basis. The choice of this basis is mainly motivated by two points. Firstly, the residual forces must be precisely represented with only a short number of coefficients. This point is very important because the shorter the number of kept coefficients is, the shorter the computational time for each global stage is. However, the solution must be correctly represented and a minimum number of coefficients is necessary. Secondly, the basis function support should be, if it is possible, compact. This point enables the correction of the solution only on parts of the interval where the residue is important. • a local stage: from the solution obtained at the global stage, a PA solution is determined on the whole path. An incremental calculation in the time domain, done at all the instants t k , projects the SKA solution onto the elastic convex using the return mapping algorithm [START_REF] Nguyen | On the elastic plastic initial-boundary value problem and its numerical integration[END_REF][START_REF] Simo | Consistent tangent operators for rate-independent elastoplasticity[END_REF]. • a residue discounting: once the solution σ PA carried out, the residue is computed.

The algorithm initiation is classically determined by solving the purely elastic problem. The convergence criterion refers to the plastic admissibility condition of the SKA solution determined during the global stage.

Wavelet analysis of a signal

Wavelet analysis [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] is based on the decomposition of a signal s into an approximation sub-space and details subspaces at different scales. From an iterative point of view, from the sub-space V j , approximation V j +1 and details W j +1 are built. Mother wavelet ϕ, that has a null mean value, is used to define the details and, at the opposite, the scaling function ψ enables to calculate the approximations. In several cases both functions are reals and have compact supports. At the different scales, the functions ψ j,n (t) = 2 -j/2 ψ(2 -j tn) and ϕ j,n (t) = 2 -j/2 ϕ(2 -j tn) with j, n ∈ Z are used to calculate the wavelet and scaling coefficients:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ a n j = s, ϕ j,n = 1 √ 2 +∞ -∞ s(t) φ(2 -j t -n) d n j = s, ψ j,n = 1 √ 2 +∞ -∞ s(t) ψ(2 -j t -n) (1)
2 -j is the scaling factor of the wavelet, and n is the translation, and so localisation, parameter. Therefore a n j and d n j describe the behavior of s around the instant 2 -j n.

The signal can therefore be decomposed into an approximation at a coarse level j 0 , completed by details at all the intermediate ones:

s(t) = j -1 i=j 0 k d i k ψ ik (t) + k a j 0 k ϕ j 0 k (t) (2)
Scaling and wavelet coefficients, at the low levels, describe the function globally. On the other hand, at the fine levels, a local description is carried out. Thus, at those scales, wavelet coefficients are small when the signal is smooth, and the detail coefficients are significant only in discontinuous areas. An efficient approximation can be carried out keeping only the important ones.

For a fixed number N of coefficients, the approximation error will be minimized considering only the N basis functions ψ j,n and ϕ j,n for which the absolute values of coefficients a n j or d n j are the largest. Therefore, from the exact relation (2), the following approximation (3) is obtained.

s(t) j -1 i=j 0 k ∈[1,m i ] d i k ψ ik (t) + k ∈[1,l] a j 0 k ϕ j 0 k (t) with N = l + i m i < N instants (3)
In the following, without distinguishing wavelet and scaling functions, the notation will be

s(t) N i=1 α i q i (t) or N i=1 β i q i (t)

New resolution of the global stage

The residue decomposition R eq (t) on a wavelet basis enables to represent its evolution with only a few coefficients. Moreover, because the support length and the localisation of the basis functions are variable, the correction of the solution can be carried out only on one part of the studied interval [0, T ]. In that case, on the other part where the obtained solution has already converged, no basis function, and so no correction, is necessary.

From K(t) u(t) = R eq (t), decomposing u(t) and R eq (t) into their basis functions, the equation

N j =1 K(t)α j (M)q j (t) = N j =1 β j (M)q j (t)
is obtained and, multiplying by the test function q i (t) leads to

N j =1 K(t)α j (M)q i (t)q j (t) = N j =1
β j (M)q i (t)q j (t).

By integrating over the whole interval [0, T ] we obtain

N j =1 T 0 K(t)q i (t)q j (t) dt α j (M) = N j =1 T 0 q i (t)q j (t) dt β j (M)
that can be rewritten

K α(M) = β(M) (4)
where

K ij = Supp q i (t) ∩ Supp q j (t) K(t)q i (t)q j (t) dt, α = ⎡ ⎣ α 1 (M) . . . α N (M) ⎤ ⎦ and β = ⎡ ⎣ β 1 (M) . . . β N (M) ⎤ ⎦
In general, the size of this system is important because the size of the matrix [ K] is N × N dof over N × N dof . However, in most cases, [ K] presents an important number of elementary matrices [ K ij ] equal to zero. One reason is linked to the multi-scale approach in the wavelet decomposition. Because of the compact length of the wavelet basis functions, q i (t) and q j (t) supports can be disjoined. In this case, the intersection of the supports being empty, [ K ij ] equals zero. Otherwise the {q k (t)} k∈{1,N} basis is, in most cases, orthonormal. So, if K(t) is constant over Supp q i (t) ∩ Supp q j (t) then

[ K ij ] = Kδ ij . Thus, if i = j , the relation [ K ij ] = [0]
is still satisfied. Another important fact to reduce the calculation time is that [ K ij ] is a real symmetric matrix. Some interesting solution methods can therefore be applied to systems of considerable size to solve the global stage such as MINRES [START_REF] Greenbaum | Iterative Methods for Solving Linear Systems[END_REF].

When plastic strains are small or when the tangent modulus is quite important compared to the Young's modulus, a modified Newton method may be a good alternative to the method previously presented. The stiffness matrix is considered constant and equal to its initial value. In that particular case

[ K ij ] = Kδ ij and K = ⎡ ⎢ ⎢ ⎣ K K . . . 0 K ⎤ ⎥ ⎥ ⎦
Only the diagonal blocs remain and, for each wavelet coefficient, the system becomes

Kα i = β i ∀i ∈ {1; N } (5)

Numerical results

In order to illustrate the method, the case of a holed plate submitted to a thermomechanical loading is treated. The constitutive law is von Mises elastoplastic with a linear kinematic hardening depending on the temperature (Fig. 1(a)). Because of the problem symmetry, only a quarter of the plate has been meshed and adequate boundary conditions have been applied (Fig. 1

(b)).

A first uncoupled thermal computation enables the determination of the temperature field T (x, t) in the whole structure. In the mechanical calculation, carried out with Matlab, a traction is applied on the surface in x = L/2 in the y-axis direction and a thermal gradient is present in the thickness. In Fig. 2, similar levels of plastic strains on the contours of the plate at the instant t = T can be observed for the incremental method and the direct method. The calculation result, presented on the [0, T ] interval, at the most loaded point (the surrounded point on Fig. 1(b)), shows again a good agreement between the two methods. A superposition of the two curves after convergence is shown in Fig. 1(c). Moreover, compared to an incremental simulation, the calculation time is reduced.

Conclusion

The direct method presented in this paper enables the determination of the solution (u, ε, ε p and σ ) of structural evolution problems over the entire time interval of interest. It belongs to the class of large time increment methods and uses the wavelet transform. During the global stage, because of the wavelet adaptability, such a basis is used to expand the different fields. Hence the different signals can be well represented, even in complex cases, with a few coefficients. The example presented shows a good agreement between the incremental method and the direct method.

Fig. 1 .

 1 Fig. 1. Calculation of a holed plate submitted to a thermomechanical loading.

Fig. 1 .

 1 Fig. 1. Calcul d'une plaque trouée soumise à un chargement thermomécanique.

Fig. 2 .

 2 Fig. 2. Residual plastic strains with the direct (left) and the incremental method (right).

Fig. 2 .

 2 Fig. 2. Déformations plastiques résiduelles avec les méthodes directe (gauche) et incrémentale (droite).


K Dang Van

Z Moumni

Evaluation of fatigue-ratcheting damage of a pressurised elbow undergoing damage seismic inputs

We present a simplified method to calculate the plastic ratchet of elbow-shaped pipes submitted to seismic loading and an internal pressure. This method is simplified in the sense that the value of the ratchet is obtained without the use of finite element method (FEM) calculations. Here we derive a formula and use it to evaluate the fatigue-ratcheting damage of an elbow. This approach is applicable to complex plastic response appropriately described by non-linear kinematics hardening, which is more realistic for stainless steel such as 316-L

Introduction

Industrial structures, designed for transportation of pressurised fluids often requires high safety levels against accidental seismic loading. For example, in the case of nuclear structures, which must fulfil a criterion against fracturing, the most vulnerable components are the curved pipes because they are more flexible than straight pipes. They suffer much more severe stresses under the same bending moment. It is well known that this is due to the ovaling of the curved pipe cross section during bending. To predict their fatigue life, under seismic loading, it is necessary to calculate the local mechanical parameters that can initiate the general collapse of a structure by low-cycle plastic fatigue.

It is very difficult to perform direct experimentation and to interpret the obtained results for many reasons.

First, it is not easy to predict the locus of crack initiation. The strain gages are rapidly de-bonded, and it becomes impossible to measure local strains. To avoid these difficulties numerical simulations are performed to evaluate the strain and stress histories within the vicinity of points of stress concentrations. These calculations present a number of difficulties since we are dealing with a shell problem presenting severe geometrical (large displacement) and material (plasticity) nonlinearities. [START_REF] Boussa | Finite pure bending of curved pipes[END_REF] presented an original treatment of the finite bending of curved pipes of arbitrary cross sections. Their model has been used to treat the previous problem [START_REF] Boussa | Flexion plane du tuyau coude ´et endommagement sous se ´isme[END_REF][START_REF] Boussa | Fatigue-Seismic Ratcheting Interactions in Pressurized Elbows[END_REF]). Nevertheless, only two simple constitutive laws have been used: perfect plasticity, and linear kinematics hardening law. The first one underestimates the ratcheting rate threshold, whereas the second is inadequate to describe this phenomenon. In addition, the calculations are rather CPU intensive, limiting our ability to perform the large number of simulations necessary to derive a design rule for a code of practice.

In this paper, we propose a simplified model that relates the strain histories at some critical points to structure characteristics and seismic loading sequences. Our method is applicable to more realistic material behaviour taking into the account non-linear kinematics hardening. It is known that this model is more realistic for 316-L stainless steel.

The paper is organised as follows. In the second section, we give an example of calculation of the finite pure bending of a curved pipes based on [START_REF] Boussa | Finite pure bending of curved pipes[END_REF] model. In this formulation, it is shown that the problem can be reduced to an axisymmetric problem augmented with one degree of freedom. The third section is devoted to the simplified method. First we derive a formula that gives the ratchet deformation in the case of perfect plasticity. It is then generalised for non-linear kinematic hardening laws. In the fourth section, we present the fatigue analysis based on the modified Masson -Coffin equation [START_REF] Coffin | The deformation and fracture of a ductile metal under superimposed cyclic and monotonic strain[END_REF] to take account of ratcheting.

Numerical finite elements calculation

We give in this sequel an example of calculation of the finite pure bending of curved pipes using Boussa-Dang Van model.

Notation

Throughout the paper, the symbols and ( ) denotes the standard tensor products of two vectors or two second-order tensors and the differentiation with respect to time.

2.2. D. Boussa-Dang Van model Boussa et al. (1996), proposed a new and general formulation for solving by finite elements method (FEM) the finite pure bending of a curved pipes. In the following, the main hypotheses of their model are first recalled. Using the virtual work principle, the authors' show that the problem, after discretisation by the FEM, is reduced to an axisymmetric problem augmented by a single degree of freedom.

Kinematics

In its initial configuration assumed to be stress free, the curved pipe of interest, shown in Fig. 1 with the main notations, is supposed to be a three-dimensional body occupying a sector V 0 of an axis symmetric body. Let _ 0 denote a typical cross section, _ 0 sol being the structure of the pipe itself and _ 0 int being the cavity formed by this structure. In addition, we have S 0 int and S 0 ext presenting the inner and outer skins, c 0 the bend angle of the curved pipe and R 0 the radius of curvature defined by: Fig. 1. General view and typical cross section of the curved pipe.

R 0 = & _ 0 r d_ 0 & _ 0 d_ 0 (1)
Let M be the position of a typical particle with cylindrical co-ordinates (R, [, Z) and a normalised natural basis (E R , E [ , E Z ):

M =O +RE R (U)+ ZE Z (2)
Under external loads, the body deforms and a particle with initial position M occupies the new position m(r, q, z) with a normalised natural basis (e r , e q , e z ):

m= O +re r (q)+ ze z (3)
The pure in plane bending of the curved pipe can then be conveniently described with the following relations:

Á Ã Í Ã Ä r = r(R, Z, t) q =(1+a(t))U z=z(R, Z, t) (4)
where a(t) is a scalar function of time. The kinematics given by Eq. (4) includes the usual assumptions for the in-plane bending problem. The assumption that cross sections remain plane is accounted for by the dependence of q on U solely, and the assumption that all the cross sections deform likewise in their plane is enforced by the independence of r and z from U. Also, the assumption of uniform stretching of parallels to the axis of symmetry is the consequence of the linearity between q and U. In addition, fixing a(t) to 0 in Eq. ( 4) yields an axisymmetric evolution. Accordingly, arbitrary values of a(t) adds only one degree of freedom to the axisymmetric problem.

External loading and boundary conditions

In the sequel the Cauchy stress tensor will be denoted by |. In its initial state, the pipe is supposed stress free. The applied loads are of two types: those (typical in three-dimensional continuum descriptions) given pointwise, and those (less common) prescribed on average. In the current problem, these are the bending moment M due to the seismic input and the normal force N due to the effect of the pipe plug. They are defined on the current configuration by:

Á Ã Í Ã Ä M= & _tsol GM |e q n e z = & _tsol (r-R t )| qq d_ t N = pA int = & _tsol | qq d_ t (5)
In Eq. ( 5) p represents the internal pressure and A int the area of _ int . The boundary conditions are the following:

The outer skin S ext is stress free:

T d = |n= 0 (6)
The inner skin S int is submitted to uniform pressure:

T d = |n= -pn (7)
At the end cross sections _ 0 sol and _ U sol the tangential component of the applied shear force vanish:

T r d = T z d = 0 (8)
Finally the end cross section _ U sol is submitted to a prescribed bending moment M and to a normal force N. They are given by Eq. ( 5). In the case of a curvature controlled loading, the moment is an unknown variable that needs to be calculated and the rotation prescribed by the condition:

q=(1+a(t))U (9)

Virtual work equation

In the case of quasi-static evolution, the principle of Virtual Work is:

P i + P e = 0 ( 10
)
where P i and P e are respectively the power of internal and external forces. Denoting the rate of deformation by d. and the virtual velocity by 7 leads to:

-& d |:d. dd + & S int 7 ˆ|n dS + & _ U sol 7 ˆ|e q d_ -& _ 0 sol 7 ˆ|e q d_ (11)
where d is the region of space occupied by the curved pipe and 7 ˆ= r ˆer +rUa ˆeq +z ˆez

Taking into account the boundary conditions Eqs. ( 5)-( 7) and Eq. ( 8), yields the following equation:

-& d |:d. dd -& S int p7 ˆn dS +(M + pRA int )Ua ˆ= 0 (12)
The constitutive law considered is the perfect plasticity given by the following equation:

| = C:(d-d p ) ( 13
)
where C is the standard constant elastic modulus tensor of the material, which is assumed to be homogeneous and isotropic, and d p the plastic deformation rate.

The total rate of deformation can be decomposed in two terms d = d _ +a; e q e q . The first one, d _ , corresponds to the axisymmetric evolution and the second, a; e q e q corresponds to the deformation of the parallels.

Differentiating Eq. ( 12) with respect to time and taking into account, the constitutive law (Eq. ( 13)) and the decomposition of the total deformation the rate, yields the following equation:

-& d d _ :C:d. _ dd -p; & S int n7 ˆ_ + & d d p :C:d. _ dd -a ˆ&V d _ :C:(e q e q ) dd -a; & d d. _ :C:(e q e q ) dd
-a; a ˆ&d (e q e q ):C:(e q e q ) dd +(M : +p; RA int )Ua ˆ+a ˆ&d d P :C:(e q e q ) dd= 0 (14)

Further information concerning the model outlined above and its finite element discretisation is found in [START_REF] Boussa | Finite pure bending of curved pipes[END_REF].

An example of calculation

As an example, the bending of a pressurised curved pipe with perfectly plastic behaviour is presented. The computations were performed with Castem 2000 code. The geometry, material parameters and loading amplitudes are given in Table 1.

Because of symmetry, the mesh was defined only for half of the section as presented on Fig. 2. The imposed rotation history is given in Fig. 3. The outcome of the finite elements calculation is the evolution of the moment versus the rotation, given in Fig. 4. One can see that the stabilised state is reached after a few cycles only. This is due to the imposed displacement condition. Fig. 5 presents the history of the hoop strain at the critical point of the structure (internal flank) versus the imposed rotation. We observe that plastic deformation increases with the number of loading cycles. This ratcheting effect is due to the internal pressure that prevents ovalisation of the section, so that the rate of deformation is not the same during the opening and the closure sequences of the bending. The FEM computation is time consuming, so that a simplified method, which permits a greater number of simulations for a sensitivity analysis, is necessary. Such simplified model is discussed in the next section.

Simplified method

The simplified method is first presented for the same problem considered above with an elastic perfectly plastic model. It is then generalised for non-linear kinematic hardening laws.

E6aluation of the hoop strain in the inner skin

The simplified method is base on the analysis of the (Fig. 5). We approximate the hoop strain in the inner skin (critical point of the structure) by the bilinear curve given in the (Fig. 6). Fig. 6 shows that our model is based on four parameters h, i, k and l that can be determined by a parametric method in the following way. Let a y and M y is, respectively, the limit rotation and moment that lead to the first plastification. [START_REF] Leconte | Chaudronnerie et tuyauteries industrielles, tome 1 Conception, Technique et Documentation[END_REF] estimates the limit moment by the formula:

M y = S y I r e C (16)
where, S y the elastic limit of the material, I the moment of inertia of the section of the pipe, r e the external radius and C a stress concentration coefficient. This last coefficient is given by [START_REF] Clarck | Bending of curved tubes[END_REF] as:

C = 1.8 h 2/3 , h= tR r m 2 (17)
Where r m the mean radius, and t the thickness of the elbow.

Considering a linear relation between a y and M y yields the following equation:

a y = k h 2/3 R 1.8r e E S y (18)
k is a flexibility coefficient given by [START_REF] Clarck | Bending of curved tubes[END_REF] as:

k = 1.73 1-6 2 h (19)
Eq. ( 18) and Eq. ( 19) lead to:

A y o = -A y c = 90× a y =86.5 R 1 -6 2 h 2/3 r e E S y (20)
The internal pressure induces initial stresses having the tendency to decrease the value of the limit rotation. These stresses are given by:

| qq = pr m t (21)
Considering the difference between the behaviour of the curved pipe in opening and closing (because of the internal pressure) a parametric calculation gives:

A y o = 112.4 R 1-6 2 h 2/3 r e E S y -0.8 pr m t A y c = -86.5 R 1-6 2 h 2/3 r e E S y -0.8 pr m t ( 22
)
Finally l is given by:

l =199 R 1-6 2 h 2/3 r e E S y -0.8 pr m t (23)
The three left parameters h, i, k are determined following the same procedure using a parametric method based on some finite elements calculation.

Parameter h

This parameter presents the slope of the elastic part of the curve. It depends only on the geometry of the pipe. The parametric calculation gives for a:

h= - 1 98.3 tr m R 2 1/3 (24)
For the curved pipe at hand one finds: h= -0.00226.

Parameters i, k

They correspond to the slopes of the plastic parts of the curve for closing (i) and opening (k). They depend on the geometry and on the elastic limit S y and the pressure p. Their evaluation is more difficult than h because it is necessary to calculate the deformation response of the structure for some cycles.

It is clear that in the absence of pressure, the behaviour of the curved pipe is identical in opening and closing. The parametric calculation gives in this case:

i= k= 3.5h (25)
Taking into account the internal pressure leads to different values for i and k. To evaluate them, we fixed the geometry and the elastic limit S y and we calculated their variations with respect of p. The results are represented in (Fig. 7):

The curve given in the figure (Fig. 7) can be approximated by the following formulas: if

if pr m S y t B 1 3 Á Ã Í Ã Ä i= 3.5+10.5 pr m S y t h k= 3.5h
1 3 B pr m S y t B 2 3 Á Ã Í Ã Ä i =7h y = 6 -7.5 pr m S y t h if 2 3 B pr m S y t B 1 Á Ã Í Ã Ä i = 11 -6 pr m S y t h k = 3 -3 pr m S y t h (26)
In (Eq. ( 26)), the adimensional number (pr m )/(S y t) measures the internal pressure effect.

Our method permits the calculation of the response of the structure (hoop strain) in the critical point (the inner skin of the flank). This response can be decompsed into two parts:

A cyclic elastic part:

cir =hA(t) ( 27
)
where cir is the hoop strain.

A cyclic plastic part given by the derivative with respect of time of the plastic ratchet r; :

r; = Á Ã Í Ã Ä 0 if elastic evolution (i-h)A : (t) if plastic evolution (closing) (k -h)A : (t) if plastic evolution (opening) (28)
To obtain the global response, one must write the algorithm corresponding to the previous evolution and solve the problem numerically. Our method gives very good results with respect to the ones obtained by the FEM method. Furthermore, it is less time consuming.

E6aluation of the plastic ratchet in the case of a seismic inputs

Let us begin by evaluating of the plastic ratchet in the case of a centred cyclic loading opening/ closing as shown in (Fig. 8). Following our simplified model the evolution of the ratchet takes place when the amplitude of the imposed rotation is \ 2l as shown in Fig. 8. For every cycle, the plastic deformation (hoop strain) cycle is given as follows:

Á Í Ä (i -h)(-2A max +2l) for closing (k-h)(2A max +2l) for opening (29)
Where A max denotes the maximum of the absolute value of A(t).

Considering Eq. ( 29) and taking into account the number of cycles N yields the following expression of the final ratchet:

r = 2Ni -k(A max -l) (30)
For the validation of our model, we calculated the ratchet due to a ten cycle loading with a maximum amplitude A max =m -0.554 (m = (maxa(t)/a y ) using both the simplified formula and the FEM. In both cases we used the geometry and material parameters given in Table 1. The pressure was taken: 11.7 MPa. The results are given in Table 2. One can observe a good correlation between the two methods.

For a realistic seismic input like the one given in (Fig. 9), we can use the same procedure. In addition, approximating the loading history as Gaussien process leads the following formula:

r= 2y| 0 A max i-k 1-erf 1 2| 0 A max (31)
where | 0 is the standard deviation of the signal with normalised maximum amplitude, and erf is the error function defined by:

erf(z)= 2 y & z 0 e -y 2 dy (32)
Table 3 summarises the results of a series of calculations of the ratchet by using either the simplified method (Eq. ( 31)) and the direct FEM. The predictions of the simplified method are very satisfactory.

Extension of the simplified method for more realistic material beha6iour

Generally, steels present a non-linear plastic behaviour, which is not well described by the perfect plasticity. Non linear kinematic hardening effects are observed. They are often coupled with isotropic hardening. Such complex memory effects are observed in 316-L steel. To account for these effects, internal hardening parameters are necessary [START_REF] Chaboche | Modelisation of the strain memory effect on the cycle hardening of 316-L stainless steel Proceedings[END_REF]. The evolution of these parameters is very sensitive to the sequences of the applied loading so that the previous version of the simplified method cannot be applied directly. An extension is then necessary. Our main idea is based on the use of an equivalent-elastic perfectly plastic model in the sense that the final ratchet and fatigue damage should be the same as the ones predicted by the FEM calculation according to non linear hardening law.

Calculation of the equi6alent elasticity limit

Suppose that we have the resulting calculations of the FEM using a non-linear kinematic hardening law. The outcome of this calculation is the internal hoop strain in the inner flank of the pipe. We then look for an equivalent perfect plastic law, which gives, for the same level of loading, similar values of the final ratchet and fatigue damage as those given by the non-linear kinematic hardening behaviour. This equivalent law is subsequently characterised by a fictitious elastic limit S* y . This limit is obtained by a parametric procedure based on a FEM calculations. The equivalent law must satisfy the following loading conditions:

Á Ã Í Ã Ä S y *\S y a y *\ a y m*Bm P*BP (33)
In Eq. ( 33) the symbol (*) denotes fictitious quantities. We recall that, p represents the internal pressure, p y the pressure level that gives the first plasticity without bending moment, S y the elastic limit, a y the angle corresponding to S y , P and m are adimensional constants that characterise the level of loading such that:

Á Í Ä p= P× p y a= m ×a y
The level of loading is determined in such a way that one obtains the same maximal opening amplitude of the elbow as with the FEM calculations. This parameter m* is then given by: m*= m× a y a y * (34)

Results of the parametric analysis

The elastic limit S y * is calculated using an iterative process: we look for the S y * that permits to obtain similar values of the ratchets given by the FEM and the simplified method by minimising their difference. Examples of results from the parametric calculations are presented in Table 4. Note that for the FEM calculations we submitted the pipe to a seismic input.

Using results from

Fatigue analysis

The last step in structural integrity assessment is the fatigue limit analysis. Let N f0 denote the number of cycles to failure under cyclic strain only and N f the number of cyclic to failure under the combination of cyclic strain and ratchet. The modified Masson-Coffin equation to take account of ratcheting is then: fatigue-ratcheting analysis is proposed. It is applicable to materials with complex behaviour and does not necessitate lengthy FEM calculations. Our model is suitable to perform a great number of simulations necessary for a parametric sensitivity analysis that could be required to derive a design rule for a code of practice.

N f N f 0 = 1-r max f -1/c (36)
where r max is the value of the ratchet given by the previous simplified calculation, f is the material ductility and c is a material constant. The cumulative damage is given by

D = S cycles 1 N f ( 37
)
where the sum is performed on the number of cycles of the seismic input. Table 5 gives some predictions of the induced damage either by using non-linear kinematic hardening law and the simplified model. These results are satisfactory.

Conclusion

In this paper, a complete simplified method for .

Fig. 2 .

 2 Fig. 2. Mesh over the half of the section.

Fig. 3 .

 3 Fig. 3. Imposed rotation versus time.

Fig. 5 .

 5 Fig. 5. Hoop strain (%) versus imposed rotation.

  The elastic range is given by the interval of rotation [2A y o , 2A y c ] where A y o and A y c are the elastic limits of the elbow during opening and closing, respectively. The difference between A y o and A y c is denoted l. The hoop strain varies linearly with the imposed rotation angle A= 90× a(t). The slope of the curve is denoted by h. The plastic range is characterised by two different curves with different slopes: i for opening and k for closing.

Fig. 6 .

 6 Fig. 6. Parameters of the simplified model. Fig. 4. Moment (Nm × 10 4 ) versus imposed rotation.

Fig. 7 .

 7 Fig. 7. Evolution of i and k with respect of the internal pressure.

Fig. 8 .

 8 Fig. 8. Plastic evolution takes place for the bold parts of the curve.

Fig. 9 .

 9 Fig. 9. Example of seismic input.

Table 2

 2 Comparaison of the simplified method and FEM results

				Table 3
				Comparaison of the simplified method and FEM results for a
				seismic input
	m	Simplified ratchet (%)	FEM ratchet (%)
				m	Simplified ratchet (%)	FEM ratchet (%)
	1	0	0
	2	3	2.9	1.85	0.56	0.61
	3	13.5	12.9	2.59	3.8	3.73
	4	24.5	21.7	3	7.09	7.04
				3.5	12.28	12.23

Table 4

 4

	one can determine a

Table 4

 4 Results of the parametric calculation

	p (MPa)

Acknowledgements

The authors would like to acknowledge Electricite ´De France (EDF) for their financial support. Thanks are due to V. Louche who carried out some Finite element calculations and to D.J. Boussa for fruitful discussions.