Optimization and response surfaces Hydrofoil optimization can be performed in different ways. Hydrofoils can be either optimized independently from the rest of the ship they will be fitted to or, in a more robust way, they can be optimized by taking into account the whole ship and simulating its behavior in a third party simulating code (a Velocity 2 Chapter 1 Introduction

Introduction

Shape optimization is the search for designs which are best suited for their intended role. In ship design, shape optimization has been performed from the beginning; through experience, shipbuilders over the generations developed types of ships that were ideal for the work and the local conditions which they encountered, given the construction methods that were available. However, it was the towing tank which first made it possible to compare different designs and select the best one, before an actual ship is built. Today, through simulation-based design, this selection of design candidates is more and more performed with computer simulation. And as computer power increases, naval architects start turning to simulationbased optimization, where not only the simulations, but also the choice of the design candidates and the geometry generation is handled in a connected, fully automatic way.

This thesis presents a simulation-based optimization procedure for ships fitted with lifting hydrofoils. It shows how all parts of the procedure have been developed to meet the demands of efficiency and flexibility which this application demands.

Design of ships and hydrofoils

The objective of most ship optimizations is to minimize the power needed to attain a certain velocity, by reducing the hydrodynamic drag. Additional objectives are to ensure the stability of the vessel, to limit its movements in waves, or to provide sufficient capacity for maneuvers. Since the costs associated with running a ship are large, even small increases of the ship performance can lead to significant gains. However, contrary to for example cars or aircraft, most ships are one-off designs or are built in series of two or three units maximum. Thus, the development budget available for design optimization is limited. Some of the world's fastest ships are fitted with hydrofoils. These are lifting surfaces located below the water surface, which partially or completely lift the ship hull out of the water in order to reduce its drag at high speeds. Motor-driven hydrofoil ships have been developed during the second half of the last century, to provide fast military patrol vessels and passenger transport.

Recently, the interest for hydrofoils has rekindled with the advance of hydrofoil-borne sailing yachts. The 34 th edition of the America's Cup, held in San Francisco in September 2013, saw 72-foot racing catamarans 'flying' on foils to attain speeds well in excess of what traditional racing yachts can obtain. It was this race which convinced the general sailing public of the feasibility of foil-borne sailing and which subsequently led to the introduction of hydrofoils for various smaller yachts.

In terms of optimization, racing yachts are similar to cargo ships since on one hand, far-going optimization is required since minute differences in performance can make the difference between winning and losing a race. On the other hand, the budgets available for the design of a sailing boat are small. In comparison, while America's Cup budgets are at least ten times more than those for any other sailing competition, the total cost of an America's Cup campaign is less than the price of a single large airliner (of which hundreds are built). Thus, yacht shape optimization must be both efficient and cost-effective.

Surrogate-based optimization

For simulation-based optimization, the geometries are expressed in terms of design parameters. In the design space consisting of all the possible combinations of these parameters, the optimization then involves the search for the parameter combination which gives the best performance. When the design space is large or when the ship behavior is complex, many evaluations of parameter combinations are required. However, each evaluation requires a complete numerical simulation and the total number of simulations that can be run is usually limited by time or financial costs.

One way to address this limitation is to build a simplified model of the ship behavior from simulations in a limited number of design points. Meta-modeling, also called surrogate modeling, is the process of generating such models of models or metamodels. The optimization is then performed over the metamodel, instead of using the real simulations. Since evaluations of a surrogate model are cheap, the computational cost for the optimization is much reduced.

Building a surrogate requires evaluating the original model at specified points and gathering the corresponding responses. The amount of sample points needed to approximate the behavior of a numerical model depends on the complexity of the response, which is not necessarily known beforehand. The chosen sampling strategy plays an important role in the performance of the surrogate model: under-sampling might not allow capturing the complexity of the phenomena and over-sampling leads to prohibitive computation times.

Adaptive sampling uses the responses of the simulations for adjusting the sampling in a sequential way while the surrogate is being constructed. Before each simulation, the design point is searched where a new simulation would be most beneficial for the metamodel; this point is simulated and added to the metamodel. Thus, the required number of simulations is minimized for the efficient construction of the surrogate.

Prediction Program (VPP) for instance). Designers usually switch from one approach to another according to the design stage.

Surrogate models can be used for both approaches but, if built in adaptive way, will be sampled quite differently. VPPs are based on separate metamodels for each component of a ship, which are combined through the VPP. Thus, metamodels are required which are accurate throughout the input space. To create such response surfaces, sampling from areas of high uncertainty will be favored, in order to ensure their reliability everywhere.

Automatic geometry optimization on the other hand requires sampling from areas likely to offer improvement over the current best design. This means that, once the global behavior of the metamodel is known, it can be beneficial to concentrate the sampling area close to the optimum.

The optimization chain Metamodel-based optimization requires the automation of three processes:

• The optimization over the design space, which includes the metamodel construction and the adaptive sampler that decides which design points will be simulated.

• The geometric modeler which translates a set of design parameters into a threedimensional geometry, suitable for simulation or for manufacturing.

• The numerical fluid simulation which evaluates the flow around candidate geometries and provides the responses on which the metamodel is based.

These three items are called in a loop which runs until the quality of the metamodel is deemed sufficient. Their coupling is performed through a book-keeping script which stores the results of each step and provides relevant data to subsequent processes.

Contents of the thesis As indicated, the objective of this thesis is to develop a robust and efficient adaptive metamodel-based optimization framework for the design of lifting hydrofoils. All the steps in the optimization chain are addressed.

First, chapter 2 introduces the lifting hydrofoil and describes different foil configurations. A short historical overview traces the roots of the hydrofoil shapes in use today. The chapter ends with a study of hydrofoil flow physics, which presents the hydrofoil behavior that is important for shape optimization and numerical simulation.

Chapter 3 then addresses the numerical simulation of the flow around hydrofoils. After an introduction of flow simulation in marine applications, the chapter presents the two-phase Navier-Stokes solver ISIS-CFD that is used in this thesis. The governing equations are introduced, as well as the specific aspects of the discretization which deal with free-surface flow. The last part of the chapter describes techniques to modify the computational mesh, such as grid deformation and adaptive refinement.

Constructing a metamodel requires interpolating in the design space between the sampled points. The metamodels in this thesis are created using Gaussian Process regression (chapter 4), which treats the sampled results as if they are a realization of a stochastic process. The prediction is then based on the most likely outcome in each point, given that the result in the sampled points is known. Different Gaussian processes are described and tested to find the one which performs best for the applications considered here.

In chapter 5, the adaptive sampling is discussed. The selection of a point to sample is formulated as a maximization problem, over an acquisition function which combines the metamodel values and an estimation of the metamodel uncertainty. A custom acquisition function is presented with an uncertainty estimate based on the variance of the Gaussian process, weighted with a cross-validation error estimation in the sampled points. Tests show that this custom function is more reliable than existing acquisition functions, since it gives the same performance for simple test cases, while retaining its good performance in more difficult situations where existing approaches fail.

Chapter 6 presents the development of a geometric modeler which creates lifting hydrofoils for sailing yachts. The modeler is based on a parametrized central line, along which an airfoil profile with varying chord, twist and shape is swept. Within the scope of this thesis, the geometric modeler has been applied in industrial and educational projects. As an illustration of the modeler's capacities, two of these applications are presented.

The elements described in all the previous chapters are brought together in chapter 7, which presents the surrogate-based shape optimization of lifting hydrofoils based on Navier-Stokes simulation. First, the tests of chapter 5 are repeated for response surfaces and optimization of two-dimensional hydrofoils. Then, optimal two-dimensional foil shapes are studied and the dependence of the optimum on the operating conditions is discussed. Finally, a three-dimensional hydrofoil created with the geometric modeler is optimized and the optimal shape is analyzed in terms of stability and cavitation risk.

Chapter 8 concludes this thesis with an analysis of the further developments which are concievable.

Hydrofoils and foilers

This chapter is dedicated to the presentation of hydrofoils and foiler vessels. Foiler vessels are ships fitted with hydrofoils. After a definition of the term hydrofoil (section 2.1), some applications of hydrofoils are presented in section 2.3 through relevant examples of foiler vessels, both powered and sailing. The last part of this chapter (section 2.4) is dedicated to the flow physics of hydrofoils.

This chapter is mainly bibliographic and is used as a basis for the choices of simulation tools, optimization processes and geometrical modeling methods, as emphasized in the conclusion of the chapter.

Hydrofoils

The term foil refers to a solid object with a shape such that, when placed in a fluid with a relative velocity and at a suitable angle, it generates a force that is not aligned with the flow (figure 2.1). This force can be decomposed in two components:

• Lift, is the force generated perpendicular to the mean fluid flow.

• Drag, is the force generated parallel the mean fluid flow. If the surrounding fluid is air or more generally a gas, the foil is called an airfoil. In the context of this thesis, foils evolving in water are considered, in which case the foil is referred to as a hydrofoil. Foils can be either used to transfer energy to the surrounding fluid and to set this fluid into motion (pumps) or to create forces and moments from a moving fluid for control and/or sustention purposes. The shape of foil is designed to generate a lift greater than its drag.

Lift

Hydrofoils in the marine context

In the marine domain, hydrofoils are lift-generating surfaces operating below the free surface. Marine vehicles, unlike their aerial counterparts, are mainly using lift-generating surfaces for control purposes. By creating forces and moments, lifting surfaces are used for steering ships (yaw control); in that case the foils are called rudders. Stabilizing fins are a type of hydrofoil used for roll control whereas pitch damping fins are used to control pitch movements. Control surfaces for submarines are called hydroplanes. A more extensive presentation of marine controle surfaces can be found in [START_REF] Molland | Marine rudders and control surfaces: principles, data, design and applications[END_REF].

This thesis will focus on a special type of hydrofoils whose role is to create a vertical lift able to support the weight of a ship in order to lift its hull above the free surface. The goal is to significantly reduce the ship drag by decreasing the wetted area of the ship, but also by decreasing the ship's wave added resistance. Seakeeping capabilities can also be improved as the ship's response to waves is reduced.

This section will present how hydrofoils are used in the marine context, with the different possible design options, the consequences in terms of flow physics and the parameters that need to be considered in a design and optimization process.

Hydrofoil lift fraction and longitudinal distribution of lift

Lift fraction Advanced marine vehicles and their hybrid derivations can be understood through the use of the sustention triangle (figure 2.2). This graphical method is used to classify marine platform concepts according to the type of lift used to sustain the water craft. This method is presented in [START_REF] Jewell | Hybrid fluid-borne vehicles[END_REF] and [START_REF] Jewell | Possible naval vehicles[END_REF]. Each vertex of the triangle corresponds to a particular way to lift a ship: unpowered static lift (buoyancy) for the top vertex, dynamic lift for the right lower vertex and powered static lift for the left lower vertex. Hovercraft are examples of vehicles making use of powered static lift. The lift is generated by a power system and is active even if the vehicle is stationary.

Hydrofoils can be designed to fully or partially lift the hull out of the water. When the set of hydrofoils only generates a fraction of the force needed to lift the ship, the rest is being generated by buoyancy and/or dynamic lift issued from the hull. This case is called a hydrofoil assisted ship.

The percentage of lift generated by the set of hydrofoils relative to the displacement of the ship is referred as the hydrofoil lift fraction. The lift fraction is equal to 100% in the case of a fully flying ship. On the sustention triangle, this category of watercraft are located on the right lower vertex (figure 2.2). A partial lift fraction or semi-flying mode is often referred to as a skimming mode and still allows for interesting drag reduction. This kind of displacement mode corresponds to the right edge of the sustention triangle.

In this thesis, only fully flying vessels are considered. As a consequence, interactions between hulls and hydrofoils will not be taken into account.

Longitudinal distribution of lifting surfaces

The lifting surface of the hydrofoils must be longitudinally fractioned for flight stability. The distribution of foil areas along the longitudinal axis of the hull represents a major classification criterion for hydrofoil vessels (Johnston, 1985).

Three main categories can be distinguished. In the classic configuration, also referred to as plane configuration, (figure 2.3a) the main part of the ship weight is supported by the 2.2 Hydrofoils in the marine context forward foil. The canard configuration (figure 2.3c) is the opposite of the conventional configuration and the main part of the ship's weight is supported by the aft hydrofoil. In the tandem configuration (figure 2.3b) both foils have the same area and support an equal part of the weight.

In this thesis multiple element hydrofoils will not be considered, the foils are only studied individually. However, the longitudinal distribution can induce different loadings on the fore and aft foil and as a consequence, the foils operate at different lift coefficients. This has to be taken into account for the design and optimization of each foil.

Fully submerged and surface piercing hydrofoils

Two main categories of hydrofoils can be distinguished based on how they interact with the free surface. They have either a fully submerged or a surface piercing foil configuration.

Fully submerged: the lifting surface is fully below the free surface, the lift is therefore relatively unaffected by the free surface. The ship requires a control system to maintain flying height and attitude.

Surface piercing: the lifting surface intercepts the free surface. The lift varies in relation with the foil submergence. Ships equipped with this type of hydrofoil can be made intrinsically stable in pitch, heave and roll so they do not require active ride control. This distinction between fully submerged and surface piercing hydrofoils is important as it implies different flows and hydrodynamic performances. Surface piercing hydrofoils strongly interact with the free surface inducing phenomena such as spray or ventilation (see section 2.4.3).

Structural, possibly streamlined elements can be connected to the lifting parts to support them. They are referred to as struts. For example, the fully submerged hydrofoil in figure 2.4b is supported by a vertical strut.

8

Chapter 2 Hydrofoils and foilers

Hydrodynamic control effectors

A vehicle interacting with a fluid, such as a plane, a submarine or a foiler vessel, can be dynamically controlled by using lifting surfaces and adjusting the force they generate. Those lifting surfaces used for control are called control effectors.

The aero/hydrodynamic force of lifting surfaces can be adjusted in various ways [START_REF] Chudoba | Stability and control of conventional and unconventional aircraft configurations: a generic approach[END_REF]. The most common hydrodynamic control effectors fitted on foiler vessels are using incidence control (figure 2.5a) and camber control through the use of a trailing edge flap (figure 2.5b). When deflected, a trailing edge flap usually results in a profile shape which is not optimal. The camber is not continuous and the hinge used to connect the flap represents a source of drag. Incidence control effectors offer a higher hydrodynamic efficiency [START_REF] Chudoba | Stability and control of conventional and unconventional aircraft configurations: a generic approach[END_REF] but require more powerful actuators. Incidence control can also be obtained through the movement of both the lifting part of the foil and the supporting strut. In that last case the actuators need to be even more powerful.

In this thesis, camber control will not be considered. Incidence control will be used to match either a required lift coefficient or a force, in the case of a constrained optimization.

Hydrofoils in use

This section illustrates hydrofoils in use on foiler vessels that are either powered or using sails, with different examples of foil layouts and types.

Motor foilers

Hydrofoils have been used on different watercraft since the end of the 19th century. The first known experiments with hydrofoils were performed in 1861 by a British engineer, Sir Thomas Moy, who was trying to understand the physics of lifting surfaces for the development of airplanes. This work led in 1897 to the construction of the first fully foiling vessel by French engineer Charles de Lambert. His boat, a catamaran powered by a steam engine (figure 2.6), was capable of flying at 9 knots of speed. In 1906, the hydrofoil boat built by Italian inventor Enrico Forlanini reached a top speed of 36 knots during testing. Foilers remained experimental craft until the end of the Second World War. Then, knowledge acquired in fluid dynamics and structures thanks to the development of aviation, enabled the creation of new foiler craft with improved performances. Later, in the 1960s, aviation turbines and electronic control systems allowed to design high performance foiler craft with fully submerged hydrofoil configurations and offshore capabilities. At that time, attack submarines in both Soviet and NATO navies were able to reach speeds in excess of 40 knots. Many navies expressed a need for ocean-going patrol craft that could outrun such high-speed submarines. Companies commissioned by the US government chose to investigate the fully submerged hydrofoil configuration. Stabilization systems were developed, relying on gyroscopes, accelerometers and sonar or radar height sensors to operate control effectors in a way similar to modern airliners or jet fighters.

Foiler vessels designed during that period represent the state of start in the domain. Later with the oil crisis most of those projects were abandoned and the ships decommissioned.

More technical and historical details can be found in Chapter 5: Hydrofoil Craft of [START_REF] Yun | High performance marine vessels[END_REF]. 2.7a). These ships, first tested in 1975, were capable of operating at speeds of 40 knots in seas up to 4 meters and could reach a maximum speed of 48 knots. They were fitted with fully submerged foils in a canard configuration. Control effectors are trailing edge flaps on both the fore and aft foils.

Based on the success obtained with the PHM, Boeing designed the Jetfoil 929 (figure 2.7b), a passenger version with a revised planing hull having a wider beam and superstructure, capable of transporting up to 250 passengers. The Boeing Jetfoil 929 and its derivatives built under licence are still operated, mostly in the Hong Kong bay and in Japan.

Surface piercing foiler vessels Even though the most technologically advanced foiler vessels were using fully submerged hydrofoils, some highly effective craft were designed with surface piercing hydrofoil. Among them was the FHE-400 "Bras d'Or", a hydrofoil patrol craft prototype built by De Havilland Aircraft Company for the Canadian navy (figure 2.8a). This 180 tonne vessel, intended as the first unit of a fast coastal patrol class, was powered by a gas turbine and able to run at 50 knots in 3.5 m seas or close to 60 knots on calm water. "Bras d'Or" used a canard layout, which is more conventional with surface piercing hydrofoils than with fully submerged foils.

In figure 2.8b is presented an example of a foiler fitted with surface piercing foils in a classical plane layout. "Aldebaran" is a unit from the RHS160 series designed and built by the Italian shipyard Rodriquez. This foiler was carrying up to 200 passengers and capable of reaching speeds up to 36 knots. The RHS160 serie has the interesting particularity of combining surface piercing hydrofoils (which are intrinsically stable) with active ride control to enhance stability. Foilers of this type still operate in the Mediterranean sea.

Sailing foilers

Sailboats are interacting simultaneously with two fluids, air and water. A sailboat is able to generate its own motion by using the energy from the wind and creating an aerodynamic force, thanks to a sail or a wing. Through its hydrodynamic components (hull and appendages), the sailboat creates a hydrodynamic force which, combined with the aerodynamic force, sets the boat into motion before reaching an equilibrium state and a steady state motion.

Sail force

Hull force Hydro Aero Wind Boat motion Figure 2.9: Schematic representation of the aero-hydro equilibrium of a sailing yacht in steady conditions (after [START_REF] Garrett | The symmetry of sailing: the physics of sailing for yachtsmen[END_REF]).

Figure 2.9 shows a schematic representation of the equilibrium reached by a sailboat in steady motion. Under the influence of the wind, the sail creates an aerodynamic force. Unless the boat is moving in the direction of the wind, this sail force is not aligned with the direction of motion: it has a lateral component which may even be larger than the forward traction. To prevent the boat from drifting sideways, this lateral force must be counterbalanced by the hull or by a hydrofoil. Hydrofoils on sailboats are thus designed to generate a horizontal sideforce to balance the force generated by the sails in addition to producing a vertical force to lift the boat out of the water (figure 2.10). For instance in the 35th America's Cup, AC50 main foils were typically required to generate a vertical force twice the magnitude of the horizontal side force (respectively 30 and 15 kN ).

There is a long history of hydrofoil assisted sailing yachts, for which the hydrofoils generate a vertical force that is not sufficient to fully lift the hull out of the water (as explained in section 2.2.1). However, until recently full flight has been restricted to experimental sailing yachts or high-end record breaking yachts.

The first prototype of a foiling sailing yacht was developed by Alexander Bell before 1920, in parallel to his work on motor foilers. In 1920, Malcolm and Thomas McIntyre were the first to patent the concept of a flying yacht equipped with surface-piercing hydrofoils. In 1954 the American engineer Gordon Baker demonstrated the "Monitor" foiling yacht, originally designed with two rigid wing sails but built with a conventional rig. The "Hydroptère" (figure 2.11), based on Eric Tabarly's 1976 prototype, is an example of a fully foiling high performance sailing yacht which uses surface piercing foils in a conventional layout. A submerged T-foil rudder, located at the transom of the main hull, is used as an elevator to control the pitch motion of the boat. The America's Cup, first disputed in 1851, is a major international sailing competition focused on design innovation. The AC50 yachts (figure 2.12) were developed and raced during the 35th edition in June 2017. These catamarans were allowed to use L-shaped foils that can almost be considered as fully submerged foils, as the part generating the vertical lift is most of the time immersed and weakly interacting with the water surface. This weak interaction with the free surface induces less passive stability but is compensated by a manually controlled, yet highly effective, hydraulic actuating system for the foil position.

Small dinghies also make use of hydrofoils. The Moth IMCA (figure 2.13) is an interesting example. For these, a mechanical control system allows to use a T-foil configuration which Similarly to motor yachts, sailing yachts could increase their performance by using automatic control systems and fully submerged hydrofoil configurations but most class rules prevent the use of such systems. One big difference between sailing and powered foilers is the speed operating range. Sailing foilers are relying on the wind to generate motion and the boat speed is directly related to the wind speed. The consequence is that hydrofoils fitted on sailing yachts have to operate in a wide range of speeds, while motor foilers can be optimized for one cruise speed.

The 35th America's Cup provides an interesting example of the complexity of designing hydrofoils for sailing yachts. During that edition of the Cup, the AC50 were sailing in wind speeds ranging from 5 to 25 knots, which meant in term of boat speed, 15 to almost 50 knots. Teams were allowed 2 sets of foils, helping them to better cover the range.

Physics of hydrofoils

This section is dedicated to the physical flow phenomena that can be encountered when operating hydrofoils. More detailed theoretical explanations can be found in chapter 6 of [START_REF] Faltinsen | Hydrodynamics of high-speed marine vehicles[END_REF].

Lift and drag of a 2D foil

Foils are generally designed so that the lift is substantially larger than the drag (figure 2.1). Figure 2.14 shows typical steady characteristics of lift and drag for a foil with a symmetric profile similar to the one shown in figure 2.15b. This figure shows that the lift grows in an approximately linear manner with an increasing angle of attack. At some point the angle of attack becomes too large and a phenomenon called stall occurs. A stall is caused by boundary layer separation and induces a reduction of the lift (see section 2.4.2). Depending on their usage, foils can have different section shapes. The usual distinction is between asymmetrical (figure 2.15a) and symmetrical (figure 2.15b) foil sections. A foil section can be described by various parameters. Thickness and camber are the parameters which have the most significant effect on the airfoil performance.

Mean camber line

Chord line Camber

Boundary layers

When a viscous fluid flows past a body, a thin layer of the flow close to the surface of the body is affected by the viscous friction. The theory for this boundary layer flow relies on the hypothesis of an adhesion of the fluid to the walls, that is, on the hypothesis of a zero relative velocity between fluid and wall. From a zero relative flow velocity at the wall, the velocity in the boundary layer increases until it reaches the outer flow velocity U , as shown in figure 2.16. The boundary layer is characterized by a velocity gradient ∂u/∂y. The distance required for the flow to nearly reach flow stream velocity (99 %), is defined as the boundary layer thickness, δ(x).

Shear stresses are large in the boundary layer and in the viscous wake formed by fluid within the boundary layer being swept downstream of the body. This shear stress is a major source of drag but it means that the viscosity only matters in the boundary layer and the viscous wake. Boundary layer theory [START_REF] Schlichting | Boundary-layer theory[END_REF] describes these effects of the fluid viscosity and the main characteristics of the viscous flow. Boundary layer separation Because of the viscous effects, the flow velocity is reduced in the boundary layer, resulting in a low flow momentum. The boundary layer flow is therefore more sensitive to a pressure derivative ∂p/∂x, compared to the main flow.

In case of a negative pressure derivative (decreasing pressure following the flow direction ∂p/∂x < 0, the fluid is accelerated and the boundary layer thins and stays attached to the surface. A negative pressure gradient is termed favorable. Similarly, if the pressure is constant in the direction of the flow, ∂p/∂x = 0, the boundary layer remains attached and it only increases in thickness.

In contrast, in case of a positive (adverse) pressure derivative ∂p/∂x > 0, the flow in the boundary layer, with its low momentum, might not be able to overcome the pressure to reach a downstream location. In that case the velocity decreases until the velocity gradient

16

Chapter 2 Hydrofoils and foilers becomes zero (figure 2.17). This is where the separation occurs. Further downstream, the velocity gradient becomes negative with the appearance of a flow reversal.

The position of the separation point depends on the pressure derivative ∂p/∂x along the body surface. It also depends on the flow conditions (laminar or turbulent) in the boundary layer ahead of the separation point. Separation is more likely to occur in a laminar boundary layer as in that case the flow has less momentum near the wall compared to a turbulent boundary layer. Laminar separation will be triggered by weaker adverse gradients compared to turbulent layers meaning that turbulent boundary layers are more resistant to separation than laminar ones.

U Suction face (-)

Pressure face (+)

In cr ea si n g p re ss ur e Adverse pressure gradients can be encountered in divergent flows or on curved walls, such as the low pressure face of a hydrofoil. As the angle of attack increases there is a further increase in the negative pressure generated by the hydrofoil, with a pressure peak located after the leading edge. The dashed line in figure 2.18 represents a pressure profile typically observed on the low pressure face of a foil. This profile has an adverse pressure gradient in the flow direction and can lead to boundary layer separation. This is the stall phenomenon

F nh CL(h/c)/CL(h/c = ∞) h/c = 0.5 h/c = 1.0 h/c = 2.0 h/c = 5.0

Figure 2.19:

Influence of the free surface on the lift coefficient (reproduced from [START_REF] Hough | Froude number effects on two-dimensional hydrofoils[END_REF]).

that leads to a loss of lift and an increase of drag. The stall phenomenon can be observed in the evolution of the lift and drag with the angle of attack, figure 2.14.

Free-surface related effects

The flow around a submerged hydrofoil is strongly influenced by the presence of the free surface. This section presents four major physical effects of the water surface which have an impact on the functioning of the hydrofoil.

Pressure field

The pressure field generated by the hydrofoil is affected by the vicinity of the free surface. The lift coefficient decreases as the hydrofoil gets closer to the free surface, because the air above the surface offers less resistance to the foil than the much heavier water. This is confirmed by figure 2.19, from [START_REF] Hough | Froude number effects on two-dimensional hydrofoils[END_REF], which shows the numerically predicted lift coefficient C l as a function of the submergence Froude number

F n h = U∞ √ gh
, for different values of the non-dimensional immersion h/c, where U ∞ is the inflow velocity, h the immersion depth (figure 2.20), g the gravity and c the chord length. These calculations are based on a linear body boundary condition and linear free-surface condition.

Waves When a hydrofoil evolves in a fluid, it creates a pressure field. In the vicinity of the free surface, the pressure drop located on the suction face is likely to deform the free surface and create a wave field in its wake (figure 2.20). This deformation of the free surface as well as the resulting wave field results in an increase in drag.

Spray Spray is the projection of a water film behind the trailing edge of a foil or strut, that eventually turns into droplets. It results from the direct interaction of the hydrofoil with the free surface and is a source of drag. Surface piercing foils usually have a pressure difference over the part which traverses the surface; this increases the spray generated by the hydrofoil.

18

Chapter 2 Hydrofoils and foilers Ventilation The low pressure generated on the upper surface of the hydrofoil may reach a value lower than the atmospheric pressure, with a resulting suction of air from the surface if the hydrofoil is close enough to the free surface (figure 2.21). This suction of air causes a loss of lift due to the increase in pressure. It also results in an increase of drag. Ventilation can also start in the tip vortex when it gets close enough to the free surface, in the wake behind the hydrofoil. The entrapped air may then travel upstream through the tip vortex core and reach the hydrofoil surface, creating loss of lift and an increase of drag.

Cavitation

Cavitation is a physical phenomenon occurring when a liquid is subjected to a drop in pressure which, if the drop is large enough, may become less than the saturated vapor pressure. In that case, the liquid undergoes a phase change and passes from a liquid state to a gaseous state (figure 2.22a).

For a propeller or hydrofoil operating in water, cavitation leads to the creation of vapor cavities (figure 2.23), resulting in a sudden loss of lift and an increase in drag. Inertial cavitation is the process where a void or bubble in a liquid rapidly collapses, producing a shock wave. This phenomenon can induce significant damage to lifting surfaces. The onset of cavitation can be explained using the Bernoulli equation and thus neglecting the viscous effects. Considering an immersed hydrofoil close to the free surface, the local pressure at a point on the face of the hydrofoil is given by the following equation, where p a is the atmospheric pressure and U the local flow velocity:

Physics of hydrofoils

p = p a + ρgh + ρ 2 U 2 ∞ 1 - U U ∞ 2 . (2.1)
For a given travel speed of the hydrofoil, the risk of cavitation increases with the local speed on the suction face of the hydrofoil. For a given foil section the maximum local speed increases with the angle of attack, meaning that cavitation is more likely to occur when the hydrofoil is operating at a high lift coefficient.

Although not caused by the proximity of the hydrofoil to the free surface, the occurrence of the cavitation phenomenon is favored by this proximity due to the surrounding pressure drop experienced by the hydrofoil as it approaches the surface. Boundary layer separation will also affect cavitation inception.

20

Chapter 2 Hydrofoils and foilers

Conclusion

This chapter gave a wide view of what hydrofoils are and of the potential issues that can be encountered when designing them.

One of the main points of interest for foiler vessels is flight stability. The choice between intrinsically stable or regulated foil configurations has direct consequences for the operating conditions of the hydrofoils. Within this two categories, hydrofoils have different configurations and shapes dictated by propulsion or structural considerations. Sailing foilers face additional design constraints with the need to generate a horizontal sideforce to balance the aerodynamic force originating from the sails.

This variety of shapes and operating conditions has consequences in the three domains of the optimization process introduced in chapter 1. The geometric modeler has to be versatile enough to generate complex shapes capable of approaching an "absolute" optimum.

Unfortunately the more richness one wants to use, the more parameters are needed, which increases the number of computations required to reach an optimum. A trade-off has to be found, taking into consideration the type of optimization algorithm used and the duration of a each simulation point which depends on the simulation tool used.

Apart from the concern of minimizing the required number of simulations, the optimization algorithm itself has to be able to optimize a hydrofoil according to how it will be operated.

The diversity of use can imply that a hydrofoil must be optimized for the narrow range of operation that is characteristic for a motor foiler, or the extended range of a sailing foiler.

In the remainder of this thesis, some limiting hypotheses regarding the optimization are used. Hydrofoils will be considered as isolated and not interacting with other devices (propellers or hydrofoils) of the ship, through their wash. Furthermore, only steady foilborne operating conditions are considered. No transition phases such as take-off or landing are studied.

Hydrofoils generate complex flows, which in return strongly influence their global performance. Simulation tools capable of capturing the physics of the operating hydrofoil are required. Since hydrofoil behavior is influenced by the free surface, tools capable of simulating and predicting the impact of the free surface are required. Boundary layer separation is also of great interest as the performance of the hydrofoil is dramatically influenced when the flow separates. Therefore, accurate modeling of viscous effects is required. The next chapter will deal with the fluid simulation aspect of a hydrofoil optimization and especially with Reynolds-averaged Navier-Stokes simulations.

Conclusion

Computational Fluid Dynamics

This chapter presents the methods and tools used for fluid simulation and the evaluation of hydrofoil performance. The main simulation software is the FINE™/Marine computing suite. The ISIS-CFD flow solver used in FINE™/Marine is presented in this chapter along with the different governing equations it solves. The chapter is provided as background information and may be skipped by a reader familiar with the subject.

Introducing marine CFD

Computational fluid dynamics (CFD) refers to the numerical resolution of the governing equations of fluid flows, with the help of computers. This technique, although initially developed in the field of aeronautics, is also extensively used in the marine industry for ship design and optimization as maritime transportation represents one of the most important components of transportation technology today.

Marine CFD is used in a wide variety of ways and at all design stages: from preliminary design to high-end calculations where accurate determination of the ship's performance is required. The physics of flows around a ship's hull like wave patterns or wake fields can be simulated with CFD, which can also be used for predicting a ship's behavior such as the interaction with waves (seakeeping) and maneuvering. Those simulations and the information collected are used for hull shape improvement but also for the design of appendages (rudders, stabilisers. . . ) or propellers.

In the field of sailing, the use of CFD is usually restricted to racing yacht or high-end superyacht design because of the costs implied by this type of simulation. The reason for this is, that sailboat design requires large numbers of simulations due to the wide range of operating points and sailing conditions that sailing yachts are subjected to.

CFD has allowed a dramatic increase in racing yacht performance over the last two decades.

Recent racing yachts make extensive use of large and flat planing hulls and/or hydrofoils for reducing drag. This quest to lift the boat out of the water results in a strong interaction between the lifting elements and the free surface (section 2.4). This represents a special challenge for CFD codes as they need to accurately simulate phenomena such as ventilation, spray or wave breaking. Foilers running at high speed are also likely to encounter cavitation on their hydrofoils which represents another challenge for CFD codes.

Presentation of the ISIS-CFD flow solver

The ISIS-CFD flow solver is developed by the METHRIC team (Modélisation des Écoulements Turbulents à Haut Reynolds Incompressibles et Couplages) in the LHEEA lab of Centrale Nantes and CNRS and is distributed by NUMECA Int. as a part of the FINE™/Marine computing suite, which is mainly devoted to marine hydrodynamics. ISIS-CFD uses the incompressible unsteady Navier-Stokes equations.

The solver is based on the finite volume method to build the spatial discretization of the transport equations. The unstructured discretization is face-based and the method is generalized to two-dimensional, rotationally-symmetric, or three-dimensional unstructured meshes for which non-overlapping control volumes are bounded by an arbitrary number of constitutive faces. The velocity field is obtained from the momentum conservation equations and the pressure field is extracted from the mass conservation constraint, or continuity equation, transformed into a pressure-equation. In the case of turbulent flows, additional transport equations for modeled variables are formulated similar to the momentum equations and they can be discretized and solved using the same principles. Free-surface flows are simulated by incompressible and non-miscible flow phases, modeled through the use of conservation equations for the volume fraction of each phase.

For the Reynolds-averaged Navier Stokes equations (RANSE), the method features several sophisticated turbulence models: apart from the classical two-equation k-and k-ω models, the anisotropic two-equation Explicit Algebraic Reynolds Stress Model (EARSM), as well as Reynolds Stress Transport Models, are available, see [START_REF] Deng | Comparison of explicit algebraic stress models and second-order turbulence closures for steady flows around ships[END_REF] and [START_REF] Duvigneau | On the role played by turbulence closures in hull shape optimization at model and full scale[END_REF], with or without rotation corrections. All models are available with wall-function or low-Reynolds near wall formulations. Hybrid LES turbulence models based on Detached Eddy Simulation (DES) are also implemented and have been validated on automotive flows with large separations, see [START_REF] Guilmineau | Numerical simulation with a DES approach for automotive flows[END_REF].

Additionally, several cavitation models based on transport equations for the vapor fraction with source terms have been implemented in the code. The models by Sauer, Merkle and Kunz are available and can be combined with the capturing of the water surface. See for example [START_REF] Dauby | Simulations d'écoulements cavitants par résolution numérique des équations de Navier-Stokes en moyenne de Reynolds[END_REF] for a study of these cavitation models.

Governing equations

This section presents the governing equations of the ISIS-CFD solver and shows how they are discretized and solved.

Conservation equations

The flow solver can deal with multi-phase flows and moving grids. In the multi-phase continuum for incompressible flow of viscous fluid under isothermal conditions, using the generalized form of Gauss' theorem, conservation equations for mass, momentum and volume fraction can be written as:

∂ ∂t V ρdV + S ρ( - → U - - → U d ) • - → n dS = 0, (3.1a) ∂ ∂t V ρU i dV + S ρU i ( - → U - - → U d ) • - → n dS = S (τ ij I j -pI i ) • - → n dS + V ρg i dV , (3.1b) ∂ ∂t V c i dV + S c i ( - → U - - → U d ) • - → n dS = 0, (3.1c)
where V is the domain of interest, or control volume, bounded by the closed surface S moving at the velocity -→ U d with a unit normal vector -→ n directed outward.

-→ U and p represent, respectively, the velocity and pressure fields. τ ij and g i are the components of the viscous stress tensor and the gravity vector, whereas I j is a vector whose components vanish, except for the component j which is equal to unity. c i is the volume fraction for fluid i and is used to distinguish the presence (c i = 1) or the absence (c i = 0) of fluid i. Since a volume fraction between 0 and 1 indicates a mixture of two fluids, the value of 0.5 is selected as a definition of the interface between the fluids.

The effective flow physical properties (viscosity and density) are obtained from the properties of each phase (µ i and ρ i ) with the following constitutive relations:

ρ = i c i ρ i ; µ = i c i µ i ; 1 = i c i . (3.2)
When the grid is moving, the so-called space conservation law must also be satisfied:

∂ ∂t V dV - S - → U d • - → n dS = 0. (3.3)
A simplified form of the general mass conservation equation (3.1a) can be obtained when considering incompressible phases with constant densities ρ i . From the constitutive relations (3.2) it is possible to isolate one arbitrary phase j with ρ j = 0:

c j = 1 - i =j c i , (3.4a) ρ = c j ρ j + i =j c i ρ i = ρ j + i =j c i (ρ i -ρ j ).
(3.4b)

Governing equations

Substituting the previous relations (3.3) and (3.4b), as well as the volume fraction equation (3.1c) into the global mass conservation equation (3.1a) yields:

0 = ∂ ∂t V ρ j + i =j c i (ρ i -ρ j ) dV + S ρ j + i =j c i (ρ i -ρ j ) ( - → U - - → U d ) • - → n dS, = ρ j ∂ ∂t V dV - S ( - → U - - → U d ) • - → n dS + i =j (ρ i -ρ j ) ∂ ∂t V c i dV + S c i ( - → U - - → U d ) • - → n dS , = ρ j S - → U • - → n dS .
(3.5)

Thus, mass conservation simplifies to:

S - → U • - → n dS = 0, (3.6)
or, in a non-integral form using the divergence operator D () ,

D - → U = 0.
(3.7)

Turbulence closure equations

In Reynolds-averaged form, the Navier-Stokes equations contain a contribution of the turbulent velocity fluctuations which is usually assimilated into the stress tensor:

τ ij = τ t ij + τ l ij , (3.8a
)

τ l ij = 2µ(S ij -S nn δ ij /3), (3.8b) τ t ij = -ρu i u j . (3.8c)
The closure of the Reynolds-averaged equations requires the definition of the turbulent Reynolds stresses τ t ij in (3.8a), in terms of known quantities in a physically consistent way.

For turbulent eddy viscosity closures, the Reynolds stress tensor is considered proportional to the mean strain-rate tensor S. All of the commonly used Reynolds stress models in this class are based on the Boussinesq hypothesis in terms of the eddy viscosity µ t where it is assumed that:

τ t ij = 2µ t (S ij -S nn δ ij /3) -2ρKδ ij /3. (3.9)
Considering incompressible flows, previous relations lead to (3.10) and, if the isotropic part of the Reynolds stress tensor (2K/3) is not explicitly needed, it can be simply absorbed into the mean pressure in (3.1b). In this case, the equations for the stress tensors reduce to the simplified form:

τ ij = τ t ij + τ l ij , τ l ij = 2µS ij , τ t ij = 2µ t S ij -2ρKδ ij /3.
(3.10)

Numerical framework

This section shows the basic approach used to transform the flow equations of the preceding section into a discretized system of equations for the velocity, pressure, and volume fractions. A pressure equation can be derived from the conservation equation (3.6) where the volumetric flux reconstruction -→ U f .

-→ S f involves a coupling between the pressure gradient and the velocity at the faces. While the methodology is based on the Rhie and Chow SIMPLE [START_REF] Rhie | A numerical study of the turbulent flow past an isolated airfoil with trailing edge separation[END_REF] algorithm, special attention is given (i) to the pressure equation formulation in presence of a density discontinuity, (ii) the way unsteady terms ∂/∂t, and (iii) pseudo-unsteady terms ∂/∂τ are interpolated so that, when an overall steady solution is expected the solution does not depend on the time step ∆t and on the local fictitious time step ∆τ . The originality of the method proposed by [START_REF] Queutey | An interface capturing method for free-surface hydrodynamic flows[END_REF] is to take into account the pressure gradient discontinuity by using a pressure equation based on -→ ∇p/ρ rather than -→ ∇p itself in the context of the SIMPLE algorithm.

Semi-continuous form of the momentum equations

After isolating the pressure gradient together with the gravity acceleration, the semidiscretized momentum equation for any Cartesian velocity component reads: (3.11) where e c , e p and e q are the weights of the current and the two previous time steps in the time scheme; {a C , a nb } are the matrix coefficients from the implicit part of the diffusive and convective terms; -→ S is a source term containing all explicit remaining contributions and 3.4 Numerical framework external force fields except gravity and pressure. A compact and useful form is introduced to express the velocity at the cell-center: (3.12) where the discretized vector -→ Û , homogeneous to gravity acceleration, includes part of the diffusion, convection and source terms. Unsteady and pseudo-unsteady contributions are explicitly kept and not transfered into this intermediate vector.

(e c + 1/∆τ C )(V ρ - → U ) c C + (eV ρ - → U ) p C + (eV ρ - → U ) q C + a C - → U c C + nb a nb - → U c nb + - → S C + (V - → ∇p) c C = (ρV ) c C - → g + (ρV ) c C - → U c0 C /∆τ C ,
- → U c C = -Cp C - → Û C + ( - → ∇p/ρ) c C -- → g + Cp C - → U c0 C /∆τ C -Cp C (eV ρ - → U ) p C + (eV ρ - → U ) q C /(ρV ) c C , with Cp C = (e c + 1/∆τ C + a C /(ρV ) c C ) -1 ,
- → Û C = nb a nb - → U c nb + - → S C /(ρV ) c C .
(3.13)

Reconstruction of volumetric fluxes

The velocity vector on the cell faces is expressed as the face interpolate of equation (3.12):

- → U c f = -Cp f - → Û f + ( - → ∇p/ρ) c f -- → g + Cp f - → U c0 f /∆τ f -Cp f e p (ρV ) p f - → U p f + e q (ρV ) q f - → U q f /(ρV ) c f , with Cp f = e c + 1/∆τ f + (a C ) f /(ρV ) c f -1 . (3.14)
Except the pressure gradient, all terms in the previous equation are interpolated from the available cell quantities (L, R) on both sides of the face with a central operator CLR() . For instance, Cp f = (e c + 1/CLR(∆τ ) + CLR(a C )/CLR(ρV )) -1 .

Then, the reconstruction of volumetric fluxes F( -→ U ) reads:

F( - → U ) = -Cp f F( - → Û ) -Cp f F( - → ∇p/ρ -- → g ) + Cp f F(U c0 )/∆τ f -Cp f e p (ρV ) p f /(ρV ) c f F( -→ U p ) -Cp f e q (ρV ) q f /(ρV ) c f F( -→ U q ).
(3.15)

Pressure equation

When the velocity reconstruction (3.14) is substituted into the continuity equation (3.7), the pressure equation is obtained as

-D Cp - → ∇p/ρ =D Cp - → Û -D Cp --→ U c0 /∆τ -D (Cp - → g ) +D Cp [(eρV ) p /(ρV ) c ] - → U p + D Cp [(eρV ) q /(ρV ) c ] - → U q .
(3.16)

Consequently, provided that pressure equation (3.16) is satisfied, the volumetric flux F( -→ U ) defined by the interpolation (3.15) is guaranteed to be conservative. Using classical Gaussian integration with adequate differencing for the normal gradient to the face, as worked out in the following section (equation (3.20)), yields the discretization over the current control volume: (3.17) where the term involving E p represents the explicit part of the pressure gradient flux. Unsteady and pseudo-unsteady fluxes are gathered into the flux F i . The matrix assembled from all control volumes is sparse, symmetric, and positive definite so that conjugate gradient based iterative solvers can be used.

- f C p f S f p R -p L hρ = f C p f F( - → Û ) + S f E p hρ -F( - → g ) + F i ,

Algorithm

The discretization of mass and momentum conservation equations yields a set of algebraic equations: one for each control volume and for each transport/conservation equation. These nonlinear and coupled equations are solved by the following segregated algorithm:

1. Initialize flow field quantities Q 0 at t = t 0 .

2 While this point-wise solver works well for linear systems arising from discretized transport equations, its efficiency is dramatically reduced for the pressure operator. This operator is elliptic and requires the complete solution of a Laplace-like equation for each step. Moreover, the corresponding linear system is ill-conditioned for highly stretched grids. This is why, depending on the behavior of the pressure equation, the flow solver uses either an Algebraic Multigrid solver or a PGMRES algorithm with Incomplete LU(k) preconditioning.

Reconstruction at the cell faces

The major difficulty when solving both air and water in the same continuum is to obtain a perfect equilibrium between the pressure gradient and the gravity term to prevent the growth of parasitic currents due to gravity, even when neglecting surface tension and viscosity effects in jump conditions. This section first explains the reconstruction of pressure and pressure gradients at the cell faces, in order to deal correctly with gravity and large density variations. Then, it details the procedures for the face reconstruction of the convective terms, such as the velocity, turbulence intensity, and the volume fractions.

Pressure equation

In the spirit of [START_REF] Rhie | A numerical study of the turbulent flow past an isolated airfoil with trailing edge separation[END_REF], the pressure equation (3.17 values with central interpolation. The pressure itself is then solved (see section 3.4.4) from the matrix assembled from all control volumes using the specific discretization of the normalized pressure gradient through the face, as explained below. The gravity term must be kept along with the pressure gradient so that a pure hydrostatic equilibrium will be exactly satisfied if the continuous hydrostatic equilibrium -→ ∇p/ρ = -→ g is satisfied.

Basis of pressure reconstruction In continuous flow, the water surface appears as a density discontinuity. Over such a discontinuity, the following jump conditions hold:

[p] = 0, (3.18a) Face reconstruction of p From Taylor series expansion on both sides of the face, combined with the jump conditions (3.18), a reconstruction of the pressure on the face can be established [START_REF] Queutey | An interface capturing method for free-surface hydrodynamic flows[END_REF] in the following compact form involving left and right side cell-centered data only (see figure 3.1 for notations). This pressure reconstruction will be used in the momentum equations.

- → ∇p ρ = - → 0 . (3.18b) - → ∇p itself,
p f = h + ρ + p L + h -ρ -p R h + ρ + + h -ρ - + ρ + ρ - ρ h -- → E + -h + - → E - h . h + h - → ∇p ρ L + h - h - → ∇p ρ R . (3.19a)
The framed term is kept explicit in each solution step, while the non-framed term is implicited in the solver. Geometrical vectors -→ E ± are introduced so that the framed term contribution goes to zero when the grid becomes orthogonal (

-→ Lf . - → n = -→ f R. - → n = 0): - → E - -→ Lf . - → n - → n - -→ Lf , - → E + -→ f R. - → n - → n - -→ f R. (3.19b)
Distances used are the projected distances to the face h ± and the projected distance h between the L and R cell centers:

h -= -→ Lf . - → n , h + = -→ f R. - → n , h = h -+ h + = -→ LR. - → n . (3.19c)
The quantity ρ homogeneous with ρ is defined by: .19d) 3.5 Reconstruction at the cell faces

ρ = h -ρ -+ h + ρ + h . ( 3
Normal gradient to the face The discretization of the pressure gradient in equation (3.17) is obtained with a reconstruction following the same rules as for the quantity on the face. The continuous term through the face is the gradient normal to the face normalized by ρ:

( - → ∇p. - → n /ρ) f : - → ∇p. - → n ρ f = 1 ρ p R -p L h + - → ∇p L . - → E -+ - → ∇p R . - → E + ρh . (3.20)
Here again, the framed (explicit) term contribution goes to zero when the grid becomes orthogonal. The non-framed term is the implicit part that goes into the matrix for the pressure equation.

Discretization of convection equations

Discretizations of the convective terms in the momentum, turbulence, and volume fraction equations require special reconstructions of the cell-centered values to the faces. To guarantee accuracy, stability and boundedness of the solutions, these schemes are developed in the Normalized Variable Diagram (NVD).

Q C Q f Q U Q D f C D U Flow direction (CD) (UD) (DD)

NVD diagram and boundedness considerations

When constructing a face reconstruction scheme, in order to avoid unrealistic oscillations, especially in the volume fraction c i which is discontinuous in nature, the search for an acceptable compromise between accuracy and boundedness is a key point [START_REF] Jasak | High resolution NVD differencing scheme for arbitrarily unstructured meshes[END_REF][START_REF] Pržulj | Bounded convection schemes for unstructured grids[END_REF][START_REF] Darwish | The x-schemes: A new consistent high-resolution formulation based on the normalized variable methodology[END_REF]. A practical way to introduce a new numerical implementation is to consider the Normalized Variable Diagram (NVD) analysis [START_REF] Leonard | Simple high-accuracy resolution program for convective modelling of discontinuities[END_REF][START_REF] Darwish | Normalized variable and space formulation methodology for high-resolution schemes[END_REF] and to follow the rules that enforce local monotonicity and the Convection Boundedness Criterion (CBC) [START_REF] Gaskell | Curvature compensated convective transport: SMART , a new boundedness preserving transport algorithme[END_REF]. The NVD is introduced in general terms here. In the following, the generic quantity Q can be thought of as the volume fraction c i , a velocity component for convection, or a turbulence quantity.

The NVD diagram was originally introduced on structured, one-dimensional grids. On such a grid, in the neighborhood of a face f , points U , C and D are selected according to the flow direction on the face and represent the downwind, central, and upwind cell centers, respectively (see figure 3.2). Then a normalized variable Q is defined as:

Q = Q -Q U Q D -Q U . (3.21)
The idea of the NVD is to represent a reconstruction scheme for the face value by Qf , the normalized value on the face, as a function of QC .

The CBC criterion [START_REF] Gaskell | Curvature compensated convective transport: SMART , a new boundedness preserving transport algorithme[END_REF] corresponds to an area in the NVD in which a scheme must lie in order to be stable and monotone; this area is shaded in figure 3.3a. The shape of this area implies that the first-order upwind differencing scheme (UDS) is the only scheme which unconditionally satisfies the boundedness criterion. The second-order centered differencing scheme (CDS) is only useful in the range 0 QC 1, as is the first-order downwind differencing scheme (DDS). Practical schemes are often designed as blendings of these elementary schemes, to obtain certain desired properties while remaining inside the CBC area. NVD on unstructured grids On arbitrary unstructured grids, the far upstream node U is not known explicitly (C and D are still chosen as the centers of the two cells next to the face, see figure 3.3b). It is not even certain that a node exists in the position where U should lie. Therefore, an alternative is chosen: an imaginary nodal quantity Q U is defined by the use of the gradient projection method [START_REF] Queutey | An interface capturing method for free-surface hydrodynamic flows[END_REF] in such a way that:

Q U = Q C - --→ CU • - → ∇Q| C , with --→ CU - --→ CD. (3.22)
The location of the imaginary point U is found as the point that is the mirror image of D relative to point C.

Examples of NVD schemes

The first NVD scheme implemented in ISIS-CFD was the Gamma Differencing Scheme (GDS) [START_REF] Jasak | Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows[END_REF], which was used in all equations.

The second scheme implemented for the discretization of the convective fluxes in both 3.5 Reconstruction at the cell faces the momentum equations and the equations for turbulence modelling is the AVLSMART [START_REF] Pržulj | Bounded convection schemes for unstructured grids[END_REF]. It has been implemented following the χ -Scheme methodology [START_REF] Darwish | The x-schemes: A new consistent high-resolution formulation based on the normalized variable methodology[END_REF], and it was demonstrated to have an improved convergence behavior without loss of accuracy in many situations. If the base scheme for the GDS scheme is the second-order Central Differencing scheme (CD), the base scheme for AVLSMART is the third-order QUICK scheme [START_REF] Leonard | A stable and accurate convective modelling procedure based on quadratic upstream interpolation[END_REF].

For the face reconstruction of the volume fraction c i in the conservation equation (3.1c), a compressive scheme BICS is used based on the robust Gamma Differencing Scheme (GDS) and the Inter-Gamma scheme [START_REF] Jasak | Interface tracking capabilities of the inter-gamma differencing scheme[END_REF] for its compressive properties.

The special aspect of BICS and similar schemes is that they preserve the sharpness of the discontinuity in the volume fraction. A later scheme BRICS [START_REF] Wackers | Free-surface viscous flow solution methods for ship hydrodynamics[END_REF] is derived from BICS but with improved computation of the upwind cell value.

Dynamic meshing

For the simulation of complex problems, it is often necessary to dynamically adjust the mesh during a computation. ISIS-CFD includes for example a mesh deformation technique to accommodate the motion of bodies and adaptive mesh refinement to dynamically adjust the local mesh resolution. To allow the relative motion of multiple bodies, the mesh can be divided in different domains which are connected through sliding interfaces (for rotation) and overset meshes (for arbitrary motions). In the context of this thesis, the first two techniques are the most important so they are described in some detail here.

Mesh deformation

To accomodate rigid-body translation and rotation of ships and other bodies, a mesh deformation technique is introduced by [START_REF] Leroyer | Numerical methods for RANSE simulations of a self-propelled fish-like body[END_REF]. This technique displaces the nodes of an existing mesh using analytical weighted mesh deformation. A weighting coefficient is computed which is equal to 1 on the moving body and 0 on the other boundaries, while it varies smoothly in the interior of the domain. This weighting coefficient is obtained as a solution of the Laplace equation, discretized and solved with the same techniques as the pressure equation (section 3.4.3).

The rigid-body displacement associated with the moving body is then computed in each node and multiplied with the weighting coefficient to obtain the node displacement. Given the boundary conditions for the weighting, this implies that the mesh on the body follows its movement, while the mesh on the other boundaries remains stationary. Usually, this weighted deformation is applied for certain degrees of freedom only, and combined with block motion for other degrees of freedom (for example, horizontal translations).

3.6.2 Grid refinement procedure and refinement criteria

An anisotropic automatic grid refinement procedure has been developed which is controlled by various flow-related criteria (Wackers et al., 2014a). The procedure is integrated completely in the flow solver. During a computation, the refinement procedure is called repeatedly. In such a call, first the refinement criterion is calculated, then in a separate step of the procedure the grid is refined based on this criterion.

For anisotropic refinement, a metric tensor is used as refinement criteria to specify different cell sizes in different directions. The refinement criterion in each cell is a 3 × 3 symmetric positive definite matrix, computed from the flow solution. This matrix is interpreted as a geometric transformation: the modified dimensions of a cell are computed by multiplying the vectors between the opposing face centers in the three cell directions with this matrix. Grid refinement is then applied in order to make the modified dimensions of all cells equal. The result in the real space is an anisotropically refined, adapted mesh.

Adaptive refinement can be used in different ways. For example, locally fine cells can reveal physical details of a flow which are hard to capture on non-adapted meshes. Furthermore, adaptive refinement can simplify the computational setup, because a part of the final mesh is created automatically. For the simulation of hydrofoils, a typical use of adaptive refinement is to capture the free surface with an adaptively generated fine grid (figure 3.4). This ensures that the grid at the surface is fine enough, independent of the exact shape and size of the waves. Furthermore, it becomes possible to simulate a foil at different immersion depths using the same original mesh, which reduces the perturbations of the results due to a remeshing for each position [START_REF] Wackers | Hessian-based grid refinement for the simulation of surface-piercing hydrofoils[END_REF]. 3.6 Dynamic meshing

ISIS-CFD parallelization

Parallelization is based on domain decomposition. The grid is divided into different partitions, which contain the cells. This is done with the help of the Metis partitioning algorithm [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF][START_REF] Schloegel | A unified algorithm for load-balancing adaptive scientific simulations[END_REF]. The interface faces on the boundaries between the partitions are shared between the partitions; information on these faces is exchanged with the MPI (Message Passing Interface) protocol. This method works with the sliding grid and overset approaches and the different sub-domains can be distributed arbitrarily over the processors without any loss of generality. Moreover, the automatic grid refinement procedure is fully parallelized with a dynamic load balancing working transparently with or without sliding grids.

This chapter introduces the Gaussian process regression model. A Gaussian process can be used to reconstruct a function over a domain, based only on the function values in a limited number of points. For this, the key assumption is to consider the (deterministic) function as if it is a realization of a stochastic process. Thanks to this probabilistic context, the Gaussian process not only provides a reconstruction of the function, but also an estimate of the data spread in each point, which indicates the uncertainty in the reconstruction. In subsequent chapters, the Gaussian process will be used to construct metamodels based on CFD simulations. Here, the approach is presented and several variations are introduced.

The chapter ends with a test to find the most suitable Gaussian process to fit the type of data encountered in the following chapters.

Historical overview

Gaussian processes (GP) are a general class of probability distributions over functions. As such they are a generalization of probability distributions, which describe finite-dimensional random variables.

Gaussian processes have been widely studied and used for many purposes. Being a type of stochastic process, they were first used for time series prediction. Work in this area dates back to the 1940s [START_REF] Wiener | Extrapolation, interpolation, and smoothing of stationary time series[END_REF]. Spatial data reconstruction based on Gaussian processes has been widely used since the 1970s in the field of geostatistics where it is termed kriging, named after the South African mining engineer D. G. Krige [START_REF] Matheron | The intrinsic random functions and their applications[END_REF]. The idea of global optimization based on Gaussian process models also dates back to the 1970s [START_REF] Močkus | On bayesian methods for seeking the extremum[END_REF]. Gaussian process models were then extended by statisticians to more general multivariate input regression problems [START_REF] O'hagan | Curve fitting and optimal design for prediction[END_REF].

Before the 1990s the use of Gaussian process models in statistics was largely confined to the spatial statistics sub-field. Since the 1990s however, increasing work with kernel methods and Bayesian inference applied to machine learning has been widening the use of Gaussian process regression and [START_REF] Williams | Gaussian processes for regression[END_REF] first described Gaussian process regression in a machine learning context. The definitive book on GPs in the context of machine learning is [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF]. Most of the material covered in this chapter can be found there.

Definition of Gaussian processes

Stochastic processes are sets of (potentially infinitely many) random variables f (x). As such, they are the stochastic equivalent of functions. An alternative way of seeing them is as distributions over infinitely long vectors. The single outcome of a stochastic process is called a sample function or realization. Gaussian processes are a particular type of stochastic processes.

Before introducing the Gaussian process, the finite-dimensional context is presented, where the univariate Gaussian (normal) distribution is given by the probability density function:

p(y|µ, σ 2 ) = 1 √ 2πσ 2 e (y-µ) 2 2σ 2 . (4.1)
This distribution is fully characterized by the mean µ and the standard deviation σ 2 , thus one often writes a stochastic variable f with a Gaussian distribution as:

f ∼ N (µ, σ 2 ). (4.2)
A multivariate Gaussian distribution is a finite set of coupled stochastic variables f for which any linear combination of the variables has a univariate Gaussian distribution. It is denoted as:

f ∼ N (µ µ µ, Σ), (4.3)
with the mean vector and the covariance matrix given by:

µ i = E[f i ], (4.4a) Σ i,j = E[(f i -µ i )(f j -µ j )].
(4.4b)

Here E denotes the expected value of a stochastic variable.

A Gaussian process is defined as a stochastic process for which any finite subset has a joint multivariate Gaussian distribution. Gaussian processes can be considered as an infinitedimensional generalization of multivariate normal distributions. This definition implies that a Gaussian process is fully specified by its mean function:

µ(x) = E[f (x)], (4.5a)
and its (positive definite) covariance function, also called kernel:

k(x, x ) = E[(f (x) -µ(x))(f (x ) -µ(x ))]. (4.5b)
A Gaussian process will be denoted as:

f (x) ∼ GP(µ, k(x, x )). (4.6)
In Gaussian process modeling, the mean function µ(x) is often assumed to be zero. With this assumption, what relates one observation to another is just the covariance function k(x, x ). The type of data that can be captured by a GP model is determined by this function, which expresses the assumed similarity between observations (similar observations should have similar target values). A zero mean prior is not a requirement -non-zero mean functions µ(x) can be used if a priori knowledge of the response structure to be fitted is available. Using a zero mean just reflects prior knowledge that the function is equally likely to be positive or negative.

Gaussian processes are useful in statistical modeling because of their properties inherited from the normal distribution. The distributions of various derived quantities can be obtained explicitly and the consistency of the Gaussian process distribution implies that the usual rules of probability apply to the collection of random variables.

Gaussian process modeling

A Gaussian process can be used to interpolate between data known only in some points of a domain. Even if these point data are issued from a deterministic source (as is often the case with a computer code), the Gaussian process regression (GPR) model treats them as a realization of a random function. Predictions in the rest of the domain are then based on the conditional probability of the process, given that the values in these points are fixed. Let us consider n realizations of a computer code. To each output realization y i corresponds a d-dimensional input vector of geometrical and flow parameters x i = (x 1 , ..., x d ) i . Those n observations are constituting the data set D n = {(x i , y i )} n i=1 used to train the GP model.

Gaussian process modeling

The basis of GP modeling is the assumption that the response of the computer code over the entire domain can be represented as a realization of a Gaussian process, which means by definition that the responses in a finite number of points form a sample from a multivariate Gaussian distribution. Consider now the probabilistic distribution in the points of D n plus one arbitrary test point x * :

f n f * ∼ N µ µ µ n µ * , K n K T * K * K * * , (4.7)
with the covariance matrix for the points of D n :

K n =     k(x 1 , x 1 ) • • • k(x 1 , x n ) . . . . . . . . . k(x n , x 1 ) • • • k(x n , x n )     (4.8a)
and the covariance and autocorrelation for the point x * :

K * = k(x * , x 1 ) k(x * , x 2 ) • • • k(x * , x n ), K * * = k(x * , x * ). (4.8b)
The objective of the Gaussian process modeling is then to provide an estimate of f * based on f n . With the Gaussian process prior on f , this requires to compute the conditional probability (the posterior in Bayesian statistics) over the value f * at the query input x * , given that f n is fixed and equal to the y i in D n . [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF] show that the conditional probability p(f * |f n ) follows the Gaussian distribution:

f * |f n ∼ N (µ * + K * K -1 n (f n -µ µ µ n ), K * * -K * K -1 n K T * ).
(4.9)

The best estimate for f * is the mean of this distribution:

f * = µ * + K * K -1 n (f n -µ µ µ n ), (4.10)
which, under the assumption that µ = 0, becomes:

f * = K * K -1 n f n , (4.11)
and the uncertainty of the estimate is given by its variance:

var(f * ) = K * * -K * K -1 n K T * . (4.12)
For these estimates, the covariance function k(x, x ) is not known in general and must be estimated, which influences the accuracy of the prediction. The choice of the covariance kernel, a key aspect of GP reconstruction, is further discussed in section 4.5. Most of the time, the response of a computer code is deterministic, i.e. re-running the code will give the same output provided the inputs were identical. However, simulations always contain numerical errors coming from the discretization and the mesh. Since these numerical errors have an unpredictable behavior, they can be considered as noise in the training set D n .

While Gaussian process modeling is initially an interpolating model, it is also possible to make it regressive in a fairly simple way. Suppose that the observation f n consists of f (x) plus random Gaussian noise ∼ N (0, σ 2 N ). The usual assumption is to consider the noise as uncorrelated in space1 . As a consequence, the covariance matrix of the noise is diagonal and each non-zero element is the variance of the error. The probability distribution for the data set and the test point becomes:

f n f * ∼ N µ µ µ n µ * , K n + σ 2 N I K T * K * K * * . (4.13)
The variance σ N is not added to the autocorrelation K * * , since the goal is to predict f * unpolluted by the noise. Under these assumptions, the prediction in x * becomes:

f * = K * K n + σ 2 N I -1 f n , (4.14)
and the variance is:

var(f * ) = K * * -K * K n + σ 2 N I -1 K T * . (4.15)
The main difference with the noise-free case is that best predictor is not an interpolator anymore: the predictions (4.14) in the observation points do not correspond to f n and the variance (4.15) at the observation points becomes non-null. The noise term acts as a smoothing effect (figure 4.2).

Covariance functions

The covariance function k(x, x ) = cov(f (x), f (x )), also called kernel function or covariance kernel, is the center part of the Gaussian process model. It expresses the similarity or dependencies between data points and, as such, contains the assumptions about the modeled function. Covariance functions are generally not fully fixed beforehand, but specified in terms of unknown hyperparameters θ, which are estimated from the data. This estimation corresponds to the learning part of the Gaussian model (see section 4.6). A kernel is a function of two inputs x and x . Two categories of kernels can be distinguished: stationary kernels depend only on the distance of two datapoints and not on their absolute values, i.e. k(x, x ) = k(d(x, x )) and are thus invariant to translations in the input space.

Non-stationary kernels depend also on the specific values of the datapoints.

The choice of the covariance function should be made in agreement with the known properties of the source from which the data points come. The particular choice of covariance function determines the properties of sample functions drawn from the GP prior (e.g. smoothness, lengthscales, amplitude etc). Thus, it is an important part of GP modeling to select an appropriate covariance function for a particular problem. Chapter 4 of [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF] gives more details and a list of correlation functions with their advantages and drawbacks.

Kernel parameters

Each kernel has a number of parameters which specify the shape of the covariance function.

They are referred to as hyper-parameters since they specify a distribution over functions, instead of being parameters which directly define a function. Apart from the parameters that are specific to each kernel, some generic parameters appear in most kernel definitions.

Length scale Most correlation functions contain a length scale, which gives an order of magnitude for the distance beyond which two points become uncorrelated. For short length scales, points have a local influence only while for large length scales, many points are correlated and the GP realizations will be smooth. A kernel can use different length scales in different directions, to take into account that a function to be reconstructed varies less in some directions than in others. This is an example of an asymmetric kernel.

Signal variance

The signal variance σ 2 is a scaling factor of the kernel k(x, x ). It determines the variation of function values from their mean: small values of σ 2 characterize functions that stay close to their mean value, while larger values allow more variation.

In the absence of noise, the prediction of a GP model does not depend on this scaling. Since both K * and K n in equation (4.11) scale with σ 2 , this parameter disappears from the expression for f * .

Noise level While the noise covariance σ 2 N I is treated separately from the kernel k(x, x ) in equations (4. 14) and(4.15), it has a major influence on the prediction. In the presence of noise, the ratio between σ and σ N determines the treatment of oscillations in the function: if the noise variance is too high with respect to the signal variance, then actual oscillations in the function will be filtered as noise. On the other hand, if the signal variance is too large, the modeled function will chase outliers created by the noise.

Examples of kernels

Squared exponential (SE) kernel

This kernel is also known as the radial-basis function (RBF) kernel, the Gaussian kernel, or the exponentiated quadratic. It is defined by:

k(x, x ) = σ 2 exp(- x -x 2 2 2 ), (4.16)
with hyperparameters θ = {σ, }. Since the SE kernel is infinitely differentiable, it is very smooth (figure 4.3). Therefore, other classes of kernels are needed to fit less regular functions.

The Matérn kernel

The Matérn kernels are a more flexible class of kernels [START_REF] Matérn | Spatial variation[END_REF]. The Matérn kernel is defined by the general form:

k(x, x ) = σ 2 2 ν-1 Γ(ν) 2 √ ν x -x ν H ν 2 √ ν x -x , (4.17)
where Γ is the gamma function and H v the modified Bessel function of the second kind.

The hyperparameters are θ = {σ, , ν}. The smoothness of this kernel increases with the parameter ν and for ν → ∞, the SE kernel is recovered. For half-integer values of ν, the expression is simplified. Notably for ν = 3 2 the Matérn kernel becomes:

k(x, x ) = σ 2 1 + √ 3 x -x exp - √ 3 x -x . (4.18)
Like the SE kernel, the Matérn kernel is stationary and in the form presented above, with the same length scale for all dimensions, the kernels are isotropic. The shape of the Matérn kernel is pointier than the SE kernel (figure 4.3).

Covariance functions

Learning the kernel hyperparameters

For GP modeling, the kernel functions are chosen a priori based on assumptions about the behavior of the data. Thus, the covariance function k(x, x ) does not express the real covariance between points in the domain, which limits the accuracy of the predictions (4.14) and especially of the variance estimation (4.15). However, the kernel hyperparameters can be freely adjusted. To obtain the best predictions, the hyperparameters need to be optimized such that the Gaussian process mimics the behavior of the real function as closely as possible.

The most commonly advocated approach for hyperparameter optimization is to maximize the marginal likelihood (see [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF], chapter 5). Simplified, this approach implies that p(f n |x 1...n ) is maximized, i.e. the Gaussian process is sought which is most likely to produce the outcome f n in the datapoints x i of D n . For example, the actual data may be oscillating. A Gaussian process with a short length scale is likely to have realizations that oscillate in the same way, while such an outcome is almost impossible for long length scales. Thus, the short length scale would be preferable.

If the chosen kernel has a different length scale for each direction, then estimating the lengthscale parameters 1 , 2 , . . . , d , implicitly determines the relevance of each dimension.

Input dimensions with relatively large lengthscales imply relatively little variation along those dimensions in the function being modeled. Such a procedure is known as Automatic Relevance Determination (ARD), since it automatically determines the dimensions which have the most influence on the outcome.

Implementation and performance evaluation

In subsequent chapters, Gaussian processes will be used to construct response surfaces based on simulations of fluid flow. This section first describes the GP setup that has been chosen for these applications. To evaluate the performance of the Gaussian process regression model, tests are then performed on analytical functions that ressemble the behavior expected from the flow simulations.

Gaussian process setup

In all applications performed here, the Python code is used and Gaussian processes are created with the library GaussianProcessRegressor [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. Both the SE and the Matérn kernel are tested. The regressions are performed without noise and the variance of the kernels is fixed at σ = 1. For the Matérn kernel, the parameter ν = 3 2 is also kept fixed.

A single length scale is used, so the kernels are isotropic. The length scale is optimized for the data being fitted. Since no ARD is performed, the tests here are performed on square domains and on functions which vary in all directions. In the applications of chapter 5 and 7, the input variables are always scaled such that the domain for the Gaussian process has a unit dimension in each direction. The cases where one dimension has little influence on the output (figure 7.3 for example) are handled relatively well by the isotropic kernels.

Test cases

The tests are performed on two-dimensional analytic functions, given below. Evaluations of these functions on regular grids of points are used as datasets for GP regression. The quality of the regression as a function of the number of points is studied.

Three hump camel This smooth test function is characterized by two minimum troughs with a central ridge; it increases sharply away from the origin (figure 4.4 left). The domain for the input is (x, y) ∈ [-2.5, 1.5] and the Three hump camel function is given by:

f (x, y) = 2x 2 -1.05x 4 + 1 6 x 6 + xy + y 2 . (4.19)

Implementation and performance evaluation

Gramacy-Lee function The Gramacy-Lee test function [START_REF] Gramacy | Gaussian processes and limiting linear models[END_REF] is more complex than the previous one. First, it consists of two sharp peaks with the minimum located close to the maximum. And second, the function close to zero on most of the domain (figure 4.4 right). This function is evaluated on the square (x, y) ∈ [-2, 6] and its expression is:

f (x, y) = x exp(-x 2 -y 2 ). (4.20)

Results

To assess the performance of the different candidate kernels, the GP models are used to predict function values on a regular grid that is much finer than the grids used to construct the model. These predictions are compared with the real values. Two measures for the error are studied, the maximum error and the statistical metric R-squared:

R 2 = 1 - n i=1 (y i -y i ) 2 n i=1 (y -y i ) 2 , (4.21)
where y i denotes the true observations on the fine grid, y i the values predicted by the GP model in these points, and y their empirical mean. The R-squared, or coefficient of determination, is a statistical measure of how close the data are to the fitted regression model. R-squared with a value of 1 indicates that the model perfectly fits the data.

Gramacy-Lee Figure 4.5 shows the convergence of the interpolation for the Gramacy-Lee function. This function has two close, sharp peaks which are only visible on fine enough grids; the coarser grids do not capture the peaks at all. A minimum of nine points per dimension is needed to represent the shape of the peaks correctly. The performance of the two kernels is similar, although there is a notable difference on the finest grid. The SE kernel produces incorrect oscillations next to the peaks around the line y = 0. These wiggles are nearly absent for the Matérn kernel.

The convergence of the error norms (figure 4.7) confirms these observations. The errors are high on coarse grids, with a maximum at 4 points per dimension (this is an unlucky situation where one point lies in the maximum, while the minimum is not sampled -see also section 5.8.2). While the maximum errors for the two kernels are similar, the R 2 measure is systematically lower for Matérn. This is caused by the lower level of oscillations away from the peaks.

Three hump camel In figure 4.6, the convergence is shown for the Three hump camel function. This function is smoother than Gramacy-Lee which means that it is better represented on coarse grids. On the other hand, contrary to the peaks of Gramacy-Lee, this function is not shaped like GP kernels; it has a polynomial behavior with maximum values on the borders. This creates problems for the SE kernel, which has marked oscillations that remain present on fine grids. The Matérn kernel does not have these wiggles at all.

48

Chapter 4 Gaussian process modeling Thus, the difference in the errors (figure 4.8) between the kernels is much more marked. While both kernels produce smooth convergence, the Matérn kernel is more precise. The difference is most noticeable on the finer grids, where both the maximum error and the R 2 measure converge faster for Matérn than for SE due to the oscillations. Thus, the ratio between the errors increases on finer grids.

Choice of the kernel In these tests, the Matérn kernel with ν = 3 2 (or at least, its implementation in GaussianProcessRegressor) appears as the best choice. The difference with SE for the Three hump camel test is important, since this function ressembles the behavior of the drag for an airfoil (see for example figure 7.3). Thus, the Matérn kernel can be expected to reduce oscillations also for hydrofoil optimization. If peaks appear in the function (like in figure 7.2), the performance of the kernel remains good.

The Matérn kernel is recommended by [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF] as the most generally applicable kernel. Based on the current test results, it is selected for all further GP regressions in this thesis.

Implementation and performance evaluation

Conclusion

Constructing simplified models from simulation data requires the capacity to interpolate between the data points. This chapter presented Gaussian Process regression as a way to perform such interpolations.

GP regression treats the function to be reconstructed as if it is a realization of a stochastic process and bases its estimation on the most likely values of this process, given that the function values in the data points are known. This approach may seem counter-intuitive if the actual data come from a deterministic source such as numerical simulation, but it has the advantage of providing an uncertainty estimation for the interpolation based on the variance of the statistical process. As confirmed by the examples in this chapter, the maxima for this uncertainty appear in the middle of the intervals between points.

Furthermore, the GP paradigm offers a straightforward way of treating noise in the input data, by adding a diagonal component to the correlation matrix to model noise that is uncorrelated between the data points. Although the noise modeling is not used in the remainder of this thesis, it is potentially interesting as a way to reduce perturbations due to numerical errors in the simulations.

Finally, tests show that the implementation of the GP regression adopted here is successful in reconstructing functions with either smooth variations or local peaks. However, the results depend on the number of data points and fine grids are needed to capture peaks. Since Gaussian Process regressions do not require a regular distribution of the points, it would be possible to concentrate the data near the most important features of the function to be interpolated. This raises the question of the optimal point placement for a given function. The next chapter addresses this question.

50

Chapter 4 Gaussian process modeling 5 Surrogate based optimization Ship design has become increasingly complex and demanding in terms of fluid dynamics requirements. Since small improvements in hydrodynamic efficiency lead to major cost savings, a far-going optimization of ship designs is required nowadays. Therefore, the use of empirically derived methods to evaluate new designs is most of the time irrelevant as these do not predict the performance accurately enough to allow precise optimization. Fortunately, the massive development of computational power has made numerical simulation widely available with computational tools such as CFD or FEM, see chapter 3. Ship design is now heavily based on numerical simulations.

Building on that trend of simulation based design and helped by the rapid innovation in algorithms and optimization techniques, engineers have begun to take a keen interest in the possibilities offered by automatic simulation-based optimization. However, accurate highfidelity numerical simulations are usually time-consuming and computationally expensive. The resolution of optimization cases with many design variables where the response of the system is moreover non-linear, leading to a large number of design evaluations, can be prohibitive. Furthermore, many ships are designed as one-off units, which means that the budget for their design is limited. Thus, ship design requires efficient optimization methods.

This chapter is dedicated to surrogate-based optimization, which is a combination of surrogate modeling and optimization. A surrogate model (chapter 4) is constructed based on a limited number of CFD simulations and the optimum is sought using this surrogate model instead of the expensive CFD computations. After an introduction of simulation-and surrogate-based optimization, the chapter presents an adaptive technique to efficiently and reliably construct surrogates for different optimization purposes. Furthermore, constrained optimization is considered. The chapter ends with tests of the adaptive algorithm on analytical functions.

Optimization: overview

In general, design optimization implies searching for design parameters which give the best performance in terms of a certain measure. From a mathematical viewpoint, this is written as follows. Let x ∈ X be a set of parameters which define a design, the so-called input parameters. The domain X contains all feasible designs. Furthermore, let f (x) be a function of x. Then the optimization implies finding the value of x which minimizes f :

x opt = arg min x∈X f (x).
(5.1) Some restrictions may apply on which values of x ∈ X are acceptable. This case is called constrained optimization, it is discussed in section 5.7.

Local minimum

f(x)
x

Local minimum

Global minimum Figure 5.1 shows the optimum of a one-dimensional function. Furthermore, it illustrates the possibility of suboptimal local minima, which form a practical difficulty for optimization algorithms that may confuse them with the global minimum.

Simulation-based design and optimization

Simulation-based design makes use of numerical simulation tools to evaluate the performance of a design or design alternatives. While these numerical simulation methods were initially only used to analyse and validate the final design, they are applied more and more to determine the performances of design alternatives in order to select the best one. Integrated in the design process from the earliest conceptual design phases, they help engineers to make informed design decisions; this design process is said to be simulation-based. However, simulation-based design remains limited to a few design evaluations and even though it allows for design improvements, it may not lead to a true optimum.

Simulation-based optimization combines numerical simulations with automatic optimization techniques to support virtual product development and is now playing an increasingly prominent role in the industry. Simulation-based optimization nevertheless remains challenging, for the following reasons:

• The computational cost of the simulations is high, if they are run with high fidelity codes.

• Most design problems are high-dimensional, which increases the number of simulations required to reach optimum.

Thus, evaluating all the design points needed for an automatic optimization with highfidelity simulation may be prohibitively expensive. The next section presents a way to address this limitation using an approximate model of the simulation results.

Surrogate based optimization

One way to address the issue of the computational cost for simulation-based optimization and to reduce the number of required simulations is to approximate the response of the simulation code with an additional model that will be cheaper to evaluate. The optimization itself will be run on this additional model. In the case of simulation-based optimization, the response is already coming from a model. Approximating that response implies building a model of a model, hence the term meta-model. (The Greek metameaning "after", or "beyond", is a prefix used to indicate a concept that is the abstraction of another concept, used to complete or supplement it). Meta-modeling is also referred to as surrogate modeling or response surface modeling.

Surrogate modeling: mathematical definition

Surrogate modeling is a technique that uses sampled data (obtained by running the numerical simulation) to construct approximation models, which are sufficient to predict the output of an expensive computer code in untested points of the design space.

Let a function f : X → R be defined on a domain X ⊆ R d , like in section 5.1. In the context of this definition, f represents the full model and is considered expensive to compute. The number of input parameters is d, these can contain either geometry parameters (for shape optimization) or operational conditions (such as the velocity or attitude of the system).

Surrogate modeling aims at generating an approximation to f from points where f is known exactly. Thus, consider a dataset

D n of n observations, D n = {(x i , y i )} n i=1
, where x i is an input vector of dimension d and y i = f (x i ) is an evaluation of the function f . Then, for a surrogate-creating algorithm A: [START_REF] Queipo | Surrogate-based analysis and optimization[END_REF] presented an overview of surrogate-based analysis and optimization methods in the aerospace science and engineering field. Such methods can be characterized as either parametric or non-parametric: parametric approaches are based on a priori assumptions about the relationship between the response variable and the design variables, while non-parametric methods rely only on the data itself while making as few assumptions as possible. This can be achieved by using simple local models in different regions of the design space to construct a global model.

f (x; D n ) = A(D n ), ( 5
Among the many methods available, Gaussian process regression (chapter 4) will be used in the following. Gaussian Process regression is a semi-parametric regression method since it contains a sub-optimization step to adjust the hyperparameters.

Sequential versus non-sequential sampling

Any surrogate-based optimization involves three separate steps:

• Performing simulations to obtain the dataset D n of design points and responses,

• Creation of the metamodel through the use of surrogate modeling techniques,

• Optimization over the metamodel.

Here, surrogate-based optimization is different from data-mining based on a fixed data-set: since the responses are obtained from numerical simulation, the design points to evaluate are chosen as a part of the procedure. Furthermore, the performance of the optimization depends on this choice of the dataset. In a one-shot or non-sequential approach, the three steps above are performed separately, in this order. This means that the design points are planned in advance, without knowledge about the response of the simulation.

An alternative approach is the sequential or adaptive sampling (figure 5.2). Here, a few points are simulated in the beginning to create an initial sample plan (see the following section 5.4) and a metamodel is constructed with these points. The quality of this metamodel is then evaluated to find the design point where adding a data point would be most beneficial (figure 5.3). This point is simulated, a new improved metamodel is constructed, etc. This sequence of model-optimizing is repeated until an ending criteria is reached, either a convergence criterion based on the metamodel quality or a maximum number of simulations.

The choice of the points to be added is based on an acquisition function (section 5.5), which usually takes into account both the estimated uncertainty in the metamodel and the proximity to the metamodel optimum. uncertain, then the adaptive sampling will cover the entire design space. This approach is called exploration, since it investigates parts of the design space where the exact response is yet unknown. In the end, this approach leads to a metamodel that is reliable everywhere and suitable as a response surface in a VPP for example (see chapter 1). For optimization on the other hand, if the metamodel quality is sufficient it may be advantageous to concentrate further points close to the optimum, since this is the region of greatest interest. Since this approach exploits the existing metamodel to guide points towards the optimum, it is known as exploitation. The right balance between exploration and exploitation is crucial for an acquisition function.

In all cases, the procedure ends with an optimization over the metamodel in order to find the best design parameters. To check the metamodel, this point can be simulated: if the true response is too far from the metamodel value, the optimum is unreliable and further iterations may be required.

Initial sampling

The first step in any surrogate-based optimization is to generate an initial set of data so as to gain a prior knowledge of the response and to start sampling the following points efficiently. In a one-shot approach, the final metamodel is based on these data only. Many different ways exist to create these initial sample plans.

Full factorial A full factorial sampling method consists in a Cartesian division of the design space. A regular distribution of points is created along each dimension and the sampling plan is the tensor product of these distributions (figure 5.4a). This approach guarantees that the effect of all the design parameters is sampled, as well as the coupling between parameters. A disadvantage is that the number of points grows exponentially with the dimension d of the design space. Furthermore, these distributions suffer from aliasing since oscillations in the function f (x) whose wavelength coincides with the distance between points are not captured by the sampling plan.

Random sequence To prevent aliasing and to reduce the number of points, a random distribution over the design space can be adopted (figure 5.4b). The disadvantage of this choice is that random points tend to be clustered, since there is nothing which forces a minimum separation between the points. Thus, a pure random sample (Monte Carlo) badly fills the space [START_REF] Giunta | Overview of modern design of experiments methods for computational simulations[END_REF]. the discrepancy, i.e. the maximal deviation between the distribution of the sample's points to an uniform distribution. The Hammersley sequence [START_REF] Hammersley | Monte Carlo methods for solving multivariable problems[END_REF] is a lowdiscrepancy sampling method where the coordinates for a point i are found by expressing the number i in various prime bases (base-2, base-3, base-5, etc.), and then reversing the order of the digits. This sequence (figure 5.4c) fills the design space without clustering.

Many more initial sampling techniques exist, like the well-known Latin hypercube and Sobol approaches. However, the importance of the initial sampling plan is limited in a sequential approach, since most of the points are placed by the adaptive sampler. In the following, full factorial plans with a low number of points per direction (2 to 4) are used as start sets. The Hammersley sequence plays a role in the searching of the acquisition function for the adaptive sampling (section 5.6).

Adaptive sampling with acquisition functions

A sequential surrogate-based optimization is in fact a nesting of two optimization processes. Its goal is to find the optimum of the objective function and, in order to do so, it must first optimize the quality of the surrogate. This section discusses the choice of the sampling strategy and introduces a robust and flexible acquisition function which can create any type of response surface.

Acquisition function: definitions

Given a set of observations D n , if we want to add a new point x n+1 to D n in order to enhance the accuracy of the model, how to select the next query point x n+1 ? The Bayesian approach consists in designing an acquisition function or infill criterion a(x). The acquisition function is an inexpensive function that can be easily evaluated throughout the design space; the next sampling point is selected as the maximum of a:

x n+1 = max x∈X a(x).
(5.3) Thus, the optimal distribution of points in D n is obtained by choosing a correctly. In general, two sampling strategies can be identified (section 5.3.3): exploration which is the sampling of points x where f (x) has a high uncertainty (typically, this means areas with few existing points), and exploitation which is sampling around regions of interest in f (typically, its minimum).

A universal formulation for a is introduced here that is able to perform both response surface creation for a VPP and design optimization. This formulation requires that an estimator U (x; D n ) for the uncertainty in f (x; D n ) is available, and that the region of 5.5 Adaptive sampling with acquisition functions interest for exploitation is the minimum of f as in equation ( 5.1). Then, the proposed acquisition function is given by:

a(x; D n ) = U (x; D n ) -β f (x; D n ), (5.4)
where β ≥ 0 is a trade-off parameter to balance exploration and exploitation. If the surrogate model is to be used as a response surface, exploration is always to be prefered so β is set to zero.

For shape optimization on the other hand, a fixed non-zero value is chosen for β. Thus, in the beginning of the optimization U is high so a will be dominated by U , favoring exploration. Later on, when more points are added and the overall shape of f is better known, U will diminish and a becomes dominated by f . As a result, the acquisition function automatically switches to exploitation.

Existing acquisition functions

Many different acquisition strategies exist. To provide a point of comparison with the custom acquisition function described below, two common strategies are described here.

Variance The Gaussian process regression method provides at any point the statistical prediction error, or variance σ(x) = var(f (x)), see equation 4.12. This variable can be used to estimate the uncertainty and can therefore be used as an acquisition function in a purely exploration-type response surface creation. In the formulation (5.4), this choice corresponds to setting U = σ and β = 0. However, the GP variance in a point x depends mostly on the distance to other points; the actual shape of the function f does not have a major influence on σ(x), see section 4.3. Thus, this choice will lead to an equidistribution of the points in space, similar to a full factorial sampling.

Lower Confidence Bound In an optimization context, the acquisition function obtained by choosing U = σ and β > 0 is known as the Lower Confidence Bound approach (LCB) [START_REF] Cox | A statistical method for global optimization[END_REF][START_REF] Cox | SDO: A statistical method for global optimization[END_REF]. The idea is that in the beginning, the variance dominates so an exploration strategy is chosen similar to the variance-based sampling described above. As the metamodel improves, the variance diminishes so the minimum of f starts to attract the new points, which means that the region around the optimum is sampled more than the rest of the parameter domain. The difference with the variancebased sampling is illustrated in figure 5.5.

Custom acquisition function

To compute the estimated uncertainty U in equation (5.4), as an alternative to the variance which mainly depends on the distance between sampling points, it would be preferable to take into account an actual error estimation for f . One way to estimate this error is to perform cross validation. Leave one out cross-validation (LOO-CV) computes the model errors in the points x i of D n , by comparing f (x i ) with the prediction of a GPR based on all the points in D n except i itself:

e i = y i -f (x i ; D n\i ), (5.5)
For the dataset E n = {(x i , e i )} n i=1 , a surrogate error model is then constructed:

ê(x; E n ) = A(E n ). (5.6)
To concentrate the sampled points in the regions where the error in f is high, the variance is weighted with this error model:

U (x; D n ) = σ(x)ê(x; E n ).
To obtain the final form for the uncertainty estimator, this expression is further modified.

First, the error model (5.6) is not perfect, especially for few sample points: the actual error may be high in zones with no sample points and the model cannot see this. Therefore, it may be necessary to sample large regions without any points, which is obtained by including the unweighted variance in the uncertainty estimator, based on a balance parameter. Furthermore, the cross-validation error is divided by the constant σ e which is the average value of the variance in the points i for the evaluation of (5.5):

U (x; D n ) = σ(x) α ê(x; E n ) σ e + (1 -α) .
(5.7)

Thus, α = 1 produces pure cross-validation based estimations, while α = 0 recovers the variance-based or LCB sampling.

Adaptive sampling with acquisition functions

One of the difficulties in using this criterion lies in the choice of the value assigned to the parameter β. This is the reason for dividing by σ e , which has the effect of normalising σ(x).

In fact, the magnitude of σ relative to f is unknown but ê is an error estimator for f so it has the same dimension. Thus, since σ(x) σe is of O(1), the custom criterion with α = 1 makes sense when β is chosen around 1. The same reasoning does not apply to LCB sampling, where β has to be determined by trial and error. This aspect is illustrated in sections 5.8.1 and 7.1.3.

Searching the acquisition function maximum

For adaptive sampling, selecting a new point to sample involves determining the maximum over the acquisition function according to equation (5.3). This search is simple on one hand, because the acquisition function is cheap to evaluate. On the other hand, the function is irregular and has many local maxima, which complicates the optimization (figure 5.7). This section presents two classical optimization techniques and shows how they have been combined to search the maximum of the acquisition function.

Generally speaking, all search strategies can be classified either as complete (and therefore global) or local. The difference between them is that complete strategies perform a systematic examination of all possible solutions in the search space whereas the other strategies only concentrate on a part of the solutions, following a known algorithm.

Gradient based algorithms This is a class of algorithms that trace a single path from a given starting point, following the direction of the local gradient everywhere on the path. By following the gradient, one makes sure that the function value is monotonically increasing along the path, which implies that the path ends up at a maximum.

The interest of gradient based search algorithms is that they can be very efficient, requiring only a few function evaluations to reach the maximum. Their main disadvantage is that they provide no guarantee whatsoever that the maximum found is global; depending on the starting point, the search may end up in a suboptimal local maximum.

Brute-force methods or grid search Typical acquisition functions are highly nonlinear with many local optima (figure 5.7), so there is a high risk of being trapped in a local optimum for a non robust search algorithm. Therefore, an exploration phase is essential.

The search algorithm needs to explore extensively the entire design space in order to find the global optimum.

This type of search can be performed with different algorithms. Evolutionary algorithms for example are known to be robust optimizers that are well suited for discontinuous and irregular objective functions. However, the surest way of exploring the design space is to evaluate the acquisition function in a fixed grid of points, similar to an initial sample plan (section 5.4). Full factorial sampling plans can be used, but even with a function that is cheap to evaluate, these samplings become too expensive when the dimension of the design space increases.

For this reason, the Hammersley sequence is used in our case. Independent of the design space dimension, the number of points in the grid is fixed at 10000.

Hybrid search method While the grid search method is sure to come close to the global maximum, it is not very precise. Therefore, the chosen method for searching the acquisition function space is a two step method, which combines a grid search with a gradient method. A Hammersley point set is generated and the acquisition function is evaluated. Then, the best n points are selected and used as start points for the gradient search method. In the following, n = 10 is chosen. Performing the optimization this way is safe and accurate but computationally expensive, and is only possible because the acquisition function is cheap to evaluate.

Exclusion zones and constrained optimization

For some optimizations, designs in a part of the design space X may be unacceptable. This can be the case if the design is subject to constraints coming from structural considerations, rules and regulations, etc. Another possibility is that the flow solver used to evaluate f does not return a solution for some points. These points must be excluded from the optimization.

Constrained optimization: definition

Bound constrained optimization problems consider the problem of optimizing an objective function subject to constraints on the design space. These can be explicit specifications of the admissible part of the domain:

minimize f (x), subjected to x ∈ X c ⊂ X.
(5.8) Figure 5.8 shows such a type of constraint. A practical example is a hydrofoil for which a minimum thickness is imposed, to satisfy structural requirements. Many optimization problems also have constraints which are unknown a priori. These constraints require the evaluation of a separate constraint function g(x) and the optimization takes the form: minimize f (x), subjected to g(x) > 0.

(5.9)

Often in fluid dynamics, computing the flow in a design point x will give both f (x) and g(x). For example, a foil section can be optimized for drag with a minimum lift coefficient of 0.6, both indicated by the simulation. C l > 0.6 is the constraint function here.

Exclusion zones for constrained optimization

Constrained optimization with adaptive sampling is handled through exclusion zones, which are the parts of the design space where sampling is not authorized. The adaptive sampler is forced to choose new points outside of these exclusion zones. For explicit constraints (equation 5.8) respecting the exclusion zone is simple: the points outside X c are not searched during the maximization of the acquisition function.

For implicit constraints (equation (5.9)) the constraint function g is modeled by a GP regression like the objective function f . This GP model is evaluated in all the points tested during the optimization of the acquisition function and the parts of the design space not satisfying the constraint are assigned an arbitrary low value for their acquisition function.

As a consequence, the sampler is not picking points in those zones. This method is known as a penalty function.

Instead of creating a supplementary GP model for the constraint function, one could assign an arbitrary high value to the objective function in points that do not satisfy the constraint and, in this way, include the constraint in the objective function. The disadvantage of this approach is that the model for f becomes discontinuous, which reduces the quality of the GP fit and may lead to erroneous optima. Therefore, a separate metamodel for the constraint function is prefered.

Exclusion zones for unconverged simulations

For certain reasons, the simulation code included into the optimization loop may not return a usable result. This can happen when the simulation ends up unconverged for instance, which is the case with the X-foil panel code which returns an empty array for lift and drag coefficient if no solution is found. For FINE™/Marine in similar situations, the flow may not converge to a steady state. A solution needs to be found to deal with these scenarios.

Two solutions were initially envisaged: assigning an arbitrary value to the considered simulation point or removing the point from the dataset D n used to build the surrogate. None of these were satisfying. Removing the point from the data set implies that the acquisition 5.7 Exclusion zones and constrained optimization function remains the same as before the point was selected. This leads to sampling the exact same point at the next iteration and results in performing the same simulation over and over again. Assigning an arbitrary value (an average value was envisaged) would corrupt the surrogate and decrease the convergence rate of the optimization or even worse, could lead the algorithm into an "artificial" local optimum.

The solution finally retained consists in preventing the sampler from exploring any undesired zones by modeling the "no return" event occurrence as a function of input parameters and get a map of go/nogo zones, similar to an implicit constraint. However, the constraint function has only two values, indicating an admissible or inadmissible point. This scenario corresponds to a classification problem that can be found in supervised machine learning.

Gaussian processes can deal with this type of problem. The modeling technique is no longer a Gaussian Process regression but a Gaussian Process classification, which uses a modified probability density function (see [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF], chapter 3. The points returning no response from the simulation are used to model the exclusion zones (figure 5.9) but are excluded from the dataset used for modeling the response f . Therefore, they are called ghost points.

Tests of the adaptive sampling

This section tests the different adaptive sampling approaches presented here on analytical test functions. For optimization, the effect of the objective function weight β in (5.4) is tested and the custom sampler is compared with LCB. Finally, response surface creation is tested. The test functions used are the Three hump camel and Gramacy-Lee functions described in section 4.7.2.

Optimization: effects of the function weight

As noted in section 5.5.1, the parameter β in (5.4) dictates the balance between exploration and exploitation. The effect of this parameter is tested here for the Gramacy-Lee function, using LCB acquisition (α = 0) and a full factorial startset of 5 × 5 points.

Figure 5.10 shows the evolution of the sampled points and the acquisition function for different settings of β, during the adaptive sampling iterative process. It also indicates the convergence of the iterations by showing the distance of each sampled point to the true optimum, in the parameter space. For the highest value of β, to the left, the acquisition functions are dominated by the solution f so the main peak is coming from the single trough in the Gramacy-Lee function. As β is reduced, the uncertainty estimation (here, the variance σ) gains in importance so stronger local maxima of the acquisition function appear in between all the sampled points.

This has the effect of spreading the points over a wider area. For the highest setting of β the added points form a nearly straight line towards the optimum and the convergence is rapid, as shown by the graph of the distance to the optimum. The middle setting of β results in a wider spread of the points and a somewhat slower convergence, although all points remain concentrated around the optimum.

The lowest value for β produces first a series of points around the optimum, followed by an exploration point in the tenth iteration. This shows that the method does not necessarily start with exploration and then change to exploitation. On the contrary, the maximum implied by the startset attracts some exploitation points first, but when the method detects that this maximum is represented sufficiently well, it switches back to exploitation. This is actually a procedure which is both efficient and safe, since it gives a good indication of the optimum early on, while still making sure through exploration that this optimum is not local. Optimizations for actual RANS simulations follow the same pattern, see for example section 7.3.3.

Thus, the highest values for β lead to fast convergence, while the lower values are safer; the right choice depends on the function being optimized. Since it is not easy for LCB to give guidelines for β, the values here were chosen through trial and error. In section 7.1.3, this problem is revisited for RANS simulations. In an optimization setting, the two acquisition functions have a similar behaviour, since the main action of the acquisition process is determined by the balance between the uncertainty U and the function f in equation (5.4). Whenever the uncertainty is reduced enough, the points will be placed around the optimum of the function, independent of the exact way in which the uncertainty is estimated. This is confirmed by figure 5.11 which shows the optimization for the Three hump camel. Both acquisition functions produce similar placements of the points and have the same type of balance between exploration and exploitation. Points are placed alternately close to the optimum and further away. The true optimum, indicated by a star in the figures, is found by both acquisition functions within 10 iterations.

However, in certain situations the LCB acquisition may be unable to reach the optimum. This is the case for example in figure 5.12 which shows the optimization for Gramacy-Lee. For this optimization, by coincidence, one of the points in the initial sample plan is placed close to the maximum which lies right next to the minimum for Gramacy-Lee. The high value of f in this point reduces the acquisition function around it (this produces the large blue spot in the left images of figure 5.12), which implies that the entire domain is sampled except the lower left corner. The minimum, which lies in this corner, is therefore not found. This same problem was observed during tests with RANS simulation, see section 7.1.3.

The custom function does not suffer from the same problem. Since the point near the maximum is the only initial point where f is nonzero, removing this point will alter the response surface completely. Thus, the cross-validation error in this point is high, which increases the weight on the variance in the lower left corner (two brighter spots can be seen around the blue spot in the top right image of figure 5.12). Thus, this corner attracts sampling points, the minimum is detected in the seventh iteration, and the optimization succeeds from there on.

Therefore, although the uncertainty estimators are not very critical for an optimization acquisition function, the custom sampler is still a safer choice since it allows to find the optimum even in certain pathological situations.

Tests of the adaptive sampling

Iteration 10

Distance to opt. For this final test, the custom and LCB acquisition functions are used to construct response surfaces for the Three hump camel and Gramacy-Lee functions. To obtain response surfaces, the parameter β is fixed at zero. The initial sample plan is a 4 × 4 full factorial startset.

The main criteria for evaluating the acquisition functions is the statistical metric R-squared R 2 of equation (4.21), evaluated on a regular 30 × 30 grid over the design space. R-squared with a value of 1 indicates that the model perfectly fits the data.

On the Three hump camel test function (figure 5.13) the performances of the two acquisition functions are similar with the variance-based acquisition function performing slightly better.

As indicated, the variance sampling creates a regular full factorial-like distribution which is well suited for this smooth test function. Thus, the R-squared metric converges faster for the variance-based sampling and it attains a higher final value. However, the quality for both surfaces is good.

On the Gramacy-Lee test function (figure 5.14) the custom acquisition function is giving notably better results than the variance-based acquisition function. The variance leads to a uniform sampling of the design space which is nearly identical to the sampling for the Three hump camel. This proves that for response surface creation, the variance-based sampling is not really an adaptive technique since it does not react to the function being reconstructed. The custom acquisition function on the other hand samples points in areas with large variability, close to the peaks of the test function. Therefore, the R-squared quality of the custom surface is higher.

Thus, the custom acquisition function is the safest choice: it performs almost as well as the variance for smooth problems and remains effective when the function has local peaks.

Conclusion

In the maritime industry today, simulation-based optimization is gaining interest. The purpose of this approach is to search for optimal geometries automatically, using numerical simulation to evaluate the performance of the geometries considered. The main limitation of simulation-based design is that performing high-fidelity simulations using nonlinear physical models, for the large number of design points required to search a multidimensional design space, is prohibitively expensive.

A possible solution to this problem comes from surrogate-based optimization. In this technique, the results for a limited number of expensive simulations are used to construct a metamodel of the simulation response, which is inexpensive to evaluate. The geometric optimization is then performed over the metamodel. In this work, Gaussian process regression is adopted to construct the metamodels.

A key aspect of simulation-based surrogates is that the design points to be simulated can be chosen. In a sequential or adaptive sampling approach, an initial metamodel is constructed from a limited number of points. Then, in an iterative process, the design point is searched where this metamodel is most lacking in quality. This point is simulated and added to the dataset, after which a new metamodel is constructed and the iteration recommences.

This approach implies a double optimization: the metamodel is used to optimize the geometry, but first the quality of the metamodel is optimized by adding points through the adaptive sampling. The choice of the points to sample is made by defining an acquisition function, different from the function to optimize, and to place each new point in the maximum of the acquisition function. The acquisition proposed here combines the function to optimize with an estimation of the uncertainty in each simulation point. Thus, points will be distributed throughout the design space to reduce the uncertainty, the so-called exploration. When the uncertainty is low enough, the new points will be concentrated around the optimum of the function, which is called exploitation.

The uncertainty in the metamodel is computed using cross-validation: the simulated value in each design point is compared with the value produced by a metamodel from which that data point has been removed. The variance of the Gaussian process is weighted with this cross-validation error to produce the final uncertainty estimate.

The difficulty in the search for the maximum of the acquisition function is that this function has many local maxima. For this reason, a hybrid search strategy is chosen: first, the entire design space is sampled in a brute-force approach where the acquisition function is computed in the points of a Hammersley distribution. Starting from the best points found, a gradient search is then initiated. This brute-force approach is acceptable since the acquisition function is inexpensive to compute.

In some situations, it may be necessary to exclude a part of the design space, for example in the case of constrained optimization or if the simulations fail for certain design parameters.

In this case, a second metamodel is constructed whose value indicates the exclusion zone: in this zone, the response of the acquisition function is set to a low value so these points are not selected.

Tests of the acquisition function show that a high weight on the function, with respect to the uncertainty, leads to a rapid convergence of the optimization. A lower weight reduces the convergence speed but leads to more exploration, which minimizes the risk of missing the global optimum. With respect to a pure variance-based uncertainty estimation, the custom acquisition function performs better for the creation of response surfaces with local

Geometric modeling

One of the core issues that an engineer must deal with when designing lifting surfaces in a computer environment, is how to describe their shape. It is impossible to treat shapes on a point by point basis, since this requires an infinite number of degrees of freedom. The typical practice is to rely on a series of curves, such as polynomials or Bezier curves, to describe the surfaces. Thus, each surface is defined by a few parameters or control points only, which reduces the total degrees of freedom to a manageable number.

In an optimization context, two main approaches can be distinguished. Deformation-based modeling uses a single base shape for the geometry, which is deformed to create a family of shapes. This deformation is parametrized with control points or deformation parameters.

Constructive modeling on the other hand, uses parameter-based curves to construct each geometry in the family from scratch. Thus, the shapes are generated independently from each other.

Within this thesis, a geometric modeler was created to generate hydrofoil shapes. The design brief for this tool was, that it had to generate the full three-dimensional shape of the hydrofoil and to be versatile enough for use throughout the foil design process, from pre-design to manufacturing. This means that the generated shapes and the data related to the foil have to be suitable for numerical simulation, for analysis by naval architects, and for CNC mold milling. These needs originate from the requirements of Groupama Sailing Team, one of Streamline's clients.

The modeler uses a mix of constructive and deformation-based methods. It is able to generate multiple design output formats like IGES, STEP, Parasolid, or structured point clouds. Furthermore, it can be used as a standalone tool or integrated in an automatic optimization loop. This chapter shows the principles of the modeler and presents two projects in which it has been applied.

Principle of the shape generation

In the modeler, the geometric design of the hydrofoil is achieved in several distinct steps (see figure 6.1). The first step, described in section 6.2, consists in defining the general outline of the foil by creating a geometric object called the spine, which serves as the centerline for the airfoil shapes (figure 6.2). This step is constructive in nature. The spine is then divided in sections of equal size and a foil profile is placed on each section.

The second construction step is the generation of these profile sections; this step is itself divided in two. First, smoothly varying profile shapes are created by blending one or more basic profiles (section 6.4). Then, these shapes are deformed with the possibility of varying several parameters such as chord, camber and thickness (section 6.5). The final step is to position the required number of sections along the spine with a regular curvilinear spacing (section 6.6).

The result of these steps is a structured point cloud which defines the hydrofoil surface. If desired, the foil envelope can be generated by fitting a NURBS surface through the cloud or by triangulating between the points. This step is described in section 6.8. is a circle arc with high radius. TEL or Tip Equivalent Length is the length of this tip element. On the right side is a graphical illustration of the tangency continuity used to create the foil spine. Arc 1 and Arc 2 are joined at tangency point A.

Spine generation

The spine of the foil is the three-dimensional curve used to support the hydrofoil sections. This spine can either be parametrized within the modeler, or imported as a point cloud with a fixed shape. For example, fixed spine shapes were used in the design process of the hydrofoils fitted on the USV co-developed by Streamline and Centrale Nantes Hydroproject (see section 6.10.2).

The parametrization option for the spine allows to optimize its shape when the modeler is included in an optimization loop. In this mode, the spine is created from basic geometric elements combined together with continuity requirements. The basic spine type used is composed of four circle arcs (see figure 6.3), which are defined by their radius and chord length. The circle arcs are joined together based on a tangency constraint. The tip of the foil is parametrized with an additional parameter which is the angle formed by the chord and the reference horizontal axis. Furthermore, one reference point is necessary to fix the position of the entire spine.

The choice of the generation technique was made because its parameters have a physical meaning and are intuitively clear for a designer. Such architectural parameters allow the users to better understand the effect of modifying them, and give meaning to the results of an automatic optimization. This aspect of the hydrofoil modeler is further discussed in section 7.3.

Spine generation

Once the parameters are known, each primary element is discretized into a number of points. The sets of points from each arc are then gathered into a common set and reinterpolated with a single B-spline. This interpolation makes it possible to work with a continuous spine and it serves to define a normalized curvilinear coordinate system (figure 6.4) with coordinate t. The origin for this coordinate is located on the head of the foil, the tip of the foil having an abscissa value of t = 1.0. Thus, the spine has two definitions:

• The one-dimensional curvilinear coordinate system t ∈ [0, 1],

• The 3D space coordinate system, where the spine is defined as x S (t).

The curvilinear coordinate is used as a reference to generate the airfoil sections, which are placed on the 3D spine to create the shape of the hydrofoil. This is described in the following sections.

Parameter distribution laws for the sections

The spine gives the foil outline but its envelope originates from the foil sections. In general, the shape of these sections is not constant; the profile chord, twist and shape vary according to the position on the spine. These parameters do not vary independently from section to section but are related by parameters distribution laws to insure the smoothness of the final shape. Those distribution laws are continuous functions that extend over the entire length of the spine and set the value of the considered parameter at each curvilinear abscissa t. Distribution laws are obtained by defining control points corresponding to specific zones of the spine. Those control points are then interpolated in order to obtain a continuous description of the parameter (see figure 6.5 for an example). To interpolate through the control points, the spline method is prefered over polynomial interpolation. Instead of a single high-order polynomial, the spline method uses piecewise polynomials with lower degrees, thus avoiding the Runge phenomenon (the non-convergence of the interpolation polynomials towards the function, due to oscillations). However, interpolation by cubic splines can also present waves of high amplitude (overshooting) which makes it inappropriate for the interpolation of section parameters (see figure 6.6). To overcome this problem, the interpolations are performed using cubic Hermite splines (Piecewise Cubic Hermite Interpolating Polynomial or PCHIP), which are similar to cubic B-splines but which have been developed to eliminate overshooting [START_REF] Fritsch | Monotone piecewise cubic interpolation[END_REF]. Mathematically, a cubic Hermite spline is a spline of which each piece is a Hermite polynomial, which is itself a linear combination of 4 unitary polynomials of degree three.

Parameter distribution laws for the sections

A Hermite polynomial allows the function and its first derivatives to be imposed at the control points. This makes it possible to ensure not only the continuity of the derivative at the control points, but also the monotonicity of the function between the points.

Basic profiles

To create the profiles to be placed at each section of the spine, it would be possible to parametrize the profile shape and to vary the parameters using distribution laws as defined above. However, the 2D hydrodynamic profiles used to construct a hydrofoil are often designed separately and in this case, they are mostly defined by point clouds and not by parametric variables. Thus, the characteristics of the profile are unknown and it is impossible to act directly on them to build the hydrofoil. Therefore, a deformation technique is chosen, where the hydrofoil shape comes from fixed profiles, with a variable transformation applied to adjust for example the scaling or camber. This transformation is described in section 6.5; the current section concerns the definition of the basic profiles.

Sometimes it may be useful to use different basic profiles at different areas of the foil. When more than one profile type is used, they are each assigned a separate zone of the foil and between these zones, a smooth transition is necessary. This is performed with the blending process presented in section 6.4.2.

NACA four-digit sections

While the basic profile shapes usually come from point-cloud libraries, an alternative is available in the geometric modeler: the NACA four-digit series. This series of profiles is parametrized, but that option is not used to vary the shape. Instead, constant choices for the parameters are applied to create fixed profiles similar to point cloud-based ones and the resulting profiles are then deformed as described later.

The four-digit series of airfoil shapes was developed by the NACA (National Advisory Committee for Aeronautics). Their creation is fairly simple and consist in generating a symmetric profile and a camber line using analytic expressions. The camber line is used to transform the symmetric profile into a cambered profile [START_REF] Abbott | Theory of Wing Sections, Including a Summary of Airfoil Data[END_REF].

NACA four-digit sections are parametrized by three non-dimensional variables:

• maximum camber C,

• longitudinal position of the maximum camber X,

• maximum thickness T ,

80

Chapter 6 Geometric modeling each given as a percentage of the chord length c. The half thickness of the symmetric profile is given by:

t c = 5T 0.2969 x c -0.1260 x c -0.3516 x c 2 + 0.2843 x c 3 -0.1015 x c 4 , (6.1)
where x is the horizontal coordinate. The mean camber line is divided in two parts and given by:

z c c =    C X 2 2X x c -x c 2 0 ≤ x c ≤ X, C (1-X) 2 (1 -2X) + 2X x c -x c 2 X < x c ≤ 1, (6.2)
To find the upper and lower sides of the airfoil, the half-thickness is applied normal to this mean camber line. The resulting profile has a trailing-edge thickness of 0.021cT ; if a zero-thickness trailing edge is desired, the profiles closest to the original ones are obtained by changing the last coefficient to 0.1036.

Maximum Camber

Maximum Thickness position of maximum camber Figure 6.7: NACA four digits foil sections. The top drawing is the symmetric section obtained with equation 6.1. The bottom drawing is the mean camber line.

Blending hydrofoil profiles of different type

By the blending of profile shapes, different existing 2D profiles or shapes can be used at different locations along the spine. All foil sections will be defined by 2D point clouds having m points each. Let the foil be constructed with N basic shapes. Once loaded from the database or generated with the 6.4 Basic profiles NACA4 modeler, the N sections are re-interpolated with a B-spline and a new point cloud is generated with an identical distribution of points along the x-axis for all section types. This re-interpolation ensures that the x-position of the section points is similar for the different sections to blend, which allows to perform the blending only for the z-position of the points. The distribution along the x-axis is adapted to have a concentration of points close to the leading edge so this zone is described accurately. The result of the re-interpolation is the structured point cloud:

x K i 1 ≤ K ≤ N 1 ≤ i ≤ m, ( 6.3)
which defines each basic shape. For the moment, the sections are vertical so y K i = 0 everywhere.

To create the n shapes which will be placed along the spine, a simple weighted average is performed on the coordinates of the base sections:

x i,j = N K=1 ψ K (t j ) x K i 1 ≤ i ≤ m 1 ≤ j ≤ n. (6.4)
The weights ψ K are obtained by performing a partition of unity, such that N K=1 ψ K (t) = 1 for every curvilinear coordinate t.

The basis functions ψ K are obtained by interpolating control points with PCHIP functions to ensure a smooth transition between shapes. The position of the transition is obtained by translating the set of control points along the t-axis. The relative transition length can be adjusted by moving the control points closer or further away from each other (see the example in figure 6.9). .9: The partition of unity used for section blending in the case of two basis functions. The grey area corresponds to the transition length where the two base profile types are mixed. Outside of this grey area, the profiles are not mixed and therefore kept in their original form.

82

Chapter 6 Geometric modeling

Section modification

Before the sections are rotated and placed on the spine, their shape is adjusted. These modifications concern the camber, the chord, and the trailing edge thickness.

Camber modification of 2D profiles

The camber of the basic profiles coming from equation ( 6.4) can be increased through conformal mapping. A conformal map is function that preserves the angles locally, which means for instance that orthogonal vectors stay orthogonal. The mapping used here keeps vertical lines straight, while transforming horizontal lines into circle segments (see figure 6.10). The amount of added camber depends on the magnitude of the transformation.

A separate procedure was developed to extract the thickness distribution and the mean camber line from an arbitrary profile. The intention was, to iteratively adjust the amount of conformal deformation until this extracted mean camber line had the desired maximum camber. However, this coupling between the two procedures was never completed.

Camber + Initial Camber Figure 6.10: Comformal mapping. The left figure shows a conformal mapping applied on a regular grid. On the right figure, the same type of transformation is applied to a hydrofoil section to increase its camber.

Scaling

Once the base profile mix is performed and the camber is adjusted, the sections are scaled to match the required chord length distribution. The scaling is performed by applying an affine scaling transformation. All the following transformations are represented by matrices, including translation (see 6.6) thanks to the use of homogeneous coordinates1 .

Chord and thickness The transformation to obtain the desired chord distribution is the following scaling:

       x y z w        i,j =        α x (t j ) 0 0 0 0 1 0 0 0 0 α z (t j ) 0 0 0 0 1               x y z w        i,j , (6.5)
where α x (t j ) and α z (t j ) represent the scaling factors in x and z-direction, respectively. If only the chord length is to be adjusted, the factors are chosen the same so the scaling is uniform. By choosing different factors for x and z, the profile thickness is also adjusted.

Trailing edge thickness Many base profile shapes do not include a blunt trailing edge. However, sharp trailing edges are difficult to manufacture and too fragile for practical use. Therefore, hydrofoils are constructed with blunt trailing edges. For the current modeler, it was desired to have a trailing edge thickness that is uniform over the hydrofoil, independent of the local chord. Therefore, a sub-optimization is performed at each section: the profiles are enlarged adaptively until the trailing edge thickness at the true chord length position equals the desired value. Then the remaining part of the profile is cut off.

Spatial positioning of the sections

After adjusting the camber and chord, the sections have their final shape. To finish the hydrofoil construction, the sections are placed and oriented on the spine by using a series of affine transformations (rotation and translation). The following operations are performed on the sections to position them correctly in the 3D space:

• Rotation around the y-axis to apply twist,

• Rotation around the x-axis to be in the plane normal to the spine,

• Translation to position the section along the spine.

Twist The application of twist is included in the section positioning because this operation does not change the shape of the profiles, it is only a rotation around the axis normal to the profile. The twist distribution is given by a function θ(t). The rotation matrix is: .6) direction of the spine. The angle for this rotation around the x-axis is defined by the shape of the spine:

       x y z w        i,j =        cos θ(t j ) 0 -sin θ(t j ) 0 0 1 0 0 sin θ(t j ) 0 cos θ(t j ) 0 0 0 0 1               x y z w        i,j . ( 6
φ j = arctan ∂y S ∂x S j . (6.7)
The resulting transformation matrix is:

       x y z w        i,j =        1 0 0 0 0 cos φ j -sin φ j 0 0 sin φ j cos φ j 0 0 0 0 1               x y z w        i,j
. (6.8)

Translation Finally, each section is translated to its correct position on the spine. This translation is performed such that the reference point of the section j moves to the coordinate x S (t j ) on the spine. Thanks to the homogeneous coordinates, this translation is treated as a final matrix multiplication:

       x y z w        i,j =        1 0 0 x S (t j ) 0 1 0 y S (t j ) 0 0 1 z S (t j ) 0 0 0 1               x y z w        i,j
. (6.9)

Attitude setting

Before generating the hydrofoil envelope from the point cloud, the foil position relative to the boat reference frame can be modified. This feature allows to investigate the influence of the hydrofoil positioning. Two parameters can be adjusted (figure 6.11): the foil extension and the cant angle, which is the foil rotation around the x-axis.

These motions mimic the way a real hydrofoil is attached to a sailing boat, i.e. by two sets of bearings. The kinematic constraints induced by these bearings are taken into account: the cant angle is changed by moving the upper bearings while the lower ones remain in place and the extension is varied with both bearings fixed in position. Thus, the cant angle and vertical extension reproduce the possible settings of the platform the hydrofoil is fitted to. Both movements require an iterative procedure to compute the position of the hydrofoil.

Surface generation

After transforming the sections and placing them at their corresponding position, the foil exists in a discretized form, described by a structured cloud of points. The final step in the 6.7 Attitude setting geometry generation consists in creating the envelope of the hydrofoil, either in continuous form or as a triangulated mesh (discrete form).

CAD format export using the Open Cascade library

The foil envelope can be generated in continuous form by fitting NURBS surfaces between the points. This fitting operation is performed with the Open CASCADE Technology2 (OCCT).

OCCT is an open source 3D modeling C++ library used by the modeler as a CAD kernel.

OCCT is used together with pythonOCC3 , an open source software framework for 3D modeling based on the Python language. The core of pythonOCC is a Python wrapper, performing a binding of the OCCT library written in C++ to the Python language.

PythonOCC offers the possibility to perform advanced geometric and topological operations as well as data export (STEP, IGES, STL import/export). The use of pythonOCC allows for a greater portability of the modeler and the possibility to fully automate the generation of 3D foil shapes. This automation is essential in order to integrate the modeler into an optimization loop. In addition, by using pythonOCC and OCCT the modeler only depends on elements under free licenses.

Triangulated mesh generation

The geometry definition used by HEXPRESS, the mesh generator of FINE™/Marine, is based on triangulated surfaces. A custom Python tool has been written to create this triangulated mesh directly from the hydrofoil point cloud.

The triangulation of the surface is made easy by the structured nature of the point cloud describing the hydrofoil. Since each section contains the same number of points, the surface between two sections can be considered as a series of rectangles formed by two subsequent points in both of the sections. A triangulation of this surface is obtained by dividing each rectangle in two triangles. This triangulated mesh is saved as a ".dom" file and can be directly used by HEXPRESS to create a volumic mesh.

This approach has two major advantages. First, the geometry description used for the simulations is identical to the shape generated by the modeler, so the geometry is exact.

And second, the procedure ensures that the topology of the hydrofoil is perfect, without any holes in the surface. Since this is a requirement for HEXPRESS, the procedure guarantees that the geometry is suitable for meshing, which is important for the robustness of an automatic optimization loop.

Foil data display

To facilitate the analysis by naval architects once a hydrofoil geometry is created, the main data elements regarding the foil characteristics are gathered and provided to the user through graphical elements or foil "identity cards" (figure 6.12). Other parameters describing the spine geometry are presented in a table. The full set of these output parameters is presented in figure 6.13. This identity card was requested by Groupama Sailing Team to improve the efficiency of the discussions during design meetings. More information about the context of this collaboration can be found in the next section.

Industrial applications of the 3D modeler

The hydrofoil modeler has been applied in different industrial and educational projects. This section presents two of these projects.

6.10.1 Hydrofoil design for the Groupama Sailing Team C-Class Groupama Sailing Team is a professional sailing team created in 1998 by Franck Cammas. The team has been involved in different races and sailing classes, like the Route du Rhum and the Trophée Jules Verne. After winning the Volvo Ocean Race in 2011/2012, the team decided to get involved in the America's Cup. Groupama Team France was founded in 2013 for that purpose and was largely made up of the existing Groupama Sailing Team. With the intention to prepare for the America's Cup, Groupama Sailing Team decided to participate in the 2013 International C-Class Catamaran Championship. This competition, also dubbed "Little America's Cup" is raced on C-Class catamarans. The C-Class is a highperformance development class of sailing catamarans which are only constrained by a maximum length, beam and sail area. The resulting boats are very light and can sail at twice the speed of the wind. They are allowed to use rigid wing sails and to fly, which explains their high performance.

At the end of 2014, Streamline was commissioned by Groupama Sailing Team to take part in the design process of the main foils for the C-Class Groupama C, shown in figure 6.14.

Industrial applications of the 3D modeler

For this design, the foil modeler was not included in an automatic optimization loop but used to generate systematic families of foils. Within those families only a single parameter was varied (figure 6.15). The purpose was to evaluate and understand the influence of the geometric parameters on the performance criteria, which were the heave stability and the drag. The hydrofoil section was taken from the 2013 hydrofoil designed through the collaboration of naval architects Guillaume Verdier and Martin Fischer. A total number of 110 foils were generated and evaluated during this campaign. The foil geometry eventually chosen was the one having the best heave stability and minimizing the drag, while matching the required F y and F z to lift the boat out of the water and counteract the side force of the sails (section 2.3.2). An example of a vertical stability analysis is shown in figure 6.16. Vertical stability is characterized by the computation of the derivative of the vertical force generated by the hydrofoil with respect to the flying altitude: S z = ∂F z ∂z .

The performances were evaluated using FOILWIZARD, a software suite created by Streamline to design and analyze lifting hydrofoils. FOILWIZARD combines the two-dimensional panel In September 2015, Franck Cammas and his crew Louis Viat won the 27 th edition of Little Cup with their C-Class catamaran fitted with a hydrofoil designed and studied by Streamline and whose 3D-shape was generated using the hydrofoil modeler presented in this chapter.

Centrale Nantes Hydrocontest USV

Centrale Nantes Hydroproject is a pedagogical project started by a teaching team of Centrale Nantes to train engineering students with practical design challenges. The students who join the team are required to design, build and test small-scale ships with the goal of participating in the Hydrocontest competition.

The Hydrocontest is an international student competition dedicated to naval energy efficiency. For this contest, students are invited to create the most energy-efficient boat. The Hydrocontest organization provides each team with the same electric engine and battery. Thus, the amount of available energy is imposed to travel a set distance and the ranking is made according to the best running time.

Teams compete in two categories: the Mass Transport category in which each vessel must hold 200 kg of ballast to simulate the displacement of a cargo vessel, and the Private All the hydrofoil geometries were generated using the 3D modeler.

Over the length of the hydrofoil, three separate basic sections were used:

• A squared top section.

• An ogival section in the intent of minimizing wave and spray induced by the foil/free surface interaction.

• The H105 developed by Tom Speers4 .

The latter two sections are shown in figure 6.18. The ogival section was generated using the CST method [START_REF] Kulfan | A universal parametric geometry representation method -CST[END_REF] implemented as a 2D profile modeler similar to the NACA four-digit generator. The blending of the sections is performed as in section 6.4.2. The boat is under continuous development and is also available as a test platform for researchers of Centrale Nantes and for industrial partners. The team won the competition in 2016 in the lightweight category and will compete again in 2018.

Conclusion

The challenge in designing a geometric modeler is to create geometries based on as few parameters as possible, while retaining the possibility to generate a wide range of designs and providing intuitive feedback to designers about the effects of each parameter. Two basic approaches were identified: deformation modeling based on the modification of a base geometry, and constructive modeling using architectural parameters. The former has the advantage of being generally applicable; the latter provides parameters which are natural for a designer, but it is limited to a smaller class of geometries.

Ultimately, a hybrid of the two approaches proved to be the best choice. In the modeler presented here, the hydrofoils are generated in two parts. First, the spine is created which specifies the outline of the foil. This spine is defined using architectural parameters such as segment radii and lengths. Then, profile sections are placed along this spine to form the hydrofoil geometry. Since airfoil shapes are often given in terms of point clouds instead of parametric expressions, it was preferable to create these profile sections by the deformation of basic airfoil shapes, rather than by direct parametric generation. Due to the way the deformation is set up, the controlling parameters such as the chord, camber and twist distribution still have an intuitive meaning.

The export to various CAD formats is an essential part of the modeler. First, it provides a straightforward way of ensuring that the analysis, the simulation, and the manufacturing of the hydrofoils are all performed on the same geometries. Thus, uncertainties about the geometry are eliminated. And furthermore, the modeler produces domains for the fluid simulation which are topologically perfect. This guarantees the reliability of the simulation chain, which is essential if the modeler and the simulations are combined in an automatic optimization loop.

The development of the hydrofoil modeler was an industrial project, which explains the emphasis that was placed on the exchange of information with the user. The modeler has been extensively tested in an industrial environment and has proven its worth for practical applications.

Conclusion 7 Applications

This chapter combines the developments described in all the previous chapters, in order to perform geometric optimization based on RANS simulations. Section 7.1 analyses the performance of the adaptive sampler from chapter 5 when it is applied to response surface creation and shape optimization of a two-dimensional airfoil. In section 7.2, the procedure is used to optimize two-dimensional airfoils under different conditions and to study the influence of these conditions, such as the distance to the water surface, on the optimum airfoil shape. Finally, in section 7.3 the sampler and the geometric modeler from chapter 6 are combined to perform a geometric optimization of a three-dimensional lifting hydrofoil.

The choice of the geometric parameters for hydrofoil optimizations is discussed, as well as the robustness of the resulting optimum.

Adaptive sampling applied to 2D profiles

This section concerns the creation of response surfaces and the geometric optimization for a two-dimensional lifting airfoil, based on RANS simulation. The tests of section 5.8 are repeated to show how the adaptive sampler behaves for such typical problems in fluid dynamics. The test case is based on the NACA four-digit airfoils.

Test case and simulations

As a simple parametric series of airfoil profiles, the NACA four-digit series is chosen (section 6.4.1). From a modeling point of view, the 4-digit profiles form a parametric series based on three input variables: the maximum camber C, its position X, and the thickness T . Thus, the input parameters represent the physical characteristics of the airfoil, like this is the case for the full three-dimensional hydrofoil modeler from chapter 6. As such, the 4-digit series is suitable as a 2D equivalent of the geometric modeler.

Test case

The ranges allowed for the geometric parameters are [0.03, 0.12] for the thickness, [0.0, 0.08] for the camber, and [0.25, 0.7] for the x-position of maximum camber. The reason for imposing a non-zero minimum T is that the meshing scripts are not adapted to zero-thickness foils and also that such a foil would not make sense from a structural point of view. The minimum and maximum of X have to be limited because the analytic expressions (6.1) and ( 6.2) create cusps in the geometry for too extreme positions of the maximum camber.

The velocity is V = 10m/s, the chord c = 1m, the density ρ = 1026kg/m 3 , and the Reynolds number Re = 8.41 • 10 6 . In this section, all test cases have monofluid flow, but in the following section the profile is placed below a free surface to see how this affects the optimum. The objective for the optimization is the minimization of the drag; to obtain a fair comparison of the different geometries, all simulations are performed at the same lift coefficient C l = 0.6. The angle of attack is adjusted dynamically during the simulations in order to obtain this constant lift coefficient (see below). Simulations and meshes All simulations of the NACA profiles presented here are performed with FINE™/Marine. To fit in the overall optimization framework (section 5.3), the simulations are completely automated: the computational setup is performed using Python scripts, the geometries are created with a simple geometric modeler which implements (6.2) and ( 6.1) and the resulting point cloud is automatically transformed into a computational domain. The mesh is automatically generated in HEXPRESS and a global parameter in the mesh generation scripts indicates the required mesh size, from 'coarse' to 'very fine'. Finally, the simulations with ISIS-CFD (chapter 3) are launched and post-treated by the Python scripts. The computations for the individual geometries are separate FINE™/Marine projects and a new HEXPRESS mesh is generated for each geometry.

In order to simulate each airfoil at a constant lift coefficient of C l = 0.6, the foils are rotated during the simulation. After a given number of time steps, the difference between the actual and the target lift is computed, which is transformed into a required change in angle of attack using the lift slope from 2D profile theory: ∆C l = 2π∆α (see for example [START_REF] Anderson | Fundamentals of Aerodynamics[END_REF]). The airfoil is slowly rotated to its new position and the forces are allowed to stabilize before a new correction is computed. Here, the correction is performed every 200 time steps and the profile is rotated over 100 time steps, for a total simulation duration of 4000 time steps. This procedure is described in more detail for the 3D optimization, see section 7.3.2.

To speed up the optimization process, the simulations are performed using coarse meshes. These meshes have about 15k cells, with 500 faces on the airfoil section. Wall function boundary conditions are applied with an average y+ of about 60. To assess the accuracy of these computations, the NACA 5408 profile is simulated on four different meshes. When free-surface flows are simulated, adaptive grid refinement (section 3.6.2) is used in ISIS-CFD to create the fine mesh around the position of the water surface. The meshes generated with HEXPRESS have no particular refinement around the free surface; this has the advantage that these meshes can be generated without taking into account the relative position of the profile and the surface, which simplifies the Python scripting of the mesh generation. Furthermore, the adaptive refinement ensures that the mesh is always fine at the free-surface location, independent of the exact shape for the airfoil-generated waves. This ensures the quality of the free-surface simulation. Figure 7.1 shows an example of a mesh with free-surface refinement.

Response surface generation

The first test in this section is the creation of a response surface for two input parameters: the thickness T and the maximum camber C. The position of the maximum camber is fixed at X = 0.4. Two series of computations are performed, one with variance-based sampling (α = 0 in equation (5.7)) and one with the custom criterion which combines the variance and the cross-validation error estimation (α = 1). Since the goal here is not an optimization, β in equation (5.4) is set to 0. All series consist of 50 samples, with a 2 × 2 full -factorial startset. the thickness, the response is characterized by a low drag for moderate C, while the drag increases rapidly when the camber is either increased or decreased. For thicker profiles, the drag varies less with the camber. This behavior will be explained in section 7.2.1.

As expected, the variance-based sampling leads to a uniform distribution of the points over the parameter space, while the custom criterion focuses more on the peaks in the drag, especially the onset of the drag increase at T = 0.06 and C = 0. As a result, the variance is higher for the custom criterion, but the maximum cross-validation is lower. Near the drag peaks, the variance-based sampling creates large cross-validation errors, which are absent for the custom criterion due to its denser sampling in these regions. Since the cross-validation is an error estimator for the response surface, this indicates that the custom criterion performs better. A contrasting view is provided by figure 7.3 which shows another response surface, C d as a function of C and X with T = 0.05. The dependence of the drag on X is minimal, so the drag peaks are replaced by high-drag ridges and the response surface is smoother than in the previous example. The regular sampling of the variance criterion performs as well as the custom criterion here, or even slightly better. Still, the custom criterion has the advantage of being able to react to local peaks, should they occur. Therefore, it is a safer choice than the variance-based sampling and it is more likely to work well.

Prediction

Optimization

For the next tests, minimizations are performed on the same two-parameter problems. In equation (5.7), α is set to either 0 (LCB) or 1 (custom sampling), while β in (5.4) is varied. Figure 7.4 shows the response and other parameters for these optimizations. Compared with the response surface creation, the points are more concentrated around the optimum. This has an effect on the variance: for the custom sampling, the kernel length scales are adapted (section 4.6) to take into account the small distance between most of the sampled points. The LCB sampling places its points further apart and covers most of the low-drag region. The reason for this is that a low value β = 0.1 had to be chosen in order to even reach the optimum.

This test with LCB exhibits the same problems that were detected for Gramacy-Lee in section 5.8.2: the initial sampling plan contains only the four corners of the domain, which correspond to the maxima of the drag for small thicknesses. Therefore, the left side of the domain repels the sampler, which never finds the minimum between the two maxima if β is chosen too high (figure 7.5). The custom sampler does not have this problem, since the points on the left side have a high cross-validation error which compensates for their high drag. Thus, these regions get sampled and the minimum is found.

Another difficulty with the LCB sampling is that no clear guideline exists to choose β, since the infill criterion mixes the prediction and the variance, which have different units and may have different orders of magnitude. The value of β = 0.1 was found by trial and error. The custom sampler on the other hand is based on the cross-validation, which has the same unit and the same order of magnitude as the prediction (since the variance is non-dimensionalized, equation (5.7)). Thus, the value β = 1 can be chosen for any problem. The effect of reducing β is given in figure 7.6 for two optimization problems: varying C and T , and varying C and X. In both cases, choosing β = 0.5 leads to more exploration and a wider spread of the points. This is confirmed in figure 7.7 which shows the convergence of the optimizations in terms of the best C d found, and the distance of each sampled point to the optimum. For β = 1, this figure shows that the distances to the optimum reduce quickly, so the method switches to exploitation early and finds the optimum rapidly. However, the method gets back to exploration every now and then. For β = 0.5, there is more exploration and the convergence is slower; this is a safer option which makes it more likely that the method finds the global optimum if there are local minima. The last image gives the convergence for α = 0 and β = 0.1: this optimization concentrates too much on exploration and the true optimum is not found, for the reasons mentioned above.

Thus, even for optimization where both the LCB and custom sampler concentrate their points around the optimum, the custom sampler is the safest option: it is less likely to miss minima that are located close to maxima, and it is easier to set the sampler parameters.

Shape optimization of 2D profiles

While the previous section tested the shape optimization of a two-dimensional airfoil, the current section studies these optimum geometries. First, a full optimization of the three NACA four-digit parameters is presented. Then this optimization is repeated for different operating conditions, to see how these affect the optimum shape: the profile is placed below a free surface, the blunt trailing edge is replaced by a sharp one, and the velocity is changed. Finally, the results are used to analyse the robustness of the optima found.

Monofluid optimization

The first optimization concerns the profile in monofluid conditions, like in section 7.1.3. All three geometry parameters are optimized. The custom acquisition function is used with α = 1 and β = 1, the number of simulation points is increased to 80 in order to thoroughly search the larger parameter space.

The convergence of this optimization is given in figure 7.8. Compared with the twoparameter case (C and T , figure 7.7), the convergence is slower and there is more exploration, even at the end of the optimization. A possible explanation for this is that the drag is not very sensitive to the position of the maximum camber, so many camber positions give results close to the optimum and must be searched. The first row of table 7.2 provides the optimal shape parameters and the drag coefficient. The optimum profile is thin, its thickness corresponds to the minimum value allowed for T . In section 7.1.2, it was shown that at small values of T , the drag is sensitive to the camber C and that the optimum lies close to geometries with high drag. This is explained in figure 7.9, which compares the velocity field at the three-parameter optimum with the extreme-value cases for C and T at X = 0.4. The sharp noses of the thin profiles (C = 0.03) are sensitive to separation (section 2.4.2) at the leading edge, so the thin-flat profile has a large separation bubble on its upper side, while the thin-cambered profile has separated flow in the hollow below the leading edge. Both situations lead to a high drag. The thick profiles (C = 0.12) on the contrary are more tolerant. Although the flow is not perfect (for example, the thick cambered profile has its stagnation point on the top surface), the thick nose allows the flow to remain attached.

Then why is the optimum profile thin? For the optimum geometry, the camber line is perfectly aligned with the incoming flow. Thus, there is no need for the flow to move around the leading edge, so the suction peak is eliminated and no separation appears. In this situation, the thinnest profile has the smallest pressure drag so it is the most advantageous. Notably, the stagnation region is small compared to the other geometries. This is confirmed in figure 7.10 which shows the pressure coefficients. The thick profiles have the sharp suction peaks for which the NACA 4-digit series is infamous; on the cambered profile the suction peak is on the bottom side so the nose generates a negative lift. This is also the case for the thin cambered profile, but the suction peaks for the thin profiles are cut off due to separation. Finally, the optimum profile has no suction peak at all and the stagnation-point flow is limited to the leading edge itself, contrary to the other geometries.

Influence of the water surface

A lifting hydrofoil always operates in the proximity of the water surface (section 2.4.3). Nevertheless, geometrical optimizations are often performed without taking into account the free-surface effects, since this reduces the computational costs and allows the use of simpler simulation tools like the panel code Xfoil. Therefore, it is important to know 7.2 Shape optimization of 2D profiles whether the presence of a free surface has an influence on the shape of the optimum profile. Furthermore, the distance of a foil to the surface may vary for different operating points, which raises the question whether the optimum shape is dependent on the immersion. For these reasons, the optimization of section 7.2.1 is repeated with the profiles placed at a depth h below the undisturbed position of a free surface. Two cases are computed, h/c = 1 and h/c = 0.25. The results of the two optimizations, the geometry parameters and force coefficients, are provided in table 7.2. Figure 7.11 gives the optimized geometries at their operational angle of attack and the pressure coefficient on the profiles. Finally, the free-surface deformations and the hydrodynamic pressures (pressure minus the hydrostatic pressure) are provided in figure 7.12. Figure 7.12 shows that the free surface reduces the suction on the top surface of the hydrofoil. Since a part of the water above the foil is replaced by the much lighter air, the flow can move out of the way of the foil more easily, which explains this reduction (see also section 2.4.3). As a consequence, the hydrofoil needs to deflect the flow more than in the monofluid case in order to create the same lift coefficient C l = 0.6, so the optimal camber C increases when the water depth is reduced. Furthermore, the lift is created more by the high pressure below the profile than by the suction above it, since the free surface reduces mostly the effect of the suction side.

Case T C X C d mono C d h/c = 1.0 C d h/c = 0
For a moderate immersion of h/c = 1, the profile shape and pressure distribution are similar to the monofluid optimum (figure 7.11). However, taking into account the water surface for the optimization creates a performance gain: if the monofluid-optimized hydrofoil is placed at h/c = 1.0 and C l = 0.6, its drag is 1.2% higher than the hydrofoil optimized at h/c = 1.0 (table 7.2). The profile optimized for h/c = 0.25 has an unusual shape, which is probably optimal for creating compression below the hydrofoil, rather than suction above.

Finally, the cross-validation computations in table 7.2 show that each shape is truly optimal: in its optimized condition, it performs better than the other two profiles. This is an indication that the optimization procedure works as expected. Also, the table shows that C d increases significantly when h is reduced, due to wave drag. This effect is stronger than the difference between geometries, which means that a straightforward drag optimization for 2D hydrofoils is to place the foil as far below the surface as possible.

Influence of the trailing edge thickness

All realistic hydrofoil shapes have a blunt trailing edge, since a sharp edge is difficult to construct and would be too fragile for practical use. From a simulation point of view however, sharp trailing edges are preferred since blunt trailing edges require small cells to capture them which makes the simulations more expensive. For HEXPRESS, capturing a blunt trailing edge typically leads to a 30% increase in the total number of cells. Thus, it is tempting to simulate foils with a sharp trailing edge and then to replace these by blunt trailing edges for the actual foil construction.

An objection to this procedure is, that a blunt trailing edge adds considerably to the drag; discussions with teams during the last America's Cup indicated that a trailing edge of normal thickness is responsible for about 5% of the total drag for a hydrofoil. This was considered as a reason to prefer simulations with blunt trailing edges. However, for shape optimization the absolute drag is less important than the difference between geometries: if the best sharp-edged airfoil is still the best when the trailing edges are cut off, then the shape optimization can be performed with sharp trailing edges. Therefore, the optimization using sharp trailing edges is tested here. A three-parameter optimization is performed for the monofluid setting, but with the last coefficient in equation ( 6.1) changed to 0.1036 which gives a zero-thickness trailing edge. The optimal parameters and drag coefficients are compared in figure 7.13 and table 7.3. The table also contains a cross-validation exercise where both blunt and sharp profiles are created with each set of coefficients.

First of all, the table shows how accurate the optimization process is. The optimal coefficients are almost identical and it is hard to see the difference between the two airfoil shapes (figure 7.13). And yet, each set of parameters performs better than the other set in the case for which it was optimized. This shows that the optimization is sensitive enough to capture the small effect of the trailing edge shape. Although the blunt trailing edge thickness of the original profile is only 0.0006c, it increases the drag by about 1.5%. Therefore, it is indeed necessary to simulate a blunt trailing edge if the exact drag is needed. However, the difference in drag between the blunt profile optimized with sharp-edge computations (C d = 7.771) and the optimal blunt profile (C d = 7.764) is less than 0.1%. This means that for this particular family of airfoil shapes, it is justified to perform the optimization using sharp trailing edges, which saves computational time. It has to be tested if this conclusion still holds in more general cases and if it remains valid for thicker trailing edges.

Influence of the velocity

As shown above, the profiles optimized for C l = 0.6 are thin and the drag increases rapidly when the camber of the optimal profile is changed. This implies that an optimal thin profile may perform badly if the lift coefficient is changed. To test this, a final optimization is performed where the inflow velocity is increased from 10 to 12 m/s. To keep the same lift force, which is realistic for a lifting hydrofoil where the lift equals the weight of the boat, the target lift coefficient is decreased to C l = 0.417.

Table 7.4 and figure 7.14 show the result. As expected, for the lower lift coefficient the optimal camber is reduced, to keep the nose aligned with the incoming flow. The flatter profile reduces the pressure coefficient on the surface (although the lift remains the same since the velocity has increased). Finally, the cross-validation shows that the two profiles perform badly in each others operating conditions. This is in agreement with figure 7.2, which shows that at C l = 0.6 a profile with C = 0.0298 is near the limit of a leading-edge stall. Thus, for the optimal profiles obtained here, a 20% change in the velocity leads to serious performance loss.

Robustness of the optima

Especially in sailing, a hydrofoil is expected to operate in a range of different conditions (section 2.3.2), so a good geometry has a low drag which is insensitive to the operating conditions. The tests in this section served among others to assess this sensitivity.

Section 7.2.2 shows that the optimal shape is not very sensitive to the free surface proximity, unless the airfoil is very close to the surface: h/c < 1. Such conditions are relevant for surface-piercing foils, which could benefit from a varying profile section depending on the local distance to the surface. For fully immersed hydrofoils, it is preferable to place the airfoil deep below the surface to minimize wave drag. In this case, some performance can be gained by optimizing with simulations that take into account the free surface, but the 7.2 Shape optimization of 2D profiles optimum profile shapes are similar to the monofluid optimum. If desired, a way to optimize a profile for free-surface conditions without actually performing free-surface simulations would be to perform a monofluid optimization at a slightly higher lift coefficient than the target value.

At least for thin trailing edges, the use of a blunt trailing edge (section 7.2.3) imposes a drag penalty which is more or less independent of the profile shape. Therefore, it appears possible to perform shape optimizations with sharp trailing edges to reduce the mesh size, while constructing the actual foil with a blunt trailing edge. It has to be tested if this conclusion holds in more general cases.

The test of the dependence on the velocity and lift coefficient (section 7.2.4) shows that the very thin profiles obtained are optimized for one angle of attack. Small changes in incidence lead to separation at the leading edge and a significant increase in drag. Furthermore, the shapes obtained are impractical from a construction point of view. Thus, hydrofoil optimization for practical applications requires either geometrical constraints or multi-point optimization. Both could be performed with our technique. For example, a minimum thickness can be specified easily by limiting the range of the design parameters, while the objective function could be changed from the drag at one lift coefficient to a weighted average of the drag at two or more operating points. These points will be discussed in more detail in the conclusion, chapter 8.

3D hydrofoil optimization

This section shows the application of the adaptive metamodel-based optimization to three-dimensional hydrofoils created with the geometric modeler from chapter 6. After a description of the test case in section 7.3.1, the dynamic positioning approach used for the hydrofoil simulation is described in section 7.3.2. Section 7.3.3 presents a three-parameter hydrofoil optimization. The last part analyses this optimization: section 7.3.4 discusses the possible choices for the optimization parameters, while section 7.3.5 studies the robustness of the optimum.

Test case description

The aim of this test case is to optimize a lifting hydrofoil with respect to drag. The optimized parameters are the tip angle (TDA) and the twist of the tip, as well as the immersion depth. The tip length, the chord, and the shaft shape are kept fixed. With a chord of 35 cm and a span of 1.3 m, this foil is somewhat larger than a C-Class hydrofoil (section 6.10.1). The foil has a Tom Speer H105 profile like the Hydrocontest light boat, section 6.10.2.

The hydrofoil is optimized at a velocity of 10 m/s, for a constant vertical force of F z = 8000 N and lateral force F y = 6000 N . The density is ρ = 1026 kg/m 3 and the viscosity µ = 0.00122 kg/ms, for a chord-based Reynolds number Re = 2.94•10 5 . Wall-law boundary conditions are used. To match F y and F z , the foil's rake (pitch) and leeway (yaw, drift) angle are adjusted in FINE™/Marine using the procedure outlined below.

To perform this optimization test quickly, coarse meshes were used. These have about 2.5M cells with 200k faces on the hydrofoil and y + = 60. Earlier studies indicate that the numerical uncertainty in the drag for these meshes is of the order of 10%. Like for the 2D case, the free surface is captured with adaptive grid refinement. Each simulation takes about 8 hours on a 20-core Xeon workstation, which means that the total time for the 50-point optimization was around 2.5 weeks.

Dynamic positioning approach

For a foiling sailing yacht, the hydrofoil produces a vertical force, to lift the boat out of the water, as well as a lateral force to counteract the side force of the sails (section 2.3.2). For a given operating point, these forces are constant: the lift equals the weight of the boat, while the side force is given by the wind strength and the sail configuration. Thus, to compare different foil designs they should be simulated at the same vertical and lateral forces (rather than in the same attitude).

Like for the 2D case (section 7.1.1), this requires adjusting the orientation of the foil with respect to the flow. However, to control both the vertical and the lateral force (section 2.3.2), two angles will be adjusted: the rake and leeway angles. The dynamic positioning is performed by deforming the mesh (section 3.6.1) to change the rake, while the leeway is modified by rotating the entire mesh around the Z-axis. A difficulty for computing the angle corrections is that their effect on the forces is coupled: a change in the rake usually modifies the lateral force as well as the lift force, and vice versa. Therefore, the dependence of the forces on changes in the angles is computed as follows.

Let F y and F z be the lateral and vertical forces, φ y the rake, and φ z the leeway. Consider a 7.3 3D hydrofoil optimization small section of the hydrofoil surface dS with normal vector n (figure 7.15). Its surface parallel to the x-axis is:

d S = dS n 2 y + n 2 z n 2 x + n 2 y + n 2 z . (7.1)
Let θ be the angle of d S with respect to the horizontal plane: θ = arctan(n y /n z ). Then the change in angle of attack for d S due to changes in the main angles dφ y and dφ z is:

dα = dφ y cos θ + dφ z sin θ. (7.2)
Using the theoretical 2D lift slope dC l = 2π dα [START_REF] Anderson | Fundamentals of Aerodynamics[END_REF], the force change on the surface becomes: df = 2πq dα d S, (7.3) where the dynamic pressure is q = 1 2 ρV 2 . Finally, the force change is decomposed as df y = df sin θ and df z = df cos θ. Integrating these incremental forces over the entire wetted surface S leads to the total change in the forces:

dF y dF z = 2πq 1 2 S sin 2 θ d S 1 2 S sin θ cos θ d S 1 2 S sin θ cos θ d S 1 2 S cos 2 θ d S dφ z dφ y . (7.4)
The coefficients 1 2 before the integrals are added because the hydrofoil has an upper and a lower side, so the projected surface of the foil is integrated twice.

To simulate a hydrofoil for a given combination of side and lift forces, the computation is started with a fixed position. After a given number of time steps, the actual forces are evaluated and corrections for the angles are computed by solving 7.4, substituting for [dF y , dF z ] T the difference between the desided and the actual forces. The angle corrections are applied over a few time steps, then the forces are allowed to converge and a new correction is applied. Since hydrofoil forces converge rapidly, the corrections can be applied often: for the current simulations, the first correction comes after 30 time steps and each subsequent one after 20 time steps, with 10 time steps for the angle modification. No explicit underrelaxation of the computed angle corrections is applied: the 2D lift slope used in (7.2) is always an overestimation in 3D so the angle corrections are systematically too small, which stabilizes the process.

Depending on the geometry of the hydrofoil, it may not be possible to solve equation (7.4). For example, if the entire foil is vertical (i.e. a rudder), cos θ = 0 everywhere so dφ y cannot be computed. This is logical since a rudder cannot generate vertical forces. A fortiori, for any straight foil, cos θ and sin θ are constant and can be moved out of the integrals; the determinant of the resulting matrix is zero. Physically speaking, since the lift force on a straight foil is normal to its projected surface, such a foil cannot generate arbitrary forces in vertical and lateral directions independently.

The dynamic positioning procedure was started by MSc student Doriane Causeur in the context of the current project [START_REF] Causeur | Shape optimisation of hydrofoils using RANS computations[END_REF] and it was completed after her graduation.

Optimization of tip angle and twist

The optimization of the tip angle, the twist, and the immersion depth is performed with the custom sampler and parameters α = 1, β = 1. With the chosen sampling parameters, a safe and rapid convergence can be expected if the response surface does not have multiple optima. This is confirmed by the convergence graph 7.16 which shows an early transition to exploitation with a fast reduction of the minimum drag. The optimization then switches back to exploration, in order to confirm that the best point found is indeed optimal. Figure 7.17 shows the response surface and the sampled points. With respect to the two-dimensional airfoil, this surface is simpler for optimization: there is a clear optimum with a smooth variation around it and no large peaks in the response. This explains why 7.3 3D hydrofoil optimization the optimum is found quickly. The point cloud shows that the exploration is efficient, since the points are all situated in the region where the drag is lowest.

The geometries, the hydrodynamic pressure and the water surface are shown in figure 7.18 for the eight corners of the parameter space and in figure 7.19 for the optimum. Due to the way the simulation scripts work, the immersion depth D is actually an elevation, it is positive when the foil is lifted out of the water (like in section 6.10.1). Furthermore, the indicated drag coefficients are not fully non-dimensional, the drag has only been divided by q and not by the immersed surface. This allows to compare hydrofoils with different immersed surfaces in a fair way.

One of the eight extreme geometries caused a divergence of the dynamic positioning (right, second row in figure 7.18). This point is removed from the metamodel dataset with an exclusion zone as described in section 5.7.3.

A comparison with the eight extreme geometries makes it possible to understand why the optimal geometry has the lowest drag. The following factors play a role:

• The charge on the hydrofoil should be neither too small nor too large. This is illustrated by the first two images of figure 7.18. In the left one, the shaft area is large so to obtain the target F y , the pressure difference over the shaft is nearly zero. Thus, this surface creates friction drag for nothing. The right image shows the other extreme, where the same geometry is lifted mostly out of the water. To generate the same side force with the diminished shaft area, the leeway has increased and the pressure difference over the shaft is large. This contributes to the lift-induced drag; furthermore, such a difference in pressure close to the surface creates large waves which produce wave drag.

• No parts of the hydrofoil should work against each other. The clearest counter-example is the left figure on the second row of figure 7.18. Here, the tip is oriented such that it creates both a vertical force and a positive side force. If the shaft were to create a positive side force as well, the total side force would be too big. Therefore, the foil operates at a large rake to create the vertical force, while adopting a negative leeway so that the shaft creates a negative side force which compensates for the tip force. This is clearly not optimal, since these counteracting parts both produce induced drag. An implication of this principle is that any tip angle TDA < 90 o is inefficient, since this makes the tip produce a negative side force that the shaft has to compensate.

• The optimal shape has little twist. The geometric modeler uses an elliptic planform for the tip, which theoretically has the perfect lift distribution, without any twist. Thus, large twist angles only perturb the pressure distribution.

Figure 7.20 tries to clarify these aspects by showing isosurfaces of the second invariant of the velocity gradient tensor. The second invariant is an indicator for the presence of vortices. The figures show three aspects of the flow: the strength of the tip vortex, the presence or absence of a vortex at the radius between the tip and the shaft (which indicates an unequal loading on the tip and the shaft), and the amount of perturbation at the free surface. This last aspect is probably less important than it looks, since most perturbations occur in the spray region rather than below the surface. Still, the figure shows that the foils which create the least vorticity have the lowest drag. However, the figure is not precise since the coarse meshes used in the simulations were not made to capture the wake of the hydrofoil. Thus, the second invariant is not well represented in the wake.

In the response surface of figure 7.17, the lowest drag values lie in a disk-shaped region which is oriented diagonally with respect to the TDA and immersion axes. This disk represents the geometries where F z and F y are well balanced without parts of the foil working against each other. It is oriented diagonally because for lower tip angles, less F y is generated by the tip, so the shaft has to sink further in the water to compensate (i.e. the 'immersion' is reduced). Within the disk, the optimum is the point where the immersed surface is ideal, i.e. the pressure difference over the foil is neither too large nor too small.

Choosing the geometry parameters to be optimized

For hydrofoil optimizations, the choice of the design parameters to be optimized is not straightforward. On one hand, the optimization will be cheaper and faster if the number of parameters is reduced. On the other hand, if certain parameters are missing, then the design space may not contain the most efficient shapes. For the design exercise in section 7.3.3, the parameter combination was chosen more or less arbitrarily. Now, the results of the exercise can be used to determine other sensible combinations of parameters.

As seen above during the analysis of the optimized geometry, the drag of a hydrofoil depends on four items:

1. The immersed area, which determines the pressure difference required on the foil: too much area causes needless friction drag, while a small area increases the induced drag (which is proportional to the lift coefficient squared [START_REF] Anderson | Fundamentals of Aerodynamics[END_REF]).

2. The lift distribution, which should vary smoothly over the length of the foil. Any large variation and especially an area with negative lift increases the induced drag.

3. The aspect ratio; the higher the aspect ratio, the lower the induced drag.

4. The proximity to the surface, which increases the wave drag.

3D hydrofoil optimization

Of these items, the first two are the most important. The aspect ratio is determined by the span and the chord, its maximum is limited by the structure of the hydrofoil. This parameter may have an influence for foils operating near the water surface, where a larger aspect ratio could bring the tip closer to the surface and increase wave drag. But for immersed foils, since the drag decreases for larger aspect ratios, this parameter should be as big as possible. Thus, there is no reason to choose it for the optimization. Reducing the wave drag requires that the tip be placed far below the water surface. This may be in direct contradiction to the first two requirements and, in any case, the surface proximity is influenced by the same design parameters as the first two items.

The immersed area is determined by the size of the foil and by its immersion depth. For any optimization, at least one parameter which changes the immersed area is needed. If both the size and the immersion depth are fixed, suboptimal shapes are found. For example, consider the optimization of the TDA with the depth fixed. The two top left images in figure 7.18 show that if the foil is completely in the water, TDA 60 o gives less drag than TDA 120 o so the first shape would be optimal, even though the tip produces a negative side force.

In the same way, a good lift distribution starts with a correct orientation of the foil in the Y Z plane such that the right ratio of F y to F z is achieved when the foil is uniformly loaded.

While the immersion plays a role here, the main influences are the TDA and the cant angle of the entire foil. These parameters play more or less the same role and at least one must be present in any optimization.

The final optimum is then found by fine-tuning the lift distribution over the hydrofoil. Parameters which have an influence are for example the twist (as seen above), or the chord variation. Furthermore, the tip-shaft junction radius influences the low-pressure peak in the joint area (see figure 7.19) and it may be useful to vary the shaft angle independently of the tip, either by optimizing the cant and TDA together or by modifying the shaft radii (section 6.2).

Why are the immersion and the cant included as geometric optimization parameters, if they do not modify the geometry? It may seem more natural to optimize the immersion and the cant separately for each geometry, in an internal optimization loop during the computation. However, such an inner optimization requires several complete simulations, which have to be repeated for each geometry. If the effects of cant and immersion are comparable for similar geometries, then it may not be necessary to evaluate them completely for each geometry. So if the cant and immersion are included in the design parameters, an accurate response surface can be created by adding only a few points where the cant and immersion are varied. The adaptive sampler can figure this out automatically. Therefore, a cheaper optimization is obtained when the cant and immersion are optimized as design parameters.

7.3.5 Robustness, stability, and cavitation risks Single-point optimization as performed above may be too limited as an approach to designing good hydrofoils, for several reasons. First, as already noted for the two-dimensional case above, the speed of the boat will vary depending on the wind so the foil cannot be optimized for one speed only. Furthermore, a sailing yacht has usually at least two or three different operating points (for example, an America's Cup catamaran performs upwind, downwind, and reaching legs), for which the loads and speed ranges are different. Finally, a good hydrofoil has other properties besides low drag, such as good stability and limited risk of cavitation.

In robust optimization, the sensitivity to the velocity is usually taken into account by adding the velocity as a design dimension and then computing the sensitivity to the velocity for each geometry from the metamodel. This gives a second performance indicator, besides the drag. Thus, a Pareto front can be drawn based on the drag and its sensitivity to the velocity. However, in this case it may become impossible to perform adaptive sampling which targets the optimum. Instead, it is necessary to create a full response surface as a basis for the optimization.

A way to treat multiple operating points could be to make separate response surfaces for each point. Of course, the foil position for each operating point can be different so the merits of each geometry are determined by looking up the best position for that geometry in each response surface, and either combine the results in one objective function or construct a Pareto front. Once again, this requires full response-surface metamodels.

To assess the risk of cavitation on the hydrofoil, computations are normally performed with a cavitation model added to the Navier-Stokes equations (see for example [START_REF] Dauby | Simulations d'écoulements cavitants par résolution numérique des équations de Navier-Stokes en moyenne de Reynolds[END_REF] for a study of cavitation models in ISIS-CFD). However, a way to estimate the cavitation risk as a simple post-treatment of non-cavitating simulations is to plot isosurfaces of the pressure where it is equal to the water vapor pressure [START_REF] Bibliography Yvin | Tip vortex cavitation inception without a cavitation model[END_REF]. In the regions where these isosurfaces appear, there is a risk of cavitation. This exercise was performed for the foils of section 7.3.3 but the results are not presented here since they show nothing: the foils at 10 m/s move relatively slowly so the risk for cavitation is nil.

Finally, if no active control system is used, the ride height stability becomes important for the hydrofoil performance: when the foil goes up, the vertical force should reduce so the boat moves back down, attaining a stable flying height. The reason that many foiling sailing boats use upward-pointing tips (TDA < 90 o ) is that these provide excellent ride height stability.

The most accurate way to compute ride height stability is to perform multiple simulations for each point, with different immersion depths but the same rake and leeway (section 6.10.1). However, a series of computations like in section 7.3.3 can be used to assess the stability qualitatively, by plotting the rake obtained as a function of the immersion. If the rake for equilibrium increases when the foil moves out of the water, then the foil at a constant rake angle will lose lift when moving up. This means that it is stable.

Figure 7.21 shows this analysis for the current optimization. The figure shows that the stable rake around the optimum decreases slightly when the hydrofoil moves up (i.e. the 'immersion' increases), so the optimum foil is unstable. However, the trend is inversed and the foils become stable when TDA < 90 o . This could be a reason to prefer foils with a lower tip angle. On the other hand, it shows the interest of active ride height control which could allow to use foils with the lowest drag.

Conclusion

In this chapter, the adaptive sampling and optimization procedure has been applied to RANS simulations of hydrofoil geometries. The first series of tests on two-dimensional hydrofoils confirms what was found in chapter 5.8: the custom acquisition function is a safer choice than variance-based or LCB acquisition. For response surface creation, the custom sampler concentrates points around peaks in the response surface if this is needed, ensuring that the surface is well represented. For optimization, it is better than LCB for finding minima which are located close to maxima of the function. Finally, it is easy to configure the custom criterion since for optimization, the weight of the solution β can always be set between 0.5 and 1.0.

The optimization procedure is then used to study optimal shapes of two-dimensional hydrofoils. The optimum shape for minimum drag at a given lift coefficient is a thin profile whose camber line is aligned with the incoming flow. Thus, both separation and suction peaks on the nose are eliminated. Close to a free surface, a part of the low pressure above the foil is lost so the optimum profile becomes more cambered. Still, a profile optimized without a free surface performs well when the surface is not too close. The effect of including a blunt trailing edge has a non-negligible effect on the drag, however a profile optimized with a sharp trailing edge is nearly optimal when the trailing edge is cut off, so shape optimization can be performed with sharp trailing edges, saving grid cells. Finally, as expected, thin profiles are limited to operations at one speed, since their drag increases rapidly when the lift coefficient is changed.

To conclude the chapter, a three-dimensional hydrofoil is optimized. In order to simulate foils at a specified vertical and lateral force, a dynamic positioning procedure is used which adjusts the rake and leeway angles in order to obtain the correct forces. The optimal foil, obtained by varying the tip angle and twist as well as the immersion depth, has a pressure distribution which varies smoothly over the hydrofoil and which does not have parts of the foil pulling against each other. A general study of the design parameters for hydrofoil optimization shows that it should at least be possible to vary the immersed surface and to orient the lift vector laterally. The chapter ends with remarks about robust optimization, multipoint optimization for upwind and downwind sailing, and straightforward analysis of a foil's stability and cavitation risk. This analysis shows the potential conflict between the optimization for minimum drag and for stability.

8

Conclusion and future work

Conclusion

In this thesis a practical hydrodynamic optimization framework for hydrofoil shape design is developed. Foil hydrodynamics depend on the free-surface proximity, boundary layer effects including separation and stall which limit the maximum lift, and risks of ventilation or cavitation. Furthermore, hydrofoil ships in foilborne conditions must maintain a stable ride attitude and minimize motions, which creates a potential conflict with the design optimization for minimal drag. The use on sailing yachts requires hydrofoils which are efficient over a large range of operating conditions. And finally, there is a tension between the far-going optimization desired for high-performance foiling craft, and the limited budgets available for ship design. Thus, optimizing the shape of lifting hydrofoils is a complex process in which many factors play a role.

In simulation-based optimization, the hydrofoil designs are evaluated through numerical simulation of the flow. The geometries variations are controlled by parameters which together form the design space of all available shapes. Simulation-based optimization requires three separate techniques:

• The actual optimization, that involves the searching of the design space and the choice of the candidate geometries to be tested.

• The geometric modeler which transforms a set of design parameters into a threedimensional geometry.

• The fluid solver, used to assess the performance of each candidate geometry.

For an optimization process, performing each of these steps automatically is necessary. Furthermore, each operation must be suitable for the specific needs of the hydrofoil.

Optimization The main difficulty for the optimization is that simulations of hydrofoil flow are expensive, so the total number of designs that can be simulated have to be limited. Surrogate-based optimization provides an elegant solution to address this problem: a small number of data points is used to construct a metamodel of the simulation response and the optimization is performed over this metamodel, which is cheap to evaluate.

Gaussian Process regression is shown to work well for constructing the metamodel, since it is general and can handle different types of reponses. Through hyperparameter optimization, the interpolation kernels are adapted to the data, which makes the method flexible. And finally, the stochastic nature of the GP provides a basis for estimating the local uncertainty of the metamodel.

Metamodels based on numerical simulation have the advantage that the datapoints can be chosen freely. In sequential sampling, these points are selected adaptively: in each step, the point in the design space is sought where new data would be most useful for the metamodel. Then a simulation is performed in this point and added to the metamodel.

The best new point is defined as the maximum of an acquisition function, which is a weighted sum of the function to be optimized and the estimated uncertainty in the metamodel. Depending on the weights, the resulting metamodel is different. If the weight on the function is zero, only points with high uncertainty will be sampled, which gives a response surface that is reliable everywhere. With a high weight on the function, points are concentrated near its optimum. This is efficient for determining the optimal point, but dangerous since the search may get stuck in a local optimum. For reliable optimization, a compromise must be chosen.

To estimate the uncertainty of the metamodel, the classical approach is to use the variance of the GP regression. However, this parameter does not account for the local behavior of the modeled function and cannot distinguish regions where the function varies rapidly. Therefore, an uncertainty estimate is proposed which multiplies the variance with a crossvalidation error, computed by comparing the true simulation result in each datapoint with a GP prediction based on all the datapoints except the one being tested. This cross-validation detects the zones where the metamodel is sensitive to new points and therefore, indicates where the uncertainty is high.

Tests on analytical functions and on hydrofoil optimizations show that this custom acquisition function is a safer choice than classical acquisition functions. For smoothly varying functions it performs equally well, while it continues to perform for difficult functions with strong peaks, where the classical approach fails to approach the optimum in some cases. Furthermore, the cross-validation error has the same unit and the same order of magnitude as the function to be optimized itself, which makes it easier to choose the weights in the acquisition function. Thus, the custom function ensures a reliable optimization.

Geometric modeler To create hydrofoil geometries, two approaches are deformationbased modeling where a base geometry is subjected to parametrized deformation, and constructive modeling where the geometry is formed by parametrized basic elements. The choice between these two approaches was based on the user: design parameters must have meaning for naval architects, since practical experience is crucial for the design optimization. Experience helps designers choose the most relevant parameters to be optimized and provides a background for analysing the results. Studying the metamodel and its dependence on the design parameters may even help naval architects learn more about hydrofoil behavior, provided that the design parameters have a physical meaning.

124

Chapter 8 Conclusion and future work

These arguments resulted in the choice to use primarily constructive modeling based on parameters which are natural to a designer.

The hydrofoil modeler presented here is based on separate descriptions of the foil 'spine' and the airfoil sections placed along this spine. The general outline of the hydrofoil is given by the spine while the set of sections generates the hydrofoil outer surface. This principle assures the flexibility of the modeler, since the spine and the sections can be generated with different methods and exchanged independently, if needed. The main spine parametrization is based on tangent circle arcs. This parametrization proved to be adapted to the requirements for industrial applications and since the circle parameters have intuitive meanings, it produces results which are well-suited for analysis.

Since most airfoil shapes are defined as point clouds rather than parametric shapes, the profile sections are constructed with deformation techniques. However, the deformation is based on natural design parameters: the camber, chord and twist. The transformation of the sections and the blending of different profile shapes work well but do not give precise control over the camber. Therefore, it would be useful to have a more general parametric 2D profile modeler, to complement the deformation methods available now. The hydrofoil modeler has proven to be robust and to generate geometries of high quality. It was tested extensively in industrial projects and has shown to be reliable in such applications.

Simulation The flow simulations are performed with the ISIS-CFD unstructured finitevolume Navier-Stokes solver. Although these simulations are costly compared with more approximate methods like panel codes, they provide high-fidelity modeling of the flow physics that influence the optimization. For example, tests showed that free-surface effects have a non-negligible influence on the optimum hydrofoil shape, so these must be taken into account. Also, modeling flow separation and ventilation requires the accurate representation of viscous effects. Thus, for reliable optimization it is important to simulate these effects correctly. And efficient surrogate modeling can reduce the required number of simulations and compensate for the costs.

The ability to rotate the hydrofoils to achieve a desired lift is important for the efficiency of the optimization. The adaptive sampler could also use fixed-position simulations to obtain an optimum at a given lift, by using constrained optimization. In this case, design points with the wrong lift are accepted but gradually excluded from the sampling zone. The dynamic rotation however ensures that every simulation has the correct lift and removes the need to add the angles of incidence (rake, leeway) as design parameters.

Finally, in an automatization loop it is essential to perform the numerical simulations automatically. Several factors contribute to this capacity. First, the Python interface provided for FINE™/Marine facilitates the scripting of the computations. Furthermore, meshing the free surface with adaptive grid refinement ensures that the water surface is well captured in all cases. And last but not least, the geometrical modeler produces the

Conclusion

125

geometry definition for the computations directly, which ensures that each geometry is flawless. This is a major help in ensuring the reliability of the computations.

Optimization results

The procedure developed here has been used successfully to optimize two-and three-dimensional hydrofoil shapes. This shows first of all that the optimizer, geometric modeler, and flow simulation function together in a fully automatic way. The computation time required is significant but not prohibitive: on a large workstation, an optimization in 2D takes a day and in 3D it requires about two weeks. If needed, these computation times can be reduced by running the simulations on clusters with more processors.

The accuracy of the optimization procedure is demonstrated by the cross-validation of profile shapes optimized for different conditions. For all the parameter variations studied, i.e. the water surface proximity, the trailing edge thickness and the velocity, the shape which is optimized for a given condition performs better in that condition than the other shapes. In the trailing edge study, where the differences in drag between the shapes are around 0.1%, this is a strong indication of the optimizer performance.

The results of the optimization show the design parameters which are important for lifting hydrofoils. If a 3D foil must generate a side force as well as a lift, then a parameter is needed which sets the balance between these two forces (such as the tip angle) in order to obtain an even load on the foil. Furthermore, the immersed area must be adjustable so the foil is not too big (which leads to needless friction drag) nor too small (which increases the induced drag).

In 2D, the camber is an important parameter, since the drag is lowest when the camber line is aligned with the incoming flow. This could be a reason to consider the camber distribution as an optimization parameter also in 3D. Finally, the optimal profile shape depends on the distance to the free surface, especially close to the surface. Thus, for surface-piercing foils it is important to allow the profile to vary along the spine.

All the geometries were optimized for a single point of operation. The results show that this has its limitations: the thin 2D profiles obtained are ill suited for velocities far from the optimal one, while the optimal tip angle in 3D depends on the ratio between the required vertical and horizontal forces that was imposed. Thus, more robust optimization would be of interest, as discussed below.

Hyrofoils cannot be optimized for low drag only. Compromises are required to account for flight stability, cavitation and ventilation risks, etc. The current procedure can assess these risks in a quantitative way: cavitation is studied with isosurface plots of the vapor pressure and the stability is inferred from the stable incidence at different immersions. Including these data in the optimization is a subject for further study.

126

Chapter 8 Conclusion and future work

Perspectives

The optimization procedure presented here has proven its worth in tests, but it is still in the initial phases of development. The framework that was designed is flexible and offers rich possibilities for further research. Regarding the use in an industrial context, two main areas of development merit attention:

• Decreasing the cost of the optimization even further,

• Making the optimization results more robust and the procedure more versatile.

These two items are discussed in this final section.

Efficiency of the optimization

Even with the help of surrogate models, simulation-based optimization remains costly since all high-fidelity computations are time-consuming. To further accelerate the optimization, the surrogate creation can be enhanced to perform optimizations which either require less simulations, or less costly ones. A number of such approches is suggested here.

Feature selection

The tests in section 7.1 show that some design parameters have a much greater influence on the design objective than others. It is wasteful to use many sample points in order to determine the influence of such nearly-irrelevant parameters.

Anisotropic kernels for the Gaussian Process regression can detect less relevant directions, provided that the length scale hyperparameters in each direction are optimized separately (section 4.6). Such Automatic Relevance Determination allows to capture the influence of the less-relevant directions accurately with a few sample points. Once the sequential sampling detects that these directions are captured accurately, it will not oversample them.

Furthermore, if multiple optimizations of similar hydrofoils are planned, then the result of the ARD for the first optimization can be used to select the most important design parameters. Subsequent optimizations could then be limited to these parameters only.

Multifidelity To reduce the average cost of the simulations, it would be interesting to combine computations with different accuracies. The idea of multifidelity optimization is to capture most of the metamodel behavior with cheap, low-fidelity simulations and to correct these results using a few high-fidelity simulations dispersed over the design space. [START_REF] Pellegrini | Resistance and payload optimization of a sea vehicle by adaptive multi-fidelity metamodeling[END_REF] introduce two separate metamodels, one for the low-fidelity response and a separate model with fewer points for the difference between the high-and low-fidelity results. Thus, a combined metamodel is obtained which can be as accurate as a high-fidelity model, using mostly low-fidelity simulation.

Perspectives

Multifidelity results are mostly obtained by using two separate simulation tools, which increases the complexity of the procedure. An alternative for Navier-Stokes simulation is to perform the low-fidelity simulations using coarse meshes. Thus, all the flow physics captured by the high-fidelity model are present in the low-fidelity one, while the simulation costs are still reduced significantly. Creating the meshes with adaptive grid refinement (section 3.6.2) allows to adjust the grid density as desired, by varying a single global threshold parameter. This approach is under development [START_REF] Wackers | Cfd-based shape optimization under limited computational resources -a study on adaptive multi-fidelity metamodeling[END_REF] and it would be interesting to combine it with the current optimization method.

Noise filtering for numerical errors In the framework of this thesis, the GP regressions for the metamodels currently do not use noise models, since the simulation results are considered exact. However, as mentioned briefly in section 4.4, the effect of numerical errors could be interpreted as noise, provided that the errors in different design points are uncorrelated.

In simulation-based optimization, this is an unorthodox requirement. Generally, efforts are made to have the same numerical errors in design points which are close, such that these errors cancel when two nearby results are subtracted. This ensures that the gradients of the response are correct; it is achieved by performing computations as much as possible on the same mesh, deformed to suit each geometry. However, if a different unstructured mesh is generated with HEXPRESS for each geometry, the numerical errors in the design points will be more or less uncorrelated. Thus, it may be possible to filter them out as noise and to end up with a metamodel that is more accurate than the simulations on which it is based.

Tests have to show if this approach is actually possible.

Efficient response surface creation based on VPP information Apart from automatic geometrical optimization, metamodels also play a role in the performance evaluation of the entire foiling ship through a Velocity Prediction Program (chapter 1), which requires that the behavior of the hydrofoil in different operating conditions is captured as a response surface. Chapter 5 shows that the sequential sampling framework is efficient for creating response surfaces by simply switching the acquisition functions. These response surfaces are built based on assumed operating ranges of the considered hydrofoil. However, the precise operating conditions of the boat result from an equilibrium and are not known beforehand. Therefore, when building the response surface, the design space explored is chosen large, which means that zones are sampled where the foiling vessel will never operate.

Inserting a VPP in the optimization loop could provide an efficient way of reducing the size of the design space. After performing an initial sampling on the entire space, the VPP would use the resulting metamodel and compute the operating conditions of the hydrofoil. Translated into a bound constraint, this information can serve to limit the space that is sampled, increasing the efficiency of the adaptive sampling. Obviously only the dimensions related to the operating conditions are involved. Before sampling each new point, the VPP is run to refine the predicted operating space. Thus, by creating a direct interaction with the VPP, the response surface is optimized for the actual operating conditions.

Robust and versatile optimization

The optimization framework presented here is powerful and flexible, so it can be adapted to other types of optimizations than the ones presented in this thesis. For example:

Robust optimization The main limitation of the results found in chapter 7 is that every geometry is optimized for a single operating point. For instance, the shape of the optimal 2D profiles depends strongly on the operating speed and the optimal profiles perform badly when the speed is changed. For practical applications, geometries are desired whose performance is more robust, i.e. geometries which work well over a range of conditions.

This requires multipoint optimization, with an objective function that is a summation or an integration over the operational conditions. Let f be a function of x = x d ∪ x v , where x d are the design parameters that are to be optimized, while x v contains the operating conditions which are variable. The optimization problem (5.1) then becomes:

x opt d = arg min x d ∈X d Xv w(x v )f (x d , x v ) dx v , (8.1)
where w(x v ) is a weighting of the operational conditions that determines the relative importance of each condition (for example, it could indicate the effective speed range).

Such an optimization can be performed with the current procedure. Since the dependence of f on both x d and x v must be determined, a metamodel is constructed that has x as input parameters. Then f in the acquisition function (5.4) is replaced by the integral from equation (8.1):

a(x; D n ) = U (x; D n ) -β Xv w(ξ ξ ξ v ) f (x d , ξ ξ ξ v ; D n ) dξ ξ ξ v . (8.2)
Thus, the contribution from f depends only on x d , it is independent of x v . This means that sample points will be concentrated around the optimal x d , where the (unmodified) uncertainty U (x; D n ) ensures that f is well represented for all values of x v .

The price to be paid for the robust optimization is, that the dimension of the metamodel increases through the addition of the degrees of freedom in x v . Furthermore, the metamodel evaluation for the acquisition function now requires a (numerical) integration instead of an evaluation in a single point. However, the capacity to create designs which are guaranteed to function over a range of velocities, immersion depths, etc. offers such an important practical advantage that this extra cost is acceptable. The robust optimization will be implemented as soon as possible.

Perspectives

Multiobjective optimization A hydrofoil cannot be optimized for minimum drag only; the optimization has to be placed in a wider context where stability, cavitation and ventilation risks are also taken into account. This creates conflicts for the optimization, since for example the profile with the lowest drag is not necessarily stable (section 7.3.5).

Assessing the drag and the stability requires a multiobjective optimization, where a Pareto front is constructed of geometries whose drag cannot be improved without reducing the stability, and vice versa. The Pareto front can be extracted from the metamodel. However, this requires a metamodel which is reliable throughout the domain, sampled with an acquisition function based on the uncertainty only. In principle it would be possible to construct the Pareto front during each sampling and to place the points close to the front, but since the stability is evaluated from a derivative of the metamodel, it is difficult to compute it accurately without a generally reliable metamodel.

The cavitation risk is a constraint, rather than an optimization parameter. The risk can be evaluated with an automatic computation of the volume where the pressure is below the vapor pressure (through the summation of cell volumes), to be included in ISIS-CFD. Then, this parameter is included in the optimization to define an exclusion zone (section 5.7).

For sailing foilers, the operation in different regimes (upwind or downwind sailing) can be treated through multipoint optimization as described above, by combining the performances at different speeds in one objective function. An alternative is to treat the performance in each operating points as a separate objective. The Pareto front provides the design teams with options for the foil geometry, having different performances in each regime.

French extended summary Cependant, ce sont les essais en bassin de carène qui ont permis de comparer de manière formelle les performances des différents candidats en vue de sélectionner le meilleur avant sa construction. Aujourd'hui, la sélection de candidats à la conception est de plus en plus souvent réalisée par simulation numérique.

A mesure que la puissance des ordinateurs augmente, les architectes navals commencent à se tourner vers l'optimisation basée sur la simulation, où non seulement les simulations, mais aussi le choix des candidats à la conception et la génération de la géométrie sont traités de manière connectée et entièrement automatisée. Cette thèse présente une procédure d'optimisation basée sur la simulation, pour les navires équipés de surfaces portantes ou hydrofoils. Elle décrit comment les éléments de la procédure ont été développés pour répondre aux exigences d'efficacité et de flexibilité que cette application requiert.

séries de deux ou trois unités. Cela impose de limiter les coûts liés à la conception et l'optimisation, ces couts ne pouvant être répartis entre les différentes unités.

Certains des navires les plus rapides au monde sont équipés d'hydrofoils. Il s'agit de surfaces portantes similaires à des ailes d'avion, situées sous la surface de l'eau, qui soulèvent de manière partielle ou totale la coque du navire hors de l'eau afin de réduire sa traînée. Des foilers à moteur ont été mis au point au cours de la seconde moitié du siècle dernier pour assurer des missions de patrouille maritime rapide ou de transport de passagers.

Plus récemment, l'intérêt pour les hydroptères s'est ravivé avec l'avancée des voiliers à foils. La 34 e édition de l'America's Cup, qui s'est tenue à San Francisco en septembre 2013, a vu des catamarans de course de 72 pieds "voler" sur foils et atteindre des vitesses bien supérieures à celle des voiliers de course traditionnels. Cette édition de l'America's Cup a particulièrement marqué les esprits et a convaincu le grand public mais également beaucoup de professionnels de l'intérêt des hydrofoils appliqués aux voiliers et a conduit par la suite à équiper d'hydrofoils divers yachts de plus petite taille.

Par beaucoup d'aspects et en particulier pour ce qui concerne leur optimisation, les yachts de course sont similaires aux navires de marine marchande. D'une part, une optimisation poussée est nécessaire, car de minuscules différences de performances peuvent faire la différence entre gagner et perdre une course. D'autre part, les budgets disponibles pour la conception d'un voilier sont faibles. En comparaison, alors que les budgets de l'America's Cup sont au moins dix fois plus élevés que ceux disponibles dans toute autre compétition de voile, le coût total d'une campagne de l'America's Cup est inférieur au prix d'un grand avion de ligne (dont des centaines sont construits). Ainsi, l'optimisation de la forme des yachts doit être à la fois efficace et rentable.

Optimisation par modèle de substitution L'optimisation de formes nécessite de définir celles-ci à l'aide de grandeurs géométriques qui constituent les paramètres de l'optimisation. La combinaison de ces paramètres constitue l'espace d'optimisation. Il convient pour réaliser une optimisation de chercher dans cet espace la combinaison de paramètres permettant les meilleures performances. Lorsque l'espace d'optimisation est grand ou lorsque le comportement du navire est complexe, de nombreuses évaluations de combinaisons de paramètres sont alors nécessaires. Cependant, chaque évaluation nécessite une simulation numérique complète et le nombre total de simulations qui peuvent être exécutées est généralement limité par le temps ou les coûts financiers. Une manière de palier à cette limitation consiste à construire un modèle simplifié du comportement du navire à partir d'un nombre limité de simulations correspondant à certaines combinaisons choisies de l'espace d'optimisation. La méta-modélisation, aussi appelée modélisation de substitution, est le processus de génération de tels "modèles de modèles" ou "meta-modèles". L'optimisation proprement dite est ensuite effectuée en interrogeant le méta-modèle, plutôt que d'exécuter

138

Extended summary des simulations réelles. Puisque les évaluations d'un modèle de substitution sont rapides et peu couteuses, le coût de calcul de l'optimisation s'en trouve réduit.

L'échantillonnage est le processus consistant à sélectionner des points de l'espace d'optimisation pour ensuite effectuer les simulations correspondantes et recueillir les réponses qui serviront à construire le méta-modèle. Le nombre de points d'échantillonnage nécessaires à l'approximation du comportement d'un modèle numérique dépend de la complexité de la réponse, qui n'est pas nécessairement connue à l'avance. La stratégie d'échantillonnage choisie joue un rôle important dans la performance d'une optimisation par méta-modèle: le sous-échantillonnage peut ne pas permettre de saisir la complexité des phénomènes et le sur-échantillonnage entraîne des temps de calcul prohibitifs.

L'échantillonnage peut être fixe, c'est-à-dire que points de simulations sont définis lors d'une étape initiale unique, ou adaptatif. Lorsque l'échantillonnage est adaptatif un nombre limité de points est initialement sélectionné puis simulé et les réponses des simulations sont utilisées afin de définir le ou les points à échantillonner ultérieurement. La construction du méta-modèle se déroule donc de manière séquentielle. Avant chaque nouvelle étape de simulation, l'espace d'optimisation est évalué et un point est sélectionné à l'emplacement où une nouvelle simulation serait la plus bénéfique pour le méta-modèle ; ce point est simulé puis ajouté au méta-modèle. Ainsi, le nombre requis de simulations est minimisé pour la construction efficace de méta-modèle.

Optimisation et surfaces de réponse L'optimisation d'un hydrofoil peut être réalisée de différentes manières. Les hydrofoils peuvent être optimisés indépendamment du reste du navire sur lequel ils seront installés ou, de manière plus robuste, ils peuvent être optimisés en prenant en compte l'ensemble du navire et en simulant le comportement de ce navire dans un code de simulation tiers tel qu'un programme de prédiction de vitesse (VPP). Les concepteurs passent généralement d'une approche à une autre en fonction de l'avancée du projet de conception.

Les méta-modèles peuvent être utilisés pour les deux approches mais, s'ils sont construits de manière adaptative, ils seront échantillonnés de manière tout à fait différente. Les VPP sont basés sur des méta-modèles séparés pour chaque composant d'un navire, qui sont combinés par l'intermédiaire du VPP. Il faut donc des méta-modèles précis dans tout l'espace de paramètres. Pour créer de telles surfaces de réponse, on privilégiera l'échantillonnage dans les zones à forte incertitude, afin d'assurer leur fiabilité partout.

D'autre part, l'optimisation automatique de la géométrie exige l'échantillonnage des zones susceptibles d'offrir une amélioration par rapport à la meilleure conception actuelle. Cela signifie qu'une fois que le comportement global du méta-modèle est connu, il peut être bénéfique de concentrer la zone d'échantillonnage près de l'optimum.

La chaîne d'optimisation L'optimisation basée sur les méta-modèles nécessite l'automatisation de trois processus :

• L'optimisation sur l'espace de conception, qui comprend la construction du métamodèle et l'échantillonneur adaptatif qui choisit les points qui seront simulés.

• Le modeleur géométrique qui traduit un ensemble de paramètres de conception en une géométrie tridimensionnelle, nécessaire à la simulation ou à la fabrication.

• La simulation numérique des fluides qui évalue le flux autour des géométries candidate et fournit les réponses sur lesquelles le méta-modèle est basé.

Ces trois éléments sont appelés dans une boucle qui tourne jusqu'à ce qu'un critère de convergence soit atteint (par exemple lorsque la qualité du méta-modèle est jugée suffisante).

Leur couplage est effectué par l'intermédiaire d'un script qui organise la boucle et stocke les résultats à chaque étape et fournit des données nécessaire pour poursuivre le processus.

Hydrofoils

Le chapitre 2 de cette thèse présente les hydrofoils et décrit les différentes configurations qu'ils peuvent adopter. Un bref aperçu historique retrace les origines des formes d'hydroptères utilisées actuellement. Le chapitre se conclut par une présentation de la physique des écoulements auxquels sont soumis les hydrofoils dont la nature impacte le comportement global. Ces éléments sont à prendre en compte pour les simulations fluides numériques et ainsi que l'optimisation.

Simulation numérique des fluides

Le chapitre 3 traite ensuite de la simulation numérique de l'écoulement (CFD) autour des hydrofoils. Après une introduction à la simulation des fluides dans les applications marines, le chapitre présente le solveur fluide ISIS-CFD, utilisé dans cette thèse. Les équations de la dynamique des fluides sont introduites, ainsi que les aspects spécifiques de la discrétisation qui concernent les écoulements avec surface libre. La dernière partie du chapitre décrit les techniques de modification du maillage telles que la déformation de maillage et le raffinement adaptatif.

Régression par processus gaussien

Les simulations fluides numériques sont souvent trop longues a réaliser pour être utilisés directement dans la résolution des problèmes d'optimisation. La solution retenue pour contourner cet obstacle consiste à remplacer le code numérique complexe par un modèle simplifié ou méta-modèle, reproduisant les réponses du code mais dont le temps de calcul 140 Extended summary est très inférieur. La methode des meta-modeles nécessite de pouvoir interpoler des points échantillonnés sur l'espace d'optimisation. Le chapitre 4 présente la régression par processus gaussien (RPG) qui est la méthode utilisée dans cette thèse pour construire les méta-modèles.

La RPG traite la fonction à reconstruire comme s'il s'agissait de la réalisation d'un processus stochastique et fonde son estimation sur les valeurs les plus probables de ce processus, étant donné que les valeurs de la fonction dans les points de données sont connues. Cette approche peut sembler contre-intuitive si les données réelles proviennent d'une source déterministe comme la simulation numérique, mais elle offre l'avantage de fournir une estimation de l'incertitude de la prédiction basée sur la variance du processus statistique. En effet cette méthode tient compte de la structure statistique spatiale de la variable estimée en considérant la corrélation entre les réponses du code, en fonction de la distance entre les variables d'entrée.

Le modèle de processus gaussien est donc caractérisé par sa fonction de moyenne mais surtout de covariance. L'utilisation du processus gaussien comme métamodèle implique un choix a priori des fonctions de régression et de covariance dont les parametres seront optimisés pour reproduire le plus fidèlement possible la réponse du code.

La régression par processus gaussien offre de nombreux avantages pour approximer un code de calcul complexe. C'est un interpolateur exact tout en étant un modèle souple, il permet une formulation analytique du prédicteur et de la variance. Le paradigme GP offre également un moyen simple de traiter le bruit dans les données d'entrée, en ajoutant une composante diagonale à la matrice de corrélation pour modéliser du bruit qui n'est pas corrélé entre les points de données. Bien que la modélisation du bruit ne soit pas utilisée dans cette thèse, elle est potentiellement intéressante comme moyen de réduire les perturbations dues à des erreurs numériques dans les simulations.

Dans ce chapitre, différents fonctions de covariance sont présentées et testées dans le but de déterminer la plus adaptée aux applications considérées ici. Les tests montrent que la mise en oeuvre de la régression GP adoptée ici réussit à reconstruire des fonctions avec des variations lisses ou des pics locaux. Cependant, les résultats dépendent du nombre de points de données et des grilles fines sont nécessaires pour capturer les pics.

Puisque les régressions du processus de Gauss n'exigent pas une distribution régulière des points, il serait possible de concentrer les données près des caractéristiques les plus importantes de la fonction à interpoler. Cela soulève la question du placement optimal des points pour une fonction donnée. Cette question est traité dans le chapitre suivant.

Optimisation par méta-modèles

Le chapitre 5 est consacré à l'optimisation basée sur les modèles de substitution qui consiste en une combinaison de création de méta-modèle et d'optimisation. Un méta-modèle (chapitre 4) est construit sur la base d'un nombre limité de simulations CFD et l'optimum est recherché en explorant ce modèle de substitution plutôt qu'en réalisant des calculs CFD coûteux.

Après une introduction à l'optimisation basée sur la simulation et les méta-modèles, le chapitre présente une technique adaptative destinée à construire efficacement et de manière fiable les méta-modèles à différentes fins d'optimisation. L'échantillonnage adaptatif y est abordé comme un élément essentiel de la méthode. La sélection d'un point à échantillonner est formulée comme un problème de maximisation d'une fonction d'acquisition qui combine les valeurs du méta-modèle et une estimation de l'incertitude du méta-modèle. Une fonction d'acquisition custom est présentée avec une estimation de l'incertitude basée sur la variance du processus gaussien, pondérée par une estimation de l'erreur de validation croisée dans les points échantillonnés.

Le chapitre se termine par des tests de l'algorithme adaptatif sur les fonctions analytiques.

Les tests montrent que cette nouvelle fonction est plus fiable que les fonctions d'acquisition existantes. Elle donne les mêmes performances pour des cas test simples, tout en offrant de bonnes performances dans des situations plus difficiles où les approches existantes échouent.

Modeleur géométrique

Le chapitre 6 présente le développement d'un modeleur géométrique dont le rôle est de générer des géométries d'hydrofoils destinés à équiper des voiliers. Le principe du modeleur consiste à définir une courbe génératrice qui va donner la forme générale de l'hydrofoil puis à venir disposer le long de cette génératrice des sections permettant de définir l'enveloppe extérieur du foil. Ces sections sont paramétrées et leur corde, twist, . . . peuvent être variés en fonction de leur position le long de la génératrice. Au cours de cette thèse, le modeleur géométrique a été utilisé dans le cadre de projets industriels et pédagogiques. Pour la création de la surface de réponse, l'échantillonnage basée sur la fonction custom concentre les points autour des pics de la surface si nécessaire, en s'assurant que la surface est bien représentée. Pour l'optimisation, la fonction custom est préférable pour trouver des minima qui sont situés à proximité des maxima de la fonction. Enfin, il est aisé de configurer le critère custom puisque pour l'optimisation, le poids de la solution β dans la fonction d'acquisition peut toujours être compris entre 0,5 et 1,0.

La procédure d'optimisation est ensuite utilisée pour étudier les formes optimales des hydrofoils 2D. La forme optimale obtenue pour une traînée minimale à coefficient de portance imposé, est un profil mince dont la ligne de cambrure est alignée avec le flux entrant. De cette manière, les pics de séparation et de succion sur le nez sont éliminés. Quand un hydrofoil est placé à proximité de la surface libre, une partie de la dépression générée sur l'extrados est perdue de sorte que le profil optimal devient plus bombé. Cependant, un profile optimisé sans surface libre donne de bons résultats lorsque la surface n'est pas trop proche. L'effet de l'inclusion d'un bord de fuite épais a un effet non négligeable sur la traînée, mais un profil optimisé avec un bord de fuite fin est presque optimal lorsque le bord de fuite est coupé, de sorte que l'optimisation de la forme peut être effectuée avec des bords de fuite fin, ce qui permet d'économiser des cellules de maillage. Enfin, comme attendu, l'utilisation des profils minces est limité a des plages de vitesse très réduites, en raison de l'augmentation rapide de leur traînée lorsque le coefficient de portance est modifié. Abstract : This thesis presents a practical hydrodynamic optimization framework for hydrofoil shape design. Automated simulation based optimization of hydrofoil is a challenging process. It may involve conflicting optimization objectives, but also impose a trade-off between the cost of numerical simulations and the limited budgets available for ship design. The optimization framework is based on sequential sampling and surrogate modeling. Gaussian Process Regression (GPR) is used to build a predictive model based on data issued from fluid simulations of selected hydrofoil geometries. The GPR model is then combined with other criteria into an acquisition function that is evaluated over the design space, to define new query points that are added to the data set in order to improve the model. A custom acquisition function is developed, based on GPR variance and cross validation of the data.

A hydrofoil geometric modeler is also developed to automatically create the hydrofoil shapes based on the parameters determined by the optimizer. To complete the optimization loop, FINE/Marine, a RANS flow solver, is embedded into the framework to perform the fluid simulations. Optimization capabilities are tested on analytical test cases. The results show that the custom function is more robust than other existing acquisition functions when tested on difficult functions. The entire optimization framework is then tested on 2D hydrofoil sections and 3D hydrofoil optimization cases with free surface. In both cases, the optimization process performs well, resulting in optimized hydrofoil shapes and confirming the results obtained from the analytical test cases. However, the optimum is shown to be sensitive to operating conditions.

DragαFigure 2 . 1 :

 21 Figure 2.1: Lift and Drag forces generated by a foil, with an angle of attack α.

Figure 2 . 2 :

 22 Figure 2.2: Sustention triangle.

Figure 2 . 3 :

 23 Figure 2.3: Longitudinal distribution of hydrofoils (bottom view; the grey rectangles represent the lifting surfaces).

Figure 2 . 4 :

 24 Figure 2.4: Hydrofoil configurations.

Figure 2 . 5 :

 25 Figure 2.5: Types of hydrodynamic control effectors.

Figure 2 . 6 :

 26 Figure 2.6: The "Hydroplane" of Charles de Lambert was fitted with four hydrofoils.

FullyFigure 2 . 7 :

 27 Figure 2.7: Fully submerged hydrofoil vessels.

Figure 2 . 8 :

 28 Figure 2.8: Surface piercing hydrofoil vessels.

Figure 2 . 10 :

 210 Figure 2.10: Outline of a typical hydrofoil used on America's Cup class yachts.

Figure 2 . 11 :

 211 Figure 2.11: Hydroptère during a speed record attempt.

Figure 2 . 12 :

 212 Figure2.12: AC50 boat of Emirates Team New Zealand, racing during the 35th America's Cup.

Figure 2 .

 2 Figure 2.13: Moth IMCA fitted with fully submerged hydrofoils.

Figure 2 . 14 :

 214 Figure 2.14: Typical lift and drag characteristics of a hydrofoil section.

Thickness ( a )

 a Figure 2.15: Hydrofoil section geometries.

Figure 2 . 16 :

 216 Figure 2.16: Boundary layer velocity profile with thickness δ(x).

Figure 2 . 17 :

 217 Figure 2.17: Development of the boundary layer separation: the velocity gradient becomes zero, then negative.

Figure 2 . 18 :

 218 Figure 2.18: Adverse pressure gradient on the low pressure face of a hydrofoil (the dashed line represents the pressure level on the surfaces of the hydrofoil).

Figure 2 . 20 :

 220 Figure 2.20: Free surface waves generated by a lifting surface.

  Figure 2.22: Cavitation.

Figure 2 . 23 :

 223 Figure 2.23: Hydrofoil tested in a cavitation tunnel. Cavitation can be observed on the suction face of the hydrofoil.

  ) is based on the partial discretization of the momentum equations (3.1b) in the face centers. The vector -→ Û and the coefficient C p contain all contributions to the momentum equations except the pressure gradient and the gravity contribution. All terms in -→ Û and C p are rebuilt from their cell

Figure 3 . 1 :

 31 Figure 3.1: Cell face notations.

Figure 3 . 2 :

 32 Figure 3.2: 1D variation of a convected quantity Q in physical space.

  Figure 3.3: Basic schemes in NVD diagram and upstream node.

Figure 3 . 4 :

 34 Figure 3.4: Grid refinement at the water surface, for a surface-piercing hydrofoil.

Figure 4 . 1 :

 41 Figure 4.1: Prediction and variance for a Gaussian process model based on a five-point data set.

Chapter 4

 4 Figure 4.1 gives an example of a GP reconstruction and shows how the variance indicates the regions where this prediction is uncertain. 40 Gaussian process modeling 4.4 Modeling with noisy training sets

Figure 4 . 2 :

 42 Figure 4.2: Interpolation of a noisy dataset by Gaussian process regression.

Figure 4 . 3 :

 43 Figure 4.3: Shape comparison of the SE and Matérn ν = 3 2 kernels for unit and σ.

Figure 4 . 4 :

 44 Figure 4.4: Analytical test functions: Three hump camel (left) and Gramacy-Lee (right).

Figure 4 . 5 :Figure 4 . 6 :

 4546 Figure 4.5: Evolution of the response surface obtained with Gaussian Process Regression with a squared exponential kernel (left) and a Matérn kernel (right). Regular grids with 3, 5, 7 and 9 points per dimension. The test function is Gramacy-Lee.

Figure 4 . 7 :

 47 Figure 4.7: Comparison of the performances of the squared exponential and Matérn kernel in modeling the Gramacy-Lee test function.

Figure 4 . 8 :

 48 Figure 4.8: Comparison of the performances of the squared exponential and Matérn kernel in modeling the Three hump camel test function.

Figure 5 . 1 :

 51 Figure 5.1: A function with multiple minima.

Figure 5 . 2 :

 52 Figure 5.2: Surrogate based optimization with sequential sampling.

Figure 5 . 3 :

 53 Figure 5.3: Evolution of the approximation model with sequential sampling.

Figure 5 . 4 :Chapter 5

 545 Figure 5.4: Sampling methods.

Figure 5 . 5 :

 55 Figure 5.5: Examples of acquisition functions (below, in green) and their corresponding function and variance (in blue).

Figure 5 .

 5 Figure 5.6 shows an example of the difference between variance-based (α = 0) and custom (α = 1) sampling for the Gramacy-Lee test function of section 4.7.2. To obtain the clearest illustration of their difference, response surface creation (β = 0) is performed. While the variance-based sampling distributes its points evenly, the custom sampling detects the two peaks in the lower left-hand corner and concentrates its points there; as a result, the function is represented much better.

Figure 5 . 6 :

 56 Figure 5.6: Distribution of sampling points for Gramacy-Lee (55 points). The colouring shows the acquisition function, the lightest colour represents the highest value.

Figure 5 . 7 :

 57 Figure 5.7: An acquisition function, which presents many local maxima.

Figure 5 . 8 :

 58 Figure 5.8: An explicitly constrained design space.

Figure 5 . 9 :

 59 Figure 5.9: Exclusion zone due to unconverged simulations. This plot represents Xfoil simulations of 2D hydrofoils, where thin foils at high angle of attack fail due to stall (section 2.4).

Figure 5 . 10 :

 510 Figure 5.10: Gramacy-Lee test function, LCB acquisition function. Beta values from left to right : 2.0, 0.5 and 0.25.

Figure 5 . 11 :Figure 5 . 12 : 5 . 8

 51151258 Figure 5.11: Comparison of acquisition function performances for optimization. LCB acquisition function is on the left, custom acquisition function on the right. Test function: Three hump camel. The star represents the global minimum of the function.

Figure 5 . 13 :Figure 5 . 14 :

 513514 Figure 5.13: Comparison of acquisition function performances for response surface creation. Variance based acquisition function is on the left, Custom acquisition function on the right. Test function: Three hump camel.

Figure 6 . 1 :

 61 Figure 6.1: Sequence of operations used to generate a 3-dimensional hydrofoil shape.

Figure 6 . 2 :

 62 Figure 6.2: Spine (red curve) and sections (black) positioned along the spine.

Figure 6 . 3 :

 63 Figure 6.3: The parametric foil spine (left) is composed of four circle arcs. The tip (green segment) is a circle arc with high radius. TEL or Tip Equivalent Length is the length of this tip element. On the right side is a graphical illustration of the tangency continuity used to create the foil spine. Arc 1 and Arc 2 are joined at tangency point A.

Figure 6 . 4 :

 64 Figure 6.4: Normalized curvilinear coordinate system and matching between object space and parameter space.

Figure 6 . 5 :

 65 Figure 6.5: Example of a parameter curve for the chord distribution. The blue dots represent the control points.

Figure 6 . 6 :

 66 Figure 6.6: Interpolation using a B-spline (dashed red line) and a PCHIP interpolator (solid blue line).

Figure 6 . 8 :

 68 Figure 6.8: Illustration of the mixing principle for two different section types.

Figure 6

 6 Figure6.9: The partition of unity used for section blending in the case of two basis functions. The grey area corresponds to the transition length where the two base profile types are mixed. Outside of this grey area, the profiles are not mixed and therefore kept in their original form.

Figure 6 . 11 :

 611 Figure 6.11: Foil attitude variation. On the left, extension variation, on the right cant angle variation.

Figure 6 .Figure 6 . 13 :

 6613 Figure 6.12: Foil "identity card". From top to bottom is presented the following information: chord distribution, camber distribution, twist distribution and section type.

Figure 6 .

 6 Figure 6.14: Groupama C foiling during the 2015 Little Cup on Lake Geneva. (photo: Pierrick Contin).

Figure 6 . 15 :

 615 Figure 6.15: Examples of foil geometry variation.

FzFigure 6 . 16 :

 616 Figure 6.16: Comparative study of the influence of different tip lengths T EL on the vertical stability as a function of the flying altitude. (The "Hull shift" indicates the altitude with respect to a reference point, positive values mean that the foil moves out of the water). The second figure shows the relative stability with respect to a reference foil.

Figure 6 .

 6 Figure 6.17: Views of the hydrofoil designed for Groupama C and installed on the boat. In these pictures, the hydrofoil is inserted in the foil case, in its maximum extended position (photo: Romain Lanos).

Figure 6 . 18 :

 618 Figure 6.18: Hydrofoil sections used for designing the hydrofoil of the Centrale Nantes Hydroproject lightweight boat.

Figure 6 . 19 :

 619 Figure 6.19: The lightweight boat of Centrale Nantes Hydroproject performing a sharp turn maneuvrability test during sea trials. The maneuver induces a strong starboard heel, raising the forward port hydrofoil out of the water (photo: Pierrick Contin).

Figure 7 . 1 :

 71 Figure 7.1: Examples of meshes for monofluid (a) and free-surface (b) simulation. The fine cells at the water surface are created with adaptive refinement. Both meshes have been deformed to rotate the airfoils to the angle of attack which produces the target C l .

Figure 7 . 2 :

 72 Figure 7.2: Response surface creation: drag coefficient as a function of T and C, for constant X = 0.4. Left: variance-based α = 0, right: custom α = 1.

Figure 7 .

 7 Figure 7.2 analyses these response surfaces, showing the predicted drag as a function of T and C, the variance, the cross-validation error and the infill criterion. At low values of

Figure 7 . 3 :

 73 Figure 7.3: Response surface creation, for constant thickness: T = 0.05. Left: α = 0, right: α = 1.

Figure 7 . 4 :Figure 7 . 5 :Figure 7 . 6 :

 747576 Figure 7.4: Optimization: drag coefficient as a function of T and C, for constant X = 0.4. Left: LCB α = 0, β = 0.1. Right: custom α = 1, β = 1.

Figure 7 . 7 :

 77 Figure 7.7: Convergence of the optimizations: minimum C d found and distance of each sampled point to the optimum. Top: custom, left: α = 1, β = 1, right: α = 1, β = 0.5. Bottom: LCB, α = 0, β = 0.1.

Figure 7 . 8 :

 78 Figure 7.8: Convergence of the three-parameter monofluid optimization: minimum C d found and distance of each sampled point to the optimum.

Figure 7 . 9 :Figure 7 . 10 :

 79710 Figure 7.9: Horizontal velocity for X = 0.4 and the extreme T and C, as well as the optimum. Top row: T = 0.12, center row: T = 0.03. Left column: C = 0, right column: C = 0.08. Bottom row: optimum.

Figure 7 .

 7 Figure 7.11: Free-surface effect: geometries of optimized profiles (left) and distribution of the pressure coefficient C p on the optimized profiles (right).

Figure 7 . 12 :

 712 Figure 7.12: Hydrodynamic pressure and free-surface position for optimized profiles with and without free surface. Left to right: monofluid, h/c = 1.0, h/c = 0.25.

Figure 7 . 13 :

 713 Figure 7.13: Influence of the trailing edge: geometries of optimized profiles (left) and distribution of the pressure coefficient C p (right).

Figure 7 . 14 :

 714 Figure 7.14: Influence of the velocity and the lift coefficient: geometries of optimized profiles (left) and distribution of the pressure coefficient C p (right).

Figure 7 . 15 :

 715 Figure 7.15: Surface element for the determination of the vertical and lateral lift slope.

Figure 7 . 16 :

 716 Figure 7.16: Convergence of the 3D hydrofoil optimization: minimum C d found and distance of each sampled point to the optimum.

Figure 7 . 17 :

 717 Figure 7.17: Prediction of drag by the 3D hydrofoil metamodel. Left: isosurfaces (0.115, 0.13, 0.15), right: sampled points.

Figure 7 .

 7 Figure 7.18: Geometry, hydrodynamic pressure and free-surface position (isolines 0.05 and 0.95 of the volume fraction) for the hydrofoils at the corners of the three-parameter design space. Each image shows (left to right) the bottom surface, the top surface, the front. The stable rake, drift (leeway) and C d are given in the figures.

Figure 7 .

 7 Figure 7.19: Geometry, hydrodynamic pressure and free-surface position (isolines 0.05 and 0.95 of the volume fraction) for the optimal hydrofoil.

Figure 7 .

 7 Figure 7.20: Isosurface Q = 20 of the second invariant of the velocity gradient tensor, colored with helicity, for the eight extreme geometries and the optimum.

Figure 7 . 21 :

 721 Figure 7.21: Rake as a function of immersion for twist=0.6 o . The black dot indicates the optimum.

  Pour conclure le chapitre, un hydrofoil tridimensionnel est optimisé. Afin de simuler des hydrofoils à portance verticale et latérale imposée, une procédure de positionnement dynamique est utilisée, qui ajuste les angles d'inclinaison (autour de l'axe y) et de dérive afin d'obtenir les forces voulues. L'hydrofoil optimal, obtenu en faisant varier le dièdre du tip et son vrillage ainsi que la profondeur d'immersion, possède une distribution de pression variant de manière progressive sur l'hydrofoil et dont les deux parties princiaples (tip et shaft) de l'hydrofoil ne génèrent pas de portances opposées. Une étude générale des paramètres de conception pour l'optimisation de l'hydrofoil montre qu'il devrait au moins être possible de faire varier la surface immergée et d'orienter latéralement la portance. Le chapitre se termine par des remarques sur l'optimisation robuste, l'optimisation multipoint pour la navigation au près et sous le vent, et l'analyse directe de la stabilité d'un foil et des risques de cavitation. Cette analyse montre le conflit potentiel entre l'optimisation de la traînée minimale et la stabilité.Le chapitre 8 conclut cette thèse par une analyse des développements ultérieurs envisagés. Optimisation de géométries d'hydrofoils par modèles de substitution construits à partir de simulations RANS Mots clés : Optimisation par modèle de substitution, régression par processus gaussien, simulations RANS, modeleur géométrique, hydrofoils, architecture navale Résumé : Cette thèse présente un framework d'optimisation pour la conception hydrodynamique de forme d'hydrofoils. L'optimisation d'hydrofoil par simulation implique des objectifs d'optimisation divergents et impose des compromis contraignants en raison du coût des simulations numériques et des budgets limités généralement alloués à la conception des navires. Le framework fait appel à l'échantillonnage séquentiel et aux modèles de substitution. Un modèle prédictif est construit en utilisant la Régression par Processus Gaussien (RPG) à partir des données issues de simulations fluides effectuées sur différentes géométries d'hydrofoils. Le modèle est ensuite combiné à d'autres critères dans une fonction d'acquisition qui est évaluée sur l'espace de conception afin de définir une nouvelle géométrie qui est testée et dont les paramètres et la réponse sont ajoutés au jeu de données, améliorant ainsi le modèle. Une nouvelle fonction d'acquisition a été développée, basée sur la variance RPG et la validation croisée des données. Un modeleur géométrique a également été développé afin de créer automatiquement les géométries d'hydrofoil a partir des paramètres déterminés par l'optimiseur. Pour compléter la boucle d'optimisation, FINE/Marine, un solveur fluide RANS, a été intégré dans le framework pour exécuter les simulations fluides. Les capacités d'optimisation ont été testées sur des cas tests analytiques montrant que la nouvelle fonction d'acquisition offre plus de robustesse que d'autres fonctions d'acquisition existantes. L'ensemble du framework a ensuite été testé sur des optimisations de sections 2D d'hydrofoil ainsi que d'hydrofoil 3D avec surface libre. Dans les deux cas, le processus d'optimisation fonctionne, permettant d'optimiser les géométries d'hydrofoils et confirmant les performances obtenues sur les cas test analytiques. Les optima se révèlent cependant être assez sensibles aux conditions opérationnelles. Title : Surrogate-based optimization of hydrofoil shapes using RANS simulations Keywords : Surrogate-based optimization, Gaussian process regression, RANS simulations, geometric modeling, hydrofoils, ship design

  . New time step t = t + ∆t,

	enough to converge when the diagonal dominance is increased by 50% with the help of
	the local time stepping artifact.
	3. Start the iterative procedure with Q = Q 0 ,
	4. If needed, compute the volume fraction for each fluid phase and update the global
	fluid properties,
	5. If needed, compute the turbulent quantities from the field of step 3,
	6. Solve the momentum equations to obtain a new prediction of the velocities,
	7. Solve the pressure equation (3.16) to obtain a new pressure field,
	8. Update the velocity face fluxes (3.15) and correct the velocity components (3.12)
	with the new pressure field,
	9. If the nonlinear residuals are not low enough, go to step 3 and update the iteration
	counter within the time step,
	10. Go to step 2 and update the time, t.
	Concerning the linear solver used in steps 4, 5, 6 for the phase concentration, turbulent
	quantities, and velocity components respectively, about 20 Gauss-Seidel iterations are
	3.4 Numerical framework

  Conclusionpeaks and is a safer choice for optimization, since it is able to detect minima which lie close to maxima.

5.9 72 Chapter 5 Surrogate based optimization

  Boats category in which the prototypes, loaded with 20 kg, represent leisure boats. Both categories have maximum dimensions of 2.50 x 2.50 x 2.00 m.Centrale Nantes Hydroproject is taking part in the contest since 2015 and collaborates with academic and industrial partners for the challenge. Since the boat competing in the lightweight category is a catamaran fitted with hydrofoils, Streamline was invited to participate in the study and design of the hydrofoil set. The chosen configuration consists in one rear fully-submerged foil and two forward surface-piercing hydrofoils (figure6.19).

6.10 Industrial applications of the 3D modeler

Table 7 .

 7 The resulting drag coefficients C d are given in table7.1. Estimating the error on the coarsest grid from such a convergence study is difficult, but the numerical uncertainty in the drag on the coarse meshes is probably at least 15%. Since the main objective of this section is to test the optimization procedure, this is acceptable. Each computation takes about 15 minutes on a 16-core Xeon workstation.

	Coarse Medium Fine	Veryfine
	α (deg.) C d • 10 3 Difference veryfine 6.3% 0.0156 0.159 8.6375 8.3448 2.7%	0.185 8.1475 8.1272 0.221 0.3% 0.0

1: Convergence of the computed drag coefficient C d

Table 7 . 2 :

 72 Optimized profile parameters and drag coefficients, depending on the distance h to the free surface. Drag coefficients C d • 10 3 are shown for the optimized condition, as well as for the other immersions.

	.25

Table 7 .

 7

	Case	T	C	X	C d blunt C d sharp
	Blunt TE (baseline) 0.0300 0.0458 0.378 7.764	7.657
	Sharp TE	0.0302 0.0450 0.401 7.771	7.640

3: Influence of the trailing edge thickness. Drag coefficients C d • 10 3 are shown for the optimized condition, as well as for the other TE shape. 7.2 Shape optimization of 2D profiles 107

Table 7 . 4 :

 74 Influence of the velocity and the lift coefficient. Drag coefficients C d • 10 3 are shown for the optimized condition, as well as for the other velocity.

  Présentation de la thèse et de son contexteEn ingénierie, l'optimisation de forme consiste à chercher les variantes géométriques les mieux adaptées à leur rôle. En architecture navale, l'optimisation de forme et notamment des carènes, est ancienne. Elle fut d'abord réalisée par empirisme, grâce à l'expérience acquise par les constructeurs de navires au fil des générations, développant des types de navires adaptés à leurs utilisations et contraintes d'exploitation ainsi qu'aux conditions météorologiques rencontrée. Ces optimisations tenaient également compte des contraintes liées aux méthodes de construction alors disponibles.

  Afin d'illustrer les capacités du modeleur, deux de ces applications sont présentées: la conception des foils principaux d'un catamaran de course de type Classe C pour Groupama Sailing Team et la conception d'un jeu de foil complet pour un drone a propulsion éléctrique pour un projet pédagogique développé au sein de Centrale Nantes.Les éléments constitutifs de la boucle d'optimisation et décrit dans les chapitres précédents sont assemblés dans le chapitre 7. Dans ce chapitre, la procédure adaptative d'échantillonnage et d'optimisation a été appliquée à des cas de conception de géométries d'hydrofoils avec simulations RANS. Tout d'abord, les tests du chapitre 5 sont répétés pour la création de surfaces de réponse et l'optimisation de sections de foil 2D. Par la suite, les sections 2D optimales sont étudiées et la dépendance de l'optimum par rapport aux conditions de fonctionnement est discutée. Pour terminer, un hydrofoil tridimensionnel généré par le modeleur géométrique est optimisé et la forme optimale est analysée en termes de stabilité et de risque de cavitation.La première série de tests sur des profile d'hydrofoils (en deux dimensions) confirme les résultats obtenus dans le chapitre 5 : il est plus sûr de baser l'échantillonnage sur la fonction d'acquisition custom que sur une fonction d'acquisition basée sur la variance.

		Applications
	142	Extended summary

2.3 Hydrofoils in use

Chapter 2 Hydrofoils and foilers

Chapter 3 Computational Fluid Dynamics

If the meshes for the different computations are made with an unstructured mesh generator, the placement of the cells in each grid is semi-random so the errors due to the grid are uncorrelated between computations. For structured grids, or if the different simulations are performed on the same grid, the numerical errors cannot be seen as uncorrelated noise.4.4 Modeling with noisy training sets

5.6 Searching the acquisition function maximum

Chapter 5 Surrogate based optimization

These coordinates express a Cartesian point (x, y, z) as (wx, wy, wz, w). For the operations performed here, w is kept fixed at 1 so the homogeneous coordinates are simplified as (x, y, z, w).6.5 Section modification

Spine alignmentNext, a rotation is applied to get the section in a position normal to the spine. Figure6.2 shows how the resulting sections are all rotated differently, to follow the

Chapter 6 Geometric modeling

https://www.opencascade.com/

http://www.pythonocc.org/ 86 Chapter 6 Geometric modeling

6.8 Surface generation

http://www.tspeer.com/Hydrofoils/h105/h105.htm 92 Chapter 6 Geometric modeling

7.2 Shape optimization of 2D profiles

Chapter 7 Applications

7.3 3D hydrofoil optimization

7.4 Conclusion

Conception de navires et d'hydrofoilsLes navires sont le plus souvent optimisés dans le but de minimiser la puissance nécessaire pour atteindre une vitesse donnée, et ce en réduisant la traînée hydrodynamique. Les autres objectifs peuvent consister à s'assurer de la stabilité du navire, limiter ses mouvements dans les vagues ou encore offrir des capacités de manoeuvrabilité adaptées. Les coûts d'exploitation d'un navire étant souvent élevés (relativement à son coût de fabrication), de faibles gains de performance peuvent permettre des gains financiers conséquents. Cependant, contrairement aux voitures ou aux avions par exemple, la plupart des navires sont des modèles uniques ou construits en très petites