A posteriori error estimates for sequential laminates in shape optimization

A posteriori error estimates are derived in the context of two-dimensional structural elastic shape optimization under the compliance objective. It is known that the optimal shape features are microstructures that can be constructed using sequential lamination. The descriptive parameters explicitly depend on the stress. To derive error estimates the dual weighted residual approach for control problems in PDE constrained optimization is employed, involving the elastic solution and the microstructure parameters. Rigorous estimation of interpolation errors ensures robustness of the estimates while local approximations are used to obtain fully practical error indicators. Numerical results show sharply resolved interfaces between regions of full and intermediate material density.


Introduction
In this paper we will address error estimation for shape optimization problems and study an adaptive finite element method based on an optimal microstructured elastic material composed of sequential laminates. A classical shape optimization problem amounts to distribute a fixed amount of rigid material throughout a domain in order to optimize a particular cost functional when confronted with surface and volume loads. However it is known that the problem is ill-posed if only scale invariant cost functionals are taken into account. As the formation of microstructures comes along with this deficiency most approaches enforce regularity by an additional regularizing cost functional such as the shape perimeter. Alternatively the problem can be relaxed to allow for more general shapes characterized by local effective material properties and densities. They stem from the shape structure on the microscale influencing the effective macroscopic description via homogenization. Optimal microstructures exist but are not unique. A particular flexible microstructure is constructed via sequential lamination. This allows for a convenient algorithmic treatment as a small set of parameters characterizing the local microstructure can be computed explicitly from the local stress. Numerical computations suggest that optimal shapes are characterized by regions with substantially different mircostructure and sharp transitions between these regions. Thus an adaptive mesh strategy based on an a posteriori error control is evidently the appropriate strategy. To this end, the error in the achieved cost has to be controlled. Review of related work. Shape optimization is a classical field in PDE constrained optimization and has been covered extensively in the literature, see e. g. the textbooks [Ben95,All02]. For the link to homogenization theory we refer to [JKO94,BD98,CD99,BD91]. Optimal microstructures where first derived by Hashin in 1962 via the concentric sphere construction for hydrostatic loads [Has62]. The sequential lamination construction dates back to the 1980s [Tar85, MT85, FM86, LC86, GC87,Ave87] and was later used in a practical numerical scheme for topology optimization in [ABFJ97]. In all these cases proofs of optimality rely on the Hashin-Shtrikman bounds on the attainable sets of effective elastic properties [HS63]. Related to the homogenization approach is the so called free material optimization, see e. g. [HKLS10]. Here the coefficients of the elasticity are the degrees of freedom to be optimized. Possibly yielding materials without a physical equivalent, best approximating realizations are searched in a post processing step, sometimes referred to as inverse homogenization. Applying the classical homogenization paradigm one asks for the optimal (geometric) microstructures and the effective material properties such that a macroscopic cost functional is optimized. In [BT10a,BT10b] a displacement approach on the microscale is investigated which consists of affine plus periodic functions on generalized periodic domains. For the numerical optimization of the microscopic geometry a boundary tracking method was used. The incorporation of topological derivatives allowed for the nucleation of small holes and additional remeshing steps ensured mesh quality. For specific geometric microstructures the effective stress response of basic affine strains on the microscopic cell boundary are investigated in [BT12] and used to describe the effective macroscopic behavior. The geometry of the locally periodic microstructures are then optimized. The same methodology was picked up in [CGRS14] while investigating microstructured materials with geometrically simple perforations. The resulting shapes offered compliance values close to optimal ones generated by twofold sequential lamination. Moreover they indicate the need to resolve sharp interfaces between regions with different types of local microstructure. To be able to both capture microscopic effects and macroscopic material behavior while maintaining feasible computational complexity in PDE simulation tailored numerical methods have to be proposed. The heterogeneous multi-scale method (HMM) presented in [EEH03, EE03, EE05, EMZ05] depicts a very general paradigm using independent macroscopic and microscopic schemes. In [Ohl05,HO09] a posteriori error estimates for elliptic homogenization problems involving a fine scale diffusion term were derived. Here the heterogeneous multi-scale method is reformulated as a direct finite element discretization of the two-scale homogenized equation in variational form. The obtained error indicator provides separate terms for errors in the macroscopic and microscopic discretization. Concerning the a posteriori control it is often not desirable to measure errors in classical function space norms but the error for a prescribed cost functionals. From a practical point of view quantities of interest for a real composite work piece like e. g. stresses in critical sections, averaged surface tensions or (mollified) pointwise stresses can be considered as cost functionals and corresponding a posteriori error estimates are derived [PO99,OV00a,OV00b,Vem04]. The approach we pick up in this contribution is the dual weighted residual (DWR) method introduced in [BR97] and applied to optimal control problems for instance in [BKR00]. Here residual type error estimates get weighted using duality methods. In [BET11] a special marking strategy was employed allowing to show quasi-optimality of the associated adaptive finite element method while maintaining sound convergence. Control and state constraints were taken into account in [BV09,VW08,LMV13]. Matrix valued L 1 optimal control in coefficients is studied in [KL13]. Shape optimization was addressed in [KV13] using a tracking type cost functional, a Helmholtz state equation, and a graph representation of the boundary. For this control function higher regularity could be shown and a priori estimates were derived. An adaptive finite element method for shape optimization via boundary tracking, remeshing and perimeter penalization was presented in [MNPV10,MNPV12]. Here a DWR type of approach is used to assess the PDE error and the actual geometric error is estimated differently. Furthermore the DWR method was used in [KDLRG14,Kal13] in the context of one shot methods for fuel ignition problems, the viscous Burgers equation and aerodynamic shape optimization. Optimal design in Navier Stokes flow is considered in [BLUU12]. Here the DWR method for optimal control problems leads to separate terms for the primal, the adjoint and the control residual. Similarly [Wol10] addresses a problem in free material optimization based on the variable thickness sheet model. The paper is organized as follows. In section 2 we will describe in condensed form the necessary background on linearized elasticity, shape optimization, the sequential laminates construction in shape optimization, and its numerical discretization. A posteriori error estimates based on the dual weighted residual approach are derived in section 3. In section 4 we discuss the estimation of weighting term using a priori regularity while section 5 focuses on their numerical approximation. Details concerning the implementation are given in section 6 and numerical results are finally presented in section 7.

Elastic shape optimization and microscopic lamination
We consider an elastic object given as a simple connected domain Ω ⊆ D ⊂ R 2 where D is a circumjacent working domain and we at first suppose that D is filled with elastic material (hard and soft) and we aim at the optimization of the distribution of the two phases. A Dirichlet boundary part Γ D ⊂ ∂D ∩ ∂Ω and a Neumann boundary Γ N ⊂ ∂D ∩ ∂Ω subject to a sufficiently regular surface load g are both supposed to be fixed and not subject to optimization. The remaining boundary ∂Ω \ (Γ D ∪ Γ N ) can be varied. The applied forcing causes stresses inside the material which maintain an inner equilibrium. It can be described by the system of partial differential equations of linearized elasticity. Let Σ denote the set of feasible stresses in accordance to the PDE, i.e.
Here u : D → R 2 denotes the displacement, n the outward pointing normal on Γ N , ε[u] = 1 2 Du + Du the symmetrized strain tensor, and C[q] : D → R 2 4 the elasticity tensor. The elasticity tensor is a forth order tensor with the symmetry properties (C[q]) ijkl = (C[q]) jikl = (C[q]) ijlk = (C[q]) klij and the ellipticity condition C[q] ξ : ξ ≥ c ξ·ξ. Usually it can be characterized locally by a small parameter vector q. We consider here material which behaves isotropic and is described by Lamé parameters λ and µ and the strain stress relation σ = 2µε[u] + λ tr ε[u]1. For later purposes we take into account the weak displacement formulation of (1), i. e. u ∈ H 1 Γ D , the space of L 2 vector valued functions with weak derivatives in L 2 (Ω) and vanishing trace on Γ D , solves the variational problem A classical shape optimization problem now amounts to the task of distributing a hard material whose elastic behavior is described by an elasticity tensor A. To this end, we consider a characteristic function χ ∈ L ∞ (D, {0, 1}) of the hard phase with D χ dx = Θ reflecting the volume constraint. The elasticity tensor on D is then defined as is an elasticity tensor for a soft material filling the remaining space. To assess the performance of a given shape for a prescribed loading a suitable cost functional has to be taken into account. We restrict ourselves to the compliance cost as a global measure of rigidity. For a displacement u = u[χ] being the unique solution of (2) for given χ it is given by ( The cost functional J[u[χ], χ] should thereby be minimized with respect to χ subject to a constraint on the volume of the hard material phase. Ultimately one is interested in the degenerate case B = 0, where there is void outside of Ω = {x ∈ D | χ(x) = 1}. This minimizing problem is ill-posed and on a minimizing sequence the onset of microstructures can be observed. Thus one investigates a relaxed formulation of this shape optimization problem. We refer to [ABFJ97, All02] for a concise summary and references therein for further details. Relaxation leads to the following reformulation: Reading from left to right, for fixed stresses σ on D solving (1), a fixed point x ∈ D and a fixed local density θ(x), an elasticity tensor C * [q] minimizing the elastic energy density is to be found among all effective homogenized elasticity tensors resulting from a mixture of materials A and B.
Additionally the Lagrangian multiplier l takes into account the amount of total material used. It can be shown that it is a decreasing function of l which facilitates an adaptation to enforce a fixed prescribed volume at any time. For the set G B θ of homogenized tensors, also referred to as G-closure, no algebraic characterization is known. However there exist sharp bounds on the elastic energy density, called Hashin-Shtrikman bounds. Theses bounds are attained within the subclass of sequentially laminated materials. The construction starts on a scale 1 by layering rigid material A and soft material B with a certain ratio and along a certain direction. By means of homogenization effective material properties can be computed and used in the next iteration in place of material B. Mechanically this takes place on a substantially larger scale 2 1 . In two dimensions and for the compliance cost it is known that the construction on these two iterations is sufficient, cf. [KL88,AK93b]. It is moreover possible to pass from the weak material B to void [AK93a,ABFJ97]. The sequential laminates microstructure is locally characterized by three parameters (cf. Fig. 1): the inclination of the main lamination direction α(x), the ratio of material spent in the second lamination stage m(x), and the overall local material density θ(x). The second lamination direction is always orthogonal to the first one and the second ratio parameter is 1 − m(x). α m 1 2 Figure 1: Sketch of twofold sequential lamination. Numerically optimized result for an L-shaped domain fixed at the bottom with downwards pointing load at the center region on the right hand side. An adaptively refined grid is shown along with the density lamination parameter θ plotted and the von Mises stress using the color coding .
All three parameters directly depend on the eigenvalues λ 1 , λ 2 and the eigenvectors τ (1) , τ (2) of the local stress tensor σ and can be computed explicitly via the mapping , see [ABFJ97]: Here λ and µ denote the Lamé parameters of the rigid isotropic material A and τ y the first and second component of the two dimensional vector τ (i) . Given these parameters the effective homogenized elasticity tensor C * [q](x) can likewise be computed explicitly: Thereby R denotes the linear mapping which transforms the tensorC given in reference configuration into the appropriate coordinate frame spanned by the eigenvectors of the stress tensor. In its definition Q are 2 × 2 rotation matrices and the Einstein summation convention is used. The bulk modulus is defined as κ = λ + µ. The tensorC is complemented using the symmetry relations and filling the remaining entries with 0. This yields a singular elasticity tensor in the sense that a strain with just off diagonal entries, corresponding to shearing deformations, will lead to a vanishing stress. In the following we will assume that the density θ is always bounded away from 0 by a small threshold. Thus we deal with the hardsoft configuration. For the ease of notation we extend the Neumann boundary to Γ N = ∂D \ Γ D and prescribe homogeneous boundary conditions g = 0 on the added parts. For the discretization we use a finite element ansatz and assume that we are given a triangulation T of D with elements T ∈ T and denote by h the piecewise constant mesh size function. Let V (2) h be the space of piecewise bi-quadratic and continuous vector valued functions on which we ask for a solution u h of the discrete weak problem a(q h ; u h , ϕ h ) = l(ϕ h ) for all ϕ h ∈ V (2) h , or more explicitly Associated to u h are stresses for which lamination parameters are computed via formula (5). For brevity of notation we write q h = (α h , m h , θ h ) := q[σ h ]. Let us emphasize that the parameter function varies from point to point. More details on the numerical realization will be given in section 6.

A posteriori error estimates based on the dual weighted residual approach
In this section we will derive estimates of the error in the compliance objective (3) when using optimal sequential laminated microstructures in shape optimization. Let u be the solution of the continuous problem (2) involving the elasticity tensor C * [q] with optimal lamination parameters q defined by (5) and u h be a discrete solution satisfying (7). We are interested in the a posteriori control of the numerical approximation error in the cost functional. To this end we pick up the dual weighted residual approach [BKR00]. In that respect we will treat the rotation parameter α differently from the other two parameters m and θ, which will from now on be denoted by ϑ := (m, θ). Let us begin with a useful observation for the elastic energy density, namely that it is invariant with respect to the rotation of the microstructure. Intuitively this follows from the fact that the lamination construction yields an orthotropic material which is always aligned with the main stress directions.
Lemma 3.1 (Invariance of the elastic energy density w. r. t. rotations). Let x ∈ D be fixed and C * [α, ϑ] be an effective elasticity tensor corresponding to an optimal rank 2 sequential laminate at x. Let u[α, ϑ] be the solution of (2). Then the local elastic energy density is invariant w. r. t. the rotation parameter α, i. e. Here the eigenvalues λ 1 , λ 2 of σ do not depend on the rotation. In fact, the rotation parameter α is derived from the orientation of the eigenvectors. Thus α only appears via σ in the first term of the right hand side. This term also represents an elastic energy density but this time involving the isotropic constituent A. Because of the isotropy of A altering the alignment of σ, e. g. by rotating to reference configuration, leaves this term unchanged. Altogether this means that the right hand side does not depend on the rotation parameter. This proofs the claim.
We now present the main result.
Theorem 3.2 (Weighted a posteriori error estimate). For a given continuous solution u to (2) with a microstructured material given by the optimal lamination parameter function q and a numerical approximation u h the following error estimate holds for the compliance objective (3): Here R is a higher order remainder term involving higher order derivatives and the local error indicators η T of first order are defined by h andm h ,θ h are chosen arbitrarily. Proof. The first step is to use the Lagrangian L(q; u, p) := J[q; u] + a(q; u, p) − l(p) for decoupled q, u, and p and take into account the PDE constraint, namely that the displacement u has to be a solution of the linearized elasticity system (2). The first order optimality conditions are ρ(q; u)(v) := L ,p (q; u, p)(v) = a(q; u, v) − l(v) , where v is a universal variable for an arbitrary test function from the corresponding space. Equation (8) is just the definition of the elastic solution u. Equation (9) defines the dual solution, which in case of compliance is p = −u. As we assumed u and u h to be solutions of the continuous and the discrete problem, respectively, (q, u, −u) and (q h , u h , −u h ) are stationary points of the Lagrangian. Thus we have Estimating the error can therefore be done by considering the difference in the Lagrangian.
Next we consider a suitable first order expansion of the Lagrangian. Here we will take special care of the rotation parameter α to be able to use Lemma 3.1 above. The Lagrangian thus only depends on u and the controls q = (α, ϑ) because the dual solution coincides with the negative primal solution. Hence, in the following we will skip the dual solution as a parameter. We now study the dependence of the elastic solution u[α] on the rotation parameter α.
As a shortcut notation we write e u (s) = u[α] − u h [α h + se α ], e α = α − α h , and e ϑ = ϑ − ϑ h for the difference between the exact and the discrete displacement and controls. Here u h is given for the continuous spatially varying parameter function α h + se α as the solution to the corresponding problem (7  (1 − s)ds . Next, we expand the lower order term and achieve

By definition we have
. This gives Furthermore, we consider the weak formulation (2) and the compliance cost functional (3). For the primal residual we get: This leads to the postulated residual terms with ρ (u) For the control residual we get using (6) Concerning the last expression we know that the controls m h , θ h being the laminate parameters are bounded. Differentiation of the entries of the effective tensor (6) w. r. t. m and θ again yields bounded expressions within the attainable range of the parameters. The tensorsC ,m andC ,θ are therefore bounded and the whole integrand is bounded in L 1 . We finally obtain the desired residual estimate with This proofs the claim.

Estimation of weights using a priori regularity
For the error estimates of Theorem 3.2 the residual terms ρ T can be computed straightforwardly from the discrete solution u h . The weights ω T , however, depend on the exact (unknown) solution u. Using a priori known regularity those terms can be estimated. For the primal solution classical interpolation estimates exist. For the laminate parameters it is possible to estimate the error by exploiting the explicit relation to the primal solution (5). The following estimates especially ensure robustness of the error estimates for a sufficiently smooth continuous solution u: Corollary 1 (A priori estimates of the weights). Under the assumption that u ∈ W 2,∞ the following estimates for the weighting terms hold: Here, h(T ) is the edge length of the square element T and denotes the smaller or equal inequality up to a constant factor depending solely on the triangulation.
Proof. For the primal error weights we chooseũ h = I h u]] and its eigenvalues and finally inserting them into formulae 5. Then, sup x∈T |m(λ 1 (σ), λ 2 (σ))− m(λ 1 (σ h ), λ 2 (σ h ))| is attained for somex ∈ T with stressesσ = σ(x), σ h = σ h (x). From the Lipschitz continuity of the function m we deduce Furthermore, the eigenvalues are Lipschitz continuous functions of the underlying matrices. In particular the Wielandt-Hoffmann inequality [HW53] for p = 2 yields The elasticity tensor C * represents a bounded linear operator, thus with |||C||| representing the associated operator norm of the tensor C. Passing from the symmetrized strain tensor to the full Jacobian yields Finally, by standard L ∞ estimates [Cia78, Theorem 3.3.7], we achieve Remark 1. The presented argument would not work for a weighting term like α −α h 0,∞,T involving the rotation parameter as it depends on the eigenvectors of the local stress and the eigenvectors do not depend continuously on the matrix. In fact this motivates the elimination of the term involving the rotation parameter in Theorem 3.2.

Numerical treatment
For a practical numerical scheme for the weighting terms ω T that were a priori estimated in section 4 we ask for an approximation based on the numerical solution u h . One approach would be to exploit higher regularity of the solution by computing a second discrete solution using a higher order scheme. However this is usually computationally demanding. To keep the cost for the evaluation of the error estimates low [BR97] suggests to merely interpolate the discrete solution to a higher order polynomial. As we use bi-quadratic finite elements to overcome checkerboard instabilities, cf. [JB98], we therefore construct polynomials of degree four on patches of four neighboring elements from the local degrees of freedom of the discrete elastic displacement u h . In the discrete problem (7) the coefficient functions q and thus C * vary from point to point. However it is known, see [Cia78, Theorem 4.1.6], that the convergence order of the finite element scheme is unaltered if the chosen quadrature order is high enough, which is the case in our computations.
To this end the parameters α, m can easily be evaluated at the quadrature points. The same holds true for the density parameter as well. However it needs to be piecewise constant on each element to avoid checkerboard instabilities, see [JB98]. This is ensured by an element wise averaging of all values obtained on quadrature points. In detail we proceed as follows. In case of the ratio parameter m (and the rotation parameter α which however was eliminated from the error estimate) the interpolated discrete deformation u h can be used again to locally obtain parameters α and m for fixed θ. For the density parameter the piecewise constant approximation can be used to construct a bilinear profile on a patch of four neighboring elements. Let us summarize the approximation of the weights in the following remark: Remark 2 (Approximation of local weights). To approximate the weights for the primal and the control error indicator on every element T and to derive a suitable error indicator we consider a patch S of the four neighboring elements including T and proceed as follows: -For the primal error indicator 5 × 5 degrees of freedom on each element of the patch are used to construct a bi-quartic polynomial on the patch. This interpolation I (4) h u h is then compared to the original approximation leading to ω (u) -For the control error indicator the ratio parameter m is computed from the interpolated elastic solution above and compared to the original value; for the density parameter θ four piecewise constant values are used to construct a bilinear profile on the patch which is then compared to the original values. In explicit

Implementation
For our computations we use an adaptive regular quadrilateral mesh provided by the QuocMesh library 1 . Refinement is done via uniform refinement of cells and handling of constrained hanging nodes. To obtain the optimal sequential lamination microstructure we reimplemented the alternating algorithm suggested in [ABFJ97]. The effective tensors are regularized by setting C 44 = 10 −2 .
To ensure coercivity we restrict m to the range [ε, 1 − ε] and θ to [ε, 1] with ε = 10 −3 . The algorithm is terminated once the change in the compliance cost value is below 10 −7 . As already mentioned we use bi-quadratic elements to overcome checkerboard instabilities. For the numerical integration it turned out to be sufficient to use a Gauss quadrature rule of order 5. The lamination parameters for the interpolated solution cannot be directly computed from the pointwise stress h u h ] as it in turn depends on the unknown parameters. Here we use a Newton method to solve the equation for the rotation parameter α and the eigenvalues λ i of the stress (from which m can be computed). The densityθ is fixed to the original value found at the current point.
For the adaptive scheme we use a Dörfler marking strategy [Dör96] with a fraction of 40% of the total error to be associated with the set of elements to be refined.

Numerical results
As the first example we consider a carrier plate under shearing computed on the unit square., cf. Fig. 2. The volume constraint was set to 33%. First we ran the alternating algorithm until convergence on uniform meshes of level l ∈ {2, . . . , 10}. The obtained values can be used to extrapolate the asymptotic value of the compliance cost by assuming the following dependence on the grid width: Using a least squares fit we found the parameters J * = 1.8399, c = 1.7645, and p = 1.0484. J * is used to replace the exact cost value when assessing error reduction in the following. We ran the discussed adaptive scheme for the carrier plate scenario, see Figures 2 and 3. Next, we take into account a cantilever scenario. Here we consider the domain D = [0, 2] × [0, 1] with Dirichlet boundary condition on the left hand side and a downward pointing force in the middle on the right hand side. The volume constraint is set to 50%, cf.   Figure 3: Difference to extrapolated value J * plotted over the number of elements, once for a uniform refinement (green) and the adaptive refinement strategy (red).
corners, i. e. the displacement in y-direction is enforced to be 0, compare for instance [All02]. Additionally the node in the lower left corner is subject to full homogeneous Dirichlet boundary conditions to ensure a unique solution. The lower boundary in between is loaded in vertical direction. The volume constraint is set to 33%.

Conclusion
In this paper we have applied the dual weighted residual approach to an elastic multi scale shape optimization problem using sequential laminates as an underlying optimal microstructure. The microstructure enters the PDE via effective material properties whose explicit formulae allow to compute sensitivities of the cost functional with respect to the describing parameters of the microstructure. The derived goal-oriented error estimate enables to appropriately steer the adaptive refinement leading to a substantial reduction of the number of degrees of freedom compared to uniform meshes and the same accuracy. In particular the goal oriented error estimation outperforms a residual type estimator solely for the elastic displacement problem. The main ingredient is the error term incorporating the sensitivity of the elastic energy density w. r. t. to the local material density. This is in accordance with the observation that the rigidity of a shape is mostly improved by an optimal distribution of the available material. The adaptive scheme leads to sharply resolved interfaces, separating regions with almost no material from regions with bulk material or an intermediate material density.