K Bhattacharya

P M Suquet

A model problem concerning recoverable strains of shape-memory polycrystals

This paper addresses a model problem of nonlinear homogenization motivated by the study of the shape-memory effect in polycrystalline media. Specifically, it numerically computes the set of recoverable strains in a polycrystal given the set of recoverable strains of a single crystal in the two dimensional scalar (or antiplane shear) setting. This problem shares a direct analogy with crystal plasticity. The paper considers typical or random polycrystals where the grains are generated by a Voronoi tesselation of random points and are randomly oriented. The numerical results show that for such microstructures, the Taylor bound appears to be the most accurate (though pessimistic) bound when the anisotropy is moderate and that recent Kohn-Little-Goldsztein outer bounds overestimate the recoverable strains when the anisotropy is large. The results also show that the stress tends to localize on tortuous paths that traverse (poorly oriented) grains as the polycrystal reaches its limit of recoverable strain.

Introduction

This paper addresses a model problem of nonlinear homogenization in the plane motivated by the shape-memory effect in polycrystals. The shape-memory effect is a phenomenon observed in some metallic alloys wherein apparently permanent deformations suffered below a critical temperature can be recovered by heating above the critical temperatures. However, the extent of strain that a material can recover is limited, and subjecting the material to larger strains can lead to irreversible deformations. So the extent of the recoverable strain is an important attribute in evaluating shape-memory alloys. Interestingly, there are materials that are good shape-memory alloys in single crystals but lose all recoverable strain as polycrystals, and there are other materials whose polycrystalline shape-memory behavior is comparable to their single-crystalline behavior. Understanding the relation between single-crystalline and polycrystalline shape-memory behavior motivates this work.

The basis of the shape-memory effect is a martensitic phase transformation, a diffusionless solid to solid phase transformation involving a change of crystal structure. In shape-memory alloys, the high temperature austenite phase has greater symmetry than the low temperature martensite phase; so the martensite phase can exist in more than one symmetry-related variants. In fact, the different variants co-exist by forming fine-scale coherent microstructures. When subjected to some deformation, the material accommodates it to the extent it can by coherently rearranging the microstructure of martensitic variants. Such deformations are recoverable by heating above the transformation temperature to since every variant returns to the unique austenite phase. Larger deformations are perforce accommodated by slip and fracture, and are unrecoverable. Thus, the recoverable strains -those that can be accommodated by the coherent rearrangement of martensitic variants -are dictated by crystallography. However, most commercial materials are polycrystals, or agglomerations of a large number of grains of the same material whose (crystal lattice) orientations and consequently set of recoverable strains differ from one another. Thus, the recoverable strains of a polycrystal are determined by the collective response of a large number of grains, and we seek to understand this. Bhattacharya and Kohn (5) modelled this problem using energy minimization. Since the martensite comes in multiple variants, the free energy as a function of strain has multiple wells with one well for each variant. With suitable normalization, the energy is zero when the strain is equal to the transformation or stress-free strain of the variant and grows away from it. Energy minimization with such multi-well or nonconvex energies leads to the formation of microstructure. Thus, at the level of a grain, one sees an effective energy of the single crystal which is the energy of a grain after it has formed suitable microstructure. Mathematically, this is obtained from the original multiwell energy by relaxation. This relaxed energy is zero on the set of recoverable strains of a single crystal -the set of strains that a single crystal can accommodate by the formation of an (almost) stress-free microstructure of martensitic variants -and grows away from it. In a polycrystal, each grain has its own effective energy and its own set of recoverable strains. Therefore, its behavior is determined by the effective energy of a polycrystal which is the energy of a polycrystalline specimen after taking into account all the interactions between the grains. Mathematically, it is obtained from the underlying energy of the single crystal by homogenization. The effective energy of a polycrystal is zero on the set of recoverable strains of the polycrystal -the set of strains that a polycrystal can accommodate by forming an (almost) stress-free and possibly spatially varying microstructure of martensitic variants in each grain -and grows away from it.

In the setting of infinitesimal strain theory, one can calculate the set of recoverable strains of a single crystal in various examples of interest and estimate them in others (3; 4). However, calculating the set of recoverable strains in a polycrystal remains an outstanding problem. Bhattacharya and Kohn (4; 5) proceeded using the so-called Taylor bound which neglects cooperative effects between grains and found that the recoverable strains depend on the transformation strains, the texture (the shapes, the orientations and arrangements of the grains) and critically on the change of symmetry during transformation (4). This allowed them to explain various observations. They also discussed the accuracy of the Taylor bound in some model problems, but this remains an open issue in general. This paper addresses the numerical computation of the set of recoverable strains for a polycrystal given the set of recoverable strains for single crystal in the two dimensional scalar (or antiplane shear) setting. In this setting, relaxation is equivalent to convexification and thus the effective energy of a single crystal is readily obtained. So we start by directly writing down the effective energy of a single crystal. In fact, we consider "locking materials", where the energy is zero on a given set of recoverable strains of a single crystal and infinite otherwise. We then investigate numerically the effective domain of recoverable strains for polycrystals. Chenchiah and Bhattacharya (11) considered a similar model in the setting of plane strain.

This problem is analogous to the problem of rigid-plastic polycrystals where

Article submitted to Royal Society the stress is constrained to lie in a given set in each grain. Indeed, the analogy is rigorous in the current scalar two-dimensional setting since divergence-free (stress) fields can be obtained by a rotation of curl-free (strain) fields. Thus one can draw on that literature. Importantly, the homogenization theory of Bouchitte ( 9) and Bouchitte and Suquet [START_REF] Bouchitte | Homogenization, Plasticity and Yield design[END_REF] is directly applicable to the problem at hand. One can also recall the various bounds and estimates for the effective properties. deBotton and Ponte Castañeda [START_REF] Debotton | Variational estimates for the creep behavior of polycrystals[END_REF] derived an upper bound using Ponte Castañeda's variational method (the linear comparison method), but this was found to be close to the classical Taylor bound for rigid-plastic polycrystals. It was improved upon by the self-consistent estimate of Ponte Castañeda and Nebozhyn [START_REF] Castañeda | Variational estimates of the selfconsistent type for the effective behaviour of some model nonlinear polycrystals[END_REF], also based on the same variational principle. The problem of polycrystals made of grains with a deficient number of slip systems motivated Kohn and Little [START_REF] Kohn | Some model problems of polycrystal plasticity with deficient basic crystals[END_REF] to consider a model problem in the scalar two-dimensional setting where each grain has two orthogonal slip systems with critical resolved shear stresses 1 and M . They showed that the yield strength of isotropic polycrystal made of such grains is bounded by √ 2M . The improvements by Nesi, Smyshlaev and Willis [START_REF] Nesi | Improved bounds for the yield stress of a model polycrystalline material[END_REF] and Goldsztein (17) also scale as √ M for large M . In fact, Goldsztein [START_REF] Goldsztein | Two-dimensional rigid polycrystals whose grains have one ductile direction[END_REF] has shown that his (4/π) √ M bound is sharp even under the assumption that the texture has square symmetry.

This paper considers both "model" microstructures like checkerboards and "typical" random microstructures close to the ones that one would expect in actual materials. Specifically, the grains in the latter are obtained by a Voronoi tesselation of the plane and their crystallographic orientations are chosen randomly. These polycrystalline microstructures are asymptotically (in the limit of an infinitely large volume element) isotropic and statistically homogeneous. They form a class of microstructures which is narrower than the entire class of isotropic microstructures.

The local fields in and the average or effective response of these polycrystalline aggregates are investigated numerically. The computational method used here is based on the fast Fourier transform and is an extension of that developed in Plasticity by Moulinec and Suquet (24) as well as Michel, et al. [START_REF] Michel | A computational method for linear and nonlinear composites with arbitrary phase contrast[END_REF]. Similar methods have been used to analyse diffusional phase transformations in solids by Dreyer and Müller (14) and viscous polycrystals by Lebensohn (21).

Our numerical approach has similarities with the recent studies of shape-memory polycrystals using the phase field models by Artemev et al. [START_REF] Artemev | Three dimensional phase field model and simulation of cubit to tetragonal martensitic transformation in polycrytals[END_REF] and Ahluwalia et al.

(2) amongst others. There are also important differences. First, we use the relaxed energy and they use the unrelaxed energy with a higher gradient regularization. So they have to explicitly resolve the internal twinning while we implicitly account for them. Second, we solve the equilibrium equation while they study the evolution. Finally, Ahluwalia et al. (2) also allow the orientation of the grains to evolve.

Our computational results lead to two interesting conclusions. The first concerns the accuracy of existing bounds and estimates. For weakly anisotropic single crystals, the computed effective set of recoverable strains for "random" polycrystals is found to be closest to the Taylor (inner or pessimistic) bound and significantly smaller than most available (outer) bounds. For large anisotropy ratios M of the single crystals, the computational results indicate that the effective recoverable strains in the polycrystal scale as M α , with α of the order of 1/4 whereas the scaling predicted by most bounds is α = 1/2. This does not contradict the results of Goldsztein [START_REF] Goldsztein | Two-dimensional rigid polycrystals whose grains have one ductile direction[END_REF] because that author considers a broader class of microstructures.

Second concerns the nature of the stress fields. Although the variable of primary interest is the overall strain (and the effective set of recoverable strains), the stress field plays an important role at the microscopic level. When the overall strain is on the boundary of the set of recoverable strain, the components of the microscopic stress field are likely to be bounded measures. In other words the stress accumulate along certain lines travelling through the poorly oriented grains in the polycrystal and carrying the force from one side of the polycrystal to the opposite side.

Energy minimization

(a) Single crystals and polycrystals

In the setting of anti-plane shear, the displacement u is a scalar field, the strain ε is a curl-free vector field, the stress σ is a divergence-free vector field,

ε = ∂u ∂x 1 , ∂u ∂x 2 , curl(ε) = 0, σ = (σ 1 , σ 2 ) , div (σ) = 0. (2.1)
The single crystal behavior is characterized by a convex strain energy w s (ε) (s stands for single crystal) or by its Legendre transform w * s (σ),

σ = ∂w s ∂ε (ε), ε = ∂w * s ∂σ (σ), w s (ε) + w * s (σ) = σ : ε. (2.2)
A typical representative volume element V of the polycrystal is comprised of domains or grains, V r , with characteristic functions χ (r) (x), which achieve a partition of V . All grains are perfectly bonded across their boundaries. The grains' shape differ but their constitutive relations, or strain energy, are identical to that of a unique single crystal up to a rotation by an angle θ (r) (for grain V r ). Therefore the strain energy at point x in the polycrystal is given by w(x, ε) = w(θ(x), ε), θ(x) = N r=1 θ (r) χ (r) (x), w(θ, ε) = w s (R T (θ).ε), (2.3) where R(θ) is the rotation by an angle θ in the plane and R T is its transpose.

For periodic polycrystals considered here, the average strain ε and the average stress σ are related by

σ = ∂ w ∂ε (ε), or equivalently ε = ∂ w * ∂σ (σ). (2.4)
The effective strain energy w and the complementary energy w * of the polycrystal are characterized by dual average variational principles :

w(ε) = inf u ∈ K(ε) w(x, ε(u)) , and w * (σ) = inf τ ∈ S(σ) w * (x, τ ) , (2.5)
where

K(ε) = {u such that u = ε.x + v, v # on ∂V }, (2.6)
and

S(σ) = {τ such that div (τ ) = 0 in V, τ = σ, τ .n -# on ∂V }. (2.7)
The symbols # and -# denote periodic and anti-periodic fields respectively.

Article submitted to Royal Society

(b) Linear elasticity

In linear elasticity, the strain energy w s and the complementary energy w * s read 2) .

w s (ε) = 1 2 µ (1) ε 2 1 + µ (2) ε 2 2 , w * (σ) = 1 2 σ 2 1 µ (1) + σ 2 2 µ (
(2.8)

µ (1) and µ (2) are the two anti-plane shear moduli of the single crystal. The ratio M = µ (2) /µ (1) measures the anisotropy of the single crystal, and M = 1 corresponds to an isotropic crystal.

(c) Domain of recoverable strains

The domain of recoverable strains in a single crystal of a shape-memory material is a closed convex set K s . The energy density of the single crystal is taken to be its indicator function I Ks :

w s (ε) = 0 if ε ∈ K s , w s (ε) = +∞ otherwise.
(2.9)

The corresponding dual function w * s is positively homogeneous of degree one in σ,

w * s (σ) = sup ε ∈ K s σ.ε. (2.10)
When the grain is rotated by an angle θ the domain of recoverable strains is rotated by the same angle,

w(x, ε) = I K(x) (ε) where K(x) = R(x).K s = {ε, R T (x)ε ∈ K s }. (2.11)
A typical form for K s is a rectangle corresponding to the four variant case in the terminology of Bhattacharya and Kohn ( 5)

K s = [-ε (1) 0 , +ε (1) 0 ] × [-ε (2) 0 , +ε (2) 0 ].
(2.12)

The aspect ratio M = ε

(1)

0 /ε (2)
0 is a measure of the anisotropy of the single crystal. Note that, unlike in the linear case where the energy is given by (2.8) , a nonlinear crystal with energy (2.9) is anisotropic even when M = 1, since the domain K s is a square (and not a ball). The corresponding complementary energy is given by

w * s (σ) = ε (2) 0 (M |σ 1 | + |σ 2 |) . (2.13)
Since the stress energy w * s (x, .) of all grains is positively homogeneous of degree 1, the effective stress energy w * inherits this feature. By duality the effective strain energy w is the indicator function of a closed convex set K which is the effective domain of recoverable strains for the polycrystal:

w(ε) = I e K (ε) where K = {ε, ε.σ ≤ w * (σ) ∀ σ ∈ R 2 }. (2.14)
The first variational problem in (2.5) provides an alternative, but rigorously equivalent (see section e), characterization of K. K is the set of macroscopic strains which are averages of microscopic strain fields lying in K(x),

K = {ε, ∃u ∈ K(ε), ε(u(x)) ∈ K(x) ∀x ∈ V, ε = ε }. (2.15)
The physical interpretation of this result is that a strain ε is recoverable at the macroscopic level whenever it is possible to find a local strain field which is recoverable at every point x and which averages to ε.

(d ) Connection with locking materials. Stress fields

The problem considered here, where the strain is assumed to remain in a closed convex set, is exactly that of locking materials. Demengel and Suquet (13) provide a mathematical discussion of variational problems associated with locking materials. It is worth recalling two features of the stress fields from their work, since they have obvious implications on the choice of numerical method and the interpretation of the results. First, the components of the stress field are expected to be bounded measures. Although it is not possible to handle measures numerically, the numerical method should be able to deal with the concentration of stress in thin bands approaching lines. Second, there is no uniqueness in the stress field.

(e) Relation with plasticity. Homogenization

The problem of recoverable strains shares several features with those of plasticity. In the theory of rigid-plastic materials, the stress σ(x) is constrained to remain in a closed convex domain P , and the stress field is subject to equilibrium equations and boundary conditions. The effective strength domain P is the set of average stresses σ for which a local stress field σ can be found which meets all the above requirements and whose average is σ (Bouchitté and Suquet [START_REF] Bouchitte | Homogenization, Plasticity and Yield design[END_REF]). This problem is equivalent to the problem at hand, upon exchange of stress and strains, of equilibrium equations and compatibility equations, of domains P and K. Further, in plasticity the space of strains with finite energy is a space of measures and the corresponding displacements lie in BD(Ω) (Temam [START_REF] Temam | Mathematical problems in Plasticity[END_REF] and Suquet ( 29)), and analogously, the stress fields are unbounded in the problem of recoverable strains.

Indeed, in the two-dimensional scalar setting considered here, the problems of recoverable strains and of plasticity are rigorously equivalent. This follows from the well-known observation of Dykhne ((15)) according to which a divergence-free vector field can be transformed into a gradient vector field (and conversely) by means of an appropriate rotation. Let R ⊥ denote the rotation by π/2 in the plane, R ⊥ . (f 1 , f 2 ) = (-f 2 , f 1 ). Whenever σ is a stress field in S(σ) then R ⊥ .σ is the gradient of a displacement field in K(R ⊥ .σ). Conversely, whenever ε derives from a displacement field in K(ε) then R ⊥ .ε is in S(R ⊥ .ε). Therefore the second variational problem (2.5) is strictly equivalent to the problem of finding the effective energy (or dissipation potential when ε is interpreted as strain-rate) for a rigid-plastic polycrystal. Consequently it is possible draw on the existing homogenization theorems of Bouchitté [START_REF] Bouchitte | Convergence et relaxation de fonctionnelles du calcul des variations à croissance linéaire. Application à l'homogénéisation en Plasticité[END_REF] and Bouchitté and Suquet [START_REF] Bouchitte | Homogenization, Plasticity and Yield design[END_REF] in the setting of rigid-perfect plasticity to study the problems (2.5). This theory alerts us to two technical difficulties. First, in plasticity the strains can be measures and therefore the stress can be measures as discussed earlier. Second, although the energy is convex, part of it pertaining to the boundary conditions may not be lower semi-continuous in the topology that is natural for Γ-convergence. Therefore the computation of the Γ-limit involves a relaxation of non lower semi-continuous terms.

While the analogy between recoverable strain and plasticity problems is exact in the two-dimensional scalar setting, it is not exact for either the three-dimensional scalar or the multi-dimensional vector setting, and one will have to revisit the homogenization theorems following Bouchitté and Suquet [START_REF] Bouchitte | Homogenization, Plasticity and Yield design[END_REF].

Bounds and estimates (a) Taylor and Sachs bounds

Taking uniform fields ε = ε and σ = σ in the variational principles (2.5) yields the Taylor and Sachs energy bounds respectively:

Taylor : w(ε) ≤ w T (ε) = w (ε), Sachs : w * (σ) ≤ w * S (σ) = w * (σ). (3.1)
When w s is the indicator function of K, the Taylor and Sachs energy bounds are also indicator functions of the sets K T and K S respectively. Further, they satisfy

K S ⊃ K ⊃ K T = x∈V K(x). (3.2)
So we refer to these sets K T and K S as the Taylor and Sachs bounds respectively. They depend in general on the texture of the polycrystal. In an equiaxed polycrystals (where all grain orientations occur with equal probability), the Taylor and Sachs bounds are balls with the following radii:

Taylor :

ε T 0 = min(ε (1) 0 , ε (2)
0 ), Sachs :

ε S 0 = 2 π ε (1) 0 + ε (2) 0 . (3.3) (b)
Bounds based on a linear comparison polycrystal deBotton and Ponte Castañeda [START_REF] Debotton | Variational estimates for the creep behavior of polycrystals[END_REF] have derived bounds for the effective properties of nonlinear polycrystals which are relevant to the present problem. We provide an alternate and simple derivation of their result. Consider first a single grain and note that for all symmetric positive definite second-order tensors M , σ.ε ≤ (σ.M .σ) 1/2 (ε.L.ε) 1/2 , with L = (M ) -1 .

(3.4)

Then it follows from (2.10) that,

w * s (σ) ≤ (σ.M .σ) 1/2 Sup ε ∈ K s (ε.L.ε) 1/2 . (3.5)
Choosing M in diagonal form, with eigenvalues 1/2µ 1 and 1/2µ 2 , yields:

w * s (σ) ≤ 1 2µ 1 σ 2 1 + 1 2µ 2 σ 2 2 1/2 2µ 1 (ε (1) 0 ) 2 + 2µ 2 (ε (2) 0 ) 2 1/2 . ( 3.6)
Similar inequalities can be obtained for all individual grains V r with different elastic moduli M (r) (or different shear moduli µ (r)

1 and µ (r)
2 ) in each grain. After averaging over all grains (and use of the Cauchy-Schwarz inequality), an upper bound is obtained for the effective potential w * s :

w * (σ) ≤ Inf µ (r) 1 , µ (r) 2 ≥ 0, r = 1, ..., N ( w * 0 (σ)) 1/2 2µ (r) 1 (ε (1) 0 ) 2 + 2µ (r) 2 (ε (2) 0 ) 2 1/2 , (3.7)
where w * 0 (σ) is the effective energy of a linear comparison composite (LCC) with the same geometry (same grains) as the actual nonlinear polycrystal but with linear Article submitted to Royal Society elastic behaviour characterized by shear moduli µ (r) 1 and µ (r) 2 in the r-th grain V r . Note that, since the shear moduli's differ from one grain to another, the different grains in the LCC are not obtained by rotation from the same single crystal.

Rigorous upper bounds for w * are obtained by bounding w * 0 from above. But only estimates for w * are obtained when estimates (such as the self-consistent scheme) are used to estimate the effective properties of the LCC.

The best bound in (3.7) is obtained by solving an optimization problem which involves 2 × N variables (the µ (r) i 's). This optimization problem can be quite large for actual polycrystals, even for untextured polycrystals. So a cruder bound can be obtained for untextured polycrystal by choosing the shear moduli in each grain to be identical (µ (r) 1 and µ (r) 2 independent of r). Then it can be proved that the optimal choice is µ 1 /µ 2 = 1/M 2 where M is the anisotropy ratio of the single crystal. The resulting bound is the Kohn-Little bound [START_REF] Kohn | Some model problems of polycrystal plasticity with deficient basic crystals[END_REF].

The variational bound (3.7) has been improved by Nesi et al (26) using the translation method. It is similar to the Taylor bound for small M and the Kohn-Little bound for large M .

Computational method

The Euler-Lagrange equations corresponding to the variational problems (2.5) read

ε(u(x)) = ∂w * ∂σ (x, σ(x)), div (σ(x)) = 0 ∀x ∈ V, σ = σ, u -ε(u) .x #, σ • n -# (4.1)
The computational method used to solve (4.1) follows closely that developed by Moulinec and Suquet [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] for plasticity problems, except that it is written here in dual form. It consists in re-writing (4.1) as a nonlinear integral equation which is then solved by successive iterations.

(a) Reduction to a nonlinear integral equation

Introducing a homogeneous linear elastic reference material with compliance M 0 , the local problem (4.1) can be re-written as

ε(u(x)) = M 0 • σ(x) + η(x), div (σ(x)) = 0 ∀x ∈ V, σ = σ, u -ε(u) .x #, σ • n -#. (4.2)
The macroscopic stress σ is prescribed, the macroscopic strain ε = ε is unknown and the eigenstrain η(x) reads

η(x) = ∂w * ∂σ (x, σ(x)) -M 0 • σ(x). (4.3)
The solution σ of (4.2) can be expressed as σ(x) = σ -∆ 0 * η(x) where ∆ 0 is the periodic Green's operator associated with the reference medium M 0 and can be easily expressed in Fourier space ( f denotes the Fourier transform of f ):

∀ξ = 0 : σ(ξ) = -∆ 0 (ξ) • η(ξ), σ(0) = σ, where ∆ 0 (ξ) = (L 0 .ξ) ⊗ (L 0 .ξ) ξ.L 0 .ξ -L 0 , L 0 = (M 0 ) -1 .      (4.4)
Article submitted to Royal Society Substituting the expression (4.3) of η into (4.4) leads to the following nonlinear integral equation for σ:

in real space σ(x) = -∆ 0 (x) * (δM (x, σ(x))) + σ, in Fourier space σ(ξ) = -∆ 0 (ξ) • δM (σ)(ξ) ∀ξ = 0, σ(0) = σ, (4.5)
where

δM (x, σ) = ∂w * ∂σ (x, σ) -M 0 • σ. (4.6) (b) Iterative algorithm
The nonlinear integral equation (4.5) can be solved iteratively by a fixed-point algorithm similar to the one used in Moulinec and Suquet [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]:

real space: σ i+1 = σ -∆ 0 * δM (σ i ) , Fourier space: σ i+1 (ξ) = -∆ 0 (ξ) • δM (σ i )(ξ) ∀ξ = 0, σ i+1 (0) = σ (4.7)
This fixed point algorithm is implemented in Fourier space. However δM (σ i+1 )(ξ) is obtained via the computation of (σ i+1 ) in real space. This requires the use of a direct and an inverse Fourier transform to pass from one space to the other. A further simplification is obtained by noting that ∆ 0 * M 0 • τ = ττ , for all divergence-free field τ and therefore -∆ 0 * δM (x,

σ i ) = σ i -∆ 0 * ∂w * ∂σ (x, σ i ).
The final form of the algorithm is, at iterate i+1 :

σ i (x) and ε i (x) being known, perform : a) εi = F(ε i ), b) convergence test, c) σi+1 (ξ) = σi (ξ) - ∆0 (ξ) : εi (ξ) ∀ξ = 0, σi+1 (0) = σ, d) σ i+1 = F -1 ( σi+1 ) e) ε i+1 (x) = ∂w * ∂σ (x, σ i+1 (x)), ∀ x,                    (4.8)
where F and F -1 denote the Fourier transform and the inverse Fourier transform. Iterations are initialized with σ 0 (x) = σ.

(c) Space discretization.

The unit-cell is discretized into pixels. The Fourier transform is replaced by the discrete Fourier Transform, which can be computed using the Fast Fourier Transform (see [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] for details). The unknowns u, ε and σ are sampled by taking their values at these discrete pixels or voxels. The spatial discretization of the image induces a corresponding spatial discretization for the frequency ξ in Fourier space.

The grid size (number of pixels) has a significant influence on the local fields (also on the effective properties but to a lesser extent) since, as will be seen later, the stress fields might have a very fine microscopic structure. Convergence of the fields is not ensured even with a fine mesh, since stress fields can accumulate on lines with zero width and are not unique. In most of our computations, grids with a resolution of 2048×2048 pixels were used and for some of them the resolution was even increased to 4096×4096.

(d ) Convergence test.

Convergence is reached when the field ε i+1 is a compatible strain field. In the current setting, curl(ε) = ∂ε 1 ∂x 2 -∂ε 2 ∂x 1 provides a measure of the incompatibility of ε. The iterative algorithm is stopped when curl(ε i+1 ) / ε is smaller than a tolerance threshold. The accuracy of the strain and stress fields depend significantly on the tolerance parameter which should be small. In our computations it was systematically chosen equal to 10 -6 . Larger values gave much faster convergence but led in certain circumstances to significant error in the local fields.

(e) Miscellaneous Remarks 1. In the problem of recoverable strains, the energy w * s is singular (non differentiable) at the origin. A regularization is introduced by adding a small (artifical or real) elasticity at the origin. For simplicity this elastic term was chosen to be isotropic and the regularized energies corresponding to (2.9) and (2.13) read as:

w s (ε) = 1 2η ε 2 1 + ε 2 2 + I Ks (ε). (4.9)
The corresponding dual energy w * s is differentiable. In our computations the parameter η was taken equal to 10 -2 ε

(2) 0 . With the regularization (4.9) the grains are elastically homogeneous and isotropic. The homogeneous reference medium is also taken isotropic and characterized by the same elastic parameter η.

2. All our results are obtained by solving the equations (4.1) under a prescribed macroscopic stress or load. Since one may have high contrast in local properties, the fixed point algorithm (4.8) may have slow or no convergence if this load is applied in one step. Therefore, we applied the load incrementally, i.e. by increasing the macroscopic stress gradually along the direction of prescribed stress. We used approximately 3000 steps for each applied load and need 2 or 3 iterations for most time steps to reach convergence (with a maximum of 20 iterations for a very few steps). We had to exercise further care since large increments may result in significant oscillations of the stress-strain relations once the boundary of the domain of recoverable strains is attained. An alternative would have been to follow the improved algorithm of Michel et al [START_REF] Michel | A computational method for linear and nonlinear composites with arbitrary phase contrast[END_REF].

3. As the loading is increased incrementally, the average strain is computed at each time-step to monitor when it reaches a limit which lies on the boundary of K. All stress and strain snapshots shown correspond to these limit strain states. We also note that the macroscopic stress is an outer normal vector to K at this point.

Test examples: checkerboards

We begin with model polycrystals in the form of checkerboards. They provide a test of the computational method since much is known analytically about them. the checkerboard and one rotated by 45 0 (see Figure 1(a)), then the Taylor bound is the intersection of the original square domain with the same square rotated by 45 0 (see Figure 1(b)). It is also known ( 5) that the Taylor bound is in fact equal to the actual domain of recoverable strains in this example. The domain K is determined computationally in the space of macroscopic strains (ε 1 , ε 2 ) by applying a macroscopic stress in different directions (every 15 0 ) :

σ = σ (cos α e 1 + sin α e 2 ) .
(5.1)

The computed limit strains are shown in Figure 1(b) as crosses, and the corresponding directions of stress are shown by arrows. All computed points lie (with an error less than 2%) on the exact surface (solid line). A few of them are concentrated at the vertices, even though they correspond to different macroscopic stresses. At these points, K has an outer normal cone instead of a single outer normal vector. Snapshots of the intensity or magnitude of the stress field (defined as σ 2 1 + σ 2 2 1/2 ) at the strain limit are shown in Figure 2. The stress localises on thin lines, consisting of only one row of pixels. When the stress inclination is 0 0 or 45 0 , the inclination of the lines where the stress localise are 0 0 and 45 0 respectively. At other inclinations of the macroscopic stress (15 0 and 30 0 ) the stress localises on these two same lines, but with different weights for each line to balance the applied average stress. This "stress localisation' is the numerical manifestation of the fact that the mathematical solutions for the stress can be measures.

(b) Fine scale solutions in highly anisotropic checkerboards We now provide an example where the stress is chanelled into bands with a multiscale structure. This is motivated by the work of Dykne and Kaganova [START_REF] Dykne | The electrodynamics of polycrystals[END_REF] and Bhattacharya, Kohn and Kozlov [START_REF] Bhattacharya | Some examples of nonlinear homogenization involving nearly degenerate energies[END_REF] who observed in the linear setting that when the single crystals are sufficiently anisotropic, the stress can even be "trapped" at multiple points (intersection of grain boundaries).

The unit cell is a checkerboard with individual crystals inclined at θ and π/2 + θ to the horizontal axis (Figure 3 with aspect ratio M . The length of the arrows in the figure correspond (schematically) to the size of the recoverable strain in the principal directions of the crystals.

The computed stress field in the limit strain state for θ = π/6, M = 100, computational resolution of 2048 × 2048 pixels and the applied stress parallel to {1, 1} is shown in Figure 3(b). It follows bands which are aligned along the directions of the smallest possible recoverable strain. The crystals "lock-up" along these directions and the stress builds up in reaction. This stress field has apparently two arms, one winding towards the center of the cell and another one unwinding from this center.

Despite the refinement of the grid, it was impossible to decide whether the stress fields spirals indefinitely or with a finite number of arms. A theoretical argument is helpful to answer, at least qualitatively, this question. A stress field constructed by Bhattacharya et al (6) (called flow and denoted by g in ( 6) ) can be used to derive an upper bound for the energy w * . It consists of a spiral stress flow winding in towards the center, crossing a bridge and a second spiral stress flow winding away from the center of the cell. The stress is uniform in each arm of the spiral, with intensity σ n in the n-th arm and parallel to the direction of the minimal recoverable strain. Its complementary energy in the n-th arm is w obtain the following bound :

w * (σ) ≤ σ 3 -t + t N - 1 2 (1 + t -2 )(1 + t) + M t N -4 (1 -t 3 )(1 -t) (cos θ) 2 (5.2)
where t = tan θ and N is the number of arms in the inflowing spiral. The right-hand side of (5.2) can be schematically written as C -at N + M t 2N where a > 0. The minimum of this function (as a function of N ) is attained when t N = a/2M , that is for a finite value of N . Therefore there is a finite number of arms but this number becomes infinite in the limit as M goes to +∞. Another question of interest is to determine whether the recoverable strain in the diagonal direction remains finite even when the domain of recoverable strains in the crystals become infinite in one direction, i.e., in the limit when M tends to +∞. The trial field described above has average σ = 1 2 σ(cos θ -sin θ)(1, 1). Since every ε in K should satisfy ε.σ ≤ w * (σ), one has :

ε 1 +ε 2 ≤ 1 cos θ -sin θ 3 -t + t N - 1 2 (1 + t -2 )(1 + t) + M t N -4 (1 -t 3 )(1 -t) cos θ 2 .
Taking the limit of this expression as N tends to +∞ gives an upper bound for the strain in the diagonal direction which is independent of M :

√ 2 2 (ε 1 + ε 2 ) ≤ √ 2 2 3 -t cos θ -sin θ .
This estimate shows that the recoverable strain in the diagonal direction is uniformly bounded, even as the recoverable strain of one direction of the single crystal becomes unbounded.

Random polycrystals

We now turn to typical or random polycrystals. By "random polycrystals" we mean random statistically homogeneous aggregates of perfectly bonded single-crystals with different orientations. We focus on untextured polycrystals where no particular orientation is priviledged. (a) About the number and the shape of the grains It is important to carefully examine the shape, orientation and number of grains in a unit cell to ensure a good approximation to a statistically homogeneous structure (see Kanit et al. (19) for similar considerations in linear problems). In this paper, two different classes of polycrystalline microstructures are considered, those with identical hexagonal grains (Figure 4(a)) and those with irregular polygonal grains obtained by a Voronoi tesselation of a set of random points (Figure 4(b)). In both cases the orientation of each individual grain is generated randomly with no priviledged orientation (untextured polycrystals). We consider 10 different realizations in each case containing a large number of grains. The microstructures with hexagonal grains contained 1024 hexagons and the computational spatial resolution was 2048 × 2048 pixels. This corresponds to a spatial resolution of 64 × 64 pixels per hexagon which turned out to be sufficient to capture correctly the high stress gradients inside the individual grains. The same overall computational spatial resolution of 2048 × 2048 pixels was used for the Voronoi tesselations but with only 512 grains in each unit-cell (to allow a sufficient resolution for the smallest grains).

Article submitted to Royal Society

An alternative approach that considered a large number of realizations with a small number of grains (typically 100 realizations with 16 grains each) was also explored. The ensemble averages of the results were found to be in agreement with the results obtained with fewer large configurations. All results presented below are with large numbers of grains.

The overall isotropy of the different configurations was checked by applying a stress (5.1) along different orientations and plotting the norm of the limit strain . Two points are worth noting. First the two sets of results for the two grain shapes show almost no difference. The shape of the grains does not seem to have an influence on the effective properties, at least for isotropic textures considered here. Second, the limit strain does not depend on the orientation of the applied stress, up to a few percent. We conclude that the domain of recoverable strains can reasonably be considered isotropic. The bounds and estimates which are also shown in this figure will be discussed later.

ε 0 = ε 2 1 + ε 2
(b) Structure of the stress and strain fields in polycrystals

The stress localisation observed in the checkerboards persists in more complex microstructures as shown in Figure 4(d). This is explored further in the 16 grain microstructure shown in Figure 5(a). The magnitude of the stress at the limit strain state for an applied macroscopic stress (5.1) with different α and M = 1 are shown in Figure 5(b,c). Most of the stress follows lines running across each crystal with sufficiently high magnitude so that they can be regarded as measures. These lines are roughly straight and parallel to the directions of lowest recoverable strain (the two principal axes of the crystals when M = 1) in each grain. This observation is consistent with the minimization of M |σ 1 |+|σ 2 | (with M = 1 in the present case) in the individual grains. The average orientation of the lines is imposed by the average stress through the condition σ = σ. The stress enters the computational unit-cell at A and is transferred to the adjacent unit cell at point A ′ in all snapshots.

For the angle α = 30 • , the apparently different lines in fact correspond to the same line extending over several unit-cells: the stress (shown by an arrow) enters the unit-cell at A, is transferred to the neighbouring cell at B, which is equivalent (modulo periodicity) to B ′ , then is transferred to a neighbouring unit-cell at C which is equivalent to C ′ and finally exits the cell at A ′ . Similar effects were also observed at α = 15 • and α = 45 • . This long-range effect -the fact that the stress lines run several unit cells before closing -is closely connected to the fact that our complementary energy w * s has a linear growth at infinity and is on the border between the convex and the nonconvex case. In the nonconvex case, it is known (7; 25) that the effective energy is defined by the minimum of the average energy, not only among admissible stress fields which are periodic on a single unit-cell but also among admissible stress fields which are periodic over multiples of the unit cell. The problem at hand is convex but not strictly convex and we find a hint of the features of nonconvex problems. 8)). The strain field is more homogeneous than the stress fields, but still not fully homogeneous. The strain magnitude ranges from 0.48 to 1.414 while the stress magnitude ranges from 0 to over 100. However, one has to be careful about the interpretation. Since the set K is not isotropic, a strain magnitude of 1 can indicate locking if only one component is non-zero while a strain magnitude as large as 1.414 may indicate no locking when both components are equal. Further, the connection to the Taylor bound is also tricky. While uniform field gives the Taylor bound, the converse -the attainment of the Taylor bound implies uniform field -is not true.

(c) Effective properties

(i) Effective modulus of linear elastic polycrystals

We begin by studying linear elastic polycrystal where the energy of a single crystal is given by (2.8) and the anisotropy ratio is given as M = µ (2) /µ (1) . It can be shown using arguments like those of Dykne (15) that the the effective shear modulus of an isotropic polycrystal is µ = µ (1) √ M . The self-consistent estimate reproduces

Article submitted to Royal Society this result. Computations were run on 10 different Voronoi tesselations for the anisotropy ratio M ranging from 1 to 10 3 . The reference medium was taken isotropic characterized by a single shear modulus µ 0 = 2M/(M + 1)µ (1) . Figure 6(a) shows that the agreement between the theory and computation is satisfactory as noted independently by Lebensohn et al. [START_REF] Lebensohn | Macroscopic properties and field fluctuations in model power-law polycrystals: full-field solutions versus self-consistent estimates[END_REF]. The scatter of the computational results is small for small M (indicating that statistical homogeneity is almost reached), but increases with the anisotropy ratio M (indicating that the deviations from statistical homogeneity increase with anisotropy). In particular the number of grains necessary to approach statistical homogeneity with a given precision depends on the anisotropy of the phases.

(ii) Effective domains of recoverable strains with small anisotropy, M = 1

We now turn to the problem of effective recoverable strains for non-textured polycrystals, when the single crystal has equal recoverable strains in both directions. A comparison between computational results and a few available bounds and estimates is presented in Figure 4 ≃ 1.21. The average of all computational results, obtained with regular hexagonal grains or with Voronoi tesselations is ε 0 /ε (2) 0 ≃ 1.08, and this is significantly smaller than the preditions of all upper bounds and estimate. The Taylor lower bound ε T 0 /ε

(2) 0 = 1, although a bit pessimistic, is indeed the closest bound/estimate to the computational results.

The fact that the variational bounds and estimate based on the linear comparison medium significantly overpredict the computational result despite the agreement observed earlier in the linear elastic polycrystals implies that the inequality (3.7) is far from being sharp. This is consistent with the fact that the stresses localise and the inequality based on taking the second moment of the stress over the grain leads to a significant loss of information.

(iii) Effective domains of recoverable strains with large anisotropy M ≫ 1.

Figure 6(b) shows a similar comparison between various bounds and estimates on the one hand and computational results on the other regarding the dependance of the effective recoverable strain ε 0 on the anisotropy ratio M . All bounds and estimates (except the Taylor and Sachs bounds) predict that the effective recoverable strain ε 0 increases as √ M for large M . Therefore, only the Kohn-Little bound is shown in the figure; in the log-log coordinates used in the figure all others have the same slope (1/2) but with slightly different initial values at M = 1. In contrast, the trend observed in the computations is a straight-line with a slope of 1/4 when M ranges from 1 to 100. The deviation from the average is significant for large M indicating that achieving statistical homogeneity for strongly anisotropic (and nonlinear) materials is difficult. Despite this observation, all computations results for all 10 Voronoi tesselations lie well below the Kohn-Little bound, and indicate that the scaling law for ε 0 is closer to M 1/4 which is significantly smaller that the scaling law M 1/2 predicted by most bounds and estimates. The Taylor bound, which is equal to 1 and independent of M , also misses the observed scaling. Interestingly, Lebensohn et al. [START_REF] Lebensohn | Macroscopic properties and field fluctuations in model power-law polycrystals: full-field solutions versus self-consistent estimates[END_REF], also using a FFT method, observed that most estimates of the self-consistent type overestimate the effective properties of viscous nonlinear polycrystals.

( a )Figure 1 .

 a1 Figure 1. Checkerboard. (a) Microstructure and orientation of the grains. (b) Effective domain of recoverable strains, numerical computations and exact results.

  (a)). The domain of recoverable strains is a rectangle Article submitted to Royal Society α = 0 0 α = 15 0 α = 30 0 α = 45 0

Figure 2 .

 2 Figure 2. Snapshots of the magnitude of the local stress field. Dark color corresponds to small stress and white corresponds to large stress.

Figure 3 .

 3 Figure 3. (a) Checkerboard with anisotropic grains. (b) Snapshot of the stress intensity; brigher (hotter) shades represents greater intensity.

Figure 4 .

 4 Figure 4. a) Hexagonal polycrystal and b) Typical Voronoi tesselation polycrystal. Different shades indicate different orientations. c) Dependence of the recoverable strain e ε0 as a function of the orientation of the applied stress (M = 1). d) The stress field in the Voronoi tesselation polycrystal. Brighter (hotter) shades represent greater intensity.

2 1/ 2

 22 as a function of the orientation α. The corresponding results are shown in Figure 4(c) for the case M = 1. Crosses correspond to 10 different configurations with hexagonal grains. The two irregular solid lines are sample results obtained with Voronoi tesselations (10 tesselations were considered but only two are shown for the clarity of the figure)

Figure 5 .

 5 Figure 5. The structure of the stress and strain fields in a polycrystal. Brighter (hotter) shades represent greater intensity.

Figure 5 (+ ε 2 2 1 / 2

 5212 Figure 5(d) shows the magnitude of the strain (defined as ε 2 1 + ε 2 2 1/2 ) for α = 30 • . It is very reminiscent of experimental observations in shape-memory polycrystals (see for example, (8)). The strain field is more homogeneous than the stress fields, but still not fully homogeneous. The strain magnitude ranges from 0.48 to 1.414 while the stress magnitude ranges from 0 to over 100. However, one has to be careful about the interpretation. Since the set K is not isotropic, a strain magnitude of 1 can indicate locking if only one component is non-zero while a strain magnitude as large as 1.414 may indicate no locking when both components are equal. Further, the connection to the Taylor bound is also tricky. While uniform field gives the Taylor bound, the converse -the attainment of the Taylor bound implies uniform field -is not true.

  (b). The highest bound is the Sachs upper bound. A variational upper bound can be obtained by using the inequality (3.7) with the Hashin-Shtrikman upper bound (it does not have an explicit expression) and is shown as VB(HS+) in Fgure 4(b). The sharpest upper bound is the translation bound of Nesi, Smyshlaev and Willis (NSW) ε NSW 0 /ε (2) 0 = 3/2 when M = 1. The estimate based on the linear comparison polycrystal and the self-consistent estimate is ε SC 0

Figure 6 .

 6 Figure 6. Random polycrystals. a) Effective modulus of a linear elastic polycrystals as a function of M = µ (2) /µ (1) and b) Effective recoverable strain as a function of M = ε (1) 0 /ε (2) 0

Article submitted to Royal Society

Acknowledgement. This research was initiated when PS spent a sabbatical in Caltech during 2000-2001. We gratefully acknowledge the partial financial support of the US National Science Foundation, US Air Force Office of Scientific Research, French Délégation Générale à l'Armement during this work.