High Performance Yacht Design Conference Auckland, 46 December, 2002 An Advanced finite element method for fluid dynamic analysis of America’s Cup boatS J. GarcíaEspinosa^{1}, julio@compassis.com R. LucoSalman, M. Salas^{2}, rluco@uach.cl M. LópezRodríguez^{3}, natutatec@nautatec.com E. Oñate^{4}, onate@cimne.upc.es
Abstract. This paper shows a step forward in the development of a stabilized method for solving the Reynolds equations (RANSE) including free surface effects. The starting point of this method is the modified governing differential equations for an incompressible turbulent viscous flow and the free surface condition, incorporating necessary stabilization terms via a Finite Increment Calculus (FIC) procedure. Time integration scheme of the method is based on an implicit monolithic second order method. Implementation of the algorithm, based on the Finite Element Method (FEM) using unstructured grids of linear tetrahedra, allows us to take into account dynamic sinkage and trim, being especially adequate for hydrodynamic analysis of racing boats. This paper also presents a validation study of numerical results obtained for America’s Cup boats.
1. INTRODUCTION In recent years, more traditional tools for designing racings boats, especially America’s Cup or IMS boats, have been augmented by the advent of advanced numerical techniques. Even though CFD techniques are still under development they are becoming essential for designing boats with the highest performance [^{1},^{2}].
This paper shows a step forward in the development of stabilized Finite Element Method (FEM) to solve the Reynolds equations (RANSE) including free surface effects. The starting point of the method is the modified governing differential equations for an incompressible turbulent viscous flow and the free surface condition, incorporating necessary stabilization terms via a Finite Increment Calculus (FIC) procedure [Error: Reference source not found,^{3},^{4},^{5}]. Time integration scheme of the method is based on implicit monolithic second order method originally proposed by Soto et al. [^{6}]. This scheme is derived by splitting the momentum equation in a similar way as in an implicit fractional step method.
The application of the method, based on unstructured grid to enhance geometry flexibility, has been specially designed for its application in naval hydrodynamics. The implicit formulation is very efficient for free surface flows, where critical time step is some order of magnitude smaller than time step required to obtain time accurate results. The implementation presented here also allows us to take into account real free surface deformations and dynamic sinkage and trim of a boat by simply updating the analysis domain.
This paper also presents a validation study of numerical results obtained for an America’s Cup racing boat. The results comparison with a towing tank testing data includes drag and lift forces and wave profiles. 2. FIC FORMULATION We consider the motion around a body of a viscous incompressible fluid including a free surface. The stabilized Finite Increment Calculus (FIC) form of governing differential equations [^{7}] for a three dimensional problem on domain Ω, can be written as:
(1a)
(1b)
(1c)
, where
In the above, u_{i} is the velocity along the i^{th} global reference axis, p is the dynamic pressure, is the wave elevation defined on boundary _{} where free surface condition is applied, and _{ij} are the viscous stress tensor components. Let n_{i} be the unit outward normal to the boundary and denoting prescribed values by an overbar, boundary conditions for the stabilized problem are:
(2a)
(2b)
(2c)
, for t (t_{0},t_{f}). The boundary Γ has been considered split into three sets of disjoint components Γ_{u}, Γ_{p} and Γ_{} , the last being the part where mixed conditions are prescribed. Vectors g_{i} and s_{i} span the space tangent to Γ_{}. Finally, Γ_{u} and Γ_{p} are the two disjoint components of Γ, where Dirichlet and Neumann boundary conditions for the velocity are prescribed. Initial conditions have to be appended to problem (1)(2).
Underlined terms in eqs. (1)(2) introduce the necessary stabilization for the numerical solution, as shown in [Error: Reference source not found]. Additional time stabilization terms can be accounted for in eqs. (1)(2), although they have been found unnecessary for the type of problem solved here [Error: Reference source not found,Error: Reference source not found].
The characteristic length distances h_{j} represent dimensions of the finite domain, where balance of momentum and mass are enforced. The characteristic distances h_{}_{j} in eq. (1c) represent dimensions of a finite domain surrounding a point, where the velocity is constrained to be tangent to the free surface. Details of the derivation of eqs. (1)(2) and recommendations for computation of the stabilization parameters can be found in [^{8},^{9}].
Eqs. (1)(2) are the starting point for deriving a variety of stabilized numerical methods for solving the incompressible NavierStokes equations. It can be shown, that a number of stabilized methods allowing equal order interpolations for velocity and pressure fields can be derived from a modified form of the momentum and mass balance equations given above [Error: Reference source not found,^{10},^{11},^{12}]. 3. IMPLICIT FRACTIONAL STEP SCHEME Let us discretize in time the stabilized momentum equation (1a) using the trapezoidal rule (or θ method) as
(3)
, where superscripts n and θ refer to the time step and to the trapezoidal rule discretization parameter, respectively. For θ = 1 the standard backward Euler scheme is obtained, which has a temporal error of 0(t). The value θ = 0.5 gives a standard Crank Nicholson scheme, which is second order accurate in time 0(t^{2}). An implicit fractional step method can be simply derived by splitting eq. (3). The resulting continuous problem, omitting the boundary and initial conditions for brevity, is as follows:
(4a)
(4b)
(4c)
Terms denoted by an overbar are those calculated with an intermediate velocity ū_{i}, which is introduced to allow the momentum splitting. The scheme (4) is stable for any time increment t as shown in [^{13}]. The error due to taking implicit advective and viscous terms in (4a) can be shown [Error: Reference source not found] to be of the same order than the error of the stabilizing terms and therefore, it has the same order of approximation than the original time discretization (3). 4. MONOLITHIC STABILISED FORMULATION At this point, it is important to introduce an associated matrix structure corresponding to the variational discrete FEM form of (4) (see [Error: Reference source not found] for details of this derivation):
(5a)
(5b)
(5c)
, where U^{n+k}, P^{n+k} are the vectors of the velocity and pressure, evaluated at time step n+k. Terms denoted by overbar identify the intermediate velocity.
By taking Ū^{n+1} from (5c) and inserting the result in (5a)(5b), the following system of equations is obtained:
(6a)
(6b)
The term E(U^{n+θ}) in (6a) is the error coming from the implicit treatment of advective and viscous terms, which is of order 0(t^{2}). However, such term can be eliminated as in (6a)(6b) by writing the following analog monolithic scheme:
(7a)
(7b)
, where index i refers to the iteration of the iterative solution. Basically, in this final formulation the convergence of the block uncoupled solution is enforced by the first term of (7b), while the pressure stability is attained by the second term of the same equation.
Finally, the free surface elevation (B) solution is included in the monolithic scheme, obtaining the following system of equations:
(8a)
(8b)
(8c)
Remark A: The effect of changes in the wave elevation are introduced in eq. (8a)(8b) as prescribed viscous traction and dynamic pressure acting on the free surface, respectively [Error: Reference source not found]. Remark B: If necessary, turbulence equations are included in the scheme in a similar way as the free surface equation. The resulting scheme is a stable, second order algorithm to solve RANSE equations. 4. NUMERICAL ASPECTS The abovepresented algorithm has been implemented within the CFD environment Tdyn [^{14}]. This software has been used to solve the examples shown in the following section. Tdyn is based on the FEM and may use almost any type of element, but it is optimized for using linear tetrahedra. An unstructured grid is used in Tdyn to enhance geometry flexibility and to speed up the initial modeling time. 4.1 Unstructured grid generation Tdyn includes a pre and postprocessor completely integrated, based on the GiD system [^{15}]. This system permits the easiest geometry definition, mesh generation, and inserting necessary data definition for an analysis, as well as further results post processing.
The discretization of a general threedimensional computational domain into an unstructured assembly of tetrahedra is accomplished by means of an advancing front grid generation procedure.
This procedure requires the geometry of the computational domain to be defined in terms of an assembly of surface patches. In this case the surface definition of a complete computational domain, consisting of hull and appendages surfaces, free surface, inflow plane, exit plane, and bottom and lateral surfaces is based on NURBS patches. 4.2 Boundary conditions Tdyn preprocessor allows boundary conditions to be assigned directly on the geometrical entities and automatically transferred to the grid. This utility permits not to reassign boundary conditions every time a new grid is generated. Furthermore, boundary conditions may be defined by analytic functions. This fact allows performing different drift angle analyses using the same grid, by changing the inflow condition. Other kind of phenomena like nonuniform flows may also be simulated in a similar way. 4.3 Control Volume Finite element solution of eqs. (1)(2) on a fixed control volume is in most cases accurate enough for design purposes. However, a great quantitative performance of results may be obtained in some cases by updating the domain taking into account free surface deformations and a dynamic sinkage and trim of the boat. The procedure used in Tdyn to update the domain is now summarized.
1) The hull model is first considered in an estimated steady state position. A planar or quasiplanar surface is used as reference for free surface calculations. An automatic unstructured grid generator based on the advancing front technique is used to generate a surface mesh and a volume mesh. A steady state simulation is performed with the hull in this initial position. 2) The net heave force and trim moment acting on the hull are calculated from the previous converged solution. The sinkage and trim corrections required by the equilibrium of this force and moment are evaluated. Free surface is updated accordingly to previous results. This is made at CAD level, generating a new NURBS surface based on the previous triangulation of the free surface. A new domain is then created within Tdyn preprocessor, repositioning the hull and using the new free surface. This process is automatic, but if necessary may be controlled by the user by means of a wizardtype tool. Finally, a new mesh is automatically generated. 3) A steady state simulation is performed again in the new domain. A converged free surface is obtained for this given hull position at the end of the present step. This process should be repeated until a convergence of the results. Experience shows that one iteration is enough in most of the cases to obtain forces results within the uncertainty band.
T Figure 1. Rioja de España. Geometry definition used in the simulation and experimental testing.
he sinkage and trim corrections are expressed in terms of the net heave force and trim moment using the following relations:
z is a correction of the sinkage at a center of gravity, is a trim angle correction, F_{z} and M_{y} are a net heave force and a trim moment. A_{wp} is the water plane area, and Iy is the corresponding moment of inertia about the y axis. The heave force and the trim moment are defined in terms of the pressure p and the viscous stress tensor components _{ij}, which can be obtained directly from the flow solver. New free surface NURBS definition, taking the resulting deformation into account, is generated in three steps: NURBS Cartesian support grid of MxN points is created. M and N are calculated in terms of the number of local maximum in the X and Y axis directions. Z coordinate of the points, representing the wave elevation, is interpolated into the grid. This interpolation is based on a weighted function of the nearest points. The nearest points are easily located by using a quadtree structure. Finally, the NURBS surface based on the support grid is generated [^{16}]. Boundaries of the NURBS are defined by projecting the original ones in the Z direction. 6. APPLICATION EXAMPLE. RIOJA DE ESPAÑA. Following, some examples of the presented algorithm application are shown. The example presented in this section have been generated using the Tdyn wizard tool, starting from the standard NURBS definition of hull and appendages. This tool allows automatic control volume creation, data insertion and mesh generation in an user controlled manner.
The analysed example is the Spanish America’s Cup boat Rioja de España, participant in the races of 1995. The geometry of the boat is presented in Figures 1 and 2. It is based on a standard NURBS surface definition of the ship, used in the design process. Numerical results are compared with an extrapolation of experimental data obtained in CEHIPAR model basin (Spain). The model was produced at scale 1/3.5 by CAM over foam and faired by hand with final laser dimensional control.
 Figure 2. NURBSbased geometry used in the analysis of the Rioja de España boat. 
The resistance tests were performed with the model free to sink and trim. The dynamometer employed was a Kemp & Remmers R50 model with six components. Experimental data available include drag, lift, sinkage, trim angles, and wave profiles at 4.27m (real scale) from symmetry plane. Wave elevation measurements were performed with an ASM wire powermeter of 1250mm length and 0.8mV/V/mm sensibility. Figure 3. Detail of the final (boundary) mesh used in the E0D0 case (adapted, taking into account sinkage, trim and free surface deformation)
A list of the studied cases is shown next.
Test  Geometry  Heel  Drift  C0D0  Hull, no appendages  0º  0º  E0D0  Hull, bulb and keel  0º  0º  E15D2  Hull, bulb, keel and rudder  15º  2º  E25D4  Hull, bulb, keel and rudder  15º  4º  E25D2  Hull, bulb, keel and rudder  25º  2º 
Every case was towed at equivalent velocities of 10, 9, 8.5, 8.0, 7.5 and 7.0 knots.
Numerical analysis of cases were carried out at real scale. Characteristics of unstructured grids of linear tetrahedra used in those analyses are shown in the following table:
Test  Symmetry  # Elements  # Nodes  C0D0  Yes  300 000  70 000  E0D0  Yes  700 000  175 000  E15D2  No  1 500 000  380 000  E25D4  No  1 500 000  380 000  E25D2  No  1 500 000  380 000 
All grids used for this validation work have been generated with the same quality criteria and using element sizes from 5mm to 2000 mm. Some details of the grid used in the E0D0 case are shown in Figures 3 and 4.
The different analyses have been performed with a two layer k turbulence model, in combination with an extended law of the wall.
 Figure 4. Detail of the mesh used in the analysis of the E0D0 case, around keelbulb union 
Time step used in the calculations was t = 0.1s and 500 time steps were enough in all the cases to achieve a steady state. Additional calculations were carried out in some of the cases with t = 0.025s and t = 0.01s, in order to verify the influence of this parameter in the accuracy of the results. No significant influence was detected for the used values.
Several graphs of comparison of simulated and experimental wave profiles are presented in Figures 5
to 9. The origin of the x axis has been taken at the fore perpendicular and the x+ sense is afore. In the non symmetric cases, the measurements were performed at the opposite side of the heeled board. Ratio of maximum amplitudes of fluctuations (noise) and waves in the experimental measurements varies from 4% to 12%. It is observed that the calculated results agree well with the experimental data. However, a hollow of the wave elevation around x = 12 is slightly underpredicted in every case. This effect is observed in the wave profiles comparisons for cases above 8 Kn. The reason for the discrepancy in this region is not clear.
 Figure 5. C0D0 10kn. Wave elevation profile. 
 Figure 6. E0D0 10kn. Wave elevation profile. 
 Figure 7. E15D2 10kn. Wave elevation profile. 
 Figure 8. E15D4 10kn. Wave elevation profile. 
 Figure 9. E25D4 10kn. Wave elevation profile. 
Unfortunately, no wave elevation experimental data is available out of those profiles and no further study could not to be performed.
Figures 10 to 13 show pressure contours on bulb and keel in the different cases, corresponding to a velocity of 8 kn. These pressure contours show typical characteristics of this kind of analysis. They display higher pressure at stagnation area close to the leading edge and a significant lowpressure area around keel laterals. Position of the stagnation area, and lowpressure zone peak value are obviously dependant on heel and drift angles of the case.
Figures 14 to 17 show, in some of the simulated cases, obtained wave pattern, corresponding to a velocity of 9 kn. Systems of divergent and transverse waves are clearly observed, as well as the effect of the stern on the wave field. The asymmetry effect induced by the drift angle of the analyses is also noted.
Finally, Figures 18 to 22 show some perspective views, including pressure and velocity contours and streamlines, corresponding to the different analyzed cases.
 Figure 10. E0D0 8kn. Pressure contours on bulb. 
 Figure 11. E15D2 8kn. Pressure contours on bulb. 
 Figure 12. E15D4 8kn. Pressure contours on bulb. 
 Figure 13. E25D2 8kn. Pressure contours on bulb. 
 Figure 14. E0D0 9kn Wave map. 
 Figure 15. E15D2 9kn Wave map. 
 Figure 16. E15D4 9kn Wave map. 
 Figure 17. E25D2 9kn Wave map. 
 Figure 18. E15D2 7.5 kn. Pressure map on appendages and streamlines. Perspective view. 
 Figure 19. E15D4 7.5 kn. Velocity modulus contours. Perspective view. 
 Figure 20. E0D0 10kn. Pressure map on free surface and appendages. Perspective view.  Resistance results are also presented in Figures 22 to 25. Values obtained in the simulation have been compared to those resulting of a standard extrapolation technique using experimental data. Resistance has been computed from forces data evaluated in the global reference axes by R = Fx·cos() + Fy·sin(), where is the drift angle of the case. It is observed that the calculated results agree well with the experimental data, both in the total values and in the shape of the curve.
It is noted that the extrapolated values include an uncertainty factor, and the exact differences of the simulation results with real values is unknown.
 Figure 21. E0D0. Resistance graph. Comparison with extrapolation from experimental data. 
 Figure 22. E15D2. Resistance graph. Comparison with extrapolation from experimental data. 
 Figure 23. E15D4. Resistance graph. Comparison with extrapolation from experimental data. 
 Figure 24. E25D4. Resistance graph. Comparison with extrapolation from experimental data. 
7. CONCLUSIONS An implicit secondorder accurate monolithic scheme, based on the FIC formulation has been presented to solve incompressible free surface flow problems. The scheme presented is especially adequate for analysis of naval problems including dynamic sinkage and trim effects. In the present work the method has been applied in a validation study of America’s Cup boat. Results obtained in the numerical study show a great agreement with experimental data available in every condition. An academic version of the software Tdyn, using the formulation presented here, can be freely used and downloaded from http://www.compassis.com. Acknowledgements The authors thank the Spanish Ministry of Education and Science represented by the National Institute of Aerospace Technique (INTA) for permitting the publication of the towing tank tests of Rioja de España, and the model basin El Pardo (CEHIPAR) for sending the full documentation of the tests. Special gratefulness to IZAR shipbuilders for allowing the publication of the Rioja de España hull forms. The authors also thank Dr. Ramon Ribó and Mr. Enrique Escolano for many useful discussions and for helping in the development of CAECAM communication algorithms.
References 1 J. García, R. Ribó y M. López (2001), “A finite element method for hydrodynamic testing of ship forms”. 1^{st} International Congress on Maritime Transport. Barcelona (Spain). 2 C. Yang, R. Löhner and O. Soto (2001), “Optimization of a Wave Cancellation Multihull Ship using CFD Tools”; Proc. 8th Int. Symp. on Practical Design of Ships and Other Floating Structures  PRADS 2001, Shangai (China). 3 E. Oñate and J. García (2001), “Finite Element Analysis of Incompressible Flows with Free Surface Waves using a Finite Calculus Formulation”. ECCOMAS 2001. Swansea (UK). 4 E. Oñate (2000), “A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation”. Comp. Meth. Appl. Mech. Engng. 182, 12, 355370. 5 E. Oñate and J. García (1999), “A methodology for analysis of ﬂuidstructure interaction accounting for free surface waves”. European Conference on Computational Mechanics (ECCM99). Munich, (Germany). 6 O. Soto, R. Löhner, J. Cebral and R. Codina (2001), “A Time Accurate ImplicitMonolithic Finite Element Scheme for Incompressible Flow”. ECCOMAS 2001. Swansea UK. 7 E. Oñate, J. García and S. Idelsohn (1998), “An alpha adaptive approach for stabilized finite element solution of advectivediffusive problems with sharp gradients”, New Adv. In adaptive Comp. Met. In Mech., P. Ladeveze and J.T. Oden (Eds.), Elsevier. 8 E. Oñate, J García (2001), “A finite element method for fluidstructure interaction with surface waves using a finite calculus formulation”, Comput. Meth. Appl. Mech. Engng., Vol. 191, 635660. 9 E. Oñate, J. García and S. Idelsohn (1997), “Computation of the stabilization parameter for the finite element solution of advectivediffusive problems with sharp gradients”, Int. J. Num. Meth. Fluids, Vol. 25, pp. 13851407. 10 E. Oñate and J. García. “A stabilized finite element method for analysis of fluid structure interaction problems involving free surface waves”. Proceedings of Fluid Structure Interaction Conference, Trondheim (Norway) 1999. 11 E. Oñate, C. Sacco, S. Idelsohn (2000), “A finite point method for incompressible flow problems”, Comp. Visual. Sci. 2, 6775. 12 R. Löhner, C. Sacco, E. Onate and S. Idelsohn (2002), “A Finite Point Method for Compressible Flow”; Int. J. Num. Meth. Eng. 53, 17651779. 13 R. Codina (2001), “Pressure stability in fractional step finite element methods for incompressible flows” Jnl. Comp. Phys., Vol. 170, 112140. 14 Tdyn Theoretical Background and Reference Manual. Available to download at http:// www.compassis.com. 15 GiD. The personal pre/postprocessor. User Manual. Available to download at http:// www.gid.cimne.upc.es. 16 G. Farin (1993), Curves and surfaces for CADG. Academic Press Inc., third edition. 17 L. Larsson, and R. E. Eliasson (1994), Principles of Yacht Design. International Marine. Camden, Maine. 18 Gutelle, P. (1994), Design of Sailing Yachts. Warsah Publising, Southampton. ^{1 }^{Compass Ingeniería y Sistemas S.A. C/ Manuel Girona 61, bajos. 08034 Barcelona, }^{Spain}^{, }^{web page: http://www.compassis.com} ^{2 }^{Universidad Austral de Chile. Casilla 567. Valdivia, Chile, web page: http://www..uach.cl} ^{3 Nautatec, Ingeniería Naval Diseño de Yates S.L. C\ Arquitecto Gaudí,11 Bajo Exterior. }^{28016 Madrid, Spain, web page: http://www.nautatec.com} ^{4 International Center for Numerical Methods in Engineering (CIMNE). }^{C/ Gran Capitan s/n. Campus Norte UPC. 08034 } ^{Barcelona, Spain, web page: http://www.cimne.com}
