H Moulinec

P Suquet

A numerical method for computing the overall response of nonlinear composites with complex microstructure

The local and overall responses of nonlinear composites are classically investigated by the Finite Element Method. We propose an alternate method based on Fourier series which avoids meshing and which makes direct use of microstructure images. It is based on the exact expression of the Green function of a linear elastic and homogeneous comparison material. First the case of elastic nonhomogeneous constituents is considered and an iterative procedure is proposed to solve the Lippman-Schwinger equation which naturally arises in the problem. Then, the method is extended to nonlinear constituents by a stepby-step integration in time. The accuracy of the method is assessed by varying the spatial resolution of the microstructures. The flexibility of the method allows it to serve for a large variety of microstructures.

Introduction

This study is devoted to a numerical method introduced by Moulinec and Suquet [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear properties of composites[END_REF], [START_REF] Moulinec | A FFT-based numerical method for computing the mechanical properties of composites from images of their microstructure[END_REF] to determine the local and overall responses of nonlinear composites. Numerous studies dealt with nonlinear cell calculations by the Finite Element Method (FEM) (see for example Adams and Donner [START_REF] Adams | Transverse normal loading of a unidirectional composite[END_REF], Christman et al [START_REF] Christman | An experimental and numerical study of deformation in metal-ceramic composites[END_REF], Tvergaard [START_REF] Tvergaard | Analysis of tensile properties for a whisker-reinforced metal-matrix composite[END_REF], Michel and Suquet [START_REF] Michel | On the strength of composite materials: variational bounds and numerical aspects[END_REF]). Most of them are limited to "simple" microstructures, one or two inclusions embedded in a volume of matrix. The need to incorporate more detailed information on the microstructure is clearly recognized. Recently, several studies have considered "complex" microstructures involving a significant number of inclusions with irregular shape. Brockenborough et al [START_REF] Brockenborough | Deformation of metal-matrix composites with continuous fibers: geometrical effects of fiber distribution and shape[END_REF], Böhm et al [START_REF] Böhm | Some simple models for micromechanical investigations of fiber arrangements in MMCs[END_REF], Nakamura and Suresh [START_REF] Nakamura | Effects of thermal residual stresses and fiber packing on deformation of metal-matrix composites[END_REF], Dietrich et al [START_REF] Ch | Stress and strain partitioning in Ag-Ni fibre composite under transverse loading. Finite element modelling and experimental study[END_REF], Becker and Richmond [START_REF] Becker | Incorporation of microstructural geometry in material modelling[END_REF] are some of the contributions to this recently developed subject. All were based on the FEM. The difficulties due to meshing and to the large number of degrees of freedom required by the analysis limit the complexity of the microstructures which can be investigated by this method.

A typical example of a complex microstructure which is difficult to mesh and therefore to handle by means of the FEM is shown in Figure 10 taken from the work of Bornert [START_REF] Bornert | Morphologie structurale et comportement mécanique; caractérisations expérimentales, approches par bornes et estimations autocohérentes généralisées[END_REF]. The digital image of this Iron/Silver blend was obtained by Scanning Electron Microscopy (SEM). The initial idea of the method proposed in [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear properties of composites[END_REF] was to make direct use of these digital images of the real microstructure in the numerical simulation. A similar idea can be found in Garboczi and Day [START_REF] Garboczi | An algorithm for computing the effective linear properties of heteregeneous materials: three-dimensional results for composites with equal phase Poisson ratios[END_REF] who used a spring network technique.

The proposed method avoids the difficulty due to meshing. It makes use of Fast Fourier Transforms (FFT) to solve the unit cell problem1 , even when the constituents have a nonlinear behavior. FFT algorithms require data sampled in a grid of regular spacing, allowing the direct use of digital images of the microstructure. The second difficulty (size of the problem) is partially overcome by an iterative method not requiring the formation of a stiffness matrix.

The interest in numerical simulations of the nonlinear response of composites has recently been strengthened by the development of theoretical methods which analytically predict the nonlinear overall behavior of composites ( Willis [15], Ponte Castañeda [START_REF] Castañeda | New variational principles in plasticity and their application to composite materials[END_REF], Suquet [START_REF] Suquet | Overall potentials and flow stresses of ideally plastic or power law materials[END_REF]). Part of the present study provides precise numerical results for uniaxial loadings which could serve as guidelines for theoretical predictions.

The body of the method and the resulting algorithms are presented in section 2. In section 3, the accuracy of the method and several numerical points are discussed (choice of the reference medium, spatial resolution ....). In section 4 the method is applied to determine the local and overall responses of composites with "random" microstructures. In all the cases considered in this study the models have been limited to two dimensional approximations. The first reason for this approximation is the limitation on current computational capability. The second reason is that many microstructural observations are two dimensional.

The numerical method 2.1 Cell problem and boundary conditions

The overall behavior of a composite is governed by the individual behavior of its constituents and by its microstructure. Its effective response to a prescribed path of macroscopic strains or stresses may be determined numerically via the resolution of the so-called "local problem" on a representative volume element (r.v.e.) V . In this study, the "representative" information on the microstructure is provided by an image (micrograph) of the microstructure with arbitrary complexity. The image contains N pixels, and independent mechanical properties are assigned individually to each pixel. Most applications involve only a limited number of phases, although in principle each pixel could be considered as an individual constituent.

The local problem consists of equilibrium equations, constitutive equations, and boundary and interface conditions. All different phases are assumed to be perfectly bonded (displacements and tractions are continuous across interfaces). Displacements and tractions along the boundary of the r.v.e. are left undetermined and the local problem is ill-posed. We choose to close the problem with periodic boundary conditions which can be expressed as follows. The local strain field ε(u(x)) is split into its average E and a fluctuation term ε(u * (x)):

ε(u(x)) = ε(u * (x)) + E or equivalently u(x) = u * (x) + E.x.

By assuming periodic boundary conditions it is assumed that the fluctuating term u * is periodic (notation: u * #), and that the traction σ.n is anti-periodic in order to meet the equilibrium equations on the boundary between two neighboring cells (notation: σ.n -#). This local problem could be solved by means of the FEM ( Suquet [START_REF] Suquet | Elements of homogenization for inelastic solid mechanics[END_REF], Guedes and Kikuchi [START_REF] Guedes | Preprocessing and postprocessing for materials based on the homogenization method with adaptative finite element methods[END_REF]). We propose an alternate method of resolution.

An auxiliary problem

First we consider the preliminary problem of a homogeneous linear elastic body with stiffness c 0 subjected to a polarization field τ (x).

σ(x) = c 0 : ε(u * (x)) + τ (x) ∀x ∈ V div σ(x) = 0 ∀x ∈ V, u * #, σ.n -#    (1)
The solution of (1) can be expressed in real and Fourier spaces, respectively, by means of the periodic Green operator Γ 0 associated with c 0 :

ε(u * (x)) = -Γ 0 * τ (x) ∀x ∈ V, (2)
or ε(ξ) = -Γ0 (ξ) : τ (ξ) ∀ξ = 0, ε(0) = 0 (3)
The operator Γ 0 is explicitly known in Fourier space (see appendix A). When the reference material is isotropic (with Lamé coefficients λ 0 et µ 0 ) it takes the form :

Γ0 ijkh (ξ) = 1 4µ 0 |ξ| 2 (δ ki ξ h ξ j + δ hi ξ k ξ j + δ kj ξ h ξ i + δ hj ξ k ξ i ) - λ 0 + µ 0 µ 0 (λ 0 + 2µ 0 ) ξ i ξ j ξ k ξ h |ξ| 4 . (4)

The periodic Lippman-Schwinger equation

The auxiliary problem can be used to solve the problem of an inhomogeneous elastic composite material with stiffness c(x) at point x under prescribed strain E :

σ(x) = c(x) : ε(u * (x)) + E ∀x ∈ V divσ(x) = 0 ∀x ∈ V, u * #, σ.n -#    (5)
For simplicity E is assumed to be prescribed, although other average conditions could be considered as well (see appendix B for prescribed stresses). A homogeneous reference material with elastic stiffness c 0 is introduced and a polarization tensor τ (x), which is unknown a priori, is defined as :

τ (x) = δc(x) : ε(u(x)), δc(x) = c(x) -c 0 . (6)
Thus, the problem reduces to the periodic Lippmann-Schwinger equation ( Kröner [START_REF] Kröner | Statistical Continuum Mechanics[END_REF]), which reads, in real space and Fourier space respectively:

ε(u(x)) = -Γ 0 (x) * τ (x) + E, ε(ξ) = -Γ 0 (ξ) : τ (ξ) ∀ξ = 0, ε(0) = E    ( 7
)
where τ is given by [START_REF] Michel | On the strength of composite materials: variational bounds and numerical aspects[END_REF]. The Lippman-Schwinger equation is an integral equation for ε(u * ).

The algorithm

Continuous algorithm

The principle of the algorithm is to use alternately ( 6) and [START_REF] Brockenborough | Deformation of metal-matrix composites with continuous fibers: geometrical effects of fiber distribution and shape[END_REF], in real space and Fourier space, respectively, in an iterative scheme, to solve (5):

Initialization : ε 0 (x) = E, ∀ x ∈ V, σ 0 (x) = c(x) : ε 0 (x), ∀ x ∈ V, Iterate i + 1 : ε i and σ i being known a) τ i (x) = σ i (x) -c 0 : ε i (x), b) τ i = F(τ i ), c) Convergence test, d) ε i+1 (ξ) = -Γ 0 (ξ) : τ i (ξ) ∀ξ = 0 and ε i+1 (0) = E, e) ε i+1 = F -1 ( ε i ), f ) σ i+1 (x) = c(x) : ε i+1 (x).                                (8)
F and F -1 denote the Fourier transform and the inverse Fourier transform. This algorithm can be further simplified by noting that

Γ 0 * (c 0 : ε) = ε.
The modified algorithm reads :

Initialization : ε 0 (x) = E, ∀ x ∈ V, σ 0 (x) = c(x) : ε 0 (x), ∀ x ∈ V, Iterate i + 1 : ε i and σ i being known a) σi = F(σ i ), b) Convergence test, c) εi+1 (ξ) = εi (ξ) -Γ0 (ξ) : σi (ξ) ∀ξ = 0 and εi+1 (0) = E, d) ε i+1 = F -1 (ε i+1 ) e) σ i+1 (x) = c(x) : ε i+1 (x), ∀ x ∈ V,                            (9)
Convergence is reached when σ i+1 is in equilibrium. The error serving to check convergence is :

e i = < ||div(σ i )|| 2 > 1/2 || < σ i > || = < ||ξ. σi (ξ)|| 2 > 1/2 || σi (0)|| .
The iterative procedure is stopped when the error e is smaller than a prescribed value (typically 10 -4 in our calculations).

Discrete algorithm

The unit cell is discretized into a regular grid consisting of

N 1 × N 2 pixels (two-dimensional problem), or N 1 × N 2 × N 3 "voxels" (tri-dimensional problem).
The data and the unknowns used in the numerical calculations are images sampled on this grid (N 1 × N 2 or N 1 × N 2 × N 3 arrays). In two dimensions, the coordinates of the pixel labeled by i 1 , i 2 are

x d (i 1 , i 2 ) = (i 1 -1) • T 1 N 1 , (i 2 -1) • T 2 N 2 , i 1 = 1, ...N 1 , i 2 = 1, ...N 2 ,
where T j is the period of the unit cell in j th direction (j = 1, 2). This discretization is classical in image processing. Images of microstructures, obtained for instance by S.E.M. (scanning electron microscopy), can therefore be directly used in calculations without any operation by the user (meshing or interpolation). This discretization is also appropriate for using Fast Fourier Transforms (FFT) packages, which contribute significantly to the performances of the method. The continuous algorithm [START_REF] Nakamura | Effects of thermal residual stresses and fiber packing on deformation of metal-matrix composites[END_REF] has been implemented in the following discrete form :

Initialization : ε 0 (x d ) = E, ∀ x d ∈ V, σ 0 (x d ) = c(x d ) : ε 0 (x d ), ∀ x d ∈ V, Iterate i + 1 : ε i and σ i known at every x d a) σi = FF T (σ i ), b) Convergence test, c) εi+1 (ξ d ) = εi (ξ d ) -Γ0 (ξ d ) : σi (ξ d ) ∀ξ d = 0 , εi+1 (0) = E, d) ε i+1 = FF T -1 (ε i ) e) σ i+1 (x d ) = c(x d ) : ε i+1 (x d ), ∀ x d ∈ V                            (10)
where x d denote the coordinates of pixels in real space, and ξ d denote the N 1 × N 2 corresponding frequencies in Fourier space. To be more specific, the discrete frequencies are (in dimension 2) when N j is even :

ξ j = (- N j 2 + 1) 1 T j , (- N j 2 + 2) 1 T j , ..., - 1 T j , 0, 1 T j , ..., ( N j 2 -1) 1 T j , N j 2 1 T j ,
and when N j is odd :

ξ j = - N j -1 2 1 T j , ..., - 1 T j , 0, 1 T j , ..., N j -1 2 1 T j .
The discrete error serving to check convergence is :

e i = 1 N d ||ξ d • σi (ξ d )|| 2 1/2 || σi (0)|| (where N = N 1 × N 2 is the total number of pixels).
When the spatial resolution is low and when the number N j of discretization point is even, a special attention must be paid to the highest frequency2 ξ j = ± Nj 2 1 Tj , j = 1 or 2. In most FFT packages, the Fourier expansion at these frequencies consists of either cos(ξ j x j ) or exp(-iξ j x j ), instead of the correct expression consisting of the two terms exp(-iξ j x j ) and exp(iξ j x j ). Therefore, even when the stress σ is correctly approached by its Fourier expansion in step a) of the algorithm [START_REF] Ch | Stress and strain partitioning in Ag-Ni fibre composite under transverse loading. Finite element modelling and experimental study[END_REF], the result of step d) may not approach accurately the Fourier expansion of the strain ε at these particular frequencies. This is because Γ0 is neither even nor odd with respect to each individual component ξ j . Oscillations were observed when (4) was used with relatively small values of N j (lower than 128). This problem was fixed by using a different expression of Γ0 in algorithm [START_REF] Ch | Stress and strain partitioning in Ag-Ni fibre composite under transverse loading. Finite element modelling and experimental study[END_REF] at these frequencies

Γ0 = c 0 -1 .
In other terms, the stress σ is forced to 0 by the algorithm at these frequencies when convergence is reached.

Nonlinear Behavior.

The algorithm can be extended to the case in which the individual constituents obey a nonlinear law, written either in terms of stresses and strains (nonlinear elasticity at infinitesimal strain) or in incremental form relating strain-rates and stress-rates (flow theory). The nonlinearity requires an appropriate modification of step e) in algorithm [START_REF] Ch | Stress and strain partitioning in Ag-Ni fibre composite under transverse loading. Finite element modelling and experimental study[END_REF]. In the present study, special attention will be paid to phases exhibiting an incremental elastic-plastic behavior at small strains governed by a J 2 -flow theory with isotropic hardening (although more general constitutive laws can be considered) :

σ = c : ( ε -εp ), εp = ṗ 3 2 s σ eq , σ eq -σ 0 (p) ≤ 0, ṗ ≥ 0. ( 11
)
ε p denotes the plastic strain, s denotes the stress deviator and p denotes the hardening parameter, which coincides with the cumulated plastic strain

ṗ(t) = 2 3 εp ij (t) εp ij (t) , p(t) = t 0 ṗ(s) ds σ eq = 3 2 s ij s ij 1 2
.

The integration in time of the constitutive law [START_REF] Becker | Incorporation of microstructural geometry in material modelling[END_REF] is achieved by means of an implicit scheme which is classical in the analysis of elastic-plastic structures by the FEM method. The time interval (or, alternatively, the loading path) is discretized into subintervals [t n , t n+1 ]. The field equations are solved for (ε n , σ n , p n ), which denote strain, stress and hardening parameter at time t n . Assuming that these fields are known at step n (time t n ), the principal unknown at step n + 1 is ε n+1 . The incremental equations [START_REF] Becker | Incorporation of microstructural geometry in material modelling[END_REF] are discretized by an implicit scheme. The unknown ε n+1 is a compatible strain field such that the associated stress field (by the constitutive law) is in equilibrium. The resulting system of equations to be solved for ε n+1 is nonlinear. The algorithm for the determination of ε n+1 reads (for simplicity the lowerscript (n + 1) is omitted below; superscripts i and i+1 refer to the iterative loop within the step) :

Initialization : ε 0 (x d ) given by [START_REF] Garboczi | An algorithm for computing the effective linear properties of heteregeneous materials: three-dimensional results for composites with equal phase Poisson ratios[END_REF], Compute σ 0 and p 0 from (ε 0 , σ n , ε n , p n ),

Iterate i + 1 : ε i and σ i are known a) σi = FF T (σ i ), b) Convergence test, c) εi+1 (ξ d ) = εi (ξ d ) -Γ0 (ξ d ) : σi (ξ d ) ∀ξ d = 0, εi+1 (0) = E n+1 , d) ε i+1 = FF T -1 (ε i+1 ) e) Compute σ i+1 and p i+1 from (ε i+1 , σ n , ε n , p n )                            (12)
More specifically a) The initial strain ε 0 at time t n+1 is extrapolated (linearly) from ε n and ε n-1 at the two previous time steps t n and t n-1 :

ε 0 (x d ) = ε n (x d ) + t n+1 -t n t n -t n-1 (ε n (x d ) -ε n-1 (x d )), ∀ x d ∈ V. ( 13
)
This choice significantly improves the convergence of the iterative process within the time step.

b) σ i and p i are computed from (ε i , σ n , ε n , p n ) (step e) in algorithm [START_REF] Bornert | Morphologie structurale et comportement mécanique; caractérisations expérimentales, approches par bornes et estimations autocohérentes généralisées[END_REF]) by a radial return method (see appendix C).

3 Convergence and accuracy of the method

Reference medium

The rate of convergence of the algorithm depends drastically on the Lamé coefficients λ 0 and µ 0 of the reference material. After several tests, the best rate of convergence was observed with

λ 0 = 1 2 inf x ∈ V λ(x) + sup x ∈ V λ(x) µ 0 = 1 2 inf x ∈ V µ(x) + sup x ∈ V µ(x)            (14)
The number of iterations at convergence is significantly influenced by several other parameters. First, as shown in Figure 1, it increases with the contrast between the phases (typically the ratio between the elastic moduli of the phases). When the contrast is infinite (rigid inclusions or voids in an elastic matrix), the algorithm no longer converges. Second, the number of iterations at convergence also depends on the complexity of the solution itself. In the example of an elastic ideally plastic matrix reinforced by stiff inclusions, the computing time increases with the tortuosity of the bands where the strain tends to localize (see below).

Implementation of the method on a vector or a parallel computer

The constitutive law acts locally in real space (i.e. applies separately to each individual point x). Similarly, Green's function Γ 0 acts locally in Fourier space, (i.e. applies separately to each individual frequency ξ). From a computational standpoint, the corresponding steps (c and e in the algorithms [START_REF] Ch | Stress and strain partitioning in Ag-Ni fibre composite under transverse loading. Finite element modelling and experimental study[END_REF] or ( 12) ) are performed by independent loops on each individual pixel in real or Fourier space. These steps can consequently be vectorized or parallelized. In addition, optimized FFT packages are available on most vector or parallel computers. The whole algorithm can therefore be efficiently implemented on these machines. It follows from the same argument that the time spent in the steps corresponding to the constitutive law and to the Lippman Schwinger equation varies linearly with the number N of pixels. The CPU time for a FFT varies as N • log 2 N . The time required by the other steps of the algorithm are comparable to the time required by the FFTs. The CPU time t for one iteration can be estimated by

k 1 × N ≤ t ≤ k 2 × N log 2 N,
where k 1 and k 2 are expected to be independent of the size N of the problem. The dependence of the CPU time on the size of the problem is shown in Figure 2. The square unit cell shown in Figure 5 is subjected to uniaxial transverse tension at 0 0 . The volume fraction of fibers is 47.5%. Both the fibers and the matrix are assumed to be elastic with elastic constants given by ( 17) and [START_REF] Suquet | Elements of homogenization for inelastic solid mechanics[END_REF]. The dependence of the CPU time on the size of the problem is approximately linear.

Optimizing the memory occupancy. The Fourier transform of a real valued function has the symmetry property f (-ξ) = f (ξ).

Since all quantities under consideration in our computation are real, this symmetry property allows us to restrict our attention to positive frequencies (the values of the fields for negative frequencies being immediately deduced). The size of the arrays can therefore be divided by 2, provided the FFT package allows for the storage of real numbers as complex numbers with the same memory occupancy.

Performances. Most computations were run on a Cray YMP with peak performance of 333 M F lops. The performance observed with our algorithm was 210 M F lops on the elastic-plastic problem described in section 4 with unit cells discretized into 1024 × 1024 pixels. The typical CPU time on one processor of this computer is less than 30 seconds for an elastic problem (with a spatial resolution of 1024 × 1024 pixels, the ratio between the Young moduli being approximately 6). When the matrix is elastic plastic, the typical CPU time for a run as described in section 4 is 4000 seconds.

Comparison with analytical solutions

To assess the accuracy and the stability of the method we examined two cases for which analytical solutions are available. Laminates. The first example concerns layered materials. As is well-known, the strain field is then uniform within each individual layer and takes different values from one layer to another. The example shown in Figure 3 corresponds to a two-phase material, both phases having equal volume fraction. The layers are parallel to the plane (x 2 , x 3 ). The constitutive materials of the layers were linear elastic with elastic characteristic given by ( 17) and [START_REF] Suquet | Elements of homogenization for inelastic solid mechanics[END_REF]. The applied loading was pure shear parallel to the layers

Σ 12 arbitrary, Σ 11 = Σ 22 = Σ 33 = Σ 13 = Σ 23 = 0.
The image was discretized into 32 × 32 pixels (good results were obtained with an even cruder resolution). The computed local strain field ε 12 is plotted in Figure 3 and shows no oscillation. In addition the numerical solution coincides with the exact solution.

Circular fiber at dilute concentration. The second example concerns the elastic strain field generated by stiff circular fibers placed at the nodes of a square lattice in a more compliant matrix. The exact solution to this problem (with periodic boundary conditions) is not known in closed form (to the authors' knowledge). However when the volume fraction of fibers is small this solution can be accurately approximated by the solution of a simpler problem, where a circular fiber (with radius a) is surrounded by a circular shell of matrix (with radius b) and subject to the boundary condition

u(x) = E.x when r = b,
where the overall strain E is the same as in the original periodic problem. When the imposed loading is an in-plane shear E 12 = 0, other E ij = 0, the displacement field has the form

u r (r, θ) = Ar 3 + Br + C r + D r 3 sin(2θ), u θ (r, θ) = 2λ + 3µ λ Ar 3 + Br + µ λ + 2µ C r -D r 3 cos(2θ),       
where r and θ are the polar coordinates in the plane. A, B, C, D, take different values in the matrix and in the fiber. They solve a system of linear equations expressing the boundary condition at r = b, the absence of singularity at r = 0, the continuity of tractions and displacements at r = a.

According to Saint Venant's principle, the local strain fields in the two problems coincide far from the boundary of the cell. Therefore at low volume fraction of fibers (a 2 /b 2 1), the solutions of the two problems are expected to coincide except in the vicinity of the boundary of the cell. The example presented in Figure 4 corresponds to a/b = 1/16. The spatial discretization used in the numerical calculation was 1024 × 1024. The component ε 12 of the strain field in a square window of width c = 4a is shown in Figure 4 (note that the unit cell itself with width 2b is much larger than the window shown). There is almost no difference between the analytical and the numerical solutions shown in (a) and (b) respectively. A more explicit comparison is made in Figure 4 (c) which shows an horizontal cut through the field ε 12 at x 2 = 0. Except from little undulations inside the inclusion, there is no significant oscillations at the fiber boundary where the field ε 12 is discontinuous. In addition the accuracy of the numerical solution is observed to increase with the spatial resolution. The discrepancy between the numerical and the analytical solutions depends on the spatial resolution and should not be attributed to a Gibbs phenomenon, i.e. to an oscillation of the Fourier series of a function in the vicinity of a discontinuity point. This oscillation is attached to the summation of the Fourier series which is not what the discrete inverse Fourier transform performs.

Discrete Fourier transform. The discrete Fourier transform, when applied to an image discretized into N 1 × N 2 pixels, is the exact Fourier transform of the image when two requirements are met : ( Brault and White [START_REF] Brault | The analysis and restoration of astronomical data via the Fast Fourier Transform[END_REF]) C1 the image is periodic with the same period (T 1 , T 2 ) as the unit cell, C2 the image cut-off frequency f c ( i.e. the frequency above which the Fourier transform of the image vanishes identically) is less than half of the sampling frequency (Shannon's theorem):

f c j < 1 2 N j T j j = 1, 2
The periodic boundary conditions which have been assumed from the true beginning of this study ensure that condition (C1) is met. However, condition (C2) is not met in general. In particular a discontinuous field has no cut-off frequency and there is no discretization able to capture this discontinuity. It is however expected that the solution of the discrete problem approaches the solution of the continuous problem when the image sampling (number of pixels) increases. A high resolution will therefore be required for problems in which high strain or stress gradients are likely to occur.

Influence of spatial resolution

As already stated the influence of the spatial resolution depends on the stress and strain gradients within the phases and therefore on the strength of the phases nonlinearities. The following examples illustrate these general considerations. The method has been applied to simulate the local and overall response of composites reinforced by unidirectional long fibers aligned along the e 3 direction. The geometry of these composites is described by a two-dimensional image of their cross section. Generalized plane strains were assumed :

u 1 (x) = u 1 (x 1 , x 2 ), u 2 (x) = u 2 (x 1 , x 2 ), u 3 (x) = E 33 x 3 . (15)
The overall strain E has four independent components E 11 , E 22 , E 12 , E 33 (the other two are equal to 0). The overall stress Σ also has four independent components. It is possible to prescribe either a path in the space of strains, or a path in the space of stresses, or alternatively some components of the strain and the other components of the stress. Classical plane strains are a particular case of the more general setting considered in [START_REF] Willis | On methods for bounding the overall properties of nonlinear composites[END_REF]. It corresponds to a path in the space of strains along which E 33 is identically 0. The need to introduce generalized plane strain is illustrated by uniaxial tension in the 0 0 direction, which corresponds to a path in the space of stresses along which

Σ 11 arbitrary, Σ 22 = Σ 12 = Σ 33 = 0. ( 16
)
The axial component E 33 of the strain is unknown and determined a posteriori by the condition Σ 33 = 0.

The assumption of generalized plane strains reduces (5) to a two-dimensional problem for the two unknowns

(u * 1 , u * 2 ).
Two classical configurations were investigated in which the fibers were placed at the nodes of a square or hexagonal lattice. The fibers were assumed to be elastic, isotropic, and characterized by a Young modulus and a Poisson ratio :

E f = 400 GPa, ν f = 0.23. ( 17
)
The fiber volume fraction was 47.5 % (for comparison, we chose the same volume fraction as in [START_REF] Böhm | Some simple models for micromechanical investigations of fiber arrangements in MMCs[END_REF]). The behavior of the matrix was varied from linear elasticity to elasto-plasticity with hardening so as to study the effect of the nonlinearity on the accuracy of the method. All the constitutive laws of the matrix which were considered can be put in the incremental form [START_REF] Becker | Incorporation of microstructural geometry in material modelling[END_REF]. Its isotropic elastic properties were characterized by a Young's modulus and Poisson coefficient

E m = 68.9 GPa, ν m = 0.35. ( 18
)
The plastic properties of the matrix were governed by the Von Mises criterion

σ eq ≤ σ 0 + Hp. (19)
The initial yield stress σ 0 was either infinite (pure linear elasticity) or given by σ 0 = 68.9 MPa. The hardening modulus H was either 0 (perfectly plastic behavior) or H = 1 171 MPa (isotropic linear hardening).

The influence of spatial resolution on the accuracy of the results was studied. The spatial resolution of the image is determined here through the square root of the total number of pixels contained in the image divided by the number of fibers in the image. For the square array, with N 1 × N 1 pixels and a single fiber in the unit cell, the spatial resolution is exactly N 1 . The hexagonal array can be viewed as a rectangular array, thus allowing the use of the Fourier technique in orthogonal coordinates, instead of the natural nonorthogonal coordinates defined by the two unit vectors of the hexagonal lattice (see Figure 5). The rectangular unit cell contains 1 + 4 × 1 4 = 2 fibers. The number of pixels along the first direction x 1 is 2 times larger than the number of pixels in the second direction x 2 . The spatial step in x 2 is 2 √ 3/3 times larger than the step in x 1 . Therefore in the hexagonal array, the spatial definition as defined above is again N 1 for an image containing 2N 1 × N 1 pixels.

Both unit cells were submitted to uniaxial tension at 0 0 and 45 0 in the sense of [START_REF] Castañeda | New variational principles in plasticity and their application to composite materials[END_REF]. The results of the overall response of the composite are shown in Tables 1 to 6. The initial response of the composite is linear and its slope defines the overall Young's modulus of the composite. When the matrix is elastic ideally plastic the overall stress applied to the composite in the direction of tension reaches (asymptotically) a limit which defines the overall flow stress of the composite. When the matrix is governed by a linear hardening, the stress-strain curve of the composite exhibits a nonlinear transition to an asymptotically linear (affine) response. The slope of this limit response is the overall hardening modulus of the composite.

Each table gives an overall material constant as a function of the spatial resolution of the image. The "error" was estimated as the relative difference between the result at a given resolution and the result at the finest resolution.

These results suggest the following remarks.

1. When both constituents are linearly elastic, the overall stiffness is not very sensitive to spatial resolution. Even at the lowest resolution (32 × 32 pixels/fiber), the estimated error was under 1% in all cases.

2. When the matrix is elastic plastic, the local and overall responses are sensitive to spatial resolution. The strain fields exhibit a strong tendency to concentrate in thin bands. The higher the nonlinearity, the thinner the bands. These stiff gradients in strain require high spatial resolution to be correctly captured.

3. The solutions may even be discontinuous when the matrix is elastic-perfectly plastic. This explains the relatively high errors at low resolution: about 15% for the square array of fibers in an elastic-perfectly plastic matrix under tension at 0 • , with a resolution of 32 × 32 pixels/fiber. Shear bands can form in the matrix under tension at 45 • . These shear bands correspond to a mode of deformation of the r.v.e. in plane strains. Therefore, for this particular loading, the effective behavior of the composite depends only on the behavior of the matrix. The overall flow stress of the composite coincides with the flow stress of the matrix under plane strain conditions, i.e. 2σ0 √ 3 . The formation of a slip plane through the matrix is well captured by the numerical method and explains the precision of the numerical result for this particular loading.

4. When the matrix has linear hardening, the strain fields are more regular than in the perfectly plastic case. The local and overall responses of the composite are less sensitive to spatial resolution. The error on the hardening modulus is about 7.5% with a resolution of 32 × 32 pixels/fiber.

This study of the influence of spatial resolution led us to use a resolution of 128 × 128 pixels/fiber in most of the examples presented in the next section.

In this section we investigate the influence of the geometrical arrangement of the fibers on the local and overall responses of nonlinear composites. Attention is again restricted to two-dimensional problems, i.e. to composites reinforced by aligned fibers. The fiber arrangement is determined by a two-dimensional image of the composite cross section.

Configurations

Two classes of fiber arrangement, regular and random, were considered. The fibers were identical circular disks and they were not allowed to overlap (impenetrability condition) except in section 4.3. In most simulations the fiber volume fraction was prescribed to 47.5%, except in section 4.3.

Standard fiber distribution. The "standard" configurations consist of a single fiber placed at the nodes of a square or an hexagonal lattice (see preceeding section). Most F.E.M. cell calculations reported in the literature are based on these standard configurations with the exceptions of Brockenborough et al (1991) and [START_REF] Böhm | Some simple models for micromechanical investigations of fiber arrangements in MMCs[END_REF] who investigated the effect of disorder in the fiber arrangement on the overall transverse properties of composites.

Random fiber distribution In the "random" configurations, the centers of the fibers were placed at random in the unit cell, subject only to the constraints of impenetrability and periodicity. The latter constraint implies that, when a fiber overlaps the boundary of the unit cell, it is split into two parts I and II (see Figure 6 ) to fit in the unit cell. The size of the images was the largest one allowed by the memory on our computer and compatible with a resolution of 128 × 128 pixels per fiber. These two constraints led to unit cells discretized into 1024 × 1024 pixels and containing up to 64 fibers.

Impenetrable fibers

Twenty three different configurations of 64 impenetrable fibers were generated randomly in the unit cell. The fibers were assumed to be elastic with material properties given by [START_REF] Suquet | Overall potentials and flow stresses of ideally plastic or power law materials[END_REF]. The matrix was an elastic plastic material governed by a J 2 flow theory [START_REF] Becker | Incorporation of microstructural geometry in material modelling[END_REF] with material properties given by ( 18) [START_REF] Guedes | Preprocessing and postprocessing for materials based on the homogenization method with adaptative finite element methods[END_REF]. The local and global responses of each configuration to a transverse uniaxial tension in the 0 0 direction (according to ( 16)) were computed with the above described method. The square array and hexagonal array were also subjected to transverse tension in the 0 0 and 45 0 directions.

Local and overall responses

The stress-strain curves predicted by the simulation are shown in Figure 7. The solid line corresponds to the mean response (average of the stress-strain curves over the 23 configurations).

These results call for the following comments :

1. The fibers were stiff and perfectly bonded to the matrix. Therefore, although the strain E 33 in the axial direction was not imposed a priori (Σ 33 was prescribed to 0), it was relatively small along the whole loading path. The strain state was consequently close to the plane strain state, explaining the strain concentrations observed in the perfectly plastic matrices. As is well known, plane strain is more favorable to these strain concentrations than is pure uniaxial tension.

2. The square lattice has a marked transverse anisotropy which is strengthened by the nonlinear behavior, which gives raise to different responses when the direction of tension makes an angle of 0 0 or 45 0 with one of the axes of the square lattice. The low value of the flow stress in the diagonal direction (45 0 ) is due to a shear plane passing through the matrix. Indeed, when a plane of shear can be passed through the weakest phase of a composite, the shear strength of the composite is exactly the strength of the weakest phase ( Drucker (1959)). In tension (under plane strains) in a direction inclined at 45 0 on this plane, the transverse flow stress of the composite is 2σ m 0 / √ 3. This is the flow stress observed in Figure 7 and Table 4 (2σ m 0 / √ 3 79.56 MPa). In conclusion, except at low volume fractions, the square array should not be used to investigate the transverse properties of transversely isotropic nonlinear composites.

the flow stress in the diagonal direction (45 0 ) is again due to a shear plane passing through the matrix.

In conclusion, the hexagonal lattice should be used with care to predict the transverse properties of nonlinear composite systems, even for hardening matrices.

4. The deviation from the average of the transverse Young's moduli computed on the different configurations is small. By contrast, the deviations in the other properties (flow stress, hardening modulus) are higher and may be attributed to the combined effects of nonlinearity and incompressibility.

5. The local plastic strains showed significant differences between the ideally plastic case and the hardening case. For the former, the strain concentrates in thin bands in the matrix. In most configurations, only a small percentage of the matrix contributes to the plastic dissipation. The overall flow stress of the composite is observed to be directly related to the "tortuosity" of these bands. Two different configurations with the corresponding zones of strain concentration are shown in Figure 8. In the first configuration slip bands inclined at approximately 45 0 on the direction of traction can be passed through the matrix, resulting in a low flow stress. Conversely, the fiber arrangement in the second direction inhibits long-range slip bands and causes these bands to deviate or the plastic deformation to spread into wider zones. The plastic dissipation and the flow stress are higher in the second configuration than in the first one. Adding more fibers in the undeformed zones would not change the plastic dissipation, or in other terms, would not affect the flow stress of the composite. These results lead us to think that, when the matrix is perfectly plastic, the geometrical parameter which governs (at first order) the flow stress of the composite is not the volume fraction of the fibers but, instead, the length of the shortest path passing through the matrix at an angle of approximately 45 0 in tension, or 0 0 in shear.

6. When the matrix is a hardening material, the plastic strain spreads all over the matrix (see Figure 8). The whole matrix contributes (although non homogeneously) to the plastic dissipation and, consequently, to the overall strengthening of the composite. In this case, the volume fraction of the fibers seems to be the relevant geometrical information (at least to first order) to predict the overall hardening of the composite.

7. In spite of the differences in the maps of plastic strains in the ideally plastic material and in the hardening matrix, the "stiffest" (respectively the "weakest") configurations in the ideally plastic case remain the stiffest (respectively the weakest) configurations in the hardening case.

Model size

The present section deals with the "representativity" of a unit cell in two aspects. First, does the unit cell contain enough heterogeneities so that the computed effective properties no longer depend on the cell size? Second, how much do different unit cells randomly generated with the same volume fraction and number of heterogeneties differ from each other? Several series of microstructures containing 4, 9, 16, 36, 64 or 256 impenetrable fibers randomly placed in the unit cell were generated. The volume fraction of fibers was identical in all simulations (47.5%) and the spatial resolution was also fixed (128 × 128 pixels/fiber). The total number of pixels in each image was therefore the number of fibers multiplied by ×128 × 128. The fibers and the matrix were respectively assumed to be elastic and elastic-perfectly plastic with materials properties given by ( 17) and ( 18) [START_REF] Guedes | Preprocessing and postprocessing for materials based on the homogenization method with adaptative finite element methods[END_REF]. The loading was uniaxial transverse tension at 0 0 (see [START_REF] Castañeda | New variational principles in plasticity and their application to composite materials[END_REF]). Statistical data on the computed Young's moduli as a function of the number of fibers in the unit cell are reported in Table 7. The mean Young's modulus and its standard deviation are defined as

Ē = 1 N s i=1,Ns E i , σ(E) = 1 N s -1 i=1,Ns (E i -Ē) 2
where E i is the Young's modulus of the i th microstructure and N s is the number of different microstructures. The error on the mean is classically estimated by the ratio

σ(E) Ē√ N s
Similar data on the overall flow stress of the composite are given in Table 8. The number of fibers in the unit cell does not significantly influence the mean overall properties, provided a lower number of fibers is compensated by a higher number of configurations. The mean Young's modulus and the mean flow stress of configurations with four fibers differ from those of configurations with 256 fibers by 0.56% and 0.74% respectively. These differences are comparable to the error on the mean itself (0.13% and 0.23% for the Young's modulus and the flow stress for configurations with 256 fibers). This is an illustration of the ergodic property : spatial averaging on one large sample is equivalent to ensemble averaging on many small samples. A related observation is that the standard deviations of the overall properties decrease as the number of fibers increases.

Spacing between fibers.

In the above analyses, the fibers were placed randomly in the unit cell with impenetrability as the only restriction. The effects of imposing a minimal space between fibers are of interest for at least two reasons. First, when the minimal spacing between the centers of the fibers increases, the ordering of the microstructure increases. As a limit case, when this minimal spacing reaches 2 √ 3 . S N (S is the surface of the unit cell, N is the number of fibers), the microstructure is completely determined and coincides with the centered hexagonal arrangement. Second, numerical difficulties could be expected when two neighboring fibers are nearly touching. Indeed, when the spatial resolution is not fine enough, the method cannot capture the high strain gradients in the necks between the two fibers.

Ten configurations with 64 fibers were generated, and a minimal space of 4 pixels between two neighboring fibers was imposed. This distance seemed sufficient to correctly describe strain concentration. The results of this study suggest the following comments :

1. When the matrix is elastic ideally plastic, the mean overall flow stress is Σ 0 = 86.9 MPa (with an estimation error of 0.44 MPa). This value is 2.0% smaller than the value obtained with no restriction on the space between fibers. It lies slightly below the flow stress of the hexagonal array subjected to tension at 0 0 (Σ 0 = 87.9 MPa). However it lies above the flow stress of the square array under tension at 0 0 or 45 0 (Σ 0 = 79.6 MPa) and of the hexagonal array under tension at 45 0 .

2. When the matrix is elastic-plastic with linear hardening, the effective hardening modulus drops significantly : H = 9382 MPa (estimation error = 123.8 MPa), instead of 10002 MPa. But it is still much higher than the hardening modulus predicted with the hexagonal array case (H = 7100 MPa at 0 0 , H = 7420 MPa at 45 0 ).

In conclusion, it seems that the "safety coating" around the fibers leads to a decrease in the overall mechanical properties of the composite, at least at the volume fraction which has been investigated.

Influence of the shape of the fibers

The above analyses show that the overall flow stress of the composite and, to a lesser extent, its overall hardening depend primarily on the tortuosity of shear bands passing through the matrix. Obviously, the volume fraction of the reinforcing phase plays a role in the possibility that such bands are formed, but for a fixed volume fraction, significant differences arise from the differences in the patterning of bands. These shear bands are locked or deviated by the fibers. The overall flow stress of the composite can (empirically) be related to the length of the shortest path passing through the matrix and making an angle of approximately 45 0 with the tensile direction. It can be expected that the shape of the fibers, which act as "shear bands barriers", is important in their capacity to inhibit shear bands. The shape of fibers is important at two levels. First it affects the arrangement of fibers in the unit cell. For instance, it can be favorable to clustering of particles, leaving large areas of inclusions-free matrix where plastic strain is likely to localize. At a smaller scale an elongated particle perpendicular to a shear band will form an effective barrier.

Random microstructures were generated with three shapes of fibers : circular, elliptical (aspect ratio= 3.333), equilateral triangles. The volume fraction was 47.5%. The unit cells contained 64 fibers and were discretized into 1024 × 1024 pixels. The center of the fibers and their orientation were chosen randomly, subject to the contraints of periodicity, impenetrability and given volume fraction. A minimal space of four pixels between two fibers was imposed to correctly capture the high strain gradients in the matrix between two neighboring reinforcements. For each fiber shape, 10 different configurations were tested. The results of the numerical simulations are given in Table 9. The Young's modulus is not significantly affected by the shape of the inclusions, at least for this particular volume fraction and for the contrast of elastic properties which was investigated (investigation of the percolation threshold for highly contrasted phases would probably lead to different conclusions). The mean flow stress of the composite with elliptical inclusions is close to that of the composite with circular inclusions (0.9% higher). However the flow stress is significantly higher for the composite with triangular inclusions (5.2% higher). This "hardening" effect can be attributed to the fact that at a given volume fraction triangles form more efficient barriers to shear band formation. This efficiency can be related to the length of the projection of the fiber orthogonally to the shear bands. The minimal length, the maximal length, and the average length over all possible orientations are reported in Table 10 for each shape of fibers at a given area s. For circular fibers these three quantities are equal to the radius of the fiber (2 s/π).

Penetrable fibers

When the matrix is elastic ideally plastic, the overall response of the composite is strongly influenced by the existence of continuous paths in the matrix, connected from one cell to the other. The contiguity of the matrix obviously plays a crucial role in the formation of these paths, which are ruled out when the matrix is not contiguous.

In order to study this effect, different configurations at different volume fractions were generated with penetrable fibers. The centers of the fibers were first chosen at random. Then the volume fraction of the reinforcing phase was controlled by increasing the radius of the fibers (all fibers at a given volume fraction had identical radius). The matrix was assumed to be elastic perfectly plastic. The results of the simulations can be analyzed as follows :

1. When the fiber volume fraction is small, shear bands can be passed through the matrix. According to Drucker's remark, the resulting overall flow stress of the composite coincides with the flow stress of the matrix under plane strains, 2σ 0 / √ 3. However, when the fiber volume fraction is very small, a nearly homogeneous deformation of the matrix is more favorable (less energy is dissipated in the plastic deformation) and no strain concentration is observed. Then the overall flow stress of the composite stands between the flow stress of the matrix σ 0 and the flow stress of the matrix under plane strains (2σ 0 / √ 3).

2. Over a certain radius, straight shear bands cannot be passed through the matrix. For a given geometrical distribution of fibers, this radius is half of the maximal distance between adjacent parallel lines passing through the centers of the fibers and inclined at ±45 o on the tensile direction. Periodic continuous paths can again be passed through the matrix but they are tortuous. The bands where the plastic strain concentrates have a nonvanishing width. The stress-strain response of the composite again reaches a limit value, one higher than the flow stress of the matrix under plane strains. The overall flow stress increases with the volume fraction of fibers, and the increase is closely related to the tortuosity of the "shear" bands.

3. When the fibers percolate and form a contiguous phase, the matrix loses contiguity. No periodic continuous path can be passed through the matrix. This leads to a drastic modification of the stressstrain curve of the composite, which is no longer limited. The composite behaves asymptotically as an elastic plastic material with linear hardening .

Complex microstructures

To illustrate the capability of the method to deal with complex microstructures, we have considered a real microstructure taken from the work of Bornert [START_REF] Bornert | Morphologie structurale et comportement mécanique; caractérisations expérimentales, approches par bornes et estimations autocohérentes généralisées[END_REF] (see also [START_REF] Bornert | Self-consistent approaches and strain heterogeneities in two-phase elastoplastic materials[END_REF]). The materials studied in [START_REF] Bornert | Self-consistent approaches and strain heterogeneities in two-phase elastoplastic materials[END_REF] were twophase iron/silver blends, manufactured with powder metallurgy techniques. The digital image was obtained by Scanning Electron Microscopy. The microstructure is shown in Figure 10 (a). Clearly meshing this microstructure for application of the FEM would be a considerable task. The present numerical method can handle such a microstructure as easily as the simpler ones shown in previous examples. In the numerical simulation each phase is considered elastic-plastic following a J 2 -flow theory with isotropic hardening of the Von-Mises type. The stress/strain curves for each constituent under uniaxial tension are shown in Figure 10 (c). The applied loading is uniaxial tension in the horizontal direction. The map of equivalent strain is shown in Figure 10 (b) at an overall strain E 11 = 3.3%. In the soft phase (silver in white) the strain is organized in bands which cannot develop over long distances due to the presence of the hard phase (iron in black). A full comparison between simulated and experimental strain maps is difficult to perform essentially because the numerical calculations are two-dimensional whereas the real material is three-dimensional in nature. Only the surface of the specimen is observed and it is in a state of plane stress, whereas the calculations are performed assuming a state of generalized plane strains. In addition the material below the surface plays a significant role on the deformation of the surface itself. The variations between the arrangement of the phases at the surface and below the surface is not taken into account by the numerical model.

Concluding remarks

A new numerical technique has been developed to investigate the local and overall response of nonlinear composites. The advantages of the method are the following :

1. Images of microstructures can be directly used in the analysis, which avoids meshing the microstructure. Complex microstructures can be investigated. Part of the efficiency of the method is due to the use of FFT packages.

2. The iterative procedure does not require the formation or inversion of a stiffness matrix.

3. Convergence is fast.

However the method has some limitations.

1. Convergence is not ensured for materials containing voids or rigid inclusions.

2. The number of degrees of freedom is high by comparison with the FEM (typically an image with 1024 × 1024 pixels is required to deal with 64 fibers). The method can be implemented only on computers with high memory capabilities.

A Green's operator of a linear elastic material

The auxiliary problem of a homogeneous material with stiffness c 0 subject to a periodic polarization field τ plays an important role in the method which has been proposed. Its solution, which can be found in several textbooks (e.g. Mura [START_REF] Mura | Micromechanics of Defects in Solids[END_REF]), can be expressed in terms of the Fourier transform of the polarization field by means of the Fourier transform of the Green's operator of the following systems of equations

σ(x) = c 0 : ε(u * (x)) + τ (x) ∀x ∈ V divσ(x) = 0 ∀x ∈ V, σ.n -#, u * #    (20)
In Fourier space, these equations take the form

σij (ξ) = i c 0 ijkh ξ h û * k (ξ) + τij (ξ), i σij (ξ) ξ j = 0. ( 21
)
(It is hoped that the index i will not be confused with the complex number i = √ -1). Eliminating σij between the two equations in [START_REF] Bornert | Self-consistent approaches and strain heterogeneities in two-phase elastoplastic materials[END_REF] yields

K 0 ik (ξ).u * k = τij (ξ) ξ j ,
where K 0 (ξ) denotes the acoustic tensor of the homogeneous material,

K 0 ik (ξ) = c 0 ijkh ξ h ξ j . Then û * k (ξ) = i N 0 ki (ξ) τij (ξ) ξ j = i 2 (N 0 ki (ξ) ξ j + N 0 kj (ξ) ξ i ) τij (ξ),
where the symmetry of τ has been used and where N 0 (ξ) denotes the inverse of K 0 (ξ). Therefore

εkh (u * ) = i 2 (ξ h û * k (ξ) + ξ k û * h (ξ)) = Γ0 khij (ξ) τij (ξ), ( 22
) with Γ0 khij = 1 4 N 0 hi (ξ) ξ j ξ k + N 0 ki (ξ) ξ j ξ h + N 0 hj (ξ) ξ i ξ k + N 0 kj (ξ) ξ i ξ h , (23)
and τij (ξ) =< τ ij (x)e -iξ.x > .

The strain field induced at each point x of the unit cell V by an initial stress τ can be determined from ( 22), ( 23) and (24). These formulas give the explicit form of the operator Γ 0 and of the operation * considered in section 2: (u * ) = -Γ 0 * τ .

Deatailed expressions of Γ 0 can be found in Mura [START_REF] Mura | Micromechanics of Defects in Solids[END_REF] for different types of anisotropy for the reference medium. Its expression is particularly simple when the material is isotropic with Lamé coefficients λ 0 and µ 0 ; the above expression becomes :

c 0 ijkh = λ 0 δ ij δ kh + µ 0 (δ ik δ jh + δ ih δ jk ) . K 0 ij (ξ) = λ 0 + µ 0 ξ i ξ j + µ 0 |ξ| 2 δ ij N 0 ij (ξ) = 1 µ 0 |ξ| 2 δ ij - ξ i ξ j |ξ| 2 λ 0 + µ 0 λ 0 + 2µ 0 . Therefore: Γkhij (ξ) = 1 4µ 0 |ξ| 2 (δ ki ξ h ξ j + δ hi ξ k ξ j + δ kj ξ h ξ i + δ hj ξ k ξ i ) - λ 0 + µ 0 µ 0 (λ 0 + 2µ 0 ) ξ i ξ j ξ k ξ h |ξ| 4
B Imposing a macroscopic stress direction.

In the above described algorithm the overall strain is prescribed by assessing the value of the Fourier transform of the strain field at the zero frequency :

ε(0) = E.
It is often convenient (or necessary) to impose the overall stress Σ, rather than the overall strain E. A typical example is provided by uniaxial tension in the transverse direction as described by [START_REF] Castañeda | New variational principles in plasticity and their application to composite materials[END_REF]. In strongly nonlinear problems it is even necessary to impose only the direction of the overall stress and to drive the loading by means of an auxiliary parameter (arc length method). The algorithm can be modified to account for loadings in the form Σ = k S 0 and E :

S 0 = t, (25)
where S 0 is the prescribed direction of overall stress (by direction of stress we refer to a direction in the 6-dimensional space of stresses), k is the unknown level of overall stress and t, which serves as a loading parameter, is the component of the overall strain in this direction. Then, the overall strain and stress E i et Σ i have to be determined by means of (25). For this purpose, at iterate i, σ i-1 and ε i-1 being known, the loading level t i being known but k i being unknown, E i and Σ i are subject to :

Σ i -c 0 : E i =< σ i-1 > -c 0 :< ε i-1 > Σ i = k i S 0 , E i : S 0 = t i    (26)
Elimination of Σ i yields

E i = k i c 0 -1 : S 0 -c 0 -1 :< σ i-1 > + < ε i-1 > (27)
and

k i = t i + (c 0 -1 :< σ i-1 > -< ε i-1 >) : S 0 c 0 -1 : S 0 : S 0
Therefore the modification brought into the algorithm ( 12) is an additional step to determine E i according to (27), which is then prescribed as the overall strain through :

εi (0) = E i .
It is worth noting that the condition < ε i >= E i is met at each step of the iterative procedure, whereas the equality < σ i >= Σ i is met only at convergence. The difference arises from the fact that σ i is deduced from the constitutive law, whereas Σ i is deduced from (26). Indeed, once convergence is reached, one has

E i = E i-1 =< ε i-1 >,
and, according to (26), Σ i =< σ i-1 >=< σ i > .

C Radial return algorithm

The equations governing a plastic material obeying a J 2 flow theory with isotropic hardening read :

σ = c : ( ε -εp ), εp = 3 2 ṗ s σ eq , (28)
ṗ = 0 when σ eq -σ 0 (p) < 0, ṗ > 0 when σ eq -σ 0 (p) = 0.

[

B    (29)
ε p is the plastic strain, p is the equivalent plastic strain ṗ = 2 3 εp : εp 1/2 . c is the stiffness tensor, assumed to be isotropic and characterized by a bulk modulus k and a shear modulus µ.

Time is discretized into intervals [t n , t n+1 ]. F n denotes the value of a function F at time t n . ε n , σ n and p n denote the strain, stress and equivalent plastic strain at time t n . Given the mechanical fields at step n, and given the strain field ε n+1 at step n + 1, the constitutive law amounts to finding the stress field σ n+1 and the equivalent plastic strain field p n+1 . Replacing time differentiation by a finite difference in (28) provides

σ n+1 -σ n = c : ε n+1 -ε n -εpn+1 × (t n+1 -t n ) .
The elastic prediction is

σ n+1 T = σ n + c : ε n+1 -ε n . (30)
After due account of plastic incompressibility, (30) gives

σ n+1 = σ n+1 T -2(t n+1 -t n ) µ εpn+1 .
Alternatively, making use of the flow rule (28) and of the decomposition of σ n+1 into a spherical stress and deviator stress tr

(σ n+1 ) = tr(σ n+1 T ) = tr(σ n ) + 3k tr(ε n+1 -ε n ) (31) s n+1 = s n+1 T - 3(t n+1 -t n ) µ ṗn+1 σ n+1 eq s n+1 (32)
(31) can be re-written, assuming that there are no initial stresses or strains at time t 0 , tr(σ n+1 ) = 3k tr(ε n+1 )

The radial return method is based on the observation that, according to (32), the deviators s n+1 and s n+1 T are proportional. The Von Mises stresses associated with σ n+1 and σ n+1 T are therefore related through :

σ n+1 eq = (σ n+1 T ) eq -3µ (p n+1 -p n ) (33)
-If (σ n+1 T ) eq < σ 0 (p n ), the step is purely elastic,

σ n+1 = σ n+1 T , p n+1 = p n .
-If (σ n+1 T ) eq ≥ σ 0 (p n ), the material plastifies at step n + 1, σ n+1 eq = σ 0 (p n+1 ) and (33) reduces to :

σ 0 (p n+1 ) + 3µ p n+1 = (σ n+1 T ) eq + 3µ p n
Assuming that hardening is positive (no softening), the function h(p) = σ 0 (p) + 3µp can be inverted to give

p n+1 = h -1 (σ n+1 T ) eq + 3µp n (34)
The case of linear hardening leads to simple inversion. Indeed, in this case, σ 0 (p) = σ 0 + Hp and (34) reduces to

p n+1 = 3µ H + 3µ p n + (σ n+1 T
) eq -σ 0 H + 3µ .

The case of a perfectly plastic material, corresponding to H = 0, is covered by the above relation. When h -1 is not available in a closed form, it can be approximated by linear interpolation. When k ∈ [h(p l ), h(p l+1 )], p = h -1 (k) is approximated by p l + (k -h(p l )) p l+1 -p l h(p l+1 ) -h(p l )

Figure 1 :

 1 Figure 1: Transverse Young's modulus of the composite. Dependence of the number of iterations at convergence on the contrast of the elastic moduli of phases (e ≤ 10 -4 ). Square array as shown in Figure 5. Spatial resolution 128 × 128 pixels. Fiber volume fraction 47.5 %. Poisson coefficients ν f = ν m = 0.35. The stiff phase is the fiber.

Figure 2 :

 2 Figure2: CPU time on one processor of a CRAY YMP as a function of the size N of the problem. The solid line is obtained by linear regression on all points and passes through the origin.

Figure 3 :Figure 4 :Figure 5 :

 345 Figure 3: Two-phase laminate. (E 1 = 68.9 GPa, ν 1 = 0.35, E 2 = 400 GPa, ν 2 = 0.23 ). Volume fraction of both phases 50%. Spatial resolution 32 × 32 pixels. Applied loading: pure shear in the plane (x 1 , x 2 ). (a): map of the local strain field ε 12 . (b) : cut through ε 12 along an arbitrary horizontal line.

Figure 6 :

 6 Figure 6: Periodic unit cell containing 16 circular fibers randomly placed.

Figure 7 :Figure 8 :

 78 Figure 7: Overall stress-strain response computed with the present method. Volume fraction of fibers: 47.5%. (a) Matrix with linear hardening. (b) Ideally plastic matrix Dotted lines: 23 configurations of 64 identical circular fibers placed randomly in the r.v.e. Thick solid line: average of the responses of the random configurations. Square 0 0 (resp: Square 45 0 ): fibers placed at the nodes of a square lattice, tension at 0 0 (resp. 45 0 ). Hexag. 0 0 (resp: Hexag. 45 0 ): fibers placed at the nodes of a hexagonal lattice, tension at 0 0 (resp. 45 0 ).

Figure 9 :Figure 10 :

 910 Figure9: Overall stress-strain response of fiber reinforced materials at different fiber volume fraction f . 100 penetrable circular fibers with increasing radius in an elastic ideally plastic matrix. (E m = 68.9 GPa, ν m = 0.35 σ 0 = 68.9 M P a and E f = 400 GPa, ν f = 0.23 ). Tension at 0 0 .

During the revision of this paper, the attention of the authors was called on a similar work by Müller[START_REF] Müller | Mathematical versus Experimental Stress Analysis of Inhomogeneities in Solids[END_REF] concerning phase transformation.

An error had crept into the expression of the highest frequency in the original paper published in Comput Methods Appl Mech Eng. The authors thank Anthony Rollett for pointing it out. HM 11/12/2020

The hexagonal lattice approaches transverse isotropy. When the matrix is a hardening material, the predictions obtained with the hexagonal lattice underestimate the stiffness of the composite, or at least are located below the average of the predictions for the random configurations in the range of overall deformations considered. Another computation, not reported here, was performed up to 30% of transverse strain, with no modification in the conclusions. A similar observation was made by Brockenborough et al (1991) for another system. When the matrix is ideally plastic, the low value of

Acknowledgements. Most computations were carried out at the Institut Méditerranéen de Technologie in Marseille; the funds being provided by the PACA region. The other computations were carried out at the Institut du Développement et des Ressources en Informatique Scientifique funded by CNRS. The authors are indebted to Michel Bornert for fruitful discussions and for providing the image of the microstructure shown in 10 (a).

End of test

Update

tr(σ n+1 ) = tr(σ n ) + 3k tr(ε n+1 -ε n )


Lucian Florin

Constantin
email: fl.constantin@univ-lille1.fr

Towards Probing a Variation of Fundamental Constants with Optical Clock Transitions of 127 I 2

Keywords: Laser Spectroscopy, Clocks, Rovibronic States, Fundamental Constants

Precision measurements of the molecular iodine optical transitions can be exploited for constraining a possible time variation of the fundamental constants. The sensitivities of the molecular frequencies to a variation of the proton-toelectron mass ratio and of the fine structure constant are calculated. Compact molecular iodine clocks, designed for space applications with improved stability performances, enable fractional frequency reproducibility at the 10 -15 level. The comparison of the optical iodine clocks based on transitions at 532 nm, 514 nm and 502 nm with the Cs frequency standard can constrain a fractional time variation of the proton-to-electron mass ratio and of the fine structure constant at the 10 -14 yr -1 level.

I. INTRODUCTION

A space-time variability of fundamental constants, allowed in theories beyond the Standard Model and General Relativity, can infirm the validity of the local position invariance hypothesis and the equivalence principle [START_REF] Uzan | Varying constants, gravitation and cosmology[END_REF]. In unified theories, the variation of a constant is associated to the variation of other constants. A possible time variation of fundamental constants was tested in laboratory using comparisons between different atomic clocks. The most precise results allowed to constrain a fractional time variation of the fine structure constant α at the 10 -17 yr -1 level and of the proton-to-electron mass ratio µ at the 10 -16 yr -1 level [START_REF] Godun | Frequency ratio of two optical clock transitions in 171 Yb + and constraints on the time variation of fundamental constants[END_REF][START_REF] Huntemann | Improved limit on a temporal variation of mp/me from comparisons of Yb + and Cs atomic clocks[END_REF]. The optical atomic clocks are not sensitive to the variation of µ. In this case one should address the case of the molecular transitions which are intrinsically sensitive to the µ variation. The absolute frequency measurements of a rovibrational transition of SF6 allowed to determine a modelfree constraint of the fractional time variation of µ at 5.6 × 10 -14 yr -1 [START_REF] Shelkovnikov | Stability of the proton-to-electron mass ratio[END_REF]. The strategy for improving the constraint with laboratory tests is to increase the accuracy of the clock transitions, to exploit transitions that are highly sensitive to the variations of fundamental constants, and to perform measurements during a long period of time.

This contribution addresses the sensitivity to variations of the fundamental constants from the molecular iodine optical transitions spanning the spectral domain between 500 nm and 700 nm. The Doppler-free lines of molecular iodine form a dense grid of optical frequency references that were exploited to build frequency standards [START_REF] Riehle | The CIPM list of recommended frequency standard values: guidelines and procedures[END_REF].

II. SENSITIVITY TO VARIATIONS OF α AND µ FROM THE OPTICAL TRANSITIONS OF THE MOLECULAR IODINE

The rovibronic energy levels of 127 I2 in the ground X electronic state and in the excited B electronic state are described with the Dunham series expansion :

( ) ( ) [ ] ( ) ( ) [ ]   + ′ ′       + ′ + = + ′ ′ ′ ′       + ′ ′ = ′ ′ ′ ′ l , k l k B l , k e B l , k l k X l , k X J J v Y T ' J , ' v E J J v Y J , v E 1 2 1 1 2 1 (1)
in function of the vibrational quantum numbers v' and v'', the rotational quantum numbers J' and J'', the electronic energy Te, and the Dunham parameters

. Y B , X l , k
The sensitivity of the energy levels to a variation of µ arises solely from the dependence of the Dunham parameters on the nuclear mass :

0 2 =       + - = µ ln d T ln d l k µ ln d Y ln d e X , B l , k (2)
The rovibronic energy levels of the X and B electronic states of 127 I2 are issued from the fine structure splitting of the iodine atom, that depends quadratically on the atomic number and the fine structure constant. The parameters of interest in this contribution are the adimensional sensitivity coefficients of the 127 I2 transitions to a variation of the constants. The sensitivity coefficients to a variation of α are assumed here to be the same as for the case of the iodine atom Kα=2. The sensitivity coefficients to a variation of µ are calculated using eq. (1,2) and the molecular parameters from [START_REF] Gerstenkom | Description of the absorption spectrum of iodine recorded by means of Fourier transform spectroscopy : the (B-X) system[END_REF]. The results are plotted in Fig. 1 for a selection of molecular iodine transitions.

III. METROLOGICAL PERFORMANCES WITH COMPACT CLOCKS BASED ON THE MOLECULAR IODINE TRANSITIONS

Compact optical setups with lasers stabilized on the molecular iodine transitions have been developed in view of several space applications [START_REF] Barbarat | Compact and transportable iodine frequency stabilized laser[END_REF][START_REF] Gürlebeck | BOOST: A satellite mission to test Lorentz invariance using high-performance optical frequency references[END_REF], including gravitational wave detection and test of the equivalence principle. The approaches are based on lasers with low intrinsic phase noise which are frequency doubled or tripled to be locked onto molecular iodine lines with high signal to noise ratio and high quality factor that are probed with the frequency modulation transfer technique. The transitions to the dissociation limit of the B electronic state have the advantage to display narrow natural linewidths [START_REF] Cheng | Sub-Doppler molecular-iodine transitions near the dissociation limit (523-498 nm)[END_REF]. A setup developed at Laboratoire SYRTE (Fig. 2) exploits a narrow linewidth Telecom laser source that is frequency tripled using the fibered waveguide periodically-poled lithium niobate (PPLN) technology and locked to a narrow line at 514 nm [START_REF] Barbarat | Compact and transportable iodine frequency stabilized laser[END_REF]. The molecular iodine stabilized laser is compared with a similar Telecom laser locked onto an ultrastable optical cavity [START_REF] Argence | Prototype of an ultra-stable optical cavity for space applications[END_REF], for which the frequency is compared to the SYRTE's frequency standards. The stabilization of the pump and probe laser intensities, the strict control of the residual amplitude modulation and the stabilization of the cold finger temperature of the cell allowed to reach an Allan deviation of 2.9 × 10 -14 τ -1/2 with a minimum of 5.7 × 10 -15 for 100 s of integration time [START_REF] Barbarat | Compact and transportable iodine frequency stabilized laser[END_REF]. The reproducibility of laser frequency was estimated fractionally at the 10 -15 level [START_REF] Philippe | Sources laser à 1,5 µm stabilisées en fréquence sur l'iode moléculaire[END_REF]. Comparable metrological performances have been obtained with a semiconductor laser that was frequency doubled and locked to a molecular iodine line at 532 nm [START_REF] Gürlebeck | BOOST: A satellite mission to test Lorentz invariance using high-performance optical frequency references[END_REF].

IV. ESTIMATED CONSTRAINTS FOR THE TIME VARIATIONS OF α AND µ FROM MOLECULAR IODINE CLOCK COMPARISONS

The time variation of the ratio between the iodine molecular clock frequency and the Cs clock frequency may be expressed in function of the time variations of µ and α as : This contribution evaluates the constraints to the variations of α and µ that can be inferred from absolute frequency measurements of the three molecular iodine transitions that are indicated in Table 1. The time variations of the absolute molecular iodine frequencies are described with three eq. ( 3) with the relevant sensitivity coefficients. The estimation is based on 100 frequency measurements performed with the same periodicity during one year. The fractional frequency uncertainty of the iodine frequency, assumed at 10 -15 , translates into a fractional uncertainty of the time derivative estimated at 4 × 10 -16 yr -1 . For this choice of clock comparisons, the uncertainties of the time derivatives of the fundamental constants are :

( ) ( ) ( ) dt d K K dt dµ µ K dt f f d f f Cs I I µ Cs I Cs I α α - + + = α α 1 1 1 1
1 14 1 14 10 8 1 1 10 6 1 1 - - - - × =       α α × =         yr . dt d u yr . dt dµ µ u (4)
The proposed approach will therefore enable a threefold improvement of the constraint on the variation of µ that has been determined by absolute frequency measurements of a rovibrational transition [START_REF] Shelkovnikov | Stability of the proton-to-electron mass ratio[END_REF]. Moreover, the molecular iodine frequency measurements can allow constraining for the first time the time variation of α with a molecular clock, independently on previous optical atomic clock comparisons.

V. CONCLUSION

This contribution proposes a new method to probe the variability of fundamental constants with precision measurements in molecular physics. The optical clocks based on the molecular iodine transitions are sensitive to the variations of the proton-to-electron mass ratio and the fine structure constant and may allow 10 -15 reproducibility. The absolute frequency measurements for a set of molecular iodine transitions may constrain conjointly the time variations of the proton-to-electron mass ratio and the fine structure constant at the 10 -14 yr -1 level.

Fig. 1 .

 1 Fig.1. Sensitivity coefficients to a variation of µ from selected optical transitions of 127 I2.

Fig. 2 .

 2 Fig. 2. An experimental setup for the molecular iodine clock. AOM acousto-optical modulator, EOM electro-optical modulator, x3 frequency tripling stage based on waveguide PPLN crystals.

  by introducing the sensitivity coefficients of the Cs clock to a variation of α, last equation neglects the contributions arising from the time variation of the relevant nuclear g-factors.

TABLE I

 I

	.	MOLECULAR IODINE TRANSITIONS SELECTED FOR
		ABSOLUTE FREQUENCY MEASUREMENTS
	127 I2 clock transition	Exp. freq. (THz)	Kµ	Kα
	127 I2 (1) : R(26) 62-0	597.366	-0.019	2
	127 I2 (2) : R(35) 44-0	583.254	-0.045	2
	127 I2 (3) : R(56) 32-0	563.260	-0.054	2

ACKNOWLEDGEMENT

This work was supported by the LABEX Cluster of Excellence FIRST-TF (ANR-10-LABX-48-01) within the Program « Investissements d'Avenir » operated by the French National Research Agency. Laboratoire PhLAM of the University of Lille is UMR 8523 of the CNRS.