Shape optimization for composite materials in linear elasticity

This article is devoted to the optimal design of the microstructure in composite materials, which are governed by the equations of linear elasticity. To this end, we combine homogenization with shape optimization. In particular, we determine the sensitivity of the homogenized coefficients of the elasticity tensor with respect to the shape of the periodic microstructure also in case of spatially varying material coefficients. We compute the respective Hadamard shape gradient and demonstrate the applicability and feasibility of our approach by numerical experiments for different problem settings.


Introduction
The equations describing the behaviour of materials with elastic properties are widely studied in mechanics. Especially, the simplified case of linear elasticity is of interest. The system of equations considered in this article describes the behaviour of an elastic material subject to surface loads. Our material is considered to be periodically heterogeneous. Specifically, on a microscopic scale, one substance envelops another substance in a periodically repeating pattern. These kind of materials can be produced amongst other techniques by additive manufacturing and find applications in many fields, for example as light weight structures in aeronautics or as synthetic 1 3 bone material in medical applications, see e.g. Hopkinson et al. (2006), Wang et al. (2016) and the references therein.
We aim in this article at optimizing the shapes in the periodically repeating pattern on the microscopic scale. Therefore, an objective function depending on the shapes in the pattern is defined. It describes a desired macroscopic property, which is achieved when the function's minimal solution is found. In order to solve this minimization problem, methods of shape sensitivity analysis are employed. Thus, we are able to compute Hadamard's shape gradient, which allows for the application of gradient based optimization methods to find the optimal shape, compare with Delfour and Zolesio (2001), Murat and Simon (1976), Pironneau (1984), and Sokolowski and Zolesio (1992) for example. In Fig. 1, a schematic of an optimal design problem is displayed. It aims at finding a microstructure of the material that minimizes the compliance of a pillar with a force acting on top of it.
Note that the optimal design of periodic structures has already been considered in many works, see Adachi et al. (2006), Ferrer et al. (2018), Hollister et al. (2002), Lin et al. (2004), Luo et al. (2017), Sigmund (1995), Wang and Kang (2019), and Wormser et al. (2017) for some of the respective results and applications. The methodology used therein is mainly based on the forward simulation of the material properties of a given microstructure. In contrast, in Allaire et al. (2019), Barbarosie (2003), Dambrine and Harbrecht (2020), Faure et al. (2017), Geoffroy-Donders et al. (2020), Haslinger and Dvořák (1995), Hübner et al. (2019), Wang et al. (2004), shape sensitivity analysis has been used. The shape derivative of the coefficients of the homogenized tensor, expressed as volume integrals, have been computed in Allaire et al. (2019), Barbarosie (2003) in case of piecewise constant material properties. Note that these shape derivatives have also been derived in the context of a level set representation of the inclusion in Nika and Constantinescu (2019). The novelty of this work is the use of Reynolds transport theorem to compute the shape derivative. Especially, we characterize the local shape derivative and allow for locally varying material coefficients, which was not done before. Note that the local shape derivative is an important tool in shape uncertainty quantification, compare Harbrecht and Li (2013) or Harbrecht et al. (2008) for example. Therefore, we follow the ideas and computations presented in Dambrine and Harbrecht (2020) and Harbrecht et al. (2022) for an elliptic diffusion problem. We emphasize that the numerical results presented in this article are in two spatial dimensions. They serve as a proof Shape optimization for composite materials in linear… of concept and foundation for future works, where we will consider the optimization of three-dimensional problems.
This article is structured as follows. In Sect. 2, we introduce the system of equations and objects we treat in this article. Especially, we briefly discuss the notation of the computations with tensors. Moreover, homogenization through a multiscale ansatz is performed. The result of the homogenization procedure is a system of equations describing a homogeneous approximation to our material. This new system contains the homogenized elasticity tensor, also called the effective elasticity tensor. Next, the objective functions under consideration are introduced. From then on, we consider a composite material, which consists of two intertwined materials with different elastic properties.
In Sect. 3, the shape derivative of the coefficients of the effective tensor is computed. To this end, we first compute the local shape derivative of the cell functions, then the shape derivative of the material coefficients, and finally the shape derivative of the least squares and compliance functional under consideration. These steps are repeated for the special case of a partly hollow material. In the two-dimensional setting, this corresponds to the special case of a perforated plate, while in the threedimensional setting, we arrive at a scaffold structure, a partly hollow three-dimensional structure.
The implementation in the two-dimensional setting is presented in Sect. 4. We use a parametrization by a finite Fourier series to represent the sought interface between the two materials in the microcell. Then, we apply parametric finite elements to discretize the state equation under consideration. Finite elements are also used in the macroscopic domain. By means of finite differences we verify our theoretical findings and the implementation.
The following two sections are dedicated to numerical experiments. In Sect. 5, we optimize for various material parameters the shape of the interface between the two given materials. Especially, we consider nonhomogeneous material parameters, different initial guesses, and also the situation where there is a void instead of an interior material. In Sect. 6, the compliance of a porous material is minimized by finding an optimal void within its microstructure. Finally, in Sect. 7, we state concluding remarks.

Linear elasticity system
Consider a fixed, bounded domain Ω ⊂ ℝ n , which is the reference configuration of an elastic body, subject to a body force f . The structure Ω is clamped on Γ D ⊂ Ω and subject to surface loads g on Γ N ⊂ Ω . Here, Γ D and Γ N are assumed to be fixed such that Γ D ∩ Γ N = � . The displacement u is then the solution of the system where the Cauchy stress tensor for the displacement u is defined as with the two Lamé constants > 0 and > 0 and the identity matrix id. This system of equations corresponds to the system introduced in Braess (2001).
Next, we shall introduce the fourth order elasticity tensor, which allows to express the Cauchy stress tensor in terms of a tensor and of the displacement. To this end, we first introduce some notation for tensors.

Tensor notation
The tensor notation applied in this article follows in wide parts Itskov (2007). Let n ∈ ℕ , a, b, c, d be column vectors in ℝ n , B, C, D matrices in ℝ n×n and id the respective identity matrix. The Frobenius product of two matrices is defined as The tensor product for vectors is defined by the relation Whereas, in case of matrices, we have the two tensor products There hold the following two relations for matrices which are given by tensor products of vectors: Moreover, we use the classical transposition operation for matrices, while for fourth order tensors there is the transposition operation The Frobenius product of a fourth-order tensor A and a matrix B is computed as

Shape optimization for composite materials in linear…
where e i , i = 1, … , n , denotes the canonical basis of vectors.
With the Lamé constants and , the fourth-order elasticity tensor A , satisfying the relation A ∶ ∇u = (u) , can be expressed as Especially, with the tensor computation rules above, we find which for B∶=∇u corresponds to the Cauchy stress tensor in (2.2).

Remark 1
In some works on linear elasticity, for instance Braess (2001) and Geoffroy-Donders (2018), the elasticity tensor is defined differently. The notation from Itskov (2007) was chosen in this article since it contains all the objects used in the following pages in a single work. Furthermore, the super-symmetric identity tensor which is introduced in Itskov (2007) allows to combine the elasticity and strain tensor, often denoted by A and , respectively, resulting in the elasticity tensor expression (2.3).
Instead of defining the elasticity tensor by means of the Lamé constants and , we can define it through the elasticity modulus E > 0 and the Poisson ratio 0 < < 1∕2 . The elasticity modulus is proportional to the resistance a material opposes to its deformation. The Poisson ratio is a description of the influence of stresses on displacements perpendicular to the specific direction of loading. For soft materials such as rubber it is near 0.5 and for typical solids around 0.2-0.3. The Lamé constants are derived from these two numbers by the formulas A fourth-order tensor A is called to be uniformly elliptic if there exists a constant > 0 such that This property is always satisfied for the elasticity tensor A defined by E and . Thus, we denote the admissible elasticity tensor by the image of the function

The multiscale ansatz and composite materials
The homogenization procedure considered in this article follows closely Allaire (2002). Homogenization is concerned with the situation where the elasticity tensor A is oscillatory. To this end, we denote the unit cell by Y = [0, 1] n and write y∶=x∕ for any x ∈ Ω and > 0 . We assume A(y) to be Y-periodic, satisfying A(y) = M E(y), (y) for all y ∈ Y with E(y) > 0 and 0 < (y) < 1∕2 . Then, in order to approximate the displacement in (2.1), one makes the multiscale ansatz where u j (x, y) are Y-periodic functions in y for all j = 0, 1, … . Note that this ansatz does not take into account the boundary conditions of (2.1) and extra terms would be needed to capture boundary layers, see e.g. Cioranescu and Saint Jean Paulin (1999).
Following Allaire (2002), this ansatz shows that the solution u of the equation can be approximated for 0 < ≪ 1 by the solution u of the homogenized equation To characterize the effective tensor A * , we introduce the cell problems where e ij ∶=e i ⊗ e j , i, j = 1, … , n , denotes the basis of matrices and is the space of Y-perioric H 1 -functions, equipped with the norm ‖ ⋅ ‖ H 1 (Y) . Then, for all i, j, k, l = 1, … , n , the entries of the effective tensor are given by A composite material consists of two materials with different elastic properties, see Fig. 1 for an illustration. We model such a material by defining the elastic tensor where is an open and connected subset of Y. Under the assumption that A 1 (y) = M E 1 (y), 1 (y) and A 2 (y) = M E 2 (y), 2 (y) for all y ∈ Y and E 1 , E 2 , 1 , 2 ∈ C 1 (Y) , this is a piecewise smooth, fourth-order elasticity tensor for all choices of .
The normal vector n on indicates the direction from the interior of to the exterior Y⧵ . The jump of a function f through the interface at a point y ∈ is thus given by y) ∶ e il + ∇w il (y) ∶ e jk + ∇w jk (y) dy.

3
Shape optimization for composite materials in linear… For composite materials, the cell problem (2.7) is augmented with the two transmission conditions on the interface With these conditions, the variational formulation of the cell problem (2.7) for v ∈ [H 1 per (Y)] n becomes Formula (2.11) together with (2.8) shows the dependency of the homogenized elasticity tensor A * ( ) on the choice of the interface . Also, A * ( ) inherits the symmetry properties and the uniform ellipticity from the elasticity tensor A(y) . Therefore, the homogenized system on Ω , written in the form with the (effective) Cauchy stress tensor is * ( , u)∶=A * ( ) ∶ ∇u , has a unique solution for all sufficiently smooth subsets .

Goal of optimization
In this article, we consider two objective functions. For a given material, we aim to minimize one of the two objective functions, which depend on the inclusions present within the repeating pattern in the microstructure of those materials.
Firstly, we ask for the shape ⊂ Y , such that the least squares matching between the effective tensor A * ( ) and a desired elasticity tensor B ∈ ℝ n 4 with coefficients b ijkl is minimized with respect to the Frobenius norm ‖ ⋅ ‖ F . Note that not all desired elasticity tensors B can be attained, since they might not satisfy the Hashin-Shtrikman bounds, see Hashin and Shtrikman (1963). Secondly, we consider an elastic body described by Ω . In the microstructure of Ω , we ask for the shape ⊂ Y , which minimizes the compliance functional where f and g are forces given in (2.12).
In a porous material, the compliance is minimized by simply filling the holes. In order to exclude this trivial solution, we impose a volume constraint | | = 0 on the (2.10) w ij (y) = 0 and A(y) ∶ e ij + ∇w ij (y) n(y) = 0 for y ∈ .
(2.11) ∫ Y A(y) ∶ (e ij + ∇w ij ) ∶ ∇v dy = 0. (2.12) void. To this end, the objective function considered in order to minimize the compliance of the elastic body Ω under the given volume constraint is where ∈ ℝ is the Langrange multiplier to be maximized and ∈ ℝ + is an appropriately chosen constant. This augmented Lagrangian functional allows to enforce the volume constraint without taking → ∞ throughout the optimization process.
Only the Lagrange multiplier is updated by whenever J C ( , ) has been minimized with respect to , cf. Nocedal and Wright (2006). In order to optimize the objective functions above by a gradient based optimization method, we shall compute the shape derivative of J LS ( ) and J C ( , ) . We begin by finding an expression for the shape derivative of the coefficients of the effective tensor A * ( ).

Local shape derivative
We introduce a vector field h ∶ Y → Y that vanishes on the boundary Y and consider the perturbation of identity operator T t ∶=I + th . T t is a diffeomorphism for t small enough, which preserves Y. We define the function where w ij is the unique solution of the cell problem (2.7) and y ij is a column vector with entries 0 except in the i-th row the entry y j . Then, ∇y ij = e ij and the cell problem (2.7) becomes (2.15)

Shape optimization for composite materials in linear…
We present next the following theorem, which specifies the system of equations the local shape derivative ij satisfies on Y. As the local shape derivative has a jump at , we introduce the broken H 1 -space in accordance with Theorem 1 (Shape derivative of the cell problem) The shape derivative ij in X of the function ij = w ij + y ij is the solution in Y of the transmission problem

Proof
The proof is carried out in the usual way, compare with Dambrine and Harbrecht (2020). We prove that the material derivative exists, characterize it, and finally express the shape derivative.
1. Computation of the material derivative. Let t ij (y) be the unique solution of the cell problem (2.7) on the perturbed domain Y t ∶={y + th(y) | | y ∈ Y} and set Then, the transported function ̃ ij (t, y)∶= t ij T t (y) , satisfies the variational equation Using the property A ∶ (BC) = (AC ⊺ ) ∶ B of the Frobenius product for square matrices A, B, C and subtracting the Eq. (2.11) from the Eq. (3.4) above, we find for any t > 0 n . Now, to this equation we add and subtract equal terms in order to obtain the equivalent equation Replacing The details of these arguments can be found by straightforward adaption of the proof for the case of conductivity, see Henrot and Privat (2010).
With the reasoning above, the material derivative exists and satisfies the variational problem: find ̇ ij ∈ [H 1 per (Y)] n such that In this equation, for m = 1, 2 , h(∇ ⊗ A m ) denotes the fourth order tensor with entries ⟨h(y), ∇a m ijkl (y)⟩ , where a m ijkl are the entries of A m (y). 2. Recovering the shape derivative. Employing integration by parts in (3.6), we can rewrite the weak problem solved by the material derivative as Shape optimization for composite materials in linear… From Delfour and Zolesio (2001), we have the relation ̇ ij = ij + ∇ ij h between the shape derivative and the material derivative. Thus, we can replace the material derivative in (3.6) by the shape derivative to arrive at Since ̇ ij has no jump at as a function in [H 1 per (Y)] n , we conclude the first transmission condition for ij at the interface : Next, with (2.10), (3.2), and the symmetry of the matrix A ∶ ∇ ij , we rewrite (3.7) as We proceed with integration by parts over and Y⧵ on the left-hand side and with integration by parts over the closed surface on the right-hand side of the equation. By using (2.10) again, this yields We deduce that ij is harmonic in both subdomains and Y⧵ , while the second transmission condition is also deduced from this equation. ◻

Sensitivity of the effective tensor
From Delfour and Zolesio (2001) or Sokolowski and Zolesio (1992), given a (3.7) This formula appears in the literature in Hadamard's shape calculus and also as the Reynolds transport theorem in continuum mechanics.
Theorem 2 (Shape derivative of the coefficients of the effective tensor) The shape derivatives of the entries a * ijkl of the effective tensor given by (2.8) are where ij (y) = w ij (y) + y ij with w ij ∈ [H 1 per (Y)] n are the unique solutions of the cell problem (2.7) for all i, j = 1, … , n.
Proof We set f (t, y) = A(y) ∶ ∇ t ij (y) ∶ ∇ t kl (y) . Since the value of the function ij has a jump through , we want to consider ∇ ij inside and outside the interface separately. Therefore, we introduce the set Θ , which may designate either the set Y⧵ or . Then, for the first term in the Reynolds transport theorem (3.9), we find Hence, we arrive at where we added the equation with Θ = Y⧵ to the one with Θ = and used the results from Theorem 1.
We compute next the second term in the Reynolds transport theorem (3.9) by using the divergence theorem, the identity P ∶ � Q(n ⊗ n) � = tr � Pnn ⊺ Q ⊺ � = ⟨Pn, Qn⟩ for matrices P, Q ∈ ℝ n×n , and the second transmission condition in (2.10) Putting the two terms together yields the desired result of the theorem. ◻

Remark 2
1. By substituting back ij = w ij + y ij , we arrive at the computable expression This condition is used in the proofs of Theorems 1 and 2. Without using the second transmission condition, the shape derivative of the coefficients of the homogenized tensor can be written as Formula (3.11) has the advantage of explicitly preserving the symmetry of the homogenized tensor. In contrast, using formula (3.10) together the symmetrization ( a * iklj + a * kijl )∕2 reestablishes the symmetry in the implementation while omitting terms being equal to 0. In our numerical examples, this variant turned out to be more accurate.

Shape derivative of the objective functions
For the least squares matching, with the help of Theorem 2 and the chain rule we can easily determine the shape derivative of the objective J LS ( ) given by (2.13).

Proposition 1 The shape derivative of the objective function J LS ( ) from (2.13) reads
For the computation of the shape derivative of the objective function J C ( ) , we use the notation (3.10) (3.12) and where t ≥ 0 and u( , ⋅) is the displacement function in case of the interface , being the solution of (2.12).
The following two lemmata are required to express the shape derivative of the compliance functional (2.14).
Lemma 1 For t > 0 , a vector field h and ⊂ Y , it holds for all v ∈ [H 1 (Ω)] n Proof Expressing the variational formulation of (2.12) once for and once for t , then subtracting one expression from the other and dividing by t, we obtain for any v ∈ [H 1 (Ω)] n Taking the limit t → 0 leads to and the claim follows. ◻ Lemma 2 For a vector field h , the shape derivative of the compliance (2.14) reads

Proof We find
This expression converges for t → 0 to

3
Shape optimization for composite materials in linear… The first two terms in this result are equal due to the symmetry properties of A * . By using Lemma 1, we thus conclude the desired result.
Combining the shape derivative of the volume of the inclusion with Lemma 2 and the chain rule, we obtain the following result.

Perforated plates and scaffold structures
We consider now the n-dimensional unit cube Y = [0, 1] n with a hole such that ⋐ Y . This corresponds to the limiting case A 2 → 0 in (2.9) for a composite material. In two spatial dimensions, we get perforated plates, while in three spatial dimensions, we get scaffold structures.
The previous results remain similar when considering Y⧵ instead of Y. This means that the cell problem becomes: seek w ij ∈ [H 1 per (Y⧵ )] n with (w ij , 1) [L 2 (Y⧵ )] n such that The homogenized equation in the present situation is while the coefficients of the effective tensor are Moreover, the system for the shape derivative reads This implies that the shape derivative of the coefficients of the homogenized tensor is given by Consequently, the shape derivative of the objective function from (2.13) reads Furthermore, the shape derivative J C ( , )[h] of the objective function J C ( , ) from (2.15) remains identical to (3.14), except for the entries of the fourth order tensor A * ( )[h] which are now given by the formula (3.15).

Implementation using finite elements
Our implementation is for the two-dimensional case, i.e. n = 2 , and was implemented in Matlab.
In the microstructure of the material, we assume that the sought domain is starlike with respect to the midpoint of the unit cell. The boundary of the interface is represented in polar coordinates by the function where the radial function r( ) is given by the finite Fourier series where N ∈ ℕ and a −N , … , a 0 , … , a N ∈ ℝ 2N+1 .
The domain Y = [0, 1] × [0, 1] is discretized by parametric finite elements with a macro triangulation which is composed of 8 triangles and 8 squares in order to resolve the interface exactly. Each of these elements is transformed from a reference geometry by the help of a diffeomorphism into Y. The squares are generated by Coons patches and the triangles with the Zenisek's formula from Zenisek (1990). The mesh on the reference triangle is composed of 4 triangular elements while the mesh of the reference square is composed 2 × 4 triangular elements. Therefore, on Y, we place 24 × 4 For this modified problem, we use piecewise linear finite elements to find an approximation to the solution of the cell problem (2.7). Since the interface is resolved by the triangulation, the convergence order of our implementation is of second order in the mesh size h with respect to the [L 2 (Y)] 2 -norm. At the macroscopic scale, we consider a rectangular elastic body Ω , i.e. a pillar of a certain height and width. At its bottom, homogeneous Dirichlet boundaries are imposed, which means that the pillar is fixed at its bottom edge. On its top edge, a constant force is applied, expressed through a nonhomogeneous Neumann boundary condition and the function g . The remaining two boundaries can move freely in space which corresponds to homogeneous Neumann boundary conditions at these boundaries. Moreover, there is no volume force, i.e., f = . This setting is implemented analogously to Alberty et al. (2002). Piecewise linear finite elements are used to solve the respective variational formulation.

Verification of the shape derivative through finite differences
Choose N ∈ ℕ and let the coefficients a −N , … , a 0 , … , a N ∈ ℝ 2N+1 be describing the interface in our numerical experiments in accordance with (4.1) and (4.2). We use the computable and symmetric expression ( a * ijkl ( )[h] + a * jilk ( )[h])∕2 from Remark 2 to evaluate the shape derivative of the coefficients of the homogenized tensor. We evaluate this expression in the midpoint of the edge lying on the interface of the triangular elements. As deformation functions, we use the functions Thus, we obtain the shape derivative of the coefficients of the homegenized elasticity tensor with respect to the 2N + 1 directions q k . We compare these 2N + 1 results with an approximation of the shape derivative obtained by first order finite differences. We compute these from the expression (2.8) of the coefficients of the homogenized tensor a * ijkl , once evaluated on Y with the interface given by the function r( ) from (4.2) and once evaluated on Y with the interface given by the function r( ) + q k ( ) . For our test runs, we set = 10 −8 .
The material properties are chosen as follows. To realize a composite material, we define the elasticity tensor outside of the interface by the values E = 1 Pa and = 0.3 and inside of the interface by E = 2 Pa and = 0.2 . Whereas, to obtain a perforated plate, we use E = 1 Pa and = 0.3 outside of the interface and no material inside of the interface, cf. (2.4). The boundary under consideration is a circle of radius 1/4 on one hand and a respective perturbed version of this circle on the other hand, the same perturbed circle as shown in the mesh plots in Fig. 2.
In Fig. 3, the convergence of the evaluation of the formula (2.8) towards the respective finite difference can be observed in all of the four cases. Therefore, we can conclude that our implementation is correct.

Computation of the volume of the inclusion
In our implementation, the definition of the inclusion by the functions and r from (4.1) and (4.2), respectively, allows for a simple computation of the volume of the inclusion by the formula

Shape optimization for composite materials in linear…
With this formula, the computation (3.13) and the functions q k from (4.4), the shape derivative of the volume of the inclusion is Therefore, the shape derivative of the function (2.15) is expressed by The entries of the fourth order tensor A * ( )[q k ] are given by (3.15) in case of a porous material. (4.5)

Fig. 3
Sum of the 2 -norms of the absolute and the relative differences between the computable formula (3.10) and the finite difference approximation of the shape derivative in four different settings. The interface is either circular or a slightly perturbed circle, composite and hollow materials are considered, and the mesh refinement levels of the finite element method range from 0 to 7. As a means of comparison, the plots of linear and quadratic decay are shown where needed 1 3

Numerical results: least squares matching
We use the expressions (3.12) and (3.16) in order to compute the shape derivative of the objective function (2.13) inside an implementation of the gradient descent method with quadratic line search. The aim of the procedure is to find optimal interface shapes which minimize the objective function (2.13) for different settings.
In the examples below, we choose N = 32 in (4.2). The algorithm was always stopped after 75 gradient steps. The results of five settings within a composite material and one within a perforated plate are presented. The computation of these examples is based on 98 304 and 65 536 finite elements for composite materials and for perforated plates, respectively, which corresponds to the mesh refinement level = 6 in our implementation. In the following five settings, the materials under consideration are defined by E 1 ∶=2 Pa and 1 ∶=0.3 outside of the interface, and in the case of composite materials as E 2 ∶=5 Pa and 2 ∶=0.2 inside the interface. The other material parameters are then specified in the respective setting.

First example: anisotropic tensor for a composite material
In our first experiment, we compute the homogenized tensor A for a composite material, where the interface is described by the circle of radius 1/4. Then, for the desired homogenized elastictiy tensor B , we introduce an anisotropy into A by adapting the two values b 1111 ∶=a 1111 + ∕2 and b 2222 ∶=a 2222 − ∕2 . With the starting profile chosen to be the circle of radius 1/4 , we apply the gradient descent method for different values of . The optimized shapes we found are shown in Fig. 4. The histories of the objective function and of the 2 -norm of the shape gradient during the optimization process are displayed in Fig. 5.

Second example: anisotropic tensor for a perforated plate
In this example, we consider perforated plates. We proceed similarly as in the first example to construct the desired homogenized elasticity tensor B . We compute the homogenized elasticity tensor A for the circle with radius 1/4 and introduce an anisotropy by setting b 1111 = a 1111 + ∕2 and b 2222 = a 2222 − ∕2 . Again, starting at the circle with radius 1/4, we then apply the gradient descent method to compute the optimal shapes for different values of . The optimized shapes we found are shown in Fig. 6. The histories of the objective function and of the 2 -norm of the shape gradient during the optimization process are displayed in Fig. 7.

Third example: various starting profiles for one anisotropic tensor
In our third experiment, we are concerned with the uniqueness of optimal shapes for the shape optimization problem under consideration. We consider the desired homogenized elasticity tensor to be anisotropic, namely the desired homogenized tensor from Sect. 5.1 with = −0.04 . We choose our initial guess to be the 1 3 Shape optimization for composite materials in linear… circle of radius 1/4 for the first run, an axes-aligned ellipse for the second run, its rotated version for the third run, and finally we randomly perturb the circle for the last run, compare with the first row in Fig. 8. We apply the gradient descent method to optimize the shape for these various initial guesses. The results are shown in the second row of Fig. 8. The histories of the objective function and of the 2 -norm of the shape gradient during the optimization process are displayed in Fig. 9. This example indicates strongly that the solution is not unique as all optimized shapes exhibit both, small values of the objective function and small norms of the discrete shape gradient. Indeed, the same has already been observed in Dambrine and Harbrecht (2020) in case of a second order diffusion equation.

Fourth example: inhomogeneous material inside a composite material
In the fourth example, we consider the elasticity modulus E(y) and the Poisson ratio (y) to be nonconstant functions in the interior of the interface, namely

Shape optimization for composite materials in linear…
In the exterior of the interface, we keep the same homogeneous material properties as in the previous examples, given by E 1 and 1 . We moreover use the anisotropic desired homogenized tensor from Sect. 5.1 with = −0.04 as desired homogenized elasticity tensor B . We take the four different initial guesses from Sect. 5.3 and apply the gradient descent method. We obtain for these initial guesses the optimized shapes shown in Fig. 10. The histories of the objective function and of the 2 -norm of the shape gradient during the optimization process are displayed in Fig. 11. Again, the optimal shape seems not to be unique.

Fifth example: convergence towards an inverted composite material
In this example, the desired homogenized tensor represents a composite material which has material properties E 1 and 1 inside and E 2 and 2 outside of the interface. This correponds to a setting where the two materials of the composite materials we considered yet have switched their position. In order to define the desired homogenized elasticity tensor B , we consider a circle with radius of 0.48. This ensures that an attainable desired tensor is derived. We use the same initial guesses as in Sect. 5.3, which yields the optimized shapes found in Fig. 12. The histories of the objective function and of the 2 -norm of the shape gradient during the optimization process are displayed in Fig. 13.
(5.1) E(y)∶= E 2 2 1 + sin( y 1 ) + sin( y 2 ) and (y)∶= 2 e −(2y 1 −1) 2 −(2y 2 −1) 2 . We shall next study cases for which we want to minimize the compliance of an object. In our first example, we consider a pillar with holes in its microstructure, subject to a force on its top. In our second example, we consider a beam, attached to the wall on one side and with a force acting from atop on the opposite side. The beam is composed of inhomogeneous materials in its microstructure.
In our implementation, we use the energy functional (4.5) in order to compute the shape derivative of the objective function (2.15) of the gradient descent method. In each step of the optimization process, we need to approximate the macroscopic solution u of the system of equations (2.12), describing the displacement of the material within the pillar or the beam, respectively, subject to the force g , as well as the solutions w ij of the cell problems (2.7) at the microscopic scale. For both examples, we choose N = 32 in the Fourier series (4.2). The macrostructure, meaning the pillar in  Example 5.5: Convergence history curves of the objective function and the 2 -norm of the shape gradient during the optimization runs the first example and the beam in the second example, is meshed with 8 192 triangular elements. The microstructure is computed on refinement level = 5 , which yields 16 384 triangular elements in the case of the hollow material and 24 576 in the composite material. In both cases, the Lagrange multiplier in the functional (2.15) is initially set to 0 and is updated in accordance with (2.16) whenever the 2 -norm of the shape gradient has been halved.

First example: compliance minimization in a partly hollow pillar
We consider a two-dimensional pillar of 10 m height and 5 m width. It is composed of a material with elasticity modulus E = 50 Pa and Poisson ratio 0.2. In this pillar, we aim at finding an optimal design for the shape of the holes in its microstructure, such that the compliance of the pillar is minimized while 10% of the pillar remain hollow.
The parameter in the objective function (2.15) is set to = 2 × 10 4 . We want 10% of the material to be hollow in its microstructure, thus we set 0 = 0.1 . A downward pointing force of 10 N is applied on each point on the top of the pillar, i.e. g(x) = (0, −10) ⊺ for x ∈ [0, 5] × 10.
We start the optimization run with a circle of radius 0.25 as initial guess. The optimization run is halted after 20 steps since the interface elongates so far that it would leave the cell. Hereby, the value of is updated 6 times during the optimization run and the value of the 2 -norm of the shape gradient evolves from approximately 35 to 49. We have J C ( 0 , 0) ≈ 247 at the beginning of the optimization run when evaluating the objective function (2.15) with = 0 , and at the end in step 20 we have J C ( 20 , 0) ≈ 107.
The results of the optimization run are displayed in Fig. 14 while the histories of the objective function and of the 2 -norm of the shape gradient during the compliance minimization process are displayed in Fig. 15. The result corresponds to the intuitive solution of the inclusion getting thin and stretching vertically to reduce the compliance of the porous material whilst retaining a certain volume.

Second example: compliance minimization in an inhomogeneous beam
We consider a two-dimensional beam of 1 m height and 10 m width. The beam is made of a material with elasticity modulus E = 100 Pa and Poisson ratio 0.3 within the inclusions of its microstructure. Outside of the inclusions, the material is inhomogeneous with an elasticity modulus Ẽ (r) ≤ 20 Pa and a Poisson ratio ̃(r) ≥ 0.2 , where for r > 0 material in the point y are given by the evaluation of the functions Ẽ (r) and ̃(r) at r ∶= ‖y − y 0 ‖ 2 ∕‖y M − y 0 ‖ 2 . This setting accounts for a material that is increasingly compressed when getting closer to the inclusion. We thereby followed the description of the evolution of the elasticity modulus and the Poisson ratio of an increasingly compressed material in Alliche (2016). Note that the inhomogeneity in the material around the inclusion can be interpreted to be caused by the manufacturing process. The material around the inclusion is compressed, thus creating space for the material to be added inside the inclusion. In this way, the density of the material increases around the inclusion which results in this inhomogeneous material considered in this example. The beam is fixed at its left edge and a downward force of 0.01 N acts at the upper right corner of the beam, i.e. g(x) = (0, −0.01) ⊺ for x = (10, 1) and g(x) = 0 elsewhere. We then aim at finding an optimal design for the shape of the inclusions in its microstructure, such that the beam's compliance is minimized while the material contained in the inclusion is 25% of the material in the beam, i.e. 0 = 0.25 . To this end, we set the parameter in the objective function (2.15) to  . 15 Convergence history curves of the objective function and the 2 -norm of the shape gradient during the compliance minimization run in a pillar = 2 × 10 −1 . This value differs from the value of in Sect. 6.1 because in both examples has been chosen proportionally to the value of the compliance (2.14) evaluated at the initial shape 0 . The initial compliance is ∼10 2 in Sect. 6.1 while here it is ∼10 −2 .
We start the optimization run with a circle of radius 0.1 as initial guess. The optimization was halted after 16 optimization steps, since the quadratic line search suggested to use a step size smaller than 10 −7 in five consecutive optimization steps. The outcome of the optimization run is found in Fig. 16. During the optimization run, the value of was updated once after the tenth optimization step, which is visible in Fig. 17. The value of the 2 -norm of the shape gradient evolves from approximately 3.19 × 10 −2 to 9.79 × 10 −3 . We have J C ( 0 , 0) ≈ 1.66 × 10 −2 at the beginning of the optimization run when evaluating the objective function (2.15) with = 0 , and at the end in step 16 we have J C ( 16 , 0) ≈ 7.14 × 10 −3 .

Fig. 16
On the right, the initial and optimized inclusions in the microstructure of the composite material are shown. The left-hand side shows three plots of a beam of height 1 m and width 10 m experiencing a downward force of 0.01 N on its right upper corner. Three settings are seen, one with the initial microstructure (top), one with the optimized microstructure (middle), and finally one with the material composed solely of the more elastic material (bottom) Fig. 17 Convergence history curves of the objective function and the 2 -norm of the shape gradient during the compliance minimization run in a beam

Conclusion
In this article, shape sensitivity analysis of the effective tensor in case of periodic materials has been performed. In particular, we computed the Hadamard representation of the shape derivative of the least squares matching of a desired material tensor. Numerical experiments by means of a finite element implementation in the two-dimensional setting have been presented for composite materials and for perforated plates. Likewise to Dambrine and Harbrecht (2020), it has been observed that the computed optimal shapes depend on the initial guess, which means that the solution of the problem under consideration seems to be in general not unique. We emphasize, however, that a volume constraint as in the case of compliance optimization is not necessary when considering the least squares matching-the optimization algorithm stagnates in the (nonunique) minimum.
Funding Open access funding provided by University of Basel.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.