Reconstructing displacements from the solution to the periodic Lippmann-Schwinger equation discretized on a uniform grid

displacements from the solution to the periodic Lippmann-Schwinger equation discretized on a uniform grid


INTRODUCTION
The motivation for the present paper comes from the numerical homogenization of periodic media.The reader is referred to standard textbooks for a detailed account of this theory, from the mechanical point of view [1,2,3] and the mathematical point of view [4,5,6].In what follows, only the essential steps of the homogenization process are highlighted.

S. BRISARD
Let Ω = (0, L 1 ) × • • • × (0, L d ) be the unit cell of a heterogeneous, Ω-periodic material, where L i denotes the length of the unit cell Ω in the i-th direction for (i = 1, . . ., d; d: dimension of the physical space).Within the framework of linear elasticity, homogenization of periodic media requires the solution to the so-called local (or corrector) problem.This problem amounts to finding the Ω-periodic displacement u such that in the whole space R d .In the above equations, C denotes the (Ω-periodic) local stiffness and E denotes the macroscopic (imposed) strain; finally, ∇ s u denotes the symmetric part of the gradient of the displacement u.It should be noted that the solution to Eq. ( 1) is unique up to an additive constant, which need not be specified for the time being.In Sec.4.4, we will explicitely impose the average displacement to be null.
From the solution to problem (1), the effective stiffness C eff is derived as the linear operator which relates the macroscopic stress to the macroscopic strain E Solving problem (1) can be a daunting task if the microstructure is highly heterogeneous (that is, if C is highly oscillatory within the unit cell Ω).In a seminal paper, Moulinec and Suquet [7] (see also [8]) introduced a new, efficient numerical scheme on a cartesian grid.It relies on the discretization of the so-called Lippmann-Schwinger equation [9,10,11] which we now proceed to introduce.It is first observed that problem (1) is equivalent to finding the Ω-periodic displacement u such that where C 0 denotes the elastic moduli of some arbitrary, homogeneous, reference material.In the above equations, ε denotes the total strain, while τ is the so-called stress polarization.By definition of the Green operator for strains Γ 0 (see Appendix A), it is readily found that [compare Eqs.(3a) and (62)] where ' * ' denotes the standard (periodic) convolution product.Eliminating the stress polarization τ between Eqs. (3b), (3c) and (4) leads to the following integral equation which is known as the Lippmann-Schwinger equation.It should be noted that in Eq. ( 5), the main unknown is no longer the periodic part u of the displacement, but the total strain ε defined by Eq. (3b).Solving Eq. ( 5) numerically proceeds in two steps.First, the continuous equation must be discretized.Then, a suitable strategy must be adopted to find the solution to the resulting linear system.
Regarding the discretization, it is convenient to introduce the concept of discrete Green operator [12], which provides a unified framework to a wide range of discretization strategies: truncation of high frequencies [8], exact evaluation for cell-wise constant stress polarizations [13], low-pass filtering [12], finite elements [14] or finite differences [15,16].Although a comparison of these approaches is out of the scope of the present paper, the recent work by Willot [15,16] deserves a special mention as the most promising strategy.
As for the solution to the resulting discretized Lippmann-Schwinger equation, various iterative schemes have been proposed, namely: the original fixed-point iterations of Moulinec and Suquet [7,8], the accelerated scheme of Eyre and Milton [17], the augmented Lagrangian solver of Michel et al. [18], Krylov subspace methods [13,12], or the polarization-based scheme of Monchiet and Bonnet [19] (see also [20] for a comparison of some of these schemes).All these schemes share a common matrix-free approach relying on the fast Fourier transform to efficiently compute matrix-vector products, as was first suggested by Moulinec and Suquet [7,8].
The combination of a discrete Green operator and an iterative linear solver constitutes what we will call a Uniform Grid Periodic Lippmann-Schwinger solver in the remainder of this paper (in short, UGPLS solver).The exact discretization and solving strategies are not relevant to the present work, as long as it is understood that UGPLS solvers produce an estimate of the strain (or, equivalently, the stress-polarization), not the displacement, on a uniform cartesian grid.While in most use-cases, only the local stresses and strains are required, it is sometimes desirable to produce kinematically admissible displacement fields, as illustrated below with two examples.
As a first example, we consider the limit analysis of periodic, heterogeneous materials, which requires kinematically admissible displacement fields in order to produce rigorous upper bounds on the macroscopic strength criterion.This can be achieved by means of dedicated finite element techniques [21,22].Alternatively, the problem can be reformulated as the homogenization of a non-linear, viscous material [2], for which UGPLS solvers can be used iteratively (see [7] for the first application of UGPLS solvers to non-linear materials).The numerical solution is then post-processed, delivering a rigorous upper bound on the collapse load.It should be noted that this approach requires the local strength criterion to be smooth or regularized [23].
As a second example, we mention the a posteriori evaluation of the error in constitutive relation as first introduced by Ladevèze and Leguillon [24] within the framework of finite element analysis.From the Prager-Synge theorem, it can be shown that the error in constitutive relation provides an upper bound on the approximation error in the energy norm (see [25] for a review, or the textbook by Ladevèze [26]).Evaluation of the error in constitutive law requires both a kinematically admissible displacement field (which is addressed in the present paper), and a statically admissible stress field (which will be presented in a paper to come).
To sum up the above discussion, UGPLS schemes are attractive techniques to solve the corrector problem in an efficient way, but fail to deliver a kinematically admissible displacement field.In the present paper, we seek to overcome this shortcoming and propose an a posteriori construction of a kinematically admissible displacement field.The proposed method applies to any UGPLS solver (regardless of the discretization strategy and the linear solver).The ideas underlying the construction are described now.
Let ε N denote the discrete strain field resulting from a UGPLS analysis, where the supscript 'N' is a discretization parameter which relates to the fineness of the grid (see Sec. 3.1).In a Galerkin setting, ε N can be considered as a cell-wise constant estimate of the true strain ε [12].From Eq. (3c), the discrete stress polarization τ N is evaluated.Likewise, it is seen as a cell-wise constant estimate of the true stress polarization τ.
If the true stress polarization τ were known exactly, then the displacement u would be the periodic solution to Eq. (3a).We therefore seek the approximate displacement field u N as the periodic solution to Eq. (3a) where the true stress polarization τ is replaced with the approximate, cell-wise constant, In other words, reconstructing the displacement field resulting from a UGPLS analysis reduces to solving the problem of the linear elastic equilibrium of a prestressed, homogeneous material with periodic boundary conditions.Of course, the exact solution to this problem is not known; rather, a finite element estimate of this solution is seeked.The remainder of this paper is therefore devoted to the numerical solution of the following problem (formulated on the whole space R d ) where the unknown displacement u is Ω-periodic.In the above equation, is the (Ω-periodic) prestress.Problem (7) is slightly more general than problem (6) due to the body forces with Ωperiodic volume density b.
Problem (7) is solved numerically by means of a finite element technique.The displacement is discretized over a d-dimensional cartesian grid; interpolation is achieved through P1 shape functions.Periodicity of the discretized displacement is enforced by construction (using periodic shape functions).Due to the periodic boundary conditions and the homogeneity of the material, the finite element approximation of Eq. ( 7) results in a block-circulant stiffness matrix.This suggests to formulate this approximation in the Fourier space, where all frequencies are uncoupled.More precisely, while the nodal displacements are the solution to a (N nodes d) × (N nodes d) system (where N nodes is the total number of nodes), the modal displacements (that is, the discrete Fourier transform of the nodal displacements) are the solution to N nodes systems of size d × d, which can be solved at virtually no cost.Of course, this approach requires the computation of direct and inverse discrete Fourier transforms: using the fast Fourier transform guarantees the overall efficiency of the proposed method.
At this point, the alternative approach of Vondřejc and collaborators (see [27,28] and references therein) should be mentioned.These authors propose a Galerkin discretization of the corrector problem with trigonometric polynomials: their approach therefore delivers non-local, kinematically admissible, approximations of the displacement field.Local approximations might sometimes be more convenient, and the approach proposed in the present paper should be seen as complementary to that of Vondřejc and collaborators.It should also be noted that, since high frequencies are discarded from the approximation space, the displacement and strain fields delivered by Vondřejc and collaborators might be subject to the Gibbs phenomenon.
The remainder of this paper is organized as follows.Sec. 2 introduces some convenient elementby-element operations on vectors, which allow for more compact expressions; basic results relating to Fourier series and the discrete Fourier transform are also recalled.We then define in Sec. 3 the discretization grid and the discretized displacement and strain.In particular, we relate (without approximation) the Fourier series expansion of the piecewise affine displacement field to the discrete Fourier transform of the nodal displacements.Sec. 4 constitutes the crux of the paper: we show that the global stiffness matrix and the vector of nodal forces have simple expressions in Fourier space.We therefore introduce the modal stiffness (d × d) matrices and force (d × 1) vectors.A summary of the method is provided in Sec.4.5.This paper closes on a 2D (plane strain elasticity) example (see Sec. 5).The reconstructed displacement field is used to produce rigorous upper bounds on the effective elastic moduli.

Element-by-element operations on the components of vectors
In the present paper, it will be convenient to consider vectors of R d (N d , Z d ) as tuples of size d, and define non-linear operations on the coordinates of the vectors.We first define element-by-element product and division of two vectors x and y as follows x j y j e j (8a) Then, f denoting a scalar function of a scalar variable, we define f (x) for any vector x ∈ R d as the product The above definitions lead to very compact expressions, which can be read naturally, as though they were written in a one-dimensional space.This allows the reader to concentrate on the overall structure of the equations, rather than their (dimension dependent) details.As a first example, where the sinc function (cardinal sine) is defined as follows As further examples, the following identities resulting from rules ( 8) and ( 9) will be extensively used below in particular, is the total number of elements in the cartesian grid described in Sec.3.1.
It should be noted that rules ( 8) and ( 9) have a shortcoming: they are not intrinsic operations.In other words, the results of these operations depend on the basis in which they are performed.However, it is observed that in the present setting, only the basis which is attached to the sides of the rectangular box Ω is relevant, and no change of basis will be performed or required.Therefore, we will continue to use bold letters to denote vectors and tensors, even if they should really be understood as their arrays of components in the global basis.

Fourier series and the discrete Fourier transform
In the present section, we briefly recall the essential formulas regarding Fourier series and the discrete Fourier transform.In particular, we clarify the normalization constants that we adopt.
Fourier series Under some regularity assumptions which are not recalled here, any Ω-periodic function f can be expressed as the sum of a Fourier series where the basis functions Ψ k are defined as follows for k ∈ Z d and the Fourier coeffcients fk are related to f Parseval's theorem then relates the L 2 Hermitian product of two Ω-periodic functions f and g to the 2 Hermitian product of their Fourier coefficients where overlined quantities denote complex conjugates.
The discrete Fourier transform We now consider a sequence of numbers (y n ) n∈Z d , which we assume N-periodic (N ∈ N d ), that is y n+pN = y n for all n, p ∈ Z d .Clearly, this infinite sequence is fully defined by the finite sequence (y n ) n∈N , where N denotes the following set of multi-indices S. BRISARD The discrete Fourier transform of (y n and we have the inversion formula The Plancherel theorem relates the 2 Hermitian product of two N-periodic sequences (y n ) n∈Z d and (z n ) n∈Z d to the 2 Hermitian product of their discrete Fourier transforms n∈N

P1 DISCRETIZATION ON A PERIODIC, UNIFORM GRID
In the present section, we define the discretized displacement and strain fields.Both are Ω-periodic; as such, they can be expanded into Fourier series.Since the displacement is piecewise affine, we show that the Fourier series expansions take very simple forms.It will be seen that the discrete Fourier transform of the nodal values of the displacement arises naturally in the derivation, without any approximation.The fast Fourier transform can therefore be used to switch efficiently to and from Fourier space.In Sec.3.1, we define the discretization grid, and show that cell-averages can be computed in Fourier space.Then, we define in Sec.3.2 the discretized displacement, and show that its Fourier series expansion is completely defined by the discrete Fourier transform of the nodal displacements.Finally, spatial derivation of the discretized displacement leads in Sec.3.3 to the discretized strain.

The discretization grid
In the present section, we define the cartesian grid which is to be used to discretize the mechanical fields resulting from problem (7).In the i-th direction (i = 1, . . ., d), the side of the unit cell Ω is divided into N i equal segments of size L i /N i (see Fig. 1).As a result, the vertices of the grid are located at where L = L i e i and N = N i e i .In the above formula, the element-by-element operations introduced in Sec.2.1 have been used.We further define the grid-cells Ω N n as follows where the set N previously defined by Eq. ( 18) should be understood as the set of node multi-indices.We will need in what follows expressions of the cell-averages of a Ω-periodic function f , as well as their discrete Fourier transform.We introduce the following notations for these quantities It should be noted that the the cell-averages of a Ω-periodic function are N-periodic; therefore, it is meaningful to compute the discrete Fourier transform of the cell-averages of a Ω-periodic function (see Sec. 2.2).Definitions (25a) and (25b) are applied below to the Fourier basis functions Ψ k defined in Eq. (15).We first find by straightforward integration where 1 denotes the vector 1 = e 1 + • • • + e d , and the sinc function has been defined above [see Eq. ( 11)].Taking the discrete Fourier transform of both sides of Eq. ( 26), we find, for h ∈ N The above formulas can be used to compute the cell-averages of any Ω-periodic function f with Fourier coefficients fk .In particular, the following expression is found for the discrete Fourier transform of the cell-averages of f with

The discretized displacement and its Fourier series expansion
Let u N denote the P1 finite element approximation of the solution to problem (7), discretized over the cartesian grid (Ω N n ) n∈N defined in Sec.3.1.It is expressed as follows S. BRISARD Figure 2. The periodic shape function Φ N in the 2D case, with where Φ N denotes the Ω-periodic P1 shape function (to be defined below), and u N n are the nodal values of the discretized displacement which are N-periodic owing to the Ω-periodicity of the discretized displacement.As already argued in Sec.2.2, the discretized displacement is therefore fully defined by its values at nodes n ∈ N.
The Ω-periodic shape function it is then extended to R d by Ω-periodicity (see Fig. 2).Thus defined, Φ N is the periodic counterpart of the usual P1 shape function.Its Fourier series expansion reads (see Appendix B) and, upon substitution in Eq. ( 29) where Eq. ( 22) has been used.In the above expression, the inner sum is recognized as the discrete Fourier transform of the nodal displacements (u N n ) n∈N , so that we finally find the compact expression of the Fourier series expansion of the discretized displacement u where the Fourier coefficients ũN k are given by the following formula 3.3.The discretized strain and its Fourier series expansion Let ε N = ∇ s u N be the strain associated to the discretized displacement u N .Its Fourier series expansion is readily obtained from the differentiation of Eq. ( 34) where the following expression is found for the Fourier coefficients εN Combining Eqs. ( 28), ( 37) and (38) then leads to the following compact expression of the discrete Fourier transform of the cell averages of the strain where the modal strain-displacement vector BN k has been introduced A lengthy but straightforward (see Appendix C.2) calculation shows that and the other components of the modal strain-displacement vector are deduced by circular permutations of the indices.

FINITE ELEMENT APPROXIMATION OF THE PROBLEM
The approximate nodal displacements for the solution to problem (7) are found by minimization of the total potential energy Π = U − V.The strain energy U = U e + U p is the sum of an elastic contribution U e and a term arising from the prestress U p , and we have where V denotes the potential of the external loads.The three above contributions are evaluated in Secs.4.1, 4.2 and 4.3, respectively, by substitution of the discretized displacement u N (defined in Sec.3.2) or strain ε N (derived in Sec.3.3) into Eqs.(42a), (42b) and (42c).The total potential energy is then computed and optimized in Sec.4.4, leading to a series of low-dimensional linear systems in the Fourier space.A summary of the method is provided in Sec.4.5.

Evaluation of the strain energy (elastic contribution)
Owing to the homogeneity of the material, the elastic part U e of the strain energy [see Eq. (42a)] can readily be computed in Fourier space by means of Parseval's identity ( 17) where the Fourier coefficients εN k of the strain are given by Eq. ( 38) and A 0 denotes the acoustic tensor, A 0 (q) = q • C 0 • q.Replacing in the above series k ∈ Z d with k + pN (with k ∈ N and p ∈ Z d ), and taking advantage of the N-periodicity of the discrete Fourier transform ûN k , we find We then introduce the so-called modal stiffness matrices KN k defined as follows and expression (44) of the strain energy reduces to As expected, it is observed that for k = 0, KN 0 = 0. Indeed, the strain energy associated with constant displacements (associated with the k = 0 mode) is zero.For isotropic materials with shear modulus µ 0 and Poisson ratio ν 0 , a lengthy but straightforward calculation (see Appendix C.3) leads to the following closed-form expression of the modal stiffness matrices and the diagonal and off-diagonal coefficients of ĤN k are fully defined (up to a circular permutation of the indices) by the following identities

Evaluation of the strain energy (contribution of the prestress)
We now turn to the second term U p of the strain energy [see Eq. (42b)].We assume that the prescribed prestress is cell-wise constant, and denote N n its value in cell Ω N n .Then, from the Plancherel theorem ( 21) and, upon substitution of Eq. ( 39) Copyright c 0000 John Wiley & Sons, Ltd.Int.J. Numer.Meth.Engng (0000) Prepared using nmeauth.clsDOI: 10.1002/nme

Potential energy of the external loads
The volume density of body forces b is assumed cell-wise constant in Eq. (42c): b N n denotes its value in cell Ω N n .Then, using again the Plancherel theorem ( 21) which, upon substitution of Eq. ( 36), leads to the following expression where the modal forces

Optimization of the total potential energy
Gathering Eqs. ( 46), ( 50) and ( 52), leads to the following expression of the total potential energy the optimization of which (with respect to the modal displacements ûN k ) results in a set of d × d linear systems of equations (one such system for each value of k ∈ N) For k = 0, the linear system (55) is ill-posed, since KN 0 = 0 [see Eq. ( 45)].This is consistent with the fact that the solution to Eq. ( 7) is unique up to an additive constant.It is customary to impose the average displacement to be zero.Then, Eq. (55) must be replaced for k = 0 with ûN 0 = 0. (56)

Summary of the method
The procedure proposed in this paper for the reconstruction of the displacement field results from the above developments (see Secs.  56)] to find the modal displacement ûN k .4. Compute the inverse Fourier transform of the modal displacements ( ûN k ) k∈N to obtain the nodal displacements (u N n ) n∈N .In the absence of body forces, it is observed that the overall cost of Algorithm 1 is dominated by that of d + s discrete Fourier transforms, where s = d (d + 1) /2.This should be compared to the cost of one iteration of the UGPLS solver, which in general requires 2s discrete Fourier transforms.Therefore, the cost of the proposed reconstruction is commensurable with that of one iteration of the UGPLS solver.Of course, the fast Fourier transform is used to compute discrete Fourier transforms; the cost of the present computation therefore grows as |N| log |N|, where |N| is the total number of cells.
S. BRISARD Figure 3.The model microstructure studied in Sec. 5. Table I.Summary of the simulations carried out on the microstructure depicted in Fig. 3.The first column is the grid size, the second column is the number of CG iterations.The third column is the direct estimate from Eq. ( 59) of the 1212 macroscopic modulus and the fourth column is the bound on this modulus, computed by means of the displacement reconstructed from Algorithm 1 and the minimum potential energy principle.

NUMERICAL RESULTS
This section provides an application of Algorithm 1.The method proposed in this paper applies equally well in two and three dimensions.However, for the sake of simplicity, we selected a plane strain elasticity problem.The geometry of the biphasic microstructure is defined in Fig. 3.This kind of model microstructures has already been used for illustration purposes by many authors [13,12,15,27,16].The material constants of the two phases are and the unit cell is subjected to a macroscopic shear strain The corrector problem (1) is solved by means of a UGPLS solver combining a conjugate gradient (CG) solver [13] and the so-called filtered Green operator [12].The elastic constants of the reference material are µ 0 = 0.9µ m and ν 0 = 0.3.The tolerance for the stopping criterion of the CG solver was set to 10 −10 , which is extremely low.The computation was carried out on square grids of increasing size (64 × 64, 128 × 128, . . ., 4096 × 4096), and was distributed over 8 cores for the largest simulations.It is observed (see Table I, second column) that the total number of iterations is fairly stable regardless of the grid size.Fig. 4 shows maps of the numerical estimates of the stress polarization, τ N for For each grid-size N, the piecewise constant estimate of the stress polarization τ N provides an estimate of the 1212 effective elastic modulus.Indeed, taking the volume average of Eq. (3c), we have Copyright c 0000 John Wiley & Sons, Ltd.Int.J. Numer.Meth.Engng (0000) Prepared using nmeauth.clsDOI: 10.1002/nme Owing to the choice of the reference material, the stress polarization is almost null in the "matrix"; therefore, only the (0, L/2) × (0, L/2) region ("inclusion") is shown.The same color scale was used for all three maps.The present images correspond to the 32 × 32 simulation.On the left image are also represented the lines along which the displacement profiles are plotted in Fig. 6.
and the resulting estimates are shown in Table I (third column).Then, Algorithm 1 is used to reconstruct the displacement u N .Fig. 5 shows the resulting maps for N 1 = N 2 = 32.For validation purposes, Fig. 6 compares the 32 × 32 UGPLS results with a standard, 256 × 256 displacement-based finite element approximation of the displacement, which serves as a reference.A very good agreement is observed, despite the fact that the grid used for the Lippmann-Schwinger based approach is rather coarse.It should be noted that the reference solution was computed with the deal.II finite element library [29,30].
To close this section, the minimum energy principle is used to produce guaranteed upper bounds on the 1212 elastic modulus with the reconstructed displacements u N It is therefore necessary to evaluate the strain energy of the heterogeneous microstructure for the displacement u N .Unlike the homogeneous case, this is carried out in the real space, using the P1 element stiffness matrix.More precisely, where K n denotes the element stiffness matrix of element n, and q n is the column vector of the full nodal displacements (including the affine part: The resulting bounds are tabulated in Table I (fourth column), and plotted in Fig. 7, alongside the estimates of the same quantity derived   I, third column) with the upper bounds (triangles) based on Eq. (60) (see also Table I, fourth column).As expected, both estimates and bounds seem to converge to the same value; however, the estimates are more accurate than the rigorous bounds.
from the average stress polarization [see Eq. ( 59)].It is observed that both bounds and estimate seem to converge to the same value; this was an expected result, which provides a quantitative validation of the procedure presented in this paper.It should be noted that, for the same grid size, the estimates are slightly more accurate than the rigorous bounds (the maximum relative distance between estimate and bound is about 1 %).

CONCLUSION
A procedure is presented to reconstruct displacements from the output of a Lippmann-Schwinger solver on a uniform, rectangular grid with periodic boundary conditions.The reconstruction is formulated as a standard problem of elastic equilibrium of a homogeneous material (also with periodic boundary conditions).This auxiliary problem is solved by means of displacement-based P1 finite elements.Owing to the uniformity of the grid and the homogeneity of the material, the global stiffness matrix exhibits a block-circulant structure, which lends itself to an efficient implementation in Fourier space.As a Copyright c 0000 John Wiley & Sons, Ltd.
Int. J. Numer.Meth.Engng (0000) Prepared using nmeauth.clsDOI: 10.1002/nme result, the overall cost of the reconstruction is similar to the cost of one iteration of the Lippmann-Schwinger solver.It is emphasized that the procedure presented here ought to be regarded as a mere post-processing of the output of a Uniform Grid Periodic Lippmann-Schwinger (UGPLS) solver.As such, it applies to any UGPLS solver, regardless of the actual discrete Green operator that is used.The method applies to both two-and three-dimensional problems.A simple, two-dimensional (plane strain) illustration is provided.The reconstructed displacements are compared to a direct calculation with displacement-based finite elements; both approaches were found to be in excellent agreement.The reconstructed displacement field is then combined with the minimum potential energy principle to produce rigorous upper bounds on the effective elastic moduli.
This work constitutes the first step towards a posteriori error estimates for UGPLS solvers, based on the error in constitutive relation.To achieve this goal, both kinematically admissible displacement field and statically admissible stress field are required.The present paper addresses the former; as for the latter, we are currently investigating how standard techniques commonly used for the aposteriori error analysis of finite element solutions can be adapted to the present case.hat function is first defined over this range as follows it is then extended to the whole real line R by 1-periodicity (see Fig. 8).Λ N thus defined is the periodic counterpart of the standard one-dimensional Lagrangian shape function of degree one.The Fourier coefficients of Λ N are found by means of Eq. ( 16) and the Fourier series expansion of Λ N therefore reads Eq. ( 31) clearly shows that the d-dimensional periodic shape function Φ N is the tensor product of d one-dimensional periodic hat functions and upon substitution of Eq. (66) where the tensorization rules ( 8), ( 9), ( 12) and ( 13) have been used.Thus Eq. ( 32) is retrieved.

Copyright c 0000
John Wiley & Sons, Ltd.Int.J. Numer.Meth.Engng (0000) Prepared using nmeauth.clsDOI: 10.1002/nme Evaluation of the potential energy of the external loads (see Sec. 4.3) further requires the DFT of the cell-averages of the displacement, which can be expressed as follows (see Appendix C.1) 4.1, 4.2 and 4.3).It proceeds in four steps.Algorithm 1 (Reconstruction of the displacement) 1. Compute the discrete Fourier transform ( ˆ N k ) k∈N of the cell values ( N n ) n∈N of the prestress.2. Compute the discrete Fourier transform ( bN k ) k∈N of the cell values (b N n ) n∈N of the body forces.3.For each discrete frequency k ∈ N (a) Use Eq. (41) to compute the modal strain-displacement vector BN k .(b) Use Eqs.(47), (48a) and (48b) to compute the modal stiffness matrix KN k .(c) Use Eq. (53) to compute the modal force FN k .(d) Solve the d × d linear system (55) [or use Eq. (

Figure 4 .
Figure 4. Maps of the normalized components of the stress polarization: τ 11 /µ m (left), τ 22 /µ m (middle) and τ 12 /µ m (right).Owing to the choice of the reference material, the stress polarization is almost null in the "matrix"; therefore, only the (0, L/2) × (0, L/2) region ("inclusion") is shown.The same color scale was used for all three maps.The present images correspond to the 32 × 32 simulation.

Figure 5 .
Figure 5. Maps of the normalized components of the displacement: u 1 /(E 12 L) (left) and u 2 /(E 12 L) (right).Unlike Fig.4, the whole unit cell (0, L) × (0, L) is represented here.The same color scale was used for the two maps (mutually symmetrical about the first bisector); the present images correspond to the 32 × 32 simulation.On the left image are also represented the lines along which the displacement profiles are plotted in Fig.6.

1 = L/4 x 1 = 3L/ 4 Figure 6 .
Figure 6.Profile of the first component of the normalized displacement, u 1 /(E 12 L), along the lines x 1 = L/4 (squares) and x 1 = 3L/4 (triangles).Solid lines correspond to a reference, displacement-based FEM solution computed on a 256 × 256 cartesian grid with P1 elements.Symbols correspond to the displacement reconstructed from the stress polarization τ N , with N 1 = N 2 = 32.Both approaches are in excellent agreement.

Figure 7 .
Figure 7.Comparison of the estimates (squares) of C eff 1212 based on Eq. (58) (see also TableI, third column) with the upper bounds (triangles) based on Eq. (60) (see also TableI, fourth column).As expected, both estimates and bounds seem to converge to the same value; however, the estimates are more accurate than the rigorous bounds.

Figure 8 .
Figure 8.The periodic hat function Λ N (ξ) defined in Appendix B. In the above figure, N = 5, and h = 1/N.