Get Complete Project Material File(s) Now! »

## Shape optimization and shape identification problems

We introduce the abstract framework of a shape optimization problem, that is an optimiza-tion problem in which the objective functional depends on the shape of the domain in which a Partial Diﬀerential Equation (PDE) is formulated and on the solution of the PDE itself. We review the most common approaches described in the literature to tackle this class of problems and we focus on gradient-based methods which rely on the computation of the so-called shape gradient, that is the diﬀerential form of the objective functional with respect to the shape. Eventually, a review of the results on a posteriori estimators of the discretization error for this class of problems is presented.

**Abstract framework**

We consider an open domain Ω ⊂ Rd (d ≥ 2) with Lipschitz boundary ∂Ω. Let VΩ be a separable Hilbert space depending on Ω, we define uΩ ∈ VΩ to be the solution of a state equation which is a linear PDE in the domain Ω:

aΩ(uΩ, δu) = FΩ(δu) ∀δu ∈ VΩ (1.1)

where aΩ(•, •) : VΩ × VΩ → R is a continuous bilinear form satisfying the inf-sup condition

inf sup aΩ(v, w) = inf sup aΩ(v, w) > 0

w∈VΩ v∈VΩ v∈VΩ w∈VΩ

and FΩ(•) is a continuous linear form on VΩ, both of them depending on Ω. Under these assumptions, problem (1.1) has a unique solution uΩ.

We introduce a cost functional J(Ω) which depends on the domain Ω. We consider the following shape optimization problem inf J(Ω) (1.2) Ω Uad

where Uad is the set of admissible domains in Rd. Several examples of this class of problems arise from applications. A very classical case of shape optimization problem is the minimal surface problem [116] in which we aim to minimize the (d − 1)-dimensional Hausdorﬀ measure Hd−1 with respect to all variations with compact support: J(Ω) := Hd−1(Ω) , Uad := {Ω ⊂ Rd : K ⊂ ∂Ω}.

Another well-known example [95, 161] of shape optimization problem is represented by the minimiza-tion of the Willmore functional J(Ω) := |H|2dS , Uad := {Ω ⊂ Rd : Hd−1(Ω) = A}

where H is the mean curvature and A ∈ R+ is given. In a more general context, the criterion J depends both on the shape of the domain Ω and on the solution uΩ of the state equation and the objective functional and the set of the admissible shapes read as follows: J(Ω) := j(Ω, uΩ), (1.3)

Uad := {Ω ⊂ Rd s.t. uΩ ∈ VΩ and aΩ(uΩ, δu) = FΩ(δu) , ∀δu ∈ VΩ}. (1.4)

Within this framework, problem (1.2) may be viewed as a PDE-constrained optimization problem, in which we aim to minimize the functional j(Ω, u) under the constraint u = uΩ, that is the minimizer u is solution of the state equation (1.1).

Several industrial problems may be formulated as shape optimization problems in which the objec-tive functional and the set of admissible shapes respectively have the forms (1.3) and (1.4). A non-exhaustive list of applications spans the optimal design of elastic structures [17], the drag minimization in aerodynamics [200] and the identification of physical inhomogeneities in Electrical Impedance To-mography [128]. A vast literature has been developed in this field during the last forty years and for a complete introduction to the subject we refer the reader to [53, 60, 178, 207, 229, 238, 242, 256] and references therein.

### A note on the ill-posedness of shape optimization problems

It is well-known in the literature that there may not exist an optimal solution in shape optimiza-tion problems if additional constraints on the nature of the shape are not imposed [88, 106, 162]. For example, in a wide number of shape optimization problems porous structures are known to be more eﬃcient than bulky ones [17] but this class of domains does not belong to the set of admissible shapes Uad. From a practical point of view, this theoretical issue is responsible for many limitations experi-enced by the strategies proposed in the literature to numerically solve shape optimization problems. Among them, we highlight the high sensitivity of the algorithms to the initialization of the guessed shape – and the consequent risk of getting trapped into local minima – and their mesh dependency. This latter aspect is especially critical since the finer the computational mesh under analysis is, the closer to a porous shape whose micro-structure is of infinitesimal size the optimization procedure is allowed to go.

A first strategy to circumvent the issue of non-existence of optimal solutions relies on relaxing the original problem (1.2) by enlarging the set of admissible shapes to include homogenized – that is, porous – structures. This approach finds a mathematically-sound justification in the homogenization theory and for more details we refer the interested reader to [16, 208]. An alternative idea is based on the introduction of additional constraints in the set of admissible shapes in order to impose restrictions on the topology of the shape and avoid highly oscillatory phenomena near the boundaries. Here, we briefly recall some possibilities investigated in the literature and we refer to the corresponding works for additional details. In [31], Ambrosio and Buttazzo proved the existence of optimal shapes by adding a constraint on the perimeter of the shape. In [106], Chenais proved that the shape optimization problem (1.2) is well-posed owing the admissible shapes are uniformly Lipschitz. The well-posedness of the shape optimization problem may be alternatively retrieved by adding a constraint on the topology of the shapes, e.g. by setting an upper bound on the number of connected components in the domain. We refer to [105, 266] for additional details on the subject.

For the rest of this thesis, we will neglect the problem of the existence of an optimal shape for the shape optimization problem (1.2). On the contrary, assuming the existence of a local minimizer, (1.2) reduces to min J(Ω). (1.5) Ω Uad

In this thesis, we investigate some numerical strategies to eﬃciently solve (1.5) by tackling the crucial aspect of defining a reliable stopping criterion for the overall optimization procedure.

**Geometrical description of the shape**

A key aspect in shape optimization problems is the eﬃcient description of the geometry of the domain under analysis. In particular, we recall that the aforementioned shape acts both as the domain on which the state equation (1.1) is posed and as the variable in the optimization problem (1.5). There exist two main kinds of framework that one can use to describe the geometry of the shape: an explicit one and an implicit one.

**Explicit treatment of the geometry**

The most intuitive way to describe the shape of the domain is by means of an explicit treatment of the geometry, e.g. through a polynomial approximation of the boundaries [67], using B´ezier curves or Non-Uniform Rational B-Splines (NURBS) [79, 244] or via a FreeForm Deformation (FFD) map [51]. In all these cases, a set of Degrees of Freedom explicitly provides all the information for the construction of the initial domain and the deformation of the shape during the evolution of the algorithm. When considering a polygonal representation of the domain, its Degrees of Freedom are the physical coordinates of the nodes in the triangulation, whereas for the case of B´ezier curves, NURBS and FFD maps their control points act as Degrees of Freedom for the description of the geometry. We remark that the aforementioned dependency is said to be explicit, since a change in the set of the Degrees of Freedom is automatically responsible for a modification of the corresponding shape.

**Implicit treatment of the geometry**

In their seminal work [217], Osher and Sethian introduced a novel paradigm for the description of the motion of a domain. Basic idea relies on using an implicit representation of the domain Ω ⊂ Rd as the subdomain identified by the negative values of an auxiliary function φ : Rd → R – also known as the level-set function – defined as follows:

φ(x) < 0

φ(x) = 0

φ(x) > 0

, x ∈ Ω

, x ∈ ∂Ω (1.6) , x ∈ Ωc := Rd \ Ω.

Owing to (1.6), the boundary of the design domain is identified by the zero-level contour of the function φ(x) and not by an explicit set of parameters as in the previously described framework. Moreover, we may express local geometric quantities of the shape in terms of the level-set function [283]. Let Ω be of class C2, for any point x on the boundary ∂Ω such that ∇φ(x) = 0 we may write the outward normal vector n(x) to ∂Ω and its mean curvature κ(x) as follows n(x) := ∇φ(x) , κ(x) := div ∇φ(x) . (1.7) | φ(x)| | φ(x)|

#### Classification of shape optimization problems

The goal of shape optimization is to find the shape of a domain that optimizes a given criterion which may depend on the solution of a Boundary Variation Problem as discussed in section 1.1. When dealing with shape optimization problems, three main aspects have to be accounted for: the number of connected regions inside the domain, their shape and their size. Owing to this observation, the following three classes of optimization problems may be identified.

**Parametric optimization**

Also known as size optimization, parametric optimization aims to improve a given objective func-tional by modifying a set of parameters describing the domain under analysis [110]. In particular, the number of connected regions and their shape is set a priori and only minor variations of the existing configuration are allowed [17], e.g. via the modification of the size of the diﬀerent regions. Despite the limited improvements it oﬀers, this optimization strategy is widely used in industrial applications owing to its simplicity and to the possibility of exploiting physical quantities as the size and the thickness of the structural elements as design variables in the problem.

**Geometric optimization**

By performing geometric optimization, one is interested in changing the shape of the domain under analysis to improve the objective functional [17]. Within this framework, the optimization variable of the problem is the boundary of the domain which may be represented either using an explicit or an implicit description as discussed in section 1.2. No a priori restrictions on the shape or the size of the domain is posed but the number of connected regions is not allowed to change. In many published works, geometric optimization is also referred to as shape optimization, whereas the latter term is more accurate to describe the ensemble of the problems discussed in this section.

**Topology optimization**

A key aspect in the engineering community is the identification of the optimal topology (or layout) of the domain under analysis in order to improve a given criterion. This problem is known as topology optimization or layout optimization and aims to identify the best configuration of the domain without imposing any constraint on the number of connected regions inside it [17].

Remark 1.1. Though one may think that topology and geometric optimization problems are extremely diﬀerent, they actually represent complementary aspects of a unique problem. On the one hand, topology optimization focuses on the identification of the optimal number of connected regions inside the domain. On the other hand, geometric optimization is concerned with the optimization of the shape of the existing regions, whose number is a priori set. With an abuse of notation, in the rest of this thesis the term shape optimization will be used to indicate both the geometric optimization paradigm and the ensemble of parametric, geometric and topology optimization. In case of ambiguity, a remark will be provided, e.g. when a given technique requires an additional restriction as the knowledge of the number of connected regions inside the domain.

**A review of the numerical techniques to handle shape optimiza-tion problems**

The interest in the development of techniques to solve shape optimization problems originally arose in the field of engineering and computational mechanics to answer the question of where to place material within a given domain in order to obtain the best structure for a specific target objective. In this section, we provide an overview of the most important methods developed over the years to handle shape optimization problems. For a more detailed comparison of these techniques, we refer to the review paper [253] by Sigmund and Maute and to the works [239] and [140, 269] which focus on more specific aspects, respectively methods that have achieved the stage of application in industrial environments and level-set approaches to structural topology optimization.

The main idea of shape optimization goes back to the seminal work of Hadamard [153] and has been later developed by many authors over several decades [53, 106, 207, 229, 254, 256]. Starting from the original idea of deforming the computational domain according to a given velocity field to improve the value of the objective functional (cf. subsection 1.4.1), we discuss the subsequent developments using the homogenization theory (cf. subsection 1.4.2), the SIMP method (cf. subsection 1.4.3) until the most recent results obtained via the level-set approach (cf. subsection 1.4.4) and the phase-field model (cf. subsection 1.4.5). Eventually, in subsection 1.4.6, we provide a short commentary on the evolutionary strategies that part of the engineering community has been extensively investigating for problems of shape optimization, highlighting the reasons why in our opinion they do not represent a suitable approach for this kind of problems.

**A moving mesh approach**

Given an initial guess for the shape, a first idea to handle shape optimization problems consists in deforming the boundary of the shape along its normal direction. In a more general framework, the deformation may be driven by any smooth velocity field which motivates the name of this method, known in the literature as the velocity approach [118]. From a practical point of view, the most intuitive implementation of the velocity approach is the moving mesh method which exploits the geometrical information of the computational mesh to deform the shape by moving the nodes of the triangulation [83].

The main advantages of this approach are the convenient framework based on the calculus of variations that can be used to handle its theoretical analysis and its straightforward numerical imple-mentation. Nevertheless, the method suﬀers from two main drawbacks. On the one hand, it requires a smooth parameterization of the boundary. On the other one, no topological changes are allowed, that is the number of connected regions inside the domain is a priori set. As pointed out in [16], the latter aspect represents a serious issue when dealing with applications in structural design since it is well-known that for a given weight, porous structures are extremely more eﬃcient than bulky ones.

In order to remedy the issues due to the smoothness properties required for the boundary and the constraint on topological changes, an approach to shape optimization based on the homogenization theory was developed (cf. e.g. [16]). In next subsections, we will discuss the theoretical improvements due to the use of homogenization as a tool to solve shape optimization problems and we will describe the numerical strategies arising from this approach.

#### The homogenization method

The application of homogenization theory to the problem of optimal design of structures was developed by several authors starting from the late 1970’s. In particular, we refer the interested reader to the works of Murat and Tartar [208, 209], Cherkaev, Lurie and co-workers [192, 193, 194] and Kohn and Strang [181, 182, 183]. During the same years, only few works [144, 146, 147] have dealt with the construction of numerical methods based on the homogenization theory and the authors mainly focused on academic problems with limited impact on real-life applications. It is with the seminal paper [61] by Bendsøe and Kikuchi that homogenization methods were used for the first time to tackle real-life problems in the field of structural shape optimization. After this pioneering contribution, several works have proved the eﬃciency of this method in shape optimization [22, 27, 122, 150]. For a complete and rigorous introduction to the subject, we refer the reader to [16]. Starting from the notion of relaxed or homogenized formulation, two main strategies may be pursed to solve shape optimization problems (cf. [16]). On the one hand, the so-called optimality criteria method relies on iteratively using the optimality conditions for the relaxed formulation to update the design parameters. On the other hand, by deriving the objective functional with respect to the design variable, gradient-based strategies may be developed. In both cases, a key aspect of the application of the homogenization method to shape optimization is represented by the choice of a distribution of intermediate densities going from 0 to 1 as design variable. Hence, the geometrical description of the structure as a shape with sharp boundaries – one may visualize it as a black-and-white image in which the white regions are void whereas the black ones feature a given material – is replaced by a graded structure representation, that is a grayscale image. Within this framework, the regions in which the density is small are filled with a soft material mimicking void, whereas a hard material is used where the density tends to 1.

From a practical point of view, the optimal shapes obtained via the homogenization method have the major drawback of not being manufacturable using the machining and molding procedures classically employed in industrial applications. The recent developments in 3D printing technologies may overcome this issue but since the original work [61] of Bendsøe and Kikuchi in 1988, the scientific community has been investigating several techniques to retrieve classical shapes starting from their homogenized optima. A widely used approach relies on the penalization of the intermediate densities to produce 0 − 1 structures [16] and led to the development of SIMP method [60]. Recently, an alternative post-processing technique for the optimal homogenized shapes by means of a projection strategy has been proposed in [221].

**SIMP method: Simplified Isotropic Material with Penalization**

Inspired by the homogenization method, in [59] Bendsøe proposed the strategy known as Simplified Isotropic Material with Penalization (SIMP) or power-law approach. The name is due to the relation-ship posed between the density design variable ρ and the mechanical properties of the material, i.e. the Young’s modulus E: E(ρ) = ρqE0. (1.8)

In (1.8), E0 stands for the Young’s modulus of the solid material and q is the penalization parameter. The SIMP was driven by the will of reducing the complexity of the homogenization method while improving at the same time the convergence of the algorithm towards 0 − 1 shapes. Several other density approaches have been proposed in the literature over the years. Among them, we cite RAMP (Rational Approximation of Material Properties) [264] but we restrict our description to SIMP since all these strategies are based on the idea of providing a continuous interpolation between solid and void through a penalization of the intermediate densities.

Though extremely eﬃcient and easy to implement, SIMP had no rigorous justification when it first appeared in the literature, especially concerning the choice of the penalization parameter q. Prompted by the promising results obtained with the power-law approach, several authors investigated it in the following years. For q = 1, the minimization of the compliance using SIMP corresponds to the well-known Variable Thickness Sheet problem [110] and has been studied in [63, 227]. For q > 1, a penalization of the intermediate densities is introduced: within this framework, small values of q are responsible for grayscale areas whereas big q’s cause the algorithm to quickly converge towards a local minimum. A physical justification of SIMP was provided in [62], where the authors showed that choosing q = 3 the method converges to almost 0 − 1 shapes and the physical realizability of the elements featuring intermediate densities is ensured. Eventually, in [34] the author proved that q = 3 is the exponent that makes density gradients equal to topological derivatives for the case of linear elasticity.

It is well-known in the literature that without smoothing or filtering, SIMP typically converges to local minima with poor performance. Moreover, a major issue of the computed solutions is their mesh dependency, as remarked by the extensive documentation on the checkerboard problem, that is the phenomenon of solutions featuring patches of alternating black and white elements [228]. To circumvent the issue of rapid oscillations of the density distribution, several techniques have been studied over the years. In the rest of this subsection, we briefly recall some ideas to improve SIMP and we refer to [253] for a detailed review of these techniques. A first improvement of SIMP may be obtained by introducing additional constraints on the length of the perimeter or explicit penalizations which depend on the gradient of the density in order to regularize the solution. A possible alternative is represented by sensitivity filters which modify the value of the sensitivity in a mesh element using a weighted average of the values in the neighboring elements within a mesh-independent radius [249]. In [74], Bourdin proposed a density filter which defines the density inside a given mesh element as a weighted average of the design variable in the neighboring elements contained within a mesh-independent radius. More recently, novel techniques have been proposed to overcome the issue of regions featuring grayscale transitions. Starting from the aforementioned filters, in [151] the authors proposed to project the resulting shape in a binary space by means of a smoothed Heaviside function.

Due to its simplicity and the good quality of the resulting shapes, SIMP has encountered a huge success within the engineering community and is currently the basic strategy implemented in most commercial softwares dedicated to topology optimization. Nevertheless, the dependency of the method on multiple parameters and the need of several post-processing techniques to retrieve a 0 − 1 shape starting from the optimum obtained through SIMP limit the robustness of the overall approach. As a matter of fact, a precise tuning of the involved parameters is required to obtain an eﬃcient optimization algorithm. Moreover, the problem-dependent nature of the aforementioned quantities to be estimated makes the resulting strategy extremely heuristic.

**Table of contents :**

Introduction (en francais)

Introduction

The role of shape optimization in the aerospace and aeronautic industry

Non-destructive testing

Optimal design

Thesis outline

**I Optimization of shape-dependent functionals and a posteriori error estimators **

**1 Shape optimization and shape identification problems **

1.1 Abstract framework

1.1.1 A note on the ill-posedness of shape optimization problems

1.2 Geometrical description of the shape

1.3 Classification of shape optimization problems

1.4 A review of the numerical techniques to handle shape optimization problems

1.4.1 A moving mesh approach

1.4.2 The homogenization method

1.4.3 SIMP method: Simplified Isotropic Material with Penalization

1.4.4 A level-set approach

1.4.5 A phase-field model

1.4.6 Gradient-free strategies

1.5 Differentiation with respect to the shape

1.5.1 Hadamard’s boundary variation method

1.5.2 Volumetric and surface expressions of the shape gradient

1.6 A posteriori error estimators for shape optimization

**2 Certification procedure for a genuine descent direction **

2.1 Optimize-then-Discretize and Discretize-then-Optimize

2.1.1 The Boundary Variation Algorithm

2.1.2 The discretized Boundary Variation Algorithm

2.1.3 Certification procedure for a genuine descent direction

2.2 Numerical error in the shape gradient

2.2.1 Bound for the error due to the approximation of a linear functional

2.2.2 Variational formulation of the error in the shape gradient

2.3 Estimate of the error in the shape gradient based on energy-norm error estimates

2.3.1 Complementary energy principle

2.3.2 Explicit residual estimates

2.3.3 Implicit residual estimates

2.3.4 Equilibrated residual estimates

2.4 Improving the estimate of the error in the shape gradient

2.4.1 Dual Weighted Residual method

2.4.2 Energy-norm estimate via the parallelogram identity

2.4.3 Equilibrated fluxes and flux-free approaches

2.5 The Certified Descent Algorithm

2.5.1 CDA based on the complementary energy principle

2.5.2 CDA based on an equilibrated fluxes approach

2.5.3 The mesh adaptation procedure

**II Shape identification: an inverse problem in Electrical Impedance Tomography **

**3 Minimization of the Kohn-Vogelius functional **

3.1 The problem of Electrical Impedance Tomography

3.2 A shape optimization approach

3.2.1 Shape gradient of the Kohn-Vogelius functional

**4 Certified Descent Algorithm based on the complementary energy principle **

4.1 Conforming Finite Element approximation

4.1.1 The state problems

4.1.2 The adjoint problems

4.2 Estimate of the error in the shape gradient via the complementary energy principle

4.2.1 Energy-norm error estimates for the state equations

4.2.2 Energy-norm error estimates for the adjoint equations

4.3 Numerical results

4.3.1 Numerical assessment of the goal-oriented estimator

4.3.2 1-mesh and 2-mesh reconstruction strategies

4.3.3 A more involved test case

4.3.4 The case of multiple boundary measurements

**5 Certified Descent Algorithm based on an equilibrated fluxes approach **

5.1 Estimate of the error in the shape gradient via the equilibrated fluxes

5.2 Equilibrated fluxes for a conforming Finite Element approximation

5.2.1 Equilibrated fluxes for the state equations

5.2.2 Equilibrated fluxes for the adjoint equations

5.2.3 Goal-oriented equilibrated fluxes error estimator

5.3 Discontinuous Galerkin approximation

5.3.1 Weak imposition of the essential boundary conditions

5.3.2 The state problems

5.3.3 The adjoint problems

5.4 Equilibrated fluxes for a Discontinuous Galerkin approximation

5.4.1 Equilibrated fluxes for the state equations

5.4.2 Equilibrated fluxes for the adjoint equations

5.4.3 Goal-oriented equilibrated fluxes error estimator

5.5 Numerical results

5.5.1 Numerical assessment of the goal-oriented estimator

5.5.2 Reconstruction of a single inclusion

5.5.3 The case of two inclusions featuring multiple boundary measurements

Summary of the results and prospective developments

**III Shape optimization: a forward problem in the design of elastic **

**6 Minimization of the compliance under a volume constraint **

6.1 The linear elasticity problem

6.1.1 The pure displacement variational formulation

6.1.2 Mixed variational formulations via the Hellinger-Reissner principle

6.1.3 A dual mixed variational formulation with weakly enforced symmetry of the stress tensor

6.2 Minimization of the compliance under a volume constraint

6.3 Volumetric shape gradient of the compliance via the pure displacement formulation

6.4 Volumetric shape gradient of the compliance via the dual mixed formulation

6.4.1 The case of strongly enforced symmetry of the stress tensor

6.4.2 The case of weakly enforced symmetry of the stress tensor

6.5 Qualitative assessment of the discretized shape gradients via numerical simulations

6.5.1 Experimental analysis of the convergence of the error in the shape gradient

6.5.2 Boundary Variation Algorithm using the pure displacement and the dual mixed formulations

6.5.3 A note on the moving mesh approach

**7 Complementary energy-based estimate of the error in the shape gradient of the compliance**

7.1 Discretization of the pure displacement formulation of the linear elasticity problem

7.1.1 The state problem

7.1.2 The adjoint problem

7.2 Estimate of the error in the shape gradient via the complementary energy principle

7.2.1 Energy-norm error estimates for the state equation

7.2.2 Energy-norm error estimates for the adjoint equation

7.3 A note on the Finite Element approximation of the complementary energy problem

Summary of the results and prospective developments

**Conclusion **

Conclusion (en francais)

A Identification of optimal grasping points for an Unmanned Aerial Vehicle

B Modeling and simulation of a bioreactor landfill

Published and preprint articles

**Bibliography**