The homogenization method for topology optimization of structures: old and new

TL;DR

This paper reviews the homogenization method for topology optimization, focusing on mathematical tools, applications to lattice structures, and practical exercises using FreeFem++, highlighting both classical and recent developments.

Contribution

It provides a comprehensive overview of homogenization techniques in topology optimization, including new theoretical insights and practical implementation methods.

Findings

Effective homogenization techniques for lattice structure optimization

Numerical methods demonstrated with FreeFem++ software

Enhanced understanding of mathematical tools for topology optimization

Abstract

These are the lecture notes of a short course on the homogenization method for topology optimization of structures, given by Gr\'egoire Allaire, during the "GSIS International Summer School 2018" at Tohoku University (Sendai, Japan). The goal of this course is to review the necessary mathematical tools of homogenization theory and apply them to topology optimization of mechanical structures. The ultimate application, targeted in this course, is the topology optimization of structures built with lattice materials. Practical and numerical exercises are given, based on the finite element free software FreeFem++.

Figures

Equations

$$\{\begin{array}{rlrl}-\mathrm{div}(h\nabla u)& =f& & \text{ in }\Omega ,\\ u& =0& & \text{ on }\partial \Omega ,\end{array}$$

$$0<h_{\min }\le h(x)\le h_{\max }<\infty .$$

$${\mathcal{U}}_{\mathrm{ad}}=\left\{h\in L^{\infty }(\Omega )\text{mathchar}580<h_{\min }\le h(x)\le h_{\max }\text{ a.e. in }\Omega ,{\int }_{\Omega }h(x)dx=h_0|\Omega |\right\},$$

$$J(h)={\int }_{\Omega }j(u)dx,$$

$$j(u)=fu.$$

$$j(u)={\left|u-u_0\right|}^2.$$

$$\underset{h\in {\mathcal{U}}_{\mathrm{ad}}}{\inf }J(h).$$

$$J(h)=\underset{x\in \Omega }{\sup }\left|\sigma (x)\right|$$

$$J(h)={\left({\int }_{\Omega }{\left|\sigma (x)\right|}^{p}dx\right)}^{1/p}.$$

$$\{\begin{array}{rlrl}-\mathrm{div}(h\nabla u)& ={\omega }^{2}u& & \text{ in }\Omega ,\\ u& =0& & \text{ on }\partial \Omega .\end{array}$$

$$\{\begin{array}{rlrl}-\mathrm{div}(h\nabla u_i)& =f_i& & \text{ in }\Omega ,\\ u_i& =0& & \text{ on }\partial \Omega .\end{array}$$

$$J(h)=\sum _{i=1}^n{c}_i{\int }_{\Omega }j(u_i)dx,$$

$$J(h)=\underset{1\le i\le n}{\max }\left\{{\int }_{\Omega }j(u_i)dx\right\}.$$

$$\partial \Omega =\Gamma \cup {\Gamma }_D\cup {\Gamma }_N,$$

$$\{\begin{array}{rlrl}-\Delta u& =0& & \text{ in }\Omega ,\\ u& =0& & \text{ on }{\Gamma }_D,\\ \frac{\partial u}{\partial n}& =g& & \text{ on }{\Gamma }_N,\\ \frac{\partial u}{\partial n}& =0& & \text{ on }\Gamma .\end{array}$$

$${\mathcal{U}}_{\mathrm{ad}}=\left\{\Omega \subset {\mathbb{R}}^N\text{mathchar}58{\Gamma }_D\cup {\Gamma }_N\subset \partial \Omega \text{and}\left|\Omega \right|=V_0\right\},$$

$$\underset{\Omega \in {\mathcal{U}}_{\mathrm{ad}}}{\inf }J(\Omega ),$$

$$J(\Omega )={\int }_{{\Gamma }_N}guds,$$

$$J(\Omega )={\int }_{\Omega }{\left|u-u_0\right|}^{2}dx.$$

$$\{\begin{array}{rlrl}-\mathrm{div}(Ae(u))& =0& & \text{ in }\Omega ,\\ u& =0& & \text{ on }{\Gamma }_D,\\ (Ae(u))\cdot n& =g& & \text{ on }{\Gamma }_N,\\ (Ae(u))\cdot n& =0& & \text{ on }\Gamma ,\end{array}$$

$$\partial \Omega =\Gamma \cup {\Gamma }_D\cup {\Gamma }_N,$$

$${\mathcal{U}}_{\mathrm{ad}}=\left\{\Omega \subset {\mathbb{R}}^N\text{mathchar}58{\Gamma }_D\cup {\Gamma }_N\subset \partial \Omega \text{and}\left|\Omega \right|=V_0\right\},$$

$$J(\Omega )={\int }_{{\Gamma }_N}g\cdot uds,$$

$$J(\Omega )={\int }_{\Omega }{\left|u-u_0\right|}^{2}dx.$$

$$\underset{\Omega \in {\mathcal{U}}_{\mathrm{ad}}}{\inf }J(\Omega ).$$

$$\underset{v\in K}{\inf }J(v).$$

$$u\in K\text{and}\exists \delta >0,\forall v\in K,\Vert v-u\Vert <\delta ⟹J(v)\ge J(u).$$

$$u\in K\text{and}J(v)\ge J(u)\forall v\in K.$$

$$\underset{n\to +\infty }{\lim }J(u^n)=\underset{v\in K}{\inf }J(v).$$

$$\forall (u^n{)}_{n\in \mathbb{N}}\text{sequence in}K,\underset{n\to +\infty }{\lim }\Vert u^n\Vert =+\infty ⟹\underset{n\to +\infty }{\lim }J(u^n)=+\infty .$$

Full text

1. Preface and acknowledgments

2. 1 Introduction

3. 1.1 Optimal design of structures

4. 1.2 Example of sizing or parametric optimization

5. 1.3 Example of shape optimization

6. 1.4 Topology optimization and the homogenization method

7. 1.5 Lattice materials in additive manufacturing

8. 1.6 Goals of these lecture notes

9. 1.7 Exercises

10. 2 Some tools in optimization

11. 2.1 Generalities

12. 2.2 Convex analysis

13. 2.3 Optimality conditions

14. 2.3.1 Minimization with equality constraints

15. 2.3.2 Minimization with inequality constraints

16. 2.3.3 Interpreting the Lagrange multipliers

17. 2.4 Dual energy

18. 2.5 Numerical algorithms

19. 2.5.1 A gradient-type algorithm (non-constrained case)

20. 2.5.2 A gradient-type algorithm (constrained case)

21. 2.6 Exercises

22. 3 Parametric optimal design

23. 3.1 Optimization of a membrane’s thickness

24. 3.2 Existence theories

25. 3.2.1 Definition of a counter-example

26. 3.2.2 Existence for a discretized model

27. 3.2.3 Existence with a regularity constraint

28. 3.3 Computation of a continuous gradient

29. 3.3.1 Adjoint state

30. 3.3.2 A simple formula for the derivative

31. 3.4 A discrete approach

32. 3.4.1 Adjoint state

33. 3.5 Numerical algorithms

34. 3.5.1 Another numerical algorithm for the compliance

35. 3.6 Thickness optimization of an elastic plate

36. 3.6.1 Regularization

37. 3.7 Exercises

38. 4 Homogenization theory

39. 4.1 Homogenization based on two-scale asymptotic expansions

40. 4.2 Composite materials

41. 4.2.1 Lamination for two phase composites

42. 4.2.2 Characterization of $G_{\theta}$

43. 4.3 The elasticity setting

44. 4.4 Numerical applications

45. 4.5 Exercises

46. 5 Topology optimization by the homogenization method

47. 5.1 Why topology optimization?

48. 5.2 Homogenization method in the conductivity setting

49. 5.2.1 Relaxed problem for the conductivity model

50. 5.2.2 Optimality conditions

51. 5.2.3 Numerical algorithm

52. 5.3 Homogenization method in the elasticity setting

53. 5.3.1 Introduction of the model of the elasticity and relaxed problem

54. 5.3.2 An explicit optimal bound $HS(\sigma)$

55. 5.3.3 Optimality conditions

56. 5.3.4 Numerical algorithm

57. 5.4 Convexification, “fictitious materials” and SIMP

58. 5.4.1 Existence of solutions

59. 5.4.2 Optimality condition

60. 5.4.3 Numerical algorithm

61. 5.4.4 Concluding remarks on the SIMP method

62. 5.5 Generalizations of the homogenization method

63. 5.6 Exercises

64. 6 Resurrection of the homogenization method: lattice materials in additive manufacturing

65. 6.1 Introduction

66. 6.2 Setting of the problem

67. 6.3 A three steps approach

68. 6.4 1st step: pre-computing the homogenized properties

69. 6.5 2nd step: parametric optimization of the homogenized problem

70. 6.6 3rd step: reconstruction of an optimal periodic structure

71. 6.6.1 Projection in the simple case without varying orientation, $\alpha\equiv 0$

72. 6.6.2 Projection in the general case with orientation, $\alpha\neq 0$

73. 6.6.3 The new approach of Allaire-Geoffroy-Pantz

74. 6.6.4 Regularization of the double angle $\beta=2\alpha$

75. 6.6.5 Computation of the map $\varphi$

76. 6.6.6 An abstract manifold setting

77. 6.7 A final post-processing/cleaning of the lattice reconstruction

78. 6.8 Other numerical examples

79. 6.9 Further issues

80. 6.10 Exercises

Preface and acknowledgments

These are the lecture notes of a short course on the homogenization method for topology optimization of structures, given by one of us, Grégoire Allaire, during the “GSIS International Summer School 2018” at Tohoku University (Sendai, Japan). Based on the slides of this course, the four other authors, Lorenzo Cavallina, Nobuhito Miyake, Tomoyuki Oka, Toshiaki Yachimura, have written the present lecture notes, which have been proofread by Grégoire Allaire. Each chapter of these lecture notes corresponds to one class, except the two first ones which were taught together.

Topology optimization of structures is nowadays a well developed field with many different approaches and a wealth of applications. One of the earliest method of topology optimization was the homogenization method, introduced in the early eighties. It became extremely popular in its over-simplified version, called SIMP (Solid Isotropic Material with Penalization), which retains only the notion of material density and forgets about true composite materials with optimal (possibly non isotropic) microstructures. However, the appearance of mature additive manufacturing technologies which are able to build finely graded microstructures (sometimes called lattice materials) drastically changed the picture and one can see a resurrection of the homogenization method for such applications. Indeed, homogenization is the right technique to deal with microstructured materials where anisotropy plays a key role, a feature which is absent from SIMP. Homogenization theory allows to replace the microscopic details of the structure (typically a complex networks of bars, trusses and plates) by a simpler effective elasticity tensor describing the mesoscopic properties of the structure.

The goal of this course is to review the necessary mathematical tools of homogenization theory and apply them to topology optimization of mechanical structures. The ultimate application, targeted in this course, is the topology optimization of structures built with lattice materials. Practical and numerical exercises are given, based on the finite element free software FreeFem++.

Finally, the authors would like to express their gratitude to the organizers of the Summer School: Reika Fukuizumi, Kei Funano, Jun Masamune, Jinhae Park, Ruo Li, Shigeru Sakaguchi, Kenjiro Terada, Takayuki Yamada and Lei Zhang. The meeting was partially supported by a grant from the JSPS A3 Foresight Program, JSPS KAKENHI Grant Numbers $26287020$ and $26400062$, GSIS and RCPAM.

Chapter 1 Introduction

1.1 Optimal design of structures

A problem of optimal design (material, shape and topology optimization) of structures is defined by three ingredients (see [Al2007-1, BS2003, HM2003, HP2018, KPTZ2000, SK1992]):

1. (a)

a model (typically a partial differential equation) to evaluate (or analyze) the mechanical behavior of a structure, 2. (b)

an objective function which has to be minimized or maximized, or sometimes several objectives (also called cost functions or criteria), 3. (c)

a set of admissible designs which precisely defines the optimization variables, including possible constraints.

The kind of optimal design problems which we focus on in these lecture notes can be roughly divided into three categories, from the “easiest” to the “most difficult” one:

Parametric or sizing optimization, for which designs are parametrized by a few variables (for example, thickness or member sizes), implying that the set of admissible designs is considerably simplified (see Figure 1-1, where the variable parameters, the thickness of the two boxes in this case, are symbolized by arrows), 2. 2.

Shape (geometric) optimization, for which all designs are obtained from an initial guess by moving its boundary without change of its topology due to the generation of new boundaries (see Figure 1-2, where an admissible shape is drawn with a broken line), 3. 3.

Topology optimization where both the shape and the topology of the admissible designs can vary without any explicit or implicit restrictions (see Figure 1-3, where the broken lines show removable holes).

The last category in the above is, of course, the most general but also the most difficult. We recall that two shapes share the same topology if there exists a continuous deformation from one to the other. In dimension $2$, topology is completely characterized by the number of holes (or, equivalently, of connected components of the boundary). In dimension $3$ it is quite more complicated. Indeed, the topology of a set in dimension $3$ is not only determined by the number of holes, but it also depends on the number and intricacy of “handles” or “loops”.

First of all, one could ask theoretical questions concerning existence, uniqueness, and qualitative properties of the solutions of these shape optimization problems. One could also study the necessary and/or sufficient conditions satisfied by the optimal shapes. Such “optimality conditions” are very important both from a theoretical and a numerical point of view. They are often the basis for numerical algorithms of gradient method type. Furthermore one can investigate the numerical computation of approximate optimal shapes. All these questions will be addressed in the following chapters.

1.2 Example of sizing or parametric optimization

First of all, we show some examples of sizing or parametric optimization. Let us consider the thickness optimization of a membrane, where $\Omega$ is a mean surface of a (plane) membrane and $h$ is the thickness in the normal direction to the mean surface $\Omega$ (see Figure 2).

In what follows, we consider our membrane to be pre-stressed at its boundary and subject to some vertical force $f$. Moreover, for small displacements, small deformations and negligible bending effects in the elasticity, the membrane deformation can be modeled by its vertical displacement $u\mathrel{\mathop{\mathchar 58\relax}}\Omega\to\mathbb{R}$, solution of the following partial differential equation, the so-called membrane model (see also [Al2007-1, K2016]),

[TABLE]

where the thickness $h$ is bounded by some given minimum and maximum values:

[TABLE]

The thickness $h$ is the optimization variable. Notice that we are dealing with a sizing or parametric optimal design problem here, because the computational domain $\Omega$ does not change.

Let us define the set of admissible thickness as follows:

[TABLE]

where $h_{0}$ is an imposed average thickness.

Remark 1.1** (Possible additional “feasibility” constraints).**

According to the production process of membranes, the thickness $h(x)$ can be discontinuous, or on the contrary continuous. A uniform bound can be imposed on its first derivative $h^{\prime}(x)$ (molding-type constraint) or on its second order derivative $h^{\prime\prime}(x)$, linked to the curvature radius (milling-type constraint).

The optimization criterion is linked to some mechanical property of the membrane, evaluated through its displacement $u$, solution of the PDE,

[TABLE]

where, of course, $u$ depends on $h$. For example, the global rigidity of a structure is often measured by its compliance, or work done by the load $f$: the smaller the work, the larger the rigidity (compliance = $-$ rigidity). In such a case, we set

[TABLE]

Another example amounts to achieve (at least approximately) a target displacement $u_{0}(x)$, which is modeled by taking

[TABLE]

Those two criteria are the typical examples studied in this course. Then, a parametric optimization problem is

[TABLE]

Other examples of objective functions are the following:

• •

Introducing the stress vector $\sigma(x)=h(x)\nabla u(x)$, we can minimize the maximum stress norm

[TABLE]

or more generally, for any $p\geq 1$, the following $p$-norm

[TABLE]

• •

For a vibrating structure, introducing the first eigenfrequency $\omega$, defined by

[TABLE]

We consider $J(h)=-\omega$ to maximize it.

• •

Multiple loads optimization: for $n$ given loads $(f_{i})_{1\leq i\leq n}$ the independent displacements $u_{i}$ are solutions of

[TABLE]

We then introduce an aggregated criterion

[TABLE]

with given coefficients $c_{i}$, or

[TABLE]

1.3 Example of shape optimization

In this section, we show two examples of shape optimization. At first let us consider a shape optimization of a membrane’s shape. A reference domain for the membrane is denoted by $\Omega$, with a boundary made of three disjoint parts

[TABLE]

where $\Gamma$ is the variable part, $\Gamma_{D}$ is the Dirichlet (clamped) part and $\Gamma_{N}$ is the Neumann part (loaded by $g$).

The vertical displacement $u$ is the solution of the following membrane model

[TABLE]

From now on the membrane thickness is fixed, equal to $1$. Moreover, we consider the parts $\Gamma_{D}$ and $\Gamma_{N}$ to be given. Thus the set of admissible shapes is

[TABLE]

where $V_{0}>0$ is a given volume. The shape optimization problem reads

[TABLE]

with, as a criterion, the compliance

[TABLE]

or a least-square functional to achieve a target displacement $u_{0}(x)$

[TABLE]

Notice that the true optimization variable is only the free boundary $\Gamma$, and therefore the topology of the shape does not change.

Another example is a shape optimization in the elasticity setting. The model of linearized elasticity gives the displacement vector field $u\mathrel{\mathop{\mathchar 58\relax}}\Omega\to\mathbb{R}^{N}$ as the solution of the system of equations

[TABLE]

with $e(u)=(\nabla u+(\nabla u)^{t})/2$ and $A\xi=2\mu\xi+\lambda(\mathop{\mathrm{tr}}\xi)\mathrm{Id}$, where $\mu$ and $\lambda$ are the Lamé coefficients, and $n$ is the outer unit normal to $\Omega$. The boundary $\partial\Omega$ is again divided into three disjoint parts

[TABLE]

where $\Gamma$ is the free boundary, the true optimization variable. The set of admissible shapes is again

[TABLE]

where $V_{0}$ is a given imposed volume. The objective function chosen is either the compliance

[TABLE]

or a least-square criterion for the target displacement $u_{0}(x)$

[TABLE]

As before, the shape optimization problem reads

[TABLE]

1.4 Topology optimization and the homogenization method

In topology optimization, not only the connected components of the boundary $\Gamma$ are allowed to move but also new connected components (holes in $2$-d) of $\Gamma$ can appear or disappear. Topology is now optimized too. In order to solve this task, we introduce the homogenization method. The homogenization method is a kind of averaging methods for partial differential equations, and is commonly used to determine the averaged (or effective, or homogenized, or equivalent, or macroscopic) parameters of a heterogeneous medium [Al2002, BLP1978, ch1000, CD1999, JKO1995, MT1997, TA2000].

How does homogenization apply to optimal design ? The homogenization method is based on the concept of “relaxation”: it makes ill-posed problems well-posed by enlarging the space of admissible “shapes”. It is crucial to introduce “generalized” shapes, that are “not too generalized”. In the homogenization method, we think of generalized shapes as “limits” of minimizing sequences of classical shapes. We can then say that homogenization allows, as admissible shapes, composite materials obtained by micro-perforation of the original material (fine mixtures of material and void).

1.5 Lattice materials in additive manufacturing

Additive manufacturing, also known as 3D printing, is a process that creates physical structures built layer by layer from a digital design by using metallic powder melted by a laser or an electron beam [GRS2015]. One of the main advantages of additive manufacturing is that, a priori, there are no limitations on the structures that can be built (unfortunately, in practice there are some limitations of manufacturability, like overhangs or the possibility of thermal residual stresses). Moreover, one can even build microstructures or lattice materials.

However, it is impossible to describe all the fine details of a lattice structure in a finite element model for optimization purposes. Therefore, homogenization theory is the right tool for dealing with lattice materials and related optimal design problems. In chapter 6, we will tackle the problem of optimizing lattice structures by using the homogenization method.

1.6 Goals of these lecture notes

The main goal of these lecture notes is to introduce the homogenization method for topology optimization of structures. The rest of these lecture notes are organized as follows. In chapter 2, we show some tools in optimization and describe numerical algorithms for computing optimal designs. In chapter 3, we consider parametric optimization problem and compute gradients of objective functions (by an optimal control approach) for further use in gradient-type algorithms. A representative example of parametric optimization is that of a membrane’s thickness. In chapter 4, we provide a brief survey on homogenization theory. In chapter 5, we apply the homogenization method to topology optimization. In chapter 6, we present a resurrection of the homogenization method for the design of lattice materials in additive manufacturing.

Through these lecture notes, numerical exercizes are proposed with the FreeFem++ code (http://www.freefem.org). FreeFem++ is a free software for solving partial differential equations by the finite element method [He2012]. Moreover, you can find some scripts of FreeFem++ for shape optimization in the web site (http://www.cmap.polytechnique.fr/~allaire/freefem_en.html) and in the corresponding educational paper [AP2006].

1.7 Exercises

Problem 1.7.1**.**

Solve (with FreeFem++) the elasticity equations for the following test cases: cantilever, bridge, MBB beam and L-beam (see Figure 7).

Chapter 2 Some tools in optimization

We review some classical result in optimization theory. More details can be found in textbooks like [Al2007-2, BGLS2006, ET1999, NW1999].

2.1 Generalities

Let $V$ be a Banach space and $K\subset V$ be a non-empty subset. Let $J\mathrel{\mathop{\mathchar 58\relax}}V\to\mathbb{R}$. We consider the following minimization problem

[TABLE]

Let us specify some basic definitions.

Definition 2.1**.**

An element $u$ is called a local minimizer of $J$ on $K$ if

[TABLE]

Moreover, an element $u$ is called a global minimizer of $J$ on $K$ if

[TABLE]

Definition 2.2**.**

A minimizing sequence of a function $J$ on the set $K$ is a sequence $(u^{n})_{n\in\mathbb{N}}\subset K$ such that

[TABLE]

By definition of the infimum value of $J$ on $K$ there always exists at least one minimizing sequence for $J$ on $K$.

Let us consider the existence of minima for optimization problems in finite dimension. The following result guarantees the existence of a minimum.

Theorem 2.3**.**

Let $K$ be a non-empty closed subset of $\mathbb{R}^{N}$ and $J$ a continuous function from $K$ to $\mathbb{R}$ satisfying the so-called “infinite at infinity” property, i.e.,

[TABLE]

Then there exists at least one minimizer of $J$ on $K$. Furthermore, from each minimizing sequence of $J$ over $K$ one can extract a subsequence which converges to a minimum of $J$ on $K$.

Proof.

Let $(u^{n})_{n\in\mathbb{N}}$ be a minimizing sequence for $J$ over $K$. In particular, since $J$ is infinite at infinity and the sequence $\left(J(u^{n})\right)_{n\in\mathbb{N}}$ is bounded, we conclude that $(u^{n})_{n\in\mathbb{N}}$ must be bounded as well. Therefore, since closed bounded sets are compact in finite dimension, there exists a subsequence $\left(u^{n_{k}}\right)_{k\in\mathbb{N}}$ that converges to a point $u\in{{\mathbb{R}}^{N}}$. Now, $u\in K$ because $K$ is closed, and $J(u^{n_{k}})$ converges to $J(u)$ by continuity. We conclude that $\displaystyle J(u)=\lim_{k\to\infty}J(u^{n_{k}})=\inf_{K}J$.

∎

Remark 2.4**.**

In an infinite dimensional vector space, a continuous function on a closed bounded set does not necessarily attain its minimum. For example, let $H^{1}(0,1)$ be the usual Sobolev space with the norm $\displaystyle\mathinner{\!\left\lVert v\right\rVert}=\left(\int^{1}_{0}(v^{\prime}(x)^{2}+v(x)^{2})dx\right)^{\frac{1}{2}}$. Let

[TABLE]

One can check that $J$ is continuous and “infinite at infinity”. Nevertheless the minimization problem

[TABLE]

does not admit a minimizer. Indeed, there exists no $v\in H^{1}(0,1)$ such that $J(v)=0$ but, still,

[TABLE]

To obtain it, we construct a minimizing sequence $(u^{n})_{n\in\mathbb{N}}$ defined for, $n\geq 1$, by

[TABLE]

as Figure 8.

We can easily check that $u^{n}\in H^{1}(0,1)$ and $(u^{n})^{\prime}=\pm 1$. Consequently,

[TABLE]

We clearly see in this example that the minimizing sequence $(u^{n})_{n\in\mathbb{N}}$ is “oscillating” more and more and it is not compact in $H^{1}(0,1)$ despite being bounded in the same space.

2.2 Convex analysis

As we have seen in Remark 2.4, continuous functions do not necessarily attain their minimum on a bounded closed set. In order to extend the result of Theorem 2.3 to the case of an infinite dimensional Hilbert space, we shall work in a convex framework.

Definition 2.5**.**

A set $K\subset V$ is said to be convex if, for any $x,y\in K$ and for any $\theta\in[0,1]$, the linear combination $\theta x+(1-\theta)y$ belongs to $K$.

Definition 2.6**.**

A function $J$, defined from a non-empty convex set $K\subset V$ into $\mathbb{R}$ is convex on $K$ if

[TABLE]

Furthermore, $J$ is said to be strictly convex if the inequality above is strict whenever $u\neq v$ and $\theta\in(0,1)$.

Theorem 2.7**.**

Let $K$ be a non-empty closed convex set in a reflexive Banach space $V$ (i.e. the dual of $V^{\prime}$ is $V$ itself), and $J$ be a convex continuous function on $K$, which is “infinite at infinity”, i.e.,

[TABLE]

Then, there exists a minimizer of $J$ in $K$.

Remark 2.8**.**

Theorem 2.7 remains true if $V$ is just the dual of some separable Banach space. In particular it holds true when $V=L^{p}(\Omega)$ with $1<p\leq\infty$.

The proof will follow along the same lines as that of Theorem 2.3. However, the infinite dimensional case is much more delicate, since it relies on weak convergence and its relations with convexity. We refer to [Al2007-2, Theorem 9.2.7 and Remark 9.2.9] for a complete proof.

Proposition 2.9**.**

Under the hypotheses of Theorem 2.7, suppose that $J$ is strictly convex. Then $J$ has at most one minimizer.

Proof.

Suppose by contradiction that $u_{1}\neq u_{2}$ are two distinct minimizers of $J$ over the closed strictly convex set $K$. If we take $\theta=1/2$ in (2.1), then we get

[TABLE]

which is contradicts the definition of minimum. ∎

Proposition 2.10**.**

If $J$ is convex on the convex set $K$, then any local minimizer of $J$ on $K$ is a global minimizer.

Proof.

Let $u$ be a local minimizer of $J$ on $K$. Thus there exists $\delta>0$ such that $J(v)\geq J(u)$ for any $v\in K\cap B(u,\delta)$. Let $w\in K\setminus B(u,\delta)$ . Our aim is to show that $J(w)\geq J(u)$. Let $\theta\in(0,1]$ be such that $u+\theta(w-u)\in B(u,\delta)$. For example we can take $\theta=\frac{\delta}{\|w-u\|_{K}}$ . Since $u+\theta(w-u)\in B(u,\delta)$, it follows that $J(u)\leq J(u+\theta(w-u))$, and by Jensen’s inequality

[TABLE]

Thus, $J(u)\leq J(w)$ follows. ∎

Convexity is not the only tool to prove existence of minimizers. Another method is, for example, compactness.

2.3 Optimality conditions

In this section, we discuss optimality conditions for objective functions.

Definition 2.11**.**

Let $V$ be a Banach space. A function $J$, defined from a neighborhood of $u\in V$ into $\mathbb{R}$, is said to be differentiable in the sense of Fréchet at $u$ if there exists $L\in V^{\prime}$ such that

[TABLE]

We call $L$ the differential (or derivative, or gradient) of $J$ at $u$ and we denote it by

[TABLE]

Remark 2.12**.**

If $V$ is a Hilbert space, its dual $V^{\prime}$ can be identified with $V$ itself thanks to the Riesz representation theorem. Thus, there exists a unique $p\in V$ such that $\langle p,w\rangle=L(w)$. We also write $p=J^{\prime}(u)$. We use this identification $V=V^{\prime}$ if $V=\mathbb{R}^{n}$ or $V=L^{2}(\Omega)$. In practice, it is often easier to compute the directional derivative $j^{\prime}(0)=\langle J^{\prime}(u),w\rangle_{V^{\prime},V}$ with $j(t)=J(u+tw)$.

Consider the variational formulation

[TABLE]

where $a$ is a symmetric coercive continuous bilinear form and $L$ is a continuous linear form. By the Lax–Milgram theorem we know that the variational formulation (2.3) admits a unique solution. Let us now define the energy

[TABLE]

The following lemma tells us the relationship between the energy $J$ and the variational formulation (2.3).

Lemma 2.13**.**

Let $u\in V$ be the unique solution of the variational formulation (2.3). Then $u$ is the unique minimizer of $J$, that is,

[TABLE]

Conversely, if $u\in V$ is a point giving an energy minimum of $J(v)$, then $u$ is the unique solution of the variational formulation (2.3).

Proof.

If $u$ is the solution of the variational formulation (2.3), then thanks to the symmetry of $a$ we have

[TABLE]

As $u+v$ is arbitrary in $V$, $u$ minimizes the energy $J$ in $V$.

Conversely, let $u\in V$ be such that

[TABLE]

For $v\in V$ we define $j(t)=J(u+tv)$. Then

[TABLE]

We differentiate $t\mapsto j(t)$,

[TABLE]

By definition, $j^{\prime}(0)=\langle J^{\prime}(u),v\rangle_{V^{\prime},V}$, thus

[TABLE]

Since $t=0$ is a minimum point of $j$, we have $a(u,v)=L(v)$ for all $v\in V$. ∎

Remark 2.14**.**

When computing the Fréchet differential of a given functional $J$ at $u$ (see the definition of $L$ and $w\mapsto L(w)$ in (2.2)), there is not always an obvious way to deduce a formula for $J^{\prime}(u)$, nevertheless most of the time it is enough to know the mapping $w\mapsto\langle J^{\prime}(u),w\rangle$.

Example 2.15**.**

For fixed $f\in L^{2}(\Omega)$, define

[TABLE]

We have

[TABLE]

Thus

[TABLE]

Notice that here we identified $L^{2}(\Omega)$ with its dual. 2. 2.

For fixed $f\in L^{2}(\Omega)$ define

[TABLE]

We have

[TABLE]

Therefore, by the usual definition of the duality pairing between $H_{0}^{1}(\Omega)$ and $H^{-1}(\Omega)$ (that comes from a formal integration by parts) we get

[TABLE]

Notice that here the space $H_{0}^{1}(\Omega)$ is not identified with its dual.

Remark 2.16**.**

If instead of the “usual” scalar product in $L^{2}$ we rather use the $H^{1}$ scalar product in the second part of Example 2.15, then we have to identify $J^{\prime}(u)$ with a different function (in other words, the definition of $J^{\prime}(u)$ depends on the scalar product used). From the directional derivative

[TABLE]

using the $H^{1}$ scalar product $\displaystyle\langle\phi,w\rangle=\int_{\Omega}({\nabla}\phi\cdot{\nabla}w+\phi w)\,dx$, we deduce that

[TABLE]

in the distributional sense. Here we identify $H_{0}^{1}(\Omega)$ with its dual.

Theorem 2.17** (Euler inequality).**

Let $K$ be a convex Banach space. Take $u\in K$ and let $J\mathrel{\mathop{\mathchar 58\relax}}K\to\mathbb{R}$ be differentiable at $u$. If $u$ is a local minimizer of $J$ in $K$, then

[TABLE]

On the other hand, if $u\in K$ satisfies this inequality and $J$ is convex, then $u$ is a global minimizer of $J$ in $K$.

Proof.

For $v\in K$ and $\delta\in(0,1]$, we have $u+\delta(v-u)\in K$. Thus, if $\delta$ is sufficiently small, since $u$ is a local minimizer of $J$ in $K$, we have

[TABLE]

We obtain inequality (2.4) by letting $\delta\to 0$ in the above.

We will now prove the second claim of the theorem. Since, by hypothesis $J$ is convex on $K$, then, the graph of $J$ always lies above its tangent plane at any point $w\in K$. In other words, the following inequality holds true for all $v\in K$:

[TABLE]

The conclusion follows by taking $w=u$. ∎

Remark 2.18**.**

If $u$ belongs to the interior of $K$, then we deduce the Euler equation $J^{\prime}(u)=0$.

Remark 2.19**.**

The Euler inequality is usually just a necessary condition (for instance, it is verified also if $u$ is a local maximizer). It becomes a necessary and sufficient condition under the further assumption that the functional $J$ is also convex.

2.3.1 Minimization with equality constraints

We consider the following problem

[TABLE]

where $F=\left(F_{1},\dots,F_{M}\right)$ is a differentiable function from $V$ into $\mathbb{R}^{M}$.

Notice that the set $K=\left\{v\in V\mathrel{\mathop{\mathchar 58\relax}}\;F(v)=0\right\}$ is not necessarily convex. We will therefore need a generalized version of the Euler inequality as stated in Theorem 2.17. To this end we introduce the set of admissible directions for our constrained optimization problem.

Definition 2.20**.**

At every point $v\in K$, the set

[TABLE]

is called the cone of admissible directions at the point $v$.

In other words, $K(v)$ is the set of all vectors that are tangent at $v$ to a curve in $K$ that passes through $v$ (hence, if $K$ is a regular manifold, $K(v)$ coincides with the tangent space to $K$ at $v$). Moreover, notice that, as the name suggests, the set $K(v)$ is a cone in the sense of convex analysis: namely, for all $\lambda\geq 0$ and $w\in K(v)$, then also $\lambda w\in K(v)$.

Proposition 2.21** (Euler inequality, general case).**

Let $u$ be a local minimum of $J$ over $K$. Then, if $J$ is differentiable at $u$, we have

[TABLE]

Proof.

With the same notations of Definition 2.20, set $w^{n}=(v^{n}-v)/\varepsilon^{n}$. By definition, we have that $w\in K(v)$ if and only if there exists a sequence $(w^{n})_{n\in\mathbb{N}}$ in $V$ and a sequence of positive real numbers $(\varepsilon^{n})_{n\in\mathbb{N}}$ such that

[TABLE]

Now, since $u$ is a local minimum of $J$ over $K$, we get

[TABLE]

Passing to the limit as $n\to\infty$ yields

[TABLE]

as claimed. ∎

Definition 2.22**.**

We call Lagrangian of problem (2.5), the function

[TABLE]

The new variable $\mu\in\mathbb{R}^{M}$ is called Lagrange multiplier for the constraint $F(v)=0$.

Lemma 2.23**.**

The constrained minimization problem (2.5) admits the following equivalent formulation using the Lagrangian:

[TABLE]

Proof.

The proof is done by cases. Notice that, if $F(v)=0$, then $J(v)=\mathcal{L}(v,\mu)$ for all $\mu\in\mathbb{R}^{M}$. On the other hand, if $F(v)\neq 0$, then $\displaystyle\sup_{\mu\in\mathbb{R}^{M}}\mathcal{L}(v,\mu)=+\infty$. Putting the two together yields

[TABLE]

∎

Theorem 2.24** (Stationarity of the Lagrangian).**

With the same notation of (2.5), assume that $J$ and $F$ are continuously differentiable in a neighborhood of $u\in V$ such that $F(u)=0$. If $u$ is a local minimizer and if the vectors $\left(F_{i}^{\prime}(u)\right)_{1\leq i\leq M}$ are linearly independent, then there exist a Lagrange multiplier $\lambda\in\mathbb{R}^{M}$ such that

[TABLE]

Proof.

First define $K=\left\{v\in V\mathrel{\mathop{\mathchar 58\relax}}\;F(v)=0\right\}$ and then the corresponding cone of admissible directions $K(u)$ by Definition 2.20. Now, since the vectors $\left(F_{i}^{\prime}(u)\right)_{1\leq i\leq M}$ are linearly independent by hypothesis, we can use the implicit function theorem in a standard way to deduce that

[TABLE]

or equivalently

[TABLE]

In particular $K(u)$ is a vector space (it is indeed the tangent space to the variety $K$ at the point $u$). Thus we can successively take $w$ and $-w$ in Proposition 2.21 to get

[TABLE]

This implies that $J^{\prime}(u)$ is generated by $\left(F_{i}^{\prime}(u)\right)_{1\leq i\leq M}$ (moreover, since the $F_{i}^{\prime}(u)$ are linearly independent, the Lagrange multipliers $\mu_{i}$ are uniquely defined). ∎

2.3.2 Minimization with inequality constraints

We consider the following minimization problem with inequality constraints

[TABLE]

where $F(v)\leq 0$ here means that $F_{i}(v)\leq 0$ for $1\leq i\leq M$, with $F=(F_{1},\dots,F_{M})\mathrel{\mathop{\mathchar 58\relax}}V\to\mathbb{R}^{M}$ differentiable.

Definition 2.25**.**

Let $u$ be such that $F(u)\leq 0$. The set

[TABLE]

is called the set of active constraints at $u$. The inequality constraints are said to be qualified at $u\in K$ if the vectors $\left(F_{i}^{\prime}(u)\right)_{i\in I(u)}$ are linearly independent.

There are other (more general) definitions of constraints qualification [BGLS2006].

Definition 2.26**.**

We call Lagrangian of the previous problem the function

[TABLE]

The new non negative variable $\mu\in\left(\mathbb{R}_{\geq 0}\right)^{M}$ is called Lagrange multiplier for the constraint $F(v)\leq 0$.

The proof of the result below is analogous to that of Lemma 2.23 and thus will be omitted.

Lemma 2.27**.**

The constrained minimization problem (2.7) is equivalent to

[TABLE]

The existence of (non negative) Lagrange multipliers, analogous to Theorem 2.24, can be proved also for a minimization problem subject to inequality constraints. We refer to [Al2007-2, Theorem 10.2.15] for a proof.

Theorem 2.28** (Stationarity of the Lagrangian for the inequality constraint).**

We assume that the constraints are qualified at $u$ satisfying $F(u)\leq 0$. If $u$ is a local minimizer, there exist Lagrange multipliers $\lambda_{1},\dots,\lambda_{M}\geq 0$ such that

[TABLE]

The condition (2.8) is indeed the stationarity of the Lagrangian since

[TABLE]

and the condition $F(u)\leq 0$ and $\lambda\cdot F(u)=0$ for $\lambda\geq 0$, is equivalent to the Euler inequality (Theorem 2.17) associated to the maximization problem $\sup\mathcal{L}(u,\mu)$ with respect to the variable $\mu$ in the closed convex set $\left(\mathbb{R}_{\geq 0}\right)^{M}$. Indeed

[TABLE]

and thus $F(u)\cdot\mu\leq F(u)\cdot\lambda=0$ for all $\mu\in(\mathbb{R}_{\geq 0})^{M}$ as claimed.

2.3.3 Interpreting the Lagrange multipliers

Define the Lagrangian for the minimization of $J(v)$ under the constraint $F(v)=c$ as follows:

[TABLE]

We claim that the value of the Lagrange multiplier represents the sensitivity of the minimal value with respect to variations of the constraint $c$. To this end, let $u(c)$ and $\lambda(c)$ denote the minimizer and the optimal Lagrange multiplier respectively. Moreover, we assume that they are differentiable with respect to $c$. Then

[TABLE]

In other words, $\lambda$ gives the derivative of the minimal value with respect to $c$ without any further calculation. Indeed

[TABLE]

where, in the last equality we used

[TABLE]

which are a consequence of Theorem 2.24 and the constraint $F(u(c))=c$ respectively.

2.4 Dual energy

In this section, we shall associate to a minimizing problem with a maximizing problem, so called dual problem. To simplify the argument, we will assume that $V$ and $Y$ are two Banach spaces. Let $V^{\prime}$ and $Y^{\prime}$ be the corresponding dual spaces. The following argument is according to [ET1999] and see the book for the more general setting. For $J\colon V\to\mathbb{R}\cup\{\infty\}$, we consider the following minimizing problem:

[TABLE]

For given problem (2.9), we are now able to define a dual problem. We shall consider a function $\Phi\colon V\times Y\to\mathbb{R}\cup\{\infty\}$ such that

[TABLE]

We define the conjugate function $\Phi^{*}\colon V^{\prime}\times Y^{\prime}\to\mathbb{R}\cup\{\infty\}$ as

[TABLE]

We call the problem

[TABLE]

the dual problem of (2.9).

In the following, we will mention the relationship between (2.9) and (2.10) in a special case. Let $\Lambda\colon V\to Y$ be a continuous linear operator. Assume that $J$ can be rewritten as

[TABLE]

where $\tilde{J}$ is a function of $V\times Y$ into $\mathbb{R}\cup\{\infty\}$. In this case, the function $\Phi$ will be

[TABLE]

Then the conjugate function $\Phi^{*}$ becomes

[TABLE]

For this case, we can see the following relationship.

Theorem 2.29**.**

Assume that $\tilde{J}$ is convex and (2.9) is finite. We also assume that there exists $v_{0}\in V$ such that $\tilde{J}(v_{0},\Lambda v_{0})<\infty$ and the function $p\mapsto\tilde{J}(v_{0},p)$ is continuous at $\Lambda v_{0}$. Then

[TABLE]

and maximizing problem (2.10) has at least one solution.

To show Theorem 2.29, we will use convex analysis. For the details of the proof, see [ET1999, Chapter III, Theorem 4.1].

Example 2.30**.**

We show an application of Theorem 2.29. Let $\Omega\subset\mathbb{R}^{N}$ be a smooth domain. We consider the Dirichlet problem

[TABLE]

where $f\in L^{2}(\Omega)$ and $h\colon\Omega\to\mathbb{R}$ is a positive given function. The solution $u$ of the problem above is the minimizer of

[TABLE]

We can apply Theorem 2.29 with

[TABLE]

In the case, we see that

[TABLE]

and hence

[TABLE]

2.5 Numerical algorithms

In this section we present some numerical algorithms in order to solve the kind of minimization problems that were treated in this chapter. All these algorithms are of iterative nature: starting from a give initial value $u^{0}$, we construct a sequence $(u^{n})_{n\in\mathbb{N}}$, which can be shown to converge to the solution $u$ of the given minimization problem under some hypotheses.

2.5.1 A gradient-type algorithm (non-constrained case)

Suppose that $V=\mathbb{R}^{N}$ (or, more generally, a Hilbert space, that we will identify with its dual $V^{\prime}$). We consider the following minimization problem without constraints:

[TABLE]

We initialize the algorithm by choosing some initial value $u^{0}\in V$ and iteratively update it as follows:

[TABLE]

where $\mu$ is a positive parameter that we choose in advance (a more sophisticate algorithm involving the optimal choice of $\mu=\mu^{n}$ for each iteration is discussed in [Al2007-1, Theorem 3.38]).

Theorem 2.31**.**

Let $V$ be a Hilbert space and suppose that the functional $J\mathrel{\mathop{\mathchar 58\relax}}V\to\mathbb{R}$ is strongly convex, that is, for some $\alpha>0$

[TABLE]

Moreover, assume that $J$ is differentiable with Lipschitz continuous derivative $J^{\prime}$. Then, if $\mu$ is small enough (depending on $\alpha$ and on the Lipschitz constant of $J^{\prime}$), the gradient-type algorithm described above converges. In other words, for all $u^{0}$, the sequence $(u^{n})_{n\in\mathbb{N}}$ defined in (2.12) converges to the solution $u$ of (2.11).

For a proof, see [Al2007-2].

Remark 2.32**.**

Choosing the right step length is not an easy task. Let us use the line search strategy as follows: start with a given step $\mu^{0}>0$. Now, at each iteration, increase the current step, $\mu_{n+1}=1.1\times\mu_{n}$, if $J$ decreases, and reduce it, $\mu_{n+1}=0.5\times\mu_{n}$ if $J$ increases.

2.5.2 A gradient-type algorithm (constrained case)

Suppose that $J$ is a real valued strictly convex differentiable functional defined on a nonempty closed convex subset $K$ of the Hilbert space $V$. The set $K$ represents the imposed constraints. We consider the following minimization problem

[TABLE]

Theorem 2.3 ensures the existence of a minimizer $u$ for (2.13) (which is unique by Proposition (2.9)). Moreover, according to Theorem 2.17, the minimizer $u$ is characterized by the condition

[TABLE]

Notice that the condition above can be rephrased as follows. For all $\mu>0$

[TABLE]

Let $P_{K}\mathrel{\mathop{\mathchar 58\relax}}V\to K$ denote the projection operator onto the convex subset $K$. Then (2.14) just states that $u$ is the orthogonal projection of $u-\mu J^{\prime}(u)$ onto $K$. In other words

[TABLE]

Therefore we devise a (projected) gradient-type algorithm, defined by the following iteration

[TABLE]

where $\mu$ is a fixed positive parameter.

Theorem 2.33**.**

Let $J$ be a differentiable strongly convex functional, with derivative $J^{\prime}$ Lipschitz continuous on $V$. Then, if $\mu$ is small enough, the projected gradient algorithm with fixed step defined above converges. In other words, for all initial values $u^{0}\in K$, the sequence $(u^{n})_{n\in\mathbb{N}}$ defined by (2.15) converges to the solution $u$ of (2.13).

We refer to [Al2007-2, Theorem 10.5.8] for a proof.

Remark 2.34**.**

Another possibility is to penalize the constraints, i.e., for small $\varepsilon$ we replace the problem

[TABLE]

Example 2.35** (Some projection operators $P_{K}$).**

Here we present some projection operators that can be computed explicitly.

• •

If $V=\mathbb{R}^{M}$ and $\displaystyle K=\prod_{i=1}^{M}\,[a_{i},b_{i}]$, then for $x=(x_{1},\dots,x_{M})\in\mathbb{R}^{M}$ we have

[TABLE]

• •

If $V=\mathbb{R}^{M}$ and \displaystyle K=\Big{\{}x\in\mathbb{R}^{M}\mathrel{\mathop{\mathchar 58\relax}}\;\sum_{i=1}^{M}x_{i}=c_{0}\Big{\}}, then

[TABLE]

• •

Similarly, if $V=L^{2}(\Omega)$ and $K=\{\phi\in V\mathrel{\mathop{\mathchar 58\relax}}\;a(x)\leq\phi(x)\leq b(x)\}$ or $\displaystyle K=\left\{\phi\in V\mathrel{\mathop{\mathchar 58\relax}}\;\int_{\Omega}\phi\,dx=c_{0}\right\}$ the corresponding projection operators $P_{K}$ can be obtained by replacing finite sums with integrals in the two examples above.

For more general closed convex sets $K$, the corresponding projection operator $P_{K}$ can be very hard to determine. In such cases one can use the so called Uzawa algorithm [Al2007-2] which looks for a saddle point of the Lagrangian.

2.6 Exercises

Problem 2.6.1**.**

For a given $f\in L^{2}(\Omega)$, $\Omega$ being a rectangle in 2D, solve the following optimization problem numerically under the constraints $0\leq u(x)\leq 1$ and $\displaystyle\int_{\Omega}u\,dx=|\Omega|/2$:**

[TABLE]

Problem 2.6.2**.**

For a given $f\in L^{2}(\Omega)$, $\Omega$ being a rectangle in 2D, and $\varepsilon>0$, solve the following optimization problem numerically under the constraints $0\leq u(x)\leq 1$:**

[TABLE]

Chapter 3 Parametric optimal design

3.1 Optimization of a membrane’s thickness

In this section, we consider a parametric optimal design problem of a membrane. Let $\Omega$ be a bounded domain in ${{\mathbb{R}}^{N}}$ ($N\geq 2$) and $f\in L^{2}(\Omega)$ be external forces. Let us consider the displacement $u\in H^{1}_{0}(\Omega)$, defined as the solution of

[TABLE]

where $h=h(x)$ is the thickness of membrane. Note that Lax–Milgram theorem ensures that there exists a unique solution $u\in H^{1}_{0}(\Omega)$ of (3.1) if $f\in L^{2}(\Omega)$. For some given constants $0<h_{\rm min}<h_{0}<h_{\rm max}$, we seek to optimize the membrane by varying its thickness $h(x)$ in the admissible set defined by

[TABLE]

We consider the following parametric shape optimization problem:

[TABLE]

where $u$ depends on $h$ through the state equation (3.1), and $j$ is a $C^{1}$ function from $\mathbb{R}$ to $\mathbb{R}$ such that $|j(u)|\leq C(u^{2}+1)$ and $|j^{\prime}(u)|\leq C(|u|+1)$. As examples of function $j$, we can take $j(u)=fu$ if we want to minimize the compliance (maximize the rigidity of the membrane), or $j(u)=|u-u_{0}|^{2}$ if we want to minimize the least-square criterion to reach a target displacement $u_{0}\in L^{2}(\Omega)$.

Before studying the existence of an optimal thickness, we show the continuity of the cost function.

Proposition 3.1**.**

The application

[TABLE]

is a continuous mapping from $\mathcal{U}_{\rm ad}$ into $\mathbb{R}$.

Proof.

The result follows immediately by composition of the two continuous functions that appear in the following lemmas: Lemma 3.2 and Lemma 3.3. ∎

Lemma 3.2**.**

The map $\displaystyle v\mapsto\int_{\Omega}j(v)\,dx$ is continuous from $L^{2}(\Omega)$ into $\mathbb{R}$.

Proof.

The result follows by the Lebesgue dominated convergence theorem. ∎

Lemma 3.3**.**

The map $h\mapsto u$, where $u\in H_{0}^{1}(\Omega)$ is the solution of (3.1), is a continuous function from $\mathcal{U}_{\rm ab}$ into $H_{0}^{1}(\Omega)$.

Proof.

Let $(h_{n})_{n\in\mathbb{N}}\subset\mathcal{U}_{\rm ab}$ be a sequence converging in the $L^{\infty}$-norm to some $h_{\infty}\in L^{\infty}(\Omega)$. Let $u_{n}\in H_{0}^{1}(\Omega)$ denote the unique solution of the membrane equation with associated thickness $h_{n}$:

[TABLE]

or the equivalent weak formulation

[TABLE]

We will prove that $u_{n}$ is a Cauchy sequence in $H_{0}^{1}(\Omega)$ and thus it converges. Take $n,m\in\mathbb{N}$ and subtract the variational formulation for $u_{n}$ (3.2) from that of $u_{m}$ for fixed $\phi\in H_{0}^{1}(\Omega)$ to be chosen later. We get

[TABLE]

Choosing $\phi=u_{m}-u_{n}$ we deduce

[TABLE]

which proves the claim. ∎

3.2 Existence theories

The question of the existence of optimal shapes is far from simple. We cannot apply the results of Chapter 2 directly since $J(h)$ is not generally convex function. In fact, there exists no optimal shape in general. General counter-examples have been found by Murat [Mu1977]. It is an important issue because this non-existence phenomenon has dramatic consequences for the numerical computations. Thus the definition of the set $\mathcal{U}_{\rm ab}$ of admissible designs has to be modified in order to obtain existence of optimal shapes. The main strategies employed to gain the existence of optimal shapes are discretization (when the admissible set is made finite dimensional), regularization (when the admissible set is made compact), and sometimes a miracle (when the given optimization problem happens to be convex).

3.2.1 Definition of a counter-example

First, let us show a counter-example to the existence of optimal design for the membrane problem. For simplicity, let $N=2$ and $\Omega=(0,1)\times(0,1)$.

We want to minimize the following objective function for $h\in\mathcal{U}_{\rm ad}$:

[TABLE]

where $e_{1}$, $e_{2}$ are the horizontal and vertical directions $(1,0)$, $(0,1)$ respectively and $u_{1}$, $u_{2}$ are the solutions of the following membrane problems:

[TABLE]

When we minimize (3.3), we want the membrane to be strong for horizontal loading (we minimize compliance in the $e_{1}$ direction), and at the same time weak for vertical loading (we maximize the compliance in the direction $e_{2}$). This property of the objective function makes the problem ill-posed in the following sense.

Theorem 3.4**.**

The infimum of (3.3) is not attained by any $h\in\mathcal{U}_{\rm ad}$.

Since the rigorous proof of Theorem 3.4 is a little bit technical, here we will only explain the main ideas by means of a “hand-waving argument”. First of all, notice that if $h$ is uniform (i.e. $h$ is a constant function), then by definition the membrane is isotropic. Therefore, also the domain $\Omega$ is isotropic, that is to say that it shows the same mechanical behavior in all direction. However, it is better to build horizontal layers of alternating small and large thicknesses in order to minimize the objective function (3.3) (see Figure 10). In other words, we are building a laminated structure that is horizontally strong but vertically weak. In order to intuitively justify this statement consider the following. Vertically, the lines of forces must cross the layers of minimal thickness: this means that the structure is thus weak with respect to vertical stress. On the other hand, horizontally, the lines of forces follow the layers of maximal thickness: this means that the structure is thus strong with respect to horizontal stress. However, since the boundary conditions are uniform, the membrane is horizontally stronger if the layers are finer, as the lines of forces are deviating from the horizontal to a lesser extent. If $h$ oscillates at a small scale, we obtain an anisotropic composite material. To reach the minimum, the oscillation scale must go to [math]. Therefore, there does not exist any real optimal design that does not involve a microstructure at an infinitely small scale. We refer the interested reader to Section 5.2 in [Al2007-1] for the details.

3.2.2 Existence for a discretized model

One way to avoid non-existence due to a loss of compactness consists in working with a discretized (and hence finite-dimensional) model. Let $(\omega_{i})_{1\leq i\leq n}$ be a partition of $\Omega$ such that

[TABLE]

We introduce the subset $\mathcal{U}_{\rm ad}^{n}$ of $\mathcal{U}_{\rm ab}$ defined by

[TABLE]

In other words, any function $h\in\mathcal{U}_{\rm ab}^{n}$ is uniquely determined by the choice of the vector $(h_{i})_{1\leq i\leq n}\in\mathbb{R}^{n}$ and thus $\mathcal{U}_{\rm ad}^{n}$ is identified with a closed subset of $\mathbb{R}^{n}$.

Theorem 3.5** (Existence in finite dimension).**

The discretized optimization problem

[TABLE]

admits at least one minimizer.

Proof.

Since $\mathcal{U}_{\rm ab}^{n}$ is a compact subset of ${{\mathbb{R}}^{N}}$ and $J(h)$ is a continuous function on $\mathcal{U}_{\rm ab}^{n}$, the existence of a minimizer of $J$ in $\mathcal{U}_{\rm ab}^{n}$ follows from Theorem 2.3. ∎

3.2.3 Existence with a regularity constraint

Another classical way of ensuring the existence of minimizers relies in imposing additional regularity. For example, consider the space $C^{1}(\overline{\Omega})$ which is a Banach space with the norm

[TABLE]

Take a given constant $R>0$ and introduce the subspace $\mathcal{U}_{\rm ad}^{\rm reg}$:

[TABLE]

The upper bound on the $C^{1}$-norm of $h$ in the definition above can be interpreted as a “feasibility” (or “manufacturability”) constraint, as, in practice, the thickness cannot vary too rapidly. Then the following theorem holds:

Theorem 3.6**.**

The regularized optimization problem

[TABLE]

admits at least one minimizer.

Proof.

Consider a minimizing sequence $(h_{n})_{n\in\mathbb{N}}\subset\mathcal{U}_{\rm ad}^{\rm reg}$ such that

[TABLE]

By definition, the sequence $(h_{n})_{n\in\mathbb{N}}$ is bounded uniformly in $n$ in the space $C^{1}(\overline{\Omega})$. We then apply a variant of Rellich theorem which states that one can extract a subsequence (still denoted by $h_{n}$ for simplicity) that converges in $C^{0}(\overline{\Omega})$ to a limit function $h_{\infty}$ (furthermore, we know that $h_{\infty}\in C^{1}(\overline{\Omega})$). We already know that $h\mapsto J(h)$ is a continuous mapping from $\mathcal{U}_{\rm ad}$ into $\mathbb{R}$ by Proposition 3.1, therefore

[TABLE]

which proves that $h_{\infty}$ is a global minimizer of $J$ in $\mathcal{U}_{\rm ad}^{\rm reg}$ as claimed. ∎

Remark 3.7**.**

Theorem 3.6 is actually a theorem of limited practical interest for the following reasons.

• •

In the practical cases, it is not clear how to choose the upper bound $R$ in the definition of $\mathcal{U}_{\rm ad}^{\rm reg}$.

• •

Usually we do not have convergence as $R$ goes to infinity.

• •

It is not clear whether, numerically, we have global or local minimizers.

• •

Numerically, an upper bound on the $H^{1}$-norm is preferred instead:

[TABLE]

3.3 Computation of a continuous gradient

In this section, we will calculate the gradient of the objective function $J(h)$. This tells us the necessary conditions for optimality of the optimal shape and allows us to establish a numerical algorithm for calculating the optimal shape.

First, we consider the boundary value problem

[TABLE]

where $h$ belong to the following convex set which is larger than \mathcal{U}_{\rm ad}$$\colon

[TABLE]

Lemma 3.8**.**

The application $h\mapsto u(h)$, which gives the solution $u(h)\in H^{1}_{0}(\Omega)$ of (3.4) for $h\in\mathcal{U}$, is differentiable and its directional derivative at $h$ in the direction $k\in L^{\infty}(\Omega)$ is given by

[TABLE]

where $v$ is the unique solution in $H^{1}_{0}(\Omega)$ of

[TABLE]

Proof.

Formally, one simply differentiates equation (3.4) with respect to $h$. However, to be mathematically rigorous one should rather work at the level of the variational formulation. To compute the directional derivative with respect to $k\in L^{\infty}(\Omega)$, we define $h(t)=h+tk$ for $t>0$. For $t>0$, let $u(t)$ be the solution for the thickness $h(t)$. Differentiating with respect to $t$ leads to

[TABLE]

and, since $h^{\prime}(0)=k$, we deduce $u^{\prime}(0)=v$.

Let us justify the above calculation by showing that the map $h\mapsto u(h)$ is differentiable in the sense of Fréchet. First, there exists a unique solution $v$ of (3.5) in $H^{1}_{0}(\Omega)$ thanks to the Lax–Milgram Theorem applied to the variational formulation

[TABLE]

We combine (3.6) with the following variational formulation for $u(t)$

[TABLE]

Since $u(1)=u(h+k)$ and $u(0)=u(h)$, we obtain by difference

[TABLE]

Taking $\phi=u(h+k)-u(h)-v$ as a test function in the above yields

[TABLE]

which implies

[TABLE]

where we used Cauchy–Schwarz’s inequality and the $H^{1}_{0}$ boundedness of $v$. Furthermore, by (3.7) we have

[TABLE]

Taking the test function as $\phi=u(h+k)-u(h)$ in (3.10), we obtain the following estimate:

[TABLE]

Combining (3.8) with (3.11), we have

[TABLE]

Therefore we obtain $u(h+k)=u(h)+v+o(k)$ as $\mathinner{\!\left\lVert k\right\rVert}_{L^{\infty}(\Omega)}\to 0$, which proves the claim. ∎

Lemma 3.9**.**

For $h\in\mathcal{U}$, let $u(h)\in H^{1}_{0}(\Omega)$ be the solution to (3.4) and

[TABLE]

where $j$ is a $C^{1}$ function from $\mathbb{R}$ into $\mathbb{R}$ such that $|j(u)|\leq C(u^{2}+1)$ and $|j^{\prime}(u)|\leq C(|u|+1)$ for any $u\in\mathbb{R}$. The application $J(h)$, from $\mathcal{U}$ into $\mathbb{R}$, is differentiable and its directional derivative at $h$ in the direction $k\in L^{\infty}(\Omega)$ is given by

[TABLE]

where $v=\langle u^{\prime}(h),k\rangle$ is the unique solution in $H^{1}_{0}(\Omega)$ of

[TABLE]

Proof.

By simple composition of differentiable applications. To justify it, one only has to check that all the terms are well defined. We omit the details of the proof. ∎

3.3.1 Adjoint state

In order to treat the derivative of the objective function $J(h)$, we introduce the adjoint state $p$, defined as the unique solution in $H^{1}_{0}(\Omega)$ of

[TABLE]

Theorem 3.10**.**

The cost function $J(h)$ is differentiable on $\mathcal{U}$ and

[TABLE]

If $h\in\mathcal{U}_{\rm ad}$ is a local minimizer of $J$ in $\mathcal{U}_{\rm ad}$, then it satisfies the necessary optimality condition

[TABLE]

for any $k\in\mathcal{U}_{\rm ad}$.

Proof.

To make explicit $J^{\prime}(h)$ from Lemma 3.9, we must eliminate $v=\langle u^{\prime}(h),k\rangle$. To this end, we employ the use of the adjoint state, solution of (3.12). multiplying the equation for $v$ by $p$ and that for $p$ by $v$, we integrate by parts

[TABLE]

[TABLE]

Comparing these two equalities we deduce

[TABLE]

for any $k\in L^{\infty}(\Omega)$. Since ${\nabla}u\cdot{\nabla}p$ belongs to $L^{1}(\Omega)$, we check that $J^{\prime}(h)$ is continuous on $L^{\infty}(\Omega)$. To obtain the condition of optimality, it suffices to apply Theorem 2.17 since $\mathcal{U}_{\rm ad}$ is a closed non-empty convex subset of $L^{\infty}(\Omega)$. ∎

Remark 3.11** (How to find the adjoint state).**

For independent variable $(\hat{h},\hat{u},\hat{p})\in L^{\infty}(\Omega)\times H^{1}_{0}(\Omega)\times H^{1}_{0}(\Omega)$, we introduce the Lagrangian

[TABLE]

where $\hat{p}$ is a Lagrange multiplier (a function) for the constraint which connects $u$ to $h$. By integration by parts we get

[TABLE]

The partial derivative of $\mathcal{L}$ at $\hat{u}=u$ in the direction $\phi\in H^{1}_{0}(\Omega)$ is given by

[TABLE]

Notice that, requiring that $\left\langle\dfrac{\partial\mathcal{L}}{\partial\hat{u}}({h},{u},{p}),\phi\right\rangle=0$ for all directions $\phi$ is nothing else than the variational formulation of the adjoint equation (3.12).

3.3.2 A simple formula for the derivative

It is possible to compute the derivative of $J$ by means of the Lagrangian in the following way:

[TABLE]

where $u$ is the state function (solution to (3.4)) and $p$ is the adjoint state (solution to problem (3.12)). Indeed, we have

[TABLE]

by definition of the state function $u$. Thus, if the map $h\mapsto u(h)$ is differentiable, we get for $k\in L^{\infty}(\Omega)$

[TABLE]

Then, taking $\hat{p}=p$, the adjoint we obtain

[TABLE]

By the above discussion, we obtain the following theorem.

Theorem 3.12**.**

Let $\mathcal{L}(\hat{h},\hat{u},\hat{p})$ be the Lagrangian defined as the sum of the objective function and the variational formulation of the state equation, i.e.,

[TABLE]

Let $p$ be the solution of the adjoint equation

[TABLE]

Assume that the solution $u=u(h)$ of the state equation (3.4) is differentiable with respect to $h$. Then the objective function $J$ is differentiable and

[TABLE]

This theorem is the practical method for computing $J^{\prime}(h)$. Once the gradient of the cost function has been obtained, it is natural and quite easy to implement a gradient method to minimize $J(h)$ numerically. In Section 3.5, we provide numerical algorithms to compute the optimal thickness.

3.4 A discrete approach

One can wonder whether the such optimal design problems get simpler after discretization. Unfortunately, the answer is “no”. In this section, we consider a discrete approach to the problems. Applying a finite element method, the equation becomes a linear system of order $n$

[TABLE]

where $K(h)$ is the rigidity matrix of the membrane (which depends on $h$), $b$ is a vector representing the forces $f$, and $y(h)$ the vector of the coordinates of the solution $u$ in the finite element basis (of dimension $n$). We also discretize the admissible set as follows:

[TABLE]

where the finite sum

[TABLE]

is an approximation of

[TABLE]

Approximating the cost function, the discrete problem becomes

[TABLE]

where $j^{\rm disc}$ is a smooth approximation of $j$ from $\mathbb{R}^{N}$ into $\mathbb{R}$. In the case of the compliance we have:

[TABLE]

In the case of a least-square criterion for a target displacement we have:

[TABLE]

where $B$ is a mass matrix. In practice, we need a way to compute the gradient of $J^{\rm disc}(h)$. This can be applied to both finding the optimality condition and the implementation of a numerical method of minimization.

First, we consider the following “naive idea”. Since $y(h)=K(h)^{-1}b$, we have

[TABLE]

where we used the notation $f^{\prime}(h)=\left(\partial f(h)/\partial h_{i}\right)_{1\leq i\leq n}$ and the second identity in (3.13) is a direct application of the formula for the derivative of a matrix. We remark that this method is not practically useful because one must solve $n+1$ linear systems with respect to the matrix $K(h)$ in order to obtain all components of $y^{\prime}(h)$. Recall that $K(h)$ is a very large matrix (of size $n\times n$) and its inverse is never explicitly computed as it would take too long. As a consequence, we do not use the explicit formula $y(h)=K(h)^{-1}b$. We rather use an adjoint method.

3.4.1 Adjoint state

Definition 3.13**.**

We define the adjoint state $p\in\mathbb{R}^{N}$ as the solution of

[TABLE]

By rearranging the second equality of (3.13) we get

[TABLE]

Now, taking the scalar product of (3.15) with $p(h)$ and that of (3.14) with $y^{\prime}(h)$, we obtain, for each component $i=1,\dots,n$:

[TABLE]

from which we deduce

[TABLE]

In practice, this is the very formula that we use for evaluating the gradient $(J^{\rm disc})^{\prime}(h)$ since it requires only to solve two linear systems.

There is no simplification in using a discrete approach rather than a continuous one. Some authors prefer to discretize first and optimize afterwards. This approach guarantees a perfect compatibility between the gradient and the cost function, but it requires a deep knowledge of the numerical solver. Here, we follow another philosophy, “first optimize in a continuous framework, then discretize”. It is much simpler, and no precision is lost if the finite element spaces are adequately chosen.

3.5 Numerical algorithms

In this section, we show numerical algorithms to seek the optimal thickness of $h$. First, we consider the following projected gradient algorithm.

To make the algorithm fully explicit, we have to specify how to compute the projection operator $P_{\mathcal{U}_{\rm ad}}$.

We define the projection operator $P_{\mathcal{U}_{\rm ad}}$ as follows:

[TABLE]

where $\ell$ is the unique Lagrange multiplier such that

[TABLE]

The determination of the constant $\ell$ is not explicit but based on an iterative algorithm. First, notice that the function

[TABLE]

is strictly increasing on the interval $[\ell^{-},\ell^{+}]$, the inverse image of the closed interval $[h_{\min}|\Omega|,h_{\max}|\Omega|]$. Thanks to this monotonicity property, we propose a simple iterative algorithm$\colon$ we first bracket the root by an interval $[\ell^{1},\ell^{2}]$ such that

[TABLE]

then we proceed by dichotomy to find the root $\ell$.

Remark 3.14**.**

**

In practice, we rather use a projected gradient algorithm with a variable step (not optimal) which guarantees the decrease of the functional $J(h_{n+1})<J(h_{n})$.

• 2.

The algorithm is rather slow. A possible acceleration is based on the quasi-Newton algorithm.

• 3.

The overhead generated by the adjoint computation is very modest$\colon$one has to build a new right-hand-side (using the state) and solve the corresponding linear system (with the same rigidity matrix).

• 4.

Convergence is detected when the optimality condition is satisfied with a threshold $\varepsilon>0$

[TABLE]

3.5.1 Another numerical algorithm for the compliance

When $j(u)=fu$, we find $p=-u$ since $j^{\prime}(u)=f$. This particular case is said to be self-adjoint. We use the dual or complementary energy (see Section 2.4)

[TABLE]

in order to rewrite the original optimization problem as a double minimization problem:

[TABLE]

and the order of minimization is irrelevant. This problem is convex and therefore it admits a minimizer.

By elementary calculation, we can show that the following lemma holds.

Lemma 3.15**.**

The function $\phi(a,\sigma)=a^{-1}|\sigma|^{2}$, defined from $\mathbb{R}_{\geq 0}\times\mathbb{R}^{N}$ into $\mathbb{R}$, satisfies

[TABLE]

where the derivative is given by

[TABLE]

In particular, since by (3.16), the graph of $\phi(a,\sigma)$ lies above its linear approximation at each point $(a_{0},\sigma_{0})$, then $\phi$ is convex.

As a result, we obtain the following.

Lemma 3.16** (Optimality conditions).**

For a given $\tau\in L^{2}(\Omega)^{N}$, the problem

[TABLE]

admits a minimizer $h(\tau)$ in $\mathcal{U}_{\rm ad}$ given by

[TABLE]

where $\ell$ is the Lagrange multiplier such that

[TABLE]

Sketch of the proof.

By Lemma 3.15 we obtain that the map $h\mapsto\int_{\Omega}h^{-1}|\tau|^{2}dx$ is convex in $\mathcal{U}_{\rm ad}$. Therefore, Theorem 2.7 ensures the existence of a minimum point $h$. This point is then characterized by the optimality condition given by Theorem 2.17. We refer to [Al2007-1, Lemma 5.2.25] for more details. ∎

Lemma 3.16 tells us the following numerical algorithm for the compliance:

Remark that, by the dual energy approach introduced in Section 2.4, minimizing (3.18) in $\tau$ is equivalent to solving the equation

[TABLE]

and then recovering $\tau_{n}$ by the formula

[TABLE]

The algorithm can be interpreted as an alternate minimization in $\tau$ and $h$ of the objective function. In particular, we deduce that the objective function always decreases through the iterations. Indeed, for all $n\geq 0$,

[TABLE]

where, for fixed $h_{n+1}$ we minimized in $\tau$ and then, for fixed $\tau_{n}$ we minimized in $h$. This algorithm can also be interpreted as an optimality criteria method.

3.6 Thickness optimization of an elastic plate

We consider the following elasticity problem for an elastic plate $\Omega$

[TABLE]

with strain tensor $e(u)=({\nabla}u+({\nabla}u)^{t})/2$. The set of admissible thicknesses is

[TABLE]

The compliance optimization reads

[TABLE]

The theoretical results are the same of previous sections. We apply the optimality criteria method for the compliance optimization (3.19). In order to compute (3.19), we use FreeFem++. You can see its scripts on the web page http://www.cmap.polytechnique.fr/~allaire/freefem_en.html.

In Figure 12, we used finite elements $P2$ for $u$ and $P0$ for $h$. However, numerical instabilities like checkerboards occur if we use finite elements $P1$ for $u$ and $P0$ for $h$ (see Figure 14). Therefore we consider a “regularization” in order to avoid the instabilities.

3.6.1 Regularization

In what follows, let us consider the “regularized” framework to avoid numerical instabilities. The main idea, similar to that introduced in Remark 2.16 is as follows.

We are going to replace the scalar product

[TABLE]

with a different one. Previously we identified $\mathcal{U}_{\rm ad}$ with a subspace of $L^{2}(\Omega)$, thus

[TABLE]

Now, we identify a “regularized” admissible set $\mathcal{U}^{\rm reg}_{\rm ad}$ to a subspace $H^{1}(\Omega)$, thus

[TABLE]

where $\varepsilon>0$ is a regularization parameter (which can be interpreted as a length scale). Therefore, we deduce a new formula for the gradient

[TABLE]

Solving (3.20) and using a gradient algorithm such as projected gradient method, we obtain regularized optimal shape (see Figure 14).

3.7 Exercises

Problem 3.7.1**.**

Check the numerical instabilities (Figure 14, left) by using FreeFem++.

Problem 3.7.2**.**

Solve (3.20) and see the regularized optimal shape (Figure 14, right) by using FreeFem++.

Chapter 4 Homogenization theory

In this section, we explain the homogenization method in order to apply shape optimization problems in Section 5. Homogenization method is one of the averaging methods for partial differential equations. It is often concerned with the derivation of (macroscopic) equations whose solutions are defined as limits of solutions to (microscopic) equations with rapidly varying coefficients. A particular case of homogenization is obtained when the coefficients of the partial differential equation are periodically and rapidly oscillating. Indeed, in many fields of science and technology one has to solve boundary value problems in periodic media. In such a case, homogenization is simpler and can be achieved, at least formally, by using asymptotic expansions. This chapter is devoted to an elementary introduction of periodic homogenization, without providing a fully rigorous justification. Of course, homogenization methods using functional analysis method were considered for mathematical justification. The interested reader is referred to the classical books [BLP1978], [CD1999], [H1996], [JKO1995], for further details.

Note that, for applications in shape optimization, one should rely, in full rigor, on a more general homogenization method, called H-convergence, introduced in [MT1997] (see the textbook [Al2002] for more details). For simplicity, we restrict ourselves to the setting of periodic homogenization which is enough for a formal understanding.

4.1 Homogenization based on two-scale asymptotic expansions

In what follows, we consider an elastic membrane made of a composite material with a fine periodic structure and apply the periodic homogenization method. We assume that ratio between the period and the characteristic size of the structure equals to $\varepsilon\ll 1$. We will find the “true” problem by the limit problem obtained as $\varepsilon\to 0$.

Let $\Omega$ be a bounded domain in ${{\mathbb{R}}^{N}}$ ($N\geq 1$) and $f=f(x)$ be a load. Then we consider the displacement $u_{\varepsilon}$, which is defined as the solution of the following boundary value problem:

[TABLE]

where the coefficient $A(y)$ satisfies the variable Hooke’s law, that is, $A(y)$ is a Y-periodic function with $Y=(0,1)^{N}$. Thus for any i-th vector of the canonical basis $e_{i}$, the coefficient $A(y)$ satisfies

[TABLE]

If we replace $y$ by $x/\varepsilon$, then we obtain that the map $x\mapsto A(x/\varepsilon)$ is a periodic of period $\varepsilon$ in all the coordinate directions $e_{1},\cdots,e_{N}$. A direct computation of $u_{\varepsilon}$ can be very expensive (since the mesh size $h$ should satisfy $h\ll\varepsilon$), thus we seek only the averaged values of $u_{\varepsilon}$. We assume that the solution $u_{\varepsilon}$ can be expanded as follows:

[TABLE]

with $u_{i}(x,y)$ function of the two variables $x$ and $y$, periodic in $y$, with periodicity cell given by $Y=(0,1)^{N}$. Plugging the series (4.2) in the equation (4.1), we use the derivation rule

[TABLE]

Then we get

[TABLE]

Substituting (4.4) into (4.2), the equation becomes a series in $\varepsilon$

[TABLE]

In order to find the solution of the limit equation as $\varepsilon\to 0$, we identify each power of $\varepsilon$. The most important terms are only the first three terms of the series. We start by a technical lemma:

Lemma 4.1**.**

Let $g\in L^{2}(Y)$ and suppose that $A(y)$ is a $Y$-periodic $N\times N$ matrices satisfying

[TABLE]

for some $\lambda>0$. Moreover, let $H^{1}_{\#}(Y)/\mathbb{R}$ denote the quotient space, $H^{1}_{\#}(Y)$ up to an additive constant, equipped with the norm $\|\nabla\cdot\|_{L^{2}(Y)}$. Then the problem

[TABLE]

admits a unique solution $v\in H^{1}_{\#}(Y)/\mathbb{R}$ if and only if

[TABLE]

Proof.

Let us check that $g$ being of zero mean over $Y$ is a necessary condition for existence. As a matter of fact, integrating the equation over $Y$, we get

[TABLE]

because of the periodic boundary condition. Indeed $A(y)\nabla_{y}v(y)$ is periodic, but the normal $n$ changes its sign on opposite faces of $Y$.

The sufficient condition is obtained by applying Lax–Milgram theorem with respect to $H^{1}_{\#}(Y)/\mathbb{R}$. Indeed, $\displaystyle a(u,v)=\int_{Y}A(y)\nabla u(y)\cdot\nabla v(y)\,dy$ is a coercive continuous bilinear form on $H^{1}_{\#}(Y)/\mathbb{R}$ by uniform ellipticity. Furthermore, the map $F\mathrel{\mathop{\mathchar 58\relax}}H^{1}_{\#}(Y)/\mathbb{R}\to\mathbb{R}$, defined by $\displaystyle F(\phi)=\int_{Y}g(y)\phi(y)\,dy,$ is a well defined bounded linear functional on $H^{1}_{\#}(Y)/\mathbb{R}$ because $g$ is a function of zero mean over $Y$. Indeed, for all $\phi\in H^{1}_{\#}(Y)$, if we let $\displaystyle\bar{\phi}=\fint_{Y}\phi(y)\,dy$ denote the mean value of $\phi$ over $Y$, then, we get

[TABLE]

where we used the Poincaré–Wirtinger inequality in the last inequality. This implies that the map $F\mathrel{\mathop{\mathchar 58\relax}}H^{1}_{\#}(Y)/\mathbb{R}\to\mathbb{R}$ defined above is bounded in the norm $\mathinner{\!\left\lVert\nabla\cdot\right\rVert}_{L^{2}(Y)}$ as claimed. Hence, by Lax–Milgram’s theorem, there exists a unique solution $v\in H_{\#}^{1}(Y)/\mathbb{R}$ such that

[TABLE]

∎

By using Lemma 4.1, we can find the solution of the limit equation. Let us consider the equations that arise when we consider the first three terms of the series in (4.5).

• $\varepsilon^{-2}$:

[TABLE]

It is a partial differential equation with respect to $y$ in $Y$ (here $x$ is just a parameter). By the uniqueness of the solution up to an additive constant, we deduce that

[TABLE]

• $\varepsilon^{-1}$:

[TABLE]

The necessary and sufficient condition of existence is satisfied. Thus, by (4.7), $u_{1}$ (seen as an element of $H^{1}_{\#}(Y)/\mathbb{R}$) depends linearly on $\nabla_{x}u(x)$. In particular, if we let $(e_{i})_{1\leq i\leq N}$ denote the canonical basis of ${{\mathbb{R}}^{N}}$, then it is easy to check that

[TABLE]

where $w_{i}$ is the solutions of the following auxiliary problems (cell problems) for $i=1,\cdots,N$:

[TABLE]

The functions $w_{i}$ are usually called the correctors.

• \varepsilon^{0}\:

[TABLE]

By using Lemma 4.1, the necessary and sufficient condition of existence of the solution $u_{2}$ is

[TABLE]

By employing the use of the representation formula (4.9), we can rewrite $u_{1}$ in terms of $\nabla_{x}u(x)$:

[TABLE]

In other words, we have succeeded in identifying the the homogenized problem

[TABLE]

where the homogenized tensor $A^{*}$ is defined by

[TABLE]

or, integrating by parts

[TABLE]

Indeed, the cell problems (4.10) yield

[TABLE]

Remark 4.2**.**

The formula for $A^{\ast}$ is not fully explicit because cell problems (4.10) must be solved. However $A^{\ast}$ does not depend on $\Omega$, nor $f$, nor the boundary conditions. It only characterizes the microstructure. Later, we shall compute explicitly some examples of $A^{\ast}$.

Under mild smoothness assumptions on the data, one can justify the expansion in $H^{1}(\Omega)$ [BLP1978, JKO1995].

Theorem 4.3**.**

Assume that the homogenized solution $u$ is smooth. Then the following expansion holds in $H^{1}(\Omega)$:

[TABLE]

In particular

[TABLE]

Remark 4.4** (Rigorous justification).**

Employing a formal asymptotic expansion is a very useful method. However we don’t know a priori whether the solution of the microscopic equation can be expanded as (4.2). We refer the interested reader to Tartar’s method [MT1997] and the two-scale convergence method[Ng1989, Al1992] for a rigorous mathematical justification.

Remark 4.5** (Homogenized coefficients $A^{\ast}$).**

In dimension $N=1$, the explicit formula for $A^{\ast}$ is the so-called harmonic mean. In dimension $N\geq 2$, there is no explicit formula for $A^{\ast}$, which has to be computed numerically. Nevertheless, one can obtain explicit bounds on $A^{\ast}$.

Remark 4.6**.**

Homogenization works for non-periodic media too (H-convergence or G-convergence).

Remark 4.7** (Asymptotic expansions for the stress).**

We assume that

[TABLE]

where $\sigma_{i}(x,y)$ is a function of the two variables $x$ and $y$, periodic in $y$ with period $Y=(0,1)^{N}$. Plugging this series in the equation (4.1), we find

[TABLE]

On the other hand,

[TABLE]

and

[TABLE]

One can prove that $\tau$ is the solution of the dual cell problem.

4.2 Composite materials

Composite materials are ubiquitous in engineering, mechanics and physics and their effective properties can be understood through homogenization theory [Al2002, ch1000, MI2001]. In what follows, we identify a composite material by its homogenized tensor $A^{\ast}$. We restrict ourselves to two-phase composites. We mix two isotropic constituents $A(y)=\alpha\chi(y)+\beta(1-\chi(y))$, where $\chi\mathrel{\mathop{\mathchar 58\relax}}Y\to\{0,1\}$ is a characteristic function. Let $\displaystyle\theta=\int_{Y}\chi(y)dy$ be the volume fraction of phase $\alpha$ and $(1-\theta)$ be that of phase $\beta$.

We focus on the characterization of $G_{\theta}$ defined as follows:

Definition 4.8** (The set of all homogenized tensors $G_{\theta}$).**

Let $G_{\theta}$ be the set of all homogenized tensors $A^{\ast}$ obtained by homogenization of the two phases $\alpha$ and $\beta$ in proportions $\theta$ and $(1-\theta)$.

Remark 4.9**.**

Of course, we have $G_{0}=\{\beta{\rm Id}\}$ and $G_{1}=\{\alpha{\rm Id}\}$. However, $G_{\theta}$ is usually a (very) large set of tensors (corresponding to different choices of $\chi(y)$).

4.2.1 Lamination for two phase composites

For two phase composites, the density $\theta(x)$, as well as the homogenized tensor $A^{\ast}(x)$, depends on the position $x$. For two-phase mixtures, an explicit characterization of $G_{\theta}$ is possible by the variational principle of Hashin and Shtrikman [HS1963]. We make the following assumptions:

• (i)

Linear model of conduction or membrane stiffness (it is more delicate for linearized elasticity and very few results are known in the non-linear case).

• (ii)

Perfect interfaces between the phases (continuity of both displacement and normal stress), no possible effects of delamination or debonding.

In dimension one, the cell problem (4.10) reads:

[TABLE]

The solution computed explicitly as follows:

[TABLE]

By (4.13), we know that $\displaystyle A^{\ast}=\int_{0}^{1}A(y)(1+w^{\prime}(y))^{2}\,dy$, which yields the harmonic mean of $A(y)$:

[TABLE]

Therefore, if we choose $A(y)=\alpha\chi(y)+\beta(1-\chi(y))$, then homogenized tensor of any two-phase material is just

[TABLE]

This formula tells us that, in one dimension, the homogenized tensor depends on the characteristic function $\chi$ by means of its volume fraction $\theta$ only.

In dimension $N\geq 2$, we cannot express $A^{\ast}$ explicitly in general as mentioned in Remark 4.5. However it is possible under the following special case. We consider parallel layers of two isotropic phases $\alpha$ and $\beta$, orthogonal to the direction $e_{1}$. Assume that $A^{\varepsilon}$ depends only on $y_{1}$. Let

[TABLE]

We denote by $A^{\ast}$ the homogenized tensor of $A(y)=\left(\alpha\chi(y_{1})+\beta(1-\chi(y_{1}))\right)I$. Then we obtain the following lemma. This lemma is a simple case of the more general Lemma 4.12.

Lemma 4.10**.**

Define $\displaystyle\lambda_{\theta}^{-}=\left(\frac{\theta}{\alpha}+\frac{1-\theta}{\beta}\right)^{-1}$ and $\lambda_{\theta}^{+}=\theta\alpha+(1-\theta)\beta$. Then we have

[TABLE]

Remark 4.11** (Interpretation (resistance $=$ inverse of conductivity)).**

In the context of electrical conductivity, the harmonic mean is the effective conductivity of a mixture of conductors placed in series (in the direction $e_{1}$), while the arithmetic mean is the effective conductivity of a mixture of conductors placed in parallel ( in any direction orthogonal to $e_{1}$).

Lemma 4.12** (Simple laminate of two non-isotropic phases).**

The homogenized tensor $A^{\ast}$ of a simple laminate made of $A$ and $B$ in proportions $\theta$ and $(1-\theta)$ in the direction $e_{1}$ is

[TABLE]

Moreover, if we assume that $(A-B)$ is invertible, then this formula is equivalent to

[TABLE]

Proof.

Recall that by definition (4.13)

[TABLE]

namely

[TABLE]

Consequently, for any $\xi\in\mathbb{R}^{N}$, we have

[TABLE]

where $\displaystyle w_{\xi}(y)=\sum_{i=1}^{N}\xi_{i}w_{i}(y)$ is the solution of

[TABLE]

Defining $u(y)=\xi\cdot y+w_{\xi}(y)$, we seek a solution $u$ such that the gradient of $u$ is constant in each phase,

[TABLE]

Thus, we have

[TABLE]

where $c_{a}$ and $c_{b}$ are constant vectors.

Let $\Gamma$ be the interface between the two phases. By continuity of (4.18) through the interface $\Gamma$, we have

[TABLE]

Since $c_{a}$ and $c_{b}$ are constant vectors, by (4.19) we have

[TABLE]

Since $(x-y)$ is orthogonal to $e_{1}$, there exists a real number $t\in\mathbb{R}$ such that $b-a=te_{1}$.

Moreover, by continuity of the flux $A(y)\nabla u\cdot n$ through the interface $\Gamma$, we have

[TABLE]

In particular, it implies -${\mathop{\mathrm{div}}}(A(y)\nabla u)=0$ in the weak sense.

Since $b-a=te_{1}$, (4.20) yields the following value for $t$:

[TABLE]

Since $w_{\xi}$ is periodic, it satisfies $\displaystyle\int_{Y}\nabla w_{\xi}\ dy=0$, thus by the definition of $u$ we have

[TABLE]

On the other hand, by (4.17) and the definition of $u$ we have

[TABLE]

Thus we obtain

[TABLE]

Since $b=a+te_{1}$ with $\displaystyle t=\frac{(A-B)a\cdot e_{1}}{Be_{1}\cdot e_{1}}$, we find

[TABLE]

Then, a simple computation gives

[TABLE]

The other formula is a consequence of the following fact: if $M$ is invertible, then

[TABLE]

where $c\in\mathbb{R}$ and $e$ is a unit vector in $\mathbb{R}^{N}$ which determines the direction of the lamination. ∎

The composite $A^{\ast}$ is said to be a single lamination in the direction $e_{1}$ of the two phases $A$ and $B$ in proportions $\theta$ and $(1-\theta)$ (see Figure 18). By varying the proportion $\theta$ and the direction $e_{1}$, we obtain a whole family of composite materials. This family can still be enlarged by laminating again these simple laminates. Then we laminate again the preceding composite with always the same phase $B$.

A sequential laminate is obtained by an iterative process of lamination where the previous laminate is laminated again with a single pure phase (always the same one). By using the special form of (4.16) (which does not deliver directly the value of $A^{\ast}$, contrary to (4.15)), the iterative or sequential laminate can be explicitly characterized. Let $(e_{i})_{1\leq i\leq p}$ be a collection of unit vectors and $(\theta_{i})_{1\leq i\leq p}$ be proportions in $[0,1]$. By (4.16) a simple laminate $A_{1}^{\ast}$ of $A$ and $B$ in proportions $\theta$, $(1-\theta)$ is

[TABLE]

This simple laminate $A_{1}^{\ast}$ can again be laminated with phase $B$, in direction $e_{2}$ and in proportions $\theta_{2}$, $(1-\theta_{2})$ respectively, to obtain a new laminate denoted by $A_{2}^{\ast}$. By induction, we obtain $A_{p}^{\ast}$ by lamination of $A_{p-1}^{\ast}$ and $B$, in direction $e_{p}$ and in proportions $\theta_{p}$, $(1-\theta_{p})$, respectively. Then the homogenized tensor $A_{p}^{\ast}$ is

[TABLE]

Replacing $(A_{p-1}^{\ast}-B)^{-1}$ in (4.21) by the similar formula defining $(A_{p-2}^{\ast}-B)^{-1}$, and so on up to $A_{0}^{\ast}\equiv A$, we obtain a formula of the same type as (4.16), namely,

[TABLE]

We remark that we always laminate an intermediate laminate with the same phase $B$. In other words, the other phase $A$ is coated by several layers of $B$. One can say that $B$ plays the role of a matrix phase, and $A$ plays the role of a core phase. Globally, $A^{\ast}$ can be seen as a mixture of $A$ and $B$ in different layers having a large separation of scales (see Figure 19).

Let us define rank-$p$ sequential laminate with matrix $B$ and inclusion $A$.

Lemma 4.13** (rank-$p$ sequential laminate).**

If we laminate $p$ times with $B$, we obtain a rank-$p$ sequential laminate with matrix $B$ and inclusion $A$, in proportions $(1-\theta)$ and $\theta$, is defined by

[TABLE]

with $\displaystyle\sum_{i=1}^{p}m_{i}=1$ and $m_{i}\geq 0$, $1\leq i\leq p$.

Proof.

By (4.22) we already have

[TABLE]

We make the change of variables

[TABLE]

which is indeed one-to-one with the constraints on the $m_{i}$’s and the $\theta_{i}$’s. ∎

Of course the same can be done when exchanging the roles of $A$ and $B$.

Lemma 4.14**.**

A rank-$p$ sequential laminate with matrix $A$ and inclusion $B$, in proportions $\theta$ and $(1-\theta)$, is defined by

[TABLE]

with $\displaystyle\sum_{i=1}^{p}m_{i}=1$ and $m_{i}\geq 0$, $1\leq i\leq p$.

Remark 4.15**.**

Sequential laminates form a very rich and explicit class of composite materials which, as we shall see, completely describes the boundaries of the set $G_{\theta}$.

4.2.2 Characterization of $G_{\theta}$

From now on, we assume that the microscopic tensor $A(y)$ is symmetric. Then $A^{\ast}$ is also symmetric. Furthermore, $A^{\ast}$ is characterized by the following variational principle:

[TABLE]

Indeed, if $w_{\xi}$ is a minimizer of (4.23), then it satisfies the Euler optimality condition

[TABLE]

By linearity, we have $w_{\xi}=\sum_{i=1}^{N}\xi_{i}w_{i}$, where $w_{i}$ ($i=1,\cdots,N$) denotes the solution of (4.10), and thus, by (4.13) we get

[TABLE]

By using the variational principle of $A^{\ast}$ (4.23), we can obtain arithmetric and harmonic mean bounds for $A^{\ast}$.

Lemma 4.16** (Arithmetic and harmonic mean bounds).**

Any homogenized tensor $A^{\ast}$ satisfies the arithmetic mean bound

[TABLE]

and the harmonic mean bound

[TABLE]

Proof.

Taking $w=0$ in the variational principle (4.23), we deduce the arithmetic mean bound. For the harmonic mean bound we enlarge the minimization space as follows. Indeed, since $\displaystyle\int_{Y}\nabla w\,dy=0$, we replace $\nabla w$ with any vector field $\zeta(y)$ with zero-average on $Y$

[TABLE]

The Euler equation for the minimizer $\zeta_{\xi}(y)$ of this convex problem is

[TABLE]

where $\lambda\in\mathbb{R}$ is the Lagrange multiplier for the constraint $\displaystyle\int_{Y}\zeta\,dy=0$. Thus

[TABLE]

and

[TABLE]

∎

Lemma 4.16 can be improved for two-phase composites. Next, we consider two isotropic phases $A=\alpha{\rm Id}$ and $B=\beta{\rm Id}$ with $0<\alpha<\beta$.

Theorem 4.17** (Hashin and Shtrikman bounds [HS1963, TA2000]).**

The set $G_{\theta}$ of all homogenized tensors obtained by mixing $\alpha$ and $\beta$ in proportions $\theta$ and $(1-\theta)$ is the set of all symmetric matrices $A^{\ast}$ with eigenvalues $\lambda_{1},\cdots,\lambda_{N}$ such that

[TABLE]

[TABLE]

[TABLE]

Furthermore, these so-called Hashin and Shtrikman bounds are optimal and attained by rank-$N$ sequential laminates.

Proof.

We first show that all matrices satisfying these inequalities (Hashin-Shtrikman bounds) belong to $G_{\theta}$. Let us start by showing that the upper bound (4.26) is attained by sequential laminates. Take a matrix $A^{\ast}$ such that

[TABLE]

Define a rank-$N$ sequential laminate $A_{L}^{\ast}$ of matrix $\beta$ and inclusion $\alpha$, with lamination directions being the (orthogonal) eigenvectors of $A^{\ast}$. By Lemma 4.13 we have

[TABLE]

We obtain $A^{\ast}=A^{\ast}_{L}$ if we can choose the $m_{i}$’s such that

[TABLE]

that is,

[TABLE]

We check that $0<m_{i}<1$ is equivalent to $\lambda_{\theta}^{-}<\lambda<\lambda_{\theta}^{+}$ and that

[TABLE]

Thus any matrix on the upper bound (4.26) is a rank-$N$ sequential laminate with matrix $\beta$ and inclusion $\alpha$. The same proof works for the lower bound (4.25) upon exchanging the role of $\alpha$ (now the matrix) and $\beta$ (now the inclusion).

A simple but lengthy computation shows that all the matrices satisfying the inequalities (4.24), (4.25) and (4.26) can be obtained as a rank-$N$ sequential laminate of two suitable matrices, one realizing the equality in the upper bound (4.26) and the other realizing the equality in the lower bound (4.25) (see the full proof of [Al2002, Theorem 2.2.13] for the details). It remains to prove that the lower and upper Hashin–Shtrikman bounds hold true. To establish the lower bound (4.25) we introduce the so-called Hashin and Shtrikman variational principle. Main idea is to use Fourier analysis and Plancherel theorem.

By definition of $A^{\ast}$, for $\xi\in\mathbb{R}^{N}$, we have

[TABLE]

Subtracting a reference material $\alpha$,

[TABLE]

We use convex duality (or Legendre transform): for any symmetric positive definite matrix $K$, the following holds

[TABLE]

Since $0<\alpha<\beta$, we apply the formula (4.27) at each point in $Y$. Then we get

[TABLE]

which becomes an inequality if we restrict the minimization to constant $\eta$ in $Y$

[TABLE]

On the other hand, because of periodicity, $\int_{Y}\nabla w\,dy=0$ which implies

[TABLE]

Overall, we obtain that, for any $\eta\in\mathbb{R}^{N}$,

[TABLE]

where $g(\chi,\eta)$ is a so-called non-local term, defined by

[TABLE]

We can now use Fourier analysis to compute $g(\chi,\eta)$. By periodicity, both $\chi$ and the test function $w$ can be written as Fourier series:

[TABLE]

Since $\chi$ and $w$ are real-valued, their Fourier coefficients satisfy

[TABLE]

The gradient of $w$ at $y\in Y$ is given by

[TABLE]

Then, Plancherel formula yields

[TABLE]

Notice that minimizing in $w(y)\in H^{1}_{\#}(Y)$ is equivalent to minimizing in $\hat{w}(k)\in\mathbb{C}$. For $k\neq 0$ the minimum is achieved by

[TABLE]

and we deduce that

[TABLE]

where $M$ is a symmetric non-negative matrix defined by

[TABLE]

Since, by Plancherel theorem, we have

[TABLE]

we deduce that the trace of $M$ is equal to $1$.

Substituting (4.29) to (4.28), for any $\xi,\eta\in\mathbb{R}^{N}$,

[TABLE]

The minimum (in $\xi$) of the inequality (4.30) is obtained when

[TABLE]

Then we deduce

[TABLE]

Thus, we have

[TABLE]

Taking the trace of this matrix inequality (4.31), and recalling that $\mathop{\mathrm{tr}}M=1$, we obtain the lower Hashin–Shtrikman bound. The proof of the upper bound is similar. ∎

4.3 The elasticity setting

In what follows, let us consider the elasticity setting. The homogenization method can be generalized to the elasticity setting. However, an explicit characterization of $G_{\theta}$ is still lacking in the elasticity setting.

We set

[TABLE]

with the identity matrix $I_{2}$, and $\kappa_{A,B}=\lambda_{A,B}+2\mu_{A,B}/N$. We assume $B$ to be weaker than $A$:

[TABLE]

We work with stresses rather than strains, thus we use inverse elasticity tensors. The similar results of the two-phase composites in the elasticity setting as follows (in details, see [Al2002, Section 2.3]):

Lemma 4.18** (Sequential laminates in elasticity).**

The Hooke’s law of a simple laminate of $A$ and $B$, in proportions $\theta$ and $(1-\theta)$, respectively, in the direction $e$, is

[TABLE]

where $f^{c}_{A}(e)$ is the tensor, defined, for any symmetric matrix $\xi$, by

[TABLE]

Proposition 4.19** (Reiterated lamination formula).**

A rank-$p$ sequential laminate with matrix $A$ and inclusions $B$, in proportions $\theta$ and $(1-\theta)$, respectively, in the directions $(e_{i})_{1\leq i\leq p}$ with parameter $(m_{i})_{1\leq i\leq p}$ such that $0\leq m_{i}\leq 1$ and $\displaystyle\sum_{i=1}^{p}m_{i}=1$, is given by

[TABLE]

Theorem 4.20** (Hashin–Shtrikman bounds in elasticity).**

Let $A^{\ast}$ be a homogenized elasticity tensor in $G_{\theta}$ which is assumed isotropic

[TABLE]

Its bulk $\kappa_{\ast}$ and shear $\mu_{\ast}$ moduli satisfy

[TABLE]

[TABLE]

[TABLE]

Furthermore, the two lower bounds, as well as the two upper bounds are simultaneously attained by a rank-$p$ sequential laminate with $p=3$ if $N=2$, and $p=6$ if $N=3$.

Proof.

We refer to [Al2002, Theorem 2.3.13] ∎

Remark 4.21**.**

These bounds do not characterize all possible isotropic homogenized tensors $A^{\ast}$ in $G_{\theta}$. In other words, there exist isotropic elasticity tensors with moduli satisfying these bounds that are not composite materials obtained by mixing phases $A$ and $B$ in proportions $\theta$, $(1-\theta)$, respectively.

Proposition 4.22** (Hashin–Shtrikman optimal energy bound).**

Let $G_{\theta}$ be the set of all homogenized elasticity tensors obtained by mixing the two phases $A$ and $B$ in proportions $\theta$ and $(1-\theta)$. Let $L_{\theta}$ be the subset of $G_{\theta}$ made of sequential laminated composites. For any stress $\sigma$,

[TABLE]

Furthermore, the minimum is attained by a rank-$N$ sequential laminate with lamination directions given by the eigendirections of $\sigma$.

Remark 4.23**.**

An optimal tensor $A^{\ast}$ can be interpreted as the most rigid composite material in $G_{\theta}$ able to sustain the stress $\sigma$. $HS(\sigma)$ is called Hashin–Shtrikman optimal energy bound. In practical conclusion, $G_{\theta}$ can be replaced by $L_{\theta}$ for compliance minimization.

4.4 Numerical applications

Let us consider the case of parametrized periodicity cells. For example, the square cell with a rectangular hole (as used in the seminal work of Bendsøe and Kikuchi [BK1988]), parametrized by $m_{1}$, $m_{2}$, and denoted by $Y(m)$ (see Figure 21).

We compute the so-called correctors or cell solutions:

[TABLE]

where $e_{ij}=(e_{i}\otimes e_{j}+e_{j}\otimes e_{i})/2$ is a basis of the symmetric tensors of order $2$, and $n$ is the normal to the hole’s boundary $\Gamma_{\rm int}$ in $Y(m)$. Hence we find a unique solution (up to an additive translation) $w_{ij}\in H^{1}_{\#}(Y(m),\mathbb{R}^{2})$ to the variational formulation:

[TABLE]

The tensor $A^{\ast}$ is then given by

[TABLE]

4.5 Exercises

Problem 4.5.1**.**

Compute numerically $A^{\ast}$ for various parameters $m_{1}$, $m_{2}$.

Problem 4.5.2**.**

Orthotropic composite: check numerically that $A^{\ast}_{1112}=A^{\ast}_{2212}=0$, then prove it theoretically.

Problem 4.5.3**.**

Check that if $m_{1}\rightarrow 1$, then $A^{\ast}$ is close to the formula of a rank-$1$ laminate.

Problem 4.5.4**.**

What happens if $m_{1}\rightarrow 0$? Is $A^{\ast}$ close to $A$?

Problem 4.5.5**.**

If $m_{1}=m_{2}$, is $A^{\ast}$ isotropic ?

Problem 4.5.6**.**

Check numerically that $A^{\ast}$ is isotropic for a honeycomb structure with hexagonal holes.

Chapter 5 Topology optimization by the homogenization method

5.1 Why topology optimization?

Shape optimization consists in “shape tracking” algorithms, hence the method cannot change the topology, such as the number of holes in the case of $2$-dimension. On the other hand, topology optimization consists in “shape capturing” algorithms, that allow us to consider the optimization in a wider class which includes the different topological properties. There are several methods of topology optimization but we focus on just one, called the homogenization method (see [Al2002, BS2003] and references therein). In the following of this chapter, we introduce a model problem for topology optimization with a constraint on the volume of holes and study it using the method of homogenization as mentioned in Chapter 4.

5.2 Homogenization method in the conductivity setting

In this section, we apply the homogenization method in the conductivity setting. As we mentioned in Section 3.2, there is no minimizer for the corresponding minimizing problem (we will explain the detail below) in general. To solve the problem, we introduce a set of generalized shapes as the limit of classical shapes. More precisely, goals of the homogenization method for topology optimization are following:

• •

To introduce the notion of generalized shapes made of composite material,

• •

To show that those generalized shapes are limits of sequences of classical shapes (in the sense of homogenization),

• •

To compute the generalized objective function and its gradient,

• •

To prove an existence theorem for optimal generalized shapes,

• •

To deduce new numerical algorithms for topology optimization.

In order to consider the limit of classical shapes, we recall one of the main results of homogenization theory. For the details of the proof, see [Al2002, Theorem 1.2.16 and 2.1.2].

Theorem 5.1**.**

Let $\Omega$ be a bounded domain in $\mathbb{R}^{N}$ and $\left(\chi_{\varepsilon}(x)\right)_{\varepsilon>0}$ be a sequence of characteristic functions in $\Omega$. Set $A_{\varepsilon}(x)=\alpha\chi_{\varepsilon}(x)+\beta(1-\chi_{\varepsilon}(x))$ for $x\in\Omega$. Then there exists a subsequence, still denoted by $\chi_{\varepsilon}$, a density $0\leq\theta(x)\leq 1$ and a homogenized tensor $A^{*}$ such that

[TABLE]

and $A_{\varepsilon}$ converges in the sense of homogenization to $A^{*}$, i.e., for all $f\in L^{2}(\Omega)$, the solution $u_{\varepsilon}$ of the problem

[TABLE]

converges strongly in $L^{2}(\Omega)$ to the solution $u$ of the homogenized problem

[TABLE]

Furthermore, for almost all $x\in\Omega$, $A^{*}(x)$ belongs to the set $G_{\theta(x)}$ defined in Definition  4.8.

5.2.1 Relaxed problem for the conductivity model

In this subsection, we introduce the conductivity model and derive generalized shapes as the limits of classical shapes by the homogenization method. We impose a simplifying assumption: the “holes” (with a Neumann, free boundary condition) are filled with a weak (“ersatz”) material $\alpha$, while the other space is filled with a material $\beta$, that is, $\alpha<\beta$. We consider a membrane we introduced in Chapter 3 with two thicknesses $\alpha$ and $\beta$, whose distribution is give by $h_{\chi}=\alpha\chi+\beta(1-\chi)$, where $\chi$ is a characteristic function which denotes the position of the holes. If $f\in L^{2}(\Omega)$ is the applied load, the displacement $u_{\chi}$ satisfies

[TABLE]

We will optimize the membrane’s shape amounts by considering the minimizing problem

[TABLE]

with

[TABLE]

where $V_{\alpha}$ is a given positive constant, denoting the volume of the holes and

[TABLE]

Since there is no minimizer of (5.2) in general, we apply Theorem 5.1 to introduce a set of generalized shapes. Let $\chi_{\varepsilon}$ be a sequence (minimizing sequence of (5.2) or not) of characteristic functions. Applying Theorem 5.1, we see that there exist $\theta$ and $A^{*}$ such that

[TABLE]

and

[TABLE]

where $u$ is the solution of (5.1) and $j$ is defined by (5.3). From this convergence result, we define the set of admissible homogenized shapes

[TABLE]

and consider the following relaxed or homogenized optimization problem:

[TABLE]

We easily check that $\mathcal{U}_{\rm ad}\subset\mathcal{U}^{*}_{\rm ad}$ if we identify $\chi\in\mathcal{U}_{\rm ad}$ with the pair $(\chi,\alpha\chi{\rm Id}+\beta(1-\chi){\rm Id})\in\mathcal{U}^{*}_{\rm ad}$. The inclusion implies that we have enlarged the set of admissible shapes. Moreover, one can prove that the relaxed problem (5.4) always admits an optimal solution, and the homogenized formulation is a relaxation of the original topology optimization problem as follows.

Theorem 5.2**.**

The homogenized formulation is a relaxation of the original topology optimization problem in the sense that:

• •

there exists, at least, one optimal composite shape $(\theta,A^{*})\in\mathcal{U}_{\rm ad}^{*}$, i.e. a minimizer of (5.4),

• •

for any minimizing sequence $\left(\chi_{n}\right)_{n\in\mathbb{N}}$ of (5.2), there exists a minimizer $(\theta,A^{*})\in\mathcal{U}^{*}_{\rm ad}$ of (5.4) such that, up to subsequence, $\chi_{n}$ converges weakly *** in $L^{\infty}(\Omega;[0,1])$ to $\theta$ and $A_{n}\mathrel{\mathop{\mathchar 58\relax}}=\alpha\chi_{n}+\beta(1-\chi_{n})$ converges to $A^{*}$ in the sense of homogenization,

• •

any composite optimal solution $(\theta,A^{*})\in\mathcal{U}_{\rm ad}^{*}$ of (5.4) is the limit of a minimizing sequence of (5.2).

Moreover, the infima of the original and homogenized objective functions coincide

[TABLE]

For the proof of Theorem 5.2, see [Al2002, Theorem 3.2.1]. Theorem 5.2 means that the topology optimization problem is not changed by relaxation. Moreover, close to any optimal composite shape, we are sure to find a quasi-optimal classical shape. This theorem is at the root of new numerical algorithms.

Remark 5.3**.**

The homogenized formulation is similar to a parametric or sizing optimization problem. This is the main reason why the homogenization method is computationally cheap and works like a shape capturing algorithm. Moreover, computing gradients or optimality conditions are thus very simple. On the other hand, the design parameters $(\theta,A^{\ast})$ are quite complicated. Another further (drastic and unjustified) simplification is to suppress the parameter $A^{\ast}$ and to keep only the material density $\theta$. This is the main idea of the SIMP method [BS2003], we will mention in Section 5.4.

5.2.2 Optimality conditions

In Section 5.2.1, we introduced the relaxed formulation (5.4) of the original optimization problem (5.2). One of the advantages of the relaxed problem is that we always have the existence of a minimizer of (5.4). There is also another advantage: we can get the optimality condition for the relaxed problem since it is possible to perform variations of composite designs. In this subsection, we will consider the optimality condition for problem (5.4).

We now compute the gradient of the following objective function

[TABLE]

where $u$ is the solution to (5.1) and $u_{0}$ is a given function in $L^{2}(\Omega)$. We introduce the adjoint state $p$ of $u$ as the unique solution in $H^{1}_{0}(\Omega)$ of

[TABLE]

By the use of the adjoint state $p$, we can obtain the derivative of the functional $J$.

Proposition 5.4**.**

Let $\alpha>0$ and $\mathcal{M}_{\alpha}$ be the set of symmetric positive definite matrices $M$ such that $M\geq\alpha{\rm Id}$. The functional $J$ is differentiable with respect to $A^{*}$ in $L^{\infty}(\Omega;\mathcal{M}_{\alpha})$, and its derivative is

[TABLE]

i.e.,

[TABLE]

where $u\in H^{1}_{0}(\Omega)$ is the unique solution of (5.1) and $p$ is the adjoint state (5.5) of $u$.

Proof.

We can prove the proposition by the same strategy that we used in Lemma 3.8, 3.9 and Theorem 3.10. Thus we omit the proof. ∎

Remark that the partial derivative with respect to $\theta$ vanishes because $\theta$ appears only in the constraint of $A^{*}$. Moreover, we can also consider the Lagrangian

[TABLE]

where $(A^{*},v,q)\in L^{\infty}(\Omega;\mathcal{M}_{\alpha})\times H^{1}_{0}(\Omega)\times H^{1}_{0}(\Omega)$. The partial derivatives of $\mathcal{L}$ with respect to $q$ and $v$ yield the state and adjoint state respectively. Furthermore, the functional $J$ is also differentiable with respect to $A^{*}$ with derivative

[TABLE]

The essential consequence of this section is the following optimality condition.

Theorem 5.5**.**

Let $(\theta,A^{\ast})$ be a global minimizer of $J$ in $\mathcal{U}_{\rm ad}^{\ast}$ which admits $u$ and $p$ as state and adjoint. Then there exists $(\tilde{\theta},\tilde{A}^{\ast})$, another global minimizer of $J$ in $\mathcal{U}_{\rm ad}^{\ast}$, which admits the same state and adjoint $u$ and $p$, and such that $\tilde{A}^{\ast}$ is a rank-$1$ simple laminate.

Proof.

We fix $\theta\in L^{\infty}(\Omega;[0,1])$ with

[TABLE]

We remark that by Theorem 4.17, the set

[TABLE]

is a convex set. Then, by Theorem 2.17 and Proposition 5.4, we have

[TABLE]

We easily check that the above inequality is equivalent to the point-wise constraint

[TABLE]

for almost all $x\in\Omega$. Fix $x\in\Omega$ which satisfies (5.6). If $\nabla u(x)$ or $\nabla p(x)$ vanishes, then any $A^{\ast}\in G_{\theta(x)}$ is optimal in (5.6) (the fact that $u$ and $p$ does not change at these points is more delicate to establish and we refer to [Al2002, TA2000] for details). Otherwise, we define two unit vectors

[TABLE]

Then (5.6) is equivalent to finding a minimizer $A^{*}(x)$ of

[TABLE]

A lower bound is easily seen to be

[TABLE]

where $\lambda^{\pm}_{\theta}$ is defined in Theorem 4.17. We can see that the lower bound is attained by a matrix $A^{1}$ corresponding to a rank-1 laminate (see [Al2002, Remark 2.2.14]). More precisely, $A^{1}$ satisfies

[TABLE]

Thus we obtain

[TABLE]

Moreover, if $A^{\ast}$ is any optimal tensor, then, $A^{\ast}$ must also satisfy

[TABLE]

Indeed, if (5.9) does not hold true, the bounds on eigenvalues in Lemma 4.17 imply the strict inequality

[TABLE]

which is a contradiction with (5.6) and (5.7). By (5.8) and (5.9), we see that

[TABLE]

for almost all $x\in\Omega$. Therefore any optimal tensor $A^{*}$ can be replaced by the rank-1 simple laminate $A^{1}$ without changing $u$ and $p$. ∎

Remark 5.6**.**

Theorem 5.5 implies that in the definition of $\mathcal{U}_{\rm ad}^{\ast}$, the set $G_{\theta}$ can be replaced by its simpler subset of rank-$1$ simple laminates. We actually use this simplification in the numerical algorithms. We remark that this simplification holds true for other objective functions as well. However, it does not hold for multiple loads optimization in general. For the details concerning the conditions that the objective function $J$ must satisfy in order to obtain the optimality conditions, see [Al2002, Chapter 3.2.2].

As stated in Remark 5.6, we can simplify the admissible set of the minimization problem (5.4). We consider the parametrization of rank-1 laminates. For simplicity, we consider the case $N=2$. A rank-1 laminate is defined by

[TABLE]

where the angle $\phi\in[0,\pi]$ determines the orientation of the unit cell. Hence the admissible set is rewritten as

[TABLE]

If we set $J(\theta,\phi)\mathrel{\mathop{\mathchar 58\relax}}=J(\theta,A^{*}(\theta,\phi))$, then the derivative of the objective function $J(\theta,\phi)$ follows by Proposition 5.4 immediately.

Proposition 5.7**.**

The objective function $J(\theta,\phi)$ is differentiable with respect to $(\theta,\phi)$ in $\mathcal{U}_{\rm ad}^{L}$, and its partial derivatives are

[TABLE]

5.2.3 Numerical algorithm

In this subsection, we show the numerical algorithm to seek the optimal shape $\theta$ of (5.4). As stated in Section 5.2.2, we can treat $J(\theta,\phi)$ instead of $J(\theta,A^{*})$ if $N=2$. We explain the projected gradient algorithm for the minimization of $J(\theta,\phi)$ by the use of Proposition 5.7.

For the details about the multiplier $\ell_{k}$, we refer to Section 3.5. In the following, we will consider two simpler self-adjoint cases. Finally, we will show the algorithm to obtain the optimal shape which is close to the classical shapes.

First example of a self-adjoint case

A first example is the minimization of the torsional rigidity (maximization of compliance)

[TABLE]

where $u$ is the solution of

[TABLE]

In this case, the adjoint state $p$ is the solution of

[TABLE]

i.e., the adjoint state is just $p=u$.

Before applying a numerical algorithm, we simplify the minimization problem by using the argument in Section 5.2.2. By the similar argument as in Proposition 5.4, we get

[TABLE]

Hence we have to decrease $A^{\ast}$ to minimize $J$. Indeed, by (5.9) any minimizer $(\theta,A^{\ast})$ satisfies

[TABLE]

i.e., the optimal composite is the worst possible conductor. This condition allows us to eliminate the angle $\phi$ and it remains to optimize with respect to $\theta$ only, i.e., instead of (5.11) we have to consider the state of the problem

[TABLE]

Moreover, we recall that the state $u$, which is the solution of (5.12), is characterized as the minimizer of the corresponding energy, i.e.,

[TABLE]

On the other hand, by the weak form of the state $u$, we have

[TABLE]

and thus

[TABLE]

Therefore, we can rewrite the minimizing problem (5.10) as

[TABLE]

where the minimum in the right hand side of the above equation is taken over the set

[TABLE]

Furthermore, since the function $(\theta,v)\mapsto\lambda_{\theta}^{-}|\nabla v|^{2}$ is convex, there are only global minima by Proposition 2.10.

Numerically, we use an algorithm based on alternate direction minimization (see Section 3.5). We solve in the domain $\Omega=(0,1)^{2}$ with phases $\alpha=1$ and $\beta=2$ respectively. We work with a volume constraint $50\%$ of phase $\alpha$. We initialize with a constant value of $\theta=0.5$ and a constant zero lamination angle $\phi=0$. We perform $30$ iterations. We show the numerical result how the objective function $J$ convergent to some value (Figure 22) and the volume fraction $\theta$ at some iteration numbers (Figure 23). Here, the horizontal axis in Figure 22 means the iteration number.

Second example of a self-adjoint case

A second self-adjoint example is a compliance minimization

[TABLE]

where $u$ is the solution of (5.11). In this case, the adjoint state is $p=-u$.

In order to apply the numerical algorithm, we will simplify the minimizing problem as we did in the first example. By the same argument as in the first example, we see that

[TABLE]

and

[TABLE]

if $(\theta,A^{*})$ is a minimizer of (5.13). Hence the optimal composite is the best possible conductor and we have only to consider the following problem instead of (5.11):

[TABLE]

Therefore, as in the previous section, we can eliminate the dependency on the angle $\phi$ and then optimize with respect to $\theta$ only. We rewrite the optimization problem thanks to the dual energy

[TABLE]

We can obtain the dual energy by the similar calculation in Example 2.30. Thus instead of (5.13) we obtain a double minimization

[TABLE]

where the minimum is taken over the set

[TABLE]

Furthermore, since the function $(\theta,\tau)\mapsto(\lambda^{+}_{\theta})^{-1}|\tau|^{2}$ is convex, there are only global minima by Proposition 2.10.

We apply the same algorithm as in the first example. The setting of the problem is also same as in the first example. Figure 24 shows the numerical result of the volume fraction $\theta$ at some iteration numbers.

Remark 5.8**.**

Thanks to the convexity properties of the functionals, the convergence to a global minimum is guaranteed. In practice, it can be checked by numerical experiments with various initializations converging to the same solution.

If one is interested by shape optimization rather than two-phase optimization, then, in numerical practice, holes can be mimicked by a very weak phase $\alpha$, such as $10^{-3}\beta$. Mathematically, when $\alpha\to 0$ we obtain Neumann boundary conditions on the holes boundaries.

Penalization

By the algorithm stated in the two examples, we obtain optimal shapes in the wider class of composite shapes. Since, in practice, we are rather interested in classical shapes, we choose to use a penalization process to force the density to take values close to [math] or $1$.

In the second algorithm, we note that if $0<\theta_{\rm opt}<1/2$, then $\theta_{\rm pen}<\theta_{\rm opt}$, while $1/2<\theta_{\rm opt}<1$, then $\theta_{\rm pen}>\theta_{\rm opt}$. Hence we see that the density $\theta$ goes to [math] or $1$ when we apply the algorithm.

Example 5.9**.**

We consider the optimal radiator (Figure 25-(a))

[TABLE]

We minimize the temperature where heating takes place

[TABLE]

This is another case of compliance minimization. Thus, the problem is self-adjoint and $p=-u$. We solve in the domain $\Omega=(0,1)^{2}$ with phases $\alpha=0.01$ and $\beta=1$. We work with a volume constraint $50\%$ of phase $\alpha$. We initialize with a value $\theta$ as Figure 25-(b). Figure 25-(c), (d) and (e) show the numerical results of the volume density at some iteration numbers under the above penalization algorithm.

5.3 Homogenization method in the elasticity setting

In this section, we will apply the homogenization method in the elasticity setting. We remark that it is very similar to the conductivity setting but there are some additional hurdles. We shall review the results without proofs, however, the basic ingredients of the homogenization method which we will consider in this section are the sames$\colon$

• •

Introduction of composite designs characterized by $(\theta,A^{*})$,

• •

Hashin-Shtrikman bounds for composites,

• •

Sequential laminates are optimal microstructures for compliance minimization, which we will consider the following.

We remark that, unfortunately, the full set of composites $G_{\theta}$ is unknown as stated in Section 4.3, unlike the case we considered in the Section 5.2.

5.3.1 Introduction of the model of the elasticity and relaxed problem

We introduce the model compliance minimization problem (Figure 26). Let $N=2$ or $3$ and $D\subset\mathbb{R}^{N}$ be a bounded domain. If the Hooke’s law $A$ is isotropic, with positive bulk and shear moduli $\kappa$ and $\mu$, we have

[TABLE]

Let $\Gamma_{D}\subset\partial D$ be the Dirichlet part and $\Gamma_{N}\subset\partial D$ be the Neumann part loaded by $g$. For any domain $\Omega\subset D$ with $\Gamma_{D}$, $\Gamma_{N}\subset\partial\Omega$, the displacement vector field $u\colon\Omega\to\mathbb{R}^{N}$ is defined as the solution of the problem

[TABLE]

where $\Gamma\mathrel{\mathop{\mathchar 58\relax}}=\partial\Omega\setminus(\Gamma_{D}\cup\Gamma_{N})$, $e(u)\mathrel{\mathop{\mathchar 58\relax}}=(\nabla u+(\nabla u)^{t})/2$ and $\lambda\mathrel{\mathop{\mathchar 58\relax}}=\kappa-2\mu/N$. We consider the following minimization problem to obtain the optimal shape such that the weight is minimized and the rigidity is maximized:

[TABLE]

where $\ell$ is a Lagrange multiplier and the infimum is taken over the all subset $\Omega$ of $D$ with $\Gamma_{D}$, $\Gamma_{N}\subset\partial\Omega$.

The shape optimization problem (5.15) can be approximated by a two-phase optimization problem: the original material $A$ and the holes of rigidity $B\approx 0$. Then the Hooke’s law of the mixture in $D$ is rewritten as

[TABLE]

Hence the admissible set becomes

[TABLE]

As in conductivity (membrane) case, we can apply the relaxation approach based on homogenization theory.

We introduce composite structures characterized by a local volume fraction $\theta(x)$ of the phase $A$ and a homogenized tensor $A^{*}(x)$, corresponding to its microstructure. The set of admissible homogenized designs is

[TABLE]

where, for fixed $x$ in $D$, $G_{\theta}(x)$ denotes the set of all possible two-phase composite materials at fixed volume fraction $\theta(x)$. In this case, the homogenized state equation is

[TABLE]

and the homogenized compliance is defined by

[TABLE]

By the above setting, the relaxed or homogenized optimization problem is derived as follows:

[TABLE]

In the elasticity setting, an explicit characterization of $G_{\theta}$ is still lacking, it is a major inconvenience of the problem (5.17). Fortunately, for compliance one can replace $G_{\theta}$ by its explicit subset of laminated composites. The key argument to avoid the knowledge of $G_{\theta}$ is that, thanks to the complementary energy minimization, the compliance can be rewritten as

[TABLE]

The complementary energy is followed by a similar argument in Example 2.30. The shape optimization problem (5.17) thus becomes a double minimization problem

[TABLE]

We will exchange the order of the minimization (5.18). Since the order of minimization is irrelevant, (5.18) can be rewritten as

[TABLE]

The minimization with respect to the design parameters $(\theta,A^{*})$ is local. Hence the above minimization becomes

[TABLE]

For a given stress tensor $\sigma$, the minimization of complementary energy

[TABLE]

is a classical problem in homogenization of finding optimal bounds on the effective properties of composite materials. It turns out that we can restrict ourselves to sequential laminates which form an explicit subset $L_{\theta}$ of $G_{\theta}$ by Proposition 4.22. Recall that in the conductivity setting, it was enough to consider only the case of rank-$1$ laminates. On the other hand, Proposition 4.22 tells us that, in the elasticity setting, rank-$1$ laminates are not enough and $G_{\theta}$ has to be replaced by the set of rank-$N$ sequential laminates instead.

As in the case of the conductivity setting in Section 5.2.1, one can prove that the problem (5.17) is a relaxation of the original shape optimization in the following sense (for the details of the proof, see [Al2002, Theorem 4.1.12]).

Theorem 5.10**.**

The homogenized formulation (5.17) is the relaxation of the original problem (5.15) in the sense where

there exists, at least, one optimal composite $(\theta,A^{*})$ minimizing (5.17),

• 2.

any minimizing sequence of classical shapes $\Omega$ for (5.15) converges, in the sense of homogenization, to a minimizer $(\theta,A^{*})$ of (5.17),

• 3.

the minimal values of the original and homogenized objective functions coincide.

5.3.2 An explicit optimal bound $HS(\sigma)$

As we stated in Section 5.3.1, we can restrict the set $G_{\theta}$ to the sequential laminates $L_{\theta}$. In this subsection, we will show the explicit computation of $HS(\sigma)$ for a special case. We will consider the case $B=0$ for the simplification of algebra. Note that the case $B=0$ is natural since the weak material is actually degenerate.

In two dimension case, we can obtain an explicit formula for the bound (see [Al2002, Theorem 2.3.35]):

Theorem 5.11**.**

Assume that $N=2$, $B=0$ and $\theta\neq 0$. Then for any stress tensor $\sigma$, the optimal bound $HS(\sigma)$ is rewritten as

[TABLE]

where

[TABLE]

and $\sigma_{1}$, $\sigma_{2}$ are the eigenvalues of $\sigma$. Furthermore, an optimal rank-$2$ sequential laminate is given by

[TABLE]

where $m_{i}$ is the parameters appeared in Lemma 4.13.

In the three dimension case, we can also obtain the explicit formula for the bound (see [Al2002, Theorem 2.3.36]), however, it is more complicated than the two dimension case. Hence we restrict ourselves to the simple case of zero Poisson ratio, i.e. $\kappa=2\mu/3$.

Theorem 5.12**.**

Assume that $N=3$, $B=0$ and $\theta\neq 0$. We also assume that the constants $\kappa$ and $\mu$ satisfies $\kappa=2\mu/3$. For any $\sigma$, we label the eigenvalues of $\sigma$ as $|\sigma_{1}|\leq|\sigma_{2}|\leq|\sigma_{3}|$. Then, we obtain (5.20) with

[TABLE]

Furthermore, in the first regime, an optimal rank-$3$ sequential laminate is given by

[TABLE]

and in the second regime, an optimal rank-$2$ sequential laminate is

[TABLE]

where $m_{i}$ is the parameters appeared in Lemma 4.13.

5.3.3 Optimality conditions

We consider the optimality condition for the minimization problem (5.19). In this subsection, we assume the same condition as in Section 5.3.2, i.e., $B=0$ and $\theta\neq 0$. If $(\theta,A^{*},\sigma)$ is a minimizer, then by Proposition 4.22, $A^{*}$ is a rank-$N$ sequential laminate aligned with $\sigma$. Moreover, Proposition 4.19 leads the explicit proportions

[TABLE]

If we consider the case $N=2$, then, by (5.21), we can rewrite the minimization appearing in the integrand (5.19) as

[TABLE]

Hence we obtain the explicit optimality formula for $\theta$ as follows:

[TABLE]

The explicit formula for $\theta$ in the case $N=3$ follows by a similar argument.

5.3.4 Numerical algorithm

In this subsection, we will introduce a numerical algorithm for the minimization problem (5.17). We use the following double “alternating” minimization in $\sigma$ and in $(\theta,A^{*})$:

Since the problem is self-adjoint, we can exchange the problem which we can consider the explicit optimality formula and the fact allows us to consider the above numerical algorithm. The algorithm uses a local microstructure $A^{*}$ and a global density $\theta$. Such algorithm is called micro-macro method.

Remark 5.13**.**

The objective function always decreases. Indeed, since $(\theta_{n},A^{*}_{n})$ minimizes the compliance under the stress $\sigma_{n}$, we have

[TABLE]

On the other hand, $\sigma_{n+1}$ minimizes the elastic complementary energy corresponding to the Hooke’s law $A^{*}_{n}$, we see that

[TABLE]

Combining the above inequalities, we obtain the claim

[TABLE]

We show the numerical results for the case of the cantilever (Figure 27). The optimal shape of the short cantilever, i.e., the domain size 10$\times$20, is displayed on the left of Figure 28. Moreover, the left figure of Figure 29 shows the convergence history of the objective function $J$. Here, the horizontal axis means the iteration number. The right figure of Figure 29 shows the transition of the quantity

[TABLE]

where the index $i$ refers to the cell number. We also show the optimal shape of the medium cantilever, i.e., the domain size is 20$\times$10 in the left of Figure 30.

Penalization

The algorithm, we considered in this subsection, compute composite shapes instead of classical shapes. We thus use a penalization technique to force the density to take values close to [math] or $1$ as in Section 5.2.3. The algorithm is the following$\colon$

Note that if $0<\theta_{\rm opt}<1/2$, then $\theta_{\rm pen}<\theta_{\rm opt}$, while, if $1/2<\theta_{\rm opt}<1$, then $\theta_{\rm pen}>\theta_{\rm opt}$. By using this algorithm, we can obtain the penalized resulting shape of the short cantilever and the medium cantilever (see the right figure of Figure 28 and 30 respectively). We also show the numerical result for the bridge problem (Figure 31 and Figure 32).

5.4 Convexification, “fictitious materials” and SIMP

In the homogenization method, composite materials are introduced, however, discarded at the end by penalization. In this section, we will consider whether we can simplify the approach by introducing merely a density $\theta$. We will use the following ”convexification” approach. A classical shape is parametrized by characteristic functions $\chi(x)$. If we convexify this admissible set, we obtain $\theta(x)\in[0,1]$. Replacing the admissible set by the convexified set, the Hooke’s law, which was $\chi(x)A$, becomes $\theta(x)A$. We will call these $\theta(x)A$ “fictitious materials”, because one can not realize them by a true homogenization process in general. Combined with a penalization scheme, this method is called SIMP method.

We consider the elasticity setting. For $\theta\in L^{\infty}(D;[0,1])$, the convexified formulation for the problem reads as follows:

[TABLE]

where $e(u)\mathrel{\mathop{\mathchar 58\relax}}=(\nabla u+(\nabla u)^{t})/2$. Moreover, the compliance minimization becomes

[TABLE]

with

[TABLE]

There is only one single design parameter, the material density $\theta$. In other words, any information concerning the microstructure $A^{*}$ has disappeared.

5.4.1 Existence of solutions

In this section, we will show the existence of the minimizer of (5.23).

Theorem 5.14**.**

The convexified formulation

[TABLE]

admits at least one solution.

Proof.

Let $\mathcal{M}^{s}_{N}$ be the set of symmetric squared matrices of order $N$. The function

[TABLE]

is convex because

[TABLE]

where the derivative $D\phi$ is given by

[TABLE]

Then, by Theorem 2.7 we can see that there exists a minimizer of (5.23). ∎

5.4.2 Optimality condition

If we exchange the minimizations in $\tau$ and in $\theta$, we can compute the optimal $\theta$ which is

[TABLE]

By using this explicit optimality formula, we can use again an “alternating” double minimization algorithm which will be shown in the following subsection.

5.4.3 Numerical algorithm

By the use of the explicit optimality formula (5.25) for $\theta$, we can apply the following double minimization algorithm.

In practice, it is extremely simple, however, the numerical results are not as good. This can be explained as follows: since the SIMP method uses very little information on the given composites (in particular, it neglects its microstructure all together), we have a lack of a relaxation theorem. In other words, applying the SIMP method changes the problem, and, as a consequence, we have no guarantee that the optimal solution obtained is really an approximation of the optimal solution of the original problem. Moreover, we have to be careful, as it could be very delicate to monitor the penalization. We show the numerical result for the case of the bridge problem (Figure 31). In Figure 33, we show the convergence history of the objective function. Here, the horizontal axis means the iteration number. In Figure 34, we show the numerical result of the optimal shape of the bridge by the algorithm which is mentioned in this subsection.

5.4.4 Concluding remarks on the SIMP method

SIMP (or convexification, or “fictitious materials”) is very simple and very popular. Actually, many commercial codes are using the method. However, since SIMP uses very little information on composites, its simplicity comes at a cost. In particular, contrary to the homogenization method, SIMP is a convexification method and not a relaxation method, i.e., it changes the problem. Hence there is a gap between the true minimal value of the objective function and that of SIMP.

5.5 Generalizations of the homogenization method

As possible generalizations of the homogenization method, we give the following three examples:

• •

Multiple loads,

• •

Vibration eigenfrequency,

• •

General criterion of the least-square type.

The first two cases are self-adjoint and we have a complete understanding and justification of the relaxation process. The third case is, however, not self-adjoint and only a partial relaxation is known. For the detail of these cases, see [Al2002].

5.6 Exercises

Problem 5.6.1**.**

Implement the optimal radiator test case (with penalization).

Problem 5.6.2**.**

Implement the SIMP method (start with exponent $p=1$ and increase to $p=3$) for compliance minimization$\colon$ cantilever, bridge, MBB beam, L-beam.

Problem 5.6.3**.**

Implement SIMP method with an adjoint for the objective function

[TABLE]

Problem 5.6.4**.**

Minimize compliance for the cantilever problem with a parametrized cell (rectangular hole in a square cell).

Chapter 6 Resurrection of the homogenization method: lattice materials in additive manufacturing

6.1 Introduction

As we mentioned in the previous chapters, the homogenization method uses true composite materials, possibly anisotropic, and this makes it complicated to implement. Therefore, in practice, the method was replaced by its much simplified version, the so-called SIMP method, which uses only fictitious isotropic materials. (For the detail of the SIMP method, see [BS2003].) Since intermediate densities (between full material and void) are penalized in the end, there is indeed no need to have a detailed knowledge and optimization of microstructures.

Nevertheless, the recent progress of additive manufacturing techniques has revived the interest for the use of graded or microstructured materials since they are now manufacturable (see Figure 6, page 6). Since the homogenization method is the right technique to deal with microstructured materials, where anisotropy plays a key role (a feature which is absent from SIMP), we could well see a resurrection of the homogenization method for such applications. There is however one final hurdle to overcome, once an optimal composite structure has been obtained, that is the projection of the optimal microstructure at a chosen finite lengthscale to get a global and detailed picture of the optimal microstructure. This is the most delicate part of this homogenization approach and the one where this chapter is most contributing.

Often (but not always) lattice materials are periodic structures, with macroscopically varying parameters. Hence we will restrict to periodic homogenization and macroscopically modulated periodic structures, i.e., the material parameters are of the type

[TABLE]

where $y\mapsto A(x,y)$ is $Y$-periodic and $x\mapsto A(x,y)$ describes the macroscopic variations. (We discuss how to choose the period and the holes in Section 6.2.) The orientation of the microstructures is rarely taken into account and optimized, although it is well-known that their orientation is a crucial and determining parameter in topology optimization (see [Al2002] and [Pe1989]). Actually, even if optimizing the microstructure orientation is not difficult, reconstructing the oriented periodic structure is a challenging issue. We propose a method to settle the difficulty by projecting the optimal microstructure on a fine mesh of the overall structure in a smoothly varying way. This idea was first introduced by [PT2008].

In this chapter, we consider the post-treatment of 2-d compliance minimization, which is based on [AGP2018]. We will improve the pioneer work [PT2008] in several aspects, which we will note below. See also [GS2018] for another homogenization method in the spirit of [PT2008]. Note that this methods can be extended for the case of 3-d compliance minimization as well [GAP2018]. One of the difficulties of this extension is to treat the rotation of the materials. Indeed, compared to the 2-d case, where rotations are parametrized by a single angle (see Section 6.4), the 3-d case is more involved and requires new ingredients.

6.2 Setting of the problem

As in the previous chapters, we will consider compliance minimization problems in $N=2$. Let $D\subset\mathbb{R}^{N}$ be a smooth bounded domain and $\Omega\subset D$ be the reference configuration of a homogeneous isotropic linear elastic body whose Hooke’s law $A$ is defined by (4.32), with Lamé coefficients $\lambda$ and $\mu$. We assume that $\Omega$ is clamped on $\Gamma_{D}\subset\partial\Omega$ and subject to surface loads $g$ on $\Gamma_{N}\subset\partial\Omega$. Also, for simplicity, these parts $\Gamma_{D}$ and $\Gamma_{N}$ of the boundary $\partial\Omega$ are assumed to be fixed and subsets of $\partial D$. The displacement vector field $u$ and the stress tensor $\sigma$ are given by the following system

[TABLE]

where $e(u)\mathrel{\mathop{\mathchar 58\relax}}=(\nabla u+(\nabla u)^{t})/2$ is the strain tensor. Now, let us consider the following compliance minimization problem:

[TABLE]

where $V>0$ is the maximum admissible volume and the objective function $J$ is the compliance

[TABLE]

As we have already seen in Chapter 5, the compliance minimization problem (6.1) does not admit a classical solution. This is why we consider the homogenized problem. We introduce composite structures characterized by the local volume density $\theta(x)$ of the material and a homogenized elasticity tensor $A^{*}(x)$, corresponding to its microstructure. Then, the homogenized or macroscopic displacement $u^{*}$ is the solution of the system

[TABLE]

By the above setting, the relaxed or homogenized optimization problem is obtained as follows:

[TABLE]

where $u^{*}$ is the solution of (6.2) and $P_{\theta}(x)$ is a given subset of effective or homogenized Hooke’s laws for some well-chosen microstructures of density $\theta(x)$.

Our aim in Chapter 6 is to propose a specific subset $P_{\theta}$ of periodic composites and to construct a minimizing sequence for (6.3). A typical example is that of a rectangular hole in a square cell (Figure 21, page 21), where the cell parameters are the lengths $m_{1},m_{2}>0$ and the rotation angle $\alpha$ (acting either on the hole or the whole cell). Also, we let the homogenized tensors be denoted by $A^{*}(m_{1},m_{2},\alpha)$. We note that the same ideas in the following are applicable to other geometries as well.

6.3 A three steps approach

To achieve our purpose, we take a three steps approach for the optimization. The first step is to pre-compute the homogenized properties $A^{*}(m_{1},m_{2},\alpha)$ for all values of the parameters. The second step is to apply a simple parametric optimization process to the homogenized problem. Compared to Chapter 3, we replace the thickness field $h$ by new parameter fields $m_{1}$, $m_{2}$ and $\alpha$ which vary in space. The third step is to choose a length scale $\varepsilon$ and reconstruct a periodic domain $A(x,\frac{x}{\varepsilon})$ approximating the optimal $A^{*}$. We remark that the third step is the trickiest part of the whole process.

The most delicate point is the combined problem of the orientation of the microstructure and the reconstruction of a macroscopically varying periodic lattice. In order to avoid these difficulties, there are two possible approaches. The first one is a “naive” approach, that is, we assume that the periodic grid is never deformed and the holes are simply rotated. The main advantage of the “naive” approach is that the reconstruction of the periodic perforated structure is very easy (see e.g. [GLRS2013]). We remark that this approach is naive because the “skeleton” of the reconstructed structure is fixed and thus it does not follow the supported stresses or forces.

The second one is a deeper approach initiated by Pantz and Trabelsi [PT2008]. The main advantage of this new approach is that the reconstructed structure adapts its geometry to the supported stresses or forces (in some sense, it looks like Michell trusses [ROZ1989] in the case of dimension $2$).

The main difficulty is then to reconstruct such a macroscopically deformed periodic perforated structure. There may be issues on the regularity of the orientation field for stresses or forces. This is the approach that we follow in the sequel.

6.4 1st step: pre-computing the homogenized properties

For the square cell with a rectangular hole, with a fixed orientation, we compute the homogenized properties $A^{*}(m_{1},m_{2})$ for a discrete sampling of $0\leq m_{1},m_{2}\leq 1$. We also compute the derivatives of $A^{*}(m_{1},m_{2})$ with respect to $(m_{1},m_{2})$ by using a shape derivative and an adjoint approach in the same manner in Section 3.3 (see [AGP2018, Section $3.3$ and $3.4$] for the details and [HP2018, SK1992] for an introduction on shape derivatives). We note that we can also compute the derivatives of $A^{*}(m_{1},m_{2})$ with respect to $(m_{1},m_{2})$ numerically by means of a finite difference method.

In the following, for simplicity, we will only consider the case of dimension $N=2$. Assume that the cell is rotated by an angle $\alpha$. We define $Y_{\alpha}(m)$ as the periodic cells with hole $m=(m_{1},m_{2})$, together with the orientation $\alpha$ of the cell. Then the dependency of the homogenized properties $A^{*}(m_{1},m_{2},\alpha)$ with respect to the angle $\alpha$ is given by

[TABLE]

where $R(\alpha)$ is the fourth-order tensor defined by

[TABLE]

and $Q(\alpha)$ is the rotation matrix of angle $\alpha$. Indeed, the variational formulation (4.33) and the characterization of $A^{*}$ (4.34) imply that the quadratic form $A^{*}(m_{1},m_{2},\alpha)\xi\mathrel{\mathop{\mathchar 58\relax}}\xi$ satisfies

[TABLE]

for any $\xi\in\mathcal{M}^{s}_{2}$. On the other hand, we see that, for $\xi\in\mathcal{M}^{s}_{2}$,

[TABLE]

Thus, the numerical computation of the homogenized properties $A^{*}(m_{1},m_{2},\alpha)$ can be restricted to the case $\alpha=0$. Note that (6.4) implies that a rotation of the cell by an angle $\pi$ does not change its Hooke’s law as $R(\pi)=-{\rm Id}$. Hence the optimal orientation can only be defined modulo $\pi$.

We show the numerical results for the entires of the homogenized tensor $A^{*}(m_{1},m_{2},0)$ and their derivatives as functions of $m$ (see Figure 37). When $m=0$, then the homogenized tensor $A^{*}(m_{1},m_{2},0)$ is equal to $A$. On the other hand, if $m$ is close to $(1,1)$, then the homogenized tensor is converging to the null tensor. Moreover, one can see easily check, that the entries of $A^{*}(m_{1},m_{2},0)$ decrease, when $m_{1}$ is fixed and $m_{2}$ is increasing (and vice versa). In other words, the cell is globally weaker when its hole is widening in one direction or the other. However, the sensitivity of the component $A^{*}(m_{1},m_{2},0)_{1111}$ to the parameter $m_{2}$ is greater than the one to the parameter $m_{1}$ (see Figure 37). Figure 37 shows the numerical results in the case $A_{1111}=A_{2222}=24.07$, $A_{1122}=12.96$ and $A_{1212}=11.11$. That is explained by the fact that, along the $y_{1}$ axis, the strength of the cell is mainly ensured by the material in the areas above and below the hole, whose sizes depend on $m_{2}$. As one could expect, the homogenized elasticity tensor is quite smooth with respect to the parameter $m$, so it is amenable to a gradient based optimization method.

6.5 2nd step: parametric optimization of the homogenized problem

Let us recall that the homogenized equation in a box $D$ is

[TABLE]

and the compliance minimization problem is

[TABLE]

Here, the minimum is taken over all functions $m_{1}$, $m_{2}\in L^{\infty}(D;[0,1])$ and $\alpha\in L^{\infty}(D)$. One can rewrite the compliance as the minimum of complementary energy (see Section 2.4) as follows

[TABLE]

where

[TABLE]

This is interesting for algorithmic purposes because we can apply the optimality criteria or alternate minimization algorithm of Chapter 3. We note that the orientation optimization with respect to $\alpha$ is very simple, due to a result of Pedersen [Pe1989]. Pedersen proved that the optimal orientation of an orthotropic cell for a given displacement field is the one where the cell is aligned with the principal eigen-direction of the strain tensor. We can easily show a similar result for stress field in the same way.

We compute the volume fraction of material in a single unit cell as

[TABLE]

also, the total volume of the lattice structure is thus

[TABLE]

To implement a volume constraint, we rely on a Lagrangian algorithm

[TABLE]

where $\ell$ is the Lagrange multiplier associated to the volume constraint ${\rm Vol}=V_{0}$.

Now let us consider a numerical algorithm. To minimize with respect to the microstructure $m$, we use the following algorithm of alternate minimization (or optimality criteria). Moreover, to minimize with respect to the orientation $\alpha$, we could use the same method as for the minimization with respect to the microstructure $m$, but Pedersen’s result [Pe1989] is a better (more efficient) algorithm than the gradient descent method to compute the optimal orientation because it is a global minimization method, proving an optimal orientation at each iteration. However, this method can usually not be generalized to other objective functions.

Let us focus on the technical details of the algorithm. The partial derivative of the Lagrangian $\mathcal{L}$ with respect to the parameter $m_{i}$ ($i=1,2$), is given by

[TABLE]

We remark that the derivative of $A^{*}$ with respect to $m_{i}$ can be obtained by employing the use of shape derivatives. For the details of the computation, see [AGP2018].

Finally, $\mathcal{P}$ is the projection operator onto the interval $[0,1]$. In the process of this projection, we have to update the Lagrange multiplier $\ell$, which is constant in $D$, by a dichotomy process designed to respect the volume constraint.

Note that except when $\sigma$ is proportional to the identity, the optimal orientation angle $\alpha$ is unique up to the addition of a multiple of $\pi$. As shown in Figure 38, this non-uniqueness creates a regularity issue for $\alpha$. For example,

• •

$\alpha$ or $\alpha+\pi$ correspond to the same orientation,

• •

where the material density is close to [math] or $1$, the orientation does not play any role (cf. the corners in Figure 38),

• •

there are real singularities of the orientation, like a fan (cf. the bottom middle in Figure 38),

• •

if the values of $m_{1}$ and $m_{2}$ are exchanged, then the optimal orientation switches from $\alpha$ to $\alpha+\pi/2$, but it does not seem to appear in our results.

These issues create some numerical difficulties, that we will explain in the next section.

6.6 3rd step: reconstruction of an optimal periodic structure

In the previous sections we computed an optimal homogenized design (with an underlying modulated periodic structure). In this section, let us enter a post-processing step. We choose a length scale $\varepsilon$ for this projection step and reconstruct a periodic shape with length scale $\varepsilon$, approximating the optimal one.

6.6.1 Projection in the simple case without varying orientation, $\alpha\equiv 0$

First, we consider the case where there is no varying orientation, that is, $\alpha\equiv 0$. The unit cells (rectangular hole in a square) are defined by

[TABLE]

The domain $D$ is paved with cells $\varepsilon Y(m)$. Since the hole size $m(x)$ is varying in $D$, the periodicity cell is macroscopically modulated and we define a projected lattice shape $\Omega_{\varepsilon}(m)$ as

[TABLE]

where $m_{1}(x),m_{2}(x)$ are functions defined on $D$ with values in $[0,1]$.

Remark 6.1**.**

The values of $m_{1}$ $m_{2}$ are not necessarily constant in each cell of the structure. Hence, the holes in the cellular structure $\Omega_{\varepsilon}(m)$ are not exactly rectangles. But, when $\varepsilon$ goes to [math], the sequence of cellular structures converges to the composite with local Hooke’s law equal to $A^{*}_{0}(m)$.

The cellular structures can be defined using level-sets. We introduce two functions $f^{m}_{\varepsilon,i}$, one for each direction

[TABLE]

and the level-set function

[TABLE]

The final structure $\Omega_{\varepsilon}(m)$ is then defined by

[TABLE]

The construction of a minimizing sequence of shapes is immediate: we just have to update the size $\varepsilon$ in the previous level-set function.

6.6.2 Projection in the general case with orientation, $\alpha\neq 0$

This section is based on the paper of Pantz and Trabelsi [PT2008]. The main idea of Pantz and Trabelsi is to find a map $\varphi=(\varphi_{1},\varphi_{2})$ from $D$ into $\mathbb{R}^{2}$ which distorts the regular square grid in order to orientate each square at the optimal angle $\alpha$ (or $\alpha+\pi$). Once $\varphi$ is found, we can proceed as before.

The final shape, now denoted $\Omega_{\varepsilon}(\varphi,m)$, is still defined by a level-set function:

[TABLE]

with $F^{\varphi,m}_{\varepsilon}=\min(f^{\varphi,m}_{\varepsilon,1},f^{\varphi,m}_{\varepsilon,2})$ and

[TABLE]

Geometrically (in $2$-d), the Jacobian $\nabla\varphi$ should be proportional to the rotation matrix defined by

[TABLE]

One possibility (Pantz and Trabelsi) is to find $\varphi$ which (roughly) minimizes

[TABLE]

However, it is not straightforward because $Q(\alpha)$ is not smooth (the orientation is not coherent a priori). Pantz and Trabelsi proposed a complicated trick to avoid this coherent orientation issue. We also mention that Groen and Sigmund [GS2018] suggested another trick from image processing to obtain a coherent orientation.

6.6.3 The new approach of Allaire-Geoffroy-Pantz

In this section, we propose a new approach, based on the paper of Allaire-Geoffroy-Pantz [AGP2018].

Let us recall that geometrically (in $2$-d), at every point $x\in D$, the Jacobian $\nabla\varphi$ should be proportional to the rotation matrix $Q(\alpha)$ defined by (6.6). Moreover, the proportions of the cell have to be preserved in order to converge to a true square and not simply to a rectangle. For this purpose, we impose $|\nabla\varphi_{1}|=|\nabla\varphi_{2}|=e^{r}$, where $r\in H^{1}(D)$ is a scalar dilation field. Then the Jacobian $\nabla\varphi$ should be

[TABLE]

This equation can be satisfied only if $\alpha$ satisfies the following conformality condition.

Lemma 6.2**.**

Let $\alpha$ be a regular orientation field and $D$ be a simply connected domain. There exists a mapping function $\varphi$ and a dilation field $r$ satisfying $\nabla\varphi=e^{r}Q(\alpha)$ if and only if

[TABLE]

We recall that, for a vector field $u=(u_{1},u_{2})$, its curl is defined as $\text{curl}\,u=\nabla\wedge u=\dfrac{\partial u_{2}}{\partial x_{1}}-\dfrac{\partial u_{1}}{\partial x_{2}}$, where $\wedge$ is the $2$-d cross product of vectors. Of course, $\text{curl}\,\nabla\varphi=\left(\text{curl}\,\nabla\varphi_{1},\text{curl}\,\nabla\varphi_{2}\right)=0$.

Proof of Lemma 6.2.

Since the domain $D\subset\mathbb{R}^{2}$ is simply connected, By Poincaré’s lemma, the map $\varphi$ exists if and only if the right hand side of (6.7) is curl-free. That is,

[TABLE]

Let $a_{1},a_{2}$ be the columns of $Q(\alpha)$. Then

[TABLE]

Since, for fixed $\alpha$, $(a_{1},a_{2})$ is an orthonormal basis of $\mathbb{R}^{2}$, by (6.9) we see that the vector ${\nabla}r$ can be decomposed as

[TABLE]

We compute

[TABLE]

It leads to

[TABLE]

Thus, again by Poincaré’s lemma, the dilation factor $r$ exists if and only if the right hand side in the above is curl-free, which leads to the harmonic condition on $\alpha$. ∎

The following proposition is a very useful property of conformal orientations.

Proposition 6.3**.**

If there exists a map $\varphi=(\varphi_{1},\varphi_{2})$ from $D\subset\mathbb{R}^{2}$ to $\mathbb{R}^{2}$ such that

[TABLE]

then all angles are preserved by the map $\varphi$. In particular, small square cells are deformed into almost square cells locally.

Sketch of the proof.

Let $x=(x_{1},x_{2})$ be the origin of a small square $S$ of side length $\varepsilon>0$ and edges $\varepsilon e_{1}$, $\varepsilon e_{2}$ (in other words, the vertices of $S$ are $x,x+\varepsilon e_{1},x+\varepsilon e_{2}$ and $x+\varepsilon e_{1}+\varepsilon e_{2}$). The map $\varphi$ then transforms $S$ into an “almost” square $\varphi(S)$ of vertices $\varphi(x),\varphi(x+\varepsilon e_{1}),\varphi(x+\varepsilon e_{2})$ and $\varphi(x+\varepsilon e_{1}+\varepsilon e_{2})$. By a Taylor expansion, we see that $\varphi(S)$ is, up to terms of order $\varepsilon^{2}$, equal to a parallelogram of origin $\phi(x)$ and edges $\displaystyle\varepsilon\frac{\partial\varphi}{\partial x_{1}}$ and $\displaystyle\varepsilon\frac{\partial\varphi}{\partial x_{2}}$. Since ${\nabla}\varphi=e^{r}Q(\alpha)$, the two edges are orthogonal, so $\varphi(S)$ is “almost” a square with side length $\varepsilon e^{r}(x)$ (possibly rotated with respect to the initial square $S$).

∎

Nevertheless, in the applications we face a problem. Since $\alpha$ is a stress eigen-direction, it has no reason to be a harmonic function in general. Even worse, $\alpha$ might not even be smooth at some places (for example, at corners or at the junction point of different boundary conditions, but at other places as well). A more profound reason lies in the following observation: both $\alpha$ and $\alpha+\pi$ give rise to the same orientation. By Pedersen’s result, we can show that the rotated Hooke’s law $A^{*}(m_{1},m_{2},\alpha)=R(\alpha)^{T}A^{*}(m_{1},m_{2},0)R(\alpha)$ depends only on the double angle $2\alpha$.

From now on, we are going to be working with the double angle $\beta=2\alpha$, thus, removing the indeterminate additive constant $\pi$. In what follows, we shall regularize the double angle $\beta=2\alpha$ and make it harmonic.

6.6.4 Regularization of the double angle $\beta=2\alpha$

As we mentioned before, the orientation $\alpha$ given by the optimization does not necessarily satisfy the conformality condition $\Delta\alpha=0$. Thus, at each iteration of the algorithm, instead of minimizing locally (by using Pedersen’s result) the following quantity

[TABLE]

let us consider to minimize globally

[TABLE]

for a small parameter $\eta>0$, under the harmonic constraint

[TABLE]

This is a non-linear (and non-quadratic) constrained optimization problem. It turns out that working with the angle $\beta$ is not so easy since the Hooke’s law $A^{*}(m_{1},m_{2},\beta)$ is highly nonlinear in terms of $\beta$. It is however quadratic with respect to the vector $b$, defined by

[TABLE]

where the matrix $S(b)$ is defined by $S(b)=R(\alpha)$. It can be shown that $S(b)$ is affine with respect to $b$ by a careful computation based on Pedersen’s result. Therefore, from now on we shall work with $b$ as the main unknown. We remark that ${\nabla}\beta=b\land{\nabla}b$. Indeed,

[TABLE]

Therefore, the objective function becomes

[TABLE]

which in turn has to be minimized under the harmonic constraint

[TABLE]

We iteratively solve the non-linear problem (6.10) under the minimization constraint (6.11) by a Newton-type approximation with an increment $\delta b$ as follows.

Find a step $\delta b^{n}\in H^{1}(D;\mathbb{R}^{2})$ and a Lagrange multiplier $p^{n+1}\in H_{0}^{1}(D)$ such that, for any $\delta c\in H^{1}(D;\mathbb{R}^{2})$ and $q\in H_{0}^{1}(D)$,

[TABLE]

and

[TABLE]

where $S^{\prime}(\delta c)$ is the directional derivative of $S(b)$ in the direction $\delta c$. Notice that, since $S$ is an affine function, then $S^{\prime}(\delta c)=S(c+\delta c)-S(c)$. At each iteration, we update the vector field $b$ as follows

[TABLE]

We now apply this regularization process together with the alternate minimization algorithm.

The above mentioned algorithm is structured as follows:

Remark 6.4**.**

In numerical practice, starting from the optimal (but not necessarily smooth) orientation, a few tens of iterations of this regularization process are enough.

Remark 6.5**.**

One advantage of this $b$-formulation is that it is insensitive to the $2\pi$-modulo of $\beta$.

As we can see by Figure 40, the regularization occurs mainly in areas where the density is close to [math] or $1$, i.e. where the homogenized material is almost isotropic and the orientation has no significant impact.

Remark 6.6**.**

Unfortunately this process does not always work in general. In other words, it does not always yield a smooth harmonic angle $\beta$. This is because of “true” singularities, i.e. singularities of the orientation that remain and thus do not allow the angle to be harmonic. The vector field is not coherently orientable in these cases, as the vector rotates by an angle of $\pm\pi$ along circles which are enclosing the singularities (see Figure 41. Such cases might be handled by means of other regularization processes (e.g. minimizing a Ginzburg-Landau energy [GE2018]).

6.6.5 Computation of the map $\varphi$

Once a harmonic angle $\alpha=\beta/2$ has been found, one needs to compute $r$ and $\varphi$ such that

[TABLE]

Since the dilation field $r$ satisfies

[TABLE]

one computes $r$ as the minimizer in $H^{1}(D)$ of

[TABLE]

Once $r$ has been computed, a naive idea would be to compute $\varphi$ as a minimizer in $H^{1}(D;\mathbb{R}^{2})$ of

[TABLE]

However, we know that, even if $\beta$ is smooth, $\alpha$ may have jumps of the type $\pm\pi$ and thus $Q(\alpha)$ may have jumps of its sign (recall that $Q(\alpha+\pi)=-Q(\alpha)$).

To compute $\varphi$ there are two possibilities.

Find a coherent orientation of $\alpha$ (i.e. choose between $\alpha$ and $\alpha+\pi$ at every point)$\colon$ this is possible only if there are no singularities (this is the approach of Groen and Sigmund [GS2018]).

• 2.

Leave the angle $\alpha$ as it is and extend $\varphi$ to be defined in an abstract manifold. This is the approach of Allaire–Geoffroy–Pantz [AGP2018] and it works also in the presence of singularities.

6.6.6 An abstract manifold setting

Let us introduce the cover space of $D$.

Definition 6.7**.**

Denote by $T$ a rotation matrix field which is a candidate for being $Q(\alpha)$. Then we define

[TABLE]

where ${\rm SO}(2)$ is the set of rotations in $\mathbb{R}^{2}$.

We note that at every point $x\in D$ the rotation satisfies $T(x)^{2}=Q(\beta)(x)$. If the angle $\alpha$ is globally orientable, then $T(x)=Q(\alpha)(x)$ or $T(x)=-Q(\alpha)(x)$, and that $\mathcal{D}$ is simply the union of two copies of $D$, consisting of the two possible signs of $Q(\alpha)$. We assume the simple case where $\alpha$ could be globally oriented (no singularity) but extend it to the singular cases.

We change our working space from $D$ to $\mathcal{D}$. The map $\varphi(x,T)$ is now defined on the manifold $\mathcal{D}$ by

[TABLE]

and the gradient operator in (6.15) defined by

[TABLE]

where $U$ is an orientable open subset of $D$,

[TABLE]

and $g_{U}$ is one of the charts

[TABLE]

with $T_{U}^{2}=Q(\beta)$ and $T_{U}\in C(U,{\rm SO}(2))$. Moreover, without loss of generality, we can assume the antisymmetric property

[TABLE]

Indeed, if $\varphi$ satisfies ${\nabla}\varphi=e^{r}T$, then the map $(\varphi(x,T)-\varphi(x,-T))/2$ still satisfies (6.15) and is antisymmetric. Thus, if the orientation $\alpha$ satisfies the conformality condition (6.8), the map $\varphi$ can be defined as a minimizer of

[TABLE]

over all maps $\varphi$ in

[TABLE]

In practice we face the problem of making the actual computations on the abstract manifold $\mathcal{D}$. In order to solve this problem we use a new idea, namely, non-conformal finite elements on $D$.

On each triangle $K$ of the mesh we compute one continuous orientation $T_{K}$ such that $T_{K}^{2}=Q(\beta)$. Then, we glue together (with $P_{1}$ discontinuous finite elements) these orientations. We compute

[TABLE]

with $g^{\pm}_{K}=\mathrm{Id}\times(\pm T_{K})$. By the antisymmetry of $\varphi$, we obtain

[TABLE]

Then, we minimize with respect to $\varphi$ in the space of $P_{1}$ discontinuous finite elements.

6.7 A final post-processing/cleaning of the lattice reconstruction

The shapes we obtained in Section 6.6 are not straightforwardly manufacturable. Indeed, there are disconnected components of the lattice structure and/or too thin members that should be removed. A final post-processing is made to cure these defects. Note that there is room for improvement in the process.

Let $h_{\rm min}$ be the minimal manufacturable lengthscale or feature size, meaning the smallest possible width of bars and diameter of holes which can be effectively built. Recall that $\varepsilon$ is our choice of a global size of cells. After deformation, the cell size is $h_{c}(x)=\varepsilon e^{-r(x)}$. Hence the local widths of the bars and holes are respectively given by $(1-m_{i}(x))h_{c}(x)$ and $m_{i}(x)h_{c}(x)$.

In the following, we distinguish two regimes, depending of the local size of the cell $h_{c}(x)$. First, if the cell size is too small, a hole and a bar of minimal width cannot coexist and then we have to choose a completely full or void cell. Hence, if $h_{c}<2h_{\rm min}$, a thresholding is applied separately to each field $m_{i}$: it is assigned the value [math] if $m_{i}<0.5$ and $1$ otherwise.

Second, when $h_{c}\geq 2h_{\rm min}$, our post-processing criterion is satisfied if

[TABLE]

Otherwise, we simply threshold the values of $m_{1}$ and $m_{2}$, according to Figure 46, in order to reach void or full materials. The thresholded $m$ is then denoted by $\tilde{m}$.

Let $O_{\varepsilon}(\varphi,\tilde{m})$ be the shape obtained from $\Omega_{\varepsilon}(\varphi,\tilde{m})$ by filling its closed holes. Numerically, the complement of $\Omega_{\varepsilon}(\varphi,\tilde{m})$ is computed step by step, by evaluating the sign of $F^{\varphi,\tilde{m}}_{\varepsilon}$. If it is positive, the current vertex belongs to the complement $\Omega^{c}_{\varepsilon}(\varphi,\tilde{m})$ and then its neighbors, which are not already visited, are added to the list of vertices which should be tested. Otherwise, the current vertex does not belong to $\Omega^{c}_{\varepsilon}(\varphi,\tilde{m})$.

We will regularize the subset $O_{\varepsilon}(\varphi,\tilde{m})$ in order to remove the disconnected bars or the bars that have one free endpoint. Numerically, we explore all the vertices of the complement as follows. For any given vertex, we check each other vertex not further away than a distance $h_{\rm min}$: if this vertex belongs to the complement too, all vertices between them are added to the complement. In this way, we suppress all disconnected bars and all bars of $O_{\varepsilon}(\varphi,\tilde{m})$ that have one free endpoint, which are not too wide. This new subset is denoted by $\tilde{O}_{\varepsilon}(\varphi,\tilde{m})$.

Finally, the post-processed structure is given by the intersection $\tilde{\Omega}_{\varepsilon}(\varphi,\tilde{m})\mathrel{\mathop{\mathchar 58\relax}}=\Omega_{\varepsilon}(\varphi,\tilde{m})\cap\tilde{O}_{\varepsilon}(\varphi,\tilde{m})$. Several post-processed structures $\tilde{\Omega}_{\varepsilon}(\varphi,\tilde{m})$ for the bridge case are displayed in Figure 47.

6.8 Other numerical examples

In this section, we show some numerical examples of the application of the whole method of this chapter. Here, we mention that the method is not always applicable because some singularities cannot be eliminated, as we will see later.

We will show the numerical results for the optimization of a cantilever (Figure 48), an MBB beam (Figure 49) and an L-beam (Figure 50). For each case, we have represented:

• (a)

the optimal orientation of the periodicity cells before regularization,

• (b)

the optimal orientation of the periodicity cells after regularization,

• (c)

the underlying lattice on which the optimal composite is built, i.e., the projection by $\varphi$,

• (d)-(f)

the sequence of shapes after post-processing for the case of $\varepsilon=0.2$, $0.1$ and $0.05$ respectively.

As shown in the case of L-beam (Figure 50), the singularities which appear in Figure 50-(a) are removed during the regularization step (Figure 50-(b)). This is a necessary condition in order to apply our method. Indeed, as we mentioned, there is a case where the singularity cannot be removed, the so-called electrical mast (see Figure 51). Figure 51-(a) shows that two negative singularities, located inside the domain, cannot neither be removed nor pushed toward the boundary during the regularization step.

Hence, the computed grid shown in Figure 51-(b) is clearly not correctly aligned with the optimal orientation of the cells in the vicinity of the singularities. To overcome this problem, at least two different strategies can be considered. One consists in modifying the regularization process in a way that forces more effectively the singularities to be eliminated (see Figure 51-(c) and (d)). Another approach is to adapt the projection step so that it is able to take singularities into account, and will be discussed in the next section.

6.9 Further issues

As we mentioned in Section 6.6, the problem of removing the singularity might persist. To overcome this problem, some preliminary remedies are proposed in the PhD thesis of Perle Geoffroy [GE2018] by either trying to eliminate them by a Ginzburg-Landau approach or compute a map $\varphi$ with the previous approach and an enriched discontinuous finite element space. In [GE2018], one may find others extensions to cases, such as different objective functions, multiple loads and 3-d problems.

6.10 Exercises

Problem 6.10.1**.**

Consider a compliance minimization problem (for one test case like cantilever, bridge, MBB beam, etc.) for the homogenized formulation. Choose an homogenized tensor corresponding to a non-isotropic microstructure (for example a square cell with a rectangular hole with fixed size). Then minimize the compliance with respect to the sole orientation of the microstructure (using Pedersen result).

Problem 6.10.2**.**

For the same compliance minimization problem as in the previous exercise, fix now the orientation and let the parameters of the microstructure (for example, the lengths $m_{1}$ and $m_{2}$ of the rectangular hole, see Figure 21) become the optimization variables. Then minimize the compliance with respect to these parameters (with fixed orientation). Compare the optimal designs and the attained minimal compliances with the previous exercise. In particular, notice that the absence of orientation optimization yields a self-penalizing effect, namely the obtained designs feature almost no intermediate densities (somehow similarly to the SIMP method).

Problem 6.10.3**.**

Combine the two previous optimization (with respect to the orientation and the size parameters) and minimize the compliance for the test case of the previous exercises.