==Ship hydrodynamics==

​

'''Eugenio  Oñate¹, Julio García¹ and Sergio R. Idelsohn²'''

​

{|

|-

|(¹)  International Center for Numerical Methods in Engineering (CIMNE)

|-

| Universitat Politécnica de Catalunya (UPC), Gran Capitán, s/n, 

|-

| 08034 Barcelona, Spain

|-

| (²)  CIMNE and International Center for Computational Methods in Engineering (CIMEC)

|-

| Universidad Nacional del Litoral and CONICET, Santa Fe, Argentina

|}

-->

​

==Abstract==

​

This chapter presents an overview of some computational methods for analysis of ship hydrodynamics problems. Attention is focused on the  description of a stabilized finite element formulation derived via a finite calculus procedure. Both arbitrary Lagrangian-Eulerian (ALE) and fully Lagrangian forms are presented. Details of the treatment of the free-surface waves and the interaction between the ship structure and the sea water are also given. Examples of application to a variety of ship hydrodynamics problems are shown.

​

'''Keywords''' Ship hydrodynamics, finite element method, finite calculus, fluid-structure interaction

​

==1 INTRODUCTION==

​

Accurate prediction of the hydrodynamic forces on a ship in motion is of paramount importance in ship design. The water resistance at a certain speed determines the required engine power and thereby the fuel consumption. Minimization of the hydrodynamic forces is therefore an important issue in ship hull design. Further, excitation of a wave pattern by ship motion not only induces wave resistance but may also limit the speed in the vicinity of the shore for environmental reasons, which must also be taken into account in ship design.

​

The usual simplification in ship hydrodynamics design is to separately consider the performance of the ship in still water and its behaviour in open sea. Hydrodynamic optimisation of a hull primarily requires the calculation of the resistance in a calm sea and the open sea effects  are generally taken into account as a wave added resistance.

​

The resistance of a ship in still water can be considered as the sum of several contributions: a viscous resistance associated with the generation of boundary layers, the wave resistance, the air resistance on the superstructure and the induced resistance related to the generation of lift forces.

​

Wave resistance in practical cases amounts to <math display="inline">10</math> to <math display="inline">60%</math> of the total resistance of a ship in still water (Raven, 1996 <span id="citeF-84"></span>[[#cite-84|[84]]]). It increases very rapidly at high speeds dominating the viscous component for high speed ships.  Furthermore, wave resistance is very sensitive to the hull form design and easily affected by small shape modifications. For all these reasons, the  possibility to predict and reduce the wave resistance is an important target.

​

The prediction of the wave pattern and the wave resistance of a ship has  challenged mathematicians and hydrodynamicists for over a century. The Boundary Element Method (BEM) is the basis of many computational algorithms developed in  past years. Here the flow problem is solved using a simple potential model.  BEM methods, termed  by  hydrodynamicists as Panel Methods may be classified into two categories. The first one uses the Kelvin wave source as the elementary singularity. The main advantage of such scheme is the automatic satisfaction of the radiation  condition. The theoretical background of this method was reviewed by Wehausen (1970) <span id="citeF-103"></span>[[#cite-103|[103]]], while computational aspects can be found in Soding (1996) <span id="citeF-87"></span>[[#cite-87|[87]]] and Jenson and Soding (1989) <span id="citeF-58"></span>[[#cite-58|[58]]]. The second class of BEM schemes uses the Rankine source as the elementary singularity. This procedure, first presented by Dawson (1977) <span id="citeF-29"></span>[[#cite-29|[29]]], has been widely applied in practice and many improvements have been addressed to account for the nonlinear wave effects. Among these, a succesful example is the Rankine Panel Method (Xia, 1986 <span id="citeF-106"></span>[[#cite-106|[106]]]; Jenson and Soding, 1989 <span id="citeF-58"></span>[[#cite-58|[58]]]; Nakos and Sclavounos, 1990 <span id="citeF-67"></span>[[#cite-67|[67]]]).

​

In addition to the important developments in potential flow panel methods for practical ship hydrodynamics analysis during the period 1960-1980, much research in the second half of the twentieth century was oriented towards the introduction of viscosity in the CFD analysis. In the 1960's the viscous flow research was mainly focused in 2D boundary layer theory and by the end of the decade several methods for arbitrary pressure gradients were available. This research continued to solve the 3D case during the following decade and an evaluation of the capability of the new methods to predict ship wave resistance was carried out at different workshops (Bai and McCarthy, 1979 <span id="citeF-7"></span>[[#cite-7|[7]]]; Larsson, 1981 <span id="citeF-61"></span>[[#cite-61|[61]]]; Noblesse and McCarthy, 1983 <span id="citeF-69"></span>[[#cite-69|[69]]]). Here application to some well specified test cases were reported and numerical and experimental results compared acceptable well for most part of the boundary layer along the hull, while wrong results were obtained near the stern. This prompted additional research and by the end of the 1980's a number of numerical procedures for solving the full viscous flow equation accounting for simple turbulence modes based on Reynolds averaged Navier-Stokes (RANS) equations were available. Considerable improvements for predicting the stern flow were reported in subsequent workshops organized in the 1990's (Kim and Lucas, 1990 <span id="citeF-59"></span>[[#cite-59|[59]]]; Reed ''et al.'' , 1990<span id="citeF-85"></span>[[#cite-85|[85]]]; Beck ''et al.'', 1993 <span id="citeF-8"></span>[[#cite-8|[8]]]; Raven, 1996 <span id="citeF-84"></span>[[#cite-84|[84]]]; Soding, 1996 <span id="citeF-87"></span>[[#cite-87|[87]]]; Janson and Larsson, 1996 <span id="citeF-57"></span>[[#cite-57|[57]]]; Alessandrini and Delhommeau, 1996 <span id="citeF-1"></span>[[#cite-1|[1]]]; Miyata, 1996 <span id="citeF-67"></span>[[#cite-67|[67]]], Löhner ''et al.'', 1998 <span id="citeF-64"></span>[[#cite-64|[64]]]). A good review of the status of CFD in ship hydrodynamics in the last part of the 20th century can be found in Larsson ''et al.'' (1998) <span id="citeF-62"></span>[[#cite-62|[62]]].

​

Independently of the flow equations used, the free-surface boundary condition has been solved in different manners. The exact free surface condition is nonlinear and several linearizations have been proposed (Baba and Takekuma, 1975 <span id="citeF-6"></span>[[#cite-6|[6]]]; Newmann, 1976 <span id="citeF-68"></span>[[#cite-68|[68]]]; Idelsohn ''et al.'', 1999 <span id="citeF-52"></span>[[#cite-52|[52]]]). Some of them use a fixed domain and others a moving one. An alternative is to solve the full nonlinear free surface equation on a reference surface which does not necessarily coincides with the free surface itself. In this way the updating of the surface mesh is minimized and sometimes is not even necessary.

​

The solution of the free-surface equation in a bounded domain brings in the necessity of a radiation condition to eliminate spurious waves. A way to introduce this condition was proposed by Dawson (1977) <span id="citeF-29"></span>[[#cite-29|[29]]] who used a finite difference (FD) formula based in four upwind points to evaluate the first derivatives appearing in the free-surface equation. This method became very popular and this is probably the main reason why a large majority of codes predicting the wave resistance of ships use FD methods on structured meshes (Larsson ''et al.'', 1998 <span id="citeF-62"></span>[[#cite-62|[62]]]).

​

Indeed the 1990's were a decade of considerable progress in CFD methods for ship hydrodynamics and the most important breakthrough was perhaps the coupled solution of the free-surface equation with the fluid flow equations. Here a number of viscous and inviscid solutions for the surface ship wave problem using finite element (FE) and finite volume (FV) methods with non structured grids were reported (Farmer ''et al.'', 1993 <span id="citeF-35"></span>[[#cite-35|[35]]], Hino ''et al.'', 1993 <span id="citeF-41"></span>[[#cite-41|[41]]];   Luo ''et al.'', 1995 <span id="citeF-66"></span>[[#cite-66|[66]]]; Storti ''et al.'', 1998a,b <span id="citeF-90"></span>[[#cite-90|[90]]-<span id="citeF-91"></span>[[#cite-91|91]]]; García ''et al.'', 1998 <span id="citeF-36"></span>[[#cite-36|[36]]]; García, 1999 <span id="citeF-37"></span>[[#cite-37|[37]]]; Alessandrini and Delhommeau, 1999 <span id="citeF-2"></span>[[#cite-2|[2]]]; Idelsohn ''et al.'', 1999 <span id="citeF-52"></span>[[#cite-52|[52]]]; Löhner ''et al.'', 1998 <span id="citeF-64"></span>[[#cite-64|[64]]], 1999 <span id="citeF-65"></span>[[#cite-65|[65]]]).

​

The current challenges in CFD research for ship hydrodynamics focus in the development of robust (stable) and computationaly efficient numerical methods able to capture the different scales involved in the analysis of practical ship hydrodynamics situations. Wave resistance coefficients for modern ship design are needed for a wide range of speeds and here the accurate prediction of the wave pattern and the hull pressure distribution at low speed (say below Froude number (<math display="inline">Fn = 0,2</math>) ) are still major challenges. Great difficulties also exist in the computation of the viscous resistance which requires very fine grids in the vicinity of the hull, resulting in overall meshes involving (at least) some <math display="inline">10^7-10^9</math> elements. Other relevant problems are the prediction of the wake details and the propeller-hull interaction. Fine meshing and advanced turbulence models are crucial for the realistic solution of these problems. Indeed the use of unstructured meshes is essential for problems involving complex shapes.

​

A different class of ship hydrodynamic problems require the modelling of breaking waves or the prediction of water inside the hull (green water) due to large amplitude waves typical of  sea keeping problems. Here Lagrangian flow methods where the motion of each flow particle is individually tracked using techniques developed for (incompressible) solid mechanics are a promising trend for solving a wide class of ship hydrodynamics problems.

​

The content of the chapter is structured as follows. In the next section the standard Navier-Stokes equations for an incompressible viscous flow are presented. The equations are formulated in an arbitrary Lagrangian-Eulerian (ALE) description allowing the independent motion of the mesh nodes from that of the fluid particles. Details of the problems posed  by the free-surface wave boundary condition are given. The difficulties encountered in the numerical solution of the fluid flow and the free-surface equations, namely the unstabilities induced by the convective terms and the limits in the approximation introduced by the incompressibility constraint are explained. A new procedure for deriving stabilized numerical methods for this type of problems based on the so called ''finite calculus'' (FIC) formulation is presented. The FIC method is based in redefining the standard governing equations in fluid mechanics by expressing the balance laws in a domain of ''finite size''. This introduces additional terms in the differential equations of the infinitessimal theory which are essential to derive stabilized numerical schemes (Oñate, 1998 <span id="citeF-70"></span>[[#cite-70|[70]]], 2000 <span id="citeF-71"></span>[[#cite-71|[71]]], 2004 <span id="citeF-72"></span>[[#cite-72|[72]]]). We present here a stabilized finite element method using equal-order linear interpolation for the velocity and the pressure variables. Both monolithic and fractional step time integration procedures are  described. A method for solving the coupled fluid-structure interaction problem induced by the motion of the ship due to the hydrodynamic forces is also presented.  A mesh moving algorithm  for updating the position of the free-surface nodes during the ship motion is given.

​

In the last part of the chapter the Lagrangian formulation for fluid flow analysis is presented as a particular case of the ALE form. As above mentioned the Lagrangian description has particular advantages for tracking the displacement of  the fluid particles in flows where large motions of the fluid surface occur. One of the advantages of the Lagrangian approach is that the convective terms dissapear in the governing equations of the fluid. In return, the updating of the mesh at almost every time step is now a necessity and efficient mesh generation algorithms must be used (Idelsohn ''et al.'', 2002 <span id="citeF-53"></span>[[#cite-53|[53]]], 2003a,b,c <span id="citeF-54"></span>[[#cite-54|[54]]-<span id="citeF-56"></span>[[#cite-56|56]]]).

​

The examples show the efficiency of the ALE and fully Lagrangian formulations to solve a variety of ship hydrodynamics problems.

​

==2 THE NAVIER-STOKES EQUATIONS FOR INCOMPRESSIBLE FLOWS. ALE FORMULATION==

​

===2.1 Momentum and mass conservation equations===

​

The  Navier-Stokes equations for an incompressible fluid in a domain <math display="inline">\Omega </math> can be written in an arbitrary Lagrangian-Eulerian form as

​

''Momentum''

<span id="eq-1"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\displaystyle{\rho \left({\partial u_i \over \partial t}+v_j {\partial u_i \over \partial x_j}+ u_i{\partial u_j \over \partial x_j} \right)+ {\partial p \over \partial x_i} - {\partial s_{ij} \over \partial x_j}-b_i =0 } \qquad \hbox{in }\Omega  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (1)

|}

​

''Mass conservation''

<span id="eq-2"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\displaystyle{{\partial u_i \over \partial x_i}=0} \qquad \hbox{in }\Omega  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (2)

|}

​

In Eq.([[#eq-1|1]]) <math display="inline">u_i</math> is the velocity along the ith global reference axis, <math display="inline">u_i^m</math> is the velocity of  the moving mesh nodes, <math display="inline">v_i=u_i-u_i^m</math> is the relative velocity between the fluid and the moving mesh nodes, <math display="inline">\rho </math> is the (constant) density of the fluid, <math display="inline">b_i</math> are body forces, <math display="inline">t</math> is the time, <math display="inline">p</math> is the pressure and <math display="inline">s_{ij}</math> are the viscous stresses related to the viscosity <math display="inline">\mu </math> by the standard expression

<span id="eq-3"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>s_{ij}=2\mu \left(\varepsilon _{ij} - {1\over 3}\delta _{ij}{\partial u_k \over \partial x_k}\right) </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (3)

|}

​

where <math display="inline">\delta _{ij}</math> is the Kronecker delta and the strain rates <math display="inline">\varepsilon _{ij}</math> are

<span id="eq-4"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\varepsilon _{ij}= {1\over 2}\left({\partial u_i \over \partial x_j}+ {\partial u_j \over \partial x_i}\right) </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (4)

|}

​

The volumetric strain rate terms in the Eqs.([[#eq-1|1]]) and ([[#eq-3|3]]) can be neglected following Eq.([[#eq-2|2]]). We will however retain these terms as they are useful for the derivation of the stabilized formulation in Section [[#5 FIC EQUATIONS FOR VISCOUS INCOMPRESSIBLE FLOW. ALE FORMULATION|5]].

​

Eqs.([[#eq-1|1]]) and ([[#eq-3|3]]) are completed with the boundary conditions

<span id="eq-5a"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>n_j\sigma _{ij}-t_i =0 \quad \hbox{on }\Gamma _t</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (5a)

|}

<span id="eq-5b"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>u_j - u_j^p =0 \quad \hbox{on }\Gamma _u</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (5b)

|}

​

where <math display="inline">\sigma _{ij}=s_{ij}- p\delta _{ij}</math> are the total stresses, <math display="inline">n_j</math> are the component of the unit normal vector to the boundary and <math display="inline">t_i</math> and <math display="inline">u_j^p</math> are prescribed tractions and velocities on the Neumann and Dirichlet boundaries <math display="inline">\Gamma _t</math> and <math display="inline">\Gamma _u</math>, respectively where <math display="inline">\Gamma = \Gamma _t \cup \Gamma _u</math> is the boundary of the analysis domain <math display="inline">\Omega </math>. The boundary conditions are completed with the initial condition <math display="inline">u_j =u_j^0</math> for <math display="inline">t=t_0</math>.

​

In above equations <math display="inline">i,j=1,n_d</math> where <math display="inline">n_d</math> is the number of space dimension (i.e. <math display="inline">n_d =3</math> for 3D). Also, throughout the text the summation convention for repeated indexes is assumed unless specified otherwise.

​

===2.2 Free-surface boundary conditions===

​

The boundary conditions ([[#eq-5a|5a]]) on the surface tractions can be written in local normal and tangential axes as

<span id="eq-6"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\begin{array}{l}\sigma _n -t_n = 0\\ \sigma _{t_j}- t_{g_j}=0  \qquad j=1,2 \quad \hbox{for 3D} \end{array} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (6)

|}

​

where <math display="inline">\sigma _n</math> and <math display="inline">\sigma _{t_j}</math> are the normal and tangential stresses, respectively and <math display="inline">t_n</math> and <math display="inline">t_{g_j}</math> are the prescribed normal and tangential tractions, respectively.

​

On the free boundary we have to ensure at all times that: 1) the pressure (which approximates the normal stres) equals the atmospheric pressure and the tangential tractions are zero unless specified otherwise, and 2) material particles of the fluid belong to the free-surface.

​

The condition on the pressure is simply especified as

<span id="eq-7"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>p=p_a </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (7)

|}

​

where <math display="inline">p_a</math> is the atmospheric pressure (usually given a zero value).

​

The condition on the tangential tractions is satisfied by setting <math display="inline">t_{g_i} = 0</math> in Eq.([[#eq-6|6]]). This is automatically accounted for by the ''natural boundary conditions'' in the weak form of the momentum equations (Zienkiewicz and Taylor Vol. 3, 2000 <span id="citeF-107"></span>[[#cite-107|[107]]]).

​

The condition on the material particles is expressed (for steady state conditions) as

<span id="eq-8"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>u_in_i =0 \qquad \hbox{or } \qquad {\boldsymbol u}^T {\boldsymbol n}=0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (8)

|}

​

i.e. the velocity vector is tangent to the free-surface. An alternative form of Eq.([[#eq-8|8]]) can be found by noting that the normal vector has the following components

<span id="eq-9"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\boldsymbol n}=\left[-{\partial \beta  \over \partial x_1},1\right]^T\qquad \hbox{for 2D }\quad  \hbox{ and }\quad {\boldsymbol n}=\left[-{\partial \beta  \over \partial x_1}, -{\partial \beta  \over \partial x_2},1\right]^T \qquad \hbox{for 3D}  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (9)

|}

​

In Eq.([[#eq-9|9]]) <math display="inline">\beta = x_3 - x_3^{ref}</math> is the free-surface elevation measured in the direction of the vertical coordinate <math display="inline">x_3</math> relative to some previously known surface, which we shall refer to as the ''reference surface'' (Figure [[#img-2|2]]). This surface may  be horizontal (i.e. the undisturbed water surface) or may be simply a previously computed surface.

​

Introducing Eq.([[#eq-9|9]]) into ([[#eq-8|8]]) gives (for 3D)

<span id="eq-10"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>u_i {\partial \beta  \over \partial x_i}-u_3 =0 \quad , \quad i=1,2 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (10)

|}

​

Eq.([[#eq-10|10]]) is generalized for the transient 3D case as

<span id="eq-11"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\displaystyle{{\partial \beta  \over \partial t}+u_i {\partial \beta  \over \partial x_i}-u_3 =0}\quad , \quad i=1,2 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (11)

|}<br/>

​

We observe that <math display="inline">\beta </math> obeys ''a pure convection'' equation with <math display="inline">u_3</math> playing the role of a (non linear) source term. The solution of Eq.([[#eq-11|11]]) with standard Galerkin FEM, or centred FD or FV methods will therefore suffer from numerical instabilities and some kind of stabilization is needed in order to obtain a physically meaningful solution. A  method to solve Eq.([[#eq-11|11]]), very popular in the context of the potential flow formulation, was introduced by Dawson (1977) <span id="citeF-29"></span>[[#cite-29|[29]]] using a four point FD upwind operator to evaluate the first derivatives of <math display="inline">\beta </math> on regular grids. Dawson's method has been extended by different authors to solve a many ship hydrodynamics problems (Raven, 1996 <span id="citeF-84"></span>[[#cite-84|[84]]]; Larsson ''et al.'', 1998 <span id="citeF-62"></span>[[#cite-62|[62]]]; Idelsohn ''et al.'', 1999 <span id="citeF-52"></span>[[#cite-52|[52]]]).

​

Solution of Eq.([[#eq-11|11]]) is strongly coupled with that of the fluid flow equations. The solution of the whole problem is highly non linear due to the pressence of the unknown velocities in Eq.([[#eq-11|11]]) and also to the fact that the free-surface position defining the new boundary conditions is also unknown. A number of iterative schemes have been developed for the solution of the non linear surface wave problem (Idelsohn ''et al.'', 1999 <span id="citeF-52"></span>[[#cite-52|[52]]]). They all basically involve solving Eq.([[#eq-11|11]]) for the new free-surface height <math display="inline">\beta </math>, for fixed values of the velocity field computed from the fluid solver in a previous iteration within each time increment. At this stage two procedures are possible, either the position of the free-surface is updated after each iteration and this becomes the new reference surface, or else an equivalent pressure of value <math display="inline">p=p_a + g (\beta -  \beta _{ref})</math>, where <math display="inline">g</math> is the gravity constant, is applied at the current reference surface as a boundary condition in the next flow iteration. The first option might require the regeneration of a new mesh, whereas the second one is less accurate but computationaly cheaper. Hence a compromise between the two alternatives is usually chosen in practice. The iterative process continues until a converged solution is found for the velocity, the pressure and the free-surface height at each time step. Details of the computational process are described in a next section.

​

An alternative method to treat the free-surface equation is based in the volume of fluid (VOF) technique (Hirt and Nichols, 1981 <span id="citeF-43"></span>[[#cite-43|[43]]]). In the VOF method the free-surface position is defined as the interface between two fluids interacting with each other, where the efect of one fluid on the other is very small (i.e. the water and the surrounding air). An interface function which takes the values 0 and 1 for each of the two fluids is transported with the fluid velocity using a time dependent advection equation. Early applications of the VOF in the context of the stabilized FEM were reported by Codina ''et al.'', 1994 <span id="citeF-26"></span>[[#cite-26|[26]]] and Tezduyar ''et al.'', 1998 <span id="citeF-102"></span>[[#cite-102|[102]]].  Examples of application of the VOF to ship hydrodynamics problems can be found in Tajima and Yabe, 1999 <span id="citeF-93"></span>[[#cite-93|[93]]] and Azcueta ''et al.'', 1999 <span id="citeF-5"></span>[[#cite-5|[5]]].

​

==3 ABOUT THE FINITE ELEMENT SOLUTION OF THE NAVIER-STOKES EQUATIONS==

​

The development of efficient and robust numerical methods for incompressible flow problems has been a subject of intensive research in last decades. Much effort has been spent in developing the so called ''stabilized'' numerical methods overcoming the two main sources of instability in incompressible flow analysis, namely those originated by the high values of the convective terms  and those induced by the difficulty in satisfying the incompressibility constraint.

​

The solution of above problems in the context of the finite element method (FEM) has been attempted in a number of ways. The underdiffusive character of the Galerkin FEM for high convection flows (which incidentaly also occurs for centred FD and FV methods) has been corrected by adding some kind of artificial viscosity terms to the standard Galerkin equations. A good review of such approach can be found in  Zienkiewicz and Taylor, Vol. 3 2000 <span id="citeF-107"></span>[[#cite-107|[107]]] and Donea and Huerta, 2003 <span id="citeF-31"></span>[[#cite-31|[31]]].

​

A popular way to overcome the problems with the incompressibility constraint is by introducing a pseudo-compressibility in the flow and using implicit and explicit algorithms developed for this kind of problems, such as artificial compressibility schemes (Chorin, 1967A <span id="citeF-14"></span>[[#cite-14|[14]]] popular way to overcome the problems with the incompressibility constraint is by introducing a pseudo-compressibility in the flow and using implicit and explicit algorithms developed for this kind of problems, such as artificial compressibility schemes (Chorin, 1967 <span id="citeF-14"></span>[[#cite-14|[14]]]; Farmer ''et al.'', 1993 <span id="citeF-35"></span>[[#cite-35|[35]]]; Peraire ''et al.'', 1994 <span id="citeF-82"></span>[[#cite-82|[82]]]; Briley ''et al.'', 1995 <span id="citeF-10"></span>[[#cite-10|[10]]]; Sheng ''et al.'', 1996 <span id="citeF-86"></span>[[#cite-86|[86]]])  and preconditioning techniques (Idelsohn ''et al.'', 1995 <span id="citeF-51"></span>[[#cite-51|[51]]]). Other FEM schemes with good stabilization properties for the convective and incompressibility terms  are based in Petrov-Galerkin (PG) techniques. The background of PG methods are the non-centred (upwind) schemes for computing the first derivatives of the convective operator in FD and FV methods. More recently a general class of Galerkin FEM has been developed where the standard Galerkin variational form is extended with adequate residual-based terms in order to achieve a stabilized numerical scheme (Codina, 1998 <span id="citeF-16"></span>[[#cite-16|[16]]], 2000 <span id="citeF-17"></span>[[#cite-17|[17]]]). Among the many FEM of this kind  we can name the Streamline Upwind Petrov Galerkin (SUPG) method (Hughes and Brooks, 1979 <span id="citeF-45"></span>[[#cite-45|[45]]]; Brooks and Hughes, 1982 <span id="citeF-11"></span>[[#cite-11|[11]]]; Tezduyar and Hughes, 1983 <span id="citeF-97"></span>[[#cite-97|[97]]]; Hughes and Tezduyar, 1984 <span id="citeF-46"></span>[[#cite-46|[46]]]; Hughes and Mallet, 1986 <span id="citeF-48"></span>[[#cite-48|[48]]]; Idelsohn ''et al.'', 1995 <span id="citeF-51"></span>[[#cite-51|[51]]]; Storti ''et al.'', 1995 <span id="citeF-88"></span>[[#cite-88|[88]]], 1997 <span id="citeF-89"></span>[[#cite-89|[89]]]; Cruchaga and Oñate, 1997 <span id="citeF-27"></span>[[#cite-27|[27]]], 1999 <span id="citeF-28"></span>[[#cite-28|[28]]]), the Galerkin Least Square (GLS) method (Hughes ''et al.'', 1989 <span id="citeF-50"></span>[[#cite-50|[50]]]; Tezduyar, 1991 <span id="citeF-94"></span>[[#cite-94|[94]]]; Tezduyar ''et al.'', 1992a <span id="citeF-99"></span>[[#cite-99|[99]]]), the Taylor-Galerkin method (Donea, 1984 <span id="citeF-30"></span>[[#cite-30|[30]]]), the Characteristic Galerkin method (Douglas and Russell, 1982 <span id="citeF-32"></span>[[#cite-32|[32]]]; Pironneau, 1982 <span id="citeF-83"></span>[[#cite-83|[83]]]; Löhner ''et al.'', 1984 <span id="citeF-63"></span>[[#cite-63|[63]]]) and its variant the characteristic Based Split (CBS) method (Zienkiewicz and Codina, 1995 <span id="citeF-108"></span>[[#cite-108|[108]]]; Codina ''et al.'', 1998 <span id="citeF-23"></span>[[#cite-23|[23]]]; Codina and Zienkiewicz, 2002 <span id="citeF-25"></span>[[#cite-25|[25]]]), pressure gradient operator methods (Codina and Blasco, 1997 <span id="citeF-22"></span>[[#cite-22|[22]]], 2000 <span id="citeF-24"></span>[[#cite-24|[24]]])  and the Subgrid Scale (SGS) method (Hughes, 1995 <span id="citeF-44"></span>[[#cite-44|[44]]]; Brezzi ''et al.'', 1997 <span id="citeF-9"></span>[[#cite-9|[9]]]; Codina, 2000 <span id="citeF-17"></span>[[#cite-17|[17]]], 2002 <span id="citeF-21"></span>[[#cite-21|[21]]]).

​

In this work a stabilized FEM for incompressible flows is derived taking as the starting point the modified governing equations of the flow problem formulated via a finite calculus (FIC) approach. The FIC method is based in invoking the balance of fluxes in a domain of finite size. This introduces naturally additional terms in the classical differential equations of infinitessimal fluid mechanics which are a function of the balance domain dimensions. The new terms in the modified governing equations provide naturally the necessary stabilization to the standard Galerkin finite element method.

​

In the next section, the main concepts of the FIC approch are introduced via a simple 1D convection-diffusion model problem. Then the  FIC formulation of the fluid flow equations and the free-surface  wave equations using the finite element method  are presented.

​

==4 BASIC CONCEPTS OF THE FINITE CALCULUS (FIC) METHOD==

​

We will consider a convection-diffusion problem in a 1D domain <math display="inline">\Omega </math> of length <math display="inline">L</math>. The equation of balance of fluxes in a subdomain of size <math display="inline">d</math> belonging to <math display="inline">\Omega </math> (Figure [[#img-1|1]]) is written as

<span id="eq-12"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>q_A - q_B=0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (12)

|}

​

where <math display="inline">q_A</math> and <math display="inline">q_B</math> are the incoming and outgoing fluxes at points <math display="inline">A</math> and <math display="inline">B</math>, respectively. The flux <math display="inline">q</math> includes both convective and diffusive terms; i.e. <math display="inline">q=u\phi - k{d\phi \over dx}</math>, where <math display="inline">\phi </math> is the transported variable, <math display="inline">u</math> is the velocity and <math display="inline">k</math> is the diffusitivity of the material.

​

We express now the fluxes <math display="inline">q_A</math> and <math display="inline">q_B</math> in terms of the flux at an arbitrary point <math display="inline">C</math> within the balance domain (Figure  [[#img-1|1]]). Expanding <math display="inline">q_A</math> and <math display="inline">q_B</math> in Taylor series around point <math display="inline">C</math> up to second order terms gives

<span id="eq-13"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>q_A= q_C - d_1 \frac{dq}{d x}\vert _C+ \frac{d^2_1}{2}\frac{d^2q}{dx^2}\vert _C + O(d^3_1)\quad ,\quad  q_B= q_C + d_2 \frac{dq}{d x}\vert _C+\frac{d^2_2}{2}\frac{d^2q}{dx^2}\vert _C + O(d^3_2) </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (13)

|}

​

Substituting Eq.([[#eq-13|13]]) into Eq.([[#eq-12|12]]) gives after simplification

<span id="eq-14"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\frac{dq}{dx}-\underline{\frac{h}{2} \frac{d^2q}{dx^2}}=0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (14)

|}

​

where <math display="inline">h=d_1-d_2</math> and all derivatives are computed at point <math display="inline">C</math>.

​

Standard calculus theory assumes  that the domain <math display="inline">d</math> is of infinitessimal size and the resulting balance equation is simply <math display="inline">{dq\over dx}=0</math>. We will relax this assumption and allow the balance domain to have a ''finite size''. The new balance equation ([[#eq-14|14]]) incorporates now the underlined term which introduces the ''characteristic length'' <math display="inline">h</math>. Obviously, accounting for higher order terms in Eq.([[#eq-13|13]]) would lead to new terms in Eq.([[#eq-14|14]]) involving higher powers of <math display="inline">h</math>.

​

Distance <math display="inline">h</math> in Eq.([[#eq-14|14]]) can be interpreted as a free parameter depending on the location of point <math display="inline">C</math> (note that <math display="inline">h=0</math> for <math display="inline">d_1=d_2</math>). However, the fact that Eq.([[#eq-14|14]]) is the exact balance equation (up to second order terms) for any 1D domain of finite size and that the position of point <math display="inline">C</math> is arbitrary, can be used to derive numerical schemes with enhanced properties, simply by computing the characteristic length parameter using an adequate “optimality” rule.

​

Consider, for instance, the modified equation ([[#eq-14|14]]) applied to the convection-diffusion problem. Neglecting third order derivatives of <math display="inline">\phi </math>, Eq.([[#eq-14|14]]) can be written as

<span id="eq-15"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>-u \frac{d\phi }{dx}+\left(k+\frac{u h}{2}\right)\frac{d^2\phi }{dx^2}=0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (15)

|}

​

We see that the FIC procedure introduces ''naturally'' an additional diffusion term into the standard convection-diffusion equation. This is the basis of the popular “artificial diffusion” method (Hirsch, 1988, 1990 <span id="citeF-42"></span>[[#cite-42|[42]]]). The characteristic length <math display="inline">h</math> is typically expressed as a function of the cell or element dimensions. The optimal or critical value of <math display="inline">h</math> can be computed from numerical stability conditions such as obtaining a physically meaningful solution, or even obtaining “exact” nodal values (Zienkiewicz and Taylor, 2000 <span id="citeF-107"></span>[[#cite-107|[107]]]; Oñate and Manzan, 1999 <span id="citeF-76"></span>[[#cite-76|[76]]], 2000 <span id="citeF-77"></span>[[#cite-77|[77]]]; Oñate, 2004 <span id="citeF-72"></span>[[#cite-72|[72]]]).

​

Equation ([[#eq-13|13]]) can be extended to account for source and time effects. The full FIC equation for the transient convection-diffusion problem can be written in compact form as

<span id="eq-16"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>r- \underline{\frac{h}{2} \frac{dr}{dx}}- \underline{{\delta \over 2} {\partial r \over \partial t}}=0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (16)

|}

<span id="eq-17"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>r = -\left(\frac{d\phi }{dt}+u \frac{d\phi }{dx}\right)+ \frac{d}{dx}\left(k \frac{d\phi }{dx}\right)+ Q </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (17)

|}

​

where <math>Q</math> is the external source and <math>\delta </math> is a time stabilization parameter (Oñate, 1998 <span id="citeF-70"></span>[[#cite-70|[70]]]; Oñate and Manzan, 1999 <span id="citeF-76"></span>[[#cite-76|[76]]]). For consistency a FIC form of the Neumann boundary condition should be used. This is obtained by invoking balance of fluxes in a domain of finite size next to the boundary <math>\Gamma _q</math> where the flux is prescribed to a value <math>\bar q</math>. The modified FIC boundary condition is

<span id="eq-18"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>k \frac{d\phi }{dx}+\bar q- \underline{\frac{h}{2}r}=0 \quad \hbox{at } \Gamma _q </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (18)

|}

​

The definition of the problem is completed with the standard Dirichlet condition prescribing the value of <math display="inline">\phi </math> at the boundary <math display="inline">\Gamma _\phi </math> and the initial conditions.

​

==5 FIC EQUATIONS FOR VISCOUS INCOMPRESSIBLE FLOW. ALE FORMULATION==

​

The starting point are the FIC equations for a viscous incompressible fluid. For simplicity we will neglect the time stabilization term, as this is not relevant for the purposes of this work. The equations are written as (Oñate, 1998 <span id="citeF-70"></span>[[#cite-70|[70]]], 2000 <span id="citeF-71"></span>[[#cite-71|[71]]]; Oñate ''et al.'', 2002 <span id="citeF-79"></span>[[#cite-79|[79]]])

​

''Momentum''

<span id="eq-19"></span>

{| 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{{1\over 2} h_j{\partial r_{m_i} \over \partial x_j}}=0  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (19)

|}

​

''Mass balance''

<span id="eq-20"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\left({\partial u_k \over \partial x_k}\right)- \underline{{h_j\over 2} {\partial  \over \partial x_j} \left({\partial u_k \over \partial x_k}\right)}=0  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (20)

|}

​

where

<span id="eq-21"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>r_{m_i} = \rho \left({\partial u_i \over \partial t}+v_j{\partial u_i \over \partial x_j}+u_i {\partial u_j \over \partial x_j}\right)+ {\partial p \over \partial x_i}- {\partial s_{ij} \over \partial x_j}-b_i \qquad i,j = 1, n_d </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (21)

|}

​

and all the terms have been defined in Section [[#2.1 Momentum and mass conservation equations|2.1.]]

​

The Neumann boundary conditions for the FIC formulation are (Oñate, 1998 <span id="citeF-70"></span>[[#cite-70|[70]]], 2000 <span id="citeF-71"></span>[[#cite-71|[71]]])

<span id="eq-22"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>n_j \sigma _{ij} -t_i + \underline{{1\over 2} h_j n_j r_{m_i}}=0 \quad \hbox{on }\Gamma _t </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (22)

|}

​

The Dirichlet and initial boundary conditions are the standard ones given in Section [[#2.1 Momentum and mass conservation equations|2.1.]]

​

The <math display="inline">h_i's</math> in above equations are characteristic lengths of the domain where balance of momentum and mass is enforced. In Eq.([[#eq-22|22]]) these lengths define the domain where equilibrium of boundary tractions is established. The sign in front the <math display="inline">h_i</math> term in Eq.([[#eq-22|22]]) is consistent with the definition of <math display="inline">r_{m_i}</math> in Eq.([[#eq-21|21]]).

​

Eqs.([[#eq-19|19]])-([[#eq-22|22]]) are the starting  point for deriving stabilized FEM for solving the incompressible Navier-Stokes equations using equal-order interpolation for all variables.

​

===5.1 Stabilized integral forms===

​

From Eqs.([[#eq-19|19]]) and ([[#eq-3|3]]) it can be obtained

<span id="eq-23"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{h_j\over 2} {\partial  \over \partial x_j} \left({\partial u_k \over \partial x_k}\right)\simeq   \sum \limits _{i=1}^{n_d}\tau _i     {\partial r_{m_i} \over \partial x_i} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (23)

|}

​

where

<span id="eq-24"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\tau _i = \left({8\mu \over 3h_i^2}+{2\rho u_i\over h_i}\right)^{-1} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (24)

|}

​

Substituting Eq.([[#eq-23|23]]) into Eq.([[#eq-20|20]]) leads to the following stabilized mass balance equation

<span id="eq-25></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\partial u_k \over \partial x_k}- \sum \limits _{i=1}^{n_d} \tau _i {\partial r_{m_i} \over \partial x_i}=0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (25)

|}

​

The <math display="inline">\tau _i</math>'s in Eq.([[#eq-23|23]]) are ''intrinsic time parameters per unit mass''. Similar parameters  also appear in other stabilized formulations (Hughes and Mallet, 1986a <span id="citeF-48"></span>[[#cite-48|[48]]]; Tezduyar, 1991 <span id="citeF-94"></span>[[#cite-94|[94]]], 2001 <span id="citeF-95"></span>[[#cite-95|[95]]]; Codina, 2002 <span id="citeF-21"></span>[[#cite-21|[21]]]). Note that  the parameters of Eq.([[#eq-24|24]]) emerge naturally form the FIC formation and take the values of <math display="inline">\tau _i =\displaystyle{3h_i^2\over 8\mu }</math> and <math display="inline">\tau _i =\displaystyle{h_i\over 2\rho u_i}</math> for the viscous limit (Stokes flow) and the inviscid  limit (Euler flow), respectively.

​

The weighted residual form of the governing equations (Eqs.([[#eq-19|19]]), ([[#eq-22|22]]) and ([[#eq-25|25]])) is

<span id="eq-26"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\int _\Omega \delta u_i \left[r_{m_i} - {h_j\over 2} {\partial r_{m_i} \over \partial x_j}\right]+ \int _{\Gamma _t} \delta u_i (\sigma _{ij} n_j - t_i +{ h_j\over 2} n_j r_{m_i}) d\Gamma =0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (26)

|}

<span id="eq-27"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\int _\Omega q \left[{\partial u_k \over \partial x_k}- \sum \limits _{i=1}^{n_d} \tau _i {\partial r_{m_i} \over \partial x_i}\right]d\Omega =0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (27)

|}

​

where <math display="inline">\delta u_i</math> and <math display="inline">q</math> are arbitrary weighting functions representing virtual velocity and virtual pressure fields, respectively. Integrating by parts above equations leads to

<span id="eq-28"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\int _\Omega \left[\delta u_i\rho \left({\partial u_i \over \partial t}+v_j {\partial {u_i} \over \partial x_j}\right)+ \delta \varepsilon _{ij}(\tau _{ij}- \delta _{ij}p)\right]d\Omega -  \int _{\Omega } \delta u_i b_i d\Omega -\int _{\Gamma _t} \delta u_i t_id\Omega + </math>

|-

| style="text-align: center;" | <math> + \sum \limits _e \int _{\Omega ^e} {h_j\over 2}{\partial \delta u_i \over \partial x_j} r_{m_i} d\Omega =0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (28)

|}

<span id="eq-29"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\int _\Omega q r_d \,d\Omega + \int _\Omega \left[\sum \limits _{i=1}^{n_d} \tau _i {\partial q \over \partial x_i} r_{m_i}\right]d\Omega =0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (29)

|}

​

The fourth integral in Eq.([[#eq-28|28]]) is computed as a sum of the element contributions to allow for discontinuities in the derivatives of <math display="inline">r_{m_i}</math> along the element interfaces. As usual <math display="inline">\delta \varepsilon _{ij} ={1\over 2} \left({\partial \delta u_i \over \partial x_j}+{\partial \delta u_j \over \partial x_i}\right)</math>. In Eq.([[#eq-28|28]]) we have neglected the volumetric strain rate term, whereas in the derivation of Eq.([[#eq-29|29]]) we have assumed that <math display="inline">r_{m_i}</math> is negligible on the boundaries.

​

===5.2 Convective and pressure gradient projections===

​

The convective and pressure gradient projections <math display="inline">c_i</math> and <math display="inline">\pi _i</math> are defined as

<span id="eq-30"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>c_i = r_{m_i} -\rho v_j {\partial u_i \over \partial x_j}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (30)

|}

<span id="eq-31"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math> \pi _i = r_{m_i} - {\partial p \over \partial x_i} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (31)

|}

​

We can now express <math display="inline">r_{m_i}</math> in Eqs.([[#eq-28|28]]) and ([[#eq-29|29]]) in terms of <math display="inline">c_i</math> and <math display="inline">\pi _i</math>, respectively which become additional variables. The system of integral equations is now augmented in the necessary number of equations by imposing that the residuals <math display="inline">r_{m_i}</math> vanish (in average sense) for both forms given by Eqs.([[#eq-30|30]]) and ([[#eq-31|31]]). The final system of integral equations is:

<span id="eq-32"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\int _\Omega \left[\delta u_i\rho \left({\partial u_i \over \partial t}+v_j {\partial {u_i} \over \partial x_j}\right)+ \delta \varepsilon _{ij}(\tau _{ij}- \delta _{ij}p)\right] d\Omega -  \int _{\Omega } \delta u_i b_i d\Omega - \int _{\Gamma _t} \delta u_i t_id\Omega +</math>

|-

| style="text-align: center;" | <math>  + \sum \limits _e \int _{\Omega ^e} {h_k\over 2}{\partial (\delta u_i) \over \partial x_k} \left(\rho v_j {\partial {u_i} \over \partial x_j} + c_i\right)d\Omega =0</math>

|-

| style="text-align: center;" | <math> \int _\Omega q \left({\partial u_k \over \partial x_k}\right) d\Omega +  \int _\Omega \sum \limits _{i=1}^{n_d} \tau _i {\partial q \over \partial x_i} \left({\partial p \over \partial x_i}+\pi _i\right)d\Omega =0</math>

|-

| style="text-align: center;" | <math> \int _\Omega \delta c_i \rho \left(\rho v_j {\partial {u_i} \over \partial x_j} + c_i\right)d\Omega =0 \qquad \hbox{ no sum in }i</math>

|-

| style="text-align: center;" | <math> \int _\Omega \delta \pi _i \tau _i \left({\partial p \over \partial x_i}+\pi _i\right)d\Omega =0\qquad \hbox{ no sum in }i </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (32)

|}

​

where <math display="inline">\delta c_i</math> and <math display="inline">\delta \pi _i</math> are appropriate weighting functions and the <math display="inline">\rho </math> and <math display="inline">\tau _i</math> terms are introduced in the last two equations for convenience. As usual <math display="inline">i,j,k=1,n_d</math>.

​

===5.3 Stabilized FIC equations for the free-surface wave condition===

​

We will derive next the FIC equation for the water wave surface.

​

Let us consider a 2D free-surface wave problem. Figure 3 shows a typical free-surface segment line <math display="inline">AB</math>. The average vertical velocity for the segment <math display="inline">\bar u_2</math> is defined as

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\bar u_2 ={u_2^A + u_2^B\over 2}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (33a)

|}

​

where <math display="inline">u_2^A</math> and <math display="inline">u_2^B</math> are the vertical velocities of the end points of the segment <math display="inline">A</math> and <math display="inline">B</math>.

​

The average vertical velocity <math display="inline">\bar u_2</math> can be computed from the wave heights at points <math display="inline">A</math> and <math display="inline">B</math> as

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\bar u_2 ={x_2^B -x_2^A\over \bar \delta }</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (33b)

|}

​

where <math display="inline">\bar \delta </math> is the time which a material particle takes to travel from point <math display="inline">A</math> to point <math display="inline">B</math> at the average velocity <math display="inline">\bar u_2</math> and <math display="inline">x_2 \equiv \beta </math> is the free-surface wave height. Equaling Eq.(33a) and (33b) gives

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{u_2^A + u_2^B\over 2}= {x_2^B -x_2^A\over \bar \delta }</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (33c)

|}

​

We can now express the vertical velocities and the wave height at points <math display="inline">A</math> and <math display="inline">B</math> in terms of values at an arbitrary point <math display="inline">C</math> as

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>u_2^A = u_2 (x_1^C -d_1, t^C-\delta _1) = u_2^C -d_1 {\partial \beta  \over \partial x_1}-\delta _1 {\partial \beta  \over \partial t}+ O (d_1^2,t_1^2) </math>

|-

| style="text-align: center;" | <math> u_2^B = u_2 (x_1^C +d_2, t^C +\delta _2) = u_2^C +d_2 {\partial \beta  \over \partial x_2}+\delta _2 {\partial \beta  \over \partial t}+ O (d_2^2, t_2^2)</math>

|-

| style="text-align: center;" | <math> x_2^A= x_2 (x_1^C -d_1, t^C-\delta _1) = x_2^C -d_1 {\partial \beta  \over \partial x_1}-\delta _1 {\partial \beta  \over \partial t}+{d_1^2\over 2} {\partial ^2 \beta \over \partial x_1^2}+ </math>

|-

| style="text-align: center;" | <math> +\, d_1\delta _1 {\partial ^2 \beta \over \partial x_1\partial t}+{\delta _1^2\over 2}{\partial ^2 \beta \over \partial t^2} + O (d_1^3, \delta _1^3)</math>

|-

| style="text-align: center;" | <math> x_2^B= x_2 (x_1^C +d_2, t^C +\delta _2) = x_2^C +d_2 {\partial \beta  \over \partial x_1}+\delta _2 {\partial \beta  \over \partial t}+{d_2^2\over 2} {\partial ^2 \beta \over \partial x_1^2}+ </math>

|-

| style="text-align: center;" | <math> +\, d_2\delta _2 {\partial ^2 \beta \over \partial x_1\partial t}+{\delta _2^2\over 2}{\partial ^2 \beta \over \partial t^2} + O (d_2^3, \delta _2^3) </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (34)

|}

​

where all the derivatives are computed at the arbitrary point <math display="inline">C</math>.

​

Substituting Eqs.(34) into (33c) and noting that <math display="inline">d\simeq u_1^C\bar \delta </math>, with <math display="inline">\bar \delta =(\delta _1+\delta _2)</math>, <math display="inline">d_1 \simeq u_1^C\delta _1</math> and <math display="inline">d_2 \simeq u_1^C\delta _2</math> (Figure 3)  and that the position of point <math display="inline">C</math> is arbitrary, gives the FIC equation for  the free-surface height (neglecting high order terms) as

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>r_\beta - \underline{{h\over 2}{\partial r_\beta  \over \partial x_1}}{h\over 2}{\partial r_\beta  \over \partial x_1} - \underline{{\delta \over 2}{\partial r_\beta  \over \partial t}}{\delta \over 2}{\partial r_\beta  \over \partial t}=0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (35)

|}

​

with

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>r_\beta = {\partial \beta  \over \partial t}+u_1 {\partial \beta  \over \partial x_1}-u_2 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (36)

|}

​

where <math display="inline">h=(d_1 -d_2)</math> and <math display="inline">\delta = (\delta _1 - \delta _2)</math> are space and time stabilization parameters. The standard infinitesimal form of the free-surface wave condition is obtained by making <math display="inline">h=\delta =0</math> in Eq.(35) giving

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\partial \beta  \over \partial t}+u_1 {\partial \beta  \over \partial x_1}-u_2=0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (37)

|}

​

which coincides with equation (11) for the 2D case.

​

A simpler FIC expression can be derived from Eq.(35) by retaining the second order space term only. This gives

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\partial \beta  \over \partial t}+u_1 {\partial \beta  \over \partial x_1}-{u_1h\over 2} {\partial ^2\beta \over \partial x_1^2}-u_2 = 0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (38)

|}

​

This can be interpreted as the addition of an artificial diffusion term where <math display="inline">{u_1h\over 2}</math> plays the role of the new balancing diffusion coefficient.

​

In the following we will use the 3D ALE form of Eq.(35) neglecting the time stabilization term, given by

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>r_\beta - \underline{{h_j\over 2}{\partial r_\beta  \over \partial x_j}}{h_j\over 2}{\partial r_\beta  \over \partial x_j}=0\qquad \hbox{with }\quad  r_\beta = {\partial \beta  \over \partial t}+ v_j {\partial \beta  \over \partial x_j}-v_3 \quad , \quad j=1,2 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (39)

|}

​

where, as usual, <math display="inline">v_i</math> are the relative velocities between the fluid and the moving mesh nodes.

​

==6 FINITE ELEMENT DISCRETIZATION==

​

===6.1 Discretization of the fluid flow equations===

​

We will choose <math display="inline">C^\circ </math> continuous linear interpolations of the velocities, the pressure, the convection projections <math display="inline">c_i</math> and the pressure gradient projections <math display="inline">\pi _i</math> over three node triangles (2D) and four node tetrahedra (3D). The interpolations are written as

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\begin{array}{l}u_i = \displaystyle \sum \limits _{j=1}^n N_j \bar u_i^j \quad , \quad p = \sum \limits _{j=1}^n N_j \bar p^j\\ c_i = \displaystyle \sum \limits _{j=1}^n N_j \bar c_i^j \quad , \quad \pi _i = \sum \limits _{j=1}^n N_j \bar \pi _i^j \end{array} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (40)

|}

​

where <math display="inline">n=3</math> (4) for triangles (tetrahedra), <math display="inline">\bar {(\cdot )}</math> denotes nodal variables and <math display="inline">N_j</math> are the linear shape functions (Zienkiewicz and Taylor, Vol 1 2000).

​

Substituting the approximations (40) into Eqs.(32) and choosing the Galerking form with <math display="inline">\delta u_i =q=\delta c_i=\delta \pi _i =N_i</math> leads to following system of discretized equations

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\boldsymbol M}\dot{\bar {\boldsymbol u}} + ({\boldsymbol A}+{\boldsymbol K}+\hat{\boldsymbol K}) \bar {\boldsymbol u} - {\boldsymbol G}\bar {\boldsymbol p}+{\boldsymbol C}\bar {\boldsymbol c}={\boldsymbol f}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (41a)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\boldsymbol G}^T \bar {\boldsymbol u} + {\boldsymbol L}\bar {\boldsymbol p}+{\boldsymbol Q}\bar {\boldsymbol \pi }={\boldsymbol 0}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (41b)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\qquad \qquad \hat{\boldsymbol C}\bar {\boldsymbol u}+ {\boldsymbol M}\bar {\boldsymbol c}={\boldsymbol 0} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (41c)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\boldsymbol Q}^T \bar {\boldsymbol p} + \hat{\boldsymbol M}\bar {\boldsymbol \pi }={\boldsymbol 0}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (41d)

|}

​

The element contributions are given by (for 2D problems)

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\displaystyle {\boldsymbol M}_{ij}\!\!=\!\! \int _{\Omega ^e} \rho {\boldsymbol N}_i{\boldsymbol N}_j d\Omega \quad ,\quad {\boldsymbol A}_{ij}= \!\int _{\Omega ^e}\! {\boldsymbol N}_i \rho v_k {\partial {\boldsymbol N}_j \over \partial x_k} d\Omega \quad ,\quad \displaystyle {\boldsymbol K}_{ij}= \! \int _{\Omega ^e}\! {\boldsymbol B}_i^T {\boldsymbol D}{\boldsymbol B}_j d\Omega </math>

|-

| style="text-align: center;" |