Free Material Optimization

Michal Kočvara and Jochem Zowe

Abstract. Free material design deals with the question of finding the stiffest structure with respect to one or more given loads which can be made when both the distribution of material and the material itself can be freely varied. We consider here the general multiple-load situation. After a series of transformation steps we reach a problem formulation for which we can prove existence of a solution; a suitable discretization leads to a semidefinite programming problem for which modern polynomial time algorithms of interior-point type are available. Two numerical examples demonstrates the efficiency of our approach.

1991 Mathematics Subject Classification: 73K40; 90C90; 90C25

Keywords and Phrases: Structural optimization; Material optimization; Topology optimization; Semidefinite programming

1 Problem Formulation

In this section we introduce the problem of free material optimization. Only basic description of the problem is given; for more details the reader is refered to [2, 5, 1]. We study the optimization of the design of an elastic continuum structure that is loaded by multiple independent forces. The material properties at each point are the design variables. We start from the infinite-dimensional problem setting, show the existence of a solution after a reformulation of the problem and, after discretization, reach a finite-dimensional formulation expressed as a semidefinite program, and as such accessible to modern numerical interior-point methods.

First we sketch the single-load model in the two-dimensional space. Let Ω ⊂ ℝ^{2} be a bounded domain (the elastic body) with Lipschitz boundary Γ. The standard notation [H^{1}(Ω)]^{2} and [H_{0}^{1}(Ω)]^{2} for Sobolev spaces of functions v : Ω → ℝ^{2} is used. By u(x) = (u_{1}(x), u_{2}(x)) with u ∈ [H^{1}(Ω)]^{2} we denote the displacement vector at point x of the body under load. Further, let

e_{ij}(u(x)) = \frac{1}{2} (\frac{∂u_{i}(x)}{∂x_{j}} + \frac{∂u_{j}(x)}{∂x_{i}}) for i, j = 1, 2

denote the (small-)strain tensor, and σ_{ij}(x), i, j = 1, 2, the stress tensor. To simplify the notation we will often skip the space variable x in u, e, etc.

Our system is governed by linear Hooke’s law, i.e., the stress is a linear function of the strain

σ_{ij}(x) = E_{ijkl}(x)e_{kl}(u(x)) (in tensor notation),  (1)

where E(x) is the (plain-stress) elasticity tensor of order 4; this tensor characterizes the elastic behaviour of material at point x. The strain and stress tensors are symmetric and also E is symmetric in the following sense:

E_{ijkl} = E_{jikl} = E_{ijlk} = E_{klij} for i, j, k, l = 1, 2.

These symmetries allow us to and interpret the 2-tensors e and σ as vectors

e = (e_{11}, e_{22}, \sqrt{2}e_{12})^{T} ∈ ℝ^{3}, σ = (σ_{11}, σ_{22}, \sqrt{2}σ_{12})^{T} ∈ ℝ^{3}.

Correspondingly, the 4-tensor E can be written as a symmetric 3 × 3 matrix

E = ( \begin{matrix}E_{1111} & E_{1122} & \sqrt{2}E_{1112} \\ & E_{2222} & \sqrt{2}E_{2212} \\ sym. & & 2E_{1212} \end{matrix} ) .  (2)

In this notation, equation (1) reads as σ(x) = E(x)e(u(x)). Henceforth, E will be understood as a matrix and we will use double indices for its elements. To allow switches from material to no-material, we work with E ∈ [L^{∞}(Ω)]^{3×3}.

We consider a partitioning of the boundary Γ into two parts: Γ = \overline{Γ}_{1} ∪ \overline{Γ}_{2}, where Γ_{1} and Γ_{2} are open in Γ and Γ_{1} ∩ Γ_{2} = ∅. Further we put

H = {u ∈ [H^{1}(Ω)]^{2} | u_{i} = 0 on Γ_{1} for i = 1 or 2 or any combination},

i.e., [H_{0}^{1}(Ω)]^{2} ⊂ H ⊂ [H^{1}(Ω)]^{2}. To exclude rigid-body movements, we assume throughout that

{v ∈ H | v_{i} = a_{i} + bx_{i}, a_{i} ∈ ℝ, i = 1, 2, b ∈ ℝ arbitrary} = ∅.

For the elasticity tensor E and a given external load f ∈ [L_{2}(Γ_{2})]^{dim} the potential energy of an elastic body as a function of the displacement u ∈ H is given by

−\frac{1}{2} ∫_{Ω}⟨Ee(u), e(u)⟩ dx + F (u) with F (u) := ∫_{Γ_{2}} f · u dx.  (3)

The system is in equilibrium for u^{∗} which maximizes (3), i.e., u^{∗} which solves

sup_{u∈H} \{−\frac{1}{2} ∫_{Ω}⟨Ee(u), e(u)⟩ dx + F (u)\} .  (4)

Under our assumptions, the supremum in (4) is equal to \frac{1}{2}F (u^{∗}); this value is known as compliance. Now the role of the designer is to choose the material function E such that the “sup” in (4) becomes as small as possible, that is, the body responds with minimal displacements in the direction of the load f. We assume E(x) to be a symmetric and positive semidefinite matrix for almost all x ∈ Ω (recall E ∈ L^{∞}(Ω)), what we write as

E(x) = E(x)^{T} ⪰ 0  a.e. in Ω.  (5)

To introduce a resource (cost) constraint for E, we use the (invariant) trace of E

tr(E(x)) := ∑_{i=1}^{3} E_{ii}(x)  (6)

and require with some given positive α that

∫_{Ω} tr(E(x)) dx ≤ α.  (7)

Further, to exclude singularities at isolated points (e.g., at boundary points of Γ_{2}) we demand that, with some fixed 0 < r^{+} ∈ L^{∞}(Ω),

tr(E(x)) ≤ r^{+}(x)  a.e. on Ω.  (8)

The feasible design functions are collected in a set

E := \{E ∈ [L^{∞}]^{3×3}(Ω) | \begin{matrix}E is of form (2) and \\ satisfies (5), (7) and (8) \end{matrix} \} .  (9)

With this definition, the single-load problem becomes

inf_{E∈E} sup_{u∈H} \{−\frac{1}{2} ∫_{Ω}⟨Ee(u), e(u)⟩ dx + F (u)\} .  (10)

Let us now assume that the structure must withstand a whole collection of independent loads f^{1}, …, f^{L} from L^{2}(Γ_{2}), acting at different times; further, the design should be the “best possible” one in this framework. This leads to the following multiple-load design (MLD) problem, in which we seek the design function E which yields the smallest possible worst-case compliance

inf_{E∈E} sup_{ℓ=1,...,L} sup_{u∈H} \{−\frac{1}{2} ∫_{Ω}⟨Ee(u), e(u)⟩ dx + F^{ℓ}(u)\};  (11)

here

F^{ℓ}(u) := ∫_{Γ_{2}} f^{ℓ} · u dx  for ℓ = 1, …, L.  (12)

2 Existence of a solution

We first eliminate the discrete character of the “ sup_{ℓ=1,...,L}” in (11). With a weight

vector λ for the loads, which runs over the unit simplex

Λ := {λ ∈ ℝ^{L} | ∑_{ℓ=1}^{L} λ_{ℓ} = 1, λ_{ℓ} ≥ 0 for ℓ = 1, …, L} ,

we get from a standard LP-argument as reformulation of (11):

inf_{E∈E} sup_{λ∈Λ} sup_{(u^{1},...,u^{L})∈H×···×H} ∑_{ℓ=1}^{L} \{−\frac{1}{2} ∫_{Ω} λ_{ℓ}⟨Ee(u^{ℓ}), e(u^{ℓ})⟩ dx + λ_{ℓ}F^{ℓ}(u^{ℓ})\} . (13)

The objective function in (13) is linear (thus convex) in the inf-variable E; it is, however, not concave in the sup-argument (λ; u^{1}, …, u^{L}). We will show that a simple change of variable yields a convex-concave version of the problem.

First note that the inf-sup value in (13) remains the same when restricting λ to the half-open set

Λ^{0} := {λ ∈ Λ | λ_{ℓ} > 0 for ℓ = 1, …, L}

and then pass from the variable (λ; u^{1}, …, u^{L}) to

(λ; v^{1} := λ_{1}u^{1}, …, v^{L} := λ_{L}u^{L}).

This converts (13) to

inf_{E∈E} sup_{(v;λ)∈V} ∑_{ℓ=1}^{L} \{−\frac{1}{2} ∫_{Ω} λ_{ℓ}^{−1}⟨Ee(v^{ℓ}), e(v^{ℓ})⟩ dx + F^{ℓ}(v^{ℓ})\} ,  (14)

where we put v := (v^{1}, …, v^{L}) and

V := {(v; λ) | v ∈ [H]^{L}, λ ∈ Λ^{0}} .

The objective function in (14)

F(E; (v; λ)) := ∑_{ℓ=1}^{L} \{−\frac{1}{2} ∫_{Ω} λ_{ℓ}^{−1}⟨Ee(v^{ℓ}), e(v^{ℓ})⟩ dx + F^{ℓ}(v^{ℓ})\}  (15)

is now concave in (v; λ) = (v^{1}, …, v^{L}; λ) ∈ V and a result due to Moreau ([4]) yileds the following existence theorem.

Theorem 1 There exists E^{∗} ∈ E such that

sup_{(v;λ)∈V} F(E^{∗}; (v; λ)) = min_{E∈E} sup_{(v;λ)∈V} F(E; (v; λ)).

Further

inf_{E∈E} sup_{(v;λ)∈V} F(E; (v; λ)) = sup_{(v;λ)∈V} inf_{E∈E} F(E; (v; λ)).

3 Discretization and Semidefinite Reformulation

Using the well-known identity for the trace of the product of a d × d matrix A and the rank-one matrix aa^{T} with a ∈ ℝ^{d}:

tr(A · aa^{T}) = ⟨Aa, a⟩  (16)

we can rewrite the objective function (15) in (14) as

F(E; (v; λ)) = −\frac{1}{2} ∫_{Ω} tr (E · ∑_{ℓ=1}^{L} λ_{ℓ}^{−1}e(v^{ℓ})e(v^{ℓ})^{T} ) dx + ∑_{ℓ=1}^{L} F^{ℓ}(v^{ℓ}).

Due to Theorem 1, we may switch the order of “inf” and “sup” in (14); further, in order to simplify, let us multiply (14) by −2 to get

inf_{(v;λ)∈U} sup_{E∈E} {∫_{Ω} tr (E · ∑_{ℓ=1}^{L} λ_{ℓ}^{−1}e(v^{ℓ})e(v^{ℓ})^{T} ) dx − 2 ∑_{ℓ=1}^{L} F^{ℓ}(v^{ℓ})} .  (17)

For convenience, we will use the same symbols for the “discrete” objects (vectors) as for the “continuum” ones (functions). Assume that Ω is partitioned into M polygonal elements Ω_{m} of volume ω_{m} and let N be the number of nodes (vertices of the elements). We approximate E by a function that is constant on each element Ω_{m}, i.e., E becomes a vector (E_{1}, …, E_{M}) of 3 × 3 matrices E_{m}— the values of E on the elements. The feasible set E is replaced by its discrete counterpart

E := \{E ∈ ℝ^{3×3M} | \begin{matrix}E_{m} = E_{m}^{T} ⪰ 0 and tr(E_{m}) ≤ r_{m}^{+} for m = 1, …, M \\ ∑_{m=1}^{M} tr(E_{m})ω_{m} ≤ α \end{matrix} \} .

To avoid merely technical details we neglect in the following the constraint

tr(E_{m}) ≤ r_{m}^{+} for m = 1, …, M.

Further assume that the displacement vector u^{ℓ} corresponding to the load-case ℓ is approximated by a continuous function that is bi-linear (linear in each coordinate) on every element. Such a function can be written as

u^{ℓ}(x) = ∑_{n=1}^{N} u_{n}^{ℓ}ϑ_{n}(x)

where u_{n}^{ℓ} is the value of u^{ℓ} at n^{th} node and ϑ_{n} is the basis function associated with this node (for details, see [3]). Recall that, at each node, the displacement has 2 components, hence u ∈ ℝ^{D}, D ≤ 2N (D could be less than 2N because of boundary conditions which enforce the displacements of certain nodes to lie in given subspaces of ℝ^{2}).

For basis functions ϑ_{n}, n = 1, …, N, we define matrices

B_{n}(x) = ( \begin{matrix}\frac{∂ϑ_{n}}{∂x_{1}} & 0 \\ 0 & \frac{∂ϑ_{n}}{∂x_{2}} \\ \frac{1}{2} \frac{∂ϑ_{n}}{∂x_{2}} & \frac{1}{2} \frac{∂ϑ_{n}}{∂x_{1}} \end{matrix} ) .

For an element Ω_{m}, let D_{m} be an index set of nodes belonging to this element. The value of the approximate strain tensor e on element Ω_{m} is then (we add the space variable x as a subscript to indicate that e_{x}(u^{ℓ}) is a function of x)

e_{x}(u^{ℓ}) = ∑_{n∈D_{m}} B_{n}(x)u_{n}^{ℓ} on Ω_{m}.

Finally, the linear functional F^{ℓ}(u^{ℓ}) reduces to (f^{ℓ})^{T} u^{ℓ} with some f^{ℓ} ∈ ℝ^{D}. As discrete version of (17) we thus obtain, after a simple manipulation,

inf_{(v;λ)∈V} sup_{E∈E} { ∑_{m=1}^{M} tr (E_{m} · ∑_{ℓ=1}^{L} λ_{ℓ}^{−1} ∫_{Ω_{m}} e_{x}(v^{ℓ})e_{x}(v^{ℓ})^{T} dx) − 2 ∑_{ℓ=1}^{L} F^{ℓ}v^{ℓ}} . (18)

Note that for each element Ω_{m} the d × d matrices ∫_{Ω_{m}} e_{x}(v^{ℓ})e_{x}(v^{ℓ})^{T} dx can be

computed explicitly using the Gaussian integration rule; namely, there exist points x_{ms} ∈ Ω_{m} and weights γ_{ms}^{2} for s = 1, …, S such that

∫_{Ω_{m}} e_{x}(v^{ℓ})e_{x}(v^{ℓ})^{T} dx = ω_{m} ∑_{s=1}^{S} γ_{ms}^{2}e_{x_{ms}}(v^{ℓ})e_{x_{ms}}(v^{ℓ})^{T} .  (19)

For instance, for linear B_{n}(.) (i.e. bilinear ϑ_{n}) one takes S = 4. Hence (18) becomes

inf_{(v;λ)∈V} sup_{E∈E} { ∑_{m=1}^{M} ω_{m}tr(E_{m}A_{m}(v, λ)) − 2 ∑_{ℓ=1}^{L} F^{ℓ}v^{ℓ}}  (20)

where

A_{m} := A_{m}(v; λ) := ∑_{ℓ=1}^{L} λ_{ℓ}^{−1} ∑_{s=1}^{S} γ_{ms}^{2}e_{x_{ms}}(v^{ℓ})e_{x_{ms}}(v^{ℓ})^{T} .  (21)

We now make one further step and introduce a dummy variable ρ_{m} for tr(E_{m}) and m = 1, …, M. Then the constraint E ∈ E in (20) splits into a global part (the global material distribution)

ρ ∈ ℝ_{+}^{M}, ∑_{m=1}^{M} ρ_{m}ω_{m} ≤ α

and a local one (the local material properties)

E_{m} = E_{m}^{T} ⪰ 0, tr(E_{m}) = ρ_{m}, for m = 1, …, M.

The “sup” over the local part can be now put under the sum:

inf_{(v;λ)∈V} sup_{ρ∈ℝ_{+}^{M} ; ∑ ρ_{m}ω_{m}≤α} \{ ∑_{m=1}^{M} ω_{m} sup_{{E_{m}=E_{m}^{T}⪰0 \atop tr(E_{m})=ρ_{m}}} tr(E_{m} · A_{m}(v, λ)) − 2 ∑_{ℓ=1}^{L} F^{ℓ}v^{ℓ}\} . (22)

Now we will analytically perform the inner “sup”, thus finally reaching a semidefinite programming formulation of the multiple-load problem.

Fix m ∈ {1, …, M} and consider the inner “sup” in (22):

sup_{{E_{m}=E_{m}^{T}⪰0 \atop tr(E_{m})=ρ_{m}}} tr(E_{m}A_{m}).  (23)

We use Lagrange theory to write this as

inf_{τ∈ℝ}{τρ_{m} + sup_{E_{m}=E_{m}^{T}⪰0} tr(E_{m}(A_{m} − τI_{d}))}  (24)

with the d × d identity matrix I_{d}. The only τ for which the inenr “sup” is finite are those with A_{m} − τI_{d} ⪰̸ 0. Hence we get for (24)

sup_{{E_{m}=E_{m}^{T}⪰0 \atop tr(E_{m})=ρ_{m}}} tr(E_{m}A_{m}) = ρ_{m} inf_{τI_{d}−A_{m}⪰0} τ.  (25)

With

τ_{m} := inf_{τI_{d}−A_{m}⪰0} τ

our discretized problem (22) becomes (note that A_{m} and thus τ_{m} depends on (v; λ))

inf_{(v;λ)∈V} sup_{ρ∈ℝ_{+}^{M} ; ∑ ρ_{m}ω_{m}≤α} { ∑_{m=1}^{M} ρ_{m}ω_{m}τ_{m} − 2 ∑_{ℓ=1}^{L} F^{ℓ}v^{ℓ}} .

The inner “sup” over ρ is a linear program for each fixed outer variable (v; λ). Hence the “sup” is attained at an extreme point of the feasible ρ−set and we can continue

inf_{(v;λ)∈V}{ max_{m=1,...,M} ατ_{m} − 2 ∑_{ℓ=1}^{L} F^{ℓ}v^{ℓ}},

which in view of (3) is the same as

inf_{{(v;λ)∈V \atop τ∈ℝ}} ατ − 2 ∑_{ℓ=1}^{L} F^{ℓ}v^{ℓ}

s.t.

τI_{d} − A_{m}(v; λ) ⪰ 0 for m = 1, …, M.

(26)

To emphasize the dependence of A_{m} on (v; λ), we have again inserted the variables. With the (d × LS)-matrix

Z_{m} := [γ_{m1}e_{x_{m1}}(v^{1}), …, γ_{ms}e_{x_{ms}}(v^{1}), . . . . . . , γ_{m1}e_{x_{m1}}(v^{L}), …, γ_{ms}e_{x_{ms}}(v^{L})]

and the (LS × LS)-matrix

Λ(λ) := diag(λ_{1}, …, λ_{1}, . . . . . . , λ_{L}, …, λ_{L})

the constraints in (26) become

τI_{d} − Z_{m}(v)Λ(λ)^{−1}Z_{m}(v)^{T} ⪰ 0

which, using a standard result on Shur complement, is equivalent to

( \begin{matrix}τI_{d} & Z_{m}(v) \\ Z_{m}(v)^{T} & Λ(λ) \end{matrix} ) ⪰ 0.

We end up with the announced semidefinite program for the discretization of (14)

inf_{{(v;λ)∈V \atop τ∈ℝ}} ατ − 2 ∑_{ℓ=1}^{L} F^{ℓ}v^{ℓ}

s.t.

( \begin{matrix}τI_{d} & Z_{m}(v) \\ Z_{m}(v)^{T} & Λ(λ) \end{matrix} ) ⪰ 0 for m = 1, …, M.

(27)

The semidefinite program (27) can be efficiently solved by modern interiorpoint polynomial time methods. The question of recovering the optimal elasticity matrices E_{1}^{∗}, …, E_{M}^{∗} from the solution of (27) is a bit technical; again we refer the reader to [1].

4 Examples

Results of two numerical examples are presented in this section. The values of the “density” function ρ are depicted by gradations of grey: full black corresponds to high density, white to zero density (no material), etc. Example 1. We consider a typical example of structural design: The two forces (or force and fixed boundary) are opposite to each other and there is a hole in between because of technological reasons. The geometry of domain Ω and the forces are depicted in Figure 1. The body can be loaded either by the forces on the left or on the right-hand side. Therefore this example has to be considered as MLD (two-load case). Symmetry allows us to compute only one half of the original domain. The resulting values of the “density” function ρ for 37 × 25 mesh are also presented in Figure 1. Again, the figure is composed from two computational domains to get the full body. Example 2. In this example we try to model a wrench. The geometry of domain Ω is depicted in Figure 2. The nut (depicted in full black in Figure 2) is considered to present a rigid obstacle for the wrench. Hence the wrench is in unilateral contact with the nut and there are no other boundary conditions. The loads are also shown in Figure 2. Note that the problem is nonlinear because of the unilateral contact conditions and that for positive vertical force we get a different design than for a negative one; hence we have to consider these two forces as two independent loads. The resulting values of the “density” function ρ for 37 × 22 discretization are shown in Figure 3.

Free Material Optimization 715 Figure 1: Example 1 Figure 2: Example 2 Figure 3: Example 2

Documenta Mathematica · Extra Volume ICM 1998 · III · 707–716

Acknowledgements

The work on this paper was partly supported by GIF-contract I0455-214.06/95, BMBF-project 03ZO7BAY and grant A1075707 of the Czech Academy of Sciences. The first author is on leave from the Czech Academy of Sciences.

References

[1] A. Ben-Tal, M. Kočvara, A. Nemirovski, and J. Zowe. Free material design via semidefinite programming. The multi-load case with contact conditions. Preprint 219, Inst. Appl. Math., Univ. of Erlangen, 1997.

[2] M. P. Bendsøe, J. M. Guades, S. Plaxton, and J. E. Taylor. Optimization of structure and material properties for solids composed of softening material. Int. J. Solids Struct., 33:1179–1813, 1995.

[3] P. G. Ciarlet. The Finite Element Method for Elliptic Problems. North-Holland, Amsterdam, New York, Oxford, 1978.

[4] J.-J. Moreau. Théorèmes “inf-sup”. C. R. Acad. Sc. Paris, t.258, Groupe 1:2720–2722, 1964.

[5] J. Zowe, M. Kočvara, and M. Bendsøe. Free material optimization via mathematical programming. Preprint 213, Inst. Appl. Math., Univ. of Erlangen, 1997.

Michal Kočvara and Jochem Zowe Institute of Applied Mathematics University of Erlangen Martensstr. 3 91058 Erlangen Germany