Julioges

Abstract

A stabilized semi-implicit fractional step finite element method for solving coupled fluid-structure interaction problems involving free surface waves is presented. The stabilized governing equations for the viscous incompressible fluid and the free surface are derived at a differential level via a finite calculus procedure. A mesh updating technique based on solving a fictitious elastic problem on the moving mesh is described. Examples of the efficiency of the stabilized semi-implicit algorithm for the analysis of fluid-structure interaction problems in totally or partially submerged bodies is presented.

Introduction

Accurate prediction of the fluid-structure interaction effects for a totally or partially submerged body in a flowing liquid including a free surface is a problem of great relevance in civil and offshore engineering and naval architecture among many other fields.

The difficulties in accurately solving the coupled fluid-structure interaction problem in this case are mainly due to the following reasons:

1. The difficulty of solving numerically the incompressible fluid dynamic equations which typically include intrinsic nonlinearities except for the simplest and limited potential flow model.
2. The obstacles in solving the constraint equation stating that the fluid particles remain on the free surface boundary which position is in turn unknown.
3. The difficulties in predicting the motion of the submerged body due to the interaction forces while minimizing the distortion of the finite elements discretizing the fluid domain, thus reducing the need for remeshing.

This paper extends recent work of the authors [1–5] to derive a stabilized finite element method which overcomes the above three obstacles. The starting point are the modified governing differential equations for the incompressible viscous flow and the free surface condition incorporating the necessary stabilization terms via a finite calculus (FIC) procedure developed by the authors [6–11]. The FIC approach has been successfully applied to the finite element and meshless solution of a range of advective-diffusive transport and fluid flow problems [1–13].

The stabilized governing equations are written in an arbitrary Lagrangian-Eulerian (ALE) form to account for the effect of relative movement between the mesh and the fluid points. These equations are solved in space-time using a semi-implicit fractional step approach and the finite element method (FEM). Free surface wave boundary effects are accounted for in the flow solution either by moving the free surface nodes in a Lagrangian manner, or else via the introduction of a prescribed pressure at the free surface computed from the wave height.

The movement of a fully or partially submerged body within the fluid due to the interaction forces is treated by solving a structural dynamic problem using the fluid forces as input loads. A method to update the mesh in the fluid domain following the movement of the submerged body with minimum element distortion is presented. The mesh update procedure is based on the finite element solution of a linear elastic problem on the mesh domain, where fictitious elastic properties are assigned so that elements suffering a larger straining are stiffer [14].

The content of the paper is structured as follows. First, details of the stabilized form of the governing equations for a viscous flow and the free surface using a finite calculus procedure are given. The semi-implicit fractional step approach using the FEM is then described. Details of the computation of the stabilization parameters are also given. Next the mesh updating procedure is presented. Finally, the efficiency of the method proposed is shown in the 3D analysis of the standard square cavity problem and of several fluid-structure interaction problems with free surface waves.

Finite calculus formulation of fluid-flow and free surface equations

The finite element solution of the incompressible Navier-Stokes equations with the classical Galerkin method may suffer from numerical instabilities from two main sources. The first is due to the advective-diffusive character of the equations which induces oscillations for high values of the velocity. The second source has to do with the mixed character of the equations which limits the choice of finite element interpolations for the velocity and pressure fields.

Solutions for these two problems have been extensively sought in the last years. Compatible velocity-pressure interpolations satisfying the inf-sup condition emanating from the second problem above mentioned have been used [15–17]. In addition, the advective operator has been modified to include some “upwinding” effects [18–25]. Recent procedures based on Galerkin Least Square (GLS) [26,27], Characteristic Galerkin [28,29], Variational Multiscale [30–32] and Residual Free Bubbles [33–35] techniques allow equal order interpolation for velocities and pressure by introducing a Laplacian of pressure term in the mass balance equation, while preserving the upwinding stabilization of the momentum equations. Most of these methods lack enough stability in the presence of sharp layers transversal to the velocity. This deficiency is usually corrected by adding new “shock capturing” stabilization terms to the already stabilized equations [36–37]. The computation of the stabilization parameters in all these methods is typically based on “ad hoc” generalizations of the parameters for the 1D linear advective-diffusive-reactive problem [38,39].

Applications of stabilized GLS FEM to fluid-structure interaction problems, mainly of aerolastic type, have been reported in [41–50].

Introduction of the free surface boundary condition in the flow equations increases considerably the difficulty of solving fluid-structure interation problems using FEM. A review of these difficulties and some solution procedures can be found in [51]. Another successful application of stabilized FEM to free surface wave problems was reported in [52].

This paper presents a different approach for deriving stabilized finite element methods for incompressible flow problems with a free surface. The starting point is the stabilized form of the governing differential equations derived via a finite calculus (FIC) procedure. This technique first presented in [6,7] is based on writting the different balance equations over a domain of finite size and retaining higher order terms. These terms incorporate the ingredients for the necessary stabilization of any transient and steady state numerical solution already at the differential equations level. Application of the standard Galerkin formulation to the consistently modified differential equations for the fluid flow problem leads to a stabilized system of discretized equations which overcomes the two problems above mentioned (i.e. the advective type instability and that due to lack of compatibility between the velocity and pressure fields). Application of the FIC method to the free surface wave problem leads to a new stabilized governing equation for the free surface which again can be solved numerically by standard Galerkin FEM. In addition, the modified differential equations can be used to derive a numerical scheme for iteratively computing the stabilization parameters [7–9].

Ilinca et al. [40] have recently proposed a stabilized FEM for incompressible advective-diffusive transport and fluid flow problems based on applying the standard Galerkin technique to the modified governing differential equations obtained by expanding the residuals around a known finite element solution using Taylor series. The set of modified equations ressembles those obtained by the FIC method using a conceptually different procedure.

Initial applications of the FIC method to solve free surface ship wave problems were reported in [1–5]. Idelsohn et al. [51] have shown that starting from the stabilized FIC form of the free surface equation allows the identification of a number of stabilized upwinding finite difference schemes traditionally used for solving free surface problems in naval architecture.

The FIC formulation presented in this paper for incompressible flows with a free surface can be considered an extension of that recently developed in [10] for finite element analysis of incompressible Navier-Stokes flows. A new formulation of the stabilized governing differential equations via the FIC method is here presented which holds for the viscous (Stokes) and zero viscosity (Euler) limit cases. The stabilized fluid flow equations are completed with the FIC form of the free surface wave equation following the ideas first presented in [2]. The set of stabilized governing equations is first discretized in time and then solved in space using a Galerkin finite element method.

A semi-implicit fractional step procedure is used for the momentum and mass balance equations allowing for equal order linear interpolations of the velocity and pressure variables over tetrahedral elements. Examples of application of the new stabilized finite element formulation to the standard square cavity flow problem and to a number of free surface ship-wave problems, including coupled fluid-structure interaction situations, are presented.

For the sake of preciseness, the basic ideas of the FIC method are given next.

Basic concepts of the finite calculus (FIC) method

Let us consider a sourceless transient problem over a one dimensional domain $AB$ of length $L$ (Figure 1). The balance of flux $q$ over a domain of finite size belonging to $L$ can be written as

$q_{A}-q_{B}=0$

where $A$ and $B$ are the end points of the finite size domain of length $h$. As usual $q_{A}$ and $q_{B}$ represent the values of the flux $q$ at points $A$ and $B$, respectively.

For instance, in an 1D advective-diffusive problem the flux $q=-cu\phi +k{d\phi \over dx}$, where $\phi$ is the transported variable (i.e. the temperature in a thermal problem), $u$ is the advective velocity and $c$ and $k$ are the advective and diffusive material parameters, respectively.

The flux $q_{A}$ can be expressed in terms of the values at point $B$ by the following Taylor series expansion

$q_{A}=q_{B}-h{\partial q \over \partial x}\vert _{B}+{h^{2} \over 2}{d^{2}q \over dx^{2}}\vert _{B}+Oh^{3}$

Substituting (2) into (1) gives after simplification and neglecting cubic terms in $h$ ${dq \over dx}-{\underline {{h \over 2}{d^{2}q \over dx^{2}}}}=0$

where all terms are evaluated at the arbitrary point $B$.

Eq. (3) is the finite form of the balance equation over the domain $AB$. The underlined term in eq.(3) introduces the necessary stabilization for the discrete solution of eq.(3) using any numerical technique. Distance $h$ is the characteristic length of the discrete problem and its value depends on the parameters of the discretization method chosen (such as the grid size [6–10]). Note that for $h\to 0$ the standard infinitesimal form of the balance equation $\left({dq \over dx}=0\right)$ is recovered.

The above process can be extended to derive the stabilized balance differencial equations for any problem in fluid or solid mechanics as

$r_{d}-{\underline {{h_{j} \over 2}{\partial r_{i} \over \partial x_{j}}}}=0$

where $r_{i}$ is the standard form of the $i$th differential equation for the infinitesimal problem, $h_{j}$ are the dimensions of the domain where balance of fluxes, forces, etc. is enforced, and $j=1,2,3$ for 3D problems. It is important to note that the numerical solution of eq.(4) (together with the appropriate stabilized boundary conditions) using Galerkin FE or central finite difference schemes leads to stable results [6–11]. Details of the derivation of eq.(4) for steady-state and transient advective-diffusive and fluid flow problems can be found in [6]. Applications of the FIC approach to the solution of these problems using Galerkin finite element and meshless procedures are reported in [1–13].

The underlined stabilization terms in eqs.(3) and (4) are a consequence of accepting that the infinitesimal form of the balance equations is an unreachable limit within the framework of a discrete numerical solution. Indeed eqs.(3) or (4) are not longer valid for obtaining an analytical solution following traditional integration methods from infinitesimal calculus theory. The meaning of the new stabilized equations makes sense only in the context of a discrete numerical method yielding approximate values of the solution at a finite set of points within the analysis domain. Convergence to the exact analytical value at the points will occur only for the limit case of zero grid size (except for some simple 1D problems [6,11]) which also implies naturally a zero value of the characteristic length parameters.

FIC formulation of viscous flow and free surface equations

We consider the motion around a body of a viscous incompressible fluid including a free surface.

The stabilized FIC form of the governing differential equations for the three dimensional (3D) problem can be written in Arbitrary Lagrangian-Eulerian (ALE) form as [2,3,10]

.3cm Momentum $r_{m_{i}}-{\underline {{1 \over 2}h_{mj}{\partial r_{m_{i}} \over \partial x_{j}}}}-{\underline {{1 \over 2}\delta {\partial r_{m_{i}} \over \partial t}}}=0\qquad {\hbox{on }}\Omega \quad i,j=1,2,3$

Mass balance -.3cm $r_{d}+{\underline {{1 \over 2}h_{dj}{\partial r_{d} \over \partial x_{j}}}}=0\qquad {\hbox{on }}\Omega \quad j=1,2,3$

Free surface -.3cm

$r_{\beta }-{\underline {{1 \over 2}h_{\beta _{j}}{\partial r_{\beta } \over \partial x_{j}}}}-{\underline {{1 \over 2}\gamma {\partial r_{\beta } \over \partial t}}}=0\qquad {\hbox{on }}\Gamma _{\beta }\quad j=1,2$

where Error al representar (función desconocida «\noalign»): {\displaystyle \begin{aligned} r_{m_i}& =&\rho \left[{\partial u_i\over \partial t} + v_j {\partial u_i\over \partial x_j} \right] + {\partial p\over \partial x_i} - {\partial \tau_{ij}\over \partial x_j} \\ r_d&=&{\partial u_i\over \partial x_i} \qquad i=1,2,3 \\ r_\beta &=&{\partial \beta \over \partial t} + v_i {\partial \beta \over \partial x_i}-v_3 \qquad i=1,2 \\ \noalign{\hbox{and}} v_i&=&u_i-u_i^m\end{aligned}}

Above, $u_{i}$ is the velocity along the $i$th global reference axis, $u_{i}^{m}$ is the velocity of the mesh nodes and $v_{i}$ is the relative velocity between the moving mesh and the fluid point $i$, $\rho$ is the (constant) density of the fluid, $p$ is the dynamic pressure defined as $p={1 \over \rho }p_{a}-gx_{3}$ where $p_{a}$ is the absolute pressure and $x_{3}$ is the vertical coordinate, $\beta$ is the wave elevation (measured with respect to a reference flat surface) and $\tau _{ij}$ are the viscous stresses related to the viscosity $\mu$ by the standard expression $\tau _{ij}=\mu \left({\partial u_{i} \over \partial x_{j}}+{\partial u_{j} \over \partial x_{i}}-\delta _{ij}{1 \over 3}{\partial u_{k} \over \partial x_{k}}\right)$ where $\delta _{ij}$ is the Kronecker delta.

The boundary conditions for the stabilized problem are written as

$n_{j}\tau _{ij}+t_{i}+{\underline {{1 \over 2}h_{mj}n_{j}r_{m_{i}}}}=0\qquad {\hbox{on }}\Gamma _{t}$

$u_{j}-u_{j}^{p}=0\qquad {\hbox{on }}\Gamma _{u}$

where $n_{j}$ are the components of the unit normal vector to the boundary and $t_{i}$ and $u_{j}^{p}$ are prescribed tractions and displacements on the boundaries $\Gamma _{t}$ and $\Gamma _{u}$, respectively.

The underlined terms in eqs.(5)–(7) introduce the necessary stabilization for the approximated numerical solution.

The characteristic length distances $h_{mj}$ and $h_{dj}$ represent the dimensions of the finite domain where balance of momentum and mass is enforced. On the other hand, the characteristic distances $h_{\beta j}$ in eq.(7) represent the dimensions of a finite domain surrounding a point where the velocity is constrained to be tangent to the free surface. The signs before the stabilization terms in eqs.(5)–(7) and (13) ensure a positive value of the characteristic length distances. The parameters $\delta$ and $\gamma$ in eqs.(5) and (7) have dimensions of time. Details of the derivation of eqs. (5)–(7) can be found in [2,6,10]. As an example, the stabilized equation for the free surface (eq.(7)) is derived in the Appendix.

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

Remark 1

In reference [10] a modified version of the Dirichlet condition (14) is used including an additional stabilization term. This term is not strictly necessary for the subsequent derivation and will be neglected here.

Alternative form of the mass balance equation

Taking the first derivative of eq.(12) gives (assuming the viscosity $\mu$ to be constant) ${\partial \tau _{ij} \over \partial x_{j}}=\mu \Delta u_{i}+{\mu \over 3}{\partial r_{d} \over \partial x_{i}}$ where $\Delta ={\partial ^{2} \over \partial x_{i}\partial x_{i}}$ is the Laplacian operator. Substituting eq.(15) into (5) gives after algebraic rearrangement, ${\partial r_{d} \over \partial x_{i}}=\left({\mu \over 3}+{\rho u_{i}h_{m_{i}} \over 2}\right)^{-1}{\bigg [}{\bar {r}}_{m_{i}}-{h_{m_{k}} \over 2}{\partial r_{m_{i}} \over \partial x_{k}}+{\rho u_{i}h_{m_{i}} \over 2}{\partial r_{d} \over \partial x_{i}}-{\delta \over 2}{\partial r_{m_{i}} \over \partial t}{\bigg ]}\quad {\hbox{no sum in}}\,\,i$ where ${\bar {r}}_{m_{i}}=r_{m_{i}}+{\mu \over 3}{\partial r_{d} \over \partial x_{i}}$ and $r_{m_{i}}$ is given by eq.(8).

Inserting eq.(16) into eq.(6) gives $r_{d}+c_{i}\left({\bar {r}}_{m_{i}}-{h_{m_{k}} \over 2}{\partial r_{m_{i}} \over \partial x_{k}}+{\rho u_{i}h_{m_{i}} \over 2}{\partial r_{d} \over \partial x_{i}}-{\delta \over 2}{\partial r_{m_{i}} \over \partial t}\right)=0\quad {\hbox{no sum in}}\,\,i$ with $c_{i}=\left({2\mu \over 3h_{d_{i}}}+{\rho u_{i}h_{m_{i}} \over h_{d_{i}}}\right)^{-1}\qquad {\hbox{no sum in}}\,\,i$

Extracting the pressure terms from the brackets in (18) gives $r_{d}-g_{ii}{\partial ^{2}p \over \partial x_{i}\partial x_{i}}+r_{p}=0$ with

$r_{p}=c_{i}{\bar {r}}_{m_{i}}-g_{ij}{\partial \over \partial x_{j}}\left(r_{m_{i}}-\delta _{ij}{\partial p \over \partial x_{i}}\right)+{\rho u_{i}h_{m_{i}} \over 2}{\partial r_{d} \over \partial x_{i}}-{\delta \over 2}{\partial r_{m_{i}} \over \partial t}\quad {\hbox{no sum in}}\,\,i$

where -.6cm $g_{ij}=\left({4\mu \over 3h_{d_{i}}h_{m_{j}}}+{2\rho u_{i}h_{m_{i}} \over h_{d_{i}}h_{m_{j}}}\right)^{-1}\quad {\hbox{no sum in}}\,\,i$

Note that for $h_{m_{i}}=h_{d_{i}}=h$ where $h$ is a typical grid dimension (i.e. the average element size), the value of $g_{ii}$ is simply

$g_{ii}=\left({4\mu \over 3h^{2}}+{2\rho u_{i} \over h}\right)^{-1}$

The stabilization parameter $g_{ii}$ has now the form traditionally used in the GLS formulation for the viscous (Stokes) limit ($u_{i}=0$) and the inviscid (Euler) limit ($\mu =0$) and deduced from ad-hoc extensions of the 1D advective-diffusive problems [18–28]. Note, however, that the general form of the stabilization parameter $g_{ij}$ is deduced here from the general FIC formulation without further extrinsic assumptions.

Indeed, the precise computation of the characteristic length values is crucial for the practical application of above stabilized expressions. This problem is dealt with in a later section.

Fractional step approach

The momentum equations (5) are first discretized in time using the following scheme

$u_{i}^{n+1}=u_{i}^{n}-{\Delta t \over \rho }{\bigg [}\rho v_{j}^{n}{\partial u_{i}^{n} \over \partial x_{j}}+{\partial p^{n} \over \partial x_{i}}-{\partial \tau _{ij}^{n} \over \partial x_{j}}-{h_{m_{k}}^{n} \over 2}{\partial r_{m_{i}}^{n} \over \partial x_{k}}-{\delta ^{n} \over 2}{\partial r_{m_{i}}^{n} \over \partial t}{\bigg ]}$

Eq.(23) is now split into the two following equations

${\begin{aligned}u_{i}^{*}&=&u_{i}^{n}-{\Delta t \over \rho }{\bigg [}\rho v_{j}{\partial u_{i} \over \partial x_{j}}-{\partial \tau _{ij} \over \partial x_{j}}-{h_{m_{k}} \over 2}{\partial r_{m_{i}} \over \partial x_{k}}-{\delta \over 2}{\partial r_{m_{i}} \over \partial t}{\bigg ]}^{n}\\u_{i}^{n+1}&=&u_{i}^{*}-{\Delta t \over \rho }{\partial p^{n} \over \partial x_{i}}\end{aligned}}$

Note that the sum of eqs.(24) and (25) gives the original form of eq.(23).

Substituting eq.(25) into the stabilized mass balance equation (20) gives the standard Laplacian of pressure form

Error al representar (función desconocida «\eqno»): {\displaystyle \left({\Delta t \over \rho} + g_{ii}^n\right) {\partial^2 p^n\over \partial x_i \partial x_i}=r_d^*+r_p^n\eqno (26a)} where Error al representar (función desconocida «\eqno»): {\displaystyle r_d^*={\partial u_i^*\over \partial x_i}\eqno (26b)}

Standard fractional step procedures neglect the contribution from the terms involving $g_{ii}$ in eq. (26a). These terms have an additional stabilization effect which improves the numerical solution when the values of $\Delta t$ are small. Also the influence of the stabilization term $g_{ii}$ has proven to be essential for obtaining a fully converged solution in steady state problems (see the square cavity example in a next section). Indeed accounting for this additional stabilization term has lead to improved numerical solutions in all problems solved. Similar conclusions have been reached in a recent work by Codina [59].

Note that the cross-derivative terms $\displaystyle {\partial ^{2}p \over \partial x_{i}\partial x_{j}}$ have been kept within the term $r_{p}$ in the r.h.s. of eq.(26a). The influence of these terms should be studied in more detail in the future.

The stabilized free surface wave equation (7) is discretized in time to give

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

A typical solution in time includes the following steps. .3cm

Step 1. Solve explicitely for the so called fractional velocities $u_{i}^{*}$ using eq. (24).

.3cm

Step 2. Solve for the dynamic pressure field $p^{n}$ solving the Laplacian equation (26a). The dynamic pressures at the free surface computed from step 6 below, in the previous time step, are used as boundary conditions for solution of eq.(26a).

.3cm

Step 3. Compute the velocity field $u_{i}^{n+1}$ at the updated configuration for each mesh node using eq.(25)

.3cm Step 4. Compute the new position of the free surface elevation ${\beta }^{n+1}$ in the fluid domain by using eq.(27).

.3cm Step 5. Compute the movement of the submerged body by solving the dynamic equations of motion in the body subjected to the pressure field $p^{n}$ and the viscous stresses $\tau _{ij}^{n+1}$.

.3cm

Step 6. Compute the new position of mesh nodes in the fluid domain at time $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.

Assuming air is at rest, the absolute pressure at the free surface at time $n+1$ obtained from the stress equilibrium condition (neglecting surface tension effects) as Error al representar (función desconocida «\eqno»): {\displaystyle p_a=\tau_{33}\eqno (28a)}

The dynamic pressure at the free surface is computed by Error al representar (función desconocida «\eqno»): {\displaystyle {p}^{n+1}=\left[{\tau_{33}\over \rho} - g \beta\right]^{n+1}\eqno (28b)} where $g$ is the gravity constant.

As already mentioned, the effect of changes in the free surface elevation are introduced in step 2 of the flow solution as a prescribed dynamic pressure acting on the free surface. Note that eq.(28b) allows to take into account the changes in the free surface without the need of updating the free surface nodes. A higher accuracy in the solution of the flow problem can however be obtained if the free surface nodes are updated after a number of time steps.

Finite element discretization

Spatial discretization is carried out using the finite element method [15]. The stabilized formulation described allows an equal order interpolation of velocities and pressure [10,15]. A linear interpolation over four node tetrahedra for both $u_{i}$ and $p$ is chosen in the examples shown in the paper. Similarly, linear triangles are chosen to interpolate $\beta$ on the free surface mesh. The velocity and pressure fields are interpolated within each element in the standard finite element manner as Error al representar (función desconocida «\eqno»): {\displaystyle u_i=\sum\limits_j N_{j} \bar{(u_i)}_j\eqno (29a)} Error al representar (función desconocida «\eqno»): {\displaystyle p_i=\sum\limits_j N_{j}(\bar p_j)\eqno (29b)} where $N_{j}$ are the linear shape functions interpolating the velocity and pressure fields, respectively, and ${\bar {(\cdot )}}$ denote nodal values [15].

Similarly the wave height is discretized as

$\beta =\sum \limits _{j}N_{\beta _{j}}({\bar {\beta }}_{j})$

where $N_{\beta _{j}}$ are linear shape functions defined over the three node triangles discretizing the free surface.

The discretized integral form in space is obtained by applying the standard Galerkin procedure to eqs.(24),(25),(26a) and (27) and the boundary conditions (13). Solution of the discretized problem follows the pattern given below.

Step 1. Solve for the nodal fractional velocities

${\bar {\bf {u}}}^{*}={\bf {M}}^{-1}{\bf {f}}_{1}^{n}$

with ${M}_{ij}=\int _{\Omega }N_{i}N_{j}d\Omega$

${\bf {f}}^{n}=\left\{{\begin{array}{c}{\bf {f}}_{1}\\{\bf {f}}_{2}\\\vdots \\{\bf {f}}_{n}\end{array}}\right\}^{n}\qquad ,\qquad {\bf {f}}_{k}^{n}=\left\{{\begin{array}{c}f_{k_{1}}\\f_{k_{2}}\\f_{k_{3}}\end{array}}\right\}^{n}$

Error al representar (función desconocida «\nonumber»): {\displaystyle \begin{aligned} f_{k_i}^n&=&\int_\Omega N_{k} \bigg[u_i -{\Delta t\over \rho} \left( \rho v_j{\partial u_i\over \partial x_j} -{h_{m_k}\over 2} {\partial r_{m_i}\over \partial x_k}-{\delta\over 2}{\partial r_{m_i}\over \partial t}\right)\bigg]^nd\Omega +\nonumber\\ &&+\int_\Omega {\Delta t\over \rho} {\partial N_{k}\over \partial x_j} \tau^n_{ij} d\Omega -\int_{\Gamma_t} {\Delta t\over \rho} N_{k}t_i^n d\Gamma\quad,\quad i=1,2,3\end{aligned}}

The solution of eq.(31) can be speeded up by diagonalizing matrix ${\bf {M}}$. Alternatively a simple Jacobi iteration procedure can be used and this has proved to converge in very few iterations.

No boundary condition is applied when solving for the fractional velocities $u_{i}^{*}$ in eq.(31) as these velocities can be interpreted as a predicted value of the actual velocities. The kinematic boundary conditions (14) are applied in step 3 as shown below.

Step 2. Solve for the nodal pressures at time n

${\bf {H}}{\bar {\bf {p}}}^{n}={\bf {q}}^{n}$

$H_{kl}=\int _{\Omega }{\partial N_{k} \over \partial x_{i}}\left({\Delta t \over \rho }+g_{ii}^{n}\right){\partial N_{l} \over \partial x_{i}}d\Omega$

$q_{k}^{n}=\int _{\Omega }{\partial N_{k} \over \partial x_{i}}u_{i}^{*}d\Omega -\int _{\Omega }N_{k}r_{p}^{n}d\Omega +\int _{\Gamma }N_{k}{\rho \over \Delta t}\left({\Delta t \over \rho }+g_{ii}^{n}\right)u_{i}^{n}n_{i}d\Gamma$

The last integral in eq.(37) can be neglected in solid walls and stationary free surfaces where the normal velocity is zero.

Recall that the dynamic pressures computed from step 6 are used as a boundary condition for solution of eq.(35).

Step 3. Solve for the nodal velocities at time n+1

${\bar {\bf {u}}}^{n+1}={\bf {M}}^{-1}{\bar {\bf {f}}}^{n}$

where ${\bf {M}}$ is given by eq.(35) and

${\bar {f}}_{k_{i}}^{n+1}=\int _{\Omega }N_{k}{\bigg [}u_{i}^{*}-{\Delta t \over \rho }{\partial p^{n} \over \partial x_{i}}{\bigg ]}d\Omega$

The kinematic boundary conditions on the nodal velocities (eq.(14)) are imposed when solving eq.(38).

Step 4. Solve for the new free surface height at the time n+1

The new free surface elevation $\beta ^{n+1}$ in the fluid domain is computed as

Error al representar (función desconocida «\pmb»): {\displaystyle \bar{\pmb \beta}^{n+1}={\bf M}_\beta^{-1} {\bf s}^{n}}

with ${\bf {M}}_{\beta }=\int _{\Gamma _{\beta }}{\bf {N}}_{\beta }^{T}{\bf {N}}_{\beta }d\Gamma$ $s_{i}^{n+1}=\int _{\Gamma _{\beta }}N_{\beta _{i}}{\bigg [}\beta ^{n}-\Delta t\left(v_{k}^{n}{\partial \beta ^{n} \over \partial x_{k}}-v_{3}^{n}-{\gamma \over 2}{\partial r_{\beta }^{n} \over \partial t}\right){\bigg ]}d\Gamma +\int _{\Gamma _{\beta }}{h_{\beta _{j}} \over 2}{\partial N_{\beta _{i}} \over \partial x_{j}}r_{\beta }d\Gamma$

In the derivation of eq.(42) the assumption that $r_{\beta }=0$ at the boundary line of the free surface domain has been made.

Steps 5 and 6 follow the process described in the previous section.

Computation of the stabilization parameters

Accurate evaluation of the stabilization parameters is one of the crucial issues in stabilized methods. Most of 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 ${\bf {h}}_{m}$ has the direction of the velocity field [6,10]. This unnecessary restriction leads to instabilities when sharp layers transversal to the velocity direction are present. This additional deficiency is usually corrected by adding a shock capturing or crosswind stabilization term [36–38].

Let us first assume for simplicity that the stabilization parameters for the mass balance equations are the same as those for the momentum equations. This implies ${\bf {h}}_{m}={\bf {h}}_{d}$ -.2cm The problem remains now finding the value of the characteristic length vectors ${\bf {h}}_{m}$. Indeed, the components of ${\bf {h}}_{m}$ can introduce the necessary stabilization along both the streamline and transversal directions to the flow.

Excellent results have been obtained in all problems solved using linear tetrahedra with the same value of the characteristic length vector for the three momentum equations defined by Error al representar (función desconocida «\pmb»): {\displaystyle {\bf h}_{m_i}=h_s {{\bf u}\over {u}}+h_{c} {{\pmb \nabla} u\over \vert{\pmb \nabla}u\vert}\qquad i=1,2,3} where $u=\vert {\bf {u}}\vert$ and $h_{s}$ and $h_{c}$ are the “streamline” and “cross wind” contributions given by Error al representar (función desconocida «\pmb»): {\displaystyle \begin{aligned} h_s&=&\max ({\bf l}^T_j {\bf u})/{u} \\ h_{c}&=&\max ({\bf l}^T_j {\pmb \nabla}u)/ \vert {\pmb \nabla}u\vert \quad , \quad j=1,n_s\end{aligned}} where ${\bf {l}}_{j}$ are the vectors defining the element sides ($n_{s}=6$ for tetrahedra).

An alternative method for computing vector ${\bf {h}}_{m}$ in a more consistent manner is explained in the next section.

As for the free surface equation the following value of the characteristic length vector ${\bf {h}}_{\beta }$ has been taken Error al representar (función desconocida «\pmb»): {\displaystyle {\bf h}_\beta =\bar h_s {{\bf u}\over {u}}+\bar h_c {{\pmb\nabla}\beta \over \vert {\pmb\nabla}\beta\vert}}

The streamline parameter has been obtained by eq.(45) using the value of the velocity vector ${\bf {u}}$ over the 3 node triangles discretizing the free surface and $n_{s}=3$.

The cross wind parameter has been computed by

Error al representar (función desconocida «\pmb»): {\displaystyle \bar h_c = \max [{\bf l}_j^T {\pmb \nabla}\beta ] {1\over \vert {\pmb \nabla}\beta \vert} \quad ,\quad j=1,2,3}

Note that the cross-wind terms in eqs.(44) and (47) account for the effect of the gradient of the solution in the stabilization parameters. This is a standard assumption in most “shock-capturing” stabilization procedures [36–39].

Regarding the time stabilization parameters $\delta$ and $\gamma$ in eqs.(5) and (7) the value $\delta =\gamma =\Delta t$ has been taken for solution of the problems presented in the paper. A more consistent evaluation following the diminishing residual technique described next is described in [9] for transient advective-diffusive problems.

Computation of the characteristic length parameters via a diminishing residual procedure

The idea of this technique first presented in [6] and tested in [7–9,11] for advective-diffusive problems is the following. Let us assume that a finite element solution for the velocity and pressure fields has been found for a given mesh. The point wise residual of the momentum equation corresponding to this particular solution is (assuming $\delta =0$ in eq.(5))

$^{1}r_{m_{i}}=r_{m_{i}}-{1 \over 2}h_{m_{j}}{\partial r_{m_{i}} \over \partial x_{j}}$

The average residual over an element can be defined as

$^{1}r_{m_{i}}^{(e)}={1 \over \Omega ^{(e)}}\int _{\Omega ^{(e)}}{}^{1}r_{m_{i}}d\Omega$

Let us assume now that an enhanced numerical solution has been found for the same mesh and the same approximation (i.e., neither the number of elements nor the element type have been changed). This enhanced solution could be based, for instance, in a superconvergent recovery of derivatives [15,53,54]. The element residual for the enhanced solution is denoted by $^{2}r_{m_{i}}^{(e)}$. The element residuals must obviously tend to zero as the solution improves and the following condition must be satisfied

$^{1}r_{m_{i}}^{(e)}-^{2}r_{m_{i}}^{(e)}\geq 0$

The above equation applies for $^{1}r_{m_{i}}^{(e)}>0$. Clearly for $^{1}r_{m_{i}}^{(e)}<0$ the inequality in eq.(51) should be changed to $\leq 0$.

Substituting eq.(49) into (51) and applying the identity condition in eq.(51) gives the following system of equations for each element which unknowns are the characteristic length parameters for the element

${\bf {A}}{\bf {h}}_{m}^{(e)}={\bf {f}}$

with Error al representar (función desconocida «\nonumber»): {\displaystyle \begin{aligned} A_{ij} &=& 2\left[{{}^2\partial r_{m_i}^{(e)}\over \partial x_j} - {{}^1\partial r_{m_i}^{(e)}\over \partial x_j}\right]\\ \nonumber\\ f_i&=& ^2r_{m_i}^{(e)}-^1r_{m_i}^{(e)}\end{aligned}}

The following “adaptive” algorithm can be proposed for obtaining a stabilized solution:

1. Solve for the numerical values of nodal velocities and pressure, for an initial value ${\bf {h}}_{m}^{(e)}={}^{o}{\bf {h}}_{m}^{(e)}$. Compute $^{1}r_{m_{i}}^{(e)}$.
2. Evaluate the enhanced velocity and pressure fields. Compute $^{2}r_{m_{i}}^{(e)}$.
3. Compute the updated value of ${\bf {h}}_{m}^{(e)}$ solving eq.(52).
4. Repeat steps (1)–(3) until a stable solution is found.

The above strategy can be naturally incorporated into the transient solution scheme previously described by simply updating the value of ${\bf {h}}_{m}^{(e)}$ after the solution for each time step has been found.

The assumption ${\bf {h}}_{d}={\bf {h}}_{m}$ can be relaxed and an independent value of the characteristic length vector ${\bf {h}}_{d}$ for the mass balance equation can be found following a similar approach as described above for computing ${\bf {h}}_{m}$.

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 [41–51].

Chiandussi, Bugeda and Oñate [14] have recently proposed a simple method for the movement of mesh nodes ensuring minimum element distortion. The method is based on the iterative solution of a fictious linear elastic problem on the mesh domain. In order to minimize mesh deformation the “elastic” properties of each mesh element are appropiately selected so that elements suffering greater movements are stiffer. The basis of the method is given below.

Let us consider an elastic domain with homogeneous isotropic elastic properties characterized by the Young modulus ${\bar {E}}$ and the Poisson coefficient $\nu$. Once a discretized finite element problem has been solved using, for instance, standard $C^{o}$ linear triangles (in 2D) or linear tetraedra (in 3D), the principal stresses ${}^{1}\sigma _{i}$ at the center of each element are obtained as

${}^{1}\sigma _{i}={\bar {E}}[\varepsilon _{i}-\nu (\varepsilon _{j}+\varepsilon _{k})]\qquad i,j=1,2,3~~{\hbox{for 3D }}$

where $\varepsilon _{i}$ are the principal strains.

Let us assume now that a uniform strain field $\varepsilon _{i}={\bar {\varepsilon }}$ throughout the mesh is sought. The principal stresses are then given by

${}^{2}\sigma _{i}=E{\bar {\varepsilon }}(1-2\nu )\qquad i=1,2,3~~{\hbox{for 3D }}$

where $E$ is the unknown Young modulus for the element.

A number of criteria can be now used to find the value of $E$. The most effective approach found in [14] is to equate the element strain energies in both analysis. Thus

${\begin{aligned}U_{1}&=&{}^{1}\sigma _{i}\varepsilon _{i}={\bar {E}}[(\varepsilon _{1}^{2}+\varepsilon _{2}^{2}+\varepsilon _{3}^{2})-2\nu (\varepsilon _{1}\varepsilon _{2}+\varepsilon _{2}\varepsilon _{3}+\varepsilon _{1}\varepsilon _{3})]\\U_{2}&=&{}^{2}\sigma _{i}\varepsilon _{i}=3E{\bar {\varepsilon }}^{2}(1-2\nu )\end{aligned}}$

Equaling eqs.(57) and (58) gives the sought Young modulus $E$ as $E={{\bar {E}} \over 3{\bar {\varepsilon }}^{2}(1-2\nu )}[(\varepsilon _{1}^{2}+\varepsilon _{2}^{2}+\varepsilon _{3}^{2})-2\nu (\varepsilon _{1}\varepsilon _{2}+\varepsilon _{2}\varepsilon _{3}+\varepsilon _{1}\varepsilon _{3})]$

Note that the element Young modulus is proportional to the element deformation as desired. Also recall that both ${\bar {E}}$ and ${\bar {\varepsilon }}$ are constant for all elements in the mesh.

The solution process includes the following two steps. .2cm

Step 1. Consider the finite element mesh as a linear elastic solid with homogeneous material properties characterized by ${\bar {E}}$ and $\nu$. Solve the corresponding elastic problem with imposed displacements at the mesh boundary.

.2cm

Step 2. Compute the principal strains and the values of the new Young modulus in each element using eq.(59) for a given value of ${\bar {\varepsilon }}$. Repeat the finite element solution of the linear elastic problem with prescribed boundary displacements using the new values of $E$ for each element. .2cm

The movement of the mesh nodes obtained in the second step ensures a quasi uniform mesh distortion. Further details on this method including other alternatives for evaluating the Young modulus $E$ can be found in [14].

The previous algorithm for movement of mesh nodes is able to treat the movement of the mesh due to changes in position of fully submerged and semi-submerged bodies. Note however that if the floating body intersects the free surface, the changes in the analysis domain geometry can be very important. From one time step to other emersion or inmersion of significant parts of the body can occur.

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 (Figure 2).

The surface mapping technique used in this work is based on transforming the 3D curved surfaces into reference planes. 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 [5].

Examples

All the examples shown next have been solved in a standard PC Pentium II 450 Mhz with a memory of 128 Mb.

Example 1. Square cavity problem

The purpose of this example is to test the stabilized formulation presented in the solution of a standard benchmark problem solved by a number of authors [22,23,40,59]. Figure 3 shows the definition of the problem solved with an unstructured mesh of 7395 linear tetrahedra for a Reynolds number value of 1.

The steady-state solution was sought using the stabilized fractional step algorithm previously described. Results in Figure 4a,b are tabulated for the horizontal velocity along the vertical centerline of the mid-section and for vertical velocity and pressure along the horizontal centerline of the same section. Numerical results are fully stable and agree well with similar solutions reported in the mentioned reference. The effect of the stabilization term $g_{ii}$ in the pressure equation (see eq.(26a)) is seen clearly in Figure 4c. The curves in this figure show the convergence towards steady state of the $L\infty$ norm of the nodal pressures with time. The curve listed as “standard” is obtained neglecting the stabilization term $g_{ii}$ in eq.(26a), whereas the second curve shows the convergence when this term is taken into account. The difference between the two curves is noticeable as the error obtained with the fully stabilized solution is several orders of magnitude smaller than that obtained neglecting the term $g_{ii}$.

Example 2. Submerged NACA 0012 profile

A 2D submerged NACA0012 profile at $\alpha =5^{\circ }$ angle of attack is studied. This configuration was tested experimentally by Duncan [55] for high Reynolds numbers (Re=400000) and modelled numerically using the Euler equations by several authors [50,51,52,56]. The submerged depth of the airfoil is equal to the chord and this was used as the length (L) for normalizing the problem. The Froude number for all the cases tested was set to $Fr={{u} \over {\sqrt {gL}}}=0.5672$ where $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 of [55].

Example 3. Sphere falling in a tube filled with liquid

The movement of a sphere falling by gravity in a cylindrical tube filled with liquid is studied. The relationship between the diameters of the sphere and the tube is 1:4. The Reynolds number for the stationary speed is 100. The mesh has 85765 elements with 13946 nodes (Figure 6).

Figures 6 and 7 show the mesh deformation and contours of the mesh deformation and of the velocity in the domain for different times, respectively. The evolution of the falling speed is shown in Figure 7c. Note the good agreement with the so called Stokes velocity computed by equaling the weight of the sphere with the resistance to the movement of the sphere expressed in terms of the velocity. Obviously, this value is slightly greater than the actual one as frictional effects are neglected.

A similar problem for a much greater number of spheres has been solved by Johnson and Tezduyar [47].

Example 4. Movement of a submerged sphere in an open channel

Figure 8 shows the geometry of the channel and the position of the sphere of 2m diameter with a weight of 1000 N and a rotational inertia of 1000 kgm$^{2}$. A mesh of 19870 linear tetrahedra with 4973 nodes has been used for the analysis.

The problem has been analyzed for values of Reynolds number = 200 and Froude number = 0.71, corresponding to a velocity of 1m/s at the inlet.

It is assumed that the sphere can only move vertically and rotate around the global $y$ axes due to the forces induced by the fluid. The vertical displacement is constrained by a spring linking the sphere to the ground. An initial vertical velocity of 1m/s for the sphere has been taken.

Figure 8 shows a plot of the time evolution of the vertical displacement of the sphere. The contours of the velocity module in the fluid on two perpendicular planes at different times is shown in Figure 8b. The deformation of the free surface at $t=0.47$ s. and 3.16 s. is shown in Figure 8c.

Example 5. Interaction of a rigid vertical cylinder with a moving stream

The definition of the problem is clearly seen in Figure 9a. The cylinder diameter is 2 m and the stream speed is 1 m/s. The Froude and Reynolds numbers are 1.0 and 200, respectively. The walls of the cylinder are assumed to be rigid in this case. A mesh of 35567 tetrahedra and 4670 nodes is used for the analysis.

Figure 9b shows the contours of the velocity module and the vertical displacement in the mesh for a time $t=4.57$ s. Note the important deformation of the free surface in this problem.

Example 6. Wigley hull

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

The same configuration was tested experimentally in [57] and modelled numerically by several authors [50,51,52,58]. We use here an unstructured 3D finite element mesh of 65434 linear tetrahedra, with a reference surface of 7800 triangles, partially represented in Figure 10.

Figure 10 also shows the results of the viscous analysis of the Wigley model in three different cases ($L_{pp}=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. Secondly, the volume mesh was updated due to free surface movement, considering the model fixed. The third case corresponds to the analysis of a real free model including the mesh updating due to free surface evaluation and ship movement (sinkage and trim). A Smagorinsky turbulence model was used in all the cases.

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

Wigley Hull. Total resistance coefficient

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

In the study of the free model the numerical values of sinkage and trim were -0.1% and 0.035, respectively, while experiment gave -0.15% and 0.04.

Figure 10a shows the pressure distribution obtained near the Wigley hull for the free model. A number of streamlines have also been plotted in the figure. The obtained mesh deformation in this case is also presented in Figure 10b.

Comparisons of the obtained body wave profile with the experimental data for the free and fixed models are shown in Figure 10b. Significant differences are found close to stern in the case of the fixed model.

The free surface contours for the truly free ship motion are shown in Figure 10c.

Conclusions

The finite calculus method makes it possible to derive stabilized forms of the governing differential equations for a viscous fluid with a free surface. Solution of the new stabilized equations written in ALE form with a semi-implicit fractional step finite element method provides a straight-forward and stable algorithm for fluid-structure interaction analysis.

The mesh-moving scheme presented ensures minimum mesh distortion for large mesh displacements. The stabilized finite element method developed is adequate for solving large scale fluid-structure interaction problems in naval architecture and offshore engineering.

An academic version of the software developed using the formulation presented can be freely downloaded from www.cimne.upc.es/shyne.

Acknowledgements

Support for this work was provided by the European Community through projects Brite-Euram BR 967-4342 SHEAKS and Esprit 24903 FLASH.

The support of Empresa Nacional BAZAN de Construcciones Navales y Militares, S.A. is also gratefully acknowledged.

Thanks are given to Mr. J. Royo for his help in computing some of the examples presented.

The authors also thank Prof. S. Idelsohn, Prof. R. Lohner and Dr. C. Sacco for many useful discussions.