Jean-Claude Michel
email: michel@lma.cnrs-mrs.fr

Pierre Suquet

NONUNIFORM TRANSFORMATION FIELD ANALYSIS: A REDUCED MODEL FOR MULTISCALE NONLINEAR PROBLEMS IN SOLID MECHANICS

published or not. The documents may come     Nonuniform transformation field analysis: a reduced model for multiscale nonlinear problems in solid mechanics

Introduction

A common engineering practice in the analysis of composite structures is to use effective or homogenized material properties instead of taking into account all details of the individual phase properties and geometrical arrangement (fiber and matrix in the case of a fiber-reinforced composite). These effective properties are sometimes difficult to measure and this difficulty has motivated the development of mathematical homogenization techniques which provide a rational way of deriving effective material properties directly from those of the individual constituents and from their arrangement or microstructure. A further interest of such predictive schemes is that material or geometrical parameters can be varied easily which opens the way for tayloring of new materials for a given application. Although homogenization has been developed for both periodic (Sanchez-Palencia [START_REF] Sanchez-Palencia | Nonhomogeneous media and vibration theory[END_REF] ) or random composites (Milton 20 ), the present study is focused on periodic composites.

Periodic homogenization of linear properties of composites is now wellestablished and the reader is referred to Bensoussan et al [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF] or Sanchez-Palencia [START_REF] Sanchez-Palencia | Nonhomogeneous media and vibration theory[END_REF] for the general theory, and to Suquet [START_REF] Suquet | Elements of Homogenization for Inelastic Solid Mechanics[END_REF] or Guedes and Kikuchi 10 (among others) for computational aspects. The central theoretical result for linear properties is that, provided that the scales are well separated, the linear effective properties of a composite are completely determined by solving a finite number of unit-cell problems. These unit-cell problems are solved once for all and their resolution yields the effective properties of the composite. Then the analysis of a structure comprised of such a composite material can be performed using these effective linear properties. In summary, for linear problems, the analysis consists of two completely independent steps, an homogenization step at the unit-cell level only, and a standard structural analysis performed at the structure level only.

In comparison, the situation for nonlinear composites is more complicated. For composites governed by a single nonquadratic but strictly convex potential (elastic potential or dissipation potential) homogenization results can be established to define an effective behavior, deriving from an effective potential (provided that the scales are well separated). However, except in very specific cases, this effective potential cannot be found by solving a small, or even a finite, number of unit-cell problems. To each macroscopic stress or strain state corresponds a unit-cell problem which has to be solved independently of the unit-cell problem for a different macroscopic state. Therefore, although there exists a homogenized behavior for the composite, the rigorous analysis of a composite structure consists of two coupled computational problems: 1) a structural problem where the (unknown) effective constitutive relations express the relations between the microscopic stress and strain fields solution of the second problem, 2) a unit-cell problem whose loading conditions are imposed by the (unknown) macroscopic stress or strain (or their rates).

Exactly the same type of complication occurs when the composite is made of individual constituents governed by two potentials, free-energy and dissipation potential, accounting for reversible and irreversible processes respectively. The most common examples of such materials are elasto-viscoplastic or elastoplastic materials. It has long been recognized by Rice, [START_REF] Rice | On the structure of stress-strain relations for time-dependent plastic deformation in metals[END_REF] Mandel [START_REF] Mandel | Plasticité classique et Viscoplasticité[END_REF] or Suquet, [START_REF] Suquet | Local and global aspects in the mathematical theory of plasticity[END_REF] that the exact description of the effective constitutive relations of such composites requires the determination of all microscopic plastic strains at the unit-cell level. For structural computations, the consequence of this theoretical result is that the number of integration points required in the analysis is equal to the product of the number of integration points at all scales, which is prohibitively large. With the increase in computational power, numerical strategies for solving these coupled problems have been proposed (see Feyel and Chaboche [START_REF] Feyel | FE2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials[END_REF] or Terada and Kikuchi 32 for instance) but are so far limited by the formidable size of the corresponding problems.

In order to derive constitutive models of the effective behavior of composites which are both useable and reasonably accurate, one has to resort to approximations. The Transformation Field Analysis (TFA) originally proposed by Dvorak and Benveniste [START_REF] Dvorak | On transformation strains and uniform fields in multiphase elastic media[END_REF] is an elegant way of reducing the number of macroscopic internal variables by assuming the microscopic fields of internal variables to be piecewise uniform. It has been extended by Fish et al [START_REF] Fish | Computational plasticity for composite structures based on mathematical homogenization: Theory and practice[END_REF] to periodic composites using asymptotic expansions. Assuming the eigenstrains to be uniform within each individual constituent, Fish et al [START_REF] Fish | Computational plasticity for composite structures based on mathematical homogenization: Theory and practice[END_REF] derived an approximate scheme which they called, for a two-phase material, the "two-point homogenization scheme". The original scheme and this extended scheme have been incorporated successfully in structural computations ( [START_REF] Dvorak | The modeling of inelastic composite materials with the transformation field analysis[END_REF][START_REF] Fish | Computational Mechanics of Fatigue and Life Predictions for Composite Materials and Structures[END_REF][START_REF] Kattan | Overall damage and elastoplastic deformation in fibrous metal matrix composites[END_REF] ). However, it has been noticed ( [START_REF] Chaboche | Towards a micromechanics based inelastic and damage modeling of composites[END_REF][START_REF] Michel | Constitutive relations involving internal variables based on a micromechanical analysis[END_REF][START_REF] Suquet | Effective properties of nonlinear composites[END_REF] ) that the application of the TFA to two-phase systems may require, in certain circumstances, a subdivision of each individual phase into several (and sometimes numerous) sub-domains to obtain a satisfactory description of the effective behavior of the composite. The need for a finer subdivision of the phases stems from the intrinsic nonuniformity of the plastic strain field which can be highly heterogeneous even within a single material phase. As a consequence, the number of internal variables needed to achieve a reasonable accuracy in the effective constitutive relations, although finite, is prohibitively high.

In order to reproduce accurately the actual effective behavior of the composite, it is important to capture correctly the heterogeneity of the plastic strain field. This observation has motivated the introduction in [START_REF] Michel | Constitutive relations involving internal variables based on a micromechanical analysis[END_REF][START_REF] Michel | Nonuniform Transformation Field Analysis[END_REF] of nonuniform transformation fields. More specifically the (visco)plastic strain within each phase is decom-posed on a finite set of plastic modes which can present large deviations from uniformity. An approximate effective model for the composite can be derived from this decomposition where the internal variables are the components of the (visco)plastic strain field on the (visco)plastic modes. This theory is called the Nonuniform Transformation Field Analysis (NTFA). For two-phase composites (nonlinear matrix and elastic fibers), comparison of the classical TFA, and of the NTFA with numerical simulations of the response of a unit-cell under monotone or cyclic loadings, has shown the accuracy of the NTFA [START_REF] Michel | Nonuniform Transformation Field Analysis[END_REF] . The present study is devoted to presentation of the NTFA and to its implementation into a macroscopic structural Finite Element analysis. It will be shown that the NTFA not only provides accurate predictions for the effective behavior of composite materials, which is its initial goal, but also provides an accurate approximation of the local fields which are the quantities of interest in predicting the lifetime of structures.

Structural problems with multiple scales

Homogenization and two-scale expansions

Structures made of composite materials naturally involve two very different length-scales. The largest scale (the macroscopic scale) is related to the structure itself and is characterized by a length L (Figure 1.1). The second and smallest scale (the microscopic scale) is related to the size of the heterogeneities in the composite material (typically the fiber scale in fiber-reinforced structures). The typical length at this scale is denoted by d. In fiber-reinforced laminates, d is of the order of the fiber diameter, whereas L is typically related to the thickness, or length, of the layered structure. When the scales are "well-separated", i.e. when the ratio η = d/L is small (η ≪ 1), one can expect all details about the microstructure to be "smeared out". In other words, the response of the structure at the macro-scale can be computed by replacing the very contrasted physical properties of the individual constituents by effective or homogenized properties (at the macro-scale).

The aim of the mathematical theories of homogenization is to determine exactly or to bound these effective properties from the information available, often partially, on the individual constituents themselves and on their arrangement (microstructure). However, if effective properties are sufficient for analysis performed in the linear range (stiffness of a composite structure, few first eigenfrequencies...) where the structure responds macroscopically as a whole, in many problems of engineering interest it is essential to take into consideration not only averaged fields, or effective properties, but also full local fields. Damage or fracture for instance are dramatically dependent on the local details of the strain or stress fields. The procedure by which the local fluctuations of fields are reconstructed from their macroscopic average is sometimes called localization and one important objective of the present approach is to propose a simple localization rule for strain and stress fields. The microstructure of periodic composites is completely known as soon as the geometry of a single unit-cell V is prescribed. For such composites, homogenization results can be obtained heuristically by means of two-scale expansions making use of the fact that the parameter η = d/L is small and that the geometry (and therefore the fields) are periodic at the microscopic scale (Sanchez-Palencia 24 , Bensoussan et al 2 ). Rigorous mathematical techniques have been developed to establish convergence theorems which usually confirm that homogenization results obtained by asymptotic expansions usually hold true (see for instance Tartar 31 ).

A brief reminder (by no means exhaustive) about two-scale expansions is given now. A function f defined on the macroscopic structure has variations at the two different spatial scales and can be denoted as f (X, x) to highlight this dependence on both variables, where X denotes the macroscopic spatial variable (structural scale) whereas x denotes the microscopic variable (at the unit-cell level). A dependence of a function on the microscopic variable x corresponds to fast oscillations of this function at the macroscopic scale, whereas a dependence on the macroscopic variable X corresponds to slower variations at the structural level.

The scale ratio η is finite and different from 0 in the actual structure (even though it is convenient mathematically to consider that it tends to 0). Therefore all mechanical fields (stress, strain, displacement...) in the actual structure depend on this ratio. For instance the displacement field and the stress field in the actual structure will be denoted by u η and σ η . The homogenized relations are obtained by taking the limit of u η and σ η as η goes to 0 and by studying the set of equations satisfied by these limits. These limits can be determined by means of two-scale expansions. For any function f η defined on the composite structure with finite scale ratio η, its two-scale expansion is defined as:

f η (X) = +∞ j=0 η j f j X, X η , (1.1)
where, by virtue of the periodicity of the microstructure, all functions f k (X, x)'s are periodic with respect to the microscopic variable x. Therefore, for a macroscopic point X, the argument X/η of the functions f j , denotes the location of X in the unit-cell at the microscopic scale. Let us recall that, setting g η (X) = g X η where g is periodic over the unitcell, the limit of g η as η goes to 0 is the average of g over the unit-cell. The convergence is weak (only in average) and not pointwise. Consequently the limit of f η as η goes to 0 is the average with respect to x of the zero-th order term in the expansion (1.1):

lim η→0 f η (X) = f 0 (X) = 1 |V | V f 0 (X, x) dx.
The homogenized (or effective) relations for the composite are therefore the relations between the limits as η goes to 0 of the fields σ η and ε η , or equivalently between the averages of the zero-th order terms in the expansion of the stress field and strain field (or strain-rate field), and additional internal variables α, depending on the constitutive relations of the individual constituents which remain to be specified (see section 1.2.2).

To understand how these zero-th order terms behave, one has to expand the unknown displacement, strain and stress fields u η , ε η and σ η in powers of η, after due account of the equations satisfied by these fields. In addition to the constitutive equations (to be specified), these equations are the compatibility equations and the equilibrium equations:

ε η ij = 1 2 du η i dX j + du η j dX i , dσ η ij dX j + F i = 0, (1.2)
where F denote the body forces applied to the structure. The derivation of a twoscale function f (X, x) which is periodic with respect to x with x = X/η is performed according to the chain-rule:

d dX = ∂ ∂X + 1 η ∂ ∂x .
Applying this derivation rule to the double-scale expansion of u η , ε η and σ η :

u η (X) = u X, X η = ∞ k=0 η k u k (X, x), ε η (X) = ε X, X η = ∞ k=0 η k ε k (X, x), σ η (X) = σ X, X η = ∞ k=0 η k σ k (X, x),                    (1.3)
one obtains the expansion of the compatibility and equilibrium equations in powers of η:

Order -1 : ε x (u 0 ) = 0, div x σ 0 = 0, Order 0 : ε 0 = ε X (u 0 ) + ε x (u 1 ), div X σ 0 + div x σ 1 + F = 0,
σ 0 , ε 0 and α 0 satisfy the constitutive relations.

                     (1.4)
Similar equations corresponding to higher-order terms in the expansions can be obtained in the same way. The operators ε x and div x in (1.4) stand for the deformation and divergence operators with respect to the microscopic variable x (with similar conventions for these operators with respect to the macroscopic variable X). The constitutive equations of the phases may involve internal variables, in which case the zero-th order terms of these internal variables also enter the relations between σ 0 and ε 0 . It follows from the first equation of the first line in (1.4) that u 0 (X, x) = u 0 (X). u 0 has no dependence on the microscopic variable (no fast oscillations in the displacement field). In addition, taking the average over the unit-cell of the first equation at order 0 (second line in (1.4)), and taking into account the fact that the average of the gradient of a periodic function vanishes identically, one obtains that:

ε X (u 0 )(X) = ε 0 (X), (1.5)
where an overlined letter denotes an averaged quantity:

ε 0 (X) = ε 0 (X, .) with f = 1 |V | V f (x) dx.
In other words the macroscopic strain ε X (u 0 ) is the average over the unit-cell of the zero-th order term in the expansion of the strain field ε η .

Unlike the displacement field, the zero-th order terms σ 0 and ε 0 of the stress and strain fields have microscopic fluctuations (i.e. they depend on both the macroscopic and the microscopic variables). It follows from the second equation in the first line of (1.4) that σ 0 is self-equilibrated at the microscopic scale, whichever body forces F are applied to the structure at the macroscopic scale. Taking the average over the unit-cell of the third line in (1.4), and noting that the average of the divergence of a periodic field vanishes identically, one finds that the average stress σ 0 = σ 0 satisfies the macroscopic equilibrium equations: div X σ 0 + F = 0.

(1.6)

The two equations (1.5) and (1.6) are valid irrespective of the constitutive behavior of the phases. The homogenized, or effective, constitutive relations relate the average stress σ 0 and the average strain ε 0 . The determination of these relations requires, in principle, a complete knowledge of the fields σ 0 and ε 0 with all their microscopic fluctuations. The dependence of these fields on the macroscopic variable X is known by solving the equilibrium problem for the structure subjected to the imposed macroscopic loading and where the effective constitutive relations are used for the composite material. Their dependence on the microscopic variable is known by solving the so-called local problem (or unit-cell problem), where the macroscopic variable X is only a parameter and will be omitted for clarity:

ε 0 (x) = ε 0 + ε x (u 1 (x)
) in V where u 1 is periodic, div x σ 0 = 0 in V , σ 0 .n anti-periodic on ∂V, σ 0 , ε 0 and α 0 are related by the constitutive equations of the phases.

  

(1.7)

The anti-periodicity condition for the traction σ 0 .n on ∂V derives from the periodicity of σ 0 and the anti-periodicity of n on opposite sides of the unit-cell V . The first line in (1.7) can be replaced by ε 0 = ε 0 and periodicity conditions. The constitutive relations of the phases have to be specified in order to further exploit these relations. For simplicity, the zero-th order terms ε 0 and σ 0 will simply be denoted by ε and σ in the rest of the paper and the dependence on the variable X will be omitted in the rest of this section.

Individual constituents

As already noted, the microstructure of periodic composites is completely specified by the knowledge of a unit-cell V , which plays, for periodic media, a role parallel to that of a representative volume element (r.v.e) in homogenization theories for random media. The unit-cell V is occupied by N homogeneous phases r) with characteristic function χ (r) (x) and volume fraction c (r) :

V (
χ (r) (x) = 1 if x ∈ V (r) , 0 otherwise, c (r) = χ (r) .
The average of a field f over the unit-cell V and over each individual phase V (r) are denoted by overlined letters f and f (r) :

f = f = N r=1 c (r) f (r) , f (r) = f r = 1 |V (r) | V (r) f (x) dx.
The composite structures of interest for this study may be subjected to thermomechanical loadings. Therefore the validity of the constitutive relations of the individual constituents must cover a wide range of temperature and strain-rates. For simplicity, attention will be restricted here to isotropic materials. We shall adopt in the sequel a viscoplastic model with nonlinear kinematic hardening proposed by Chaboche, [START_REF] Chaboche | Cyclic Viscoplastic Constitutive Equations, Part I: A thermodynamically consistent formulation[END_REF] generalizing the Armstrong-Fredericks constitutive relations:

σ = L : (ε -ε vp ), εvp = 3 2 ṗ s -X (σ -X) eq , ṗ = ε0 ((σ -X) eq -R) + σ 0 n , Ẋ = 2 3 H εvp -ξ X ṗ, R = R(p),            (1.8)
where (.) + denotes the Mc Cauley bracket (positive part):

A + = A if A ≥ 0, A + = 0 if A ≤ 0.
When the phases are isotropic, their elastic properties are characterized by a bulk modulus k and a shear modulus G. Kinematic hardening effects are characterized by the back-stress X whereas isotropic hardening manifests itself through the dependence of the yield stress R(p) on the cumulated viscoplastic strain p defined as ṗ = (2/3 εvp : εvp ) 1/2 . To simplify notations it is useful to introduce the viscoplastic potential:

ψ(A, R) = σ 0 ε0 n + 1 (A eq -R) + σ 0 n+1 , (1.9)
by means of which the second line of the constitutive relations (1.8) can be written as

εvp = ∂ψ ∂A (σ -X, R), ṗ = - ∂ψ ∂R (σ -X, R).
(1.10)

The model (1.8) (and subsequent refinements which will not be considered here) is commonly used in the analysis of the lifetime of metallic or polymeric structures under repeated thermomechanical loadings (see Samrout et al [START_REF] Samrout | Model for 28CrMoV5-8 steel undergoing thermomechanical cyclic loadings[END_REF] and Amiable et al 1 among others). The material parameters of the model, namely the elastic moduli L, the rate-sensitivity exponent n, the flow-stress σ 0 , the isotropic hardening function R(p), the kinematic hardening modulus H and the spring-back coefficient ξ, are strongly temperature-dependent. For simplicity, thermal loadings and thermal strains will not be considered in the present analysis, but the strong temperature-dependence of the material parameters will be accounted for. For instance the rate-sensitivity exponent n can vary from 5 to 20 for Aluminum alloys when the temperature varies from 20 0 C to 500 0 C. The method proposed here will make use of certain objects, called plastic modes, identified at a given temperature but used over the whole range of temperature with the appropriate material parameters. In other words, these plastic modes do not need to be identified at each temperature.

Unit-cell problem. Effective response of heterogeneous materials

As seen in section 1.2.1, the first order terms of the stress and strain field solve a unit-cell problem (also called the local problem) consisting of the equilibrium and compatibility equations (1.7) and the constitutive relations (1.10). All material properties are assumed to be uniform in each individual phases:

L(x) = N r=1 L (r) χ (r) (x), ψ(x, A, R) = N r=1 ψ (r) (A, R) χ (r) (x).
The overall stress σ and the overall strain ε are the averages of their local counterparts σ and ε (for simplicity the dependence on the macroscopic variable X of all fields will be omitted):

σ = σ , ε = ε . (1.11)
The homogenized effective relations are the relations between the macroscopic stress σ (and its time-derivatives) and the overall strain ε (and its timederivatives).

To find these relations, an history of macroscopic strain ε(t) is prescribed on a time interval [0, T ] generating a time-dependent local stress field σ(x, t). Its average σ(t) is the macroscopic stress whose history is therefore related to the history of ε(t).

The local problem to be solved to determine σ(t) reads:

σ(x, t) = L(x) : (ε(x, t) -ε vp (x, t)), εvp (x, t) = ∂ψ ∂A (σ(x, t) -X(x, t), R(x, t)), div x σ(x, t) = 0, ε(t) = ε(t), boundary conditions.                (1.12)
In view of the local periodicity of the structure, periodic boundary conditions are assumed on the boundary of the unit-cell.

The average of the local stress field σ(x, t) is the macroscopic stress response of the composite to a prescribed history of macroscopic strain ε(t). Unfortunately, except in very specific situations (e.g. laminates), these effective relations for nonlinear materials cannot be given in closed form. They are accessible only numerically, along a prescribed path. An important consequence of this observation for the computational analysis of a composite structure, is that the macroscopic and microscopic levels are intimately coupled. At the structural level, the macroscopic strain ε(X, t) is a function of position and a problem similar to (1.12) has to be solved at every macroscopic point X or, in a computational analysis, at every macroscopic integration point. As pointed out by Fish and Shek 8 , history data has to be updated at a number of integration points equal to the product of the numbers of integration points at all scales at each time increment.

To avoid the computational difficulty associated with the coupling of scales, approximations are introduced to render the resolution of the local problem (1.12) possible in closed form or amenable to simple algebra.

An auxiliary elasticity problem

Before introducing approximate resolution schemes for the local problem (1.12), it is important to emphasize that the stress and strain fields are solution of a linear elasticity problem on the unit-cell when the fields of internal variables are known. Indeed, assuming that the viscoplastic part of the strain is prescribed, the stress and strain fields in the r.v.e. solve the following linear elastic problem, with appropriate boundary conditions (for simplicity the time dependence of the fields is omitted):

σ(x) = L(x) : (ε(x) -ε vp (x)), divσ(x) = 0, ε = ε. (1.13)
Assume that ε vp (x) is known. It plays the role of a thermal strain in thermoelasticity when the temperature is prescribed, or that of a transformation strain in phase transformation problems.

The solution of (1.13) can be expressed in terms of ε vp and ε by a straightforward application of the superposition principle. Consider first the case where ε vp is identically 0. Problem (1.13) is then a standard elasticity problem and its solution can be expressed by means of the elastic strain-localization tensor A(x) as:

ε(x) = A(x) : ε. (1.14)
Consider next the case where ε = 0 and ε vp (x) is arbitrary. Problem (1.13) can then be written as an elasticity problem with eigenstress (sometimes called polarization stress)

τ (x) = -L(x) : ε vp (x) σ(x) = L(x) : ε(x) + τ (x), divσ(x) = 0, ε = 0. (1.15)
Introducing the nonlocal elastic Green operator Γ(x, x ′ ) of the nonhomogeneous elastic medium, the solution of (1.15) can be expressed as:

ε(x) = -Γ * τ (x) where Γ * τ (x) def = 1 |V | V Γ(x, x ′ ) : τ (x ′ ) dx ′ .
(1.16) The superposition principle applied to (1.14) and (1.16) gives that the solution of (1.13) reads as:

ε(x) = A(x) : ε + 1 |V | V D(x, x ′ ) : ε vp (x ′ ) dx ′ = A(x) : ε + D * ε vp (x),
(1.17) where the nonlocal operator D(x, x ′ ) = Γ(x, x ′ ) : L(x ′ ) gives the strain at point x created by a transformation strain at point x ′ .

Nonuniform transformation field analysis (NTFA)

Motivation: approximate resolution of the local problem

The Transformation Field Analysis (TFA), originally developed by Dvorak [START_REF] Dvorak | The modeling of inelastic composite materials with the transformation field analysis[END_REF] (see also references herein), is based on the assumption that the viscoplastic strains are uniform within each individual domain V (r) :

ε vp (x, t) = N r=1 ε (r) vp (t) χ (r) (x).
(1.18)

The determination of the field ε vp (x) is therefore reduced to the determination of the tensorial variables ε (r)

vp , r = 1, ...N . Using this decomposition, the macro-scopic stress reads as:

σ = N r=1 c (r) σ (r) , σ (r) = σ r = L (r) : (ε (r) -ε (r) vp ), (1.19)
where

ε (r) = ε r = A (r) : ε + N s=1 D (rs) : ε (s) vp , r = 1, ..., N, (1.20)
and

A (r) = A r , (1.21)
D (rs) = 1 c (r) 1 |V | 1 |V | V V χ (r) (x)Γ(x, x ′ ) : L(x ′ )χ (s) (x ′ ) dx ′ dx. (1.22)
The evolution of ε

vp is governed by the constitutive relations of the individual phases applied to the average stresses and thermodynamic forces on the phases. Assuming that these constitutive relations take the form (1.8) (or (1.10)), with material properties labelled by the phase r, the evolution equations for the generalized variables ε (r) vp read as:

ε(r) vp = ∂ψ (r) ∂A (σ (r) -X (r) , R (r) ), ṗ(r) = - ∂ψ (r) ∂R (σ (r) -X (r) , R (r) ), Ẋ(r) = 2 3 H (r) ε(r) vp -ξ (r) X (r) ṗ(r) , R (r) = R (r) (p (r) ).          (1.23) When a prescribed path ε(t), t ∈ [0, T ]
is prescribed in the space of macroscopic strains, the corresponding history of the average strains ε (r) (t) and viscoplastic strains ε A nice feature of the TFA is that its implementation is relatively easy. However applying the TFA to two-phase systems using plastic strains which are uniform in each phase yields predictions of the overall behavior of the composite which can be unreasonably stiff (Suquet [START_REF] Suquet | Effective properties of nonlinear composites[END_REF] , Chaboche et al [START_REF] Chaboche | Towards a micromechanics based inelastic and damage modeling of composites[END_REF] ). The origin of this excessive stiffness is to be seeked in the intrinsic nonuniformity (in space) of the actual plastic strain field which can be highly heterogeneous even within a single material phase, a feature which is disregarded by the TFA. Dvorak et al [START_REF] Dvorak | The modeling of inelastic composite materials with the transformation field analysis[END_REF] have obtained better results by subdividing each phase into several subdomains. Unfortunately, although the refinement does improve the predictions, a rather fine subdivision of the phases is often necessary to achieve a satisfactory agreement (see Michel, Galvanetto and Suquet [START_REF] Michel | Constitutive relations involving internal variables based on a micromechanical analysis[END_REF] ), resulting in a prohibitive increase of the number of internal variables entering the effective constitutive relations. These observations have motivated the development of alternative approximate schemes (Michel and Suquet [START_REF] Michel | Nonuniform Transformation Field Analysis[END_REF] ).

Nonuniform transformation fields

The aim of the NTFA is to account for the nonuniformity of the plastic strain field. The field of anelastic strains is decomposed on a set of fields, called plastic modes, µ (k) :

ε vp (x, t) = M k=1 ε (k) vp (t) µ (k) (x).
(1.24)

Unlike in the classical Transformation Field Analysis, the modes µ (k) are nonuniform (not even piecewise uniform) and depend on the spatial variable x. The idea is that their spatial variations capture the salient features of the plastic flow in the unit-cell. They can be determined either analytically or numerically. Their total number, M , can be different (larger or smaller) from the number N of phases. The µ (k) 's are tensorial fields whereas the corresponding variables ε

(k)
vp are scalar variables.

Further assumptions will be made to simplify the theory:

H1: The support of each mode is entirely contained in a single material phase. It follows from this assumption that one can attach to each mode a characteristic function χ (k) , elastic moduli L (k) and a dissipation potential ψ (k) which are those of the phase supporting this mode. M (r) will denote the number of modes with support in a given phase V (r) . H2: The modes are incompressible:

tr µ (k) = 0. (1.25)
This assumption stems from the fact that the µ (k) are meant to represent (visco)plastic strain fields. As a consequence of this assumption, the field ε vp given by the decomposition (1.24) is incompressible, expected, with no restriction on the components ε (k)

vp . H3: The modes are orthogonal:

µ (k) : µ (ℓ) = 0 when k = ℓ. (1.26)
This condition is obviously met when the modes have their support in different material phases but has to be imposed to the modes when their support are in the same material phase. H4: The modes are normalized:

µ (k) eq = 1.
(1.27)

Reduced variables and influence factors

Using the decomposition (1.24) into (1.17), one obtains that:

ε(x) = A(x) : ε + M ℓ=1 η (ℓ) (x) ε (ℓ) vp , (1.28)
where

η (ℓ) (x) = D * µ (ℓ) (x)
is the strain at point x due to the presence of an eigenstrain µ (ℓ) (x ′ ) at point x ′ , the average strain ε being zero. Upon multiplication of equation (1.28) by µ (k) and averaging over V , one obtains

e (k) = a (k) : ε + M ℓ=1 D (kℓ) N ε (ℓ) vp , (1.29)
where the reduced strains e (k) , the reduced localization tensors a (k) and the influence factors D (kℓ) N

(N stands for NTFA) are defined as

e (k) = µ (k) : ε , a (k) = µ (k) : A , D (kℓ) N = µ (k) : η (ℓ) . (1.30)
By analogy with the equation defining the reduced strain e (k) in (1.30), one can define:

e (k) vp = µ (k) : ε vp = µ (k) : µ (k) ε (k) vp (no summation over k). (1.31)
Reduced stresses can be associated by duality to the generalized viscoplastic strains ε

(k)
vp (the notations are chosen so as to highlight the analogy between the reduced stress τ (k) and the resolved shear stress on the k-th system in crystal plasticity):

τ (k) = µ (k) : σ , x (k) = µ (k) : X .
(1.32)

Constitutive relations for the reduced variables

It remains to specify the reduced constitutive relations relating the reduced strains and stresses. A first set of equations is obtained upon substitution of the stress-strain relation (1.12) 1 into the definition (1.32) 1 of the reduced stresses τ (k) :

τ (k) = µ (k) : L : (ε -ε vp ) .
Elastic isotropy of the phases and assumptions H1 and H2 for the modes µ (k) lead to:

τ (k) = 2G (k) (e (k) -e (k) vp ), (1.33)
where G (k) denotes the shear modulus of phase r containing the support of mode k.

The second set of equations concerns the evolution of the generalized variables e (k) vp and x (k) . Using the definition (1.31) of e (k) vp and equations (1.9) and (1.10) for the evolution of the viscoplastic strain field ε vp (x), one obtains that:

ė(k) vp = µ (k) : εvp = 3 2 ṗ µ (k) : A A eq , A = σ -X, ṗ = - ∂ψ ∂R (A eq , R)
.

(1.34) At this stage an additional approximation must be introduced to derive a relation between the ė(k) vp 's, the τ (ℓ) 's and x (ℓ) 's. Different approximations are discussed in Michel and Suquet [START_REF] Michel | Nonuniform Transformation Field Analysis[END_REF] (uncoupled and coupled models) to which the reader is referred for further details. It follows from this work that the most accurate model is the so-called coupled model where the force acting on a mode is the quadratic average of all the generalized forces acting on all modes contained in the same phase. For a given phase r, the generalized force A (r) is defined as

A (r) =   M (r) k=1 τ (k) -x (k) 2   1/2 .
(1.35)

In this relation M (r) denotes the number of modes having their support in phase r. Then the relation (1.34) is modified by replacing A eq by A (r) and R by R (r) :

ė(k) vp = 3 2 ṗ(r) τ (k) -x (k) A (r) , ṗ(r) = - ∂ψ (r) ∂R A (r) , R (r) , R (r) = R (r) (p (r) ),
(1.36) where, again, r is the phase containing the support of µ (k) . The plastic multiplier ṗ(r) is the same for all modes having their support in the same phase r.

Finally, in order to obtain an evolution equation for the x (k) 's the last equation in (1.8) is multiplied by µ (k) and averaged over V :

ẋ(k) = µ (k) : Ẋ = 2 3 H (k) ė(k) vp -ṗ ξ µ (k) : X . (1.37)
Then, replacing as previously A eq by A (r) and R by R (r) in the expression of the plastic multiplier ṗ, one obtains:

ẋ(k) = 2 3 H (k) ė(k) vp -ṗ(r) ξ (k) x (k) . (1.38)
In summary the constitutive relations for the model are:

e (k) = a (k) : ε + M ℓ=1 D (kℓ) N ε (ℓ) vp , τ (k) = 2G (k) (e (k) -e (k)
vp ),

A (r) =   M (r) k=1 τ (k) -x (k) 2   1/2 , R (r) = R (r) (p (r) ), ė(k) vp = 3 2 ṗ(r) τ (k) -x (k) A (r) , ṗ(r) = ∂ψ (r) ∂A eq A (r) , R (r) , ẋ(k) = 2 3 H (k) ė(k) vp -ṗ(r) ξ (k) x (k) .                                                    (1.39)
The systems of differential equations (1.39) is to be solved at each integration point of the structure (macroscopic level). At each time increment, knowing the increment in macroscopic strain, the resolution of the system yields the e

(k)
vp 's from which the ε (k) vp 's can be obtained by inversion of (1.31). Once the internal variables ε (k) vp are determined, the local stress field in the composite resulting from (1.13) and (1.28) reads as:

σ(x, t) = L(x) : A(x) : ε(t) + M k=1 ρ (k) (x)ε (k) vp (t), where ρ (k) (x) = L(x) : η (k) (x) -µ (k) (x) .          (1.40)
The effective constitutive relations for the composite are obtained by averaging this stress field:

σ(t) = L : ε(t) + M k=1 ρ (k) ε (k) vp (t). (1.41)
The localization tensors a (k) , the influence factors

D (kℓ)
N , the effective stiffness L and the tensors ρ (k) are computed once for all.

Choice of the plastic modes

The plastic modes are essential for the accuracy of the method. However there is no universal choice for these modes and they should rather be chosen according to the type of loading which the structure is likely to be subjected to. This implies that the user has an a priori idea of the triaxiality of the macroscopic stress field, as well as of its intensity and its time history. For instance when the structure schematically depicted in Figure 1.1 is subjected to pure bending, the macroscopic stress is expected to have a strong uniaxial component. Therefore the plastic modes should incorporate information about the response of the unit-cell under uniaxial tension (and compression if the response is not symmetric in tension/compression). But close to points where the plate is supported, the macroscopic stress will likely exhibit a non negligible amount of transverse shear and transverse normal stress so that plastic modes accounting for the unit-cell response under transverse shear and transverse tension-compression should also be present in the set of modes. Similarly if one is interested in the response of the structure under monotone loading with limited amplitude, the information about the response of the unit-cell will be limited to certain monotone loading paths in stress space up to a limited amount of deformation.

Given the complexity of the microstructures under consideration, the plastic modes are not determined analytically but numerically from actual viscoplastic strain fields in the unit-cell. Different unit-cell responses along the different loading paths of macroscopic stresses stemming from the above qualitative analysis are determined numerically. Second, the plastic modes are extracted from the microscopic viscoplastic strain fields at a given macroscopic strain, which depends on the range of macroscopic strains which is expected in the structural computation. Different or additional loadings can be considered, depending on the problem and keeping in mind that it is desirable to approach as closely as possible the macroscopic loading paths expected at the different integration points of the composite structure.

One of the building assumption of the NTFA is the mode orthogonality (hypothesis H3). If this prerequisite is obviously met when the modes have their support in different material phases, it has to be imposed to the modes which have their support in the same material phase. Let θ (k) (x), k = 1, ..., M T (r) be potential candidates to be plastic modes in phase r. The procedure used to obtain these fields will be detailed in due time but they will not satisfy assumption H3 in general. The Karhunen-Loève decomposition (also known as the proper orthogonal decomposition or as the principal component analysis) is used to construct a set of (visco)plastic modes µ (k) (x), k = 1, ..., M T (r) from these fields θ (k) (x):

µ (k) (x) = MT (r) ℓ=1 v (k) ℓ θ (ℓ) (x), (1.42)
where v (k) and λ (k) are the eigenvectors and eigenvalues of the correlation matrix g:

MT (r) j=1 g ij v (k) j = λ (k) v (k) i , g ij = θ (i) : θ (j) . (1.43)
It is straightforward to check that the resulting modes are orthogonal (as any set of eigenvectors of symmetric matrices):

µ (k) : µ (ℓ) = λ (k) if k = ℓ, otherwise 0. (1.44)
Another advantage of the Karhunen-Loève decomposition is that the NTFA model is almost insensitive to modes with small intensity, or in other words to modes µ (k) corresponding to small eigenvalues λ (k) . Therefore, in practice, among the total M T (r) modes, it is sufficient to consider in the model the first M (r) modes corresponding to the largest eigenvalues (see Roussette et al [START_REF] Roussette | Nonuniform transformation field analysis of elastic-viscoplastic composites[END_REF] for more details).

Reduced localization tensors and influence factors.

Once the plastic modes are chosen, the localization and influence tensors can be determined by solving only linear problems. The strain localization tensor A is obtained by solving successively 6 linear elasticity problems a :

σ(x) = L(x) : ε(u(x)), div (σ(x)) = 0, ε = ε, (1.45)
where ε is taken to be equal successively to one of the second-order tensors i (ij) with components

i (ij) mn = 1 2 (δ im δ jn + δ in δ jm ) .
a 6 problems in dimension 3, but only 3 problems in plane strain problems and 4 problems in generalized plane strain problems, see Michel et al [START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF] for further details.

Let u (ij) and σ (ij) denote the displacement field and the stress field solution of (1.45) with ε = i (ij) . The components of the fourth-order strain-localization tensor A, of the fourth-order effective stiffness tensor L and of the second-order reduced strain-localization tensor a (k) read as:

A ijmn (x) = ε ij (u (mn) (x)), L ijmn = σ (mn) ij , a (k) ij = µ (k) : ε(u (ij)
) .

(1.46) To obtain the influence factors D (kℓ) N and the second-order tensors ρ (k) , M linear elasticity problems have to be solved: k) . Let u (ℓ) denote the displacement field solution of (1.47) with µ = µ (ℓ) . Note that ρ (ℓ) is the stress field solution of (1.47). Then:

σ(x) = L(x) : (ε(u(x)) -µ(x))) , div (σ(x)) = 0, ε = 0, (1.47) with µ = µ (
D (kℓ) N = µ (k) : ε(u (ℓ) ) . (1.48)
The Finite Element Method (FEM) was used in the two examples presented in section 1.3.9 and 1.4.4 to solve the linear elasticity problems (1.45) and (1.47).

Time-integration of the NTFA model. Strain control

This section is devoted to the time-integration of the NTFA model at the level of a single macroscopic material point when the individual constituents are elastoviscoplastic (the reader is referred to Michel and Suquet [START_REF] Michel | Computational analysis of nonlinear composite structures using the Nonuniform Transformation Field Analysis[END_REF] for rate-independent elasto-plasticity). The history of macroscopic strain ε(t) is prescribed on the time interval [0, T ].

The equations (1.39) to be solved form a system of nonlinear differential equations. Its time-integration requires a time-discretization and an iterative procedure within each time-step. The time interval [0, T ] is decomposed into a finite number of time-steps [t, t + ∆t]. All reduced variables at time t are assumed to be known. The reduced variables and the macroscopic stress at time t + ∆t are obtained as follows.

Time step t + ∆t, iterate i + 1:

The reduced strains (e (k) ) i t+∆t , k = 1, .., M being known,

• Step 1: Compute the plastic multipliers (p (r) ) i t+∆t , r = 1, ..., N , the reduced stresses (τ (k) ) i t+∆t and the reduced back-stresses (x (k) ) i t+∆t , k = 1, .., M (see following paragraph).

• Step 2: Compute the reduced viscoplastic strains (e

(k) vp ) i t+∆t and (ε (k) vp ) i t+∆t . For k = 1, .., M : (e (k) vp ) i t+∆t = (e (k) ) i t+∆t - (τ (k) ) i t+∆t 2G (k) , (ε (k) vp ) i t+∆t = (e (k) vp ) i t+∆t µ (k) : µ (k) .
• Step 3: Compute the macroscopic stress σ i t+∆t :

σ i t+∆t = L : ε t+∆t + M k=1 ρ (k) (ε (k) vp ) i t+∆t .
• step 4: Update the reduced strains (e (k) ) i t+∆t . For k = 1, .., M :

(e (k) ) i+1 t+∆t = a (k) : ε t+∆t + M ℓ=1 D (kℓ) N (ε (ℓ) vp ) i t+∆t .
Go to 1.

The convergence test reads:

Max k = 1, ..., M (e (k) ) i+1 t+∆t -(e (k) ) i t+∆t < δ ε t+∆t -ε t .
A typical value for δ is δ = 10 -6 . The norm for second-order tensors used in the right-hand-side of the convergence test is a = max i,j |a ij | .

Step 1 in details.

In order to determine the plastic multipliers (p (r) ) i t+∆t , the reduced stresses (τ (k) ) i t+∆t and the reduced back-stresses (x (k) ) i t+∆t at step 1 of the above described procedure, the last four equations (1.39) are re-written in the form of a first-order differential equation for these three unknowns. This is done for each phase separately. For a given phase r, the differential system to be solved in the time interval [t, t + ∆t] can be written as:

ẏ = f (y), (1.49)
where the initial data at time t (beginning of the time interval) is known from the previous time step, and where:

y = {y s } s=1,2M (r)+1 , f = {f s } s=1,2M (r)+1 , {y s } s=1,M (r) = {τ (k) } k=1,M (r) , {f s } s=1,M (r) = 2G (k) ė(k) - 3 2 ṗ(r) τ (k) -x (k) A (r) k=1,M (r)
,

{y s } s=M (r)+1,2M (r) = {x (k) } k=1,M (r) , {f s } s=M (r)+1,2M (r) = H (k) τ (k) -x (k) A (r) -ξ (k) x (k) ṗ(r) k=1,M (r)
,

{y s } s=2M (r)+1 = {p (r) }, {f s } s=2M (r)+1 = ṗ(r) .                                               
(1.50) In (1.50), the generalized force A (r) is known according to the τ (k) 's and the x (k) 's, k = 1, .., M (r), by (1.35), the plastic multiplier ṗ(r) according to A (r) and p (r) by (1.36), and the strain-rates ė(k) , k = 1, .., M (r), are given by:

ė(k) = (e (k) ) i t+∆t -(e (k) ) t ∆t . ( 1

Time-integration of the NTFA model. Stress control

It is often convenient (or necessary) to impose the direction of the macroscopic stress σ:

σ t = λ(t)Σ 0 , (1.52)
where Σ 0 is the imposed direction of stress. This is typically the situation which is met in the simulation of the response of a composite to uniaxial tension.

In rate-independent plasticity, especially in ideal plasticity, or in viscoplasticity with power-law materials such as those considered in section 1.3.9, it is not appropriate to control directly the level of stress λ(t). An arc-length method is preferable (Michel et al [START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF][START_REF] Michel | A computational method for linear and nonlinear composites with arbitrary phase contrast[END_REF] ) and the loading is applied by imposing:

ε t : Σ 0 = Ė0 t.
where Ė0 is the imposed strain-rate (in the direction of the applied stress). As in the strain-controlled method, all reduced variables at time t are assumed to be known. At time t + ∆t, the condition ε t+∆t : Σ 0 = Ė0 (t + ∆t) is imposed. The macroscopic stress λ(t + ∆t) is to be determined in addition to the reduced variables. An iterative procedure is used to impose the direction of stress (1.52) as follows:

Time step t + ∆t, iterate i + 1:
The reduced strain (e (k) ) i t+∆t , k = 1, .., M being known and a macroscopic strain ε i t+∆t meeting the condition ε i t+∆t : Σ 0 = Ė0 (t + ∆t) being known: • 1, 2: Perform steps 1, 2 of the procedure described in section 1.3.7.

• 3, 4: Perform steps 3 and 4 of the same algorithm using ε i t+∆t in place of ε t+∆t .

• 5: Compute the level of macroscopic stress:

λ i t+∆t =
Σ 0 : L -1 : σ i t+∆t Σ 0 : L -1 : Σ 0 , and update ε i t+∆t :

ε i+1 t+∆t = ε i t+∆t + L -1 : λ i t+∆t Σ 0 -σ i t+∆t . Go to 1.
The test used to check convergence now reads:

Max k = 1, ..., M (e (k) ) i+1 t+∆t -(e (k) ) i t+∆t < δ ε i+1 t+∆t -ε t , ε i+1 t+∆t -ε i t+∆t < δ ε i+1 t+∆t . A typical value for δ is δ = 10 -6 .

Example 1: Effective response of a dual-phase inelastic composite

The composite materials under consideration in this section are two-phase composites where the two phases play similar (interchangeable) roles in the microstructure.

Material data

Both phases are elastoviscoplastic with linear elasticity and a power-law viscous behavior (corresponding to the dissipation potential (1.9) with R = 0). The material characteristics of phase 1 and 2 read respectively:

E (1) = 100 GPa, ν (1) = 0.3, σ (1) 0 = 250 MPa, ε0 = 10 -5 s -1 , n 1 = 1, and
E (2) = 180 GPa, ν (2) = 0.3, σ (2) 0 = 50 MPa, ε0 = 10 -5 s -1 .
The rate-sensitivity exponent n 2 of phase 2 is varied from 1 to 8. This variation corresponds to the fact that the rate-sensitivity exponent varies significantly with temperature and our objective here is to assess the accuracy of the NTFA used with a single set of plastic modes determined independently for an intermediate value of the rate-sensitivity exponent.

Microstructure

The two-dimensional unit-cell consists of 80 "grains" in the form of regular and identical hexagons. The material properties of these hexagons are prescribed randomly to be either that of phase 1 or that of phase 2 under the constraint that both phases have equal volume fraction (c 1 = c 2 = 0.5). 25 different configurations have been generated (same configurations as in Michel et al [START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF] ). One configuration has been selected among these 25 realizations, namely the one which gives, when n 1 = n 2 = 1 and when the phases are incompressible, an effective response which the closest to the exact result for interchangeable microstructures (given by the self-consistent scheme). This configuration is shown in Figure 1.2 (phase 1 is the darkest phase). Each hexagon is discretized into 64 eight-node quadratic finite elements with 4 Gauss points (5120 quadratic elements and 15649 nodes in total). The unit-cell is subjected to an in-plane simple shear loading with a uniform strain-rate γ = √ 3 ε0 :

ε(t) = γ(t) 2 (e 1 ⊗ e 2 + e 2 ⊗ e 1 ), γ(t) = γ t. (1.53)
Periodic boundary conditions are applied on the boundary of the unit-cell. All computations are performed with the same global time-step ∆t = 2/ √ 3 s.

Plastic modes

The plastic modes retained for application of the NTFA are given by the Karhunen-Loève procedure from two initial fields θ (k) (x) in each phase corresponding respectively to viscoplastic strain fields determined numerically at small strains (γ = 0.03%) and at large strains (γ = 12%). The procedure delivers four orthogonal modes, two modes with support in phase 1 and two modes with support in phase 2. Snapshots of the equivalent strain µ

Discussion of the results

The macroscopic stress-strain response (σ 12 versus γ) is shown in Figure 1.5 when n 2 = 1 and n 2 = 8. The full-field computation which serves as the reference is shown as a solid line. NTFA (1) refers to the NTFA model with a single mode in each phase (the viscoplastic strain field at large strains) whereas NTFA(2) refers to the NTFA model with two modes per phase. If the model NTFA(1) predicts accurately the asymptotic stress response at large strains, the model NTFA( 2) is required for a better agreement in the transient regime where elastic and viscous effects are of comparable order, since, as can be seen in Figures 1.3 and 1.4, the features of the modes for small and large strains are rather different.

The variation of the macroscopic asymptotic stress (creep stress at constant strain-rate) is shown in Figure 1.6 as a function of the rate-sensitivity exponent m 2 = 1/n 2 of phase 2. The full-field computations are shown as stars. The solid line corresponds to the NTFA model implemented with plastic modes which vary with n 2 (viscoplastic strain fields are computed for each value of n 2 and the corresponding plastic modes are deduced by means of the Karhunen-Loève procedure). The results shown as NTFA(n 2 = n) were obtained by considering a single set of modes identified once for all with n 2 = n. NTFA(n 2 = 1) overestimates the macroscopic creep stress of the composite for large values of n 2 . The snapshot of the modes for n 2 = 8 shows a rather significant amount of strain localization in phase 2. NTFA(n 2 = 8) overestimates the creep stress for small nonlinearity (which is consistent with the property of minimization of the dissipation potential). NTFA(n 2 = 3) is a reasonable compromise.

Application of the NTFA to structural problems 1.4.1. Implementation of the NTFA method

The implementation of the NTFA method consists of four different steps. The first two steps are "material steps" in the sense that they are concerned only with computations performed at the unit-cell level, independently of any macroscopic structural problem, except for the choice of the modes which is influenced by the type of macroscopic stress that the material is likely to sustain (as explained in section 1.3.5). These two first steps can be performed once for all. The two last steps are the structural computation itself and a localization step which is essential in the prediction of more local phenomena (such as the lifetime of the structure in fatigue). The four steps are as follows:

Step A: Prior to the resolution of any structural problem, choices and preliminary computations have to be made following sections 1. ρ (k) entering the expression of the macroscopic stress (1.41). This is done once for all by solving linear elasticity problems on the unit-cell (see section 1.3.6).

Step B: Set up a time-integration scheme to integrate the constitutive relations (1.39) along a prescribed path of macroscopic strain ε(t) or macroscopic stress σ(t). This can be done using the schemes proposed in sections 1.3.7 and 1.3.8.

Step C: Incorporate the NTFA model (or more specifically the time-integration scheme of step B) into a Finite Element code. Find the history of macroscopic stresses σ(X, t) and strains ε(X, t) at every macroscopic material point X in the structure.

Step D: It is often useful to determine the local strains and stresses ε(X, x) and σ(X, x) in the actual composite structure and not only the macroscopic strain and stress ε(X) and σ(X) (which are a smoother fields, being averages of the corresponding local fields over a volume element). This localization step is greatly simplified by the NTFA.

Unlike in the exact homogenized problem where the microscopic and macroscopic variables are closely coupled, all steps can be performed independently in the present approach. Steps A and B have already been discussed in section 1.3 and we shall concentrate the discussion on steps C and D.

Implementation of the NTFA in a Finite Element code (step C)

This section deals with the incorporation of the NTFA in a Finite Element code to solve a structural problem. After discretization of the structure into macroscopic finite elements, the unknowns pertaining to the structural (i.e. macroscopic) problem are denoted by overlined letters e.g. u(X), σ(X).... Arrays of discrete unknowns are denoted with braces, e.g. {u} denotes the array of discrete unknowns associated with the displacement field u, matrices are denoted with brackets, e.g.

[K] denotes the assembled stiffness matrix associated with the effective stiffness of the composite L.

The structural problem is solved incrementally after time discretization of the interval of study. All significant variables (displacement, stresses) being known at time t, the unknowns at time t + ∆t are determined by the equilibrium equations and the macroscopic (or homogenized) constitutive relations.

Equilibrium of the structure implies Go to step a.

T {v -u t+∆t } {R} t+∆t = 0, {R} t+∆t = -

The global stiffness matrix [K]

I t+∆t can be chosen among many different possibilities. One of the simplest one, which was used in the example presented in section 1.4.4, is the initial elastic stiffness: No particular convergence problem was observed with this elementary method.

In the convergence test used in step b, the norm of the equilibrium residues is checked:

max j {|R j |} I t+∆t < δ max j ′ {|R j ′ |} I t+∆t , (1.56)
where {R} I t+∆t denotes the array of reactions at nodal points on the boundary of the structure. A typical value for δ is δ = 10 -6 .

Computation of {σ}

I t+∆t

Consider an integration point X of a finite element e in the structure. The strain at X reads {ε(X)}

I t+∆t = [B(X)] {u e } I t+∆t , {u e } I t+∆t = {u e } t + {∆u e } I t+δt .
Then the iterative procedure of section 1.3.7, applied with ε t+∆t = {ε(X)} I t+∆t , delivers the stress {σ(X)} I t+∆t at point X.

Note that the procedure of section 1.3.7 requires the knowledge of the reduced variables at time t. These variables are (p r ) t , r = 1, ...N and (e (k) ) t , (τ (k) ) t , (x (k) ) t , k = 1, ...M . It is therefore necessary to store these scalar variables at each integration point of each finite element of the structure.

Localization rules

The strain and stress fields ε(X) and σ(X) delivered by the structural analysis are averaged fields. Their value at a macroscopic point X is the average over the microscopic variable x of the zero-th order terms ε 0 (X, x) and σ 0 (X, x) in the expansion of the strain and stress fields , when x varies in the unit-cell. The averaged fields do not capture the rapid oscillations (and most importantly the peaks) of the actual strain and stress fields at the microscopic scale.

Mathematical analysis shows that these zero-th order terms in the asymptotic expansion (1.3) provide, after rescaling, a better approximation of ε η (X) and σ η (X) than ε(X) and σ(X) by setting:

εη (X) = ε 0 X, X η , ση (X) = σ 0 X, X η . (1.57)
In linear elasticity it has been shown theoretically (Suquet 27 ) and observed numerically (Feyel and Chaboche [START_REF] Feyel | FE2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials[END_REF] ) that εη and ση are pointwise approximations of ε η and σ η and not only weak approximations (as are ε and σ), except in a boundary layer close to the boundary of the structure where the periodicity conditions can be in contradiction with the actual boundary conditions (boundary layer terms must be added to have a good approximation up to the boundary).

In linear elasticity the zero-th order terms ε 0 and σ 0 in the expansion of the strain and stress fields are nothing else than the local fields ε and σ and are therefore related to their average by means of the localization tensors A and B:

ε 0 (X, x) = A(x) : ε(X), σ 0 (X, x) = B(x) : σ(X).
(1.58) Therefore a good approximation of the actual strain and stress fields can be obtained by solving independently the structural problem to find the macroscopic fields ε(X) and σ(X) and six unit-cell problems to find the stress-localization tensors A and B. Then the two results are combined by means of (1.58) to give a good approximation of the actual strain and stress fields in the composite structure (with a possible exception at the boundary, as discussed in section 1.2.1). In other words, the local fields ε(X, x) and σ(X, x) (or good approximations of them) can be obtained by post-processing the macroscopic strain and stress fields ε(X) and σ(X). A full decoupling of scales can be achieved.

In nonlinear problems, and in particular in presence of viscoplasticity or plasticity, no simple relation such as (1.58) exists. Rigorously speaking, there is no explicit decoupling of scales. If no approximation is made, the microscopic fields ε 0 (X, x) and σ 0 (X, x) are intimately coupled to the macroscopic fields ε(X) and σ(X) and all microscopic and macroscopic fields must be determined in the course of a coupled computation. The field localization is not performed as a postprocessing step but as part of the structural analysis. As already underlined, the cost of this computational procedure can be formidable.

The NTFA avoids this complication, thanks to the relations (1.28) and (1.40), admittedly at the expense of the approximation (1.24). First, as shown in section 1.4.2 the structural problem is solved independently of the unit-cell calculations (performed once for all). Second, the microscopic fields are deduced from their macroscopic counterpart by means of the explicit and linear relations (1.28) and (1.40):

ε(X, x, t) = A(x) : ε(X, t) + M k=1 η (k) (x)ε (k) vp (X, t), σ(X, x, t) = L(x) : A(x) : ε(X, t) + M k=1 ρ (k) (x)ε (k) vp (X, t).            (1.59)
The macroscopic state variables (ε(X), ε

vp (X)) are outputs of the structural computation performed with the homogenized NTFA model. The relation (1.59) can be used to post-process these fields and obtain an accurate approximation of the actual strain and stress fields ε η and σ η by setting:

εη (X, t) = A X η : ε(X, t) + M k=1 η (k) X η ε (k) vp (X, t), ση (X, t) = L X η : A X η : ε(X, t) + M k=1 ρ (k) X η ε (k) vp (X, t).            (1.60)

Example 2: Fatigue of a metal-matrix composite structure

In this section the NTFA model is applied to a structural problem. A plate composed of a inner core (thickness 4 mm), made of a metal-matrix composite, surrounded by two outer layers of pure matrix (thickness 0.5 mm each) is subjected to a cyclic four-point bending test. By symmetry only half of the plate is considered as shown in Figure 1.1 where the unit-cell generating the core of the plate by periodicity is also shown.

The matrix is elasto-viscoplastic with purely nonlinear kinematic hardening (the isotropic hardening is negligible R(p) = σ y ):

E m = 60 GPa, ν m = 0.3, σ y = 20 MPa, n = 5, η = σ 0 ε-1 n 0 = 100 MPa s 1 n , H = 10 GPa, ξ = 1000 MPa.      (1.61)
The metal matrix is reinforced by long circular fibers arranged at the nodes of a square array. The fiber volume fraction is 25%. The fibers are linear elastic with Young modulus and Poisson ratio:

E f = 300 GPa, ν f = 0.25. (1.62)
The plate is simply supported at points B and B ′ and periodic (in time) displacements at points A and A ′ are prescribed. Depending on the amplitude of the displacement, the structure is likely to undergo viscoplastic deformations leading to fatigue failure. There are three possible failure mechanisms at the microscopic scale: fiber breakage, fiber-matrix debonding and matrix failure. At high temperature, when the matrix is viscoplastic as considered in this study, matrix damage is the dominant mechanism (Llorca 13 ). Therefore a first modeling assumption is that failure of the composite occurs by matrix failure. To predict matrix failure, a model due to Skelton [START_REF] Skelton | Energy criterion for high temperature low cycle fatigue failure[END_REF] for low-cycle fatigue is used (a comparison of different lifetime prediction methods including Skelton's model can be found in Amiable et al [START_REF] Amiable | A comparison of lifeprediction methods for a thermal fatigue experiment[END_REF] ). The model is based on the energy dissipated by viscoplasticity during the stabilized cycle:

w = cycle σ : εvp dt.
(1.63) Skelton's model is based on the assumption (confirmed experimentally) that the number of cycles to failure N f for a material under cyclic thermomechanical fatigue tests in the low-cycle regime is related to the energy dissipated w by:

w N β f = C, (1.64)
where C and β are material constants independent of the thermomechanical loading.

In the framework of these two working assumptions (failure of the composite governed by matrix failure, and matrix failure governed by the criterion (1.64)), one can predict the lifetime of the composite structure subjected to a cyclic thermomechanical loading at the expense of resolving the stress and strain fields at the smallest scale in order to apply the criterion (1.64). This procedure is extremely heavy and the aim of this section is to demonstrate that an accurate prediction can be obtained by means of the NTFA at a much reduced cost, involving only a purely macroscopic computation, followed by a proper postprocessing of the macroscopic fields.

Meshes

The fine mesh accounting for all microstructural details of the heterogeneous structure is shown in Figure 1.7a. The mesh of the inner core is obtained by repeating the mesh of the unit-cell shown in Figure 1.7c which consists of 80 six-node triangular elements (3 Hammer points) in the fiber and 128 eight-node quadrilateral elements (4 Gauss points) in the matrix, for a total of 208 elements and 577 nodes. The same unit-cell mesh was used for the unit-cell preliminary computations (effective properties, plastic modes, influence factors, localization fields A and η (k) ). The resulting mesh for the heterogeneous structure consists of 26880 quadratic elements (6 or 8 nodes) and 71601 nodes in total. The mesh used in the homogenized computations is shown in Figure 1.7b and consists of only 600 eight-node quadrilateral elements and 1941 nodes.

Loading

The boundary conditions applied to the right half of the cross-section of the plate are (refer to Figure 1.1 for the location of points A, A ′ , B and B ′ ):

X 1 = 0 : u 1 (0, X 2 ) = 0, t 2 (0, X 2 ) = 0, - h 2 ≤ X 2 ≤ h 2 , Point A : t 1 (X A 1 , h 2 ) = 0, u 2 (X A 1 , h 2 ) = u, Point A' : t 1 (X A 1 , - h 2 ) = 0, u 2 (X A 1 , - h 2 ) = u, Point B : t 1 (X B 1 , - h 2 ) = 0, u 2 (X B 1 , - h 2) = 0, Point B' : t 1 (X B 1 , h 2 ) = 0, u 2 (X B 1 , h 2 ) = 0, X 1 = L : t 1 (L, X 2 ) = 0, t 2 (L, X 2 ) = 0, - h 2 ≤ X 2 ≤ h 2 ,                                 
(1.65) with h = 5 mm, L = 30 mm, X A 1 = 10 mm, X B 1 = 25 mm. The traction on the boundary of the structure is denoted by t = σ.N . The vertical displacement u imposed at points A and A ′ is periodic in time with period T . It is a piecewise linear function of time, varying linearly between u max and -u max as shown in Figure 1.8. The loading frequency f = 1/T is prescribed f = 0.1 Hz, whereas the maximal displacement at points A and A ′ is varied u max = 0.15, 0.2, 0.25, 0.35 and 0.5 mm. The loading frequency being kept constant in the different loading cases, varying the maximal displacement prescribed to A and A ′ results in different velocities for these points and therefore in different strain-rates in the structure. All computations were performed with a global time-step ∆t = 0.25 u max 10 -2 s. The choice of the modes depend in general on the type of loading that the structure is likely to undergo. Although it is expected that the dominant stress will be uniaxial tension-compression in the horizontal direction, transverse shear and even transverse normal stress cannot be excluded. So the three types of stress (horizontal, vertical and shear) will be considered in the analysis leading to the choice of the modes. The NTFA model is implemented with five plastic modes in the matrix, and the macroscopic model has therefore five internal variables . These modes were obtained by subjecting the unit-cell to cyclic loading along three different directions of macroscopic stress:

Σ (1) = e 1 ⊗ e 1 + Σ (1)
33 e 3 ⊗ e 3 , ε 33 = 0,

Σ (2) = e 1 ⊗ e 2 + e 2 ⊗ e 1 + Σ (2)
33 e 3 ⊗ e 3 , ε 33 = 0,

Σ (3) = e 2 ⊗ e 2 + Σ (3) 33 e 3 ⊗ e 3 , ε 33 = 0.              (1.66)
The components

Σ (i)
33 are a priori left free and determined a posteriori as the reactions to the constraint ε 33 = 0. The computations at the unit-cell level are performed in plane strains, in concordance with the plane strain conditions which prevail at the structural level.

The unit-cell is subjected to a cyclic loading along each of the three stress directions (1.66). The problem is strain-controlled (as described in section 1.3.8). The macroscopic strain in the imposed stress direction varies between ε max : Σ (i) and -ε max : Σ (i) , with ε max : Σ (i) = 0.0025, i = 1, 2, 3. The variation of the macroscopic strain in time is a triangular profile similar to that shown in Figure 1.8 where the prescribed strain-rate is ε : Σ (i) = 10 -3 s -1 , i = 1, 2, 3. All computations are performed with the same global time-step ∆t = 10 -2 s until the response of the material point undergoing the largest viscoplastic dissipation (as defined through the scalar quantity (1.63)) reaches a stabilized cycle.

For each of the three loading cases (1.66) the viscoplastic strain fields at each quarter of all cycles are stored. In other words, for a given cycle c beginning at time t c and with period T , the viscoplastic strain fields at time t c , t c + T /4, t c + T /2 and t c + 3T /4 are stored. This is done for all cycles until the "hottest" point in the unit-cell reaches a stabilized cycle. Let θ

(k) i (x), k = 1, ..., M (i)
T , i = 1, 2, 3 denote the whole set of fields stored according to this procedure. M (i) T denotes the total number of fields stored along the i-th loading direction Σ (i) . The number of modes is first reduced for each loading direction by applying the Karhunen-Loève procedure described in section 1.3.5 separately to the three family of fields θ (k) i (x), i = 1, 2, 3. The modes with the highest intensity (corresponding to the highest eigenvalue of the correlation matrix) are extracted for each loading case. The five modes finally retained for further use in the NTFA are the shear mode (macroscopic stress being a pure shear) with the highest intensity and the two modes with highest intensity for the two other loading cases (tension-compression in the horizontal and vertical direction respectively). Taken separately, these modes are sufficient for the NTFA to reproduce accurately the response of the unit-cell along the loading direction from which they were extracted. Lastly, since these five modes were selected independently, they do not necessarily meet the orthogonality condition (1.44). Another application of the Karhunen-Loève procedure is made, leading finally to five modes satisfying all the desirable requirements. Snapshots of the equivalent strain of the five modes are given in Figure 1.9. A more local comparison can be performed by examining the stress-strain response, not of the whole unit-cell as was done in Figure 1.10, but at the material point in the unit-cell undergoing the largest dissipated energy (1.63). This is done for uniaxial tension-compression in the horizontal direction in Figure 1.11. Again the agreement is seen to be excellent.

Finally, it is also of interest to compare the prediction of the model for the energy dissipated along the stabilized cycle with full-field simulations. This is done in Figure 1.12. The model makes use of the localization rules (1.60) to estimate the energy (1.63). The NTFA model captures well the local distribution of the energy dissipated in the unit-cell. The energy turns out to be maximal at the fiber-matrix interface. The reference FEM simulation gives w max = 2.134 MPa, whereas the NTFA model predicts w max = 2.231 MPa. ergy dissipated along the stabilized cycle, for two different structural simulations:

a) The first simulation is performed with a very fine mesh of the heterogeneous structure 1.7a) and accounts for all detailed heterogeneities. It is considered as the exact response of the composite structure with a small but non-vanishing value of η.

b) The second simulation is performed on a coarse mesh, using at each integration point of the mesh the homogenized NTFA model.

A first element of comparison is provided in Figure 1.13 where the forcedisplacement (the force is the sum of the reactions at points A and A ′ ) response of the structure predicted by the homogenized NTFA model (grey dashed line) is compared to the detailed simulation with full account of the heterogeneities (black solid line). The two graphs correspond to three different values of the maximal displacement u max = 0.25 mm (at left) and u max = 0.15 and 0.5 mm (at right). The agreement is good in all cases. A more local comparison can be made by examining the response of the most severely loaded unit-cell in the structure (where the energy dissipated is maximal). The stress and strain fields for the NTFA model are obtained by means of relations (1.60). The quantities used for comparison in Figure 1.14 are the stress and strain averaged on this particular unit-cell. The agreement in the stress level is rather good, but the NTFA seems to slightly overestimate the amount of local strain.

Finally, as exposed in the introduction of this section, the quantity of interest here is the lifetime of the structure which is directly related to the energy dissipated at the "hottest" point in the structure through the relation (1.64). The use of the NTFA model raises two questions: 2) Is the amount of energy dissipated correctly estimated by the model?

To answer these questions, the heterogeneous Finite Element analysis and the macroscopic structural simulation using the homogenized NTFA model are run until the structure reaches a stabilized cycle. The energy dissipated along this stabilized cycle is directly available in the heterogeneous simulation. In the NTFA model it can be directly deduced from the macroscopic results by means of the localization rules (1.60). To answer the first question, the two snapshots (fullfield computation and NTFA model) of the energy w over the whole structure are shown in Figure 1.16 (u max = 0.25 mm). This very local quantity is reasonably well predicted by the NTFA model. A close-up of the same energy distribution in the region where w is maximal is shown in Figure 1.17. As can be seen from these figures, the location of the hottest point is well predicted by the NTFA model. To answer the second question, the stabilized cycles at the hottest point in the structure are shown in Figure 1.15. Given the very local character of this information, the agreement of the model's prediction with the detailed computation can be considered as good, the model overestimating the amount of strain at this hottest point. A further pointwise comparison of the maximum w max of the energy is shown in Figure 1.18. Independently of maximal displacement prescribed to the structure, the NTFA overestimates by about 25% the maximum of the dissipated energy (this estimation is related to the overestimation of the strain at the hottest point). Therefore the lifetime of the structure will be underestimated by a similar amount, which is a quite reasonable error (on the safe side), given the fact that no coupled multiscale computation is required by the NTFA model but only a postprocessing of a purely macroscopic simulation.

Conclusion

The Nonuniform Transformation Field Analysis is a newly proposed micromechanical scheme for multiscale problems with nonlinear phases. This model is based on a drastic reduction of the number of variables describing the microscopic (visco)plastic strain field performed by means of the Karhunen-Loève procedure (proper orthogonal decomposition). It delivers effective constitutive relations for nonlinear composites expressed in terms of a small number of internal variables which are the components of the microscopic plastic field over a finite set of plas-tic modes. This reduced model can be easily incorporated in a structural computation. A numerical scheme is proposed to integrate in time the homogenized constitutive relations at each integration point of the structure. The predictions of the model compare well to results of large-scale computations over the whole composite structure, accounting for all detailed information. The agreement is good not only in terms of global quantities (force/displacement) but also in terms of local quantities. For instance the lifetime of a structure subjected to cyclic loading has been predicted with a fatigue criterion based on the energy dissipated along a cycle in the matrix. The agreement between the model and the large-scale heterogeneous computation is very good.

Fig. 1 . 1 .

 11 Fig. 1.1. Composite structure (left) and unit-cell (right).

  vp (t) in each phase can be obtained by integrating in time the systems of differential equations (1.19) 2 , (1.20) and (1.23).

Fig. 1 . 2 .

 12 Fig. 1.2. Covering of the unit-cell by regular hexagons of phase 1 (dark) and phase 2. Realization used for the implementation of the NTFA.

  eq of the four modes are shown in Figures 1.3 and 1.4.

Fig. 1 . 3 .

 13 Fig. 1.3. Dual-phase material. Plastic modes for phase 1. Snapshot of the equivalent strain µ (k) eq , k = 1, 2. At top: n 2 = 1. At bottom: n 2 = 8. From left to right: modes for small and large strains. The look-up table is the same for all four snapshots.

Fig. 1 . 4 .

 14 Fig. 1.4. Dual-phase material. Plastic modes for phase 2. Snapshot of the equivalent strain µ (k) eq , k = 3, 4. At top: n 2 = 1. At bottom: n 2 = 8. From left to right: modes for small and large strains. The look-up table is the same for all four snapshots.

12 Fig. 1 . 5 .Fig. 1 . 6 .

 121516 Fig. 1.5. Dual-phase material. Response of the unit-cell under macroscopic shear deformation (1.53). At left: n 2 = 1. At right: n 2 = 8.

T

  [B] {σ} t+∆t dX , (1.54) where v is an arbitrary kinematically admissible displacement field, [B] is the classical FE matrix relating displacements and strains, i.e. {ε e } = [B]{u e }, and e denotes a finite element. Equation (1.54) is a nonlinear equation which can be solved by an iterative Newton scheme as follows. Time step t + ∆t, iterate I + 1: {∆u} I t+∆t being known at each nodal point of the structure, • Step a: Compute the stresses {σ} I t+∆t at each integration point of each finite element of the structure (see paragraph below). • Step b: Check convergence. If convergence is not reached, solve the linear system:

  ] , where [k e ] = e T [B][ L] [B] dX. (1.55)

Fig. 1 . 7 .

 17 Fig. 1.7. Meshes used in the analysis of the composite plate shown in figure 1.1. (a): fine mesh of the heterogeneous structure. (b): Coarse mesh used for the analysis of the homogenized structure by means of the NTFA model. (c): mesh of the unit-cell generating, by periodicity, the mesh of the inner core of the plate as shown in (a).

Fig. 1 . 8 .

 18 Fig. 1.8. History of the prescribed displacements at points A and A ′ .

Fig. 1 . 9 .Fig. 1 .

 191 Fig. 1.9. Snapshots of the equivalent strain µ (k)eq , k = 1, ..., 5 for the five orthogonal plastic modes in the matrix. The look-up table is the same for all five snapshots.

1. 4 . 4 . 5 .Fig. 1 .

 4451 Fig. 1.11. Unit-cell response. Comparison between full-field FEM computations (black solid line) and the NTFA model with the five modes shown in Figure 1.9 (grey dashed line). Stress-strain response at the hottest point in the unit-cell. Tension-compression in the horizontal direction (loading case i = 1 in (1.66)).

Fig. 1 .

 1 Fig. 1.12. Unit-cell response. Snapshot of the energy w dissipated in the unit-cell by viscoplasticity along the stabilized cycle. Uniaxial horizontal tension-compression. At left: Reference full-field FEM simulation. At right: Prediction of the NTFA model. The look-up table is the same for both snapshots.

15 Fig. 1 .

 151 Fig. 1.13. Four-point bending. Comparison between the heterogeneous Finite Element analysis (black solid line) and the NTFA homogenized model (grey dashed line). Global force-displacement response. At left: umax= 0.25 mm. At right: umax= 0.15 and 0.5 mm.

1 ) 15 Fig. 1 .

 1151 Fig. 1.14. Four-point bending. Comparison between the heterogeneous Finite Element analysis (black solid line) and the NTFA homogenized model (grey dashed line). Average-stress/average-strain response of the most solicited unit-cell in the structure. At left: umax= 0.25 mm. At right: umax= 0.15 and 0.5 mm.

15 Fig. 1 .

 151 Fig. 1.15. Four-point bending. Stress/strain response at the hottest point in the structure. Comparison between the heterogeneous Finite Element analysis (black solid line) and the NTFA homogenized model (grey dashed line). At left: umax= 0.25 mm. At right: umax= 0.15 and 0.5 mm.

Fig. 1 .Fig. 1 . 18 .

 1118 Fig. 1.17. Four-point bending. Distribution of the dissipated energy w (normalized by its maximum). Stabilized cycle. umax= 0.25 mm. Close-up in the most severely loaded region. At left: Full heterogeneous simulation (reference). At right: Prediction of the NTFA model using the localization rules.The look-up table is the same for both snapshots.

  uum Micromechanics, volume 377 of CISM Lecture Notes, pages 197-264. Springer Verlag, New York, 1997. 31. L. Tartar. Estimations de coefficients homogénéisés. In R. Glowinski and J.L. Lions, editors, Computing methods in applied sciences and engineering, volume 704 of Lecture Notes in Mathematics, pages 364-373. Springer Verlag, Berlin, 1977. 32. K. Terada and N. Kikuchi. A class of general algorithms for multi-scale analyses of heterogeneous media. Comp. Meth. Appl. Mech. Engng, 190:5427-5464, 2001.