A finite element method for hydrodynamic testing of ship forms

Introduction

The prediction of the resistance and wave pattern, joint to the study of the flow around a ship are of major relevance in naval architecture. The analytical and numerical solutions of this problem has challenged mathematicians and hydrodynamicists for over a century.

In recent years, the advent of advanced numerical schemes for the Navier Stokes (NS) equations and the rapid development of computer hardware has enabled a realistic prediction of this problems.

The troubles in accurately solving this problem are mainly due to the difficulty of solving numerically the incompressible fluid dynamic equations, which include significant non linearities and the obstacles in solving the constraint equation stating that at the free surface boundary the fluid particles remain on that surface which position is in turn unknown.

Among the schemes developed over the last decade for the solution of the incompressible NS equations the fractional step schemes [1,2,3,4] yield highly accurate, pressure-stable [2] results by integrating in an explicit manner the advective terms of the NS equations. However, in most of the cases of interest for the naval architecture [5] the use of these techniques requires tens of thousands of time steps per simulation, rendering the schemes impractical. Most of the artificial compressibility and preconditioned schemes suffer from the same shortcoming.

On the other hand, the monolithic schemes treat, in general, the advective term in an implicit manner, which avoids the mentioned disadvantages. Nevertheless, these methods are very expensive from a computational point of view: the velocity and pressure discrete equations are coupled.

This paper presents advances in recent work of the authors to derive a fractional step scheme based on the stabilized finite element method that allows overcoming the above mentioned problem, resulting in a efficient time accurate scheme.

The starting point is the modified governing differential equations for the incompressible turbulent viscous flow and the free surface condition incorporating the necessary stabilization terms via a finite calculus (FIC) procedure developed by the authors [6,7,8,9]. This technique is based on writting the different balance equations over a domain of finite size and retaining higher order terms. These terms incorporate the ingredients for the necessary stabilization of any transient and steady state numerical solution already at the differential equations level.

The resulting stabilized equations are integrated in space using the standard finite element method, and in time using an implicit and uncoupled second order fractional step method, based on the scheme originally proposed by Soto [10].

FIC formulation

We consider the motion around a body of a viscous incompressible fluid including a free surface. The stabilized finite calculus (FIC) form of the governing differential equations [11] for the three dimensional (3D) problem can be written as

$r_{m_{i}}{\underline {-{\frac {1}{2}}h_{j}{\frac {\partial r_{m_{i}}}{\partial x_{j}}}}}=0{\begin{array}{cccc}&on&\Omega &i,j=1,2,3\end{array}}$

(1a)

$r_{d}{\underline {-{\frac {1}{2}}h_{j}{\frac {\partial r_{d}}{\partial x_{j}}}}}=0{\begin{array}{cccc}&on&\Omega &j=1,2,3\end{array}}$

(1b)

$r_{\beta }-{\underline {{\frac {1}{2}}h_{{\beta }_{j}}{\frac {\partial r_{\beta }}{\partial x_{j}}}}}=0{\begin{array}{cccc}&on&\Omega &j=1,2\end{array}}$

(1c)

where

${\begin{array}{c}r_{m_{i}}={\frac {\partial u_{i}}{\partial t}}+{\frac {\partial }{\partial x_{j}}}\left(u_{i}u_{j}\right)+{\frac {\partial p}{\partial x_{i}}}-{\frac {\partial {\tau }_{ij}}{\partial x_{j}}}\\r_{d}={\frac {\partial u_{i}}{\partial x_{i}}}{\begin{array}{cc}&i=1,2,3\end{array}}\\r_{\beta }={\frac {\partial \beta }{\partial t}}+u_{i}{\frac {\partial \beta }{\partial x_{i}}}{\begin{array}{cc}&i=1,2\end{array}},3\end{array}}$

In above u_i is the velocity along the i-th global reference axis, p is the dynamic pressure, β is the wave elevation and τ_ij are the viscous stress tensor components.

Let n_i be the unit outward normal to the boundary Ω and denoting by an overbar prescribed values, the boundary conditions for the stabilized problem are:

${\begin{array}{c}u={\overline {u}}{\begin{array}{ccc}&on&{\Gamma }_{u}\end{array}},\\{\begin{array}{ccc}p={\overline {p}}&and&\end{array}}n_{j}{\tau }_{ij}{\underline {-{\frac {1}{2}}h_{j}n_{j}r_{m_{i}}}}={\overline {t}}_{i}{\begin{array}{ccc}&on&{\Gamma }_{p}\end{array}}\\u_{j}n_{j}={\overline {u}}_{n}{\begin{array}{cc},&n_{j}{\tau }_{ij}g_{i}{\underline {-{\frac {1}{2}}h_{j}n_{j}r_{m_{i}}g_{i}}}={\overline {t}}_{1}\end{array}}{\begin{array}{ccc}{\begin{array}{cc}{\begin{array}{cc}&and\end{array}}&n_{j}{\tau }_{ij}s_{i}{\underline {-{\frac {1}{2}}h_{j}n_{j}r_{m_{i}}s_{i}}}={\overline {t}}_{2}\end{array}}&on&{\Gamma }_{\tau }\end{array}}\end{array}}$

(2)

for t ∈ (t₀,t_f). The boundary Γ has been considered split into three sets of disjoint components Γ_u, Γ_p and Γ_τ, the latter 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).

The underlined terms in eqs. (1)-(2) introduce the necessary stabilization for the numerical solution. Additional time stabilization terms can be accounted for in eqs. (1)-(2) [6,7] although they have been found unnecessary for the type of problem solved here.

The characteristic length distances h_j represent the dimensions of the finite domain where balance of momentum and mass is enforced. The characteristic distances h_βj in eq. (1c) represent the 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 [6,12,13].

Eqs. (1)-(2) are the starting point for deriving a variety of stabilized numerical methods for solving the incompressible Navier-Stokes equations with a free surface. It can be shown that a number of stabilized finite element methods and new meshless methods allowing equal order interpolations for the velocity and pressure fields can be from the modified form of the momentum and mass balance equations given above [6-9].

Implicit Fractional Step Formulation

Let us discretize in time the stabilized momentum equation (1a) using the trapezoidal rule (or θ method) as

${\frac {u_{i}^{n+1}-u_{i}^{n}}{\delta t}}+{\frac {\partial }{\partial x_{j}}}\left(u_{i}^{n+\theta }u_{j}^{n+\theta }\right)+$${\frac {\partial p^{n+\theta }}{\partial x_{i}}}-{\frac {\partial {\tau }_{ij}^{n+\theta }}{\partial x_{j}}}-$${\frac {1}{2}}h_{j}{\frac {\partial r_{m_{i}}^{n+\theta }}{\partial x_{j}}}=0$

(3)

the 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 the standard Crank Nicholson scheme, which is second order accurate in time 0(t²).

A classical 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:

${\frac {{\overline {u}}_{i}^{n+1}-u_{i}^{n}}{\delta t}}+{\frac {\partial }{\partial x_{j}}}\left({\overline {u}}_{i}^{n+\theta }{\overline {u}}_{j}^{n+\theta }\right)+{\frac {\partial p^{n}}{\partial x_{i}}}-{\frac {\partial {\overline {\tau }}_{ij}^{n+\theta }}{\partial x_{j}}}-{\frac {1}{2}}h_{j}{\frac {\partial {\overline {r}}_{m_{i}}^{n+\theta }}{\partial x_{j}}}=0$

(4a)

$\delta t{\frac {{\partial }^{2}}{\partial x_{i}\partial x_{i}}}\left(p^{n+1}-p^{n}\right)={\frac {\partial u_{i}^{\ast }}{\partial x_{i}}}+{\frac {1}{2}}h_{j}{\frac {\partial r_{d}}{\partial x_{j}}}$

(4b)

$u_{i}^{n+1}={\overline {u}}_{i}-\delta t{\frac {\partial }{\partial x_{i}}}\left(p^{n+1}-p^{n}\right)$

(4c)

where the terms denoted by an overbar are those calculated with the intermediate velocity ū_i, which is introduced to allow the momentum splitting.

The error due to taking implicit advective and viscous terms in (4a) can be shown [10] to be of the same order than the error of the stabilizing term and therefore, it has the same order of approximation than the original time discretization (3).

Monolithic Stabilized Scheme

At this point, it is important to introduce the associated matrix structure corresponding to the variational discrete form of (4) (see [6] for details of this derivation):

$M{\frac {1}{\delta t}}\left({\overline {U}}^{n+1}-U^{n}\right)+K\left({\overline {U}}^{n+\theta }\right){\overline {U}}^{n+\theta }-GP^{n}=0$

(5a)

$\left(\delta t+\tau \right)L\left(P^{n+1}-P^{n}\right)-HU^{n}+D{\overline {U}}^{n+1}=0$

(5b)

${\frac {1}{\delta t}}\left(U^{n+1}-{\overline {U}}^{n+1}\right)-\left(P^{n+1}-P^{n}\right)=0$

(5c)

By taking Ūⁿ⁺¹ from (5c) and inserting the result in (5a)-(5b), the following system of equations is obtained:

$M{\frac {1}{\delta t}}\left(U^{n+1}-U^{n}\right)+K\left(U^{n+\theta }\right)U^{n+\theta }-GP^{n}=E\left(U^{n+\theta }\right)$

(6a)

$\delta t\left(L-DM^{-1}G\right)\left(P^{n+1}-P^{n}\right)-H\left(U^{n},P^{n}\right)=-DU^{n+1}$

(6b)

The term E(U^n+θ) in (6a) is the error coming from the implicit treatment of the advective and viscous terms, which is of order 0(t²). However, such term can be eliminated as in (6a)-(6b) by writing the following analog monolithic scheme:

$M{\frac {1}{\delta t}}\left(U^{n+1,i}-U^{n}\right)+K\left(U^{n+\theta ,i}\right)U^{n+\theta ,i}-GP^{n+1,i-1}=0$

(7a)

$\delta t\left(L-DM^{-1}G\right)\left(P^{n+1,i}-P^{n}\right)-H\left(U^{n+1,i},P^{n+1,i}\right)=-DU^{n+1,i}$

(7b)

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 wave elevation coupling effect is included in the monolithic scheme, obtaining the following system of equations:

$M{\frac {1}{\delta t}}\left(U^{n+1,i}-U^{n}\right)+K\left(U^{n+\theta ,i}\right)U^{n+\theta ,i}-GP^{n+1,i-1}=0$

(8a)

$\delta t\left(L-DM^{-1}G\right)\left(P^{n+1,i}-P^{n}\right)-H\left(U^{n+1,i},P^{n+1,i}\right)=-DU^{n+1,i}$

(8b)

$M{\frac {1}{\delta t}}\left(B^{n+1,i}-B^{n}\right)+A\left(U^{n+\theta ,i}\right)B^{n+\theta ,i}-R(U^{n+\theta ,i},B^{n})=MU_{3}^{n+\theta ,i}$

(8b)

where the U_i^n+θ is the vector of the velocity along the i-th global reference axis.

Figure shows the flow chart of the solution algorithm of eqs. (8).

Figure 1. Flow chart of the solution algorithm of eqs. (8).

Application Examples

The algorithm here presented have been implemented within the CFD environment Tdyn [14,15]. This software has been used to solve the examples shown next.

Example 1. DTMB 5415 model

The first case here analysed is a benchmark test case, the David Taylor Model Basin (DTMB) 5415 model. The geometry used in the analysis was obtained from the Gothemburg 2000 Workshop database [16] and is based on the NURBS definition shown in Figure 2.

Figure 2. DTMB 5415. Hull geometry based on the NURBS entities.

Figure 3. DTMB 5415. Surface mesh on the hull.

The obtained results are compared with the experimental data available in the EFD database developed by the Iowa Institute of Hydraulic Research, the Instituto Nazionale per Studi ed Esperienze di Architettura Navale and the Naval Surface Warfare Center [16].

Figure 4. DTMB 5415. Wave Profile on hull.

The main characteristics of this analysis are listed below:

Length:	5.72 m	Beam:	0.5 m
Draught:	0.248 m	Wetted Surface:	4.861m2
Velocity:	2.1 m/s	Froude Number:	0.28
Viscosity:	0.001 Kg/ms	Density:	1000 Kg/m3
Reynolds number:	12.3 106

The analysis was carried out for three different grids (from 150.000 to 600.000 linear tetrahedra, corresponding to 25.000 and 115.000 nodes) in order to qualitatively analyse the influence of the element size in the solution. Here, only the results corresponding to the finest grid are shown. The smallest element size used was 0.002 m and the maximum 0.750 m. The surface mesh of the DTMB 5415 used in the last analysis is shown in Figure 3. The Smagorinsky turbulence model with the extended law of the wall was used.

Figures 4 and 5 show the wave profile on the hull and in a cut at y/L =0.082, respectively. Numerical results obtained with the numerical model are compared to the experimental data.

Figure 5. DTMB 5415. Wave Profile y/L = 0.082.

Example 2. KVLCC2 Model

The next example is another benchmark test case, the KVLCC2 model. The geometry used in the example, presented in Figure 6, is based in a NURBS surface definition obtained from the Hydrodynamic Performance Research team of Korea (KRISO) 17. The obtained results are compared with the experimental data available in the KRISO database [19].

Figure 6. KVLCC2. Hull geometry based on the NURBS entities.

Figure 7. KVLCC2. Surface mesh on the hull.

Several simulations of this model were carried out, in order to complete the comparison with the experimental data available. The mesh used in all cases has 380.000 linear tetrahedra elements and 70.000 nodes. The smallest element size used was 0.001 m and the maximum was 0.50 m. The surface mesh of the KVLCC2 used is shown in Figure 7.

Figure 8. KVLCC2. Wave Profile on hull.

Test 1.- Wave pattern calculation. The main characteristics of this analysis are listed below:

Length:	5.52 m	Beam:	0.82 m
Draught:	0.18 m	Wetted Surface:	8.08 m2
Velocity:	1.05 m/s	Froude Number:	0.142
Viscosity:	0.00126 Kg/ms	Density:	1000 Kg/m3
Reynolds number:	4.63 106

The turbulence model used in this case was the K model. Figures 8 and 9 show the wave profiles on the hull and in a cut at y/L = 0.082 obtained in test 1, compared to the experimental data. The obtained results are quantitatively good close to the hull, but a lost of accuracy can be observed in the profiles further from the hull. This is probably due to the fact that the element sizes are not small enough in this area, for an accurate enough capture of the fluid perturbation.

Figure 9. KVLCC2. Wave Profile y/L = 0.0964.

Test 2.- Wake analysis at different planes. Several turbulence models were used in this phase (Smagorinsky, K and K − ε model) in order to verify the quality of the results. Here, only the results corresponding to the last case are shown. However, it is important to note that the velocity maps obtained even for the simplest Smagorinsky model were qualitatively good, showing the accuracy of the fluid flow solver scheme used.

Figure 10. Contours of the X component of the velocity on a plane at 2.71 m from the aft perpendicular. Comparison with experimental data (right).

Figure 11. Contours of the X component of the velocity on a plane at 2.82 m from the aft perpendicular. Comparison with experimental data (right).

The main characteristics of this analysis are listed below:

Length:	2.76 m	Beam:	0.41 m
Draught:	0.09 m	Wetted Surface:	2.02 m2
Velocity:	25 m/s	Froude Number:	0.0
Viscosity:	3.05 10-5 Kg/ms	Density:	1.01 Kg/m3
Reynolds number:	4.63 106

Figure 10 shows the contours of the X component of the velocity on a plane at 2.71 m from the aft perpendicular compared with the experimental data available.

Figure 11 shows the map of the X component of the velocity on a plane at 2.82 m from the aft perpendicular compared with the experimental data available.

Example 3. 79’ Catamaran

The application example here presented is the parametric analysis of a 79´ catamaran designed by the company NAUTATEC. The main characteristics of this boat are listed next (see Figure 12):

Length overall	23.95 m
Length between perpendiculars	23.35 m
Total breadth	11.75 m
Extreme breadth of every hull	1.65 m
Draught with appendages	1.31 m
Hull depth	2.40 m
Displacement	26 730 Kg
Maximum sail area	222 m2
Passengers	120

This study started from a fixed well-known geometry of the hulls, being the object of the study, the effect of the distance between hulls in the resistance and lateral forces.

Three different alternatives of distance between hulls studied were:

Case		Name
Distance between hulls centerlines	10.1 m	79cata5_0
Distance between hulls centerlines	9.1 m	79cata4_5
Distance between hulls centerlines	8.1 m	79cata4_0

Figure 12. NURBS based geometrical definition of the 79´catamaran. Different alternatives studied.

The analyses were carried out for a range of ship speed from 8 kn to 28 kn, being the real range from 8 kn to 16 kn. Regarding the drift angles, four cases were considered: 0.0º, 1.0º, 2.5º and 5.0º, being the project range from 0.0º to 2.5º.

The grid used for the analyses consisted of 450.000 linear tetrahedra. The fluid viscosity taken was 1.3 Kg/ms and the density 1024 Kg/m³. The k turbulence model with the extended law of the wall was used.

The basic results of these analyses are presented next:

Figure 13. Resistance and lift curves for the 79cata4_5 case.

Figure 14. Resistance curves for the different cases analyzed.

The results shown in Figure 14 allowed to select as optimal configuration, regarding minimum resistance, the 79cata5_0 geometry, which was finally built.

Conclusions

An implicit second-order accurate monolithic scheme, based on the FIC formulation was presented to solve incompressible free surface flow problems. The final system of equations resulting from the time and space discretization is solved in each time step in an uncoupled manner.

The numerical experience indicates that the formulation is very efficient for free surfaces flows, when the critical time step of the problem is some orders of magnitude smaller than the time step required to obtain time-accurate results (physical time step), and that its time accuracy is excellent.

Numerical results obtained in the 3D viscous analysis of complex ship geometries indicate that the proposed numerical method can be used for confidence for practical hydrodynamic design purposes in naval architecture.

Acknowledgements

The authors want to thank Dr. O. Soto for many useful discussions.

References

[1] Chorin. On the convergence of discrete approximation to the Navier-Stokes equations. Math. Comput., 23, 1969.

[2] R. Codina. Pressure stability in fractional step finite element methods for incompressible flows. J. Comp. Phys., Submitted for publication, 2000.

[3] J. García, E. Oñate, H. Sierra, C. Sacco y S. Idelsohn. A stabilised numerical method for analysis of ship hydrodynamics. ECCOMAS 98 (Vol. II). Atenas 1998.

[4] R. Löhner, C. Yang, E. Oñate, and S. Idelsohn. An unstructured grid-based, parallel

free surface solver. AIAA-97-1830, 1997.

[5] Löhner, C. Yang, E. Oñate and S. Idelsohn, An unstructured grid-based parallel free surface solver, Appl. Num. Math., 31, 271--293, 1999.

[6] E. Oñate and J. García. Finite Element Analysis of Incompressible Flows with Free Surface Waves using a Finite Calculus Formulation. ECCOMAS 2001. Swansea UK 2001.

[7] E. Oñate. A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation. Comp. Meth. Appl. Mech. Engng. 182, 1-2, 355-370, 2000.

[8] E. Oñate and J. García. A methodology for analysis of fluid-structure interaction accounting for free surface waves. European Conference on Computational Mechanics (ECCM99). Munich, Germany, September 1999.

[9] 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 (1999).

[10] O. Soto, R. Löhner, J. Cebral and R. Codina. A Time Accurate Implicit-Monolithic Finite Element Scheme for Incompressible Flow. ECCOMAS 2001. Swansea UK 2001.

[11] E. Oñate, J. García and S. Idelsohn, An alpha adaptive approach for stabilized finite element solution of advective-diffusive problems with sharp gradients, New Adv. In adaptive Comp. Met. In Mech., P. Ladeveze and J.T. Oden (Eds.), Elsevier (1998).

[12] E. Oñate, Derivation of stabilized equations for advective-diffusive transport and fluid flow problems, Comput. Meth. Appl. Mech. Engng., Vol. 151, 1-2, pp. 233-267 (1998).

[13] E. Oñate, J. García and S. Idelsohn, Computation of the stabilization parameter for the finite element solution of advective-diffusive problems with sharp gradients, Int. J. Num. Meth. Fluids, Vol. 25, pp. 1385-1407 (1997).

[14] Tdyn: Theoretical Background. Available at http://www.compassis.com. 2001.

[15] GiD. The personal pre/postprocessor. User Manual. Available at [http:// http://]www.gid.cimne.upc.es. 2001.

[16] David Taylor Model Basin 5415 Model Database. http://www.iihr.uiowa.edu/gothenburg2000/5415/combatant.html.

[17] Korea Research Institute of Ships and Ocean Engineering (KRISO). http://www.iihr.uiowa.edu/gothenburg2000/KVLCC/tanker.html.