Abstract

We consider the optimization of the topology and geometry of an elastic structure O⊂R2 subjected to a fixed boundary load, i.e. we aim to minimize a weighted sum of material volume Vol(O), structure perimeter Per(O) and structure compliance Comp(O) (which is the work done by the load). As a first simple and instructive case, this paper treats the situation of an imposed uniform uniaxial tension load in two dimensions. If the weight ε of the perimeter is small, optimal geometries exhibit very fine-scale structure which cannot be resolved by numerical optimization. Instead, we prove how the minimum energy scales in ε, which involves the construction of a family of near-optimal geometries and thus provides qualitative insights. The construction is based on a classical branching procedure with some features unique to compliance minimization. The proof of the energy scaling also requires an ansatz-independent lower bound, which we derive once via a classical convex duality argument (which is restricted to two dimensions and the uniaxial load) and once via a Fourier-based refinement of the Hashin–Shtrikman bounds for the effective elastic moduli of composite materials. We also highlight the close relation to and the differences from shape optimization with a scalar PDE-constraint and a link to the pattern formation observed in intermediate states of type-I superconductors.

1. Problem formulation and discussion

The aim of this paper is to gain some insight into optimal geometries of elastically loaded structures, that is, geometries which achieve an optimal balance between some measure of structural rigidity, material weight and structural complexity. We here consider the case where structural rigidity and complexity are expressed in terms of the elastic compliance and the geometry perimeter, respectively. If a low structural complexity or equivalently a small perimeter is of minor importance, optimal geometries often exhibit very fine microstructure, which is difficult to capture by numerical optimization. However, following the lead by Kohn & Müller [1], one can achieve some understanding by exploring how the energy of optimal geometries scales in the model parameters. Such an analysis involves constructing a family of geometries with a certain scaling and proving that no other construction can do better. As a first instructive case, we focus here on two spatial dimensions and a simple geometry with an imposed uniaxial load.

(a) Introduction to compliance minimization

The elastic compliance of a structure O⊂R2 subjected to a surface load F:Γ→R2 on a part of the boundary Γ⊂∂O is defined as
Comp(O)=12∫ΓF⋅uda
for the equilibrium displacement u:O→R2, which is the minimizer of the free energy [2]
E[u]=∫O12Cϵ(u):ϵ(u)dx−∫ΓF⋅udawithϵ(u)=12(∇u+∇uT),
potentially subjected to some Dirichlet boundary conditions on part of ∂O∖Γ. Here, C denotes the elasticity tensor of the material, and we have abbreviated A:B=tr(ATB) for A,B∈R2×2. For the sake of simplicity, we have neglected any body forces. The compliance measures the mechanical work of the load and thus can be interpreted as the inverse structural rigidity.

Fixing a domain of interest Ω⊂R2 with Γ⊂∂Ω, one now seeks the geometry O⊂Ω which minimizes the compliance under a constraint or a penalty on the volume Vol(O)=∫O1dx, which may be interpreted as the material weight or material costs. To prevent geometries with infinitely fine microstructure and thereby ensure existence of minimizers [3,4], additional regularization is required such as a penalization of the geometry perimeter Per(O)=H1(∂O) (H1 denotes the one-dimensional Hausdorff measure). Note that because we penalize Per(O), we always work with geometries of finite perimeter, and throughout this paper, ∂O shall refer to the so-called essential boundary of O [5]. The perimeter regularization may be regarded as a measure of structural complexity or production costs of the design. Altogether, we are led to the optimization problem of finding the optimal shape O⊂Ω with respect to the objective functional
J[O]=αComp(O)+βVol(O)+εPer(O)
for weighting parameters α,β,ε>0.

Compliance minimization represents by now a classical field in optimal design and has been mathematically analysed and numerically implemented in different variants, using level set formulations [6], phase field approaches [7–9], multiple materials [7,8], design-dependent loads [6,7], nonlinear elasticity [9] or other regularizations such as constraining the topological genus [10]. Alternative objectives to the compliance such as the L2-norm of the internal stresses or the expected excess compliance for a probability distribution of loads are much less understood but also possible [11,12].

(b) Energy scaling law for uniaxial load geometry

In this paper, we restrict ourselves to the simple case of a rectangular domain Ω=[0,ℓ]×[0,L] and a boundary load σ¯n imposed on ∂Ω (n denotes the unit outward normal) which is consistent with a uniform uniaxial stress field σ¯=(000F). In other words, a normal tension of magnitude F is applied at the top and bottom boundary Γ=[0,ℓ]×{0,L}⊂∂Ω (figure 1). Note that this load configuration implies the implicit geometry constraint ∂O⊃Γ. Furthermore, for the sake of simplicity, we shall consider an isotropic material with shear modulus μ and zero Poisson's ratio, so that the compliance is finally given by Compμ,F,ℓ,L(O)=∞ if ∂O⊅Γ and otherwise
Compμ,F,ℓ,L(O)=12∫∂Ω(σ¯n)⋅udawithσ¯=(000F),
in which the equilibrium displacement u minimizes the free energy
Eμ,F,ℓ,L[u]=∫Oμ|ϵ(u)|2dx−∫∂Ω(σ¯n)⋅uda
(the superscripts here refer to the model parameters). As can be readily verified, the objective functional
Jα,β,ε,μ,F,ℓ,L[O]=αCompμ,F,ℓ,L(O)+βVol(O)+εPer(O)
satisfies the scale invariance
Jα,β,ε,μ,F,ℓ,L[LO]=βL2J1,1,ε/βL,1/4,Fα/4μβ,ℓ/L,1[O]
so that it suffices to consider the case α=β=L=1, μ=14 which we assume throughout. Instead of J1,1,ε,1/4,F,ℓ,1, we, in the following, simply write Jε,F,ℓ. We show that the optimal geometry satisfies the following scaling law.

(a) Load geometry considered in this paper with a uniform normal tension F at the top and bottom. The optimal design O is sought inside Ω. (b) Sketch of optimal construction (here with three branching levels), which is composed of several unit cells.

Theorem 1.1 (Optimal energy scaling for uniaxial normal load)

In theorem 1.1, we require |F| to stay bounded away from 1 by some δ, where the choice δ=116 just reduces the level of technicality in the proof of the upper inequality. For |F| close to 1, results for related problems (cf. [13,14] and the later discussion in §1d) suggest a different energy scaling involving a power of (1−|F|). In that case, one requires additional constraints on the smallness of ε relative to (1−|F|). Let us also remark that the factor |F|1/3 in the scaling law is only relevant for |F| near zero, because it is of order 1 otherwise.

The upper bound is proved in §2a by constructing a family of designs with the claimed scaling. Those designs follow a classical branching pattern with slight adaptations to compliance optimization (cf. figure 1). Section 2b then puts forward two proofs of how the lower bound scales in ε, one based on convex duality and the other based on a Fourier argument. Alternative proofs of lower bounds are of interest, because they enrich our collection of methods for the analysis of material microstructure. Furthermore, the convex duality proof (in the provided, very elementary and intuitive version) only holds for |F| smaller than and bounded away from 12 in which case it yields the optimal prefactor cℓ|F|1/3, whereas the Fourier argument seems more generalizable but only provides the optimal prefactor for |F| bounded away from 0 and 1.

In the limit ε=0, minimizing geometries cease to exist, because the functional J0,F,ℓ is not lower semi-continuous along sequences of geometries with weakly converging characteristic functions. However, the infimum of J0,F,ℓ can be expressed as the minimum of the relaxation rel(J0,F,ℓ) (the lower semi-continuous envelope), which is given by J0∗,F,ℓ and is well understood (cf. §2 and [15–17]). In the case of our uniaxial load, the minimizers of the relaxed problem are known to be rank-1 laminates, i.e. near-optimal geometries are almost one-dimensional and consist of very many very thin parallel strands along the loading direction. This is the borderline case between two very different situations:

— In one situation, the boundary load is consistent with a uniform stress field σ¯ with two eigenvalues of opposite sign (e.g. a shear load). In that case, the only optimal geometries are known to be rank-2 laminates aligned with the two orthogonal principal stress directions [18]. A one-dimensional branching as in §2a then no longer suffices for non-zero ε, and the energy scaling law is quite different [19].

— In the other situation, σ¯ has two eigenvalues of equal sign (e.g. a hydrostatic pressure load). In that case, there are many different optimal microstructures such as rank-2 and higher-rank laminates or confocal ellipse constructions. Depending on the shape of Ω, there might even be minimizers of J0,F,ℓ without microstructure, so that minO⊂ΩJε,F,ℓ[O] is independent of ε: a simple example is the unit disc Ω={x∈R2||x|<1} under hydrostatic pressure—the optimal geometry O will be an annulus.

The uniaxial load case treated here is relatively simple and thus serves as a good starting point for the understanding of compliance minimization. Let us finally mention that our restriction to zero Poisson's ratio considerably simplifies the analysis, but is not expected to have a strong influence, because the optimal geometries are composed of truss-like structures for which the lateral contraction is less relevant. Future work will deal with the case of non-zero Poisson's ratio.

(c) Dual reformulation via stress field

It is a classical result in linearized elasticity that the compliance can alternatively be expressed in terms of the equilibrium stress σ instead of the equilibrium displacement u [20]. For our choice of parameters, we have
CompF,ℓ(O)=∫O|σ|2dx=minσ∼∈ΣadO∫O|σ∼|2dx,
where the set ΣadO of statically admissible stress fields is given by divergence-free symmetric tensor fields satisfying the prescribed stress boundary conditions,
ΣadO={σ:Ω→Rsym2×2|divσ=0inΩ,σ=0inΩ∖O,σn=σ¯non∂Ω}.
Indeed, testing the Euler–Lagrange equation for the equilibrium displacement, ∂uEF,ℓ[u]=0, with u itself, we obtain
CompF,ℓ(O)=12∫∂Ω(σ¯n)⋅uda=∫O14|ϵ(u)|2dx=−EF,ℓ[u].
Furthermore, by Fenchel duality,
∫O14|ϵ(u∼)|2+|σ∼|2dx≥∫Oσ∼:ϵ(u∼)dx=∫∂O(σ∼n)⋅u∼da−∫Odivσ∼⋅u∼dx
with equality if and only if σ∼=12ϵ(u∼), i.e. σ∼ is the stress induced by the displacement u∼. Hence, for any u∼ and σ∼∈ΣadO, we have
∫O|σ∼|2dx≥∫∂Ω(σ∼n)⋅u∼da−∫O14|ϵ(u∼)|2dx=−EF,ℓ[u∼],
which implies the desired expression for the compliance. Summarizing, an alternative formulation of the compliance minimization problem reads
minO⊂ΩJε,F,ℓ[O]forJε,F,ℓ[O]=minσ∼∈ΣadO∫O|σ∼|2dx+Vol(O)+εPer(O).

(d) Two related optimal design problems

Compliance minimization can be viewed as PDE-constrained optimization in which the PDE constraint is vector-valued. The corresponding scalar version is the optimization of a heat or electricity conductor for prescribed normal flux f¯⋅n:∂Ω→R with f¯=(0F) on the boundary, i.e.
minimizeJscalε,F,ℓ[O]=DissF,ℓ(O)+Vol(O)+εPer(O),
where the energy dissipation is defined as
DissF,ℓ(O)=12∫∂Ω(f¯⋅n)udxforu=argminu∼:O→R⁡∫O14|∇u∼|2dx−∫∂Ω(f¯⋅n)u∼dx
with u:O→R the equilibrium temperature or electric potential. The corresponding dual reformulation via fluxes f=∇u reads
minO⊂ΩJscalε,F,ℓ[O]forJscalε,F,ℓ[O]=minf∼∈ΦadO∫O|f∼|2dx+Vol(O)+εPer(O),
where the set of admissible fluxes is given by
ΦadO={f:Ω→R2|divf=0inO,f=0inΩ∖O,f⋅n=f¯⋅non∂Ω}.

In the above scalar problem, a pointwise flux constraint ∇u⋅n=f⋅n=f¯⋅n on ∂Ω is imposed. A less restrictive problem formulation is to extend the flux to R2 and just penalize the difference f−f¯ on R2∖Ω,
minO⊂ΩJscal,relε,F,ℓ[O]forJscal,relε,F,ℓ[O]=minf∼∈Φad,relO∫O|f∼|2dx+∫R2∖Ω|f∼−f¯|2dx+Vol(O)+εPer(O)
with
Φad,relO={f:R2→R2|divf=0inR2,f=0onΩ∖O}.
This can also be interpreted as penalizing the H−1/2-norm of f−f¯ on ∂Ω instead of strictly requiring f⋅n=f¯⋅n on ∂Ω [21]. This latter optimization problem turns out to be equivalent to the variational pattern formation problem in intermediate states of type-I superconductors: if an external magnetic field f¯=(0F) is applied to a sample Ω⊂R2 of superconducting material, then the material spontaneously partitions into superconducting and normal regions, and an induced magnetic field f:R2→R2 develops, which is divergence-free and zero in the superconducting region. Denoting by O the normal region, the free energy can be modelled as being composed of a magnetic term and a small interfacial energy between both regions [13],
JSCε,F,ℓ[O]=minf∈Φad,relO∫O∪(R2∖Ω)|f−f¯|2dx+∫Ω∖O|f¯|2−1dx+εPer(O).
The relaxed energy for ε=0 has the minimum JSC,0∗,F,ℓ=−(|F|−1)2ℓ. Using ∫0ℓf2(x1,x2)dx1=Fℓ for any f∈Φad,relO and every cross section x2, we obtain
JSCε,F,ℓ[O]−JSC,0∗,F,ℓ=minf∈Φad,relO∫O|f|2−2f⋅f¯dx+∫Ω|f¯|2−1dx+∫R2∖Ω|f−f¯|2dx+Vol(O)+εPer(O)+(|F|−1)2ℓ=Jscal,relε,F,ℓ[O]−J0∗,F,ℓ.

By a straightforward adaptation of the proofs presented in the next section, the energy scaling theorem 1.1 is also valid for the above presented scalar and relaxed scalar problem. In fact, by choosing the second row of σ∼:Ω→Rsym2×2 as a test flux f∼:Ω→R2, we find
Jε,F,ℓ[O]≥Jscalε,F,ℓ[O]≥Jscal,relε,F,ℓ[O]
so that the upper bound would only have to be proved for the compliance minimization and the lower bound for the type-I superconductor problem. Nevertheless, in the next section, we shall present the proof of the upper as well as the lower bound in the context of compliance minimization, because the convex duality proof of the lower bound is very instructive and intuitive, and the Fourier-type argument for the lower bound will form the basis for the analysis of a compliance minimization scaling law in the more complicated case of an imposed shear load.

Note that the scaling law of the type-I superconductor problem has already been proven in [13] (improved estimates for extreme regimes of small or large applied field in three dimensions are given in [14]). The lower bound is obtained by simply observing that the energy is larger than the energy of a particular micromagnetic pattern formation problem with known energy scaling law. This micromagnetic scaling law was derived in two dimensions in [22], following the classic approach due to Kohn & Müller [1] which we also present for compliance minimization in §2b(i), and in three dimensions in [21], where the argument is based on interpolation inequalities whose derivation is not unlike the Fourier-based proof of the lower bound presented in §2b(ii).

2. Energy scaling of optimal elastic design for uniaxial load

Here, we prove the energy scaling law theorem 1.1 for the load case from figure 1 and the compliance minimization energy
Jε,F,ℓ[O]=minσ∈ΣadO∫O|σ|2dx+Vol(O)+εPer(O).
As previously discussed, for non-zero ε, one pays an excess energy over J0∗,F,ℓ=infO⊂ΩJ0,F,ℓ[O]. By identifying O with all those points in Ω, where the stress field is non-zero, we can express this infimum as
J0∗,F,ℓ=infσ∈ΣadΩ∫Ωg(σ)dxwithg(σ)={0ifσ=0,|σ|2+1else.
The integral is not lower semi-continuous with respect to a weakly converging sequence of stress fields σ so that the infimum is not achieved; however, one may replace the infimum by the minimum of the lower semi-continuous envelope, which is obtained by quasi-convexifying the integrand g and which in [15–17] is shown to be
J0∗,F,ℓ=minσ∈ΣadΩ∫Ωg∼(σ)dxwithg∼(σ)={2(|σ1|+|σ2|−|σ1σ2|)if|σ1|+|σ2|≤1,1+σ12+σ22else
with σ1,σ2 the two eigenvalues of σ.

The minimum value J0∗,F,ℓ=2ℓ|F| for |F|≤1 is achieved by σ=σ¯.

Intuitively, owing to the uniaxial load geometry only one row of the stress field σ contributes to the relaxed energy. For integrands with vector-valued arguments, quasi-convexification is equivalent to convexification, and indeed, we also have J0∗,F,ℓ=minσ∈ΣadΩ∫Ωg∗∗(σ)dx for the convexification g** of g with g**(σ)=2|σ| if |σ|<1 and g**(σ)=1+|σ|2 else. This gives a hint as to why convex duality techniques suffice for the proof of a lower bound in §2b(i).

We next prove the upper bound and then provide two different arguments for the lower bound. In the following, for two expressions A and B, we shall abbreviate by A≲B that there exists an independent constant C for which A≤CB. A≳B is used equivalently to B≲A, and A∼B shall stand for A≲B and B≲A.

(a) Upper bound: a classical branching construction

Close to the centre x2∼12, near-optimal structures must be very coarse to save interfacial energy, whereas the material must be distributed quite evenly along the boundary x2∼0, x2∼1 in order to support the load. Hence, one may expect a stepwise structural refinement from the centre to the boundary via branching (cf. figure 1), similar to constructions known from type-I superconductors [13]. As is most useful for branching constructions, we shall proceed in the standard steps; in addition, we just consider the upper half x2∈[12,1] for the construction, because the lower one will be symmetric.

Step 1. Specify the geometry of a unit cell and find its excess energy and optimal aspect ratio. The typical tree-like branching structure from type-I superconductors or similar problems (figure 2) does not suffice since all branches will bend towards each other under load, producing lots of compliance. However, by introducing a cross-truss in between all branching pairs, the tree turns into a truss structure in which each truss is either compressed or dilated, which costs much less than bending.

(a) Classical branching structure describing the pattern emerging in intermediate states of type-I superconductors. (b) Deformation of the same structure under a tension load. All branches bend, which results in a high compliance.

The unit cell of width w and height l is depicted in figure 3. Note that the provided stress field is divergence-free in the distributional sense to conform with the admissibility condition in §1c. The total normal load acting on its top or bottom is given by |F|w. In order to have uniaxial stress of magnitude 1 in each truss (which from the convexification of g we know to be a preferred value), we choose a=(|F|w/2)tan⁡α and b=|F|w/2cos⁡α, where tan⁡α=w/4l and we will ensure w/l≤1. The excess energy over the relaxed energy 2|F|wl per unit cell for ε=0 can then be estimated as
ΔJcell=Compcell+Volcell+εPercell−2|F|wl∼2(2lcos⁡αb+aw2)+4εlcos⁡α+2εw2−2|F|wl∼2|F|wl(w216+l2)+|F|w38l+4εl2+w216+εw−2|F|wl∼|F|w3l+ε(l+w),
where Volcell and Percell denote the volume and perimeter of a unit cell, respectively, and where the cell compliance Compcell is given by integrating the squared stress |σ|2 over the unit cell. The above expression for ΔJcell is minimized by l=|F|w3/ε. We thus choose the dimensions from figure 3 as
l=|F|w3ε,α=arctan(14ε|F|w),a=18|F|wϵ,b=|F|w21+ε16|F|w.
Using the stress field given in figure 3, a detailed but straightforward calculation of the total excess energy of the unit cell reveals
ΔJcell(w)≲|F|w3ε
under the condition w<l.

Sketch of a unit cell; the domains of constant stress are numbered. The right graph serves to indicate the geometrical parameters. The geometry as well as the stresses are symmetric with respect to the middle vertical axis. The 3-zones are such that at the top of the unit cell they match exactly the 3-zones at the bottom of the two attached unit cells on the next finer level. The non-horizontal triangle sides of zones 4 and 5 in the right half of the unit cell are either parallel or normal to (sin⁡αcos⁡α).

Step 2. Compute coarsest unit cell width and estimate total bulk energy. The unit cells will be arranged in different levels as shown in figure 1. Note that the horizontal scaling of a unit cell by a factor of 2 from one level to the next is compatible with the perfect connection of the horizontal top and bottom bases of the cells. At x2=12, the width wcoarse of the unit cells will be largest, and it will be halved with each new layer of unit cells, i.e. the kth layer has unit cell width wk=wcoarse/2k and height lk=|F|wk3/ε=2−3k/2|F|wcoarse3/ε. The total height of Ω equals the accumulated height of the unit cell layers. Hence, using n layers,
1=2∑k=0nlk=2|F|wcoarse3ε∑k=0n2−3k/2∼|F|wcoarse3ε,
which implies
wcoarse∼ε|F|3.
The width of the coarsest unit cell has to be smaller than the width ℓ of Ω from which we obtain the construction constraint
ε≲ℓ3|F|.
The total bulk excess energy can now be estimated from above by
ΔJbulk≤∑i=0nℓwiΔJcell(wi)∼ℓwcoarseΔJcell(wcoarse)∼ℓ|F|1/3ε2/3.

Step 3. Introduce boundary layer. We have to stop layering unit cells as soon as wk reaches the same size as lk. This implies the maximum number of layers n=log2⁡(|F|wcoarse/ε)∼23log2⁡(|F|/ε) from which we obtain the constraint ε≲|F| for our construction. Furthermore, we have to introduce a boundary layer between this nth layer and x2=1. The boundary layer and a corresponding test stress field are shown in figure 4. The top of the finest unit cell layer exactly matches the bottom of the boundary layer from figure 4. Note that wn∼ε/|F|, so that the volume of the boundary layer scales like ℓwn∼ℓε|F|, which is less than ℓ|F|1/3ε2/3 as long as
ε≤|F|4.
The contribution to the compliance is of the same order. Finally, the perimeter contribution is of order ℓε, which is even smaller, so that altogether the boundary layer does not interfere with the bulk energy scaling.

Boundary layer at x2=1 with test stress field. In the definition of the stress σ9 inside each of the periodic circular sectors, (ξ1,ξ2) are the coordinates of the local coordinate system with origin in the sector apex.

(b) Lower bound: a real space and a Fourier approach

We now prove the inequality
cFℓε2/3≤minO⊂ΩJε,F,ℓ[O]−J0∗,F,ℓ,
where the constant cF scales like |F|1/3 for |F| close to zero and is of order 1 for |F| bounded away from zero. The first proof imitates the classical approach due to Kohn & Müller [1] and provides a very intuitive idea of how the deformation of near-optimal structures will be distributed. To avoid technicalities, the proof assumes that |F| is smaller than and bounded away from 12. The second approach to prove the lower bound is based on a refinement of the Hashin–Shtrikman bounds for the effective elastic moduli of composite materials. This proof attains the correct scaling in ε, but not quite the correct expression in |F| for the prefactor in the case that |F| is small. Nevertheless, we present this proof, because it is valid also for |F|≥12 and furthermore forms the basis for the analysis of compliance minimization under shear loads [19]. It is an open question whether its deficiency of the suboptimal scaling in |F| can be remedied by a refinement of the argument.

(i) Lower bound by convex duality

Abbreviate minOJF,ε,ℓ[O] by J^ and set ΔJ=J^−J0F,∗,ℓ=J^−2|F|ℓ. From ΔJ≥εPer(O) for the optimal O, we obtain that there is a generic cross section x2=x^2∼12 on which (up to a constant factor), there are at most ΔJ/ε interfaces. Furthermore, J^≥Vol(O) implies that the material volume fraction on x^2 may be assumed to be less than a constant times J^/ℓ, which by the previous section is ≲|F|. We shall try to compute the excess energy of such a configuration over the relaxed energy J0F,∗,ℓ. Denote by N the number of material gaps on x2=x^2 and denote their union by Γ∼=((a1,b1)∪⋯∪(aN,bN))×{x^2}. Then
JF,ε,ℓ[O]=minσ∈ΣadOσn=0onΓ∼∫O|σ|2dx+Vol(O)+εPer(O)≥infσ∈ΣadΩσn=0onΓ∼∫Ωg(σ)dx.
Note the additional boundary condition σn=0 on Γ∼, which expresses that no normal stress can be transmitted across the material gaps. By convexification, the above is greater than or equal to
minσ∈ΣadΩσn=0onΓ∼∫Ωg∗∗(σ)dx=minσ:Ω→Rsym2×2σn=σ¯non∂Ωσn=0onΓ∼supu:Ω→R2∫Ωg∗∗(σ)+u⋅divσdx=minσ:Ω→Rsym2×2supu:Ω→R2∫Ω∖Γ∼g∗∗(σ)−ϵ(u):σdx+∫∂Ω(σ¯n)⋅uda≥supu:Ω→R2∫Ω∖Γ∼minσ∈Rsym2×2g∗∗(σ)−ϵ(u):σdx+∫∂Ω(σ¯n)⋅uda=supu:Ω→R2∫Ω∖Γ∼−g∗(ϵ(u))dx+∫∂Ω(σ¯n)⋅uda,
where the Lagrange multiplier u∈H1(Ω∖Γ∼;R2) has the interpretation of a displacement and the Legendre–Fenchel dual g* to g is given by
g∗(ϵ)=(g∗∗)∗(ϵ)={0if|ϵ|≤2,|ϵ|24−1else.
Using this formula, the excess energy ΔJ is bounded below by
supu:Ω→R2|ϵ(u)|≤2onΩ∖Γ∼∫∂Ω(σ¯n)⋅uda−2|F|ℓ.

We now construct an appropriate test displacement u. Without loss of generality, assume the case of σ¯ being a tension load (for σ¯ a compression, the direction of the displacement is just reversed). The test displacement u is chosen to pull the fissures Γ∼ open (cf. figure 5): take the ansatz
u=(02ηx2)+2f(x1)(0sgn(x2−x^2)),Du=(00±2f′(x1)2η),ϵ(u)=(0±f′(x1)±f′(x1)2η)
for some η∈R and f:R→R so that |ϵ(u)|2=4η2+2f′(x1)2. We choose
f(x)=2(1−η2){x−aiifx∈[ai,ai+bi2],bi−xifx∈[ai+bi2,bi],0else.
Introducing di=bi−ai and d=(1/N)∑i=1Ndi, we directly obtain
∫∂Ω(σ¯n)⋅uda−2|F|ℓ=2|F|[ℓ(η−1)+1−η22∑i=1Ndi2]≥2|F|[ℓ(η−1)+1−η22Nd2]≥2|F|ℓ(1+(Nd)42N2ℓ2−1),
where the last step comes from maximization over η, yielding η2=1/(1+(Nd)4/2N2ℓ2). Using Nd≳(ℓ−J^)∼ℓ (here we exploit that |F| lies well below 12) as well as N≲ΔJ/ε≲ℓ|F|1/3ε−1/3 (from the previous section), we arrive at the bound
ΔJ≳2|F|ℓ(1+(ε|F|)2/3−1)≳ℓ|F|1/3ε2/3,
where we used that ε<|F|.

Illustration of the domain Ω with the material gaps Γ∼ and of the test displacement, a uniform vertical dilation by the factor 2(1+η) superimposed with a piecewise linear deformation pulling the fissures Γ∼ open.

(ii) Lower bound via refined Hashin–Shtrikman bounds

The lower Hashin–Shtrikman bound [23] provides lower estimates on the elastic moduli of composite materials. In particular, given an applied stress field, the Hashin–Shtrikman bound tells us the least compliance with respect to that stress that a microstructured material with prescribed material fraction can have. A neat derivation in discrete Fourier space is given in [18]. We refine that calculation for our setting in order to extract

(i) the energy costs associated with a deviation of the material trusses from the vertical direction and

(ii) the costs associated with a deviation from the optimal material volume.

These two pieces of information will be combined with

(iii) a Fourier-estimate that encodes the strict pointwise boundary condition σn=σ¯n for the material stress σ and

(iv) a Fourier-estimate of Per(O)

to yield the desired bound.

Let χ: Ω→{0,1} be the characteristic function of the optimal geometry O and θ=(1/Vol(Ω))∫Ωχ(x)dx be its volume fraction. Furthermore, define γ:R2→R by γ=χ−θ on Ω, extended by zero outside Ω. For k∈R2∖{0}, abbreviate k¯=k/|k| and denote the Fourier transform of a function f:R2→R by
f¯(k)=∫R2f(x)e−2πik⋅xdx
(the inverse Fourier transform is given by gˇ(x)=∫R2g(k)e2πik⋅xdk). Finally, let Σ0ad denote the symmetric divergence-free tensor fields on R2 vanishing outside Ω,
Σad0={η:R2→Rsym2×2|divη=0inR2,η=0inR2∖Ω}.
Note that those tensor fields have a zero normal component on ∂Ω, ηn=0. Repeating and adapting the calculation in [18], the refined Hashin–Shtrikman bound is obtained as follows:
Comp(O)+Vol(O)=minη∈Σad0(σ¯+η)(1−χ)=0onΩ∫Ω|σ¯+η|2+χdx≥limsupK→0⁡minη∈Σad0∫Ωχ|σ¯+η|2+χ+(1−χ)K−1|σ¯+η|2dx=limsupK→0⁡minη∈Σad0∫Ω|σ¯+η|2+χ+(1−χ)(K−1−1)|σ¯+η|2dx=limsupK→0⁡minη∈Σad0∫Ω|σ¯+η|2+χ+(1−χ)maxτ[2(σ¯+η):τ−(K−1−1)−1|τ|2]dx,
using Fenchel duality in the last step. Assuming the strain τ to be spatially constant and exploiting ∫Ωηdx=0, the previous expression is greater than or equal to
minη∈Σad0∫Ω|σ¯+η|2+χ+(1−χ)2(σ¯+η):τdx=Vol(Ω)(|σ¯|2+2(1−θ)σ¯:τ+θ)+minη∈Σad0∫R2|η|2−2γη:τdx≥ℓ(|σ¯|2+2(1−θ)σ¯:τ+θ)+minη^(k)∈Rsym2×2andη^(k)k=0forallk∈R2∫R2|η^|2−2γ^¯η^:τdk
using Parseval's identity (recall γ=χ−θ on Ω and γ=0 outside Ω). Letting superscript ⊥ denote a counterclockwise rotation by π/2, this expression is minimized by η^=γ^(k¯⊥⋅τk¯⊥)k¯⊥⊗k¯⊥ with k¯=k/|k|, yielding
ℓ(|σ¯|2+2(1−θ)σ¯:τ+θ)−∫R2|γ^|2|k¯⊥⋅τk¯⊥|2dk=ℓ(|σ¯|2+2(1−θ)σ¯:τ+θ−θ(1−θ)max(τ12,τ22))+∫R2|γ^|2[max(τ12,τ22)−|k¯⊥⋅τk¯⊥|2]dk,
using ∫R2|γ^|2dx=∫Ω|γ|2dx=ℓθ(1−θ). Here, τ1 and τ2 denote the two singular values of τ. The integral is non-negative, and the remainder is maximized by τ=σ¯/θ, which we shall assume from now on. Given this τ, the integrand equals (F2/θ2)|γ^|2(1−k¯14)=(F2/θ2)|γ^|2k¯22(1+k¯12)≥(F2/θ2)|γ^|2(k22/(k12+k22)), so that altogether the excess compliance and volume are estimated below by
ΔJelast=Comp(O)+Vol(O)−J0∗,F,ℓ≥ℓ(1θF2+θ−2|F|)+F2θ2∫R2|γ^|2k22k12+k22dk=ℓ(|F|−θ)2θ+F2θ2∫R2|γ^|2k22k12+k22dk.
This estimate quantifies the penalty for deviating from the material fraction θ=|F| and for the Fourier coefficients of χ or rather γ lying away from the horizontal axis.

Because from the upper bound we know ΔJelast≲ℓ|F|1/3ε2/3, we directly see |F|1/3ε2/3≳(|F|−θ)2/θ=|F|(|F|/θ−1)(1−θ/|F|), which can be solved for |F|/θ to yield |F|/θ=1+O((ε/|F|)1/3)∼1. Hence, we may reduce the above estimate to
ΔJelast≳∫R2|γ^|2k22k12+k22dk.

The fact that the structure O stops abruptly at x2=0,1 and has to bear a load there enters via the continuous Fourier space analogue of [21], lemma 2.4, which essentially says that the major L2-mass of a (characteristic) function with support in Ω lies beyond the k2-frequency |k2|=1 in Fourier space. For the sake of completeness and readability, let us briefly reproduce the short argument:
∫R|γ^|2k22k12+k22dk2≥∫|k2|>1/4|γ^|2k22k12+k22dk2≥116k12+1(∫R|γ^|2dk2−∫|k2|≤1/4|γ^|2dk2)≥12116k12+1∫R|γ^|2dk2,
where in the last step, we have used |f¯|2≤∥f∥L2([0,1])2=∥f¯∥L2(R)2 for any f:[0,1]→R by Hölder's inequality so that ∫|k2|≤1/4|γ^|2dk2≤12∫R|γ^|2dk2. By Fubini's theorem, we thus arrive at
ΔJelast≳∫R2|γ^|211+k12dk.

Finally, we estimate the perimeter Per(O) in Fourier space by a technique analogous to the proof of [24], lemma 3, where the authors exploit the simple fact Per(O)≥(1/|s|)∥χ−χ(⋅+s)∥L2(R2)2 for all s∈R2. We perform the same calculation, only we work in continuous instead of discrete Fourier space, use γ instead of χ, and solely consider shifts along the x1-axis:
Per(O)≳1L∫L/22L1s∥γ−γ(⋅+s,⋅)∥L2(R2)2ds∼1L2∫L/22L∫R2|γ^(k)(1−e2πisk1)|2dkds≥1L2∫L|k1|≥1|γ^|2∫L/22L|1−e2πisk1|2dsdk≳1L∫L|k1|≥1|γ^|2dk
for any (small) length scale L>0, using Parseval's and Fubini's theorem.

Combining the estimates for ΔJelast and Per(O) and assuming L≲1,
ΔJelastL2+LPer(O)∼ΔJelast(1+1L2)+LPer(O)≳(1+1L2)∫L|k1|≤111+k12|γ^|2dk+∫L|k1|≥1|γ^|2dk≥∫R2|γ^|2dk=ℓθ(1−θ).
The choice L=ε1/3 now yields
ΔJelast+εPer(O)≳ℓθ(1−θ)ε2/3≳ℓ|F|ε2/3.
For |F| bounded away from 0, this bound behaves like ℓε2/3, which was to be shown. For |F| close to zero, however, note that the obtained scaling in |F| is suboptimal.

3. Discussion

Let us here briefly comment on the differences between the structural design task treated here and related problems.

Compared with the geometry optimization problem under a scalar PDE constraint, minO⊂ΩJscalε,F,ℓ, the energy scaling is the same (cf. §1d). Furthermore, as explained in §1d, the two rows of a stress field can be interpreted as two decoupled fluxes of a scalar quantity so that any test geometry and associated test stress for Jε,F,ℓ yield a test geometry and test flux for Jscalε,F,ℓ. However, the opposite is not true, so that both optimization problems are not equivalent. Indeed, a geometry with optimal energy scaling for Jscalε,F,ℓ is typically not sufficient for Jε,F,ℓ. The reason lies in the fact that elastic deformations satisfy two balance equations instead of a single one: the conservation of linear momentum, which expresses that the stress field is divergence-free (a condition also present in the scalar problem), as well as the conservation of angular momentum, which can be expressed as a symmetry constraint on the stress and is special to the elasticity setting. The symmetry of the stress tensor implies that its rows are actually not two independent fluxes. Choosing a near-optimal flux for Jscalε,F,ℓ as one row of the stress field thus induces a constraint for the other row, which typically produces very unfavourable energy. Mechanically, this additional conservation of angular momentum manifests as energy-costly bending of the geometry (cf. figure 2). As a remedy, cross-trusses have to be introduced as in figure 1.

In the case of an applied shear load σ¯=(0FF0), the energy scaling is different (for details, we refer to [19]). A shear load can be interpreted as a superposition of a uniaxial tension and an orthogonal uniaxial compression. Hence, in a near-optimal geometry, we expect to see branching trees similar to the construction from §2a, only in two orthogonal directions. It turns out that the branching trees occur on two different scales: along one direction, they extend over the whole domain (just as for the uniaxial load), along the other direction they only connect the different larger branching trees. In [19], the proof of the lower bound is based on a refinement of the Hashin–Shtrikman bounds analogous to §2b(ii). The ingredients (i)–(iv) are essentially the same, only they have to be adapted and augmented by two further estimates: estimate (i) now has to express the cost associated with a deviation of the geometry from the two preferred orthogonal directions, and the two additional estimates penalize an unbalanced material distribution and an unbalanced spatial distribution between the two directions.

Finally, let us note that slight variations of the load geometry will not alter the energy scaling. For instance, the proofs of lower and upper bound still apply if at x1=0 and x1=ℓ we consider periodic instead of homogeneous Neumann boundary conditions. In addition, if the upper and lower boundary are slightly perturbed normally into a smooth but, for example, no longer straight curve, the proof of the lower bound remains valid with the obvious modifications, whereas the construction has to be adapted a little as indicated in figure 6: unit cells are shifted up or down and connected to unit cells on the next hierarchy level via vertical struts. As a final note on the corresponding three-dimensional problem, in analogy to the pattern formation in type-I superconductors [14], we expect a richer behaviour with a number of different optimal geometries for different regimes. Future work will deal with this and other three-dimensional cases.