A FEM fluid-structure interaction algorithm for analysis of the seal dynamics of a Surface-Effect Ship

Abstract

This paper shows the recent work of the authors in the development of a time-domain FEM model for evaluation of the seal dynamics of a surface effect ship. The fluid solver developed for this purpose, uses a potential flow approach along with a stream-line integration of the free surface. The paper focuses on the free surface-structure algorithm that has been developed to allow the simulation of the complex and highly dynamic behavior of the seals in the interface between the air cushion, and the water.

The developed fluid-structure interaction solver is based, on one side, on an implicit iteration algorithm, communicating pressure forces and displacements of the seals at memory level and, on the other side, on an innovative wetting and drying scheme able to predict the water action on the seals. The stability of the iterative scheme is improved by means of relaxation, and the convergence is accelerated using Aitken’s method.

Several validations against experimental results have been carried out to demonstrate the developed algorithm.

1. Introduction

The Transformable Craft (T-Craft) is a novel ship concept of the US Office of Naval Research, operative in multiple modes. T-Craft can deploy in an unloaded condition from the intermediate support base to the seabase, and then be used as a high speed connector to the shore, transporting wheeled and tracked vehicles through the surf zone and onto the beach.

T-Craft has been conceived as a Surface-Effect Ship (SES). A SES is a non-amphibious vehicle supported by an air cushion, with flexible seals at the bow and stern, and twin hulls, like a catamaran, at the sides. Due to the lack of air leakage at the craft sides, lift power can be reduced significantly compared with other type of Air-Cushion Vehicles (ACV). Also, it is possible to install conventional water propellers or waterjet propulsion, with rather smaller machinery space requirements compared to that for air propellers or fans used on ACVs. Furthermore, the SES can operate in modes of full displacement, partial air-cushion support, and full aircushion support.

Predicting the overall performance of a SES is of paramount importance to support the design phase, as the motion of the ship can be affected by the interaction between the air, the cushion, the ship structure, the seals, the sea waves and the sea bottom in the shallow water region. Different approaches with different types of complexity and accuracy have been taken to cope with this type of analyses.

In the last decade, there have been extensive applications of Navier-Stokes models to naval hydrodynamics problems. For example, Oñate and García-Espinosa [1] presented a stabilized FEM for fluid structure interaction with free surface. In [2] Löhner et al. developed a FEM capable of tracking violent free surface flows interacting with objects. Also García-Espinosa et al. [3] developed a new technique to track complex free surface shapes. More recently, in [4], an application for the calculation of the flow about a SES in still water, using a commercial Volume of Fluid model, has been presented. While, in [5], Mousaviraad et al. uses an URANS solver for evaluating the manoeuvring performance of a SES. While the outcome of the analyses is outstanding, the CPU-time reported in this paper, makes this model quite unaffordable for being used during design stages.

Actually, it is a common consensus that solvers based on the Navier-Stokes equations are too expensive computationally speaking when it comes to simulate unsteady naval hydrodynamics problems. These sorts of problems can be more efficiently calculated using potential flow theory. This approach, jointly with the Stokes perturbation approximation, is widely used for analysis of seakeeping problems [6]. In [7], Connell et al., uses a boundary-element time-domain potential flow solver to calculate the multi-body seakeeping behaviour of a T-Craft SES and a LMSR in different scenarios. While, in [8], the same computational solver is adapted to calculate the manoeuvre of a SES.

Despite the complexity of the above referred SES computational models, none of them takes into account the seal dynamics, or the effect of free surface-seal interaction. However, it is well known the relevance of this interaction in the unsteady dynamics of a SES [9][10]. The complexity of this phenomenon makes impossible to develop a theoretical background, and prompts many design parameters to be traditionally decided by empirical formulas [9]. Actually, only limited theoretical and computational models have been developed to analyze seal dynamics [11][12][13].

This paper shows an extension of the work presented by Serván-Camas and García-Espinosa [6] in the development of an efficient seakeeping solver. In particular, it is focused in the recent work of the authors in the development of a computational model for the analysis of the complex and highly dynamic behavior of the seals in the interface between the air cushion, and the water of a T-Craft [14]. The fluid solver developed for this purpose, uses a potential flow approach along with a stream-line integration of the free surface. While this approximation is much simpler than using RANS computations, significant outcomes can be obtained as well, allowing to significantly reducing computational time by 2 or 3 orders of magnitude even when computing on a regular desktop or laptop.

The developed fluid-structure interaction solver is based, on one side, on an implicit iteration algorithm, using a TCP/IP sockets link, able to communicate pressure forces and displacements of the seals at memory level and, on the other side, on an innovative wetting and drying scheme able to predict the water action on the seals.

The outline of this paper is as follows. In section 2, we summarize the mathematical modeling and the details of the integration of the fluid equations. Section 3 describes de structural solver. Section 4 introduces the wetting and drying scheme and the fluid-structure interaction solver. In section 5 different validation cases are presented. Finally, we arrive at some conclusions in section 6.

2. Hydrodynamic model

2.1 Governing equations

We consider the first order diffraction-radiation problem of a ship moving with forward speed ${\textstyle {\mathit {\boldsymbol {u}}}}$. We assume that the flow is incompressible and irrotational, so that the velocity field ${\textstyle {\boldsymbol {v}}}$, is derived from a potential function ${\textstyle \phi }$ as: ${\textstyle {\boldsymbol {\,v=}}\nabla \phi }$ . The following assumptions are made on the order of magnitude of the velocity components and free surface elevation ${\textstyle \eta }$ :

${\mathit {\boldsymbol {u}}}=O\left(1\right);\,{\phi }_{x}=O\left(1\right);\,{\phi }_{y}=O\left(1\right);\,{\phi }_{z}=O\left(\epsilon \right)$ ${\phi }_{\alpha \beta }=O\left(\epsilon \right);\,\,{\eta }_{x}=O\left(\epsilon \right);\,\,{\eta }_{y}=O\left(\epsilon \right)$

Based on the previous assumption, the governing equations for the first order diffraction-radiation wave problem on a frame of reference fixed on the ship are:

${\nabla }^{2}\phi =0$ ${\partial }_{t}\phi -{\mathit {\boldsymbol {u}}}\cdot {\nabla }_{h}\phi +{{\frac {1}{2}}\nabla }_{h}\phi \cdot {\nabla }_{h}\phi +$$\,{\frac {P}{\rho }}+g\eta =0$ ${\partial }_{t}\eta +\left({\mathit {\boldsymbol {-u}}}+{\nabla }_{h}\phi \right)\cdot {\nabla }_{h}\eta -$${\partial }_{z}\phi =0$ $\nabla \phi \cdot {\mathit {\boldsymbol {n}}}_{b}=-{\mathit {\boldsymbol {u}}}\cdot {\mathit {\boldsymbol {n}}}_{b}$ ${\partial }_{z}\phi =0$

in Ω in z = 0 in z = 0 in Γb in z = -H

(1)

where ${\textstyle {\nabla }_{h}\,}$ is the gradient in the horizontal plane, ${\textstyle \Omega \,}$ is the fluid domain, ${\textstyle {\Gamma }_{b}\,}$ represents the wetted surface of the ship, $H$ is the water depth, ${\textstyle {\mathit {\boldsymbol {n}}}_{b}}$ is the normal of the ship wetted surface pointing outwards the ship, and $P$ is the free surface pressure.

2.2 FEM solver

The finite element method, FEM, is adopted in this work to solve the system of equations (11-12). The FEM statement of the problem is as follows:

Let ${\textstyle {Q}_{h}^{\ast }}$ be the finite element space to interpolate functions, constructed in the usual manner. From this space, we can construct the subspace ${\textstyle {Q}_{h,\phi }}$, that incorporates the Dirichlet conditions for the potential ${\textstyle \,\phi }$ . The space of test functions, denoted by ${\textstyle {\,Q}_{h}}$, is constructed as ${\textstyle {\,Q}_{h,\phi }}$, but with functions vanishing on the Dirichlet boundary. The weak form of the problem can be formulated as:

Find ${\textstyle {\,[{\phi }_{h}]\epsilon Q}_{h}}$, by solving the discrete variational problem:

$\int _{\Omega }^{}\nabla {\upsilon }_{h}\cdot \nabla {\phi }_{h}d\Omega =$ $=\int _{{\Gamma }^{B}}^{}{\upsilon }_{h}\cdot {\hat {\phi }}_{n}^{B}d\Gamma +\int _{{\Gamma }^{{Z}_{0}}}^{}{\upsilon }_{h}\cdot {\hat {\phi }}_{n}^{{Z}_{0}}d\Gamma +\int _{{\Gamma }^{{Z}_{-H}}}^{}{\upsilon }_{h}\cdot {\hat {\phi }}_{n}^{{Z}_{-H}}d\Gamma \quad \,\forall {\upsilon }_{h}\in {Q}_{h}$

(2)

Where ${\textstyle {\hat {\phi }}_{n}^{B},\,{\hat {\phi }}_{n}^{{Z}_{0}}\,}$ and ${\textstyle {\hat {\phi }}_{n}^{{Z}_{-H}}}$ are the potential normal gradients corresponding to the Neumann boundary conditions on body, radiation boundary, free surface and bottom, respectively. At this point, it is useful to introduce the associated matrix form

${\overline {\overline {\mathit {\boldsymbol {L}}}}}{\mathit {\boldsymbol {\phi }}}={\mathit {\boldsymbol {b}}}^{B}+$${\mathit {\boldsymbol {b}}}^{\Gamma o}+{\mathit {\boldsymbol {b}}}^{{\Gamma }_{-H}}$

(3)

Where ${\textstyle {\overline {\overline {\mathit {\boldsymbol {L}}}}}}$ is the standard laplacian matrix, and ${\textstyle \,{\mathit {\boldsymbol {b}}}^{B}}$, ${\textstyle {\mathit {\boldsymbol {b}}}^{\Gamma o}\,}$ and ${\textstyle {\mathit {\boldsymbol {b}}}^{{\Gamma }_{-H}}}$ are the vector resulting of integrating the corresponding boundary condition term. Further details on the FEM solver can be found in [15].

2.3 Free surface boundary condition

Solving the free surface boundary condition efficiently is essential when dealing with free surface-structure interaction (FSSI) problems. The free surface conditions can be rewritten as:

${\partial }_{t}\phi +{\mathit {\boldsymbol {U}}}\cdot {\nabla }_{h}\phi -{{\frac {1}{2}}\nabla }_{h}\phi \cdot {\nabla }_{h}\phi +$$\,{\frac {P}{\rho }}+g\eta =0$ ${\partial }_{t}\eta +{\mathit {\boldsymbol {U}}}\cdot {\nabla }_{h}\eta -{\partial }_{z}\phi =$$0$

(4)

where ${\textstyle {\mathit {\boldsymbol {U=-u+}}}{\nabla }_{h}\phi }$ is the base flow. In this work no linearization is assumed since it is considered that ${\textstyle {\nabla }_{h}\phi =O(1)}$. Therefore, the convective velocity ${\textstyle {\mathit {\boldsymbol {U}}}}$ will be different in each time step, and the numerical scheme will be adapted accordingly. Three main issues must be kept in mind:

The numerical schemes adopted for solving the kinematic-dynamic free surface boundary conditions are based on an implicit midpoint rule for the kinematic condition, and backward Euler for the dynamic condition. Then ${\textstyle {\phi }^{n+1}}$ can be imposed as a Dirichlet boundary condition in the system (2). The schemes read as follows:

${\phi }^{n+1}={\phi }^{n}-\Delta t\,{\mathit {\boldsymbol {U}}}^{n+1}\cdot {\left({\nabla }_{h}\phi \right)}^{n+1}-$$\Delta t\cdot {\left({\nabla }_{h}\phi \right)}^{n+1}\cdot {\left({\nabla }_{h}\phi \right)}^{n+1}-$$\Delta t\,\left({P}^{n+1}/\rho +g{\eta }^{n+1}\right)=0$ ${\eta }^{n+1}={\eta }^{n}-\Delta t\,{\mathit {\boldsymbol {U}}}^{n+1/2}\cdot {\left({\nabla }_{h}\eta \right)}^{n+1/2}+$$\Delta t{\left({\partial }_{z}\phi \right)}^{n+1/2}=0$

(5)

It is well known that the convective terms appearing in the above equations can result in solutions corrupted by spurious node-to-node oscillations, if no special care is taken when evaluating them. In this work, the convective term is obtained by differentiating along streamlines:

${\mathit {\boldsymbol {U}}}^{n+1}\cdot {\left({\nabla }_{h}\phi \right)}^{n+1}=\left|{\mathit {\boldsymbol {U}}}^{n+1}\right|{\left({\partial }_{L}\phi \right)}^{n+1}$ ${\mathit {\boldsymbol {U}}}^{n}\cdot {\left({\nabla }_{h}\eta \right)}^{n}=\left|{\mathit {\boldsymbol {U}}}^{n}\right|{\left({\partial }_{L}\eta \right)}^{n}$

(6)

where ${\textstyle {\partial }_{L}}$ denotes the derivative along the streamline. The streamline derivatives are estimated using a three points upstream and one point downstream differential operator inspired by the work of [16]. The scheme has been conceived to work with unstructured meshes. The use of unstructured meshes enhances geometry flexibility and speed ups the initial modelling time. Furthermore it helps to reduce the number of elements by refining only in those areas of interest. An automatic unstructured grid generator based on the advancing front method is used to generate triangular surface grids and tetrahedral volumetric meshes.

Figure 1: Streamline discretization.

Figure 1 shows the tracing of the streamline at node C. The left (L) and forward left (FL) points are the upstream points, while the right (R) point corresponds to the downstream point. The values of the velocity potential ${\textstyle \phi }$ and free surface elevation ${\textstyle \eta }$ at L, FL and R points are obtained by linear interpolation between the nodes of the edges where they lie on. The streamline differential operator reads as:

${\partial }_{L}\phi \approx {\delta }_{L}\phi ={\alpha }_{R}{\phi }_{R}+{\alpha }_{C}{\phi }_{C}+$${\alpha }_{L}{\phi }_{L}+{\alpha }_{FL}{\phi }_{FL}$ ${\partial }_{L}\eta \approx {\delta }_{L}\eta ={\alpha }_{R}{\eta }_{R}+{\alpha }_{C}{\eta }_{C}+$${\alpha }_{L}{\eta }_{L}+{\alpha }_{FL}{\eta }_{FL}$

(7)

where the stencils are obtained as a one dimensional third order finite difference scheme:

$\left[{\begin{matrix}{\begin{matrix}1&\quad \,\,1\\{\Delta x}_{R}&\quad \,\,0\end{matrix}}&{\begin{matrix}1&1\\{-\Delta x}_{L}&\quad \,\,-{\Delta x}_{L}-{\Delta x}_{FL}\end{matrix}}\\{\begin{matrix}{\left({\Delta x}_{R}\right)}^{2}&0\\{\left({\Delta x}_{R}\right)}^{3}&0\end{matrix}}&\quad {\begin{matrix}{\left({\Delta x}_{L}\right)}^{2}&\,{\left({\Delta x}_{L}+{\Delta x}_{FL}\right)}^{2}\\{-\left({\Delta x}_{L}\right)}^{3}&-{\left({\Delta x}_{L}+{\Delta x}_{FL}\right)}^{3}\end{matrix}}\end{matrix}}\right]\left[{\begin{matrix}{\begin{matrix}{\alpha }_{R}\\{\alpha }_{C}\end{matrix}}\\{\begin{matrix}{\alpha }_{L}\\{\alpha }_{FL}\end{matrix}}\end{matrix}}\right]=$$\left[{\begin{matrix}{\begin{matrix}0\\1\end{matrix}}\\{\begin{matrix}0\\0\end{matrix}}\end{matrix}}\right]$

(8)

3. Structural solver

For the structural model of the seals it has been considered membrane and shell finite elements.

The membrane finite element used in this work is a 3-noded membrane triangle. This element uses a total Lagrangian formulation and assumes a linear-elastic behaviour of the material [17]. The kinematic description is formulated in large displacements and large deformations, making use of the right Cauchy-Green deformation tensor.

The shell finite element used in this paper is a corotational 3-noded triangle shell element [18]. The core formulation of the bending part corresponds to a classical Discrete Kirchhoff plate Triangle (DKT). The core formulation of the membrane part corresponds to the Assumed Natural Deviatoric Strain (ANDES) approach [19]. The ANDES formulation represents a linear strain triangle with 3 corner nodes and 3 DOFs per node (two in-plane displacements and a drilling rotation).

An energy conserving/decaying time-stepping algorithm has been used for the dynamic analysis of the problem [20] [21] [22]. This approach ensures the stability in the non-linear geometric range.

In general, shell elements are the preferred option to model the seals behaviour, because they can consider the actual flexural rigidity of the seals, despite being small. However, since the thickness of the seals is very small, numerical instabilities can appear. This happens because the tangent matrix is poorly conditioned by the large ratio between membrane and flexural rigidity of the seals. In these cases, the only solution is to neglect the flexural rigidity by using membrane elements. In a dynamic analysis the membrane elements no longer have these problems since the tangent stiffness matrix is properly conditioned by the mass matrix.

4. Coupled fluid-structure solver

4.1 Free surface-seal interaction (FSSI)

Modelling the behaviour of the seals of aircushion vehicles is not a trivial issue, due to the complex interaction of the seals with free surface. In this work, a new algorithm for handling the free surface-seal interaction problem is formulated. It is based on finding an equivalent pressure field to be applied over the free surface such that the elevation of this one is limited by the location of the seal. That is to say, the seal act as an upper limit for the free surface elevation.

The free surface boundary conditions are applied in different ways depending on whether the free surface is in contact with the seal or not. If the free surface is not in contact, the boundary conditions are applied as if there is no seal, but if there is contact, the implementation will be different in order to ensure that the free surface does not penetrate the seal and the necessary pressure to fulfil this condition is calculated, and imposed as boundary condition on the free surface.

It will be said that the free surface node where the algorithm is to be applied is dry if the seal is not in contact with the free surface at that location, and wet if it is. Figure 2 illustrates the wet and dry concepts.

The main challenge for an algorithm like this is to be capable of capturing when a node goes from dry to wet and vice versa, as well as estimating the pressure field on the wet nodes. Figure 3 shows the wetting and drying iterative algorithm. For a dry node, the implementation of both, the kinetic and dynamic boundary condition is the same as for any other node not interacting with the seal. However, for wet nodes, the free surface boundary condition is imposed via imposing that the free surface elevation matches the seal elevation, ${\textstyle \,{H}_{seal}}$, and ensuring that there is no flow across the seal. These two conditions are represented by the following equations:

${\eta }_{wet}^{n+1}={H}_{seal}^{n+1}$ ${\phi }_{z}^{n+1}=2\left({\frac {{\eta }^{n+1}{-\eta }^{n}}{\Delta t}}+\,\left|{\mathit {\boldsymbol {U}}}^{n+1/2}\right|{\left({\delta }_{L}\eta \right)}^{n+1/2}\right)-$${\phi }_{z}^{n}$

(9)

where ${\textstyle {\left({\partial }_{z}\phi \right)}^{n+1}}$ is to be imposed as Neumann condition instead of imposing ${\textstyle {\phi }^{n+1}}$ as Dirichlet condition.

Figure 2: Wet and dry seal regions for free surface boundary condition implementations.

The change from being a dry node to become a wet node is identified via the kinematic BC though the condition ${\textstyle {\eta }_{wet}^{n+1}>{H}_{seal}^{n+1}\,}$ . On the other hand, the switch from being a wet node to become a dry one is carried out by comparing the dynamic pressure with the reference pressure. The dynamic pressure on wet nodes is obtained by the dynamic BC:

${P}_{wet}^{n+1}=-\rho \left[{\frac {{\phi }^{n+1}{-\phi }^{n}}{\Delta t}}+\,\left|{\boldsymbol {U}}^{n+1}\right|\cdot {\left({\delta }_{L}\phi \right)}^{n+1}-\right.$$\left.{\frac {1}{2}}{\left({\nabla }_{h}\phi \right)}^{n+1}\cdot {\left({\nabla }_{h}\phi \right)}^{n+1}+\right.$$\left.{\frac {1}{2}}{\phi }_{z}^{n+1}{\phi }_{z}^{n+1}+\,g{\eta }^{n+1}\right]$

(10)

We found that when the seal is interacting with the free surface, the term ${\textstyle {\phi }_{z}}$ might become of the same order of magnitude than ${\textstyle {\boldsymbol {U}}}$. Hence the term ${\textstyle {\frac {1}{2}}{\phi }_{z}^{2}}$has been kept when calculating the pressure on a free surface node in contact with the seal.

Figure 3: Wetting and drying algorithm.

4.2 Free surface-flexible seals coupling algorithm

The fluid-structure coupling is performed by an implicit interfield iteration algorithm [23], based on a Block Gauss-Seidel method [24]. For the i-th iteration of the current time step n, let ${\textstyle {\boldsymbol {\phi }}^{n,i}}$, ${\textstyle {\boldsymbol {\eta }}^{n,i}}$, ${\textstyle {\boldsymbol {p}}^{n,i}}$ be the vectors of the velocity potential, free surface elevation and pressure in the corresponding nodes of the fluid, ${\textstyle {\boldsymbol {x}}^{n,i}}$ the vector containing the three components of the displacement of the nodes of the seals structure, with ${\textstyle {\boldsymbol {x}}_{k}^{n,i}}$ the vector of the k-th component of the displacement field, and ${\textstyle {\boldsymbol {H}}_{seal}^{n,i-1}}$ the vector of the elevation field of the seal nodes, then the fluid-structure interaction algorithm, can be written as:

${\begin{matrix}\left({\boldsymbol {p}}^{n,i},\,{\boldsymbol {\phi }}^{n,i}\right)={\boldsymbol {F}}^{p}\left({\boldsymbol {\phi }}^{n,i-1},{\boldsymbol {\eta }}^{n,i-1},{\boldsymbol {H}}_{seal}^{n,i-1},{{\boldsymbol {\eta }}^{n-1,i-1}{\boldsymbol {,\phi }}}^{n-1}\right)\\{\boldsymbol {x}}^{n,i}={\boldsymbol {S}}^{q}\left({\boldsymbol {x}}^{n,i-1},{\boldsymbol {p}}^{n,i},{\boldsymbol {x}}^{n-1},{\boldsymbol {x}}^{n-2}\right)\\{\boldsymbol {H}}_{seal}^{n,i-1}={\boldsymbol {H}}_{seal}^{0}+{\boldsymbol {x}}_{k}^{n,i}\end{matrix}}$

(11)

Where ${\textstyle {\boldsymbol {F,\,S}}}$ represent the fluid and solid solver, and p and q are the number of internal iterations done within the non-linear loop of each solver. This way, the pressure field computed in the hydrodynamic solver is sent to the structural solver to compute the seals deformation. The resulting displacements are used to compute the new seal elevation on the fluid interface, and then calculate the pressure field for the following iteration.

Figure 4: Fluid-structure interaction algorithm.

Since the hydrodynamics and structural solvers used in this work are independent, the strategy developed to communicate both is based on the interchange of information at memory level by means of TCP-IP sockets. This way, only minor adaptations of the solvers were required. The required functionality has been implemented in C++, in form of library. The communication protocol developed for the library is quite simple. Actually, it just requires sending and receiving functions for vectorial fields, including an implicit waiting instruction in the receiving function for synchronization purposes. Furthermore, this library has been extended to be able to interpolate data between the boundary meshes of the fluid and solid interfaces (see Figure 4). The interpolating algorithm is based on an octtree search to find the interpolating element of every node of the other mesh.

Let ${\textstyle {\boldsymbol {I}}^{FS}{\boldsymbol {,\,I}}^{SF}}$ be the interpolator operators from the fluid boundary mesh to the solid interface, and vice-versa, then the fluid-structure interaction algorithm reads

${\begin{matrix}\left({\boldsymbol {p}}^{n,i},\,{\boldsymbol {\phi }}^{n,i}\right)={\boldsymbol {F}}^{p}\left({\boldsymbol {\phi }}^{n,i-1},{\boldsymbol {\eta }}^{n,i-1},{\boldsymbol {H}}_{seal}^{n,i-1},{{\boldsymbol {\eta }}^{n-1,i-1}{\boldsymbol {,\phi }}}^{n-1}\right)\\{\check {\boldsymbol {p}}}^{n,i}={\boldsymbol {I}}^{FS}\left({\boldsymbol {p}}^{n,i}\right)\\{\boldsymbol {x}}^{n,i}={\boldsymbol {S}}^{q}\left({\boldsymbol {x}}^{n,i-1},{\check {\boldsymbol {p}}}^{n,i},{\boldsymbol {x}}^{n-1},{\boldsymbol {x}}^{n-2}\right)\\{\check {\boldsymbol {x}}}^{n,i}={\boldsymbol {I}}^{SF}\left({\boldsymbol {x}}^{n,i}\right)\\{\boldsymbol {H}}_{seal}^{n,i-1}={\boldsymbol {H}}_{seal}^{0}+{\check {\boldsymbol {x}}}_{k}^{n,i}\end{matrix}}$

(12)

Even though the iterations within each time step lead to the monolithic solution of the coupled problem, sometimes relaxation techniques are combined with this method to ensure the convergence of the coupling scheme. In this work, a relaxation method based on Aitken’s acceleration has been adopted [25]. Let ${\textstyle {\omega }_{i}}$ be the optimal relaxation parameter, obtained via Aitken’s method, the resulting algorithm can be written as:

${\begin{matrix}\left({\boldsymbol {p}}^{n,i},\,{\boldsymbol {\phi }}^{n,i}\right)={\boldsymbol {F}}^{p}\left({\boldsymbol {\phi }}^{n,i-1},{\boldsymbol {\eta }}^{n,i-1},{\boldsymbol {H}}_{seal}^{n,i-1},{{\boldsymbol {\eta }}^{n-1,i-1}{\boldsymbol {,\phi }}}^{n-1}\right)\\{\hat {\boldsymbol {p}}}^{n,i}=\left(1-{\omega }_{i}\right){\hat {\boldsymbol {p}}}^{n,i-1}+{\omega }_{i}{\boldsymbol {p}}^{n,i}\\{\check {\boldsymbol {p}}}^{n,i}={\boldsymbol {I}}^{FS}\left({\hat {\boldsymbol {p}}}^{n,i}\right)\\{\boldsymbol {x}}^{n,i}={\boldsymbol {S}}^{q}\left({\boldsymbol {x}}^{n,i-1},{\check {\boldsymbol {p}}}^{n,i},{\boldsymbol {x}}^{n-1},{\boldsymbol {x}}^{n-2}\right)\\{\check {\boldsymbol {x}}}^{n,i}={\boldsymbol {I}}^{SF}\left({\boldsymbol {x}}^{n,i}\right)\\{\boldsymbol {H}}_{seal}^{n,i-1}={\boldsymbol {H}}_{seal}^{0}+{\check {\boldsymbol {x}}}_{k}^{n,i}\end{matrix}}$

(13)

As stated above, the free surface-seals interaction algorithm is based on finding an equivalent pressure field to be applied over the free surface such that the elevation of this one is limited by the location of the seal. For this purpose, the elevation of the seal over the reference free surface is required. The algorithm (3) assumes that the variation of this elevation is only due to the vertical displacements of the structure. This is a good enough approximation to predict the forces of the water on the seals and the free surface elevation in the seals region.

5. Validation of FSSI algorithm with flexible bow seals

5.1 Validation against experiments at University of Michigan

Zalek and Doctors [26] carried out at University of Michigan an extensive experimental study of free surface and flexible seal interaction. Some of this data has been already used for validation purpose against numerical results based on a smoothed particle hydrodynamics approach [27]. Figure 5 sketches the basic arrangements of the experiments, such as seal length, location of stiffeners and cables and undeflected seal angle. Table 1 provides the seal properties used for the experiments [27].

Figure 5: Experiments setup at University of Michigan

Thickness	3.175 mm
Young Modulus	12.74 MPa
Poisson coefficient	0.25
Density	1107 kg/m3

Table 1: Seal properties

In this work, the experimental cases shown in Table 2 were analysed with the proposed fluid-structure interaction model. The test runs were simulated modelling the left half part of seal and tank, using the actual width and depth of the channel of the University of Michigan. The seals were modelled with shells elements, with the properties provided in Table 1.

Test run [26]	Seal Immersion		Velocity		Internal Pressure
	(in)	(m)	(ft/s)	(m/s)	(in. H2O )	(Pa)
1	9	0.2286	6	1.8288	1.42	353.3
10	9	0.2286	8	2.4384	1.01	251.3
11	9	0.2286	8	2.4384	1.91	475.3
16	9	0.2286	8	2.4384	4.24	1055
32	7.5	0.1905	6	1.8288	1.28	318.5
45	7.5	0.1905	9	2.7432	4.01	997.8

Table 2: Test run particulars

Table 3 provides the main characteristics of the computational mesh used for the simulation, as well as the length of the numerical channel used. For test run 45, a larger channel was necessary to avoid upstream and downstream interactions with the inlet and outlet boundaries respectively.

Test run	Seal		Water volume		Time step (s)	Domain length (m)
	Number of Nodes	Number of Elements	Number of Nodes	Number of Tetrahedra
1	234	424	58136	297188	0.002	5.3
10	234	424	58315	297888	0.001	5.3
11	234	424	58315	297888	0.001	5.3
16	234	424	58240	297643	0.002	5.3
32	234	424	56112	287493	0.002	5.3
45	234	424	80325	419569	0.001	15.3

Table 3: Computational domain and numerical particulars

In the experiments, three stiffeners (rods) were fixed to the seal to reduce the transversal deformation of the seals. Since we had no data regarding the properties of such stiffeners, we simulated them by introducing very stiff beams at the same locations, but weightless, assuming that the ones used in the experiments were light enough to avoid significant perturbation in the seal behaviour. Furthermore several cables were attached to the stiffeners rods to avoid the forward displacement of the seal due to the internal pressure. Since there is no available information about the length and exact disposition of these cables, they were simulated in the computational model by a boundary condition, constraining the longitudinal displacement of the points of the seal where the cables were attached.

Figure 6 shows the quasi-steady deformed front seal and free surface elevation respectively for test run number 10. Figure 7 compares de average seal deformation obtained numerically and experimentally. Although the computational results do not match exactly the experimental ones, they agree well in terms of the main trends under different conditions of velocity, seal immersion and internal pressure. Anyhow, the differences are probably within the uncertainty range of the data, since the unsteady and three-dimensional effects are relevant in most of the studied cases.

Figure 6: Seal deformation and structural mesh for test run 10.

Figure 7: Seal deformation: experimental results (squares) versus numerical (solid line).

Figure 8 provides the top view of the hydrodynamic pressure distribution acting on the seal due to the free surface seal interaction. It can be observed that while in some cases the distribution is almost uniform in the transverse direction it is not in others, resulting in relevant three dimensional effects, as mentioned before.

Figure 8: Top view of pressure distribution over seal

5.2 Validation against experiments at the U.S. Navy´s Large Cavitation Channel (LCC)

Wiggins et al. [28] developed a large scale testing platform for finger bow seals at the U.S. Navy´s Large Cavitation Channel. The purpose of those test were to provide experimental results to be used for validation purposes of analytical/numerical approaches to the problem. Figure 9 sketches the model setup at the LC facilities

Figure 9: LCC model setup sketch

The seal material used is the same one that the one used in the University of Michigan experiments and described in the previous section. In this case, a membrane model was used for the computational analysis of the structure, since preliminary tests with a shell model showed an unstable behavior, probably due to the poor condition of the tangent matrix.

With regards to the test condition, we used a free stream flow of 6 m/s, seal immersion of 15 cm, and internal pressure of 1.5 kPa [28]. Other data needed for the setup of the case under study are provided in Table 4.

Breadth of cushion (m)	1.52
Sidewall depth (m)	1.68
Sidewall draft at zero immersion (m)	0.42
Channel width (m)	3.05
Number of Fingers 5	5
Finger height (m)	1.00
Finger width (m)	0.31
Bow seal angle (º)	50.00

Table 4: Tests particulars [28]

Table 5 provides the main characteristics of the computational mesh used for the simulation, as well as the length of the numerical channel used. Figure 10 shows the computational mesh used for solving the fluid problem.

Fingers	Number of Nodes	7023
	Number of Elements	13436
Water domain	Number of Nodes	120942
	Number of Tetrahedra	659687
Time step (s)		0.0005
Domain length (m)		21

Table 5: Computational domain and numerical particulars

Figure 10: Computational mesh for potential flow solver

Figure 11 (a) shows the steady seal deformation obtained in the experiments under the abovementioned test conditions. Figure 11 (b), (c) and (d) show seal deformation obtained by the numerical approach proposed in this work under the same test conditions. Free surface elevation at the seal and aircushion areas are shown in Figure 12 (a). A top view of the pressure distribution over finger seals is shown in Figure 12 (b).

Figure 11: Experimental results (a) versus numerical results (b), (c) and (d). Picture obtained in the experimental tests reported in [28] (a); Front view of deformed fingers (b); deformed finger detail (c); top view (d).

Figure 12: Numerical results: free surface elevation (a) and top view of pressure distribution on fingers (b)

Next, a video of the T-Craft analysis using the computational model developed in this work is presented.

Video 1. T-Craft analysis using the computational model developed in this work.

As can be seen, the bow seals deformation obtained in this analysis agrees with the qualitative experimental information available in reference [28]. In fact, the deformation mode obtained in the computational model is quite similar to that found experimentally.

6. Conclusions

This paper reports the advances in the development of a fluid-structure interaction algorithm for evaluation of a SES seal dynamics.

The fluid solver developed for this purpose, uses a potential flow approach along with a stream-line integration of the free surface. The work has been focused on the FSSI algorithm that has been developed to allow simulating complex and highly dynamic behavior of the seals in the interface between the air cushion, and the water.

The developed fluid-structure interaction solver is based, on one side, on an implicit iteration algorithm, communicating pressure forces and displacements of the seals at memory level and, on the other side, on an innovative wetting and drying scheme able to predict the water action on the seals. The seal dynamics are solved using a Lagrangian FEM formulation for a 3-noded membrane/shell triangle.

Special care has to be taken to ensure the convergence of the fluid-structure interaction algorithm due to the highly flexible and non-linear nature of the seal. We found that relaxation method based on Aitkens acceleration methods worked well under all the case studies analysed.

Different validation and demonstration cases have confirmed the validity of the approach to study the complex dynamic behaviour of the seals.

7. Acknowledgements

This work relates to Department of the Navy Grant N62909-10-1-7053 issued by Office of Naval Research Global. The United States Government has a royalty-free license throughout the world in all copyrightable material contained herein.

8. References

[1] E. Oñate, J. García-Espinosa, A finite element method for fluid-structure interaction with surface waves using a finite calculus formulation, Comp. Methods Appl. Mech. and Eng. 2001; 191: 635-660.

[2] R. Löhner, C. Yang, E. Oñate, On the simulation of flows with violent free surface motion and moving objects using unstructured meshes, Comp. Methods Appl. Mech. Engng. 2007; 53: 1315-1338.

[3] J. García-Espinosa, A. Valls, E. Oñate, ODDLS: A new unstructured mesh finite element method for the analysis of free surface flow problems, Int. J. Numer.  Meth.  Fluids 2008; 76 (9): 1297-1327.

[4] D. J. Donnelly, W. L. Neu, Numerical Simulation of Flow About a Surface-Effect Ship, 11th International Conference on Fast Sea Transportation FAST 2011, Honolulu, Hawai. September 2001.

[5] S. M. Mousaviraad, S. Bhushan, F. Stern, CFD Prediction of Free-Running SES/ACV Deep and Shallow Water Maneuvering in Calm Water and Waves, MARSIM 2012. Singapore, April 23-27, 2012.

[6] B. Serván-Camas, J. García-Espinosa, Accelerated 3D multi-body seakeeping simulations using unstructured finite elements, J. Comp. Phys. 252 (2013) 382–403.

[7] B. S. H. Connell, W. M. Milewski, B. Goldman, D. C. Kring, Single and Multi-Body Surface Effect Ship Simulation for T-Craft Design Evaluation, 11th International Conference on Fast Sea Transportation FAST 2011, Honolulu, Hawai. September 2001.

[8] D. C. Kring, M. K. Parish, W. M. Milewski, B. S. H. Connell, Simulation of Maneuvering in Waves for a High-Speed Surface Effect Ship, 11th International Conference on Fast Sea Transportation FAST 2011, Honolulu, Hawai. September 2001.

[9] L. Yun, A. Bliault, Theory &Design of Air Cushion Craft, Elsevier 2005.

[10] N. Hirata, O. M. Faltinsen, Computation of Cobblestone effect with unsteady viscous flow under a stern seal bag of a SES, Journal of Fluids and Structures, 2000: 14, 1053–1069.

[11] L. J. Doctors, Nonlinear motion of an air-cushion vehicle over waves, Journal of Hydronautics 1975: 9 (2), 44–57.

[12] P. A. Sullivan, P. A. Charest, T. Ma, Heave stiffness of an air cushion vehicle bag and finger skirt, Journal of Ship Research 1994: 38 (4), 302–307.

[13] Q. Yang, V. Jones, L. McCue, Investigation of Skirt Dynamics of Air Cushion Vehicles under Non-linear Wave Impact Using a SPH-FEM Model, 11th International Conference on Fast Sea Transportation FAST 2011, Honolulu, Hawaii. September 2001.

[14] Evaluating Manoeuvering and Seakeeping Performance of a Surface Effect Ship, Technical Report, CIMNE IT-630, October 2012.

[15] O. C. Zienkiewicz, R. L. Taylor, J. Z. Zhu, The Finite Element Method: Its Basis and Fundamentals, Butterworth–Heinemann, ISBN 0750663200, 2005.

[16] S. R. Idelsohn, E. Oñate, C. Sacco, Finite element solution of free-surface ship-wave problems, Int. J. Numer. Meth. Engng. 45, 50-528 (1999)

[17] R. L. Taylor, E. Oñate, P.-A. Ubach, Finite element analysis of membrane structures, Textile composites and inflatable structures. Springer (2005) 47-68.

[18] C. Felippa, B. Haugen, A unified formulation of small-strain corotational finite elements: I. Theory, Computer Methods in Applied Mechanics and Engineering 194 (2005) 2285–2335

[19] C. A. Felippa, C. Militello, Membrane triangles with corner drilling freedoms: II. The ANDES element, Finite Elem Anal Des, (1992) 12:189–201.

[20] M. A. Crisfield, Non-linear Finite Element Analysis of Solid and Structures, Advanced Topics, vol. 2, Wiley, 1997.

[21] F. S. Almeida, A. M. Awruch, Corotational nonlinear dynamic analysis of laminated composite shells, Finite Elements in Analysis and Design 47 (2011) 1131–1145.

[22] A. Ibrahimbegovic, S. Mamouri, Energy conserving/decaying implicit time-stepping scheme for nonlinear dynamics of three-dimensional beams undergoing finite rotations, Computer Methods in Applied Mechanics and Engineering 191 (2002) 4241–4258.

[23] C. A. Felippa, K. C. Park, C. Farhat, “Partitioned analysis of coupled mechanical systems,” Comput. Methods Appl. Mech. Eng., vol. 190, pp. 3247–3270, 2001.

[24] J. G. Valdés, “Nonlinear analysis of orthotropic membrane and shell structures including fluid-structure interaction,” PhD thesis, Barcelona Tech., 2007.

[25] B. M. Irons, R. C. Tuck, “A Version of the Aitken Accelerator for Computer Iteration,” Int. J. Numer. Methods Eng., vol. 1, pp. 275–277, 1969.

[26] S. F. Zalek , L. J. Doctors, Experimental Study of the Resistance of Surface-Effect-Ship Seals, Proceeding of the 28th Symposium on Naval Hydrodynamics, 2010.

[27] R. J. Paredes, Smoothed particle hydrodynamics applied to fluid structure interaction problems involving hydroelastic response, Ph.D. thesis, Stevens Institute of Technology, United States of America, 2013.

[28] A. D. Wiggins, S. F. Zalek, M. Perlin, S. L. Ceccio, Development of large scale surface effect ship bow seal testing plataform. 11th International Conference on Fast Sea Transportation FAST 2011, Honolulu, Hawai. September 2001.