Fourier-based schemes for computing the mechanical response of composites with accurate local fields
FFourier-based schemes for computing the mechanicalresponse of composites with accurate local fields
Franc¸ois Willot
Mines ParisTech, PSL Research University, Centre for Mathematical Morphology,35, rue St Honor´e, 77300 Fontainebleau, France
Abstract
We modify the Green operator involved in Fourier-based computational schemesin elasticity, in 2D and 3D. The new operator is derived by expressing continuummechanics in terms of centered di ff erences on a rotated grid. Use of the modifiedGreen operator leads, in all systems investigated, to more accurate strain and stressfields than using the discretizations proposed by Moulinec and Suquet (1994) orWillot and Pellegrini (2008). Moreover, we compared the convergence rates ofthe “direct” and “accelerated” FFT schemes with the di ff erent discretizations. Thediscretization method proposed in this work allows for much faster FFT schemeswith respect to two criteria: stress equilibrium and e ff ective elastic moduli. Keywords:
FFT methods, Homogenization, Heterogeneous media, Linearelasticity, Computational mechanics, Spectral methods
1. Introduction
Fourier-based algorithms, or “FFT” methods for short, are an e ffi cient ap-proach for computing the mechanical response of composites. Initially restrictedto linear-elastic media, FFT tools are nowadays employed to treat more involvedproblems, ranging from viscoplasticity [1] to crack propagation [2]. In FFT meth-ods, the microstructure is defined by 2D or 3D images and the local stress andstrain tensors are computed along each pixel or “voxel” in the image. Coupledwith automatic or semi-automatic image segmentation techniques [3], this al-lows for the computation of the mechanical response of a material directly from Email address: [email protected] (Franc¸ois Willot)Tel.: +
33 1 64 69 48 07.
Preprint submitted to CRAS December 30, 2014 a r X i v : . [ m a t h . NA ] D ec xperimental acquisitions, like focused ion beam or 3D microtomography tech-niques [4]. The latter often deliver images containing billions of voxels, for whichFFT methods are e ffi cient [5, 6]. This allows one to take into account representa-tive volume elements of materials which are multiscale by nature such as concreteor mortar [7]. From a practical viewpoint, the simplicity of FFT methods is at-tractive to researchers and engineers who need not be experts in the underlyingnumerical methods to use them. Nowadays, FFT tools are available not only asacademic or free softwares [8, 9] but also as commercial ones [10].In the past years, progress has been made in the understanding of FFT algo-rithms. Vondˇrejc and co-workers have recently shown that the original methodof Moulinec and Suquet [11] corresponds, under one technical assumption, toa particular choice of approximation space and optimization method [12] (seealso [13]). This property allows one to derive other FFT schemes that use standardoptimization algorithms, such as the conjugate gradient method. In this regard,making use of variational formulations, e ffi cient numerical methods have beenproposed that combine FFTs with an underlying gradient descent algorithm [14,15].Di ff erent approximation spaces or discretization methods have also been pro-posed, where, contrary to the original scheme, the fields are not trigonometricpolynomials anymore. Brisard and Dormieux introduced “energy-based” FFTschemes that rely on Galerkin approximations of Hashin and Shtrikman’s varia-tional principle [15, 13] and derived modified Green operators consistent with thenew formulation. They obtained improved convergence properties and local fieldsdevoid of the spurious oscillations observed in the original FFT scheme [13, 16].In the context of conductivity, accurate local fields and improved convergencerates have also been obtained from modified Green operators based on finite-di ff erences [17]. These results follow earlier works where continuum mechan-ics are expressed by centered [18, 19] or “forward and backward” finite di ff er-ences [20].This work focuses on the e ff ect of discretization in FFT methods. It is or-ganized as follows. We first recall the equations of elasticity in the continuum(Sec. 2). We give the Lippmann-Schwinger equations and the “direct” and “accel-erated” FFT schemes in Sec. (3). In Sec. (4), a general formulation of the Greenoperator is derived that incorporates methods in [20], and a new discretizationscheme is proposed. The accuracy of the local stress and strain fields are exam-ined in Sec. (5) whereas the convergence rates of the various FFT methods areinvestigated in Sec. (6). We conclude in Sec. (7).2 . Microstructure and material elastic response We are concerned with solving the equations of linear elasticity in a square orcubic domain
Ω = [ − /
2; 1 / d in dimension d ( d = σ i j ( x ) = C i j , kl ( x ) ε ( x ) , ∂ i σ i j ( x ) ≡ , ε i j ( x ) = (1 / (cid:104) ∂ i u j ( x ) + ∂ j u i ( x ) (cid:105) , (1)where ε ( x ) is the strain field, σ ( x ) the stress field, u ( x ) the displacement vectorfield, C ( x ) the local elasticity tensor and x is a point in Ω . Tensorial componentsrefer to a system of Cartesian coordinates ( e ; e ) in 2D and ( e ; e ; e ) in 3D. Thematerial has an isotropic local elastic response that reads: C i j , kl ( x ) = λ ( x ) δ i j δ kl + µ ( x )( δ ik δ jl + δ il δ jk ) , (2)where δ is the Kronecker symbol and λ ( x ) and µ ( x ) are constant-per-phase Lam´e’sfirst and second coe ffi cients. The local bulk modulus κ = λ + (2 / d ) µ and the elasticmoduli take on values: λ ( x ) = λ α , κ ( x ) = κ α , µ ( x ) = µ α , in phase α . For simplicity, we restrict ourselves to binary media and, by conven-tion, α = α = ν = ν = .
25 so that, in 3D and 2D [21], we have µ α /κ α = .
6. The contrast of properties χ reads: χ = κ κ = µ µ = λ λ , (3)where 0 ≤ χ ≤ ∞ . In the matrix, we also fix κ = d = µ = . d = λ = . d = λ = . d = χ . In 3D, theYoung modulus is E = / E = χ/ χ = χ = ∞ .Periodic boundary conditions are applied with the material subjected to anoverall strain loading ε : σ i j ( x ) n j ( x ) − x ∈ ∂ Ω ) , (cid:104) ε kl ( x ) (cid:105) = ε kl , (4)where n is the normal at the boundary ∂ Ω of the domain Ω , oriented outward, − (cid:104)·(cid:105) denotes the spatial mean over Ω . The resultinge ff ective elastic tensor (cid:101) C is computed from: (cid:104) σ i j ( x ) (cid:105) = (cid:101) C i j , kl ε kl . (5)3 . Lippmann-Schwinger equation and FFT methods Fourier methods are by principle based on the Lippmann-Schwinger equa-tions. The latter follow from (1) and (4) as [22]: τ i j ( x ) = σ i j ( x ) − C i j , kl ε kl ( x ) , ε i j ( x ) = ε i j − (cid:90) x (cid:48) d d x (cid:48) G i j , kl ( x (cid:48) − x ) τ kl ( x (cid:48) ) , (6)where we have introduced a homogeneous “reference” elasticity tensor C and itsassociated polarization field τ and Green operator G . In the above we assume (cid:104) G (cid:105) = ε = (cid:104) ε ( x ) (cid:105) holds. The Green operator has, in the Fourier domain,the closed form [23]: G i j , kl ( q ) = (cid:26) q i (cid:104) q m C m j , kn q n (cid:105) − q l (cid:27) sym , (7)where q (cid:44) sym indicates minorsymmetrization with respect to the variables ( i , j ) and ( k , l ). Hereafter, we assumethat C is a symmetric, positive-definite, isotropic tensor defined by its bulk ( κ )and shear ( µ ) moduli, or Lam´e coe ffi cient ( λ ). Accordingly, when q (cid:44)
0, theGreen operator G is also symmetric definite and we have: G i j , kl ( q ) = µ (cid:32) q i q l | q | δ jk (cid:33) sym − λ + µ λ + µ q i q j q k q l | q | . (8)The “direct scheme” [11] consists in applying Eqs. (6) iteratively as: ε k = ≡ ε , ε k + = ε − G ∗ ( σ − C : ε k ) ( k ≥ . (9)In Moulinec and Suquet’s method, the convolution product ( ∗ ) above is computedas an algebraic product in the Fourier domain, making use of (8). Discrete Fouriertransforms are used to switch between the space Ω and Fourier domain F . Thisamounts to representing the strain field as a trigonometric polynomial [12] of theform: ε ( x ) = L d (cid:88) q ∈F ε ( q ) e i q · x , (10)where ε ( q ) is the discrete Fourier transform of ε ( x ). Accordingly: ε ( q ) = (cid:88) x ∈ Ω ε ( x ) e − i q · x . (11)4imilar forms are used for the stress and displacement fields. In practice, thedomain Ω is discretized on a square or cubic grid of L d voxels and the operator G ( q ) in (7) is evaluated along equispaced Fourier modes: q i = π m i L , m i = (cid:40) − ( L / , ..., L / , ( L even) , − ( L − / , ..., ( L − / , ( L odd) . (12)As noted in [24], when L is even, the relation: G ( q ) ∗ = G ( − q ) , (13)where G ( q ) ∗ is the complex conjugate of G ( q ), is not verified when one of thecomponents of q is equal to the highest frequency q i = π (i.e. m i = L / G ( q ) τ ( q ) used to compute thestrain field has non-zero imaginary part even if τ ( x ) is purely real. To fix thisproblem, we follow [24] and set: G ( q ) = (cid:16) C (cid:17) − , if q i = π for some i . (14)This choice enforces σ ( q ) = ε in (10) is not strictly-speaking irrotational, because of the lack ofconstraint at high Fourier modes for the strain field. We briefly mention anotheroption that we explored in this work. It consists in forcing the symmetry by re-placing G ( q ) with: G ( q ) ∗ + G ( − q )2 , (15)when one of the components of q equals π , which enforces ε ( q ) = G ( q ) at the highestfrequencies has little influence on the convergence rate, except at small resolution.When L <
128 pixels, faster convergence was achieved with the choice G ( q ) = (cid:16) C (cid:17) − . Furthermore, in the 2D example studied in this work, the choice G ( q ) = (cid:16) C (cid:17) − led to smaller oscillations, consistently with observations in [24]. The useof (15) was therefore not pursued further. We emphasize that, when L is odd, thisdiscrepancy disappears and no special treatment is needed.Refined FFT algorithms have been introduced to overcome the slow conver-gence rate of the direct scheme, observed for highly-contrasted composites, mostnotably the “accelerated scheme” [25] and “augmented Lagrangian” [26] meth-ods. In this work, we use the extension of the accelerated scheme to elastic-ity [26, 27]: ε k + = ε k + (cid:16) C + C (cid:17) − : C : (cid:104) ε − ε k − G ∗ (cid:16) C : ε k − C : ε k (cid:17)(cid:105) . (16)5he convergence rates of the accelerated and direct schemes depend on the choiceof the reference tensor C . For the accelerated scheme the optimal choice is [25,27]: κ = √ κ κ , µ = (cid:112) µ µ . (17)For the direct scheme, upper bounds on the eigenvalues of the Green operatorsuggest the choice [27]: κ = β (cid:16) κ + κ (cid:17) , µ = β (cid:16) µ + µ (cid:17) , (18)with β = /
4. Discretization and approximation space
In this section, we derive the expression of a modified Green operator G (cid:48) thatreplaces G defined in (8) and (12). We give it in a form that includes previously-proposed modified operators [20] and also introduce a new one. In the following, we assume that the strain and stress fields are defined on agrid of points in 2D, one per pixel. Eqs. (10) and (11) are used to apply discreteFourier transforms, but we do not postulate a representation in the continuum any-more. The equilibrium and strain admissibility conditions (1) are approximatedby means of finite di ff erences on which we apply the discrete transforms (10)and (11). In [20], this results in the following form: k ∗ i ( q ) σ i j ( q ) = , ε i j ( q ) = (1 / (cid:104) k i ( q ) u j ( q ) + k j ( q ) u i ( q ) (cid:105) , (19)where k and k ∗ represent “discrete” gradient and divergence operators, respec-tively. In the centered scheme, one takes k equal to: k C i ( q ) = i sin ( q i ) , (20)whereas in scheme [20], one chooses, for k : k W i ( q ) = e i q i − . (21)These expressions correspond, respectively, to the centered scheme: ∂ j σ i j ( x ) ≈ σ i j ( x + e j ) − σ i j ( x − e j )2 , ∂ j u i ( x ) ≈ u i ( x + e j ) − u i ( x − e j )2 , (22)6nd to the forward-and-backward di ff erence scheme: ∂ j σ i j ( x ) ≈ σ i j ( x ) − σ i j ( x − e j ) , ∂ j u i ( x ) ≈ u i ( x + e j ) − u i ( x ) . (23)Using (19), the resulting Green operator reads: G (cid:48) i j , kl ( q ) = (cid:26) k i ( q ) (cid:104) k m ( q ) C m j , kn k ∗ n ( q ) (cid:105) − k ∗ l ( q ) (cid:27) sym , (24)which is homogeneous in k , so that, with r i = k i / | k | : G (cid:48) i j , kl ( q ) = ( λ + µ ) (cid:16) r i r ∗ l δ jk (cid:17) sym − λ Re (cid:16) r i r ∗ j (cid:17) Re (cid:16) r k r ∗ l (cid:17) − µ r i r j ( r k r l ) ∗ µ (cid:20) λ + µ ) − λ (cid:12)(cid:12)(cid:12) r + r (cid:12)(cid:12)(cid:12) (cid:21) , (25)where Re( · ) denotes the real part of the enclosed quantity. The denominator onthe right-hand side is strictly positive due to the triangle inequality. Hereafter, theoperator G (cid:48) is denoted by G C when k = k C and by G W when k = k W . Note thatthe operator G in (8) is recovered from (25) by setting k i ( q ) = i q i . Now the Greenoperator G (cid:48) is complex and follows the minor and major symmetries: G (cid:48) i j , kl ( q ) = G (cid:48) ji , kl ( q ) = G (cid:48) i j , lk ( q ) , G (cid:48) i j , kl ( q ) = (cid:104) G (cid:48) kl , i j ( q ) (cid:105) ∗ . (26)Furthermore we have: G C,W ( q ) ∗ = G C,W ( − q ) , (27)for the schemes (22) and (23), including when q i = π and when L is even, because k C , Wi is real at the frequency q i = π . Therefore, the fix (14) in Sec. (3) is notnecessary. However, a problem of a di ff erent nature arises when using the centeredscheme (20) when L is even. Eq. (25) does not define the Green operator G C atthe three frequencies q = ( π ; 0), (0; π ) and ( π ; π ), for which k C ( q ) =
0. This isbecause the second equation in (22) has in general non-unique solutions for thedisplacement field u ( x ). Indeed, when L is even, the displacement is defined up toa linear combination of 2-voxels periodic fields. They are given by the following8 independent fields: v m ( x ) = δ mn , v m ( x ) = ( − i δ mn , v m ( x ) = ( − j δ mn , v m ( x ) = ( − i + j δ mn , x = (cid:16) iL ; jL (cid:17) , i , j = − L , ..., L , n = , . (28)The operator G C remains finite when q approaches one of the modes ( π ; 0), (0; π )or ( π ; π ), but can not be continuously extended at these modes. To fix this problem,we set, for the centered scheme: G C i j , kl ( q ) = , if L is even and ( q i = q i = π ) for all i , (29)7hich enforces ε ( q ) = ε ( x ) is accord-ingly admissible, and the stress field σ ( x ) is divergence-free, in the sense of (22).We explored the alternate choice G C ( q ) = (cid:16) C (cid:17) − in (29). Almost identical conver-gence rates and oscillations were observed for the two options, however the choice G C ( q ) = (cid:16) C (cid:17) − does not produce an irrotational strain field and is not consideredfurther. We emphasize that no special treatment is required for the operator G W athigh frequencies since k W (cid:44) q (cid:44) G with G C or G W in (9) and (16), we derive “direct” and “ac-celerated” schemes that solve (22) or (23). In the limit of very fine resolution, wehave G C,W ( q ) ≈ G ( q ) when q →
0, which guarantees that the strain and stressfields do not depend on the employed discretization. This property holds for anychoice of k such as k ∼ i q when q → e + e and e − e [20]. In a domain con-taining a single centered disc, the scheme produces fields that break the axialsymmetries of the problem. In fact, the discretization (23) is actually one of fourpossible choices, all of them breaking the symmetries. Attempts to force the sym-metry by averaging over the four Green operators or over the fields themselves, asproposed in [20], are not explored in this work. The fomer method indeed leadsto less accurate “di ff use” local fields. The latter necessitates to run four di ff erentcomputations, in 2D, instead of one, which is cumbersome.In the rest of this section, we derive a discrete scheme in 2D di ff erent from (22)and (23). In this scheme, the displacement field is evaluated at the 4 corners of thepixels and the strain and stress fields are evaluated at the centers of the pixels. Wefirst express these fields in the 45 o -rotated basis: f = e + e √ , f = e − e √ , (30)by: u i = R iI u I , ε i j = R iI ε IJ R (cid:48) J j , σ i j = R iI σ IJ R (cid:48) J j , R iJ = − δ i δ J √ , (31)where uppercase indices refer to components in the rotated grid. We discretize (1)8 igure 1: A pixel with edges parallel to the Cartesian axis ( e ; e ). Superimposed: 45 o -rotatedbasis ( f ; f ). The strain and stress fields are evaluated at the pixel center x (square). The displace-ment and the divergence of the stress field lie along the pixel corners (disks). in the rotated basis by the centered di ff erences (see Fig. 1): σ IJ ( x ) = C IJ , KL ( x ) ε KL ( x ) , (32a) σ I ( x ) − σ I (cid:16) x − √ f (cid:17) + σ I (cid:16) x + f − f √ (cid:17) − σ I (cid:16) x − f + f √ (cid:17) = , (32b) ε KL ( x ) = √ (cid:104) u K (cid:16) x + f L √ (cid:17) − u K (cid:16) x − f L √ (cid:17) + u L (cid:16) x + f K √ (cid:17) − u L (cid:16) x − f K √ (cid:17)(cid:105) (32c)where x lie at the centers of the pixels and x ± f I / √ e ; e ) and applying the backward discreteFourier transform (10) we arrive again at (19) with the following expression for k : k R i ( q ) = i2 tan (cid:18) q i (cid:19) (cid:16) + e i q (cid:17) (cid:16) + e i q (cid:17) . (33)We denote by G R the corresponding Green operator, derived by substituting k = k R in (25). The operator G R is real and also verifies: G R ( q ) = G R ( − q ) . (34)However, when L is even, k R = q = ( π ; π ) and G R is not defined by (25)at this frequency. Again, this is because (32c) gives the displacement field up tolinear combinations of the 4 independent fields v , , (see 28). Accordingly we set: G R ( q ) = , when L is even, d = q = q = π, (35)which enforces strain compatibility and stress equilibrium, in the sense of (32).9 .2. Three dimensions We follow the same methodology in 3D. The equilibrium and strain admissi-bility conditions (19) are unchanged, as well as the expression for the vectors k C,W in (20) and (21) resulting from (22) and (23). In 3D, we also extend (33) as: k R i ( q ) = i4 tan (cid:18) q i (cid:19) (cid:16) + e i q (cid:17) (cid:16) + e i q (cid:17) (cid:16) + e i q (cid:17) , (36)for the rotated scheme. The strain and stress fields are now evaluated at the cen-ters of the voxels and the displacement field at their corners. Derivatives of thedisplacement are estimated by di ff erences at opposite corners. For the strain com-ponents ε and ε : ε ( x ) ≈ (cid:88) m = ± n = ± (cid:20) u (cid:18) x + e + m e + n e (cid:19) − u (cid:18) x − e − m e − n e (cid:19)(cid:21) , (37a) ε ( x ) ≈ (cid:88) m = ± n = ± (cid:20) u (cid:18) x + e + m e + n e (cid:19) − u (cid:18) x − e − m e − n e (cid:19) + u (cid:18) x + e + m e + n e (cid:19) − u (cid:18) x − e − m e − n e (cid:19)(cid:21) , (37b)where x lie at the center of a voxel. The expression for the strain component ε (resp. ε ) is obtained after exchanging the indicia 1 and 2 (resp. 1 and 3) in (37a).The component ε (resp. ε ) is derived from (37b) by exchanging the indicia 3and 1 (resp. 3 and 2). Stress divergence is discretized in a similar manner. Its firstcomponent reads: ∂ i σ i ( x ) ≈ (cid:88) m = ± n = ± (cid:20) σ (cid:18) x + e + m e + n e (cid:19) − σ (cid:18) x − e − m e − n e (cid:19) + σ (cid:18) x + e + m e + n e (cid:19) − σ (cid:18) x − e − m e − n e (cid:19) + σ (cid:18) x + e + m e + n e (cid:19) − σ (cid:18) x − e − m e − n e (cid:19)(cid:21) , (38)where x lie at one of the edges of a voxel. The components ∂ i σ i and ∂ i σ i areobtained from (38) by circular permutations of the indicia 1, 2 and 3. We notethat (37) and (38) are the natural generalization of (32) to d =
3, expressed in theCartesian basis ( e ; e ; e ). 10n 3D, Eq. (24) yields, for the Green operator: G (cid:48) i j , kl ( q ) = (cid:16) λ + µ (cid:17) (cid:16) r i r ∗ l δ jk (cid:17) sym + λ (cid:20)(cid:16) r i r ∗ l s jk (cid:17) sym − Re (cid:16) r i r ∗ j (cid:17) Re (cid:16) r k r ∗ l (cid:17)(cid:21) − µ r i r j r ∗ k r ∗ l µ (cid:20) λ + µ ) − λ (cid:12)(cid:12)(cid:12) r + r + r (cid:12)(cid:12)(cid:12) (cid:21) (39)where again r i = k i / | k | and s is the symmetric second-order tensor: s j j = r i r ∗ k ) , s jk = − r k r ∗ j ) Im( r k r ∗ i ) , i (cid:44) j (cid:44) k (cid:44) i , (40)with Im( · ) the imaginary part of the enclosed complex quantity. Like in 2D, theoperator G (cid:48) follows minor and major symmetries (26).Again, the operators G C , G W and G R are derived using the expression for G (cid:48) in (39) with k = k C , k W and k R , respectively. The symmetries (27) and (34)are verified in 3D as well. But again, a special treatment is needed for G C and G R when L is even, at the modes q for which k C,R ( q ) =
0. Like in 2D, thedisplacement is undefined at these frequencies and the Fourier coe ffi cients of thestrain field are zero and so we set G C,R = d = G R ( q ) = , if L is even, d = q i = q j = π with i (cid:44) j . (41)The operators G C , G W and G R are, in 2D and 3D, periodic functions where,contrary to G , high frequencies are cut. Accordingly, we expect faster conver-gence rates for schemes using operators derived from finite di ff erences and moreexact local fields, as was previously observed in the conductivity problem [17].We also expect higher accuracy for the local fields when employing G R ratherthan the other discrete operators G W and G C . First, the operator G R is based oncentered di ff erences which are more precise than forward and backward di ff er-ences, used in G W . Second, derivatives are evaluated more locally when using G R rather than G C . Indeed, the latter are computed at points separated by 2 voxels for G C instead of √ √ G R .The above considerations guided the choice for the discretization schemes (37)and (38), leading to k R and G R . Clearly, many other choices are possible, andEq. (39) gives a general class of Green operators based on finite-di ff erences. Thelatter depend on the choice for the complex vector k . However, a systematicinvestigation of such discrete schemes is beyond the present study.In the the rest of this study, we estimate the accuracy of the local fields and ofthe e ff ective properties predicted by the various schemes, as well as their conver-gence rates. We denote by DS and AS the direct and accelerated schemes defined11y (9) and (16) respectively, when G is used. We denote by DS C , DS W , DS R and AS C , AS W and AS R , the same algorithms obtained by substituting G with G C,W,R respectively. We emphasize that, for a given Green operator, the direct and accel-erated schemes produce the same strain and stress fields, up to round-o ff errors.
5. Local strain and stress fields accuracy
Hereafter we consider the 2D ‘four-cell’ microstructure, where the periodicdomain Ω is divided into 4 identical squares of surface fraction 25%. Its non-trivial solution with singular fields at the corners makes it a good benchmark fornumerical schemes. Furthermore, the microstructure is discretized exactly at anyresolution, provided L is even. In the following, we make use of a simplified ver-sion of the four-cell microstructure made of a single quasi-rigid square inclusionembedded in a matrix (Fig. 2). We set the contrast to χ = . The material issubjected to the macroscopic strain loading: ε i j = (cid:16) δ i δ j + δ j δ i (cid:17) . (42)We determine the strain and stress fields predicted by FFT schemes when usingthe Green operators G or G C,W,R . The fields are computed using the acceleratedscheme (16) at discretizations L = − . These variations are the e ff ect of round-o ff errors in dou-ble precision floating point numbers. These computations allow us to compare thee ff ect of the discretization, independently of the algorithm used for convergence.We focus on the stress component σ ( x ) parallel to the applied loading in asmall region [ − .
04; 0 . around the corner of the inclusion (Fig. 3). At lowresolution L = . .
5, which amount to 0 .
24% ofthe pixels. Using the same color scale for all images, the smallest stress value,equal to 1 .
5, is shown in dark blue whereas the highest, equal to 3 .
5, is in darkred. Green, yellow and orange lie in-between.As expected, in the limit of very fine resolution, all methods tend to the samelocal stress field, as shown by the similar field maps obtained at resolution L = G leads to spurious oscillations alongthe interfaces of the inclusion, up to resolutions as big as 2048 pixels, a side-e ff ect12 , κ μ , κ e e Figure 2: Elementary periodic domain
Ω = [1 /
2; 1 / containing a square inclusion with elasticmoduli µ , κ (top-left, shown in white) embedded in a matrix (shown in gray) with elastic moduli µ , κ . noticed in [17] in conductivity. The oscillations do not disappear after computinglocal averages of the fields (not shown).Strong oscillations are produced by schemes using G C as well, not only in thequasi-rigid inclusion, but also in the matrix. We observe checkerboard patternsin the former, and vertical and horizontal alignments in the latter, at resolution1024 . These oscillations are greatly reduced by the use of G W . Still, due to thenon-symmetric nature of G W , the stress is not correctly estimated along a lineof width 1 pixel oriented upward from the inclusion corner. Similar patterns areobserved, in other directions, along the three other corners of the inclusion (notshown). These issues are solved when using G R which produces a stress field thatrespects the symmetries of the problem. Furthermore, use of G R greatly reducesoscillations compared to G and G C . In this section, we consider a 3D material analogous to the four-cell mi-crostructure in 2D. We divide the periodic domain into 8 identical cubes of volumefraction 12 . Ω and contained in the region [ − /
4; 1 / . Again, we apply a macroscopic strainloading of the form (42). The inclusion is quasi-rigid compared to the matrix withcontrast of properties χ = . We compute the strain and stress fields predictedby each Green operator using the accelerated scheme. As in Sec. (5.1) we let theiterations converge up to round-o ff errors in double precision.A 2D section of the stress component σ ( x ) is represented in Fig. (4). The13 = L = L = GG C G W G R Figure 3: Stress component σ ( x ) predicted by the various FFT schemes at the three resolutions L = − .
04; 0 . . The center of the region is thebottom-left corner of the square inclusion in Fig. (2). e , of equa-tion x = − . .
5, this time less than 0 .
04% of the voxels, andrepresent all field using the same color scale.At high resolution L = G and G C,W,R are close to one another. However, stress patterns near the corners of the inclusionare less pronounced with G than with the other methods. At smaller resolutions L =
256 and L = G are notably di ff erent fromthe others, suggesting slower size-convergence with this operator. Furthermore,the field maps computed at resolution L =
512 confirms the results obtained in2D: strong oscillations are observed inside the inclusion when using G and G C .The two methods produce artificial patterns directed vertically and horizontally,close to the interface. Conversely, the fields produced by G W and G R have thesmallest oscillations, but that of G W are not symmetric. When L = ff ect only slowly disappears when L is increased. Thesolution resulting from the use of G R does not su ff er from this problem. As in 2D,it produces symmetric fields. Furthermore, the latter are close to one another atall resolutions and contain almost no oscillations. Contrary to the previous sections, we now consider a microstructure with-out singularities (edges or corners) and focus on the e ff ect of the Green operatordiscretization on the e ff ective elastic properties. In the rest of this section, theelementary domain Ω contains one spherical inclusion of volume fraction 20%,so that the material is a periodic array of spheres. The spheres are very soft withcontrast of properties χ = − . We compute the e ff ective elastic modulus (cid:101) C , produced by either G or G C,W,R at increasing resolutions L =
32, 64, 128, 256 and512. Again, we use the accelerated scheme and iterations are stopped when thestress field maximum variation over two iterations in any pixel is less than 2 10 − .Results are shown in Fig. (5) and are compared with the analytical estimate in [28].When the resolution increases, the e ff ective elastic modulus (cid:101) C , increases up toa limit value that we estimate to about 1 . ± . G which has the slowest conver-gence with respect to the system size. At fixed resolution, the error on the pre-15 = L = L = GG C G W G R Figure 4: 2D section of the stress component σ ( x ) along the plane x = − . L = L =
512 and L = L C Sphere 20%, =10 - 4
Cohen (2004) ~ L= ¥ (estimate) Figure 5: Apparent elastic modulus (cid:101) C , estimated by FFT methods using the Green operators G and G C,W,R (black and red), at increasing resolution L . Orange: estimate in [28]; violet: estimateof the asymptotic e ff ective modulus (cid:101) C , using FFT data. dictions given by G is about 2 times larger than the one provided by G R , which,among all methods, gives the best estimate. The operators G C and G W stand in-between. This is another indication of the benefits of the operator G R .
6. Convergence rate
In this section, we estimate the rates of convergence of the direct and accel-erated schemes DS, DS
C,W,R , AS and AS
C,W,R , that use the various Green oper-ators. All schemes enforce stress equilibrium at convergence only, therefore wefollow [27] and consider a criterion based on the L -norm: η = || div( σ ) || ||(cid:104) σ (cid:105)|| = ||(cid:104) σ (cid:105)|| (cid:115) | Ω | (cid:90) Ω d d x | div( σ ) | = ||(cid:104) σ (cid:105)|| (cid:115)(cid:88) q | k ( q ) · σ ( q ) | , (43)where η (cid:28) ||(cid:104) σ (cid:105)|| is the Frobeniusnorm: ||(cid:104) σ (cid:105)|| = (cid:88) i , j (cid:104) σ i j ( x ) (cid:105) .
17n (43) we set k = k C,W,R for the schemes using G C,W,R and k = i q when usingthe Green operator G , so that k · σ ( q ) is the divergence of the stress field in theFourier domain, estimated according to the various discretization schemes.We now estimate the convergence rates on a random microstructure. In thefollowing, the domain Ω is a (periodized) Boolean model of spheres of resolu-tion L =
64 and volume fraction 17%, below the percolation threshold of thespheres — of about 29% [29]. To obtain meaningful comparisons, we use thesame randomly-generated microstructure for all schemes. This particular config-uration contains 743 spheres of diameter 5 voxels, about 13 times smaller than thesize of Ω .Taking ν = ν = ν = .
25 for the reference Poisson ratio, we computenumerically the number of iterations N ( E ) required to reach the precision η ≤ − , for varying reference Young moduli E , in the range 0 < E <
1. Weconsider the Boolean model of spheres with contrast χ = − and the variousaccelerated schemes AS and AS C,W,R (Fig. 6). Within the range 0 < E (cid:46) .
03, thenumber of iterations N ( E ) is about the same for all accelerated schemes. When E > .
03, however, N ( E ) is a strongly increasing function of E for scheme AS,contrary to the other schemes AS C,W,R . For the latter, N ( E ) decreases with E upto a local minimum, beyond which variations are much less sensitive to E . Oneunique local minimum around E ≈ .
09 is found for scheme AS R , whereas theschemes AS W and AS C exhibit two local minima.The e ff ect of the Poisson ratio is also investigated numerically. We let ν = . ± .
01 and 0 . ± .
05 for various values of E with χ = − and ob-serve a strong increase of the number of iterations N ( E ), for the schemes AS andAS C,W,R . The same behavior is observed for the direct scheme DS and DS
C,W,R with χ = − . Therefore, in the following, we fix the Poisson ratio to ν = . χ . Thisleaves one parameter, E , to optimize on. We use the gradient descent methodto determine a local minimum of N ( E ) for arbitrary contrast and scheme DS,DS C,W,R , AS and AS
C,W,R . As above, N ( E ) is the number of iterations necessaryto reach η ≤ − . We choose E = . E + E ) for schemes DS, DS C,W,R and E = √ E E for schemes AS and AS C,W,R as initial guess for E , suggestedby (18) and (17). At each step, we determine if E is to be increased or decreased,by comparing N ( E ) with N ( E + δ E ) where δ E = . E . It frequently happensthat N ( E ) = N ( E + δ E ). In that case, we compare the values of the precision η after N ( E ) iterations and follow the direction that minimizes η . Iterations arestopped whenever N ( E ) is unchanged after two descent steps.18 .0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 E N( E ) ASAS C AS W AS R χ =10 -5 Figure 6: Number of iterations N ( E ) required to achieve convergence, as a function of thereference Young modulus E , for the accelerated schemes AS and AS C,W,R using various Greenoperators. Convergence is achieved when the precision η = − is reached. The microsctructureis a Boolean model of quasi-porous spheres with χ = − . The gradient descent method determines a local minimum rather than theglobal minimum, which is sub-optimal. To check the validity of the results, fur-ther numerical investigations are carried out for χ = − , 10 and schemes DSand DS C,W,R . The method predicts the global minimum in these cases. This alsoholds for schemes AS and AS
W,R with χ = − , 10 , but not for scheme AS C with χ = − . However, in this case the number of iterations N ( E ) are very similar atthe two local minima, as shown in Fig. (6). In the following, the results given bythe gradient descent method are used as-is.Results for the optimal reference E are indicated in Fig. (7). For the directscheme, the optimal reference follows (18) with 0 . ≤ β ≤ . β smaller than 1 / ff ers depending on the Green operator in use. For the scheme AS withGreen operator G , the choice (17) is optimal except in the region χ ≤ − wherethe value of E tend to a small constant of about 0 .
01. Similar behavior is foundfor the schemes AS
C,R for which: E ≈ E + √ χ, (44)19ith E = .
07 for AS C and E = .
12 for AS R . Similar behavior has beenobserved numerically in [17], in the context of conductivity. For the scheme AS W that uses G W , the optimal choice for E follows the same pattern as above with E = . − ≤ χ ≤ − . This behavior is an e ff ect of thepresence of two local minima, similar to that shown in Fig. (6) for χ = − .Convergence rates, computed with optimized reference, are represented inFig. (8) as a function of the contrast, in log-log scale. Results for χ = χ when χ (cid:29) /χ when χ (cid:28)
1. For the accelerated scheme AS, the numberof iterations is smaller and follows √ χ when χ (cid:29) / √ χ when χ (cid:28)
1, withone exception. At very high contrast of properties χ < − , including at χ = , η = − .When χ <
1, the schemes DS
C,W,R and AS
C,W,R , that use G C,W,R , converge aftera number of iterations not exceeding 430. As shown in Fig. (8), the number ofiterations is nearly constant in the range 0 ≤ χ ≤ − . As expected, the accel-erated schemes AS C,W,R are faster than the direct schemes DS
C,W,R , with schemeAS R proving the fastest. With this scheme and when χ <
1, the number of itera-tions is at most 168. Again, these results are qualitatively similar with that givenin [17] in the context of conductivity.For rigidly-reinforced media ( χ >
C,W,R and AS
C,W,R follow the same powerlaw behaviors, with respect to χ , asthat of DS or AS. In all considered schemes, the number of iterations continu-ously increases with the contrast. Di ff erences are observed between the variousaccelerated schemes AS and AS C,W,R , with AS R the fastest. The use of Green op-erators associated to the problem for the strain field, as undertaken here, results inconvergence properties that are worse in the region χ > χ <
1. Inthis respect, benefits are to be expected from the use of dual Green operators [17],associated to the problem for the stress fields. ff ective elastic moduli In this section, we focus on the accelerated schemes AS and AS
C,W,R and ex-amine the rate of convergence of the various schemes with respect to the e ff ectiveelastic moduli. We consider the same Boolean microstructure as given in Sec. (6)but discretized on higher resolution grids of 256 and 512 voxels. The volumefraction of the spheres are respectively 16 .
82% and 16 . -6 -4 -2 χ -2 -1 E DSDS C DS W DS R ASAS C AS W AS R Figure 7: Optimal reference Young modulus E as a function of the contrast of properties χ , forthe various FFT methods, in log-log scale. Direct schemes: black and red (nearly superimposedto one another); accelerated schemes: blue and green. Results for a porous material ( χ =
0) areindicated at the left of the graph. The material is a Boolean model of spheres with volume fraction17%. -6 -4 -2 χ iterations DSDS C DS W DS R ASAS C AS W AS R Figure 8: Number of iterations as a function of the contrast of properties χ , for the various FFTmethods, in log-log scale. Direct schemes: black and red; accelerated schemes: blue and green.Results for a porous material ( χ =
0) are indicated at the left of the graph. The material is aBoolean model of spheres with volume fraction 17%. contrast of properties take on two values χ = − and χ = , so that the spheresare quasi-porous or quasi-rigid.We perform iterations of the schemes AS and AS C,W,R using the optimizedreference moduli found in the previous section, on the lower resolution grid. Weapply the macroscopic strain loading: ε i j = δ i δ j . At each iteration and for each scheme, we compute the elastic modulus (cid:101) C , ,derived from the mean (cid:104) σ (cid:105) of the stress component σ . The convergence ratetoward the elastic modulus is represented in Fig. (9) for quasi-porous spheres with L =
512 and in Fig. (10) for quasi-rigid spheres with L = R which are allrepresented. Dotted lines are guide to the eyes. In the porous case, much betterconvergence is obtained with scheme AS R than with schemes AS and AS C,W , asshown in Fig. (9). The estimate for (cid:101) C , predicted by AS and AS C,W presentstrong oscillations that are much reduced with AS R . After about 7 iterations, theestimate given by AS R is valid to a relative precision of 10 − . To achieve the sameprecision, more than 50 iterations are needed for schemes AS and AS C,W .22
10 20 30 40 50 60 70 80 90 iterations C ASAS C AS W AS R L=512, χ=10 −4 HS ~ Figure 9: Estimate of the elastic modulus (cid:101) C , as a function of the number iterations performed,for a 3D Boolean model of quasi-porous spheres. Black symbols: accelerated schemes AS, AS C ,AS W ; red: scheme AS R (orange: Hashin and Shtrikman upper bound). The situation is notably di ff erent for quasi-rigid spheres (Fig. 10). For allschemes, a much higher number of iterations is required to determine the elasticmodulus (cid:101) C , with a precision of 10 − . The slower convergence rate followsthat observed in Sec. (6), where convergence is much poorer for χ > χ <
1, and where AS R is less advantageous compared to the other schemes .Nevertheless, in this case also, as shown in Fig. (10), smaller oscillations areobserved in the estimate for (cid:101) C , when using AS R rather than schemes AS orAS C,W .
7. Conclusion
In this work, a novel discretization method has been proposed in 2D and 3D foruse in Fourier-based schemes. The core of the proposed scheme is a simple modi-fication of the Green operator in the Fourier domain. The results obtained confirmthose achieved in the context of conductivity [17]. Compared to schemes usingtrigonometric polynomials as approximation space, or to other finite-di ff erencesmethods, superior convergence rates have been observed in terms of local stressequilibrium, but also in terms of e ff ective properties. More importantly, the solu-tion for the local fields, predicted by the new discretization scheme is found to be23 iterations C ASAS C AS W AS R L=256, χ=10 ~ Figure 10: Estimate of the elastic modulus (cid:101) C , as a function of the number iterations performed,for a 3D Boolean model of quasi-rigid spheres. Black lines: accelerated schemes AS, AS C , AS W ;red: scheme AS R . more accurate than that of other methods, especially at the vicinity of interfaces.This property is important when applying FFT methods to solve more complexproblems like large strain deformation [31]. The new method also provides betterestimates for the e ff ective elastic moduli. Furthermore, its estimates does not de-pend on the reference medium, because the scheme is based on a finite-di ff erencesdiscretization of continuum mechanics. Although not explored in this work, themodified Green operator can be used with most other FFT iterative solvers, likethe “augmented Lagrangian” [32] or with FFT algorithms that are less sensitive tothe reference [14, 15], leading to the same local fields. Acknowledgements
The author thanks Carnot M.I.N.E.S for support through grant 20531.
References [1] S.-B. Lee, R. A. Lebensohn, and A. D. Rollett. Modeling the viscoplasticmicromechanical response of two-phase materials using fast fourier trans-forms.
International Journal of Plasticity , 27(5):707–727, 2011.242] J. Li, S. Meng, X. Tian, F. Song, and C. Jiang. A non-local fracture modelfor composite laminates and numerical simulations by using the ff t method. Composites Part B: Engineering , 43(3):961–971, 2011.[3] M. Faessel and D. Jeulin. Segmentation of 3D microtomographic imagesof granular materials with the stochastic watershed.
Journal of Microscopy ,239(1):17–31, 2010.[4] F. Willot, L. Gillibert, and D. Jeulin. Microstructure-induced hotspots in thethermal and elastic responses of granular media.
International Journal ofSolids and Structures , 50(10):1699–1709, 2013.[5] C. F. Dunant, B. Bary, A. B. Giorla, C. Pniguel, J. Sanahuja, C. Toulemonde,A. B. Tran, F. Willot, and J. Yvonnet. A critical comparison of several nu-merical methods for computing e ff ective properties of highly heterogeneousmaterials. Advances in Engineering Software , 58:1–12, 2013.[6] F. Willot and D. Jeulin. Elastic and electrical behavior of some random mul-tiscale highly-contrasted composites.
International Journal for MultiscaleComputational Enginneering: special issue on Multiscale modeling and un-certainty quantification of heterogeneous materials , 9(3):305–326, 2011.[7] J. Escoda, F. Willot, D. Jeulin, J. Sanahuja, and C. Toulemonde. Estimationof local stresses and elastic properties of a mortar sample by FFT computa-tion of fields on a 3D image.
Cement and Concrete Research , 41(5):542–556, 2011.[8] Morphhom software, http://cmm.ensmp.fr/morphhom , accessed Octo-ber 4, 2014.[9] CraFT software, http://craft.lma.cnrs-mrs.fr , accessed October 4,2014.[10] GeoDict software, , accessed October 4, 2014.[11] H. Moulinec and P. Suquet. A fast numerical method for computing thelinear and non linear mechanical properties of the composites.
Comptesrendus de l’Academie des sciences, S´erie II , 318(11):1417–1423, 1994.[12] J. Vondˇrejc, J. Zeman, and I. Marek. An FFT-based Galerkin method forhomogenization of periodic media.
Computers & Mathematics with Appli-cations , 68(3):156–173, 2014. 2513] S. Brisard and L. Dormieux. Combining Galerkin approximation techniqueswith the principle of Hashin and Shtrikman to derive a new FFT-basednumerical method for the homogenization of composites.
ComputationalMethods for Applied Mechanical Engineering , 217(220):197–212, 2012.[14] J. Zeman, J. Vodrejc, J. Novak, and I. Marek. Accelerating a FFT-basedsolver for numerical homogenization of a periodic media by conjugate gra-dients.
Journal of Computational Physics , 229(21):8065–8071, 2010.[15] S. Brisard and L. Dormieux. FFT-based methods for the mechanics of com-posites: A general variational framework.
Computational Materials Science ,49(3):663–671, 2010.[16] S. Brisard and F. Legoll. Periodic homogenization using the lippmann–schwinger formalism, 2014. http://arxiv.org/abs/1411.0330 .[17] F. Willot, B. Abdallah, and Y.-P. Pellegrini. Fourier-based schemes withmodified green operator for computing the electrical response of heteroge-neous media with accurate local fields.
International Journal for NumericalMethods in Engineering , 98(7):518–533, 2014.[18] W. H. M ¨uller. Mathematical vs. experimental stress analysis of inhomo-geneities in solids.
Journal de Physique , 6(C1):139–148, 1996.[19] C.M. Brown, W. Dreyer, and W.H. M¨uller. Discrete fourier transforms andtheir application to stress–strain problems in composite mechanics: a con-vergence study.
Proceedings of the Royal Society of London. Series A: Math-ematical, Physical and Engineering Sciences , 458:1967–1987, 2002.[20] F. Willot and Y.-P. Pellegrini. Fast Fourier transform computations andbuild-up of plastic deformation in 2D, elastic-perfectly plastic, pixelwise-disordered porous media. In
D. Jeulin, S. Forest (eds), “Continuum Modelsand Discrete Systems CMDS 11” , pages 443–449, Paris, 2008. ´Ecole desMines.[21] M. F. Thorpe and P. N. Sen. Elastic moduli of two-dimensional compos-ite continua with elliptic inclusions.
Journal of the Acoustical Society ofAmerica , 77(5):1674–1680, 1985.[22] G. W. Milton.
The Theory of Composites . Cambridge Univ. Press, Cam-bridge, 2002. 2623] S. Kanaun and V. Levin.
Self-consistent methods for composites . Springer,Dordrecht, The Netherlands, 2008.[24] H. Moulinec and P. Suquet. A numerical method for computing the overallresponse of nonlinear composites with complex microstructure.
ComputerMethods in Applied Mechanics and Engineering , 157(1):69–94, 1998.[25] D. J. Eyre and G. W Milton. A fast numerical scheme for computing the re-sponse of composites using grid refinement.
The European Physical JournalApplied Physics , 6(1):41–47, 1999.[26] J.-C. Michel, H. Moulinec, and P. Suquet. A computational scheme for lin-ear and non-linear composites with arbitrary phase contrast.
InternationalJournal for Numerical Methods in Engineering , 52(1-2):139–160, 2001.[27] H. Moulinec and F. Silva. Comparison of three accelerated FFT-basedschemes for computing the mechanical response of composite materials.
International Journal for Numerical Methods in Engineering , 97(13):960–985, 2014.[28] I. Cohen. Simple algebraic approximations for the e ff ective elastic moduliof cubic arrays of spheres. Journal of the Mechanics and physics of Solids ,52(9):2167–2183, 2004.[29] F. Willot and D. Jeulin. Elastic behavior of composites containing booleanrandom sets of inhomogeneities.
International Journal of Engineering Sci-ence , 47(2):313–324, 2009.[30] H. Moulinec and P. Suquet. Comparison of FFT-based methods for comput-ing the response of composites with highly contrasted mechanical properties.
Physica B: Condensed Matter , 338(1–4):58–60, 2003.[31] N. Lahellec, J.-C. Michel, H. Moulinec, and P. Suquet. Analysis of inhomo-geneous materials at large strains using fast fourier transforms. In
Proc. IU-TAM Symposium on Computational Mechanics of Solid Materials at LargeStrains , pages 247–258. Kluwer Academic Publishers, 2001.[32] J.-C. Michel, H. Moulinec, and P. Suquet. A computational method basedon augmented lagrangians and fast fourier transforms for composites withhigh contrast.