Structural optimization under overhang constraints imposed by additive manufacturing processes: an overview of some recent results

Grégoire Allaire, Charles Dapogny, Rafael Estevez, Alexis Faure, Georgios Michailidis

► To cite this version:

Grégoire Allaire, Charles Dapogny, Rafael Estevez, Alexis Faure, Georgios Michailidis. Structural optimization under overhang constraints imposed by additive manufacturing processes: an overview of some recent results. Cedya + CMA 2017, Jun 2017, Cartagena, Spain. hal-01589906

HAL Id: hal-01589906

https://hal.science/hal-01589906

Submitted on 21 Sep 2017

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Applied Mathematics and Nonlinear Sciences

http://journals.up4sciences.org

# Structural optimization under overhang constraints imposed by additive manufacturing processes: an overview of some recent results

Grégoire Allaire<fnref n="1" />, Charles Dapogny 2†, Rafael Estevez<fnref n="3" />, Alexis Faure<fnref n="3" /> and Georgios Michailidis<fnref n="3" />.

**Submission Info**

Communicated by Referees

Received Day Month Year

Accepted Day Month Year

Available online Day Month Year

**Abstract**

The purpose of this note is to report on the recent work [2, 3], where new structural optimization strategies are proposed so that the optimized designs are free of *overhang* regions, which jeopardize their constructibility by additive manufacturing technologies. After showing numerical evidence that the intuitive angle-based criteria alone are insufficient to overcome this difficulty, a new constraint functional of the domain is introduced, which aggregates the self-weights of all the intermediate structures appearing in the course of the layer by layer assembly of the total structure. The mathematical analysis of this constraint is outlined and an algorithm is proposed to accelerate the significant computational effort entailed by the implementation of these ideas. Eventually, a numerical validation and several concrete examples are discussed.

**Keywords:** Shape optimization, additive manufacturing, level set method, shape derivative.

**AMS 2010 codes:** 65K10 - 93B40 - 49M30 - 74P05 - 74P20

# 1 Introduction

The additive manufacturing technologies have recently demonstrated a unique potential in constructing structures with a high degree of complexity, thereby allowing to process almost directly the designs predicted by topology optimization algorithms. These breakthroughs however come along with new challenges. One of them is to overcome the difficulty of building shapes showing large *overhangs*, i.e. regions hanging over void without sufficient support from the lower structure.

To give a hint of why such regions are problematic from the point of view of additive manufacturing, a few generalities about these technologies are in order; see [10] or the introduction in [6] for further details. Additive manufacturing is a common label for a whole range of processes, which all begin by a decomposition (or *slicing*) of the produced shape into a series of two-dimensional layers; these are then constructed successively, one on top of the other; see Figure 1.

<footnote n="1">¹ Centre de Mathématiques Appliquées, Ecole Polytechnique, CNRS, Université Paris-Saclay, 91128 Palaiseau, France</footnote>

<footnote n="2">² Laboratoire Jean Kuntzmann, Université Grenoble-Alpes, CNRS, Grenoble INP, F-38000, France</footnote>

<footnote n="3">³ SIMaP, Université Grenoble-Alpes, CNRS, Grenoble INP, F-38000, France</footnote>

<footnote>† Corresponding author.
Email address: charles.dapogny@univ-grenoble-alpes.fr</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Sketch of the slicing procedure, at the beginning of all additive manufacturing processes.</figcaption></figure>

The construction of each individual layer may be achieved owing to quite different technologies; for instance:

* *Material extrusion methods*, such as the well-known *Fused deposition modeling* (FDM) technique, proceed by extruding a molten filament from a nozzle, which is deposed under the form of rasters;

* *Powder bed fusion methods*, such as Electron Beam Melting (EBM) or Selective Laster Melting (SLM), rely on a laser to bind the grains of a (metallic) powder together.

The physical origin of the problems caused by overhangs during the construction process depends on the particular technology: in the case of material extrusion methods, assembling a large overhang feature implies that the machine tool has to depose material over void; in the case of powder-bed fusion technologies, overhangs induce large residual stresses, which eventually entail warpage.

The prevailing remedy in the literature to the issue of overhangs consists in erecting a scaffold structure alongside with the assembly of the desired shape; see e.g. [8]. This auxiliary, sacrificial structure has to be removed at the end of the construction process, which is costly and cumbersome.

A different idea is to impose since the design optimization investigations that the shape should be self-supporting - that is, free of overhangs. In this direction, ad hoc criteria, based on a minimum angle between the structural boundary and the horizontal directions, have been used hitherto in order to tackle this issue [11, 13, 14, 19].

In the recent journal articles [2, 3], we have introduced new formulations to deal with overhang constraints in shape and topology optimization algorithms. The purpose of the present note is to give an introduction to this model; as such, it does not contain any original result. The reader is referred to the aforementioned articles for a complete mathematical exposition, and for the discussion of efficient numerical strategies dedicated to the treatment of overhang constraints.

The key idea of our work is to introduce a new constraint functional $P_{sw}(\Omega)$ for shape optimization problems, which appraises the constructibility of the shape $\Omega$ at each stage of its layer by layer assembly; in particular, overhang constraints are naturally addressed by this formulation. To achieve this, in the setting of the optimization problem, we distinguish the mechanical situation where the *final* (completed) shape $\Omega$ is utilized, on which the optimization criterion is based, and that where $\Omega$, and all the associated *intermediate shapes* $\Omega_h$, are under construction, which guides the definition of our constraint functional.

This note is organized as follows. In Section 2, we introduce the shape optimization problem at stake. In Section 3, we describe the intuitive attempt to rely on geometric functionals to deal with overhang constraints, and explain why it proves insufficient for this purpose. Then, in Section 4, we present our mechanical model of the context in which $\Omega$ is constructed, and we formulate our manufacturing constraint functional $P_{sw}(\Omega)$ accordingly. Several details of its mathematical analysis are outlined, notably the calculation of its shape derivative, and an algorithm taking advantage of its intrinsic structure is proposed, which allows to accelerate significantly the (costly) calculations entailed by the evaluation of $P_{sw}(\Omega)$ and its derivative. Eventually, a numerical validation and several examples are provided in Section 5.
## 2 Presentation of the shape optimization problem

A shape is a bounded, regular domain $\Omega \subset \mathbb{R}^d$, $d = 2, 3$, filled with a linear elastic material with Hooke's law $A$. In the context of its final utilization, $\Omega$ is clamped on a subset $\Gamma_D \subset \partial\Omega$, and it is submitted to surface loads $f \in L^2(\Gamma_N)^d$ applied on a region $\Gamma_N$ of $\partial\Omega$ disjoint from $\Gamma_D$. The remaining part $\Gamma := \partial\Omega \setminus \overline{\Gamma_D \cup \Gamma_N}$ is traction-free. The elastic displacement $u_\Omega$ is the unique solution in $H^1_{\Gamma_D}(\Omega)^d := \{u \in H^1(\Omega)^d, u = 0 \text{ on } \Gamma_D\}$ to the mechanical system:

$$ \left\{ \begin{array}{ll} -\mathrm{div}(Ae(u_\Omega)) = 0 & \text{in } \Omega, \\ u_\Omega = 0 & \text{on } \Gamma_D, \\ Ae(u_\Omega)n = 0 & \text{on } \Gamma, \\ Ae(u_\Omega)n = f & \text{on } \Gamma_N. \end{array} \right. \qquad (1) $$

For simplicity, the criterion $J(\Omega)$ measuring the performance of shapes is the compliance:

$$ J(\Omega) = \int_{\Omega} Ae(u_{\Omega}) : e(u_{\Omega}) dx = \int_{\Gamma_N} f \cdot u_{\Omega} ds. \qquad (2) $$

Our optimization problem then reads:

$$ \min_{\Omega \in \mathcal{U}_{ad}} J(\Omega), \quad \text{such that} \quad P(\Omega) \le \alpha, \qquad (3) $$

in which

* $\mathcal{U}_{ad}$ is a set of (smooth) admissible shapes $\Omega$, whose boundaries enclose the non optimizable regions $\Gamma_D$, $\Gamma_N$, and another one, $\Gamma_0$, which may overlap $\Gamma_D$ or $\Gamma_N$, and which is defined in Section 4 below;

* $P(\Omega)$ is a constraint functional, meant to enforce the constructibility of shapes by additive processes, whose device is the central issue of the present note;

* $\alpha$ is a tolerance threshold.

Obviously, other constraints could be added to (3), e.g. constraints on the volume $\mathrm{Vol}(\Omega) := \int_{\Omega} dx$, or the perimeter $\mathrm{Per}(\Omega) := \int_{\partial\Omega} ds$ of shapes.

Most popular optimization algorithms for the numerical resolution of (3) rely on the derivatives of $J(\Omega)$ and $P(\Omega)$ with respect to the domain; these are understood in the framework of Hadamard's method (see e.g. [1, 15, 18, 21]): a function $F(\Omega)$ of the domain is shape differentiable if the underlying mapping

$$ \theta \mapsto F(\Omega_{\theta}), \quad \text{where} \quad \Omega_{\theta} := (\mathrm{Id} + \theta)(\Omega), \qquad (4) $$

from $W^{1,\infty}(\mathbb{R}^d, \mathbb{R}^d)$ into $\mathbb{R}$, is Fréchet differentiable at $\theta = 0$; the corresponding derivative is denoted by $F'(\Omega)(\theta)$. Often, the deformations $\theta$ in (4) are restrained to a subset of $W^{1,\infty}(\mathbb{R}^d, \mathbb{R}^d)$; in the following, we shall consider the sets

$$ \Theta^k = \{ \theta \in \mathcal{C}^{k,\infty}(\mathbb{R}^d, \mathbb{R}^d), \theta = 0 \text{ on } \Gamma_D \cup \Gamma_N \cup \Gamma_0 \}, $$

where $k \ge 1$ and $\mathcal{C}^{k,\infty}(\mathbb{R}^d, \mathbb{R}^d)$ is the set of $k$ times continuously differentiable functions from $\mathbb{R}^d$ into itself, whose derivatives up to order $k$ are uniformly bounded.

For example, the shape derivative of (2) reads, for deformations $\theta \in \Theta^k$, $k \ge 1$ (see e.g. [5]):

$$ J'(\Omega)(\theta) = - \int_{\Gamma} Ae(u_{\Omega}) : e(u_{\Omega}) \theta \cdot n ds. $$

# 3 The geometric attempt to deal with overhang constraints

## 3.1 Definition of the geometric functionals

The perhaps most intuitive way to account for overhang constraints - and the most common one in the literature - relies on criteria involving the angle between the normal vector $n \equiv n_\Omega$ to the structural boundary and the build
direction of the additive machine tool. In this section, for simplicity and without loss of generality, this direction is chosen to be the d-th vector $e_d$ of the canonical basis of $\mathbb{R}^d$. In other terms, this geometric description of overhang features brings into play anisotropic perimeter functionals of the form:

$$P_g(\Omega) = \int_{\partial\Omega} \varphi(n_{\Omega}) ds, \qquad (5)$$

where $\varphi: \mathbb{R}^d \to \mathbb{R}$ is a given function of class $\mathscr{C}^1$. Two particular instances of such functions $\varphi$ may be thought off:

* The choice

$$\varphi_a(n) := (n \cdot e_d + \cos \nu)^2, \text{ where } (s)_- := \min(s, 0), \qquad (6)$$

and $\nu$ is a threshold angle, penalizes the regions of $\partial\Omega$ where the angle between the normal vector $n$ and the negative vertical direction $-e_d$ is smaller than $\nu$.

* The choice

$$\varphi_p(n) = \prod_{i=1}^{m} (n - n_{\psi_i})^2, \qquad (7)$$

where the $\psi_i: \mathbb{R}^d \to \mathbb{R}$, $i = 1, ..., m$ are given pattern functions, and $n_{\psi_i} := \frac{\nabla \psi_i}{|\nabla \psi_i|}$ are the corresponding normalized gradients, compels $n$ to be close to at least one of the directions $n_{\psi_i}$.

When it comes to the shape derivative of (5), the result of interest is the following [7].

**Proposition 1.** The function $P_g(\Omega)$ defined by (5) is shape differentiable at any admissible shape $\Omega \in \mathcal{U}_{\text{ad}}$ when deformations $\theta$ are in $\Theta^k$, $k \ge 1$. Its shape derivative reads:

$$P'_g(\Omega)(\theta) = \int_{\Gamma} \kappa \varphi(n) \theta \cdot n ds - \int_{\Gamma} \nabla_{\partial\Omega}(\varphi(n)) \cdot \nabla_{\partial\Omega}(\theta \cdot n) ds,$$

where $\kappa: \partial\Omega \to \mathbb{R}$ is the mean curvature of $\partial\Omega$ and $\nabla_{\partial\Omega}\psi := \nabla\psi - (\nabla\psi \cdot n)n$ is the tangential gradient of a smooth enough function $\psi: \partial\Omega \to \mathbb{R}$.

## 3.2 Insufficiency of the geometric constraint functionals: the ‘dripping effect’

In spite of their simplicity, purely geometric criteria of the form (5) generally fail to prevent the appearance of overhang features.

To illustrate this point, let us anticipate a little on Section 5 where our numerical setting is presented in more details, and consider the benchmark two-dimensional MBB Beam example, as depicted in Figure 2. In a computational domain $D$ with size $6 \times 1$, the structure is anchored at its bottom-right corner, and the vertical displacement is set to 0 at its bottom-left corner. A unit vertical load $f = (0, -1)$ is applied at the middle of its upper side. From the manufacturing point of view, the shape is assembled from bottom to top, i.e. $\Gamma_0$ coincides with the lower side of $D$. Taking advantage of the symmetry of the mechanical problem, only half the working domain $D$ is considered during the optimization process; it is meshed by using $300 \times 100 \mathbb{Q}_1$ elements.

Starting from the initial shape $\Omega_0$ of Figure 3 (top), we solve the compliance minimization problem

$$\begin{array}{ll} \min & J(\Omega) \\ \Omega & \text{s.t. } \mathrm{Vol}(\Omega) \le \alpha_v \mathrm{Vol}(D) \end{array} \qquad (8)$$

with the threshold $\alpha_v = 0.3$ for the volume constraint. The resulting optimized design $\Omega^*$ is depicted in Figure 3 (bottom) and contains large overhanging regions. These overhangs are of great physical significance for the performance of $\Omega^*$; hence, it is expected that their removal will prove difficult.

The results displayed in Figure 4 are typical of the ‘optimized’ shapes resulting from the use of geometric functionals such as $P_a(\Omega)$ or $P_p(\Omega)$ to penalize the overhangs formed by members of such great structural significance. They are obtained by solving the new problem

$$\begin{array}{ll} \min & (1 - \alpha_g) \frac{J(\Omega)}{J(\Omega^*)} + \alpha_g \frac{P_g(\Omega)}{P_g(\Omega^*)}, \\ \Omega & \text{s.t. } \mathrm{Vol}(\Omega) \le \alpha_v \mathrm{Vol}(D) \end{array} \qquad (9)$$

with the parameter $\alpha_g = 0.50$, using
<figure><img src="image_4.png" /><figcaption>Fig. 2 Setting of the two-dimensional MBB beam example.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3 (Top) initial and (bottom) optimized shapes for Problem (8) in the two-dimensional MBB Beam test-case of Section 3.2.</figcaption></figure>

* the function $\varphi_a$ given by (6) and the threshold angle $\nu = 45^\circ$ as for Figure 4 (top),

* the function $\varphi_p$ given by (7) and the pattern functions $\psi_i : \mathbb{R}^2 \to \mathbb{R}$ defined as follows, as for Figure 4 (bottom):

$$ \psi_i(x) := n_{\psi_i} \cdot x, \text{ where } n_{\psi_i} = (\cos v_i, \sin v_i) \text{ and } v_i = -\frac{\pi}{4} + \frac{6i\pi}{40}, \quad i = 0, \dots, 10; \qquad (10) $$

in other terms, the constraint functional $P_p(\Omega)$ imposes that the orientation of the boundary $\partial\Omega$ should comply with that of one of the straight lines with normal vector $n_{\psi_i}$. These directions $n_{\psi_i}$ are uniformly sampled among the set of those making an angle with the negative vertical direction comprised in $[\frac{\pi}{4}, \frac{5\pi}{4}]$.

One first observes that several parts in the resulting designs do comply with the desired orientation, but these are the parts whose structural significance is negligible.

More importantly, oscillations arose on the boundaries of some members, in particular those which are close to horizontal and bear a great part of the loading. This dripping effect has already been observed in the literature; see for instance [19]. In the present situation, the structure contains regions where the minimization of the compliance $J(\Omega)$ urges the formation of large horizontal features; the oscillatory patterns have little impact on the mechanical performance (at least when it is measured in terms of the compliance), but they lead to a dramatic decrease in the values of the geometric functionals $P_a(\Omega)$ and $P_p(\Omega)$. Obviously, the algorithm prefers creating oscillating boundaries than rearranging the large horizontal bars, which would undermine significantly the structural performance of the shape. In this sense, the geometric criteria are not 'strong enough' to steer the algorithm to a different optimization path.
<figure><img src="image_6.png" /><figcaption>Fig. 4 Optimized shapes resulting from Problem (9) in the two-dimensional MBB Beam example, using (top) $\varphi \equiv \varphi_a$ and the threshold angle $\nu = 45^\circ$, and (bottom) $\varphi \equiv \varphi_p$ and the pattern functions $\psi_i$ defined in (10).</figcaption></figure>

# 4 Description and analysis of the mechanical constraint functionals

Assessing the drawbacks of the geometric functionals of Section 3.2, we now introduce a new mechanical constraint functional $P_{sw}(\Omega)$, which relies on a simplified model for the additive construction process.

## 4.1 Definition of the self-weight manufacturing compliance $P_{sw}(\Omega)$

The constraint $P_{sw}(\Omega)$ relies on the mechanical situation of $\Omega$ in the course of the manufacturing process: $\Omega$ is enclosed in a box $D = S \times (0, H)$, $S \subset \mathbb{R}^{d-1}$, representing the build chamber with a vertical build direction $e_d$. For $h \in (0, H)$,

$$ \Omega_h := \Omega \cap \{x = (x_1, \dots, x_d) \in \mathbb{R}^d, 0 < x_d < h\} $$

is the intermediate shape describing the stage where $\Omega$ is assembled up to height $h$. The boundary $\partial\Omega_h$ is here decomposed in a somewhat different fashion from that of Section 2:

$$ \partial\Omega_h = \Gamma_0 \cup \Gamma_h^l \cup \Gamma_h^u, $$

where:

* $\Gamma_0 = \{f \in \partial\Omega_h, x_d = 0\}$ is the contact region between $\Omega$ and the build table,

* $\Gamma_h^u = \{f \in \partial\Omega_h, x_d = h\}$ is the upper side of the intermediate structure,

* $\Gamma_h^l = \partial\Omega_h \setminus \Gamma_0 \cup \Gamma_h^u$ is the lateral surface.

Eventually, we define $l_h := \{x \in \partial\Omega, x_d = h\}$, the part of the boundary $\partial\Omega$ that lies at height $h$; see Figure 5.

Each intermediate shape $\Omega_h$ is clamped on $\Gamma_0$, and is subjected to gravity effects, accounted for by a body force $g \in H^1(\mathbb{R}^d)^d$. Its elastic displacement $u_{\Omega_h}^c \in H_{\Gamma_0}^1(\Omega_h)^d$ satisfies:

$$ \left\{ \begin{array}{ll} -\mathrm{div}(Ae(u_{\Omega_h}^c)) = g & \text{in } \Omega_h, \\ u_{\Omega_h}^c = 0 & \text{on } \Gamma_0, \\ Ae(u_{\Omega_h}^c)n = 0 & \text{on } \Gamma_h^l \cup \Gamma_h^u, \end{array} \right. \qquad (11) $$

so that the self-weight $c_{\Omega_h}$ of $\Omega_h$ reads:

$$ c_{\Omega_h} = \int_{\Omega_h} Ae(u_{\Omega_h}^c) : e(u_{\Omega_h}^c) dx = \int_{\Omega_h} g \cdot u_{\Omega_h}^c dx. \qquad (12) $$

Our constraint $P_{sw}(\Omega)$ of the total structure $\Omega$ aggregates the self-weights of all the intermediate shapes, and therefore deserves the name of self-weight manufacturing compliance:

$$ P_{sw}(\Omega) = \int_{0}^{H} j(c_{\Omega_h}) dh, \qquad (13) $$

where $j: \mathbb{R} \to \mathbb{R}$ is a smooth function (in practice, we use $j(s) = s$).
<figure><img src="image_7.png" /><figcaption>Fig. 5 Intermediate shape $\Omega_h$ at height $h$ during the construction of the final structure $\Omega$: the red zone is the lower boundary $\Gamma_0$ and the blue zone is the upper boundary $\Gamma_h^u$.</figcaption></figure>

**Remark 1.** Let us point out an important bias in this model for the manufacturing process: it implicitly relies on the assumption that every layer of material is constructed instantaneously (and so, it does not take into account the way each such layer is assembled). This feature has important consequences on the interpretation of the numerical results, as we shall see in Section 5.

## 4.2 Shape derivative of $P_{sw}(\Omega)$

The statement of our main result about the shape differentiability of $P_{sw}(\Omega)$ requires additional notations. In this section, $\Omega \in \mathcal{U}_{ad}$ is a given admissible shape. We introduce two open sets $\mathcal{O}_1 \in \mathcal{O}_2$ in $\mathbb{R}^d$ and a smooth function $\chi: \mathbb{R}^d \to \mathbb{R}$ such that:

$$ \{x \in \partial\Omega \setminus \bar{\Gamma}_0, n(x) \cdot e_d = \pm 1\} \subset \mathcal{O}_1, 0 \le \chi \le 1, \chi \equiv 0 \text{ on } \mathcal{O}_1, \text{ and } \chi \equiv 1 \text{ on } \mathbb{R}^d \setminus \bar{\mathcal{O}}_2. $$

In other words, $\mathcal{O}_1$ is a neighborhood of the 'flat regions' of $\partial\Omega \setminus \bar{\Gamma}_0$, and $\chi$ is a cutoff function whereby these regions will be ignored in the analysis.

In this context, the relevant sets for perturbations of $\Omega$ are the Banach spaces

$$ X^k = \{ \theta = \chi \tilde{\theta}, \tilde{\theta} \in \Theta^k \}, \text{ equipped with the quotient norm } ||\theta||_{X^k} = \inf \{ ||\tilde{\theta}||_{\mathcal{C}^{k,\infty}(\mathbb{R}^d, \mathbb{R}^d)}, \theta = \chi \tilde{\theta} \}. \quad (14) $$

In particular, vector fields $\theta \in X^k$ vanish at the points of $\partial\Omega \setminus \bar{\Gamma}_0$ where the normal vector $n$ is parallel to $e_d$.

**Theorem 2.** The functional $P_{sw}(\Omega)$ given by (13) is shape differentiable at $\Omega$, in the sense that the mapping $\theta \mapsto P_{sw}(\Omega_\theta)$, from $X^k$ into $\mathbb{R}$ is differentiable for $k \ge 1$. Its derivative is:

$$ \forall \theta \in X^k, P'_{sw}(\Omega)(\theta) = \int_{\partial\Omega \setminus \bar{\Gamma}_0} \mathcal{D}_{\Omega} \theta \cdot n \, ds, \quad (15) $$

where the integrand $\mathcal{D}_{\Omega}$ is defined by, for a.e. $x \in \partial\Omega \setminus \bar{\Gamma}_0$:

$$ \mathcal{D}_{\Omega}(x) = \int_{x_d}^{H} j'(c_{\Omega_h}) (2g \cdot u_{\Omega_h}^c - Ae(u_{\Omega_h}^c) : e(u_{\Omega_h}^c)) (x) \, dh. \quad (16) $$
Sketch of the proof: The study of the dependence of $P_{sw}(\Omega)$ with respect to the domain $\Omega$ is not standard since it involves a continuum of domains $\Omega_h$, $h \in (0, H)$, obtained from $\Omega$ by truncation at fixed heights $h$. The proof proceeds in three steps:

(i) At first, we prove that $\theta \mapsto P_{sw}(\Omega_\theta)$ is differentiable when variations $\theta$ are horizontal, i.e.

$$ \theta \in X_H^k, X_H^k := \{ \theta \in X^k, \theta \cdot e_d = 0 \}, $$

and we prove that (15) - (16) hold in this case. This is possible since, in this case, the truncation at level $h$ $(\Omega_\theta)_h$ of the deformed shape $\Omega_\theta$ coincides exactly with the deformation $(\mathrm{Id} + \theta)(\Omega_h)$ of the truncated domain $\Omega_h$. Hence, in this particular case, Theorem 2 follows almost directly from classical techniques from shape optimization, as in e.g. [12, 21].

(ii) We prove that an arbitrary deformation $\theta \in \Theta^k$ can be equivalently represented by a horizontal one $\xi(\theta) \in X_H^k$, in the sense that

$$ (\mathrm{Id} + \theta)(\Omega) = (\mathrm{Id} + \xi(\theta))(\Omega), $$

and we analyze the properties of the mapping $\theta \mapsto \xi(\theta)$ (in particular, we calculate its derivative). This relies on the implicit function theorem after rewriting an identity of the form $\Omega_\theta = \Omega_\xi$ under the form $\mathcal{F}(\theta, \xi) = 0$, for $\theta, \xi \in X^k$ and an appropriate mapping $\mathcal{F}$.

(iii) Theorem 2 follows from the conclusions of Steps (i) and (ii) by an application of the chain rule.

<figure><img src="image_8.png" /></figure>

Remark 2. Formulas (15) and (16) have a quite intuitive structure: the shape gradient $\mathcal{D}_\Omega(x)$ of $P_{sw}(\Omega)$ at some point $x \in \partial\Omega \setminus \bar{\Gamma}_0$ involves the shape gradients of all the self-weights $c_{\Omega_h}$ of the intermediate structures $\Omega_h$ lying above $x$, i.e. such that $h > x_d$.

## 4.3 Practical calculation of the mechanical constraint and its derivative

The numerical evaluations of $P_{sw}(\Omega)$ and $P'_{sw}(\Omega)(\theta)$, or equivalently $\mathcal{D}_\Omega$, rely on a discretization of the height interval $(0, H)$ with a sequence $0 = h_0 < h_1 < \dots < h_N = H$. The intuitive, ‘0th-order’ method to calculate approximations $P_N^0$ and $\mathcal{D}_N^0$ of $P_{sw}(\Omega)$ and $\mathcal{D}_\Omega$ consists in replacing $c_{\Omega_h}$ and $u_{\Omega_h}^c$ by piecewise constant quantities on each interval $I_i := (h_i, h_{i+1})$ before applying (13) or (15):

$$ c_{\Omega_h} \approx c_{\Omega_{h_{i+1}}} \quad \text{and} \quad u_{\Omega_h}^c \approx u_{\Omega_{h_{i+1}}}^c \quad \text{on } \Omega_h, \text{ for } h \in I_i. $$

Because this procedure is low-order, the subdivision $\{h_i\}_{i=0,\dots,N}$ of $(0, H)$ has to be quite fine so that the accuracy of this approximation process is guaranteed; this brings about many numerical resolutions of the elasticity system (11) for the $u_{\Omega_{h_i}}^c$, which is very costly.

The computational efficiency of this ‘0th-order’ method can be improved by constructing a higher-order approximation of the mappings $h \mapsto c_{\Omega_h}$ and $h \mapsto u_{\Omega_h}^c$ on each interval $I_i$; this is achieved by an interpolation procedure based on the values and the derivatives of these mappings at the $h_i$, which have then to be calculated.

### 4.3.1 Derivatives of $h \mapsto c_{\Omega_h}$ and $h \mapsto u_{\Omega_h}^c$

Consider a fixed shape $\Omega \in \mathcal{W}_{ad}$, and a height $h \in (0, H)$ satisfying:

$$ \text{For any } x \in \ell_h, \text{ the normal vector } n(x) \text{ is not parallel to } e_d. \qquad (17) $$

Our first result is about the derivative of $h \mapsto c_{\Omega_h}$; see [2] for the proof.

**Proposition 3.** In the above context, the mapping $h \mapsto c_{\Omega_h}$ is differentiable at h and:

$$ \left. \frac{d}{dh} (c_{\Omega_h}) \right|_h = \int_{\Gamma_h^u} (2g \cdot u_{\Omega_h}^c - Ae(u_{\Omega_h}^c) : e(u_{\Omega_h}^c)) ds. \qquad (18) $$

We now turn to the issue of differentiating the mapping $h \mapsto u_{\Omega_h}^c$, an operation which has yet to be given an adequate meaning. In this direction, we prove in [2] that (see also [1, 16] for related notions):
• There exists $t_0 > 0$ and a mapping $(-t_0, t_0) \ni t \mapsto \phi_t$ satisfying the following properties:

(i) For $t \in (-t_0, t_0)$, $\phi_t$ is a diffeomorphism of $\mathbb{R}^d$, mapping $\Omega_h$ onto $\Omega_{h-t}$ such that
$$ \phi_t(\Gamma_0) = \Gamma_0, $$

(ii) The mapping $t \mapsto (\phi_t - \text{Id})$, from $(-t_0, t_0)$ into $W^{1,\infty}(\mathbb{R}^d, \mathbb{R}^d)$ is of class $\mathcal{C}^1$, (19)

(iii) Introducing $V(x) := \frac{d\phi_t(x)}{dt}|_{t=0} \in W^{1,\infty}(\mathbb{R}^d, \mathbb{R}^d)$, one has,

For $x \in \Gamma_h^u$, $V(x) \cdot e_d = -1$, and for $x \in \Gamma_h^l$, $V(x) \cdot n(x) = 0$.

• The mapping $t \mapsto u_{\Omega_h}^c \circ \phi_t$ is differentiable from $(-t_0, t_0)$ into $H_{\Gamma_0}^1(\Omega_h)^d$. Its derivative at $t=0$ is called the *Lagrangian derivative* $Y_{\Omega_h}$ of $h \mapsto u_{\Omega_h}^c$. In general, $Y_{\Omega_h}$ depends on the mapping $\phi_t$ used in its definition.

• The quantity $U_{\Omega_h} := Y_{\Omega_h} - \nabla u_{\Omega_h}^c V$, which we identify as the *Eulerian derivative* of $h \mapsto u_{\Omega_h}^c$ is the natural candidate for defining its 'derivative'. $U_{\Omega_h}$ is the solution in $H_{\Gamma_0}^1(\Omega_h)^d$ to the system:

$$ \left\{ \begin{array}{ll} -\mathrm{div}(Ae(U_{\Omega_h})) = 0 & \text{in } \Omega_h, \\ U_{\Omega_h} = 0 & \text{on } \Gamma_0, \\ Ae(U_{\Omega_h})n = 0 & \text{on } \Gamma_h^l, \\ Ae(U_{\Omega_h})n = \frac{\partial}{\partial n}((Ae(u_{\Omega_h}^c)n) & \text{on } \Gamma_h^u. \end{array} \right. \qquad (20) $$

In particular, $U_{\Omega_h}$ is independent of the mapping $\phi_t$ as long as it fulfills (19).

### 4.3.2 Practical interpolation algorithm for the calculation of $P_{sw}(\Omega)$ and $\mathcal{D}_{\Omega}$

The considerations of Section 4.3.1 suggest the following interpolation procedure for calculating first-order approximations, say $P_N^1$ and $\mathcal{D}_N^1$, of $P_{sw}(\Omega)$ and $\mathcal{D}_{\Omega}$ respectively. This allows for an accurate and computationally efficient calculation of these quantities, which uses a coarser subdivision $\{h_i\}_{i=1,...,N}$ of $(0,H)$ than what is needed for the calculation of the 0th-order approximate values $P_N^0$ and $\mathcal{D}_N^0$.

1. For $i = 0, ..., N$, calculate the compliances $c_{\Omega_{h_i}}$ as (12) and the displacements $u_{\Omega_{h_i}}^c$ by solving (11).

2. For $i = 0, ..., N$, calculate the derivative $\frac{d}{dh}(c_{\Omega_h})|_{h=h_i}$ of the compliance owing to Proposition 3.

3. For $i = 1, ..., N$, calculate the Eulerian derivative $U_{\Omega_{h_i}}$ at $h_i$ by using (20).

4. On each interval $I_i$, $i = 0, ..., N-1$, the compliance $c_{\Omega_h}$ is approximated by a cubic spline $\tilde{c}_i(h)$ which is uniquely (and explicitly) determined by the data:

$$ \tilde{c}_i(h_i) = c_{\Omega_{h_i}}, \quad \tilde{c}_i(h_{i+1}) = c_{\Omega_{h_{i+1}}}, \quad \tilde{c}'_i(h_i) = \left. \frac{d}{dh}(c_{\Omega_h}) \right|_{h_i}, \quad \text{and} \quad \tilde{c}'_i(h_{i+1}) = \left. \frac{d}{dh}(c_{\Omega_h}) \right|_{h_{i+1}}. \qquad (21) $$

5. For $i = 0, ..., N-1$ and $h \in I_i$, $u_{\Omega_h}^c$ is approximated by the function $\tilde{u}_h$ defined for a.e. $x \in \Omega_h$ by:

$$ \tilde{u}_h(x) = u_{\Omega_{h_{i+1}}}^c(x) + (h_{i+1} - h) U_{\Omega_{h_{i+1}}} (x); \qquad (22) $$

notice that the above relation does make sense for $x \in \Omega_h$ regardless of the height $h \in (h_i, h_{i+1})$ since $u_{\Omega_{h_{i+1}}}^c$ and $U_{\Omega_{h_{i+1}}}$ are both well-defined on $\Omega_h \subset \Omega_{h_{i+1}}$.

## 4.4 Alternative mechanical models for the behavior of intermediate shapes

Our definition of the constraint $P_{sw}(\Omega)$ in Section 4.1 is guided by the description (11) of the physical behavior of the intermediate shapes $\Omega_h$ during the construction process: they are only subjected to gravity effects.
Interestingly, different mechanical models could be considered instead of (11), including fictitious ones, thereby giving different focuses to the constraint functional. For instance, one may imagine replacing the uniform gravity load $g$ in (11) by an artificial force, applied only on the upper region of each intermediate shape $\Omega_h$: $\Omega_h$ is now subjected to a body force $g_h$ defined by:

$$g_h(x) = \begin{cases} g & \text{if } x_d \in (h - \delta, h), \\ 0 & \text{otherwise,} \end{cases} \qquad (23)$$

where $\delta > 0$ is a small parameter, and $g \in H^1(\mathbb{R}^d)^d$ is a given function. The elastic displacement $u_{\Omega_h}^a$ of $\Omega_h$ in this context is the unique solution in $H_{\Gamma_0}^1(\Omega_h)^d$ to the system:

$$\left\{ \begin{array}{ll} -\mathrm{div}(Ae(u_{\Omega_h}^a)) = g_h & \text{in } \Omega_h, \\ u_{\Omega_h}^a = 0 & \text{on } \Gamma_0, \\ Ae(u_{\Omega_h}^a)n = 0 & \text{on } \Gamma_h^l \cup \Gamma_h^u. \end{array} \right. \qquad (24)$$

The related upper-weight of $\Omega_h$ then reads:

$$c_{\Omega_h}^a = \int_{\Omega_h} Ae(u_{\Omega_h}^a) : e(u_{\Omega_h}^a) dx = \int_{\Gamma_h^u} g \cdot u_{\Omega_h}^a ds,$$

and the corresponding upper-weight manufacturing compliance is defined by:

$$P_{uw}(\Omega) = \int_{0}^{H} j(c_{\Omega_h}^{a}) dh. \qquad (25)$$

As we shall see in Section 5, this formulation is well-suited when it comes to penalizing more specifically the upper region of each intermediate shape $\Omega_h$.

As far as the shape derivative of $P_{uw}(\Omega)$ is concerned, the exact same proof as that of Theorem 2 can be worked out, taking advantage of the definition (23) of $g_h$, and the conclusions of this Theorem extend verbatim to this new case. Also, using a similar analysis to that of Section 4.3 (and working out similar calculations as in [2]), it is possible to calculate the ‘derivatives’ of the mappings $h \mapsto c_{\Omega_h}^a$ and $h \mapsto u_{\Omega_h}^a$, which results in expressions similar to (18) and (20), up to additional terms accounting for the dependence of the load $g_h$ on $h$.

Let us eventually emphasize that our choice of the linearized elasticity systems (11) or (24) for the description of the behavior of the intermediate shapes is merely incidental. One could for instance rely on a similar construction, involving instead the heat equation, so as to model cooling effects within the intermediate shapes, and thereby residual stresses; see [4] about this idea.

# 5 Numerical illustrations

We end this note with several illustrations of the above theoretical considerations; in Section 5.2, we test the algorithms of Section 4.3 for the evaluation of the constraint $P_{sw}(\Omega)$ and its derivative; the subsequent examples (Section 5.3 and below) discuss several numerical strategies for preventing the appearance of overhang features during the shape optimization process. Before entering the core of the matter, we briefly outline our numerical setting.

## 5.1 A few words about the numerical implementation

In all our examples, we rely on the level set method on a fixed computational mesh to represent shapes and their deformations; see [17,20] about the level set method, and [5,22] about its use in shape and topology optimization. In a nutshell, one shape $\Omega$, enclosed in a larger, fixed working domain $D$, is regarded as the negative subdomain of a scalar ‘level set’ function $\phi: D \to \mathbb{R}$, that is:

$$\left\{ \begin{array}{ll} \phi(x) < 0 & \text{if } x \in \Omega, \\ \phi(x) = 0 & \text{if } x \in \partial\Omega, \\ \phi(x) > 0 & \text{if } x \in D \setminus \bar{\Omega}. \end{array} \right. \qquad (26)$$
The evolution in time of a shape $\Omega(t)$, driven by a velocity field with normal component $V(t,x)$, can be modeled by the following Hamilton-Jacobi equation, the solution $\phi(t,x)$ of which, is a level set function for $\Omega(t)$:

$$ \frac{\partial \phi}{\partial t}(t,x) + V(t,x) |\nabla \phi(t,x)| = 0, \quad t > 0, \quad x \in D. \qquad (27) $$

In our applications, the (scalar) velocity field $V(t,x)$ stems from the resolution of Problem (3) by means of an SLP-type optimization algorithm similar to that presented in [9], and the (pseudo-) time $t$ represents the descent step.

From a numerical point of view, the working domain $D$ is a box in two or three space dimensions; it is discretized by means of a Cartesian mesh $\mathcal{G}$, i.e. $\mathcal{G}$ is composed of square elements in $2d$ and cubes in $3d$. The level set function $\phi$ is discretized at the vertices of $\mathcal{G}$ and the Hamilton-Jacobi equation (27) is solved by using an explicit second-order upwind scheme on $\mathcal{G}$, as presented e.g. in [20].

Since the computational mesh $\mathcal{G}$ is fixed throughout the optimization process of the shape $\Omega$, no mesh of $\Omega$ is available for the Finite Element resolution of linearized elasticity systems of the form (1). To circumvent this difficulty, we rely on the so-called ‘ersatz-material’ approximation which consists in filling the void $D \setminus \bar{\Omega}$ with a very soft material with Hooke’s tensor $\varepsilon A$ (in practice $\varepsilon = 10^{-3}$), thus transferring a system posed on $\Omega$ into an approximate one posed on $D$; see for instance [5].

In all examples, the Young’s modulus of the considered elastic material is normalized as $E = 1$ and the Poisson’s ratio is set to $\nu = 0.33$.

## 5.2 Validation of the approximations of Section 4.3

Our first numerical example aims to assess the computational efficiency of the first-order interpolation algorithm proposed in Section 4.3 for the calculation of the constraint $P_{sw}(\Omega)$ and its derivative $\mathcal{D}_{\Omega}$.

Let us consider again the two-dimensional MBB Beam described in Section 3.2. At first, we calculate the functional $P_{sw}(\Omega)$ and its shape derivative $\mathcal{D}_{\Omega}$ in the particular case where $\Omega = \Omega_0$ (the shape displayed in Figure 3 (top)), by using a uniform subdivision of $(0, H)$ made of 100 layers and the $0^{th}$-order approximation scheme, i.e. we evaluate $P_{100}^0$ and $\mathcal{D}_{100}^0$, which serve as reference values for the comparisons in this section. We then calculate the $0^{th}$- and $1^{st}$-order approximations $P_N^i$ and $\mathcal{D}_N^i$, $i = 0, 1$ associated to several subdivisions of $(0, H)$ made of $N$ intervals with equal length. We are interested in the behavior of the relative errors:

$$ \text{err}(P,N,i) = \frac{|P_N^i - P_{100}^0|}{P_{100}^0}, \quad \text{and} \quad \text{err}(\mathcal{D},N,i) = \frac{||\mathcal{D}_N^i - \mathcal{D}_{100}^0||_{L^2(\partial\Omega\setminus\Gamma_0)}}{||\mathcal{D}_{100}^0||_{L^2(\partial\Omega\setminus\Gamma_0)}}. $$

The results are displayed on Figure 6 (bottom): while the $1^{st}$-order approximation method does not bring a lot of improvement when it comes to evaluating the constraint functional $P(\Omega)$, it allows for a substantial gain (i.e. a faster convergence with respect to the number $N$ of subdivisions) in the evaluation of its derivative.

<figure><img src="image_9.png" /><figcaption>Fig. 6 Relative errors of the $0^{th}$- and $1^{st}$-order approximations of (top) $P_{sw}(\Omega_0)$ and (bottom) its derivative $\mathcal{D}_{\Omega_0}$.</figcaption></figure>
## 5.3 Optimization of the two-dimensional MBB Beam using a constraint on the self-weight manufacturing compliance $P_{sw}(\Omega)$

Still in the context of the two-dimensional MBB Beam described in Section 3.2, we turn to the compliance minimization problem under volume constraint:

$$ \begin{array}{ll} \displaystyle \min_{\Omega} & J(\Omega) \\ \text{s.t.} & \mathrm{Vol}(\Omega) \le \alpha_v \mathrm{Vol}(D), \quad \text{where } J(\Omega) \text{ is the compliance (2).} \end{array} \qquad (28) $$

We first solve (28), starting from the initial design $\Omega_0$ of Figure 3 (top), with threshold value $\alpha_v = 0.3$. For the reader's convenience, the resulting design $\Omega^*$ is reprinted from Figure 3 (bottom) in Figure 7 (a).

We then add our mechanical constraint $P_{sw}(\Omega)$ to the problem (28), and now solve:

$$ \begin{array}{ll} \displaystyle \min_{\Omega} & J(\Omega) \\ \text{s.t.} & \mathrm{Vol}(\Omega) \le \alpha_v \mathrm{Vol}(D), \\ & P_{sw}(\Omega) \le \alpha_c P_{sw}(\Omega^*), \end{array} \qquad (29) $$

for different values of the parameter $\alpha_c$; the results are represented on Figure 7 (b)-(d). Understandingly, decreasing values of $\alpha_c$ result in structures which are more rigid from a manufacturing point of view (i.e. with lower values of $P_{sw}(\Omega)$) and more compliant from a structural perspective (i.e. with higher values of $J(\Omega)$).

Let us observe however that the optimized shapes still contain overhangs, which are mainly concentrated on their upper regions. This phenomenon arises for mainly two reasons:

1. As mentioned in Remark 1, the fact that each layer is assumed to be deposited instantaneously has an impact on the values of the self-weights $c_{\Omega_h}$ given by (12): the rigidity of completely flat overhangs (as those appearing in figures 7 (b)-(d)) is actually overestimated; roughly speaking, our modeling only involves the intermediate stages of the construction process where each layer is completely assembled, and where these flat regions are thereby connected to the lower structure. Hence, all the situations where these regions are hanging over void (for instance, at the moments when they are only half-constructed, and when the self-weight would typically take large values) are overlooked.

2. Secondly, considering a uniform self-weight loading for each intermediate shape $\Omega_h$ leads to a concentration of the shape gradient of $J(\Omega)$ at regions closer to $\Gamma_0$. Intuitively, this comes from the fact that such regions appear in a greater number of intermediate shapes, which is reflected on the shape derivative of $P_{sw}(\Omega)$; see Theorem 2. Thus, the algorithm favours the elimination of the overhangs located in these regions, while upper parts are comparatively less influenced by this requirement.

One remedy for the former issue would consist in using a more precise, 'pixel by pixel' model of the construction process introduced in Section 4.1: one could indeed consider the trajectory of the machine tool during the assembly of each layer of the shape, and take as 'intermediate shapes' several snapshots in the assembly of each layer. Unfortunately, beyond difficult modeling issues, relying on this idea would cause the computational cost (which is already significant) of the constraint functional to increase dramatically. Therefore, this proposition is not examined in the present article.

A second remedy amounts to changing the mechanical problem (11) characterizing the behavior of the intermediate shapes. More precisely, as proposed in Section 4.4, we consider fictitious body forces applied on the upper region of these shapes, thus better penalizing thin horizontal patterns in the intermediate shapes $\Omega_h$. This is the purpose of the next section.

## 5.4 Improvement brought by the upper-weight manufacturing compliance $P_{uw}(\Omega)$

We now take on the two-dimensional MBB Beam example of Section 5.3, replacing the self-weight manufacturing compliance $P_{sw}(\Omega)$ with the upper-weight manufacturing compliance $P_{uw}(\Omega)$ in the formulation (29) of the optimization problem, with the expectation that it shows a better ability when it comes to removing overhanging features.

The thickness $\delta$ of the regions where $g$ is applied in the definition (23) of the body force $g_h$ is of the order of the mesh size: $\delta = \Delta x$.
<figure><img src="image_10.png" /><figcaption>**Fig. 7** Optimized shapes for the two-dimensional MBB Beam example of Section 5.3: (a) optimized shape $\Omega^*$ for Problem (8) (i.e. without additive manufacturing constraints), and optimized shapes for Problem (29) using parameters (b) $\alpha_c = 0.50$, (c) $\alpha_c = 0.30$, and (d) $\alpha_c = 0.10$.</figcaption></figure>

One observes that the algorithm tends to add more features oriented along the build direction and to connect them together by creating small archs, which have optimal rigidity for self-weight loadings. The results are much more satisfying than those obtained by using the functional $P_{sw}(\Omega)$ when it comes to removing overhanging features.

## 5.5 A three-dimensional example

Our last example is about the optimization of a three-dimensional bridge, as depicted in Figure 9: the dimensions of the working domain $D$ are 6 × 1 × 1. The structure is clamped at its lower-right corners, while at its lower-left corners, the vertical displacement is prevented. A uniform load $f = (0,0,-1)$ is applied on the upper side of the structure. From the manufacturing point of view, the shape is constructed from top to bottom, i.e. $\Gamma_0$ coincides with the upper side of $D$ (i.e. the deck of the bridge). Due to the double symmetry of the mechanical problem, only one quarter of $D$ is meshed and the cubic grid $\mathcal{G}$ is composed of 90 × 15 × 30 elements.

We start by minimizing the structural compliance $J(\Omega)$ of the bridge under volume constraint, taking the full working domain $D$ as initial shape; i.e. we solve (8) with $\alpha_v = 0.10$. Figure 10 (left column) shows the resulting
<figure><img src="image_11.png" /><figcaption>**Fig. 8** Optimized shapes for the two-dimensional MBB Beam example of Section 5.4, solving Problem (29) with the upper-weight manufacturing compliance $P_{uw}(\Omega)$ and parameters (a) $\alpha_c = 0.30$, (b) $\alpha_c = 0.10$, (c) $\alpha_c = 0.05$, and (d) $\alpha_c = 0.03$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9** Setting of the three-dimensional bridge test-case.</figcaption></figure>
optimized shape $\Omega^*$, which presents several overhangs of significant structural importance.

We then add a constraint on the upper-weight manufacturing compliance $P_{uw}(\Omega)$, but instead of solving (29) we rather consider the following optimization problem:

$$ \begin{array}{ll} \min & \mathrm{Vol}(\Omega), \\ \text{s.t.} & J(\Omega) \le J(\Omega^*), \\ & P_{uw}(\Omega) \le \alpha_c P_{uw}(\Omega^*), \end{array} \qquad (30) $$

where we look for a shape having at least the same rigidity as $\Omega^*$, but which is also more rigid from a manufacturing perspective. Figures 10 (right column) and 11 (left column) show the optimized shapes for $\alpha_c = 0.70$ and $\alpha_c = 0.10$ respectively. For $\alpha_c = 0.70$ only slight changes are observed: the upper bar becomes thinner and the central bars become thicker and are relocated closer to $\Gamma_0$. For $\alpha_c = 0.10$ the changes are drastic and the same trend as in the two-dimensional MBB beam example is observed.

## Acknowledgements

G. A. is a member of the DEFI project at INRIA Saclay Ile-de-France; the work of C.D. is partially supported by the IRS-CAOS grant from Université Grenoble-Alpes.
<figure><img src="image_13.png" /><figcaption>**Fig. 10** Optimized designs for the three-dimensional bridge example of Section 5.5, (left) without manufacturing constraints, (right) solving Problem (30) with $\alpha_c = 0.7$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. 11** (Left) Different views of the optimized shape for the three-dimensional bridge example of Section 5.5, solving Problem (30) with $\alpha_c = 0.1$; (right) another view on the three-dimensional bridges for Problem (30) with (top) no manufacturing constraint, (middle) $\alpha_c = 0.7$ and (bottom) $\alpha_c = 0.1$.</figcaption></figure>
References

[1] G. ALLAIRE, *Conception optimale de structures*, vol. 58 of Mathématiques & Applications, Springer-Verlag, Berlin, (2007), doi: 10.1007/978-3-540-36856-4.

[2] G. ALLAIRE, C. DAPOGNY, A. FAURE, AND G. MICHAILIDIS, *Shape optimization of a layer by layer me- chanical constraint for additive manufacturing*, C. R. Acad. Sci. Paris, Ser. I, 355, (2017), pp. 699–717, https://doi.org/10.1016/j.crsa.2017.04.008.

[3] G. ALLAIRE, C. DAPOGNY, A. FAURE, G. MICHAILIDIS AND R. ESTEVEZ, *Structural optimization under overhang constraints imposed by additive manufacturing technologies*, submitted, (2017).

[4] G. ALLAIRE, L. JAKABCIN, *Taking into account thermal residual stresses in topology optimization of structures built by additive manufacturing*, (in preparation).

[5] G. ALLAIRE F. JOUVE A.M. TOADER, *Structural optimization using shape sensitivity analysis and a level-set method*, J. Comput. Phys., 194 (2004) pp. 363–393, doi:10.1016/j.jcp.2003.09.032.

[6] A. CLAUSEN, *Topology optimization for additive manufacturing*, doctoral thesis of the Technical University of Denmark, (2016).

[7] C. DAPOGNY, A. FAURE, G. MICHAILIDIS, G. ALLAIRE, A. COUVELAS, AND R. ESTEVEZ, *Geometric constraints for shape and topology optimization in architectural design*, Computational Mechanics, 59, (2017), pp. 933–965, https://doi.org/10.1007/s00466-017-1383-6.

[8] J. DUMAS, J. HERGEL, AND S. LEFEBVRE, *Bridging the gap: automated steady scaffoldings for 3d printing*, ACM Transactions on Graphics (TOG), 33 (2014), p. 98, doi:10.1145/2601097.2601153.

[9] P. DUNNING AND H. KIM, *Introducing the sequential linear programming level-set method for topology optimization*, Structural and Multidisciplinary Optimization, 51 (2015), pp. 631–643, https://doi.org/10.1007/s00158-014-1174-z.

[10] I. GIBSON, D.W. ROSEN AND B. STUCKER, *Additive manufacturing technology: rapid prototyping to direct digital manufacturing*, Springer Science Business Media, Inc, (2010), doi:10.1007/978-1-4419-1120-9.

[11] A.T. GAYNOR AND J.K. GUEST, *Topology optimization considering overhang constraints: Eliminating sacrificial support material in additive manufacturing through design*, Struct. Multidisc. Optim., (2016), doi:10.1007/s00158-016-1551-x.

[12] A. HENROT AND M. PIERRE, *Variation et optimisation de formes: une analyse géométrique*, vol. 48, Springer, 2005, doi:10.1007/3-540-37689-5.

[13] M. LANGELAAR, *Topology optimization of 3D self-supporting structures for additive manufacturing*, Additive Manufacturing, 12, (2016), pp. 60–70, https://doi.org/10.1016/j.addma.2016.06.010.

[14] A. M. MIRZENDEHDEL AND K. SURESH, *Support structure constrained topology optimization for additive manufacturing*, Computer-Aided Design, (2016), http://dx.doi.org/10.1016/j.cad.2016.08.006.

[15] F. MURAT AND J. SIMON, *Sur le contrôle par un domaine géométrique*, Technical Report RR-76015, Laboratoire d'Analyse Numérique (1976).

[16] A. NOVOTNY AND J. SOKOLOWSKI, *Topological derivatives in shape optimization*, Springer, Heidelberg (2012), doi:10.1007/978-3-642-35245-4.

[17] S. OSHER AND J. SETHIAN, *Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations*, Journal of computational physics, 79 (1988), pp. 12–49, https://doi.org/10.1016/0021-9991(88)90002-2.

[18] O. PIRONNEAU, *Optimal shape design for elliptic systems*, Springer-Verlag, New York (1984), https://doi.org/10.1007/BFb0006123.

[19] X. QIAN, *Undercut and overhang angle control in topology optimization: a density gradient based integral approach*, International Journal for Numerical Methods in Engineering, 111, (2017), pp. 247–272, doi: 10.1002/nme.5461.

[20] J. SETHIAN, *Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science*, Cambridge university press, 1999, https://doi.org/10.1108/k.2000.29.2.239.3.

[21] J. SOKOŁOWSKI, J.-P. ZALESIO, *Introduction to shape optimization: shape sensitivity analysis*, Springer Series in Computational Mathematics, Vol. 10, Springer, Berlin (1992), https://doi.org/10.1007/978-3-642-58106-9.

[22] M. WANG, X. WANG, AND D. GUO, *A level set method for structural topology optimization*, Computer methods in applied mechanics and engineering, 192 (2003), pp. 227–246, https://doi.org/10.1016/S0045-7825(02)00559-5.