A non-linear finite element method on unstructured meshes for added resistance in waves

<!-- metadata commented in wiki content


<span id='_Toc411871622'></span><span id='_Toc418534736'></span><span id='_GoBack'></span><big>'''A NON-LINEAR FINITE ELEMENT MEHOD ON UNSTRUCTURED MESHES FOR ADDED RESISTANCE IN WAVES'''</big>

Julio Garcia-Espinosa <sup>a,b</sup>, Borja Servan-Camas <sup>a,*</sup>

<sup>a</sup> ''Centre Internacional de Mètodes Numèrics en Enginyeria (CIMNE), Gran Capitán s/n, 08034 Barcelona, Spain.''

<sup>b</sup> ''Universitat Politècnica de Catalunya - BarcelonaTech (UPC), Barcelona, Spain.''

''<sup>*</sup> Corresponding author. Tel.: +34 93 405 4627.''

''E-mail address: [mailto:bservan@cimne.upc.edu bservan@cimne.upc.edu] (B. Serván-Camas).''

''URL: [http://www.cimne.com/ http://www.cimne.com/] (B. Serván-Camas).''

-->===ABSTRACT===

In this work a finite element method (FEM) is proposed to solve the problem of estimating the added resistance of a ship in waves in the time domain and using unstructured meshes. Two different schemes are used to integrate the corresponding free surface kinematic and dynamic boundary conditions: the first one based on streamlines integration (SLI); and the second one based on the Streamline-Upwind Petrov-Galerking (SUPG) stabilization. The proposed numerical schemes have been validated in different test cases, including towing tank tests with monochromatic waves. The results obtained in this work show the suitability of the present method to estimate added resistance in waves in a computationally affordable manner.

'''Keywords: '''Added resistance in waves, FEM, unstructured mesh, potential flow.

===1. Introduction===

In order to get a good prediction of the total resistance of a ship in a seaway the added resistance in waves has to be taken into account. And although the added resistance in waves is a second order component, it may reach up to 40% of the total average resistance for some specific wave frequencies where wave radiation increases.

Despite of the importance of the added resistance in waves and its impact on the long term performance of a ship, there are not many works in the literature coping with this problem using computational methods. Some formulation have been derived by different authors as an extension of the mean drift force induced by waves in the presence of a water current (Maruo 1957; Journée 1992). These approaches basically modify damping and added mass coefficients, as well as the mean drift forces, to account for the current effect. However, convective terms in the free surface boundary conditions are neglected, limiting this approach to low Froude numbers. A more sophisticated approach is to solve the three dimensional potential flow problem in the frequency domain using the boundary element method (BEM), and obtain the added resistance using Maruo’s approach thereafter (Maruo 1960; Maruo 1963; Liu at al. 2011). However, Maruo’s approach is still limited to low Froude numbers, since it does not account for the free surface convective terms either. Few works can be found trying to solve the added resistance in waves in the time domain considering the convective terms in the free surface boundary condition and using potential flow (Joncquez 2009; Seo et al. 2013; Afshar 2014; Park et al. 2016). The latter approaches are based on the BEM, consider linearized convective terms, and are constrained to structured meshes.

Non-linear viscous CFD solvers have also been applied to analyze details of the added resistance in waves (Sadat-Hosseini et al. 2013; el Moctar et al. 2015;Tezdogan et al. 2015; Kim et al. 2017a; kim et al. 2017b; Kim et al. 2017c; Ozdemir and Barlas 2017; Hizir et al. 2018). Recently Sigmund and el Moctar (2018) carried out state of the art numerical and experimental investigations on four types of ships using the finite volume method. The effect on added resistance of the ship's speed, wave steepness, skin friction and diffraction-radiation was analyzed in detail. Among other findings, it was found that:  calm water and added resistance must be computed on the same numerical grid; radiation forces are more influenced by the ship's speed than by the diffraction forces; in short waves and at model scale frictional resistance accounts for some 20% of the total added resistance, while this percentage is not as high for full scale; and, finally, that the assumption of a quadratic dependence of added resistance on wave amplitude holds for moderate and long waves, while it does not in short waves because wave diffraction is dominant. However, since the use of CFD tools is computationally expensive, its application is limited to research analysis. On the other hand, potential flow solver are more suitable to evaluate systematically a large number of navigation conditions.

The main motivation of this work is the development of a new time-domain computational model capable of solving the wave making and added resistance problem efficiently. This can be achieved by offering am accurate approach not subject to the use of structured meshes.

The proposed method is based on potential flow and applies the free surface boundary conditions on the mean water level for computational efficiency. The main novelties of the model are the use of the Finite Element Method (FEM) on unstructured meshes (Servan-Camas and Garcia-Espinosa 2013, Garcia-Espinosa et al. 2015; Servan-Camas 2016; Oñate et al. 2018) in the context of potential flow. The proposed numerical scheme is derived using a non-linear (NL) formulation for the convective terms, although the Neumann-Kelvin (NK) and double-body (DB) flow linearizations have been analyzed as well.

The first scheme proposed for integrating the free surface boundary conditions is based on computing the convective term as a derivative along the streamlines. To this end a third order finite difference (FD) scheme using three points upwind and one point downwind has been tailored. The second scheme proposed is based on the FEM with Streamline-Upwind Petrov-Galerkin (SUPG) stabilization (Brooks and Hughes 1982).

<span id='_Toc434595907'></span>

===2. Mathematical models for wave resistance problems===

<span id='_Toc434595910'></span>

====2.1 Problem statement====

We consider the first-order diffraction-radiation problem of a body moving on the horizontal plane. Potential flow is assumed and the problem is formulated in terms of the velocity potential and free surface elevation. The following assumptions are made on the order of magnitude of the velocity components and free surface elevation:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">\left| {\boldsymbol{V}}_{b}\sim O(1)\right|</math> , <math display="inline">\left| \boldsymbol{u}\sim O(1)\right|</math> , <math display="inline">{\varphi }_{x}\sim O(1)</math>, <math display="inline">{\varphi }_{y}\sim O(1)</math>,
|}
|  rowspan='2' style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(1)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">{\varphi }_{z}\sim O(\epsilon )</math>, <math display="inline">{\varphi }_{\alpha \beta }\sim O(\epsilon )</math>, <math display="inline">{\xi }_{x}\sim O(\epsilon )</math>, <math display="inline">{\xi }_{y}\sim O(\epsilon )</math> ,
|}
|}


<span id='_Toc434595911'></span>where <math display="inline">\varphi</math>  the velocity potential, ξ is is the free surface elevation,  <math display="inline">{\boldsymbol{V}}_{b}</math> is the forward speed of the moving body and <math display="inline">\boldsymbol{u}</math>''' '''is the uniform water current vector. The velocity potential and free surface elevation are split into incident and diffraction-radiation components:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">\varphi =\psi +\phi</math> ,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(2)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">\xi =\zeta +\eta</math> ,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(3)
|}


where <math display="inline">\phi</math>  and <math display="inline">\eta</math>   are the diffraction-radiation velocity potential and free surface elevation, and <math display="inline">\, \psi</math>  and <math display="inline">\zeta</math>   are the incident wave potential and free surface elevation (modeled as a set of Airy waves).

<span id='_Toc411871626'></span><span id='_Ref433611076'></span><span id='_Toc434595912'></span>

====2.2 Governing equations in a moving frame of reference====

The governing equations are solved in a frame of reference fixed to the moving body. Let the two dimensional movement of the local frame be defined by the horizontal velocity <math display="inline">{\boldsymbol{V}}_{b}</math> (see <span id='cite-_Ref411871029'></span>[[#_Ref411871029|Figure 1]]). This frame of reference is assumed to match the global frame at time zero. For an observer sitting on the ship the flow field around is given by the relative motion <math display="inline">{\boldsymbol{U}}_{b}\boldsymbol{=u-}{\boldsymbol{V}}_{b}</math>.  Then the first-order diffraction-radiation equations in the local frame of reference are (see Servan-Camas 2016):

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\Delta \phi =0</math> <math display="inline">in\, \Omega</math> ,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(4)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\partial \phi }{\partial t}+{\boldsymbol{U}}_{b}\cdot {\nabla }_{h}\phi +\frac{1}{2}{\nabla }_{h}\phi {\cdot \nabla }_{h}\phi +</math><math>{\nabla }_{h}\psi {\cdot \nabla }_{h}\phi +g\eta =0</math> <math display="inline">on\, z=0</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(5)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\partial \eta }{\partial t}+\left( {\boldsymbol{U}}_{b}+\right. </math><math>\left. {\nabla }_{h}\phi \right) \cdot {\nabla }_{h}\eta +{\nabla }_{h}\phi \cdot {\nabla }_{h}\zeta -</math><math>\frac{\partial \phi }{\partial z}=0</math> <math display="inline">on\, z=0</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(6)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\nabla \phi \cdot {\boldsymbol{n}}_{p}=-\left( {\mathit{\boldsymbol{v}}}_{p}+\right. </math><math>\left. \boldsymbol{u}+{\mathit{\boldsymbol{v}}}_{\psi }\right) \cdot {\boldsymbol{n}}_{p}</math> <math display="inline">on\, {S}_{B}^{0}</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(7)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\partial \phi }{\partial z}=0</math> <math display="inline">on\, z=-H</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(8)
|}


where <math display="inline">{\mathit{\boldsymbol{v}}}_{\psi }=\nabla \psi</math> , <math display="inline">{S}_{B}^{0}</math> is the body surface below <math display="inline">z=</math><math>0</math>,  <math display="inline">{\mathit{\boldsymbol{v}}}_{p}</math> is the velocity of a point P on <math display="inline">{S}_{B}^{0}</math>, and <math display="inline">H</math> is the water depth The pressure induced by diffracted and radiated waves is

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">{P}_{p}=-\rho \left( \frac{\partial \phi }{\partial t}+\right. </math><math>\left. \left( {\boldsymbol{U}}_{b}+{\nabla }_{h}\phi +{\nabla }_{h}\psi \right) \cdot {\nabla }_{h}\phi \mathit{\boldsymbol{\, }}-\right. </math><math>\left. \frac{1}{2}{\nabla }_{h}\phi \cdot {\nabla }_{h}\phi \right)</math>          <math display="inline">in\, \Omega</math> .
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(9)
|}


{| style="width: 73%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;width: 0%;"| [[Image:Draft_Garcia-Espinosa_464519160-image1.jpeg|432px]] 
|}


<span id='_Ref411871029'></span><div id="_Toc429661856" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 1: Global and local frame of reference.</span></div>

<span id='_Ref412102917'></span><span id='_Toc434595913'></span>

====2.3 Flow linearization====

The previous governing equations have been obtained under the assumption <math display="inline">\nabla \phi \sim O(1)</math>. Then the free surface boundary conditions remain non-linear and can be written as follows

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\partial \phi }{\partial t}+\left( {\boldsymbol{U}}_{b}+\right. </math><math>\left. {\nabla }_{h}\phi \right) \cdot {\nabla }_{h}\phi -\frac{1}{2}{\nabla }_{h}\phi {\cdot \nabla }_{h}\phi +</math><math>{\nabla }_{h}\psi {\cdot \nabla }_{h}\phi +g\eta =0</math> <math display="inline">on\, z=0</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref510773805'></span>(10)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\partial \eta }{\partial t}+\left( {\boldsymbol{U}}_{b}+\right. </math><math>\left. {\nabla }_{h}\phi \right) \cdot {\nabla }_{h}\eta +{\nabla }_{h}\phi \cdot {\nabla }_{h}\zeta -</math><math>\frac{\partial \phi }{\partial z}=0</math> <math display="inline">on\, z=0</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref510773816'></span>(11)
|}


where <math display="inline">{\boldsymbol{U}}_{b}+{\nabla }_{h}\phi</math>  represents the convective velocity and must be updated every iteration to account for the variations of <math display="inline">{\boldsymbol{U}}_{b}</math>''' '''and <math display="inline">{\nabla }_{h}\phi</math> . Retaining this assumption allows for simulating flows where <math display="inline">{\boldsymbol{U}}_{\boldsymbol{b}}</math>''' '''is not steady. However, linearization of the convective velocity is a common practice when <math display="inline">{\boldsymbol{U}}_{b}</math>''' '''is steady and under assumptions of low speed or slender bodies. Two commonly used linearizations are:

The NK linearization: assumes <math display="inline">{\nabla }_{h}\phi \sim O(\epsilon )</math>. Then the convective velocity becomes the apparent velocity <math display="inline">{\boldsymbol{U}}_{b}\boldsymbol{=u-}{\boldsymbol{V}}_{b}</math>.

The DB linearization: assumes <math display="inline">{\nabla }_{h}\phi ={\nabla }_{h}{\phi }^{DB}+</math><math>{\nabla }_{h}{\phi }^{\ast }</math> with <math display="inline">{\nabla }_{h}{\phi }^{DB}\sim O(1)</math> and <math display="inline">{\nabla }_{h}{\phi }^{\ast }\sim O(\epsilon )</math> . The <math display="inline">{\nabla }_{h}{\phi }^{DB}</math>  flow field is obtained from the double-body approximation (see Servan-Camas 2016).

<span id='_Toc424744183'></span><span id='_Toc434595914'></span>

====2.4 Second-order pressure correction over body====

A second-order correction on the pressure calculation is proposed in order to better estimate hydrodynamic loads. The idea is to retain second-order terms depending on the first-order solution. Although this will not provide a full second-order solution, it will retain the terms contributing to mean drift loads ( <math display="inline">{\boldsymbol{U}}_{b}=</math><math>0</math>), and added resistance in waves ( <math display="inline">{\boldsymbol{U}}_{b}\mathit{\boldsymbol{\not =}}0</math>). The second-order pressure correction is:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{P}_{p}^{c\left( 2\right) }\, =-\rho {\mathit{\boldsymbol{r}}}_{p}^{1}\cdot \nabla \left( {P}_{p}^{I}+\right. </math><math>\left. {P}_{p}^{DR}\right) -\rho \left( \frac{1}{2}{\phi }_{z}{\phi }_{z}+{\phi }_{z}{\psi }_{z}+\right. </math><math>\left. \frac{1}{2}\nabla \psi \cdot \nabla \psi \right) ,</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(12)
|}


where <math display="inline">{P}_{p}^{I}</math> and <math display="inline">{P}_{p}^{DR}</math> are pressures induced by the incident and diffracted-radiated waves respectively,  <math display="inline">{\mathit{\boldsymbol{r}}}_{p}^{1}</math> is the first-order displacement of a point P over the body surface, and <math display="inline">{P}_{p}^{c\left( 2\right) }</math> is the pressure correction term. The first term of the right hand side account for the pressure variation due to first order movements, and the second term accounts for second-order terms depending on first-order terms.

<span id='_Toc424744186'></span><span id='_Toc434595915'></span><span id='_Ref509562039'></span>

====2.5 Second-order pressure correction on waterline====

The movements of a ship and the free surface elevation modify the wetted surface. In order to take into account this change, a linear pressure is assumed from the waterline to the free surface (see <span id='cite-_Ref423598660'></span>[[#_Ref423598660|Figure 2]]):

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <br/> <math display="inline">{P}_{p}^{wl}={P}_{p}^{I}+{P}_{p}^{DR}-\rho g{r}_{pz}</math>, 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(13)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{P}_{p}^{wl\left( 2\right) }\, \left( z\right) ={P}_{p}^{wl}-\left( {P}_{p}^{wl}-\right. </math><math>\left. {P}_{p}^{fs}\right) \frac{z}{\xi -{r}_{pz}},</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(14)
|}


where <math display="inline">{P}_{p}^{I}</math> and <math display="inline">{P}_{p}^{DR}</math>  are evaluated at z=0, <math display="inline">{P}_{p}^{wl}</math> is the first-order pressure at z=0, and <math display="inline">{P}_{p}^{wl\left( 2\right) }</math> represents the vertical pressure distribution over the water line (see Servan-Camas 2016).

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_464519160-image2.jpeg|432px]] </div>

<span id='_Ref423598660'></span><div id="_Toc429661858" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 2: Pressure estimation at the waterline area.''' '''</span></div>

<span id='_Ref429553979'></span><span id='_Ref430774547'></span><span id='_Toc434595916'></span><span id='_Toc411494943'></span>

===3. Numerical schemes===

====3.1 Problem statement====

Let <math display="inline">{Q}_{h}^{\ast }</math> be the finite element (FE) space of interpolating functions, defined in the usual manner. From this space we can construct the subspace <math display="inline">{Q}_{h,\phi }</math> that incorporates the Dirichlet conditions for the velocity potential <math display="inline">\phi</math> . The space of test functions, denoted by <math display="inline">{Q}_{h}</math>, is constructed as <math display="inline">{Q}_{h,\phi }</math>, but with functions vanishing on the boundary where the Dirichlet conditions are imposed. Then the discrete variational finite element problem associated with the Eqs. (27)-(31) can be written as follows (Belibassakis et al. 2013): Find  <math display="inline">\left[ {\phi }_{h}\right] \in {Q}_{h,\phi }</math>, by solving

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">\int_{\Omega }^{}\nabla {\nu }_{h}\cdot \nabla {\phi }_{h}d\Omega =</math><math>\int_{{\Gamma }^{B}}^{}{\nu }_{h}\cdot {\hat{\phi }}_{n}^{B}d\Gamma +\int_{{\Gamma }^{R}}^{}{\nu }_{h}\cdot {\hat{\phi }}_{n}^{R}d\Gamma +\int_{{\Gamma }^{{Z}_{o}}}^{}{\nu }_{h}\cdot {\hat{\phi }}_{n}^{{Z}_{o}}d\Gamma \quad \, \forall {\nu }_{h}\in {Q}_{h}</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(15)
|}


where <math display="inline">{\hat{\phi }}_{n}^{B}</math>, <math display="inline">{\hat{\phi }}_{n}^{R}</math>, <math display="inline">{\hat{\phi }}_{n}^{{Z}_{o}}</math> are the potential normal gradients corresponding to the Neumann boundary conditions on bodies, radiation boundary and free surface, respectively. At this point it is useful to introduce the associated matrix form (Zienkiewicz et al. 2005a)

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">\overline{\overline{\boldsymbol{L}}}\boldsymbol{\phi }=</math><math>{\boldsymbol{b}}^{{Z}_{0}}+{\boldsymbol{b}}^{B}+{\boldsymbol{b}}^{R}</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(16)
|}


where <math display="inline">\overline{\overline{\boldsymbol{L}}}</math> is the standard FE Laplacian matrix, and <math display="inline">{\boldsymbol{b}}^{{Z}_{0}}</math> is a vector accounting for the right hand side of the dynamic condition of the free surface, and <math display="inline">{\boldsymbol{b}}^{B}</math> and <math display="inline">{\boldsymbol{b}}^{R}</math> are the vectors resulting from integrating the other boundary condition terms.

Two different schemes have been developed to integrate the free surface boundary conditions (Eq. <span id='cite-_Ref510773805'></span>[[#_Ref510773805|(10)]] and Eq. <span id='cite-_Ref510773816'></span>[[#_Ref510773816|(11)]]): one based on FD streamline integration (FD-SLI), and the other based on FEM with SUPG stabilization scheme (FEM-SUPG). Both schemes have been tailored to account for the non-linear convective terms and to be used on unstructured meshes.

<span id='_Ref511204222'></span>

====3.2 FD-SLI====

The first-order free surface boundary conditions can be written in a general form as follows:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\partial \phi }{\partial t}+\boldsymbol{U}\cdot {\nabla }_{h}\phi +\frac{{P}_{fs}}{\rho }+g\eta +R=</math><math>0</math> <math display="inline">on\, z=0</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref458067729'></span>(17)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\partial \eta }{\partial t}+\boldsymbol{U}\cdot {\nabla }_{h}\eta -\frac{\partial \phi }{\partial z}+S=</math><math>0</math> <math display="inline">on\, z=0</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref458067730'></span>(18)
|}


where <math display="inline">R</math> and <math display="inline">S</math> represent remaining terms depending on the linearization used. The numerical schemes adopted for solving the kinematic-dynamic free surface boundary conditions are based on Adams-Bashforth-Moulton schemes, using an explicit scheme for the kinematic condition, and an implicit one for the dynamic condition. Then <math display="inline">{\phi }^{n+1}</math> is imposed as a Dirichlet boundary condition. The schemes read as follows:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">{\phi }^{n+1}+\Delta t{\left( \mathit{\boldsymbol{U}}\cdot {\nabla }_{h}\phi \right) }^{n+1}=</math><math>{\phi }^{n}-\Delta t\left( \frac{1}{\rho }{P}_{fs}^{n+1}-g{\eta }^{n+1}-{R}^{n+1}\right)</math> <math display="inline">on\, z=0</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref412020700'></span>(19)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\eta }^{n+1}={\eta }^{n}-\Delta t{\left( \boldsymbol{U}\cdot {\nabla }_{h}\eta \right) }^{n}+</math><math>\Delta t\left( {\phi }_{z}^{n}-{S}^{n}\right)</math> <math display="inline">on\, z=0</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref424802532'></span>(20)
|}


where <math display="inline">\boldsymbol{U}</math>''' '''is the convective velocity and the convective terms are evaluated by differentiating along streamlines as:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\left( \boldsymbol{U}\cdot {\nabla }_{h}\phi \right) }^{n+1}={\left| \mathit{\boldsymbol{U}}\right| }^{n+1}{\partial }_{L}{\phi }^{n+1}\, ,</math>
|}
|  rowspan='2' style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(21)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\left( \boldsymbol{U}\cdot {\nabla }_{h}\eta \right) }^{n}={\left| \mathit{\boldsymbol{U}}\right| }^{n}{\partial }_{L}{\eta }^{n}\, ,</math>
|}
|}


where <math display="inline">{\partial }_{L}</math> denotes the derivative along the streamline. This streamline derivative is estimated using a two points upstream and one point downstream FD operator. <span id='cite-_Ref326165045'></span>[[#_Ref326165045|Figure 3]] shows the tracing of the streamline at node 0. The left (-1) and forward left (-2) points are upstream points, while the right (1) point corresponds to the downstream point. Then the streamline differential operator reads as:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\partial }_{L}{\phi }_{0}={w}_{1}{\phi }_{1}+{w}_{0}{\phi }_{0}+{w}_{-1}{\phi }_{-1}+</math><math>{w}_{-2}{\phi }_{-2},</math>
|}
|  rowspan='2' style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(22)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\partial }_{L}{\eta }_{0}={w}_{1}{\eta }_{1}+{w}_{0}{\eta }_{0}+{w}_{-1}{\eta }_{-1}+</math><math>{w}_{-2}{\eta }_{-2},</math>
|}
|}


where <math display="inline">{\phi }_{1}</math>, <math display="inline">{\phi }_{-1}</math>, <math display="inline">{\phi }_{-2}</math> are interpolated between <math display="inline">\left( {\phi }_{1a},{\phi }_{1b}\right)</math> ,  <math display="inline">\left( {\phi }_{-1a},{\phi }_{-1b}\right)</math> , and <math display="inline">\left( {\phi }_{-2a},{\phi }_{-2b}\right)</math>  respectively. In matrix form, the resulting scheme can be written as:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">\left( \overline{\overline{\mathit{\boldsymbol{I}}}}+\right. </math><math>\left. \Delta t{\overline{\overline{\boldsymbol{W}}}}^{\mathit{\boldsymbol{n+1}}}\right) {\mathit{\boldsymbol{\phi }}}^{n+1}=</math><math>\overline{\overline{\mathit{\boldsymbol{I}}}}\left( {\mathit{\boldsymbol{\phi }}}^{n}-\right. </math><math>\left. \Delta t\left( \frac{1}{\rho }{\mathit{\boldsymbol{P}}}_{fs}^{n+1}-g{\mathit{\boldsymbol{\eta }}}^{n+1}-\right. \right. </math><math>\left. \left. {\mathit{\boldsymbol{R}}}^{n+1}\right) \right)</math> <math display="inline">on\, z=0</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref412022091'></span>(23)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\overline{\overline{\mathit{\boldsymbol{I}}}}{\mathit{\boldsymbol{\, \eta }}}^{n+1}=</math><math>\left( \overline{\overline{\mathit{\boldsymbol{I}}}}-\Delta t{\overline{\overline{\boldsymbol{W}}}}^{\mathit{\boldsymbol{n}}}\right) {\mathit{\boldsymbol{\, \eta }}}^{n}+</math><math>\Delta t\left( {\mathit{\boldsymbol{\phi }}}_{z}^{n}-{\boldsymbol{S}}^{n}\right)</math> <math display="inline">on\, z=0</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(24)
|}


where <math display="inline">\overline{\overline{\boldsymbol{W}}}</math> and <math display="inline">\overline{\overline{\mathit{\boldsymbol{I}}}}</math> are the streamline convective and identity matrices respectively. The stencils are tailored using Taylor series expansion to get a third-order FD scheme along the streamline.

The associated matrix form to the FE formulation for the governing equations is:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">\overline{\overline{{\boldsymbol{L}}^{SL}}}\boldsymbol{\phi }=</math><math>{\boldsymbol{b}}^{SL}+{\boldsymbol{b}}^{B}+{\boldsymbol{b}}^{R}</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(25)
|}


where <math display="inline">\overline{\overline{{\boldsymbol{L}}^{SL}}}</math> is the standard Laplacian matrix modified to account for the left hand side of Eq.<span id='cite-_Ref412022091'></span>[[#_Ref412022091|(23)]], <math display="inline">{\boldsymbol{b}}^{SL}</math> is a vector containing the right hand side of Eq. <span id='cite-_Ref412022091'></span>[[#_Ref412022091|(23)]], and <math display="inline">{\boldsymbol{b}}^{B}</math> and <math display="inline">{\boldsymbol{b}}^{R}</math> are the vectors resulting from integrating the other boundary condition terms.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_464519160-image3.jpeg|432px]] </div>

<span id='_Ref326165045'></span><span id='_Ref320697593'></span><div id="_Toc429661859" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 3: SLI discretization</span></div>

<span id='_Ref429554348'></span><span id='_Toc434595920'></span>

====3.3 FEM-SUPG====

Alternatively a FEM-SUPG (Hughes and Mallet 1986) has been also developed for integrating the free surface boundary conditions. The dynamic and kinematic boundary conditions can be seen as  convective-transport equations of the form:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\partial \chi }{\partial t}+\boldsymbol{U}\cdot {\nabla }_{h}\chi +Q=0.</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(26)
|}


The standard FE problem associated to the above equation is

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">\forall i\, \int_{\Omega }^{}{N}^{i}\, \left( {N}^{j}\frac{\partial {\chi }_{j}}{\partial t}+\boldsymbol{U}\cdot {\nabla }_{h}{N}^{j}{\chi }_{j}+{N}^{j}{Q}_{j}\right) d\sigma =</math><math>0,</math> 
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(27)
|}


where <math display="inline">{N}^{i}\,</math>  is the corresponding shape function associated to the mesh node <math display="inline">i</math> and <math display="inline">{\chi }_{j}</math> and <math display="inline">{Q}_{j}</math> represent the values of <math display="inline">{\chi }_{h}</math> and <math display="inline">{Q}_{h}</math> at node <math display="inline">j</math>, respectively. The above equation is discretized in time as follows:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\forall i\, \int_{\Omega }^{}{N}^{i}{N}^{j}\frac{{{\chi }_{j}}^{n+1}-{{\chi }_{j}}^{n}}{\Delta t}d\sigma +</math><math>\, \int_{\Omega }^{}{N}^{i}\left( {\boldsymbol{U}}^{n+\theta }\cdot {\nabla }_{h}{N}^{j}\right) {{\chi }_{j}}^{n+\theta }d\sigma +</math><math>\, \int_{\Omega }^{}{{N}^{i}N}^{j}{Q}_{j}^{n+\theta }d\sigma =0,</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(28)
|}


where <math display="inline">{A}^{n+\theta }=\left( 1-\theta \right) {A}^{n}+</math><math>\theta {A}^{n+1}</math>. It is well known that the numerical solution of the previous equation may show spurious oscillations (Zienkiewicz et al. 2005b; Donea and Huerta 2003). Hence the SUPG formulation is used to introduce the necessary stabilization. The resulting scheme is

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\forall i\, \, \left[ \int_{\Omega }^{}{N}^{i}{N}^{j}d\sigma +\sum _{e}^{}\frac{{h}_{e}}{2\left| {\left( \hat{{\boldsymbol{U}}^{e}}\right) }^{n+\theta }\right| }\int_{{\Omega }^{e}\, }^{}\left( {\left( {\boldsymbol{U}}^{e}\right) }^{n+\theta }\cdot {\nabla }_{h}{N}^{i}\right) {N}^{j}d\sigma \right] \left( {\chi }_{j}^{n+1}-\right. </math><math>\left. {\chi }_{j}^{n}\right)</math>

<math>+\Delta t\, \left[ \int_{\Omega }^{}{N}^{i}\left( {\left( {\boldsymbol{U}}^{\boldsymbol{e}}\right) }^{n+\theta }\cdot {\nabla }_{h}{N}^{j}\right) d\sigma \right] {\chi }_{j}^{n+\theta }</math>

<math>+\Delta t\, \left[ \sum _{e}^{}\frac{{h}_{e}}{2\left| {(\hat{{\boldsymbol{U}}^{e}})}^{n+\theta }\right| }\int_{{\Omega }^{e}}^{}\left( {\left( {\boldsymbol{U}}^{e}\right) }^{n+\theta }\cdot {\nabla }_{h}{N}^{i}\right) \left( {\left( {\boldsymbol{U}}^{e}\right) }^{n+\theta }\cdot {\nabla }_{h}{N}^{j}\right) d\sigma \right] {\chi }_{j}^{n+\theta }</math>

<math>+\Delta t\left[ \int_{\Omega }^{}{N}^{i}{N}^{j}d\sigma +\sum _{e}^{}\frac{{h}_{e}}{2\left| {\mathit{\boldsymbol{(}}\hat{{\boldsymbol{U}}^{e}}\mathit{\boldsymbol{)}}}^{n+\theta }\right| }\int_{{\Omega }^{e}\, }^{}\left( {\left( {\boldsymbol{U}}^{e}\right) }^{n+\theta }\cdot {\nabla }_{h}{N}^{i}\right) {N}^{j}d\sigma \right] {Q}_{j}^{n+\theta }=</math><math>0.</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(29)
|}


In matrix form:

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\left( \overline{\overline{\mathit{\boldsymbol{M}}}}+{\overline{\overline{\mathit{\boldsymbol{M}}}}}_{supg}^{n\mathit{\boldsymbol{+}}\theta }\right) {\mathit{\boldsymbol{\chi }}}^{n+1}\boldsymbol{+\Delta }t\left( {\overline{\overline{\mathit{\boldsymbol{C}}}}}^{n\mathit{\boldsymbol{+}}\theta }+\right. </math><math>\left. {\overline{\overline{\mathit{\boldsymbol{C}}}}}_{supg}^{n\mathit{\boldsymbol{+}}\theta }\right) {\theta \mathit{\boldsymbol{\chi }}}^{n+1}=</math>

<math display="inline">\left( \overline{\overline{\mathit{\boldsymbol{M}}}}+\right. </math><math>\left. {\overline{\overline{\mathit{\boldsymbol{M}}}}}_{supg}^{n+\theta }\right) {\mathit{\boldsymbol{\chi }}}^{n}\boldsymbol{-\Delta }t\left( {\overline{\overline{\mathit{\boldsymbol{C}}}}}^{n\mathit{\boldsymbol{+}}\theta }+\right. </math><math>\left. {\overline{\overline{\mathit{\boldsymbol{C}}}}}_{supg}^{n\mathit{\boldsymbol{+}}\theta }\right) {\left( 1-\theta \right) \mathit{\boldsymbol{\chi }}}^{n}-</math><math>\left( \overline{\overline{\mathit{\boldsymbol{M}}}}+{\overline{\overline{\mathit{\boldsymbol{M}}}}}_{supg}^{n+\theta }\right) {\mathit{\boldsymbol{Q}}}^{n+\theta },</math>''' '''
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(30)
|}


Where

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\left( {U}^{e}\right) }^{n+\theta }=\sum _{{k}_{e}}^{}{N}^{{k}_{e}}{U}_{{k}_{e}}^{n+\theta }</math>

<math>\overline{\overline{\mathit{\boldsymbol{M}}}}=\int_{\Omega }^{}{N}^{i}{N}^{j}d\sigma ,</math>

<math>{\overline{\overline{\mathit{\boldsymbol{M}}}}}_{supg}^{n+\theta }=\sum _{e}^{}\left[ \frac{{h}_{e}}{2\left| {\left( \hat{{\mathit{\boldsymbol{U}}}^{e}}\right) }^{n+\theta }\right| }\int_{{\Omega }^{e}}^{}\left( {\left( {\mathit{\boldsymbol{U}}}^{e}\right) }^{n+\theta }\cdot {\nabla }_{h}{N}^{i}\right) {N}^{j}d\sigma \right] ,</math>

<math>{\overline{\overline{\mathit{\boldsymbol{C}}}}}^{n+\theta }=\sum _{e}^{}\int_{{\Omega }^{e}}^{}{N}^{i}\left( {\mathit{\boldsymbol{U}}}_{e}^{n+\theta }\cdot {\nabla }_{h}{N}^{j}\right) ,</math>

<math>{\overline{\overline{\mathit{\boldsymbol{C}}}}}_{supg}^{n+\theta }=\sum _{e}^{}\left[ \frac{{h}_{e}}{2\left| {\left( \hat{{\mathit{\boldsymbol{U}}}^{e}}\right) }^{n+\theta }\right| }\int_{{\Omega }^{e}}^{}\left( {\left( {\mathit{\boldsymbol{U}}}^{e}\right) }^{n+\theta }\cdot {\nabla }_{h}{N}^{i}\right) \left( {\mathit{\boldsymbol{U}}}_{e}^{n+\theta }\cdot {\nabla }_{h}{N}^{j}\right) d\sigma \right] .</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(31)
|}


The FEM-SUPG dynamic and kinematic boundary conditions become

{| class="formulaSCP" style="width: 100%;border-collapse: collapse;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\left( \overline{\overline{\mathit{\boldsymbol{M}}}}+{\overline{\overline{\mathit{\boldsymbol{M}}}}}_{supg}^{n\mathit{\boldsymbol{+}}\theta }\right) {\mathit{\boldsymbol{\phi }}}^{n+1}\boldsymbol{+\Delta }t\left( {\overline{\overline{\mathit{\boldsymbol{C}}}}}^{n\mathit{\boldsymbol{+}}\theta }+\right. </math><math>\left. {\overline{\overline{\mathit{\boldsymbol{C}}}}}_{supg}^{n\mathit{\boldsymbol{+}}\theta }\right) {\theta \mathit{\boldsymbol{\phi }}}^{n+1}=</math>

'''            ''' <math display="inline">\left( \overline{\overline{\mathit{\boldsymbol{M}}}}+\right. </math><math>\left. {\overline{\overline{\mathit{\boldsymbol{M}}}}}_{supg}^{n+\theta }\right) {\mathit{\boldsymbol{\phi }}}^{n}\boldsymbol{-\Delta }t\left( {\overline{\overline{\mathit{\boldsymbol{C}}}}}^{n\mathit{\boldsymbol{+}}\theta }+\right. </math><math>\left. {\overline{\overline{\mathit{\boldsymbol{C}}}}}_{supg}^{n\mathit{\boldsymbol{+}}\theta }\right) {\left( 1-\theta \right) \mathit{\boldsymbol{\phi }}}^{n}</math>

<math display="inline">-\left( \overline{\overline{\mathit{\boldsymbol{M}}}}+\right. </math><math>\left. {\overline{\overline{\mathit{\boldsymbol{M}}}}}_{supg}^{n+\theta }\right) \left( \frac{1}{\rho }{\mathit{\boldsymbol{P}}}^{n+\theta }+\right. </math><math>\left. g{\mathit{\boldsymbol{\eta }}}^{n+\theta }+{\mathit{\boldsymbol{R}}}^{n+\theta }\right) \mathit{\boldsymbol{.}}</math>''' '''
|}
|  colspan='2'  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref412022628'></span>(32)
|-
|  colspan='2'  style="text-align: center;width: 5px;text-align: right;white-space: nowrap;"| <math>\left( \overline{\overline{\mathit{\boldsymbol{M}}}}+{\overline{\overline{\mathit{\boldsymbol{M}}}}}_{supg}^{n\mathit{\boldsymbol{+}}\theta }\right) {\boldsymbol{\eta }}^{n+1}\boldsymbol{+\Delta }t\left( {\overline{\overline{\mathit{\boldsymbol{C}}}}}^{n\mathit{\boldsymbol{+}}\theta }+\right. </math><math>\left. {\overline{\overline{\mathit{\boldsymbol{C}}}}}_{supg}^{n\mathit{\boldsymbol{+}}\theta }\right) {\theta \boldsymbol{\eta }}^{n+1}=</math>

'''            ''' <math display="inline">\left( \overline{\overline{\mathit{\boldsymbol{M}}}}+\right. </math><math>\left. {\overline{\overline{\mathit{\boldsymbol{M}}}}}_{supg}^{n+\theta }\right) {\boldsymbol{\eta }}^{n}\boldsymbol{-\Delta }t\left( {\overline{\overline{\mathit{\boldsymbol{C}}}}}^{n\mathit{\boldsymbol{+}}\theta }+\right. </math><math>\left. {\overline{\overline{\mathit{\boldsymbol{C}}}}}_{supg}^{n\mathit{\boldsymbol{+}}\theta }\right) {\left( 1-\theta \right) \boldsymbol{\eta }}^{n}</math>

<math display="inline">+\left( \overline{\overline{\mathit{\boldsymbol{M}}}}+\right. </math><math>\left. {\overline{\overline{\mathit{\boldsymbol{M}}}}}_{supg}^{n\mathit{\boldsymbol{+}}\theta }\right) \left( \frac{1}{\rho }{\mathit{\boldsymbol{\phi }}}_{z}^{n\mathit{\boldsymbol{+}}\theta }-\right. </math><math>\left. {\mathit{\boldsymbol{S}}}^{n+\theta }\right) \mathit{\boldsymbol{.}}</math>
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(33)
|}


In this work an explicit ( <math display="inline">\theta =0</math>) and implicit ( <math display="inline">\theta =</math><math>1</math>) schemes are used for the kinematic and dynamic conditions respectively. The matrix form associated to the FE formulation of the governing equations is:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">\overline{\overline{{\boldsymbol{L}}^{SUPG}}}\boldsymbol{\phi }=</math><math>{\boldsymbol{b}}^{SUPG}+{\boldsymbol{b}}^{B}+{\boldsymbol{b}}^{R}</math>,
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(34)
|}


where <math display="inline">\overline{\overline{{\boldsymbol{L}}^{SUPG}}}</math> is the standard Laplacian matrix modified to account for the left hand side of Eq. <span id='cite-_Ref412022628'></span>[[#_Ref412022628|(32)]], <math display="inline">{\boldsymbol{b}}^{SUPG}</math> is a vector accounting for the right hand side of Eq. <span id='cite-_Ref412022628'></span>[[#_Ref412022628|(32)]], and <math display="inline">{\boldsymbol{b}}^{B}</math>  and <math display="inline">{\boldsymbol{b}}^{R}</math>  are the vectors resulting of integrating the other boundary condition terms.

<span id='_Ref430775275'></span><span id='_Toc434595921'></span>

===4. Hydrodynamic loads for wave resistance problems===

<span id='_Toc434595924'></span>

====4.1 First-order loads====

Hydrodynamic forces and moments are obtained from direct pressure integration over the body surface (Bai and Teng 2013). The body gravity center ( <math display="inline">G</math>) will be used as a reference for body movements and moments acting on it. The resulting expressions for forces ( <math display="inline">\mathit{\boldsymbol{F}}</math>) and moments ( <math display="inline">\mathit{\boldsymbol{M}}</math>''')''' can be found in (Servan-Camas 2016).

<span id='_Toc424744197'></span><span id='_Ref430004083'></span><span id='_Toc434595927'></span>

====4.2 Second-order correction====

Second-order correction loads are obtained by calculating the components of the second-order loads depending on the first-order solution as described in section <span id='cite-_Ref509562039'></span>[[#_Ref509562039|2.5]] (<span id='cite-_Ref423598660'></span>[[#_Ref423598660|Figure 2]]). Integrating the corresponding second-order pressure correction:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{F}}}_{C}^{(2)}=\int_{{S}_{B}}^{}{P}_{p}^{c\left( 2\right) }{\mathit{\boldsymbol{n}}}_{p}ds+</math><math>\int_{{\Gamma }_{wl}}^{}{P}_{p}^{wl\left( 2\right) }\frac{{\mathit{\boldsymbol{n}}}_{p}^{0}}{\sqrt{1-{{n}_{pz}^{0}}^{2}}}dl,</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref509569238'></span>(35)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{M}}}_{C}^{(2)}=\int_{{S}_{B}}^{}{P}_{p}^{c\left( 2\right) }\overset{\rightarrow}{\mathit{\boldsymbol{GP}}}\mathit{\boldsymbol{\times }}{\mathit{\boldsymbol{n}}}_{p}ds+</math><math>\int_{{\Gamma }_{wl}}^{}{P}_{p}^{wl\left( 2\right) }\overset{\rightarrow}{\mathit{\boldsymbol{GP}}}\mathit{\boldsymbol{\times }}{\mathit{\boldsymbol{n}}}_{p}\frac{{\mathit{\boldsymbol{n}}}_{p}^{0}}{\sqrt{1-{{n}_{pz}^{0}}^{2}}}dl.</math>
|}
|  style="text-align: right;width: 5px;text-align: right;white-space: nowrap;"|(36)
|}


The mean value of the longitudinal component of <math display="inline">{\mathit{\boldsymbol{F}}}_{C}^{(2)}</math> has a special meaning in the presence of waves: when <math display="inline">\mathit{\boldsymbol{U}}=</math><math>0</math> it becomes the mean drift force due to the waves; and when <math display="inline">\mathit{\boldsymbol{U}}\not =0</math> it becomes the added resistance in waves. Additionaly second-order correction of dynamic and hydrostatic loads resulting from the first-order rotation (see Servan-Camas 2016) are also taken into account in this work.

===5. Verification and validation===

<span id='_Toc434595956'></span>

====5.1 Wave making resistance of an elliptic pressure distribution====

<span id='_Toc434595957'></span>'''5.1.1 Problem description'''

<span id='_Ref424625021'></span>Newman and Poole (1962) derived an analytical solution for the wave resistance of a pressure distribution moving with constant forward speed along the free surface of a canal with constant width and depth. This solution was obtained based on a first-order approximation with uniform flow (NK). In particular, an analytical expression for the elliptic pressure distribution was obtained. The main particulars of the elliptic pressure distribution used to validate the present model are given in <span id='cite-_Ref431291735'></span>[[#_Ref431291735|Table 1]]:

<div id="_Ref431291735" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 1: Particulars for elliptic pressure distribution.</span></div>

{| style="margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Length (L)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Beam (B)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.5 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Water depth (H)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">5 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Channel width (W)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">10 m</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Free surface pressure (</span> <math display="inline">{P}_{fs}</math><span style="text-align: center; font-size: 75%;">)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1 Pa</span>
|}


The wave resistance is obtained by integrating the pressure distribution over the free surface as follows:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{F}_{W}=-\int_{{S}_{FS}}^{}{P}_{fs}{n}_{x}ds.</math>
|}
|  style="text-align: right;vertical-align: top;width: 5px;text-align: right;white-space: nowrap;"|(<span style="text-align: center; font-size: 75%;">37</span>)
|}


It is assumed that the pressure distribution moves in the x direction. Then, in order to obtain the wave resistance, it is mandatory to be able to estimate the free surface deformation within the pressure patch.

<span id='_Toc434595958'></span>'''5.1.2 Convergence analysis'''

Different unstructured meshes have been used for the numerical calculations. It has been found that refinement at the edge of the pressure distribution is mandatory due to the discontinuity in the pressure and its large influence in the free surface deformation, especially at low Froude numbers where wave lengths are smaller. <span id='cite-_Ref430786777'></span>[[#_Ref430786777|Table 2]] provides the particulars for each mesh and <span id='cite-_Ref424625686'></span>[[#_Ref424625686|Figure 4]] shows a detail of the finest one.

<div id="_Ref430786777" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 2: Characteristic sizes of the used meshes.</span></div>

{| style="width: 85%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"|
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Pressure patch edge (m)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Pressure patch (m)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Free surface outside (m)</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Volume (m)</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Mesh 1</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.002</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.02</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.02</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.05</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Mesh 2</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.003</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.03</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.03</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.075</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Mesh 3</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.004</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.04</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.04</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.1</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Mesh 4</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.006</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.06</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.06</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.15</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Mesh 5</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.008</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.08</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.08</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.20</span>
|}


Wave resistance was calculated for all the meshes at Froude number Fr=0.2, using the FD-SLI and FEM-SUPG schemes presented in sections <span id='cite-_Ref511204222'></span>[[#_Ref511204222|3.2]] and <span id='cite-_Ref429554348'></span>[[#_Ref429554348|3.3]] respectively. <span id='cite-_Ref425158075'></span>[[#_Ref425158075|Table 3]] provides the values of the relative errors. <span id='cite-_Ref425169910'></span>[[#_Ref425169910|Figure 5]] shows the free surface elevation computed for Fr=0.2 with the FEM-SUPG and the finer mesh.

<div id="_Ref425158075" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 3: Dimensionless wave resistance </span> <math display="inline">{\mathit{\boldsymbol{C}}}_{\mathit{\boldsymbol{w}}}\mathit{\boldsymbol{=\rho g}}{\mathit{\boldsymbol{F}}}_{\mathit{\boldsymbol{w}}}\mathit{\boldsymbol{/(}}{\mathit{\boldsymbol{P}}}_{\mathit{\boldsymbol{fs}}}^{\mathit{\boldsymbol{2}}}\mathit{\boldsymbol{B)\, }}</math><span style="text-align: center; font-size: 75%;"> at Fr=0.2.</span></div>

{| style="width: 61%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  colspan='2'  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Analytical </span> <math display="inline">{C}_{w}</math><span style="text-align: center; font-size: 75%;">=0.1840</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Mesh 5</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Mesh 4</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Mesh 3</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Mesh 2</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Mesh 1</span>
|-
|  rowspan='2' style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Streamline</span>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"| <math>{C}_{w}</math>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.550</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.640</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.730</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.780</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.820</span>
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">Error</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">15.77 %</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">10.88 %</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">5.99 %</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">3.27 %</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.10 %</span>
|-
|  rowspan='2' style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">FEM-SUPG</span>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"| <math>{C}_{w}</math>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.643</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.710</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.775</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.795</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.820</span>
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|<span style="text-align: center; font-size: 75%;">Error</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">10.74 %</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">7.08 %</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">3.54 %</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">2.46 %</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.10 %</span>
|}


<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_464519160-image4.jpeg|432px]] </div>

<span id='_Ref424625686'></span><div id="_Toc429661881" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 4: Detail of the finest mesh for the elliptic pressure patch.</span></div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
''' [[Image:Draft_Garcia-Espinosa_464519160-image5.jpeg|438px]] '''</div>

<span id='_Ref425169910'></span><div id="_Toc429661882" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 5: Free surface elevation for elliptic pressure distribution at Fr=0.2.</span></div>

<span id='_Toc434595959'></span>'''5.1.3 Verification'''

The wave making resistance has been calculated for several Froude numbers for the finest mesh using the FD-SLI and FEM-SUPG. <span id='cite-_Ref424633190'></span>[[#_Ref424633190|Table 4]] compares analytical and numerical results. It can be observed that the numerical results are able to reproduce the analytical ones with small errors.

<div id="_Ref424633190" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 4: Dimensionless wave resistance </span> <math display="inline">{\mathit{\boldsymbol{C}}}_{\mathit{\boldsymbol{w}}}\mathit{\boldsymbol{\, }}</math><span style="text-align: center; font-size: 75%;">for elliptic pressure distribution.</span></div>

{| style="width: 48%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  rowspan='2' style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Fr</span>
|  colspan='3'  style="border: 1pt solid black;text-align: center;"| <math display="inline">{C}_{w}=</math><math>\rho g{F}_{w}/({P}_{fs}^{2}B)\,</math> 
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">FD-SLI</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">FEM-SUPG</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">Analytical</span>

<span style="text-align: center; font-size: 75%;"> (Newman and Poole 1962)</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.2</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.82</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.82</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.84</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.3</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">2.18</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">2.16</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">2.18</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.4</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.67</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.65</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">1.64</span>
|-
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">0.5</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">2.73</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">2.66</span>
|  style="border: 1pt solid black;text-align: center;"|<span style="text-align: center; font-size: 75%;">2.66</span>
|}


<span id='_Toc434595960'></span>

====5.2 Wave making resistance of a Wigley hull====

===5.2.1 Problem description===

The wave making resistance of a Wigley hull is analyzed. Different simulations have been performed considering the NL, NK, and DB free surface conditions. The FD-SLI and FEM-SUPG schemes have been used for integrating the free surface boundary conditions. Wave resistance coefficients with and without second-order correction have been obtained.

The standard Wigley hull with beam to length ratio B/L=0.1 and amidships section coefficient C<sub>m</sub>=0.6667 was used. An unstructured mesh was generated with an element size around the ship of 0.01 m. The mesh generated consists of 84351 nodes and 476083 tetrahedral elements (see <span id='cite-_Ref425172689'></span>[[#_Ref425172689|Figure 6]]). Only half domain was required due to the symmetry of the problem.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_464519160-image6.jpeg|432px]] </div>

<span id='_Ref425172689'></span><div id="_Toc429661883" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 6: Detail of the mesh generated for the Wigley hull.</span></div>

<span id='_Toc434595962'></span>'''5.2.2 Second-order correction and flow linearizations'''

<span id='cite-_Ref424636948'></span>[[#_Ref424636948|Figure 7]] compares the dimensionless first-order wave making resistance coefficients <math display="inline">{C}_{w1}=</math><math>{F}_{Dx}^{1}/(0.5\rho S{U}^{2})</math> and the dimensionless second-order-corrected wave making resistance coefficients <math display="inline">{C}_{w2}=</math><math>\left( {F}_{Dx}^{1}+{F}_{Cx}^{2}\right) /(0.5\rho S{U}^{2})</math>. They have been obtained using the NL, NK and DB assumptions for the case of a fixed Wigley hull model in still water. The wave making resistance forces <math display="inline">{F}_{Dx}^{1}</math> and <math display="inline">{F}_{Cx}^{2}</math> have been calculated as the longitudinal forces obtained by pressure integration (section <span id='cite-_Ref430775275'></span>[[#_Ref430775275|4]]). It is observed that the increase of resistance due to the second-order correction term is noticeable.

The NK linearization leads to smaller values of wave resistance when compared to the NL and DB approaches, while the NL and DB predict similar results up to Fr=0.3.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_464519160-image7.jpeg|600px]] </div>

<span id='_Ref424636948'></span><span id='_Toc429661884'></span><span style="text-align: center; font-size: 75%;">Figure 7: Wave resistance for Wigley hull using NL, DB, and NK approximations. Dimensionless first-order wave making resistance (</span> <math display="inline">{C}_{w1}</math><span style="text-align: center; font-size: 75%;">), and dimensionless second-order-corrected wave making resistance (</span> <math display="inline">{C}_{w2}</math><span style="text-align: center; font-size: 75%;">).</span>

<span id='_Toc434595963'></span>'''5.2.3 Validation'''

<span id='cite-_Ref424638994'></span>[[#_Ref424638994|Figure 8]] compares the second-order-corrected <math display="inline">{C}_{w2}</math> numerical results obtained in this work with NL flow approximation against several experimental results: ITTC (1984) (collection of experimental data), Shearer and Cross (1965), University of Iowa, University of Tokyo (UT), and the Ship Research Institute of Japan (SRI). The results from Iowa, UT and SRI can be found in Ju (1983). Experimental data shows large dispersion and this might be mainly due to the uncertainty of model testing and the assumption of a constant form factor with speed.

Taking into account the large dispersion of the experimental data a fairly good agreement between the numerical and the experimental results has been obtained. This agreement is better for the higher Froude numbers. Lower Froude numbers require a much finer discretization in order to reproduce the shorter wave pattern. <span id='cite-_Ref425416177'></span>[[#_Ref425416177|Figure 9]] compares the wave profiles at the water line obtained numerically against experimental results. A good agreement is also found in this case.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_464519160-image8.jpeg|432px]] </div>

<span id='_Ref424638994'></span><div id="_Toc429661885" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 8: FD-SLI</span> <span style="text-align: center; font-size: 75%;">and FEM-SUPG versus experimental results.</span></div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_464519160-image9.jpeg|600px]] </div>

<span id='_Ref425416177'></span><div id="_Toc429661886" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 9: Wave profile over Wigley hull.</span></div>

<span id='_Toc434595968'></span>

====5.3 Added resistance in waves of four modified Wigley hulls====

===5.3.1 Problem description===

In this section computed added resistance is compared to that obtained by Journée (1992). Journée carried out an extensive model test campaign measuring the added resistance in waves of four modified Wigley hulls.. The main particulars of the models are given in <span id='cite-_Ref432579756'></span>[[#_Ref432579756|Table 5]].

<div id="_Ref432579756" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 5: Modified Wigley hulls particulars.</span></div>

{| style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"|
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|Wigley I
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|Wigley II
|  style="border: 1pt solid black;text-align: center;"|Wigley III
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|Wigley IV
|-
|  style="border: 1pt solid black;text-align: center;"|Amidships section coefficient C<sub>m</sub>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|0.9090
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|0.9090
|  style="border: 1pt solid black;text-align: center;"|0.6667
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|0.6667
|-
|  style="border: 1pt solid black;text-align: center;"|Length to breadth ratio, L/B
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|10
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|5
|  style="border: 1pt solid black;text-align: center;"|10
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|5
|-
|  style="border: 1pt solid black;text-align: center;"|Length L (m)
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|1
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|1
|  style="border: 1pt solid black;text-align: center;"|1
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|1
|-
|  style="border: 1pt solid black;text-align: center;"|Breadth B (m)
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|0.1
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|0.2
|  style="border: 1pt solid black;text-align: center;"|0.1
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|0.2
|-
|  style="border: 1pt solid black;text-align: center;"|Draught d (m)
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|0.0625
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|0.0625
|  style="border: 1pt solid black;text-align: center;"|0.0625
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|0.0625
|-
|  style="border: 1pt solid black;text-align: center;"|Displacement (m3)
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|0.003504
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|0.007008
|  style="border: 1pt solid black;text-align: center;"|0.002889
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|0.005778
|-
|  style="border: 1pt solid black;text-align: center;"|KG (m)
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|0.05667
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|0.0625
|  style="border: 1pt solid black;text-align: center;"|0.05667
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|0.0625
|-
|  style="border: 1pt solid black;text-align: center;"|Radius of inertia for pitch, ''k''<span style="text-align: center; font-size: 75%;">''yy''</span> (m)
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|0.25
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|0.25
|  style="border: 1pt solid black;text-align: center;"|0.25
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|0.25
|-
|  style="border: 1pt solid black;text-align: center;"|Pitch damping (N/(m/s))
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"| <math>0.02\sqrt{4{I}_{yy}{K}_{55}}</math>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"| <math>0.10\sqrt{4{I}_{yy}{K}_{55}}</math>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"| <math>0.02\sqrt{4{I}_{yy}{K}_{55}}</math>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"| <math>0.02\sqrt{4{I}_{yy}{K}_{55}}</math>
|}


<span id='_Toc434595970'></span>'''5.3.2 Validation'''

The added resistance <math display="inline">{R}_{aw}</math> has been estimated by subtracting the second-order-corrected longitudinal force in still water to the average value of the second-order-corrected longitudinal force obtained in the presence of a monochromatic wave. Then the dimensionless wave resistance is calculated as:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{R}_{aw}^{\ast }=\frac{{R}_{aw}}{\rho g{A}^{2}{B}^{2}/L}\, ,</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(38)
|}


where A is the wave amplitude.

Simulations have been carried out with the streamlines and FEM-SUPG scheme using the NL flow approximations. Pitch damping was introduced to model the viscous dissipation effects, and prevent the excessive pitch movement around resonance obtained otherwise (see <span id='cite-_Ref432579756'></span>[[#_Ref432579756|Table 5]]). The amount of damping introduced was calibrated to best fit the RAOs response across all Froude numbers tested.

Journée (1992) reported that in some case studies (for instance Wigley II at Fr=0.3 and Fr=0.4) it was not possible to stabilize the experiments, and therefore no results were obtained in those cases. This gives an idea on the complexity of carrying out this type of experiments, and the uncertainty that might be involved (Park et al. 2015). Moreover the added resistance in waves is a second-order force and usually quite small when compared to first-order forces, which makes very complicate to measure and separate both quantities. In fact, quite large models are needed to be able to measure second order forces, which is not the case of the referred experiments where the model length was 3 m. Since experimental data might have significant uncertainties, such results must be analyzed qualitatively rather than quantitatively.

<span id='cite-_Ref433611942'></span>[[#_Ref433611942|Figure 10]] shows a snapshot of the wave pattern generated by the Wigley III hull at Fr=0.2 advancing against a monochromatic wave with a wave length equal to the length of the hull. <span id='cite-_Ref511119114'></span>[[#_Ref511119114|Figure 11]] provides the time evolution of the resistance in still water and in wave for the same case. The added resistance can be obtained from the difference between the still water resistance and the mean value of the resistance in wave. <span id='cite-_Ref430014337'></span>[[#_Ref430014337|Figure 12]] shows the hull geometries and mesh sizes used.

<span id='cite-_Ref430155216'></span>[[#_Ref430155216|Figure 13]]-16 show the heave and pitch RAOs, as well as the dimensionless added resistance <math display="inline">{R}_{aw}^{\ast }</math>, versus the dimensionless wavelength for the four modified Wigley hull models.

Overall the results obtained with the FD-SLI and the FEM-SUPG schemes are quite similar to each other. Only in a few cases slight deviations are observed. Regarding the experimental results general trends are well captured, although in some cases some scattering is observed. Numerical results for RAOS and added resistance in waves follow reasonably well the experimental data although the goodness of the fitting depends on each case. Resonance frequencies are recovered well although peak values of added resistance are often underestimated.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_464519160-image10.jpeg|600px]] </div>

<span id='_Ref433611942'></span><div id="_Toc429661892" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 10: Snapshot of wave pattern for Wigley I. Fr=0.2. λ<sub>w</sub>/L=1.</span></div>

<div id="_Ref511119114" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
{| style="width: 61%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;width: 61%;"| [[Image:Draft_Garcia-Espinosa_464519160-image11.jpeg|354px]]  [[Image:Draft_Garcia-Espinosa_464519160-image12.jpeg|354px]] 
|}
<span style="text-align: center; font-size: 75%;">Figure 11: Dimensionless time values of second-order-corrected longitudinal force for Wigley I. Fr=0.2. λ<sub>w</sub>/L=1.</span></div>


<div id="_Ref430014337" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_464519160-image13.jpeg|600px]]
<span style="text-align: center; font-size: 75%;">Figure 12: Mesh details and modified Wigley hulls view.</span></div>


<div id="_Ref430155216" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_464519160-image14.jpeg|600px]]
<span style="text-align: center; font-size: 75%;">Figure 13: RAOs and </span> <math display="inline">{\mathit{\boldsymbol{R}}}_{\mathit{\boldsymbol{aw}}}^{\mathit{\boldsymbol{\ast }}}</math><span style="text-align: center; font-size: 75%;">for Wigley I hull.</span></div>


<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_464519160-image15.jpeg|360px]] </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 14: RAOs and </span> <math display="inline">{\mathit{\boldsymbol{R}}}_{\mathit{\boldsymbol{aw}}}^{\mathit{\boldsymbol{\ast }}}</math><span style="text-align: center; font-size: 75%;"> for Wigley II hull.</span></div>


<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_464519160-image16.jpeg|600px]] </div>

<span id='_Ref428881649'></span><div id="_Toc429661891" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 15: RAOs and </span> <math display="inline">{\mathit{\boldsymbol{R}}}_{\mathit{\boldsymbol{aw}}}^{\mathit{\boldsymbol{\ast }}}</math><span style="text-align: center; font-size: 75%;"> for Wigley III hull.</span></div>


<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_464519160-image17.jpeg|600px]] </div>

<div id="_Ref430155204" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 16: RAOs and </span> <math display="inline">{\mathit{\boldsymbol{R}}}_{\mathit{\boldsymbol{aw}}}^{\mathit{\boldsymbol{\ast }}}</math><span style="text-align: center; font-size: 75%;"> for Wigley IV hull.</span></div>

===6. Summary and conclusion===

A FEM based model to compute the wave making resistance in still water and added resistance in waves has been proposed. The proposed method solves the governing equations in the time-domain, using unstructured meshes, and is not constrained to linearized free surface boundary conditions. The mathematical model is based on potential and Stokes’ perturbation method, which makes it computationally attractive compared to CFD solvers. Two different numerical schemes have been proposed to solve the free surface boundary conditions: the first one is based on a third order FD-SLI scheme using three upwind plus one downwind points; and the second one is based on the FEM-SUPG stabilization.

The proposed schemes have been verified comparing against analytical solutions available for a moving pressure patch, and validated against experimental results by comparing wave making resistance and wave profiles. The verification and validation analyses have proved that fairly good results are obtained by the FD-SLI and FEM-SUPG schemes.

Regarding the use of linearized convective terms, the following conclusions are made: the NK linearization might be valid for low Froude numbers (lower than 0.15), but it is not advisable for higher velocities, even for slender bodies like the Wigley hull; wave making resistance predictions using the DB linearization are quite close to the NL approach.

The numerical results of added resistance in waves have been validated against experimental data available for four modified Wigley hulls at different Froude numbers and for a number of monochromatic waves. Added resistance highly depend on the wave induced heave and pitch movements, requiring a good calibration of the viscous damping is mandatory. The added resistance in waves is a second-order force and quite small when compared to first-order forces, which makes it very complicate to measure experimentally. In fact quite large models are needed to be able to measure second order forces, which is not the case of the refereed experiments where the model length was 3 m. Hence experimental data has significant uncertainties and its results must be analyzed qualitatively rather than quantitatively. Taking this into account, it can be concluded that the present methodology has been capable of estimating fairly well heave and pitch RAOs, and the added resistance trends.

===Acknowledgements===

This research work was partially supported by the X-SHEAKS project (ENE2014-59194-C2-1-R) of the Ministerio de Economía y Competitividad (Spain).

===REFERENCES===

Afshar, M. (2014). Towards predicting the added resistance of slow ships in waves. Ph.D. Thesis, Technical University of Denmark.

Bai, W., & Teng, B. (2013). Simulations of second-order wave interaction with fixed and floating structures in time domain. ''Ocean Engineering, 74'', 168-177.

Belibassakis, K., Gerostathis, T., Kostas, K., Politis, C., Kaklis, P., Ginnis, A., & Feurer, C. (2013). A BEM isogeometric method for the ship wave-resistance problem. ''Ocean Engineering, 60'', 53-67.

Brooks, A., & Hughes, T. (1982). Streamline upwind/Petrov-Galerkin formulations for convection dominated floaws with particular emphasis on the incompressible Navier-Stokes equations. ''Computer Methods in Applied Mechanics and Engineering, 32'', 199-259.

Donea, J., & Huerta, A. (2003). ''Finite element methods for flow problems.'' John Wiley & Sons, Ltd. doi:10.1002/0470013826.

el Moctar, O., Sigmund, S., & Schelling, T. (2015). Numerical and experimentla analysis of added resistance of ships in waves. ''Proceedings of the ASME 34th International Conference on Ocean, Offshore and Artic Engineering OMAE 2015.'' St. John´s, NL, Canada.

Hizir, O., Kim, M., Turan, O., Day, A., & Incecik, A. (2018). Numerical studies on non-linearity of added resistance and ship motions of KVLCC2 in short and long waves. ''International Journal of Naval Architecture and Ocean Engineering'', 1-11.

ITTC84. (1984). ''Report of Resistance Comitee.'' Göteborg: Proceeding of the 17th International Towing Tank Conference.

Joncquez, S. (2009). Second-order forces and moments acting on ships in waves. Ph.D. Thesis, Technical University of Denmark.

Journée, J. (1992). ''Experiments and calculations on 4 Wigley hull forms in head waves.'' Delft University of Technology. Report 0909.

Ju, S. (1983). ''Study of total viscous drag for the Wigley parabolic ship form.'' IIHR Report no. 261.

Kim, M., Hizir, O., Turan, O., & Incecik, A. (2017a). Numerical studies on added resistance and motions of KVLCC2 in head seas for various ship speeds. ''Ocean Engineering, 140'', 466-476.

Kim, M., Hizir, O., Turan, O., Day, S., & Incecik, A. (2017b). Estimation of added resistance and ship speed loss in a seaway. ''Ocean Engineering, 141'', 465-476.

Kim, Y.-C., Kim, K.-S., Kim, J., Kim, Y., Park, I.-R., & Jang, Y.-H. (2017c). Analysis of added resistance and seakeeping responses in head sea conditions for low-speed full ships using URANS approach. ''International Journal of Naval Architecture and Ocean Engineering, 9'', 614-654.

Liu, S., Papanikolaou, A., & Zaraphonitis, G. (2011). Prediction of added resistance in waves. ''Ocean Engineering, 38'', 614-650.

Maruo, H. (1957). The excess resistance in rough seas. ''International Shipbuilding progress, 4''(35).

Maruo, H. (1960). The drift of a body floating on waves. ''Journal of Ship Research., 4''(3), 1-10.

Maruo, H. (1963). Resistance in waves. ''60th Anniversary Series.'' ''8'', pp. 67-102. The Society of Naval Architects of Japan.

Newman, J., & Poole, F. (1962). ''The wave resistance of a moving pressure distribution in a canal.'' David Taylor Model Basin, Report 1619.

Oñate, E., García-Espinosa, J., Idelsohn, S., & Serván-Camas, B. (2018). Ship hydrodynamics. In ''Encyclopedia of Computational Mechanics Second Edition.'' John Wiley & Sons, Ltd. doi:10.1002/9781119176817.

Ozdemir, Y., & Barlas, B. (2017). Numerical study of ship motions and added resistance in regular incident waves of KVLCC2 model. ''International Journal of Naval Architecture and Ocean Engineering, 9'', 149e159.

Park, D.-M., Kim, Y., Seo, M.-G., & Lee, J. (2016). Study on added resistance of a tanker in head waves at different drafts. ''Ocean Engineering, 111'', 569-581.

Park, D.-M., Lee, J., & Kim, Y. (2015). Uncertainty analysis for added resistance experiment o fKVLCC2 ship. ''Ocean Engineering, 95'', 143-156.

Sadat-Hosseini, H., Wu, P.-C., Carrica, P., Kim, H., Toda, Y., & Stern, F. (2013). CFD verification and validation of added resistance and motions of KVLCC2 with fixed and free surge in short and long head waves. ''Ocean Engineering, 59'', 240-273.

Seo, M.-G., Park, D.-M., Yang, K.-K., & Kim, Y. (2013). Comparative study of ship added resistance in waves. ''Ocean Engineering, 73'', 1-15.

Servan-Camas, B. (2016). A time-domian finite element method for seakeeping and wave resistance problems. Doctoral Thesis, ETSI Navales, Universidad Politecnica de Madrid. [http://oa.upm.es/39794/1/BORJA_SERVAN_CAMAS.pdf http://oa.upm.es/39794/1/BORJA_SERVAN_CAMAS.pdf]

Servan-Camas, B., & Garcia-Espinosa, J. (2013). Accelerated 3D multi-body seakeeping simulation using unstructured finite elements. ''Journal of Computational Physics, 252'', 382-403.

Sigmund, S., & el Moctar, O. (2018). Numerical and experimental investigation of added resistance of<br/>different ship types in short and long waves. ''Ocean Engineering, 147'', 51-67.

Shearer, J., & Cross, J. (1965). The experimental determination of the components of ship resistance for a mathematical model. ''Transactions of teh Royal Institution of Naval Architecs, 107'', 459-473.

Tezdogan, T., Demirel, Y., Kellett, P., Khorasanchi, M., & Incecik, A. (2015). Full-scale unsteady RANS CFD simulations of ship behaviour and performance in head seas due to slow steaming. ''Ocean Engineering, 97'', 186-206.

Zienkiewicz, O., Taylor, R., & Nithiarasu, P. (2005a). ''The finite element method for fluid dynamics.'' Elsevier butterworth-Heinemann.

Zienkiewicz, O., Taylor, R., & Zhu, J. (2005b). ''The finite element metho: its basis and fundamentals   (sixth edition).'' Elsevier butterworth-Heinemann.