A-Class VPP
Marine

Draft report.

Sail Plan and Appendages Numerical Optimization Procedure for Multihulls

A numerical optimization procedure of catamarans sail plan and appendages is descripted. The method integrates a parametric geometry model, an automatic computational domain generator and a Velocity Prediction Program (VPP) based on a combination of CFD computations and analytical models. The sailing speed and course angle are obtained, with an iterative process, solving the forces and moment equilibrium system of equations. The performance of the boat is estimated by a model that provides a very fast evaluation of the resistance at a given velocity and any appendages configuration. The closure of the equations system is assured by the sail forces and the position of the aerodynamic center of effort provided by RANS analyses. The hull forces are modelled by empirical analytical formulations whose coefficients are tuned against a matrix of known solutions of the isolated demihull. The appendages aerodynamic polars are estimated applying preliminary design criteria from the aerospace literature. The procedure permits to find the combination of appendages configuration, rudders setting, sail planform, shape and trim that maximize the VMG. The method can be applied for the design of any type of multihull of any dimension and can be adapted to face any sail plan configuration, included rigid wing sails. The test case used for its development is the A–Class catamaran.

Nomenclature

alfa = Angle of incidence dCLda = Lift curve slope
aCL0 = Zero lift incidence dCLhdb = Hull side force curve slope
beta = Leeway angle Fb = Boat aerodynamic resistance
betaa = Apparent wind angle Fh = Sail heeling force
betat = True wind angle FM = Crew aerodynamic resistance
gamma = Rudder setting angle Ft = Sail thrust force
delta = Foils dihedral angle h = Foils spanwise aerodynamic center
epsilon = Dagger board twist hb = Height of boat center of gravity
teta = Mast side setting angle hh = Arm of sail aerodynamic resultant
lambda = Aspect ratio k = Demihull shape factor
roa = Air density L = Lift
row = Water density LH = Hull side force
fi = Heeling angle lm = Arm of crew righting moment
CD = Drag coefficient r = Dagger board root stagger angle
CD0 = Drag coefficient at zero lift SH = Demihull reference surface
Cf = Skin friction coefficient SS = Sail reference surface
CL = Lift coefficient Swet = Demihull wet surface
CL0 = Lift coefficient at zero incidence V = Boat speed
CT = Hull total resistance coefficient Va = Apparent wind speed
CW = Residual resistance coefficient VT = True wind speed
D = Drag WBE = Boat empty weight
DH = Hull drag WBO = Boat operative displacement
dmin = Distance between hulls centreline WM = Crew weight

Performance prediction model

Analytical formulations for hull and appendages are coupled with the CFD based sail aerodynamic model in order to develop a function that provides the performance of the boat (in term of VMG) for a given sail plan and appendage configuration.

Fig. 1

Fig. 1: Forces acting on the boat and reference frame

Hull forces model

The hull forces are modelled by analytical formulations both taken from literature and developed by a comparison with a matrix of RANS computations solutions at several velocities, attitudes, displacements and leeway angles performed on the isolated demihull.
Two formulations were adopted for the side force and for the resistance.

Side force model

The hull is modelled as a lifting body. The side force (force laying in a plane parallel to the water plane and normal to the sailing direction) is modelled as:

LH

where the reference surface SH is the is the side projection on the symmetry plane of the submerged part of the demihull. It changes with displacement and is estimated by CAD. According to the geometry used to setup the method here developed, the analytical formulation that models its relation with displacement can be split in two regions: an exponential grow starting from zero displacement and a linear behaviour after a defined input WBO0 value of displacement:

SH

where ksh1 and ksh2 are computed knowing two surface values in the linear region and taush is estimated adding a known value in the exponential region.
From the matrix of CFD solutions obtained on the reference demihull, the lift curve slope dCLhdbeta seams to change with displacement, velocity and leeway angle. The proposed formulation developed to best match the CFD solutions at the several conditions is:

dCLHdb

Substituting it in the hull side force expression we obtain:

LH
eq. 1

Parameters kh1, kh2, tauh1 and tauh2 are determined from the known solutions of the isolated hull.
The graphs in Figure 2 compare the analytical solutions at two velocities and three leeway angles with a matrix of CFD solutions obtained on a reference demihull geometry.

Fig. 2

Fig. 2: Comparison between demihull side force analytical model and CFD solutions

Drag model

The bare demihull total resistance is expressed by:

DH

The hull wet surface Swet grows with displacement. The values are evaluated by CAD assuming the submerged volume to be determined cutting the hull with a flat plane representing the water surface. The model is similar to that used to model SH:

Swet

As for the SH model, the value of displacement WBO0 from which the curve begins to be linear depends on the geometry and has to be evaluated by CAD. The terms ksw1, ksw2 and tausw are found knowing the hull wet surface at two displacement values in the linear an one in the exponential region. The reason to model the reference surfaces starting from zero displacement is to provide the procedure the capability to analyse also configurations in which the lifting contributions of the foils are predominant to the hull’s one.

Fig. 3

Fig. 3: Demihull wet surface model

The total resistance coefficient can be modelled by a combination of a friction and a residual component [1].

CT

The form factor 1 + k accounts to the over velocity generated by the thick shape of the body [2]. Its value has to be evaluated from literature or from a known bare hull drag value.
The skin friction coefficient is estimated according to the ITTC–57 friction line expression:

Cf

where RN is the Reynolds number referred to the hull waterline length. Good agreement with CFD computations were evidenced (in figure 4 it is plotted in function of Froude number).

Fig. 4

Fig. 4: Skin friction coefficient

A significant simplification was chosen to model the residual drag coefficient. A combination of two exponential formulations, for speeds lower and higher than a critical value and function of displacement and Froude number, was adopted. The coefficients are tuned against the matrix of known demihull performance solutions. The proposed formulations are:

Cw

The parameters ww, kw1 and kw2 have to be tuned to match the set of known solutions. The value of FNcr is assumed to be the critical Froude number FNcr04 but, if considered more appropriate, its value can be chosen to best match the matrix of known drag solutions. In the case of the reference A Class catamaran, for instance, the critical Froude number of 0.4 is at 5.7 knots of speed but it was preferred to model the wave drag with a value of FNcr equal to 0.46 to corresponds a speed of 6.6 knots (Figure 5).

Fig. 5

Fig. 5: Residual drag coefficient model

Figure 6 compares the computed and modelled viscous and residual drag. The strong simplification adopted in the residual component modelling is justified by its lower influence respect to the friction component.

Fig. 6

Fig. 6: Computed and modelled drag breakdown for the reference hull

An additional factor that takes into account for the drag increment due to leeway angle was added. Such dependency was assumed to be quadratic with leeway angle. From CFD computations it also showed to be linearly dependent to displacement and exponentially to velocity. The proposed factor, to be multiplied to the drag model, is:

additional factor

The final expression of total resistance is then:

DH
eq. 2

The parameters to tune are kbeta, taubeta and wbeta.
Figure 7 compares the modelled hull drag increase due to the leeway angle with the CFD computations at two velocities.

Fig. 7

Fig. 7: Drag increase with leeway angle

Appendages forces model

Dagger boards and rudders are modelled as wings. The lift and drag formulations are defined applying preliminarily design criteria from aerospace literature.

Foils lift expression

The formulation of foils lift is:

L
eq. 3

The lift coefficient, in the linear region of the lift curve of a non–symmetric foil, can be expressed as:

CL

The effective velocity Veff is the component of the boat velocity normal to the foil leading edge (for rectangular planform) and aeff is its angle of incidence. Veff is the only velocity component responsible to the generation of lift (the friction contribution can be neglected). The spanwise component, in fact, does not influence the lift but only cause a shifting of the boundary layer [3].
From geometrical consideration, the effective velocity component Veff, of the downwind (unswept) dagger board for instance, is given by the following relation:

Veff

which, for low values of the leeway angle beta, can be approximated to the boat speed V.

Veff

The effective angle of incidence is related to the leeway angle by the relation

tanaeff

For small values of aeff and beta we can approximate

tanaeff

The effective incidence can then be written as

aeff

From the above considerations and referring to the Figure 1, the lift coefficients of dagger boards and rudders are expressed in function of beta as following:

CL
eq. 4

The foils lift curve slopes are obtained applying empirical formulations used in aeronautics in the predevelopment phase [4]. Assuming a linear twist and constant airfoil section along the full span, the 3D lift curve slope (with dimension invdeg) is modelled as:

dCLda
eq. 5

where p, b and lambda are respectively the foil planform perimeter subtracting the root chord length, the draft and the aspect ratio of the mirrored full span wing. The term dCLda2D is the 2D lift curve slope of the airfoil.
The dagger boards zero incidence lift coefficient CL0D is obtained solving the lift curve equations for CLD-0 (first two expressions of eq. 4) substituting to the factor between squared brackets the 3D zero lift incidence which, for the downwind dagger board for instance, is given by:

aCL0

in which epsilon is the foil twist and r the root stagger angle. The 2D lift curve slope of the airfoil and the zero lift incidence aCL02D, as well as its drag polar, are provided from the airfoil aerodynamic database.
The lift forces of all foils are obtained substituting the lift coefficient expressions (eq. 4) in the lift formulation (eq. 3).

Foils drag expression

The wing drag formulation is:

D

where the drag coefficient expressed as function of lift coefficient is:

CD
eq. 6

The values of J, e, vmin and w are reported in literature [5] as function of aspect and taper ratio. In order to integrate the models in the automatic procedure, analytical formulations that approximate the tabled values for rectangular and elliptical foils planforms were developed. The above expression, applied to dagger board and rudder, provides their drag coefficient formulation in function of lift coefficient.
The 2D airfoil aerodynamic data required by the model can be provided by an external database or computed (if the airfoil is a variable of design). An automatic procedure to characterize any 2D geometry, based on the XFoil panel/boundary layer solver, was implemented and integrated in the loop.

Boat global forces and moment equilibrium

The forces equilibrium equations of the complete boat, referred to a frame with the X axis aligned with the sailing direction and Z axis perpendicular to the water plane, are:

X equilibrium

Dtot
eq. 7

Y equilibrium (assuming to neglect the lift components of crew and boat)

FH
eq. 8

Z equilibrium

WM
eq. 9

The heeling moment equilibrium around the center of buoyancy of the downwind hull gives, assuming to neglect the influence of the mast setting on the center of gravity (it was introduced a variable assuming the possibility to move the mast upwind or downwind), is expressed as:

WBE
eq. 10

The yaw and pitching moment equilibrium were neglected at this stage.
Developing the eq. 7, eq. 8, eq. 9 and eq. 10, substituting the appendages lift and drag forces formulation (eq. 3, eq. 4 and eq. 6) and including the hull side force equation (eq. 1), we obtain, assuming the boat velocity V and the height of sail aerodynamic center of effort hh to be given as input, a system of five equations and five unknowns (Dtot, Fh, WBO, LH and beta). The solution of the equations system is implemented as a script function that produces as output the boat total resistance Dtot and sail heeling force Fh (which are the parameter we are interested in to be compared with the CFD sail solution) at a given speed V, center of effort height hh and set of parameters characterizing the boat configuration. Figure 8 outlines the inputs and the outputs of the function.

Fig. 8

Fig. 8: Scheme of the equilibrium equations function

Sail aerodynamic model

The sail heeling and thrust forces, in windward sailing conditions, together with the resultant aerodynamic center of pressure, are obtained by a steady fully turbulent MRF Navier–Stokes analysis. The two equation k – ω Shear Stress Transport (SST) of Menter turbulence model [6] was used and Wall Functions applied to model the wall boundary layer. In order to include the boat interference effect on the sail, the complete geometry, included the crew, is modelled. The boat is moving on a local reference frame at the given speed V and direction betat respect to the true wind. The wind boundary layer velocity profile, in the absolute reference frame, is imposed at farfield [7].
Figure A1 in appendix reports the CFD solution obtained on an A–Class catamaran sailing windward with 10 knots of true wind speed at 10 meters from the water plane. The left contour plot details the structures of the sail tip and root vortices. The wind boundary layer velocity magnitude is reported on a plane behind the boat. The right plot evidences the evolution of the separation behind the mast.

Boat performance solution

The sail aerodynamic solution provides the required closure of the boat performance prediction problem. The CFD analysis is coupled with the boat analytical model to solve the equilibrium system of equations.
The procedure starts with a first sailing speed V and direction betat guess. A CFD analysis, with the selected sail plan, shape and trim, is then run at this conditions. Sail forces and center of effort are extracted and used to verify the equilibrium system. The verification consists in checking if the boat total resistance and the sail heeling force, computed by the analytical model (for the given hull and appendage configuration), are equal respectively to the sail thrust force and the heeling force deriving from the CFD computation:

DtotFh

If not, new values of boat speed and true wind angle are selected. The CFD computation is restarted at the new conditions and the procedure is repeated until the equilibrium equations criteria are verified (with a prescribed tolerance). The VPP problem is completed with the production, as output, of the “Velocity Made Good” (VMG) value, which represents the performance of the boat [8] with the prescribed sail geometry (considered rigid) and appendages.
In order to speed up the VPP solution convergence, the procedure of sailing conditions exploration is split into two nested cycles. The principle is to use the external loop, which involves the RANS computation, to model analytical polars formulations of the sail aerodynamics that indicate the inner search algorithm the direction where to find the sailing conditions that verify the equilibrium. The estimated sail aerodynamic model is then refined every external cycle until equilibrium is verified in both loops.
The analytical polars formulations used to model the sail aerodynamics are similar to the one adopted to model the appendages. The sail lift and drag coefficients are expressed as:

CLS
eq. 11
CDS
eq. 12

The term eS is the sail induced drag factor (always lower than 1) or “Oswald efficiency factor” and is related to the shape of the spanwise load distribution. For preliminary design purpose it is reported as function of aspect and taper ratio. The coefficient kS is related to the sail camber.
The thrust and heeling forces are expressed in function of lift and drag coefficients by the equations system:

FtFh
eq. 13

The apparent wind speed Va and angle betaa are obtained as function of the true wind VT (for convention usually measured at 10 meters from the sea) and its angle betat by the relations [9]:

betaA
eq. 14
VA
eq. 15

Substituting the eq. 14, eq. 15, eq. 11 and eq. 12 in the system of eq. 13 we obtain the formulation of Ft and Fh in function of V and betat, that will be used in the verification criteria of the equilibrium equations system.
The modelled sail aerodynamic drag polar is a quadratic formulation with lift coefficient. It than requires three solutions to be estimated. Its progressive update follows different criteria during the first three iterations of the external cycle:

  • In the first iteration the sail lift curve slope dCLdaS and the induced drag factor eS are estimated from literature as function of sail aspect and taper ratio. The value of zero–lift drag coefficient CD0S is roughly imposed from experience. The sail lift and drag coefficients CLS and CDS, obtained from the CFD analysis, are used to the estimate CL0S eq. 11 and kS from eq. 12.
  • In the second iteration the additional CFD solution is used to complete the analytical lift curve formulation adjusting the values of the lift curve slope dCLdaS and zero–incidence lift coefficient CL0S. The parameters updated in the polar curve are CD0S and kS while the value of eS is still guessed.
  • In the third iteration the analytical drag polar formulation is completed with the computation of the induced drag factor eS which is last unknown parameter. The lift curve is updated connecting a quadratic formulation to the previous computed linear part.
  • In all the following iterations the sailing condition estimation are performed modelling the polars regions under investigation updating both curves by a generic quadratic formulation using the closest three solutions.

The Figure 9 reports, for a typical A–Class catamaran sail plan, an example of the evolution of the polars computation during the progress of the first three iterations and the estimation of the values to be used for the computation of the sailing conditions in the fourth iteration.

Fig. 9

Fig. 9: Example of sail analytical polars computation progress

If the sail aerodynamic conditions fall in the linear region of the lift curve three iterations are in general sufficient to converge. If not, the reported analytical polars formulation are no more valid. The quadratic formulations, with which the non–linear sail aerodynamics is modelled in the following iterations, simply constitute interpolating curves whose coefficients have no particular physical meaning. The complete workflow of the performance prediction model is described in the flow chart of Figure A3 in appendix.
The searching criterion of the aerodynamic coefficients in the inner cycle is driven by an optimization procedure, based on the Nelder–Mead Simplex algorithm [10], whose objective function is the minimization of the differences between the forces derived from the boat analytical model and CFD computation in the equilibrium verification criteria:

ObjFun

The values of V and betat that satisfy the equilibrium are used in the next external loop iteration where the sail polars are update. The iterations continue up to the satisfaction of the equilibrium system in both loops. When convergence is reached, the boat VMG is computed and produced as output.
It was experienced that with this procedure the RANS computations required to obtain a convergence rarely was higher than four of five (if sail is not stalled or, in general, if separations are not too large). Furthermore, a restarting procedure from the previous solution and a progressive reduction of the CFD solver number of iteration, was implemented. This strategy showed to be very efficient in boosting the convergence but its robustness is related to the capability to select starting sailing conditions as realistic as possible. The procedure fails in case of sudden sail separation. A check if complete stall occur has been implemented in order to reject such solutions.

Optimization environment development

The above described performance prediction procedure has been integrated in an optimization environment in which the optimal sail plan, trim and appendage configuration is searched. The method integrates, in an automatic process, the sail parametric CAD model, the computational domain generation module and the VPP model.

Parametric CAD model

Both the master sail plan geometry and the architecture of the structured computational domain is generated and made parametric in the CAD model. Spline curves and points, which will constitute the edges and the vertices of the mesh blocks, are updated together with the sail geometries. This strategy provides a great control of the grid topology shape and reduces the risk of failure in the computational domain generation.
The A–Cat sail plan consists in a single mast/main–sail configuration. The CAD parameters were selected with the aim to investigate the largest possible range of geometries. Traditional sail plans, wing masts or rigid sails with a small portion of flexible sail can be generated. The sail is built by a loft surface between the luff and leech curves which are used as guides. The foot curve, an intermediate arbitrarily positioned section and the head are used as control sections. In the similar manner, the mast is generated from three geometries at the same sail stations. The planform is controlled by reference surface, aspect ratio, taper ratio and by other parameters (as tip and root inclination, luff and leech camber in the sail plane etc…) that give the possibility to investigate any kind of shape. The examples in Figure 10 give the sense on the flexibility of the parametric model.

Fig. 10

Fig. 10: Examples of possible sail planforms

The sail sections are modelled by cubic Bezier curves. The first point of the control polygon is connected to the mast luff, the last one coincides with the leech of the sail. The four coordinates of the two intermediate control points are parameters of the geometry (red polylines in Figure 11). The mast sections are generated by spline curves controlled by three parameters: the tangent tension at the leading edge, the tangent tension at the luff point and the angle between the latter tangent and the symmetry plane of the mast. The spanner and the three sections angle are sail setting parameters. The input reference surface area is kept unchanged. After the geometry creation, the final sail area is measured and the loft surface cut in order to restore the required value.
The CAD model is linked to more than 40 parameters. All of them are available to be used as design variables of an optimization problem.

Fig. 11

Fig. 11: Examples of possible mast–sail sections geometry

Computational domain

The domain has to be updated rapidly and automatically every iteration. The grid generation using a structured hexahedral topology fulfils this requirement. Such topology is however very demanding to be created on complex geometries. A mixed strategy was then adopted. The domain was divided in several regions with common boundaries and each region meshed applying the more opportune strategy. A structured CH–grid topology was created in a limited volume around the sail and dimensioned to envelope the full range of possible geometries (Figure 12).

Fig. 12

Fig. 12: Structured hexahedral mesh around the sail rig

Due to the boat and crew shape complexity, an unstructured hybrid prism/tetra mesh was generated in a volume between the sail structured mesh boundaries and the water plane. The sail/boat mesh assembly is contained in a box, four boat lengths large and tall, in which another hybrid prism/tetra mesh portion was generated. The remaining volume was meshed with hexahedral cells growing toward the top with a progression aimed to better model the inflow air boundary layer. The full domain is 10 boat lengths wide and is extended 10 boat lengths upstream and downstream the model. In the final domain resulting from the assembly, the several parts are connected by inner zonal interface boundaries in which a simple “flow through” condition between the non–conformal mesh zones is imposed. Figure 13 shows the assembly of the parts with the interfaces in evidence.

Fig. 13

Fig. 13: Details of complete computational domain assembly

This inner structured sail grid is the only mesh part subjected to update every optimization iteration. After the CAD model update, the process is guided by a batch procedure which loads the new CAD file, opens the predefined mesh topology, projects the blocks vertices, edges and faces on the new geometry, recomputes the grid, assembles the new part with the rest of the unchanged domain and export the final mesh in the CFD solver format.

Optimization procedure

A schematization of the optimization cycle is reported in Figure A2 in appendix. A set of batch procedures are in charge of generating the CAD model, of updating the mesh and of managing the performance prediction module descripted in the flow chart in Figure A3 which provides as output the boat performance in term of VMG. The optimization algorithm selects a variables combination every iteration and the cycle progress up to a convergence to the configuration that should represents “the best” boat. The variables of design can be all the parameters that control the sail plan geometry and all the parameters that characterize the appendage configuration (surfaces, aspect ratio, dihedral, etc…)…

…………………………….

To be continued.!!.. Report under development…

 

Appendix

Fig. A1

Fig. A1: Contour plot of CFD solutions on an A–Class catamaran

 

Fig. A2

Fig. A2: Scheme of the optimization procedure

 

Fig. A3

Fig. A3: Flow chart of the boat performance prediction model

Bibliography

[1] M. Insel and A. Molland, “An Investigation into Resistance Components of High–Speed Displacement Catamarans”, Transaction of the Royal Institute of Naval Architects, 134, pp. 1 – 20, 1992.

[2] S. F. Hoerner, “Fluid&nd–amic Drag”, Bricktown, New Jersey: Hoerner Fluid Dynamics, 1965.

[3] V. Losito, “Fondamenti di Aeronautica Generale”, Accademia Aeronautica, 1991.

[4] J. Roskam and C. T. E. Lan, “Airplane Aerodynamics and Performance”, DARcorporation, 2000.

[5] I. H. Abbott and A. E. Von Doenhoff, “Theory of Wing Sections”, New York: Dover Publications, 1959.

[6] F. R. Menter, “Two–Equation Eddy–Viscosity Turbulence Models for Engineering Applications”, AIAA Journal, vol. 32, n. 8, pp. 1598 – 1605, August 1994.

[7] A. Musker, “Explicit Expression for the Smooth Wall Velocity Distribution in a Turbulent Boundary Layer”, AIAA Journal, vol. 17, pp. 655 – 657, 1979.

[8] L. Larsson and R. E. Eliansson, “Principle of Sailing Yacht Design”, London: Adlard Coles Nautical, 1997.

[9] A. Claughton, J. Wellicome and A. Shenoi, “Sailing Yacht Design – Theory”, Longman, 1998.

[10] J. A. Nelder and R. Mead, “Simplex Method for Function Minimization”, The Computer Journal, vol. 7, p. 308 – 313, 1965.