Published in Encyclopedia of Computational Mechanics, Encyclopedia of Computational Mechanics, E. Stein, R. de Borst and T.J.R. Hughes (Eds.), John Wiley & Sons Ltd, Vol. 3, Chapter 18, pp. 579 - 607, 2004
DOI: 10.1002/9781119176817

## 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 ${\textstyle 10}$ to ${\textstyle 60\%}$ of the total resistance of a ship in still water (Raven, 1996 ). 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) , while computational aspects can be found in Soding (1996)  and Jenson and Soding (1989) . The second class of BEM schemes uses the Rankine source as the elementary singularity. This procedure, first presented by Dawson (1977) , 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 ; Jenson and Soding, 1989 ; Nakos and Sclavounos, 1990 ).

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 ; Larsson, 1981 ; Noblesse and McCarthy, 1983 ). 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 ; Reed et al. , 1990; Beck et al., 1993 ; Raven, 1996 ; Soding, 1996 ; Janson and Larsson, 1996 ; Alessandrini and Delhommeau, 1996 ; Miyata, 1996 , Löhner et al., 1998 ). 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) .

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 ; Newmann, 1976 ; Idelsohn et al., 1999 ). 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)  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 ).

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 , Hino et al., 1993 ; Luo et al., 1995 ; Storti et al., 1998a,b [90-91]; García et al., 1998 ; García, 1999 ; Alessandrini and Delhommeau, 1999 ; Idelsohn et al., 1999 ; Löhner et al., 1998 , 1999 ).

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 (${\textstyle Fn=0,2}$) ) 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 ${\textstyle 10^{7}-10^{9}}$ 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 , 2000 , 2004 ). 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 , 2003a,b,c [54-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 ${\textstyle \Omega }$ can be written in an arbitrary Lagrangian-Eulerian form as

Momentum

 $\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$
(1)

Mass conservation

 $\displaystyle {{\partial u_{i} \over \partial x_{i}}=0}\qquad {\hbox{in }}\Omega$
(2)

In Eq.(1) ${\textstyle u_{i}}$ is the velocity along the ith global reference axis, ${\textstyle u_{i}^{m}}$ is the velocity of the moving mesh nodes, ${\textstyle v_{i}=u_{i}-u_{i}^{m}}$ is the relative velocity between the fluid and the moving mesh nodes, ${\textstyle \rho }$ is the (constant) density of the fluid, ${\textstyle b_{i}}$ are body forces, ${\textstyle t}$ is the time, ${\textstyle p}$ is the pressure and ${\textstyle s_{ij}}$ are the viscous stresses related to the viscosity ${\textstyle \mu }$ by the standard expression

 $s_{ij}=2\mu \left(\varepsilon _{ij}-{1 \over 3}\delta _{ij}{\partial u_{k} \over \partial x_{k}}\right)$
(3)

where ${\textstyle \delta _{ij}}$ is the Kronecker delta and the strain rates ${\textstyle \varepsilon _{ij}}$ are

 $\varepsilon _{ij}={1 \over 2}\left({\partial u_{i} \over \partial x_{j}}+{\partial u_{j} \over \partial x_{i}}\right)$
(4)

The volumetric strain rate terms in the Eqs.(1) and (3) can be neglected following Eq.(2). We will however retain these terms as they are useful for the derivation of the stabilized formulation in Section 5.

Eqs.(1) and (3) are completed with the boundary conditions

 $n_{j}\sigma _{ij}-t_{i}=0\quad {\hbox{on }}\Gamma _{t}$
(5a)

 $u_{j}-u_{j}^{p}=0\quad {\hbox{on }}\Gamma _{u}$
(5b)

where ${\textstyle \sigma _{ij}=s_{ij}-p\delta _{ij}}$ are the total stresses, ${\textstyle n_{j}}$ are the component of the unit normal vector to the boundary and ${\textstyle t_{i}}$ and ${\textstyle u_{j}^{p}}$ are prescribed tractions and velocities on the Neumann and Dirichlet boundaries ${\textstyle \Gamma _{t}}$ and ${\textstyle \Gamma _{u}}$, respectively where ${\textstyle \Gamma =\Gamma _{t}\cup \Gamma _{u}}$ is the boundary of the analysis domain ${\textstyle \Omega }$. The boundary conditions are completed with the initial condition ${\textstyle u_{j}=u_{j}^{0}}$ for ${\textstyle t=t_{0}}$.

In above equations ${\textstyle i,j=1,n_{d}}$ where ${\textstyle n_{d}}$ is the number of space dimension (i.e. ${\textstyle n_{d}=3}$ 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 (5a) on the surface tractions can be written in local normal and tangential axes as

 ${\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}}$
(6)

where ${\textstyle \sigma _{n}}$ and ${\textstyle \sigma _{t_{j}}}$ are the normal and tangential stresses, respectively and ${\textstyle t_{n}}$ and ${\textstyle t_{g_{j}}}$ 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

 $p=p_{a}$
(7)

where ${\textstyle p_{a}}$ is the atmospheric pressure (usually given a zero value).

The condition on the tangential tractions is satisfied by setting ${\textstyle t_{g_{i}}=0}$ in Eq.(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 ).

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

 $u_{i}n_{i}=0\qquad {\hbox{or }}\qquad {\boldsymbol {u}}^{T}{\boldsymbol {n}}=0$
(8)

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

 ${\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}}$
(9)

In Eq.(9) ${\textstyle \beta =x_{3}-x_{3}^{ref}}$ is the free-surface elevation measured in the direction of the vertical coordinate ${\textstyle x_{3}}$ relative to some previously known surface, which we shall refer to as the reference surface (Figure 2). This surface may be horizontal (i.e. the undisturbed water surface) or may be simply a previously computed surface.

Introducing Eq.(9) into (8) gives (for 3D)

 $u_{i}{\partial \beta \over \partial x_{i}}-u_{3}=0\quad ,\quad i=1,2$
(10)

Eq.(10) is generalized for the transient 3D case as

 $\displaystyle {{\partial \beta \over \partial t}+u_{i}{\partial \beta \over \partial x_{i}}-u_{3}=0}\quad ,\quad i=1,2$
(11)

We observe that ${\textstyle \beta }$ obeys a pure convection equation with ${\textstyle u_{3}}$ playing the role of a (non linear) source term. The solution of Eq.(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.(11), very popular in the context of the potential flow formulation, was introduced by Dawson (1977)  using a four point FD upwind operator to evaluate the first derivatives of ${\textstyle \beta }$ on regular grids. Dawson's method has been extended by different authors to solve a many ship hydrodynamics problems (Raven, 1996 ; Larsson et al., 1998 ; Idelsohn et al., 1999 ).

Solution of Eq.(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.(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 ). They all basically involve solving Eq.(11) for the new free-surface height ${\textstyle \beta }$, 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 ${\textstyle p=p_{a}+g(\beta -\beta _{ref})}$, where ${\textstyle g}$ 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 ). 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  and Tezduyar et al., 1998 . Examples of application of the VOF to ship hydrodynamics problems can be found in Tajima and Yabe, 1999  and Azcueta et al., 1999 .

## 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  and Donea and Huerta, 2003 .

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 ) 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 ; Farmer et al., 1993 ; Peraire et al., 1994 ; Briley et al., 1995 ; Sheng et al., 1996 ) and preconditioning techniques (Idelsohn et al., 1995 ). 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 , 2000 ). Among the many FEM of this kind we can name the Streamline Upwind Petrov Galerkin (SUPG) method (Hughes and Brooks, 1979 ; Brooks and Hughes, 1982 ; Tezduyar and Hughes, 1983 ; Hughes and Tezduyar, 1984 ; Hughes and Mallet, 1986 ; Idelsohn et al., 1995 ; Storti et al., 1995 , 1997 ; Cruchaga and Oñate, 1997 , 1999 ), the Galerkin Least Square (GLS) method (Hughes et al., 1989 ; Tezduyar, 1991 ; Tezduyar et al., 1992a ), the Taylor-Galerkin method (Donea, 1984 ), the Characteristic Galerkin method (Douglas and Russell, 1982 ; Pironneau, 1982 ; Löhner et al., 1984 ) and its variant the characteristic Based Split (CBS) method (Zienkiewicz and Codina, 1995 ; Codina et al., 1998 ; Codina and Zienkiewicz, 2002 ), pressure gradient operator methods (Codina and Blasco, 1997 , 2000 ) and the Subgrid Scale (SGS) method (Hughes, 1995 ; Brezzi et al., 1997 ; Codina, 2000 , 2002 ).

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 ${\textstyle \Omega }$ of length ${\textstyle L}$. The equation of balance of fluxes in a subdomain of size ${\textstyle d}$ belonging to ${\textstyle \Omega }$ (Figure 1) is written as

 $q_{A}-q_{B}=0$
(12)

where ${\textstyle q_{A}}$ and ${\textstyle q_{B}}$ are the incoming and outgoing fluxes at points ${\textstyle A}$ and ${\textstyle B}$, respectively. The flux ${\textstyle q}$ includes both convective and diffusive terms; i.e. ${\textstyle q=u\phi -k{d\phi \over dx}}$, where ${\textstyle \phi }$ is the transported variable, ${\textstyle u}$ is the velocity and ${\textstyle k}$ is the diffusitivity of the material.

We express now the fluxes ${\textstyle q_{A}}$ and ${\textstyle q_{B}}$ in terms of the flux at an arbitrary point ${\textstyle C}$ within the balance domain (Figure 1). Expanding ${\textstyle q_{A}}$ and ${\textstyle q_{B}}$ in Taylor series around point ${\textstyle C}$ up to second order terms gives

 $q_{A}=q_{C}-d_{1}{\frac {dq}{dx}}\vert _{C}+{\frac {d_{1}^{2}}{2}}{\frac {d^{2}q}{dx^{2}}}\vert _{C}+O(d_{1}^{3})\quad ,\quad q_{B}=q_{C}+d_{2}{\frac {dq}{dx}}\vert _{C}+{\frac {d_{2}^{2}}{2}}{\frac {d^{2}q}{dx^{2}}}\vert _{C}+O(d_{2}^{3})$
(13)

Substituting Eq.(13) into Eq.(12) gives after simplification

 ${\frac {dq}{dx}}-{\underline {{\frac {h}{2}}{\frac {d^{2}q}{dx^{2}}}}}=0$
(14)

where ${\textstyle h=d_{1}-d_{2}}$ and all derivatives are computed at point ${\textstyle C}$.

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

Distance ${\textstyle h}$ in Eq.(14) can be interpreted as a free parameter depending on the location of point ${\textstyle C}$ (note that ${\textstyle h=0}$ for ${\textstyle d_{1}=d_{2}}$). However, the fact that Eq.(14) is the exact balance equation (up to second order terms) for any 1D domain of finite size and that the position of point ${\textstyle C}$ 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 (14) applied to the convection-diffusion problem. Neglecting third order derivatives of ${\textstyle \phi }$, Eq.(14) can be written as

 $-u{\frac {d\phi }{dx}}+\left(k+{\frac {uh}{2}}\right){\frac {d^{2}\phi }{dx^{2}}}=0$
(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, 1990 ). The characteristic length ${\textstyle h}$ is typically expressed as a function of the cell or element dimensions. The optimal or critical value of ${\textstyle h}$ can be computed from numerical stability conditions such as obtaining a physically meaningful solution, or even obtaining “exact” nodal values (Zienkiewicz and Taylor, 2000 ; Oñate and Manzan, 1999 , 2000 ; Oñate, 2004 ).

Equation (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

 $r-{\underline {{\frac {h}{2}}{\frac {dr}{dx}}}}-{\underline {{\delta \over 2}{\partial r \over \partial t}}}=0$
(16)

 $r=-\left({\frac {d\phi }{dt}}+u{\frac {d\phi }{dx}}\right)+{\frac {d}{dx}}\left(k{\frac {d\phi }{dx}}\right)+Q$
(17)

where $Q$ is the external source and $\delta$ is a time stabilization parameter (Oñate, 1998 ; Oñate and Manzan, 1999 ). 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 $\Gamma _{q}$ where the flux is prescribed to a value ${\bar {q}}$. The modified FIC boundary condition is

 $k{\frac {d\phi }{dx}}+{\bar {q}}-{\underline {{\frac {h}{2}}r}}=0\quad {\hbox{at }}\Gamma _{q}$
(18)

The definition of the problem is completed with the standard Dirichlet condition prescribing the value of ${\textstyle \phi }$ at the boundary ${\textstyle \Gamma _{\phi }}$ 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 , 2000 ; Oñate et al., 2002 )

Momentum

 $r_{m_{i}}-{\underline {{1 \over 2}h_{j}{\partial r_{m_{i}} \over \partial x_{j}}}}=0$
(19)

Mass balance

 $\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$
(20)

where

 $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}$
(21)

and all the terms have been defined in Section 2.1.

The Neumann boundary conditions for the FIC formulation are (Oñate, 1998 , 2000 )

 $n_{j}\sigma _{ij}-t_{i}+{\underline {{1 \over 2}h_{j}n_{j}r_{m_{i}}}}=0\quad {\hbox{on }}\Gamma _{t}$
(22)

The Dirichlet and initial boundary conditions are the standard ones given in Section 2.1.

The ${\textstyle h_{i}'s}$ in above equations are characteristic lengths of the domain where balance of momentum and mass is enforced. In Eq.(22) these lengths define the domain where equilibrium of boundary tractions is established. The sign in front the ${\textstyle h_{i}}$ term in Eq.(22) is consistent with the definition of ${\textstyle r_{m_{i}}}$ in Eq.(21).

Eqs.(19)-(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.(19) and (3) it can be obtained

 ${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}}$
(23)

where

 $\tau _{i}=\left({8\mu \over 3h_{i}^{2}}+{2\rho u_{i} \over h_{i}}\right)^{-1}$
(24)

Substituting Eq.(23) into Eq.(20) leads to the following stabilized mass balance equation

 ${\partial u_{k} \over \partial x_{k}}-\sum \limits _{i=1}^{n_{d}}\tau _{i}{\partial r_{m_{i}} \over \partial x_{i}}=0$
(25)

The ${\textstyle \tau _{i}}$'s in Eq.(23) are intrinsic time parameters per unit mass. Similar parameters also appear in other stabilized formulations (Hughes and Mallet, 1986a ; Tezduyar, 1991 , 2001 ; Codina, 2002 ). Note that the parameters of Eq.(24) emerge naturally form the FIC formation and take the values of ${\textstyle \tau _{i}=\displaystyle {3h_{i}^{2} \over 8\mu }}$ and ${\textstyle \tau _{i}=\displaystyle {h_{i} \over 2\rho u_{i}}}$ for the viscous limit (Stokes flow) and the inviscid limit (Euler flow), respectively.

The weighted residual form of the governing equations (Eqs.(19), (22) and (25)) is

 $\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$
(26)

 $\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$
(27)

where ${\textstyle \delta u_{i}}$ and ${\textstyle q}$ are arbitrary weighting functions representing virtual velocity and virtual pressure fields, respectively. Integrating by parts above equations leads to

 $\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_{i}d\Omega +$ $+\sum \limits _{e}\int _{\Omega ^{e}}{h_{j} \over 2}{\partial \delta u_{i} \over \partial x_{j}}r_{m_{i}}d\Omega =0$
(28)

 $\int _{\Omega }qr_{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$
(29)

The fourth integral in Eq.(28) is computed as a sum of the element contributions to allow for discontinuities in the derivatives of ${\textstyle r_{m_{i}}}$ along the element interfaces. As usual ${\textstyle \delta \varepsilon _{ij}={1 \over 2}\left({\partial \delta u_{i} \over \partial x_{j}}+{\partial \delta u_{j} \over \partial x_{i}}\right)}$. In Eq.(28) we have neglected the volumetric strain rate term, whereas in the derivation of Eq.(29) we have assumed that ${\textstyle r_{m_{i}}}$ is negligible on the boundaries.

### 5.2 Convective and pressure gradient projections

The convective and pressure gradient projections ${\textstyle c_{i}}$ and ${\textstyle \pi _{i}}$ are defined as

 $c_{i}=r_{m_{i}}-\rho v_{j}{\partial u_{i} \over \partial x_{j}}$
(30)

 $\pi _{i}=r_{m_{i}}-{\partial p \over \partial x_{i}}$
(31)

We can now express ${\textstyle r_{m_{i}}}$ in Eqs.(28) and (29) in terms of ${\textstyle c_{i}}$ and ${\textstyle \pi _{i}}$, respectively which become additional variables. The system of integral equations is now augmented in the necessary number of equations by imposing that the residuals ${\textstyle r_{m_{i}}}$ vanish (in average sense) for both forms given by Eqs.(30) and (31). The final system of integral equations is:

 $\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_{i}d\Omega +$ $+\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$ $\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$ $\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$ $\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$
(32)

where ${\textstyle \delta c_{i}}$ and ${\textstyle \delta \pi _{i}}$ are appropriate weighting functions and the ${\textstyle \rho }$ and ${\textstyle \tau _{i}}$ terms are introduced in the last two equations for convenience. As usual ${\textstyle i,j,k=1,n_{d}}$.

### 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 ${\textstyle AB}$. The average vertical velocity for the segment ${\textstyle {\bar {u}}_{2}}$ is defined as

 ${\bar {u}}_{2}={u_{2}^{A}+u_{2}^{B} \over 2}$
(33a)

where ${\textstyle u_{2}^{A}}$ and ${\textstyle u_{2}^{B}}$ are the vertical velocities of the end points of the segment ${\textstyle A}$ and ${\textstyle B}$.

The average vertical velocity ${\textstyle {\bar {u}}_{2}}$ can be computed from the wave heights at points ${\textstyle A}$ and ${\textstyle B}$ as

 ${\bar {u}}_{2}={x_{2}^{B}-x_{2}^{A} \over {\bar {\delta }}}$
(33b)

where ${\textstyle {\bar {\delta }}}$ is the time which a material particle takes to travel from point ${\textstyle A}$ to point ${\textstyle B}$ at the average velocity ${\textstyle {\bar {u}}_{2}}$ and ${\textstyle x_{2}\equiv \beta }$ is the free-surface wave height. Equaling Eq.(33a) and (33b) gives

 ${u_{2}^{A}+u_{2}^{B} \over 2}={x_{2}^{B}-x_{2}^{A} \over {\bar {\delta }}}$
(33c)

We can now express the vertical velocities and the wave height at points ${\textstyle A}$ and ${\textstyle B}$ in terms of values at an arbitrary point ${\textstyle C}$ as

 $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})$ $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})$ $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}}+$ $+\,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})$ $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}}+$ $+\,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})$
(34)

where all the derivatives are computed at the arbitrary point ${\textstyle C}$.

Substituting Eqs.(34) into (33c) and noting that ${\textstyle d\simeq u_{1}^{C}{\bar {\delta }}}$, with ${\textstyle {\bar {\delta }}=(\delta _{1}+\delta _{2})}$, ${\textstyle d_{1}\simeq u_{1}^{C}\delta _{1}}$ and ${\textstyle d_{2}\simeq u_{1}^{C}\delta _{2}}$ (Figure 3) and that the position of point ${\textstyle C}$ is arbitrary, gives the FIC equation for the free-surface height (neglecting high order terms) as

 $r_{\beta }-{\underline {{h \over 2}{\partial r_{\beta } \over \partial x_{1}}}}-{\underline {{\delta \over 2}{\partial r_{\beta } \over \partial t}}}=0$
(35)

with

 $r_{\beta }={\partial \beta \over \partial t}+u_{1}{\partial \beta \over \partial x_{1}}-u_{2}$
(36)

where ${\textstyle h=(d_{1}-d_{2})}$ and ${\textstyle \delta =(\delta _{1}-\delta _{2})}$ are space and time stabilization parameters. The standard infinitesimal form of the free-surface wave condition is obtained by making ${\textstyle h=\delta =0}$ in Eq.(35) giving

 ${\partial \beta \over \partial t}+u_{1}{\partial \beta \over \partial x_{1}}-u_{2}=0$
(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

 ${\partial \beta \over \partial t}+u_{1}{\partial \beta \over \partial x_{1}}-{u_{1}h \over 2}{\partial ^{2}\beta \over \partial x_{1}^{2}}-u_{2}=0$
(38)

This can be interpreted as the addition of an artificial diffusion term where ${\textstyle {u_{1}h \over 2}}$ 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

 $r_{\beta }-{\underline {{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$
(39)

where, as usual, ${\textstyle v_{i}}$ 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 ${\textstyle C^{\circ }}$ continuous linear interpolations of the velocities, the pressure, the convection projections ${\textstyle c_{i}}$ and the pressure gradient projections ${\textstyle \pi _{i}}$ over three node triangles (2D) and four node tetrahedra (3D). The interpolations are written as

 ${\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}}$
(40)

where ${\textstyle n=3}$ (4) for triangles (tetrahedra), ${\textstyle {\bar {(\cdot )}}}$ denotes nodal variables and ${\textstyle N_{j}}$ are the linear shape functions (Zienkiewicz and Taylor, Vol 1 2000 ).

Substituting the approximations (40) into Eqs.(32) and choosing the Galerking form with ${\textstyle \delta u_{i}=q=\delta c_{i}=\delta \pi _{i}=N_{i}}$ leads to following system of discretized equations

 ${\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}}$
(41a)

 ${\boldsymbol {G}}^{T}{\bar {\boldsymbol {u}}}+{\boldsymbol {L}}{\bar {\boldsymbol {p}}}+{\boldsymbol {Q}}{\bar {\boldsymbol {\pi }}}={\boldsymbol {0}}$
(41b)

 $\qquad \qquad {\hat {\boldsymbol {C}}}{\bar {\boldsymbol {u}}}+{\boldsymbol {M}}{\bar {\boldsymbol {c}}}={\boldsymbol {0}}$
(41c)

 ${\boldsymbol {Q}}^{T}{\bar {\boldsymbol {p}}}+{\hat {\boldsymbol {M}}}{\bar {\boldsymbol {\pi }}}={\boldsymbol {0}}$
(41d)

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

 $\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$ ${\hat {\boldsymbol {K}}}_{ij}\!\!=\!\!\int _{\Omega ^{e}}\rho {h_{k}v_{l} \over 2}\displaystyle {\partial {\boldsymbol {N}}_{i} \over \partial x_{k}}\displaystyle {\partial {\boldsymbol {N}}_{j} \over \partial x_{l}}d\Omega \!\!\quad ,\quad \!\!\displaystyle {\boldsymbol {G}}_{ij}=\!\int _{\Omega ^{e}}\!({\boldsymbol {\nabla }}N_{i})N_{j}d\Omega \!\!\quad ,\!\!\quad {\boldsymbol {C}}_{ij}=\!\int _{\Omega ^{e}}\!{h_{k} \over 2}\displaystyle {\partial {\boldsymbol {N}}_{i} \over \partial x_{k}}{\boldsymbol {N}}_{j}d\Omega \,d\Omega$ $\displaystyle L_{ij}\!\!=\!\!\int _{\Omega ^{e}}{\boldsymbol {\nabla }}^{T}N_{i}[\tau ]{\boldsymbol {\nabla }}N_{j}d\Omega ~,~{\hat {\boldsymbol {C}}}_{ij}=\!\int _{\Omega ^{e}}\!{\boldsymbol {N}}_{i}\rho ^{2}v_{k}\displaystyle {\partial {\boldsymbol {N}}_{j} \over \partial x_{k}}d\Omega ~,~{\boldsymbol {\nabla }}=\!\left\{{\begin{array}{c}\displaystyle {\partial ~ \over \partial x_{1}}\\\displaystyle {\partial ~ \over \partial x_{2}}\end{array}}\right\}~,~[\tau ]=\!\left[{\begin{array}{cc}\tau _{1}&0\\0&\tau _{2}\end{array}}\right]$ ${\boldsymbol {Q}}\!\!=\!\![{\boldsymbol {Q}}^{1},{\boldsymbol {Q}}^{2}]\quad ,\quad Q_{ij}^{k}=\int _{\Omega ^{e}}\tau _{k}{\partial N_{i} \over \partial x_{k}}N_{j}d\Omega \quad {\hbox{ no sum in }}k$ (42) ${\hat {\boldsymbol {M}}}\!\!=\!\!\left[{\begin{array}{cc}{\hat {\boldsymbol {M}}}^{1}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\hat {\boldsymbol {M}}}^{2}\end{array}}\right]\quad ,\quad {\hat {M}}_{ij}^{k}=\int _{\Omega ^{e}}\tau _{k}N_{i}N_{j}d\Omega \quad ,\quad {\boldsymbol {N}}_{i}=N_{i}\left[{\begin{array}{cc}1&0\\0&1\end{array}}\right]$ ${\boldsymbol {f}}_{i}\!\!=\!\!\int _{\Omega ^{e}}N_{i}{\boldsymbol {b}}d\Omega +\int _{\Gamma ^{e}}N_{i}{\boldsymbol {t}}d\Gamma \quad ,\quad {\boldsymbol {b}}=[b_{1},b_{2}]^{T},\,{\boldsymbol {t}}=[t_{1},t_{2}]^{T}$

with ${\textstyle i,j=1,n}$ and ${\textstyle k,l=1,n_{d}}$.

In above ${\textstyle {\boldsymbol {B}}_{i}}$ is the standard strain rate matrix and ${\textstyle {\boldsymbol {D}}}$ the deviatoric constitutive matrix (for ${\textstyle \displaystyle {\partial u_{i} \over \partial x_{i}}=0}$). For 2D problems

 ${\boldsymbol {B}}_{i}=\left[{\begin{array}{cc}{\partial N_{i} \over \partial x_{1}}&0\\0&{\partial N_{i} \over \partial x_{2}}\\{\partial N_{i} \over \partial x_{1}}&{\partial N_{i} \over \partial x_{2}}\end{array}}\right]\quad ,\quad {\boldsymbol {D}}=\mu \left[{\begin{array}{ccc}2&0&0\\0&2&0\\0&0&1\end{array}}\right]$
(43)

Note that the stabilization matrix ${\textstyle {\hat {\boldsymbol {K}}}}$ brings in an additional orthotropic diffusivity of value ${\textstyle \rho \displaystyle {h_{k}v_{l} \over 2}}$. Matrices ${\textstyle {\boldsymbol {A}},{\hat {\boldsymbol {K}}}}$ and ${\textstyle {\hat {\boldsymbol {C}}}}$ are dependent on the velocity field. The solution process can be advanced in time in a (quasi-nearly) implicit iterative manner using the following scheme.

Step 1

 ${\bar {\boldsymbol {u}}}^{n+1,i}={\bar {\boldsymbol {u}}}^{n}-\Delta t{\boldsymbol {M}}^{-1}[({\boldsymbol {A}}^{n+\theta _{1},i-1}+{\boldsymbol {K}}+{\hat {\boldsymbol {K}}}^{n+\theta _{1},i-1}){\bar {\boldsymbol {u}}}^{n+\theta _{1},i-1}-{\boldsymbol {G}}{\boldsymbol {p}}^{n+\theta _{2},i-1}+{\boldsymbol {C}}{\bar {\boldsymbol {c}}}^{n+\theta _{3},i-1}-{\boldsymbol {f}}^{n+1}]$
(44)

Step 2

 ${\boldsymbol {p}}^{n+1,i}=-{\boldsymbol {L}}^{-1}[{\boldsymbol {G}}^{T}{\bar {\boldsymbol {u}}}^{n+1,i}+{\boldsymbol {Q}}{\bar {\boldsymbol {\pi }}}^{n+\theta _{4},i-1}]$
(45)

Step 3

 ${\bar {\boldsymbol {c}}}^{n+1,i}=-{\boldsymbol {M}}^{-1}{\hat {\boldsymbol {C}}}^{n+1,i}{\bar {\boldsymbol {u}}}^{n+1,i}$
(46)

Step 4

 ${\bar {\boldsymbol {\pi }}}^{n+1,i}=-{\hat {\boldsymbol {M}}}^{-1}{\boldsymbol {Q}}^{T}{\bar {\boldsymbol {p}}}^{n+1,i}$
(47)

In above ${\textstyle \theta _{i}}$ are time integration parameters with ${\textstyle 0\leq \theta _{i}\leq 1}$ and ${\textstyle {\bar {(\cdot )}}^{n,i}}$ denotes nodal values at the ${\textstyle n}$th time step and the ith iteration. ${\textstyle {\boldsymbol {A}}^{n+\theta _{1},i-1}\equiv {\boldsymbol {A}}({\bar {\boldsymbol {u}}}^{n+\theta _{1},i-1})}$ etc. Also ${\textstyle (\cdot )^{n+\theta _{i},0}\equiv (\cdot )^{n}}$ for the computations in step 1 at the onset of the iterations.

Steps 1, 3 and 4 can be solved explicitely by choosing a lumped (diagonal) form of matrices ${\textstyle M}$ and ${\textstyle {\hat {\boldsymbol {M}}}}$. In this manner the main computational cost is the solution of step 2 involving the inverse of a Laplacian matrix. This can be solved very effectively using an iterative method.

For ${\textstyle \theta _{i}\not =0}$ the iterative proces is unavoidable. The iterations follow until convergence is reached. This can be measured using an adequate error norm in terms of the velocity and pressure variables, or the residuals. Indeed some ot the ${\textstyle \theta _{i}}$'s can be made equal to zero. Note that for ${\textstyle \theta _{2}=0}$ the algorithm is inconditionable unstable. A simple form is obtained making ${\textstyle \theta _{1}=\theta _{3}=\theta _{4}=0}$. This elliminates the non linear dependence with the velocity of matrices ${\textstyle {\boldsymbol {A}}}$ and ${\textstyle {\hat {\boldsymbol {K}}}}$ during the iterative scheme.

Convergence of above solution scheme is however difficult for some problems. An enhanced version of the algorithm can be obtained by adding the term ${\textstyle {\hat {\boldsymbol {L}}}({\bar {\boldsymbol {p}}}^{n+1,i}-{\bar {\boldsymbol {p}}}^{n+1,i-1})}$ where ${\textstyle {\hat {L}}_{ij}=\Delta t\int _{\Omega ^{e}}{\boldsymbol {\nabla }}^{T}N_{i}{\boldsymbol {\nabla }}N_{j}d\Omega }$ to the equation for the computation of the pressure in the second step. The new term acts as a preconditioner of the pressure equation given now by

 ${\bar {\boldsymbol {p}}}^{n+1,i}=-[{\boldsymbol {L}}+{\hat {\boldsymbol {L}}}]^{-1}[{\boldsymbol {G}}^{T}{\bar {\boldsymbol {u}}}^{n+1,i}+{\hat {\boldsymbol {L}}}{\bar {\boldsymbol {p}}}^{n+1,i-1}+{\boldsymbol {Q}}{\bar {\boldsymbol {\pi }}}^{n+\theta _{4},i-1}]$
(48)

Note that the added term vanishes for the converged solution (i.e. when ${\textstyle {\bar {\boldsymbol {p}}}^{n+1,i}={\bar {\boldsymbol {p}}}^{n+1,i-1}}$).

An alternative to above algorithm is to use the fractional step method described in the nex section.

### 6.2 Fractional step method

The pressure can be split from the discretized momentum equations (see Eq.(44)) as

 ${\bar {\boldsymbol {u}}}^{*}={\bar {\boldsymbol {u}}}^{n}-\Delta t{\boldsymbol {M}}^{-1}[({\boldsymbol {A}}^{n+\theta _{1}}+{\boldsymbol {K}}+{\hat {\boldsymbol {K}}}^{n+\theta _{1}}){\bar {\boldsymbol {u}}}^{n+\theta _{1}}-\alpha {\boldsymbol {G}}{\bar {\boldsymbol {p}}}^{n}+{\boldsymbol {C}}{\bar {\boldsymbol {c}}}^{n+\theta _{3}}-{\boldsymbol {f}}^{n+1}]$
(49)

 ${\bar {\boldsymbol {u}}}^{n+1}={\bar {\boldsymbol {u}}}^{*}+\Delta t{\boldsymbol {M}}^{-1}{\boldsymbol {G}}\delta {\bar {\boldsymbol {p}}}$
(50)

In above equations ${\textstyle \alpha }$ is a variable taking values equal to zero or one. For ${\textstyle \alpha =0}$, ${\textstyle \delta p\equiv p^{n+1}}$ (first order scheme) and for ${\textstyle \alpha =1}$, ${\textstyle \delta p=\Delta p}$ (second order scheme) (Codina, 2001 ). Note that in both cases the sum of Eqs.(49) and (50) gives the time discretization of the momentum equations with the pressures computed at ${\textstyle t^{n+1}}$. The value of ${\textstyle {\bar {\boldsymbol {u}}}^{n+1}}$ from Eq.(50) is substituted now into (41b) to give

 ${\boldsymbol {G}}^{T}{\bar {\boldsymbol {u}}}^{*}+\Delta t{\boldsymbol {G}}^{T}{\boldsymbol {M}}^{-1}{\boldsymbol {G}}\delta {\bar {\boldsymbol {p}}}+{\boldsymbol {L}}{\bar {\boldsymbol {p}}}^{n+1}+{\boldsymbol {Q}}{\bar {\boldsymbol {\pi }}}^{n+\theta _{4}}=0$
(51a)

The product ${\textstyle {\boldsymbol {G}}^{T}{\boldsymbol {M}}^{-1}{\boldsymbol {G}}}$ can be approximated by a laplacian matrix, i.e.

 ${\boldsymbol {G}}^{T}{\boldsymbol {M}}^{-1}{\boldsymbol {G}}={\hat {\boldsymbol {L}}}\quad {\hbox{with }}{\hat {\boldsymbol {L}}}\simeq \int _{\Omega ^{e}}{1 \over \rho }({\boldsymbol {\nabla }}^{T}N_{i}){\boldsymbol {\nabla }}N_{j}\,d\Omega$
(51b)

A semi-implicit algorithm can be derived as follows.

Step 1 Compute the nodal fractional velocities ${\textstyle {\boldsymbol {u}}^{*}}$ explicitely from Eq.(49) with ${\textstyle {\boldsymbol {M}}={\boldsymbol {M}}_{d}}$ where subscript ${\textstyle d}$ denotes a diagonal matrix.

Step 2 Compute ${\textstyle \delta {\bar {\boldsymbol {p}}}}$ from Eq.(51a) (using Eq.(51b)) as

 $\delta {\bar {\boldsymbol {p}}}=-({\boldsymbol {L}}+\Delta t{\hat {\boldsymbol {L}}})^{-1}[{\boldsymbol {G}}^{T}{\bar {\boldsymbol {u}}}^{*}+{\boldsymbol {Q}}{\bar {\boldsymbol {\pi }}}^{n+\theta _{4}}+{\boldsymbol {L}}{\bar {\boldsymbol {p}}}^{n}]$
(52)

Step 3 Compute the nodal velocities ${\textstyle {\bar {\boldsymbol {u}}}^{n+1}}$ explicitely from Eq.(50) with ${\textstyle {\boldsymbol {M}}={\boldsymbol {M}}_{d}}$

Step 4 Compute ${\textstyle {\bar {\boldsymbol {c}}}^{n+1}}$ explicitely from Eq.(46) using ${\textstyle {\boldsymbol {M}}_{d}}$.

Step 5 Compute ${\textstyle {\bar {\boldsymbol {\pi }}}^{n+1}}$ explicitely from Eq.(47) with ${\textstyle {\hat {\boldsymbol {M}}}={\bar {\boldsymbol {M}}}_{d}}$.

This algorithm has an additional step than the iterative algorithm of Section 6.1. The advantage is that now Steps 1 and 2 can be fully linearized by choosing ${\textstyle \theta _{1}=\theta _{3}=\theta _{4}=0}$. Also the equation for the pressure variables in Step 2 has improved stabilization properties due to the additional laplacian matrix ${\textstyle {\hat {\boldsymbol {L}}}}$.

Details on the stability properties of the FIC formulation for incompressible fluid flow problems can be found in Oñate et al. (2002) .

### 6.3 Discretization of free-surface wave equation

The solution in time of Eq.(39) can be written in terms of the nodal velocities computed from the flow solution, as

 $\beta ^{n+1}=\beta ^{n}-\Delta t\left[v_{i}^{n+1,i}{\partial \beta ^{n} \over \partial x_{i}}-v_{3}^{n+1,i}-{h_{\beta _{j}} \over 2}{\partial r_{\beta }^{n} \over \partial x_{j}}\right]$
(53)

Eq.(53) can now be discretized in space using the standard Galerkin method and solved explicitely for the nodal wave heights at ${\textstyle t^{n+1}}$. Typically the general algorithm will be as follows:

1. Solve for the nodal velocities ${\textstyle {\bar {\boldsymbol {u}}}^{n+1}}$ and the pressures ${\textstyle {\bar {\boldsymbol {p}}}^{n+1}}$ in the fluid domain using any of the algorithms of Sections 6.1 and 6.2. When solving for the pressure equation impose ${\textstyle p^{n+1}=p_{a}}$ at the free-surface ${\textstyle \Gamma _{\beta }}$.
2. Solve for the free-surface elevation ${\textstyle \beta ^{n+1}}$ (viz. Eq.(53)).
3. Compute the new position of the mesh nodes in the fluid domain at time ${\textstyle t^{n+1}}$. Alternatively, regenerate a new mesh.

The mesh updating process can also include the free-surface nodes, although this is not strictly necessary. An hydrostatic adjustement can be implemented once the new free-surface elevation is computed simply by imposing the pressure at the nodes on the reference surface as

 $p^{n+1}=p_{a}+\rho g\Delta \beta \quad {\hbox{with}}\quad \Delta \beta =\beta ^{n+1}-\beta ^{ref}$
(54)

where ${\textstyle g}$ is the gravity constant.

Eq.(54) takes into account the changes in the free-surface without the need of updating the reference surface nodes. A higher accuracy in the flow solution can be obtained by updating these nodes after a number of time steps.

## 7 FLUID-SHIP INTERACTION

The algorithms of previous section can be extended to account for the ship motion due to the hydrodynamic forces. Here the ship hull structure can be modelled as a rigid solid defined by the three translations and the three rotations of its center of gravity. Alternatively, the full deformation of the ship structure can be computed by modelling the actual stiffness of the hull, the decks and the different structural members. Indeed, the first option usually suffices for the hull shape optimization.

In both cases the computation of the ship motion involves solving the dynamic equations of the ship structure written as

 ${\boldsymbol {M}}_{s}{\ddot {\boldsymbol {d}}}+{\boldsymbol {K}}_{s}{\boldsymbol {d}}={\boldsymbol {f}}_{ext}$
(55)

where ${\textstyle {\boldsymbol {d}}}$ and ${\textstyle {\ddot {\boldsymbol {d}}}}$ are the displacement and acceleration vectors of the nodes discretizing the ship structure, respectively, ${\textstyle {\boldsymbol {M}}_{s}}$ and ${\textstyle {\boldsymbol {K}}_{s}}$ are the mass and stiffness matrices of the structure and ${\textstyle {\boldsymbol {f}}_{ext}}$ is the vector of external nodal forces accounting for the fluid flow loads induced by the pressure and the viscous stresses. Clearly the main driving forces for the motion of the ship is the fluid pressure which acts in the form of a surface traction. Indeed Eq.(55) can be augmented with an appropriate damping term. The form of all the relevant matrices and vectors can be found in standard books on FEM for structural analysis (Zienkiewicz and Taylor, Vol 2 2000 ).

Solution of Eq.(55) in time can be performed using implicit or fully explicit time integration algorithms. In both cases the values of the nodal displacements, the velocities and the accelerations at time ${\textstyle t^{n+1}}$ are found.

A simple coupled fluid-ship-structure solution in time using, for instance, the fractional step method of Section 6.2 (for ${\textstyle \theta _{1}=\theta _{2}=\theta _{3}=\theta _{4}=0}$) involves the following steps.

Step 1 Solve for the fractional velocities ${\textstyle {\bar {\boldsymbol {u}}}^{*}}$ using Eq.(49). Here use of ${\textstyle \alpha =1}$ is recommended.

Step 2 Compute ${\textstyle \delta {\bar {\boldsymbol {p}}}}$ from Eq.(51a) solving a simultaneous system of equations.

Step 3 Compute explicitely the nodal velocities ${\textstyle {\bar {\boldsymbol {u}}}^{n+1}}$ from Eq.(50) with a diagonal mass matrix.

Step 4 Compute explicitely the projected convective variables ${\textstyle {\bar {\boldsymbol {c}}}^{n+1}}$ from Eq.(46) using ${\textstyle {\boldsymbol {M}}_{d}}$.

Step 5 Compute explicitely the projected pressure gradients ${\textstyle {\bar {\boldsymbol {\pi }}}^{n+1}}$ from Eq.(47) using ${\textstyle {\hat {\boldsymbol {M}}}_{d}}$.

Step 6 Compute explicitely the new position of the free-surface elevation ${\textstyle {\bar {\boldsymbol {\beta }}}^{n+1}}$ from Eq.(53).

Step 7 Compute the movement of the ship by solving the dynamic equations of motion for the ship structure under the hydrodynamic forces induced by the pressures ${\textstyle {\bar {\boldsymbol {p}}}^{n+1}}$ and the viscous stresses ${\textstyle {\boldsymbol {s}}^{n+1}}$.

Step 8 Update the position of the mesh nodes in the fluid domain at ${\textstyle t^{n+1}}$ by using the mesh update algorithm described in the next section. The updating process can also include the free-surface nodes although this is not strictly necessary to be done at every time step and the hydrostatic adjustment of the pressure acting on the free-surface (Section 6.3) can be used as an alternative.

Obviously, the use of a value different from zero for the ${\textstyle \theta _{1},\theta _{2},\theta _{3}}$ or ${\textstyle \theta _{4}}$ parameters will require an iterative process between steps 1 and 8 until a converged solution for the variables describing the fluid and ship motions is found.

A common option is to update the position of the mesh nodes only when the iterative process for computing the fluid and ship variables has converged. Clearly the regeneration of the mesh is unavoidable when the distorsion of the elements exceed a certain limit.

## 8 A SIMPLE ALGORITHM FOR UPDATING THE MESH NODES

Different techniques have been proposed for dealing with mesh updating in fluid-structure interaction problems. The general aim of all methods is to prevent element distortion during mesh deformation (Tezduyar, 2001a ; Tezduyar et al., 1992bc [100-101]).

Tezduyar et al. (1992c)  and Chiandussi et al. (2000)  have proposed simple method for moving the mesh nodes based on the iterative solution of a fictious linear elastic problem on the mesh domain. In the method introduced in Tezduyar et al. (1992c) , the mesh deformation is handled selectively based on the element sizes and deformation modes, with the objective to increase stiffening of the smaller elements, which are typically located near solid surfaces. In Chiandusi et al. (2000)  in order to minimize the mesh deformation the “elastic” properties of each mesh element are appropiately selected so that elements suffering greater movements are stiffer. A simple and effective procedure is to select the Poisson's ratio ${\textstyle \nu =0}$ and compute the “equivalent” Young modulus for each element by

 $E={{\bar {E}} \over 3{\bar {\varepsilon }}^{2}}(\varepsilon _{1}^{2}+\varepsilon _{2}^{2}+\varepsilon _{3}^{2})$
(56)

where ${\textstyle \varepsilon _{i}}$ are the principal strains, ${\textstyle {\bar {E}}}$ is an arbitrary value of the Young modulus and ${\textstyle {\bar {\varepsilon }}}$ is a prescribed uniform strain field. ${\textstyle {\bar {E}}}$ and ${\textstyle {\bar {\varepsilon }}}$ are constant for all the elements in the mesh.

In summary, the solution scheme proposed by Chiandusi et al. (2000)  includes the following two steps.

Step 1. Consider the FE mesh as a linear elastic solid with homogeneous material properties characterized by a prescribed uniform strain field ${\textstyle {\bar {\varepsilon }}}$, an arbitrary Young modulus ${\textstyle {\bar {E}}}$ and ${\textstyle \nu =0}$. Solve a linear elastic problem with imposed displacements at the mesh boundary defined by the actual movement of the boundary nodes. An approximate solution to this linear elastic problem, such as that given by the first iterations of a conjugate gradient solution scheme, suffices for practical purposes.

Step 2. Compute the principal strains in each element. Repeat the (approximate) FE solution of the linear elastic problem with prescribed boundary displacements using the values of ${\textstyle E}$ of Eq.(56).

The previous algorithm is able to treat the movement of the mesh due to changes in position of fully submerged and semi-submerged bodies such as ships. However if the floating body intersects the free-surface, the changes in the analysis domain can be very important as emersion or inmersion of significant parts of the body can occur within a time step. A possible solution to this problem is to remesh the analysis domain. However, for most problems a mapping of the moving surfaces linked to mesh updating algorithm described above can avoid remeshing. The surface mapping technique used by the authors is based on transforming the 3D curved surfaces into reference planes (Figure 4). This makes it possible to compute within each plane the local (in-plane) coordinates of the nodes for the final surface mesh accordingly to the changes in the floating line. The final step is to transform back the local coordinates of the surface mesh in the reference plane to the final curved configuration which incorporates the new floating line (García, 1999 ; Oñate and García, 2001 ).

## 9 MODELLING OF THE TRANSOM STERN FLOW

The transom stern causes a discontinuity in the domain and the solution of the free-surface equation close to this region is inconsistent with the convective nature of the equation. This leads to instability of the wave height close to the transom region. This instability is found experimentally for low speeds. The flow at a sufficient high speed is physically more stable although it still can not be reproduced by standard numerical techniques (Reed et al., 1990) .

A solution to this problem is to apply adequate free-surface boundary conditions at the transom boundary. The obvious condition is to fix both the free surface elevation ${\textstyle \beta }$ and its derivative along the corresponding streamline to values given by the transom position and the surface gradient. However, prescribing those values can influence the transition between the transom flux and the lateral flux, resulting in unaccurate wave maps.

The method proposed in García and Oñate (2003)  is to extend the free-surface below the ship. In this way the necessary Dirichlet boundary conditions imposed at the inflow domain are enough to define a well posed problem. This method is valid both for the wetted and dry transom cases and it is also applicable to ships with regular stern.

This scheme does not work for partially wetted transoms. This situation can occur for highly unsteady flows where wake vortex induces the free-surface deformation and the flow remains adhered to the transom. To favour the computation of the free-surface, an artificial viscosity term is added to the free-surface equation in the vicinity of the transom in these cases.

## 10 LAGRANGIAN FLOW FORMULATION

The Lagrangian formulation is an effective (and relatively simple) procedure for modelling the flow of fluid particles undergoing severe distorsions such as water jets, high amplitude waves, braking waves, water splashing, filling of cavities, etc. Indeed the Lagrangian formulation is very suitable for treating ship hydrodynamic problems where the ship undergoes large motions. An obvious “a priori” advantage of the Lagrangian formulation is that both the ship and the fluid motion are defined in the same frame of reference.

The Lagrangian fluid flow equations are obtained by noting that the velocity of the mesh nodes and that of the fluid particles are the same. Hence the relative velocities ${\textstyle v_{i}}$ are zero in Eq.(21) and the convective terms vanish in the momentum equations, while the rest of the fluid flow equations remain unchanged. Also, the solution of the free surface equation is not needed in the Lagrangian formulation, as the free surface motion is modelled naturally by computing the displacement of (all) the fluid particles at each time step. In any case, the new position of the free surface nodes must be properly identified as described below.

The FEM algorithms for solving the Lagrangian flow equations are very similar to those for the ALE description presented earlier. We will focus in the second order fractional step algorithm of Section 6.2 (for ${\textstyle \theta _{1}=\theta _{4}=0}$ and ${\textstyle \alpha =1}$) accounting also for fluid-ship interaction effects.

Step 1 Compute explicitely a predicted value of the velocities ${\textstyle {\boldsymbol {u}}^{*}}$ as

 ${\bar {\boldsymbol {u}}}^{*}={\bar {\boldsymbol {u}}}^{n}-\Delta t{\boldsymbol {M}}_{d}^{-1}[{\boldsymbol {K}}{\bar {\boldsymbol {u}}}^{n}-{\boldsymbol {G}}{\bar {\boldsymbol {p}}}^{n}-{\boldsymbol {f}}^{n+1}]$
(57)

Note that matrices $A$ and ${\textstyle {\bar {\boldsymbol {K}}}}$ of Eq.(48) emanating from the convective terms have been elliminated.

Step 2 Compute ${\textstyle \delta {\bar {\boldsymbol {p}}}}$ from Eq.(52).

Step 3 Compute explicitely ${\textstyle {\bar {\boldsymbol {u}}}^{n+1}}$ from Eq.(50) with ${\textstyle {\boldsymbol {M}}={\boldsymbol {M}}_{d}}$.

Step 4 Compute ${\textstyle {\bar {\boldsymbol {\pi }}}^{n+1}}$ explicitely from Eq.(46).

Step 5 Solve for the motion of the ship structure by integrating Eq.(55).

Step 6 Update the mesh nodes in a Lagrangian manner as

 ${\boldsymbol {x}}_{i}^{n+1}={\boldsymbol {x}}_{i}^{n}+{\bar {\boldsymbol {u}}}_{i}\Delta t$
(58)

Step 7 Generate a new mesh and identify the new free surface nodes.

Further details of above algorithm can be found in Idelsohn et al. (2002 , 2003a ) and Oñate et al. (2003 , 2004 ).

The mesh regeneration process can be effectively performed using the extended Delaunay Tesselation described in Idelsohn et al. (2003b,c [55-56])). This method allows the fast generation of good quality meshes combining four node tetrahedra (or three node triangles in 2D) with hexahedra and non standard polyhedra such as pentahedra (or quadrilaterals and pentagons in 2D) where linear shape functions are derived using non-Sibsonian interpolation rules. The mesh regeneration can take place at each time step (this has been the case for the examples presented in the paper), after a prescribed number of time steps, or when the nodal displacements induce significant element distorsions.

The identification of the free-surface nodes in the Lagrangian analysis can be made using the Alpha Shape method. This is based on the search of all nodes which are on an empty Voronoi sphere with a radius greater than a specified distance defined in terms of the mesh size. For details see Edelsbrunner and Mucke (1994)  and Idelsohn et al. (2002 , 2003a,b,c [54-56]). A known value of the pressure (typically zero pressure or the atmospheric pressure) is prescribed at the free surface nodes.

The conditions of prescribed velocities or pressures at the solid boundaries in the Lagrangian formulation are usually applied on a layer of nodes adjacent to the boundary. These nodes typically remain fixed during the solution process. Contact between water particles and the solid boundaries is accounted for by the incompressibility condition which naturally prevents the water nodes to penetrate into the solid boundaries. This simple way to treat the water-wall contact is another attractive feature of the Lagrangian flow formulation. For details see Idelsohn et al. (2002 , 2003a ) and Oñate et al. (2003 , 2004 ).

## 11 MODELLING THE STRUCTURE AS A VISCOUS FLUID

A simple and yet effective way to analyze the rigid motion of solid bodies in fluids with the Lagrangian flow description is to model the solid as a fluid with a viscosity much higher than that of the surrounding fluid. The fractional step scheme of Section 10 can be readily applied skipping now step 5 to solve for the motion of both fluid domains (the actual fluid and the fictitious fluid modelling the rigid body displacements of the solid). An example of this type is presented in Sections 14.6.4 and 14.6.7.

Indeed this approach can be further extended to account for the elastic deformation of the solid treated now as a visco-elastic fluid. This will however introduce some complexity in the formulation and the full coupled fluid-structure interaction scheme previously described is preferable.

## 12 COMPUTATION OF THE CHARACTERISTIC LENGTHS

The evaluation of the stabilization parameters is a crucial issue in stabilized methods. Most existing methods use expressions which are direct extensions of the values obtained for the simplest 1D case. It is also usual to accept the so called SUPG assumption, i.e. to admit that vector ${\textstyle {\boldsymbol {h}}}$ has the direction of the velocity field. This restriction leads to instabilities when sharp layers transversal to the velocity direction are present. This deficiency is usually corrected by adding a shock capturing or crosswind stabilization term (Hughes and Mallet, 1986 ; Tezduyar and Partk, 1986 ; Codina, 1993 ). Indeed, in the FIC formulation the components of ${\textstyle {\boldsymbol {h}}}$ introduce the necessary stabilization along both the streamline and transversal directions to the flow.

Excellent results have been obtained in all problems solved with the ALE formulation using linear tetrahedra with the value of the characteristic length vector defined by

 ${\boldsymbol {h}}=h_{s}{{\boldsymbol {u}} \over {u}}+h_{c}{{\boldsymbol {\nabla }}u \over \vert {\boldsymbol {\nabla }}u\vert }$
(59)

where ${\textstyle u=\vert {\boldsymbol {u}}\vert }$ and ${\textstyle h_{s}}$ and ${\textstyle h_{c}}$ are the “streamline” and “cross wind” contributions given by

 ${\boldsymbol {h}}_{s}=\max({\boldsymbol {l}}_{j}^{T}{\boldsymbol {u}})/{u}$
(60)

 ${\boldsymbol {h}}_{c}=\max({\boldsymbol {l}}_{j}^{T}{\boldsymbol {\nabla }}u)/\vert {\boldsymbol {\nabla }}u\vert \quad ,\quad j=1,n_{s}$
(61)

where ${\textstyle {\boldsymbol {l}}_{j}}$ are the vectors defining the element sides (${\textstyle n_{s}=6}$ for tetrahedra).

The form chosen for the characteristic length parameters is similar to that proposed by other authors (see references of previous paragraph and also Donea and Huerta, 2003  and Tezduyar, 2001b ).

As for the free-surface equation the following value of the characteristic length vector ${\textstyle {\boldsymbol {h}}_{\beta }}$ has been taken

 ${\boldsymbol {h}}_{\beta }={\bar {h}}_{s}{{\boldsymbol {u}} \over {u}}+{\bar {h}}_{c}{{\boldsymbol {\nabla }}\beta \over \vert {\boldsymbol {\nabla }}\beta \vert }$
(62)

The streamline parameter ${\textstyle {\bar {\boldsymbol {h}}}_{s}}$ has been obtained by Eq.(60) using the value of the velocity vector ${\textstyle {\boldsymbol {u}}}$ over the 3 node triangles discretizing the free-surface and ${\textstyle n_{s}=3}$.

The cross wind parameter ${\textstyle {\bar {h}}_{c}}$ has been computed by

 ${\bar {h}}_{c}=\max[{\boldsymbol {l}}_{j}^{T}{\boldsymbol {\nabla }}\beta ]{1 \over \vert {\boldsymbol {\nabla }}\beta \vert }\quad ,\quad j=1,2,3$
(63)

The cross-wind terms in eqs.(59) and (62) take into account the gradient of the solution in the stabilization parameters. This is a standard assumption in shock-capturing stabilization procedures.

A more consistent evaluation of ${\textstyle {\boldsymbol {h}}}$ based on a diminishing residual technique can be found in Oñate and García (2001) .

## 13 TURBULENCE MODELLING

The detailed discussion on the treatment of turbulent effects in the flow equations falls outside the scope of this chapter. Indeed any of the existing turbulence models is applicable.

In the examples presented next we have chosen a turbulence model based on the Reynolds averaged Navier-Stokes equations where the deviatoric stresses are computed as sum of the standard viscous contributions and the so called Reynold stresses. Here we have chosen the Boussinesq assumption leading to a modification of the viscosity in the standard Navier-Stokes equations as sum of the “physical” viscosity ${\textstyle \mu }$ and a turbulent viscosity ${\textstyle \mu _{T}}$.

One of the simplest and more effective choices for ${\textstyle \mu _{T}}$ is the Smagorinski LES model giving

 $\mu _{T}=C_{l}h^{e}(2\varepsilon _{ij}\varepsilon _{ij})^{1/2}$
(64)

where ${\textstyle h^{e}}$ is the element size and ${\textstyle C}$ is a constant (${\textstyle C\simeq 0.01}$). In the examples analyzed ${\textstyle h^{e}}$ has been taken equal to ${\textstyle (\Omega ^{e})^{1/2}}$ and ${\textstyle (V^{e})^{1/3}}$ for 2D and 3D situations, respectively.

Many other options are possible such as the one and two equations turbulence models (i.e. the ${\textstyle k}$ model and the ${\textstyle k-\varepsilon }$ and ${\textstyle k-w}$ models) and the algebraic stress models. For further details the reader is refered to specialized publications (Celik el al., 1982 ; Wilcox, 1994 ).

## 14 EXAMPLES

The examples chosen show the applicability of the ALE and Lagrangian formulations presented to solve ship hydrodynamics problems. The fractional step algorithm of Section 6.2 with ${\textstyle \theta _{1}=\theta _{3}=\theta _{4}=0}$ and ${\textstyle \alpha =1}$ has been used. The first ALE example is the flow past a submerged NACA 0012 profile. The next ALE examples include the analysis of a Wigley hull, a scale model of a commercial ship and two American Cup racing sail boats. Numerical results obtained with linear tetrahedral meshes are compared with experimental data in all cases.

The last series of examples show applications of the Lagrangian formulation to simple schematic 2D ship hydrodynamics situations.

### 14.1 Submerged NACA 0012 profile

A submerged NACA0012 profile at ${\textstyle \alpha =5^{\circ }}$ angle of attack is studied using a 3D finite element model. This configuration was tested experimentally by Duncan (1983)  and modelled numerically using the Euler equations by several authors (Hino et al., 1993 ; Löhner et al., 1998 , 1999 ; Idelsohn et al., 1999 ). The submerged depth of the airfoil is equal to the chord length L. The Froude number for all the cases tested was set to ${\textstyle Fr={{U} \over {\sqrt {gL}}}=0.5672}$ where ${\textstyle U}$ is the incoming flow velocity at infinity.

The stationary free surface and the pressure distribution in the domain are shown in Figure 5. The non-dimensional wave heights compare well with the experimental results.

### 14.2 Wigley hull

We study here the well known Wigley Hull, given by the analytical formula ${\textstyle y=0.5B(1-4x^{2})(1-z^{2}/D^{2})}$ where ${\textstyle B}$ and ${\textstyle D}$ are the beam and the draft of the ship hull at still water.

The same configuration was tested experimentally (Noblesse and McCarthy, 1983 ; Wigley, 1983  ) and modelled numerically by several authors (Farmer et al., 1993 ; Storti et al., 1998 ; Idelsohn et al., 1999 ; Löhner et al., 1999 ). We use an unstructured 3D finite element mesh of ${\textstyle 65434}$ linear tetrahedra, with a reference surface of ${\textstyle 7800}$ triangles, partially represented in Figure 6. A Smagorinsky turbulence model was chosen.

Three different viscous cases were studied for ${\textstyle L=6m,F_{n}=0.316,\mu =10^{-3}Kg/m.s}$. In the first case the volume mesh was considered fixed, not allowing free-surface nor ship movements. Next, the volume mesh was updated due to free-surface movement, while keeping the model fixed. The third case corresponds to the analysis of a real free model including the mesh updating due to free-surface displacement and ship movement (sinkage and trim).

Table 1 shows the obtained total resistance coefficient in the three cases studied compared with experimental data.

 Experimental Numerical Test 1 ${\textstyle 5.2\quad 10^{-3}}$ ${\textstyle 4.9\quad 10^{-3}}$ Test 2 ${\textstyle 5.2\quad 10^{-3}}$ ${\textstyle 5.3\quad 10^{-3}}$ Test 3 ${\textstyle 4.9\quad 10^{-3}}$ ${\textstyle 5.1\quad 10^{-3}}$

Numerical values obtained for sinkage and trim were ${\textstyle -0.1\%}$ and ${\textstyle 0.035}$, respectively, while experiments gave ${\textstyle -0.15\%}$ and ${\textstyle 0.04}$.

Figure 6a shows the pressure distribution obtained near the Wigley hull for the free model. Some streamlines have also been plotted. The mesh deformation in this case is shown in Figure 6a

Comparison of the obtained body wave profile with experimental data for the free and fixed models is shown in Figure 6b.Significant differences are found close to stern for the fixed model. The free-surface contours for the truly free ship motion are shown in Figure 6c.

### 14.3 KVLCC2 hull model

The example is the analysis of the KVLCC2 benchmark model. Here a partially wetted tramsom stern is expected due to the low Froude number of the test. Figure 7 shows the NURBS geometry obtained from the Hydrodynamic Performance Research team of Korea (KRISO). Numerical results are compared with the experimental data available from KRISO (2000).

The smallest element size used was ${\textstyle 0.001}$ m and the largest ${\textstyle 0.50}$ m. The surface mesh chosen is also shown in Figure 7. The 3D mesh included ${\textstyle 550.000}$ tetrahedra. The tramsom stern flow model described was used.

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

• Length: ${\textstyle 5.52m}$, Beam (at water plane): ${\textstyle 0.82m}$, Draught: ${\textstyle 0.18m}$, Wetted Surface: ${\textstyle 8.08m^{2}}$.
• Velocity: ${\textstyle 1.05m/seg}$, Froude Number: ${\textstyle 0.142}$.
• Viscosity: ${\textstyle 0.00126Kg/m\cdot seg}$, Density: ${\textstyle 1000Kg/m^{3}}$, Reynolds number: ${\textstyle 4.63\cdot 10^{6}}$.

The turbulence model chosen was the ${\textstyle k}$ model. Figures 8 and 9 show the wave profiles on the hull and in a cut at ${\textstyle y/L=0.0964}$ obtained for Test 1, compared to experimental data. The results are quantitatively good.

Test 2.- Wake analysis at different planes. Several turbulence models were used (Smagorinsky, ${\textstyle k}$ and ${\textstyle k-\epsilon }$ model) in order to verify the quality of the results. Here, only the results from the ${\textstyle k-\epsilon }$ model are shown. The velocity maps obtained even for the simplest ${\textstyle Smagorinsky}$ model were qualitatively good. The main characteristics of this analysis are:

• Length: ${\textstyle 2.76m}$, Beam (at water plane): ${\textstyle 0.41m}$, Draught: ${\textstyle 0.09m}$, Wetted Surface: ${\textstyle 2.02m^{2}}$.
• Velocity: ${\textstyle 25m/seg}$.
• Viscosity: ${\textstyle 3.05\cdot 10^{-5}Kg/m\cdot seg}$, Density: ${\textstyle 1.01Kg/m^{3}}$, Reynolds number: ${\textstyle 4.63\cdot 10^{6}}$.

Figures 10 and 11 present some results for Test 2. Figure 10 shows the contours of the axial (X) component of the velocity on a plane at ${\textstyle 2.71m}$ from the orthogonal aft. Figure 11 shows the map of the kinetic energy at this plane. Experimental results are also plotted for comparison. Further results are reported in García and Oñate (2003).

### 14.4 American Cup Rioja de España Model

The Rioja de España boat representing Spain in the American Cup's edition of 1995 was analyzed. Figure 12 shows the geometry of the boat based on standard NURBS surfaces. Numerical results are compared with an extrapolation of experimental data obtained in CEHIPAR basin (Spain) using a 1/3.5 scale model. Resistance tests were performed with the model free to sink and trim. Experimental data include drag, lift, sinkage, trim angles and wave profiles at ${\textstyle 4.27m}$ (real scale) from symmetry plane. Every model was towed at equivalent velocities of ${\textstyle 10}$, ${\textstyle 9}$, ${\textstyle 8.5}$, ${\textstyle 8.0}$, ${\textstyle 7.5}$ and ${\textstyle 7.0}$ knots.

Numerical analysis were carried out at real scale. Characteristics of unstructured grids of four node linear tetrahedra used are shown in Table 2 together with the parameters of some of the model studied.

 Test Geometry Heel Drift Symmetry Elements nodes E0D0 Hull, bulb and keel $0^{\circ }$ $0^{\circ }$ Yes 700 000 175 000 E15D2 Hull, bulb, keel and rudder $15^{\circ }$ $2^{\circ }$ No 1 500 000 380 000 E15D4 Hull, bulb, keel and rudder $15^{\circ }$ $4^{\circ }$ No 1 500 000 380 000 E25D2 Hull, bulb, keel and rudder $25^{\circ }$ $2^{\circ }$ No 1 500 000 380 000

All grids used were generated with the same quality criteria and using element sizes from ${\textstyle 5mm}$ to ${\textstyle 2000mm}$. Some details of the grid used in the E0D0 case are shown in Figures 13 and 14. A ${\textstyle k-\varepsilon }$ turbulence model in combination with an extended law of the wall was chosen for the analysis.

A time step of ${\textstyle 0.1s}$ was used and this sufficed to achieve a steady state solution in all cases. Additional calculations were carried out with ${\textstyle \Delta t=0.025s}$ and ${\textstyle 0.01s}$, in order to verify the influence of the time increment in the accuracy of the results. No significant influence was detected for the selected results.

Figure 16 shows a comparison of simulated and experimental wave profiles. The origin of the x axis has been taken at the fore perpendicular and the x+ sense is afore. In the non symmetric case, the measurements were performed at the opposite side of the heeled board. The ratio of maximum amplitudes of fluctuations (noise) and waves in the experimental measurements varies from ${\textstyle 4\%}$ to ${\textstyle 12\%}$

Figures 16 and 17 show pressure contours on the bulb and keel for different cases, corresponding to a velocity of ${\textstyle 8kn}$. Figures 18 to 19 show some of the wave patterns obtained for a velocity of ${\textstyle 9kn}$. Figures 21 and 22 show some perspective views, including pressure and velocity contours, streamlines and cuts for some cases analyzed. Finally Figures 23 and 24 show resistance graphs where numerical results are compared with values extrapolated from experimental data. For further details see García et al. (2002) .

### 14.5 American Cup BRAVO ESPAÑA Model

The finite element formulation presented was also applied to study the racing sail boat Bravo España participating in the 1999 edition of the American Cup. The mesh of linear tetrahedra used is shown in Figure 25. Results presented in Figures 2527 correspond to a non symmetrical case including appendages. Good comparison between experimental data and numerical results was again obtained Further details can be found in García and Oñate (2002) .

### 14.6 Lagrangian flow examples

A number of simple problems have been solved in order to test the capabilities of the Lagrangian flow formulation to solve ship hydrodynamics problem. All examples have been analyzed with the algorithm described in Section 10. We note again that the regeneration of the mesh has been performed at each time step.

#### 14.6.1 Rigid ship hit by wave

The first example is a very schematic representation of a ship when is hit by a big wave produced by the collapse of a water recipient (Fig. 28). The ship can not move and initially the free-surface near the ship is horizontal. Fixed nodes represent the ship as well as the domain walls. The example tests the suitability of the Lagrangian flow formulation to solve water-wall contact situations even in the presence of curved walls. Note the breaking and splashing of the waves under the ship prow and the rebound of the incoming wave. It is also interesting to see the different water-wall contact situations at the internal and external ship surfaces and the moving free-surface at the back of the ship.

#### 14.6.2 Horizontal motion of a rigid ship in a reservoir

In the next example (Fig. 29) the same ship of the previous example moves now horizontally at a fixed velocity in a water reservoir. The free-surface, which was initially horizontal, takes the correct position at the ship prow and stern. Again, the complex water-wall contact problem is naturally solved in the curved prow region.

#### 14.6.3 Semi-submerged rotating water mill

The example shown in Figure 30 is the analysis of a rotating water mill semi submerged in water. The blades of the mill are treated as a rigid body with an imposed rotating velocity, while the water is initially in a stationary flat position. Fluid structure interactions with free-surfaces and water fragmentation are well reproduced in this example.

#### 14.6.4 Floating wood piece

The next example shows an initially stationary recipient with a floating piece of wood where a wave is produced on the left side. The wood has been simulated by a liquid of higher viscosity as described in Section 11. The wave intercepts the wood piece producing a breaking wave and displacing the floating wood as shown in Figure 31.

#### 14.6.5 Ships hit by an incoming wave

In the example of Figure 32 the motion of a fictitious rigid ship hit by an incoming wave is analyzed. The dynamic motion of the ship is induced by the resultant of the pressure and the viscous forces acting on the ship boundaries. We note that the horizontal displacement of the gravity center of the ship was fixed to zero. In this way, the ship moves only vertically although it can freely rotate. The position of the ship boundary at each time step defines a moving boundary condition for the free surface particles in the fluid.

Figure 32 shows instants of the motion of the ship and the surrounding fluid. It is interesting to see the breaking of a wave on the ship prow as well as on the stern when the wave goes back. Note that some water drops slip over the ship deck.

Figure 33 shows the analysis of a similar problem using precisely the same approach. The section of the ship analyzed corresponds now to that of a real container ship. Differently to the previous case the rigid ship is free to move laterally due to the sea wave forces. The objective of the study was to asses the influence of the stabilizers in the ship roll. The figures show clearly how the lagrangian method predicts the ship and wave motions in a realistic manner.

#### 14.6.6 Tank-ship hit by a lateral wave

Figure 34 represents a transversal section of a tank-ship and a wave approaching to it. The tank-ship is modelled as a rigid body and it carries a liquid inside which moves freely within the ship domain.

Initially (${\textstyle t=0.00}$) the ship is released from a fixed position and the equilibrium configuration found is consistent with Arquimides principle. During the following time steps the external wave hits the ship and both the internal and the external fluids interact with the ship boundaries. At times ${\textstyle t=6.15}$ and ${\textstyle 8.55}$ breaking waves and some water drops slipping along the ship deck can be observed. Figure 35 shows the pressure pattern at two time steps.

#### 14.6.7 Rigid cube falling in a recipient with water

The last examples show the capacity of the Lagrangian formulation to analyze the motion of solids within water. In the first example a solid cube is initially free and falls down within a water recipient. In this example, the rigid solid is modeled first as a fictitious fluid with a higher viscosity, similarly as for the floating solid of Section 14.6.4. The results of this analysis are shown in Figure 36. Note that the method reproduces very well the interaction of the cube with the free surface as well as the overall sinking process. A small deformation of the cube is produced. This can be reduced by increasing further the fictitious viscosity of the cube particles.

The same problem is analyzed again considering now the cube as a rigid solid subjected to pressure and viscous forces acting in its boundaries. The resultant of the fluid forces and the weight of the cube are applied to the center of the cube. These forces govern the displacement of the cube which is computed by solving the dynamic equations of motion as described in the fractional step algorithm of Section 10, similarly as for the rigid ships of the three previous examples. Similarly as in those cases the moving cube contours define a boundary condition for the fluid particles at each time step.

Initially the solid falls down freely due to the gravity forces (Figure 37). Once in contact with the water surface, the alpha-shape method recognizes the different boundary contours which are shown with a thick line in the figure. The pressure and viscous forces are evaluated in all the domain and in particular on the cube contours. The fluid forces introduce a negative acceleration in the vertical motion until, once the cube is completely inside the water, the vertical velocity becomes zero. Then, Arquimides forces bring the cube up to the free-surface. It is interesting to observe that there is a rotation of the cube. The reason is that the center of the floating forces is higher in the rotated position than in the initial ones.

Figure 38 shows a repetition of the same problem showing now all the finite elements in the mesh discretizing the fluid. We recall that in all the Lagrangian problems here described the mesh in the fluid domain is regenerated at each time step combining linear triangles and quadrilateral elements as described in Section 10. Note that some fluid particles separate from the fluid domain. These particles are treated as free boundary points with zero pressure and hence fall down due to gravity.

It is interesting to see that the final position of the cube is different from that of Figure 37. This is due to the unstable character of the cube motion. A small difference in the numerical computations (for instance in the mesh generation process) shifts the movement of the cube towards the right or the left. Note however that a final rotated equilibrium position is found in both cases.

Further examples of the Lagrangian flow formulation can be found in Oñate et al. (2003 , 2004 ) and Idelsohn et al. (2002 , 2003a ).

## 15 CONCLUDING REMARKS

An overview of different finite element schemes for solving ship hydrodynamics problems has been presented. The necessary stabilization in the numerical computation has been provided by the finite calculus formulation of the fluid flow equations accounting for surface wave effects.

Stabilized finite algorithms for solving a variety of ship hydrodynamics situations using ALE and fully Lagrangian descriptions have been presented. Both monolithic and fractional step algorithms have been derived. The interaction of the motion of the ship structure with the fluid equations can be accounted for in a straight forward manner within the general flow solution schemes. The ALE formulation is particularly adequate for ship hydrodynamics problems involving free-surface waves of moderate amplitude. The Lagrangian flow formulation allows to solve in a simple and effective manner ship hydrodynamics problems involving large motions of the free-surface and complex fluid-structure interaction situations.

## ACKNOWLEDGEMENTS

The authors thank Prof. R. Löhner, Prof. P. Bergan, Dr. C. Sacco and Dr. H. Sierra for many useful discussions.

Thanks are also given to Dr. F. del Pin and Mr. N. Calvo for their help and involvement in solving the problems with the Lagrangian formulation.

The authors are grateful to Copa America Desafio Español SA and CEHIPAR for providing the geometry and experimental data for the racing boats analyzed and to Det Norske Veritas for providing the geometry of the ship with stabilizers for the analysis shown in Figure 33.

The Spanish company COMPASS Ingeniería y Sistemas SA provided the computer code Tdyn for the numerical simulation of the examples solved with the ALE formulation (Tdyn, 2004 ). This support is gratefully acknowledged.

Back to Top

### Document information Published on 29/05/19

DOI: 10.1002/9781119176817
Licence: CC BY-NC-SA license

### Document Score 0

Views 1
Recommendations 0

### Keywords 