Keywords: shape optimisation, wave-making resistance, Neumann-Kelvin problem, boundary integral equation, level-set method

First, all the praises and thanks be to Allah with His compassion and mercifulness to allow me accomplishing this PhD project.

I would like to Acknowledge the Ministry of Research, Technology and Higher Education for the BPPLN study grant from Octobre 2014 to Septembre 2018, which gave me the opportunity to obtain the PhD in University of Poitiers, France.

I would like to express my sincere gratitude to my advisor, Dr Julien Dambrine, for his enthusiastic encouragement, valuable and constructive suggestions during the planning and development of this research work. I would like to express my great appreciation to my thesis director, Dr Morgan Pierre, for his patient guidance, advice and help during my study. I will always remember the help from both of you since the first day my family and I arrived at Poitiers. I also extend my grateful thanks to the PhyDROMAT team chaired by Germain Rousseaux, as a collaboration between math and physics department in University of Poitiers.

I would like to acknowledge Dr Marion Darbas and Prof. Edouard Oudet who have been willing to become my referees, and also Prof. Marc Dambrine, Prof. Noureddine Igbida, Dr Mario Ricchiuto and Dr Germain Rousseaux as my jury.

Special gratitude goes out to Prof. Dr Thamrin Usman, DEA, the rector of the University of Tanjungpura (Untan) for all the support for the success of my PhD study. My special thanks to Prof. Gérard Mauco and his wife, Christine Mauco, for their friendliness and support given to my family and me during our stay in Poitiers.

I would also like to extend my thanks to Alexandra Sarti the head of laboratory, Samuel Boissiere the head of doctoral school SISMI, Alicia Lecesve the secretary of the doctoral school, and Silvy Perez the school service of SFA, for their hospitality and services. Special i thanks also to technicians of the laboratory of LMA University of Poitiers, Benoit Metrot for his help in offering me the resources in running the program and troubleshooting the software and hardware problems, Jocelyne Attab for her hospitality and help in printing manuscript and the other documents, Myriam Andre and Nathalie Marlet for their warmth and support in searching the bibliography and Bridgitte Brault and Nathalie Mongin for their help along my study. I would also like to thanks my labmates in LMA University of Poitiers.

I would also like to acknowledge the support and the family-like atmosphere offered by my Indonesian colleagues in Poitiers. I would also like to give my thankful to my family in Netherland.

I am grateful to my parents H. Ateng Suhendar dan Hj. Enan Asmanah, who have provided me with moral and emotional support in my life, who always say I can do it! Many thanks also to my mother-in-law Hj. Rumiati for her du'a and support and also to my fatherin-law Alm. Ahmad Dhani Ramadan for his spirit at the first step in this PhD journey; we succeed Papa!!!. I am also grateful to my lovely big brother Asep Wahyudin and his wife Hali for their du'a and encouragement when I fall into a "deep gap" and also Mbak Yeni, keep your spirit on!, Abang Fatahillah Abror, Wahyu Hidayatullah and Lintangsari and also my other family members and friends who have supported me along the way.

Finally, I give my significant praise to my husband, Dr Yoga Satria Putra for his support, encouragement, and love. My daughters, Annisa Syifa Ramadhani and Afiyah Nada Zahirah, my lovely son Ahmad Yusuf Putra, I can not reach this achievement without you. I believe you can do better than me. After all that we have faced together... we did it!!!

General Introduction

Many problems in modern industries involve the question: "what is the best design for the shape of a certain object, regarding a given objective function, and taking into account several physical constraints?". The answer to the question is indispensable in the design and construction of industrial structures. The goal of shape optimisation is to provide a solid and applicable mathematical framework to answer this type of questions. Perhaps the most notable example of shape optimisation in the context of fluid mechanics arises in aeronautical engineering, in which the shape of the wings of an aircraft needs to be optimised from the point of view of the drag, with a constraint on the lift (see, for instance [61], chap. 9). Other examples of shape optimisation for aeronautical engineering can be found in [START_REF] Epstein | Accurate CFD driven optimization of lifting surfaces for wingbody configuration[END_REF], [START_REF] Schilling | Efficient shape optimization for certain and uncertain aerodynamic design[END_REF], [START_REF] Peter | Review: Numerical sensitivity analysis for aerodynamic optimization: A survey of approaches[END_REF], [START_REF] Jameson | Aerodynamic shape optimization techniques based on control theory[END_REF], [45] and many others. In acoustic engineering, some problem examples are: designing a sound isolation room, building an acoustic horn (see [START_REF] Bängtsson | Shape optimization of an acoustic horn[END_REF]) and building sound barrier in transportation noise reduction (see [START_REF] Liu | Shape optimization of sound barrier using an isogeometric fast multipole boundary element method in two dimensions[END_REF]). Finally, structural engineering was probably the largest provider of both theoretical and practical problems in shape optimisation, and has driven the development most of the tools that we have used in this thesis. For instance, the textbook of G. Allaire [START_REF] Allaire | Conception optimale de structures[END_REF], often uses a flexibility minimisation problem for a cantilever as an application of the shape optimisation tools presented.

The wave making resistance & ship design

The motivation for the study presented in this thesis is the reduction of the fuel consumption of ships. This problem is of great importance from an economical standpoint, since, as Figure 1 -Water waves behind a ship [START_REF] Contributors | Kelvin wake pattern generated by a small boat[END_REF] and a submarine [15] specified in a 2010 report of the IMO1 [START_REF]Imo's contribution to sustainable maritime development: Capacitybuilding for safe, secure and efficient shipping on clean oceans through the integrated technical co-operation programme[END_REF] over 90% of the world's trade is carried by sea.

Beside the purely economical aspect, efficiency is also associated to lower greenhouse effect gaz emissions.

As stated by The American Bureau of Shipping [START_REF]The american bureau of shipping (abs), ship energy efficiency measures: Status and guidance[END_REF], the fuel consumption can be reduced by optimising size and capacity of the ship, by reducing the speed, and by minimising the resistance of the ship. While the size, capacity and cruise speed are mostly constrained by the desired amount of goods transported per unit of time, the resistance to motion of a ship can be reduced by small modifications on the shape of its hull (see for instance [START_REF]The american bureau of shipping (abs), ship energy efficiency measures: Status and guidance[END_REF][START_REF] Harries | Advanced hydrodynamic design of container carriers for improved transfort efficiency[END_REF]).

The resistance of a ship is defined by the sum of all the forces exerted on the hull in the direction of the motion. The wetted surface of the vessel raises a force called the viscous drag or viscous resistance which has an opposite direction to the relative flow velocity and is proportional to the wetted surface. In addition, when the ship moves, it generates a system of waves on the water/air interface which we also call the called the free surface. In figure 1, we show pictures of the famous v-shaped pattern of the wake behind a moving ship and a submarine. The energy required to generate these waves is associated to a force developing around the body. We call this force the wave-making resistance. This wave-making resistance along with the viscous resistance and appendage resistance (such us eddy-making resistance) compose the total resistance of the ship's hull [START_REF]The american bureau of shipping (abs), ship energy efficiency measures: Status and guidance[END_REF]35].

Even though in fact these resistances cannot be separated from one another, for simplification and computational purposes, the calculation is usually considered independently.

However, the contribution of the wave-making resistance to the total resistance of ship is important and varies for different types of ships. Using TEU 2 to measure a ship's cargo carrying capacity, for a container ship with a capacity of respectively 4.500 TEU and 8.000 TEU, called Panamax and Post-Panamax, the wave-making resistance are 23% and 18% respectively of the total resistance [START_REF]The american bureau of shipping (abs), ship energy efficiency measures: Status and guidance[END_REF]. The wave-making resistance increases with ship speed [START_REF]The american bureau of shipping (abs), ship energy efficiency measures: Status and guidance[END_REF] and even could give 50% of the total resistance of a high-speed ship and/or a ship with full hull form [7]. Because of the potential gain described above and the important shape sensitivity of the wave-making resistance (see [START_REF] Weinblum | Applications of wave resistance theory to the problem of ship design[END_REF], [7]), we will focus on this particular aspect of ship hydrodynamics.

Besides the benefits of efficiency and low greenhouse effect gaz emissions, the reduction of the wave-making resistance could also help prevent the damage caused by wash waves on river banks in the context of fluvial transportation. The figure 2a shows the traveling waves produced by a moving container ship in a British canal. In figure 2b, traveling waves are also generated by a boat in Kapuas 3 river. With time, these wash wave can have a destructive behaviour on the infrastructure built to facilitate navigation in the case of artificial canals, or destroy habitats and hence endanger aquatic or semi-aquatic species in the case of natural rivers.

Although the long term goal of the present research is to determine and optimise the wave-making resistance of ships, this thesis will focus on a simpler two-dimensional submerged body case. Many technical difficulties are associated to the case of ships. First, the three dimensional nature of the problem make the calculations more difficult, and computationally expansive. Second, ships are surface piercing bodies, and this adds some difficulty in the mathematical setting of the problem because it involves a triple line (solid/liquid/air interface), on which the nature of the boundary conditions is not very clear. A remarkable work on this topic can be found in [START_REF] Delhommeau | Les problèms de diffraction-radiation et de résistance de vagues: étude théorique et résolution numérique par la méthode des singularités[END_REF] in which the author, not only obtains the wavemaking resistance of a ship in the full 3D case, but is also able to recover the lift and torques. Unfortunately, for general shape optimisation, the matter is even more difficult. Although our simplified approach cannot be applied directly for ship hydrodynamics, we will build tools that we believe can be extended in the 3D surface-piercing case without too many efforts. 2 Twenty-Foot Equivalent Unit, 1 TEU is approximately 39 m 3 . 3 West Borneo, Indonesia. Kapuas river is the longest river in Indonesia with 1143 km long. This river is a fundamental waterway which can support cargo and passenger shipping in West Borneo while connecting the center of the island and its western coast (see figure 3 for examples of water-crafts in Kapuas river) (a) A container in a British canal [52] (b) A boat in Kapuas river, Indonesia 2 Mathematical modelling of wake and resistance

A brief history of the ship waves problem

The problem of ship waves has been studied for several hundreds years. One of the first study on this topic was conducted by Leonardo da Vinci (1458-1519) who examined the movement of objects in water, in a qualitative manner. His descriptions and sketches underlined the fact that, free surface flows (as we call them today) seem to follow natural laws, that, in his time, were far from reach. Later, in 1687, Isaac Newton developed the first theory of water resistance of moving ships, and determined that the drag should be a quadratic function of its velocity. His formula was obtained by deducing the drag from the impact of a large amount of particles of water on the hull, each of them carrying some momentum. His formula came 62 years before the first complete theory of perfect fluids was formulated by Euler. An important contribution was provided in 1871 by William Froude which gave a separation between viscous resistance (well described with Newton's formula) and wavemaking resistance (that arises from Euler's equations). He later discovered a criterion for the description of ship waves which depends on the gravitational force, the depth of a body and its velocity, this criterion was later coined as the Froude number. Hereinafter, Lord Kelvin, in 1886, studied the hydrodynamic forces acting on an obstacle generating a stationary wave pattern on the fluid surface. For a complete review of historical references on this topic, we refer to [START_REF] Gotman | Navigating the wake of past efforts[END_REF].

In 1898 John Henry Michell put forward a theory of water resistance for thin ships in an ideal fluid with infinite depth [START_REF] Michell | The wave resistance of a ship[END_REF]. He deduced, from the theory of linear ship waves developed by Kelvin, the first mathematical formula for the wave resistance of a ship that takes into account an actual (supposed thin) hull shape. His formula writes:

R w (f ) = 4ρgν 3 π ∞ 1 |I(λ, f )| 2 λ 4 √ λ 2 -1 dλ , (1)
where:

I(λ, f ) = ω f (x 1 , x
3 ) e -νλ 2 z e iνλx dx 1 dx 3 .

(2)

Here, f (x 1 , x 3 ) denotes the half-width of the hull at a point (x 1 , x 3 ) of the center-plane (the plane that cuts the ship in half along its length). For this contemporaries, Michell's formula was considered impractical, because it consists of triple integrals that are difficult to obtain for general values of f . With the advent of computers and numerical simulation, the calculation of the wave-making resistance through Michell's formula for general ship forms became accessible (see [44] for instance). Later, the development of sophisticated numerical methods made possible the prediction of the wave resistance with more and more accurate models, and today, Computational Fluid Dynamics (CFD) have developed to a point that it is possible to test the full behaviour of a ship in what we call "virtual towing tanks" (in analogy with virtual wind tunnels for aeronautics). Despite this, thin ship theory regained some popularity recently (see for instance [START_REF] Noblesse | The Neumann-Michell theory of ship waves[END_REF], [START_REF] Dambrine | A theoretical and numerical determination of optimal ship forms based on michell's wave resistance[END_REF] and [8]) as it provides with an insight of the physics of the wave-making resistance that the Navier-Stokes equations (on which CFD is based) fails to give.

Mathematical formulation in the 2D non-thin case

The thin ship theory is rigorous and powerful from a mathematical point of view, but it is much more difficult to deal with ships which are not slender. There is no general formula for the wave resistance in this case. Nevertheless, Havelock [35,[START_REF] Havelock | The forces on a circular cylinder submerged in a uniform stream[END_REF] obtained an analytical formula for an obstacle with a circular cylinder form. We use it in chapter III for a numerical validation of the wave-making resistance. Havelock also obtained various results on the wave-making resistance which can be found in his collection of papers [START_REF] Havelock | The collected papers of Sir Thomas Havelock[END_REF]. We consider a fully submerged obstacle which is at rest in a fluid flow (see figure 4 for the illustration of geometrical setting of the problem). The velocity of the fluid at infinity far away from the obstacle is U ∞ e 1 where e 1 is the unit horizontal vector. Under the assumption that the fluid is homogeneous, inviscid and incompressible, the Euler equations read:

∇ • U = 0 in Ω + , ( 3
)
∂ t U + (U • ∇)U = 1 ρ ∇P -ge 2 in Ω + , ( 4
)
where U is the velocity field of the fluid, Ω + is the fluid domain, P is the pressure, ρ is the density of the fluid and g is the acceleration of the gravity field. The Euler equations were introduced by Euler in 1755, and, as we will see in chapter I, it can be derived from the (later discovered) Navier-Stokes equations by neglecting the effect of viscosity.

Assuming that the initial velocity field is curl-free, it can be shown by taking the curl of (4) (and assuming that the solution is regular enough) that the velocity field stays curlfree for all times. By the Helmhotz-Hodge theorem, we can seek therefore U as the gradient of a potential flow, namely

U = ∇Ψ in Ω + . ( 5
)
S 0 S Ω - Ω + e 1 e 2
U ∞ e 1 n

Figure 4 -Setting of the problem.

Then, ( 3) and ( 4) can be written:

∆Ψ = 0 in Ω + , ( 6
)
∂ t Ψ + 1 2 |∇Ψ| 2 + gx 2 = 1 ρ (P -P atm ) in Ω + , (7)
where P atm is the pressure at the free surface. The equation ( 7) is known as Bernoulli's equation. Further, in the case of the wave-making resistance, we have to deal with the boundary constraint i.e the condition on the free surface. This problem has its own difficulty from a mathematical point of view, because it deals with a partial differential equation which has to be solved for both an unknown function U and an unknown domain Ω + . To overcome this difficulty, by supposing that the waves have small amplitude, it is possible by linearisation to transfer the conditions on the free surface to the horizontal plane (denoted S 0 ). This technique has been applied for instance by Brard in [START_REF] Brard | Le problème de Neumann-Kelvin[END_REF] which proposed the linearization of the free surface condition. By taking the steady state for the problem, we obtain the so called Neumann-Kelvin problem (NK-problem) which reads:

                 ∂ 2 11 Φ + ν ∂ 2 Φ = 0 , on S 0 , ∆Φ = 0 , in Ω + , ∂ n Φ = -U ∞ n • e 1 , on Γ , |∇Φ| → 0 , for |x| → ∞ . (8)
Here, n is the normal to the boundary Γ pointing outwards of the obstacle, and Φ = Ψ-U ∞ x 1 is the perturbed potential. As stated in [9], the NK problem begins with the work of Rayleigh [START_REF] Rayleigh | The form of standing waves on the surface of running water[END_REF] and Kelvin [START_REF] Kelvin | On stationary waves in flowing water[END_REF] and also further by Russel [START_REF] Russel | Report on waves[END_REF]. They studied analytically the waves on the free surface generated by the movement of an obstacle in the fluid. Further, the theory around the Neumann-Kelvin problem is discussed and is used to formulate the resistance of a ship (see for instance [START_REF] Delhommeau | Les problèms de diffraction-radiation et de résistance de vagues: étude théorique et résolution numérique par la méthode des singularités[END_REF], [START_REF] Brard | Le problème de Neumann-Kelvin[END_REF], [START_REF] Dern | Un problème de Neumann-Kelvin bien posé[END_REF], [START_REF] Bauer | Le problème de Neumann-Kelvin I[END_REF], [6], [START_REF] Euvrard | Sur la régularité de la fonction de Green du problème de Neumann-Kelvin[END_REF]). Inspired by these works, in this study, we will use the NK-problem to obtain the formula of the wave-making resistance which we want to minimise. The method is introduced in next section.

The boundary integral method

Several methods have been used to deal with the steady state of the irrotational, inviscid and incompressible fluid which is described by the NK-problem. Finite element methods (FEM) have been examined by several authors to solve this problem, for instance [START_REF] Doctors | An application of the finite element technique to boundary value problems of potential flow[END_REF] and [START_REF] Bai | A variational method in potential flow with a free surface[END_REF]. FEM is a numerical method to approximate the solutions of partial differential equations which is characterized by a variational formulation. Since it needs an explicit representation of the domain, hence it is less adapted to moving domains. Another method in solving the problem is to use a boundary integral equation. In this method, the solution of the wave problem is represented by an integral equation on the boundary of the obstacle which also handles the condition at infinity. Numerical methods to approximate this boundary integral equation are called boundary element methods (BEM). They are used for instance by Noblesse [START_REF] Noblesse | Velocity representation of free-surface flow and Fourier-Kochin representation of waves[END_REF] in the free surface potential flow. Brebbia [START_REF] Brebbia | Boundary element methods for potential problems[END_REF] compared the boundary element method with the finite element method in the case of a flow around a cylinder between parallel plates; it shows that the BEM has a better computational eficiency (less points and less computer time).

However, the BEM has its own difficulties. It is necessary to find the explicit fundamental solution which handles the free surface condition, and this isn't always easy. In general, the fundamental solution (also called Green's function) depends on the problem and on the contraints that have to be satisfied, hence the method gives a restriction on the generality of the problem. However in the NK-problem, the result in finding the associated Green's function has been obtained for instance by Havelock [START_REF] Havelock | The effect of speed of advance upon the damping of heave and pitch[END_REF], Ursell [72] and Wehausen and Laitone [START_REF] Wehausen | Surface waves. Encyclopedia of Physics[END_REF] in different forms. In this present work we give a calculation of the fundamental solution that uses the exponential integral, a numerically well known special function. This fundamental solution is equivalent with that in [START_REF] Kuznetsov | Linear Water Waves[END_REF]. Once we find the Green's function, we have the advantage of the method, i.e we only have to calculate the solution of the problem by solving an integral equation where the unknown function is defined on the boundary. Hence, only the boundary of the domain needs to be discretized [16]. In addition, the BEM can solve the exterior problem where the domain is unbounded (this is the case in this study), as easily as the interior problem [16].

The shape optimisation method

We want to compute the shape of an obstacle which will minimise the wave-making resistance for the NK-problem. In order to avoid a resistance going to zero, we assume that the area of the obstacle is fixed.

Our idea is to use a standard descent gradient method for the minimisation problem. However, we deal with a set of shapes which does not define a vector space. Following a standard approach in shape optimisation, we use small deformations of the shape to define a notion of gradient of the wave-resistance functional with respect to the shape. This computation of the shape gradient is based on a Lagrangian method; we also use a Lagrangian derivative to deal with the boundary integral equation. This idea of shape derivative was first introduced by Hadamard [31] in 1907 in application of elastic problems. Later the method was developed by Murat and Simon in [START_REF] Murat | Sur le contrôle par un domaine géométrique[END_REF][START_REF] Simon | Differentiation with respect to the domain in boundary value problems[END_REF], and Allaire [START_REF] Allaire | Conception optimale de structures[END_REF]. More recent, Costabel and Le Louër in [17] established the shape derivative using Gâteaux differentiability analysis with application in electromagnetic scattering problems.

For numerical purposes, there are two major difficulties that we should overcome in this work. Firstly, we need to calculate the solution of the boundary problem so that we can calculate the wave-making resistance for a given shape. Secondly, we have to deal with the minimisation problem with respect to the shape which includes a moving domain. For both steps, we need to have an appropiate representation of the shape of the obstacle (in 2D case it is a curve). To manage the first part, the calculation of the integration on the boundary when we use the explicit parametrisation seems to be relatively more efficient because it uses less data. On the other hand, when we meet a moving domain as it is the case in the second part, we will meet a difficulty in the remeshing procedure. In contrast, if we use an implicit representation of the curve, dealing with a moving domain is rather easy. Moreover in this numerical issue, nowadays, the level set method has its superority, because it can be used to solve problems where simulation of boundary evolutions is involved [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF]. The first contribution on this method started by the seminal work of Osher and Setian in [START_REF] Osher | Fronts propagating with curvature dependent speed: Algorithms based on hamilton-jacobi formulations[END_REF]. Further, the level set methodology can perform shape changes on a domain without remeshing [76] because the method adds dynamics to implicit surfaces.

As we mentioned previously, in this thesis we couple the level set method with the boundary integral formulation to obtain the solution of the minimisation wave-making resistance problem with respect to the shape. Previous works on the level-set method to solve boundary integral equations for moving boundary problems have been made for instance by Garzon et al [29] in the context of wave breaking over a sloping beach, and Chen et al [13] for the Mullins-Sekerka dynamics. Particularly, in this thesis, we deal with evolving domains to obtain the shape minimising the wave-making resistance for a given depth and velocity.

Outline of the thesis

In this thesis, we have contributed to a better knowledge of a geometric shape optimisation for the wave-making resistance problem. The study will be presented as follows.

In the first chapter, we will present the formulation that is used in the calculation of the wave-making resistance. Here, we start with giving the derivation the Neumann-Kelvin problem from the Navier-Stokes equations. We also give the notes on the existence and uniqueness of the solution to the problem from [START_REF] Kuznetsov | Linear Water Waves[END_REF]. Then we describe the fundamental solution of the problem which will be used in representing the problem in the form of an integral boundary equation. From the knowledge of the solution of the boundary integral equation, we will give the derivation of the wave-making resistance formula. Further, to make it useful in the next chapter, we also give the notes and notation of several boundary integral operators at the end of the chapter.

The second chapter is devoted to solving shape optimisation problems which involve boundary integral equations. We begin with the introduction of the shape derivative and shape gradient and with some useful lemmas that we will need on the boundary variation method. Further, we explain the Lagrangian method and its application to the shape optimisation of the wave making resistance. We give the formal calculation of the derivative of the Lagrangian that will lead us to the "gradient" of the shape. Thus, we can use the obtained shape gradient in a gradient descent method to obtain the optimal shape. Chapter 3 is dedicated to the numerical features of the level-set method. We present a general introduction on level-set methods in the first section, which begins with the representation of curves and the signed distance function. We also give the explanation about normal motion of the signed distance function. Further, we use level-sets and tubular neighbourhoods of the boundary to obtain the numerical approximation of boundary integrals. Some simple tests are also given to enhance the explanation. Finally, in the last section we give the procedure to apply the aforementioned method to a model boundary integral equation, namely the Neumann exterior problem.

In the last chapter, we present the numerical implementation on the shape optimisation problem by combining the knowledge from the previous chapters. Having introduced some useful notation in the first section, we describe the computation of the wave-making resistance in the second section. In the third section, we give a detailed presentation of the shape optimisation algorithm. The numerical results are discussed at the end of the chapter.

A general conclusion sums up this thesis and we give perspectives of future work.

Contents

Chapter I

The Neumann-Kelvin problem

In this chapter we introduce the Neumann-Kelvin problem, which will be our model for the calculation of the wave-making resistance. The goal is to predict the steady flow and wave pattern around an obstacle placed just below the free surface placed in a uniform stream (see figure I.1 below). Many results from [START_REF] Kuznetsov | Linear Water Waves[END_REF] will be cited extensively throughout the chapter as they will be particularly useful for us. In particular, we rely on [START_REF] Kuznetsov | Linear Water Waves[END_REF] for the existence results, and for the expression of the wave-making resistance.

The chapter is organised as follows: first, we will explain how the Neumann-Kelvin model can be obtained formally from the incompressible Navier-Stokes equation (section 1), then the calculation of the fundamental solution associated to our problem will be presented (section 2), which will lead to the boundary integral formulation of the problem (section) for which the existence results from [START_REF] Kuznetsov | Linear Water Waves[END_REF] and the idea begin their proof will be briefly presented (section 3). Then we will explain how the wave-making resistance can be calculated once the solution to the aforementioned boundary integral equation is known (section 4). Finally, in section 5 we will present the expression of boundary integral operators (weakly singular, strongly singular, hypersingular) related to our problem which will be useful for later purposes.

From Navier-Stokes to Neumann-Kelvin

Let us introduce the following notations (see figure I.1): Ω + will denote the fluid domain, Ω -will represent the interior of the obstacle. The boundaries of Ω + are the free surface S and the border of the obstacle, denoted Γ. The vector n = (n 1 , n 2 ) denotes the inwards unit normal vector of ∂ + Ω , and τ the tangent vector (oriented anticlockwise). We also denote S 0 as the straight horizontal line at the level of water far upstream.

Free-surface Navier-Stokes Equations

In their incompressible version, the Navier-Stokes equations read:

   ρ(∂ t V + (V • ∇)V ) = ∇ • σ -ρge 2 , ∇ • V = 0 , in Ω + (I.1)
where V : (t, x) → V (t, x) accounts for the fluid's velocity, and σ is the fluid's stress tensor, while ρ represents the fluid's density and g the acceleration of the gravity field. For a complete derivation of the incompressible Navier-Stokes system from the classical laws of conservation of mass, momentum and energy, we refer to [START_REF] Batchelor | An Introduction to fluid dynamics[END_REF]. The stress tensor contains two components:

• the Newtonian inner stress, proportional to the strain tensor,

• the hydrodynamic pressure, which value ensures the incompressibility everywhere and for all times (the pressure P : (t, x) → P (t, x) is hence an unknown of the problem).

The stress tensor hence reads :

σ = ν ∇V + t ∇V 2 + P Id , (I.2)
where ν stands for the fluid's dynamic viscosity. On the free surface S t , the normal stress vanishes:

σ n S = 0 , on S t , (I.3)
where n S is the outwards unit normal vector at the free surface. This condition is called the "dynamic condition" for the free-surface. Of course, in such water-wave phenomena, the free surface evolves in time. Let:

S t = {x(t, s) ; s ∈ R} , (I.4)
then x satisfies the following equation:

ẋ(t, •) = V (x(t, •)) , (I.5)
where ẋ is the Lagrangian derivative of x with respect to t. This equation is called the "kinematic condition" since it describes the motion of the interface, which is carried by the fluid's flow. The boundary condition on Γ have to ensure the impermeability of the obstacle, which reads:

V • n = 0 , on Γ . (I.6)
For the Navier-Stokes equation, an other condition for the tangent component of the velocity on the obstacle is required in order to avoid ill-posedness. The most common approach is to impose a so-called "no-slip" condition: V • τ = 0.

The perfect fluid hypothesis : Euler's equations

When the the viscosity is neglected (perfect fluid hypothesis), the stress tensor writes:

σ = P Id . (I.7)
By injecting the above stress tensor into the Navier-Stokes system (I.1), we obtain the freesurface Euler's equations :

   ρ(∂ t V + (V • ∇)V ) = ∇P -ρge 2 , in Ω + , ∇ • V = 0 , in Ω + , (I.8)
The dynamic boundary free surface condition yields :

σ n S = P n S = 0 , (I.9) which leads to the following boundary condition on the pressure: P = 0 , on S . (I.10)

The irrotational flow hypothesis : Bernoulli's equations

In this section we introduce the irrotational flow hypothesis. We will first show that this hypothesis is self-consistent, and then we will use this hypothesis to re-write the Euler equations and the boundary conditions in terms of a velocity potential. First, let us take the curl of the conservation of momentum equation in (I.8) :

∂ t ∇ × V + ∇ × ((V • ∇)V ) = 1 ρ ∇ × ∇(P -ρgx 2 ) = 0 , (I.11)
by denoting ω = ∇ × V , we remark that:

(V • ∇)V = 1 2 ∇|V | 2 -V × ω , (I.12) hence: ∇ × ((V • ∇)V ) = -∇ × (V × ω) , (I.13) = -((∇ • ω)V -(∇ • V )ω + (ω • ∇)V -(V • ∇)ω) . (I.14)
Since ∇ • V = 0 and ∇ • ω = 0, we have:

∇ × (V • ∇)V = -(ω • ∇)V + (V • ∇)ω . (I.15)
Substituting (I.15) into (I.11), we obtain the vorticity equation :

∂ t ω + (V • ∇)ω = (ω • ∇)V . (I.16)
From the above vorticity equation, and assuming ω and V are smooth, it is easy to prove that, if ω = 0 at a given time, then ω = 0 for all times. In the rest of the manuscript we will suppose that ω = 0. The Helmholtz-Hodge theorem states that any curl-free smooth vector field can be expressed as the gradient of a potential. Hence, since ω = 0, there exists a potential Ψ such that :

V = ∇Ψ , (I.17) since V is also divergence-free in Ω + , the potential solves the Laplace equation:

∇ • ∇Ψ = ∆Ψ = 0 , in Ω + . (I.18)
The momentum equation in (I.8), combined with (I.17), and the identity (I.12) yields:

∇ ∂ t Ψ + 1 2 |∇Ψ| 2 -1 ρ P + gx 2 = 0 , in Ω + , (I.19)
which leads to :

∂ t Ψ + 1 2 |∇Ψ| 2 -1 ρ P + gx 2 = c(t) , in Ω + , (I.20)
where c is some function of time. The equation (I.20) is often referred in literature as the equation of Bernoulli. In (I.20) the pressure is still unknown, except on S (where the dynamic free-surface condition P = 0 applies), hence taking the trace of the above expression onto S allows us to define the following dynamic free-surface condition for the potential:

∂ t Ψ + 1 2 |∇Ψ| 2 + gx 2 = c(t) , on S , (I.21)
Considering that Ψ is defined up to a function of time, we can take c independent of time without any harm:

∂ t Ψ + 1 2 |∇Ψ| 2 + gx 2 = c , on S , (I.22)
The impermeability condition on Γ writes:

∇Ψ • n = 0 , on Γ . (I.23)
In the following section, we introduce a graph representation of the free-surface in order to obtain a simpler formulation of the kinematic free-surface condition.

Graph representation of the free-surface

From now on, we will consider that S is represented with a function describing the election of the free surface with respect to a reference level (say, x 2 = 0). Hence:

S t = {(x 1 , η(t, x 1 )) ; x 1 ∈ R} . (I.24)
Let us consider an arbitrary point (x 1 (t), x 2 (t)) on the free surface, carried by the flow. From the kinematic free-surface condition we have:

d dt x 1 (t) = V 1 (t, x 1 (t), x 2 (t)) , (I.25) d dt x 2 (t) = V 2 (t, x 1 (t), x 2 (t)) . (I.26)
By the definition of η, we have x 2 (t) = η(t, x 1 (t)), hence:

V 2 (t, x 1 (t), η(t, x 1 (t))) = d dt η(t, x 1 (t)) = ∂ t η(t, x 1 (t)) + d dt x 1 (t) ∂ 1 η(t, x 1 (t)) , (I.27) = ∂ t η(t, x 1 (t)) + V 1 (t, x 1 (t), η(t, x 1 (t))) ∂ 1 η(t, x 1 (t)) .
(I.28)

This equation can be rewritten as the kinematic boundary condition:

∂ t η(t, x 1 ) + V 1 (t, x 1 , η(t, x 1 )) ∂ 1 η(t, x 1 ) = V 2 (t, x 1 , η(t, x 1 )) , (I.29)
where t ∈ R +, * and x 1 ∈ R. With this representation of the interface, the Bernoulli system of equations for water-waves writes:

                 ∂ t Ψ(t, x 1 , η(t, x 1 )) + 1 2 |∇Ψ(t, x 1 , η(t, x 1 ))| 2 = -gη(t, x 1 ) + c , for x 1 ∈ R , ∂ t η(t, x 1 ) + ∂ 1 Ψ(t, x 1 , η(t, x 1 )) ∂ 1 η(t, x 1 ) = ∂ 2 Ψ(t, x 1 , η(t, x 1 )) , for x 1 ∈ R , ∆Ψ = 0 , in Ω + , ∂ n Ψ = 0 , on Γ . (I.30)
Note that nothing ensures that the graph representation is consistent over time. The wellknown phenomenon of wave-breaking indicates that a graph representation is not always possible, at least in some conditions.

Up to now we haven't described the behaviour of the solution at infinity. Far away from the obstacle, the velocity is assumed to be uniform and horizontal:

V → U e 1 , for |x| → ∞ . (I.31)
Which can be translated as :

∇Ψ → U e 1 , for |x| → ∞ . (I.32)
For the sake of the presentation, we introduce the perturbation potential:

Φ = Ψ -U x 1 .
(I.33)

I.1 From Navier-Stokes to Neumann-Kelvin

The perturbation potential satisfies:

                         ∂ t Φ(t, x 1 , η(t, x 1 )) + 1 2 |∇Φ(t, x 1 , η(t, x 1 )) + U e 1 | 2 = -gη(t, x 1 ) + c , for x 1 ∈ R , ∂ t η(t, x 1 ) + (∂ 1 Φ(t, x 1 , η(t, x 1 )) + U ) ∂ 1 η(t, x 1 ) = ∂ 2 Φ(t, x 1 , η(t, x 1 )) , for x 1 ∈ R , ∆Φ = 0 , in Ω + , ∂ n Φ = -U n • e 1 , on Γ , |∇Φ| → 0 , for |x| → ∞ , (I.34)
with c = U 2 by consistency. In the following section we present a linearised version of the above model, for waves of small amplitude and impulsion.

The linearised equations of water-waves

Our goal is to expand formally all the terms of (I.34), at first order with respect to (η, Φ). We remark that: where S 0 = {(x 1 , x 2 ) ∈ R 2 ; x 2 = 0}. Injecting the above equation into (I.34), we obtain the following system by keeping only first order terms:

Φ(t, x 1 , η(t, x 1 )) = Φ(t, x 1 , 0) + η(t, x 1 )∂ 2 Φ(t, x 1 , 0) + o(|η|) , (I.
                         ∂ t Φ + U ∂ 1 Φ = -gη , on S 0 , ∂ t η + U ∂ 1 η = ∂ 2 Φ , on S 0 , ∆Φ = 0 , in Ω + , ∂ n Φ = -U n • e 1 , on Γ , |∇Φ| → 0 , for |x| → ∞ . (I.38)
The Neumann-Kelvin problem consists in finding a steady-state for (I.38). Formally, by setting ∂ t Φ = ∂ t η = 0 in the two first equations of (I.38), we obtain:

   U ∂ 1 Φ = -gη , on S 0 , U ∂ 1 η = ∂ 2 Φ , on S 0 , (I.39)
Taking the derivative with respect to x 1 in the first equation, and substituting the second equation, we get the free surface boundary condition as:

∂ 2 11 Φ + ν ∂ 2 Φ = 0 , on S 0 , (I.40)
where ν = g/U 2 . The others equations in (I.38) remain the same (except for the last one), as they do not contain any time derivative. However, the boundary condition |∇Φ| → 0 for |x| → ∞, contains a difficulty in the two-dimensional case because of the (unidirectional, non-damped) train of waves generated by the obstacle downstream. As explained in [START_REF] Kuznetsov | Linear Water Waves[END_REF], it is only possible to impose:

           |∇Φ| = O(1) , for |x| → ∞ , |∇Φ| → 0 , for x 1 → -∞ , |∇Φ| → 0 , for x 2 → -∞ . (I.41)
which is less than the previous condition, since we don't ensure |∇Φ| → 0 for x 2 → +∞ (downstream). This condition cannot guarantee the uniqueness of the solution, even up to a constant. A so-called radiation condition can be imposed downstream to recover uniqueness results, however, as stated in [START_REF] Kuznetsov | Linear Water Waves[END_REF] this condition cannot be completed until the problem is partly solved. One approach is to add some dissipation in the model, in order to damp the waves downstream and to keep the condition |∇Φ| → 0 for |x| → ∞:

                 ∂ 2 11 Φ + ν ∂ 2 Φ = ε ∂ 1 Φ , on S 0 , ∆Φ = 0 , in Ω + , ∂ n Φ = -U n • e 1 , on Γ , |∇Φ| → 0 , for |x| → ∞ , (I.42)
where ε accounts for the damping which is often called the "Rayleigh dissipation". One can solve this problem by defining the problem's fundamental solution, and then, taking ε → 0, we can recover the appropriate set of conditions at infinity by studying the behaviour of the solution far away from the obstacle. In the next section we present the calculation of the fundamental solution of the Neumann-Kelvin problem we have described above.

Fundamental Solution of the Neumann-Kelvin problem

In ship hydrodynamics, the fundamental solution, also known as Green's function, is an important tool both from the theoretical and from the numerical standpoint (see for instance [START_REF] Delhommeau | Les problèms de diffraction-radiation et de résistance de vagues: étude théorique et résolution numérique par la méthode des singularités[END_REF] and [START_REF] Rahman | Three dimensional Green's function for ship motion at forward speed[END_REF] for numerical simulations) as it the base ingredient for the boundary integral reformulation of the problem. For our simplified 2D problem, we rely on the same tools as those used in ship hydrodynamics since our long-term goal is to use the same methods for the actual 3D problem.

The main problem here is to obtain a Green's function that satisfies the linearised freesurface condition (I.40), as well as the right behaviour at infinity (logarithmic behaviour upstream and deep, and bounded downstream). This will prove to be very useful in the next section, since, though a single-layer representation of the solution, the only condition that will be left to impose will be (I.32), which only involves the boundary of the obstacle. The Green's function for the Neumann-Kelvin problem is already known in literature and appears in various equivalent forms: Havelock [START_REF] Havelock | The effect of speed of advance upon the damping of heave and pitch[END_REF], Wehausen and Laitone [START_REF] Wehausen | Surface waves. Encyclopedia of Physics[END_REF] presented the Green's function as a double integral, while Kuznetsov in [START_REF] Kuznetsov | Linear Water Waves[END_REF] used a single integral to express the Green's function. Rahman in [START_REF] Rahman | Three dimensional Green's function for ship motion at forward speed[END_REF] replaced the Green's function with double integral in [START_REF] Wehausen | Surface waves. Encyclopedia of Physics[END_REF] by a single integral using the complex exponential integral. On this particular topic, references in literature rarely give an insight on how this function is obtained. For this reason, in the following section, we detail the calculation of the Green's function.

Calculation of the Green's function

For the Neumann-Kelvin problem (I.38), the Green's function G y for a source y fully submerged (y ∈ R × R -, * ) should satisfy (at least):

     ∆G y = δ y , in R × R -, * ∂ 2 11 G y + ν ∂ 2 G y = 0 , on S 0 , (I.43)
where δ y is the Dirac distribution supported on {y}, and hence the first equation has to be understood in the sense of distributions. In the following, we introduce the following lifting function to get rid of the Dirac distribution on the right-hand side of the first equation of (I.43): G y (x) := where y = (y 1 , -y 2 ) is the symmetric of y with respect to the free surface. This function is known to satisfy ∆G y = δ y thanks to the first term (Green's function of the Laplace equation), while the second term ensures, by symmetry, that ∂ 2 G y = 0 on S 0 , which will be useful for later. Let us write G y using this lifting:

G y (x) = 1 2π ln(|x -y|) + 1 2π ln(|x -y|) + H y (x). (I.45)
The function H y that we need to obtain solves:

∆H y = 0, (I.46)
and:

∂ 2 11 H y (x 1 , 0) + ν∂ 2 H y (x 1 , 0) = -∂ 2 11 G y (x 1 , 0). (I.47)
since ∂ 2 G y = 0 on S 0 as we have seen earlier.

By taking the Fourier transform of (I.46) with respect to x 1 , we obtain

∆H y (k, x 2 ) = -k 2 Ĥy (k, x 2 ) + ∂ 2 22 Ĥy (k, x 2 ) = 0 , (I.48)
which is a second order differential equation with respect to the variable x 2 . This equation has solutions of the form:

Ĥy (k, x 2 ) = Ae |k|x 2 + Be -|k|x 2 . (I.49)
where A and B are real functions of k that we will determine below. First, for the Green's function to be bounded for x 2 → -∞, we require that B = 0, which leads to

Ĥy (k, x 2 ) = Ae |k|x 2 . (I.50)
In order to find A, we take the Fourier transform of (I.47) with respect to x 1 , we get:

-k 2 Ĥy (k, 0) + ν∂ 2 Ĥy (k, 0) = -F[∂ 2 11 G y (•, 0)](k), (I.51)
where F stands for the Fourier transform. We hence have to get the Fourier transform of:

∂ 2 11 G y (x 1 , 0) = 1 π -(x 1 -y 1 ) 2 + y 2 2 ((x 1 -y 1 ) 2 + y 2 2 ) 2 .
(I.52)

The Appendix B shows how this Fourier transform is obtained. We have:

F[∂ 2 11 G y (•, 0)](k) = |k|e -iky
Ĥy (k, x 2 ) = e -iky 1 +|k|(y 2 +x 2 ) |k| -ν (I.56)
In order to recover H y , we have to take the inverse Fourier transform of (I.56).

H y (x) = 1 2π R e -iky 1 +|k|(y 2 +x 2 )
|k| -ν e ikx 1 dk (I.57)

Let us denote:

u(k) = e -iky 1 +|k|(y 2 +x 2 ) |k| -ν e ikx 1 .
We remark that u(-k) = u(k), hence:

H y (x) = 1 2π R u(k)dk = 1 2π 2 Re R + u(k) dk (I.58) = 1 π Re R + e -ik(y 1 -x 1 )+k(y 2 +x 2 ) k -ν dk. (I.59)
Let us denote: ω(x, y) = (y 1 -x 1 ) + i(x 2 + y 2 ). (I.60)

We obtain

H y (x) = 1 π Re R + e -ikω(x,y) k -ν dk. (I.61)
It is clear that the integral (I.61) is an improper integral. This problem is related to the lack of decay of the solution at infinity downstream. It is however possible to give a meaning in [START_REF] Kuznetsov | Linear Water Waves[END_REF] to (I.61) by avoiding the singularity from below, with a half-circle of center ν and radius (see figure I.2), where is meant to vanish:

H y (x) = 1 π lim →0 Re l - e -ikω(x,y) k -ν dk. (I.62)
where l -is a complex path of integration described in I.2.

Remark 1. The choice we made for the path of integration might seem arbitrary, and there could be lots of other ways to give a meaning to (I.61) with a limiting process. For instance we could consider the Cauchy principal value of this integral, or avoiding the singularity from above. We will see later that the choice we have made will imply that the Green's function we obtain satisfies physically sound conditions at infinity (both upstream and downstream). Moreover, taking the "below" path can be shown to be equivalent to consider a vanishing Raleigh dissipation in the model (see Appendix A).

The formula (I.62) gives us an expression for the Green's function, but it is rather impractical for the numerical point of view, as it involves the limit in ε. Fortunately, as it will be shown in the next subsection, this limit can be related to the exponential integral function, which is well known, and for which series expansions exists.

Relation with the exponential integral function

The exponential integral function we consider is defined by:

E 1 (a) = R + + a e -z z dz, (I.63)
where the a ∈ C, and hence the path of integrations goes from a to infinity parallel to the real axis (see figure I.2). This function is well known, from its behaviour at infinity to numerical approximations (see for instance [START_REF] Olver | Asymptotics and special function[END_REF]). [START_REF] Bai | A variational method in potential flow with a free surface[END_REF] -Paths of integration for z → e -z /z in (I.65) and (I.66). On the left: Im(-iνω > 0), which implies that y 1 < x 1 (the point of observation x is downstream of the source y). On the right: Im(-iνω < 0), which implies that y 1 > x 1 (the point of observation x is upstream of the source y).

R iR -iνω -iνω + R P ε R iR -iνω -iνω + R P ε Figure I.
Let us show that (I.62) can be rewritten as the integral of e -z /z on a certain path P (see figure I.3) which is built from P = {iω(k -ν); k ∈ R + }, and avoids the singularity (at 0) from below:

P e -z z dz = R + e -iω(k-ν) iω(k -ν) iω dk = e iων R + e -iωk k -ν dk. (I.64)
Now let us rewrite the integral over P as an integral on -iνω + R by making use of the Cauchy's integral theorem. As seen on figure I.3, there are two cases:

• If (y 1 -x 1 ) < 0 (downstream, see figure I.3 left), we close the path with an arc going from P to -iνω + R, and use the decay at infinity for the positive real part of e -z /z, and we remark that the path encloses the singularity at z = 0. The Cauchy's integral theorem yields:

P e -z z dz = -iνω+R e -z z dz + 2πi Res(z → e -z z , 0) = E 1 (-iνω) + 2πi (I.65) • If (y 1 -x 1 ) > 0 (upstream, see figure I.3 right)
, by the same arguments, and noticing that this time the integral is outside of the interior of the path of integration, we have:

P e -z z dz = E 1 (-iνω) (I.66)
Summing up (I.64), (I.65), and (I.66), we have:

H y (x) = 1 π    Re{e -iνω E 1 (-iνω)}, if: y 1 > x 1 Re{e -iνω [E 1 (-iνω) + 2πi]}, if: y 1 < x 1 (I.
G y (x) = 1 2π ln(|x -y|) + 1 2π ln(|x -y|) + 1 π        Re{e -iνω E 1 (-iνω)}, if: y 1 > x 1 Re{e -iνω [E 1 (-iνω) + 2πi]}, if: y 1 < x 1 (I.68)
In the rest of the chapter, we will denote G(x, y) = G y (x) for the sake of presentation.

Asymptotic behaviour of the Green's function

In this section we study the behaviour of the Green's function (I.68) far away from the source. This will show that the choice of path of integration we have made to give sense to (I.61) ensures that waves are propagated downstream. We first recall the following asymptotic formula, which can be found in [START_REF] Olver | Asymptotics and special function[END_REF]:

e z E 1 (z) = O 1 |z| , for z ∈ R -. (I.69)
This leads, along with (I.68) to the asymptotic formula:

G y (x) = r(x, y) + 1 π ln(|x -y|) + 2        0 if: y 1 > x 1 e ν(x 2 +y 2 ) sin(ν(y 1 -x 1 )) if: y 1 < x 1 (I.70)
where :

r = O(1/|x -y|) and ∇r = O(1/|x -y| 2 )
. This formula can also be found in [START_REF] Kuznetsov | Linear Water Waves[END_REF], albeit in a different form. From (I.70), we remark that the asymptotic behaviour of the Green's function downstream (i.e. for y 1 < x 1 ) is periodic, with a period of 1/ν. The figure I.4 shows a plot of the Green's function (I.68), along with the asymptotic approximation (I.70) (neglecting the rest r), and their difference.

Boundary integral formulation of the problem

Let us recall the Neumann-Kelvin problem:

                                 ∆Φ = 0 , in Ω + , ∂ 2 11 Φ + ν ∂ 2 Φ = 0, on S 0 , ∂ n Φ = -U n • e 1 , on Γ , |∇Φ| = O(1) , for |x| → ∞ , |∇Φ| → 0 , for x 1 → -∞ (unif. w.r.t. x 2 ),
and The idea behind the boundary integral formulation of such a problem is to replace the boundary condition (I.73) on Γ by a source term in (I.71), considered on the whole half-space R × R -, consisting in a single layer distribution carried by Γ:

x 2 → -∞ (unif. w.r.t.
< sl α , ϕ > D ,D = Γ α(x s ) ϕ(x s ) ds ,
where α ∈ C(Γ) is the (yet unknown) weight function associated with the distribution of sources. Since the support of this distribution is Γ, solving ∆Φ = sl α on R × R -, will imply that ∆Φ = 0 on both Ω + and Ω -. When x ∈ Ω + ∪ Ω -, we obtain a solution of (I.71) and (I.72):

Φ(x) = Γ G(x, y) α(y)ds y , (I.76)
where G is the Green's function previously defined in section 2, designed in such a way that it solves (I.72). The unknown α should be determined in such a way that (I.73) is satisfied. Note that ∂ n Φ in (I.73) has to be understood as the trace of ∇Φ • n onto Γ from the exterior (we will denote it γ + Γ [∇Φ]). Using (I.76), the condition (I.73) becomes:

γ + Γ z → n(z) • ∇ z Γ G(z, y) α(y) ds y (x) = -U n(x) • e 1 , ∀x ∈ Γ . (I.77)
Furthermore, it is possible to show that (see [46] for example), when Γ is of class C 2 :

γ + Γ z → n(z) • ∇ z Γ G(z, y) α(y) ds y (x) = Γ ∂ n(x) G(x, y) α(y) ds y - 1 2 α(x) .
(I.78)

I.3 Boundary integral formulation of the problem

which leads to the following boundary integral equation for α:

1 2 α(x) - Γ ∂ n(x) G(x, y) α(y) ds y = U n(x) • e 1 , ∀x ∈ Γ . (I.79) Remark 2.
In principle, the singularity of the kernel ∂ n(x) G(x, y) which is of the type :

∂ n(x) G(x, y) = (x -y) • n(x) |x -y| 2 + o(1) in the limit y → x
should be a priori non integrable (like y → 1/y on R) and hence the integral should be understood as a Cauchy principal value. However, in two dimension of space (as it is the case here) it turns out that when the boundary is C 2 , we have (see [START_REF] Kublik | An implicit interface boundary integral method for Poisson's equation on arbitrary domains[END_REF]):

(x -y) • n(x) |x -y| 2 = 1 2 H(x) + o(d Γ (x, y)) , (I.80)
where H(x) denotes the curvature of Γ at the point x. Hence the kernel in (I.79) is in fact continuous in our case and the integral has to be understood in the usual sense.

Asymptotic behaviour of solutions

The following theorem, found in [START_REF] Kuznetsov | Linear Water Waves[END_REF] ensures that a solution defined by (I.76) satisfies the conditions (I.74)-(I.75), when α solves (I.79).

Theorem 1. Let Γ be a C 2 boundary, let Φ be defined given by (I.76), with α ∈ C(Γ). Then we have:

Φ(x) = c + R(x) + Q ln(|x|) +    0 if: x 1 < 0 Im C e -iνx 1 +νx 2 if: x 1 > 0 (I.81) where R(x) = O(1/|x|) and ∇R(x) = O(1/|x| 2 )
, and:

Q = 1 π Γ α(y) ds (I.82) C = 2 Γ α(y)e iνy 1 +νy 2 ds (I.83)
The proof uses the asymptotic behaviour of the Green's function. Combining (I.76) and (I.70), we have:

Φ(x) = Γ r(x, y)α(y) ds y + 1 π Γ ln(|x -y|)α(y) ds y + 2        0 if: x 1 < 0 Γ e ν(x 2 +y 2 ) sin(ν(y 1 -x 1 ))α(y) ds y if: x 1 > 0 (I.84)
From the continuity of α, and the behaviour of r at infinity, he have:

Γ r(x, y)α(y) ds y = R 1 (x) (I.85) where R 1 (x) = O(1/|x|) and ∇R 1 (x) = O(1/|x| 2 ).
The second term writes:

Γ ln(|x -y|)α(y) ds y = ln(|x|) Γ ln(|x -y|) ln(|x|) α(y) ds y (I.86) = ln(|x|) Γ α(y) ds y + Γ ln |x -y| |x| α(y) ds y (I.87)
Elementary calculations show that, for y given in the lower half plane: This comes from the following identity, that can be found in [46], which is a consequence

ln |x -y| |x| = O 1 |x| , (I.88) ∇ ln |x -y| |x| = O 1 |x| 2 , (I.89) hence: Γ ln(|x -y|)α(y) ds y = Q ln(|x|) + R 2 (x) (I.90) where R 2 (x) = O(1/|x|) and ∇R 2 (x) = O(1/|x| 2 ). Finally: Γ e ν(x 2 +y 2 ) sin(ν(y 1 -x 1 ))α(y) ds y = Im 2 Γ e iν(y 1 -x 1 )+ν(x 2 +y 2 ) α(y) ds y , (I.91) = Im C e -iνx 1 +νx 2 . (I.92) Stating that R = R 1 + R 2 completes

I.3 Boundary integral formulation of the problem

of the jump relations for solutions of Laplace equation represented with single or double layer potentials:

Γ ∂ n(x) G(x, y) ds x =            -1 if y ∈ Ω -, -1/2 if y ∈ Γ , 0 if y ∈ Ω + . (I.93)
Integrating (I.79) over Γ leads to:

Γ 1 2 α(x) - Γ ∂ n(x) G(x, y)α(y) ds y ds x = Γ U n(x) • e 1 ds x = 0 , (I.94)
from the flux-divergence theorem. The left-hand side can be rewritten as:

Γ 1 2 α(x) - Γ ∂ n(x) G(x, y)α(y) ds y ds x = 1 2 Γ α(x) ds x - Γ α(y) Γ ∂ n(x) G(x, y) ds x ds y ,
(I.95) by exchanging the order of integration on the second term. Now, using (I.93), we have:

Γ 1 2 α(x) - Γ ∂ n(x) G(x, y)α(y) ds y ds x = Γ α(x) ds x , (I.96)
hence, from (I.94), we obtain Q = 0.

Existence of solutions and uniqueness

In this section we recall the existence and uniqueness results that were obtained in [START_REF] Kuznetsov | Linear Water Waves[END_REF] for the Neumann-Kelvin problem, and give a brief account of the main ideas behind the proofs of these results. We recall that the mathematical formulation of the Neumann-Kelvin problem reads: find

Φ ∈ C 2 (Ω + ) such that                    ∆Φ = 0 , in Ω + , ∂ 2 11 Φ + ν ∂ 2 Φ = 0 , on S 0 , ∂ n Φ = -U ∞ n • e 1 , on Γ , sup Ω + |∇Φ| < ∞ and |∇Φ| → 0 as x 1 → -∞ . (I.97)
This problem is extensively studied in [START_REF] Kuznetsov | Linear Water Waves[END_REF] for a general right-hand side f instead of -U ∞ n • e 1 in the third equation above. A central idea is to relate this problem to the integral formulation (I.79). We sum up the main results and ideas.

Theorem 2. Assume that Ω -is a bounded and simply connected domain of the lower half-plane with C 3 boundary Γ, and such that Ω -is included in the (open) lower half-plane. Then the integral equation (I.79) has a unique solution α ∈ C 0 (Γ) for all ν > 0 with a possible exception for a finite number of values. Moreover, α belongs to C 1,θ (Γ) for all θ ∈ (0, 1), and the Neumann-Kelvin problem (I.97) is solvable for all ν > 0 with a possible exception for a finite number of values.

We first note that the function -U ∞ n • e 1 belongs to C 2 (Γ), and we can apply the results in [START_REF] Kuznetsov | Linear Water Waves[END_REF]. Any solution α ∈ C 1,θ (Γ) of (I.79) provides a solution Φ of the Neumann-Kelvin problem (I.97) thanks to formula (I.76). We note that at this point, this solution Φ of (I.97) may not be unique. Now, the idea is that (I.79) is a Fredholm equation of the second kind that we denote:

(Id + T ν )α = f . (I.98)
First, It can be shown that this Fredholm equation has a unique solution in C 0 (Γ) for ν = 0 and ν = ∞, as it consist in solving the Laplace equation with respectively a Dirichlet and a Neumann boundary condition at S 0 , as it can be observed formally on (I.97). By continuity, for ν > 0 small enough and ν > 0 large enough, (I.98) is uniquely solvable. Moreover, the Fredholm operator (Id + T ν ) depends analytically on ν with values in L(C 0 (Γ)), the Banach space of continuous linear operators from C 0 (Γ) into itself. By the principle of isolated zeros, combined with the Fredholm alternative, the unique solvability of (I.98) can be lost for a finite values of ν ∈ (0, ∞) at most.

Uniqueness of solution to the Neumann-Kelvin problem (I.97) can be proved in some cases. In particular, we have : Theorem 3. Let the assumptions of Theorem 2 hold, and assume moreover that x • n ≥ 0 on Γ. Then problem (I.97) and the integral equation (I.79) are equivalent in the following sense. Any solution of (I.97) has a unique representation in the form (I.76) (up to a constant term), where α satisfies (I.79). Conversely, subsituting any solution α of (I.79) into (I.76), one obtains a solution of (I.97).

The geometric condition x • n ≥ 0 on Γ is satisfied for instance if Γ is split by the axis x 1 = 0 into two parts that can be specified explicitly: x 1 = g ± (x 2 ). In the case of a ball, it has also been shown that the Neumann-Kelvin problem (I.97) has a unique solution (up to a constant term) for all values of ν.

The wave-making resistance

In this section, we explain how the wave-making resistance is determined from the knowledge of the solution of the Neumann-Kelvin problem. Moreover, by making use of the behaviour of the potential far away from the obstacle, we will introduce the so-called "energy method" for the calculation of the wave-making resistance.

Integration of the normal stress

The forces exerted by the fluid on the obstacle is the result of the integration of the normal stress on the surface of the obstacle:

F = - Γ σ n ds . (I.99)
As we have seen in section 1, we have σ = P Id, hence:

F = - Γ P n ds . (I.100)
This force has two components: the horizontal component is called the drag or resistance and the vertical component is called the lift. Our interest focus on the resistance, which is the quantity we will seek to minimise later. Since all these forces are zero in the absence of a free-surface (because of d'Alembert's paradox, see [START_REF] Batchelor | An Introduction to fluid dynamics[END_REF]), the resistance is only caused by the production of waves behind the obstacle. For this reason, the horizontal force will be called the wave-making resistance:

R w = - Γ P n 1 ds , (I.101)
where n 1 = n • e 1 . Using Bernoulli's formula (I.20), at steady-state, we have:

P = c + ρgx 2 -ρ |V | 2 2 . (I.102)
Hence:

R w = - Γ c + ρgx 2 -ρ |V | 2 2 n 1 ds . (I.103)

Chapter I. The Neumann-Kelvin problem

First, we remark that:

R w = - Γ (c + ρgx 2 )n 1 dΓ = - Γ   c + ρgx 2 0   • n ds (I.104) = - Ω - ∇ •   c + ρgx 2 0   dx = 0 . (I.105)
Finally:

R w = ρ Γ |V | 2 2 n 1 ds . (I.106)
Describing the velocity in terms of the perturbation potential V = ∇Φ + U e 1 , we have:

R w = ρ Γ |∇Φ| 2 2 + U ∂ 1 Φ n 1 ds . (I.107)
Now we will show that the wave-making resistance can be expressed as a function of the amplitude of the waves downstream rather than integrating the pressure on the boundary.

The energy method

First, we give a description of the behaviour of the solution of the Neumann-Kelvin problem at infinity in all directions. Let Φ be a solution of the Neumann-Kelvin problem, then we have, as a consequence of Theorem 1:

∇Φ(x 1 , •) = O(|x| -2 ) uniformly as x 1 → -∞ , (I.108) ∇Φ(•, x 2 ) = O(|x| -2 ) uniformly as x 2 → -∞ , (I.109) ∇Φ(x 1 , x 2 ) = ν Re   C   -i 1   e -i(νx 1 + π 2 )+νx 2   + O(|x| -2 ) uniformly as x 1 → +∞ , (I.110)
where the complex amplitude reads:

C = 2 Γ Φ∂ n E -∂ n ΦE ds , (I.111)
with:

E(x) = e ν(ix 1 +x 2 ) . (I.112)
Following [START_REF] Kuznetsov | Linear Water Waves[END_REF], we start by showing that the wave-making resistance can be calculated using the asymptotic behaviour of the solution described above. Following [START_REF] Kuznetsov | Linear Water Waves[END_REF], we intro-

I.4 The wave-making resistance

duce the rectangle:

R a,b =] -a, a[×] -b, 0[ ,
with a and b large enough for the rectangle to enclose the obstacle Ω -. On one hand, we have:

R a,b ∩Ω + ∆Φ ∂ 1 Φ dx = 0 , (I.113)
on the other hand, integrating by parts, we have:

R a,b ∩Ω + ∆Φ ∂ 1 Φ dx = ∂(R a,b ∩ Ω + ) ∂ ñΦ ∂ 1 Φ ds - R a,b ∩Ω + ∇Φ • ∂ 1 ∇Φ dx , (I.114) = ∂(R a,b ∩ Ω + ) ∂ ñΦ ∂ 1 Φ ds - R a,b ∩Ω + ∂ 1 |∇Φ| 2 2 dx , (I.115) = ∂(R a,b ∩ Ω + ) ∂ ñΦ ∂ 1 Φ - |∇Φ| 2 2 ñ1 ds , (I.116)
where ñ is the outwards unit normal vector of R a,b ∩ Ω + . The boundary ∂(R a,b ∩ Ω + ) is the union of Γ and the border ∂R a,b of the rectangle. Note that on Γ, we have ñ = -n (as it was defined before). We deduce from the above that:

Γ ∂ n Φ ∂ 1 Φ - |∇Φ| 2 2 n 1 ds = ∂R a,b ∂ ñΦ ∂ 1 Φ - |∇Φ| 2 2 ñ1 ds (I.117)
The left-hand side writes:

Γ ∂ n Φ ∂ 1 Φ - |∇Φ| 2 2 n 1 ds = Γ -U n 1 ∂ 1 Φ - |∇Φ| 2 2 n 1 ds = -R w , (I.118) hence: R w = ∂R a,b |∇Φ| 2 2 ñ1 -∂ ñΦ ∂ 1 Φ ds . (I.
R w = T + B + L + R. First, we have: B = b a,b |∇Φ| 2 2 ñ1 -∂ ñΦ ∂ 1 Φ ds = a -a ∂ 2 Φ(x 1 , -b) ∂ 1 Φ(x 1 , -b) dx 1 , (I.120)
and:

L = l a,b |∇Φ| 2 2 ñ1 -∂ ñΦ ∂ 1 Φ ds = 1 2 0 -b |∂ 2 Φ(-a, x 2 )| 2 -|∂ 1 Φ(-a, x 2 )| 2 dx 2 . (I.121)
Thanks to the asymptotic formula (I.108)-(I.109), we can see easily that the two terms above vanish as a → ∞ and b → ∞. Let us now investigate the other terms:

T = t a,b |∇Φ| 2 2 ñ1 -∂ ñΦ ∂ 1 Φ ds = a -a -∂ 2 Φ(x 1 , 0) ∂ 1 Φ(x 1 , 0) dx 1 . (I.122)
Using the free-surface boundary condition, we obtain:

T = a -a ν -1 ∂ 2 11 Φ(x 1 , 0) ∂ 1 Φ(x 1 , 0) dx 1 = 1 2ν ∂ 1 Φ(a, 0) 2 -∂ 1 Φ(-a, 0) 2 , (I.123)
Using (I.108)-(I.110), we get:

T = ν 2 Re C e -iνa 2 + O(a -1 ) . (I.124)
The remaining term writes:

R = r a,b |∇Φ| 2 2 ñ1 -∂ ñΦ ∂ 1 Φ ds = 1 2 0 -b |∂ 2 Φ(a, x 2 )| 2 -|∂ 1 Φ(a, x 2 )| 2 dx 2 . (I.125)
Again, using the asymptotic formula (I.110), we get:

R = ν 2 2 Re C e -iνa-i π 2 2 -Re C e -iνa 2 0 -b e 2νx 2 dx 2 + O(a -1
) .

(I.126)

The limit b → ∞ yields:

R = ν 4 Re C e -iνa-i π 2 2 -Re C e -iνa 2 + O(a -1 ) . (I.127)
By summing T , B, L and R we have:

R w = ν 4 Re C e -iνa-i π 2 2 + Re C e -iνa 2 + O(a -1 ) , (I.128) = ν 4 Im C e -iνa 2 + Re C e -iνa 2 + O(a -1 ) (I.129) = ν 4 C e -iνa 2 + O(a -1 ) = ν 4 |C| 2 + O(a -1 ) (I.130)
Hence, taking the limit a → ∞, we get that:

R w = ν 4 |C| 2 . (I.131)
This formula shows that the wave-making resistance is proportional to the square of the amplitudes of the waves generated by the obstacle downstream. Recalling the value of C, from (I.83), we obtain that:

R w = ρ ν 4 Γ α Γ (x) E(x)ds x 2 .
(I.132)

Traces of some boundary integral operators

In the following section we detail the expression of several boundary integral operators which will be useful for later purposes. The proofs are not given in detail here, these results being well known for the Green's function of the Laplace equation [START_REF] Hsiao | Boundary Integral Equations[END_REF]. We only show how to adapt these results in our case. Let us denote, for x ∈ Ω + ∪ Ω -:

S[u](x) = Γ G(x, y)u(y) ds y , (I.133) D 1 [u](x) = Γ ∂ n(x) G(x, y)u(y) ds y , D 2 [u](x) = Γ ∂ n(y) G(x, y)u(y) ds y , (I.134) H 1 [u](x) = Γ ∇ 2 xx G(x, y)n(x) • n(x) u(y) ds y , H 2 [u](x) = Γ ∇ 2 xy G(x, y)n(y) • n(x) u(y) ds y , (I.135)
where

∇ 2 xx = ∇ x ∇ T x and ∇ 2 xy = ∇ x ∇ T y .
Let us now define the traces from Ω + of these functions as operators on u:

S[u] = γ + Γ [z → S(u)(z)] , (I.136) D 1 [u] = γ + Γ [z → D 1 (u)(z)] , D 2 [u] = γ + Γ [z → D 2 (u)(z)] , (I.137) H 1 [u] = γ + Γ [z → H 1 (u)(z)] , H 2 [u] = γ + Γ [z → H 2 (u)(z)] . (I.138)
The kernel in S being weakly singular, we have classically:

S[u](x) = Γ G(x, y)u(y) ds y . (I.139)
The singularity of G(x, y) comes from the term (2π) -1 ln(|x -y|) in (I.68), so, exactly as in the Laplace equation (see [46] or [START_REF] Hsiao | Boundary Integral Equations[END_REF]), we get :

D 1 [u](x) = Γ ∂ n(x) G(x, y)u(y) ds y - 1 2 u(x) , (I.140) D 2 [u](x) = Γ ∂ n(x) G(x, y)u(y) ds y + 1 2 u(x) . (I.141)
Let us now calculate the traces on Γ, from Ω + of:

H 1 [u](x) = Γ ∇ 2 xx G(x, y)n(x) • n(x) u(y) ds y , (I.142) H 2 [u](x) = Γ ∇ 2 xy G(x, y)n(y) • n(x) u(y) ds y . (I.143)
Using the following identity on 2 × 2 matrices:

Au • v = [tr(A)Id -A T ]Ru • Rv (I.144)
where R is a π/2 rotation matrix and tr(A) is the trace of the matrix A, we obtain:

∇ 2 xx G(x, y)n(x) • n(x) = tr(∇ 2 xx G(x, y))Id -∇ 2 xx G(x, y) τ (x) • τ (x) , (I.145) ∇ 2 xy G(x, y)n(y) • n(x) = tr(∇ 2 xy G(x, y))Id -∇ 2 yx G(x, y) τ (y) • τ (x) . (I.146)
The Green's function (I.68) of the Neumann-Kelvin problem has the form:

G(x, y) = g(x -y) + h(x -y) , (I.147)
where ∆g = 0 and ∆h = 0 on R 2 \ {0}. Firstly, we have:

tr(∇ 2 xx G(x, y)) = ∆g(x -y) + ∆h(x -y) = 0 , (I.148)
since x = y and x = y, for x ∈ Ω + . Secondly, we have:

tr(∇ 2 xy G(x, y))Id -∇ 2 yx G(x, y) = -∇ 2 xy (g(x -y) -h(x -y)) . (I.149)
This leads to :

∇ 2 xx G(x, y)n(x) • n(x) = -∇ 2 xx G(x, y)τ (x) • τ (x) , (I.150) ∇ 2 xy G(x, y)n(y) • n(x) = -∇ 2 xy G * (x, y)τ (y) • τ (x) . (I.151)
where G * (x, y) = g(x -y) -h(x -y). On one hand, we have:

∇ 2 xy G * (x, y)τ (y) • τ (x) = ∂ τ (x) ∂ τ (y) G * (x, y) (I.152)
On the other hand:

∂ τ (x) ∂ τ (x) G(x, y) = ∇ x (∇ x G(x, y) • τ (x)) • τ (x) (I.153) = (∇ x ∇ T x G(x, y) τ (x) + ∇ x τ (x) T ∇ x G(x, y)) • τ (x) (I.154) = ∇ 2 xx G(x, y) τ (x) • τ (x) + ∇ x G(x, y) • (∇ x τ (x) T ) T τ (x) . (I.155)
We have:

(∇τ T ) T τ =   ∂ 1 τ 1 τ 1 + ∂ 2 τ 1 τ 2 ∂ 1 τ 2 τ 1 + ∂ 2 τ 2 τ 2   , (I.156) =   ∂ 1 n 2 n 2 -∂ 2 n 2 n 1 -∂ 1 n 1 n 2 + ∂ 2 n 1 n 1   , (I.157) = 1 2 ∇(|n| 2 ) -∇ • n n = -∇ • n n = -H n , (I.158)
where we recall that H denotes the curvature of Γ. This leads to :

∇ 2 xx G(x, y) τ (x) • τ (x) = ∂ τ (x) ∂ τ (x) G(x, y) + H(x) ∇ x G(x, y) • n(x) . (I.159)
Gathering (I.152) on one hand and (I.159) on the other hand, we obtain, for x ∈ Ω + :

H 1 [u](x) = -∂ τ (x) ∂ τ (x) Γ G(x, y) u(y) ds y -H(x) Γ ∂ n(x) G(x, y) u(y) ds y (I.160) H 2 [u](x) = -∂ τ (x) Γ ∂ τ (y) G * (x, y) u(y) ds y = ∂ τ (x) Γ G * (x, y) ∂ τ u(y) ds y . (I.161)
since Γ is a closed smooth curve. Taking the trace of both the expression above leads to, for x ∈ Γ:

H 1 [u](x) = -∂ τ (x) ∂ τ (x) Γ G(x, y) u(y) ds y -H(x) Γ ∂ n(x) G(x, y) u(y) ds y + 1 2 H(x)u(x) , (I.162)
and:

H 2 [u](x) = ∂ τ (x) Γ G * (x, y) ∂ τ u(y) ds y . (I.163)
Hence the form:

H 1 [u] = -∂ 2 τ,τ S[u] -H D 1 [u] , (I.164) H 2 [u] = ∂ τ S[∂ τ u] , (I.165)
where:

S[u](x) = Γ G * (x, y)u(y) ds y . (I.166)
Chapter II

Shape Optimisation

In this chapter, we will give a presentation of the boundary variation method for shape optimisation along with its applications to the wave-making resistance problem. The boundary variations method consists in studying the influence (at first order) of small deformations of a domain (bounded open set with a smooth boundary) on a quantity (e.g. the surface area, the perimeter,...) that depends on the shape of the aforementioned domain. Typically, for applications in physics and engineering, the quantity of interest depends on the solution of a partial differential equation which, in turn, is posed on the domain we wish to optimise. Our general presentation will follow the main ideas that can be found in [START_REF] Allaire | Conception optimale de structures[END_REF] and [38], and we will briefly recall some of the results that can be found in these two textbooks. As stated in [38], the first occurence of this method can be traced back to the 1907 article of Jacques Hadamard called Mémoire sur le problème d'analyse relatif á l'équilibre des plaques élastiques encastrées [31]. Remarkably, in this long article, the problem is studied using a boundary integral formulation of involving a Green's function, exactly as it will be the case here.

This chapter is organised as follows: in the first section, we will recall the principle of the boundary variation method, and introduce the main notations and the concept of "shape derivative" and "shape gradient", on which we will rely for the future algorithm. In the section 2 we will present some useful lemmas on the calculation of various shape derivatives, which we will use as building blocks for our application. The Lagrangian method which allows an easy formal calculation of the shape gradient of a quantity depending on the solution of a partial differential equation posed on the variable domain will be presented in section 3. After recalling how the method works for partial differential equations formulated as variational problems, we propose an adaptation for the case of boundary integral equations. Finally, the last section of this chapter is dedicated to the calculation of the shape gradient of the wave-making resistance as it was defined in chapter I.

Presentation of the boundary variation method

Le us now consider a problem of minimisation:

min Ω∈O J(Ω) , (II.1)
where O is the set of open, non-empty bounded domains with a smooth boundary. In general, for optimisation, it is useful to define derivatives of the objective function with respect to variable of interest. First it provides a necessary optimality condition (the derivative of the optimal function at the optimal value has to be zero), and also, in some cases, this derivative can be used in optimisation algorithms such as the gradient descent method or Newton's method. Let us recall the notion of Fréchet derivative of a function in Banach space:

Definition 1 (Fréchet differentiability, Fréchet derivative). Let E and F be two Banach spaces and V an open subset of E. Let us consider f , a continuous mapping from V to F . We say that f is Fréchet differentiable at x ∈ V , if and only if there exists a linear and continuous mapping g x : E → F , called the derivative of f at x, satisfying:

f (x + h) -f (x) -g x (h) F = o( h E ) (II.2)
where:

lim h E →0 o( h E ) h E = 0 (II.3)
The set of admissible domains involved in our problem is not adapted to this definition, since it is not even a vector space. The idea behind the boundary variations method is to restrict the admissible domains to domains that can be obtained through the image of a reference domain Ω 0 with a smooth and reversible (with a smooth inverse) deformation. The classical presentation of the method in [START_REF] Allaire | Conception optimale de structures[END_REF] makes use of the space W 1,∞ to define diffeomorphisms:

O ad = T (Ω 0 ) ; (T -Id) ∈ W 1,∞ (R 2 , R 2 ) and (T -1 -Id) ∈ W 1,∞ (R 2 , R 2 ) , (II.4) where W 1,∞ (R 2 , R 2 ) denotes the set of L ∞ (R 2 , R 2 )
functions with L ∞ first order distributional derivatives, which is a Banach space when equipped with the norm:

u 1,∞ = |k|≤1 D k u ∞ . (II.5)
In (II.4), T will be called the placement function, and θ = T -Id, will be called the displace-

II.1 Presentation of the boundary variation method

ment field (see figure II.1).

x T (x) θ(x) The set of admissible shape (II.4) can be rewritten as:

Ω 0 T (Ω 0 )
O ad = (Id + θ)(Ω 0 ) ; θ ∈ W 1,∞ (R 2 , R 2 ) and ((Id + θ) -1 -Id) ∈ W 1,∞ (R 2 , R 2 ) , (II.6)
Now that the set of domains is parametrised with the function θ, it is natural to rewrite our objective function in terms of θ:

J((Id + θ)(Ω 0 )) = JΩ 0 (θ) , (II.7)
and to compute derivatives with respect to θ. Unfortunately, the set of smooth deformation fields that induce a smooth and reversible transformation is not a vector space either. For instance, if we define θ(x 1 , x 2 ) = (0, -1 2 x 2 ), then we have θ ∈ W 1,∞ (R 2 , R 2 ), also (Id + θ) is inversible and:

((Id + θ) -1 -Id)(x 1 , x 2 ) = (0, x 2 ) (II.8) which belongs to θ ∈ W 1,∞ (R 2 , R 2 )
. Now, if we consider the transformation induced by 2θ, we obtain that:

(Id + 2θ)(x 1 , x 2 ) = (x 1 , 0) , (II.9)
which is non-inversible. Luckily, the following lemma, that can be found in [START_REF] Allaire | Conception optimale de structures[END_REF] ensures that if a given smooth deformation field is small enough, then it induces a bijective transformation with a smooth inverse:

Lemma 1. Let θ ∈ W 1,∞ , satisfying θ W 1,∞ < 1 then, (Id + θ) is bijective and ((Id + θ) -1 -Id) ∈ W 1,∞ (R 2 , R 2 ).
The consequence of this lemma is that, on a neighbourhood of 0, it is sufficient to consider θ ∈ W 1,∞ (R 2 , R 2 ), which equipped with the norm (II.5) is a Banach space. Hence it is possible to define the derivability of JΩ 0 , at least at θ = 0: Definition 2 (shape derivative). We say that a function J of the domain is shape differentiable at Ω 0 ∈ O, when JΩ 0 defined in (II.7) is Fréchet differentiable at 0. We call the shape derivative the linear and continuous form J Ω 0 of W 1,∞ satisfying:

JΩ 0 (θ) = J(Ω 0 ) + J Ω 0 (θ) + o( θ W 1,∞ ) .
(II.10)

We denote the shape derivative J (Ω 0 ) = J Ω 0 (0). Also, from now on, we will denote:

θ W 1,∞ = ||θ||.
The next proposition, that can also be found in [START_REF] Allaire | Conception optimale de structures[END_REF] shows that of shapes derivatives applied to a displacement field θ depend only on the normal component of θ on the boundary Γ 0 of the reference domain.

Proposition 1. Let Ω be an open bounded domain of R 2 . Let J : O ad → R be a shape differentiable function. If (θ 1 , θ 2 ) ∈ (W 1,∞ (R 2 , R 2 )) 2 are such that θ 1 -θ 2 ∈ C 1 (R 2 , R 2 ) and θ 1 • n = θ 2 • n, on Γ 0 ∂Ω 0 with θ • n = θ n , then J Ω (θ 1 ) = J Ω (θ 2 ).
(II.11)

Remark 4. As a consequence of the previous lemma, we will replace J (Ω 0 )(θ) by J (Ω 0 )(θ n n) in the following calculations, where θ n = θ • n, and n is a vector field of R 2 that coincides with the exterior normal vector of Ω 0 on Γ 0 .

The shape derivative we defined above as a linear form is a rather abstract object, however, in all the calculations below, it can be associated with a displacement field of the boundary, though the L 2 (Γ) dot product. Definition 3 (shape gradient). Let J : O ad → R be a shape differentiable function. Assume that:

J (Ω 0 )(θ) = Γ θ(x) • w Ω 0 (x) ds (II.12)
then, w Ω 0 is called the shape gradient, and denoted ∇ Ω J Ω 0 .

Useful lemmas

In this section we present the a few lemmas on the calculation of the shape derivative of various functions of the domain, which be useful in the next section.

Integrals on ∂Ω

Let us recall that, in chapter I, the wave-making resistance (which will be our objective function) is written as an integral on the boundary of the obstacle. Hence it is useful to get the integral of functionals of the type:

U (Ω) = ∂Ω f ds . (II.13)
The following proposition, can be seen in [START_REF] Allaire | Conception optimale de structures[END_REF] for instance, in a general form. We give a simpler form and a simple proof in 2D involving a parametrisation of the boundary:

Lemma 2 (change of variables in a boundary integral). Let Ω 0 be a smooth domain of R 2 , and

Ω = T (Ω 0 ), with T = Id + θ, and θ ∈ W 1,∞ (R 2 , R 2 ).
We have:

∂Ω f ds = ∂Ω 0 f • T J(θ) ds (II.14) with: J(θ) = |τ + ∇θ T τ |,
where we recall that τ is the (positivly oriented) tangent vector of ∂Ω 0 , and ∇θ = (∂ i θ j )

Let us introduce x : I → R 2 , a natural parametrisation of Ω 0 . Let us define:

x(s) = x(s) + θ(x(s)) .

(II.15)

From the definition of Ω, it is clear that x is a parametrisation of ∂Ω. Hence:

∂Ω f ds = I f (x(s))|x (s)|ds . (II.16)
Substituting x in (II.15) into the equation (II.16) leads to:

Ω f ds = I f (x(s) + θ(x(s))) d ds (x + θ(x))(s) ds .
(II.17)

First we recognise that x(s) + θ(x(s)) = T (x(s)), and then, we obtain by the chain rule that:

d ds (x + θ(x)) = x + ∇θ T (x)x . (II.18)
Finally, since x is a natural parametrisation of ∂Ω 0 , we have that x = τ (x) and |x | = 1.

Gathering (II.17) and (II.18), we get:

∂Ω f ds = I f (T (x(s))) τ (x(s)) + ∇θ T (x(s))τ (x(s)) ds = Ω 0 f • T τ + ∇θ T τ ds . (II.19)
The following lemma will be used a lot in the following sections:

Lemma 3 (derivative of J). Let J be defined as in lemma 2. We have:

J 0 (θ) = Hθ n (II.20)
where H denotes the curvature of ∂Ω 0 . Note that J 0 represents the Fréchet derivative of θ → J(θ) at

θ = 0.
We expand J at first order with respect to θ:

J(θ) = |τ + ∇θ T τ | , (II.21) = 1 + 2τ • ∇θ T τ + o(|∇θ|) , (II.22) = 1 + τ • ∇θ T τ + o(|∇θ|) . (II .23)
Hence: J 0 (θ) = τ • ∇θ T τ = τ • ∇θτ . Using the proposition 1, we have:

J 0 (θ) = τ • ∇(θ n n)τ , (II.24) = τ • (∇θ n n T + θ n ∇n)τ . (II.25)
On one hand we have, n T τ = 0, on the other hand, we use the following identity (only true in 2D):

τ • ∇nτ = n • (tr(∇n) -(∇n) T )n , (II.26)
to obtain (since ∇nn = 0 as we have seen before in section 5 of chapter I):

τ • ∇nτ = tr(∇n) = ∇ • n = H , (II.27)
by definition of the curvature in two dimension of space. Finally gathering (II.25) and (II.27), we get:

J 0 (θ) = θ n H .
(II.28)

(Id + θ)Ω Ω x (Id + θ)Γ Γ x x + θ(x)

Derivative of a function that depends on the domain

As mentioned before, shape optimisation problems often involve functions which depend on the domain or its boundary (for instance, the boundary potential α for the wave-making resistance). First, we explain how we can define the shape derivative of a function U Ω : Ω → R by comparing U Ω and U (Id+θ)Ω pointwise: Definition 4 (Eulerian derivative). Let x ∈ Ω, the Eulerian derivative of U Ω is defined at x if J(Ω) = U Ω (x) is shape differentiable according to the definition (2). It is hence defined by the linear and continuous mapping on W 1,∞ , which maps θ to ∂ E Ω U Ω satisfying: Let us now consider a function V Γ : Γ → R, where Γ is the boundary of a bounded open set Ω. In this case Γ is closed, hence the Eulerian derivative V Γ has no a-priori meaning (see figure II.2, left). If x ∈ Γ, no matter how small the (non-zero) deformation θ we consider, it is impossible to ensure x ∈ (Id + θ)Γ for all directions of deformation. One solution is to make the point of observation "follow" the domain as it is deformed. Let us define the pulled back on Γ version of V Γ :

U (Id+θ)Ω (x) = U Ω (x) + ∂ E Ω U Ω (θ)(x) + o( θ
Vθ (x) = V (Id+θ)Γ (x + θ(x)) (II.30)
Comparing V Γ and Vθ for θ small leads to the following definition:

Definition 5 (Lagrangian derivative). Let x ∈ Γ, we say that the Lagrangian derivative of V Γ is well defined at x when θ → Vθ (x) is Fréchet differentiable at θ = 0. Hence the Lagrangian derivative is defined by the linear and continuous mapping:

θ → ∂ L Γ V Γ (θ)(x) satisfying: Vθ (x) = V Γ (x) + ∂ L Γ V Γ (θ)(x) + o( θ ).
(II.31) Remark 6. The Lagrangian derivative of V Γ (θ) is the derivative of V (Id+θ)Γ at point (Id + θ)x pulled back to Γ at θ = 0. It writes:

∂ L V Γ (θ)(x) = ∂ θ Vθ | θ=0 (θ)(x).
(II.32)

One example of such a function defined only on the boundary of the domain is the normal vector. The next section describes the details of the calculation of the Lagrangian derivative of both normal vector n Γ and tangent vector τ Γ .

The normal and tangent vector

In this subsection, we detail the calculation of the Lagrangian derivative of both the normal and tangent vector associated to the boundary of a domain. Let us compute the Lagrangian derivative of tangent vector. Let Γ be the reference domain, the deformed domain is Γ = (Id + θ)Γ. As before, we consider x : I → R 2 , a natural parametrisation of Γ. Again, the tangent vector at the point x(s) can be written as:

τ Γ (x(s)) = x (s) x (s) = x (s). (II.33)
As we have seen before, s → x(x) + θ(x(s)) defines a parametrisation of Γ . The pulled-back on Γ tangent vector of Γ writes:

τθ (x(s)) = (x + θ(x)) (s) (x + θ(x)) (s) . (II.34)
Our goal is now to expand (II.34) at first order. In subsection 2.1, we already obtained:

(x + θ(x)) (s) = τ Γ (x(s)) + ∇θ T (x(s)) τ Γ (x(s)) , (II.35) and (x + θ(x)) (s) = 1 + τ Γ • ∇θ T (x(x))τ Γ (x(s)) + o( θ ). (II.36) Further, 1 (x + θ(x)) (s) = 1 -τ Γ • ∇θ T (x(s))τ Γ (x(s)) + o( θ ).
(II.37)

Summing up equations (II.35) and (II.37), we obtain:

τθ = τ Γ + ∇θ T τ Γ -τ Γ • ∇θ T τ Γ τ Γ + o( θ ) (II.38) = τ Γ + ∇θ T τ Γ • n Γ n Γ + o( θ ). (II.39)
From the equation (II.39), we recognise that the shape derivative of the tangent vector is:

∂ L Γ τ Γ (θ) = ∇θ T τ Γ • n Γ n Γ (II.40)
From now on, and for the sake of simplicity, we denote n = n Γ and τ = τ Γ , since there is no ambiguity possible anymore. We hence rewrite equation (II.40):

∂ L Γ τ (θ) = ∇θ T τ • n n . (II.41)
We now use the fact (the proposition (1)) that any shape differentiable quantity depends, at first order, only on the normal component of the deformation field:

∂ L Γ τ (θ) = ∂ L Γ τ (θ n n). (II.42)
Hence we obtain:

∂ L Γ τ (θ) = ∇(θ n n) T τ • n n (II.43) = n ∇θ T n + θ n ∇n T τ • n n (II.44) = n ∇θ T n τ • n n + θ n ∇n T τ • n n (II.45) Since t ∇θ n τ = ∇θ n • τ = ∂ τ θ n and ∇n T τ • n = τ • ∇nn = 0, we have: ∂ L Γ τ (θ) = ∂ τ θ n n. (II.46)
The normal vector is obtained by rotation -π 2 of the tangent, hence we obtain the Lagangian derivative of normal vector:

∂ L Γ n(θ) =R(-π 2 ) ∂ L Γ τ (θ) (II.47) =∂ τ θ n R(-π 2 )n (II.48) = -∂ τ θ n τ. (II.49)

The boundary integral operator

In this subsection we present the calculation of the shape derivative of the boundary integral operator that appears in (I.79). Let us recall the double layer operator in section (5), for x ∈ Γ:

D 1 [α Γ ](x) = γ + Γ [D 1 [α Γ ]] (x) = γ + Γ z → Γ ∇ z G(z, y) • n Γ (z)α Γ (y) ds y (x). (II.50)
Where the boundary potential α Γ is supposed to solve (I.79) (it hence depends intrinsically on Γ). Hence D 1 [α Γ ] depends on Γ in several ways: first it is the trace of an integral on Γ, involving n Γ , then it also depends on α Γ . In this section our focus will be on the first dependence. The main idea is to see D 1 [α Γ ] as a function that depends on Γ and to get the shape derivative of this function. Since D 1 [α Γ ] is defined only Γ, we will try to obtain its Lagrangian derivative. The first step is hence to introduce a deformation T = Id + θ of Γ, and to pull back (II.50) written for T (Γ) onto Γ:

D 1 [α T (Γ) ] • T = γ + T (Γ) z → T (Γ) ∇ z G(z, y) • n T (Γ) (z) α T (Γ) (y) ds y • T. (II.51) Using the fact that γ + T (Γ) (f ) • T = γ + Γ (f • T ) for f continuous in Ω +
, such that f can be continuously extended in Ω + as it is the case for D 1 [α Γ ], we obtain:

D 1 [α T (Γ) ] • T = γ + Γ z → T (Γ) ∇ z G(T (z), y) • n T (Γ) (T (z)) α T (Γ) (y) ds y (II.52)
Let us write it in the term of Γ and θ, using the change of variable formula (II.14), we get:

D 1 [α Γ ] • T =γ + Γ z → Γ ∇ z G(T (z), T (y)) • n T (Γ) (T (z)) α T (Γ) (T (y)) J(θ)(y) ds y (x) (II.53) =γ + Γ z → Γ ∇ z G(z + θ(z), y + θ(y)) • ňθ (z) αθ (y) J(θ)(y) ds y (x) , (II.54)
by using the pulled-back function notation introduced in the subsection 2.2. In the following, we consider (II.54) this as a pulled-back operator Ď1, θ applied to the pulled-back boundary potential αθ . Hence, given any continuous function v defined on Γ, Ď1, θ [v] reads:

Ď1, θ [v] = γ + Γ z → Γ ∇ z G(z + θ(z), y + θ(y)) • ňθ (z) v(y) J(θ)(y) ds y , (II.55)
Here v does not depend on θ, but as stated before, in this section, our efforts are focused on the shape-dependence of D 1 , and not on its argument α Γ . We will see in the next section that the calculation of the shape derivative of α Γ is of few importance from a practical point of view.

Let us now calculate the Fréchet derivative of Ď1, θ [v] with respect to θ at 0. To this end, we define Ď1, θ [v]:

Ď1, θ [v](x) = Γ ∇ x G(x + θ(x), y + θ(y)) • ňθ (x)v(y)J(θ)(y) ds y .
(II.56)

We then have:

Ď1, θ [v] = γ + Γ Ď1, θ [v] .
Let us now expand (II.56) at first order with respect to θ. We have (recalling that we denoted n = n Γ ):

• From Taylor's formula:

∇ x G(x + θ(x), y + θ(y)) = (∇ x G)(x, y) + ∇ 2 xx G(x, y)n(x) θ n (x) + ∇ 2 xy G(x, y)n(y) θ n (y) + o( θ ). (II.57)
• Form the subsection 2.3:

ňθ (x) = n(x) -τ (x) ∂ τ θ n (x)τ (x) + o( θ ) . (II.58)
• From the subsection 2.1:

J(θ) = 1 + H(y) θ n (y) + o( θ ) . (II.59)
Combining (II.57), (II.58) and (II.59), we get:

Ď1, θ [v](x) = Ď1,0 [v](x) + θ n (x) Γ ∇ 2 xx G(x, y)n Γ (x) • n(x) v(y) ds y + Γ ∇ 2 xy G(x, y)n Γ (y) • n Γ (x) θ n (y) v(y) ds y -∂ τ θ n (x) Γ ∇ x G(x, y) • τ Γ (x) v(y) ds y + Γ ∂ n(x) G(x, y) H(y) θ n (y) v(y) ds y + o( θ ). (II.60)
Recalling the definition of S, D 1 , H 1 , and H 2 in equations (I.133), (I.134), and (I.135), the equation (II.60) can be written as:

Ď1,θ [v] = Ď1,0 [v] + θ n H 1 [v](x) + H 2 [θ n v] -∂ τ θ n ∂ τ S[v] + D 1 [Hθ n v] + o( θ ). (II.61)
Further, let us take the trace of Ď1, θ [v](x) using the result in chapter I, section 5, equations (I.136), (I.140), (I.166), and (I.165):

Ď1,θ [v] = Ď1,0 [v] -θ n ∂ 2 τ,τ S[v] -θ n H D 1 [v] + ∂ τ S[∂ τ (θ n v)] -∂ τ θ n ∂ τ S[v] + D 1 [Hθ n v] + o( θ ). (II.62)
Hence:

∂ θ ( Ď1,θ [v])| θ=0 (θ) = -θ n ∂ 2 τ,τ S[v] -θ n H D 1 [v] + ∂ τ S[∂ τ (θ n v)] -∂ τ θ n ∂ τ S[v] + D 1 [Hθ n v] (II.63)

The Lagrangian method

In this section we detail a method to obtain the shape derivative of a quantity involving a state equation, in which the shape derivative of the state equation is not required. In particular, we adopt it to our case which involves the boundary integral equations.

The "classical" case

We consider a shape optimisation problem consisting in finding an open bounded domain Ω minimising some function U (Ω) = G(Ω, u Ω ), where u Ω is a solution in a space V Ω of functions defined on Ω to some variational problem:

a Ω (u Ω , v) = b Ω (v) , for all v ∈ V Ω .
(II.64)

II.3 The Lagrangian method

The Lagrangian method consists in defining the so-called "Lagrangian functional":

L(Ω, u, v) = G(Ω, u) + a Ω (u, v) -b Ω (v) (II.65)
Clearly, when u is replaced by u Ω in the above expression, the second term vanishes:

L(Ω, u Ω , v) = G(Ω, u Ω ) = U (Ω) (II.66)
Hence, applying the chain rule :

U (Ω) = ∂ Ω L(Ω, u Ω , v) + ∂ u L(Ω, u Ω , v) , ∂ E Ω u Ω (II.67)
where ∂ E Ω u Ω denotes the Eulerian shape derivative of u Ω at Ω, which is defined as the only linear mapping on the set of W 1,∞ deformation fields satisfying, for all x ∈ Ω:

u (Id+θ)Ω (x) = u Ω (x) + ∂ E Ω u Ω (θ)(x) + o(|θ| W 1,∞ ) . (II.68) Remark 1.
Note that the consistence of the above expression is not obvious since it consists in comparing at the same point x ∈ Ω the function u Ω with a function u (Id+θ)Ω which is defined on a different set. We could argue that for θ given, this comparison cannot be performed, at least for some points of Ω. Fortunately Ω is an open set and so, it is always possible to find a small-enough deformation field θ such that x belongs to both Ω and (Id + θ)Ω.

Remark 2. The existence of ∂ E Ω u Ω is assumed here, but it is a fact that should be proved beforehand.

The calculation ∂ E Ω u Ω is tedious, and if we could find v Ω such that ∂ u L(Ω, u Ω , v Ω ) = 0, it would not be necessary, since, then:

U (Ω) = ∂ Ω L(Ω, u Ω , v Ω ) (II.69)
The calculation of v Ω satisfying the above condition involves solving a so-called "adjoint equation" which is obtained by nulling out ∂ u L.

The case of boundary integral equations

Consider now that our objective function is in the form U (Γ) = G(Γ, u Γ ) where u Γ , defined on V Γ a set of functions defined on Γ, solves a variational problem of the form:

a Γ (u Γ , v) = b Γ (v) for all v ∈ V Γ (II.70)
Like before, we could define a Lagrangian:

L(Γ, u, v) = G(Γ, u) + a Γ (u, v) -b Γ (v) (II.71)
and the following would still be true for any v ∈ V Γ :

U (Γ) = L(Γ, u Γ , v) (II.72)
However, it is in this case not possible to formally apply the chain rule as before to obtain U (Γ), since the Eulerian derivative of u Γ is not defined at all (because Γ is a closed set).

As defined, one workaround is to pull everything back onto the domain Γ around which the shape derivative is being evaluated. Let us recall:

ǔθ = u (Id+θ)Γ • (Id + θ) (II.73)
then, ǔθ is defined on Γ and it solves the following variational problem:

ǎθ (ǔ θ , v) = bθ (v) for all v ∈ V Γ (II.74)
where :

ǎθ (u, v) = a (Id+θ)Γ (u • (Id + θ) -1 , v) and bθ = b (Id+θ)Γ . Also, if we denote Ǧ(θ, u) = G((Id + θ)Γ, u • (Id + θ) -1
), we have that :

U ((Id + θ)Γ) = Ǧ(θ, ǔθ ) (II.75)
Hence, defining the Lagrangian:

Ľ(θ, u, v) = Ǧ(θ, u) + ǎθ (u, v) -bθ (v) (II.76)
we clearly have:

U ((Id + θ)Γ) = Ǧ(θ, ǔθ ) = Ľ(θ, ǔθ , v) .
(II.77)

II.4 Application to the wave-making resistance problem

Since U (Γ) is in fact the derivative of U ((Id + θ)Γ) with respect to θ at θ = 0, we have:

U (Γ) = ∂ θ Ľ(0, ǔ0 , v) + ∂ u Ľ(0, ǔ0 , v) , ∂ θ ǔ0 . (II.78)
Remark 3. In literature (see for instance [START_REF] Allaire | Conception optimale de structures[END_REF]), ∂ θ ǔ0 is in fact oftenly referred as the "Lagrangian derivative" or "material derivative" of u Γ with respect to the shape of Γ.

Hence, if we can find the function v 0 that nulls out ∂ u Ľ(0, ǔ0 , v), that is, v is a solution to adjoint equation, then :

U (Γ) = ∂ θ Ľ(0, ǔ0 , v 0 ) (II.79)
where ǔ0 solves:

ǎ0 (ǔ 0 , v) = b0 (v) for all v ∈ V Γ (II.80)
which is in fact the same as :

a Γ (ǔ 0 , v) = b Γ (v) for all v ∈ V Γ (II.81)

Application to the wave-making resistance problem

Wrapping up (I.132) and (I.79), our shape optimisation problem reads : find the domain Γ that minimises:

R w (Γ) = ρ ν 4 Γ α Γ (x) E(x)ds x 2 , (II.82)
where α Γ solves:

1 2 α Γ (x) - Γ ∂ n Γ (x) G(x, y) α Γ (y) ds y = U ∞ n Γ (x) • e 1 , for all x ∈ Γ , (II.83)
or equivalently:

D 1 [α Γ ] = -U ∞ n Γ • e 1 , on Γ . (II.84)
To make notations easier, for the rest of the section, we drop the constant ρ ν/4. In the next section we describe the shape optimisation method that we used, based on the calculation of the shape gradient. Following the approach mentioned above, we first rewrite (II.82) and (II.84) in terms of θ. First (II.82) becomes:

R w ((Id + θ)Γ) = (Id+θ)Γ α (Id+θ)Γ (x)E(x)ds x 2 .
(II.85)

Using the change of variable formula for integral over curves (or surfaces) (lemma 2), we have:

R w ((Id + θ)Γ) = Γ αθ (x)E(x + θ(x)) J θ (x) ds x 2 .
(II.86)

Similarly, pulling back both sides of (II.83), we have:

Ď1,θ [ αθ ] = -U ∞ ňθ • e 1 , (II.87)
Let us define the Lagrangian:

Ľ(θ, α, q) = Γ α(x)E(x + θ(x)) J θ (x) ds x 2 + Γ q(x) Ď1 θ [α θ ](x) + U ∞ ňθ (x) • e 1 ds x (II.88)
As explained in subsection 3.2, the derivation of this Lagrangian with respect to α will provide with the adjoint equation.

Determination of the adjoint equation

In this section, let us calculate the adjoint equation by nulling out ∂ α Ľ(0, α0 , q). Remarking that Ď1,0 = D 1 , we write (II.88) for θ = 0:

Ľ(0, α, q) = Γ α(x)E(x) ds x 2 + Γ q(x) [D 1 [α](x) + U ∞ n(x) • e 1 ] ds x (II.89)
Let us compute the derivative of L with respect to α by expanding it with respect to α:

Ľ(0, α + α, q) = Γ (α + α)(x)E(x) ds x 2 + Γ q(x) [D 1 [α + α](x) + U ∞ n(x) • e 1 ] ds x (II.90) = Ľ(0, α, q) + 2Re Γ αE ds Γ αE ds + Γ qD 1 [ α] ds + o( α ) (II.91)
Further, we recognize:

∂ α Ľ(0, α, q) , α = 2 Γ α(y) Γ Re E(x)E(y) α(x) ds x ds y + Γ q(x)D 1 [ α](x)ds x . (II.92)
By identifying a L 2 scalar product in second term and denoting D * 1 the adjoint of D 1 , we get:

∂ α Ľ(0, α, q) , α = Γ α(y)   2 Γ Re E(x)E(y) α(x) ds x + D * 1 [q](y)   ds y , (II.93)
where:

D * 1 [q](y) = Γ ∂ n(x) G(x, y) q(x) ds x - 1 2 q(y) .
(II.94)

Considering that q 0 ensure ∂ α Ľ(0, α0 , q 0 ) , α = 0 for ∀α ∈ C 0 (Γ), means that, knowing α0 from equation (II.83):

D 1 [ α0 ] = -U ∞ n(x) • e 1 , (II.95)
q 0 solves the adjoint equation:

D * 1 [q 0 ](y) = -2 Γ Re E(x)E(y) α0 (x) ds x .
(II.96)

Calculation of the derivative of the Lagrangian with respect to θ

We calculate ∂ θ Ľ(0, u, v) by separating (II.88) into 3 terms:

Ľ1 (θ, α, q) = Γ α(x) E(x + θ(x)) J θ (x) ds x 2 , (II.97) Ľ2 (θ, α, q) = Γ q(x) Ď1,θ [α](x) ds x , (II.98) Ľ3 (θ, α, q) = Γ q(x) U ∞ ňθ (x) • e 1 ds x .
(II.99)

Noting that :

E(x + θ(x)) = E(x) + θ n (x) ∂ n E(x) + o(|θ|) ,
(II.100) and:

J θ (x) = 1 + θ n (x) H(x) + o(|θ|) , (II.101)
where H represents the curvature of Γ, we have:

Ľ1 (θ, α, q) = Γ α(x) (E(x) + θ n (x) ∂ n E(x)) (1 + θ n (x) H(x)) ds x 2 + o(|θ|) , = Γ α(x)E(x)ds x + Γ θ n (x) α(x) (H(x)E(x) + ∂ n E(x)) ds x 2 + o(|θ|) , = Ľ1 (0, α, q) + 2 Γ θ n (x) α(x) Γ α(y) Re H(x) E(x) + ∂ n E(x) E(y) ds y ds x + o(|θ|) .
Hence, the derivative with respect to θ of Ľ1 writes:

∂ θ Ľ1 (0, α, q)(θ) = 2 Γ θ n (x) α(x) Γ α(y) Re H(x) E(x) + ∂ n E(x) E(y) ds y ds x . (II.102)
Let us now move on the second term. From the derivative of Ď1, θ with respect to θ we have determined in subsection 2.4, we have:

∂ θ Ľ2 (0, α, q)(θ) = Γ q(x) (D 1 [Hθ n α](x) -θ n H(x)D 1 [α](x)) ds x - Γ q(x)∂ τ θ n (x)∂ τ S[α](x) -S[∂ τ (θ n α)](x) ds x . (II.103)
Integrating by parts the second term, we obtain:

∂ θ Ľ2 (0, α, q)(θ) = Γ q(x) (D 1 [Hθ n α](x) -θ n H(x)D 1 [α](x)) ds x + Γ ∂ τ q(x) θ n (x)∂ τ S[α](x) -S[∂ τ (θ n α)](x) ds x . (II.104)
Finally, we compute the expansion of the Ľ3 (θ, α, q)(θ) by using the shape derivative of the normal vector obtained in subsection 2.3, we obtain:

Ľ3 (0, α, q)(θ) = Ľ3 (0, α, q)(θ) -U ∞ Γ ∂ τ θ n (x) q(x) τ (x) • e 1 ds x + o( θ ).
(II.105)

Hence we have:

∂ θ Ľ3 (0, α, q)(θ) = -U ∞ Γ ∂ τ θ n (x) q(x) τ (x) • e 1 ds x , (II.106) = U ∞ Γ θ n (x) ∂ τ (q(x) τ (x) • e 1 ) ds x , (II.107) = U ∞ Γ θ n (x) ∂ τ q(x) τ (x) • e 1 ds x + U ∞ Γ θ n (x) q(x) ∂ τ (τ (x) • e 1 ) ds x .
(II.108) By using the fact that ∂ τ τ = -Hn as in equation (I.158), we obtain:

∂ θ Ľ3 (0, α, q)(θ) = U ∞ Γ θ n (x) ∂ τ q(x) τ (x) • e 1 ds x -U ∞ Γ θ n (x) q(x) H(x) n(x) • e 1 ds x .
(II.109) Further, summing (II.102), (II.104) and (II.109), we obtain the derivative of Lagrangian Ľ:

∂ θ Ľ(0, α, q)(θ) =2 Γ θ n (x) H(x) α(x) Γ α(y) Re E(x)E(y) ds y ds x (II.110) +2 Γ θ n (x) α(x) Γ α(y) Re ∂ n E(x)E(y) ds y ds x (II.111) + Γ q(x) (D 1 [Hθ n α](x) -θ n H(x)D 1 [α](x)) ds x (II.112) + Γ ∂ τ q(x) θ n (x)∂ τ S[α](x) -S[∂ τ (θ n α)](x) ds x (II.113) +U ∞ Γ θ n (x) ∂ τ q(x) τ (x) • e 1 ds x -U ∞ Γ θ n (x) q(x) H(x) n(x) • e 1 ds x .
(II.114)

We can simplify the above expression by noticing that :

Γ q(x)D 1 [α H θ n ](x) ds x = Γ α(x) H(x) θ n (x) D * 1 [q](x)ds x (II.115)
which leads to:

∂ θ Ľ(0, α, q)(θ) = Γ θ n (x) H(x) α(x) D * 1 [q](x) + 2 Γ α(y) Re E(x)E(y) ds y ds x (II.116) - Γ q(x) H(x) θ n (x) D 1 [α](x) + U ∞ n(x) • e 1 ds x (II.117) +2 Γ θ n (x) α(x)
Γ α(y) Re ∂ n E(x)E(y) ds y ds x (II.118)

+ Γ ∂ τ q(x) θ n (x)∂ τ S[α](x) -S[∂ τ (αθ n )](x) ds x (II.119) +U ∞ Γ θ n (x) ∂ τ q(x) τ (x) • e 1 ds x .
(II.120)

When α and q solve respectively the equation of state (II.95) and the adjoint equation (II.96), it is not difficult to see that the two first terms are canceled (we get rid of them already for the sake of simplicity). Moreover, we notice that

Γ ∂ τ q(x) S[∂ τ (αθ n )](x) ds x = ∂ τ q, S[∂ τ (αθ n )] L 2 (Γ) = S * [∂ τ q], ∂ τ (αθ n ) L 2 (Γ) = Γ ∂ τ (y) (α(y)θ n (y)) S * [∂ τ q](x) ds y ,
where S * [u](x) = Γ G * (x, y)u(x)ds x . By integrating by parts, we obtain:

Γ ∂ τ q(x) S[∂ τ (αθ n )](x) ds x = - Γ θ n (y)α(y)∂ τ S * [∂ τ q](x)ds y (II.121)
Hence, we have the derivative of Lagrangian with respect to θ as follows:

∂ θ Ľ(0, α, q)(θ) =2 Γ θ n (x) α(x) Γ α(y) Re ∂ n E(x)E(y) ds y ds x (II.122) + Γ θ n (x) ∂ τ q(x) ∂ τ S(α)(x)ds x + Γ θ n (x)α(x)∂ τ S * [∂ τ q](x)ds x (II.123) +U ∞ Γ θ n (x) ∂ τ q(x) τ (x) • e 1 ds x . (II.124)
From the derivative of Lagrangian with respect to θ in equation (II.124), we write the shape derivative of the wave-making resistance as:

R w (Γ)(θ) = Γ θ n (x) 2α 0 (x) Γ α0 (y) Re ∂ n E(x)E(y) ds y + ∂ τ q 0 (x) ∂ τ S(α 0 )(x) + α0 (x)∂ τ S * (∂ τ q 0 )(x) + U ∞ ∂ τ q 0 (x) τ (x) • e 1 ds x , (II.125)
where, we recall, α0 and q 0 solve:

1 2 α0 (x) - Γ ∇ x G(x, y) • n(x) α0 (y) ds y = U ∞ n Γ (x)
• e 1 , for all x ∈ Γ, (II.126)

1 2 q 0 (y) - Γ ∇ x G(x, y) • n(x) q 0 (x) ds x = 2 Γ
Re E(x)E(y) α0 (x) ds x , for all y ∈ Γ (II.127) Moreover, identifying a L 2 (Γ) scalar product in the equation (II.125), we obtain the following expression for the shape gradient:

∇ Γ R w (x) = n(x) α0 (x) Γ α0 (y) Re ∂ n E(x)E(y) ds y + ∂ τ q 0 (x) ∂ τ S(α 0 )(x) + α0 (x)∂ τ S * (∂ τ q 0 )(x) + U ∞ ∂ τ q 0 (x) τ (x) • e 1 . (II.128)

Shape gradient descent method and constraints

In this subsection, we describe the principle of the numerical method that we used to obtain a local minimizer of R w . We present also the volume constraint that we imposed to prevent a trivial solution i.e. the obstacle shrink and continue to shrink. We want to avoid the minimum wave resistance is reached whenever there is no obstacle in the fluid.

Here we use the gradient descent method to find numerically a shape which is a local minimizer of R w . The idea of the shape gradient descent method is to transform Γ for each step of the method with a deformation field θ chosen in such a way that R w ((Id + θ)Γ k ) is minimal for a given step size θ . Definition 6. Let the boundary of our shape be Γ such that

∇ Γ R w = 0. A descent direction for R w (Γ) is a vector field θ ∈ W 1,∞ ensuring: R w ((Id + θ)Γ) < R w (Γ) (II.129)
for θ is a sufficiently small in order to ensure Id + θ is a W 1,∞ diffeomorphism.

Further, let us find a descent direction of the method. Let us recall the expansion of R w :

R w ((Id + θ)Γ) = R w (Γ) + ∇ Γ R w (Γ) , θ L 2 + o( θ ) . (II.130)
From there, it is clear that taking θ = -δ r ∇ Γ R w (Γ), for a given δ r > 0, is a descent direction, provided that δ r is small enough. To show that, let us replace θ by the shape gradient in the above expression:

R w ((Id + θ)Γ) = R w (Γ) -δ r ∇ Γ R w (Γ) 2 L 2 + o(δ r ) . (II.131) Since -δ r ∇ Γ R w (Γ) 2 L 2 + o(δ r
) < 0 for δ r small enough, we have that R w ((Id + θ)Γ) < R w (Γ), hence θ is a descent direction. It is moreover clear from (II.130) that ∇ Γ R w (Γ) is the "steepest" descent direction for R w . The algorithm 1 exploits this idea by deforming step by step a given initial domain according to its local shape gradient.

Remark 7. It is not clear at all from (II.128) that ∇ Γ R w (Γ) belongs to W 1,∞ . One could argue that this depends on the regularity of α Γ and q Γ , which in turn depend on the regularity of Γ, which might change over the iterations of the algorithm 1. Moreover, we do not give here any information on how δ r should be chosen, and it is quite possible that the value δ r ensuring the condition that (Id + θ) is a W 1,∞ diffeomorphism might vanish to 0 as the algorithm iterates. These theoretical questions are rather difficult, and beyond the scope of this thesis. From a practical point of view, we will proceed empirically for the determination of δ r (it is hence considered provided here).

As many optimisation problems, the wave-making resistance problem should be endowed with some constraints. Two situations should be avoided: i) a sequence of smaller and smaller obstacles for which R w (Γ n ) → 0, ii) a sequence of deeper and deeper obstacles for which R w (Γ n ) → 0 also.

Clearly, the first situation can be addressed by adding to the problem the constraint µ(Ω -) = v where µ is the measure of the surface area of Ω -and v > 0 is given. The second situation can be avoided by setting the center of gravity of the domain: Ω -x 2 dx = d where d is a constant. However, we will see in the numerical results that it is not necessary to enforce this condition since local minima of finite depth will be found.

In order to enforce the measure constraint, our approach consists in "projecting" the shape of the obstacle after each step on the set of shapes of measure v by applying a normal displacement with a magnitude a which is determined numerically. Let κ(Γ) be a measure of the domain enclosed by Ω -. The algorithm 2 explains how this deformation can be obtained numerically.

Data: A shape Γ, a tolerance , a target surface area v and a step δ.

Result:

The "projected" Γ = (Id + an) Γ, where a is such that µ(Ω -) = v. Both the step of the gradient descent method and the projection involve a normal displacement of the curve Γ which defines the shape of the obstacle. As we will see in the next chapter, from a numerical standpoint, this operation can be difficult when we consider the discretisation of a parametrised curve. Luckily, this operation is straightforward if we choose to use the level-set method.

Chapter III

The Level-set Method

In this chapter, we give a presentation of the level-set method that we will use to implement numerically the algorithm presented at the end of chapter II. Our presentation is focused on the main ingredients that we need : deforming the boundary, and solve boundary integral equations.

As much as it is possible, this presentation is self-contained and we will often refer to [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF] (in particular in section 1), which contains all the numerical methods we will use to represent and displace the boundary (and also shows some very useful applications in computational physics and imaging). For the calculation of boundary integrals and the approximation of solutions of boundary integral equations that we present in sections 2 and 3, we follow the work of C. Kublik and R. Tsai, in particular [START_REF] Kublik | An implicit interface boundary integral method for Poisson's equation on arbitrary domains[END_REF] [49] [START_REF] Engquist | Discretization of Dirac delta functions in level set methods[END_REF].

1 General introduction to level-set methods

Explicit and implicit representation of curves

In general, there are two ways to represent curves and surfaces: explicitly or implicitly. The explicit representation considers the curve Γ as given by the image of a mapping γ :

I ⊂ R → R 2 called the parametrisation of Γ: Γ = {γ(s) | s ∈ I} .
(III.1)

As an example we have, for instance C(0, 1) = {(cos(s), sin(s)) | s ∈ [0, 2π[}. Wherever γ is differentiable, and the derivative is non-zero (we say that γ is non-stationary), the unit tangent vector is given by:

τ (γ(s)) = γ (s) |γ (s)| . (III.2)
The unit normal vector is obtained by rotating τ by -π/2 (considering that γ runs on Γ anticlockwise). This parametrisation is called explicit because the points of Γ are given explicitly by γ. From a numerical standpoint, we build parametric curves from a collection of control points (x i ) i=1..N , and the curve between these points is defined using linear, quadratic or cubic interpolation (see figure III.1). Among all possible parametrisations of Γ, one is called the natural parametrisation and satisfies |γ (s)| = 1. This type of representation of curves was already used in chapter II, for instance in the calculation of the shape derivative of the normal and tangent vectors. Now let us move on the implicit representation of curves of the plane. In the implicit representation, the curve is given by the inverse image of a single value (by convention 0) through a function φ : R 2 → R. The function φ is then called the level-set function, and Γ writes:

Γ = {x ∈ R 2 | φ(x) = 0} . (III.3)
As before, we have, for example which implies (assuming that the gradient of φ, denoted ∇φ, does not take zero as a value anywhere on Γ) that the gradient of the level-set function is orthogonal to the tangent vector on Γ, hence, co-linear to the normal vector. Furthermore, if we impose that φ < 0 in Ω -and φ > 0 in Ω + , we have that the gradient of φ points outwards. We then have:

C(0, 1) = {x ∈ R 2 | |x| -1 = 0}.
n = ∇φ |∇φ| on Γ . (III.5)
This representation is called implicit because in order to obtain the points of Γ we need to solve the equation φ(x) = 0. From a numerical stand point, since Φ is a function of R 2 , it is necessary to give a discretisation of the region of R 2 that contains the curve. This mesh is not necessarily body fitted for Γ, and the usual practice is to use cartesian meshes because of their simplicity. Let us consider the points

(x ij ) i=1..N 1 j=1..N 2
of a cartesian mesh of a square region containing Γ. Let φ ij be the values of the level-set function on these nodes, then the interpolation (bilinear, cubic, ...) of φ between the node points provides a piece wise construction of the curve on the mesh (figure III.1).

At first glance, is seems that the explicit representation of curves is more natural, since, numerically, it involves less data. As an example, let us consider that we want to represent a circle of center 0 and radius 1 with N p points for the explicit representation (2N p data in 2D). If we want to represent the same circle with a similar accuracy with a level-set function, that would require (LN p/π) 2 mesh points (where L is the length of the square enclosing the circle). Therefore, the implicit representation requires O(N 2 p ) data, while O(N p) data is required for the explicit representation. The reason behind this discrepancy is the very large amount of "useless" data considered in the implicit representation. Indeed, the only values of the level-set function that matter are the ones for which a neighbour node has a different sign. Considering only those points, we recover O(N p ) data for the implicit representation (see figure III.2).

III.1 General introduction to level-set methods

The main advantage of the implicit method lies in the numerical treatment of moving surfaces that will be described later in section 1.3. Before that, we present a particular kind of level-set functions which are analogous to the "natural" parametrisation of curves in the sense that their gradient is of norm 1.

The signed distance function

Let us now focus on the implicit representation of curves. Let us consider our curve as a boundary Γ of an open bounded domain Ω -, as before let us denote Ω + the exterior domain. We define the signed distance function d as:

d(x) =            0, for x ∈ Γ , -min x 0 ∈Γ |x -x 0 | for x ∈ Ω -, min x 0 ∈Γ |x -x 0 | for x ∈ Ω + .
(III.6)

The function d gives the closest distance between x and the boundary, with a positive sign in the exterior domain, and a negative sign in the interior domain. Let us give two examples of such a function:

• For the circle of center 0 and radius one (figure III.3, left):

d(x) = x 2 1 + x 2 2 -1 . (III.7)
• For the square of center 0 and side two see (figure III.3, right):

d(x) =      max(|x 1 |, |x 2 |) -1 , for |x 1 | ≤ 1 or |x 2 | ≤ 1 , (|x 1 | -1) 2 + (|x 2 | -1) 2 , elsewhere. (III.8)
Le us show that, in general, the signed distance function defined in (III.6) satisifies an Hamilton-Jacobi equation called the eikonal equation. Let us consider x ∈ Ω + such that the projection of x onto Γ, denoted Π Γ (x), exists and is unique 1 . Hence:

x Γ d(x) u Π Γ (x)
d(x) = |x -Π Γ (x)| (III.9)
First let us prove that ∇d is colinear to:

u = Π Γ (x) -x |Π Γ (x) -x| . (III.10)
In order to find the direction of the gradient of d we seek the direction of steepest descent of d. Let us solve the minimisation problem:

min |θ|= d(x + θ) . (III.11)
First, we prove that d(x + θ) ≥ |x -Π Γ (x)| -. We have:

d 2 (x + θ) = (x + θ -Π Γ (x + θ)) • (x + θ -Π Γ (x + θ)) , (III.12) = |x -Π Γ (x + θ)| 2 + 2 + 2θ • (x -Π Γ (x + θ)) . (III.13)
From Cauchy-Schwartz inequality we have that:

θ • (x -Π Γ (x + θ)) ≥ -|θ| |x -Π Γ (x + θ)| = -|x -Π Γ (x + θ)| , (III.14)

III.1 General introduction to level-set methods

which, combined with (III.13) leads to:

d 2 (x + θ) ≥ |x -Π Γ (x + θ)| 2 + 2 -2 |x -Π Γ (x + θ)| , (III.15) ≥ (|x -Π Γ (x + θ)| -) 2 (III.16)
If is small enough then we have:

d(x + θ) ≥ |x -Π Γ (x + θ)| -. (III.17)
By the definition of the projection we have that:

|x -Π Γ (x + θ)| ≥ |x -Π Γ (x)|
, hence, the following estimate holds for d(x + θ):

d(x + θ) ≥ |x -Π Γ (x)| -. (III.18)
Now, we give a particular value of θ for which the equality holds. A simple geometrical argument (see figure III.4) shows that,

d(x + u) = |x -Π Γ (x)| -. (III.19)
and so that the minimum of d(x + θ) is attained for θ = u. We deduce from this that the steepest descent direction is u, hence, the gradient of d at x is colinear to u. Now, let us calculate |∇d| by evaluating a variation in the direction u:

∇d(x) • u = lim →0 d(x + u) -d(x) (III.20)
From figure III.4 we can see that: Π Γ (x + u ) = Π Γ (x), hence:

d(x + u) = |x + u -Π Γ (x)| = |x -Π Γ (x)| -= d(x) -, (III.21)
from the co-linearity of x -Π Γ (x) and u. Then we have:

d(x + u) -d(x) = d(x) --d(x) = -1 . (III.22)
Then the gradient of d in the direction u writes:

∇d(x) • u = -1 . (III.23)
Recalling the ∇d(x) is co-linear to u, which is of norm 1. We finally have:

|∇d(x)| = 1 , (III.24)
for x ∈ Ω + . Doing the same in Ω -, in order to get d from Γ we need to solve the system:

     |∇d(x)| = 1 , for x ∈ Ω + ∪ Ω -, d(x) = 0 , for x ∈ Γ , (III.25)
also known as the eikonal equation. From (III.25) and III.5, the normal vector at x ∈ Γ can be written:

n(x) = ∇d(x) |∇d(x)| = ∇d(x). (III.26)
Furthermore, we define the curvature of the boundary as the divergence of the normal vector. Hence, we have:

H(x) = ∇ • n(x) = ∆d(x). (III.27)
Even though III.25 provides an equation for the signed distance function, it is difficult to solve or even prove that solutions exists. In particular, one cannot expect solutions in the classical C 1 sense, since the signed distance function is not continuously differentiable in general. As an example we can cite the case of the square defined in (III.8) for which:

lim h→0 ± ∇d(h, 0) = ±e 1 , lim h→0 ± ∇d(0, h) = ±e 2 . (III.28)
Luckily, these "accidents" of non-derivability only happen on regions of zero measure. Knowing that, we could search for generalised solutions that are only Lipschitz continuous and solve the eikonal equation almost everywhere (i.e. everywhere but for a set of zero Lebesgue measure). Looking for such generalised solutions, we obtain an infinite amount of solutions, which is not a satisfactory situation either (see for instance [66]). It is however possible to select among all these solutions the one that actually defines the signed distance function. This solution is called the viscosity solution. It is obtained by solving the modified equation:

     -∆d + |∇d (x)| -1 = 0, for x ∈ Ω + ∪ Ω -, d (x) = 0 for x ∈ Γ , (III.29)
which is solvable in the classical sense. Then, passing to the limit → 0, we obtain the viscosity solution. This limit is very difficult to obtain in the general case. A more practical approach was introduced by S.N. Kruzkov for the eikonal equation in which the "generalised solution" is only required to have all its discrete second derivatives (in all directions, for all steps) bounded from below (in addition to all the previous requirements). Such a solution has been proven to be the same as the aforementioned viscosity solution (see [START_REF] Kruzkov | Generalized solutions of the Hamilton-Jacobi equations of eikonal type. I. formulation of the problems; existence, uniqueness and stability theorems; some properties of the solutions[END_REF]). This notion of viscosity solutions has then been extended to the general case of first order Hamilton-Jacobi equations by M.C. Crandall and P.-L. Lions. in their famous article of 1981 [18]. In a later section we will show how to obtain an numerical approximation of the viscosity solution of (III.25).

Normal motion

We recall that one of the most important ingredients in the algorithms 1 and 2 described at the end of chapter II is the displacement in the normal direction of the boundary of the obstacle. Let Γ k be the initial obstacle, and Γ k+1 the deformed obstacle though the deformation field θ n n, we have:

Γ k+1 = (Id + θ n n)Γ k . (III.30)
Using a parametrisation γ k for Γ k , we obtain a new parametrisation for Γ k+1 with:

γ k+1 = γ k + θ n (γ k ) R(-π/2) γ k |γ k | , (III.31)
Now, using a level-set function φ k for Γ k , we get a new level-set function for Γ k+1 by solving for φ k+1 the equation:

φ k+1 • (Id + θ n n) = φ k . (III.32)
where, we recall n = ∇φ/|∇φ|. In principle, θ n is given only on Γ k , but it is possible to extend it in Ω + ∪ Ω -. As we will see below, when φ k is a signed distance function and θ n is extended in a constant-along-normals manner φ k+1 can be obtained easily. First, since φ k is a signed distance function, we have (far away from the caustics):

φ k (x + θ n (x)n(x)) -φ k (x) = θ n (x) , (III.33)
which can be rewritten:

φ k (x + θ n (x)n(x)) -θ n (x) = φ k (x) . (III.34)
Since θ n is constant along normals, we have θ n (x) = θ n (x + n(x)θ n (x)). Hence:

(φ k -θ n )(x + θ n (x)n(x)) = φ k (x) . (III.35)
Comparing (III.32) and (III.35), we get:

φ k+1 = φ k -θ n , (III.36)
provided φ k is a signed distance function and θ n is constant along normals. Unfortunately, φ k+1 is not a signed distance function (unless θ n is constant), and it will be necessary to apply a redistancing procedure in order to obtain a signed distance function from φ k+1 . This procedure, along with the construction of a constant-along-normals θ n function, will be detailed in subsection 1.4.

Redistancing and extension along normals

In this subsection we detail how it is possible to obtain a signed distance function d and a constant-along-normals function a from an initial level-set function d 0 and non-constantalong normal a 0 . Recall from subsection 1.2, that d and a have to solve:

|∇d| = 1 , in Ω , (III.37) ∇d • ∇a = 0 , in Ω , (III.38) d = 0 , on Γ , (III.39) a = a 0 , on Γ , (III.40)
where Γ = d -1 0 ({0}), and Ω a bounded open domain enclosing Γ. We consider that Ω = ]0, 1[×]0, 1[ and Γ ⊂ Ω, without loss of generality. There are several methods to achieve this goal (see for instance [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF], for a general review of these methods and [68] for more details on the Fast Marching Method). Our approach follows [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF] in which (III.37)-(III.40) are solved by solving their time dependent counterpart:

∂ t d + sign(d 0 )(|∇d| -1) = 0 , in R + × Ω , (III.41) ∂ t a + sign(d 0 )∇d • ∇a = 0 , in R + × Ω , (III.42) d(0, x) = d 0 (x) , for x ∈ Ω , (III.43) a(0, x) = a 0 (x) , for x ∈ Ω . (III.44)
It is clear that the steady state of (III.41)-(III.44) solves (III.37)-(III.40). The idea behind this particular formulation is that, both (III.41) and (III.42) are equations of the form:

∂ t u + H(•, u, ∇u) = 0 , in R + × Ω , (III.45) u(0, x) = u 0 (x) , for x ∈ Ω . (III.46)
We consider a discretisation of space on a Cartesian grid : (x 1 i,j , x2 i,j ) for i = 1..N 1 and j = 1..N 2 , where x 1 i,j = ih 1 , and x 2 i,j = jh 2 , with h 1 = 1/(N 1 -1) and h 2 = 1/(N 2 -1). Following [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF], we approach (III.45) with:

u n+1 i,j = u n i,j -δt H x i,j , u n i,j , (∂ + 1 u) n i,j , (∂ - 1 u) n i,j , (∂ + 2 u) n i,j , (∂ - 2 u) n i,j , (III.47)
where u n i,j represents the approximation of u(t n , x i,j ), and

(∂ + k u) n i,j (r.p. (∂ - k u) n i,j
) represents the right (r.p. left) approximation of the space derivative of u in the direction k at time t n and point x i,j . These left and right derivative approximations are obtained in our case using a WENO5 2 scheme, as it is presented in [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF] and initially in [69]. The possible choices of a discrete Hamiltonian H is also explained in details in [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF], and here we have chosen the Godunov scheme, which for (III.41) writes:

H(III.41) x, d, d + 1 , d - 1 , d + 2 , d - 2 =                            σ(x) max(-min(d + 1 , 0), max(d - 1 , 0)) 2 + max(-min(d + 2 , 0), max(d - 2 , 0)) 2 -1 if d ≥ 0 , σ(x) max(-min(d - 1 , 0), max(d + 1 , 0)) 2 + max(-min(d - 2 , 0), max(d + 2 , 0)) 2 -1 if d < 0 , (III.48)
where σ(x) is a smooth approximation of sign(d 0 (x)), given by:

σ(x) = d 0 (x) |∇d 0 | 2 d 0 (x) 2 + h 2
for h = min(h 1 , h 2 ). For (III.42) the same Godunov scheme gives the classical upwind-type discretisation:

H(III.42) x, a, a + 1 , a - 1 , a + 2 , a - 2 =    min(σ(x)d - 1 a - 1 , σ(x)d + 1 a + 1 ) if a - 1 ≤ a + 1 max(σ(x)d - 1 a - 1 , σ(x)d + 1 a + 1 ) if a - 1 > a + 1 +    min(σ(x)d - 2 a - 2 , σ(x)d + 2 a + 2 ) if a - 2 ≤ a + 2 max(σ(x)d - 2 a - 2 , σ(x)d + 2 a + 2 ) if a - 2 > a + 2 (III.49)
where, in (III.49), d ± k , for k = 1, 2 are obtained from the (supposedly already computed) signed distance function approximation d i,j .

We illustrate this method with a simple text case where Γ is the zero level-set of:

d 0 (ρ, θ) = 1 -ρ -1 2 cos(2θ) , (III.50)
and the scalar function a 0 we wish to extend in a constant-along normals fashion is:

a 0 (x) = sin(6πx 2 ) . (III.51)
The figure III.5 shows the evolution of the scheme (III.47) for the redistancing procedure and the extension procedure. We remark that the information is propagated from Γ and outwards at a rate of δt/h units of space per iteration. Note that it is required that δt/h < 1 for the stability of the scheme (Courant-Freidrichs-Lewy condition), hence the information propagates of at most h units of space per iteration. Hence, in principle, O( √ N ) iterations are required to achieve convergence (where N is the number of points in the mesh). Hence, the complexity of the algorithm is N 2 . In comparison, the Fast-Marching method has a complexity of N log(N ). However, this is of little importance, since the only values of (d, a) we need for our applications are values in a small neighbourhood of Γ. If this neighbourhood T has a thickness m h for m ∈ N, then, by the above argument, approximately m iterations are required to obtain (d, a) satisfying (III.37)-(III.38) on T.

Normal motion: implicit v.s. explicit

In order to illustrate the above method, we consider the sequence of shapes defined by: where:

Γ k+1 = (Id + θ n n)Γ k , (III .52)
θ n (x) = α ((x 2 -1 2 ) -x 2 1 ) 2 + (x 2 -1 2 ) 2 + x 2 1 -2 , (III .53)
with α = 10 -2 (see figure III.6). For this normal velocity, we expect the sequence Γ k to converge towards a U-shaped curve (see figure III.6, red). As an initial shape, we set Γ 0 = C(0, 1).

Let γ k i be an explicit discretisation of Γ k , then, using Euler's explicit scheme, we have:

γ k+1 i = γ k i + θ n (γ k i )n(γ k i ) . (III.54)
Let d k be the signed distance function of Γ k , then, using the method described in subsection 1.3, we have:

dk+1 = dk + θ * n , (III.55)
where θ * n is the constant-along-normals extension of θ n | Γ k , obtained with the method described in section 1.4. Finally, d k+1 is obtained by redistancing dk+1 (see section 1.4 again).

The figure III.7 shows the comparison between the explicit and implicit representations, at starting time of the sequence Γ k , and steady state. The mesh steps were chosen in such a way that the accuracy of the implicit and explicit numerical representations of Γ 0 are similar. We can see that the explicit representation fails to capture the steady state fully, even though all the control points seem to lie on the expected curve. The main issue here is the fact that these control points are concentrating in some areas and rarefying in others.

III.2 Boundary integrals & level-set representations

Γ Γ η η T Γ (x)
x

Boundary integrals & level-set representations

Beside the evolution of the curve though deformation, the other key ingredient of the shape optimisation algorithm is the resolution of boundary integral equations, and hence the calculation of boundary integrals. Kublik et al. in [START_REF] Kublik | An implicit interface boundary integral method for Poisson's equation on arbitrary domains[END_REF] proposed a formulation for constructing boundary integral methods to solve Poisson's equation with a level-set method. One of the perks of this method is to provide boundary potentials (and hence a shape gradient) that are readily constant along normals.

Thickening of Γ

In this section we present a method, initially presented in [START_REF] Kublik | An implicit interface boundary integral method for Poisson's equation on arbitrary domains[END_REF], which will allow us to approach an integral of the type:

I = Γ f ds , (III.56)
from the knowledge of the signed distance function d of Γ. We start by defining parallel curves Γ η :

Γ η = (Id + ηn)Γ (III.57)
where η ∈ R, and small enough for Id + ηn to be one-to-one (see figure III.8). The inverse mapping then writes Id -ηn.

Chapter III. The Level-set Method

Denoting T η = Id -ηn, we have from the change of variable formula in subsection 2.1 of chapter II:

Γ f ds = Γη f • T η J η ds , (III.58)
where, in 2D:

J η = |τ -η∇n T τ | (III.59)
From chapter I, we recall that:

∇n T τ • τ = H , (III.60) ∇n T τ • n = τ • ∇nn = 0 , (III.61)
hence, ∇n T τ = Hτ . This leads to:

J η = 1 -η H , (III.62)
which is non-zero for η smaller than the minimum radius of curvature of Γ, that we will call η * from now on. Since d -η is the signed distance function of Γ η , we have:

J η = 1 -d∆d . (III.63)
Let us now we compute the average of (III.58

) for η ∈ [-η * , η * ]. Let δ be a function R → R, supported in [-η * , η * ], satisfying the moment conditions: η * -η * δ(η)dη = 1, (III.64)
and:

η * -η * δ(x) x dx = 0. (III.65)
By the moment condition (III.64), we obtain:

I = Γ f ds η * -η * δ(η)dη = η * -η * Γη f • T η J η ds δ(η)dη
Recalling that, Γ η = d -1 (η), T η = Π Γ and J η = 1 -d∆d, and denoting f * = f • Π Γ we have:

I = η * -η * d -1 (η) f * (1 -dH) δ(d) dsdη .
(III.66)

III.2 Boundary integrals & level-set representations

Using the co-area formula, we obtain:

I = R 2 f * (1 -d∆d) δ(d(x)) |∇d| dx.
Finally, since |∇d| = 1, we obtain:

Γ f ds = R 2 f * (x) (1 -d(x)∆d(x)) δ(d(x)) dx . (III.67)
Remark 8. Our presentation, in particular the calculation of J η has been simplified by the fact that we considered smooth curves of the plane. Extensions to 3D curves and surfaces, with possibly corner and end-points exist, and can be found in [START_REF] Kublik | Integration over curves and surfaces defined by the closest point mapping[END_REF]. However, the particular case presented here is sufficient for what we need.

Numerical methods

In this section, we give the numerical methods to approximate the integral on boundary in section 2.1 i.e. equation (III.67). We denote by R a rectangular domain enclosing the obstacle Γ. As in subsection 1.4 this domain is discretised using a cartesian grid x ij = (x 1 i , x 2 j ), for i = 1..N 1 and j = 1..N 2 , where the discrete values x 1 i and x 2 j are respectively evenly spaced with the step h 1 and h 2 . For any function of R, say f , we will denote f ij an approximation of f (x ij ). We will also denote the centred finite-difference operator:

∇ h f ij = f i+1,j -f i-1,j 2h 1 e 1 + f i,j+1 -f i,j-1 2h 2 e 2 .
(III.68)

Finally, we denote the 5-points centered finite difference laplacian operator:

∆ h f ij = f i+1,j -2f ij + f i-1,j h 2 1 + f i,j+1 -2f ij + f i,j-1 h 2 2 . (III.69)
The projection of a grid point x i,j onto Γ reads:

Π Γ (x i,j ) = x i,j -d i,j ∇ h d i,j , (III.70)
which is not a grid point. Hence in order to obtain an approximation of f * (x i,j ), that we denote f * i,j , one solution could be to approximate f at Π Γ (x i,j ) by using 2D interpolation, which are conveniently implemented in Matlab. Another solution would be to extend f in a constant-along normals manner, as explained in subsection 1.4. Further, we use the averaging kernel as in [START_REF] Kublik | An implicit interface boundary integral method for Poisson's equation on arbitrary domains[END_REF]:

δ(x) =        1 2η * (1 + cos( πx η * )) if |x| ≤ η * 0, else. (III.71)
Finally the integral (III.67) is approximated using the Riemann sum as follows:

Γ f ds ∼ h 1 h 2 N 1 i=1 N 2 j=1 f * i,j δ(d i,j ) J i,j . (III.72)
where

J i,j = 1 -d i,j ∆ h d i,j . Since δ is supported in [-η * , η * ],
it is clear that the values δ(d i,j ) are non zero only for a certain amount of indices i and j. We denote

N = {(i, j) | δ(d i,j ) = 0} . (III.73)
Denoting m k for k = 1..N N the elements of N, we have:

Γ f ds ∼ h 1 h 2 N N k=1 f * m k δ(d m k ) J m k , (III.74)
which we will rewrite, with a small abuse of language:

Γ f ds ∼ h 1 h 2 N N k=1 f * k δ(d k ) J k . (III.75)

Convergence of the integration method

Let us denote:

g(x) = f * (x)δ(d(x))(1 -d(x)δd(x)) (III.76)
then, from (III.67), we have that R g = Γ f , and so, from the fact that the method of rectangles approaches the integral on R with a second order accuracy with respect to the mesh size h (we consider h = h 1 = h 2 here for the sake of simplicity), we have:

Γ f ds -h 2 N N k=1 f * k δ(d k ) J k ≤ c h 2 ||D 2 g|| ∞ , (III.77)
where c is a positive constant, and D 2 g is a matrix made of all second order derivatives of g (see appendix for the details). So, it is clear that the approximation (III. [START_REF] Weinblum | Applications of wave resistance theory to the problem of ship design[END_REF]) is (at least) of second order accuracy with respect to h. However we remark that x → δ(d(x)) depends on the thickness 2η * of the tubular neighborhood of Γ considered (see equation (III.71)). For this reason, the above statement holds only when η * is given and fixed, so the number of points involved in the calculation of (III.75) grows like N 2 . In order to reduce the computational cost of the method, it is interesting to investigate the case where η * is scaled with the mesh size h. We then have no guarantee that (III.77) will ensure even the convergence of the method, since, as we can deduce from (III.71), we have δ ∼ h -3 in this case.

The seminal article of Engquist et al. [START_REF] Engquist | Discretization of Dirac delta functions in level set methods[END_REF] provides an interesting insight of this phenomenon. The authors consider the simple test case of the calculation of the length of a straight line that is not necessary parallel to the grid axes. This length correspond to the calculation of the integral:

I = Γp,q 1 ds , (III.78)
where Γ p,q = Span((p, q)) ∩ R, and (p, q) ∈ N 2 have no other common denominator than 1.

The signed distance function writes in this case:

d(x) = -px 1 + qx 2 √ p 2 + q 2 .
(III.79)

Since ∆d = 0, the approximation of I with (III.75) writes:

Ĩ = h 2 N N k=1 δ(d(x k )) . (III.80)
In this particular case, when we set η * = βh, the authors of this article prove that the sum (III.80), converges to the correct value of I only when:

β = p + q √ p 2 + q 2 , (III.81)
in other words, when β = |∇d| 1 . The general case remains an open problem to our knowledge, but the following tests shows that this result is robust to situations where Γ is not a straight line.

A simple test

Here we describe an example of calculating a perimeter of circle with radius 1 and center 0.

I = C(0,1)
ds .

(III.82)

We have here

f (x) = 1. We take R = [- √ 2, √ 2] × [- √ 2, √ 2
] (to avoid points of Γ lying exactly on the mesh nodes), and denote N and h respectively the number of points in each direction and the mesh element size. The signed distance function reads:

d(x) = (x 2 1 + x 2 2 ) -1 , (III.83)
We follow the method in section 2.2 to compute the numerical approximation of the length of the circle. The approximate value of the integral then writes:

Ĩ = h 2 N N k=1 J k δ(d k ) , (III.84)
where k is the index of the points that lie in the chosen tubular neighbourhood of Γ with a thickness η * . Since I = 2π, the relative error of the approximation III.84 writes:

E = | Ĩ -2π| 2π . (III.85)
We present here the convergence results for different thicknesses of the tubular neighbourhood. We consider three cases:

η * = constant, η * x 1 ,x 2 = 2h|∇d(x 1 , x 2 )| 1 and η * = 2h (see Figure III.9

for an illustration of how the thickness of the tubular neighborhood changes with h).

The figure 2.3 represent the convergence results for these three different cases. We notice that the best result is obtained when we set η * constant (first case), where the convergence rate is better than second order. For fairness, this case should be set aside since is computational cost is not comparable to the two others. We remark that in the second case η * = 2h|∇d| 1 we recover a fairly clean first order convergence, while in the third case η * = 2h the trend is not very clear. This confirms the intuition of Engquist et al., and shows that one should in principle choose η * = 2h|∇d| 1 to have both efficient and accurate approximation of boundary integrals with (III.75).

Application to boundary integral equations

In this section we present an application of the method presented in section 2 to the resolution of the Laplace equation with boundary integral equations. To be consistent with the application we are aiming for (the Neumann-Kelvin problem), we choose the exterior Neumann problem as our example. Note that this application was already presented in [START_REF] Kublik | An implicit interface boundary integral method for Poisson's equation on arbitrary domains[END_REF] and [START_REF] Engquist | Discretization of Dirac delta functions in level set methods[END_REF] notably (albeit for interior problem, but the difference is minor). We consider the Neumann problem: Using the same argument as in chapter I, (see for instance [START_REF] Hsiao | Boundary Integral Equations[END_REF]46])), α has to solve the boundary integral equation:

                 ∆u x = 0 , in Ω + , ∂u(x) ∂n x = g(x) , on Γ ,
Γ α(y) ∂ n(x) G(x, y) ds y - 1 2 α(x) = g(x)
, for all x ∈ Γ, (III.92)

where:

∂ n(x) G(x, y) = x -y |x -y| • n(x) (III.93)
Further, we use the equation (III.90) to recover u in Ω.

Numerical approximation

We consider the same notations as in subsection 2.2, and introduce the points (x k ) k=1..N N of the grid that lie inside the tubular neighborhood of thickness η * . We also denote x * k the projection of x onto Γ obtained with:

x * k = x k -d k ∇ h d k (III.94)
The approximation of α(x * k ) will be denoted d k . Following (III.75), the boundary integral equation (III.92) is approached by:

h 1 h 2 N N k=1 α k (∂ n(x) G) m,k J k δ k -1 2 α m = g m , (III.95)
where g m = g(x * m ), and δ k = δ(d k ). The kernel H m,k should be, in principle given by (we recall that n = ∇d):

(∂ n(x) G) m,k = x * m -x * k |x * m -x * k | • ∇ h d k . (III.96)
However, when |x * m -x * k | becomes close to zero, the evaluation of III.96 involves ratios of quantities which are of the same order as the truncation error of the method (it is even undefined for k = m). To circumvent this issue the authors of [START_REF] Kublik | An implicit interface boundary integral method for Poisson's equation on arbitrary domains[END_REF] expanded H(x, •) at first order around x, to obtain the following approximation:

∂ n(x) G(x, y) ∼ 1 4π H(x) , for x ∼ y , (III.97)
and replace H by its approximation when |x * m -x * k | is smaller than a fraction of h:

(∂ n(x) G) m,k = 1 4π ∆ * h d k , for |x * m -x * k | < (III.98)
Now we rewrite (III.95) as a linear system:

D α = g , (III.99)
where:

D = (h 1 h 2 H m,k J k δ k )m=1..N N k=1..N N -1 2 I N N , α = (α m ) m=1..N N , g = (g m ) m=1..N N , (III.100)
where I N N is the identity matrix of size N N . The system III.99 is then solved using an iterative method called the bi-conjugate gradient method, which is a variation of the conjugate gradient method that allows the matrix of the system to be non-symmetric (see [73] for more details on iterative methods). Note that, because of the strong diagonal component of the matrix D, the system is well-conditioned. This is the reason behind the choice of a singlelayer representation of the solution in (III.90), which leads to the resolution of a Fredholm equation of the second kind (III.92). Now, once α is known, if we want to recover an approximation of the solution u at a given point x, we simply approach (III.90) with:

u(x) = h 1 h 2 N N k=1 α k G(x, x * k ) J k δ k . (III.101)
In practice, for the following section, we will recover the solutions at the nodes x i,j of the same grid that was used for the definition of the level-set function, but note that this is not mandatory, and we can request any point of the domain Ω + with (III.101).

Validation tests

We consider a particular test-case in which Γ = C(0, 1), and

g(θ) = -n cos(nθ) . (III.102)
for n ∈ N. The solution of the exterior Neumann problem for this case is known from the method of separation of variables:

u(r, θ) = r -n cos(nθ) . (III.103)
Using the method detailed in subsection 3.1, we obtain results which we present in Figure III.11. In these figures, we plot the numerical, and exact solution for points of the grid in the exterior domain, along with their difference (in logarithmic scale). We remark that we have a fairly good agreement between the exact and approached solution. Let us also note that when n = 1, the solution we get is exactly the solution to the potential flow problem on a circular cylinder (see the Figure III.11, first row), since in this case g = n • e 1 .

In order to give a validation of the method, we compare the exact value of α with the numerical values. Our error criterion reads:

E = N N k=1 |α(x * k ) -α k | N N k=1 |α k | (III.104)
The exact value comes from the definition of α using the jumps relations:

α(θ) = ∂ r u + (1, θ) -∂ r u -(1, θ) = 2n cos(n θ) , (III.105)
where u + is the solution to the exterior Neumann problem, and u -is the solution of the interior problem, which was also obtained by the method of separation of variables. Here we chose to validate the method by calculating the error on the boundary potential α rather than on the solution u as it is usually done. The reason behind this is that in the following chapter, the shape gradient that we are going to use is defined directly by α. The main drawback is that we cannot get a manufactured solution for α as easily as it is the case for u. This is the reason behind the particular choice of Γ and g we have made in (III.102) and (III.103).

In the tables III.1, III.2 and III.3, we present the error E as a function of the resolution of the grid on which the level-set function is defined, for different thicknesses of the tubular neighborhood, and different values of N . As expected, the case η * = 2h|∇d| 1 converges with a clean second order convergence rate. Surprisingly, the rate of convergence of the error for the case where η * is fixed is lower than in the other cases. Also the case η * = 2h which should not be convergent shows a second order convergence. All these unexpected results might be a consequence of the choice of the test-case, in which the solution α is an eigenvalue of the boundary integral operator (g and α are colinear). Even though the solution is not an eigenvalue for the discrete operator, this particular choice (motivated by the knowledge of the exact solution) might be the cause of the strange phenomena observed.

D * : u → Γ ∂ n(x) G(x, •)u(x) ds x -1 2 u , (IV.3)
We introduce their numerical counterparts S, D, and D * defined by the numerical method presented in section 3 of chapter III: Further, we recall that:

S = (h 1 h 2 G m,k J k δ k )m=1..N N k=1..N N , (IV.4) D = (h 1 h 2 (∂ n G) m,k J k δ k )m=1..N N k=1..N N -1 2 I N N , (IV.5) D * = (h 1 h 2 (∂ n G) k,m J k δ k )m=1..N N k=1..N N -1 2 I N N , ( IV
n k = ∇ h d k , τ k = R( π 2 ) n k , H k = ∆ * h d k . (IV.7)
Where we recall that ∆ * h u is the interpolation of ∆ h u k at point x * k . For various purpose, it might be useful to calculate the tangent derivative of a function defined on Γ. Let u be such a function, and u • Π Γ its "thickened" version on the neighbourhood of Γ. Then we approach the tangent derivative of u at x k with:

∂ τ,h v k = τ k • ∇ * h (u • Π) k , (IV.8)
where, again, we recall that ∇ * h v k is the interpolation of ∇ h v k at point x * k . Also in all the following presentation, we identify the linear mapping (v k ) k=1..N N → (∂ τ,h v) k=1..N N with its matrix by a small abuse of notations.

Calculation of the wave-making resistance

The goal of this section is to present how the aforementioned numerical method can be used in order to calculate the wave making resistance of an arbitrary (smooth) obstacle in a constant flow by using the Neumann-Kelvin model. We recall from chapter I that the wave-making resistance reads:

R w (Γ) = ρ ν 4 Γ α Γ (x) E(x)ds x 2 , (IV.9)
where the boundary potential α solves:

D[α] = -U ∞ n • e 1 , (IV.10)
and:

E(x) = e ν(ix 1 +x 2 ) . (IV.11)
Using the numerical method from chapter III, and the notation from section 1, we obtain the following approximation for the wave-making resistance:

R w (Γ) ∼ ρ ν 4 h 1 h 2 N N k=1 α k E k J k δ k 2 , (IV.12)
where E k = E(x * k ), and α = (α k ) k=1..N N solves the linear system:

D α = -U ∞ n 1 , (IV.13)
where n 1 = (n k • e 1 ) k=1..N N , and D is given by equaiton (IV.5).

Let us consider the case where the obstacle is a circular cylinder of radius a and depth f , and compute the wave-making resistance vs Froude number profile. Here the Froude number serves as a dimensionless version of the velocity, and it is defined by: Fr = U ∞ / gf , (IV.14) so the chosen reference space scale is the depth of the obstacle. We have at our disposal an exact formula for the wave-making resistance in the Neumann-Kelvin model in the case of any submerged circular cylinder. Remarkably, this solution was discovered more than a century ago by Sir Thomas Henry Havelock [START_REF] Havelock | The forces on a circular cylinder submerged in a uniform stream[END_REF]. The "recipe" for the calculation of the wave-making resistance through Havelock's solution is presented in appendix C.

Let us now compare the exact and the numerical profiles obtained with (IV.12) and (IV.13) for a circular cylinder of depth 0.7m and radius 0.5m (to be consistent with the examples given in the historical article of T.H. Havelock). To this end, we introduce the level-set function

d(x) = |x + f e 2 | -a , (IV.15)
and consider a rectangular domain:

R = - √ 5 2 a, √ 5
2 a × - √ 5 2 a -f, √ 5
2 a -f , (IV.16) enclosing the circular boundary, which will be discretised using a 128 × 128 grid. The figure IV.1 presents these R w vs Fr profiles (left exact, right numerical). We recover a classical feature of the wave-making resistance, that is, its preeminence in the regime Fr = 0.3 ∼ 1.5.

The figure IV.1 also shows the good agreement between the exact and numerical approximation with a relative error less than 0.2% at its maximum (the errors shown here are relative errors). Further, we investigate the rate of convergence of our numerical method with respect to h, for different values of the velocity. The convergence results shown in table IV For practical purposes, besides the wave-making resistance, it might be useful to recover the velocity potential Φ in the exterior domain Ω + and the free-surface deformation function η. Recall from chapter I that: and: η(x 1 ) = U g ∂ 1 Φ(x 1 , 0) , (IV.18) which, using the numerical method previously described translates as:

Φ(x) = Γ G(x,
Φ(x) ∼ h 1 h 2 N N k=1 G(x, x k ) α k J k , (IV.19)
and:

η(x 1 ) ∼ h 1 h 2 U g N N k=1 ∂ x 1 G((x 1 , 0), x k ) α k J k . (IV.20)
Note that here the "query point" x (respectively x 1 ) can be any element of Ω + (respectively . This grid is hence independent of the one chosen to compute α, and is also represented on figure IV.3. For η, we chose 300 values evenly spaced between -9 and 3. These illustrations show the characteristics of the wake, in which both the amplitude and the wave-length of the wake depend on the Froude number. Note that the free surface deformation in figure IV.3 does not represent the actual upper boundary of the fluid domain Ω + which is flat because that we used a linearised model of water-waves. The free surface η has hence to be seen just as a function that "lives" on this flat upper bound. For this reason, it can happen that the free surface represented crosses the obstacle without causing any problem in the model (both from the theoretical and numerical point of view), as long as the obstacle chosen is not piercing the top of Ω + .

R
8 -6 -4 -2 0 X (m) -2 -1 0 1 2 Y (m) -8 -6 -4 -2 0 X (m) -2 -1 0 1 2 Y (m) -8 -6 -4 -2 0 X (m) -2 -1 0 1 2 Y (m) -8 -6 -4 -2 0 X (m)

Discretisation of the problem

Let us now explain how the normal displacement w is approached from the knowledge of approximate values (d k ) k=1..N N of the signed distance function d. We start by approaching the state and adjoint equations (IV.22)-(IV.23), the same way as in section (2):

D α = -U ∞ n 1 , (IV.27) D * q = M 1 α , (IV .28)
where:

M 1 = -2 Re(E m E k ) J k δ k m=1..N N k=1..N N . (IV.29)
Again using the integration method described in chapter III, we obtain, from (IV.25)

w = α * M 2 α + ∂ τ,h q * ∂ τ,h S[α] + α * ∂ τ,h S * [∂ τ,h q] + U ∞ ∂ τ,h q * τ 1 . (IV.30)
where * denotes the element-by-element product, and:

M 2 = 2 Re[(∂ n E) m E k ] J k δ k m=1..N N k=1..N N . (IV.31)
The state α and adjoint state q are obtained by solving (IV.27) and (IV.28) and the approximation of the operator S defined in equation (I.166) in chapter I, reads:

S * = (h 1 h 2 G * k,m J k δ k )m=1..N N k=1..N N . (IV.32)
Finally we recall that ∂ τ,h is the approximation operator for the tangent derivative (see section 1).

Grid resizing

In all the above, we supposed that the rectangular domain R enclosing Γ that we use for discretisation is fixed. However, in our algorithm the shape of the obstacle evolves, and it is difficult (and a poor choice) to anticipate the size of the final obstacle and to define R such that Γ n ⊂ R for all n ∈ N.

For this reason, we will change R over the iterations of our algorithm in such a way that it always encloses Γ n as tightly as we will require. In order to avoid additional interpolation, the changes we will apply to R will consist in adding or removing mesh elements layers on the right, left, top and bottom sides of R. These element are added and removed in such a way that the set of points of the tubular neighbourhood associated with our evolving curve Γ is at a distance of at least ph and at most (p + 1)h to the border of the box, where the parameter p ∈ N defines the thickness of the "safe zone" chosen.

The figure (IV.4) shows the evolution of the computational domain though this process, on a particular case of the evolution of a curve with a normal displacement given by:

θ n (x) = 4 -(8x 2 1 + x 2 2 ) , (IV.33)
and starting from a circle. In this case the steady state (plotted in black) is expected to be an ellipse, which cannot be represented properly with the initial cartesian mesh (first picture on the left). Through the iterations of this normal motion, our method adds elements on the top and bottom and removes elements on the right and the left so that the safe zone (in magenta) is tightly respected.

Measure constraint

We recall the algorithm defined in chapter II that we will use to enforce the measure constraint on the obstacle: Data: A shape Γ, a tolerance , a target surface area v and a step δ.

Result:

The "projected" Γ = (Id + an) Γ, where a is such that µ(Ω -) = v. Let us now suppose that the signed distance function d of Γ is known. The constant normal displacement "Γ := (Id + a n) Γ" can be translated in terms of d as "d := d + a". We hence simply modify d by adding a constant that has to be chosen in such a way that the measure of the domain it encloses has the required value. The most notable feature of this is that, after this modification, d remains a signed distance function. Further, we compute µ(Ω -) by using a smooth approximation of the indicator function of the domain:

µ(Ω -) ∼ μ(d) = h 1 h 2 N 1 i=1 N 2 j=1 1 2 1 -erf d(x i,j ) h (IV .

Presentation of the results

Let us first study the influence of the velocity on the optimal shapes obtained with the algorithm that we outlined in figure IV.5, and by using the parameters defined in section 3.4. The depth is fixed at 0.7m, we recall that the radius of the initial shape is 0.5m, and the velocities are calculated in such a way that the values of the Froude number (defined by (IV.35)) are 10 values evenly spaced between 0.3 and 1.4, which is the range of values for which the wave-making resistance phenomena are the most important. The figure IV.6 presents the results of the shape optimisation algorithm for some chosen values of the Froude number : 0.4222, 0.6667, 0.9111, 1.0333, 1.2778. As an indicator of the wave-making resistance, we also plot the wake generated by these obstacles (we recall, from (I.131) that the wave-making resistance is proportional to the square of the amplitude of the wake far downstream). We compare the initial shape/wake in blue with the final shape/wake, in red. The reduction in the wake amplitude is quite remarkable: in many cases, the wake downstream is indistinguishable from a flat profile. The ratio between the wave-making resistance for the final shape and the wave-making resistance for the initial shape is presented in the table IV.2, and we notice that in almost all cases, the wave-making resistance is reduced by several orders of magnitudes (one in the worst case, more than six in the best case).

Let us now consider the Froude number to be fixed, and let us take four different values of the depth : 0.7000, 0.8444, 0.9889, 1.2778, and Fr = 0.6667. In figure IV.7 we plot our results the same way as in figure IV.6, with these new parameters. As it was the case in figure IV.6, we remark that the choice of the depth has an influence on the typical length of the obstacle, though the influence of the depth seems less important than the influence of the Froude number. To finish, we present in figure IV.8 a summary of all the optimised shape we obtained, presented on an array of parameters conveniently chosen: the depth for the y-axis, and (U 2 ∞ /g)(= Fr 2 /f ) on the x-axis, the later being the typical wave-length for the wake in a flow predicted with the Neumann-Kelvin model. Each shape with a different colour represents a run of the shape optimisation algorithm with parameters given by its position on the grid.

Fr

The triangular symbols mark parameters (f, Fr) for which the algorithm fails to converge, and the sequence of obstacles obtained is a sequence of deeper and deeper obstacles. We refer to these cases as "sinking cases". This is a situation that was anticipated in chapter II, subsection 4.3 in which we remarked that the wave-making resistance should go to 0 as the depth of the obstacle goes to ∞. Also, as we can see in the algorithm depicted in figure IV.5, the measure constraint is taken into account to avoid shrinking obstacles, but nothing is done to prevent the obstacle from sinking. It is all the more surprising to find cases where the algorithm actually converges. We will try to give a tentative explanation for the existence of these local minima that remain close to the surface in subsection 5.2. -All our shape optimisation cases. We represented each optimised shape on a grid that represent the two parameters we have chosen to study (we however give a length scale for the obstacles). Each colour represents a different simulation. The cases for which our algorithm fails to converge (the so-called sinking cases) are represented with triangular markers.

(a)

-8 -6 -4 -2 0 X (m) -2 -1 0 1 2 Y (m) (b) -8 -6 -4 -2 0 X (m) -2 -1 0 1 2 Y (m) (c) -8 -6 -4 -2 0 X (m) -2 -1 0 1 2 Y (m) (d) -8 -6 -4 -2 0 X (m) -2 -1 0 1 2 Y (m) (e) -8 -6 -4 -2 0 X (m)
(a) -8 -6 -4 -2 0 X (m) -2 -1 0 1 2 Y (m) (b) -8 -6 -4 -2 0 X (m) -2 -1 0 1 2 Y (m) (c) -8 -6 -4 -2 0 X (m) -2 -1 0 1 2 Y (m) (d) -8 -6 -4 -2 0 X (m) -2 -1 0 1 2 Y (m)

Discussion

In this section we discuss the results presented in section 4, and try go give an explanation for the optimal shapes we obtained.

The length of the obstacles

When examining figure IV.8, one striking feature of the optimal shapes is that their typical length seem to fit linearly with the typical wave-length U 2 ∞ /g. This can be seen in figure IV.9, in which each point represents the length of an optimal obstacle obtained for a different value of the Froude number. We recall U 2 ∞ /g = Fr 2 /f , hence the influence of the Froude number on the length of the obstacle we obtain is quadratic. We try to give an explanation for this by using an analogy with the idea of the Rankine body. It is a classical result of fluid mechanics that, when we consider a 2D potential flow in an infinite domain, it is possible to imitate an obstacle by considering a couple of source/sink placed on an axis parallel to the flow at infinity (see for instance [START_REF] Batchelor | An Introduction to fluid dynamics[END_REF], pp. 458-460). Further, the strength of this sink/source couple determines the volume of the body. The "virtual" obstacle thus obtained is called a Rankine Body. In this case, the velocity po-

IV.5 Discussion

tential perturbation is the sum of two Green's functions, one for each source position:

Φ(x) = s (G(x, y + 1 2 le 1 ) -G(x, y -1 2 le 1 )) , (IV.36)
where s is the strength associated with the source/sink couple, and l is the length separating the source and the sink, and of course, G is the Green's function of the Laplace equation. We plot in figure IV.10 an instance of such a flow. It is important to note that this idea does not hold in general for a flow with a free-surface. Let us however consider a Rankine "pseudobody" produced by a source and a sink in the Neumann-Kelvin problem:

Φ(x) = s (G(x, y + 1 2 le 1 ) -G(x, y -1 2 le 1 )) . (IV.37)
The figure IV.11 shows the instance of such a flow (there is no closed streamline in this case). From (IV.37), we could try to find which value of l would produces the minimal wave-making resistance. As we have seen in equation (I.131) of chapter I, the wave-making resistance is proportional to the square of the amplitude of the waves far away downstream. Luckily, for the case of a flow produced by a source (or a sink), the behaviour of the solution at infinity is known (see (I.84)), and, by linearity, the contribution of the source and the sink can be subtracted. From (I.84), we obtain that behaviour of the free surface elevation η(x) = -(U ∞ /g)∂ x Φ(x 1 , 0) for such a Rankine pseudo-body, far away downstream is:

η ∞ (x 1 ) = 1 U ∞ cos(ν(x 1 + 1 2 l)) -cos(ν(x 1 -1 2 l)) = - 2 U ∞ sin(νx 1 ) sin(ν 1 2 l) . (IV.38)
Hence, the wave-making resistance is proportional to: sin(ν 1 2 l) 2 , and the smallest length for which this resistance is equal to zero is:

l = 2 π ν = 2π U 2 ∞ g . (IV.39)
While this relation does not fit what we observe in IV.9, it could explain the apparent linear relation observed between the length of the obstacle we obtain, and the parameter U 2 ∞ /g. The discrepancy between the slopes could be explained by the fact that the length of such a Rankine pseudo-body is not equal to the distance between the source and the sink, and depends on the strength s (it is in fact difficult to even define properly).

The sinking cases

In this section we attempt to give an explanation for the fact that our algorithm finds locally optimal shape of finite depth. In order to simplify the situation, we will consider a shape optimisation problem in which the only degree of freedom available is the depth. Starting with a circle, the shape remains a circle, and so the wave-making resistance in this case can be obtained with the formula of Havelock. The figure IV.12 shows a plot of the wave-making resistance for a circle of radius 0.5m as a function of the depth an the aforementioned parameter U 2 ∞ /g. We remark that, for some values of U 2 ∞ /g, there is a depth for which the wave-making resistance is maximal. This means that starting above this maximal, the gradient descent method should move the circle closer to the surface, while starting below will lead to a "sinking" sequence of domains. This interpretation is not entirely satisfactory because the two distinct regions above do not fit with what we see in figure IV.8. This is mainly due to the fact that despite the initial shape is a circle, it quickly deforms into a very different shape for which our analogy does not hold anymore. A better solution would be to allow linear deformations such as extensions along the horizontal axis. We also note that this dichotomy between optimal shapes that remain close to the surface and "sinking" shapes has already been observed in [START_REF] Dambrine | Regularity of optimal ship forms based on michell's wave resistance[END_REF], on a 3D version of the Neumann-Kelvin with a "slender" approximation of the shape.

Wave making resistance vs Froude profiles

We start by remarking that each shape that was previously obtained is designed to be optimal for a given value of the Froude number, that we will call here Fr * . We expect the value of the wave-making resistance to be minimal for Fr * , but we have no information on how each optimal obstacle behave for other values of Fr. In particular, for ship design (the 3D surface piercing case), it is interesting to know wether the optimal regime Fr * can be reached starting from 0 without experiencing a peak of resistance that would be too difficult to overcome.

The figure IV.13 shows plots of the wave-making resistance vs Froude number profiles for different optimal obstacles that were obtained for f = 0.7 and different values of Fr * . Two regimes can be observed: below a critical value Fr * 0.8, the optimal regime shows no large peak before Fr * ; above this critical value, a very high peak of resistance has to be passed before reaching Fr * . Again, this type of behaviour is consistent with what was observed in [START_REF] Dambrine | A theoretical and numerical determination of optimal ship forms based on michell's wave resistance[END_REF] with "slender"-type models. Froude number Fr (indicated with a dot). The redline is the profile for the optimal shape, the blue line indicates the profile for the initial shape. Left column: Fr * < 0.8; right column Fr * > 0.8.

Applicability for the analogue gravity experiments

As stated in the introduction of this manuscript, despite the fact that the main motivation for the study of this problem comes from ship design, the problem we studied here is two dimensional, and cannot be transposed immediately to naval engineering. However, we hope that the methods developed here can be extended in 3D without too many efforts. Despite all this, there is an application in physics for which our work could be relevant: analogue gravity experiments. In these experiments, the behaviour of water-waves in spatially variable currents are used as an analogue for the propagation of light in a region of spacetime deformed by the presence of a black-hole. In particular, G. Rousseaux, L.-P. Euvé and collaborators (see [START_REF] Euvé | Observation of noise correlated by the hawking effect in a water tank[END_REF]) at the Institut Pprime (CNRS/Université de Poiters/ENSMA/ENSIP), have recently observed an analogue of the Hawking effect 1 in an observation of the blocking of waves generated downstream of a variable current. Here, the variation of the current is obtained by the presence of an obstacle in the flow. One of the many technical challenges is to design an obstacle that produces a large enough and well controlled change of velocity, but for which the wake is as small as possible. The reason for this is to reduce the interaction between the wake and the waves that are purposely produced downstream. The figure IV.14 shows a complete view of the flow around the obstacle previously obtained by fixing the initial depth at 0.7m, and changing the Froude number (see figure IV.14). First, we remark that the wake is indiscernable, as is was shown in the previous pictures. Also, in each case, an acceleration of the fluid around the obstacle can be observed. At the free surface, an increase of about 20% of the velocity is observed just above the obstacle. The obstacle shapes we obtained hence seem like good candidates for this application. 1 The Hawking effect is a radiation theoretically emitted by black-holes that was predicted by S. Hawking in 1974 [START_REF] Hawking | Black hole explosions ?[END_REF], and is also known by the name: black hole evaporation. This effect is so tenuous that is was never observed directly. This has motivated the design of analogue experiments using water-waves or acoustic waves. eral given depths and velocities. We found that there is a huge reduction of the wave-making resistance between the initial shape and the final shape, with a ratio varying between 1 and 6 orders of magnitude in the best case. One striking feature of the (converged) optimal shapes is that their typical length seems to fit linearly with the typical wave length of the problem.

By imposing the area of the shape, we prevented the obstacle to shrink. But we did not impose the center of gravity of the object, so that in several cases, our algorithm produced a sequence of shapes going down far away from the surface, with a wave-making resistance going to zero (these are what we called the sinking cases). It was all the more surprising to obtain optimal shapes not far away from the free surface.

A shape which is optimal for a given velocity is not necessarily optimal for all velocities. By considering the wave resistance vs velocity, we observed that sometimes, the optimal regime cannot be reached without experiencing a peak of resistance. This can be a serious limitation in application to ship hydrodynamics.

Perspectives

After the shape optimisation for the wave-making resistance problem that we presented in this thesis, we propose several tasks which could be done to enhance the theory, the numerical approach, as well as the applications related to this problem.

From a theoretical point of view, it would be interesting to prove the existence of an optimal shape for our problem. This looks like a difficult question. Indeed, for a given obstacle, existence and uniqueness of the potential holds only except for a finite number of velocities. It require smoothness of the boundary and there are topological constraints. A possibility could be to start with the case of thin obstacles, by seeking a formula similar to Michell's.

We note that we calculated a shape derivative by a Lagrangian method which is a formal computation [START_REF] Allaire | Conception optimale de structures[END_REF]. It would be interesting to give a functional setting for this computation.

From a numerical point of view, we would like to extend the method presented in this thesis to the 3D setting, in the case of a surface-piercing body. This would be closer to applications in ship design. We could start from the fundamental solution of the NK-problem given by Delhommeau in his thesis [START_REF] Delhommeau | Les problèms de diffraction-radiation et de résistance de vagues: étude théorique et résolution numérique par la méthode des singularités[END_REF]. We could also study first the easier case of a 3D fully submerged body, or the 2D problem with a surface-piercing body. For a given obstacle, the NK-problem has been studies for these last two situations in [START_REF] Kuznetsov | Linear Water Waves[END_REF].

In view of applications to analogue gravity experiments, it would also be interesting to

We introduce F y = ∂ 1 Hy,2 , then F y solves:

   ∂ 2
11 F y (x 1 , 0) + ν∂ 2 F y (x 1 , 0) -ε∂ 1 F y (x 1 , 0) = ε∂ 2 11 G y (x 1 , 0) ∆F y = 0 .

(6)

Using the same methods and notations as in chapter I subsection 2, we get:

Hy,1 (x) = 1 π Re R +
e -ikω(x,y) k -ν -iε dk (7)

F y (x) = - ε π Re R + e -ikω(x,y) k -ν -iε dk . ( 8
)
Let us denote:

I = R + e -ikω k -ν -iε dk . ( 9
)
We apply the change of variables k ← k -iε, then:

I = e -iεω R + -iε e -ikω k -ν dk (10)
Let us note in the equation [START_REF] Brard | Le problème de Neumann-Kelvin[END_REF] above, when ε → 0, then the term e -iεω tends to 1. Hence, we have:

Hy,1 (x) = 1 π Re R + e -ikω(x,y) k -ν dk (11) (12)
Recalling equation (I.62) (without writing the limit), and coupling with equations [START_REF] Brebbia | Boundary element methods for potential problems[END_REF], we obtain:

Hy,1 (x) -H y (x) = 1 π Re R + -iε e -ikω k -ν dk - l - e -ikω k -ν dk (13)
Using the Cauchy's integral theorem (see the path of integration in figure A.1), we have:

Hy,1 (x) -H y (x) = 1 π Re -iε 0 e -ikω k -ν . ( 14
)
When ε → 0, . Hence, in the limit ε → 0, we have:

-iε 0 e -ikω k -ν → 0 R iR ν l - ε R -iε -iε
Hy,1 (x) = H y (x). (15)
Moreover, from the finiteness of I when ε → 0, we obtain that F y (x) → 0 as ε → 0. Hence, Hy,2 (x) = c for any arbitrary constant is a solution for Hy,2 . The Green's function we calculate being defined up to a constant, by choosing c = 0, we have, in the limit ε → 0: Hy (x) = H y (x).

(16)

In conclusion, integrating k → e -ikω k-ω on P - ε for ε → 0 in equation (I.62) gives the same result as considering an evanescent Rayleigh dissipation.

B Fourier transform of -(x-a) 2 +b 2 ((x-a) 2 +b 2 ) 2

Let us consider:

f (x) = -(x -a) 2 + b 2 ((x -a) 2 + b 2 ) 2 . ( 17
)
The Fourier transform of f is given by the integral:

f (k) = R f (x)e -ikx dx = R -(x -a) 2 + b 2 ((x -a) 2 + b 2 ) 2 e -ikx dx . ( 18
)
First, we apply the change of variables x ∈ x + a:

f (k) = R -x 2 + b 2 (x 2 + b 2 ) 2 e -ik(x+a) dx, = e -ika R -x 2 + b 2 (x 2 + b 2 ) 2 e -ikx dx.
Let us denote:

g(x) = -x 2 + b 2 (x 2 + b 2 ) 2 e -ikx . (19)
We first notice that g, when seen as a complex-valued function, is holomorphic on C\{i |b| , -i |b|}.

Moreover, the decay of g at infinity along the imaginary axis depends on the sign of k. First let us consider k > 0. In this case, we have e -ikx ≤ 1 for Im(x) ≤ 0 Let us consider a path P - R consisting in the boundary of the intersection of a disc of center -i |b| and radius R and the negative imaginary half-plane, oriented clockwise (see figure A.2).

From the theorem of residues, we have:

P - R g(x) dx = -2πi Res(g, -i |b|) . ( 20
)
On one hand we have:

P - R g(x) dx = R cos α -R cos α g(x) dx + A - R g(x) dx, (21)
where A - R denotes the arc of the circle of center -i |b| and radius R, between the angles α and -π -α, paramatrized by: We also have:

x = Re -iθ -i |b| for θ ∈ [-α, π + α] . ( 22
g(x) = -R 2 + 2 |b| 2 R 4 e -ikx . ( 23
)
So, it is clear that |g(x)| = O( 1 R 2 ) on A - R . Hence, taking the limit R → ∞,

A - R g(x) dx ≤ c R 2 A - R ≤ c R → 0.
This leads to, since α → 0 for R → ∞:

lim R→∞ P - R g(x) dx = R g(x) dx . ( 24
)
On the other hand, we have, since g has a double pole at -i |b|:

C Havelock's solution for a circular cylinder

As we discuss in chapter I, we present the condition of the fluid flows through a submerged body under a free surface condition with a uniform velocity that is called the Neumann-Kelvin problem (NK problem). Let us recall this problem in equations (I.42):

                 ∂ 2 11 Φ + κ 0 ∂ 2 Φ = ε ∂ 1 Φ , on S 0 , ∆Φ = 0 , in Ω + , ∂ n Φ = -U ∞ n • e 1 , on Γ , |∇Φ| → 0 , for |x| → ∞ , ( 34
)
with κ 0 = g/U 2 ∞ . Sir Thomas Havelock in [START_REF] Havelock | The forces on a circular cylinder submerged in a uniform stream[END_REF] gave the analitical value of the solution for the two dimensional problem [START_REF] Havelock | The collected papers of Sir Thomas Havelock[END_REF] with Γ is a circular cylinder. He writes the wave reistance in the form of infinite series of parameters.

R w = 4π 2 κ 2 0 a 4 f (κ 0 ) f * (κ 0 )e -2κ 0 f , ( 35
)
where f is the depth of the obstacle, a is the radius of obstacle, f (κ) is given by:

f (κ) = b 0 + b 1 (κa) + b 2 2! (κa) 2 + b 3 3! (κa) 3 + • • • , ( 36
)
and b 0 , b 1 , b 2 , • • • is obtained as the solution of the system:

(1 + q 1 γ 2 )b 0 + q 2 γ 3 b 1 + q 3 γ 4 2! b 2 + q 4 γ 5 3! b 3 + • • • = 1 q 2 γ 3 2! b 0 + 1 + q 3 γ 4 2! b 1 + q 4 γ 5 2!2! b 2 + q 5 γ 6 2!3! b 3 + • • • = 0 q 3 γ 4 3! b 0 + q 4 γ 5 3! b 1 + 1 + q 5 γ 6 2!3! b 2 + q 6 γ 7 3!3! b 3 + • • • = 0 • • • = 0,

Figure 2 -

 2 Figure 2 -The waves generated by ship/boat in canal/river

Figure 3 -

 3 Figure3boats as watercrafts in Kapuas river, Indonesia[START_REF] Pontianak | Transportasi motor air melintasi sungai kapuas[END_REF]

Figure I. 1 -

 1 Figure I.1 -Main notations for the geometrical setting of the problem.

  35) = Φ(t, x 1 , 0) + o(|(η, Φ)|) , (I.36) = Φ| S 0 + o(|(η, Φ)|) , (I.37)

Figure I. 2 -

 2 Figure I.2 -Path of integration for z → e -izω z -ν in [50]

67 )

 67 Thus from (I.45) and (I.67) our Green's function of the Neumann-Kelvin problem writes:

Figure I. 4 -

 4 Figure I.4 -Plot of the Green's function, along with its asymptotic approximation, and their difference, for ν = g/U 2 = 0.5 and a source at (0, -2).

the proof of the theorem. Remark 3 .

 3 If α moreover solves (I.79), we have Q = 0, and hence we recover the conditions at infinity (I.74)-(I.75) for Φ defined as in (I.76).

  119) Let us denote t a,b , b a,b , l a,b , r a,b respectively the top, bottom, left and right parts of ∂R a,b . We denote by T , B, L and R the corresponding integrals, so that

Figure II. 1 -

 1 Figure II.1 -Illustration of the displacement field θ and the placement function T .

Figure II. 2 -

 2 Figure II.2 -On the left: for a small enough deformation, regardless of the direction of this deformation, any interior point x belongs to both Ω and (Id + θ)Ω, since Ω is open. On the right: when x ∈ Γ, then no matter how small the (non-zero) displacement field is, it is impossible to ensure x ∈ (Id + θ)Γ, for all directions of deformation.

a := 1 whileAlgorithm 2 :

 12 |κ(Γ) -v| > do Normal displacement: Γ := (Id + a n) Γ Correction of the step: a := a -δ(κ(Γ) -v)end Algorithm for the projection on the constant measure constraint

  Considering a smooth and non-stationary explicit representation γ of the curve, we clearly have φ(γ(s)) = 0 for all s ∈ I, hence: d ds φ(γ(s)) = γ (s) • ∇φ(γ(s)) = 0 , (III.4)

Figure III. 1 -

 1 Figure III.1 -Explicit and implicit representation of curves, from the point of view of numerical approximations. Top: different parametric curves given by the same set of nodes, using different interpolations (linear, spline and piecewise cubic). Bottom: different implicit curves given by the same set of discrete values for the level-set function (in red), using two different interpolations (bilinear and cubic).

Figure III. 2 -

 2 Figure III.2 -Data required for the explicit and implicit representations for a circle. In black, we represented the points of a discrete parametric curve, the mesh of the level-set function for a similar accuracy is represented in blue. In red are represented the "useful" mesh elements.

Figure III. 3 -

 3 Figure III.3 -Plotting of the two exemples of the signed distance function given in (III.7) and (III.8). In red, we represent the level set 0.

Figure III. 4 -

 4 Figure III.4 -Depiction of x and its projection Π Γ (x)

Figure III. 5 -

 5 Figure III.5 -Illustration of the redistancing and extension procedures (iterations 1, 100, 250). First row: evolution of the signed distance function according to (III.47)-(III.48). Second row: evolution of the scalar function a according to (III.47)-(III.49). In red: zero level-set of d 0 . Domain Ω =] -2, 2[×] -2, 2[, grid resolution: 128 × 128.

nFigure III. 7 -

 7 Figure III.6 -Normal velocity plot from (III.52), in red we represented the zero level of θ n .

Figure III. 8 -

 8 Figure III.8 -Parallel curves Γ η and the projection mapping

Figure III. 9 -

 9 Figure III.9 -Depiction of mesh refinements, and the tubular neighbourhoods (in red) for three approaches examined here. The columns represent three different resolutions : 30 × 30, 60 × 60 and 120 × 120. The rows correspond to η * = 0.2, η * = 2h|∇d| 1 and η * = 2h.

  where Ω + is the exterior domain which boundary is Γ which is smooth, say Γ ∈ C 2 . An additional condition is mandatory for solvability: Γ g ds = 0 . (III.89) As in chapter I section 3, equation (I.76), we introduce a boundary integral representation of the solution: u(x) = Γ α(y) G(x, y) ds y , for x ∈ Ω = Ω + ∪ Ω -, (III.90) where α is the single layer boundary potential and E is the Green's function for the Laplace equation on R 2 : G(x, y) = 1 2π log |x -y| . (III.91)

Figure III. 11 -

 11 Figure III.11 -Comparison between the exact and the computed solution for our test-case. The rows represent different values of n (respectively 1, 3, 5), the first column is the numerical solution, the second column is the exact solution, and the third is the log of their differences.

  .6) where G m,k and (∂ n G) m,k are evaluations of the Green's function associated with the Neumann-Kelvin problem (see section 2 of chapter I) and its derivative at the points (x * m , x * k ), that we regularise for |x * m -x * k | < as in chapter III section 3.

3 Figure IV. 1 -

 31 Figure IV.1 -Comparisons of the wave-making resistance vs Froude number profiles. Top left: exact value from T.H. Havelock's solution; top right: computed value with (IV.12) and (IV.13); bottom: relative error.

  ). The figure IV.2 shows the velocity potential obtained with (IV.[START_REF] Dambrine | Regularity of optimal ship forms based on michell's wave resistance[END_REF]) on a cylinder of depth 0.7 and radius 0.5, for four different values of the Froude number: 0.6, 0.8, 1, 1.2. We also represent the streamlines obtained from this potential. The figure IV.3 shows the plot of the free-surface deformation obtained with (IV.20), for the same cylinder and the same values of the Froude number. The query points chosen for Φ, are defined by a 200 × 50 Cartesian grid of [-6, 2] × [-2, 0]

Figure IV. 2 -

 2 Figure IV.2 -Plots of the flow around a circular cylinder, for four different values of the Froude number (from top to bottom 0.6, 0.8, 1, 1.2). The velocity potential was obtained with (IV.19). The colours indicate the value of the velocity potential Φ, the black lines are streamlines, and the red line is the boundary of the obstacle.

-

Figure IV. 3 -

 3 Figure IV.3 -Plots of the free surface, for four different values of the Froude number (from top to bottom 0.6, 0.8, 1, 1.2). The value of the free-surface elevation was obtained with (IV.20).

Figure IV. 4 -

 4 Figure IV.4 -Iterative normal motion of a curve using (IV.33) as a normal displacement functon. In magenta we represent the "safe zone" for p = 2, in red are the elements whose center belongs to the tubular neighbourhood of Γ.

a := 1 while

 1 |κ(Γ) -v| > do Normal displacement: Γ := (Id + a n) Γ Correction of the step: a := a -δ(κ(Γ) -v) end

Algorithm 3 :

 3 [START_REF] Havelock | The collected papers of Sir Thomas Havelock[END_REF] where erf denotes the error function. With all that, we obtain a new projection algorithm in which Γ is implicitly represented by d:Data: A signed distance function d,a tolerance , a target surface area v and a step δ. Result: The signed distance function d = d + a, where a is such that μ(d) = v. a := 1 while |μ(d) -v| > do Normal displacement: d := d + a Correction of the step: a := a -δ (μ(d) -v) end Projection algorithm for the measure constraint: level-set version

Figure IV. 6 -

 6 Figure IV.6 -Plots of the obstacles and corresponding free surface deformations. The blue curves represents the initial situation, the red curves represent the wake and obstacle at the end of the shape optimisation process. The horizontal dashed line represents the zero water level. The depth is fixed at 0.7; from top to bottom, we have: Fr = 0.4222, 0.6667, 0.9111, 1.0333, 1.2778

Figure IV. 7 -

 7 Figure IV.7 -Plots of the obstacles and corresponding free surface deformations. The blue curves represents the initial situation, the red curves represent the wake and obstacle at the end of the shape optimisation process. The horizontal dashed line represents the zero water level. The Froude number is fixed at 0.6667; from top to bottom, we have: f = 0.7, 0.8444, 0.9889,1.2778

Figure IV. 9 - 2 ∞

 92 Figure IV.9 -Plot of the lengths of the optimal obstacles shown in figure IV.8 as a function of U 2 ∞ /g.

Figure IV. 10 -

 10 Figure IV.[START_REF] Brard | Le problème de Neumann-Kelvin[END_REF] -Example of a flow given by a source/sink distribution. In black we show the streamlines, in red, a particular streamline that outlines the shape of the Rankine body, and the source/sink positions are represented in blue and magenta.

Figure IV. 11 -

 11 Figure IV.11 -Example of a flow given by a source/sink distribution in the Neumann-Kelvin model. In black we show the streamlines, and the source/sink positions are represented in blue and magenta.

Figure IV. 12 -

 12 Figure IV.12 -Wave making resistance of a circular cylinder as a function of the depth and the parameter U 2 ∞ /g. The white line joins the maxima of R w with respect to f , for different values of U 2 ∞ /g.

3 Figure

 3 Figure IV.13 -Plots of the R w vs Fr profiles, for optimal shapes with different values of the "target"

Fr 0 .

 0 4222 0.6667 0.9111 1.0333 1.2778 V max /V min 1.2535 1.2716 1.2270 1.1948 1.1650

Figure

  Figure IV.14plot of the flow for the optimal shapes previously examined in figure IV.6. The colours represent the speed. The colour axis are purposely chosen between 0 and the maximum speed at the surface.

Figure A. 1 -

 1 Figure A.1 -Path of integration for z → e -izω z -ν in equation (13)

Figure A. 2 -

 2 Figure A.2 -Path of integration P - R for g (clockwise orientation).

  z)(z + i |b|)] . (25)A simple, yet tedious calculation gives:lim z→-i|b| d dz [g(z)(z + i |b|)] = i k 2 e -k|b| .(26)

  1 e -|k||y 2 | .

				(I.53)
	Hence, since in our case y 2 < 0, we have, from (I.53) and (I.51):
	-k 2	Ĥy (k, 0) + ν∂ 2 Ĥy (k, 0) = -|k|e -iky 1 e |k|y 2 .	(I.54)
	Further, we get the following value of A by using equation (I.50) and (I.54):
		A =	e -(iky 1 )+|k|y 2 |k| -ν	(I.55)
	Further, we substitute (I.55) into equation (I.50):

Table III .

 III 1 -Error and local order of convergence in the exterior Neumann problem for n = 1.

		η * = 2h |∇d| 1		η * = 4h		η * = 0.01
	N	E (10 -3 ×)	order E (10 -3 ×) order E (10 -4 ×) order
	256	0.2917	-	0.5256	-	0.7769	-
	512	0.0806	1.8551	0.1263	2.0573	0.3239	1.2624
	1024	0.0211	1.9352	0.0323	1.9651	0.1268	1.353
	2048	0.0055	1.9397	0.0079	2.0322	0.0508	1.3203

Table III .

 III 2 -Error and local order of convergence in the exterior Neumann problem for n = 3.

		η * = 2h |∇d| 1		η * = 4h		η * = 0.01
	N	E (10 -3 ×)	order E (10 -3 ×) order E (10 -4 ×) order
	256	0.3624	-	0.4995	-	0.7449	-
	512	0.0919	1.9794	0.1232	2.0199	0.3192	1.2226
	1024	0.0232	1.9859	0.0319	1.9468	0.1256	1.3456
	2048	0.0059	1.9753	0.0079	2.0233	0.0505	1.3145

Table III .

 III 3 -Error and local order of convergence in the exterior Neumann problem for n = 5.

		η * = 2h |∇d| 1		η * = 4h		η * = 0.01
	N	E (10 -3 ×)	order E (10 -3 ×) order E (10 -4 ×) order
	256	0.3448	-	0.4739	-	0.7166	-
	512	0.0896	1.9436	0.1201	1.9809	0.3143	1.1892
	1024	0.0229	1.9667	0.0316	1.928	0.1244	1.3365
	2048	0.0058	1.9715	0.0078	2.0143	0.0502	1.3102

  .1 show a reasonable rate of convergence of the method, with a fairly good overall accuracy of around 2% in the worst case.

		Fr = 0.3	Fr = 0.9	Fr = 1.3
	N	Error Order Error Order Error Order
	32 0.0214	-	0.0176	-	0.0164	-
	64 0.0070 1.6042 0.0071 1.3086 0.0066 1.3048
	128 0.0016 2.1056 0.0004 4.0390 0.0002 5.3313
	256 0.0004 2.0862 0.0002 1.2454 0.0001 1.4025

Table IV .

 IV 1 -Convergence results for our numerical method on the wave-making resistance on the case of a circular cylinder of depth 0.7m and radius 0.5m, for three different values of the Froude number.

Table IV . 2 -

 IV2 Table showing the reduction of the wave-making resistance obtained with the shape optimisation algorithm, for f = 0.7, and for 10 evenly spaced values of the Froude number between 0.3 and 1.4.

		0.3000	0.4222	0.5444	0.6667	0.7889
	R w ratio	0.0870	1.9188e-06 8.8098e-07 5.6949e-06 3.9373e-05
	Fr	0.9111	1.0333	1.1556	1.2778	1.4000
	R w ratio 8.1079e-05 6.8562e-05 9.1518e-04 8.7941e-04 9.2692e-05

Table IV .

 IV [START_REF] Bai | A variational method in potential flow with a free surface[END_REF] speed ratio at the free surface in the cases presented in figure IV.14.

International Maritime Organisation, which is the United Nations regulatory agency for the maritime sector.

this does not happen in general, and is the reason behind the lack of smooth solution for the eikonal equation (even for smooth domains), as we will see later.

Weighted Essentially Non-Oscillatory discretisation of order 5 accuracy.

The method we used for the evolution of the explicit curve is rather naive, and for the sake of fairness we mention that some numerical methods exist to prevent this issue. For instance the ALE 3 method (see[START_REF] Donea | Arbitrary Lagrangian-Eulerian Methods[END_REF] for a review on these methods) applies a tangential motion to the points in such a way that the spacing between the control points of the curve[START_REF] Bai | A variational method in potential flow with a free surface[END_REF] Arbitrary Lagrangian-Eulerian method

Acknowledgement

Data:

The initial shape Γ 0 , a step δ r , a velocity U ∞ , a tolerance Tol.

Result:

The final shape Γ f . err := 1 R 0 := R w (Γ 0 ) Γ := Γ 0 n := 1 while err > Tol do Get α0 by solving:

Get q 0 by solving D * 1 [q 0 ] = -2 Γ Re(E(•)E(y))α 0 (y) ds y Get ∇ Γ R w (Γ) from (II.128), knowing α0 , and q 0 Deform the domain in the gradient direction: Γ := (Id -δ r ∇R W (Γ))Γ

Chapter IV

Numerical implementation and results

In this final chapter, all the mathematical and numerical methods presented previously are combined for the purpose of implementing the shape optimisation algorithm 1 presented in chapter II. We start by introducing some notation for the discretisation of the problem. Then we present a validation of the numerical method used to solve the boundary integral formulation of the Neumann-Kelvin problem in a particular case where the exact solution is known. Various technical issues related to the actual implementation of the method are presented along with a full description of the algorithm. Finally we show some results obtained with this algorithm and discuss their relevance for practical applications.

General notation

We start by recalling the notation of chapter III for the discretisation of the problem. Let R be a rectangular domain containing the boundary Γ of the obstacle. Let (x i,j ) i=1..N 1 j=1..N 2 be a Cartesian mesh of R, with spacings h 1 and h 2 (sometimes referred as h when h 1 = h 2 ). The approximation of a function f : R → R M (for M ∈ R) at point x i,j will be denoted f i,j .

Let

where d is the signed distance function of Γ, and η * is the thickness of the tubular neighbourhood considered. Typically, as explained in chapter III, we will take η * = 2h. We note that for the Neumann-Kelvin problem we did not observe a big difference between the choices η * = 2h and η * = 2h|∇d| 1 . Let (m k ) k=1..N N be the sequence of couples of indexes of N (for some ordering). For the sake of simplicity we will denote f k for f m k (no ambiguity remains here because only one index appears).

Recall from chapter I the boundary integral operators involved in the wave-making re-

Implementation of the methods

In this section we give a presentation of the shape optimisation algorithm we designed from the methods and techniques developed in all the previous chapters. Recall the expression of the shape gradient from chapter II:

where α and q solve respectively:

We also recall that the principle of the shape gradient descent method is to change, step-bystep the boundary with:

where δr is a step parameter that has to be chosen small enough to ensure the stability of this process, and the deformation "direction" θ is chosen opposite to the shape gradient: θ = -w n, where:

When d n is the (supposedly known) signed distance function of Γ n , a level-set function dn+1 of Γ n+1 can be obtained with:

where w * is the constant-along-normals extension of w around Γ k . Once dn+1 is known with (IV.26), it is possible to obtain a signed distance function by solving the eikonal equation |∇d n+1 | = 1 with the boundary condition d n+1 | ( dn ) -1 (0) = 0, as explained in the section 1.4 of chapter III.

The full algorithm

The figure IV.5 presents a complete depiction of the implementation of the shape optimisation algorithm 1, using the level-set representation of the obstacle boundary Γ. Our algorithm uses a stopping criterion Crit that compares the current signed distance function d with the one obtained at the previous iteration (denoted d old ). The idea is to compare the value of the wave-making resistance computed with d and d old . If their difference is smaller than a tolerance value TOL, then the algorithm is stopped. The value of TOL is scaled on the initial wave-making resistance. To be complete, we specify that this criterion is in fact averaged on several steps (in our case 1000 steps) to prevent the algorithm to stop just because two consecutive steps gave values for the wave-making resistance that are too close.

The input parameters for this algorithm are the grid step sizes h 1 and h 2 , the initial shape Γ 0 , the velocity U ∞ , the deformation step δr. The initial shape will be a circle of center (-f, 0) and radius a. The parameter of importance here is the ratio f /a. The reason behind this is that the Neuman-Kelvin model is invariant with respect to the Froude number:

As an illustration, let us consider a situation in which both f and a are multiplied by 2 (both are multiplied in order to keep the ratio constant), then the solution of the Neumann-Kelvin system obtained with a velocity of U ∞ / √ 2 will be the same (up to a multiplicative constant).

Hence, since the relevant parameter is f /a, in all the following simulations, we will keep a = 0.5m, and change the depth f and the upstream velocity U ∞ . The initial grid resolution is 128 × 128, and is motivated by a balance between the computational cost and the accuracy (for an idea of the accuracy, we refer to the table IV.1). Finally, the step δr is changed at each iteration, and we choose δr = 10 -2 h/ max(w), which means a displacement of at most 10 -2 h per iteration. The constant 10 -2 that we used here was obtained empirically.

General conclusion

In this study, we computed numerically several 2D shapes which minimise the wave-making resistance of a moving obstacle fully submerged in a fluid with a free surface.

We considered an object immersed in a homogeneous, inviscid and incompressible fluid with a free surface and which moves at a constant velocity. The flow around the object was described by the Neuman-Kelvin (NK) equations, a well-known problem [START_REF] Kuznetsov | Linear Water Waves[END_REF] which can be derived from the irrotational Euler equations with a free surface by a linearisation technique. We used a single layer formulation to represent the NK-problem; in this approach, the fundamental solution handles the Laplace equation and the linearised free surface boundary condition. Using an energy method, we obtained a formula which gives the wave-making resistance as an integral on the boundary of the domain. This formula allowed an accurate computation of the wave-making resistance. Moreover, it showed that the value of this resistance is proportional to the square of the amplitude of the waves generated by the obstacle downstream.

In order to find a local minimum of the wave resistance (for a shape with a given area, in order to avoid degeneracy), we used a gradient descent method. The gradient with respect to the shape is obtained by computing the derivative of the Lagrangian of the problem in a "pulled-back" form through the boundary variation method.

For the computation of the boundary integrals which appear in the formulation of the NK-problem and in the wave-making resistance, we implemented a level-set method. Such a method is very interesting when dealing with moving boundaries as we do in our shape optimisation method. Level-set methods can be greedy, but we managed to reduce the computational cost by working on a tubular neighbourhood of the boundary of the object.

The sequence of shapes given by our shape gradient descent method converged for sev-work with the full Navier-Stokes equations instead of the Neuman-Kelvin problem. For instance, we could integrate a shape descent method to the 2D level-set code developed by James (see for instance [START_REF] James | Cut-cell method: application to water waves generated by a submerged obstacle[END_REF]).

Appendix

A The Rayleigh dissipation and complex integration

The goal of this appendix is to prove that considering

We split Hy into two parts, Hy = Hy,1 + Hy,2

where Hy,1 solves

and Hy,2 solves

Hence, for k > 0, we have:

For k < 0, we proceed in a similar manner, by integration on P + R , which is the symmetric of P - R with respect to the real axis. Again, using the theorem of residues, we have:

The same way as before, we have:

Hence, for k < 0, we have: Finally, we have:

when we put

with li is the logaritmic integral. The surface elevation is given also by Havelock in [START_REF] Havelock | The forces on a circular cylinder submerged in a uniform stream[END_REF] for x far from the source:

with f * (κ 0 ) denotes a complex conjugate of f (κ 0 ).

Résumé

Abstract

In this thesis, we compute the shape of a fully immersed object with a given area which minimises the wave resistance. The smooth body moves at a constant speed under the free surface of a fluid which is assumed to be inviscid and incompressible. The wave resistance is the drag, i.e. the horizontal component of the force exerted by the fluid on the obstacle. We work with the 2D Neumann-Kelvin equations, which are obtained by linearising the irrotational Euler equations with a free surface. The Neumann-Kelvin problem is formulated as a boundary integral equation based on a fundamental solution which handles the linearised free surface condition. We use a gradient descent method to find a local minimiser of the wave resistance problem. A gradient with respect to the shape is calculated by a boundary variation method. We use a level-set approach to calculate the wave-making resistance and to deal with the displacements of the boundary of the obstacle. We obtain a great variety of optimal shapes depending on the depth of the object and its velocity.