Sports | Watersports » Paul Kerdraon - Design and models optimisation of a sailing yacht dynamic simulator

Degree Project in Naval Architecture Second Cycle Stockholm, 2017 Design and models optimisation of a sailing yacht dynamic simulator Paul KERDRAON KTH Royal Institute of Technology Centre of Naval Architecture Contents Acknowledgements 4 Abstract 5 Index of variables 6 Index of acronyms 8 Introduction 9 I Performance of a sailing yacht and dynamic simulation 10 1 The physical model of sailing 10 1.1 Apparent wind 10 1.2 Sailing boat equilibrium 10 2 Yacht design and performance prediction 11 2.1 Velocity Prediction Programmes 11 2.2 Limits of VPP studies 11 3 A dynamic simulator 3.1 Objectives 3.2 State of the art 3.3 Principle 3.4 Main hypothesis 3.5 Coordinate systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 12 12 13 14 15 4 Studied ships 17 5 Presented work 17 II 18 Sailing boats dynamic simulator 1 Equations of motion 18 1.1 General presentation 18 1.2 Linearised equations: added mass and damping 19 1.3 Symmetry considerations 20 2 A wave model for the simulator 20 2.1 Interest 20 2.2 Model and assumptions 20 3 Steady and unsteady waves 22 3.1 Ship in waves 22 3.2 Wave-form

parameter 23 3.3 High forward speed 24 4 Appendages 4.1 Model 4.2 Appendages dynamics 4.3 Added mass and damping 4.4 Dynamics in waves 4.5 Angles of attack in waves 4.6 Wave loads 4.7 Other dynamical aspects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 26 26 27 27 27 28 28 5 Hulls 5.1 Model 5.2 Restoring force

5.3 Wave loads 5.4 Added mass and damping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 29 29 30 30 6 Hulls hydrodynamic coefficients 6.1 Numerical methods 6.11 Seakeeping prediction 6.12 Panel methods 6.2 Discussion 6.3 Determination from CFD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 31

31 33 33 34 . . . . . . . . 7 Aerodynamic response 35 III 36 Stability of offshore sailing foilers 1 Foil 1.1 1.2 1.3 1.4 1.5 and stability, state of the art Foil principle . Conditions of use and class rules . The question of the stability . Pitch stability and airplane equivalence Heave stability and foil design . 2 Dynamic considerations 2.1 Initial interest 2.2 Pitch study 2.3 Heave: one-parameter study 2.4 Heave: multi-parameter study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 36 37 38 38 39 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 41 41 43 44 Conclusion 47 Bibliography 48 A Heaving foil in head sea 49 3 Acknowledgements I would like to express my deep gratitude to the whole VPLP team, especially Vannes’ agency, Vincent, Charlotte, Xavier, Anne, Quentin, Daniele, Xavier, Nicolas, Adrien, Antoine, Gabriel, Erwan, Patricia and Philippe. Thanks for having welcomed me so generously and taken time to show and explain your work and projects. Through the numerous conversations we had and the warm atmosphere which reigns there, I have been able to learn a lot, to enlarge my scientific culture as well as my general knowledge. I would like to thank in particular Xavier and Adrien for their kind supervision of my master thesis. Our discussions have been of great help and you were always able to

lend a helping hand or to propose new ideas when I had the impression to be in a dead-end. 4 Abstract This master thesis is concerned with the design and optimisation of a dynamic velocity prediction programme for high performance sailing yachts. The simulator uses response surfaces methods, maximising the computational efficiency Insights are given on dynamic theoretical aspects and models are discussed with respect to the studied ships, 100 feet offshore racing trimarans. The consistency of the currently used models is assessed, and feasible solutions and methods are proposed and implemented as solutions to the identified defaults. The simulator is subsequently used on a concrete case with the objective of assessing the effects of precise geometrical features of appendages (rudders, foils and boards) on the ship dynamical stability properties. Area, extension, cant angle and tip-shaft angle in particular are studied Tests cases are developed and multi-parameters studies carried

out. Dynamic results are compared to usual static stability criteria Simulations show some inconsistencies of dynamic responses with the behaviour expected from the static criteria. Such processes allow better understanding of the design and scantling of the studied appendages. 5 Index of variables Variables Aij Coefficient of the hydrodynamic force in ith direction in phase with body acceleration in the jth direction (Added mass coefficient) AWA Apparent wind angle AWS Apparent wind speed Bij b Cij Coefficient of the hydrodynamic force in ith direction in phase with body velocity in jth direction (Damping coefficient) Ship max beam Coefficient of the hydrostatic force in ith direction due to body motion in jth direction (restoring force) c Appendage mean chord cg Group speed cw Wave velocity (= ω/k) Fw Wave hydromechanical loads g Acceleration of gravity h Seabed depth Iij Mass moment or product of inertia i Imaginary unit k Wave number (= ω 2 /g) L

Boat length li Appendage lever-arm in pitch Mij Mass matrix m Mass of the body n Unit normal vector of body/free surface boundary, directed into the fluid p Pressure p0 Atmospheric pressure pD Dynamic pressure pH Hydrostatic pressure S Appendage area s Appendage span T Period t Time ∗ Pseudo-time variable t TWA True wind angle TWS True wind speed V Translational velocity vector Vb Boat mean speed α Wave angle in earth-fixed coordinate system αc Kelvin angle αw Wave slope β Shaft-tip angle Γ Appendage cant angle 6 ∆ Displacement (mass) δ Wave velocity amplitude to boat speed ratio (= ζ0 ω/Vb )  Phase shift ζ Surface elevation ζ0 Wave amplitude η Leeway angle ηi ith degree of freedom θ Pitch angle λ Wavelength µ Angle of wave incidence ν Fluid kinematic viscosity ρ Specific mass of water ϕ Roll angle φ Velocity potential ψ Yaw angle ; heading Ω Angular velocity vector ω Wave

frequency - Circular frequency of motion ωe Circular frequency of encounter ∇ Displacement (volume) Non dimensional numbers Fn √ Froude number (= Vb / gL) Istab Heave stability index Re Reynolds number (= Vb L/ν) Vstab Horizontal tail volume coefficient γ Wave-form parameter (= ωVb /g) p Frequency number (= ω L/g) ξ Operators ∆ Laplace operator ∇ Gradient operator Coordinate systems R0 = (x0 , y0 , z0 ) Earth-fixed coordinate system Rb = (xb , yb , zb ) Body-fixed coordinate system Re = (xe , ye , ze ) Element-fixed coordinate system Rf = (xf , yf , zf ) Fluid coordinate system Rt = (xt , yt , zt ) Linear seakeeping coordinate system 7 Index of acronyms AVL Athena Vortex Lattice AWA Apparent Wind Angle AWS Apparent Wind Speed BEM Boundary Elements Methods CE CFD Centre of Effort Computational Fluid Dynamics CG Centre of Gravity CLR Centre of Lateral Resistance DES Detached Eddy Simulation DNS Direct Numerical Simulation

DVPP Dynamic Velocity Prediction Programme EFD Experimental Fluid Dynamics FSI Fluid Structure Interactions GSM Green Function Method HSST High Speed Strip Theory IMOCA LES International Monohull Open Class Association, 60 ft monohull Large Eddy Simulation ORMA Ocean Racing Multihull Association, 60 ft trimarans RANSE Reynolds-Averaged Navier Stokes Equations (also RANS) RSM Rankine Source Method STF Salvesen - Tuck - Faltinsen method, conventional strip theory TWA True Wind Angle TWS True Wind Speed VPP Velocity Prediction Programme 8 Introduction VPLP is a young and international team of naval architects and yacht designers. It was founded in 1983 by Marc Van Peteghem and Vincent Lauriot Prévost and has since then proved successful in a number of design concepts, including multi and monohulls, power and sailing concepts, from high performance sailing yachts to production cruise boats and super yachts. Among other projects, they collaborated with Alain

Thebaud on the Hydroptère, first sail craft to ever reach more than 50 knots over 1 nm in 2009. Figure 1: The Hydroptère, breaking the speed record. Figure 2: Oracle’s trimaran, USA 17. G. Martin-Raget G. Grenier/Oracle Racing After having been contacted by Russell Coutts in 2007, and after numerous optimisations and improvements, they designed Oracle’s USA 17, winner of the 33rd America’s Cup, a light and powerful 90’ trimaran with a canting wing-sail. More recently, this winter, their 31 m trimaran Sodebo IV sailed by Thomas Coville broke the solo round the world record in 49 days, 3 h and 7 mn while team IDEC lead by Francis Joyon on a VPLPdesigned trimaran won the Jules Vernes trophy (fastest circumnavigation of the world) in only 43 days, 23 h and 30 mn with an average speed over ground of 26.9 kn In monohull design, the 2016 Vendée Globe has seen the participation of twelve 60’ IMOCA designed by the team, among which 6 new generation boats with foils imagined

in collaboration with G. Verdier Armel Le Cléac’h, winner of the edition, sailed one of them. Figure 4: Banque Populaire VIII, winner of the 2016 Vendée Globe. T Martinez / BPCE Figure 3: IDEC, training off Belle-Ile, Brittany. JM. Liot / DPPI / IDEC 9 Part I Performance of a sailing yacht and dynamic simulation 1 The physical model of sailing This first part aims at recalling the main concepts and models involved in the current understanding of sailing ships. 1.1 Apparent wind Figure 5: The wind triangle. Due to its forward motion, the wind perceived by the moving boat is not the wind that a motionless point (weather station for instance) would measure. The latter, the wind blowing in the earth coordinate system, is usually referred to as True Wind (TW), defined by its norm, the True Wind Speed (TWS), and its angle, the True Wind Angle (TWA). The interaction of the True − Wind with the boat speed (BS) creates an Apparent Wind (AW) (see Figure 5) which is received

by the sails and defined in terms of Apparent Wind Angle (AWA) and Apparent Wind Speed (AWS): −− −− − AW = T W − BS. (1) Hence, the faster a ship is sailing, the closer she is to the wind - that is with a smaller angle. 1.2 Sailing boat equilibrium As one knows, sailing boats are propelled by the wind force on the sails. As for an airplane wing, the difference of pressure resulting from unequal velocity of the air circulating on each of the sail sides creates a lift force perpendicular to the sail main chord. Figure 6: The forces equilibrium of a sailing boat. Top and rear view CE: centre of effort, CG: centre of gravity, CLR: centre of lateral resistance. However, unlike a propeller, this force is not directed along the boat axis but has a lateral component (see Figure 6) which tends to make the boat drift on the water. Hence the actual ship’s speed vector is not directed along the boat centreline, but shifted of an angle called leeway in the direction of the

aerodynamic side force. 10 To limit this phenomenon – which is stronger as the boat gets closer to the wind – a high number of sailing boats are equipped with boards or keel blades which are profiled wings as well and generate a lift force counteracting this aerodynamic side force. Those appendages, the rudder(s) and the hull(s) are also subjected to a hydrodynamic drag (viscous friction, form drag, wave making) which when an equilibrium is reached, are equal to the aerodynamic driving force. Furthermore, as the centre of the aerodynamic force – the centre of effort (CE) is usually located at a point much further from the centre of gravity (CG) of the boat than the centre of the hydrodynamic force – the centre of lateral resistance (see Figure 6), a boat does not sail upright but heeled and trimmed so that aeroand hydrodynamic forces and their respective lever-arms are balanced. For more insight, the interested reader is referred to Fossati [2009]. 2 2.1 Yacht design and

performance prediction Velocity Prediction Programmes As it has been briefly exposed, the motion of a sailing yacht is a complex issue. Numerous design parameters are strongly linked and the alteration of one calls for changes in others. As an example, an increase in sail area will lead to an increase in heeling moment, which in turn can be limited by using a longer keel blade or a heavier bulb at the end of the blade, increasing consequently the ship’s displacement. Sailing boats design is therefore a thorough study of complex compromises. With the development of computational tools, numerical methods have emerged to assess the performance gains and losses from those design choices. The first Velocity Prediction Programme (VPP) was developed at the Massachusetts Institute of Technology in the 70s (see Kerwin and Newman [1979]) Its aim was to predict the speed of a sailing boat based on hydrodynamic and aerodynamic forces models and a numerical optimisation process. Today VPPs are

widely used among sailing boat architects, in a broad range of use, for instance: – To assess design or structural choices of hull forms, – To determine the location (especially longitudinal) of appendages, or their implant angles on the ship, – To evaluate different sail plans or rigs, – To compare given sailing boats, – To compute polar curves for use in routing software. VPPs results are often presented in the form of polar curves which show for given True Wind Angle (TWA) and True Wind Speed (TWS) the resulting ship’s speed. Depending on the use done with the VPP, different parameters can be used as variable inputs to be optimised. However, the ship’s speed is generally the optimised output variable From the different geometrical characteristics of the hull(s), appendages and sail plan and a mathematical model, the equations of motion are numerically solved. The used models can be obtained from numerous methods: Computational Fluid Dynamics (CFD), panel code, model

tests, vortex lattice methods, etc. 2.2 Limits of VPP studies VPPs output data are thus corresponding to equilibrium states, in steady conditions. At sea however, it is quite improbable to reach such conditions as wind and waves have - in reality - changing speeds and directions, as the crew moves and acts on the ship. The ship is therefore permanently subjected to accelerations in all six degrees of freedom, as a consequence of external (wind, waves, currents) and internal (crew position, rudder angle) inputs variations. A VPP cannot show the yacht performance in such states In this way, designs that proves good on the paper based on steady optimisation can turn out to be ill-adapted to dynamic and variable conditions. A simple example would be a hull form that has a low drag in steady conditions but cannot make its way through waves Experience and know-how are useful qualities to face those issues 11 3 3.1 A dynamic simulator Objectives In the competitive field of high

performance sailing, constant improvements and optimisations are needed. Understanding the dynamical behaviour of a given design under unsteady conditions may well provide the advantage needed to make the difference on the water. The current development of foils (see Part III), for which stability is a fundamental question in the design process, is also very demanding on dynamical responses. A better knowledge of the effective limits of stability can lead to substantial performance progress as stability is generally paid with drag increase. VPLP Design has therefore decided to launch the development of a dynamic simulator. Instead of solving a system of equations to determine an equilibrium under a set of constraints, such a programme is a timedomain tool that solves at each time step the equations of motion by integrating the loads. Figure 7: Front view of the simulator visual interface. 3.2 State of the art Following the wide spread of VPP use, the idea of modelling unsteady

situations has made its way. The concept of solving the equations of motion in the time domain has interested many researchers Initially, many Dynamic VPP (DVPP) were based on precomputed models through extended VPP, with the aim of computing manoeuvres for instance. The recent years have shown the apparition of dynamic VPP based on RANSE (Reynolds-Averaged Navier Stokes Equations) solver coupled with optimisation algorithm (see especially Roux et al. [2008]) and Böhm [2014]) For the time being they mainly aim at computing equilibrium – as static VPP – but with more accuracy and flexibility. Such methods enable the simulation of unsteady behaviour in response of perturbations, which is of great interest for seakeeping studies for instance. However, the computational power needed to run such codes demands access to powerful clusters and remains beyond the reach of companies such as VPLP. 12 From coupling two different solvers (potential code for aerodynamic forces and RANS

solver for hydrodynamics one in the case of Roux et al. [2008]) to entire RANS methods (see Korpus [2007]), the use of computationally demanding solvers and methods widens The latter do not fully simulate body motion but find optimum configurations by interpolating various hull and sails state and meshes results. The field of dynamic yachts performance simulation is gaining a wide interest and the range of existing methods soars. The choice between Experimental Fluid Dynamics (EFD), potential codes and viscous CFD solvers is a balance of financial, computational and time costs. As an order of magnitude, the VPP polars from the RANSE-VPP of Böhm [2014] have needed 1344 h in total to compute 60 optimum boat speeds. Having already developed a wide range of methods and response surfaces for static VPP and in regards to the computational available power, VPLP’s DVPP is a coefficient-based VPP, which as a conventional VPP computes loads by calling precomputed data. 3.3 Principle The

simulator is implemented on MATLAB Simulink. Simulink is a MATLAB integrated block diagram environment providing a simulation oriented tool with a high modularity potential As previously stated, the simulator uses precomputed models: multidimensional fitting is done on a wide range of data with inputs being the state of the considered element (immersion, speed, angles, etc.) and outputs being the efforts on this element. The data are coming from: – A commercial CFD software environment (FINE Marine, NUMECA) for hulls. – Athena Vortex Lattice (AVL) for appendages. AVL is a free program developed by the Massachusetts Institute of Technology (MIT) and based on an extended vortex lattice method to compute lift and induced drag of lifting surfaces. – Wind tunnel tests for the aerodynamic efforts. To increase accuracy and numerical stability, the models are corrected in extreme situations such as a hull or appendage getting out of the water. At each time step, all the models are

interrogated by the simulator to compute the global efforts (forces and moments) on the system which is considered as a rigid body. According to Newton’s Second Law, those efforts can then lead – through successive integration – to the speeds (translation and rotation) and positions of the system at the current time step (see Figure 8). Figure 8: Working principle of the dynamic simulator. 13 3.4 Main hypothesis The hypothesis on which the whole simulator is built can be sorted in two categories: those being inherent to the models used for the precomputed data, and those proper to the simulator. In the first category, one important assumption is the no-deformation hypothesis. It is widely known now that this is not true, especially for the highly loaded appendages such as the hydrofoils. Many researchers are currently working on Fluid-Structure Interactions (FSI) experimentations and models to account for the potentially large shape changes occurring under load, and which

can have substantial consequences on the angle of attacks for example. In this respect, sail models are a particular burning issue. Their fully different behaviour up- and downwind, the large deformations they encounter and the difficulty to account accurately for all the existing adjustment make extremely difficult the modelling of the aerodynamic loads. A simple twist/flat model is used here In the second category, a direct consequence of the described method is that each element is considered independently from the others and no coupling is taken into account. In particular, no flow perturbation is computed. For instance, the fluid encountered by the rudders is supposed unperturbed by the upstream parts of the boat, which obviously is not true. This also supposes that there is no interaction between hulls. In the studied cases, the potentially affected parts are the central and leeward hulls, but not the windward one as due to the ship dihedral angle, it is usually not in the

water. Regarding unsteady generated waves, this is a reasonable hypothesis for sufficient high speed. p The transverse radiated waves are moving with the group velocity cg = 1/2 g/k (see Section II.3) They do not reach the other hull at distance b/2 during the time needed for the ship to advance of L as long as: Vb ≥ Lg 2L × cg = . b bωe (2) However, as far as the steady waves are concerned (the ship wash), the assumption is more involving. It can be shown (see Wave Resistance and Wash, in Faltinsen [2005]), that the angle of the wash, named Kelvin angle, is αc ≈ 1/3 for high speed yachts in deep water. From one hull, it does not reach the other one if the length of the ship verifies: b L≤ . (3) 2 tan αc A numerical evaluation shows that such a condition is not respected for the studied ship. Approximations are also made in order to be able to solve and integrate the loads. In this way, the attitude of the ship is considered constant throughout each time step so that

loads can be computed and integrated. This is valid as long as the time step is sufficiently small in comparison with all other characteristic times of the system. Besides, it is obvious that the fluid velocity direction and magnitude are not uniform over a given element of the boat, especially when she is in unsteady motion. To be able to use precomputed data, a reference point is fixed and the fluid velocity at this point is used to compute the loads. Technically it is possible with such a code to simulate a full 6 degrees of freedom system. However to do so, one actually needs a rudder angle command law or a manual user-simulator interaction. The simulation presented herein are therefore only made in degree of freedom, with fixed heading. 14 3.5 Coordinate systems Several coordinate systems are used in the simulator and in this report. They are defined hereafter The earth-fixed coordinate system R0 is the reference frame, assumed Galilean. The third coordinate vector z0 is

directed upward, while x0 points to the North and consequently y0 to the West Its vertical origin z0 = 0 corresponds to the free surface at rest while in the horizontal plane, it is defined as the ship starting point (x0 , y0 ) = (0, 0). The boat coordinate system Rb is the body-fixed coordinate system. xb is directed along the ship centreline , yb points to port, and zb completes the system, upwards Its origin is defined as the ship centre of gravity. Figure 9: Coordinate transformation from R0 to Rb . The ship attitude at a given time is expressed by the three angles describing the orientation of Rb in R0 . Those are taken as (z,y,x) Tait–Bryan angles with intrinsic rotations and therefore corresponds respectively to heading, trim and roll (Figure 9). Each element of the ship has also its own coordinate system Re , used to compute the loads. Those coordinates are defined as the element main directions. Similarly, the angles to orientate Re in Rb fully describe the implantation

of each appendage on the boat, and are respectively named yaw, rake and cant (Figure 10). When computing the loads, a local fluid coordinate system Rf is defined, with the first component vector xe in the direction of the fluid local velocity, ze upwards and ye completing the right-handed coordinates. The rotation angles enabling to go from Rf to Re define the effective angles of attack of the fluid on the appendage (see Section II.42) Figure 10: Appendages implantation angles. 15 Finally, in the case of the linear seakeeping model, the coordinate system described by Salvesen et al. [1970] is used. This frame, Rt , is fixed to the translating ship mean position and enables the description of the ship’s small oscillations about this position. xt is in the ship direction, yt to port and zt upwards. Its vertical origin is taken in the plane of the undisturbed free surface so that the coordinates of the centre of gravity in Rt are (0, 0, zG ). Figure 11 shows the coordinate

system and the definition of the six degrees of freedom (ηi ) describing the ship’s translatory and angular displacements. Figure 11: Linear seakeeping coordinate system. Original picture from Wikipedia. 16 4 Studied ships For the time being, the simulator has only been used on large multi-hulls. The main particulars of two of them are given in Table 1 and Figures 12 and 13. They are equipped with three T-rudders (rudders with an horizontal lifting surface), two foils and one centreboard. Ship Macif 100 Banque Populaire IX Length 30 m 32 m Beam 21 m 23 m Displacement 14.5 t 15 t Draught 4.5 m 5m Air draught 35 m 38 m Sail area upwind 430 m2 610 m2 Sail area downwind 650 m2 890 m2 Architect VPLP VPLP Shipyard CDK Technologie CDK Technologie 2015 2017 Ultime Ultime Launching Class Table 1: Macif 100 and Banque Populaire IX main particulars. Data from ultimboat.com ça m’a bien motivé pour aller pédaler Figure 12: Macif 100. JM.

Liot/DPPI/MACIF Figure 13: Banque Populaire IX. Voile Banque Populaire This kind of trimarans sails at high speed (more than 30 kn), corresponding to a Froude number of more than 0.9 Designed with very slender hulls, they have very low beam-to-length, and draft-to-length ratios 5 Presented work The work herein presented is made of two main parts. In the first one, improvements are brought to the simulator – especially with regard to dynamic behaviours and wave models – and validated In the second part, the simulator is used to study the dynamic response and stability of sailing foilers in full fly and hybrid modes (foiling with a small percentage of hull in the water). 17 Part II Sailing boats dynamic simulator 1 1.1 Equations of motion General presentation The starting point for deriving the 6-degree-of-freedom equations of motion of a sailing boat is Newton’s Second law. For translation motions, at the centre of gravity, it is expressed in the earth-fixed

coordinate system R0 as:   d mV̇ = F, (4) dt where m is the boat’s mass, V̇ her speed vector and F the total forces she is subjected to. It can be decomposed in hydrodynamic, aerodynamic and other components (propeller for instance): F = Fhydro + Faero + Fext − mgz0 . (5) Fhydro results from pressure and viscous stresses on the wetted surface of the hulls and appendages. Faero is mainly the loads exerted by the wind on the sail plan. Similarly, for angular momentum the equilibrium in the earth fixed coordinates gives: ¯ dIΩ dt = MCG . (6) R0 Ω is the body angular velocity vector and I¯ is the matrix of moment of inertia of the ship at CG with respect to Rb .   Ixx −Ixy −Ixz   I¯ = −Iyx Iyy −Iyz  (7) −Izx −Izy Izz The ship being usually symmetric with respect to the centreplane xz, all diagonal terms involving the ycoordinate are zero or negligible. As previously, the total moment on the boat around the centre of gravity M can be

expressed as: M = Mhydro + Maero + Mext . (8) As I¯ remains constant in Rb (rigid body with mass assumed constant), it is therefore of interest to rewrite (6) in this coordinate system: ¯ dIΩ ¯ =M . + Ω ∧ IΩ (9) CG dt Rb Finally, equations (4) and (9), form the 6 equations of motion of the ship: ( mV̈ = F, ¯ ¯ I Ω̇ + Ω ∧ IΩ = M . CG 18 (10) 1.2 Linearised equations: added mass and damping In this Section, the ship is considered as translating at boat speed Vb . The equations of motion are derived in the coordinate system Rt introduced in Section I.35 The centre of gravity, G, is located at (0, 0, zG ) from the origin O0 of Rt . The small linear and angular displacements relative to the translating ship mean position are (ηi )i∈{1,.,6} , respectively surge, sway, heave, roll, pitch and yaw (4) remains valid in R0 , however, as in this case the centre of gravity is not the origin of the moving coordinate system, its expression in Rt must take into

account the motion of the centre of gravity with respect to O0 . dVG dt = R0 dVG dt + Rt dVO0 dt + Rt dΩRt /R0 dt
G ∧ O0 G + ΩRt /R0 ∧ ΩRt /R0 ∧ O0 G + 2ΩRt /R0 ∧ VR , (11) t R0 where the first and second term of the right-hand side are respectively η̈ t = (η̈1 , η̈2 , η̈3 ) and 0, while the derivative of the angular velocity is: dΩRt /R0 dt = R0 dΩRt /R0 dt + ΩRt /R0 ∧ ΩRt /R0 = η̈r . (12) Rt Having assumed small motions, one can rewrite (11) neglecting all second order terms, which in (4) yields: mη̈t + mη̈r ∧ O0 G = F. (13) As far as angular momentums are concerned, a similar development can be made yielding: I¯t η̈r + mO0 G ∧ η̈t = M. (14) where I¯t and M are respectively the ship’s moments of inertia and loads about the axis of Rt . Compiling (13) and (14) yields :   ¯ η̈ = F , M̄ (15) M ¯ is the generalized mass matrix of the ship: where M̄  m 0  0 m   0 0 ¯ M̄ =   0 −mz G 

mzG 0 0 0 0 0 0 −mzG m 0 0 I44 0 0 0 −I46 mzG 0 0 0 I55 0  0 0   0   −I46   0  I66 (16) In the linear seakeeping model, the loads F and M are usually expressed as functions of η, η̇ and η̈ in the form:   F = −Aη̈ − B η̇ − Cη + Fw . (17) M A, B and C are respectively the ship’s added mass, damping and hydrostatic matrices, and Fw the vector of the wave motion hydromechanical loads. In this form the equations of motion appear fully similar to any simple mechanical system. But it hides a major difference: due to the free surface along which it sails, the added mass and damping coefficients of a surface ship are no longer constant but functions of the frequency of motions. Those linear equations of motion aim at introducing the concept of added mass and damping. In reality the dependencies of the loads on the degrees of freedom are more complex and not always linear. For instance, as one knows appendages lift force is quadratic in

speed. Hereafter are presented the methods and models used to compute the loads. Added mass corresponds to the additional inertia due to the surrounding fluid which is accelerated as the ship moves, while damping stands mainly for the energy losses due to viscosity and radiated waves. 19 1.3 Symmetry considerations For usual ships, symmetric with respect to their centreplane, there is only a weak coupling between in-plane (surge, heave, trim) and out-of-plane (sway, roll, yaw) motions. That is to say, all damping, added mass and restoring coefficients involving a coupling between surge, heave and pitch on one hand and sway, roll and yaw on the other hand are zero. This enables the simplification of Equation (15) in two uncoupled sets of three equations. For a sailing boat however, such a simplification is more questionable as it moves heeled and drifting sideways, breaking the symmetry. It is therefore expected to observe some coupling between the two groups of motion. As long

as leeway and heel remain small – as they generally ought to be for good performance – the coupling terms should however remain of small order of magnitude. To illustrate this point, restoring coupling terms are derived and compared in Section II.52 for a hull with and without a leeway angle 2 2.1 A wave model for the simulator Interest Unlike the Americas Cup, offshore sailing competitions bring sailors far from the shores, and - especially in the Southern oceans - which are wavy areas. It is therefore important that the designed boats have good seakeeping properties Using the simulator to better understand a given ship design response in waves is therefore a interesting opportunity. 2.2 Model and assumptions To account for the effect of waves on the ship’s parts, a pressure and velocity distribution must be derived. Assuming the fluid is incompressible, mass and momentum conservation laws write: ∇ · u = 0, Du = −∇p + µ∇2 u − ρgz0 , Dt where D/Dt = ∂/∂t +

u · ∇ stands for the material derivative. ρ (18) (19) At the ship scale, the viscous forces are negligible with respect to the inertia terms, this is shown by computing the Reynolds number: Inertia Vb L Re = = (20) Viscous ν with Vb and L the ship’s characteristic speed and length, and ν the fluid kinematic viscosity, 1.2 10−6 m2 /s for salted water at 15‰. For the ships presented earlier in this report, this leads to Re ≈ 108 However, viscosity remains substantial near solid boundaries and in the wake. From now on, viscosity is neglected, preventing any study of unstable breaking waves. Moreover, as such a fluid starting from rest is irrotational, we can define the fluid velocity potential φ: u = ∇φ 20 (21) In those conditions, Equations (18) and (19) leads respectively to Laplace and Bernoulli equations: ∆φ = 0, (22) ∂φ 1 p 2 + (∇φ) + + gz0 = constant. ∂t 2 ρ (23) Let’s now consider a three dimensional domain, with Cartesian coordinates as

shown in Figure 14. z is pointing upwards and has its origin at the calm water surface. The free surface at time t is given by z0 = ζ (x0 , y0 , t), while the seabed is defined by z0 = −h (x0 , y0 ). n indicates the unit normal vector to the considered surface, pointing into the water. Figure 14: Wave model coordinate system. Along the seabed, there is no fluid going through the boundary, this yields the following boundary condition: ∇φ · n = 0 at z0 = −h (x0 , y0 ) . (24) The free surface boundary condition is more involving as its position is itself an unknown. There is a similar condition on normal velocity, called kinematic free surface condition, stating that the fluid normal velocity at the boundary is equal to the surface’s own velocity:   ∂ζ z0 · n. ∇φ · n = (25) ∂t This condition can also be expressed as D/Dt (z0 − ζ(x0 , y0 , t)) = 0 on ζ(x0 , y0 , t). A Taylor expansion near z0 = 0 yields: ∂φ ∂ζ ∂ζ ∂φ = + . (26) ∂z ∂t ∂x ∂x

The location of the free surface adds an unknown to the system and an additional condition must be found. This is done through the dynamic boundary condition, which states that at the free surface, momentum is conserved. Surface tension being negligible above a few 10−2 m it is not to be considered here, and the focus will be on gravity waves, driven by kinetic and gravitational potential energies. As far as normal forces are concerned, the equilibrium leads to a pressure balance on each side of the free surface: the fluid pressure given by Bernoulli equation in (23) must equals the atmospheric pressure p0 . 1 ∂φ 2 = − (∇φ) − gζ. ∂t 2 (27) The full development of the wave model is not to be presented herein, but may be found in Garme [2011]. To simplify further the problem, the waves are supposed linear, that is to say their amplitudes are small compared to their wavelengths: ζ0  λ. This is a fairly acceptable assumption as most of the waves that affect a ship

behaviour have rather low amplitude-to-wavelength ratio, even if foaming and breaking waves are present. This allows to neglect the quadratic term in Bernoulli’s equation. Subsequently, the superposition principle applies and a sum of solutions to the problem remains solution. Assuming the seabed depth is constant and using a separation of variables method, the wave velocity potential is derived as: φ (x0 , y0 , t) = ζ0 g cos (kx0 − ωt + ) cosh (k (z0 + h)) . ω cosh (kh) (28)  stands for an eventual phase shift with a reference. In offshore sailing conditions, the water depth h is far bigger than the wavelength λ: kh  1. The potential hence yields: φ (x0 , y0 , t) = ζ0 g −kz0 e cos (kx0 − ωt + ) . ω 21 (29) Spatial and temporal periodicities are linked by the gravity waves dispersion relation: ω2 . g k= (30) And finally using respectively the velocity potential definition, the dynamic free surface condition and Bernoulli’s equation, one can derive

the waves velocity, elevation and pressure field: 3 3.1 u = ∇φ, (31) ζ (x0 , y0 , t) = ζ0 cos (kx0 − ωt + ) , (32) p (x0 , y0 , z0 , t) = −ρgz0 − ρgζ0 e−kz0 sin (kx0 − ωt + ) , = pH + pD . (33) (34) Steady and unsteady waves Ship in waves When a motionless body is put in such a wave system, it perturbs the previously computed potential as a new boundary condition is introduced along the body’s surface. This is usually dealt with by defining a diffraction potential φD which cancels out the normal component of the incident waves potential φI on the body’s surface: ∇ (φI + φD ) = 0 on Sbody . (35) The pressure field generated by the diffracted waves generates a new load on the body called diffraction force. Furthermore, due to the dynamic pressure pD from the incident waves, another load, called Froude-Krylov force is exerted on the body. It is expressed as: ZZ FFK = − pD ndS. (36) Sbody Let’s now consider a moving ship, µ will designate

the angle between the ship velocity vector and the wave main direction of propagation – 0◦ corresponding to following seas – and α the angle between the earth-fixed first coordinate and the waves, as shown in Figure 15. The relative position in the wave coordinates is: X = x0 cos(α) + y0 sin(α). (37) Figure 15: Wave angles definition. Due to its motion, the ship does not perceive the waves at the same frequency they have in the earth-fixed system. An encounter period Te is therefore defined to account for the period of the received waves: Te = λ . cw − Vb cos(µ) (38) Vb cos(µ) stands for the ship’s speed along the wave propagation axis and cw is the wave velocity, hence their difference corresponds to the waves relative velocity. Similarly, the corresponding circular frequency, the frequency of encounter is expressed as: ωe = 2π Te =ω− (39) ω2 Vb cos(µ) g 22 (40) When the ship is advancing in waves, the different waves systems add and interact. In a

first approach, it is often assumed that the potentials can be added to express the resulting potential: φ = φR + φI + φD . (41) As said, radiated waves are produced by the ship’s oscillations and therefore have the frequency of encounter ωe when the ships is oscillating in waves. On the contrary, both incident and diffracted waves keep the initial wave frequency, ω. It is to be noticed that the ship being now in motion, the boundary condition along her surface becomes: ∇ (φI + φD ) · n = Vb · n on Sship . (42) However, in reality, the different wave systems do not simply add to each other but interact, making an accurate computation of the global flow more involving. The main cause of interaction is the non-linearity of the moving body and free surface boundary conditions. The combined problem may therefore not coincide with the superposition of steady and unsteady solutions. The non-linear interaction terms may have the same order as the linear ones and - in those

situations - cannot be neglected. In this way, many models and theories restrict themselves to cases where a linearisation of the free surface conditions remain valid for both problems. Such restrictions are usually slender ship, low frequencies and slow speed or non-piercing deeply submerged bodies. 3.2 Wave-form parameter Nevertheless, the situation gets more involving as the forward speed increases. The steady flow of the translating ship becomes substantial and interacts with both the incident waves and the unsteady radiated waves generated by the oscillations. When studying the hydrodynamic coefficients of a moving ship, the wave-form parameter γ carries a significant information on the environment in which the ship proceeds. It is defined as: γ= ωVb . g (43) The importance of this parameter can be shown by considering ship radiated waves. Produced at the same frequency the ship oscillates, ωe , they verify the same equations as the incident waves (see Section II.22) and

can be written in the form: φR = Ae−kz±ikx−iωe t = φr e −iωe t (44) (45) Applying kinematic and dynamic boundary conditions on the global fluid velocity potential Φ = Vb x + φR , and combining them together while keeping only first order terms, yields:   2 ∂ 2 φR ∂φR 1 ∂p ∂ 2 φR ∂p 2 ∂ φR + 2Vb + Vb +g =− + Vb , on the free surface. (46) ∂t2 ∂t∂x ∂x2 ∂z ρ ∂t ∂x Assuming the air pressure remains constant over time and space at p0 , we have p = p0 on the free surface and the right-hand side of (46) is zero. With the given wave potential, it yields: " # 2 ∂ ∂ −iωe + Vb +g φr = 0. (47) ∂x ∂z This leads, through the definition of φr , to an equation in k. For waves propagating in the same direction as the ship is advancing (−ikx term in (45)), solutions only exist for γ < 1/4, they are:  p g  1 − 2γ ± 1 − 4γ , (48) k1,2 = 2 2Vb while for the waves in the opposite direction, there is no condition on τ ,

and the admissible wave numbers are:  p g  k3,4 = 1 + 2γ ± 1 + 4γ , (49) 2 2Vb 23 As said, k1 and k2 waves propagate in the ship’s forward direction. However, they move at a speed given by the group velocity (energy propagation): r 1 g i cg = (50) 2 ki It turns out that k1 waves are slower than the ship and remain therefore downstream, but k2 waves are faster and propagate upstream. Therefore, as long as γ > 1/4, the waves (k3,4 ) radiated by an oscillating ship moving forward at Vb are confined behind: the ship advances in still water (see Figure 16, right). On the contrary, when γ < 1/4, some waves (k2 ) are faster and the ship perturbs the environment in which she evolves (id., left) Finally for γ = 1/4, the wave front of radiations, perpendicular to the course moves at the same speed as the ship. Figure 16: Wave pattern of an oscillating source with forward speed. Reproduced from Yeung and Kim [1985]. Those differences in the encountered free surface state

lead unsurprisingly to variations of the hydrodynamic coefficients. It has indeed been observed (see Newman [1961]) that both added mass and damping coefficients present irregularities at γ = 1/4. Newman also showed in his study that at this value of γ, surge, heave and pitch damping coefficients go through a mathematical irregularity and become infinite while out-of-plane coefficients remain bounded. 3.3 High forward speed The ships which are studied here have particularly high speed. Also a suitable theory must be used when deriving the coefficients and loads in order to account for the substantial wave systems generated by the moving ships - steady and unsteady. As far as unsteady waves are concerned, it has been shown that divergent ones dominate over transverse ones at high Froude (see Yeung and Kim [1985]). The interactions between steady and unsteady wave systems are complex. When deriving unsteady wave patterns, they may have significant effects on the free surface and

body boundary conditions. One can compute, for a given speed, the minimal circular frequency which verifies γ > 1/4. With Vb = 30 , this limit frequency is ωlim = 0.159 rads−1 If the ship is oscillating at her natural frequency, this gives a maximal limit to it, corresponding to T0 = 40 s. 24 But when the ship is in waves, this conditions the encountered frequency ωe which is the frequency of the ship’s oscillations in waves: Vb cos µ g −ω 2 +ω > . (51) g 4Vb The absolute value is used as for high speed with small waves angles (following to beam seas), the encountered frequency may be negative as the ship is faster than the waves. The condition (51) is therefore equivalent to:  g  2 Vb cos µ  +ω− > 0 and ωe > 0 −ω g 4Vb (52) V cos µ g   +ω 2 b −ω− > 0 and ωe < 0 g 4Vb In either case, the left-hand part is more restrictive than the right-hand one. Details are given for the first case, the second case is dealt with

similarly. For following to beam seas, 0◦ ≤ µ ≤ 90◦ , the sign of the second order coefficient is negative and the relation (52) is true only between the roots of the polynomial. For following seas, the wave form parameter is always less or equal to 1/4. The admissible domain γ > 1/4 then grows with µ For beam to head seas, 90◦ ≤ µ ≤ 180◦ , the second order coefficient is positive and (52) is verified outside of the roots of the polynomial. As in addition the circular frequency must be positive, the only condition √ which holds is ω > g(1 − 1 − cos µ)/2Vb cos µ. For those wave angles, the lower limit for ω goes from 015 to 0.13 rads−1 at 30 kn which corresponds to maximal periods 40 s and 48 s In the case of our study, such long waves are not of interest. As said, negative frequencies of encounter are possible only in following to beam seas. For µ ∈ [90; 180], the roots of (52) are negative which is unphysical for wave circular frequency (see

Figure 17). Figure 17: Domain of validity of (52) (green area). Figure 18: Allowable ω for γ > 1/4 at 30 kn. Finally, as one can see in Figure 18, the frequency domain for which the wave form parameter is below 1/4 is rather limited, especially in head to beam seas when sailing at high speed. In those waves, only the longest periods (larger than approximatively 40 to 50 s at 30 kn) do not enable a sufficiently large value of the wave form parameter. For beam to following seas, the allowable frequencies occupy two different domains, one corresponding to positive ωe , and the other for negative ωe The size of the separation between the two domains decreases as the boat increases. In this way according to this model and for sufficient high speed, γ > 1/4 is verified on a large frequency domain, especially from 90 to 180◦ wave angles. Nonetheless, the model considered only the steady potential U x with the radiated waves. In waves, the existence of interaction would

probably change the previously derived relations. 25 4 Appendages 4.1 Model Appendages loads are computed through a model that uses as inputs: ˆ the rotation angles transforming the appendage coordinates to the fluid ones, ˆ the fluid velocity at a reference point for the appendage, ˆ the emerged length of the appendage, ˆ the extension of the appendage (which percentage of the appendage emerges from the hull), ˆ eventually, particular settings specific to some appendages, such as board trimmer, or rudder flaps. The output result is corrected so that if needed a piercing drag is added or, on the contrary, if the appendage leaves the water, the loads are cancelled out. 4.2 Appendages dynamics Considered as lifting surfaces, the appendages model is based on dynamical behaviour. Also such a model remains valid in unsteady situations as long as the inputs are based on effective values – that is to say computed and updated before each call to the model. Therefore, from

the attitudes and positions of the ship platform at the current time and the position of the appendage in the ship coordinate, the height of immersion is computed. Knowing the velocity vector of the fluid in the appendage coordinates enables also the computation of the rotation angles needed for the model. Knowing two vectors it is possible to define a rotation (axis u and angle α) and a scalar κ, so that the first vector rotated and multiplied by the scalar is equal to the second one. In the studied case, the initial vector is the xe vector (||xe || = 1) of the appendage coordinate system while the second one is the fluid velocity vector in the appendage coordinates VApp/RApp . The characteristic values of the transformation are: κ = ||VApp/RApp ||, (53) u = xe ∧ V̄App/RApp . (54) If ||ū|| = 0, then there is no rotation, α = 0. In all other cases, α is defined as:   xe · VApp/RApp α = arccos . κ (55) It is then possible to write the corresponding rotation matrix,

first using the axis/angle definition (56) and second with the coordinate transformation one (57):  u2x (1 − c) + c  R = ux uy (1 − c) + uz s ux uz (1 − c) − uy s  ux uy (1 − c) − uz s u2y (1 − c) + c uy uz (1 − c) + ux s cos θ cos ψ  = sin φ sin θ cos ψ − cos φ sin ψ cos φ sin θ cos ψ + sin φ sin ψ  ux uz (1 − c) + uy s  uy uz (1 − c) − ux s u2z (1 − c) + c cos θ sin ψ sin φ sin θ sin ψ + cos φ cos ψ cos φ sin θ sin ψ − sin φ cos ψ  − sin θ sin φ cos θ  cos φ cos θ (56) (57) (ux , uy , uz ) are the ū vector coordinates while c = cos α and s = sin α. A simple identification between the two definitions enables the computation of the effective angles of attack (ψ, θ, φ) of the fluid on the appendage. 26 4.3 Added mass and damping As explained in Schmitke [1977] and Schmitke [1978], added mass of appendages are negligible, especially because those terms are, unlike damping

ones, independent of speed. As an order of magnitude, the added mass of an appendage of mean chord c and span s in the direction perpendicular to its surface is: Aapp = π ρsc2 4 (58) To get an idea of the distribution of the hydrodynamics coefficients, one can consider an horizontal 2D foil oscillating harmonically in heave. Theodorsen [1935] has produced a complete solution expressing the foil lift as a function of its heave η3 and pitch motion η5 :  c  ρ (59) L = − πc2 (η̈3 − U η̇5 ) − ρπU cC(kf ) η̇3 − U η5 − η̇5 . 4 4 U is the foil velocity while C(kf ) is the Theodorsen function. Numerical calculations from this expressions shows that at 30 kn, damping and restoring coefficients are respectively 3 and 2 orders of magnitude larger than added mass ones. Similar results are obtained via Schmitke’s expressions, which incorporate finite span and free-surface correction factors. Added mass is therefore neglected for the appendages. Damping is accounted

for via the effective angles of attack and the magnitude of the incoming fluid velocity vector. 4.4 Dynamics in waves In waves, a method similar to the one presented in Section II.42 is used, including the value of local free surface elevation and wave velocity. Using the derived wave model, the local free surface elevation is known at the appendage location and its emerged length can thus be computed: Happ = zapp/R0 − zF S/R0 at app = zapp/R0 − ζ0 cos (kXapp − ωe t) . (60) (61) The effective velocity is also computed at the appendage through the appendage reference point velocity and the wave velocity at this point: (62) Vapp/Rf ,R0 = Vapp/R0 ,R0 –VRf /R0 ,R0 The first term is computed from the position of the appendage reference point in the boat coordinates and the boat translational and rotational speeds, the second term is given by equation (31). It is then possible to compute the effective angles of attack of the appendage with respect to the fluid using the method

given in Section II.42 A linearised expression of the loads on an L-foil heaving in head waves is derived in Appendix A. 4.5 Angles of attack in waves Although the effect of waves on appendages immersion is quite explicit, the changes of angles of attack they imply is less obvious. One can however point out that the larger the boat speed with respect to the wave induced velocity, the smaller the angles of attack variations. To get an idea of the order of magnitude of the effect, a simple model is proposed: the ship is assumed to be uniformly translating at V¯b , with neither angular nor vertical velocities. The implantations angles of the considered appendage are denoted Yawi , Rakei and Canti , while for conciseness the notation ES = −e−kzapp sin(kX − ωe t) is used. Considering the waves and ship speed component along the horizontal plane, the effective yaw can be expressed as:   ζ0 ωe ES sin(µ + η) Yaweff = Yawi + η + arctan . (63) Vb − ζ0 ωe ES cos(µ + η)

Assuming sufficiently small waves, one can write: δ = ζ0 ωe /Vb  1, which in (63) yields: Yaweff ≈ Yawi + η + δ × ES sin(µ + η). 27 (64) In the horizontal plane, the fluid velocity is then: Veff = Vb (1 − δ × ES cos(µ + η)). Let’s now consider the vertical plane rotated of the effective yaw angle. A similar development is possible, and one can write:   ζ0 ωe ES Rakeeff = Rakei + η5 + arctan , (65) Vb (1 − δ × ES cos(µ + η)) ≈ Rakei + η5 + δ × ES. (66) A similar expression may be derived for the effective cant angle. In this way, for a ship speed of Vb = 30 kn in head waves of period T = 10 s and amplitude ζ0 = 1 m, we have: Te = 5.03 s, δ = 008 and ∆Yaw = Yaweff − (Yawi + η) = ± 46◦ if leeway is neglected This model is highly simplified, especially in the way it neglects the yacht motions that those waves would imply. In particular, for the largest δ, that is high amplitudes and frequencies, the actual waves slope αw are substantial

and the consecutive ship motions are not negligible: αw ≈ ζ0 ω 2 ζ0 = . λ 2πg (67) However it shows that in waves, the angles of attack variations can be relatively important. One can fear that such a dynamic behaviour may cause the appendage to stall in the seaway. This is an important issue as in this case the lift force would suddenly drop, which in fly mode can cause important loss of speed and ride height, and possibly structural damages if the yacht hits the water when loosing height. 4.6 Wave loads As detailed previously, incident waves cause two types of loads: Froude-Krylov and diffraction. The first one is expressed according to: ZZ FFK = − pD ndS. (68) Sapp pD is the dynamic pressure due to the incident waves. In the general case, the wavelength is considerably bigger (10 − 100 m) than the appendages dimensions (≈ 1 m). The variations of pressure on the different surfaces of the appendage are therefore rather low. This allows to neglect the Froude-Krylov

force on an appendage (see Wave-induced motions in foilborne conditions in Faltinsen [2005]). As far as as diffraction loads are concerned, deriving a proper expression is more complex. In Faltinsen [2005] (see id.), a simple model is derived for wavelength assumed much longer than the foil mean chord The principle is to apply condition (42) on a particular point of the appendage. As it is supposed much smaller than the wavelength, the variation of potential along it are assumed small In the work presented by Schmitke ([1977],[1978]), only the diffraction loads on the hulls are accounted for. 4.7 Other dynamical aspects Numerous studies in the aerodynamic field (on flapping airfoil for instance) have demonstrate that lifting surfaces dynamical response do not only lie on added mass and damping coefficients but that other effects are also at stake. For the time being, they are neglected in the dynamic simulator Dynamic stall is probably the major ones. As shown in Section II45, waves

can lead to substantial changes in angles of attack. Dynamic stall happens at larger angles of attack than static stall and is due to a separation of the incoming flow leading to vortex formation Substantial variations in the hydrodynamic loads are encountered. Stall recovery – flow reattachment – happens on the contrary at an angle lower than in static This creates a hysteresis phenomenon. It has also been shown that the work done over a full oscillation cycle produces a damping effect, which may be negative (and cause flutter for instance). 28 5 Hulls 5.1 Model Coming from CFD runs, the response surfaces of the hulls use as inputs: ˆ the norm of the fluid velocity, ˆ the leeway angle, ˆ the current trim and sinkage. The provided speed accounts for a forward motion and not for any vertical part. 5.2 Restoring force The response surfaces outputs F RS therefore mainly account for the restoring forces. For special needs, it is possible to explicitly derive the Cij

terms around a certain point: Cij = dFiRS . dηj (69) Figure 19 shows the lift force function of heave for the central hull of one of the considered ship, at zero leeway and 30 kn of speed. As expected from the non verticality of the ship’s side, the response is not linear – that is to say the restoring coefficient C33 is not constant. Using a constant value is therefore only valid for small motions. Figure 19: Hull lift force depending on heave. As for a conventional boat symmetric with respect to the centre-plane, when there is no leeway, there is no coupling between in-plane and out-of-plane degrees of freedom. This is no longer true when leeway appears The response surface gives the following orders for the restoring coefficients at 30 kn of boat speed, with respectively 0 and 1◦ of leeway:     0 0 104 0 −103 0 0 0 104 0 −103 −102 0 0 0 0 −103 0 103 0 0 0 103  103          5 4 5 4 ◦ ◦    0 0 −10 0 10 0

0 0 −10 0 10 10  C0 =  C1 =  , . 0 0 0 0 0 0 0 0  0 0 0 0      0 0 106 0 −105 0 0 106 0 −105 −102  0  0 0 0 0 0 105 0 0 −105 0 104 105 Those numbers are only given to get an idea of the orders of magnitude. No restoring coefficient but the direct outputs from the response surfaces are used for the simulations. 29 5.3 Wave loads The Froude-Krylov force is computed via added volumes and only taken into account for heave and pitch. Each hull is separated in a number of longitudinal slices, which are themselves cut in vertical parts for which the vertical positions (z1 , z2 , ., zN ) and the mean horizontal areas (S1 , S2 , ., SN ) are known For each longitudinal slice, the free surface elevation is computed using the previously detailed wave model. It is compared to the positions of the boundaries of the vertical slices and to the rest-level of the free surface. For all vertical slices between the rest-level

and the wave elevation, a volume δV = ∆z × Si is added, corresponding to a lift force ρg∆zSi (∆z may be negative). An example is given in Figure 20. For the pitching moment, a lever arm corresponding to the longitudinal distance from the slice to the centre of gravity is added. Figure 20: Added volume method. No diffraction force is taken into account. This is valid for short incident waves, in which the diffraction force is negligible with regard to the Froude-Krylov one (see Gerritsma [1991]). However, for longer waves, this assumption leads to growing errors. Therefore developing a model to account for such effect is needed to be able to operate simulations on a wide range of wavelengths, see Section II.6 5.4 Added mass and damping Added mass and damping are a key aspect of hulls dynamic response. However they are involving to determine According to Faltinsen [2005], added mass in surge is roughly 5% of the ship’s mass, which can be considered negligible. Out-of-plane

coefficients are not critical as the appendages provide a consequent damping due to their main orientation. Rudders, keel blade, and foil shafts are usually almost vertical and therefore strongly counteract sway, roll and yaw perturbations. The sail plan also provides roll damping In this way, it is especially regarding heave and pitch that the determination of hulls added mass and damping is crucial. Next section offers a review and discussion of the existing methods and theories to derive those coefficients. 30 6 Hulls hydrodynamic coefficients As we have shown in the previous section, hulls dynamic behaviour is so far not addressed satisfactorily. The high speed effects make the development of a model even more complicated. In this section, a review of the state of the art on this matter is made and the methods are discussed in terms of adequacy and feasibility. Most of the presented methods enables in the same process the computation of the diffraction loads. 6.1 6.11

Numerical methods Seakeeping prediction Numerical prediction of ship seakeeping performances is – as in many cases – a compromise between computational speed and accuracy. Typical Reynolds numbers for ship hydrodynamics problems are 107 - 109 , meaning fully turbulent flow on a substantial part of the boat. DNS (Direct Numerical Simulations) are unusable in such flow cases for lack of computational power. Two of the main categories of methods are potential flow solvers and turbulent solvers such as Euler, RANS (Reynolds-Averaged Navier-Stokes), LES (Large Eddy Simulation) and DES (Detached Eddy Simulation). As they only need to solve one linear differential equation instead of four non-linear coupled ones, potential flow solvers are substantially faster. Moreover, potential flow methods are mainly using Boundary Elements Methods (BEM) which enables the discretisation only the boundaries (2 or 3D) of the domain and not the whole fluid. This provides a substantial advantage in

computational speed as grid generation is a computationally heavy operation. Nonetheless, as they need continuous and simple free surfaces potential flow methods are usually irrelevant for breaking waves and slamming for instance. The recent years have seen and increasing development and use of turbulent models. RANS solvers remain the common reference and are now widely used in hydrodynamics problems (see especially Roux et al. [2008] and Böhm [2014]). In this approach, the driving principle is to decompose velocity and pressure in two values, an average part and a fluctuation part. While the latter is used for turbulence modelling, only the average parts are used to solve Navier-Stokes equations. While enabling convincing accuracy, RANS solvers demand a turnaround time and a computational power which are currently too high to generate wide response surfaces for the simulator. In addition, it would then maybe be more profitable to directly implement a full RANS simulator. Also in a

first approach, attention was drawn to potential methods, which for an accuracy that is thought almost sufficient, allows much faster computations. The main existing potential flow methods can be sorted as follow: Strip theory: Initiated in the late 50s, today’s theory is mainly based on the work of Salvesen, Tuck and Faltinsen (see Salvesen et al. [1970]), and referred to as STF strip method The main principle is to use the ship slenderness to reduce the three dimensional problems to a number of two dimensional problems expressed in strip x = constant. The steady flow φs is neglected while the free surface condition is simplified. It also includes corrections to account for (low) speed effects. The two dimensional problem is solved either analytically or by panel method. Analytical solving is carried throughout multi-parameters conformal mapping transforming circular sections to ship sections (Lewis sections). This approach prevent the study of peculiar geometries such as submerged

bulbous bow. Panel – or close-fit – methods are presented later in Section 6.12 Strip methods carry a lot of inherent assumption, but have however proven extremely efficient (speed and price) for many conventional low Froude ships. Such methods are often completed with empirical corrections Frequency domain The consistency of conventional strip method for Froude number larger than 0.4 is often discussed Comparison with experimental data have shown in several cases good agreement up to F n ≈ 0.6, when dynamical sinkage and trim as well as steady wave profile along strips are taken into account However, as computed earlier, the studied ships have Froude with order of magnitude around 1 and conventional strip theory is thus irrelevant. 31 Unified theory To account for lower frequencies waves, the unified theory of Newman [1978] couples a flow in the near field (inner solution) found from a two-dimensional approach due to the body slenderness, and a three dimensional far field

flow (outer solution) derived from a point sources distribution along the ship centreline. The distance of separating the ship from the overlapping region is chosen so that the error is minimised. 2D1⁄2 theory Also refers to as HSST for high speed strip theory, the methods is still based on the consideration of a number of two dimensional ship strips. To account for high speed effects, the strips are no longer considered separately but with respect to the upstream one, with initial conditions at the bow. For each strip, the boundary conditions helps solving wave elevation and velocity potential while longitudinal derivatives are evaluated as numerical differences with the upstream strip. Such methods are still able to use the fully non-linear free surface equations The initial condition at the bow is taken as a zero disturbance and is therefore not valid for small wave form parameters γ. This theory neglects transverse waves but steady and unsteady diverging waves are well

represented, leading to an improved accuracy in the computed hydrodynamic coefficients Yeung and Kim [1981] introduced a pseudo-time variable t∗ = −(x − L/2)/Vb . With such a notation, ship strips are periodic in real-time (when t∗ is kept constant). The hull boundary condition is linearised about its mean position for unsteady solution computation. Compared to conventional strip theory, Yeung and Kim introduced quasi-three-dimensional effects avoiding the previous high speed restriction. In their paper, Faltinsen and Zhao [1991] have developed a different high speed theory in which the interaction between the local steady flow and the unsteady waves is accounted for. The steady flow is derived first which enables then the computation of the unsteady flow with a free surface condition linearised about the steady solution. Non-linear boundary conditions are possible when calculating the steady solution As previously, the zero disturbance initial condition at the bow implies

sufficient high speed to be reasonable (γ > 1/4). Comprehensive theory: Similarly to the unified theory of Newman, Yeung and Kim [1985] developed a prolongation of their theory for lower frequencies. The outer field of Newman [1978] which takes into account three dimensionality and forward speed can directly be added to Yeung and Kim previous inner solution. The initial condition is however not obvious for small γ and a pseudo-time reverse flow function is introduced. Forward and reverse solutions are summed and modulated to give rise to the full inner solution. 3D methods For non slender body, the previous approach are inconsistent and full 3D methods are to be used (3D panels for instance). Due to their computational cost and the considered ship’s slenderness, this report will not go any deeper into 3D approaches. Conditions Low Fn High Fn γ> 1/4 STF Salvesen et al. [1970] HSST Yeung and Kim [1981] Faltinsen and Zhao [1991] Any γ Unified theory Newman [1978]

Comprehensive Yeung and Kim [1985] Figure 21: Seakeeping numerical methods for slender hulls. 32 6.12 Panel methods The main principle of panel methods, also called boundary element methods, is to discretise the boundaries of the flow domain in a finite number of elements generating simple flow with an unknown strength. Enforcing the boundary conditions on the panels generates at set of coupled linear equation leading to the determination of the flows strength. Laplace equation being linear, all potentials are summed up to express the global flow potential. Such a method can be used in two or three dimensions, in which cases boundary elements are respectively segments or panels Hereafter, the term panel refers without distinction to the relevant one Usually, the solid body boundary condition is enforced on the hull at only one collocation point per panel - its centre in most cases. However other methods exist Three main boundary elements methods are used and will be briefly

presented here: Green-function, Rankine singularity and combined methods. Green-function methods (GSM): Panels are distributed along the averaged wetted surface of the hull. The velocity potential of each panel being a Green-function, it fulfils Laplace equation, the free surface and radiation conditions (without consideration of the steady potential). The unknown source strength is derived from the coupled equations given by the non-penetration condition on each panel. The time history of the flow appears in involving integrands named convolution terms. GSM cannot account for non-linear boundary conditions Wehausen and Laitone [1960] provides a list of elementary Green functions for diverse sources behaviour (varying strength, steady translation, pulsating, etc.) Rankine source methods (RSM): Rankine singularity methods allow to completely account for the steady potential as well as more complicated boundary conditions. However, unlike the previous method, the free surface has also

to be panelled when using Rankine sources, thus a higher computational cost. Source strengths are computed from boundary conditions (hull and free surface). Special care needs to be taken regarding numerical reflection on the borders of the numerical domain, generally through the use of numerical beaches. Combined methods: As one can guess combined GFM-RSM methods aim at avoiding the weaknesses of the previous methods by gathering the advantages of both. Different possibilities exist Initially the idea was to match near-field RSM potentials to far-field GSM solutions. To increase speed or computational stability, some techniques have demonstrated a substantial efficiency. Desingularisation has in this way shown great capability in avoiding the singularity of the potential and its derivatives on the panel. The principle is to set the source point not on the panels but at some distance, outside the fluid domain. Desingularisation is often used for free surface, and sometimes for body

boundaries - demanding special care for the narrow parts. Desingularised source must neither be too close to each others nor too far from the collocation point they depend on (usually between 0.5 to 2 grid spacing) Also laying on desingularised source, the patch method confers a highly improved accuracy without any increase in the grid generation time. Instead of being enforced at only one collocation point per panel, the non-penetration boundary condition is computed as an integral along the whole panel. Potential and velocities are calculated at the path corners while pressure uses averaged velocity as well as corners velocities. However, patch method makes the computation of potential higher derivatives more complex than for regular panel method. 6.2 Discussion Numerical methods have been presented and their major drawbacks and advantages explained. In order to be able to select a suitable method in regard of the needs and studied ships, the main data are recalled and discussed

in the present section. Examples of studied ships have been presented in Section I.4 They are high speed ships with Froude number higher than 0.9 From the data provided in Table 1, one can see that with low beam to length and draft to length ratios, the hulls can be considered slender. 33 Nonetheless, for multihulls, the interactions between hulls may play an important role, especially at low speed. Discussion of such interactions has been presented in Section I.34 Due to the high ship speed, one can expect them to be rather small. But if two dimensional methods were to be used, corrections to account for the possible reflections of diffracted and radiated waves between hulls may be needed. One of the major constraints on the choice of a numerical method is the computational power it demands. Deriving hydrodynamics coefficients for high speed boats using three dimensional methods is highly expensive in terms of computational power and time. As long as this is an issue, a 2D1/2

methods seems the most efficient way to derive hydrodynamic coefficients with enough accuracy Faltinsen and Zhao’s theory brings the possibility to account for steady interactions, while the condition γ > 1/4 is not too restrictive at high speed as discussed in Section II.33 As far as viscosity is concerned, experimental or empirical corrections should be brought especially regarding roll as it is known that roll damping is usually mainly driven by viscosity. Viscous effect corrections are also highly necessary in short waves, as orbital velocities are substantial and potential flow forces are small (see Vugts [1970]). Considering the fast evolution of numerical techniques and equipment, it is highly conceivable that within a short period, RANS methods become a feasible solution. The implementation of any of those numerical methods demands a thorough work of development and validation that falls beyond the scope of this master thesis. However, a list of possible solutions was

presented and discussed, allowing for a better understanding of the elements at stake in the choice of a seakeeping numerical prediction method. 6.3 Determination from CFD To get an idea of the orders of magnitude of the studied ships hydrodynamic coefficients at different pulsations, Computational Fluid Dynamics (CFD) is thought to be able to provide some help. The aim would be to examine heave and pitch coefficients and the related coupling terms For a given frequency of oscillations, the following process is proposed: the simulation is ran in 1 degree of freedom with a fixed forward speed, and the body exited in heave or pitch at the chosen frequency. From the pressure distribution on the hulls surface, it is possible to determine the heave and pitch loads. In 1 degree of freedom, with excitement ηj = ηi0 ejωt , we can rewrite the equations of motion (see 15): j ∈ {3, 5} i, j ∈ {3, 5}, i 6= j (Mjj + Ajj )η̈j + Bjj η̇j + Cjj ηj = Fj , (Mij + Aij )η̈j + Bij η̇j +

Cij ηj = Fi , (Direct) (Coupling) (70) (71) This yields, with Fi = Fi0 ej(ωt+δ) :   Fi0 ejδ = −ω 2 (Mij + Aij ) + jωBij + Cij ηj0 . (72) As said previously, the restoring term Cij is known from the steady motion fittings, Fi0 and δ are measured from the CFD runs and for each j ∈ {3, 5} we have two sets (direct and coupling) of two coupled equations (arguments and modulus) with two unknowns: Aij and Bij . This method is computationally very expensive – runs must last several excitation periods to have reliable results, runs must be made for each frequency, each degree of freedom – but provides the advantage of only demanding technologies that are already used in the company (unlike 2D1/2 codes). Therefore it enables the calculation of the hydrodynamic coefficients in a pre-study with low engineering cost. Obviously hull displacement and speed also affects the hydrodynamic coefficients. If possible different values should therefore also be tried. As far as

excitation frequencies is concerned, it is expected that the behaviour at low and high values differs substantially. Extremal numbers need thus to tested,as well as intermediary ones. p Usually resonance happens around ω L/g ≈ 4 − 5, using values of this order seems thus reasonable. Such runs have however not yet been ran due to a busy planning of the computational tools of the company. 34 7 Aerodynamic response The modelling of the aerodynamic loads in unsteady situation is currently a central issue in sailing yacht research that draws numerous scientists and engineers (see for instance Fossati and Muggiasca [2013]). The huge complexity variations between laminar flows in upwind conditions and turbulent in downwind conditions have been and is still largely studied. In addition to the dynamic phenomena (added mass, dynamic stall) that have previously been described for lifting surfaces, the sail plan presents the major difference of being highly deformable and to be in

reality non rigidly fixed to the rest of the boat. In this way, under the ship unsteady motions, the sails can flap, while, due to the unsteady wind loads, its shape changes along time. Aerodynamic is therefore a field of huge interest for fluid structure interactions (FSI) studies. Last, it has been shown that the sail plan brings an important contribution in added mass, especially around the longitudinal axis (see Graf et al. [2007]) And in a first approach, a simple added mass coefficient could be taken into account: π (73) Asail = − ρair cAsail zCE . 4 c being the sail main chord, Asail its area and zCE the vertical coordinate of the centre of efforts. However, for the time being, no added mass coefficient has been implemented. 35 Part III Stability of offshore sailing foilers 1 Foil and stability, state of the art 1.1 Foil principle Since Monitor in 1957 (see Figure 23), the use of foils for boats in general and sailing yachts in particular has soared. From V-shaped,

they have now browsed a substantial part of the alphabet: T, L (Fig. 22), C, J, etc. The main principle is that at a speed high enough, the vertical force those lifting surfaces generate allows to decrease the hull displacement-to-length ratio and its drag. The study of the lift to drag ratios of a hull and a foil (see Figure 24) is a good illustration of foils domain of interest. For low speeds, foils only add drag, while the hull(s) generate the whole lift. As speed increases, foils produce an increasing amount of lift which induces an additional drag cost. This initial disadvantage is broken when the produced lift is sufficient to raise the hull at a draught at which the total hull and foil drag is lower than the drag the hull alone would produce. From this speed on, the advantage grows, finally reaching a full fly mode at which there is no more hull volume in the water. This takeoff speed is part of a whole strategical choice when designing the foil. For a given ship displacement,

a low takeoff speed demands a larger foil lifting surface than a higher takeoff speed, but this has a cost, in drag in particular. Cavitation is also a driving issue. Figure 23: Monitor, in 1957. Figure 24: Indicative lift to drag ratios of a hull of sailing multihull and a L foil. Figure 22: Foil terminology. From cupinfocom 36 1.2 Conditions of use and class rules In sailing, asymmetrical ship configurations generally leads to better performance than symmetrical ones. Liquid ballast and canting keels are good examples This principle applies also when foils are set to the boat: with a symmetrical configuration (both foil generating lift), the windward foil lift generates a heeling moment that adds with the sails’ one, decreasing the ship’s power (her righting moment). But a one side configuration (leeward lift only) requires switching the foil configuration at each tack. Hence the choice of symmetrical/asymmetrical foil use must be weighed against the needed tack

frequency but also against the crew size. America’s cup five-man crew – helped with hydraulic oil – has not the same constraints than a single-handed sailor on an A-class catamaran (see Figure 25). However other solutions exist, as for instance, differential rake: unloading the windward foil to decrease heeling moment. (a) A-class catamaran Z foil. a-catorg (b) AC45 L foil. T Burlet / Groupama Team France (c) IMOCA 60 ”Dali” foil. Voile Banque Populaire (d) Kite foil. Alpine Foil Figure 25: Different foils for different needs. When designing a foil, numerous aspects are to be taken into account and precise specifications must be established. Is it a foil for foil-assisted sailing (partly easing hull drag) or for full flying mode? For inshore regatta or offshore racing, with waves and limited repairing possibilities? For a small solo yacht or a huge record-breaker trimaran? In additions, as sailing ships are often subjected to specific design rules, the foil may also have

to comply with class rules such as for instance a minimal distance from tip end to centreline in the catamaran A-class. This is a major aspect in foil design innovation, conditioning the developed concept while compelling engineers and architects to imagine new shapes that meet the current specifications. 37 1.3 The question of the stability On a conventional sailing boat, the sailor acts on the sails, his position on the boat and the rudder angle. Speed, leeway, trim, heel, sinkage and heading are consequences of those settings. In gusts upwind, the boat may suddenly heel and accelerate and in reaction the sailor can let out some sheet or luff a bit. The yacht then self-adjusts to those new settings. With foils, things are getting more complicated as sinkage and trim variations are can be much larger and swifter. Due to the strong quadratic link between speed and lift, small changes in speed – in a gust for example – can lead to substantial variations of the boat’s

attitude. An immediate counteraction when stability is lost is often needed, which is not always possible to execute for the crew, especially in single-handed sailing. This issue is often addressed through the use of control systems – mechanical or electronic. This is active stability (see Figure 26). Such solutions raise other questions, regarding robustness, durability and maintenance in particular. Besides, electronic systems demand an energy, which is an expensive resource aboard a sailing yacht, especially in offshore racing and has a cost in terms of weight or drag (hydro generators, fuel, etc). Finally class and race rules do not always allow such equipment on board. Figure 26: An example of mechanical control: the moth and its palper. Junichi Hirai / Mach 2 Moth This shows the interest of passive stability, which is the capacity of the system to go back to its previous equilibrium without corrective action. Passive stability systems can also be mixed with active control.

Finally, it is important to note that passive stability has generally a direct cost in terms of drag. Being able to design a system that presents sufficient stability characteristics while minimizing the additional resistance is therefore a central matter. Currently, foiling yacht stability is mainly studied regarding pitch and heave stability. The next two sections aim at presenting the tools and measures on which those studies are based 1.4 Pitch stability and airplane equivalence Pitch stability is reached – as for an airplane – with rear lifting surfaces (the stabilizer for planes, the T-rudders for foiling boats). The principle is that following a perturbation in pitch, the subsequent change in angle of attack on the rudder elevator provides a change in lift which with the lever-arm counteracts the perturbation moment. However, the main foil also sustains a change of angle of attack. Stability is therefore reached only if the total changes provide negative feedback, that

is: dM y < 0. (74) dη5 It can be shown that after a perturbation in trim, the changes due to the lift variations of both appendages are the main variations in My . Equation (74) therefore yields: d(−lf Fzf ) d(lr Fzr ) + < 0. dη5 dη5 (75) lf and lr are the respective lever arms of the lift force Fzf and Fzr . To account for the sails moment from the driving force, the lever arms are not directly the longitudinal distances xf or xr from the appendage to the centre of gravity, but distances corresponding to an ”effective” centre of gravity that accounts for the sails moment: Fs (76) li = xi + zCE x Mg 38 Still in a rough two dimensional approach, the lift forces are expressed as: Fzi = 1 ρSi V 2 CLi . 2 (77) Considering that CLf ≈ CLr , Equation (75) yields finally: lf Sf < lr Sr . (78) In the aeronautic field, the horizontal tail volume coefficient Vstab is therefore used to assess a plane pitch stability: lr Sr (79) Vstab = lf Sf Usually a threshold

value Vcrit > 1 is used and the rear surface is calculated as Sr = lf Sf Vcrit /lr . In the case where there are more than one rudder/foil, the previous equations remain valid by referring to the total available elevator area. Figure 27: In-plane equilibrium in fly mode. However, a major difference with airplanes is that foiling yachts are moving at the interface of two fluids, appendages can completely leave the water under strong perturbations. Moreover other forces (sail force in particular) have been neglected but they do have an impact on pitch equilibrium. The drag cost of stability is here obvious, as the larger the rudder elevator surface the larger the drag it generates. There is thus a great interest in trying to put the rudders as far as possible from the centre of gravity. But in the design process, this distance is usually conditioned by other constraints and not a variable parameter to ensure pitch stability. 1.5 Heave stability and foil design Heave stability is

a complex issue for which the optimum solution may not have been found yet. It is strongly linked to the foil geometry. The main principle is that as the boat rises above the water, the lift force must decrease for the boat to find an equilibrium position. Due to the height of the foils (a few tens of centimetres to a few meters), the stabilization must happen in a narrow range. However, it may be useful to have a foil that is initially unstable, so that it may rise quickly but then become stable so that the ride height stabilizes. 39 Two main stabilizing reactions are at stake: ˆ For lifting surfaces that pierce the water surface (V foil for instance), when the boat rises, the wetted lifting area decreases, the lift force decreases. ˆ When the piercing part is mainly vertical (L foil), the stability is provided by the leeway: as the boat rises, it loses vertical area in the water – the aerodynamic side force remaining constant – the leeway increases. This makes the angle of

attack on the tip (responsible for the major part of the lift force) decreases and therefore the lift force decreases. One could argue that a surface piercing tip also provides heave stability due to the reducing tip wetted area. However this is not a desirable situation as it implies an additional piercing drag, it is more beneficial to keep the tip fully immersed and to only lay on leeway-modulated stability. Since the 34th America’s cup, L-foils have been widely used for multi-hulls. For such a shape, two parameters are essentially at stake for heave stability: ˆ Cant: increasing the cant increases the strength in the Fy /Fz coupling. This can be seen through the dihedral angle (that is the ratio vertical/horizontal force: the larger it is, the greater is the loss of vertical force following an increase of ride height, ˆ Tip angle: angle between shaft and tip. The smaller it is the more stable the foil For a given increase in ride height and therefore in leeway due to the loss

of vertical surface, a tip with a smaller angle will encounter a larger loss of angle of attack and is therefore more stable. When the angle is less than 90◦ between shaft and tip, the side forces they generate are opposite, as the angle gets smaller, the negative influence of the tip side force increases. This is therefore a drawback for performance and a compromise must be sought after. Depending on the angle between wind and boat, the needs are not the same. Off the wind, the aerodynamic side force is lessened, and less vertical area is needed from the foil to maintain a correct course On the contrary, upwind, the side force is maximal and a substantial vertical surface is needed to avoid a too large leeway. Such variable needs can be met by a varying cant system but also by adapting the ride height to the course: fly high downwind, low upwind. This ride height may be reached by adjusting the rake angle Currently, to assess the heave stability properties of a foil geometry, a

common approach (see Fisher [2014]) is to draw the evolution of the derivative of the lift force with regard to ride height, remaining at the same value of side force (that is enabling the boat to drift more to compensate the loss of vertical wetted surface while the aerodynamic side force remains the same): dFz . (80) Istab ∝ − dz Fy With such a definition, the larger the index Istab , the larger the foil heave stability. Figure 28 shows such a stability diagram for a L-foil. One can see that as said the stability index increases with cant. Figure 28: Heave stability index. In multi-hull offshore racing, foil assisted sailing is the norm since the early 2000s (in the ORMA class for instance). Sailing in full fly mode around the world may be the next step Reliability, durability in waves and offshore conditions, and of course stability are some of the driving issues. And as always in high performance racing , the balance between reliability and performance is to be pushed forward,

as far as safety is guaranteed. Understanding and managing a good stability is a key aspect. One of the simulator’s main goals is to study stability and improve the static criteria that are currently used. 40 2 2.1 Dynamic considerations Initial interest Most solutions to reach heave or pitch stability involve drag penalties. For instance in the case of L-foils, as the tip produces drift force, the shaft needs to produce more anti-drift force leading to an increase in drag (induced drag). There is therefore a potential benefit in being able to safely decrease the stability goals by showing that static studies are too conservative. Dynamic studies of hydrofoils linear stability have already been carried out regarding full fly mode (see for instance Heppel [2015]). With the simulator, the aim is to be able to account for totally non-linear phenomena as for example boat’s elements that go out of the water (rudders elevators for instance) or on the contrary fall in the water

(the hulls). It may indeed be more profitable to sail in hybrid mode than in a complete flying mode. Besides, in some peculiar situations, non-linear response may be the only real behaviour Considering a perturbation in negative heeling is a good example: due to the boat dihedral and the heel, the windward foil and rudder are not expected to produce any lift/yaw moment and are thus retracted. Therefore only the sails – and sometimes the central T-rudder – provide heeling moment. Those moments do not vary a lot when the heel decreases, there is thus no – or little – counteraction when the ship is loosing heeling angle. In many cases this is the central hull that has to enter the water and provide the counter effect. This is a non-linear effect that the present study aims at capturing. Finally, such a study must be considered as part as an global design and optimisation process which has also other constraints. Some properties of the appendages are previously fixed in the

process such as for instance the longitudinal positions Stability is therefore to be reached only by altering some given parameters (geometry, cant, horizontal area, etc.) for which changes remain possible In the following studies, a ship similar to those presented in Section I.1 is entered in the simulator The boat is flying at a chosen equilibrium (speed, attitudes, ride height) which is adjusted via its aerodynamic properties so that all compared simulations start from a comparable situation. A perturbation in pitch or heave is then applied in the form of a Dirac function. The response of the ship to this perturbation is measured 2.2 Pitch study In this frame, it has been decided to study the link between dynamic stability and the static index for pitch stability while varying the rudders surface. This was done through the measure of the time of convergence – in the case where it exists. When the system is fully unstable and diverges, the convergence time is set to infinity, it

happened however only on the boundaries of the simulated values domain. Linear studies could have been carried out and a logarithmic decrement could have been approximated but it would not reflect the numerous non-linearities of the response. As can be seen from Figure 29, the global evolution of the convergence time is consistent with the expectations: it initially strongly decreases with the stability index, while being almost constant for intermediary values. For the highest values, the convergence time increases, the counteraction may be too large. One important observation is that the boat reaches good convergence already for Vstab ≈ 0.9 Sizing the rudders so that Vstab > 1 may therefore be far too conservative, and unnecessary to have sufficient stability. Figure 29: Time of convergence as a function of the stability index. 41 However a unique study is not enough and must be completed. In additions, it is obvious that other outputs are at stake in such a situation and

convergence time is not the only data that conditions the quality of the stability and the capacity of a system to recover from a perturbation. It is also of interest to work with stability index altered through other variables than the rudders elevators surfaces. In this way, the same simulations and measurements were carried out for a centre of gravity 2 m aft from the previous one. Results are presented in Figure 30. One can see that the results are similar, for a given stability index, the convergence time are relatively close. The minimal convergence time is the same in both situation. Having a centre of gravity closer to the bow turns out to be even better that what is shown by the stability index – which is already lower for a given rudder area when CG moves further aft. Figure 30: Stability study for different position of CG. The stability index, as described in Section III.14, is based on a two dimensional approach. It is therefore interesting to study how variations in

the third dimension affect the results. The stability index is computed using all the rear lifting surfaces, which in the current case are made of the leeward float rudder and the central rudder. For a given stability index – that is to say a given total rear lifting area –, it is however possible to change the relative part of central and leeward rudders. Figure 31 shows the result of a study in which three different area variations are studied: on both rudders, on the central one only and on the leeward one only. Figure 31: Stability study for different central and leeward rudders area. It turns out that the change in convergence time are rather small. Having a smaller central rudder than leeward seems to be weakly beneficial. The generated variations remain however small with regard to the variations due to the stability index in general. The two dimensional approach appears as satisfactory, even for three dimensional changes. In conclusion, the study has shown that the

stability index is a consistent measure of the system stability in pitch. Nonetheless, it seems to be possible to decrease the aimed reference index when designing rudder elevators. 42 2.3 Heave: one-parameter study A similar study was carried out for heave stability. The static stability index was compared to the effective dynamical response in the simulator. Different foil configurations were tested Figure 32 shows as previously the convergence time as a function of the static stability index. Each point corresponds to a given cant, 5 values spaced of 5◦ each were tested. As explained earlier on, increasing cant is expected to lead to increasing stability properties. This is shown by the evolution of the static index, which is negative for the smallest cant and then increases as the foil is canted in. A first and important observation is that the first point has converged. Although the convergence time is about 4 times larger than for the other values, this shows that even

with a foil considered as unstable, the boat simulation can converge. As the stability index increases, the convergence time decreases slowly. As previously there may therefore be an advantage in using a foil that is considered as having weak static stability properties but which in dynamical conditions proves almost as performing as more stable ones. Figure 32: Convergence time as a function of the stability index. Different foil set-up were simulated. In Figure 33, the results for different shaft lengths are shown. One can notice the almost complete similarity between the two results. This is an interesting result The imposed perturbation is a Dirac function which amplitude is fixed to ∆z = 20 cm. Therefore, as the foils are identical in all but shaft length, such a perturbation leads to the same loss in vertical area. As they do not have the same shaft area As , the side force distribution is not the same: Fy = Py(S1) × A(S1) = Py(S2) × A(S2) , s s (81) where S1 and S2 refer

to each of the foils. The loss of a given vertical area therefore do not lead to the same loss in side force and consequently to the same leeway increase. The equilibration phenomenon that occurs following the perturbation is therefore not exactly the same but leads nonetheless to highly similar convergence times. Figure 33: Stability study for different shaft lengths. One can however notice that for the minimum cant, the convergence time difference is substantial. Actually, as tconv = 300 s = tsim the simulation time, the shortest shaft length has not even reached an equilibrium yet at the end of the simulation. The substantial differences are probably due to the fact that for those low stability set-ups, the motion amplitudes and the angular velocities are important. This could lead to situations in which the shaft length variations are non-negligible consequences. 43 Foils with different tip-shaft angles were also tested, a reference value β0 and 10◦ under (β0− ) and

above (β0+ ) it. The results can be found in Figure 34 One can see that on a certain range of stability index (cants), the convergence time is relatively constant and similar, especially for β0 and β0− . The β0− case is slightly better however. One could have expected more substantial advantage. The case β0+ is quickly really unstable and some points – very small cants – have lead to complete divergence. Figure 34: Stability study for tip angles. 2.4 Heave: multi-parameter study As for the pitch stability study, such a 1-measure study cannot be completely trusted. The cost of static stability (via cant or tip angle) is known in terms of drag increase. But other aspects which have proven important in sea trials are at stake. This is especially the case of leeway From full scale trials, it was observed that the aimed static stability properties can have a high cost in terms of leeway. And there is no point in having a fully stable yacht if stability is reached by

strongly drifting – and consequently losing speed. Besides, heave stability has also a strong impact on heel stability: due to the lever arm corresponding approximately to half the ship’s beam, the foil vertical force is intimately linked to the ship’s heel as it provides righting moment. It can be seen from many footages of small foiling multihulls that negative heeling responses to acceleration for example are very common. As the speed increases, the foil lift increases quadratically while the aerodynamic heeling moment varies more slowly, which can lead the ship to negative heeling. In a similar way, the response to a heel positive or negative perturbation is mainly driven by the subsequent foil vertical force changes. Therefore, this study was extended to a two-parameter study with three perturbation cases: heave, negative heel and positive heel. The principle is to measure the convergence time as well as the leeway variation amplitude in response of those perturbations.

Different foil settings were simulated, as well as different board extensions, in line with the sea trials studied cases. As for the usual heave stability study, a stability index can be built that provides informations about the ship’s response following a heel perturbation: i Istab ∝ σi dFz dP i . (82) Fy where P i stands for the considered perturbation. σ i is ±1, depending on the sign of the heeling perturbation: in order to have a counteraction effect, it must have the same sign as the perturbation. 44 Similarly, the change in leeway following those perturbations can be computed at iso-Fy , showing the drift cost of the studied settings. A radar diagram can then be built showing both indexes for all three perturbation cases (see Figure 35). The left-hand part of the diagram shows the static leeway indexes, while the right-hand one shows the static stability indexes. Data are adimensionalised, the aim is to carry out comparisons and relative measurements, not to

obtain absolute values. For the leeway indexes, the centre of the diagram stands for the maximal changes (worst case), while the external line stands for no change. As far as the stability indexes are concerned, their values increase from the centre of the diagram to the external line. In this way, the closer the results are to the external line, the better the properties Figure 35: Static stability coefficients and leeway derivatives for different foil cants. + (blue: Γ− 0 , red : Γ0 , green: Γ0 ) From the static coefficients of Figure 35, one can see the known static distribution for which a larger cant provides better stability in heave (and also negative heel, as this case is rather similar, the foil immersion decreases). One can also observe the slight static advantage in having a more vertical foil (less cant) regarding leeway. Now that the index has been derived, simulations are launched in order to compare them to the actual dynamic response. Results are shown in Figure

36 The right-hand part of the diagram shows the convergence time – the closer to the external line the shorter the time –, the left-hand part shows the leeway variation amplitudes – as previously, the closer to the external line, the smaller the amplitude. Results are presented for three different cants: a + ◦ reference value Γ0 , 5◦ under (Γ− 0 ) and 10 above (Γ0 ). Figure 36: Dynamic stability coefficients and leeway amplitudes for different foil cants. + (blue: Γ− 0 , red : Γ0 , green: Γ0 ) 45 Although the convergence times for heave and negative heel perturbations are consistent with the stability indexes, the response is unexpected for positive heel as – unlike the indexes – the convergence times are better for larger cants. Besides, although the heave stability index are rather uniformly distributed (the gain in going from cant Γ− 0 to cant Γ0 is approximately the same as for going from Γ0 to Γ+ 0 ), there is a substantial gap between the

convergence time of cant Γ− 0 on one hand, and cants Γ0 and Γ+ 0 on the other hand. It is interesting to note that, conformally to the static leeway indexes, the variations between the different foil settings are extremely small. All cases show a leeway variation amplitude of approximately 1◦ . Finally, other settings changes can be compared. For instance, as leeway is now under close surveillance it is interesting to observe the consequences of changing the board extension. Figure 37 shows the diagrams for extensions 100% and 50%, with cant Γ0 . Unlike the previously presented cases, larger leeway changes are here observed (from 1.5◦ to 25◦ depending on the perturbations) It is also particularly striking to see that although the heave convergence time is almost unchanged, the stability properties are much poorer regarding heel perturbations. This shows how the board is a key element regarding stability. First it enables to limit the leeway increase due to the loss of

foil side force and second, as can be seen from the convergence times, it also makes the system find its equilibrium more efficiently. Figure 37: Dynamic stability coefficients and leeway amplitudes for different board extensions. (blue: 100%, red : 50%) Such study is useful to understand how stability is affected by the different couplings at stake. Linked to sea trials it also strengthen the confidence one can have in the results. It also provides a useful tool to compare different foil designs and optimise stability in an early appendage design phase As shown, the unsteady response may differ from the behaviour expected from static stability indexes. The comparison with full scale sea trials has enabled the isolation of the key parameters. One-parameter approaches are not sufficient and multi-variables tests are necessary. Finally this may, in a close future, allow a decrease of the margin that is taken when designing foils to ensure stability, certainly enabling performance

gains. 46 Conclusion Unsteady sailing yacht response is a developing field offering numerous possibilities. This report presented the current state of VPLP’s DVPP and the implemented aero- and hydrodynamic models. A list of the issues that remain to be tackled is established and solutions are discussed. A concrete example of a use of the simulator is made throughout the study of the relationship between usual static-approach stability criteria and effective dynamic response. In the near future, such numerical tools will enable to precisely simulate the yacht behaviour in key phases such as tacks and jibes. Manoeuvring issues and motion stability may also be addressed as well Unsteady motion in waves, seakeeping and routing optimisations are also possible fields of use of a dynamic simulator. Through a better understanding of the ship transient response, designers and engineers are able to imagine more and more performing sailing yachts concepts. And as often in high performance

offshore sailing, such innovations gradually spreads in the leisure industry for the benefit of all sailors. 47 References C. Böhm A Velocity Prediction Procedure for Sailing Yachts with a Hydrodynamic Model Based on Fully Coupled RANSE-Free Surface Simulations. PhD thesis, Delft University of Technology, 2014 O. M Faltinsen Hydrodynamics of high speed marine vehicles Cambridge University Press, 2005 O. M Faltinsen and R Zhao Numerical predictions of ship motions at high forward speed In Philosophical Transactions of the Royal Society of London, volume 334, pages 241–252, 1991. M. Fisher Recent foil design development In The foiling week, Malcesine, July 2014 F. Fossati Aero-hydrodynamics and the performance of sailing yachts Adlard Coles Nautical, 2009 F. Fossati and S Muggiasca Experimental investigation of sail aerodynamic behavior in dynamic In Transactions of Society of Naval Architects and Marine Engineers, volume 120, pages 327–367, 2013 K. Garme Marine hydromechanics

KTH Centre for Naval Architecture, 2011 J. Gerritsma Forced oscillation experiments In Philosophical Transactions of the Royal Society of London A, volume 334, pages 199–211, 1991 K. Graf, M Pelz, V Bertram, and H Söding Added resistance in seaways and its impact on yacht performance In 18th Chesapeake Sailing Yachts Symposium, Annapolis, 2007. P. Heppel Flight dynamics of sailing foilers In 5th High Performance Yacht Design Conference, Auckland, March 2015. J. E Kerwin and J N Newman A summary of the H Irwing Pratt ocean race handicapping project In SNAME 4th Chesapeake Sailing Yachts Symposium, 1979. R. Korpus Performance prediction without empiricsm: A RANS-based VPP and design optimization capability In 18th Chesapeake Sailing Yacht Symposium, volume 14, Annapolis, March 2007. J. N Newman The damping of an oscillating ellipsoid near a free surface Journal of Ship Research, 5(3): 44–58, 1961. J. N Newman The theory of ship motions Advances in Applied Mechanics, 18:221–283,

1978 Y. Roux, M Durand, A Leroyer, P Queutey, M Visonneau, J Raymond, J-M Finot, F Hauville, and A. Purwanto Strongly coupled VPP and CFD RANSE code for sailing yacht performance prediction In 3rd High Performance Yacht Design Conference, Auckland, December 2008. N. Salvesen, E O Tuck, and O Faltinsen Ship motions and sea loads In Transactions of Society of Naval Architects and Marine Engineers, volume 78, pages 250–287, 1970. R. T Schmitke Prediction of wave-induced motions for hullborne hydrofoils Journal of Hydronautics, 11(3): 93–99, 1977. R. T Schmitke Ship sway, roll, and yaw motions in oblique seas In Transactions of Society of Naval Architects and Marine Engineers, volume 86, pages 26–46, 1978. T. Theodorsen General theory of aerodynamic instability and mechanism of flutter Technical Report 496, NACA, 1935. J. H Vugts The hydrodynamic forces and ship motions in waves PhD thesis, TU Delft, 1970 J. V Wehausen and E V Laitone Surface waves Handbuch der Physik, volume 9 S

Flügge, 1960 R. W Yeung and S H Kim Radiation forces on ships with forward speed In Proc Third Int Conf Numerical Ship Hydrodynamic, volume 7, pages 499–516, 1981. R. W Yeung and S H Kim A new development in the theory of oscillating and translating slender ships In Fifteenth Symposium Naval Hydrodynamics, pages 195–218, 1985. 48 A Heaving foil in head sea The following part aims at proposing a small simplified model of the loads encountered by a foil heaving in waves. The oscillations and waves amplitudes are assumed small with respect to the shaft length and the case of a tip piercing the free surface is not considered. The added mass term is neglected As shown on Figure 38, the foil cant is Γ and the tip opening angle is β. The shaft immersed length at rest is l and the angle of attack at rest on shaft and elevator are respectively αs0 and αe0 . Aelev , celev , Ashaft and cshaft refer respectively to the elevator’s and shaft’s wetted area and main chord. The

considered waves have the form: φ = ζ0 e−kz cos(kx − ωt) To simplify the expressions, flow velocities are expressed at the same depth z for shaft and elevator, and ζz is defined as ζz = ζ0 e−kz . Furthermore, the notations C = cos(kx − ωt) and S = sin(kx − ωt) are used. Figure 38: Considered situation. Figure 39 shows the incident fluid speed vector in the shaft coordinate. This leads to an incident speed Vsi and an angle of attack αs : h i1/2 2 2 Vsi = (U − ωe ζz S) + (η̇3 + ωe ζz C) sin2 Γ  1/2 ωe ζz S ≈U 1− U η̇ sin Γ ωe ζz sin ΓC 3 αs = αs0 + + U U Figure 39: Dynamic angle of attack. by neglecting all second order terms. αs0 is the shaft angle of attack in steady motion. The loads in the vertical directions can be expressed as follow: z z F z = Fshaft + Felev , = Lshaft sin(Γ) + Lelev sin(β − Γ). The lift forces of each part can be expressed as:       1 η3 η̇3 sin Γ ωe ζz C dCLs ωe ζz S ζ0 C 0 Lshaft = ρcshaft l

+ − × αs + + × 1− U 2, 2 cos Γ cos Γ U U dα U     1 η̇3 sin(β − Γ) ωe ζz sin(β − Γ)C dCLe ωe ζz S 0 Lelev = ρAelev × αe + + × 1− U 2. 2 U U dα U 49 (83) (84) (85) Combining and sorting equations (83), (84) and (85) yields: 1 dC s dC e F/ ρU 2 = Ashaft L αs0 sin Γ + Aelev L αe0 cos(β − Γ) 2 dα dα dCLs 0 tan Γ + η3 × Ashaft α dα s l   dCLe sin2 (β − Γ) dCLs sin2 Γ + Aelev + η̇3 × Ashaft dα U dα U   s dCL ωe ζz 2 0 ωe ζz 0 ζ0 + Ashaft sin Γ C − αs sin Γ S + αs tan Γ C dα U U l   e dCL ωe ζz 2 0 ωe ζz + Aelev sin (β − Γ)C − αs sin(β − Γ)S dα U U Steady motion term Restoring term Damping term Wave loads term To get an order of magnitude of those values, one can assume dCLi /dα ≈ 2π and consider the following case: ζ0 = 1 m, ζz = 0.77 m, Te = 5 s, The loads are divided by 1/2ρAtotal U 2 : Load term Steady motion Restoring term Damping term Wave load 50 Coefficient 1.92 0.46

m−1 0.17 m−1 s 0.87C − 034S l = 2 m, U = 15 m/s