## A stabilized finite element method for analysis of fluid-structure interaction problems involving free surface waves

E. Oñate1,2, J. García-Espinosa1,2

1 Centre Internacional de Metodes Numerics a l'Enginyeria - CIMNE, Barcelona, Spain
2 Universitat Politècnica de Catalunya (UPC), Barcelona, Spain

## Summary

A stabilized semi-implicit fractional step finite element method for analysis of coupled fluid interaction problems involving free surface waves has been presented. A procedure for automatic movement of mesh nodes during the coupled solution process has been developed. The method is adequate for solving large scale fluid-structure interaction situations in naval architecture and offshare engineering problems.

## 1 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 if great relevance in 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 non linearities except for the simplest and limited potential flow model.

2. The obstacles in solving the constraint equation stating that at the free surface boundary the fluid particles remain on that surface which position is in turn unknown.

3. The difficulties in solving the problem of motion of the submerged body due to the interaction forces while minimizing the distorsion of the finite elements discretizing the fluid domain thus reducing the need of remeshing.

This paper presents a stabilized finite element method which allows to overcome 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 increment calculus (FIC) procedure developed by the authors [1–4]. 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–5].

The modified governing equations are solved in space-time using a semi-implicit fractional step approach and the finite element method (FEM). Free surface wave effects are accounted for via the introduction of a prescribed pressure at the free surface computed from the wave height.

The movement of the 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 for the fluid domain following the movement of the submerged body which minimizes element distorsion is presented. The mesh update procedure is based on the iterative finite element solution of a linear elastic problem on the mesh domain where fictitions elastic properties are assigned so that elements suffering higher movements are stiffer [6].

The content of the paper is structured as follows. First details of the stabilized semi-implicit fractional step approach using the FEM is described. Next the mesh updating procedure is presented. Finally some examples of a coupled fluid-interaction problem are given.

## 2 Stabilized finite element formulation for the fluid flow equations

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

The stabilized form of the governing differential equations for the three dimensional (3D) problem can be written as

Momentum

 ${\displaystyle r_{m_{i}}-{\underline {{1 \over 2}h_{mj}{\partial r_{m_{i}} \over \partial x_{j}}}}{1 \over 2}h_{mj}{\partial r_{m_{i}} \over \partial x_{j}}=0\qquad {\hbox{on }}\Omega \quad i,j=1,2,3}$
(1a)

Mass balance

 ${\displaystyle r_{d}-{\underline {{1 \over 2}h_{dj}{\partial r_{d} \over \partial x_{j}}}}{1 \over 2}h_{dj}{\partial r_{d} \over \partial x_{j}}=0\qquad {\hbox{on }}\Omega \quad j=1,2,3}$
(1b)

Free surface

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

where

 ${\displaystyle r_{m_{i}}=}$ ${\displaystyle \rho \left[{\partial u_{i} \over \partial t}+{\partial \over \partial x_{j}}(u_{i}u_{j})\right]+{\partial p \over \partial x_{i}}-{\partial \tau _{ij} \over \partial x_{j}}-b_{i}}$
(2a)
 ${\displaystyle r_{d}=}$ ${\displaystyle \rho {\partial u_{i} \over \partial x_{i}}\qquad i=1,2,3}$
(2b)
 ${\displaystyle r_{\beta }=}$ ${\displaystyle {\partial \beta \over \partial t}+u_{i}{\partial \beta \over \partial x_{i}}-u_{3}\qquad i=1,2}$
(2c)

In above ${\textstyle u_{i}}$ is the velocity along the ith global reference axis, ${\textstyle \rho }$ the (constant) density of the fluid, ${\textstyle p}$ the pressure, ${\textstyle \beta }$ the wave elevation, ${\textstyle b_{i}}$ the body forces acting in the fluid and ${\textstyle \tau _{ij}}$ the viscous stresses related to the viscosity ${\textstyle \mu }$ by the standard expression

 ${\displaystyle \tau _{ij}=\mu \left({\partial u_{i} \over \partial x_{j}}+{\partial u_{j} \over \partial x_{i}}-\delta _{ij}{2 \over 3}{\partial u_{k} \over \partial x_{k}}\right)}$
(3)

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

The distances ${\textstyle h_{mj},h_{dj}}$ and ${\textstyle h_{\beta _{j}}}$ are termed characteristic lengths and represent the dimensions of the finite domain where balance of momentum, mass and transport of fluid particles is enforced. Details of the derivation of eqs.(2) can be found in [1].

A more convenient form of equation (1b) can be written by assuming ${\textstyle h_{dj}=-2\tau _{d}u_{j}}$ where ${\textstyle \tau _{d}}$ is an intrinsic time parameter. Under this assumption and using eq.(1a) the stabilized form of the mass balance equation can be written as (neglecting high order terms) [4]

 ${\displaystyle r_{d}-\tau _{d}{\partial {\bar {r}}_{m_{i}} \over \partial x_{i}}=0}$
(4)

where

 ${\displaystyle {\bar {r}}_{m_{i}}=\rho \left({\partial u_{i} \over \partial t}+u_{j}{\partial u_{i} \over \partial x_{j}}\right)+{\partial p \over \partial x_{i}}-{\partial \tau _{ij} \over \partial x_{j}}}$
(5)

The boundary conditions for the stabilized problem are written as

 ${\displaystyle n_{j}\tau _{ij}+t_{i}+{\underline {{1 \over 2}h_{mj}n_{j}r_{m_{i}}}}{1 \over 2}h_{mj}n_{j}r_{m_{i}}=0\qquad {\hbox{on }}\Gamma _{t}}$
(6a)
 ${\displaystyle u_{j}-u_{j}^{p}=0\qquad {\hbox{on }}\Gamma _{u}}$
(6b)

where ${\textstyle n_{j}}$ are the components of the unit normal vector to the boundary and ${\textstyle t_{i}}$ and ${\textstyle u_{j}^{p}}$ are prescribed tractions and displacements on the boundaries ${\textstyle \Gamma _{t}}$ and ${\textstyle \Gamma _{u}}$, respectively. The underlined stabilized tersm appearing in the Neumann boundary condition (6a) are obtained via the FIC approach [1].

Eqs.(1–6) are the starting point for deriving a variety of stabilized numerical methods for solving the incompressible Navier-Stokes equations. 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 [4]. A semi-implicit fractional step finite element procedure for solution of eqs.(1a),(1c),(4) and (6) is presented in next section.

### 2.1 Stabilized fractional step method

Let us discretize in time the stabilized momentum equation (1a) as

 ${\displaystyle \rho \left[{u_{i}^{n+1}-u_{i}^{n} \over \Delta t}+{\partial \over \partial x_{j}}(u_{i}u_{j})^{n}\right]+{\partial p^{n+1} \over \partial x_{i}}-{\partial \tau _{ij}^{n} \over \partial x_{j}}-b_{i}^{n}-{1 \over 2}h_{mj}{\partial r_{m_{i}}^{n} \over \partial x_{j}}=0}$
(7)

A fractional step method (also termed “segregation” or “splitting” procedure) can be simply derived by splitting eq.(7) as follows

 ${\displaystyle u_{i}^{*}=}$ ${\displaystyle u_{i}^{n}-\Delta t\left[{\partial \over \partial x_{j}}(u_{i}u_{j})-{1 \over \rho }{\partial \tau _{ij} \over \partial x_{j}}-{1 \over \rho }b_{i}-{1 \over 2\rho }h_{mj}{\partial r_{m_{i}}^{n} \over \partial x_{j}}\right]^{n}}$
(8a)
 ${\displaystyle u_{i}^{n+1}=}$ ${\displaystyle u_{i}^{*}-{\Delta t \over \rho }{\partial p^{n+1} \over \partial x_{i}}}$
(8b)

Note that addition of eqs.(8) gives the original stabilized momentum equation (7).

Substitution of (8a) into eq.(4) gives after some algebra [4]

 ${\displaystyle (\Delta t+\tau _{d})\Delta p^{n+1}={\partial u_{i}^{*} \over \partial x_{i}}-\tau _{d}{\partial g_{i}^{n} \over \partial x_{i}}=0}$
(9)

where

 ${\displaystyle g_{i}=\rho \left({\partial u_{i} \over \partial t}+u_{j}{\partial u_{i} \over \partial x_{j}}\right)-{\partial \tau _{ij} \over \partial x_{j}}-b_{i}}$
(10)

Standard fractional step procedures neglect the contribution from the terms involving ${\textstyle \tau _{d}}$ in eq.(9). It can be shown that these terms have an additional stabilization effect which improves the numerical solution when the values of ${\textstyle \Delta t}$ are small.

The free surface wave equation (1c) can be also discretized in time to give

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

A typical solution in time includes the following steps.

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

Step 2. Solve for the pressure field ${\textstyle p^{n+1}}$ solving the Laplacian equation (9). The pressures at the free surface computed from step 4 below in the previous time step are used as boundary conditions for solution of eq.(9).

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

Step 4. Compute the new position of mesh nodes ${\textstyle {x}_{j}^{n+1}}$ in the fluid domain by using the mesh update algorithm described in next section.

Step 5. Compute the fractional velocity and pressure fields at the new position of the nodes. This can be simply done by the following expression

 ${\displaystyle ({\hat {u}}_{i}^{*})_{j}=}$ ${\displaystyle (u_{i}^{*})_{j}-[{d}_{j}^{n+1}]^{T}({\nabla {\hbox{ }}}{u}_{i}^{*})_{j}\quad i=1,2,3}$
(12a)
 ${\displaystyle {\hat {p}}_{j}^{n+1}=}$ ${\displaystyle p_{j}^{n+1}-[{d}_{j}^{n+1}]^{T}({\nabla {\hbox{ }}}{p}_{i}^{n})_{j}}$
(12b)

In eq.(12) ${\textstyle j}$ is the node number, ${\textstyle {\hat {(\cdot )}}_{j}}$ are values of ${\textstyle u_{i}^{*}}$ in the updated nodes after mesh movement whereas the r.h.s. includes values of ${\textstyle u_{i}^{*}}$ and ${\textstyle p}$ in the original position, ${\textstyle {d}_{j}=({x}_{j}^{n+1}-{x}_{j}^{n})}$ and ${\textstyle {\nabla {\hbox{ }}}}$ is the gradient vector.

An alternative to eq.(12b) is to compute the nodal pressures in the updated mesh configuration by solving once more eq.(9) using the values of ${\textstyle {\hat {u}}_{i}^{*}}$ obtained from eq.(12a). This will ensure better satisfaction of the incompressibility condition in the updated mesh.

Step 6. Compute the velocity field ${\textstyle u_{i}^{n+1}}$ at the updated configuration for each mesh node

 ${\displaystyle u_{i}^{n+1}={\hat {u}}_{i}^{*}-\Delta t{\partial {\hat {p}}^{n+1} \over \partial x_{i}}}$
(13)

Step 7. Solve for the updated value of the free surface elevation ${\textstyle \beta ^{n+1}}$ using eq.(11). Compute the pressure in the free surface from Benouilli equation as

 ${\displaystyle p^{n+1}=p^{o}+\rho g(\beta ^{n+1}-\beta ^{o})}$
(14)

where ${\textstyle \beta ^{o}}$ and ${\textstyle p^{o}}$ are reference values of the free surface elevation and the pressure respectively and ${\textstyle g}$ is the gravity constant.

As already mentioned the effect of changes in the free surface elevation is introduced in the step 2 of the flow solution as a prescribed pressure acting on the free surface.

The accuracy of above transient solution process depends on the time step size which should satisfy stability criteria for the coupled solution. Indeed larger time steps can be used if the values at time ${\textstyle n}$ in above equations are computed at ${\textstyle n+1/2}$. The solution process becomes now implicit and an iteration loop within each time step is then required.

### 2.2 Finite element discretization

Space discretization is carried out using the finite element method [7]. The velocity and pressure fields are interpolated within each element in the standard finite element manner as

 ${\displaystyle u_{i}=}$ ${\displaystyle \sum \limits _{j}N_{u_{j}}{\bar {(u_{i})}}_{j}}$ ${\displaystyle p_{i}=}$ ${\displaystyle \sum \limits _{j}N_{p_{j}}{\bar {p}}_{j}}$
(15)

where ${\textstyle N_{u_{j}}}$ and ${\textstyle N_{p_{j}}}$ are the shape functions interpolating the velocity and pressure fields, respectively and ${\textstyle {\bar {(\cdot )}}}$ denote nodal values.

Similarly the wave height is discretized as

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

where ${\textstyle N_{\beta _{j}}}$ are shape functions defined over the nodes discretizing the free surface.

It is worth noting that the stabilized formulation described allows an equal order interpolation of velocities and pressure [4]. A linear interpolation over triangles (2D) and tetrahedra (3D) for both ${\textstyle u_{i}}$ and ${\textstyle p}$ is chosen in the examples shown in the paper. Similarly linear elements are chosen to interpolate ${\textstyle \beta }$ on the free surface mesh.

The discretized integral form is obtained by applying the standard Galerkin procedure to eqs.(8a),(8b),(9),(11),(12) and the boundary conditions (6a). The resulting expressions follow the pattern given in [4].

The computation of the stabilization parameters ${\textstyle h_{m_{j}},h_{\beta _{j}}}$ and ${\textstyle \tau _{d}}$ can be based in the diminishing residual technique explained in [1–5]. In the example presented in Section 4 the simpler option ${\textstyle h_{m_{j}}=h_{\beta _{j}}=h^{(e)}}$ and ${\textstyle \tau _{d}={h^{(e)} \over 2\vert {u}\vert }}$ with ${\textstyle h^{(e)}=[V^{(e)}]^{1/3}}$ where ${\textstyle V^{(e)}}$ is the element volume has been taken.

## 3 A simple algorithm for stable updating of mesh nodes

Finite element solution of fluid-structure interaction problems usually requires the update of the analysis mesh as described in previous section. A typical example is the study of movement of an object within a flowing liquid where the fluid mesh needs to be continuously updated accordingly to the changes in position of the object due to the interaction forces.

Chiandussi, Bugeda and Oñate [6] have recently proposed a simple method for movement of mesh nodes ensuring minimum element distorsion. The method is based on the iterative solution of a fictitions 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 ${\textstyle {\bar {E}}}$ and the Poisson coefficient ${\textstyle \nu }$. Once a discretized finite element problem has been solved using, for instance, standard ${\textstyle C_{o}}$ linear triangles (in 2D) or linear tetraedra (in 3D), the principal stresses ${\textstyle {}^{1}\sigma _{i}}$ at the center of each element are obtained as

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

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

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

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

where ${\textstyle E}$ is the unknown Young modulus for the element.

A number of criteria can be now used to find the value of ${\textstyle E}$. The most effective approach found in [6] is to equal the element strain energies in both analysis. Thus

 ${\displaystyle U_{1}=}$ ${\displaystyle {}^{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})]}$
(19)
 ${\displaystyle U_{2}=}$ ${\displaystyle {}^{2}\sigma _{i}\varepsilon _{i}=3E{\bar {\varepsilon }}^{2}(1-2\nu )}$
(20)

Equaling eqs.(19) and (20) gives the sought Young modulus ${\textstyle E}$ as

 ${\displaystyle 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})]}$
(21)

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

The solution process includes the following two steps.

Step 1. Consider the finite element mesh as a linear elastic solid with homogeneous material properties characterized by ${\textstyle {\bar {E}}}$ and ${\textstyle \nu }$. Solve the corresponding elastic problem with imposed displacements at the mesh boundary. These displacements can be due to a prescribed motion of a body within a fluid, to changes in the shape of the domain in an optimum design problem, etc.

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

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

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 posible solution to this problem is to remesh the analysis domain. However for m ost problems, a mapping of the moving surfaces linked to mesh updating algorithm described above can avoid remeshing (Figure 1).

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

Figure 1: Changes in the fluid interface in a floating body.

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

Figure 2 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${\textstyle ^{2}}$. A mesh of 19870 linear tetrahedra with 4973 nodes has been used for the analysis.

Figure 2: Geometry of the chanel with submerged sphere.

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

It is assumed that the sphere can only move vertically and rotate a fluid around the global ${\textstyle 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 3: Time evolution of vertical displacement of sphere.

Figure 3 shows a plot of the time evolution of the vertical displacement of the sphere. The position of the sphere at different time intervals is shown in Figure 4. A plot of velocity vectors in the fluid displayed on Figure 5. Pressure contours in the fluid domain are shown in Figure 6.

Figure 4. Position of sphere and mesh at different times. a) t=0.47s, b) t=0.94s, c) t=1.22s, d) t=1.83s, e) t=2.57s, f) t=3.16s.
Figure 5: Velocity contours at different times. a) t=0.47s, b) t=0.94s, c) t=1.22s, d) t=1.83s, e) t=2.57s, f) t=3.16s.
Figure 6: Pressure distribution in fluid domain at different times. a) t=0.47s, b) t=0.94s, c) t=1.22s, d) t=1.83s, e) t=2.32s, f) t=2.85s.

The streamlines in the vecinity of the sphere at a certain time are shown in Figure 7. Figure 8 shows contours of the surface wave elevation at two particular times.

Figure 7: Streamline tracks at t=1.83s.
Figure 8: Free surface elevation at different times: a) t=0.47s, b) t=3.16s.

## 5 Conclusions

A stabilized semi-implicit fractional step finite element method for analysis of coupled fluid interaction problems involving free surface waves has been presented. A procedure for automatic movement of mesh nodes during the coupled solution process has been developed. The method is adequate for solving large scale fluid-structure interaction situations in naval architecture and offshare engineering problems.

## 6 Acknowledgements

This work was partially supported by Empresa Nacional BAZAN de Construcciones Navales y Militares S.A. This support is gratefully acknowledged.

## 7. References

1. E. Oñate - Derivation of stabilized equations for advective-diffusive transport and fluid flow problems. Comput. Meth. Appl. Mech. Engng., Vol. 151, 1-2, pp. 233–267, 1998.

2. E. Oñate, J. Garcia and S. Idelsohn - Computation of the stabilization parameter for the finite element solution of advective-diffusive problems. Int. J. Num. Meth. Fluids, Vol. 25, pp. 1385–1407, 1997.

3. E. Oñate, J. Garcia and S. Idelsohn - An alpha-adaptive approach for stabilized finite element solution of advective-diffusive problems with sharp gradients. New Adv. in Adaptive Comp. Met. in Mech., P. Ladeveze and J.T. Oden (Eds.), Elsevier, 1998.

4. E. Oñate - A finite element method for incompressible viscous flows using a finite increment calculus formulation. Research Report N. 150, CIMNE, Barcelona, January 1999.

5. E. Oñate and S. Idelsohn - A mesh free finite point method for advective-diffusive transport and fluid flow problems. Computational Mechanics, 21, 283–292, 1988.

6. Chiandusi, G., Bugeda, G. and Oñate, E. - A simple method for update of finite element meshes. Research Report 147, CIMNE, Barcelona, January 1999.

7. Zienkiewicz, O.C. and Taylor, R.C. - The finite element method, 4th Edition, Vol. 1, McGraw Hill, 1989.

Back to Top

### Document information

Published on 25/05/16

Licence: CC BY-NC-SA license

### Document Score

5

Views 244
Recommendations 0

### claim authorship

Are you one of the authors of this document?