Weijie Wang

Amato Evan

Cyrille Flamant

Christophe Lavaysse

On the decadal scale correlation between African Dust and Sahel rainfall: the role of Saharan heat Low-forced winds

On the decadal scale correlation between African dust and Sahel rainfall: the role of Saharan Heat Low-forced winds

Introduc)on

By mass, aeolian dust is the most pervasive aerosol in the atmosphere. The presence of suspended dust affects the local energy balance through direct and indirect effects, modifies the hydrological cycle via radiative forcing and modification of cloud microphysical properties, and is a vehicle for the long-range transport of nutrients to global oceans and terrestrial land surfaces.

North Africa is the world's largest dust source; accounting for more than 55% of the global dust emissions. Long-term in-situ observations at and coral reef proxy and satellite have shown that dust emission from North Africa peaked during the mid-1980s and has followed a downward trend through at least the late 2000s. Studies have found dust cover over the Atlantic is anticorrelated with precious-year Sahelian precipitation. However, there are several aspects of such a theory that are not consistent. For example, satellite imagery suggests that the vast majority of dustemitting regions lie to the north of the vegetated region of the Sahel.

Here we reconcile this contradiction by showing that both Saharan surface wind fields over the major regional dust emitting regions, and the northward propagation of the monsoon flow and thus Sahel rainfall, are forced by the thermodynamic state of a meteorological feature termed the Saharan Heat Low.

The$Saharan$heat$low$(SHL)

The$SHL$is$defined$as$the$atmospheric$thickness$between$the$700$hPa$and$925hPa levels,$using$geopotenDal$height$from$ERA$Interim$reanalysis.

PrecipitaDon

The$VersionM2$Global$PrecipitaDon$Climatology$Project$(GPCP).

Dust$aerosol$opDcal$depth

Retrieved$from$satellite$radiance$measurements$from$the$Advanced$Very$High ResoluDon$Radiaometer$(AVHRR).

Dust emission rates

Derived from Spinning Enhanced Visible and Infrared Imager infrared dust index images.

Winds

ERA Interim reanalysis

Data,and,Methods

Upper: Seasonal (June-August) time series of Zhl (red), precipitation (P, blue) and td (yellow). The monthly mean precipitation series is averaged over 5° W to 20° E, 10° to 20° N, and td is averaged over 10° to 65° W and 0° to 30°N. The thin lines are the seasonal means and the thick lines are the 5-season smoothed time series (via a 1-4-7-4-1 filter).

Weijie Wang, Amato Evan, Cyrille Flamant and Christophe Lavaysse contact: wjwang@ucsd.edu

Results

Conclusion

Climatology (A) and anomaly (B) of precipitation and winds

We show that interannual variability in Sahelian rainfall, and surface wind speeds over the Sahara, are the result of changes in lower tropospheric air temperatures over the SHL. As the SHL warms an anomalous tropospheric circulation develops that reduces windspeeds over the Sahara and displaces the monsoonal rainfall northward, thus simultaneously increasing Sahelian rainfall and reducing dust emission from the major dust "hot-spots" in the Sahara.

Left:

  Climatology (A) and anomaly (B) of precipitation and winds. (A) Map of annual dust emission, climatological precipitation and 925 hPa climatological wind vectors. Dust emission rates are derived from Spinning Enhanced Visible and Infrared Imager infrared dust index images in units of % of annual North Africa dust emission. The magenta contours are precipitation in units of mm day-1 averaged for the summertime from 1979 through 2012 and wind vectors in units of m s-1 are averaged for the same period of time. The red circle indicates the location of the SHL during the summer. The brown boxes are the major dust hotspots of the Bodélé Depression (14° to 24° E, 14° to 21° N), the depression in the lee of Aïr and Adrar Mountains (2° to 10° E, 16° to 24° N), and the Mauritania and Western Sahara source region (12°W to 2° E, 18° to 27° N). (B) Map of the coefficients of the regression of the 925 hPa climatological wind vector, the wind speed and precipitation onto Zhl for the summertime from 1979 through 2012. Shaded regions with color contours are statistically significant at the 90% level.

Left:

  Histogram counts (per 0.2m/s) of 10 m wind speed from ERA-Interim reanalysis for the 300 strongest (red) and weakest (blue) phases of the SHL in the Bodélé Depression (A), a depression in the lee of Aïr and Adrar Mountains (B), and the Mauritania and Western Sahara (C) in the summertime from 1979 through 2012. The means of the wind speed for strong and weak SHL are statistically significantly different in ABC based on the 2-tailed t-test. Shown in panels DEF are the wind speed histogram differences for the warm minus cool SHL cases. The grey shaded region represents wind speed threshold for dust emission. Right: Difference in aerosol optical depth (DAOD) for 300 warmest minus 300 coolest SHL from MODIS Aqua Version 5, A. T. Evan, C. Flamant, C. Lavaysse, On the decadal scale correlation between African dust and Sahel rainfall: The role of Saharan heat low-forced winds. Sci. Adv. 1, e1500646 (2015).


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

Pierre Suquet
email: suquet@lma.cnrs-mrs.fr

A model-reduction approach to the micromechanical analysis of polycrystalline materials

Keywords: Micromechanics, model reduction, crystal plasticity, creep

The present study is devoted to the extension to polycrystals of a model-reduction technique introduced by the authors, called the Nonuniform Transformation Field Analysis (NTFA). This new reduced model is obtained in two steps. First the local fields of internal variables are decomposed on a reduced basis of modes as in the NTFA. Second the dissipation potential of the phases is replaced by its tangent second-order (TSO) expansion. Thanks to the second approximation the reduced evolution equations of the model can be entirely expressed in terms of quantities which can be pre-computed once for all. Roughly speaking, these pre-computed quantities depend only on the average and fluctuations per phase of the modes and of the associated stress fields. The accuracy of the new NTFA-TSO model is assessed by comparison with full-field simulations on two specific applications, creep of polycrystalline ice and response of polycrystalline copper to a cyclic tension-compression test. The new reduced evolution equations is faster than the full-field computations by two orders of magnitude in the two examples.

Introduction

Micromechanics of polycrystalline materials is primarily concerned with the incorporation of crystallinity, microscopic mechanisms of deformation and micro-structural informations (such as texture, distribution of orientations...) into a continuum, or macroscopic, description of Plasticity. The subject has received a great deal of attention since the early work of Taylor [START_REF] Taylor | The mechanics of plastic deformation of crystals. Part I: Theoretical[END_REF] on the subject. The problem of determining the local stress and strain fields in polycrystals is of formidable complexity. Fortunately, the development of computing capacities permits numerical simulations of the response of polycrystalline aggregates in realistic conditions (threedimensional polycrystals, several hundreds of grains, texture evolution). A considerable body 1 of literature has been devoted to the subject of computational micromechanics of polycrystals since the seminal paper of Asaro [START_REF] Asaro | Micromechanics of crystals and polycrystals[END_REF] (a comprehensive review and more references can be found in [START_REF] Roters | Crystal Plasticity Finite Element Methods[END_REF]).

However these full-field simulations have serious limitations.

1. First, they do not provide macroscopic (or homogenized) constitutive relations. They only describe the response of a polycrystalline aggregate under specified loading paths. They can be used as a building block in a nested FE 2 approach where both microscopic problems and macroscopic problems (which remain coupled in full generality) (see [START_REF] Feyel | FE2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials[END_REF][START_REF] Terada | A class of general algorithms for multi-scale analyses of heterogeneous media[END_REF][START_REF] Kanouté | Multiscale methods for composites: A review[END_REF] for a similar approach in composite materials). But so far the formidable size of the computational problems has limited the use of FE 2 methods to rather academic problems. They can also be used to calibrate a purely phenomenological model at the macroscopic scale ( [START_REF] Pellegrino | Numerical homogenization of periodic composite materials with non-linear material components[END_REF]), provided an appropriate model can be identified. However, it is rather frustrating that one has to postulate, rather than derive, the form of the macroscopic constitutive relations. The best approximation could well be outside the classical engineering models.

2. Second, they are still very costly if they have to be repeated several times, for instance when identifying material parameters. One would prefer to have a small (admittedly approximate) model reproducing as closely as possible the response of the full-field simulations in a prescribed range of macroscopic stress states.

Reduced-order modelling aims precisely at the reduction of problems with a large size to a smaller size. It is, however, important to be specific about what is called a model. It could be a computational model of large size (such as illustrated by the FE 2 method) in which case the reduced-order model is a computational model of smaller size. But it could also be a mechanical (or physical) model, governed by differential equations, in which case the reduced-order model is a set of differential equations for another (reduced) mechanical model, involving much less variables. In this study we aim at deriving a reduced-order model for micromechanics which achieves both aims: identify a suitable set of reduced variables with the differential equations governing their evolution AND reduce the size (and therefore the cost) of computations.

Model-reduction has a long history in Fluid Mechanics (see for instance [START_REF] Sirovich | Turbulence and the dynamics of coherent structures[END_REF][START_REF] Holmes | Structures, Dynamical Systems and Symmetry[END_REF]) and in many other fields of computational physics ( [START_REF] Lucia | Reduced-order modeling: new approaches for computational physics[END_REF]). Its use in Solid Mechanics is more recent (see Ryckelinck [START_REF] Ryckelynck | Multi-level a priori Hyper-reduction of mechanical models involving internal variables[END_REF], Chinesta and Cueto [START_REF] Chinesta | PGD-Based Modeling of Materials, Structures and Processes[END_REF] and the references herein). But, to the authors' knowledge, it has been mostly applied in Solid Mechanics with the objective of reducing the cost of computations and rarely with the objective of deriving reduced constitutive relations.

The problem at hand in polycrystalline materials (the same holds composite materials), is that the macroscopic response depends intimately on the whole local microscopic fields, plastic slips in polycrystals, or plastic strain fields in composite materials (and other fields of internal variables describing the hardening of the material) which have to be computed and stored in the course of the macroscopic computations. This is extremely expensive. One of the earliest, and pioneering, attempt to reduce the complexity of the plastic strain fields in micromechanics of materials is the Transformation Field Analysis (TFA) of Dvorak [START_REF] Dvorak | Transformation field analysis of inelastic composite materials[END_REF] which assumes uniformity of the plastic strain in the phases or in subdomains. It has been further developed by Dvorak himself [START_REF] Dvorak | The modeling of inelastic composite materials with the transformation field analysis[END_REF], extended to periodic composites by Fish [START_REF] Fish | Computational plasticity for composite structures based on mathematical homogenization: Theory and practice[END_REF] and has 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, the assumption of piecewise uniform plastic strain fields is far from reality and in order to reproduce accurately the actual effective behavior of the composite, it is essential to capture correctly the heterogeneity of the plastic strain field which requires a large number of subdomains.

This last observation has motivated the introduction by Michel and Suquet [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 the Nonuniform Transformation Field Analysis (NTFA) where the (visco)plastic strain field within each phase is decomposed on a finite set of plastic modes which can present large deviations from uniformity. The reduced variables are the components of the (visco)plastic strain field on the (visco)plastic modes. Approximate evolution laws for these variables have been be proposed (Michel and Suquet [30,[START_REF] Michel | Nonuniform Transformation Field Analysis: a reduced model for multiscale nonlinear problems in Solid Mechanics[END_REF]). It has been subsequently shown by the authors [START_REF] Michel | A model-reduction approach in micromechanics of materials preserving the variational structure of constitutive relations[END_REF] that these evolution laws can be interpreted in the framework of a secant formulation of homogenization. A significant advantage of the NTFA is that it provides localization rules allowing for the reconstruction of local fields which are used to predict local phenomena such as the distribution of stresses or the plastic dissipation at the microscopic scale (Michel and Suquet [START_REF] Michel | Nonuniform Transformation Field Analysis: a reduced model for multiscale nonlinear problems in Solid Mechanics[END_REF]). This model, which will be called the original NTFA model, was first applied to two-dimensional situations by the authors ( [START_REF] Michel | Nonuniform Transformation Field Analysis[END_REF][START_REF] Michel | Computational analysis of nonlinear composite structures using the Nonuniform Transformation Field Analysis[END_REF][START_REF] Michel | Nonuniform Transformation Field Analysis: a reduced model for multiscale nonlinear problems in Solid Mechanics[END_REF]). It has subsequently been applied to three-dimensional problems by Fritzen and Bölhke [START_REF] Fritzen | Three-dimensional finite element implementation of the nonuniform transformation field analysis[END_REF] and extended to phases with transformation strains by Largenton et al [START_REF] Largenton | Extension of the nonuniform transformation field analysis to linear viscoelastic composites in the presence of aging and swelling[END_REF]. A step towards a more rational derivation of the evolution equation for the reduced variables has been achieved by Fritzen and Leuschner [START_REF] Fritzen | Reduced basis hybrid computational homogenization based on a mixed incremental formulation[END_REF], who proposed a hybrid form of the incremental variational principles for materials governed by two potentials. However their hybrid model does not decouple the microscopic and macroscopic scales and does not lead to explicit constitutive relations which can be expressed only in terms of "pre-computed" quantities (computed "off-line" once for all). To achieve this final (and necessary) step towards a macroscopic model which can be pre-computed, the authors have proposed a tangent second-order linearization of the equations incorporating averages and quadratic fluctuations of the local fields ( [START_REF] Michel | A model-reduction approach in micromechanics of materials preserving the variational structure of constitutive relations[END_REF]). These averages and fluctuations can be computed from those of a few pre-computed fields with the help of the reduced variables, resulting in a compact macroscopic model. This method, which was first proposed for composite materials, is extended here to polycrystalline materials.

The paper is organized as follows. The general class of constitutive relations for the in-dividual grains considered in this study is recalled in section 2 and illustrated by two specific examples of crystal plasticity models, the model of Méric and Cailletaud [START_REF] Méric | Single crystal modeling for structural calculations. part 2: Finite element implementation[END_REF] for metals and the model of Castelnau et al [START_REF] Castelnau | Elastoviscoplastic micromechanical modeling of the transient creep of ice[END_REF] (modified by Suquet et al [START_REF] Suquet | Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF]) for ice. Then, after recalling the notion of effective potentials in section 3, the notion of plastic modes on which the NTFA model is based, is exposed in section 4. The reduced variables are the weights of the local fields of internal variables on these modes. Then the hard problem of deriving the evolution equation for the reduced variables is addressed in section 5. The tangent-second-order (TSO) expansion is recalled and applied to polycrystals. Finally, two study cases for ice and copper are presented in section 6. The capability of the NTFA-TSO model to reproduce the full-field simulations at a much lower cost is shown. The computational time is reduced by two orders of magnitude with respect to full-field simulations for one of the models.

2 Constitutive relations for single crystals

General formulation

The constitutive relations used in the present study at the single crystal level are formulated in the general framework described by Rice [START_REF] Rice | Inelastic constitutive relations for solids: an internal-variable theory and its application to metal plasticity[END_REF] and Halphen and Nguyen [START_REF] Halphen | Sur les matériaux standard généralisés[END_REF]. The internal rearrangements (in the present case the slips along specific directions) taking place in the material under the application of a load are described by a finite set of internal variables denoted collectively as α. These internal variables, together with the strain ε, are the state variables of the material at a given material point (attention is restricted here to infinitesimal strains and changes in temperature are not taken into account). The free-energy density w(ε, α) of the material is a function of the state variables and the Cauchy stress σ is the derivative of the free-energy with respect to the strain ε, the internal variables being fixed (usual thermo-elastic relations),

σ = ∂w ∂ε (ε, α). (1)
The thermodynamic forces (also called driving forces) associated with the internal variables are defined as

A = - ∂w ∂α (ε, α), (2)
and it is assumed that the evolution of the internal variables is driven by these forces

α = F(A). (3)
When the constitutive functions F exhibit additional symmetries (Onsager's reciprocity relations),

∂ αi ∂A j = ∂ αj ∂A i , i.e. ∂F i ∂A j (A) = ∂F j ∂A i (A), (4)
there exists a potential ψ(A) (called the force potential), such that

F(A) = ∂ψ ∂A (A), i.e. α = ∂ψ ∂A (A).
When both w and ψ are convex functions of their arguments, the corresponding materials are called generalized standard materials ( [START_REF] Halphen | Sur les matériaux standard généralisés[END_REF][START_REF] Germain | Continuum Thermodynamics[END_REF]). Introducing the convex dual ϕ of ψ (ϕ is the dissipation potential), the relation ( 2) can be inverted into

A = ∂ϕ ∂ α ( α). (5)
Upon elimination of A between ( 2) and ( 5), the constitutive relations for a generalized standard material (GSM) take the compact form,

σ = ∂w ∂ε (ε, α), ∂w ∂α (ε, α) + ∂ϕ ∂ α ( α) = 0. (6)
When Onsager symmetry relations (4) are not satisfied, the material is said to be non-standard.

Assuming that the function F is invertible with inverse G, the constitutive equations are

σ = ∂w ∂ε (ε, α), ∂w ∂α (ε, α) + G( α) = 0. (7)

Two crystal plasticity models

In crystal plasticity the viscoplastic strain results from slips on different slip systems,

ε = ε e + ε v , ε e = M e : σ, ε v = S s=1 γ s m s , m s = b s ⊗ s n s , ( 8
)
where S is the total number of slip systems, m s is the Schmid tensor of system s with slip plane normal n s and slip direction b s (orthogonal to n s ), the symbol : denotes contraction over 2 indices, and ⊗ s indicates the symmetrized dyadic product

(b ⊗ s n) ij = 1/2(b i n j + b j n i ).
It remains to specify how the slips γ(s) evolve. These evolution laws are phenomenological and model dependent. We will review two crystal plasticity models which can be formulated in the form of ( 1), ( 8), (3).

The model of Méric-Cailletaud [27]

Initially proposed for cubic metals, it extends to crystal plasticity the well-known phenomenological model of Armstrong and Frederick [START_REF] Armstrong | A mathematical representation of the multiaxial Bauschinger effect[END_REF] accounting for both kinematic and isotropic hardening. The slip-rates on the individual systems are given by

γs = γ0,s (|τ s -x s | -r s ) + K s ns sgn (τ s -x s ) , (9)
where τ s is the resolved shear stress on system s

τ s = σ : m s ,
x s is a back-stress which evolves with the plastic activity according to

ẋs = c s γs -d s x s | γs | , (10)
and r s is a threshold associated with the isotropic hardening of system s and is a function of the (generalized) cumulated plastic slips p s on all systems through the relation

r s = r ini,s + Qb S s =1 h s,s p s , where ṗs = (1 -bp s ) | γs | . ( 11
)
Note that p s (denoted ρ s in [START_REF] Méric | Single crystal modeling for structural calculations. part 2: Finite element implementation[END_REF]) coincides with the cumulated plastic slip when b = 0. When b = 0, p s can be still expressed in terms of the cumulated plastic slip by integration of the differential equation for p s

p s (t) = 1 b (1 -e -b t 0 | γs(θ)|dθ ).

Ice I h

The crystal plasticity model for ice proposed by Castelnau et al [START_REF] Castelnau | Elastoviscoplastic micromechanical modeling of the transient creep of ice[END_REF] and modified in Suquet et al [START_REF] Suquet | Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF] is similar to the model of Méric-Cailletaud [START_REF] Méric | Single crystal modeling for structural calculations. part 2: Finite element implementation[END_REF], except that isotropic hardening is modelled by a multiplicative term rather than with an additive term. More specifically the slip-rates on the individual systems are given by

γs = γ0,s |τ s -x s | r s ns sgn (τ s -x s ) , (12)
where the back-stress evolves with the plastic activity according to an Armstrong-Frederick type law including static recovery :

ẋs = c s γs -d s x s | γs | -e s x s . (13)
The critical resolved shear stress associated with isotropic hardening satisfies the differential equation

ṙs = (r sta,s -r s ) ṗs , ṗs = S s =1 h s,s | γs | , (14)
where the hardening matrix h s,s is symmetric (see table 4). This differential equation can be integrated into

r s (p s ) = r sta,s + (r ini,s -r sta,s )e -ps , ṗs = S s =1 h s,s | γs | . ( 15
)

Internal variables, thermodynamic forces, free-energy

In order to recognize that the constitutive relations of sections 2.2.1 and 2.2.2 can be written in the canonical form (1) (3), the first step is to define the state variables of the system. The internal variables (in addition to the strain ε) are the viscoplastic strain1 ε v , a set of variables β s | s=1,...S analogous but different from the plastic slips to describe kinematic hardening and the (generalized) cumulated plastic slip p s to describe isotropic hardening,

α = (ε v , β, p) , β = (β s )| s=1,...,S , p = (p s )| s=1,...,S . (16)
In most crystal plasticity models, the free-energy function is the sum of three terms, corresponding respectively to the elastic (recoverable) energy, the energy stored in the kinematic hardening of the material and the energy stored in the isotropic hardening,

w(x, ε, α) = 1 2 (ε -ε v ) : L(x) : (ε -ε v ) + 1 2 S s=1 c s β 2 s + w sto (p), (17)
where the stored energy w sto (p) depends on the constitutive model under consideration. In the Méric-Cailletaud model the stored energy takes the form

w sto (p) = S s=1 r ini,s p s + 1 2 Qb S s=1 S s =1 h s,s p s p s . (18)
In the case of ice, the stored energy reads as

w sto (p) = S s=1 w sto,s (p s ), w sto,s (p s ) = ps 0 r s (q) dq. (19)
The thermodynamic forces associated with the state variables ( 16) are A = (A v , A β , A p ), with

A v = L : (ε -ε v ) = σ, A βs = -c s β s = -x s , A ps = -r s = - ∂w sto ∂p s (p). (20)
With these definitions and notations it can be recognized that the constitutive relations of sections 2.2.1 and 2.2.2 can be formulated as a differential equation in the form (3). The detailed derivation is given in appendix A.

3 Polycrystals and effective potentials

Notations

A polycrystal is regarded as an aggregate of a large number N of identical single crystals with different orientations and with a size which is small compared to that of the aggregate. Let V denote the volume occupied by a representative volume element (r.v.e.) of the polycrystal and V (g) , g = 1, .., N denote the domains occupied by the N grains. The characteristic functions of the domains V (g) are denoted by χ (g) and spatial averages over V and V (g) are denoted by . and . (g) respectively. Then, for example,

ε = 1 |V | V ε(x) dx = N g=1 c (g) ε (g) , ε (g) = 1 |V (g) | V (g) ε(x) dx,
where the c (g) = χ (g) are the volume fractions of the grains. The compact notation

ε (g) = ε (g) , σ (g) = σ (g) , (21)
will also be used to denote the average strain and stress in each grain. All grains are assumed to be homogeneous, perfectly bonded at grain boundaries and their constitutive relations and potentials are deduced from that of the single crystal by a rotation Q (g) . For instance the slip systems m (g) s and elastic moduli L (g) in grain g are obtained from the slip systems m s and stiffness L of the single crystal by the relations

m (g) s = Q (g) T .m s .Q (g) , L (g) = Q (g) T .Q (g) T .L.Q (g) .Q (g) ,
where Q (g) T denotes the transpose of Q (g) . The potentials w and ϕ (or operator G) depend on the position x inside V , in such a way that

w(x, ε, α) = N g=1 χ (g) (x) w (g) (ε, α), ϕ(x, α) = N g=1 χ (g) (x)ϕ (g) ( α), G(x, α) = N g=1 χ (g) (x) G (g) ( α).            (22)

Local problem

The r.v.e. V is subjected to a path of macroscopic strain ε(t) and periodicity conditions2 on ∂V .

The local problem to be solved to determine the local fields σ(x, t), ε(x, t) and α(x, t) consists of the generalized thermoelastic problem [START_REF] Kattan | Overall damage and elastoplastic deformation in fibrous metal matrix composites[END_REF], in which the field of internal variables α(x) is fixed, coupled with the differential equation ( 24) at every point x in the volume element,

σ(x, t) = ∂w ∂ε (x, ε(x, t), α(x, t)), div σ(x, t) = 0 in V × [0, T ], ε(x, t) = ε(t) + 1 2 (∇u * (x, t) + ∇u * T (x, t)) in V × [0, T ],
u * periodic on ∂V, σ.n anti-periodic on ∂V,

           (23)
and

∂w ∂α (x, ε(x, t), α(x, t)) + ∂ϕ ∂ α (x, α(x, t)) = 0 in V × [0, T ]. (24)
The second line in [START_REF] Kattan | Overall damage and elastoplastic deformation in fibrous metal matrix composites[END_REF] imply in particular that the path of macroscopic strain is prescribed3 ε(x, t) = ε(t). For non-standard materials the differential equation ( 24) is to be replaced by [START_REF] Duval | Rate controlling processes in the creep of polycrystalline ice[END_REF] at every material point in V .

It is assumed that when the field α is prescribed, the periodic (thermo-elastic) boundary value problem [START_REF] Kattan | Overall damage and elastoplastic deformation in fibrous metal matrix composites[END_REF] has a unique solution, the local strain field ε.

The homogenized (or effective) response of the composite along this path {ε(t), t ∈ [0, T ]} is the history of average stress {σ(t), t ∈ [0, T ]} where σ(t) = σ(x, t) .

Effective potentials

When the constitutive relations of the single crystal derive from two potentials it can be shown that the homogenized constitutive relations for the polycrystal also derive from two potentials, which are the averages over the r.v.e. V of the local potentials. The state variables at the macroscopic level consist of the average strain ε and of an infinite number of internal variables which are the values of the field of internal variables at every microscopic point x ∈ V (see Rice [START_REF] Rice | On the structure of stress-strain relations for time-dependent plastic deformation in metals[END_REF], Suquet [START_REF] Suquet | Elements of Homogenization for Inelastic Solid Mechanics[END_REF] and appendix B),

α = {α(x)} x∈V . ( 25
)
and the effective potentials are defined as

w (ε, α) = w(ε, α) = Inf ε∈K(ε) w(ε, α) , ϕ( ˙ α) = ϕ( α) , (26)
where

K(ε) = {ε = ε + 1 2 (∇u * + ∇u * T ), u * periodic on ∂V }.
Admittedly, the need for an infinite number of internal variables at the macroscopic level is in apparent contradiction with the aim homogenization, which is to express the macroscopic constitutive relations in terms of purely macroscopic variables. The existence of two macroscopic potentials is essentially a qualitative result, which can be turned into a much more practical result when the local fields of internal variables can be expressed by means of a finite number of variables. This is exactly what the NTFA does (next section).

NTFA decomposition 4.1 Reduced variables

One of the earliest model-reduction technique in micromechanics is the Transformation Field Analysis (TFA) of Dvorak [START_REF] Dvorak | Transformation field analysis of inelastic composite materials[END_REF], where the plastic strain field (or more generally the transformation strain field) is considered to be uniform in the individual phases, or in specific sub-domains within the phases. As a result of this assumption, the whole plastic strain field is characterized by its values over these domains which can be taken as a finite number of internal variables. However this is often far from reality and in addition, as noticed by Suquet [START_REF] Suquet | Effective properties of nonlinear composites[END_REF], the assumption of piecewise uniformity for the plastic strain field leads to an over-stiffening of the effective response of the volume element. By contrast, the plastic strain field is often seen, in experiments and in numerical simulations, to be organized along specific patterns. This observation has led the authors ( [START_REF] Michel | Constitutive relations involving internal variables based on a micromechanical analysis[END_REF][START_REF] Michel | Nonuniform Transformation Field Analysis[END_REF]) to relax the condition of uniformity by assuming only that the local plastic strain field can be decomposed on a basis of a few shape functions (the observed patterns), called plastic modes. The corresponding decomposition in the present context of internal variables (instead of plastic strains) reads as

α(x, t) = M k=1 ξ (k) (t) µ (k) (x), (27)
where the fields µ (k) (x) are the modes and the ξ (k) are the reduced variables. The modes have the same tensorial character as the internal variables α. To avoid a possible indeterminacy in the definition of the reduced variables ξ (k) , it is further assumed that the modes µ (k) are linearly independent fields.

The modes can be inferred from experimental observations, analytical results or from numerical simulations. Regarding computational methods, it suffices to say that several techniques, such as the Proper Orthogonal Decomposition (POD) (also known under different names, Singular Value Decomposition, or Karhunen-Loève decomposition), are available in the literature to select the modes. The reader is referred to [START_REF] Chinesta | PGD-Based Modeling of Materials, Structures and Processes[END_REF] (among others) for general references and to [START_REF] Michel | Nonuniform Transformation Field Analysis[END_REF][START_REF] Michel | Computational analysis of nonlinear composite structures using the Nonuniform Transformation Field Analysis[END_REF][START_REF] Michel | Nonuniform Transformation Field Analysis: a reduced model for multiscale nonlinear problems in Solid Mechanics[END_REF][START_REF] Largenton | Extension of the nonuniform transformation field analysis to linear viscoelastic composites in the presence of aging and swelling[END_REF] and to section 6.2.1 for specific examples of the snapshot POD applied to the NTFA. How the modes are chosen or generated is not the objective of the present study, and the modes are assumed to be known in the remainder of this section.

Reduced kinetics: a hard problem

Identifying relevant modes for the decomposition ( 27) is probably the easiest part of the model. Of higher complexity is the derivation of the reduced (or coarse) "kinetic" equations, i.e. the evolution equations for the reduced variables ξ. When the decomposition ( 27) is incorporated in the definition [START_REF] Lucia | Reduced-order modeling: new approaches for computational physics[END_REF] of the effective potentials it is readily seen that the general structure of the constitutive relations deriving from two potentials is preserved, but now for a finite number of internal variables ξ, with the following potentials

w(ε, ξ) = Inf ε∈K(ε) w(ε, α) , ϕ( ξ) = ϕ( α) , with α(x, t) = M k=1 ξ (k) (t) µ (k) (x). (28)
The differential equations [START_REF] Largenton | Extension of the nonuniform transformation field analysis to linear viscoelastic composites in the presence of aging and swelling[END_REF] are reduced to M ordinary differential equations for the finite set of reduced variables ξ,

∂ w ∂ξ (k) (ε, ξ) + ∂ ϕ ∂ ξ(k) ( ξ) = 0, k = 1, ...M. (29)
In conclusion, the NTFA decomposition results in a model with a finite number of internal variables for which M ordinary differential equation can be derived, provided the two effective potentials w and ϕ can be computed explicitly.

However two difficulties arise. The first difficulty stems from the fact that the differential equation ( 29) is a doubly nonlinear equation. It would be much easier to integrate if it were in the form ξ = f(ξ,ε(t)). Second, the model is a purely macroscopic one (the ultimate goal of homogenization) provided that the potentials w and ϕ can be pre-computed explicitly as functions of (ε, ξ) and ξ respectively. As will be seen in section (4.3.2), w can be determined explicitly at the expense of solving "off-line" (once for all) a few linear thermo-elastic unit-cell problems. This decoupling between the micro and the macro scales cannot be achieved for ϕ which must be computed "on-line". Indeed, each time that ϕ is called in the course of the integration of ( 29) one has to -compute the fields α(x) = M =1 ξ( ) µ ( ) (x) by calling the modes µ (k) which must be stored (expensive in terms of memory), -compute ϕ(x, α(x)) at every local point x ((expensive in terms of CPU), -average over the r.v.e.

The high cost of this "on-line" computation diminishes significantly the performance and the interest of model-reduction.

To overcome both difficulties Michel and Suquet [START_REF] Michel | A model-reduction approach in micromechanics of materials preserving the variational structure of constitutive relations[END_REF] have proposed a rather systematic approach recalled in section 5. It was initially applied to two-phase composites and is extended to crystal plasticity in the present study.

Application of the NTFA to crystal plasticity 4.3.1 Reduced variables

The NTFA decomposition is applied to the choice ( 16) of internal variables, common to both models of crystalline plasticity. The fields of internal variables are decomposed as

ε v (x, t) = Mv k=1 ξ (k) v (t) µ (k) v (x), β s (x, t) = M βs =1 ξ ( ) βs (t) µ ( ) βs (x), p s (x, t) = N g=1 p (g) s (t) χ (g) (x), s = 1, ..., S.            (30)
Note that the isotropic hardening variables p s are assumed to be uniform in each grain, whereas the viscoplastic strain and the kinematic hardening variables β s are decomposed on more complex (but less numerous) "shape functions". The modes µ v are second-order tensorial fields whereas the modes µ βs are scalar fields. M v and M βs denote the number of modes for the viscoplastic strain and the the kinematic hardening variables respectively. The reduced variables ξ are defined as

ξ = (ξ v , ξ β , p), ξ v = (ξ (k) v )| k=1,..,Mv , ξ βs = (ξ ( ) βs )| =1,..,M βs , p s = (p (g) s )| g=1,.
.,N , s = 1, ..., S.

Effective energy and reduced forces

When the free-energy takes the form [START_REF] Germain | Continuum Thermodynamics[END_REF], the thermoelastic problem equivalent to the minimization in ( 26) is linear and reads as

σ = L(x) : (ε -ε v ), div σ = 0, ε(t) = ε(t), periodicity conditions on ∂V. ( 31
)
After due account of the decomposition (30) for ε v , the solution of this problem can be expressed by the superposition theorem as

ε(x) = A(x) : ε + Mv =1 (D * µ ( ) v )(x)ξ ( ) v , (32)
where A(x) is the strain concentration tensor, expressing the local strain field caused by an average strain ε when the eigenstrain ε v vanishes, D(x, x ) is the nonlocal Green operator expressing the strain at point x resulting from a unit eigenstrain at point x , when the average strain vanishes and * denotes the convolution product in space. The stress field corresponding to (32) can be expressed as

σ(x, ε, ξ v ) = L(x) : A(x) : ε+ Mv =1 ρ ( ) (x)ξ ( ) v , ρ ( ) (x) = L(x) : (D * µ ( ) v ) -µ ( ) v (x), (33)
and the fields of thermodynamic forces [START_REF] Holmes | Structures, Dynamical Systems and Symmetry[END_REF] are functions of (ε, ξ) given by

A v (x, ε, ξ v ) = σ(x, ε, ξ) = L(x) : A(x) : ε + Mv k=1 ρ (k) (x)ξ (k) v , A βs (x, ξ βs ) = -c s M βs =1 µ ( ) βs (x)ξ ( ) βs A ps (x, p s ) = - N g=1 χ (g) (x)r (g) s , r (g) s = r s (p (g) ).                        (34)
The effective free-energy w(ε, ξ) = w(ε, α) reads as

w(ε, ξ) = 1 2 ε : L : ε -ε : Mv k=1 ξ (k) v a (k) + 1 2 Mv k,k =1 (L (kk ) -D (kk ) )ξ (k) v ξ (k ) v + 1 2 S s=1 M βs , =1 H ( ) s ξ ( ) βs ξ ( ) βs + N g=1 c (g) w (g) sto (p (g) ), (35)
where

L = A T : L : A , a (k) = µ (k) v : L : A , D (kk ) = µ (k) v : L : (D * µ (k ) v ) , L (kk ) = µ (k) v : L : µ (k ) v , H ( ) s = c s µ ( ) βs µ ( ) βs .    (36)
The reduced forces associated with the reduced variables are defined as

a v = - ∂ w ∂ξ (ε, ξ), a βs = - ∂ w ∂ξ βs (ε, ξ), a ps = - ∂ w ∂p s (ε, ξ). ( 37
)
Taking into account the expression (35) of the effective energy, these forces read as

a (k) v = a (k) : ε + Mv =1 (D -L) (k ) ξ ( ) v , a ( ) βs = - M βs =1 H ( ) s ξ ( ) βs , a (g) ps = -c (g) r s (p (g) ). ( 38
)
It remains to specify the evolution equations for the reduced variables ξ or equivalently for the reduced forces a. This "reduced kinetics" (also referred to as the "reduced dynamics" in [START_REF] Michel | A model-reduction approach in micromechanics of materials preserving the variational structure of constitutive relations[END_REF]) is derived in the next section.

Reduced kinetics: the hard problem!

The first difficulty, which was already identified in section 4.2, is the double nonlinearity of the differential equation ( 29) which we would prefer to write as ξ = f(ξ,ε(t)). This can be done by inverting ∂ ϕ ∂ ξ or equivalently by computing the convex dual ϕ * of ϕ. The differential equation (29) could then be rewritten as,

a = - ∂ w ∂ξ (ε, ξ), ξ = ∂ ϕ * ∂a (a), or equivalently ξ = ∂ ϕ * ∂a - ∂ w ∂ξ (ε, ξ) . ( 39
)
The differential equation [START_REF] Ryckelynck | Multi-level a priori Hyper-reduction of mechanical models involving internal variables[END_REF] for ξ is now in usual form. However computing ϕ * is a formidable task in general, since ϕ is not even known explicitly.

Hybrid model (Fritzen and Leuschner [14])

To avoid the computation of ϕ * , Fritzen and Leuschner [START_REF] Fritzen | Reduced basis hybrid computational homogenization based on a mixed incremental formulation[END_REF] make use of a potential which is

closely related to ϕ * , ψ(a(ε, ξ)) = ψ(A(ε, ξ)) , (40)
where ψ = ϕ * and A(ε, ξ) is given by [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] in the present context. Michel and Suquet [START_REF] Michel | A model-reduction approach in micromechanics of materials preserving the variational structure of constitutive relations[END_REF] have shown that the potential ψ is not the dual of ϕ in full rigour, but is an accurate approximation of it when enough modes are taken in the decomposition [START_REF] Méric | Single crystal modeling for structural calculations. part 2: Finite element implementation[END_REF]. The "reduced kinetics" (ordinary differential equations for the reduced variables ξ) then becomes

a = - ∂ w ∂ξ (ε, ξ), ξ = ∂ ψ ∂a (a). (41)
In practice, it is easier to re-write (41) as a differential equation for a. This can be done by noting that a is a function of (ε, ξ) and that its time derivative reads as

ȧ = ∂a ∂ε : ε + ∂a ∂ξ . ξ.
Using the evolution equation ( 41) for ξ and the chain rule for computing derivatives, one obtains

ȧ = ∂a ∂ε : ε + ∂ ψ ∂ξ (ε, ξ(ε, a)). (42)
Equations ( 41) and ( 42) can be used for standard materials governed by two potentials. For nonstandard materials for which the evolution of the internal variables is ruled by (3), the evolution equations are obtained by first noting that, by the definition (40) of ψ ,

∂ ψ ∂a (ε, ξ(ε, a)) = ∂ψ ∂A (A). ∂A ∂a , ∂ ψ ∂ξ (ε, ξ(ε, a)) = ∂ψ ∂A (A). ∂A ∂ξ .
Then (41) 2 and (42) can be alternatively written as

ξ = ∂ψ ∂A (A). ∂A ∂a or ȧ = ∂a ∂ε : ε + ∂ψ ∂A (A). ∂A ∂ξ .
Finally, when the crystal plasticity model is non-standard (which is the case of the two models presented in sections 2.2.1 and 2.2.2), F(A) is substituted for ∂ψ ∂A (A) in these relations and the following differential equations are obtained,

ξ = F(A). ∂A ∂a or ȧ = ∂a ∂ε : ε + F(A). ∂A ∂ξ . ( 43
)
One could use ( 42) and ( 43) as the reduced kinetics and consider that the NTFA model is now complete. Unfortunately, as Michel and Suquet [START_REF] Michel | A model-reduction approach in micromechanics of materials preserving the variational structure of constitutive relations[END_REF] noted, the potential (40) requires "on-line computations" of the local fields A(x) and of ψ(x, A(x)), similar to the "on-line computations" required to evaluate ϕ * . Therefore this potential alone does not solve the second difficulty identified in section 4.2 which is responsible for the deterioration of the speed-up that one could expect from the model reduction. So an additional step is needed to arrive at a purely homogenized model (requiring only on "off-line" computations).

TSO approximation

To achieve this additional step, the authors [START_REF] Michel | A model-reduction approach in micromechanics of materials preserving the variational structure of constitutive relations[END_REF] have proposed to use a technique previously developed in homogenization of nonlinear materials governed by a single potential. A Taylor expansion at second-order ψ T SO is substituted for the exact potential ψ in [START_REF] Sirovich | Turbulence and the dynamics of coherent structures[END_REF]. In the present context of crystal plasticity where the potential ψ is the sum of potentials on slip systems (see (60)), the second-order expansion of ψ results from the expansion of each individual potential ψ s in grain (g)

ψ s (A (g) ) + ∂ψ s ∂A (A (g) ). A -A (g) + 1 2 A -A (g) . ∂ 2 ψ s ∂A 2 (A (g) ). A -A (g) , (44)
where A (g) = A (g) is the average of A over grain (g) . Upon averaging over each grain and then over the whole volume, the following approximate effective potential is found

ψ(ε, ξ) ψ T SO (ε, ξ) = N g=1 c (g) S s=1 ψ s (A (g) ) + 1 2 ∂ 2 ψ s ∂A 2 (A (g) ) : C (g) (A) , (45)
where C (g) (A) denotes the quadratic fluctuations of A over grain (g):

C (g) (A) = (A -A (g) ) ⊗ (A -A (g) ) (g) .
Therefore for standard materials the proposed approximation consists in substituting ψ T SO for ψ in [START_REF] Suquet | Effective properties of nonlinear composites[END_REF].

For non-standard materials, the approximation consists in substituting the right-hand side terms in [START_REF] Suquet | Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF] with their second order expansion with respect to all the fields which have intragranular fluctuations:

ξ = F(A). ∂A ∂a T SO or ȧ = ∂a ( ) ∂ε : ε + F(A). ∂A ∂ξ T SO . ( 46
)
How to get these TSO expansions for each specific model is detailed in appendix C.

Application to crystal plasticity

Since in both models the evolution equation for ε v derives from a potential (see appendix A)

εv = ∂ψ ∂A v (A),
the evolution equation for the corresponding force a v is given by [START_REF] Suquet | Effective properties of nonlinear composites[END_REF],

ȧ(k) v = ∂a (k) v ∂ε : ε + ∂ ψ T SO ∂ξ (k) v (ε, ξ(ε, a)). (47)
The potential ψ has a common form in both models and can be decomposed as a sum of individual potentials on the slip systems (see appendix A)

ψ s (A) = f s (τ s -x s , r s ), τ s = A v : m s , x s = -A βs , r s = -A ps , (48)
where f s depends on the model under consideration, Méric-Cailletaud:

f s (τ, r) = K s γ0,s n s + 1 (|τ | -r) + K s ns+1 , Ice: f s (τ, r) = r γ0,s n s + 1 |τ | r ns+1 .          (49)
The approximate TSO potential in grain g reads as

ψ T SO (A(ε, ξ)) (g) = S s=1 f s (τ (g) s -x (g) s , r (g) s ) + 1 2 ∂ 2 f s ∂τ 2 (τ (g) s -x (g) s , r (g) s )C (g) (τ s -x s ) , (50)
where

τ s (x) = A v (x) : m s (x), τ (g) s = τ s (g) , C (g) (τ s ) = (τ s -τ (g) s ) 2 (g) , x (g) s = x s (g) = c s M βs =1 ξ ( ) βs µ ( ,g) βs , µ ( ,g) βs = µ ( ) βs (g) , r (g) s = r s (p (g) ).          (51)
In deriving (50) use has been made of the fact that r s has no fluctuation in grain r (therefore only the fluctuations of τ s and x s are taken into consideration). Then the evolution equation (47) for the reduced force a v (corresponding to the state variable ε v which is standard) reads as

ȧ(k) v = a (k) : ε + N g=1 c (g) S s=1 ∂f s ∂τ (τ (g) s -x (g) s , r (g) s ) ∂τ (g) s ∂ξ (k) v + 1 2 ∂ 2 f s ∂τ 2 (τ (g) s -x (g) s , r (g) s ) ∂C (g) (τ s -x s ) ∂ξ (k) v + 1 2
∂ 3 f s ∂τ 3 (τ (g) s -x (g) s , r (g) s )C (g) (τ s -x s ) ∂τ (g) s ∂ξ (k) v . ( 52
)
By contrast, the variables β s and p s are non-standard (their evolution is not governed by the derivative of a potential) and the evolution equations for the corresponding reduced forces are, in general, given by [START_REF] Weertman | Creep of ice[END_REF]. In the present context, the derivatives of a βs and a ps with respect to ε vanish. Moreover the only force depending on ξ βs is A βs and the only forces depending on p s are the A p s ( and even only A ps in the case of Méric-Cailletaud). The general equations [START_REF] Weertman | Creep of ice[END_REF] then reduce to

ȧ( ) βs = F βs (A) ∂A βs ∂ξ ( ) βs T SO , ȧ(g) ps = S s =1 F p s (A) ∂A p s ∂p (g) s T SO . ( 53
)
The corresponding equations for the two crystal plasticity models are given in details in appendix C.

Applications

In the present section, examples of the implementation of the NTFA approach for the two constitutive relations of sections 2.2.1 and 2.2.2 are given. Some of the methodological tools are common to both examples and are exposed before proceeding to more specific applications.

6.1 Full-Field simulations

Computational method

The full-field simulations can be performed by any appropriate method. The computational method used in the present study is based on fast Fourier transforms which is a fast and flexible alternative to Finite Elements for volume elements subjected to periodic boundary conditions ( [START_REF] Michel | Effective properties of composite materials with periodic microstructure: a computational approach[END_REF]). It was initially proposed for composites ( [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]) and then extended to polycrystals by Lebensohn [START_REF] Lebensohn | N-site modelling of a 3d viscoplastic polycrystal using Fast Fourier Transforms[END_REF]. All differential equations appearing in the constitutive relations were integrated using a 4th-order Runge-Kutta scheme with controlled time-stepping ( [START_REF] Michel | Nonuniform Transformation Field Analysis: a reduced model for multiscale nonlinear problems in Solid Mechanics[END_REF]).

Representative volume element

Full-field simulations are performed on artificial polycrystalline microstructures generated by Voronoi tesselation (under the constraint of periodicity). The centers of the grains are generated by a Poisson's process whereas the orientation of each grain is chosen randomly. In order to assess the sensitivity of the effective response to the microstructure, 8 different realizations with 500 grains were generated. A typical realization is shown in figure 2 the single crystal level are those of section 2.2.1 (Méric-Cailletaud) with material parameters for copper taken from Gérard et al [START_REF] Gérard | Comparison of experimental results and finite element simulation of strain localization scheme under cyclic loading[END_REF] and recalled in appendix D. The loading path is specified in (57). The responses of the eight representations of the volume element are compared in figure 1(a). As can be seen there is very little discrepancy between the different configurations. Then the local stress fields were compared. The norm of the stress field at the end of the cycle, defined as

||σ|| = √ σ ij σ ij , (54)
is plotted in figure 1(b). Again the dispersion is small and the distribution of the stress is consistent between the different configurations. Therefore it can be concluded that a volume element containing 500 grains is sufficient to represent accurately the overall response of the polycrystal.

Spatial discretization

Similarly a parametric study has been conducted to assess the sensitivity of the full-field simulations to the spatial resolution of the images used in the FFT computations. The same Voronoi tesselation of the volume element has been discretized into 75 3 , 125 3 and 175 3 voxels. The predictions of the effective response were found to be almost on top of each other and their distribution of the local fields did not show any significant deviation (for brevity these results are not shown here). Therefore all simulations presented in the remainder of this study were performed with the configuration of figure 2(a) at an intermediate spatial resolution of 125 3 voxels.

6.2 Implementation of the NTFA

Snapshot POD

The modes corresponding to all internal variables ε v , γ s and β s are generated by the Karhunen-Loève method (or POD) method which is now briefly recalled. First full-field simulations are performed along loading paths which are chosen by the user (these paths are sometimes called the learning paths). The choice of these paths depends largely on the application that the user has in mind. They are very similar to the tests an experimentalist would perform to identify the relevant material parameters of a phenomenological constitutive relation. Along a given loading path, T snapshots of the appropriate internal variable fields computed by the full-field method are stored at time steps. These snapshots θ (k) are tensorial or scalar fields (depending on the internal variable under consideration). Then, following the classical POD procedure, the correlation matrix of these snapshots is formed

g (k ) = θ (k) (x) : θ ( ) (x) ,
and its eigenvalues λ (k) and eigenvectors v (k) are computed. The actual modes are expressed as a combination of the eigenvectors as

µ (k) (x) = T =1 v (k) θ ( ) (x), k = 1, ..., T,
where T denotes the number of snapshots. Note that the modes, being the eigenvectors of the symmetric matrix g, are orthogonal to each other. Another well-known feature of the POD is that the quantity of relevant information (its correlation with the set of snapshots) contained in an eigenvector v (k) is expressed by the magnitude of the corresponding eigenvalue λ (k) . This property can be used to truncate the set of modes, retaining the most relevant ones. Ordering the eigenvalues in decreasing order, the M modes corresponding to the largest eigenvalues can be selected by applying a threshold criterion

M k=1 λ (k) / T k=1 λ (k) ≥ α. (55)
In the present study the threshold was set to α = 0.9999 = 1 -10 -4 .

Two variants of the modes for the viscoplastic strain field

The NTFA method has been implemented with modes for the viscoplastic strain ε v and for the kinematic hardening variables β s (s = 1, ..., S) (the isotropic hardening variables p s are uniform per grain). Two different strategies can be adopted for generating the modes for ε v .

-In the first strategy (denoted NTFA-TSO 1 in the following), the modes µ

(k)
v are extracted by applying the POD to snapshots of the viscoplastic strain ε v . In the examples investigated in the course of the present study, 3 to 5 modes were sufficient to meet the criterion (55) with a threshold α = 0.9999 = 1 -10 -4 . The number of modes obviously depends on many different parameters, such as the heterogeneity of the plastic strain field or the complexity of the learning loading path.

-In the second strategy (denoted NTFA-TSO 2 in the following), the modes for the individual slips γ s are first generated,

γ s (x, t) = Ms =1 ξ ( ) s (t) µ ( ) s (x), s = 1, ...S,
where M s denotes the number of modes for system s. The modes for the viscoplastic strain field ε v can then be formed as

µ (k) v = µ ( ) s m s , k = + s <s M s . (56)
In the two examples investigated here, 2 to 5 modes per slip system were necessary to achieve the criterion (55).

As will be seen later, each approach has its own merits:

-The formulation based on the decomposition of the γ s 's allows for the reconstruction of the individual slips, which is not possible with the first formulation. This is the best option for reconstructing the activity on the different groups of systems. However in this second approach, the number of modes resulting from the decomposition (56), M v = S s=1 M s , can be significantly larger that the number of modes required by the first formulation.

-For this reason, the first formulation based on the direct extraction of the modes from snapshots of ε v involves less internal variables and is significantly faster (see sections 6.3.5 and 6.4.3).

The modes µ ( ) βs are identical in both NTFA approaches. The isotropic variables p s are assumed to be piecewise uniform in both approaches.

Example 1: calibration of material parameters for polycrystalline ice 6.3.1 Reduced-order modelling and identification

Constitutive models in crystalline plasticity often involves a large number of parameters. These parameters can be extracted from experimental data by several algorithms (see an example of the use of the Levenberg-Marquardt algorithm in [START_REF] Herrera-Solaz | An inverse optimization strategy to determine single crystal mechanical behavior from polycrystal tests: Application to AZ31 Mg alloy[END_REF]). The different iterations of the identification algorithm require multiple resolutions of the direct problem and this can be excessively timeconsuming. A robust reduced-order model can save a significant part of the CPU time required to solve the direct problem. To establish its robustness, it must be proved that the full-field simulations can be properly reproduced by the NTFA model. Then, it should also be proved that the sensitivity of the response to the variation of one (or several) material parameters is properly captured by the model. Of course the NTFA modes have to be identified once for all.

The calibration of the material parameters entering the constitutive relations of section 2.2.2 has been done in [START_REF] Castelnau | Elastoviscoplastic micromechanical modeling of the transient creep of ice[END_REF] for an earlier version of the model and in [START_REF] Suquet | Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF] for the present model, in several steps (just by hand, with out resorting to any identification algorithm such as the Levenberg-Marquardt algorithm).

1. First, data compiled by Duval [START_REF] Duval | Rate controlling processes in the creep of polycrystalline ice[END_REF] can be used to determine the stationary flow stress and the stress sensitivity exponent n (s) for basal slip. Altough there is quite a large spread in these experimental results from different authors, the stress sensitivity exponent for basal slip can be directly identified from these experimental data (numerical values are reported in table 4), whereas the stationary flow stress depends on both the stationary reference stress r (s)

sta and the stationary backstress x (s) .

2. Next, data from Weertman [START_REF] Weertman | Creep of ice[END_REF] can be used for the identification of the transient creep regime of basal systems. Single crystals are subjected to a uniaxial compression test at different strain-rates, with c-axis oriented at 45 • from the loading direction. These tests shed light on the softening of basal slip in the transient regime. The static recovery term e (s) in the constitutive law (13) helps achieving the correct stationary stress at very small strain-rates (since x (k) tends to a constant value c (s) /d (s) at large shear γ (s) if static recovery is not introduced).

3. For prismatic and pyramidal systems, a first estimation of the stress sensitivity exponents n (s) and stationary values r

(s)
sta can be found in [START_REF] Castelnau | Modelling viscoplastic behavior of anisotropic polycrystalline ice with a self-consistent approach[END_REF][START_REF] Castelnau | Elastoviscoplastic micromechanical modeling of the transient creep of ice[END_REF]. These estimations are based on experimental data for secondary creep where the elastic strains are negligible. In this purely viscous regime, the material parameters n (k) and r (r) sta were tuned in [START_REF] Castelnau | Elastoviscoplastic micromechanical modeling of the transient creep of ice[END_REF] to reproduce by a mean-field homogenization technique the overall response (stress sensitivity exponent and flow stress) of untextured polycrystals. [START_REF] Castelnau | Modelling viscoplastic behavior of anisotropic polycrystalline ice with a self-consistent approach[END_REF]. The parameters which are left undetermined by these previous studies (initial reference stresses r ini , hardening parameters h s,s ), were adjusted (trial-and-error) in [START_REF] Suquet | Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF] by comparing full-field simulations with the experimental data compiled by Ashby and Duval [START_REF] Ashby | The creep of polycrystalline ice[END_REF] for the transient creep response of untextured (isotropic) polycrystals.

In the initial study of Suquet et al [START_REF] Suquet | Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF] the calibration of these parameters required numerous full-field simulations on large polycrystalline microstructure (4096 grains) discretized into 256 3 voxels. The simulations were performed with a FFT-based solver ( [START_REF] Suquet | Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF]) and one simulation took 92 hours on a 8 cores Xeon W5580 workstation (CPU clock 3.2GHz). It is immediately seen that the computational effort required for the identification of all parameters was quite substantial and that a significant part of this computational effort could be saved by using a reduced-order model such as the NTFA. A first full-field simulation could be performed with a rather crude estimate of the material parameters. Extracting suitable modes from this first simulation, the NTFA-TSO model could be used to estimate, approximately, the unknown material parameters which could be fine-tuned with full-field simulations. It remains to check that this scenario is plausible. In other words, firstly that the NTFA-TSO model with the right set of parameters can capture correctly the experiments. Secondly that the trends of the full-field simulations are well reproduced when the material parameters are varied. And thirdly that the acceleration in CPU time provided by the reduced order model is significant.

Full-field simulations

We first proceed with the full-field simulations by the FFT. The present implementation of the polycrystalline constitutive relations of section 2.2.2 follows [START_REF] Suquet | Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF]. The configuration of figure 2(a) is discretized into 125 3 voxels of equal size is used. A creep experiment is simulated in which the overall stress is taken the form

σ(t) = σ(t)e 3 ⊗ e 3 ,
where the overall stress σ(t) is linearly increased from 0 to 1 MPa in 10 constant time steps of 0.05s. Then σ(t) is kept constant at 1 MPa up to t = 120000 s with the following time steps:

∆t =     
0.5s when 0.5s < t ≤ 600s, 5s when 600s < t ≤ 6000s, 10s when t > 6000s.

Three full field simulations were performed with different material data. The first data set (denoted data1 below), given in appendix E, is the data set identified in [START_REF] Suquet | Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF] using a different FFT code and a different set of configurations (4096 grains compared to 500 grains in the present study). It is used here to check the representativity of the configurations with 500 grains used in the present study by comparison with the results with 4096 grains. The other two data sets correspond to different values of the latent hardening parameters and have been chosen to highlight the sensitivity of the response to the hardening parameters data 1:

h 3 = 125, h 4 = 0, data 2: h 3 = 0, h 4 = 0, data 3: h 3 = 1100, h 4 = 1200,
(all units are MPa). The results of the full-field simulations with data 1 (solid line) are compared with experimental results (circles) in figure 3. It is seen that the configuration with 500 grains is sufficient to correctly capture the experimental data and can be used for implementation of the NTFA method.

NTFA

The procedure outlined in section 6.2.1 was followed. The POD delivered the following number of modes for each variant of the model, Note that since the isotropic variables p s are assumed to be uniform in each grain, the number of modes for these variables is the total number of grains (multiplied by the number of systems). The three modes generated by the POD for the viscoplastic strain in the formulation NTFA-TSO1 are shown in figure 2. The modes of higher number incorporate more and more details on the local viscoplastic strains. Figure 3: Creep test on polycrystalline ice. Comparison between full-field simulations (black solid line) and the predictions of the two variants of the NTFA. Material parameters identified in [START_REF] Suquet | Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF]. Left: NTFA-TSO 1 (red dashed line). Right: NTFA-TSO 2 (blue dot-dashed line).

Effective response

The predictions of the two NTFA-TSO approaches for the overall response of the polycrystalline aggregate are compared in figures 3 and 4 with full-field simulations. The following observations can be made.

1. The two NTFA models give very close predictions of the overall response of the aggregate. The choice of one or the other model is rather motivated by the local fields that need to be reconstructed (as already mentioned). For parameter identification only, the NTFA-TSO 1 is sufficient.

2. The agreement between the NTFA predictions and the full-field simulations is seen to be slightly better in figure 3 than in figure 4. This is due to the fact that the modes used in all NTFA simulations of figure 4 are those determined with the material parameters data 1. Obviously the agreement would be better if for each set of parameters, the NTFA simulations were performed with modes identified with these particular parameters.

3. Finally, the plots in figure 4 shows that the modes determined with one set of parame- ters can be used with other parameters and still capture correctly the influence of these parameters. In this example, the NTFA models with the material parameters data 2 and data 3 capture accurately the influence of the hardening parameters h 3 and h 4 observed in the full-field simulations, even though all NTFA simulations were performed with modes determined with the set of material parameters data 1. The NTFA can predict that when these hardening parameters are set above their nominal value (data 2 in figure 4) the aggregate is too stiff and the strain-rates drops too rapidly. By contrast, when these parameters are below their nominal value (data 3) the aggregate is too soft and the drop in strain-rate is not fast enough. Therefore the NTFA model can be used in an iterative algorithm of identification of the material parameters, with modes which are determined once for all. It is fair to emphasize that, for good accuracy, the modes should however contain the main features of the local fields of internal variables at convergence.

Speed-up

The CPU times for the different methods (Full-field simulations, NTFA-TSO 1 and NTFA-TSO 2) are compared in Table 1 (Intel Xeon X5687 @ 3.6 GHz). The different speed-ups, measured by the CPU ratios shown in Table 1 can be commented as follows:

1. The acceleration provided by the NTFA varies between 10 2 (NTFA-TSO 2) and 4.10 2 (NTFA TSO 1). In other words, a single full-field simulation which took more than 3 days can be performed in less than one hour and sometimes in a few minutes. This reduction can have a dramatic impact on parameter identification.

2. The second version of the NTFA, where each individual slip is decomposed on a set of modes, is 4 times slower than the first version. This is due to the fact that 30 modes (and therefore 30 reduced variables ξ) are required for the individual slips whereas only 3 modes (and therefore only 3 reduced variables) are required for the plastic strain in the first formulation. However this second version of the method has the merit that the slip activity for each slip system, or each family of slip systems, can be reconstructed by post-processing the macroscopic results (see section 6.3.6).

By comparison with the acceleration observed by the authors ([33]

) for two-phase materials, the acceleration provided by the NTFA for polycrystals is less spectacular (10 2 compared to 10 4 ).This is due to the loop over all grains (corresponding to the summation over r in all the evolution equations) and to the loop over all slip systems (summation over s).

Full-field (FFT) Intel Xeon X5687 @ 3.6 GHz.

Local fields

A significant advantage of the NTFA method, already emphasized in [START_REF] Michel | Nonuniform Transformation Field Analysis[END_REF][START_REF] Michel | Computational analysis of nonlinear composite structures using the Nonuniform Transformation Field Analysis[END_REF][START_REF] Michel | Nonuniform Transformation Field Analysis: a reduced model for multiscale nonlinear problems in Solid Mechanics[END_REF][START_REF] Michel | A model-reduction approach in micromechanics of materials preserving the variational structure of constitutive relations[END_REF] and [START_REF] Largenton | Extension of the nonuniform transformation field analysis to linear viscoelastic composites in the presence of aging and swelling[END_REF] is that the local fields can be easily reconstructed from the knowledge of the reduced state variables (ε, ξ) and from pre-computed fields:

1. The local fields of state variables are, by construction of the NTFA, reconstructed by means of the decomposition (30).

2. The local fields of thermodynamic forces (stress, resolved shear, back-stress on individual slip systems) are reconstructed by means of the relations [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF].

The second version of the NTFA (where modes are determined first for the slips on the individual slip systems) allows for the reconstruction of the activity on each slip system or each family of slip systems. The activity of the 3 families of systems (basal, prismatic, pyramidal) as found from the full-field simulations and as obtained by post-processing the results of the NTFA-TSO2 model are compared in figure 5 at the end of the creep test. The activity at time t is defined here as

s∈F |γ s (t)| ,
where the symbol s ∈ F means that the summation is extended to all systems in a given family F . Note that the color codes are different for each family of systems (the activity on the prismatic and pyramidal systems are less than on the basal systems). This good agreement is confirmed by the probability distribution. It is also worth noting that the formulation NTFA-TSO 2 involving more modes and therefore more internal variables does not significantly improve on the first formulation NTFA-TSO 1.

Speed-up

The CPU times for the different methods (Full-field simulations, NTFA-TSO 1 and NTFA-TSO 2) are compared in table 2. All simulations were performed with the same global time-step ∆t = 1.5 10 -3 s and the same time-integration scheme (4th-order Runge Kutta with controlled time-stepping). The acceleration provided by the first formulation is of the order of 5 10 2 comparable to that observed for ice (section 6.3.5). By contrast the second formulation NTFA-TSO 2 is 10 times slower than the first formulation which is worse than in the case of ice (where it was slowed down by a factor of 4). The difference resides in the larger number of modes used in the decomposition of the viscoplastic strain when the second formulation is chosen (60 modes in copper, compared to 30 modes in ice).

Full-field (FFT) NTFA-TSO

Local stress fields

The norm of the local stress fields at the end of the cycle (t = 15s) predicted by the two NTFA formulations with the localization relation [START_REF] Michel | A model-reduction approach in micromechanics of materials preserving the variational structure of constitutive relations[END_REF] are compared in figure 7 with full-field simulations. The color scale used for the full-field simulations is different from that used for the NTFA predictions as the NTFA tends to underestimate the maximal stress. This underestimation of the maximal stress is also evidenced in figure 6(b) where the probability distribution predicted by the NTFA is narrower, indicating a more homogeneous stress map in the model than in the full-field simulations. The second formulation (NTFA-TSO 2) is in better agreement with the full-field simulations. Even though the maximal stress is underestimated by the model, figure 7 shows that the 'hot-spots' (locations where the stress is maximal) are well predicted.

Closing remarks and perspectives

The present study can be extended in different directions. The first direction is that of multiaxial loadings. At this point, it is important to make a distinction between the learning loading paths, which are used to generate the snapshots from which the modes are extracted, and the validation loading paths which should be different from the learning paths to show that the method as much wider range of applicability than the learning paths. The distinction between the two types of loading paths has been made for instance in Largenton et al [START_REF] Largenton | Extension of the nonuniform transformation field analysis to linear viscoelastic composites in the presence of aging and swelling[END_REF], where for linear (but ageing) viscoelastic materials the learning paths were uniaxial, but the validation paths were multi-axial. It was found that the predictions of the reduced-order model based on modes determined from uniaxial paths were in very good agreement with full-field simulations. The same conclusion wasreached by the authors [START_REF] Michel | Nonuniform Transformation Field Analysis[END_REF][START_REF] Michel | Computational analysis of nonlinear composite structures using the Nonuniform Transformation Field Analysis[END_REF] for two-phase composites with the earlier version of the NTFA. In all our checks on two-phase composites (not reported here) the present version NTFA-TSO of the model has been found to be more accurate than the earlier version of the NTFA. It is therefore expected that the observations made with the earlier model will apply to the new model, namely accurate predictions for multi-axial loadings with modes generated from uniaxial tests for moderate nonlinearities. Extending these observations to highly nonlin-ear polycrystals is an open question and it is expected that the accuracy of the predictions based on uniaxial modes will deteriorate as the nonlinearity increases. We still believe that for moderate nonlinearity (such as in ice where the rate-sensitivity exponents vary between 2 and 4), uniaxial learning paths can be used successfully. It should however be clearly recognized that a reduced-order model based a finite (and hopefully small) number of tests cannot be expected to be accurate for any loading path. The learning paths, their direction and their amplitude, should contain enough information and should not be too different from the validation paths, or the paths used in the application of the model. At high values of the nonlinearity exponents, it is likely that the results of the reduced-order model based only on uniaxial paths will not be as satisfactory as they would be for moderate values of these exponents. One could use more complex learning loading paths involving multi-axial stress states. The POD will still deliver a set of modes and we believe (but this remains to be checked) that the number of modes will not be significantly larger than that for uniaxial paths, when all directions are covered (6 directions in stress space).How to design properly the multi-axial loading paths is an open question. So far, this choice is left to the user and, admittedly, this is more art than science. A more rational way of designing optimal (or nearly optimal) paths is desirable.

A second direction of investigation would be to explore the possibility of determining the modes "on-the-fly", so that they follow better the actual loading paths. Alternatives to the POD do exist, such as the Proper Generalized Decomposition (see for instance Chinesta and Cueto [START_REF] Chinesta | PGD-Based Modeling of Materials, Structures and Processes[END_REF]). The benefit in terms of computational cost is clear as refined modes would be introduced only when needed. The downside of such an approach is that the gain would be mostly computational and no constitutive relations would be established, which is one of the main aims of the present study.

Finally another important question which is not addressed here pertains to the control of the error between the exact model and the reduced-order model. A related question is the convergence of the reduced-order model to the exact model in the limit as the number of modes goes to infinity. The control of the error introduced by reduced-order modeling is a difficult question in general, and in particular in the present model where two different approximations are made, first the projection on a reduced basis and second the linearization fo the potential by the TSO expansion. The error due to the first approximation only, i.e. the error due to the POD itself with the exact potential, is the error between the reduced model ( 29) and the exact model. A similar question is adressed by Homescu et al [START_REF] Homescu | Error estimation for reduced-order models of dynamical systems[END_REF] in a general framework, and the error due to small perturbations in the coefficients of the ODE are also estimated. However the examples given in this paper are either linear or nonlinear with a quadratic and explicit nonlinearity. The error estimates are based on an adjoint problem the form of which is not straightforward to anticipate in the present micromechanical context. This direction of research deserves to be pursued in a future work, but its feasibility has not been checked by the authors. The second source of error in the present NTFA-TSO model is inherent to the tangent linearization. It should be noted that the very same question of estimating the error between the exact effective potential and its TSO approximation arises in nonlinear homogenization (without any mention of reduced-order modeling). To the authors' knowledge, this question has not been adressed in the literature. These two remarks show that estimating the error due to the model reduction is an open and important question which deserves specific attention in future work.

Conclusion

The new version of the Nonuniform Transformation Field Analysis (Michel and Suquet [START_REF] Michel | A model-reduction approach in micromechanics of materials preserving the variational structure of constitutive relations[END_REF]) has been extended to polycrystals. The findings of this study can be summarized as follows:

1. It was found that the general framework of constitutive relations with internal variables is convenient to apply model-reduction techniques based on reduced basis. Specific examples of two crystal plasticity laws have been explored in details.

2. In the first version of the NTFA ( [START_REF] Michel | Nonuniform Transformation Field Analysis[END_REF]) the plastic modes were restricted to have their support in a single phase. This restriction prohibited the application of the NTFA to polycrystals. This restriction is removed in the present work.

3. The evolution equations for the reduced variables which have been proposed are based on an expansion to second order (TSO) of an underlying potential, when it exists. They can be entirely expressed in terms of quantities which can be pre-computed once for all (off-line). Roughly speaking, these pre-computed quantities depend on the average and fluctuations per grain of the modes and of the associated stress fields.

4. The accuracy of the new NTFA-TSO model has been assessed by comparison with fullfield simulations. It is excellent for mild nonlinearities (strain-rate sensitivity exponent of about 3 in ice), but deteriorates for strong nonlinearity.

5. The speed-up obtained with the new reduced evolution equations over the full-field computations is more than 2 orders of magnitude.

The first relation in [START_REF] Fish | Computational plasticity for composite structures based on mathematical homogenization: Theory and practice[END_REF] can be inverted as

p s = 1 Qb S s =1 (h -1 ) s,s (r s -r ini,s ),
and the differential equation for p s reads as

ṗs = 1 + 1 Q S s =1 (h -1 ) s,s (A p s -A ini,p s ) ∂ψ ∂A βs (A) ,
where A ini,p s = -r ini,s . Finally the constitutive relations ( 8), ( 9), ( 10) and ( 11) can be written as α = F(A) where α and A are given by ( 16) and ( 20) respectively and

F(A) =       F v (A) F βs (A) F ps (A)       =            ∂ψ ∂A v (A) ∂ψ ∂A βs (A) + d s c s ∂ψ ∂A βs (A) A βs 1 + 1 Q S s =1 (h -1 ) s,s (A p s -A ini,p s ) ∂ψ ∂A βs (A)            . (63)
I h ice. Note that (13) can be re-written as (recall that x βs = c s β s = -A βs and divide the two sides of ( Eq. ( 14) is simply

ṗs = S s =1 h s,s ∂ψ ∂A β s (A) .
Finally the constitutive relations ( 8), ( 12), ( 13) and ( 14) can be written in the form (3) with α and A given by ( 16) and ( 20) respectively and

F(A) =       F v (A) F βs (A) F ps (A)       =            ∂ψ ∂A v (A) ∂ψ ∂A βs (A) + d s c s ∂ψ ∂A βs (A) A βs + e s c s A βs S s =1 h s,s ∂ψ ∂A β s (A) .            . ( 65
)
Appendix B: effective potentials

In order to check that ( 25) and ( 26) define the effective response of the polycrystal (the correct average stress history), note first that the solution ε of the variational problem in ( 26) is precisely the solution of the local problem [START_REF] Kattan | Overall damage and elastoplastic deformation in fibrous metal matrix composites[END_REF]. It remains to prove that the force associated with ε is precisely the average of the stress solution of [START_REF] Kattan | Overall damage and elastoplastic deformation in fibrous metal matrix composites[END_REF] and that the differential equation satisfied by α is consistent with [START_REF] Largenton | Extension of the nonuniform transformation field analysis to linear viscoelastic composites in the presence of aging and swelling[END_REF]. To this end, one has to compute the thermodynamic forces associated with ε and α,

∂ w ∂ε (ε, {α(x)} x∈V ) , A = { A x } x∈V , A x = - ∂ w ∂α(x) (ε, {α(x)} x∈V ) .
The first force can be computed using Hill's lemma

∂ w ∂ε (ε, α) = ∂w ∂ε (ε, α) : ∂ε ∂ε = σ : ∂ε ∂ε = σ : ∂ε ∂ε = σ, since ∂ε ∂ε = I.
The first force is therefore the average stress associated with the solution of [START_REF] Kattan | Overall damage and elastoplastic deformation in fibrous metal matrix composites[END_REF]. The second force is computed by taking the Frechet derivative of w

-A. α * = lim t→0 w (ε, α + t α * ) -w (ε, α) t = lim t→0 w (ε, α + t α * ) -w (ε, α) t = ∂w ∂α (x).α * (x) , (66)
which implies that A x coincides with the local force A(x) at point x. Following an identical procedure it can be shown that

∂ ϕ ∂ ˙ α ( ˙ α). α * = ∂ϕ ∂ α ( α).α * (x) . (67)
The evolution equations for the standard model corresponding to the potentials (26) read as

∂ w ∂ α (ε, α) + ∂ ϕ ∂ ˙ α ( ˙ α) = 0. (68)
which, by virtue of (66) and (67) shows that the evolution equations for the field of internal variables α(x) are precisely [START_REF] Largenton | Extension of the nonuniform transformation field analysis to linear viscoelastic composites in the presence of aging and swelling[END_REF].

Appendix C: TSO approximation for the two crystal plasticity models

C.1 Evolution equation for a β

For the non-standard variable β s (and similarly for p s ) the TSO approximation of the evolution equations for the generalized forces consists in making use of a second-order expansion of F βs in (53), where is standard and the TSO approximation of this is

F βs (A) = ∂ψ ∂A βs (A) + d s c s ∂ψ ∂A βs (A) A βs + e s c s A βs , (69)
∂ ψ T SO ∂ξ ( ) βs (A) = N g=1 c (g) - ∂f s ∂τ (τ (g) s -x (g) s , r (g) s ) ∂x (g) s ∂ξ ( ) βs + 1 2 ∂ 2 f s ∂τ 2 (τ (g) s -x (g) s , r (g) s ) ∂C (g) (τ s -x s ) ∂ξ ( ) βs - 1 2
∂ 3 f s ∂τ 3 (τ (g) s -x (g) s , r (g) s )C (g) (τ s -x s ) ∂x (g) s ∂ξ ( ) βs . (70)
The TSO approximation of the second term in F βs (A) ∂A βs ∂ξ

( ) βs is N g=1 c (g) d s c s ∂ψ ∂A βs (A) x s ∂x s ∂ξ ( ) βs (g) = N g=1 c (g) d s c s ∂f s ∂τ (τ (g) s -x (g) s , r (g) s ) x (g) s ∂x (g) s ∂ξ ( ) βs + 1 2
∂C (g) (x s ) ∂ξ ( ) βs + ∂ 2 f s ∂τ 2 (τ (g) s -x (g) s , r (g) s ) C (g) (τ s -x s , x s ) ∂x (g) s ∂ξ ( ) βs - 1 2 ∂C (g) (τ s -x s ) ∂ξ ( ) βs x (g) s + 1 2 ∂ 3 f s ∂τ 3 (τ (g) s -x (g) s , r (g) s )C (g) (τ s -x s )x (g) s ∂x (g) s ∂ξ ( ) βs sign ∂f s ∂τ (τ (g) s -x (g) s , r (g) s ) , (71)
where C (g) (τ sx s , x s ) = (τ sτ (g) sx s + x (g) s )(x sx (g) s ) (g) .

The third term in F βs (A) ∂A βs ∂ξ ( ) βs

does not require any approximation and is simply

C.2 Evolution equation for the isotropic hardening variables

According to [START_REF] Weertman | Creep of ice[END_REF], the evolution equation for the reduced variable p

It results from (48) that up to second-order in the fluctuations of τ s and x s (p s has no fluctuations), ∂ψ ∂A βs (A) ∂f s ∂τ (τ (g) sx (g) s , r (g) s ) + ∂ 2 f s ∂τ 2 (τ (g) sx (g) s , r (g) s )(τ sx sτ (g) s + x (g) s )

+ 1 2

∂ 3 f s ∂τ 3 (τ (g) s -x (g) s , r (g) s )(τ s -x s -τ (g) s + x (g) s ) 2 sign ∂f s ∂τ (τ (g) s -x (g) s , r (g) s ) , (75)
Therefore ṗ(g) s = 1bp (g) s ∂f s ∂τ (τ (g) sx (g) s , r (g) s ) + 1 2 ∂ 3 f s ∂τ 3 (τ (g) sx (g) s , r (g) s )C (g) (τ sx s ) , (76) where = sign ∂f s ∂τ (τ (g) sx (g) s , r (g) s ) .

Ice. For ice, the differential equation ( 73

.

Alternatively a differential equation for r ,

where the last equality has been obtained by noting that only A ps in grain g depends on p (g) s . Using the relations a (g) ps = c (g) A (g) ps , r (g) s = -A (g) ps , r s (p) = τ sta,s + (τ ini,sτ sta,s )e -p , one finds a differential equation for r

β s , r (g) s -x (g)
Table 4: Material parameters used in the full-field simulations for single crystals of ice at -10 • C. Units are MPa and s -1 . After Suquet et al [START_REF] Suquet | Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF].

Appendix E: material data for copper

Copper is a FCC material with twelve slip systems. The list of the slip sytems can be found in table 5.

The cubic symmetries can be used to reduce the hardening matrix to 6 parameters h i : selfhardening is expressed through h 0 , whereas latent hardening is expressed by the other 5 h i 's. As suggested by [START_REF] Gérard | Field measurements and identification of crystal plasticity models[END_REF] (chapter 6), the elastic moduli are taken to be isotropic (although this is probably a crude approximation) with Young's modulus E = 120000 MPa and Poisson ratio ν = 1/3. Following [START_REF] Gérard | Comparison of experimental results and finite element simulation of strain localization scheme under cyclic loading[END_REF] the four coefficients c s , d s , r 0,s , b s were taken to be the same for all systems (no dependence on s).

Figure 1 :

 1 Figure 1: Full-field simulations of a tension-compression cyclic test on polycrystalline copper with eight different realizations. (a) Effective response. (b) Distribution of the norm of the local stress field at the end of the cycle (t = 15s).

1 . 12 s=1M s = 30 ,

 11230 NTFA-TSO 1 : M v = 3, M βs = 2 for s = 1, .., 6 (basal and prismatic systems) and M βs = 1 for s = 7, .., 12 (pyramidal systems). 2. NTFA-TSO 2 : M v = M s = 3 and M βs = 2 for s = 1, .., 6 (basal and prismatic systems) and M s = 2 and M βs = 1 for s = 7, .., 12 (pyramidal systems).

Figure 2 :

 2 Figure 2: Creep test on polycrystalline ice. (a) Configuration with 500 grains. (b), (c), (d) Snapshots of the modes for the viscoplastic strain in the formulation NTFA-TSO 1 (the equivalent strain µ eq is shown). (b) Mode 1. (c) Mode 2. (d) Mode 3.

Figure 4 :

 4 Figure 4: Creep test on polycrystalline ice. Influence of the latent hardening parameters h 3 and h 4 on the response of the polycrystalline aggregate. Comparison between full-field simulations (black solid line) and the predictions of the two variants of the NTFA. Left: NTFA-TSO 1 (red dashed line). Right: NTFA-TSO 2 (blue dot-dashed line).

Figure 5 :

 5 Figure 5: Polycrystalline ice. Snapshots of the activity on each family of slip systems at the end of the creep test. (a) basal systems, (b) prismatic systems, (c) pyramidal systems. Top row: fullfield simulations. Middle row: NTFA-TSO 2 after postprocessing. Bottom row: Probability distribution of the slip activity on each family of slip systems.

Figure 6 :

 6 Figure 6: Tension-compression test on polycrystalline copper. NTFA predictions vs full-field simulations (a) Effective response. (b) Probability distribution of the norm of the local stress field at the end of the cycle (t = 15s).

Figure 7 :

 7 Figure 7: Polycrystalline copper. Norm of the local stress field at the end of the cycle (t = 15s). (a) Full-field. (b) NTFA-TSO 1. (c) NTFA-TSO 2. Units are MPa.

  with A βs = -x s and e s = 0 in the model of Méric and Cailletaud. The first term in the expression of F βs (A) ∂A βs ∂ξ ( ) βs

  is obtained by equalling ȧ( ) βs to the sum of the right-hand sides of (70), (71) and (72).

. 1 c

 1 According to[START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] only A ps = -r s depends on a (g) a (g) ps χ (g) (x),∂A ps ∂a (g) ps = 1 c (g) χ (g) (x).

Finally

  the differential equation for p (g) s reduces to ṗ(g) s = F ps (A)(g) .(73)It remains to expand F p s (A)(g) to second-order in the fluctuations of A. This expansion depends on the constitutive relations under consideration.Méric-Cailletaud. The differential equation takes the formṗ(g) s = 1bp (g) s ∂ψ ∂A βs (A) (g) = 1bp (g) s ∂ψ∂A βs (A)(g) .

sc

  can be obtained. The evolution equation (53) for the reduced force a (g ) F p s (A)∂A p s ∂p (g) s (g ) = c (g) F ps (A) (g) ∂A

  (a) Slip systems. (b) Hardening matrix.The hardening coefficients identified in Suquet et al[START_REF] Suquet | Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF] areh 1 = 70 MPa, h 2 = 110 MPa, h 3 = 125 MPa, h 4 = 0. (79)The elastic moduli are (in MPa),L 1111 = L 2222 = 13930, L 1122 = 7082, L 1133 = 5765, L 2323 = L 1313 = 3014, L 1212 = 3424.The remaining material parameters are identical within the same family of slip systems, but differ from one family to the other.

h 0 h 1 h 2 h 3 h 4 h 5 (

 5 MPa.s 1/n ) -(MPa) -(MPa)

Table 1 :

 1 Polycrystalline ice. Comparison of CPU times for the different simulations. Processor

			NTFA-TSO 1	NTFA-TSO 2
		Reference	CPU ratio= FFT/ NTFA CPU ratio= FFT/ NTFA
	data 1	261 499 s.	614 s.	2 279 s.
	h 3 = 125, h 4 = 0		CPU ratio = 425	CPU ratio = 115
	data 2	269 913 s.	569 s.	2 303 s.
	h 3 = 0, h 4 = 0		CPU ratio = 474	CPU ratio = 117
	data 3	274 544 s.	570 s.	2 255 s.
	h 3 = 1100, h 4 = 1200		CPU ratio = 482	CPU ratio = 121

Table 2 :

 2 Polycrystalline copper. Comparison of CPU times for the different simulations. Processor Intel Xeon X5687 @ 3.6 GHz.

		1	NTFA-TSO 2
	Reference	CPU ratio= FFT/ NTFA CPU ratio= FFT/ NTFA
	581 333 s.	1 198 s.	15 964 s.
		CPU ratio = 484	cpu ratio = 36

  . Bas. Bas. Pris. Pris. Pris. Pyr. Pyr. Pyr. Pyr. Pyr.

					BasBas. h 1	h 1	h 1	h 3	h 3	h 3	h 4	h 4	h 4	h 4	h 4
					Bas. h 1	h 1	h 1	h 3	h 3	h 3	h 4	h 4	h 4	h 4	h 4
					Bas. h 1	h 1	h 1	h 3	h 3	h 3	h 4	h 4	h 4	h 4	h 4
					Pris. h 3	h 3	h 3	h 2	h 2	h 2	h 4	h 4	h 4	h 4	h 4
					Pris. h 3	h 3	h 3	h 2	h 2	h 2	h 4	h 4	h 4	h 4	h 4
					Pris. h 3	h 3	h 3	h 2	h 2	h 2	h 4	h 4	h 4	h 4	h 4
					Pyr. h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4
					Pyr. h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4
					Pyr. h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4
					Pyr. h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4
					Pyr. h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4
					Pyr. h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4	h 4
										(b)
	+	1 2	∂ 3 f s ∂τ 3 (τ	(g) s -x	(g) s , r s )sign (g)	∂f s ∂τ	(τ	(g) s -x	(g) s , r	(g)

s ) C (g) (τ sx s ) .

Table 3 :

 3 Ice single crystal.

  A 3 A 6 B 2 B 4 B 5 C 1 C 3 C 5 D 1 D 4 D 6 A 2 h 0 h 1 h 1 h 3 h 4 h 4 h 2 h 4 h 5 h 2 h 5 h 4 A 3 h 1 h 0 h 1 h 4 h 2 h 5 h 4 h 3 h 4 h 5 h 2 h 4 A 6 h 1 h 1 h 0 h 4 h 5 h 2 h 5 h 4 h 2 h 4 h 4 h 3 B 2 h 3 h 4 h 4 h 0 h 1 h 1 h 2 h 5 h 4 h 2 h 4 h 5 B 4 h 4 h 2 h 5 h 1 h 0 h 1 h 5 h 2 h 4 h 4 h 3 h 4 B 5 h 4 h 5 h 2 h 1 h 1 h 0 h 4 h 4 h 3 h 5 h 4 h 2 C 1 h 2 h 4 h 5 h 2 h 5 h 4 h 0 h 1 h 1 h 3 h 4 h 4 C 3 h 4 h 3 h 4 h 5 h 2 h 4 h 1 h 0 h 1 h 4 h 2 h 5 C 5 h 5 h 4 h 2 h 4 h 4 h 3 h 1 h 1 h 0 h 4 h 5 h 2 D 1 h 2 h 5 h 4 h 2 h 4 h 5 h 3 h 4 h 4 h 0 h 1 h 1 D 4 h 5 h 2 h 4 h 4 h 3 h 4 h 4 h 2 h 5 h 1 h 0 h 1 D 6 h 4 h 4 h 3 h 5 h 4 h 2 h 4 h 5 h 2 h 1 h 1 h 0

	System Plane Direction A 2 (111) [011] A 3 (111) [101] A 6 (111) [110] B 2 (111) [011] B 4 (111) [101] B 5 (111) [110] C 1 (111) [011] C 3 (111) [101] C 5 (111) [110] D 1 (111) [011] D 4 (111) [101] D 6 (111) [110] (a)	A 2 (b)

Table 5 :

 5 Slip systems and hardening matrix for copper.

Table 6 :

 6 Méric-Cailletaud model. Material parameters for copper. After[START_REF] Gérard | Comparison of experimental results and finite element simulation of strain localization scheme under cyclic loading[END_REF].

Another possible and frequent choice is to choose the plastic slips γ s as internal variables, at the expense of a larger number of variables.

Other boundary conditions can be considered provided the Hill-Mandel condition is satisfied(Suquet [41]).

Alternatively one could prescribe the path of macroscopic stress.

Acknowledments:

The authors acknowledge the support of the Labex MEC and of A*Midex through grants ANR-11-LABX-0092 and ANR-11-IDEX-0001-02. The authors want to thank the referees for their comments which led to the discussion in section 6.5.

Example 2: polycrystalline copper under cyclic loading

The example considered in this section is inspired by the work of Gérard [START_REF] Gérard | Field measurements and identification of crystal plasticity models[END_REF]. A cyclic tensioncompression test is performed on a polycrystalline copper specimen with σ = σ(t) e 1 ⊗ e 1 , -0.01 ≤ ε 11 ≤ 0.01, ε11 = ± 1 3 10 -2 s -1 .

(57)

The material parameters for copper at room temperature are taken from [START_REF] Gérard | Comparison of experimental results and finite element simulation of strain localization scheme under cyclic loading[END_REF] and are recalled in appendix E

NTFA modes

The response of the configuration of figure 2(a), discretized into 125 3 voxels, was simulated along the cycle (57). A snapshot of all internal variables (ε v , γ s , β s ) was stored every 1/8th of the cycle, for a total of 10 snapshots for each variable. The procedure of section 6.2.1 was applied with the usual value α = 1 -10 -4 of the tolerance. The following number of modes were found by the POD, 1. NTFA-TSO 1 :

2. NTFA-TSO 2 :

M s = 60, M = 5 and M βs = 7, s = 1, .., 12.

Effective response

The effective response predicted by the two variants of the NTFA are compared in figure 6 to full-field simulations. The agreement at the end of the first loading ramp is reasonable but not excellent, especially the hardening modulus is not well predicted. The agreement is more satisfactory at the end of the cycle. Similarly the probability distribution of the local stress field at the end the cycle is not fully satisfactory (see figure 6(b)). Two sources of discrepancy can be identified.

1. First the nonlinearity is very strong (high rate-sensitivity exponent n = 20). It is very likely that the second-order Taylor expansion is not accurate at this level of nonlinearity. It is also likely that taking the second-order expansion around the average stress (or forces) in the grains is not the best choice. Improvement on the TSO expansion is needed.

2. Second, it has been observed that the variables β s undergo large variations (in time) during the cycle. It could be that the modes are not able to follow accurately these variations, leading to inaccuracies in the determination of the back-stresses x s and in the overall hardening.

Appendix A: alternative writing of the constitutive relations for single crystals This aim of this section is to show that the two crystal plasticity models presented in section 2.2.1 and 2.2.2 can written in the form

where the thermodynamic forces associated with the state variables are

given in [START_REF] Holmes | Structures, Dynamical Systems and Symmetry[END_REF]. In both models, the evolution of the viscoplastic strain-rate is governed by a potential ψ

where the potential ψ is the combination of individual potentials on the different slip systems

which have different expressions for the two models Méric-Cailletaud:

(61) The other evolution equations for β s and p s have to be considered separately for the two models.

Méric-Cailletaud crystal plasticity model. Note that (10) can be re-written as (recall that x βs = c s β s = -A βs and divide the two sides of (10) by c s )

Eq. ( 62) can be re-written as

C.3 Useful relations

The following definitions and relations can be useful to the reader:

βsµ ( ,g) βs ) (g) .

In the above, the symbol :: is to be understood as the full contraction of two fourth-order tensors over their four indices.

Appendix D: material data for ice

Ice I h is an HCP material with 12 slip systems. The list of the slip sytems can be found in table 3. The hardening matrix is taken in the simple form shown in table 3. It depends on 4 independent parameters.