J Martínez-Frutos

G Allaire

C Dapogny

F Periago
email: periago@upct.es.

Structural optimization under internal porosity constraints using topological derivatives

Keywords: Porosity, structural optimization, compliant mechanisms, manufacturing constraints, shape and topological derivatives, level set method

Porosity is a well-known phenomenon occurring during various manufacturing processes (casting, welding, additive manufacturing) of solid structures, which undermines their reliability and mechanical performance. The main purpose of this article is to introduce a new constraint functional of the domain which controls the negative impact of porosity on elastic structures in the framework of shape and topology optimization. The main ingredient of our modelling is the notion of topological derivative, which is used in a slightly unusual way: instead of being an indicator of where to nucleate holes in the course of the optimization process, it is a component of a new constraint functional which assesses the influence of pores on the mechanical performance of structures. The shape derivative of this constraint is calculated and incorporated into a level set based shape optimization algorithm. Our approach is illustrated by several two-and three-dimensional numerical experiments of topology optimization problems constrained by a control on the porosity effect.

Introduction

Fostered by the developments of scientific computing methods and the soaring in computational resources, optimal design of structures has been a very active research field since the early seventies. Typical structural optimization problems lie at the interface between mathematics, mechanics, physics and engineering; in a nutshell, they aim at finding the optimal distribution of material within a given design domain in the sense that it minimizes a certain objective function or mechanical criterion.

Although efficient algorithms are now available to address a wide variety of structural optimization problems, the predicted optimal shapes are generally 'conceptual', insofar as some real-life issues are ignored in the physical description of the considered situation. Recently, the perspective of bridging this gap between mathematical models and concrete applications has aroused a great enthusiasm among the structural optimization community, and important progress has been made. As such, the inevitable uncertainties on the physical data of the problems -e.g. the applied loads, the mechanical properties of the constituent material, or the geometry of the optimized shape itself -have been modelled and incorporated into the design optimization process in [START_REF] Allaire | A linearized approach to worst-case design in parametric and geometric shape optimization[END_REF][START_REF] Allaire | A deterministic approximation method in shape optimization under random uncertainties[END_REF][START_REF] Dambrine | Computing quantities of interest for random domains with second order shape sensitivity analysis[END_REF][START_REF] Lazarov | Topology optimization with geometric uncertainties by perturbation techniques[END_REF][START_REF] Martínez-Frutos | Robust shape optimization of continuous structures via the level set method[END_REF][START_REF] Martínez-Frutos | Risk-averse structural topology optimization under random fields using stochastic expansion methods[END_REF], to name a few references; see also the short survey [START_REF] Maute | Topology Optimization under uncertainty[END_REF]. Another key issue is to understand and to take into account the constraints imposed by the manufacturing process. A classical example is related to the thickness of the optimized shape: many optimal designs obtained numerically contain too thin members which are likely to break during their assembly by e.g. casting or molding; see for instance [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF][START_REF] Sigmund | Manufacturing tolerant topology optimization[END_REF] and the references therein.

Porosity is another phenomenon which shows up in the course of most fabrication processes of elastic materials, and induces many unwanted effects. For instance, in the context of casting manufacturing methods, porosity is caused by bubbles of air which are trapped during the solidification stage [START_REF] Monroe | Porosity in castings[END_REF] (See Figure 1(a)). This internal porosity may lead to field failures in applications that require pressure-tight components. Porosity is also fairly common in welding [START_REF] Zhou | Porosity formation and prevention in pulsed laser welding[END_REF]. There, it is entailed by small gas holes distributed randomly all over the weld bead (See Figure 1(b)). Porosity as well as other defects closely related to it, such as blow holes, gas pockets and slag inclusions, may weaken the weld. Porosity is also a typical defect for the new additive manufacturing techniques, when the process parameters are not adequately tuned [START_REF] Liu | Current and future trends in topology optimization for additive manufacturing[END_REF], [START_REF] Milewski | Additive Manufacturing of Metals. From Fundamental Technology to Rocket Nozzles, Medical Implants, and Custom Jewelry[END_REF]. For example, porosities can appear in a process like Electron Beam Melting (EBM) [START_REF] Tapia | Prediction of porosity in metal-based additive manufacturing using spatial Gaussian process models[END_REF] (See Figure 1(c)). In this context, the density of pores and their impact on fatigue is not so clear as in standard engineering components; let us indeed mention the recent work [START_REF] Tammas-Williams | The Influence of Porosity on Fatigue Crack Initiation in Additively Manufactured Titanium Components[END_REF] which studies the influence of porosity on fatigue crack initiation in additively manufactured titanium components: the size and location of porosity defects is a crucial factor in the fatigue life of the produced part. According to [START_REF] Tammas-Williams | The Influence of Porosity on Fatigue Crack Initiation in Additively Manufactured Titanium Components[END_REF], the Examples of porosity defects: (a) Shrinkage porosity in casting (source: http://efoundry.iitb.ac.in/Defects/shrinkage.html); (b) weld metal porosity (source: [START_REF] Haboudou | Reduction of porosity content generated during Nd:YAG laser welding of A356 and AA5083 aluminium alloys[END_REF]); (c) Hydrogen pores in a 3D-printed aluminum alloy (source: [START_REF] Weingarten | Formation and reduction of hydrogen porosity during Selective Laser Melting of AlSi10Mg[END_REF]).

proximity to the surface and the pore aspect ratio make the initiating porosity within the top 3% most harmful defects for metallic structures.

Hence, porosity generally affects the reliability and quality of materials, and eliminating these effects -or at least keeping them under an acceptable level -is of utmost importance in real-life applications. To the best of our knowledge, porosity has not been accounted for in the formulation of structural optimization problems so far.

The purpose of the present paper is to propose a formulation of structural optimization problems which includes a control of porosity effects as a constraint, with the goal of computing optimal designs that are robust with respect to the emergence of small holes during the manufacturing process. The key ingredient in our modeling is the concept of topological derivative. By contrast with its traditional role in structural optimization as a means to nucleate small holes inside the shape in the course of the optimization process [START_REF] Allaire | Structural optimization using topological and shape sensitivity via a level-set method[END_REF][START_REF] Allaire | Minimum stress optimal design with the level set method[END_REF][START_REF] Amstutz | A new algorithm for topology optimization using a level-set method[END_REF][START_REF] Novotny | Topological Derivatives in Shape Optimization[END_REF], in the present paper, we use it in order to approximate the mechanical performance of a structure riddled with tiny holes. From this basis, we formulate an integral constraint functional of the domain whereby the behavior of structures (in terms of compliance or von Mises stress, for instance) should not be degraded too much by the presence of small internal holes with arbitrary location. The application of classical methods for shape and topology optimization in this context then allows for the calculation of optimal shapes which are robust to internal porosity effects.

The remainder of this article is organized as follows. Section 2 is devoted to the presentation of the considered structural optimization problems and to the mathematical construction of our constraint functionals enforcing robustness of shapes with respect to porosity effects. Section 3 then deals with the computation of the shape derivatives of our objective and constraint functionals. The main features of our numerical topology optimization algorithm are outlined in Section 4 with a particular emphasis on the augmented Lagrangian algorithm used to incorporate our constraint functionals, and on the level set method for representing shapes and their evolutions. Several twoand three-dimensional numerical results are discussed in Section 5. Eventually, a short conclusion and opening perspectives are provided in Section 6.

Setting of the problem

In this section, we present the various shape and topology optimization problems of elastic structures at stake in this work, with a particular emphasis on the construction of our porosity constraint functionals.

The linearized elasticity system

In our context, a shape Ω is a bounded, Lipschitz domain in R d (d = 2 or 3 in applications), whose boundary ∂Ω is decomposed into three disjoint regions:

∂Ω = Γ D ∪ Γ N ∪ Γ 0 , |Γ D | > 0;
shapes Ω are clamped on the part Γ D and surface loads g ∈ L 2 (Γ N ) d are applied on Γ N . They are additionally submitted to body forces

f ∈ L 2 (R d ) d .
The description of the deformation of Ω classically involves the linearized strain tensor

e(u) = 1 2 ∇u + ∇u T = 1 2 ∂u i ∂x j + ∂u j ∂x i 1≤i,j≤d
associated to the displacement field u : R d → R d ; the stress tensor σ(u) := Ae(u) is inferred via the Hooke's tensor A characterizing the properties of the constituent (isotropic, linear) elastic material of shapes. The latter is defined by, for arbitrary ξ in the set S d (R) of symmetric d × d matrices:

Aξ = 2µξ + λ(Trξ)I,
where λ and µ are the Lamé moduli of the material, satisfying µ > 0 and λ + 2µ/d > 0, and I is the d × d identity matrix. Note that the physical behavior of the material may be equivalently described in terms of its bulk modulus κ = λ + 2µ/d, and its Young's modulus E and Poisson's ratio ν, which are related to λ and µ as follows:

E = 4κµ κ + µ , ν = κ -µ κ + µ if d = 2, and E = 9κµ 3κ + µ , ν = 3κ -2µ 2(3κ + µ) if d = 3. (2.1)
In this situation, the displacement u Ω of a shape Ω belongs to the space

H 1 Γ D (Ω) d , where H 1 Γ D (Ω) = u ∈ H 1 (Ω), u = 0 on Γ D ,
and it is the unique solution in this space to the system of linearized elasticity:

           -div(Ae(u Ω )) = f in Ω, u Ω = 0 on Γ D , (Ae(u Ω ))n = g on Γ N , (Ae(u Ω ))n = 0 on Γ 0 , (2.2)
where n is the unit normal vector to ∂Ω, pointing outward Ω. Equivalently, (2.2) rewrites, under variational form, as the search for the unique

u Ω ∈ H 1 Γ D (Ω) d such that: ∀v ∈ H 1 Γ D (Ω) d , Ω Ae(u Ω ) : e(v) dx = Ω f • v dx + Γ N g • v ds. (2.3) Remark 1.
In the developments ahead, as is customary in shape optimization, we shall require more regularity from f and g so as to ensure the differentiability with respect to the domain of the various considered objective functions. Without entering too much into mathematical details, we shall simply assume henceforth that f and g are smooth enough.

Shape optimization of elastic structures

The generic form of the structural optimization problems considered in this article is: min

Ω∈U ad J(Ω) s.t. Vol(Ω) = V T , C(Ω) = 0. (2.4)
Loosely speaking, (2.4) is about finding a shape Ω which minimizes a certain objective function J (Ω) within a class U ad of admissible shapes. Two constraints are added: the first one of them imposes that the volume Vol(Ω) = Ω dx of Ω match with a user-defined target value V T ; the second one is encoded by a shape functional C(Ω) meant to impose that the mechanical performance of the shape is not too much degraded by porosity effects. The construction of C(Ω) is discussed in details in the next Section 2.3. As far as the set U ad of admissible shapes is concerned, introducing a fixed working domain D ⊂ R d , we take:

U ad = {Ω ⊂ D smooth, Γ D ∪ Γ N ⊂ ∂Ω} . (2.5)
Let us emphasize that Γ D and Γ N are kept fixed during the optimization process so that the only part of ∂Ω which is subject to optimization is the traction-free region Γ 0 . When it comes to the objective criterion J(Ω), it accounts for three possible functionals. The first one of them is the compliance

J 1 (Ω) = Ω f • u Ω dx + Γ N g • u Ω ds, (2.6)
as a measure of the global rigidity of the structure Ω. In the case where Ω stands for a compliant mechanism, the considered cost functional aims at controlling the displacement of Ω on a region Γ of its boundary which is not subject to optimization (i.e. Γ is also imposed to be a part of ∂Ω in the definition (2.5) of U ad in this case):

J 2 (Ω) = Γ F Γ (x, u Ω (x)) dx, (2.7)
where

F Γ : Γ × R d u →
R is a smooth function satisfying adequate growth conditions. Eventually, we shall be interested in the control of the negative effects of porosity on the stress of the structure. To this end, we will use the functional

J 3 (Ω) = Ω j (x, σ (u Ω )) dx, (2.8)
where j : R d x × S d (R) → R is a smooth function satisfying, again, growth conditions; from the physical point of view, J 3 (Ω) may represent the equivalent von Mises stress intensity in the shape Ω, or for the discrepancy between the actual stress σ(u Ω ) inside Ω and a target stress tensor σ 0 .

To simplify notations slightly, we shall assume in the next sections that J(Ω) is positive, but the developments extend readily to the general case.

A constraint functional for the appearance of porosity effects

In the present work, we search for shapes which are optimal with respect to one of the aforementionned criteria, but whose performance is also robust with respect to the emergence of small holes (or bubbles of void) in the course of the manufacturing process.

Let J(Ω) stand for any of the three objective functions (2.6), (2.7) or (2.8); we demand that the value J(Ω) be not too dramatically increased when small holes are nucleated inside Ω.

More precisely, let x ∈ Ω and let ω ⊂ R d be a smooth open bounded set containing the origin, standing for the rescaled pattern of the hole. For ρ > 0, we consider the hole ω x,ρ = x + ρω and the perforated domain Ω x,ρ = Ω \ ω x,ρ . We impose that

J (Ω x,ρ ) ≤ (1 + η) J (Ω) , (2.9)
where η > 0 is a user-defined tolerance parameter. (For simplicity we assume that J(Ω) is always non-negative, which ensures that (2.9) implies that the nucleation of a small hole does not increase too much the objective.) In order to turn the requirement (2.9) into a constraint functional, we rely on the notion of topological derivative.

Definition 2.1. [Topological derivative] Under the above notations and assuming a traction-free (Neumann) boundary condition on the small hole ω x,ρ , the objective function J(Ω) has a topological derivative D T J(x) ∈ R at x if the following asymptotic expansion holds in the vicinity of ρ = 0:

J (Ω x,ρ ) = J (Ω) + ρ d D T J (x) + o ρ d , where lim ρ→0 |o(ρ d )| ρ d = 0. (2.10)
Using (2.10), (2.9) is equivalent to:

D T J (x) ≤ η ρ d J (Ω) - o(ρ d ) ρ d ,
and so, up to the leading order in ρ, (2.9) rewrites: .11) We now aim to impose that (2.11) hold at any point x ∈ Ω, and to this end, we reformulate this pointwise constraint by means of an integral functional as in [START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF]: that (2.11) hold for all x ∈ Ω is equivalent to the integral equality

D T J (x) -η ρ J (Ω) ≤ 0, with η ρ = η/ρ d . ( 2
Ω [D T J (x) -η ρ J (Ω)] + 2 dx = 0, (2.12)
where [t] + := max(t, 0) is the positive part of t ∈ R.

All things considered, the considered structural optimization problems in this work are:

(P ij )      min Ω∈U ad J i (Ω) , i = 1, 2 such that Vol(Ω) = V T , C j (Ω) = 0, j = 1, 3.
where the J i (Ω), i = 1, 2 are given by (2.6), (2.7), and the constraint functionals C j (Ω) read:

C j (Ω) = 1 2 Ω c Ω,j (x) dx, with c Ω,j (x) := [D T J j (x) -η ρ J j (Ω)] + 2 (x), j = 1, 3. (2.13)
Likewise, we denote by (P i ) the problem of minimizing J i (Ω) under the volume constraint Vol(Ω) = V T -that is, the counterpart of (P ij ) where the porosity constraint is omitted.

Remark 2.

• Of course, the stress-based functional J 3 (Ω) could also be taken as the objective function of the problems (P ij ), and J 2 (Ω) could be turned into a porosity constraint of the form (2.13).

These two possibilities add no particular difficulty to the approach proposed in this paper; howver, since they are not illustrated in Section 5, we do not discuss them further.

• Strictly speaking, the pointwise constraint (2.11) (or equivalently its integral counterpart (2.12)) imposes that the performance J(Ω) of Ω be not degraded by more than a certain percentage η when a single hole ω with size ρ is inserted in Ω. The same derivation as above featuring a set ω of the form ω = ∪ i ω xi,ρ , ω xi,ρ = x i + ρω i with evident notations (and working with the natural extension of the notion of topological derivative to this case) reveals that (2.11) and (2.12) actually account for the robustness of Ω with respect to the inclusion of finitely many holes with size ρ, and arbitrary locations inside Ω.

• We have been a little elusive in our Definition 2.1 of the topological derivative. In the cases of interest in this work, J(Ω x,ρ ) involves the elastic displacement u Ωx,ρ , the solution to the counterpart of (2.2) posed on the perforated domain Ω x,ρ . Hence, the topological derivative depends on the shape of the model hole ω and on the boundary conditions imposed on the boundary ∂ω x,ρ of the hole in (2.2), which we have not yet specified. For simplicity, from now on in this paper it is assumed that ω is the unit ball in R d . Homogeneous Neumann boundary conditions are imposed on the boundary ∂ω x,ρ of the hole, which account for lack of material inside ω x,ρ . Imposing homogeneous Dirichlet conditions instead would represent a weld, or a rivet -an interesting case which is not addressed in this work. Note that one could in a similar way impose robustness of shapes with respect to inclusions ω x,ρ of another material (not void) inside Ω.

• In (2.9), (2.10) and (2.11), we have assumed that ω ρ is an "interior" hole, i.e. for small enough ρ > 0, ω x,ρ ∩ ∂Ω = ∅. The reader is referred to [START_REF] Samet | The topological asymptotic with respect to a singular boundary perturbation[END_REF] for an explicit computation of the topological derivative associated to holes intersecting the boundary of Ω (i.e. where x ∈ ∂Ω, using our notations). Alternatively, imposing shapes to be robust with respect to such holes could be achieved in the framework of robustness with respect to geometric uncertainties (see for instance [START_REF] Allaire | A deterministic approximation method in shape optimization under random uncertainties[END_REF] and the references therein).

3 Shape derivatives of the considered objective and constraint functionals

The application of a steepest-descent method for the numerical resolution of the problems (P ij ) requires the calculation of the shape derivatives of the objective and constraint functions (2.6), (2.7), (2.8) and (2.13), which is the purpose of this section. We start by recalling a few basic facts about shape derivatives. Their definition is based on the so-called Hadamard boundary variation method [START_REF] Allaire | Conception optimale de structures[END_REF][START_REF] Henrot | Variation et optimisation de formes, une analyse géométrique[END_REF][START_REF] Murat | Sur le contrôle par un domaine géométrique[END_REF][START_REF] Sokolowski | Introduction to Shape Optimization : Shape Sensitivity Analysis[END_REF], whereby variations of a shape Ω are considered under the form

(Id + θ)(Ω), θ ∈ W 1,∞ (R d , R d ), ||θ|| W 1,∞ (R d ,R d ) < 1. Definition 3.1 (Shape derivative). A function F (Ω) of the domain is shape differentiable at a particular shape Ω if the underlying mapping θ → F ((Id + θ)(Ω)), from W 1,∞ (R d , R d ) into R is Fréchet differentiable at θ = 0. The corresponding Fréchet derivative θ → F (Ω)(θ)
is the shape derivative of F (Ω), and the following expansion holds in the neighborhood of 0

∈ W 1,∞ (R d , R d ): F ((Id + θ) (Ω)) = F (Ω) + F (Ω) (θ) + o (θ) , with lim θ→0 |o (θ) | θ W 1,∞ (R d ,R d ) = 0.
As an example of the shape derivative of a particular functional -which will come in handy for later purposes -let us state the following result; see e.g. [1, Proposition 6.22] for a proof. Theorem 3.2. Let Ω be a smooth bounded open domain and let φ(x) ∈ W 1,1 (R d ). The shape functional

F (Ω) = Ω φ(x) dx
is shape differentiable at Ω and its shape derivative is given by

F (Ω)(θ) = ∂Ω θ(x) • n (x) φ(x) ds, θ ∈ W 1,∞ (R d , R d ). (3.1)
In practice, we shall work with a subset Θ ad ⊂ W 1,∞ (R d , R d ) of admissible deformations, in such a way that variations (Id+θ)(Ω) of an admissible shape Ω ∈ U ad stay admissible. In particular, since the regions Γ D and Γ N of the boundary of the admissible shapes Ω are kept fixed, we shall assume that admissible deformation vector fields θ ∈ Θ ad are 'smooth enough', and satisfy θ = 0 on Γ D ∪ Γ N .

We now study separately the shape differentiability of the considered three cost functionals (2.6), (2.7), (2.8), and of the associated porosity constraint functions (2.13). So as to keep notations simple insofar as possible, we reuse from one situation to the other the notations p Ω , q Ω and z Ω to stand for the various adjoint states arising in the expressions of the shape derivatives.

3.1 Shape derivatives of the compliance functional J 1 (Ω) and the constraint C 1 (Ω)

The shape and topological derivatives of (2.6) have been already computed in [START_REF] Allaire | Structural optimization using shape sensitivity analysis and a level-set method[END_REF]Th. 3.6] and in [START_REF] Garreau | The topological asymptotic for PDE systems: the elasticity case[END_REF][START_REF] Sokolowski | Topological derivatives of shape functionals for elasticity systems[END_REF], respectively; we quote the result:

Proposition 3.1. The functional J 1 (Ω) is shape differentiable at any admissible shape Ω ∈ U ad and its shape derivative is:

J 1 (Ω) (θ) = Γ0 (2f • u Ω -Ae (u Ω ) : e (u Ω )) θ • n ds. (3.2)
Proposition 3.2. The functional J 1 (Ω) has a topological derivative at any point x ∈ Ω, whose expression reads:

D T J 1 (x) = T d 1 (u Ω )(x)
, where T d 1 (u) is defined by:

T d 1 (u) = α d 1 Ae(u) : e(u) + β d 1 Tr(Ae(u))Tr(e(u)) + γ d 1 f • u,
and the coefficients α d 1 , β d 1 and γ d 1 depend on the dimension d:

α 2 1 = 4µ π(λ + 2µ) 2µ(λ + µ) , β 2 1 = (λ -µ) π(λ + 2µ) 2µ(λ + µ) , γ 2 1 = -2π, and
α 3 1 = 20µ π (λ + 2µ) µ (9λ + 14µ) , β 3 1 = (3λ -2µ) π (λ + 2µ) µ (9λ + 14µ) , γ 3 1 = - 8π 3 .
Based on these results, we now calculate the shape derivative of the porosity constraint functional C 1 (Ω) defined by (2.13), which amounts, grossly speaking, to calculating the shape derivative of the topological derivative of J 1 (Ω). To address at the same time the cases d = 2 and d = 3, it is convenient to introduce the functional

C 1 (Ω) := 1 2 Ω [T (u Ω ) -η ρ J 1 (Ω)] 2 + dx, (3.3)
where T (u) := α Ae(u) : e(u) + β Tr(Ae(u))Tr(e(u)) + γ f • u and α, β, γ are arbitrary, fixed real numbers.

Theorem 3.3. The functional C 1 (Ω) defined by (3.3) is shape differentiable at any admissible shape Ω ∈ U ad and its shape derivative reads:

C 1 (Ω) (θ) = 1 2 Γ0 [T (u Ω ) -η ρ J 1 (Ω)] 2 + θ • n ds -η ρ Ω [T (u Ω ) -η ρ J 1 (Ω)] + ds J 1 (Ω)(θ) + Γ0 (Ae(u Ω ) : e(p Ω ) -f • p Ω ) θ • n ds, (3.4)
where J 1 (Ω) (θ) is given by (3.2) and the adjoint state p Ω ∈ H 1 Γ D (Ω) d is the unique solution to the variational problem:

∀v ∈ H 1 Γ D (Ω) d , Ω Ae(p Ω ) : e(v) dx = - Ω [T (u Ω ) -η ρ J 1 (Ω)] + T (u Ω )(v) dx, (3.5)
with

T (u)(v) = 2α Ae(u) : e(v) + 2β Tr(Ae(u))Tr(e(v)) + γ f • v. (3.6)
Proof. The shape differentiability of C 1 (Ω) follows from classical arguments involving the implicit function theorem as presented for instance in [START_REF] Henrot | Variation et optimisation de formes, une analyse géométrique[END_REF]; here, we limit ourselves with a formal calculation of its shape derivative using the method of Céa (see for instance [START_REF] Allaire | Conception optimale de structures[END_REF][START_REF] Céa | Conception optimale ou identification des formes, calcul rapide de la dérivée directionnelle de la function coût[END_REF]).

For Ω ∈ U ad and (u,

p) ∈ H 1 Γ D (R d ) d × H 1 Γ D (R d ) d
, we consider the Lagrangian:

L (Ω, u, p) = 1 2 Ω [T (u) -η ρ J 1 (Ω)] 2 + dx + Ω Ae(u) : e(p) dx - Ω f • p dx - Γ N g • p ds, (3.7)
and for given Ω ∈ U ad , we look after the stationary points (u, p) of L(Ω, •, •).

• Equating the partial derivative ∂L ∂p (Ω, u, p) to 0 allows to infer that u = u Ω , the unique solution to (2.2).

• Equating the partial derivative ∂L ∂u (Ω, u, p) to 0 leads to:

∀ u ∈ H 1 Γ D (R d ) d , Ω Ae(p) : e( u) dx + Ω [T (u) -η ρ J 1 (Ω)] + T (u)( u) dx = 0,
where T (u)( u) is given in (3.6) (Notice that the calculation of (3.6) uses the isotropy of the Hooke's tensor A). This allows to characterize p = p Ω , the unique solution to the adjoint variational problem (3.5).

Eventually, to calculate the shape derivative of C 1 (Ω), we rely on the fact that, for arbitrary

p ∈ H 1 Γ D (R d ) d , C 1 (Ω) = L(Ω, u Ω , p);
thence, a formal differentiation with respect to Ω yields:

C 1 (Ω)(θ) = ∂L ∂Ω (Ω, u Ω , p)(θ) + ∂L ∂u (Ω, u Ω , p)(u Ω (θ)),
where u Ω (θ) is the 'Eulerian derivative' of the mapping Ω → u Ω . Now choosing p = p Ω in the above identity, the last term in the right-hand side of the above identity vanishes and:

C 1 (Ω)(θ) = ∂L ∂Ω (Ω, u Ω , p Ω )(θ),
which leaves us with the calculation of the above derivative. To this end, Theorem 3.2 yields, after a little work:

C 1 (Ω)(θ) = 1 2 Γ0 [T (u Ω ) -η ρ J 1 (Ω)] 2 + θ • n ds -η ρ Ω [T (u Ω ) -η ρ J 1 (Ω)] + ds J 1 (Ω)(θ) + Γ0 (Ae(u Ω ) : e(p Ω ) -f • p Ω ) θ • n ds,
which is the expected result.

3.2 Shape derivative of the cost functional J 2 (Ω)

The following result, proved in e.g. [START_REF] Allaire | Structural optimization using shape sensitivity analysis and a level-set method[END_REF][START_REF] Laurain | A level set-based structural optimization code using FEniCS[END_REF], supplies the expression of the shape derivative of J 2 (Ω).

Proposition 3.3. The functional J 2 (Ω) is shape differentiable at any admissible shape Ω ∈ U ad and its shape derivative reads:

J 2 (Ω) (θ) = Γ0 Ae (u Ω ) : e (p Ω ) θ • n ds, (3.8)
where the adjoint state p Ω is the unique solution in H 1 Γ D (Ω) d to the variational problem:

∀v ∈ H 1 Γ D (Ω) d , Ω Ae(p Ω ) : e(v) dx = - Γ F Γ (x, u Ω )(v) ds, (3.9)
where F Γ (x, u)(v) is the derivative of u → F Γ (x, u) in the direction v.

Remark 3.

• Let us recall that the region Γ ⊂ ∂Ω which is the support of the integral featured in J 2 (Ω) is not subject to optimization (hence, θ = 0 on Γ in (3.8)).

• As we have already mentioned, we shall not deal with the porosity constraint associated to J 2 (Ω) via (2.13), even though it would be no more difficult than in the considered cases in this article.

3.3 Shape derivative of the stress-based cost functional J 3 (Ω) and of the constraint C 3 (Ω)

The shape derivative of J 3 (Ω) has been computed in [START_REF] Allaire | Minimum stress optimal design with the level set method[END_REF]Th. 3.3]. Precisely, one has:

Proposition 3.4. The functional J 3 (Ω) defined by (2.8) is shape differentiable at any admissible shape Ω ∈ U ad and its derivative reads:

J 3 (Ω) (θ) = Γ0 (j (x, σ(u Ω )) + Ae (u Ω ) : e (p Ω ) -f • p Ω ) θ • n ds, (3.10)
where the adjoint state p Ω is now defined as the unique solution in H 1 Γ D (Ω) d to the system:

   -div(Ae(p Ω )) = div(Aj (x, σ(u Ω ))) in Ω p Ω = 0 on Γ D (Ae(p Ω ))n = -(Aj (x, σ(u Ω ))) n on Γ N ∪ Γ 0 , (3.11)
where j (x, σ) ∈ S d (R) is the partial derivative of j(x, σ) with respect to σ. This equivalently rewrites, in variational form:

∀v ∈ H 1 Γ D (Ω) d , Ω Ae(p Ω ) : e(v) dx + Ω Aj (x, σ(u Ω )) : e(v) dx = 0.
From now on in this section, we restrict ourselves to the case where d = 2 and j(x, σ) = ζ(x)||σ|| 2 , where ζ is a smooth function devoted to localizing the stress in a specific region of Ω. This is, to the best of our knowledge, the only situation where an explicit expression of the topological gradient of J 3 (Ω) is known [START_REF] Hlavacek | On Topological Derivatives for Elastic Solids with Uncertain Input Data[END_REF][START_REF] Sokolowski | Topological derivatives of shape functionals for elasticity systems[END_REF].

Proposition 3.5. The topological derivative of J 3 (Ω), in the case j(x, σ) = ζ(x)||σ|| 2 and d = 2, is given by:

D T J 3 (x) = T 2 3 (u Ω , p Ω )(x), (3.12)
where we have introduced: The shape derivative of the porosity constraint functional C 3 (Ω) defined by (2.13) is provided in the following theorem. Its proof is on all points similar to that of Theorem 3.3, and it is therefore omitted.

T 2 3 (u, p) := -2π f • p + 4Ae(u) : Ae(u) -(Tr(Ae(u))) 2 + 1 E (4Ae (
Theorem 3.4. Let Ω ∈ U ad ; the porosity constraint functional C 3 (Ω) is shape differentiable at Ω and its derivative reads:

C 3 (Ω) (θ) = 1 2 Γ0 [T 2 3 (u Ω , p Ω ) -η ρ J 3 (Ω)] 2 + θ • n ds -η ρ Ω [T 2 3 (u Ω , p Ω ) -η ρ J 3 (Ω)] + dx J 3 (Ω)(θ) + Γ0 (Ae(u Ω ) : e(z Ω ) -f • z Ω )θ • n ds + Γ0 (Ae(p Ω ) : e(q Ω ) + ζσ(u Ω ) : σ(q Ω )) θ • n ds, (3.13
) where the first adjoint state p Ω solves (3.11), and the second q Ω and third z Ω adjoint states are the unique solutions in H 1 Γ D (Ω) d to the respective variational problems:

∀v ∈ H 1 Γ D (Ω) d , Ω Ae(q Ω ) : e(v) dx = - Ω [T 2 3 (u Ω , p Ω ) -ηJ 3 (Ω)] + ∂T 2 3 ∂p (u Ω , p Ω )(v) dx, (3.14)
and Remark 4. The fact that the shape derivative of the constraint functional C 3 (Ω) depends on 3 adjoint states is not a surprise. A similar phenomenon appeared and is explained in [START_REF] Allaire | A linearized approach to worst-case design in parametric and geometric shape optimization[END_REF].

∀v ∈ H 1 Γ D (Ω) d , Ω Ae(z Ω ) : e(v) dx = - Ω [T 2 3 (u Ω , p Ω ) -ηJ 3 (Ω)] + ∂T 2 3 ∂u (u Ω , p Ω )(v) dx - Ω ζ(x)σ(q Ω ) : σ(v) dx.
4 Level-set based optimization algorithm

Augmented Lagrangian minimization algorithm

The numerical resolution of any of the constrained optimization problems (P ij ), for given i = 1, 2, j = 1, 3, rests on a simple augmented Lagrangian algorithm, as in our previous work [START_REF] Allaire | Mutli-phase structural optimization via a level-set method[END_REF]; see [START_REF] Nocedal | Numerical optimization[END_REF] for a general presentation of this method. In a nutshell, the resolution of (P ij ) is replaced by that of a series of unconstrained minimization problems (whose features at iteration n -associated shape Ω n , etc. -are denoted with a n superscript) of the augmented Lagrangian functional L ij (Ω, n j , γ n j ) defined by:

L ij (Ω, j , γ j ) = J i (Ω) -1 j C j (Ω) -2 j (Vol(Ω) -V T ) + γ 1 j 2 C j (Ω) 2 + γ 2 j 2 (Vol(Ω) -V T ) 2 , (4.1)
where j = 1 j , 2 j and γ j = γ 1 j , γ 2 j are estimates of the Lagrange multipliers and penalty parameters for the associated constraint, respectively.

The shape derivative of L ij is given by:

L ij (Ω, j , γ j ) (θ) = J i (Ω) (θ) -1 j C j (Ω) (θ) -2 j Γ0 θ • n ds + γ 1 j C j (Ω) C j (Ω) (θ) + γ 2 j (Vol(Ω) -V T ) Γ0 θ • n ds, (4.2)
where the shape derivatives J i (Ω) (θ) are given by (3.2), (3.8), and C j (Ω) (θ) are defined by (3.4), (3.13). The Lagrange multipliers 1 j and 2 j are updated at each iteration n according to the rule:

1 j n+1 = 1 j n -γ 1 j C j (Ω n ) , j = 1, 3, 2 j n+1 = 2 j n -γ 2 j (Vol(Ω n ) -V T ) , j = 1, 3, (4.3)
and the penalty parameters γ 1 j , γ 2 j are increased every 5 iteration.

A brief review of the level set method for shape optimization

Following an original idea by Osher and Sethian [START_REF] Osher | Front propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF], then introduced in the shape optimization context in [START_REF] Allaire | Structural optimization using shape sensitivity analysis and a level-set method[END_REF][START_REF] Sethian | Structural boundary design via level set and immersed interface methods[END_REF][START_REF] Wang | A level set method for structural topology optimization[END_REF], we rely on the level set method for representing our shapes and their deformations. In a nutshell, the level set method consists in representing a domain Ω(t) evolving in time as the negative subdomain of an auxiliary 'level set' function ψ ≡ ψ(t, x) defined on the working domain D in such a way that:

∀x ∈ D, ∀t ∈ (0, T ) ,    ψ(t, x) < 0 if x ∈ Ω(t), ψ(t, x) = 0 if x ∈ ∂Ω(t), ψ(t, x) > 0 if x / ∈ Ω(t). (4.4)
Denoting by V (t, x) the normal component of the velocity field driving the motion of Ω(t), the level set function ψ solves the Hamilton-Jacobi equation:

∂ψ ∂t (t, x) + V (t, x)|∇ψ(t, x)| = 0, t ∈ (0, T ) , x ∈ D. (4.5)
In practice, ψ(t, •) is discretized at the vertices of a Cartesian grid of D and (4.5) may be conveniently solved by using an explicit first-order upwind scheme; see [START_REF] Sethian | Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science[END_REF].

In our context of the resolution of any of the problems (P ij ), t stands for a fictitious time that accounts for the step-size parameter in the descent algorithm; the normal velocity V stems from the shape derivative computed in (4.2), which has the convenient structure:

L ij (Ω, j , γ j ) (θ) = Γ0 g ij θ • n ds. (4.6)
As usual, the velocity V is regularized by substituting the L 2 with the H 1 inner product for the shape derivative. This guarantees a descent direction for the augmented Lagrangian L ij (Ω, j , γ j ) (see [START_REF] Allaire | Structural optimization using shape sensitivity analysis and a level-set method[END_REF] for details).

Remark 5.

• As is customary in the practice of the level set method for shape optimization, the state system (2.2) as well as the different adjoint systems are approximately solved owing to the Finite Element method on the whole working domain D (instead of being exactly solved on Ω, which is impossible in our context since no mesh of Ω is available): the so-called ersatz material method is used, whereby the void region D \ Ω is filled with a very soft material with Hooke's tensor εA, where ε 1.

(4.7)

• The level set function ψ is periodically reinitialized, so that it does not present too steep, or too flat variations which may cause numerical instabilities; see [START_REF] Chopp | Computing minimal surfaces via level set curvature flow[END_REF],

Numerical experiments

In this section, three numerical examples are presented in order to analyze the influence of a control over porosity on the performance of the optimal design of structures and mechanisms. In all the experiments, for each of the problems (P ij ), the Lagrange multipliers and the penalty parameters of the augmented Lagrangian function (4.1) are initialized heuristically according to the following rule:

1 j = -0.1(J i (Ω) /V T ) 2 j = -0.1(C j (Ω) /V T ) γ 1 j = |0.01 1 j /V T | γ 2 j = |0.01 2 j /V T |
In all cases, the problem (P ij ) under scrutiny is solved one first time without taking into account the porosity constraint -i.e. (P i ) is solved. The resulting shape Ω * then serves as the initial design for the resolution of (P ij ), properly speaking. The near-optimal guess provided by Ω * improves the convergence efficiency of the optimization process. Moreover, as it will be shown in the numerical results below, doing so does not put much restriction on the design space.

Eventually, as far as the stopping criterion is concerned, the optimization stops at the first iteration n ∈ N for which the following two conditions are satisfied: elements. The shape is clamped at its upper side and a vertical load g = (0, -3) is applied at the middle of its right-hand side (see Figure 2(a)). The initial shape is displayed in Figure 2(b).

max 1≤k≤5 {|L n ij -L n-k ij | ≤ 0.01L n ij } and |Vol(Ω n ) -V T | ≤ 0.01, where L n ij ≡ L ij (Ω n , n j , γ n j ).

Two-dimensional optimization of L-shaped beam

In this setting, the Young's modulus of the elastic material filling shapes is E = 1 and its Poisson's ratio is ν = 0.3; the value ε = 10 -5 is chosen as for the ersatz material in (4.7).

At first, we solve the compliance minimization problem (P 1 ) under the volume constraint V T = 0.35 Vol(D) (i.e. without any porosity constraint). The resulting shape Ω * is depicted in Figure 4 (a).

Then using Ω * as an initial shape, we add to the problem a porosity constraint associated to the compliance of shapes with a tolerance parameter η = 0.01 (see (2.9)) -i.e. (P 11 ) is solved. We consider porosity holes of size ρ = 0.625, which correspond to the minimum size of one mesh cell. Of course, other values of ρ could be considered but recall that it is only the value η ρ = η/ρ d which matters from a numerical point of view.

Figure 3 shows the evolution of the compliance, the porosity constraint and the volume in the course of the optimization process: within 15 iterations, the value of the porosity constraint C 1 (Ω) is reduced from 6935 to 232, at the cost of worsening the structural compliance J 1 (Ω), which is increased by 6.57%. The optimized designs with and without taking into account the porosity constraint C 1 (Ω) are compared in Figure 4, where the contours of the nodal projection of c Ω,1 are also represented (see (2.13)). Comparing Figures 4 (a) and (b) reveals that the final shape obtained with the porosity constraint presents a smoothed reentrant corner. This modification of the local geometry of the shape entails a reduction of the peak value of c Ω,1 by two orders of magnitude (see Figures 4 (c) and (d)). A more accurate comparison between both designs is presented in Figure 5 (a): only local modifications of the optimal shape near the reentrant corner are performed for this choice η = 0.01 of the tolerance parameter.

❲ ✐ ❤ ♣ ✁ r ✁ s ✐ ② ❝ ✁ ✂ s r t ✐ ✂ ❲ ✐ ❤ ✁ ✉ ♣ ✁ r ✁ s ✐ ② ❝ ✁ ✂ s r t ✐ ✂ (a) (b)
In a second time, the topology optimization problem (P 11 ) is solved again with a lower tolerance value η = 0.001. A reduction of η implies that the regions of the domain which are sensitive to porosity are larger and so the relative importance of the porosity constraint in our problem is increased. The optimized designs in this case are presented in Figure 6. Similarly to the previous case, the contours of c Ω,1 show a peak value in the corner. However, a larger number of regions with a positive value of the constraint are present in this case. As a result, the corner radius of the optimized shape is larger than in the previous case. Furthermore, the regions connecting two or more structural members are larger and more rounded (see Figure 5 (b)).

In the second part of this experiment, the effect of the stress porosity constraint is studied: we now solve Problem (P 13 ) with the value η ρ = 0.0256 in the definition of the constraint C 3 (Ω). This value together with a hole size ρ = 0.625 leads to a tolerance parameter η = 0.01. The results are presented in Figure 7; they are very similar to those obtained previously, by using the compliance porosity constraint C 1 (Ω) with the same value for the tolerance parameter η. Again, a reduction in the maximum value of the integrand c Ω,3 of C 3 (Ω) by two orders of magnitude is achieved by rounding the reentrant corner. It is important to highlight that the modifications in the final shape entail a significant reduction of the stress level

s Ω (x) := ||σ(u Ω )|| 2 (x), x ∈ D, (5.1)
in the structure, as is shown in Figure 7 (d).

3D double hook

Our second experiment takes place in the case of three-dimensional structures. For this purpose, the double hook benchmark is considered, whose details are depicted in Figure 8 (a). The working domain D is tessellated as a regular grid, whose active (resp. inactive) elements are represented in white (resp.gray). The T-shaped active domain is simply-supported at the top of the T-leg and two symmetric vertical loads g = (0, 0, -5) are applied in the middle of the faces at the end of the T-shafts. The initial shape is represented in Figure 8 (b). The elastic material filling shapes is characterized by its Young's modulus E = 1 and Poisson's ratio ν = 0.3; the parameter ε = 10 -4 is chosen for the ersatz material. Taking advantage of the symmetry of the situation, only one-quarter of D is used for the Finite Element model associated to the numerical resolution of the linearized elasticity system (2.2) (and the corresponding adjoint equations). Since we are interested in studying the impact of porosity on shapes, and because the latter is expected to weaken especially regions of the structure with small feature sizes, a low volume constraint V T = 0.04 Vol(D) is considered so as to urge designs towards such characteristics. This raises the need for a high-resolution discretization of the quarter domain, which is meshed with 192 × 48 × 272 hexahedral eight-node linear Finite Elements. As a result the Finite Element model is composed of 2.506.752 Elements and 7.745.283 degrees of freedom. This very high-resolution discretization may cause severe difficulties in terms of memory consumption and computational time. Furthermore, the performance of classical iterative solvers typically deteriorates in situations featuring a high contrast in material properties, which is the case in the present situation, since the ersatz material approximation is used; see Remark 5. To alleviate these issues, the Finite Element system associated to the resolution of (2.2) (and the attached adjoint equation) is solved using the GPU implementation of a matrix-free geometrical multigrid PCG described in [START_REF] Martínez-Frutos | GPU acceleration for evolutionary topology optimization of continuum structures using isosurfaces[END_REF].

In a first time, we minimize the compliance J 1 (Ω) of Ω under the volume constraint Vol(Ω) = V T without adding any porosity constraint-i.e. Problem (P 1 ) is solved. This results in the optimized shape Ω * represented on Figure 10 (a). We then incorporate our porosity constraint C 1 (Ω) attached to the compliance to the problem, for a tolerance parameter η = 0.0003 and ρ = 0.5 (one cell size): starting from Ω * as new initial design, we solve Problem (P 11 ). The evolution of the compliance, the porosity constraint and the volume in the course of the optimization process is depicted in Figure 9 (a-c). From the value C 1 (Ω * ) = 66, the porosity constraint function C 1 (Ω) smoothly converges to the value 0.062 within 20 iterations. During this process, the compliance J 1 (Ω) increased by (only) 1.8%. The isolines of (the nodal projection of) the integrand c Ω,1 of this constraint are shown in Figure 10 for the initial and final designs. As expected, a significant reduction in the peak value of c Ω,1 is achieved. This reduction is obtained by removing thin bars and smoothing sharp regions. This is in agreement with the results obtained in Section 5.1. However, a detailed comparison of the initial and final designs (See Figure 11) reveals a new phenomenon in the reentrant corner of the hook, where the surface in the out-of-plane direction is inflated.

Inverter mechanism

Compliant mechanisms typically feature very small hinges that may break due to defects of the manufacturing process. Therefore, a high level of porosity in the final manufactured design can completely ruin their optimality. In this third experiment, inspired by [START_REF] Allaire | A deterministic approximation method in shape optimization under random uncertainties[END_REF][START_REF] Allaire | Thickness control in structural optimization via a level set method[END_REF][START_REF] Bendsøe | Topology optimization: Theory, methods and applications[END_REF][START_REF] Laurain | A level set-based structural optimization code using FEniCS[END_REF], we address this issue by incorporating our porosity constraint into the shape and topology optimization formulation.

The setting is that depicted in Figure 12 (a): the working domain D is a 100 × 100 square, which is discretized by 128 × 128 Q 1 elements. The considered shapes are clamped at the top and bottom corners of the left edge, and an input force g = (100, 0) is applied at the middle of the left edge. The mechanical setting of this example differs slightly from the general situation of Section 2.1 (see Remark 3): in addition to Γ D and Γ N , the boundary of shapes contains another non optimizable component Γ, located at the center of the right edge. A spring is attached at Γ, with small stiffness k s = 0.1, which has the effect of adding the term

k s Γ u Ω • v ds
to the left-hand side of the variational formulation (2.3) of the linear elasticity system (and a similar modification occurs in the variational formulation (3.9) of the adjoint state).

The material properties in the present case are characterized by E = 1, ν = 0.3 and ε = 10 -2 . The volume constraint is V T = 0.2 Vol(D). The experiment consists in maximizing the horizontal component of the displacement of the inverter mechanism in the 'output' region Γ, while minimizing it in the region Γ N where loads are applied; to this end, we consider the objective function J 2 (Ω) given by (2.7), where the integrand reads:

F Γ (x, u) = α in u 1 χ Γ N (x) + α out u 1 χ Γ (x), (5.2)
where u 1 is the horizontal component of u ∈ R 2 , χ Γ N (resp. χ Γ ) is the characteristic function of Γ N (resp. Γ), and α in = 1, α out = 2.

In the above context, we first solve the minimization problem of J 2 (Ω) under the volume constraint V T = 0.2 Vol(D), without adding any porosity constraint -i.e. we solve (P 2 ). The resulting optimized shape Ω * is represented on Figure 14 (a).

Then, like in Section 5.1, two different experiments are conducted considering porosity constraint in compliance (associated to the function C 1 (Ω)) and stress (based on C 3 (Ω)). The first case, with a tolerance parameter η = 0.05 and a hole size ρ = 0.781, shows a smooth convergence (Figure 13) that completely eliminates the influence of porosity effects: C 1 (Ω) is reduced from 8.3 10 5 to 0. This reduction comes at the expense of a degradation of the target output displacement: the value of the objective J 2 (Ω) is increased by 33%. This degradation is very sensitive to the the value of the tolerance parameter η: on the one hand, large values of η lead to optimal mechanisms with large output displacement, but very sensitive to porosity effects. On the other hand, small values of η make the resulting design more robust against porosity at the cost of degrading the output displacement. These conclusions are also supported by the modifications performed in the optimal shape of the mechanism, which is depicted in Figure 14. Figure 14(c) and(d) show that high values of the nodal projection of c Ω,1 are located in weak regions prone to create hinges. These values are eliminated by increasing the thickness of weak regions, which also increases the stiffness of the mechanism. This behavior is even more apparent when the tolerance parameter η is decreased to 0.01. In the latter case, one observes in Figure 15 and Figure 16 how a reduction of the peak value of c Ω,1 is allowed by more important changes in the final shape. In the second part of this experiment, the effect of the stress porosity constraint C 3 (Ω) on the optimized design is studied: the shape and topology optimization problem (P 23 ) is solved with η = 0.01 and a hole size ρ = 0.781 (one cell size). The numerical results, presented in Figure 17, are in the same spirit as those obtained in Section 5.1. The topology of the optimized shape is very similar to that obtained in the case of the constraint C 1 (Ω) for the same value of η. A reduction by three orders of magnitude in the peak value of c Ω,3 is obtained by strengthening the weak regions. Additionally, one observes in Figure 17 (d) that the modifications in the final shape entail a reduction of the stress level (5.1) in the mechanism.

Conclusions and perspectives

This work is a contribution to the modelling and influence in topology optimization of the porosity phenomenon -the parasitic emergence of small bubbles of void inside an elastic structure in the course of its manufacturing process. To the best of our knowledge, it is the first attempt in the literature to control this phenomenon from the early stage of the design optimization (see section 4.2 in [START_REF] Liu | Current and future trends in topology optimization for additive manufacturing[END_REF]). This work has proposed a new constraint functional, elaborating from the notion of topological derivative, which lends itself to mathematical analysis, and incorporation into 'classical' shape and topological optimization algorithms, so as to urge the optimized shapes to be less sensitive to porosity. The two-and three-dimensional obtained results reveal the very strong impact of the porosity effect on the 'optimality' of designs, and they raise interesting questions in connection with the stress level within the structure. This work could be adapted to various different physical situations, including for instance the case where small grains of another material (and no longer void; see Remark 2) pollute a bulk phase. Of course, our proposed constraint functional depends on a few parameters: the tolerance η, the porosity size ρ and the shape of the hole ω (here, a ball). These parameters could be fitted to experimental data. In particular, anisotropic hole shapes could be interesting to study for modelling micro-cracks.

  Figure 1:Examples of porosity defects: (a) Shrinkage porosity in casting (source: http://efoundry.iitb.ac.in/Defects/shrinkage.html); (b) weld metal porosity (source:[START_REF] Haboudou | Reduction of porosity content generated during Nd:YAG laser welding of A356 and AA5083 aluminium alloys[END_REF]); (c) Hydrogen pores in a 3D-printed aluminum alloy (source:[START_REF] Weingarten | Formation and reduction of hydrogen porosity during Selective Laser Melting of AlSi10Mg[END_REF]).

  u) : Ae(p) -Tr(Ae(u))Tr(Ae(p))) , E is the Young modulus of the considered material (see (2.1)), and the adjoint state p Ω is the solution to (3.11).

(3. 15 ) 3 ∂u 3 ∂p

 1533 In (3.14) and(3.15), the partial derivatives of T 2 3 (u, p) are given by the formulas:∂T 2 (u, p)(v) = -2π 8Ae(u) : Ae(v) -2Tr(Ae(u))Tr(Ae(v)) + 1 E (4Ae(p) : Ae(v) -Tr(Ae(p))Tr(Ae(v))) ,and:∂T 2 (u, p)(v) = -2π f • v + 1 E (4Ae(u) : Ae(v) -Tr(Ae(u))Tr(Ae(v))) .

1 g = 3 DFigure 2 :

 132 Figure 2: Setting of the L-shaped beam problem of Section 5.1: (a) working domain D and boundary conditions; (b) initial shape.

Figure 3 :

 3 Figure 3: Convergence histories for the objective, the porosity and volume constraint functionals in the 2d L-shaped beam test case of Section 5.1.

Figure 4 :

 4 Figure 4: Optimized designs in the L-shaped beam test case of Section 5.1; (left column) without porosity constraint: (a) optimized shape and (c) corresponding distribution of the criterion c Ω,1 ; (right column) including the porosity constraint C 1 (Ω) with ρ = 0.625 and η = 0.01: (b) optimized shape and (d) corresponding distribution of the criterion c Ω,1 .

Figure 5 :

 5 Figure 5: Comparison between the optimized L-shaped beams in Section 5.1, obtained with and without including a porosity constraint C 1 (Ω) associated to the compliance, for a hole size ρ = 0.625 and a tolerance (a) η = 0.01 and (b) η = 0.001.

Figure 6 :

 6 Figure 6: Optimized designs in the L-shaped beam test case of Section 5.1; (left column) without porosity constraint: (a) optimized shape and (c) corresponding distribution of the criterion c Ω,1 ; (right column) including the porosity constraint C 1 (Ω) with η = 0.001: (b) optimized shape and (d) corresponding distribution of the criterion c Ω,1 .

Figure 7 :

 7 Figure 7: Optimized designs in the L-shaped beam test case of Section 5.1; (left column) without porosity constraint: (a) distribution of the criterion c Ω,3 and (c) distribution of the stress (5.1); (right column) including the porosity constraint C 3 (Ω) with ρ = 0.625 and η = 0.01: (b) distribution of the criterion c Ω,3 and (d) distribution of the stress s Ω .

Figure 8 :

 8 Figure 8: Setting of the 3d double hook test case of Section 5.2: (a) computational domain D and boundary conditions, and (b) initial shape.

Figure 9 :Figure 10 :

 910 Figure 9: Convergence histories for the objective, the porosity and volume constraint functionals in the 3d double hook test case of Section 5.2.

Figure 11 :

 11 Figure 11: Comparison between the optimal designs for the 3d double hook problem of Section 5.2 obtained with and without porosity constraint; (a) cut in the plane y = 0; zoom near the corner of (b) the optimized shape Ω * without porosity constraint, and (c) the optimized shape considering the porosity constraint C 1 (Ω).

Figure 12 :

 12 Figure 12: Setting of the inverter test case of Section 5.3: (a) computational domain D and boundary conditions; (b) initial shape.

Figure 13 :

 13 Figure 13: Convergence histories for the optimization problem (P 21 ) in the inverter mechanism test case of Section 5.3 for the threshold parameter η = 0.05.

Figure 14 :Figure 15 :Figure 16 :

 141516 Figure 14: Optimized inverter mechanisms in the test case of Section 5.3, (left column) without any porosity constraint: (a) optimized design and (c) corresponding distribution of the criterion c Ω,1 ; (right column) with a porosity constraint C 1 (Ω) related to the compliance, with η = 0.05: (b) optimized design and (d) corresponding distribution of the criterion c Ω,1 .

Figure 17 :

 17 Figure 17: Optimized inverter mechanisms in the test case of Section 5.3, (left column) without any porosity constraint: (a) distribution of the criterion c Ω,3 and (c) distribution of the stress level s Ω ; (right column) with a porosity constraint C 3 (Ω) related to the stress, with η = 0.01: (b) distribution of the criterion c Ω,3 and (d) distribution of the stress level s Ω .

Acknowledgements. This research was partially supported by the AEI/FEDER and UE under the contracts DPI2016-77538-R and MTM2017-83740 and by the Fundación Séneca (Agencia de Ciencia y Tecnología de la Región de Murcia (Spain)) under the contract 19274/PI/14. The work of C.D. is partially supported by the IRS-CAOS grant from University Grenoble-Alpes. G. A. is a member of the DEFI project at INRIA Saclay Ile-de-France. The work of G. A. is partially supported by the SOFIA project funded by BPI.