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.

<span id='_Ref526828097'></span>Among the schemes developed over the last decade for the solution of the incompressible NS equations the fractional step schemes [<span id='cite-1'></span>[[#1|1]],<span id='cite-2'></span>[[#2|2]],<span id='cite-3'></span>[[#3|3]],<span id='cite-4'></span>[[#4|4]]] yield highly accurate, pressure-stable [<span id='cite-_Ref526828097'></span>[[#_Ref526828097|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 [<span id='cite-5'></span>[[#5|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.

<span id='_Ref526828328'></span><span id='_Ref526828342'></span><span id='_Ref526828461'></span>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 [<span id='cite-6'></span>[[#6|6]],<span id='cite-7'></span>[[#7|7]],<span id='cite-8'></span>[[#8|8]],<span id='cite-9'></span>[[#9|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.

<span id='_Ref526828500'></span>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 [<span id='cite-10'></span>[[#10|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 [<span id='cite-11'></span>[[#11|11]]] for the three dimensional (3D) problem can be written as

{| class="formulaSCP" style="width: 100%; text-align: left;"
|-
|
{| style="text-align: left; margin:auto;width: 100%"
|-
| style="text-align: center;" | <math>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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1a)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;"
|-
|
{| style="text-align: left; margin:auto;width: 100%"
|-
| style="text-align: center;" | <math>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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1b)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;"
|-
|
{| style="text-align: left; margin:auto;width: 100%"
|-
| style="text-align: center;" | <math>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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1c)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;"
|-
|
{| style="text-align: left; margin:auto;width: 100%"
|-
| style="text-align: center;" | <math>\begin{array}{c}
r_{m_i}=\frac{\partial u_i}{\partial t}+\frac{\partial }{\partial x_j}\left(u_iu_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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | 
|}

In above ''u<sub>i</sub>'' is the velocity along the i-th global reference axis, ''p'' is the dynamic pressure, ''&#x03b2;'' is the wave elevation and ''&#x03c4;<sub>ij</sub>'' are the viscous stress tensor components.

Let ''n<sub>i</sub>'' be the unit outward normal to the boundary Ω and denoting by an overbar prescribed values, the boundary conditions for the stabilized problem are:

{| class="formulaSCP" style="width: 100%; text-align: left;"
|-
|
{| style="text-align: left; margin:auto;width: 100%"
|-
| style="text-align: center;" | <math>\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_jn_jr_{m_i}}={\overline{t}}_i\begin{array}{ccc}
 & on & {\Gamma }_p
\end{array}\\
u_jn_j={\overline{u}}_n\begin{array}{cc}
, & n_j{\tau }_{ij}g_i\underline{-\frac{1}{2}h_jn_jr_{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_jn_jr_{m_i}s_i}={\overline{t}}_2
\end{array} & on & {\Gamma }_{\tau }
\end{array}
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

for ''t ''&#x2208; (''t<sub>0</sub>,t<sub>f</sub>''). The boundary Γ has been considered split into three sets of disjoint components Γ''<sub>u</sub>'', Γ''<sub>p</sub>'' and Γ''<sub>&#x03c4;</sub>'', the latter being the part where mixed conditions are prescribed.

Vectors ''g<sub>i</sub>'' and ''s<sub>i</sub>'' span the space tangent to Γ''<sub>&#x03c4;</sub>''. Finally, Γ''<sub>u</sub>'' and Γ''<sub>p</sub>'' 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) [<span id='cite-_Ref526828328'></span>[[#_Ref526828328|6]],<span id='cite-_Ref526828342'></span>[[#_Ref526828342|7]]] although they have been found unnecessary for the type of problem solved here.

The ''characteristic length'' distances ''h<sub>j</sub>'' represent the dimensions of the finite domain where balance of momentum and mass is enforced. The characteristic distances ''h<sub>&#x03b2;j</sub>'' 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 [<span id='cite-_Ref526828328'></span>[[#_Ref526828328|6]],<span id='cite-12'></span>[[#12|12]],<span id='cite-13'></span>[[#13|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 [<span id='cite-_Ref52682'></span>[[#_Ref52682|6]]-<span id='cite-_Ref526828461'></span>[[#_Ref526828461|9]]].

==Implicit Fractional Step Formulation==

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

{| class="formulaSCP" style="width: 100%; text-align: left;"
|-
|
{| style="text-align: left; margin:auto;width: 100%"
|-
| style="text-align: center;" | <math>\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)+</math><math>\frac{\partial p^{n+\theta }}{\partial x_i}-\frac{\partial {\tau }_{ij}^{n+\theta }}{\partial x_j}-</math><math>\frac{1}{2}h_j\frac{\partial r_{m_i}^{n+\theta }}{\partial x_j}=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (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<sup>2</sup>).

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:

{| class="formulaSCP" style="width: 100%; text-align: left;"
|-
|
{| style="text-align: left; margin:auto;width: 100%"
|-
| style="text-align: center;" | <math>\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</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4a)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;"
|-
|
{| style="text-align: left; margin:auto;width: 100%"
|-
| style="text-align: center;" | <math>\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4b)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;"
|-
|
{| style="text-align: left; margin:auto;width: 100%"
|-
| style="text-align: center;" | <math>u_i^{n+1}={\overline{u}}_i-\delta t\frac{\partial }{\partial x_i}\left(p^{n+1}-p^n\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4c)
|}

where the terms denoted by an overbar are those calculated with the intermediate velocity ''ū<sub>i</sub>'', which is introduced to allow the momentum splitting.

The error due to taking implicit advective and viscous terms in (4a) can be shown [<span id='cite-_Ref526828500'></span>[[#_Ref526828500|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 [<span id='cite-_R'></span>[[#_R|6]]] for details of this derivation):

{| class="formulaSCP" style="width: 100%; text-align: left;"
|-
|
{| style="text-align: left; margin:auto;width: 100%"
|-
| style="text-align: center;" | <math>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</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5a)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;"
|-
|
{| style="text-align: left; margin:auto;width: 100%"
|-
| style="text-align: center;" | <math>\left(\delta t+\tau \right)L\left(P^{n+1}-P^n\right)-HU^n+D{\overline{U}}^{n+1}=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5b)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;"
|-
|
{| style="text-align: left; margin:auto;width: 100%"
|-
| style="text-align: center;" | <math>\frac{1}{\delta t}\left(U^{n+1}-{\overline{U}}^{n+1}\right)-\left(P^{n+1}-P^n\right)=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5c)
|}

By taking ''Ū<sup>n+1</sup>'' from (5c) and inserting the result in (5a)-(5b), the following system of equations is obtained:

{| class="formulaSCP" style="width: 100%; text-align: left;"
|-
|
{| style="text-align: left; margin:auto;width: 100%"
|-
| style="text-align: center;" | <math>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)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6a)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;"
|-
|
{| style="text-align: left; margin:auto;width: 100%"
|-
| style="text-align: center;" | <math>\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6b)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: left;"
|-
|
{| style="text-align: left; margin:auto;width: 100%"
|-
| style="text-align: center;" | <math>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</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7a)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;"
|-
|
{| style="text-align: left; margin:auto;width: 100%"
|-
| style="text-align: center;" | <math>\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (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:

{| class="formulaSCP" style="width: 100%; text-align: left;"
|-
|
{| style="text-align: left; margin:auto;width: 100%"
|-
| style="text-align: center;" | <math>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</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8a)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;"
|-
|
{| style="text-align: left; margin:auto;width: 100%"
|-
| style="text-align: center;" | <math>\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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8b)
|}
{| class="formulaSCP" style="width: 100%; text-align: left;"
|-
|
{| style="text-align: left; margin:auto;width: 100%"
|-
| style="text-align: center;" | <math>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}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8b)
|}

where the ''U<sub>i</sub><sup>n+θ</sup>'' 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).

[[Image:draft_García-Espinosa_298267853-picture-Group 48.svg|center|600px]]

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">'''Figure 1. Flow chart of the solution algorithm of eqs. (8).'''</span></div>

==Application Examples==

The algorithm here presented have been implemented within the CFD environment Tdyn [<span id='cite-14'></span>[[#14|14]],<span id='cite-15'></span>[[#15|15]]]. This software has been used to solve the examples shown next.

===Example 1. DTMB 5415 model===

<span id='_Ref526947824'></span>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 [<span id='cite-16'></span>[[#16|16]]] and is based on the NURBS definition shown in <span id='cite-_Ref526945208'></span>[[#_Ref526945208|Figure 2]].

{| style="width: 100%;" 
|-
|  style="vertical-align: top;width: 41%;"| [[Image:draft_García-Espinosa_298267853-image15-c.png|600px]] 
|-
|  style="text-align: center;vertical-align: top;"|<span id='_Ref526945208'></span><span style="text-align: center; font-size: 75%;">'''Figure 2. DTMB 5415. Hull geometry based on the NURBS entities.'''</span>
|}


{| style="width: 100%;" 
|-
|  style="vertical-align: top;width: 100%;"| [[Image:draft_García-Espinosa_298267853-image16-c.png|600px]] 
|-
|  style="text-align: center;vertical-align: top;"|<span id='_Ref526945245'></span><span style="text-align: center; font-size: 75%;">'''Figure 3. DTMB 5415. Surface mesh on the hull.'''</span>
|}


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 [<span id='cite-_Ref526947824'></span>[[#_Ref526947824|16]]].

{| style="width: 100%;" 
|-
|  style="vertical-align: top;width: 100%;"| [[Image:draft_García-Espinosa_298267853-image17.png|600px]] 
|-
|  style="text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">'''Figure 4. DTMB 5415. Wave Profile on hull.'''</span>
|}


The main characteristics of this analysis are listed below:

{| style="width: 100%;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;vertical-align: top;"|Length: 
|  style="border: 1pt solid black;vertical-align: top;"|5.72 m
|  style="border: 1pt solid black;vertical-align: top;"|Beam: 
|  style="border: 1pt solid black;vertical-align: top;"|0.5 m
|-
|  style="border: 1pt solid black;vertical-align: top;"|Draught: 
|  style="border: 1pt solid black;vertical-align: top;"|0.248 m
|  style="border: 1pt solid black;vertical-align: top;"|Wetted Surface: 
|  style="border: 1pt solid black;vertical-align: top;"|4.861m<sup>2</sup> 
|-
|  style="border: 1pt solid black;vertical-align: top;"|Velocity: 
|  style="border: 1pt solid black;vertical-align: top;"|2.1 m/s
|  style="border: 1pt solid black;vertical-align: top;"|Froude Number:
|  style="border: 1pt solid black;vertical-align: top;"|0.28
|-
|  style="border: 1pt solid black;vertical-align: top;"|Viscosity: 
|  style="border: 1pt solid black;vertical-align: top;"|0.001 Kg/ms
|  style="border: 1pt solid black;vertical-align: top;"|Density: 
|  style="border: 1pt solid black;vertical-align: top;"|1000 Kg/m<sup>3</sup>
|-
|  style="border: 1pt solid black;vertical-align: top;"|Reynolds number: 
|  style="border: 1pt solid black;vertical-align: top;"|12.3 10<sup>6</sup>
|  style="border-top: 1pt solid black;border-left: 1pt solid black;vertical-align: top;"|
|  style="border-top: 1pt solid black;vertical-align: top;"|
|}


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 <span id='cite-_Ref526945245'></span>[[#_Ref526945245|Figure 3]]. The Smagorinsky turbulence model with the extended law of the wall was used.

Figures 4 and 5 show the wave proﬁle 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.

{| style="width: 100%;" 
|-
|  style="vertical-align: top;width: 100%;"|
[[Image:draft_García-Espinosa_298267853-image18.png|center|600px]]

|-
|  style="text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">'''Figure 5. DTMB 5415. Wave Profile y/L = 0.082.'''</span>
|}


===Example 2. KVLCC2 Model===

The next example is another benchmark test case, the KVLCC2 model. The geometry used in the example, presented in <span id='cite-_Ref526945909'></span>[[#_Ref526945909|Figure 6]], is based in a NURBS surface definition obtained from the Hydrodynamic Performance Research team of Korea (KRISO) <span id='cite-17'></span>[[#17|17]]. The obtained results are compared with the experimental data available in the KRISO database [19].

{| style="width: 100%;" 
|-
|  style="vertical-align: top;width: 41%;"| [[Image:draft_García-Espinosa_298267853-image19.png|600px]] 
|-
|  style="text-align: center;vertical-align: top;"|<span id='_Ref526945909'></span><span style="text-align: center; font-size: 75%;">'''Figure 6. KVLCC2. Hull geometry based on the NURBS entities.'''</span>
|}


{| style="width: 100%;" 
|-
|  style="vertical-align: top;width: 100%;"| [[Image:draft_García-Espinosa_298267853-image20-c.png|600px]] 
|-
|  style="text-align: center;vertical-align: top;"|<span id='_Ref526945928'></span><span style="text-align: center; font-size: 75%;">'''Figure 7. KVLCC2. Surface mesh on the hull.'''</span>
|}


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 <span id='cite-_Ref526945928'></span>[[#_Ref526945928|Figure 7]].

{| style="width: 100%;" 
|-
|  style="vertical-align: top;width: 100%;"|
[[Image:draft_García-Espinosa_298267853-image21.png|center|600px]]

|-
|  style="text-align: center;vertical-align: top;"|<span id='_Ref526946394'></span><span style="text-align: center; font-size: 75%;">'''Figure 8. KVLCC2. Wave Profile on hull.'''</span>
|}


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

{| style="width: 100%;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;vertical-align: top;"|Length: 
|  style="border: 1pt solid black;vertical-align: top;"|5.52 m
|  style="border: 1pt solid black;vertical-align: top;"|Beam: 
|  style="border: 1pt solid black;vertical-align: top;"|0.82 m
|-
|  style="border: 1pt solid black;vertical-align: top;"|Draught: 
|  style="border: 1pt solid black;vertical-align: top;"|0.18 m
|  style="border: 1pt solid black;vertical-align: top;"|Wetted Surface: 
|  style="border: 1pt solid black;vertical-align: top;"|8.08 m<sup>2</sup> 
|-
|  style="border: 1pt solid black;vertical-align: top;"|Velocity: 
|  style="border: 1pt solid black;vertical-align: top;"|1.05 m/s
|  style="border: 1pt solid black;vertical-align: top;"|Froude Number:
|  style="border: 1pt solid black;vertical-align: top;"|0.142
|-
|  style="border: 1pt solid black;vertical-align: top;"|Viscosity: 
|  style="border: 1pt solid black;vertical-align: top;"|0.00126 Kg/ms
|  style="border: 1pt solid black;vertical-align: top;"|Density: 
|  style="border: 1pt solid black;vertical-align: top;"|1000 Kg/m<sup>3</sup>
|-
|  style="border: 1pt solid black;vertical-align: top;"|Reynolds number: 
|  style="border: 1pt solid black;vertical-align: top;"|4.63 10<sup>6</sup>
|  style="border-top: 1pt solid black;border-left: 1pt solid black;vertical-align: top;"|
|  style="border-top: 1pt solid black;vertical-align: top;"|
|}


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.

{| style="width: 100%;" 
|-
|  style="vertical-align: top;width: 100%;"|
[[Image:draft_García-Espinosa_298267853-image22.png|center|600px]]

|-
|  style="text-align: center;vertical-align: top;"|<span id='_Ref526946381'></span><span style="text-align: center; font-size: 75%;">'''Figure 9. KVLCC2. Wave Profile y/L = 0.0964.'''</span>
|}


''Test 2.- ''Wake analysis at different planes. Several turbulence models were used in this phase (Smagorinsky, K and K − &#x03b5; 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.

{| style="width: 100%;" 
|-
|  style="vertical-align: top;width: 100%;"| [[Image:draft_García-Espinosa_298267853-image23-c.png|600px]] 
|-
|  style="vertical-align: top;"|<span id='_Ref526946547'></span><span style="text-align: center; font-size: 75%;">'''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).'''</span>
|}


{| style="width: 100%;" 
|-
|  style="vertical-align: top;width: 100%;"| [[Image:draft_García-Espinosa_298267853-image24-c.png|600px]] 
|-
|  style="vertical-align: top;"|<span id='_Ref526946706'></span><span style="text-align: center; font-size: 75%;">'''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).'''</span>
|}


The main characteristics of this analysis are listed below:

{| style="width: 100%;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;vertical-align: top;"|Length: 
|  style="border: 1pt solid black;vertical-align: top;"|2.76 m
|  style="border: 1pt solid black;vertical-align: top;"|Beam: 
|  style="border: 1pt solid black;vertical-align: top;"|0.41 m
|-
|  style="border: 1pt solid black;vertical-align: top;"|Draught: 
|  style="border: 1pt solid black;vertical-align: top;"|0.09 m
|  style="border: 1pt solid black;vertical-align: top;"|Wetted Surface: 
|  style="border: 1pt solid black;vertical-align: top;"|2.02 m<sup>2</sup> 
|-
|  style="border: 1pt solid black;vertical-align: top;"|Velocity: 
|  style="border: 1pt solid black;vertical-align: top;"|25 m/s
|  style="border: 1pt solid black;vertical-align: top;"|Froude Number:
|  style="border: 1pt solid black;vertical-align: top;"|0.0
|-
|  style="border: 1pt solid black;vertical-align: top;"|Viscosity: 
|  style="border: 1pt solid black;vertical-align: top;"|3.05 10<sup>-5</sup> Kg/ms
|  style="border: 1pt solid black;vertical-align: top;"|Density: 
|  style="border: 1pt solid black;vertical-align: top;"|1.01 Kg/m<sup>3</sup>
|-
|  style="border: 1pt solid black;vertical-align: top;"|Reynolds number: 
|  style="border: 1pt solid black;vertical-align: top;"|4.63 10<sup>6</sup>
|  style="border-top: 1pt solid black;border-left: 1pt solid black;vertical-align: top;"|
|  style="border-top: 1pt solid black;vertical-align: top;"|
|}


<span id='cite-_Ref526946547'></span>[[#_Ref526946547|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.

<span id='cite-_Ref526946706'></span>[[#_Ref526946706|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 <span id='cite-_Ref526827589'></span>[[#_Ref526827589|Figure 12]]):

{| style="width: 68%;" 
|-
|  style="vertical-align: top;"|Length overall 
|  style="text-align: right;vertical-align: top;"|23.95 m
|-
|  style="vertical-align: top;"|Length between perpendiculars
|  style="text-align: right;vertical-align: top;"|23.35 m
|-
|  style="vertical-align: top;"|Total breadth
|  style="text-align: right;vertical-align: top;"|11.75 m
|-
|  style="vertical-align: top;"|Extreme breadth of every hull
|  style="text-align: right;vertical-align: top;"|1.65 m
|-
|  style="vertical-align: top;"|Draught with appendages
|  style="text-align: right;vertical-align: top;"|1.31 m
|-
|  style="vertical-align: top;"|Hull depth
|  style="text-align: right;vertical-align: top;"|2.40 m
|-
|  style="vertical-align: top;"|Displacement
|  style="text-align: right;vertical-align: top;"|26 730 Kg
|-
|  style="vertical-align: top;"|Maximum sail area
|  style="text-align: right;vertical-align: top;"|222 m<sup>2</sup>
|-
|  style="vertical-align: top;"|Passengers
|  style="text-align: right;vertical-align: top;"|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:

{| style="width: 85%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|Case
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|Name
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|Distance between hulls centerlines
|  style="border: 1pt solid black;text-align: right;vertical-align: top;"|10.1 m
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|79cata5_0
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|Distance between hulls centerlines
|  style="border: 1pt solid black;text-align: right;vertical-align: top;"|9.1 m
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|79cata4_5
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|Distance between hulls centerlines
|  style="border: 1pt solid black;text-align: right;vertical-align: top;"|8.1 m
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|79cata4_0
|}


{| style="width: 100%;" 
|-
|  style="vertical-align: top;width: 31%;"|
[[Image:draft_García-Espinosa_298267853-image25.png|center|600px]]

|  style="vertical-align: top;width: 73%;"|
[[Image:draft_García-Espinosa_298267853-image26.png|center|600px]]

|-
|  colspan='2'  style="vertical-align: top;"|<span id='_Ref526827589'></span><span style="text-align: center; font-size: 75%;">'''Figure 12. NURBS based geometrical definition of the 79´catamaran. Different alternatives studied.'''</span>
|}


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<sup>3</sup>. The k turbulence model with the extended law of the wall was used.

The basic results of these analyses are presented next:


[[Image:draft_García-Espinosa_298267853-picture-Group 43.svg|center|600px]]

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">'''Figure 13. Resistance and lift curves for the 79cata4_5 case.'''</span></div>


[[Image:draft_García-Espinosa_298267853-picture-Group 17.svg|center|600px]]

<div id="_Ref526947001" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">'''Figure 14. Resistance curves for the different cases analyzed.'''</span></div>

The results shown in <span id='cite-_Ref526947001'></span>[[#_Ref526947001|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==

<div id="1"></div>
[[#cite-1|[1]]] Chorin. ''On the convergence of discrete approximation to the Navier-Stokes equations''. Math. Comput., 23, 1969.

<div id="2"></div>
[[#cite-2|[2]]] R. Codina. ''Pressure stability in fractional step finite element methods for incompressible flows''. J. Comp. Phys., Submitted for publication, 2000.

<div id="3"></div>
[[#cite-3|[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.

<div id="4"></div>
[[#cite-4|[4]]] R. Löhner, C. Yang, E. Oñate, and S. Idelsohn. ''An unstructured grid-based, parallel''

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

<div id="5"></div>
[[#cite-5|[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.

<div id="6"></div>
[[#cite-6|[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.

<div id="7"></div>
[[#cite-7|[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.

<div id="8"></div>
[[#cite-8|[8]]] E. Oñate and J. García. ''A methodology for analysis of ﬂuid-structure interaction accounting for free surface waves''. European Conference on Computational Mechanics (ECCM99). Munich, Germany, September 1999.

<div id="9"></div>
[[#cite-9|[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).

<div id="10"></div>
[[#cite-10|[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.

<div id="11"></div>
[[#cite-11|[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).

<div id="12"></div>
[[#cite-12|[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).

<div id="13"></div>
[[#cite-13|[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).

<div id="14"></div>
[[#cite-14|[14]]] ''Tdyn: Theoretical Background''. Available at [http://www http://www].compassis.com. 2001.

<div id="15"></div>
[[#cite-15|[15]]] ''GiD. The personal pre/postprocessor''. User Manual. Available at  [http:// http://][http://www.gid.cimne.upc.es www.gid.cimne.upc.es]. 2001.

<div id="16"></div>
[[#cite-16|[16]]] ''David Taylor Model Basin 5415 Model Database''. [http://www.iihr.uiowa.edu/gothenburg2000/5415/combatant.html. http://www.iihr.uiowa.edu/gothenburg2000/5415/combatant.html.]

<div id="17"></div>
[[#cite-17|[17]]] ''Korea Research Institute of Ships and Ocean Engineering (KRISO)''. [http://www.iihr.uiowa.edu/gothenburg2000/KVLCC/tanker.html. http://www.iihr.uiowa.edu/gothenburg2000/KVLCC/tanker.html.]