OPTIMIZATION OF BODIES WITH LOCALLY PERIODIC MICROSTRUCTURE BY VARYING THE PERIODICITY PATTERN

. This paper describes a numerical method to optimize elastic bodies featuring a locally periodic microscopic pattern. A new idea, of optimizing the periodicity cell itself, is considered. In previously published works, the authors have found that optimizing the shape and topology of the model hole gives a limited ﬂexibility to the microstructure for adapting to the macroscopic loads. In the present study the periodicity cell varies during the optimization process, thus allowing the microstructure to adapt freely to the given loads. Our approach makes the link between the microscopic level and the macroscopic one. Two-dimensional linearly elastic bodies are considered, however the same tech-niques can be applied to three-dimensional bodies. Homogenization theory is used to describe the macroscopic (eﬀective) elastic properties of the body. Numerical examples are presented, in which a cantilever is optimized for diﬀerent load cases, one of them being multi-load. The problem is numerically heavy, since the optimization of the macroscopic problem is performed by optimizing in simultaneous hundreds or even thousands of periodic structures, each one using its own ﬁnite element mesh on the periodicity cell. Parallel computation is used in order to alleviate the computational burden.

1. Introduction. The main motivation of the present paper comes from the study of bodies having locally periodic microstructure and optimization of their macroscopic properties, in the context of linearized elasticity. A body having locally periodic microstructure is a body whose material coefficients vary at a microscopic scale, and this variation is locally periodic in the sense that, around each point of the body, the material coefficients vary according to a periodic pattern. This pattern changes from point to point (at the macroscopic scale). Homogenization theory allows one to accurately describe the macroscopic behaviour of such a body by means of so-called cellular problems, which are elliptic PDEs subject to periodicity conditions, see e.g. [1], [10], [12]. Porous materials, that is, bodies with locally periodic infinitesimal perforations, can be described in a similar manner, see [11]. Other approaches to the optimization of locally periodic microstructures are described in [14] and [15]. In [7], the authors have presented an algorithm for optimization of locally periodic structures. This algorithm links the microscopic level with the macroscopic problem. At the microscopic level, the optimization process described in [7] focuses on the shape of the microscopic hole and on its topology (that is, the number of its connected components). However, the periodicity group (or, more precisely, the two vectors defining the periodicity cell) was kept fixed along the optimization process.
In the present paper, the authors present new results obtained by varying the vectors defining the periodicity. The structure thus obtained has much more freedom do adapt itself to the macroscopic loads, by adjusting the periodicity cell's aspect and orientation along the optimization process. In order to compute the desired variation of the vectors defining the periodicity, the authors have computed the derivatives of the homogenized elastic coefficients with respect to these two vectors. These derivatives are then linked with the derivative of the macroscopic objective functional with respect to the elastic coefficients.
The layout of the paper is as follows. Section 2 presents the mathematical background on periodic microstructures. Section 3 focuses on the computation of the derivatives of the homogenized elastic coefficients with respect to the geometric parameters defining the microstructure, with special emphasis on the derivative with respect to the vectors defining the periodicity pattern. Section 4 extends these notions and results to locally periodic microstructures. Section 5 describes the optimization process, while in Section 6 some numerical examples are presented.
2. Periodic microstructures. As a preliminary for describing mathematically the notion of body having locally periodic microstructure, we briefly recall some notations and results on periodic microstructures. For more details, see [5].
Consider a parallelogram Y in R 2 which defines the periodicity of the microstructure. Often Y is taken to be the unit square for the sake of simplicity, but in the present work Y will be a general parallelogram defined by two linearly independent vectors g 1 , g 2 ∈ R 2 : A periodic microstructure is a body whose material coefficients (i.e., for the case of linear elasticity, its rigidity tensor) varies according to a periodic pattern C, which is a Y -periodic fourth-order tensor field C : R 2 → R 16 . Typically, but not necessarily, the pattern tensor field C takes only two values, modeling a periodic mixture between two given component materials (see Figure 1). A function φ is said to be Y -periodic if φ(y + g i ) = φ(y) for i ∈ {1, 2}.
The effective (macroscopic) behaviour of the body is described by the so-called homogenized elastic tensor, denoted by C H and defined in terms of solutions of elliptic partial differential equations subject to periodicity conditions. More precisely, for a given macroscopic strain (represented by a symmetric 2×2 matrix A), consider the problem known in homogenization theory as cellular problem. Here, represents the symmetric part of the gradient. The solution w A of the above problem depends (linearly) on the matrix A. Consequently, the matrix σ defined by representing the macroscopic stress associated to w A , depends linearly on the macroscopic strain A and this dependency defines the homogenized elastic tensor C H through C H A = σ, that is, In the above, |Y | denotes the area of the periodicity cell Y . An equivalent definition of C H can be given in terms of energy type products : for two given symmetric matrices A and B, one has Note that, for a function w A of the form w A (y) = Ay + φ A (y), with φ A Y -periodic, one has (see in [5] Lemma 1 and its consequences). This can be expressed in a more concise form by introducing the set LP of linear plus periodic functions defined as Functions w ∈ LP are characterized by the property The cellular problem (2) may be written in strain formulation as follows: where A is a prescribed strain matrix.
The case of porous structures, that is, infinitesimal perforations in a given material, can be treated in a similar way (see [5,Section 4]). Consider a model hole, which is a compact set T ⊂ Y (see Figure 2), where Y is the periodicity cell. The cellular problem describing the behaviour of this perforated material is defined on the perforated space, denoted by R 2 perf , and a Neumann boundary condition is imposed on the boundary of the model hole.
The direct generalization of the cellular problem (2) for porous materials is stated below (the base material C is now considered to be constant) : LOCALLY PERIODIC MICROSTRUCTURES

437
The homogenized tensor C H can be defined as Remark 1. Note that, through careful interpretetation, problem (7) makes sense even if the model hole T exits partially the periodicity cell Y , as long as it does not touch any of its translations T + k 1 g 1 + k 2 g 2 , k ∈ Z 2 , k = 0. This is important since, in the optimization process, the model hole often crosses the boundary of Y . 3. Derivatives of the homogenized coefficients. The homogenized elastic tensor C H depends, of course, on the geometric parameters defining the periodic structure under consideration : the shape and topology of the model hole T , as well as the periodicity pattern. In order to use gradient-based optimization algorithms, the derivatives of these dependencies must be computed. The aim of this section is the computation of the derivative of the homogenized tensor with respect to the periodicity pattern (at the end of the section we recall the fomulae for the shape and topology derivatives). This periodicity pattern is given by the aspect and orientation of the periodicity cell Y , or, equivalently, by the two vectors g 1 and g 2 defining Y , see formula (1). Our goal is to compute the derivative of the homogenized elastic coefficients (components of the tensor C H ) with respect to the vectors g 1 and g 2 .
We shall analyse the case of a mixture of materials, governed by the system (2) or, equivalently, (6). The case of porous materials, governed by the system (7), can be treated in a similar manner. Consider a linear application V (see Figure 5) that transforms the periodicity cell Y inỸ generated by the vectors g 1 = V ( g 1 ) and g 2 = V ( g 2 ). The corresponding fourth-order tensor fieldC verifiesC(V (y)) = C(y).
To a certain extent, this is similar to the approach in [13]. The main difference is that, in [13], the authors try to adapt the linear transformation in order for the homogenized coefficients of the deformed microstructure to match a previously computed elastic tensor; this adaption is done as a post-processing step. In the present paper, everything is done in one optimization process.
The homogenized coefficientsC H for the deformed configuration are defined by evaluating, for two given symmetric matrices A and B, the energy-type quantity wherew A andw B are the solutions of problem (6) in the deformed configuration. By changing the integration domain to Y , the above expression writes as and having in mind thatC(V (y)) = C(y) and |Ỹ | = |det∇V ||Y | = |detV ||Y |, By taking advantage of the symmetry properties of the elastic tensor C, the quantity CH A, B can be written in the form In order to express the above quantity in terms of functions defined on the original configuration Y , we introduce the functionsw A (y) =w A (V (y)) andw B (y) = w B (V (y)). One has thenw A (z) =w A (W (z)), with W = V -1 , and thus ∂w i We now turn to the study of the (infinitesimal) variations induced in the quantity CH A, B by infinitesimal variations of V . The prefix δ will be used to denote an infinitesimal variation of the quantity following it. Remark 3. Along the above computations, we have assumed implicitly that the linear application V is invertible, that is, det V = 0. This is a natural assumption; actually, in practice, V will be chosen as the sum between the identity matrix and a small variation δV . It is also possible to impose a constraint like det V = 1 (by including a projection step in the algorithm); by linearization, this constraint becomes trδV = 0. See, however, Remark 4 below.
Recall that the matrices A and B are fixed, as well as the elastic tensor C of the (undeformed) mixture. Thus : The symbol δ commutes with the partial derivatives, hence Expressing the variation of W = V -1 in terms of the variation of V is a mere exercise of calculus : We now localize the computations, that is, we compute the desired derivative around the point V = Id (the identity matrix). By doing this, W becomes also Id, and w A =w A = w A . We shall also switch the notation from δC H to D P C H (the derivative of C H with respect to periodicity). Thus : In the first parcel of the expression above, the integral can be seen as the product between C ijkl ∂w k B ∂y l , which is a zero-divergence vector field due to the cellular problem which is a gradient. The third parcel has a similar structure. Due to a compensated compactness result (e.g. Lemma 1.3.1 in [1] adapted to the particular case of periodic fields), these parcels (first and third) can be written as products of averages : We stress that the second and fourth parcels cannot be written in the same way because of the extraneous terms δV mj and δV pl , respectively. 1 We now apply two properties already refered in Section 2, equations (4) and (5). Note that, sincew A is a linear plus periodic function whose linear part is AV , δw A is a linear plus periodic function whose linear part is AδV . Similarly, the linear part of δw B is B δV , hence The above formula is ready for implementation, since it gives the derivative of the homogenized tensor in terms of computable quantities multiplying (components of) δV . By choosing the matrices A and B in a basis of symmetric matrices, one may write equation (8) in a compact form : where P is a sixth order tensor describing the sensitivity of C H with respect to the periodicity pattern. Two animations are available at [17], illustrating the process of optimizing properties of the homogenized elastic tensor by varying the periodicity cell.  That is, if one performs a homothety on the cell Y and on the hole T , the homogenized elastic tensor of the microstructure will remain unchanged, since the formulas defining it are based on an averaging operation; note the division by the area of Y in formulas (3), (4) and so on. Attending to this property, the algorithm rescales from time to time the cell Y in order to prevent extreme increase or decrease in the size of the cell. However, numerical tests show that this operation is not necessary, because the derivative (8) of the homogenized tensor with respect to the periodicity is also scale invariant and thus the algorithm does not change significantly the size of the cell Y . Actually, it suffices to observe that, in formula (8), if one takes δV to be a multiple of the identity matrix, the derivative becomes zero. This means that the steepest descent direction produced by an optimization algorithm will be such that trδV = 0 (δV orthogonal to the identity matrix). This is roughly equivalent to imposing det V = 1 (see also Remark 3), so Y will have almost constant area. However, this does not eliminate the possibility of extreme elongation of the cell (see Remark 5 below).
Remark 5. Nothing can be done to prevent the cell from becoming a long and thin rectangle, especially when the example is such that laminates are convenient. However, one can prevent the cell from becoming sharp-angled. Note that it is possible for different cells to define the same periodicity pattern, as illustrated in Figure 6. To prevent extreme sharpening of the angles of the periodicity cell, the algorithm "normalizes" the cell from time to time, that is, it chooses, among all possible equivalent cells, the one closest to a rectangle (actually, the one having the lowest perimeter). For instance, if faced with a cell like the one drawn in Figure 6b or Figure 6c, the algorithm will choose instead the cell drawn in Figure 6a, leaving the mesh unchanged.
We state here the formulae of the shape and topology derivatives; the interested reader may find more details in [5,Sections 5 and 6].  Figure 7. Then, the corresponding variations in the homogenized tensor are described by where ds denotes the superficial measure on ∂T , while λ and µ are the Lamé coefficients of the (constant) base material tensor C. See equation (32) in [5].
By choosing the matrices A and B in a basis of symmetric matrices, one may write equation (10) in a compact form: where S is a fourth order tensor describing the shape sensitivity of C H .
The sensitivity of C H with respect to topology variations, that is, with respect to the nucleation of an infinitesimal hole at an arbitrary location y ∈ Y \ T , is described by see equation (24) in [5]. See also [16]. This can be again expressed in a compact form where T is a fourth order tensor describing the topology sensitivity of C H .

4.
Locally periodic microstructures. We call "locally periodic" a microstructure which, in the neighbourhood of each point of the macroscopic domain, has a periodic character. This periodic microgeometry can vary from point to point (at the macroscopic level), as described in [7,Section 3]. Such a microstructure is sometimes called in the literature "graded structure", or "functionally graded material".
For mixtures of materials, this amounts to letting the pattern tensor field C depend on both x (the macroscopic variable) and y (the microscopic variable). Note that in [7] the periodicity cell was kept fixed and thus it did not depend on x.
In the present work the periodicity pattern is allowed to vary from point to point (at the macroscopic level) and is also varied along the optimization process. Remark 6. We do not focus here on describing with mathematical rigor the notion of locally periodic structure. The interested reader can find in [8, Chapter 1, Section 6] a description of a microstructure with fixed periodicity cell and variable microgeometry. In [9], a rigorous description is given for a microperforation with variable periodicity cell but identical (spherical) holes. The main goal of the present paper is to deal with bodies made of pieces of periodic microstructure (these pieces will be actually the macroscopic finite elements, as explained in Section 6). In each of these pieces, the microstructure is assumed to be truly infinitesimal and thus these pieces can be glued together with no difficulty.
The periodicity cell depends on the macroscopic variable x, and is denoted by Y (x). The cellular problem (2) depends on x as a parameter : The homogenized elastic tensor C H depends now on x and is defined by C(x, y) y (w A (x, y)), y (w B (x, y)) dy , where A and B are two arbitrary strain matrices. A locally periodic porous material is described in a similar manner. The model hole T now varies with x and shall be denoted by T (x).
The corresponding perforated plane will be denoted by R 2 perf (x) The base material C is now considered to be constant. The cellular problem (7) becomes The homogenized ellastic tensor C H (x) is defined by where A and B are two arbitrary strain matrices. The behaviour of the macroscopic body is described by the homogenized elastic tensor C H (x). More precisely, suppose that the macroscopic body occupies the domain Ω ⊂ R 2 and that the boundary of Ω is split in two parts (the Neumann part Γ N and the Dirichlet part Γ D ). The body is fixed on Γ D and is subject to the superficial force g on Γ N . Then, the state of deformation of the body is described by the solution u of Remark 7. One may wonder how microstructures with different periodicity patterns can be glued together in order to manufacture a macroscopic body with the desired properties. This is a pertinent question, but it is outside the scope of the present paper. See e.g. [13] for an approach to this question. The goal of the present work is to deal with macroscopic bodies made of "patches", see Remark 6.
Remark 8. Actually, the objection risen in Remark 7 applies, although to a lesser extent, even to the context of locally periodic microstructures with constant periodicity pattern, like the ones studied in [7]. Indeed, perforations of different shapes can sometimes be difficult to glue together. See, e.g., [15] The expressions of the derivatives of the homogenized elastic coefficients, obtained in Section 3 for periodic microstructures, are directly generalized for locally periodic microstructures. The only difference is that now a parameter x appears in the formulae; thus, equations (9), (11) and (13) become : 5. Optimization process. The aim of the present paper is to optimize macroscopic properties of the body under study, by varying the details of its locally periodic microstructure, more specifically, by varying the shape and topology of the holes T (x), as well as their periodicity pattern.
Consider an objective functional Φ depending on the solution u of problem (17). A typical example is the minimization of the compliance of the body

CRISTIAN BARBAROSIE AND ANCA-MARIA TOADER
Usually, there is also a constraint on the volume of material: where Remark 9. Note that, in (21), we have chosen the objective functional Φ to be half of the functional used in [7]. The reason is that the authors have learned meanwhile that definition (21) is more consistent with the mechanical notion of elastic stored energy. See also [4].
The optimization process passes through a chain of dependencies : that is, to a family of cellular holes T (x) one associates the homogenized tensor C H (x) as defined by (16), then the macroscopic elastic deformation u as defined by (17) and finally the value of the objective functional Φ. A simpler chain of dependencies holds for the volume V of material.
Describing the variation δΦ of the objective functional in terms of δC H and δu is relatively simple. The difficult part is to eliminate δu, thus expressing δΦ in terms of δC H only, and this is the object of the adjoint method. However, for the particular case of the compliance functional defined in (21), the problem turns out to be self-adjoint and one obtains See [7,Section 4] for details, taking into account Remark 9 above. As for the first part of the chain (23), equations (18), (19) and (20) describe the variation δC H associated to an infinitesimal variation in each of the structural parameters : shape of the holes/inclusions, their topology, and their periodic arrangement, respectively.
An optimization algorithm is used which alternates shape, topology and periodicity optimization at the cellular level (the approach of alternating different directions of optimization was proposed, e.g., in [2]). The macroscopic domain Ω is divided into finite elements (rectangular ones in our approach, but this is not essential). In each macroscopic finite element, the microgeometry is supposed to be constant. Consequently, the homogenized elastic tensor C H and the material density θ are also constant in each macroscopic finite element.
The cellular geometry describing the periodic microstructure in each macroscopic finite element is discretized using a (triangular) finite element mesh on the respective periodicity cell. Some of these triangles are filled with an isotropic elastic material of Lamé constants λ and µ while others are void, thus defining the model hole. This mesh is used for solving numerically the cellular problem (15). The periodicity conditions in (15) are implemented by identifying the opposite sides of Y and by keeping track of the linear part A of w A . For more details, see [7], [6] and [3].
The computer implementation has three components : a FORTRAN program which deals with each of the cellular meshes and optimizes the microgeometry, a C++ program (using the libMesh library) which deals with the mascroscopic analysis, and a Python script which links these two programs and provides data exchange. The computation is strongly parallelized : the Python script connects to several machines and spreads cellular optimization processes, by using a regular internet link and the protocol ssh. For more details, see [7, Section 5]. 6. Numerical examples. In order to illustrate the method, results for three examples are shown. A cantilever is considered, ocuppying a rectangle of size 1.5 × 1, clamped on its left side and subject to loads on its right side. This rectangle is divided in 600 = 30 × 20 Lagrange finite elements of type Q9. 2 The base elastic material C is taken to have Young modulus E = 1 and Poisson coefficient ν = 0.3. More specifically, with µ = 0.38461538 and λ = 0.576923.
Remark 10. This value of Λ is half of the Lagrange multiplier used in [7,Section 6]. This choice facilitates the comparison of the numerical results, see Remark 9 above.
The 600 cellular meshes (one per macroscopic element), are composed of triangular Lagrange P 1 finite elements. The number of elements varies during the optimization process between 1000 and 2000 triangles per cellular mesh. Starting with this initial structure, the microgeometry is varied. Shape, topology and periodicity optimization steps are alternated during the optimization process.

No penalization of intermediate densities is applied.
Each numerical example is illustrated by two Figures. Firstly, Figures 9, 11  and 13 show the configuration of the body, the applied loads, and three graphics  [7] (these values have been rescaled by a factor of 0.5, in accordance with Remarks 9 and 10 above). Secondly, in Figures 10, 12 and 14 the optimized structure is presented, following the same conventions as in Figure  8 : in the center, the density of material, θ, is drawn as a function defined in the macroscopic domain, using levels of grey. Around this central zone, magnifications of the microstructure in 8 chosen macroscopic finite elements are shown. In the background of each zoom, the periodicity cell is drawn in grey.
Note that it is difficult to represent graphically the homogenized tensor C H of the optimized structure, since it involves 6 scalar functions. An interactive version of Figures 10, 12 and 14 is available at [18].
In the first example, a load of intensity 1 is applied in the middle point of the right side of the cantilever. Figure 9 shows the setting of the problem and the evolution of the Lagrangian Φ + ΛV , of the compliance Φ and of the volume of material V , as functions of the number of iterations.
In the graphic of the Lagrangian, the continuous line shows the evolution of the Lagrangian when varying the periodicity pattern, while the dotted line shows the evolution in the previous approach, with constant periodicity cell, as described in [7] (these values have been rescaled by a factor of 0.5, in accordance with Remarks 9 and 10 above). One can see that the present approach is more efficient, reaching a significantly lower value of the Lagrangian. Note that there are intervals where the Lagrangian is almost constant, and the volume of material is truly constant. These are the optimization steps when the periodicity pattern is varied. It is interesting Figure 10. Cantilever with one load in the middle, optimized structure that these optimization steps, which at first view do not contribute for decreasing significantly the value of the Lagrangian, actually adjust the structure's periodicity pattern and help the subsequent optimization steps (mainly shape optimization) to decrease the value of the Lagrangian more efficiently. Figure 10 shows the density of material, in grayscale, and eight zooms of the optimized microstructure in chosen points of the body. Note how very different material densities are obtained : from a full material ( Figure 10, zoom a) to an almost empty cell ( Figure 10, zoom f ). Note also the great variety of microgeometries produced by the algorithm. In Figure 10, zooms c and e, we see structures close to rank-1 laminates, like those already obtained in [7]; however, now the laminates have much more freedom to orient themselves according to the stress in the macroscopic body. For instance, the microstructures represented in Figure 10, zooms d and e, exhibit a privileged direction having a slight slope, wich would be impossible to obtain if the periodicity cell were kept constant (square).
No symmetry is imposed in either of the examples. The results depicted in Figure  10 present a rather good symmetry, as expected. The symmetry can be observed even at the microscopic level, as shown in Figure 10, zooms g and h.
In the second example, a load of intensity 1 is applied in the lower right corner of the rectangle, see Figure 11. The optimized structure is represented in Figure 12. Some of the microstructures resemble rank-1 laminates (Figure 12, zooms b and c), others seem to mimic rank-2 laminates (Figure 12, zooms e, f and g), while others are more difficult to classify ( Figure 12, zooms a, d and h).
The approach presented here is not limited to the minimization of the compliance: other, more complicated, functionals can be dealt with. The third example involves a multi-load situation. Three loads are considered, as shown in Figure 13, acting independently of each other. The objective functional is the average of the three compliances (for the three load cases) : Results are presented in Figure 14. The optimized microstructures are somewhat more elaborated than in the previous two examples. Fewer laminates are obtained.
These microstructures are required to resist well to three different states of (macroscopic) stress.
In all examples, one notices more intermediate grayscales when comparing with the respective results in [7]. Most importantly, a significantly lower value of the Lagrangian is obtained.
In some cases, the extreme flexibility of the periodicity cell becomes apparent, like in Figure 14, zoom h, which is a good illustration of Remark 5.  Figure 11. Cantilever with load in the lower corner, problem setting and convergence history Figure 12. Cantilever with load in the lower corner, optimized structure Remark 11. The topology optimization steps do not contribute significantly to the optimization process. The authors have observed that, quite often, during a topology optimization step, a small hole is created and is later destroyed during subsequent shape optimization steps. Also, in the final, optimized, structures, it is  Figure 13. Cantilever with three independent loads, problem setting and convergence history Figure 14. Cantilever with three independent loads, optimized structure very rare to see more than one hole in the periodicity cell. 3 The only exception is shown in Figure 14, zoom b, where a very small hole can be seen. If the optimization process were continued, it is likely that this hole would have been destroyed by subsequent shape optimization steps.

7.
Conclusions and future work. The present paper presents an algorithm for the optimization of bodies having locally periodic perforations. Shape, topology and periodicity optimization steps are performed following an alternate directions approach. The problem in numerically heavy (hundreds of cellullar meshes are optimized, each having thousands of finite elements). The implementation is massively parallel, in order to alleviate the computational burden.
Our approach is related to free material optimization in the sense that it uses the derivative of the objective functional with respect to the homogenized elastic coefficients. The method is general : any objective functional can be treated, as long as its derivative with respect to the macroscopic material coefficients can be computed.
The numerical results are encouraging, showing good agreement with results from the literature.
The upgrade of the algorithm to three-dimensional problems is the object of ongoing work. The main difficulty is the implementation of the finite element mesh on the cube with its opposite faces identified (which is equivalent to meshing the threedimensional torus), especially its deformation and regeneration. Robust algorithms for mesh deformation and mesh regeneration in R 3 are difficult to find.
Other directions for improving the algorithm in the future are: the implementation of a quasi-Newton algorithm in order to accelerate the convergence and the implementation of different optimization criteria. For instance, it could be interesting to impose an upper bound on the (microscopic) pointwise stress.