Salah-Eddine Zerrouq

Morgan Pierre

Robust optimal ship hulls

Keywords: Michell's wave resistance, Dirichlet energy, probability distribution, robust deterministic control, geometric shape optimization, Newton's method, lumped mass matrix. Mathematics Subject Classification: 49J55, 49M15, 65K10, 76B20

published or not. The documents may come

Introduction

In this paper, we are interested in optimizing the hydrodynamic properties of a ship. In this regard, the most famous feature for minimizing the water resistance to the motion of a ship is perhaps the bulbous bow, a protruding bulb at the front of the ship, just beneath the waterline [START_REF] Ferreiro | The social history of the bulbous bow[END_REF]. Schematically speaking, the wave created by the bulb cancels the wave created by the bow of the ship; this reduces the wake and the energy lost in creating it. However, a bulbous bow is generally optimized for a cruising speed and it is generally not optimal for other speeds. Here, we want to find a hull which is optimal for a range of speeds.

The starting point in our approach is a simplified model in which the water resistance is the sum of a wave resistance (related to the wake of the ship) and of a viscous resistance (related to friction between the water and the hull of the ship). Such a decomposition is standard in ship hydrodynamics [START_REF] Birk | Fundamentals of ship hydrodynamics: Fluid mechanics, ship resistance and propulsion[END_REF], but we focus here on a specific model where, for a given form of the hull and for a given speed, the wave resistance is computed by Michell's formula [START_REF] Michell | The wave resistance of a ship[END_REF] and the viscous resistance is proportional to the wetted surface of the hull and to the square of the speed.

In Michell's formula, which is based on a linear potential flow theory, a function represents half of the wetted hull and the other half is obtained by symmetry. In [START_REF] Dambrine | A theoretical and numerical determination of optimal ship forms based on Michell's wave resistance[END_REF], the problem of finding a half hull function which minimizes the total resistance for a given speed of the ship and a given volume of the hull was solved both theoretically and numerically. In this case, the domain of definition of the hull function was also given (for instance, a rectangle, as in the Wigley hulls [START_REF] Michelsen | Wave resistance solution of Michell's integral for polynomial ship forms[END_REF]). This unified former results from the literature. In particular, the existence and uniqueness result of Krein and Sizov [START_REF] Krein | On the form of a ship of minimum total resistance[END_REF][START_REF] Sizov ; Professor | The seminar on ship hydrodynamics[END_REF] was extended to a Sobolev setting and the famous bulbous bow was recovered numerically for moderate values of the speed, as in [START_REF] Hsiung | Optimal ship forms for minimum total resistance[END_REF][START_REF] Lian-En | Optimal ship forms for minimal total resistance in shallow water[END_REF].

In [START_REF] Dambrine | Regularity of optimal ship forms based on Michell's wave resistance[END_REF][START_REF] Dambrine | Continuity with respect to the speed for optimal ship forms based on Michell's formula[END_REF], the more modern question of optimizing also the domain of definition of the hull function, for a given area, was considered. An optimal support was proved to exist, assuming for compactness that all the admissible supports belonged to a bounded hold-all domain. Several numerical simulations showed again the presence of a bulbous bow in a specific speed regime.

The simulations above pointed out that the shape of the bulb is very sensitive to the speed of the ship. Our purpose here is to propose a more robust model in which the variations or incertainties in the cruising speed are taken into account. Thus, we assume that the speed is a random variable whose probability distribution is known and we minimize the expectation of the total resistance for a given volume of the hull. We shall use a deterministic algorithm to compute the minimizer. In control theory, this approach is known as a robust deterministic control (see, e.g., [2, p.469, case (d)]). As previously, we will consider two situations. First, we handle the case where the domain of definition of the hull is fixed, and secondly we optimize this domain of definition, for a given area.

We first explain our model based on Michell's formula and on a given probability distribution of the speed in Section 2. Then, in Section 3, we consider the situation where the domain of definition of the hull function is fixed. We prove the existence and uniqueness of an optimal hull with given volume. We also establish the regularity of this optimal hull when the domain is a rectangle. Numerical simulations are given for two different kinds of probability distributions. Because of the averaging process, the bulbous bow is no longer present in the robust hull.

In Section 5, we optimize the support of the optimal hull for a given area. We prove the existence of an optimal support and then we focus on numerical simulations. In order to reduce the computational time, we adapt to our problem Newton's method for shape optimization. This is a technical matter because of the complex form of second order shape derivatives, so we first explain our discretization on several model problems in Section 4.

In the context of shape optimization, Newton's method is scarce in the literature because of its complexity. It was first implemented in [START_REF] Novruzi | Newton's method in shape optimisation: a three-dimensional case[END_REF], where the authors give a discretization of the second order shape derivative as a full matrix. We also refer to [START_REF] Harbrecht | A Newton method for Bernoulli's free boundary problem in three dimensions[END_REF] where the author uses a conjugate gradient method to compute the inverse of the Hessian matrix at every iteration. More recently, a level-set approach was developed in [START_REF] Allaire | Second-order shape derivatives along normal trajectories, governed by Hamilton-Jacobi equations[END_REF][START_REF] Vie | Second-order derivatives for shape optimization with a level-set method[END_REF].

For our problem, we use a geometric shape optimization approach with a fitted mesh and a P 1 finite element discretization. The use of a trapezoidal rule for the numerical integration on the boundary allows us to deal only with a diagonal Hessian matrix. This diagonal matrix reminds the method of lumped masses used for the finite element approximation of parabolic problems [START_REF] Thomée | Galerkin finite element methods for parabolic problems[END_REF]Chapter 15]. It can also be regarded as a quasi-Newton method. This is detailed in Section 4 for the Dirichlet energy. A similar discretization is used for our optimal design problem in Section 5. The numerical results show that the robust hull obtained by this approach still has a bulbous bow, even if the averaging process has a smoothing effect on the geometry of the hull.

The formal optimal design problems

The total resistance of water to the motion of a ship is the force required to tow the ship in calm water a constant speed. A traditional and simplified approach (see, e.g., [START_REF] Birk | Fundamentals of ship hydrodynamics: Fluid mechanics, ship resistance and propulsion[END_REF]) uses the decomposition R total = R wave + R viscous .

(2.1)

The wave resistance reflects the energy that goes into creating the wave.

In 1898, Michell [START_REF] Michell | The wave resistance of a ship[END_REF] gave a formula for the wave resistance, valid for any shape of the hull, for a constant speed and in an infinite domain. Experiments starting with Wigley in the 1920's showed a reasonable good agreement between theory and experiment (see the review by Gotman [START_REF] Gotman | Study of Michell's integral and influence of viscosity and ship hull form on wave resistance[END_REF] and references therein). We first recall Michell's formula and then we give a simple formula for the viscous resistance. 2.1. Michell's wave resistance formula. Consider a ship moving with constant speed U on the surface of an unbounded fluid. We assume that the coordinates xyz are fixed to the ship: the xy-plane is the water surface and z is vertically downward. The immerged half hull surface is represented by a continuous nonnegative function (see Figure 1)

y = f (x, z) ≥ 0, (x, z) ∈ D.
The set D on which f is defined is split into three parts: its interior D (an open subset of the upper half-plane with a Lipschitz boundary), in which f (x, z) > 0, its boundary Γ N at the surface z = 0 and its boundary Γ 0 under the surface, on which f (x, z) = 0 (see Figure 2).

It is assumed that the fluid is incompressible, inviscid and that the flow is irrotational. The motion has lasted long enough so that a steady state has been reached.

Michell's formula [START_REF] Michell | The wave resistance of a ship[END_REF] reads

R wave = 4ρg 2 πU 2 ∞ 1 (I 1 (λ) 2 + I 2 (λ) 2 ) λ 2 √ λ 2 -1 dλ, (2.2)
x z y 0

Γ 0 U D Γ Ν Figure 2.
The domain of parameters (x, z)

with

I 1 (λ) = D ∂f (x, z) ∂x exp - λ 2 gz U 2 cos λgx U 2 dxdz, (2.3)
I 2 (λ) = D ∂f (x, z) ∂x exp - λ 2 gz U 2 sin λgx U 2 dxdz. (2.4)
In this formula, U (in m • s -1 ) is the speed of the ship, ρ (in kg • m -3 ) is the constant density of the fluid and g (in m • s -2 ) is the standard gravity. The variables x, z and f (x, z) are expressed in meters. Consequently, the integrals I 1 (λ) and I 2 (λ) are in m 2 and R wave (in Newton) has the dimension of a force. The integration parameter λ has no dimension: it can be interpreted as λ = 1/ cos θ, where θ is the angle between the ship's path and the direction of wave propagation [5, p. 310], [START_REF] Tuck | Sea wave pattern evaluation -Part 1 report: Primary code and test results (surface vessels)[END_REF]. For θ = 0 (λ = 1), waves follow the ship with their crests aligned perpendicular to the ship's course (transverse waves) and for θ = ±π/2 (λ = +∞), waves have crests parallel to the ship's path.

We stress that Michell's approach is a linear theory in which the ship is assumed to be "thin", which means that the angles made by the hull surface with the longitudinal plane of symmetry are small, i.e.

0 ≤ f < < 1, ∂f ∂x < < 1 and ∂f ∂x < < 1 in D. (2.5)
In our approach, we do not assume the conditions (2.5) a priori, but we will recover them in a weak form by penalizing these constraints thanks to the viscous resistance.

Remark 2.1. In Michell's model, the flow is irrotational and incompressible, so we seek a potential in the form -U x + φ where φ is small and satisfies

       -∆φ = 0 in R × R + × R + (a quarter space), ∂ xx φ - g U 2 ∂ z φ = 0 on z = 0, ∂ y φ + U ∂ x f = 0 on y = 0 + . (2.6)
In (2.6), the first equation is the incompressibility condition whereas the third equation is a linearization of the impermeability condition on the hull. The second equation is known as the Neumann-Kelvin condition [START_REF] Kuznetsov | Linear water waves[END_REF]. It is obtained by combining a linearized no-slip condition and a linearized Bernoulli equation at the free surface.

In [START_REF] Michell | The wave resistance of a ship[END_REF], Michell managed to compute a physical solution to problem (2.6) which satisfies the condition lim x→+∞ |∇φ| = 0. This means that in front of the ship, there is no wake. The wave resistance is the drag force in this model and it is computed as

-2 R 2 δp∂ x f (x, z)dxdz,
where δp = ρU ∂ x φ(x, 0, z) stands for the difference of pressure due to the ship.

In [START_REF] Sretensky | Sur la détermination de la résistance ondulatoire d'un navire se déplaçant à la surface de l'eau d'une profondeur finie[END_REF], Sretensky solved problem (2.6) in the case of a finite depth H by means of a Fourier transform. He recovered Michell's wave resistance formula by letting H tend to infinity.

2.2.

The viscous resistance and the total resistance. In formula (2.1), a traditionnal approach is to express the viscous resistance as (see, e.g., [START_REF] Birk | Fundamentals of ship hydrodynamics: Fluid mechanics, ship resistance and propulsion[END_REF])

R viscous = 1 2 ρU 2 C F A . (2.7)
The viscous resistance R viscous accounts for the effects of viscosity which are not present in Michell's model. In (2.7), C F is the viscous drag (or friction) dimensionless coefficient and A (in m 2 ) is the surface area of the ship's wetted hull. The coefficient C F can sometimes have a complicated form [START_REF] Birk | Fundamentals of ship hydrodynamics: Fluid mechanics, ship resistance and propulsion[END_REF] but throughout this manuscript, we assume that C F is a positive constant. Since the graph of f represents the ship's half-hull, A is given by:

A = 2 D 1 + |∇f (x, z)| 2 dxdz. (2.8)
For a thin ship, |∇f | is uniformly small (see (2.5)), and the integral above can be approximated by performing a Taylor expansion of 1 + |∇f | 2 at first order:

A = 2 D 1dxdz + D |∇f (x, z)| 2 dxdz + o(||∇f || 2 ∞ ) . (2.9)
A good approximation of the viscous drag for small ∇f reads

R viscous = 1 2 ρU 2 C F 2|D| + D |∇f (x, z)| 2 dxdz , (2.10)
where |D| is the area of D. Summing up, the total resistance for a ship hull defined by f reads

R total = R wave + R viscous , (2.11)
where R wave is defined by (2.2)-(2.4) and R viscous is defined by (2.10). In this formula, the cruising speed of the ship, U , is constant.

2.3.

The formal optimal design problems. The variables g and ρ are known physical constants and we assume for simplicity that the viscous drag coefficient C F is a given constant. Then the total resistance defined by (2.11), (2.2) and (2.10) depends only on the function f : D → R and on the speed U , so we write R total (f, U ). The half volume of the hull (a positive real number) is given and it is denoted by V .

2.3.1.

Robust optimization of the hull for a given domain. In [START_REF] Dambrine | A theoretical and numerical determination of optimal ship forms based on Michell's wave resistance[END_REF], Dambrine, Pierre and Rousseaux solved the following convex problem for a fixed domain D, by using an appropriate H 1 functional setting: for a given speed U > 0,

   Find the function f D which minimizes R total (f, U ) in the set f : D → R, f = 0 on Γ 0 and D f (x, z)dxdz = V . (2.
12)

The solution to this problem may be very sensitive to variations of U . In order to have a more robust problem, we assume now that U is a random function with realizations denoted by U ω and we consider the problem:

   Find the function f D which minimizes E(R total (f, U )) in the set f : D → R, f = 0 on Γ 0 and D f (x, z)dxdz = V . (2.
13)

It will be interesting to compare R total (f D , U ) with R total (f D , U ). The cost function R total (f, U ) can more generally be replaced by

J(f, U ) = h(U ) × R total (f, U )
where h : (0, +∞) → (0, +∞) is a continuous function of U . For instance, we may choose

• h(U ) = U T where T is a time, in which case J(f, U ) is an energy; • h(U ) = 1, in which case J(f, U ) is the total resistance.
Remark 2.2. Since Michell's wave resistance is computed for a steady state, this model means that in the term E(h(U )R total (f, U )), we neglect the transitory states between two steady states.

Remark 2.3. It would be more natural to add the condition f ≥ 0 in the set of functions for problem (2.12), in order to avoid self-crossing of the hull, as in [START_REF] Dambrine | A theoretical and numerical determination of optimal ship forms based on Michell's wave resistance[END_REF]. Here, we will check this condition numerically. That is, we first solve problem (2.12), a quadratic-linear minimization problem which reduces to a linear problem (cf. (3.16)), and we check a posteriori that f D ≥ 0 in D. The same approach holds for problem (2.13).

2.3.2.

Robust optimization of the domain and of the hull. Next, we consider the domain D as a variable of the problem. The total resistance is written R total (D, f D , U ) where for a given domain D and a given speed U > 0, f D is the solution to problem (2.12). We fix an area a > 0. In [START_REF] Dambrine | Regularity of optimal ship forms based on Michell's wave resistance[END_REF][START_REF] Dambrine | Continuity with respect to the speed for optimal ship forms based on Michell's formula[END_REF], Dambrine and Pierre studied the following shape optimization problem: for a given constant speed U , (2.15)

    
For the numerical resolution of problem 2.15 in Section 5.2, we will compute the first order and second order shape derivatives of the set functional. For its theoretical resolution in Section 5.1, we will seek for an optimal support rather than an optimal domain.

3. Robust optimal hull for a fixed domain 3.1. Functional setting. Let D be a bounded and connected open subset of the upper half plane {(x, z) ∈ R 2 : z > 0}. We assume that D has a Lipschitz boundary ∂D. For the sake of simplicity, we also assume throughout section 3 that the part of ∂D which intersects the x-axis, namely ∂D ∩ (R × {0}), is a segment of the xaxis (possibly empty). We let Γ N denote the relative interior of this segment and Γ 0 = ∂D \ Γ N .

We work with the Sobolev space

H(D) = u ∈ H 1 (D) : u = 0 on Γ 0 in the sense of traces .
Let u ∈ H(D) (u is the hull function) and V > 0 (V is the speed of the ship). For Michell's wave resistance, we set

α = g V 2 (3.1)
and (2.2) becomes

R wave (u, α) = 4ρgα π ∞ 1 |S u (λ)| 2 λ 2 √ λ 2 -1 dλ, with (cf. (2.3)-(2.4)) S u (λ, u, α) = D ∂u ∂x (x, z)e -iλαx e -λ 2 αz dxdz.
The number α is known as the Kelvin wave number. The value 1/α (in m) is the typical wavelength of the transverse waves. Transverse waves follow the ship with their crests and troughs aligned perpendicular to the ship's course and the wavelength is the distance between two successive crests. Since the speed V of the ship can be recovered from α through V = g/α, knowing α is equivalent to knowing V . We will use the variable α because it is more convenient.

Integrating S u by parts with respect to x and taking advantage of the boundary condition u = 0 on ∂Ω, we find that

R wave (u, α) = 4ρgα 3 π ∞ 1 |T (u, α, λ)| 2 λ 4 √ λ 2 -1 dλ, (3.2)
where, for all λ > 0,

T (u, α, λ) = D u(x, z)e -iλαx e -λ 2 αz dxdz.
The viscous resistance (2.10) reads

R viscous (u, α) = ρg 2α C F 2|D| + D |∇u(x, z)| 2 dxdz .
The total resistance is

R total (u, α) = R viscous (u, α) + R wave (u, α). (3.3)
We assume that the Kelvin wave number is a random variable ν : Ω → R on a complete probability space (Ω, A, P). Moreover, we assume that the probability distribution P ν of ν is a measure which has a compact support in (0, +∞). We recall that the probability distribution of ν is defined for every Lebesgue-measurable set B ⊂ R by P ν (B) = P(ν ∈ B) = P ({w ∈ Ω : ν(ω) ∈ B}) . Concerning the Kelvin wave number, ν denotes the random variable, whereas α is a positive real number.

We recall that if X : Ω → R is a random variable which is integrable with respect to P, its expectation is defined by

E(X) = Ω X(ω)dP(ω).
If ϕ : (0, +∞) → R is a continuous function, then ϕ(ν) : Ω → R is a random variable and we have the well-known formula,

E[ϕ(ν)] = Ω ϕ(ν(ω))dP(ω) = R ϕ(s)dP ν (s).
In view of (3.3), we consider the cost function

J D (u) = E h ρ,g,C F (ν) R viscous (u, ν) + R wave (u, ν) , (3.4)
where h ρ,g,C F : (0, +∞) → (0, +∞) is a continuous function which may depend on the physical constants ρ, g and on the viscous coefficient C F . We note that it could be interesting to consider C F as a random variable as well (possibly dependent on α), but we assume for simplicity that C F is constant. Following [START_REF] Dambrine | Continuity with respect to the speed for optimal ship forms based on Michell's formula[END_REF], we introduce the normalized viscous resistance functional

J 0 (u) = D |∇u| 2 dxdz
and the normalized wave resistance functional

J wave (u, α) = 8α 4 π ∞ 1 |T (u, α, λ)| 2 λ 4 √ λ 2 -1 dλ, (3.5)
where T is defined as previously for all u ∈ H(D), α > 0 and λ > 0 by

T (u, α, λ) = D u(x, z)e -iλαx e -λ 2 αz dxdz. (3.6)
The cost function (3.4) reads

J D (u) = E h ρ,g,C F (ν) ρg 2ν C F 2|D| + J 0 (u) + 1 C F J wave (u, ν) . By setting h(α) = h ρ,g,C F (α) ρg 2α C F , (3.7)
we have

J D (u) = E h(ν) 2|D| + J 0 (u) + 1 C F J wave (u, ν) . (3.8)
We note that h generally depends on the fixed parameters ρ, g and C F . We let V > 0 denote the half-volume of the hull. We consider the set

H V (D) = u ∈ H(D) : D u dxdz = V ,
which is a closed affine subspace of H(D). Our robust optimization problem reads:

Find u D ∈ H V (D) such that J D (u D ) = min {J D (u) : u ∈ H V (D)} . (3.9)
In (3.9), the set D is given and the positive parameters ρ, g, C F and V are fixed.

Theoretical results.

The following result will prove useful.

Lemma 3.1. Let q ∈ (1, +∞) and let q = q/(q -1) ∈ (1, +∞) denote the conjugate exponent of q. Assume that H ∈ L q (D × D). Then for all u, v ∈ L q (D), we have

D×D |H(x, z, x , z )u(x, z)v(x , z )|dxdzdx dz ≤ H L q (D×D) u L q (D) v L q (D) .
(3.10) Moreover, for each u ∈ L q (D), the function

(x, z) → D H(x, z, x , z )u(x , z )dx dz (3.11)
belongs to L q (D).

Proof. Estimate (3.10) is a consequence of Hölder's inequality. The claim on the function (3.11) follows from (3.10), Fubini's theorem and a duality argument.

By formally switching the integrals in the expression (3.5)-(3.6), we see that Michell's normalized wave resistance can be written

J wave (u, α) = D×D k α (x, z, x , z )u(x, z)u(x , z )dxdzdx dz (3.12)
where

k α (x, z, x , z ) = 8α 4 π K(α(x -x ), α(z + z )) (3.13)
and

K(X, Z) = ∞ 1 e -λ 2 Z cos(λX) λ 4 √ λ 2 -1 dλ. (3.14)
This formal calculation was rigorously proved in [8, Appendix A]. It was shown that Michell's kernel belongs to L 5/4-ε (D × D) and that this estimate is optimal if D contains a half-disk centered on the x-axis.

The results from [8, Appendix A] are summarized in the proposition below. We first note that K is defined and continuous on R×(0, +∞), thanks to the exponential term, so that k α is continuous on (R × (0, +∞)) 2 . Proposition 3.2. Let α > 0. Michell's normalized wave resistance kernel k α (3.13) belongs to L q (D × D) for all 1 ≤ q < 5/4. For each q > 5 and for each u ∈ L q (D), the formulations for J wave (u) given by (3.5)-(3.6) and (3.12)-(3.13)-(3.14) are equal.

Let q ∈ (1, 5/4) and let q = q/(q -1) ∈ (5, +∞) be the conjugate exponent of q. Since H 1 (D) is continuously imbedded in L q (D) for all q ∈ [1 + ∞) [START_REF] Adams | Sobolev spaces[END_REF], Proposition 3.2 and Lemma 3.1 show that for all u ∈ H(D), J wave (u) < +∞.

The following result is proved in [9, Lemma 5.2 (ii)].

Lemma 3.3. For each u ∈ H(D), the map

α → J wave (u, α)
is continuous on (0, +∞).

Recall that P ν has a compact support in (0, +∞), so we may assume that its support is included in [ν min , ν max ] with 0 < ν min ≤ ν max < +∞. We define the kernel

H ν (x, z, x , z ) = νmax ν min h(s)k s (x, z, x , z )dP ν (s),
which has finite values for all (x, z, x , z

) ∈ D × D (since z + z > 0). Proposition 3.4.
The kernel H ν belongs to L q (D × D) for all 1 ≤ q < 5/4 and for all u ∈ H(D), we have

E h(ν)J wave (u, ν) = D×D H ν (x, z, x , z )u(x, z)u(x , z )dxdzdx dz . (3.15)
Proof. We choose q ∈ (1, 5/4) and we denote by q ∈ (5, +∞) the conjugate exponent of q. Since the function t → t q is convex on [0, +∞), by Jensen's inequality, we have

H ν (x, z, x , z ) q ≤ νmax ν min hq (s)|k s (x, z, x , z )| q dP ν (s).
Thus, by (3.13) and Fubini's theorem,

D×D H ν (x, z, x , z ) q dxdzdx dz ≤ νmax ν min D×D ĥ(s)|K(s(x -x ), s(z + z ))| q dxdzdx dz dP ν (s),
where ĥ : (0, +∞) → (0, +∞) is a continuous function, namely ĥ(s) = hq (s) 8s 4 π q .

We perform the change of variable (x, z, x , z ) = (sx, sz, sx , sz ) and we find

D×D H ν (x, z, x , z ) q dxdzdx dz ≤ D×D |K(x -x , z + z )| q dxdzdx dz νmax ν min 1 s 4 ĥ(s)dP ν (s),
In the right handside above, the first integral is finite thanks to Proposition 3.2. The second integral is finite since [ν min , ν max ] is compactly embedded in (0, +∞). This proves that H ν belongs to L q (D × D). For the computation of

E h(ν)J wave (u, ν) ,
we use the expression (3.12). The estimates above combined with Hölder's inequality (as in (3.10)) show that we may apply apply Fubini's theorem. This yields (3.15).

Remark 3.1. Let α 0 > 0 and assume that P ν is the Dirac delta function at α 0 , δ α 0 .

Then the kernel H ν is simply

H ν (x, z, x , z ) = h(α 0 )k α 0 (x, z, x , z ).
We are in position to prove:

Theorem 3.5. Problem (3.9) has a unique solution u D which is also the unique solution in H 1 (D) of the boundary value problem

                   -∆u(x, z) + 1 C F E h(ν) D×D H ν (x, z, x , z )u(x , z )dx dz = C, ∀(x, z) ∈ D, D u dxdz = V , u = 0 on Γ 0 , ∂u ∂n = 0 on Γ N . (3.16)
Proof. The functional T defined by (3.6) depends linearly on its first argument u, so that for each α > 0, u → J wave (u, α) is a quadratic and convex functional on H(D). Thus, u → E( h(ν)J wave (u, ν)) defined by (3.5) is convex on H(D). On the other hand, the function J 0 is stricly convex on H(D), thanks to the homogeneous Dirichlet boundary condition on Γ 0 and the Poincaré inequality. Thus, J D is stricly convex on H(D) and since H V (D) is an affine subspace of H(D), problem (3.9) has a most one solution in H V (D). Let (u n ) be a minimizing sequence for problem (3.9) in H V (D). We have

J D (u n ) = E h(ν)2|D| + E h(ν) J 0 (u n ) + 1 C F E h(ν)J wave (u n , ν) .
Each one of the three terms in the right hand-side is nonnegative so J D (u n ) is nonnegative and the sequence (u n ) is bounded in H 1 (D) (thanks to the term J 0 (u n ) and the Poincaré inequality). Up to a subsequence, (u n ) converges weakly in H 1 (D) to some u D , which belongs to H V (D) since the latter is a closed convex subset of H 1 (D). By lower semi-continuity of J 0 ,

J 0 (u D ) ≤ lim inf n J 0 (u n ). Now, let q ∈ (5, +∞). The space H 1 (D) is compactly embedded in L q (D) [1] so (u n ) converges strongly to u D in L q (D). By Proposition 3.4, E h(ν)J wave (u n , ν) → E h(ν)J wave (u D , ν) . Thus, J D (u D ) ≤ lim inf n J D (u n ). since (u n
) is a minimizing sequence, this proves that u D is a minimizer. The Euler-Lagrange equation associated to problem (3.9) yields the boundary value problem (3.16). The constant C in the first line is the Lagrange multiplier associated to the volume constraint D udxdz = V . Conversely, each solution to (3.16) is a critical point of J D on H V (D) and by convexity, it is a minimizer. This concludes the proof. Remark 3.2. In particular, the statement of Theorem 3.5 is true if P ν = δ α 0 for some α 0 > 0. This means that the ship has a constant speed. In this case, H ν = h(α 0 )k α 0 (see Remark 3.1). The same holds for Proposition 3.6 and Theorem 3.7 below. Proposition 3.6. If D is symmetric with respect to z-axis, then u D is even with respect to x.

Proof. Since D is symmetric with respect to the z-axis, for all (x, z) ∈ D, we have (-x, z) ∈ D. Let ǔ ∈ H(D) be defined by

ǔ(x, z) = u D (-x, z), ∀(x, z) ∈ D.
Performing the change of variable x → -x in T (ǔ, α, λ) (see (3.6)), we find that

T (ǔ, α, λ) = D u D (x, z)e iλαx e -λ 2 αz dxdz. Thus, |T (ǔ, α, λ)| = |T (u D , α, λ)| and consequently, by (3.5), J wave (ǔ, α) = J wave (u D , α), ∀α > 0.
Similarly, we have J 0 (ǔ) = J 0 (u D ) and so

J D (ǔ) = J D (u D ).
Since ǔ belongs to H V (D), this shows that ǔ is a solution to problem (3.9). By uniqueness of the solution, ǔ = u D . Theorem 3.7. If D is a rectangle, then the solution u D to problem (3.9) belongs to W 2,5/4-ε (D) for all ε > 0 small enough.

Proof. Let q ∈ (1, 5/4). Since u D belongs to H 1 (D) which is continously embedded in L q (D), we deduce from (3.16), Proposition 3.4 and Lemma 3.1 that u D solves the PDE -∆u D + f = C in D, where f belongs to L q (D) and C is constant. The domain D is a rectangle and u D satisfies homogeneous Dirichlet boundary conditions on three sides and homogeneous Neumann boundary conditions on one side. By elliptic regularity on polygons [START_REF] Grisvard | Elliptic problems in nonsmooth domains[END_REF], u D belongs to W 2,q (D).

Numerical simulations.

In this section, we present numerical results for the optimal hulls which minimize the expectation of the total resistance (3.3).

3.3.1.

The two optimal design problems. The cost function is

J D (u) = E (R total (u, ν)) = E (R viscous (u, ν) + R wave (u, ν)) .
In other words, we choose h ρ,g,C F (α) = 1 or equivalently (cf. (3.7))

h(α) = ρg 2α C F (3.17)
in the cost function J D defined by (3.8). The domain D is a rectangle whose upper side is on the x-axis. We recall that the robust optimization problem (3.9) reads

Find u D ∈ H V (D) such that J D (u D ) = min {J D (u) : u ∈ H V (D)} . (3.18)
We shall compare the optimal hull u D with the solution to the following (non robust) optimization problem:

Find ūD ∈ H V (D) such that R total (ū D , α) = min {R total (u, α) : u ∈ H V (D)} , (3.19
) where the value α is set to α = E(ν). By Remark 3.2, problem (3.19) has a unique solution. We shall also compare these optimal hulls with standard Wigley hulls (see, e.g., [START_REF] Michelsen | Wave resistance solution of Michell's integral for polynomial ship forms[END_REF]).

Problem (3.19) was investigated numerically in [START_REF] Dambrine | A theoretical and numerical determination of optimal ship forms based on Michell's wave resistance[END_REF] (see also references therein). We have also computed the optimal hull ūD for different values of α. The domain is a rectangle with length L = 2.2 m and draft T = 0.3 m. The half volume of the hull is V = 0.06 m 3 . The other parameters are

ρ = 1000 kg • m -3 , g = 9.81 m • s -1 and C F = 0.01. (3.20)
For the space discretization of the problem, we used a P 1 finite element approach (see Section 3.3.2). A triangulation of the rectangular domain into 1956 triangles and 1075 vertices was used.

It is convenient to introduce the length Froude number

F r L = V √ gL = 1 √ αL , (3.21)
which is a dimensionless version of the speed V . With these parameters, we have recovered the results from [START_REF] Dambrine | A theoretical and numerical determination of optimal ship forms based on Michell's wave resistance[END_REF]. Namely, for intermediate Froude numbers F r L ∈ [0.5, 1], the hull has the famous bulbous bow, which reduces the contribution of the wave resistance in the total resistance. A bulbous bow is seen in Figure 5 (bottom) for F r L = 0.6.

For large or small Froude numbers (F r L > 1 or F r L < 0.5), the influence of Michell's wave resistance is small and the optimal hull ūD mainly minimizes the viscous resistance. Figure 3 (bottom) shows ūD for F r L = 0.277.

Numerical approximation.

The expectation of the water resistance J D (u) is given by

J D (u) = E h(ν) 2|D| + J 0 (u) + 1 C F J wave (u, ν) = 2E[ h(ν)]|D| + E[ h(ν)]J 0 (u) + 1 C F E[ h(ν)J wave (u, ν)].
By dropping the constant term, we see that minimizing J (u) is the same as minimizing To solve this problem, we adopt a finite element approach, in the sense that the optimal hull u D is sought in a finite dimensional space

JD (u) = D |∇u| 2 dxdz + 1 C F E[ h(ν)] D×D H ν (x, z, x , z )u(x, z)u(x , z )dxdzdx dz .
V h ⊂ H(D) ⊂ H 1 (D).
Given a basis of P 1 (continuous and piecewise linear) finite elements {φ 1 , φ 2 , . . . , φ n }, we can write

for all u h ∈ V h , u h = n i=1 u i φ i . (3.23)
Then we have the following discretization of JD (u h ):

JD (u h ) = U t M 0 + 1 C F E[ h(ν)] M w U, (3.24)
where U = (u 1 , . . . , u n ) t is the vector of the coordinates of u h and M 0 , M w are matrices approximating J 0 and u h → E[ h(ν)J wave (u h , ν)]. This yields the following discrete optimization problem:

u D = argmin U t M 0 + 1 C F E[ h(ν)] M w U : u h ∈ V h , D u h dxdz = V . (3.25)
The computation of the stiffness matrix

M 0 = D ∇φ i • ∇φ j dxdz 1≤i,j≤n
is standard. This matrix M 0 is symmetric and positive definite. It is also sparse and non-diagonal.

The matrix M w is symmetric and positive semi-definite, but in contrast to M 0 , M w is a full matrix because it is related to a nonlocal operator. Let us explain briefly how M w is computed. Before discretizing, we write

E[ h(ν)J wave (u, ν)] = E[R wave (u, ν)] = 4ρg π E ν 3 ∞ 1 |T (u, ν, λ)| 2 λ 4 √ λ 2 -1 dλ = 4ρg π νmax ν min α 3 f (α) ∞ 1 |T (u, α, λ)| 2 λ 4 √ λ 2 -1 dλdα
where T (u, α, λ) is given by (3.6) and f is the probability density function of ν on the interval [ν min , ν max ] ⊂ (0, +∞). The approximation of the integral is done as follows. For a given (α, λ), the integral

T (u h , α, λ) = D u h (x, z) exp(-iλαx -λ 2 αz)dxdz
is computed by exact integration over each triangle of the mesh which approximates the domain D. The integration over α is handled by a numerical integration (the trapezoidal rule).

Concerning the density f , we consider two situations:

• A continuous uniform probability distribution for the Kelvin number on [ν min , ν max ], in which case

f (α) = 1 ν max -ν min on [ν min , ν max ]; (3.26)
• A continuous uniform probability distribution for the speed V = g/α, as in (3.30).

For a given α > 0 and u h ∈ V h , the approximation of the integral

∞ 1 |T (u h , α, λ)| 2 λ 4 √ λ 2 -1 dλ
is the same as in [START_REF] Dambrine | A theoretical and numerical determination of optimal ship forms based on Michell's wave resistance[END_REF]. In particular, we use the value

2 1 1 √ λ 2 -1 = ln(2 + √ 3)
as in [START_REF] Tarafder | Computation of wave-making resistance of Wigley hull form using Michell's integral[END_REF] to remove the singularity at λ = 1. The integral is also truncated at +∞ and we use an appropriate midpoint formula for the remainder of the integral. has a compact support included in [ν min , ν max ] with 0 < ν min < ν max < +∞. The bounds ν min , ν max are chosen so that the corresponding length Froude numbers are F r L min = 0.2 and F r Lmax = 1. Thus, we have

ν min = 1 F r 2 Lmax L = 0.45 and ν max = 1 F r 2 L min L = 11.36.
The case of a continuous uniform distribution of the Kelvin wave number.

In this case, the density f is given by (3.26). The mean value of ν is

E[ν] = (ν max + ν min )/2 = 5.90. (3.27)
This corresponds to a Froude number equal to F r L = 0.277. Since F r L is a decreasing function of α (see (3.21)), the distribution of Froude numbers is concentrated near 0.2. Figure 3 gives the shapes of the computed optimal hulls u D (middle) and ūD (bottom). The top figure represents a Wigley hull w with triangular cross section and parabolic horizontal section, which is defined by [START_REF] Michelsen | Wave resistance solution of Michell's integral for polynomial ship forms[END_REF][START_REF] Tarafder | Computation of wave-making resistance of Wigley hull form using Michell's integral[END_REF]

f (x, z) = B 2 1 - z T 1 - 4x 2 L 2 , (x, z) ∈ [-L/2, L/2] × [0, T ]. (3.28)
The beam B of w is chosen such that

D w(x, z)dxdz = BLT 6 = V . (3.29)
Namely, B = 6V /(LT ) = 0.5454 m.

The three hulls have a very similar shape. We also note that the optimal hulls are symmetric back and front, in agreement with Proposition 3.6. However, they perform quite differently. This can be seen in Figure 4 which represents, for each hull, the total resistance (in Newton) on the range of Froude numbers [0.2, 1.0]. We stress that the constant term 2 h(α)|D| is not represented in this figure (see (3.22)). For each hull, the expectation of the total resistance is also given.

As expected, u D clearly has the best expectation on the total range [0.2, 1.0] (24 N for u D , 33 N for the Wigley hull w and 36 N for ūD ).

In contrast, ūD has the minimal total resistance for the mean Froude number F r L = 0.277. The hull ūD is even optimal among the three shapes for F r L ∈ [0.2, 0.31]. This is mainly due to a small contribution of the wave resistance in this interval. Interestingly, the expectation of the Wigley hull w is slightly better than the expectation of ūD : a bulbous bow is interesting only in a limited range of speeds. The mean value of the Froude number is obviously F r L = 0.6.

By an appropriate change of variables, we see that ν = g/V 2 has a probability density given by

f (α) = g/α 2α(V max -V min ) on [ν min , ν max ]. (3.30)
Figure 5 shows the computed optimal hulls u D , ūD and the Wigley hull w. There is a striking difference between ūD , which has a bulbous bow, and u D , which has no bulbous bow and looks like a Wigley hull.

The hydrodynamic properties of each hull appear clearly in Figure 6 which represents, for each hull, the total resistance on the range of Froude numbers [0.2, 1.0] (up to the constant term 2 h(α)|D|). The expectation of the total resistance is best for u D (83 N), as expected. It is much better than for the Wigley hull w (128 N), but it is surprisingly very close to the expectation of ūD (84 N). For each hull, the total resistance is the sum of the viscous resistance and of the wave resistance. Since the viscous resistance is proportional to V 2 (cf. (2.10)), it is interesting to introduce for a given hull u the ratio

C viscous (u) = 1 F r 2 L R viscous (u, α) -2 h(α)|D| ,
which does not depend on α.

We have

C viscous (u D ) = 25 N, C viscous (w) = 33 N and C viscous (ū D ) = 62 N.
Thus, ūD has the worst viscous resistance for every speed: a bulbous bow is clearly not optimal at high speed (F r L close to 1). At the mean speed F r L = 0.6, ūD has a better total resistance than u D and w because it has a much better wave resistance. For F r L ∈ [0.2, 0.4], u D has a better total resistance than ūD and this is due to the wave resistance since the difference in the viscous resistance between the two hulls is less than (62 -25) × 0, 4 2 ≈ 6 N. For each hull, we notice the well-known oscillations of Michell's wave resistance for F r L < 0.4.

Newton's method for shape optimization

In this section, we first explain Newton's method for shape optimization in the continuous version. Then, we describe its discretization along with numerical exemples in the case of a PDE-fee example, before adapting it to the case of the Dirichlet energy. This is known as an optimize-then-discretize approach. 4.1. The continuous version. For k ≥ 1, W k,∞ (R 2 ; R 2 ) denotes the usual Sobolev space based on L ∞ (R 2 ; R 2 ). In particular, W 1,∞ (R 2 ; R 2 ) is the space of bounded Lipschitz continuous functions from R 2 into R 2 . In order to work with classical derivatives, we consider for k ≥ 1 the space

C k,∞ (R 2 ; R 2 ) = C k (R 2 ; R 2 ) ∩ W k,∞ (R 2 ; R 2 ), endowed with the norm of W k,∞ (R 2 ; R 2 ). For k ∈ N, O k is
E : C k,∞ (R 2 ; R 2 ) → R θ → E ((Id + θ)(Ω)) .
The function E is said to be shape-differentiable at Ω if E is Fréchet-differentiable at 0, that is, if there exists a continuous linear map E (0;

•) : C k,∞ (R 2 ; R 2 ) → R such that E(θ) = E(0) + E (0; θ) + o θ C k,∞ (R 2 ;R 2 ) .
We denote E (Ω; θ) = E (0; θ).

The second order shape derivative is defined as follows [37, Definition 3.2.12]. Definition 4.2. The function E of Definition 4.1 is said to be twice shape-differentiable at Ω if E is Fréchet-differentiable in a neighborhood U of 0 in C k,∞ (R 2 ; R 2 ) and if the first derivative E defined by

E : U → C k,∞ (R 2 ; R 2 ) θ → E (θ; •),
is Fréchet-differentiable at 0. We denote by E (0; θ, ξ) the second Fréchet derivative at 0, θ and ξ being respectively the first and second directions of derivation. We also denote E (Ω; θ, ξ) := E (0; θ, ξ).

A set Ω ∈ O k is a critical shape if E (Ω; θ) = 0 for all θ ∈ C k,∞ (R 2 ; R 2 ). In particular, a minimizer of E in O k is a critical shape.
Formally, Newton's method for finding a critical shape reads as follows. Start from a set Ω 0 ∈ O k and for p = 0, 1, 2, . . . until convergence:

• Find θ p ∈ C k,∞ (R 2 ; R 2 ) such that E (Ω p ; θ p , ξ) = -E (Ω p ; ξ), ∀ξ ∈ C k,∞ (R 2 ; R 2 ); (4.1) • Set Ω p+1 = (Id + θ p )(Ω p ) (4.2)
If Ω ∈ O k+1 and if E : O k → R is shape differentiable at Ω, it is well-known since Hadamard [START_REF] Hadamard | OEuvres de Jacques Hadamard. Tomes I[END_REF] that E (Ω; θ) depends only on the trace of the normal component of θ on the boundary Γ = ∂Ω (see, e.g., [19, Proposition 5.9.1]). Moreover, if Ω ∈ O k+2 , then E (Ω; •, •) is not invertible because each vector field θ ∈ C k,∞ (R 2 ; R 2 ) which has compact support in R 2 \ Γ belongs to the kernel of E (Ω; •, •). This is a consequence of the structure theorem of E (Ω; •, •) [START_REF] Novruzi | Structure of shape derivatives[END_REF].

Thus, we solve equation (4.1) for vector fields θ, ξ defined on Γ. Since a small perturbation of Γ is fully defined by a normal deformation, in Newton's equation (4.1) we will seek θ p ∈ C k,∞ such that its tangential component on Γ, denoted by θ Γ , is zero. Thus, (4.1) becomes

find θ ∈ C k,∞ (R 2 ; R 2 ) such that θ Γ = 0 and E (Ω; θ, ξ) = -E (Ω; ξ), ∀ξ ∈ C k,∞ (R 2 ; R 2 ) such that ξ Γ = 0, (4.3)
where we have omitted the subscript p for simplicity. This approach simplifies E (Ω; θ, ξ) because it cancels the tangential components of θ and ξ which are otherwise present if Ω is sufficiently regular [START_REF] Novruzi | Structure of shape derivatives[END_REF]. Moreover, this is a consistent approach. Indeed, the purpose is to converge to a critical shape Ω and for a critical shape, E (Ω ; θ, ξ) depends only on the normal components of the traces of θ, ξ on Γ [19, Remark p. 246]. A related idea was used with Newton's algorithm for shape optimization in [START_REF] Allaire | Second-order shape derivatives along normal trajectories, governed by Hamilton-Jacobi equations[END_REF] in the context of a level-set discretization.

Thus, in order to solve problem (4.3), we first seek the normal component θ • n on Γ. Then, we extend θ on Ω with a Hilbertian extension-regularization technique as described in [4, Section 5.2]. In our case, we used the inner product inherited from linear elasticity. 4.2. PDE-free example. Let f ∈ C 3 (R 2 ; R). We consider the minimization problem inf

Ω∈O 3 E(Ω) (4.4)
where

E(Ω) = Ω f (x)dx
is the function to be minimized. With appropriate assumptions on f (see, e.g., [37, Proposition 3.3.5]), an optimal shape Ω is explicitly given by the set of negative values of f , Ω = {x ∈ R 2 | f (x) < 0}. (4.5) We seek to compute Ω by Newton's method, and we shall follow an optimize-thendiscretize approach.

Assume that Ω ∈ O 3 . Since f is smooth, the function E from Definition 4.1 is of class C 2 on C 1,∞ (R 2 ; R 2 ) and the shape derivatives of E(Ω) are given by [START_REF] Henrot | Shape variation and optimization[END_REF]

E (Ω; ξ) = Γ (ξ • n)f, E (Ω; θ, ξ) = Γ (θ • n)(ξ • n) Hf + ∂ n f + Γ Z θ,ξ f, (4.6)
where Γ is the boundary of Ω, H is the mean curvature of Γ, and

Z θ,ξ = θ Γ • D Γ n ξ Γ -θ Γ • ∇ Γ (ξ • n) -ξ Γ • ∇ Γ (θ • n). (4.7)
For θ, ξ ∈ C 1,∞ (R 2 ; R 2 ) such that θ Γ = 0 and ξ Γ = 0, we have Z θ,ξ = 0. Thus, Newton's equation (4.3) reads: find

θ n ∈ C 1 (Γ) such that Γ θ n ξ n Hf + ∂ n f = - Γ ξ n f, ∀ξ n ∈ C 1 (Γ). (4.8)
This problem has to be properly discretized. We explain our approach which allows to obtain a diagonal hessian matrix when working with P 1 normal Lagrange finite element deformation vectors. 4.2.1. Discretization of Newton's equation. Given a conforming triangulation T h of a polygonal domain Ω h which approximates the domain Ω, Γ h denotes the boundary of the mesh which consists of an ordered n be -tuple of ordered vertices (x 1 , x 2 , . . . , x n be ). Let

V h = u h ∈ C(Ω h , R) ∀ T i ∈ T h , u h|T i ∈ P 1
be the Lagrange finite element space of continuous functions that are piecewise polynomials of degree 1. We denote by (φ i ) 1≤i≤N h the nodal basis of V h defined by

∀ 1 ≤ i, j ≤ N h , φ i (a j ) = δ ij ,
where N h is the dimension of V h and (a j ) 1≤j≤N h are the vertices of the triangulation T h . We denote by n h (x i ) the discrete normal vector to the boundary at the vertex x i which is approximated as the rotate of the tangent at

x i , τ h (x i ) = ------→ x i-1 x i+1 ------→ x i-1 x i+1
as shown in Figure 7. We search for a descent direction θ h : R 2 → R 2 defined on the vertices of T h . Since (4.8) contains boundary integrals, we first compute the descent direction θ h on the boundary vertices x 1 , . . . , x n be before extending it to all of the triangulation T h through Hilbertian extension techniques (see [START_REF] Allaire | Shape and topology optimization[END_REF]). We assume that θ h is normal to the boundary Γ h at every vertex x i , that is

∀i = 1, . . . , n be , ∃α i ∈ R, θ h (x i ) = α i n h (x i ).
(4.9)

We use the continuous piecewise linear (P 1 ) extension of θ h on Γ h , namely

θ h (x) = n be i=1 α i ϕ i (x)n h (x i ), (4.10)
where (ϕ i ) 1≤i≤n be is the nodal basis of continuous and piecewise linear functions on Γ h defined by

∀1 ≤ i, j ≤ n be , ϕ i (x j ) = δ ij . (4.11)
Here and below, κ denotes the discrete curvature computed by a standard formula (see, e.g., [START_REF] Zerrouq | Optimisation de forme robuste pour la mécanique des solides et des fluides[END_REF]Appendix A.3]). By choosing the normal deformation ξ h (x) = ϕ i (x)n h (x i ), the discrete version of Newton's equation (4.8) is:

Find θ h in the form (4.10) such that (4.12) for i = 1, . . . , n be . Here n h (x) = n be i=1 ϕ i (x)n h (x i ) can be chosen as the P 1 interpolate of the normal.

Γ h θ h (x) • n h (x i )ϕ i (x) κ(x)f (x) + ∂ n h f (x) = - Γ h f (x)ϕ i (x)n h (x i ) • n h (x),
Another issue is the discretization of the boundary integral on Γ h . This can be a major source of error and lack of consistency as noted in [START_REF] Vie | Second-order derivatives for shape optimization with a level-set method[END_REF]. We choose a trapezoidal rule, which is known to converge rapidly for periodic smooth functions. In view of (4.11), equation (4.12) becomes

α i κ(x i )f (x i ) + ∂ n h f (x i ) = -f (x i ), for i = 1, . . . , n be , (4.13)
where we used that θ h (x i ) • n h (x i ) = α i and where One can see that the system of equations (4.13) can be easily solved, since the left-hand side has a diagonal matrix which is invertible as long as it contains no zero entries. This reminds the method of lumped masses used for the finite element approximation of parabolic problems [START_REF] Thomée | Galerkin finite element methods for parabolic problems[END_REF]Chapter 15]. In our case, we also deal with the normal deformations.

∂ n h f (x i ) = ∇f (x i ) • n h (x i ). ( 4
In conclusion, finding a discrete descent direction θ h through Newton's equation boils down to solving the system of equations (4.13) for the values of α i at each vertex, and taking ∀i = 1, . . . , n be θ(x i ) = α i n h (x i ).

(4.15)

We stress that at each step (4.2) of Newton's algorithm, a filtering of the boundary mesh is necessary in order to compute correctly the discrete curvature. We used Savitsky-Golay filters (see [START_REF] Zerrouq | Optimisation de forme robuste pour la mécanique des solides et des fluides[END_REF]Appendix A.3]).

Numerical examples.

Next, we compute the optimal shapes for problem (4.4) and for three different functions f i , namely

f 1 (x, y) = (x 2 + y 2 ) 5 -2a 5 (x 5 -10x 3 y 2 + 5xy 4 ) + a 10 -b 10 , with a = 0.95, b = 0.953, f 2 (x, y) = ((x -0.5) 2 + y 2 )((x + 0.5) 2 + y 2 ) -0.51 4 , f 3 (x, y) = max(x 2 -0.6 2 , y 2 -0.6 2 ).
These functions were inspired by [START_REF] Vie | Second-order derivatives for shape optimization with a level-set method[END_REF]Section 12.1].

The initial shape is a unit disc that we successively deform using the discretization (4.13) and a constant time step t = 1 (cf. (4.2)). The results are given in Figures 8910.

One can see that the optimal shapes are easily computed after 7 iterations for f 1 and f 2 . The convergence rate is very fast although it seems geometric rather than Figure 10 helps to illustrate an important issue with the computation of optimal shapes. The discretization (4.13) assumes that the optimal shapes are smooth, not only to be able to define the different discrete quantities necessary for the algorithm such as the normals to the boundary n h , but also the structure of the second order shape derivative adopted in (4.6) is only valid for domains Ω that are at least C 3 . Thus Newton's algorithm will fail to capture optimal shapes that contain corners or cusps. This is noticed in Figure 10 where the set of negative values f 3 is clearly a square, but the algorithm never captures the square and gets stuck oscillating around shapes close to a square. This oscillation can be seen in Figure 12 where we give the L ∞ (Ω h ) norm of the computed descent directions θ h . We see that θ h ∞ converges rapidly for E(Ω, f 1 ) and E(Ω, f 2 ) after only 7 iterations. In contrast, for E(Ω, f 3 ), where the optimal shape is a square, we can see that while the value of the norms is low, it fails to converge and keeps on oscillating.

Moreover, Figure 12 shows that the L ∞ norm is constant around the optimal shape. Our experiments suggest that this can be taken as a stopping criterion for Newton's algorithm. This contrasts with first order algorithms where θ h ∞ tends to oscillate around critical shapes. 4.3. The Dirichlet energy. We now consider a functional that depends on the domain Ω through the solution u of a PDE. We consider the case of the Dirichlet energy

E(Ω) = Ω |∇u| 2 , (4.16)
where u is solution to the Dirichlet equation

-∆u = 1 in Ω, u = 0 on Γ. (4.17)
It is well known (see for example [START_REF] Bucur | Do optimal shapes exist?[END_REF][START_REF] Pironneau | Optimal shape design for elliptic systems[END_REF][START_REF] Velichkov | Existence and regularity results for some shape optimization problems[END_REF]) that the following shape optimization problem, max We seek to compute this optimal shape Ω with our algorithm and to this end, we need boundary expressions for the first and second order derivatives E and E . These derivatives have been extensively studied in shape optimization literature (see for example [START_REF] Henrot | Shape variation and optimization[END_REF][START_REF] Laurain | Distributed and boundary expressions of first and second order shape derivatives in nonsmooth domains[END_REF] and references therein). Here, we give the expressions in a way suited to our discretization (see [START_REF] Zerrouq | Optimisation de forme robuste pour la mécanique des solides et des fluides[END_REF]Section 8.2.1] or [19, Section 5.9.6]). We have: Theorem 4.1. Let Ω ∈ O 3 . Then the Dirichlet energy E(Ω) defined in (4.16) is twice shape differentiable, and for all θ, ξ ∈ C 1,∞ (R 2 ; R 2 ) we have

Ω∈O 3 |Ω|=m E(Ω), (4.18)
E (Ω; ξ) = Γ (ξ • n)(∂ n u) 2 (4.19)
and

E (Ω; θ, ξ) = 2 Γ (ξ • n)∇u θ • ∇u + Γ (θ • n)(ξ • n) (∂ n + H)|∇u| 2 + Γ Z θ,ξ |∇u| 2 , (4.20)
where Z θ,ξ is given by (4.7). The function u θ is the shape derivative of the state function u and is solution to the following problem:

-∆u θ = 0 in Ω, u θ = -(θ • n)∂ n u on Γ. (4.21)
The expression of the second order shape derivative E given by (4.20) may seem costly for numerical applications since one needs to compute u θ at every iteration. We shall see below how to deal with this difficulty. Remark 4.2. Note that the first term in the right-hand side of (4.20) is bilinear and symmetric since ∇u = (∂ n u)n (recall that u Γ = 0) and so

Γ (ξ • n)∇u θ • ∇u = Γ (ξ • n)∂ n u∂ n u θ = - Γ u ξ ∂ n u θ = - Ω ∇u ξ • ∇u θ .
Before proceeding with the discrete optimization problem, we need the first and second order shape derivatives of the volume (area) constraint V (Ω) = Ω dx. The following corollary is easily obtained by taking f (x, y) = 1 in (4.6).

Corollary 4.2. Let Ω ∈ O 3 . Then the area functional V (Ω) is twice shape differentiable, and for all θ, ξ ∈ C 1,∞ (R 2 ; R 2 ) we have

V (Ω; ξ) = Γ (ξ • n), (4.22)
and

V (Ω; θ, ξ) = Γ H(θ • n)(ξ • n) + Γ Z θ,ξ , (4.23)
where Z θ,ξ is given by (4.7).

Discretization of Newton's equation with volume constraint.

Let us now look at the shape optimization problem (4.18). We now have to deal with the volume constraint throughout the optimization process, and to do so we follow the results in [6, chapter 14]. In general, for a given constraint c(Ω), the right approach is to minimize the objective E(Ω) while keeping the constraint satisfied, which is called a primal-dual method. To describe this, according to the first-order optimality conditions, we know that when the constraint is qualified at a solution Ω , there exists a Lagrange multiplier λ * ∈ R such that

E (Ω ; ξ) + λ * c (Ω , ξ) = 0, ∀ξ, c(Ω ) = 0. (4.24)
Thus, introducing the Lagrangian

L(Ω p , λ p ) = E(Ω p ) + λ p c(Ω p ), (4.25)
Newton's method defines a step in (Ω, µ) at (Ω p , λ p ) by linearizing the system (4.24). One finds

L p A T p A p 0 θ p µ p = - L p c p , (4.26)
where L p := L (Ω p ; θ, ξ) and L p := L (Ω p ; ξ) are the first and second order shape derivatives of the Lagrangian L at Ω p (for a fixed λ p ), c p = c(Ω p ) and A p := c (Ω p ; ξ) is the shape derivative of the constraint. Given a solution (θ p , µ p ) to (4.26), Newton's method defines the next iterate (Ω p+1 , λ p+1 ) by

Ω p+1 = (I + θ p )(Ω p ) and λ p+1 = λ p + µ p . (4.27)
We then proceed to build the left hand side of (4.26) using the discretization described earlier in Section 4.2.1. Once again, we seek a continuous piecewise linear vector field θ h as in (4.10). By choosing the discrete normal vector field ξ h (x) = ϕ i (x)n h (x i ) where ϕ i is the hat function at x i , we obtain from (4.22) that A p is the row vector with entries

A p [i] = Γ h ϕ i (x)n h (x i ) • n h (x)dx ≈ h - i + h + i 2 n h (x i ) • n h (x i ) = h - i + h + i 2 .
As previously, we used the trapezoidal rule to evaluate the integral above. Moreover, for each vertex x i , h + i and h - i are the lengths of the edges [x i-1 , x i ] and [x i , x i+1 ] arriving at x i .

The matrix L p is the Hessian of the Lagrangian L. Using the expressions (4.20) and (4.23), we obtain that L p is a matrix of size n be × n be with entries given by

L p [i, j] = Γ h 2(ϕ i (x)n h (x i ) • n h (x))∇u h n j h (x) • ∇u h (x) + Γ h (ϕ i (x)n h (x i ) • n h (x))(ϕ j (x)n h (x j ) • n h (x))(∂ n h + κ(x))|∇u h (x)| 2 + λ Γ h (ϕ i (x)n h (x i ) • n h (x))(ϕ j (x)n h (x j ) • n h (x))κ(x)dx.
In this expression, u h is the discrete P 1 solution to the state equation (4.17) and u h n j h is the discrete (P 1 ) solution to problem (4.21) with θ(x) = ϕ j (x)n h (x j ). Using a trapezoidale rule to evaluate the integrals, we obtain

L p [i, j] = h - i + h + i 2 2∇u h n j h (x i ) • ∇u h (x i ) +δ ij h - i + h + i 2 (∂ n h + κ(x i ))|∇u h (x i )| 2 + λκ(x i ) .
Owing to the term δ ij , the matrix L p is diagonal, except for the first term involving ∇u h n j h (x i ) which is nonlocal. In practice, we only keep the diagonal terms in L p , leading to a quasi-Newton method.

Going back to the continuous problem (4.21), we have u θ = (θ • n)∂ n u. In particular, tor θ = n, this gives u n (x i ) = ∂ n u(x i ). Thus, we can write

∇u n = ∇(∇u • n) = D 2 u n + Dn.∇u = D 2 u.n + ∂ n u Dn.n (∇ Γ u = 0) = D 2 u.n (Dn.n = 0 since n is unitary on Γ).
This allows us to rewrite ∇u h n i h

(x i ) = D 2 u h (x i ).n h (x i ).
This could seem problematic at first sight. Indeed, u h is only P 1 and thus ∇u h = (∂ x u h , ∂ y u h ) T is only a P 0 function on T h , i.e., it is constant on each triangle T i . We thus seek to lift this function as a P 1 function. The easiest way to do so is by interpolation, but our experiments show that the most accurate way to do this is through a least-squares matching of the gradient ∇u h on the vertices x i of Γ. This is described in detail in Appendix A. In particular, we define the operator ∇ h (see (A.6)),

∇ h : V h =⇒ V h × V h ϕ h → ∇ h ϕ h ≈ ∇ϕ h .
Not only is this an accurate way of approximating the gradient, but its (almost) local character also allows us to cut down on unnecessary computations on the interior of the domain Ω h . Indeed, L p is only defined on the boundary Γ h , so we only need the values of the gradient as a P 1 function on a neighborhood of the boundary vertices x i . Notice also that this solves the issue of computing the term ∂ n |∇u| 2 , since we can compute

∂ n h |∇u h | 2 (x i ) = ∇ h i |∇ h u h | 2 • n h (x i )
, where, for a continuous function g, ∇ h i g = ∇ h g(x i ). We stress that the operator ∇ h : C 0 (Ω h ) → V h × V h is actually defined for any continuous function. Here, it applies to the continuous

P 2 function |∇ h u h | 2 .
Ultimately, the diagonal entries of the Hessian L p and the entries of L p read

L p [i, i] = 2∇ h i ∇ h u h .n h (x i ) • ∇ h i u h + ∇ h i |∇ h u h | 2 • n h (x i ) + κ(x i ) λ + |∇ h i u h | 2 (4.28) L p [i] = λ + ∇ h i u h • n h (x i ) 2 . (4.29)
We have divided these entries by (h + i + h - i )/2 which is common to both sides. 4.3.2. Numerical results. In this section, we test our discretization to compute the optimal shape Ω for minimizing the Dirichlet energy (4.16) under volume (area, since d = 2) constraint |Ω| = 0.4. Recall that the optimal shape in this case is a disc of area 0.4. Starting from a random guess, we successively compute a descent direction by solving (4.26) that we use to move the shape with a constant step t = 1 as defined in (4.27). We also compare Newton's method with an augmented Lagrangian method [START_REF] Nocedal | Numerical optimization[END_REF].

As stated in the introduction of [START_REF] Feppon | Null space gradient flows for constrained optimization with applications to shape optimization[END_REF], the tuning of the penalty parameters for the augmented Lagrangian is strongly case-dependent. Thus, we perform a series of tests with different penalty parameters and time steps and we retain the fastest in order to compare it with Newton's algorithm. We take a constant time step Figure 13 shows the starting shape and the optimal shape. In both cases, the algorithm converges to a ball, as expected.

As seen in Figures 14 and15, Newton's algorithm is remarkably fast and converges in about 20 iterations. In contrast, the augmented Lagrangian takes 10 times longer to converge. Figure 16 shows that the volume constraint is reasonably satisfied during the whole process.

5.

Optimal hull with optimized support 5.1. Theoretical results. We introduce a "bounding box", namely a bounded and connected open subset Q of the upper half plane {(x, z) ∈ R 2 : z > 0}. We assume that Q has a Lipschitz boundary ∂Q such that ∂Q ∩ (R × {0}) is a segment of the x-axis (possibly empty). We denote by ∂Q N the (relative) interior of this segment and ∂Q 0 = ∂D \ ∂Q N .

For instance, Q can be a half disc with a large radius and centered on the x-axis or a rectangle with one side of its boundary included in the x-axis.

We introduce the Sobolev space

H(Q) = u ∈ H 1 (Q) : u = 0 on ∂Q 0 in the sense of trace ,
equipped with the H 1 norm. For each function u ∈ H 1 (Q), we denote its support by

D u = {(x, z) ∈ Q : u(x, z) = 0}.
The set D u is unique up to a set of zero Lebesgue-measure and its area |D u | does not depend on the choice of the representative of u.

Let a ∈ (0, |Q|) (an area) and V > 0 (a volume). Following [START_REF] Dambrine | Regularity of optimal ship forms based on Michell's wave resistance[END_REF], we introduce the set

C V ,a (Q) = u ∈ H(Q) : Q u dxdz = V and |D u | ≤ a . The set C V ,a is a closed subset of H(Q).
As in Section 3.1, we assume that the Kelvin wave number is a random variable ν : Ω → R whose probability distribution has a compact support in (0, +∞).

For u ∈ H(Q), the functional is

J (u) = E h(ν) 2a + J 0 (u) + 1 C F J wave (u, ν) ,
where h : (0, +∞) → (0, +∞) is a continuous function which depends on the fixed parameters ρ, g and C F . The functions J 0 and J wave are defined as previously, except that the integration on D is replaced by an integration on Q. That is, we have

J 0 (u) = Q |∇u| 2 dxdz and J wave (u, α) = 8α 4 π ∞ 1 |T (u, α, λ)| 2 λ 4 √ λ 2 -1 dλ,
where T is defined for all u ∈ H(Q), α > 0 and λ > 0 by

T (u, α, λ) = Q u(x, z)e -iλαx e -λ 2 αz dxdz.
We recall that the bounding box Q and the function h are given, and that the positive parameters ρ, g, C F , a and V are fixed. The problem reads:

Find u ∈ C V ,a (Q) such that J (u ) = min u∈C V ,a (Q) J (u).
(5.1)

The support D u of u will be called on optimal domain. We note that u is generally not unique, because J is invariant with respect to translations along the x-axis. Moreover, the set C V ,a (Q) is not convex, so that an optimal domain D u is not necessarily unique up to translations along the x-axis. We have:

Theorem 5.1. Problem (5.1) has at least one solution u .

Proof. Let (u n ) n be a minimizing sequence in C V ,a (Q). We have

J (u n ) = 2aE h(ν) + E h(ν) J 0 (u n ) + 1 C F E h(ν)J wave (u n , ν) ,
and each term in the sum above is nonnegative. This implies that the sequence (J 0 (u n )) n is bounded, and so (u n ) is bounded in H 1 (Q), by the Poincaré inequality. Let q > 5. By compactness, there is a subsequence still denoted by (u n ) and u ∈ H(Q) such that (u n ) converges to u weakly in H(Q), strongly in L q (Q) and a.e. in Q. By Fatou's lemma, we have

|D u | = Q 1 u dxdz ≤ lim inf n Q 1 un dxdz = |D un | ≤ a,
where 1 v denotes the characteristic function of a function v ∈ H(Q). Thus, u belongs to C V ,a (Q). By lower semi-continuity of J 0 , we have

J 0 (u ) ≤ lim inf n J 0 (u n ). ( 5.2)
We may apply Proposition 3.4 with the set D replaced by Q. Since (u n ) converges strongly in L q (Q), we have E h(ν)J wave (u n , ν) → E h(ν)J wave (u , ν) .

(5.

3)

The relations (5.2) and (5.3) show that

J (u ) ≤ lim inf n J (u n ).
Since (u n ) is a minimizing sequence, u is a solution to problem (5.1).

Numerical simulations.

We want to compute the optimal shape for minimizing the expectation of the total water resistance under area constraint. That is, for a given area a > 0, we wish to solve numerically the following problem, min |D|=a J (D), (5.4) among bounded open subsets D of the upper half plane which are admissible. As a shortcut, we have denoted by J the functional

D → J D (u D ),
where J D is defined by (3.8) with the choice (3.17 where

J (D) = D |∇u D | 2 + 1 C F E[ h(ν)] D×D H ν (x, z, x , z )u D (x, z)u D (x , z ),
and u D is the solution to the boundary value problem (cf. (3.16))

                   -∆u(x, z) + 1 C F E h(ν) D×D H ν (x, z, x , z )u(x , z )dx dz = C, ∀(x, z) ∈ D, D u dxdz = V , u = 0 on Γ 0 , ∂u ∂n = 0 on Γ N .
(5.6) Throughout Section 5.2, D is a simply connected and bounded open subset of the upper half plane with Lipschitz boundary. We assume that ∂D ∩ (R × {0}) is a segment and we denote by Γ N the relative interior of this segment. We also assume that the curve Γ 0 = ∂D \ Γ N is of class C 3 (see Figure 2). 5.2.1. Shape derivatives. The shape sensitivity analysis of J is very similar to the case of the Dirichlet energy, assuming enough regularity. The result is very similar to Theorem 4.1 (up to the minus sign). The calculations for a smooth domain and a smooth kernel are given in [START_REF] Zerrouq | Optimisation de forme robuste pour la mécanique des solides et des fluides[END_REF]Section 8.3] (see also [START_REF] Henrot | Shape variation and optimization[END_REF]).

In our case, we will use the following expressions.

Let R 2 + = R × [0, +∞) denote the closed upper half plane. For all θ, ξ ∈ C 1,∞ (R 2 + ; R 2 ) such that θ • n = 0 and ξ • n = 0 on the x-axis, we have J (D; ξ) = - Γ 0 (ξ • n)(∂ n u) 2 (5.7) and J (D; θ, ξ) = -2 Γ 0 (ξ • n)∇u θ • ∇u - Γ 0 (θ • n)(ξ • n) (∂ n + H)|∇u| 2 - Γ 0 Z θ,ξ |∇u| 2 (5.8)
where u solves (5.6), Z θ,ξ is given by (4.7) and u θ is solution to the following problem:

                   -∆u θ + 1 C F E h(ν) D×D H ν (•, •, x , z )u θ (x , z )dx dz = C in D, D u θ dxdz = 0, u θ = -(θ • n)∂ n u on Γ 0 , ∂u θ ∂n = 0 on Γ N .
(5.9)

The same discretization as in Section 4.3.1 can be used to apply Newton's algorithm to the set functional J , assuming enough regularity both on the domain D and the kernel H ν . Even in situations where these regularity constraints are not satisfied, for example when the domain contains corners or when the kernel is not regular enough, the discretization (4.26) provides enough flexibility to get around these issues, as explained below in Section 5.2.2. 5.2.2. Tips for Newton's algorithm. In this section, we seek to compute the optimal shape D by Newton's method discretized as in Section 4.3. There are two main obstructions to regularity in the expression of the second order shape derivative (5.8). First, Michell's kernel H ν does not belong to L 2 (D). It is only in L 5 4 -ε (D) for all ε > 0 small enough, as proved in [START_REF] Dambrine | Regularity of optimal ship forms based on Michell's wave resistance[END_REF]. Second, the boundary expression of the second order derivative given by (5.8) is only given for domains D that are at least C 3 (with Γ N = ∅ and Γ 0 = ∂D) while D here is only C 3 on each boundary Γ 0 and Γ N . We shall always have two corners at the intersection of the boundaries.

These regularity issues make it a priori difficult to attempt a second order shape optmization method. Since the boundary expression of the second order shape derivative is ill-defined, one could then try to find a different expression that takes into account the irregularity of the boundary. This is extensively studied in [START_REF] Laurain | Distributed and boundary expressions of first and second order shape derivatives in nonsmooth domains[END_REF] for the case of the Dirichlet energy and volume functional, and the interested reader can find in it expressions for the second order derivative given for domains that are only Lipschitz.

In our case, it turns out that the discretization (4.13) gives enough flexibility to get around these issues when building the discrete Hessian matrix. Indeed, the discrete Hessian shall only be constructed on the vertices of the boundary Γ 0 . Since our discretization is strictly dependent on the discrete normal n h (x i ) on a vertex x i , in order to remain consistent, we exclude the vertices x i with singularities (in this case, the first and last vertex of Γ 0 ) when building the discrete Hessian. The value θ h (x i ) on each one of these two vertices is then obtained by keeping the horizontal component of the neighboring value θ h (x i+1 ) computed on Γ 0 . In this way, the constraint on Γ N is preserved, namely

θ h (x i ) • n h (x i ) = 0 on Γ N .
Finally, while theoretically the solution u is not in H 2 (D), we still can discretize quantities such as ∇ 2 u using numerical gradients, specifically the least-squares gradient defined in Appendix A. This least-square gradient will enjoy higher regularity and ensures that the resulting discrete Hessian matrix is well conditioned. These are all the tips that were used to help Newton's algorithm to converge on an irregular domain with an irregular solution, and the results are given in the next For the sake of comparison with the results in [START_REF] Dambrine | Regularity of optimal ship forms based on Michell's wave resistance[END_REF], we introduce the area Froude number

F 2 r = 1 α √ a = V 2 g √ a ,
where a = 0.518 m 2 is the area of D. The bounds ν min and ν max are computed as For our numerical results, we focus on the Froude interval [0.6, 1.0]. This is the range of speed where it is most interesting to have a bulb in order to minimize the total ship resistance, as seen in [START_REF] Dambrine | Continuity with respect to the speed for optimal ship forms based on Michell's formula[END_REF].

We first compute the optimal domain D which minimizes the expectation of the total water resistance J . Figure 18 (bottom) shows the optimal domain computed with Newton's algorithm. The corresponding optimal hull u is represented in Figure 19 (bottom). The convergence of the algorithm is given in Figure 17. We can see that Newton's algorithm converges very fast without the need of a line search, since a constant time step t = 1 was used throughout the optimisation process.

As a comparison, we have also computed the domain which minimizes the total resistance for the mean Kelvin number α = E(ν) related to the Froude number F r = 0.73 (Figure 18, top) and the corresponding optimal hull ū (Figure 19, top). We use the same algorithm (Newton's method), except that the averaged kernel H ν is replaced by Michell's kernel h(α)k α (see Remark 3.2), so the computation is less costly.

The two optimal domains and the corresponding hulls u and ū are different: the bulbous bow is more pronounced for ū than for u . In contrast u * is longer than ū since it takes into account higher Froude numbers. We refer here to [START_REF] Dambrine | Continuity with respect to the speed for optimal ship forms based on Michell's formula[END_REF]Section 7] where it is seen that the length of the optimal domain increases with the Froude number.

Despite these remarkable geometric differences, we can see in Figure 20 that the total resistance is very similar for both hulls on the whole interval [0.6, 1.0]. The expectation of the total resistance is equal to 26.3 N for u and it is equal to 27.1 N for ū. We note however that ū is optimal for α = E(ν) and slightly better for Froude where A and c are the n × 2 matrix and n-dimensional vector given by:

A =      x 1 1 -x 1 0 x 2 1 -x 2 0
x 1 2 -x 1 0 x 2 2 -x 2 0 . . . . . . where A T is the transpose of A. In (A.4), A T c is a two-dimensional vector and A T A is the 2 × 2 matrix defined by

A T A =     n i=1 (x 1 i -x 1 0 ) 2 n i=1 (x 1 i -x 1 0 )(x 2 i -x 2 0 ) n i=1 (x 1 i -x 1 0 )(x 2 i -x 2 0 ) n i=1 (x 2 i -x 2 0 ) 2     .
The matrix A T A is invertible since det(A T A) > 0 as a consequence of the Cauchy-Schwarz inequality. By repeating the procedure above, we obtain a value for the gradient of u at every vertex of the triangulation T h , and by linear interpolation, a function (v h , w h ) in V h × V h . This approach for approximating the gradient ∇u as a continuous P 1 function is pretty robust in practice. Unfortunately, it may still be affected by numerical errors, and this is a reason why it is often recommended to smooth the resulting quantities.

Figure 1 .

 1 Figure 1. The immerged hull: half of it is represented by a positive function y = f (x, z). The other half hull is obtained by symmetry.

Find a set

  D U which minimizes R total (D, f D , U ) among all bounded regular open subsets D of the lower half-plane such that |D| = a.

(2. 14 )

 14 Here and below, |D| is the area of the domain D. Let h be a positive continuous function of U . Since U is fixed, h(U ) is a constant and we may use the cost function h(U )R total (D, f D , U ) in the problem above without changing its solution(s) D U . Now we assume again that U is a random function and we are interested in the more robust problem: Find a set D which minimizes E [h(U )R total (D, f D , U )] among all bounded regular open subsets D of the lower half-plane such that |D| = a.

(3. 22 )

 22 By Theorem 3.5, a minimizer u D of JD (u) in H V (D) is a solution to the linear boundary value problem(3.16).

3. 3 . 3 .

 33 Numerical results. The parameters L, T , V , ρ, g and C F are the same as in Section 3.3.1. The probability distribution of the Kelvin wave number ν = g/V 2

Figure 3 .

 3 Figure 3. Mass distribution (left), and 3D results (right) of the Wigley hull w (top), the optimal hull u D (middle, optimized for a uniform distribution of ν) and the optimal hull ūD (bottom, optimized for α = E[ν] with ν uniform).

Figure 4 .

 4 Figure 4. Total water resistance (in N) for the three hulls u D , ūD and w (case of a uniform distribution of ν).

Figure 5 .

 5 Figure 5. Mass distribution (left), and 3D results (right) of the Wigley hull w (top), the optimal hull u D (middle, optimized for a uniform distribution of V ) and the optimal hull ūD (bottom, optimized for α = E[ν] with V uniform).

  the set of bounded open subsets of R 2 that are of class C k . The first order shape derivative is defined classically as follows [37, Definition 3.2.11]. Definition 4.1. Let k ≥ 1 and E(Ω) be a function from O k into R. We define

Figure 6 .

 6 Figure 6. Total water resistance (in N) for the hulls u D , ūD and w (case of a uniform distribution of the speed V ).

Figure 7 .

 7 Figure 7. Discrete setting around a vertex x i

Figure 8 .Figure 9 .Figure 10 .

 8910 Figure 8. Optimal (right) and intermediate shapes (left and middle) for minimizing E(Ω, f 1 )

Figure 11 .

 11 Figure 11. The error to the best shape Ω at every iteration of the minimization process for the different functionals E(Ω, f i )

Figure 12 .

 12 Figure 12. The L ∞ (Ω h ) norm of the computed descent direction θ h at every iteration of the minimization process for the different functionals E(Ω, f i )

Figure 13 .

 13 Figure[START_REF] Frey | Mesh generation[END_REF]. The starting shape (left), a nearly optimal shape (middle) and the best shape (right), with the magnitude of the descent direction θ for the minimization of the Dirichlet energy.

Figure 14 .Figure 15 .

 1415 Figure 14. Convergence of the Dirichlet energy

Figure 16 .

 16 Figure 16. Convergence of the area constraint |V (Ω p ) -V target |.

  ) for h and u D is the solution to the minimization problem(3.18).By(3.22), solving (5.4) is the same as solving min |D|=a

Figure 17 .

 17 Figure 17. Convergence results for the expectation of the total resistance and for the volume constraint using Newton's algorithm.

Figure 18 .

 18 Figure 18. The optimal domains for minimizing the total resistance at F r = 0.73 (top) and for minimizing the expectation of the total resistance for F r ∈ [0.6, 1] with ν uniform (bottom).

Figure 19 .

 19 Figure 19. The optimal hulls for minimizing the total water resistance at F r = 0.73 (ū, top) and for minimizing the expectation of the total resistance for F r ∈ [0.6, 1] with ν uniform (u , bottom).

Figure 20 .

 20 Figure 20. The total water resistance of the two hulls

Figure 21 .

 21 Figure 21. Neighborhood of a vertex x 0 in a mesh T h .

1 )

 1 -u(x 0 ) u(x 2 ) -u(x 0 ) . . . u(x n ) -u(x 0 ) (A.2) is clearly overdetermined, and in practice we search b as a solution to the least-square problem A T Ab = A T c, (A.4)

Appendix A. Gradient of a piecewise linear function

In this appendix, we define a continuous gradient of a P 1 function. This is used in our shape optimization algorithms.

For a domain D in R 2 with a Lipschitz boundary Γ, we denote by T h a conforming triangulation of D with typical mesh size h. The triangulation T h contains n t triangles T i . The approximation of the domain D is the polygonal domain D h = ∪ nt i=1 T i and we denote by Γ h the boundary of D h .

The space of continous piecewise linear (P 1 ) functions is denoted by V h , namely

where P 1 is the space of polynomials of degree less than or equal to 1.

If u belongs to V h , then the gradient ∇u is only a discontinuous P 0 function on T h . i.e. it is equal to a constant value α i on each triangle T i . In many cases, we seek to lift this gradient to a P 1 function defined on the vertices of T h through a stable approximation.

Let x 0 be a vertex of T h and denote by x 1 , . . . , x n the neighboring vertices as given in Figure 21. For a smooth function f : R 2 → R, we have the following Taylor To this end, it is common practice to use the screened Poisson equation, i.e. for a P 1 function ψ ∈ V h , we trade ψ for the solution ψ ∈ V h to the following variational problem:

where ε is chosen of the order of mesh size, that is, ε ≈ h. Summing up, for a function u ∈ V h , we define its continuous P 1 gradient as

where ṽh ∈ V h and wh ∈ V h are obtained by smoothing v h and w h respectively, and (v h , w h ) ∈ V h × V h is obtained at each vertex of the triangulation by the local least square problem described above. We note that this definition of ∇ h u makes senses for any continuous function u ∈ C 0 (D h ).