# Added mass evaluation with a finite-volume solver for applications in fluid-structure interaction problems solved with co-simulation

Camille Yvin, Alban Leroyer, Michel Visonneau, P. Queutey

► To cite this version:

Camille Yvin, Alban Leroyer, Michel Visonneau, P. Queutey. Added mass evaluation with a finite-volume solver for applications in fluid-structure interaction problems solved with co-simulation. Journal of Fluids and Structures, 2018, 81, pp.528-546. 10.1016/j.jfluidstructs.2018.05.008. hal-02391254

HAL Id: hal-02391254

https://hal.science/hal-02391254

Submitted on 26 Apr 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Added mass evaluation with a finite-volume solver for
applications in fluid-structure interaction problems solved
with co-simulation

C. Yvin<sup>a,b</sup>, A. Leroyer<sup>b,*</sup>, M. Visonneau<sup>b</sup>, P. Queutey<sup>b</sup>

<sup>a</sup> SIREHNA/Naval Group, Technocampus Ocean, 5 rue de l'Halbrane, 44340 Bouguenais, France

<sup>b</sup> Laboratoire de recherche en Hydrodynamique, Energétique et Environnement Atmosphérique, École Centrale de Nantes, CNRS-UMR 6598, 44321 Nantes Cedex 3, France

This work is devoted to the implementation and analysis of an efficient fluid-structure in-teraction coupling algorithm used in a
cooperative simulation context for rigid bodies. This framework makes possible the study of complex fluid-structure configurations,
through the coupling between CFD and dedicated multibody dynamics solvers.

A specific focus is here laid on the characterisation of the numerical properties of the proposed coupling algorithm in terms of convergence speed and stability. In order to stabilise the segregated algorithm while keeping non-intrusive in the structure solver, a relaxation operator based on an artificial added mass technique is used. To compute an efficient relaxation operator, an original resolution of the added mass effect is implemented into the fluid finite-volume solver. Good convergence properties are observed for rigid bodies with six degrees of freedom even in case of strong destabilising added mass effects.

**1. Introduction**

Co-simulation (Co-operative simulation) is a simulation methodology which uses several numerical solvers to work together in order to solve a global problem (Jürgens, 2009). It becomes possible and interesting because current numerical methods are mature enough to solve complex models in specific science fields (structural dynamics, fluid dynamics, thermodynamics, combustion, electromagnetism, etc.). Thus, the use of existing and assessed solvers looks an appealing and pragmatic way to handle multi-physics configurations.

For example, in naval architecture, Computational Fluid Dynamics (CFD) solvers are used more and more to check and/or
improve ships performances in resistance, propulsion, manoeuvrability and seakeeping (Queutey et al., 2012). Because CFD
solvers are dedicated and optimised to the resolution of the fluid problem, they are not developed to handle complex
mechanical systems like driven appendages, flexible floating wind turbines, wave energy converter, moored ships or towed
underwater systems for instance. When the mechanical systems are complex, Computational Structural Dynamics (CSD)
solvers can be used, but like CFD solvers, they are not dedicated to solve the whole problem. Indeed, in CSD solvers, the fluid
loading is often replaced by analytical models which are most of the time too simple to accurately address the complexity of
the flow physics (Masarati, 1999). To fully and accurately solve the whole problem when both fields are complex and highly
coupled, co-simulation is attractive due to the re-use of existing software including both advanced physical models and a
very strong numerical adaptation for each elementary problem (Kassiotis, 2009).

<footnote>* Corresponding author.

E-mail addresses: camille.yvin@sirehna.fr (C. Yvin), alban.leroyer@ec-nantes.fr (A. Leroyer), michel.visonneau@ec-nantes.fr (M. Visonneau), patrick.queutey@ec-nantes.fr (P. Queutey).</footnote>
Co-simulation is not only a method used by scientists and engineers to solve a general problem but an efficient work method too. Indeed, the use of different solvers allows different development teams to collaborate. From a practical point of view, the validation and maintenance of the solvers are facilitated (Hou et al., 2012) because the problem is clearly divided. However, the communication between each solvers must be performed at specific moments to be efficient. Thus, particular attention should be paid to coupling strategies to be both accurate and robust.

In Fluid-Structure Interaction (FSI), the problem is commonly divided to three problems (the fluid, the structure and the interface) but this is not always the case. Indeed, an FSI problem can be seen as a unique global problem and the resolution can be performed with a monolithic approach (Hübner et al., 2004; Heil, 2004; Dettmer and Perić, 2007; Saksono et al., 2007; Papadakis, 2008; Wick, 2011). This approach is difficult to implement in a co-simulation context because major modifications in each solver have to be implemented. Moreover the resolution is not an easy task because there are a lot of different unknowns of different nature, which, being partially coupled, have to be solved at the same time (Cervera et al., 1996). For the classical FSI problem that we are interested in, the parametrisation of the problem should contain the pressure, the velocity, the turbulence as well as the free surface unknowns, the mesh and the structural kinematics and the structural internal forces.

The partitioned approach (or block-iterative method Cervera et al., 1996) is most often used because of its ease of implementation. It consists in solving the fluid and structure problems in a segregated way. This approach is used here because it can be easily applied in a co-simulation context. However, the convergence may be slow (Degroote et al., 2009) and the coupling algorithm stability is not guaranteed due to the added mass effect (Belanger et al., 1995; Wall et al., 2006; Söding, 2001) (cf. Section 5). In order to stabilise the coupling algorithm, the artificial added mass method is sometimes used (Leroyer, 2004). This technique, which consists in increasing artificially the inertia of the structure and adding an acceleration dependent force, must be modified in a co-simulation context to avoid being too intrusive (cf. Section 5). In this paper, it is shown that the use of a physical relaxation operator based on the knowledge of the added-mass leads to an optimal convergence of the coupling in addition to be fully integrated, since the number of non-linear iterations to solve the fluid flow with FSI is similar to the one while imposing the motion (cf. Section 6). At the same time, the present work proposes an original approach to compute the classical added-mass operator when the frequency tends to infinite which overrides the limitations of classical linear potential solvers.

In this work, two different solvers are used. The first one solves the fluid problem and the other solves the mechanical behaviour. They are ISIS-CFD and MBDyn. Because of the abilities of both solvers, complex mechanical systems which strongly interact with fluids in naval architecture can be accurately studied.

Sections 2 and 3 briefly describe the two solvers. Section 4 introduces the FSI problem and the coupling algorithm used. Section 5 describes how the added mass effect is taken into account in the coupling algorithm. Section 6 explicits the numerical evaluation of the added mass effect in the framework of a finite-volume approach. Finally, two different applications are proposed in Section 7.

## 2. ISIS-CFD solver

ISIS-CFD is available as a part of the FINE™/Marine computing suite which is dedicated to marine applications. This is an incompressible unsteady Reynolds-averaged Navier-Stokes (RANS) solver developed by ECN-CNRS (Queutey et al., 2012). This solver is based on a cell-centred unstructured finite-volume method. Pressure-velocity coupling is obtained through a Rhie & Chow SIMPLE-type method. Free-surface flow is addressed with an interface capturing method, by solving a convection equation for the volume fraction of water, which is discretised with specific compressive discretisation schemes (Queutey and Visonneau, 2007). An Arbitrary Lagrangian Eulerian (ALE) formulation is used to take into account modification of the fluid spatial domain (Leroyer et al., 2008). It is associated with robust and fast grid deformation techniques (Leroyer and Visonneau, 2005). The temporal discretisation scheme is the Backward Difference Formula of order 2 (BDF2) when dealing with unsteady configurations. For each time step, an inner loop (denoted by non-linear loop) associated to a Picard linearisation is used to solve the non-linearities of the system.

## 3. MBDyn solver

MBDyn (Multi-Body Dynamics), is an open-source solver under the GNU GPL license developed at the Dipartimento di Ingegneria Aerospaziale of the Politecnico di Milano. It is intended for the simultaneous solution of multi-discipline problems including non-linear dynamics, aero-servo-elasticity, smart piezo-structural components, electric and hydraulic components. It is aimed at the modelling of complex systems (Ghiringhelli et al., 1999).

To solve the kinematic laws of a multi-body mechanical system, the Redundant Coordinate Set (RCS) formulation is used. This means that every inertial body has six rigid body Degrees of Freedom (DOF) even if they are constrained by joints for instance. Additional holonomic or nonholonomic constraint equations are added which introduce algebraic unknowns that are analogous to the Lagrange multipliers and directly represent the reaction forces and couples (Masarati, 1999). All these equations are written in the form of a set of first order Algebraic Differential Equations (ADE). Thus every rigid body is represented by 12 equations (6 equations that represent the momenta and 6 others the Newton's law that links the rate of changes in momentum due to forces). In this work, the BDF2 scheme is used. The ability to take into account multidisciplinary complex systems and the simplicity of implementation are the main advantages of this formulation.
## 4. FSI problems with a partitioned approach

FSI problems can be seen as a particular application of decomposition domain methods which consists in solving a boundary value problem by splitting it into smaller boundary value problems on subdomains. These methods are suitable for parallel computing but also for FSI problems because interactions occur only at the interface between the fluids and the solids. This is the reason why it is relevant to express this problem with the Steklov-Poincaré operator and its inverse (Deparis et al., 2006). These non-linear and time dependent operators can also be seen as the fluid or the solid solver.

Thus, the physical domain $\Omega$ is split into a fluid $\Omega^f$ and solid domains $\Omega^s$. The Steklov-Poincaré operator $S$, and its inverse $S^{-1}$, are defined for each domain as:

$$S_d(\mathbf{x}^d) = \lambda^d \quad S_d^{-1}(\lambda^d) = \mathbf{x}^d \qquad (1)$$

where the superscript $d$ represents the considered domain, $\mathbf{x}$ is the vector position and $\lambda$ the stress which acts at the interface between the fluid and the structure. At the fluid-structure interface $\Gamma = \bar{\Omega}^f \cap \bar{\Omega}^s$, the interface continuity and the action-reaction principle have to be satisfied, namely:

$$\mathbf{x}^s = \mathbf{x}^f = \mathbf{x}^\Gamma \quad \text{and} \quad \lambda^s + \lambda^f = \mathbf{0} \quad \text{on } \Gamma \qquad (2)$$

These two equations can be expressed with the operators previously defined in different formulations like the Steklov-Poincaré, the Picard (fixed-point) or the Newton (root-finding) formulations (Deparis et al., 2006). The first one is suitable when the time simulation to solve the fluid and structure problem are of the same order because the parallelisation is naturally achieved as $S_s$ and $S_f$ or their inverse are not composed. This is not the case here because the fluid problem needs much more computation time than the solid problem. This is due to the high number of unknowns in the fluid discretisation. To use a Newton strategy, the Jacobian of the problem must be computed or approximated. It can be more expensive than solving a simple fixed-point equation because of the non-linearity of the coupling terms (Dettmer and Perić, 2008). A simple implicit Block Gauss-Seidel (BGS) approach is classically used to solve the following fixed point formulation which sums up the coupled problem:

$$\lambda^f = -S_f \circ S_s^{-1}(\lambda^f) \qquad (3)$$

Eq. (3) is a non-linear time dependent equation and normally has to be fully solved at each time step. Explicit coupling schemes in time, which perform only one iteration to solve the Eq. (3), can produce or dissipate energy at the interface if the problem is highly coupled. Moreover, they are very sensitive to the destabilising added mass effect (Badia and Codina, 2000). In this work, only implicit time integration schemes are considered. This means that the fluid, the structural and the interface problems are solved implicitly. The coupling and time iterations are denoted $i$ and $n$, respectively. A simple implicit BGS algorithm, like the one presented in Fig. 1a or in Eqs. (4) and (5), is easy to implement but needs a lot of coupling iterations before reaching convergence and can be rapidly unstable when the added-mass effect becomes important (Förster et al., 2007).

$$\mathbf{x}^\Gamma|_{n+1}^{i+1} = S_s^{-1}(\lambda^f|_{n+1}^i) \qquad (4)$$

$$\lambda^f|_{n+1}^i = -S_f(\mathbf{x}^\Gamma|_{n+1}^i) \qquad (5)$$

In order to reduce the simulation time, the Aitken $\Delta^2$ relaxation technique is often used (Küttler and Wall, 2008). This technique is based on a geometrical approximation (tangent method) of the coupled problem. To be efficient, the solvers have to produce physical results, and consequently, each operator has to be solved accurately. This means that for an implicit coupling algorithm, the fluid problem has to be solved completely several times at each time step, putting a strain on the computational time.

In this work, the classical Steklov-Poincaré operator for the fluid problem is modified to reduce the CPU time of this part. This new operator, denoted as $S_f^*$, does not represent a global fluid resolution any more but only one Picard iteration of the fluid solver (Queutey and Visonneau, 2007) as has already been proposed (Leroyer, 2004; Badia and Codina, 2000). It can be seen as a fluid linearised version of the Steklov-Poincaré operator. In a few words, at each time step, the iterative resolution of the fluid problem is embedded with the iterative resolution of the coupled problem. On the contrary, because the CPU time to solve the structural problem is weak compared to the other parts, the structure Poincaré-Steklov operator is fully resolved at each coupling iteration. This algorithm is denoted as a Block Gauss-Seidel with Internal Fluid Coupling (BGS-IFC) approach.

## 5. Added-mass effect

The added-mass effect often occurs in naval architecture because the density of the water is not negligible compare to the density of the structure. As previously mentioned, one of the main techniques to eliminate this destabilising added-mass effect is the artificial added-mass method. To illustrate this point in a co-simulation context, a one-dimensional case is presented. A rigid body of mass $m$ has one degree of freedom in translation $x$. The fluid force and the structure force (gravity,
joints, etc.) are denoted as $f_s$ and $f_s$, respectively, and are computed by the fluid and the structure solvers. The added mass effect can be seen as the part of the pressure fluid force which is proportional and opposed to the acceleration of the body. The fluid force is split accordingly and the added mass part is denoted by $-m_a \ddot{x}$ while the other part is denoted by $\bar{f}_f$. This splitting is just for the presentation because in reality a classical CFD solver cannot separate the added mass effect from the total fluid force as it is commonly done in potential flow. The equation which has to be solved is:

$$m \cdot \ddot{x} = f_s + \bar{f}_f = f_s + \bar{f}_f - m_a \cdot \ddot{x} \quad (6)$$

Even if an implicit integration scheme is used, the resolution of the previous equation is unstable when the added mass effect is close or higher than the mass (Söding, 2001). The original artificial added mass technique affects the inertia of the body by an artificial added mass coefficient $\overline{m}_a$ and thus decreases the dependency of the right hand side member with the acceleration. The corresponding equation which has to be solved by the CSD solver is:

$$(m + \overline{m}_a) \cdot \ddot{x}|_{n+1}^{i+1} = f_s|_{n+1}^i + \bar{f}_f|_{n+1}^i - m_a \cdot \ddot{x}|_{n+1}^i + \overline{m}_a \cdot \ddot{x}|_{n+1}^i \quad (7)$$

One can notice that when convergence is reached, which means $\ddot{x}|_{n+1}^{i+1} = \ddot{x}|_{n+1}^i$, Eq. (6) is recovered. A theoretical lower limit for the artificial added mass coefficient to guarantee the stability of the coupling algorithm can be found when $\overline{m}_a > (m_a/m - 1)/2$. This is valid when an implicit Euler scheme of order 1 is used (Söding, 2001). However, it can be easily seen that in order to fully eliminate the acceleration dependence of the right member of Eq. (7), the artificial added mass coefficient $\overline{m}_a$ should be equal to $m_a$. Discussion about its optimum value can be found in the next section. The last term of the right side acts as a force and must be given by the coupling algorithm. On the other side, the structure solver has to be modified due to the modification of the left member. In a co-simulation context, this is not a good solution because modification of the structural solver can be a difficult task, especially for bodies with six DOF. To avoid this, a relaxation factor $\alpha$ is introduced (Eq. (8)). Eq. (7) is multiplied by $\alpha$. The corresponding equation obtained is given in Eq. (9).

$$\alpha = \frac{1}{1 + \overline{m}_a/m} \quad (8)$$

$$m \cdot \ddot{x}|_{n+1}^{i+1} = \alpha (f_f|_{n+1}^i + f_s|_{n+1}^i) + \alpha \cdot \overline{m}_a \cdot \ddot{x}|_{n+1}^i \quad (9)$$

This last equation can be seen as a relaxation of the acceleration because the terms in brackets on the right side represent the acceleration computed by the structural solver multiplied by the mass. Thus, if $\ddot{x}$ represents this acceleration, Eq. (9) can also be written with the following equations:

$$m \cdot \ddot{\tilde{x}}|_{n+1}^{i+1} = f_f|_{n+1}^i + f_s|_{n+1}^i \quad (10)$$

$$\ddot{\tilde{x}}|_{n+1}^{i+1} = \tilde{x}|_{n+1}^i + \alpha \cdot (\tilde{\delta}|_{n+1}^{i+1} - \ddot{\delta}|_{n+1}^i) \quad (11)$$

There are two differences between the original added artificial added mass method (cf. Eq. (7)) and this new set of equations. Indeed, in the first one, the structural force $f_s$ is computed at the same time as the new acceleration. This is not the case for the proposed solution because the relaxation occurs after a CSD resolution. In the same way, the non-linearities due to the rotation motion are not solved in the same manner with this two different approaches. Despite of this, good convergence properties were observed if the artificial added mass coefficient is well chosen.

The final fluid-structure algorithm used is summarised by Eq. (12) and presented in Fig. 1b where **R** is the general relaxation operator previously defined for a one DOF problem (Eq. (8)) and $\ddot{\delta}^s$ is the interface position before the relaxation step. It is named Block Gauss-Seidel with Internal Fluid Coupling (BGC-IFC) algorithm with relaxation. It can be noticed that at the first iteration of a time step, the kinematics of the fluid-structure interface is given by a prediction of order 2 with constant acceleration.

$$\lambda^f|_{n+1}^{i+1} = -S_f^* \circ \mathbf{R} \circ S_s^{-1} (\lambda^f|_{n+1}^{i+1}) \quad (12)$$

In case of a six DOF rigid body, the interface displacement can be replaced by the generalised position (position and orientation) of the body denoted as $\delta$ and the added mass effect can be represented by a symmetric matrix of rank six denoted as **MA** (cf. Eq. (17)). The corresponding artificial added mass matrix is denoted as $\overline{MA}$. Thus, the acceleration relaxation is given by:

$$\ddot{\delta}|_{n+1}^{i+1} = \ddot{\delta}|_{n+1}^i + \mathbf{R}_a \cdot (\ddot{\delta}|_{n+1}^{i+1} - \ddot{\delta}|_{n+1}^i) \quad (13)$$

where $\mathbf{R}_a$ is the acceleration relaxation operator. The application of Eq. (7) to a linearised six DOF problem gives:

$$\mathbf{R}_a = (\mathbf{I}_d + \mathbf{M}^{-1} \cdot \overline{\mathbf{MA}})^{-1} \quad (14)$$

where $\mathbf{I}_d$ is the identity matrix. After the acceleration relaxation step, velocities and generalised positions must be reconstructed according to the time integration scheme used. The reconstruction of the orientation depends also on the parametrised rotation used. Moreover, due to the non linearities of rotation motions, the general relaxation operator **R**, which consists of a relaxation of the acceleration and a reconstruction of velocities and generalised positions, is non-linear.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Two implicit coupling algorithms (x vector position, λ stress at the fluid-structure interface, $S_s^{-1}$ structure operator, $S_f$ fluid operator, $S_i^*$ linearised fluid operator).</figcaption></figure>

# 6. Added mass computation

## 6.1. Introduction

As it was shown, the artificial added mass technique can be considered as a particular relaxation technique suitable to the FSI problem. As every relaxation method, it is important to find a correct value of the relaxation parameter in order not to deteriorate the convergence speed. The first numerical experience shows that a good convergence behaviour is obtained if the artificial added mass is close to the “physical” added mass.

Added mass effect is sometimes hard to tackle because reference values are not always available and furthermore it can change during the simulated time of the computation (presence of walls or other bodies, confinement, large free surface modification, etc.). This is the reason why an added mass pressure equation has been implemented in ISIS-CFD to compute dynamically the added mass matrix. This computation can be performed once and for all before the simulation or at different times if necessary.

If the notion of added-mass is unequivocal when dealing with a body in a monofluid inviscid flow, the situation differs when dealing with viscous flow. When periodic motion is encountered, added mass effect is sometimes defined as the fluid force which is in phase with the accelerations (Vikestad et al., 2000; Uzunoğlu et al., 2001; Brennen, 2008). It is generally used for specific applications where the rigid body has very few degrees of freedom. In this case, all effects (viscous effects, history effects, etc.) are taken into account (Brennen, 1982).

The added mass effects can also be defined as the fluid response due to a brutal variation of the rigid body velocity (impulse acceleration) (Leonard and Roshko, 2001; Söding, 2001; Guyon et al., 2001; Wakaba and Balachandar, 2007). This is the definition which is retained here.

To derive an equation, an analysis of the time scales of the terms of the Navier-Stokes equations is proposed. The Navier-Stokes equations, under the convective form for an incompressible flow, can be written as:

$$ \frac{\partial u_i}{\partial t} + u_j \frac{\partial u_i}{\partial x_j} = \frac{1}{\rho} \left[ -\frac{\partial p}{\partial x_i} + \frac{\partial}{\partial x_j} \left( \mu \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) \right) \right] + f_i^v \quad (15) $$

where u is the fluid velocity, p the pressure, f^v the volumetric forces (gravity), ρ the fluid density and μ the dynamic viscosity. Let T, L and U the characteristic time, length and velocity of the fluid. The time scales of the different terms are in order: T (temporal term), L/U (convective term), zero (pressure term), L^2ρ/μ (viscous diffusion terms) and infinity (volumetric
production which is equal to gravity) (Guyon et al., 2001). Thus, when the characteristic time scale tends to zero (impulsion), the fluid response is dominated by the pressure term, namely:

$$ \frac{\partial \tilde{u}_i}{\partial t} = - \frac{1}{\rho} \frac{\partial \tilde{p}}{\partial x_i} \qquad (16) $$

The corresponding velocity and the pressure field are denoted as $\tilde{u}$ and $\tilde{p}$ to avoid any confusion with the real fields. The equations driving the added-mass effect for a monofluid potential flow are recovered. These equations can also be derived by considering the difference between two flows which differ only at the level of the acceleration of the body (Söding, 2001). This approach has already been analysed and confirmed on a simple example (Wakaba and Balachandar, 2007).

Hence, for each degree of freedom of the body, a unit acceleration can be applied and transmitted to the fluid domain through the kinematic condition induced at the surface of the body. The resultant pressure field integrated over the body leads to the corresponding column of the added mass matrix (cf. Section 6.2) which is defined as:

$$ f_{ma} = -MA \cdot \ddot{\delta} \qquad (17) $$

In that case, MA depends only on the geometry of the body and is symmetric (Landau and Lifshitz, 1987). When the fluid domain is not infinite, this matrix depends also on the spatial configuration (confinement, other bodies, etc.). When a free-surface flow is considered, both phases are modelled by a single fluid which properties (density and viscosity) change with respect to the volume fraction. Thus, the free-surface is naturally taken into account when computing the added-mass effect through the spatial evolution of the density $\rho$ which is part of Eq. (16).

This is not the case when a classical potential flow solver is used because the kinematic and the dynamic boundary conditions at the free surface have also to be satisfied Söding (2001). In addition to that, for unsteady cases where the free surface position dramatically change with respect to the considered body, (for example in the cases of drop test where the body is at the initial stage fully surrounded by air but not after the water entry or in configurations where the free-surface has large distortion in time with possibly breaking waves, ...), the proposed approach lets the possibility to update the added-mass effects at given time steps.

By taking the divergence of Eq. (16) and because the fluid velocity is divergence free, the temporal term also vanishes (Eq. (18)). The corresponding integral form of this equation is given at Eq. (19). The Green-Ostrogradski theorem is used to get a surface equation which is integrated over all the volume controls (Eq. (20)) which are denoted as V.

$$ -\frac{\partial}{\partial x_i} \left( \frac{1}{\rho} \frac{\partial \tilde{p}}{\partial x_i} \right) = 0 \qquad (18) $$

$$ \iiint_V \nabla \cdot \left( -\frac{\nabla \tilde{p}}{\rho} \right) \cdot dv = 0 \qquad (19) $$

$$ \iiint_S -\frac{1}{\rho} \nabla \tilde{p} \cdot \mathbf{n} \cdot d\mathbf{s} = 0 \qquad (20) $$

## 6.2. Boundary conditions and matrix evaluation

Boundary conditions are identical to the one used for the classical pressure, except at the considered body where unit acceleration are used for each DOF (Söding, 2001). Because the viscous effects are not taken into account, continuity of normal velocity conditions taken from the Eq. (16) are applied. So, if the cell face is a part of the considered body, the following Neumann boundary condition is used:

$$ \frac{\partial \tilde{u}}{\partial t} \cdot \mathbf{n} = -\frac{1}{\rho} \nabla \tilde{p} \cdot \mathbf{n} = \left[ \frac{d\tilde{v}}{dt} + \frac{d\tilde{w}}{dt} \wedge (\mathbf{x} - \mathbf{x}_G) \right] \cdot \mathbf{n} \qquad (21) $$

where $\mathbf{x}_G$ is the gravity centre of the rigid body and $\tilde{v}$ and $\tilde{w}$ are its linear and angular velocities:

$$ \ddot{\delta} = \frac{d}{dt} \left[ \frac{\tilde{v}}{\tilde{w}} \right] \qquad (22) $$

For each DOF j, a pressure field $\tilde{p}_{j,i}$ is evaluated through the numerical resolution of Eq. (20) and the application of a unit acceleration for these DOF (acceleration of the other DOF are set to zero). This field is integrated over the considered body to get the three forces and the three moments at the gravity centre. The added mass matrix is defined as:

$$ -MA_{ij} = \begin{cases} - \left( \iiint_S \tilde{p}_{j,i} \mathbf{n} \, \mathrm{d}\mathbf{s} \right) \cdot \mathbf{e}_i & \text{for } i = 1, 2, 3 \\ - \left( \iiint_S \tilde{p}_{j,i} (\mathbf{x} - \mathbf{x}_G) \wedge \mathbf{n} \, \mathrm{d}\mathbf{s} \right) \cdot \mathbf{e}_{i-3} & \text{for } i = 4, 5, 6 \end{cases} $$
<figure><img src="image_3.png" /><figcaption>Fig. 2. Added mass in translation of an infinite rigid cylinder radius $r_1$ in a pipe radius $r_2$ (numerical results with the finer mesh).</figcaption></figure>

### 6.3. Discretisation

Discretisation of Eq. (20) is performed with the general methodology presented in Queutey and Visonneau (2007) which is very briefly sump up here. The finite-volume method gives:

$$ \sum_f -\left(\frac{\nabla \tilde{p}}{\rho}\right)_f \cdot \mathbf{n}_f \cdot S_f = 0 \qquad (23) $$

Pressure and volumetric mass are evaluated at the centre of the cells so the reconstruction from the centre of the cell to the centre of the face of the pressure gradient divided by the volumetric mass is needed. More precisely, it is the normal projection to the face that has to be reconstructed. In a case of a two phases flow, the pressure gradient and the volumetric mass are not continuous so special attentions have to be made. It can be shown that a general reconstruction scheme which can handle the discontinuity due to the presence of a free surface can be written as (Queutey and Visonneau, 2007):

$$ \left(\frac{\nabla \tilde{p}}{\rho} \cdot \mathbf{n}_f\right)_f = \frac{1}{\hat{\rho}} \left(\frac{\tilde{p}_R}{h} - \frac{\tilde{p}_L}{h} + \frac{E_{\nabla \tilde{p}}}{h}\right) \qquad (24) $$

where the explicit part is:

$$ E_{\nabla \tilde{p}} = (\rho_R h^+ \mathbf{e}^+ + \rho_L h^- \mathbf{e}^-) \cdot \left[ \frac{h^+}{h} \left(\frac{\nabla \tilde{p}}{\rho}\right)_L + \frac{h^-}{h} \left(\frac{\nabla \tilde{p}}{\rho}\right)_R \right] \qquad (25) $$

L and R are the centres of the two finite volumes which share the face f. $h$, $h^+$ and $h^-$ are geometrical coefficients and $\hat{\rho}$ is a coefficient which depends on the previous geometrical coefficients but also $\rho_R$ and $\rho_L$. The gradients at the cells centres of the explicit part are evaluated with a Green-Gauss theorem (Queutey and Visonneau, 2007).

### 6.4. Validation

First, the added mass in translation of an infinite rigid cylinder in a pipe is compared to the analytical solution which is given by:

$$ \mathbf{MA}_{11} = \left[ \frac{1 + \left(\frac{r_1}{r_2}\right)^2}{1 - \left(\frac{r_1}{r_2}\right)^2} \right] \cdot \rho \pi r_1^2 \qquad (26) $$

where $r_1$ is the cylinder radius, $r_2$ is the pipe radius and $\rho$ the volumetric mass. The simulation results are given in Fig. 2. The maximum error is converging with the mesh density: 1.27% for a mesh with 900 cells, 0.4% with 3600 cells and 0.096% with 14 400 cells.

Second, the case of a full scale ship Cargo Series 60 is presented. The ship length is about 125 m. Results are compared to the ones obtained with the frequency-domain seakeeping code Aqua+ (Delhommeau, 1987). It is important to note that the temporal approach is identical to the frequency approach at the limit of null period.

Potential flow solvers need specific surface meshes of the ship. Two different meshes are tested. They are made up of 1400 and 5600 cells, respectively, for the half ship. The coarser is presented in Fig. 3. For the proposed method, two typical finite-volume hex-dominant meshes are used. The coarser and the finer meshes are made up of 2.9 and 3.6 millions cells, respectively, for the entire ship. The boat is discretised with 53 000 and 161 000 of cells, respectively. The coarser is presented in Fig. 4.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Cargo Series 60 — Coarser Aqua+ mesh.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Cargo Series 60 — Coarser ISIS-CFD mesh.</figcaption></figure>

The six pressure fields obtained for finest meshes are presented in Fig. 5 for both solvers. The added mass matrix obtained with the proposed method with the finest mesh is given in Eq. (27) with S.I. units. Null values are set due to the Y-Z symmetric plane. Numerically, they are only negligible compared to the others due to discretisation and iterative errors, but they are not significant. Some differences (in percentage) are given from Eqs. (28) to (30). Only non-null terms are analysed. The reference values come from the case which is in subscript. In all cases, significant differences are only observed for some extra-diagonal terms. It is important to notice that these terms are more difficult to obtain because they are due to the small asymmetries of the boat. Nevertheless, with the proposed method, the two different meshes produce similar values (cf. Eq. (29)). Moreover, the differences between the potential solver and the proposed method (cf. Eq. (30)) are significant only for the two extra-diagonal terms which are not fully converged with the potential solver (cf. Eq. (28)). Such differences are not visible when comparing the pressure fields (cf. Fig. 5).

In this case, the loop required to solved completely the discrete problem (cf. Eq. (24)) due to explicit terms is even not mandatory to obtain reasonable estimations, since diagonal and extra-diagonal terms at the first iteration are different from the converged ones by only 2% and 3%, respectively (cf. Fig. 6). In order to completely solve the explicit terms, 4 iterations are needed.

$$ \mathbf{MA} = \begin{bmatrix} 1.65e5 & 0 & -1.72e5 & 0 & 2.04e7 & 0 \\ 0 & 5.48e6 & 0 & 2.95e6 & 0 & -4.45e6 \\ -1.72e5 & 0 & 8.85e6 & 0 & -2.61e7 & 0 \\ 0 & 2.96e6 & 0 & 3.23e7 & 0 & -1.06e8 \\ 2.05e7 & 0 & -2.61e7 & 0 & 4.33e9 & 0 \\ 0 & -4.45e6 & 0 & -1.06e8 & 0 & 6.02e9 \end{bmatrix} \qquad (27) $$

$$ \mathbf{E}_{\text{Aqua+ coarse [\%]}} = \begin{bmatrix} -0.7 & -0.3 & -0.6 & -10.5 \\ -0.8 & -3.9 & -0.0 & -3.5 \\ -0.8 & -1.2 & -1.7 & -1.9 \\ -0.1 & -0.8 & -0.8 & -0.3 \\ -0.1 & -0.8 & -0.8 & -0.3 \end{bmatrix} \qquad (28) $$

$$ \mathbf{E}_{\text{ISIS-CFD coarse [\%]}} = \begin{bmatrix} -0.2 & 0.2 & 0.1 & 1.4 \\ 0.1 & 1.6 & 0.2 & -0.3 \\ 0.1 & -0.2 & 0.8 & -1.4 \\ 0.2 & -0.1 & 0.1 & -0.3 \\ 0.3 & -0.1 & -0.3 & 0.2 \end{bmatrix} \qquad (29) $$

$$ \mathbf{E}_{\text{ISIS-CFD fine [\%]}} = \begin{bmatrix} 4.2 & 3.0 & 2.4 & -14.8 \\ 2.5 & 9.3 & -1.9 & -9.8 \\ -2.0 & -1.9 & -1.5 & -3.3 \\ -0.4 & -1.5 & -2.0 & -14.7 \end{bmatrix} \qquad (30) $$
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Cargo Series 60 – Pressure fields in Pa for each degrees of freedom – View from the bottom – Aqua+ (left) and ISIS-CFD (right).</figcaption></figure>

### 6.5. Remarks

The proposed method to evaluate the added mass matrix is interesting for several reasons. Firstly, no additional work has to be done because the fluid mesh can be the same as the one used for the fluid-structure interaction problem. When the flow around moving bodies is solved, the added-mass effects are indeed part of the pressure field. As a consequence, the fluid mesh needs to be fine enough to capture this physical feature. However, when the goal is only to compute the pressure field due to added-mass, the mesh used for the full problem contains fine cells close to the wall to capture the boundary layer which are not required for the computation of this impulsive inviscid problem. These cells could be removed, but even if they are present, they do not affect the solution. Secondly, only one linear iteration for each degree of freedom is enough to correctly estimate the added mass matrix of a rigid body which means that the relaxation operator can be dynamically evaluated during the simulation without penalising the simulation time: even if a new evaluation of the full 6 × 6 added-mass matrix is requested every 20 times step, the additional CPU time is no more than 1%. Finally, this method does not suffer from any limitations of linear potential solvers. For example, contrary to the latter, fine meshes can be used so all the details of the geometry can be entirely respected. Moreover, a large amount of space configurations, like complex free surface position (wave breaking) or important confinement (shallow water, interactions with close bodies), are naturally taken into account.
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Cargo Series 60 — Added mass as a function of iterations to converge the explicit part of the discretisation.</figcaption></figure>

<table><caption>**Table 1**
Applications — Inertial characteristics.</caption><thead><tr><th rowspan="2">Number</th><th colspan="1">Case 1</th><th colspan="4">Cases 2</th><th rowspan="2">Unit</th></tr><tr><th>1</th><th>2.1</th><th>2.2</th><th>2.3</th><th>2.5</th><th>2.7</th></tr></thead><tbody><tr><td><i>d</i></td><td>0.1</td><td>0.1</td><td>0.2</td><td>0.3</td><td>0.5</td><td>0.7</td><td>-</td></tr><tr><td><i>m</i></td><td>400</td><td>400</td><td>800</td><td>1200</td><td>2000</td><td>2800</td><td>kg</td></tr><tr><td><i>I</i><sub>xx</sub>/<i>B</i><sub>3</sub></td><td>140</td><td>140</td><td>285</td><td>430</td><td>710</td><td>990</td><td>kg m<sup>2</sup></td></tr><tr><td><i>I</i><sub>yy</sub>/<i>B</i><sub>3</sub></td><td>540</td><td>880</td><td>1800</td><td>2600</td><td>4400</td><td>6200</td><td>kg m<sup>2</sup></td></tr><tr><td><i>I</i><sub>zz</sub>/<i>B</i><sub>3</sub></td><td>670</td><td>1000</td><td>2000</td><td>3000</td><td>5000</td><td>7000</td><td>kg m<sup>2</sup></td></tr></tbody></table>

# 7. Applications

## 7.1. Generalities

The present applications are devoted to the study of light parallelepipeds with six DOF. The fluid considered is water. Its density and dynamic viscosity are equal to $\rho = 1000 \text{ kg m}^{-3}$ and $\mu = 0.001 \text{ kg m}^{-1} \text{ s}^{-1}$, respectively. The densities of the parallelepipeds are smaller than one, thus, they tend to go up due to the buoyant force. Dimensions of the parallelepipeds are $L_x = 4 \text{ m}$, $L_y = 2.0 \text{ m}$ and $L_z = 0.5 \text{ m}$. At time $t = 0$, the roll, pitch and yaw angles of the parallelepipeds are equal to 20°, 15° and 35°, respectively. The coordinate systems attached to the body and expressed at the centres of gravity G are denoted as $B_3$. The geometric dimensions and inertial properties of the parallelepipeds are chosen in order to have important anisotropic added mass effects. It is important to notice that if the added mass is not taken into account or, if the artificial added mass matrix is set equal to the generalised inertia matrix as it is often done, the coupling algorithm does not converge.

For the first application, the gravity centre G matches with the geometric centre C and one density is considered (0.1). In this case, the acceleration relaxation operator $R_a/B_3$ is diagonal because both the added mass and the inertia matrix are diagonal at the gravity centre. For the second application, several densities are considered (0.1, 0.2, 0.3, 0.5 and 0.7) and the geometric centre is shifted from the geometric centre C by 1 m in first direction of $B_3$. Consequently $R_a/B_3$ are no longer diagonal because some degrees of freedom are coupled regarding to the added mass effect. The inertial characteristics used are given in Table 1. Inertias in rotation are expressed at the gravity centre of the considered parallelepiped.

The geometry, the initial orientation and the important points are presented at the Fig. 7. The fluid domain is a 40 m cube. The mesh is always the same. It is made up of 1.0 million cells leading to 8200 boundary faces around the parallelepiped. The parallelepipeds kinematics are taken into account by a rigid motion of the mesh. The pressure field is fixed at one face of the domain. Otherwise, the velocity is fixed to zero. No-slip conditions are applied on the faces of the parallelepipeds and no turbulence model is used.

The components of the position of the gravity centres are denoted as $(x, y, z)$. The Cardan angles $(\phi, \theta, \psi)$ are used to represent the orientation of the parallelepipeds. The components of the linear and angular (instantaneous rotation vector) velocities are denoted $(u, v, w)$ and $(p, q, r)$, respectively. They are sometimes expressed in $B_3$ in order to facilitate the analysis by degrees of freedom. The angular Cardan velocities are denoted $(\phi, \dot{\theta}, \psi)$. The components of the fluid forces and moments are given as $(f_x, f_y, f_z)$ and $(m_x, m_y, m_z)$, respectively. They are always expressed at the gravity centre of the considered parallelepiped.
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Applications — Initial conditions (C geometric centre, G gravity centre).</figcaption></figure>

In order to precisely analyse the convergence properties of the algorithm, 30 coupling iterations per time step are performed. It means that the results do not depend on the definition of a stopping criteria. However, a criteria of convergence is computed during the post-treatment steps for each degree of freedom. It is defined as the necessary number of iterations in order that the difference between the relaxed acceleration and the computed one at the end of the time step is smaller than a limit. In a nutshell, more this quantity is low more the coupling algorithm is efficient. This number is evaluated at each time step and averaged along the simulation time which is always about 15 s (as a consequence, the convergence criteria is a real number). The limit is set to 0.005 m s⁻² for degrees of freedom in translation and 0.005 rad s⁻² for the ones in rotation. For example, the criteria of convergence for the surge motion is evaluated as:

$$ \text{mean}_n \left( \min_{j \in [1,30]} (j, \forall i \ge j |\tilde{u}_n^i| - |\tilde{u}_n^{30}|) \le 0.005 \right) \quad (31) $$

## 7.2. Application 1

The added mass matrix computed from the proposed method is given in Eq. (32) and the corresponding acceleration relaxation operator (cf. Eq. (14)) is given in Eq. (33). One can see that the relaxation operator is diagonal and highly anisotropic. Moreover, some relaxation coefficients are especially low which means that the added mass effects can be very strong. Small relaxation coefficients usually leads to a dramatical increase of the number of iterations to reach the convergence state. As it will be shown, it is not the case since the relaxation has a physical origin in this procedure.

$$ \mathbf{MA}_{/B_3} = \begin{bmatrix} 564 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1255 & 0 & 0 & 0 & 0 \\ 0 & 0 & 11348 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1633 & 0 & 0 \\ 0 & 0 & 0 & 0 & 9657 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1274 \end{bmatrix} \quad (32) $$

$$ \mathbf{R}_{a/B_3} = \begin{bmatrix} 0.415 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0.242 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0.034 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0.080 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0.053 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0.343 \end{bmatrix} \quad (33) $$

Several time steps $\Delta t$ from 0.002 s to 0.050 s are considered in this study. The parallelepiped behaviour is presented in Fig. 8 when $\Delta t = 0.010$ s. The order of magnitude of the Reynolds number based on the averaged ascending velocity after 2 s is about $10^6$. An ascending dead leaf motion shape is always observed for each time step but the kinematics can slightly differ after 7 s as it is shown at Fig. 9.

The differences are due to the strong coupling between the motion of the parallelepiped and the fluid flow which is complex, very unstable and unsteady as it is shown in Fig. 10. As a consequence, this case with extreme added-mass effects is really sensitive to numerical errors, so that the accumulation of time discretisation errors can lead to significant differences during large simulation run-time. However, for each configuration, involving its proper time discretisation error, the convergence of the FSI coupling was fully achieved. What is analysed next is thus the convergence of the FSI discrete problem independently of the spatial and temporal discretisation error, but not the absolute convergence of the motion, which is really dependent on the chosen time step. The number of FSI iterations to reach the adequacy between the structure
<figure><img src="image_10.png" /><figcaption>Fig. 8. Parallelepiped 1 — Trajectory and orientation — $\Delta t = 0.025$ s.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Parallelepiped 1 — Position and orientation.</figcaption></figure>

solver and the fluid solver below a given tolerance is studied as a function of the solid density, the time step and the relaxation method used, as a criteria of the efficiency of the FSI algorithm.

Convergence information are presented in Table 2. One can notice that the number of iterations to reach the convergence depends on the considered degree of freedom. Indeed, the highest numbers of iterations in translation and rotation are obtained for the sway $\dot{v}$ and the roll $\ddot{p}$ motions respectively. In order to explain the differences between the degrees of freedom, the quadratic time averages of the added mass forces and the temporal derivatives of the accelerations (cf. Eq. (17)) are proposed at Tables 3 and 4. The highest added mass effects are not observed for the sway $\dot{v}$ and the roll $\ddot{p}$ motion but
<figure><img src="image_12.png" /><figcaption>Fig. 10. Parallelepiped 1 — Trajectory, orientation and iso-surface of dynamic pressure ($p = -1000$ Pa) coloured with fluid velocity- $\Delta t = 0.010$ s.</figcaption></figure>

<table><caption>**Table 2**
Parallelepiped 1 — Averaged number of iterations to reach the convergence criteria (Eq. (31)).</caption><thead><tr><th rowspan="2">Δt [s]</th><th colspan="3">Translation [-]</th><th colspan="3">Rotation [-]</th></tr><tr><th>ū<sub>/B3</sub></th><th>v̇<sub>/B3</sub></th><th>ẇ<sub>/B3</sub></th><th>ṗ<sub>/B3</sub></th><th>q̇<sub>/B3</sub></th><th>ṙ<sub>/B3</sub></th></tr></thead><tbody><tr><td>0.002</td><td>5.4</td><td><b>6.5</b></td><td>2.8</td><td><b>2.6</b></td><td>1.2</td><td>1.8</td></tr><tr><td>0.005</td><td>7.9</td><td><b>8.6</b></td><td>6.5</td><td><b>6.3</b></td><td>2.6</td><td>3.9</td></tr><tr><td>0.010</td><td>8.7</td><td><b>10.2</b></td><td>8.2</td><td><b>8.5</b></td><td>4.6</td><td>5.3</td></tr><tr><td>0.025</td><td>9.2</td><td><b>10.6</b></td><td>8.0</td><td><b>9.4</b></td><td>6.6</td><td>7.0</td></tr><tr><td>0.050</td><td><b>10.7</b></td><td><b>10.9</b></td><td>8.8</td><td><b>12.6</b></td><td>9.1</td><td>9.8</td></tr></tbody></table>

<table><caption>**Table 3**
Parallelepiped 1 — Quadratic time average of the added mass forces.</caption><thead><tr><th rowspan="2">Δt [s]</th><th colspan="3">Translation [N]</th><th colspan="3">Rotation [N m]</th></tr><tr><th>f<sup>ma</sup><sub>x</sub>/B<sub>3</sub></th><th>f<sup>ma</sup><sub>y</sub>/B<sub>3</sub></th><th>f<sup>ma</sup><sub>z</sub>/B<sub>3</sub></th><th>m<sup>ma</sup><sub>x</sub>/B<sub>3</sub></th><th>m<sup>ma</sup><sub>y</sub>/B<sub>3</sub></th><th>m<sup>ma</sup><sub>z</sub>/B<sub>3</sub></th></tr></thead><tbody><tr><td>0.002</td><td>9.1 10<sup>2</sup></td><td>1.8 10<sup>3</sup></td><td><b>9.3 10<sup>3</sup></b></td><td>1.6 10<sup>3</sup></td><td><b>3.9 10<sup>3</sup></b></td><td>6.2 10<sup>2</sup></td></tr><tr><td>0.005</td><td>9.6 10<sup>2</sup></td><td>1.8 10<sup>3</sup></td><td><b>9.5 10<sup>3</sup></b></td><td>1.6 10<sup>3</sup></td><td><b>4.3 10<sup>3</sup></b></td><td>7.2 10<sup>2</sup></td></tr><tr><td>0.010</td><td>9.4 10<sup>2</sup></td><td>1.9 10<sup>3</sup></td><td><b>8.4 10<sup>3</sup></b></td><td>1.7 10<sup>3</sup></td><td><b>4.3 10<sup>3</sup></b></td><td>7.6 10<sup>2</sup></td></tr><tr><td>0.025</td><td>1.1 10<sup>3</sup></td><td>2.1 10<sup>3</sup></td><td><b>1.1 10<sup>3</sup></b></td><td>1.6 10<sup>3</sup></td><td><b>4.1 10<sup>3</sup></b></td><td>6.3 10<sup>2</sup></td></tr><tr><td>0.050</td><td>1.2 10<sup>3</sup></td><td>1.9 10<sup>3</sup></td><td><b>8.9 10<sup>3</sup></b></td><td>1.6 10<sup>3</sup></td><td><b>4.6 10<sup>3</sup></b></td><td>7.8 10<sup>2</sup></td></tr></tbody></table>

<table><caption>**Table 4**
Parallelepiped 1 — Quadratic time average of the temporal variation of accelerations.</caption><thead><tr><th rowspan="2">Δt [s]</th><th colspan="3">Translation [m s<sup>-3</sup>]</th><th colspan="3">Rotation [rad s<sup>-3</sup>]</th></tr><tr><th>ū<sub>/B3</sub></th><th>v̇<sub>/B3</sub></th><th>ẇ<sub>/B3</sub></th><th>ṗ<sub>/B3</sub></th><th>q̇<sub>/B3</sub></th><th>ṙ<sub>/B3</sub></th></tr></thead><tbody><tr><td>0.002</td><td><b>4.1</b></td><td><b>4.5</b></td><td>2.7</td><td><b>3.7</b></td><td>1.3</td><td>2.2</td></tr><tr><td>0.005</td><td>4.1</td><td><b>5.5</b></td><td>2.9</td><td><b>4.4</b></td><td>1.5</td><td>2.7</td></tr><tr><td>0.010</td><td>4.0</td><td><b>6.4</b></td><td>2.8</td><td><b>4.6</b></td><td>1.5</td><td>3.1</td></tr><tr><td>0.025</td><td>4.2</td><td><b>6.4</b></td><td>3.4</td><td><b>4.4</b></td><td>1.3</td><td>2.9</td></tr><tr><td>0.050</td><td><b>4.9</b></td><td><b>5.7</b></td><td>2.7</td><td><b>4.1</b></td><td>1.6</td><td>2.7</td></tr></tbody></table>

for the heave $\dot{w}_{/B3}$ and the pitch $\dot{q}_{/B3}$ motions, see Table 3. At the opposite, the highest variations of the accelerations are obtained for the sway $\dot{v}$ and the roll $\dot{p}$ motions which are the ones found for the highest numbers of iterations, see Table 4. From this analysis, it can be deduced that the convergence behaviour (driven by the proposed relaxation operator) is not sensitive to the highest added mass effects but to the highest modification of kinematics. It means that the added mass effects are well taken into account.

The accelerations before and after the relaxation step at $t = 10$ s for $\Delta t = 0.010$ s are presented in Fig. 11. Overshoots are due to the strong added mass effect and the fact that fluid forces do not represent a physical converged state because of the use of a linearised Steklov-Poincaré operator for the fluid problem (cf. Section 4). Nevertheless, it can be seen that the degrees of freedom which are strongly relaxed (heave, roll and pitch) converged in the same way as the others.

These results are also used to perform a new simulation where the motion is not anymore solved but directly imposed to the results previously obtained. In other words, the converged kinematics of the parallelepiped is directly imposed from the first fluid non-linear coupling iteration at each new time step. The FSI problem is then replaced by a fluid computation with an imposed body motion configuration. Thus, to illustrate the convergence properties of the FSI coupling algorithm, the only unknowns left, namely the fluid forces, are compared in Fig. 12. This figure shows that, when the variations of fluid forces are significant, the necessary number of iterations (around 15) to get a converged state is almost the same whether the motion is imposed or solved with an adequate relaxation method like the proposed one. In a few words, with the exception of the
<figure><img src="image_13.png" /><figcaption>Fig. 11. Parallelepiped 1 – Accelerations before (—) and after (—) the relaxation step using the relaxation operator at $t = 5.0 \text{ s} - \Delta t = 0.025 \text{ s}$ (new time step at iteration 1 and after iteration 30).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. Parallelepiped 1 – Fluid forces when motion is solved (—) and imposed (—) at $t = 5.0 \text{ s} - \Delta t = 0.025 \text{ s}$ (new time step at iteration 1 and after iteration 30).</figcaption></figure>

fluid mesh modification step which occurs only one time by time step when the motion is imposed, the computational time can be considered the same.

This application showed that the proposed algorithm is suitable and efficient for typical FSI problems even with large added mass effects. Moreover, it has been shown that the use of anisotropic and sometimes very low relaxation coefficients does not deteriorate the convergence properties of the coupling algorithm. Only diagonal relaxation operators were used in this study. The influence of the extra-diagonal coefficients is studied in the second application.

### 7.3. Application 2

In this application, several density ratios, which are denoted as $d$, are considered. The added mass matrix computed from the proposed method is given in Eq. (34). The corresponding acceleration relaxation operators $R_a^{\text{full}}_{/B_3}$ (cf. Eq. (14)) are given in Table 5. Due to the dimensions and the inertial characteristics of the considered parallelepiped, the heave is coupled to the pitch and the sway is coupled to the yaw. Moreover, the diagonal terms dominate the extra-diagonal terms for the yaw motion whereas they are of the same order for the sway, the heave and the pitch.

$$ \mathbf{MA}_{/B_3} = \begin{bmatrix} 569 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1257 & 0 & 0 & 0 & -1256 \\ 0 & 0 & 11432 & 0 & 11433 & 0 \\ 0 & 0 & 0 & 1633 & 0 & 0 \\ 0 & 0 & 11433 & 0 & 21166 & 0 \\ 0 & -1256 & 0 & 0 & 0 & 2535 \end{bmatrix} \quad (34) $$

In order to independently study the influence of the diagonal coefficients, one additional diagonalised acceleration relaxation operator is defined. The extra diagonal coefficients of the acceleration relaxation operator $R_a^{\text{full}}_{/B_3}$ are arbitrary set to zero. This additional acceleration relaxation operator is denoted as $R_a^{\text{diag}}_{/B_3}$. In a nutshell, the use of $R_a^{\text{full}}_{/B_3}$ is in
<table><caption>**Table 5**

Parallellepipeds 2 — Relaxation operator $\mathbf{R}_a^{\text{full}/\beta_3}$.</caption><thead><tr><th>d</th><th>R<sub>a</sub><sup>full/β<sub>3</sub></sup></th></tr></thead><tbody><tr><td rowspan="5">0.1</td><td>0.413 0 0 0 0 0 0.294</td></tr><tr><td>0 0.330 0 0 0 0.077 0</td></tr><tr><td>0 0 0.068 0 0 0 0</td></tr><tr><td>0 0 0 0.080 0 0 0</td></tr><tr><td>0 0 -0.035 0 0.080 0 0.388</td></tr><tr><td rowspan="5">0.2</td><td>0.584 0 0 0 0 0 0.325</td></tr><tr><td>0 0.468 0 0 0 0 -0.134 0</td></tr><tr><td>0 0 0.123 0 0 0 0</td></tr><tr><td>0 0 0 0.148 0 0 0</td></tr><tr><td>0 0 -0.061 0 0.144 0 0.532</td></tr><tr><td rowspan="5">0.3</td><td>0.678 0 0 0 0 0 0.314</td></tr><tr><td>0 0.552 0 0 0 0 -0.177 0</td></tr><tr><td>0 0 0.168 0 0 0 0</td></tr><tr><td>0 0 0 0.207 0 0 0</td></tr><tr><td>0 0 -0.081 0 0.196 0 0.614</td></tr><tr><td rowspan="5">0.5</td><td>0.779 0 0 0 0 0 0.274</td></tr><tr><td>0 0.656 0 0 0 0 -0.236 0</td></tr><tr><td>0 0 0.240 0 0 0 0</td></tr><tr><td>0 0 0 0.303 0 0 0</td></tr><tr><td>0 0 -0.108 0 0.278 0 0.710</td></tr><tr><td rowspan="5">0.7</td><td>0.831 0 0 0 0 0 0.237</td></tr><tr><td>0 0.719 0 0 0 0 -0.273 0</td></tr><tr><td>0 0 0.296 0 0 0 0</td></tr><tr><td>0 0 0 0.378 0 0 0</td></tr><tr><td>0 0 -0.124 0 0.339 0 0.766</td></tr></tbody></table>

most of the cases the optimum because all the coupling terms are taken into account and the use of $\mathbf{R}_a^{\text{diag}}/\beta_3$ can be seen as a purely numerical decoupling of degrees of freedom.

According to the parallelepiped densities, different behaviours are observed because the ratio between the gravity and the buoyant force changes in the same way as the density. The trajectories and orientations of the different parallelepipeds are given in Fig. 13. In all cases, these trajectories represent 15 s of physical time. The results can slightly depend on the time step because the parallelepipeds are totally free and the flow is unsteady, complex and unstable as it has already been stated in the first application.

The analysis of this case is divided into two parts because similar observations are made when the density is higher or equal than 0.3 and when it is smaller than 0.3 which are extreme cases rarely encountered in naval applications. Firstly, the number of iterations to reach a converged state when the density is higher or equal than 0.3 is given in Table 6. When one relaxation operator produces significantly better results, the results are in bold letters. One can notice that when the degrees of freedom are uncoupled, the numbers of iterations to reach a converged state logically do not depend on the diagonalisation of the relaxation operator. On the contrary, when they are coupled, these numbers of iterations are smaller when all coupling terms are taken into account (cases “full”) and higher when the acceleration relaxation operator is diagonalised (cases “diag1”).

The number of iterations to reach a converged state when the density is smaller than 0.3 is given in Table 7. The previous conclusions can also be made in these cases but only when very small time steps are used. When large time steps are used, the numbers of iterations are smaller when the relaxation operator is diagonalised (cases “diag”) and not any more when all coupling terms are taken into account (cases “full”). More precisely, the diagonalisation is very efficient at large time step when the diagonal terms of the relaxation operator are very low and do not dominate the extra diagonal terms. This is the case for the heave and the pitch motions (cf. Table 5). As regards the two left coupled degrees of freedom, the one whose diagonal coefficient dominates the extra-diagonal one, namely the yaw, converges faster when the full relaxation operator is used. On the contrary, the extra-diagonal coefficient of the sway motion is of the same order than its diagonal coefficient (cf. Table 5) and this motion converges faster when a diagonalisation is performed. Thus, the diagonalisation does not significantly improve the convergence of these two degrees of freedom.

A typical situation where the convergence is slow and when all coupling terms are taken into account is presented in Fig. 14. In this figure, the results are also compared to the ones obtained when a diagonalisation is performed. One can notice that when all coupling terms are taken into account, the coupled degrees of freedom mutually prevent each other from converging because the not relaxed accelerations carry on oscillating. On the contrary, when a numerical diagonalisation is performed, a classical convergence behaviour is recovered.

This convergence problem when using the full relaxation operator does not occur any more when a smaller time step is used because the more the time step becomes small the more the variations of the fluid forces are dominated by the added mass effect (Badia et al., 2008; Durand, 2012). Consequently, if a small time step is used, the variation of accelerations
<figure><img src="image_15.png" /><figcaption>Fig. 13. Parallelepipeds 2 – Trajectory and orientation - $\Delta t = 0.025$ s.</figcaption></figure>

<table><caption>**Table 6**
Parallelepipeds 2 – Criteria of convergence – $d \ge 0.3$.</caption><thead><tr><th rowspan="2">d [-]</th><th rowspan="2">Δt [s]</th><th rowspan="2">Case</th><th colspan="2">Uncoupled DOF</th><th colspan="4">Coupled DOF</th><th rowspan="2">Rating</th></tr><tr><th>Surge<br>ū<sub>/B<sub>3</sub></sub></th><th>Roll<br>p̂<sub>/B<sub>3</sub></sub></th><th>Sway<br>ū̇<sub>/B<sub>3</sub></sub></th><th>Yaw<br>f̂<sub>/B<sub>3</sub></sub></th><th>Heave<br>ū̇<sub>/B<sub>3</sub></sub></th><th>Pitch<br>q̇<sub>/B<sub>3</sub></sub></th></tr></thead><tbody><tr><td rowspan="6">0.3</td><td rowspan="2">0.010</td><td>"full"</td><td>2.4</td><td>12.1</td><td><strong>3.7</strong></td><td><strong>2.1</strong></td><td>2.2</td><td>1.5</td><td><strong>1</strong></td></tr><tr><td>"diag"</td><td>2.4</td><td>12.0</td><td>4.0</td><td>2.7</td><td>2.4</td><td>1.7</td><td>2</td></tr><tr><td rowspan="2">0.025</td><td>"full"</td><td>3.7</td><td>10.2</td><td>4.3</td><td><strong>3.0</strong></td><td><strong>4.2</strong></td><td><strong>2.9</strong></td><td><strong>1</strong></td></tr><tr><td>"diag"</td><td>3.7</td><td>10.1</td><td>4.4</td><td>3.4</td><td>4.9</td><td>3.9</td><td>2</td></tr><tr><td rowspan="2">0.050</td><td>"full"</td><td>5.8</td><td>9.2</td><td>6.8</td><td>5.2</td><td><strong>6.4</strong></td><td><strong>5.2</strong></td><td><strong>1</strong></td></tr><tr><td>"diag"</td><td>5.9</td><td>9.0</td><td>6.7</td><td>5.1</td><td>7.5</td><td>6.5</td><td>2</td></tr><tr><td rowspan="6">0.5</td><td rowspan="2">0.010</td><td>"full"</td><td>3.9</td><td>5.2</td><td>3.3</td><td>2.1</td><td><strong>6.0</strong></td><td><strong>4.4</strong></td><td><strong>1</strong></td></tr><tr><td>"diag"</td><td>3.8</td><td>5.3</td><td>3.1</td><td>2.0</td><td>6.8</td><td>6.3</td><td>2</td></tr><tr><td rowspan="2">0.025</td><td>"full"</td><td>6.1</td><td>7.3</td><td>5.2</td><td>3.6</td><td><strong>7.7</strong></td><td><strong>6.8</strong></td><td><strong>1</strong></td></tr><tr><td>"diag"</td><td>6.2</td><td>7.2</td><td><strong>4.6</strong></td><td><strong>3.2</strong></td><td>8.7</td><td>7.9</td><td>2</td></tr><tr><td rowspan="2">0.050</td><td>"full"</td><td>10.2</td><td>13.2</td><td>10.8</td><td>7.6</td><td><strong>11.9</strong></td><td><strong>10.1</strong></td><td><strong>1</strong></td></tr><tr><td>"diag"</td><td>10.0</td><td>13.0</td><td>10.7</td><td>7.8</td><td>13.7</td><td>12.1</td><td>2</td></tr><tr><td rowspan="6">0.7</td><td rowspan="2">0.010</td><td>"full"</td><td>2.4</td><td>3.1</td><td>1.4</td><td>1.2</td><td><strong>4.3</strong></td><td><strong>3.4</strong></td><td><strong>1</strong></td></tr><tr><td>"diag"</td><td>2.5</td><td>3.3</td><td>1.4</td><td>1.3</td><td>5.1</td><td>4.0</td><td>2</td></tr><tr><td rowspan="2">0.025</td><td>"full"</td><td>3.8</td><td>5.3</td><td>2.9</td><td>2.1</td><td><strong>6.2</strong></td><td><strong>5.4</strong></td><td><strong>1</strong></td></tr><tr><td>"diag"</td><td>3.9</td><td>5.1</td><td>2.7</td><td>2.0</td><td>7.4</td><td>6.8</td><td>2</td></tr><tr><td rowspan="2">0.050</td><td>"full"</td><td>6.5</td><td>8.4</td><td>5.9</td><td>3.5</td><td><strong>8.5</strong></td><td><strong>7.5</strong></td><td><strong>1</strong></td></tr><tr><td>"diag"</td><td>6.5</td><td>8.5</td><td>5.6</td><td>3.6</td><td>9.4</td><td>8.5</td><td>2</td></tr></tbody></table>

between two time steps is dominated by the added mass effect and consequently, well taken into account by the proposed full relaxation operator (cf. Fig. 15). Again, the different converged states between the "full" and "diag" cases are due to the accumulation of small convergence error in time.

This application showed that taking into account all coupling terms to compute the relaxation operator decreases the number of iterations needed to reach a converged state in most cases. A diagonalisation of the relaxation operator must be
<table><caption>**Table 7**

Parallelepipeds 2 — Criteria of convergence — *d* < 0.3.</caption><thead><tr><th rowspan="2">d [-]</th><th rowspan="2">Δt [s]</th><th rowspan="2">Case</th><th colspan="2">Uncoupled DOF</th><th colspan="4">Coupled DOF</th><th rowspan="2">Rating</th></tr><tr><th>Surge<br>ū/B<sub>3</sub></th><th>Roll<br>ρ̇/B<sub>3</sub></th><th>Sway<br>ū̇/B<sub>3</sub></th><th>Yaw<br>ρ̇/B<sub>3</sub></th><th>Heave<br>ω̇/B<sub>3</sub></th><th>Pitch<br>q̇/B<sub>3</sub></th></tr></thead><tbody><tr><td rowspan="6">0.1</td><td rowspan="2">0.005</td><td>"full"</td><td>6.2</td><td>19.5</td><td>6.8</td><td><strong>4.3</strong></td><td>4.6</td><td><strong>2.3</strong></td><td><strong>1</strong></td></tr><tr><td>"diag"</td><td>7.7</td><td>26.2</td><td>6.3</td><td>6.5</td><td>4.5</td><td>3.5</td><td>2</td></tr><tr><td rowspan="2">0.010</td><td>"full"</td><td>8.7</td><td>26.6</td><td>9.5</td><td><strong>6.2</strong></td><td>14.1</td><td>6.0</td><td>2</td></tr><tr><td>"diag"</td><td>8.7</td><td>26.6</td><td><strong>8.5</strong></td><td>7.9</td><td><strong>6.5</strong></td><td>6.0</td><td><strong>1</strong></td></tr><tr><td rowspan="2">0.025</td><td>"full"</td><td>9.4</td><td>19.2</td><td>13.2</td><td><strong>9.4</strong></td><td>19.9</td><td>12.1</td><td>2</td></tr><tr><td>"diag"</td><td>9.5</td><td>18.0</td><td><strong>11.6</strong></td><td>11.1</td><td><strong>8.9</strong></td><td><strong>8.2</strong></td><td><strong>1</strong></td></tr><tr><td rowspan="6">0.2</td><td rowspan="2">0.010</td><td>"full"</td><td>6.9</td><td>15.0</td><td>6.5</td><td><strong>4.2</strong></td><td>5.2</td><td><strong>3.2</strong></td><td><strong>1</strong></td></tr><tr><td>"diag"</td><td>6.9</td><td>14.9</td><td>5.8</td><td>5.2</td><td>5.2</td><td>4.4</td><td>2</td></tr><tr><td rowspan="2">0.025</td><td>"full"</td><td>6.8</td><td>15.8</td><td>9.2</td><td><strong>6.5</strong></td><td>10.1</td><td><strong>5.4</strong></td><td>2</td></tr><tr><td>"diag"</td><td>6.7</td><td>15.9</td><td>9.1</td><td>7.4</td><td><strong>7.2</strong></td><td>5.9</td><td><strong>1</strong></td></tr><tr><td rowspan="2">0.050</td><td>"full"</td><td>8.3</td><td>15.0</td><td>12.1</td><td><strong>8.6</strong></td><td>18.2</td><td>12.0</td><td>2</td></tr><tr><td>"diag"</td><td>7.8</td><td>14.2</td><td><strong>10.8</strong></td><td>10.0</td><td><strong>8.5</strong></td><td><strong>7.9</strong></td><td><strong>1</strong></td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Parallelepipeds 2 — Accelerations in heave and pitch before the relaxation step for "full" case (---) and "diag" case (···) and after the relaxation step for "full" case (—) and "diag" case (—) at $t = 8.25 \text{ s} - \Delta t = 0.025 \text{ s} - d = 0.1$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Parallelepipeds 2 — Accelerations in heave and pitch before the relaxation step for "full" case (---) and "diag" case (···) and after the relaxation step for "full" case (—) and "diag" case (—) at $t = 8.25 \text{ s} - \Delta t = 0.005 \text{ s} - d = 0.1$.</figcaption></figure>

preferred under the following conditions: very large added mass effects (very small relaxation coefficients), no domination of the diagonal over the extra-diagonal coefficients of the relaxation operator and large time step.

## 8. Conclusion

In this paper, an efficient coupling algorithm for fluid-structure interaction problems with large added mass effect is presented. The key-points of this coupling algorithm are the replacement of the Steklov-Poincaré operator by a linearised one and the computation of an added-mass equation to correctly evaluate the relaxation operator. When the added mass matrix changes during the simulation time, a dynamic evaluation can be performed because only one iteration for each degree of freedom is enough to correctly estimate the added mass effects.

It was shown that for extreme cases, the relaxation operator requires to be diagonalised but for typical problems encountered in naval architecture, the full relaxation operator can be used because a converged state is obtained sooner. Moreover, accessing to the added mass information within the fluid solver is very interesting because it can automatically
be performed during the simulation time when dealing with unsteady configurations. Moreover, the classical potential solver limitations do not exist any more with a finite-volume solver. Another key-point of this coupling algorithm is that the implementation can be performed in a co-simulation context, which means that the structural and the fluid solver can be used as black boxes. The only requirement is that the fluid solver can communicate with another solver during its non linear resolution process. Finally, comparison between solved and imposed motion simulations shows excellent convergence properties because the number of iterations to reach a converged state is almost the same in both cases. Regardless of the efficiency of the code to solve the fluid flow with imposed moving boundaries without feedback of the structure, the FSI algorithm presented here can be then considered as optimal in term of convergence.

This coupling algorithm has been extensively validated with different and complex maritime applications (Yvin, 2014). The extension of the proposed method to flexible bodies requires additional works. Indeed, the definition of a relaxation operator for a flexible body based on the added mass effect is not straightforward.

## Acknowledgements

This work was supported by the ANRT (the French National Association for Research and Technology) through a CIFRE convention (589/2011) and was granted access to the HPC resources of IDRIS under the allocations 2013-21308 and 2014-21308 made by GENCI (Grand Equipement National de Calcul Intensif).

## References

* Badia, S., Codina, R., 2000. On some fluid-structure iterative algorithms using pressure segregation methods. Application to aeroelasticity. Int. J. Numer. Methods Eng. 72, 46-71.

* Badia, S., Quaini, A., Quarteroni, A., 2008. Modular vs. non-modular preconditioners for fluid-structure systems with large added-mass effect. Comput. Methods Appl. Mech. Engrg. 197 (49-50), 4216-4232.

* Belanger, P., Faïdoussis, M.P., De Langre, E., 1995. Time-marching analysis of fluid-coupled systems with large added mass. Amer. Inst. Aeronaut. Astronaut. (AIAA) J. 33 (4), 752-757.

* Brennen, C.E., 1982. A Review of Added Mass and Fluid Inertial Forces, Tech. Rep. N62583-81-MR-554. Naval Civil Engineering Laboratory, Port Hueneme, California.

* Brennen, C.E., 2008. Fundamentals of Multiphase Flows. Cambridge University Press.

* Cervera, M., Codina, R., Galindo, M., 1996. On the computational efficiency and implementation of block-iterative algorithms for nonlinear coupled problems. Eng. Comput. 13 (6), 4-30.

* Degroote, J., Bathe, K.-J., Vierendeels, J., 2009. Performance of a new partitioned procedure versus a monolithic procedure in fluid-structure interaction. Comput. Struct. 87 (11-12), 793-801.

* Delhommeau, G., July 1987. Les problèmes de diffraction-radiation et de résistance de vagues: étude théorique et résolution numérique par la méthode des singularités (Ph.D. thesis), École Nationale Supérieure de Mécanique.

* Deparis, S., Discacciati, M., Fourestey, G., Quarteroni, G., 2006. Fluid-structure algorithms based on Steklov-Poincaré operators. Comput. Methods Appl. Mech. Engrg. 195 (41-43), 5797-5812.

* Dettmer, W., Perić, D., 2007. A fully implicit computational strategy for strongly coupled fluid-solid interaction. Arch. Comput. Methods Eng. 14, 205-247.

* Dettmer, W., Perić, D., 2008. On the coupling between fluid flow and mesh motion in the modelling of fluid-structure interaction. Comput. Mech. 43, 81-90.

* Durand, M., 2012. Interaction fluide-structure souple et légère, applications aux voiliers (Ph.D. thesis), École Centrale de Nantes.

* Förster, C., Wall, W., Ramm, E., 2007. Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows. Comput. Methods Appl. Mech. Engrg. 196 (7), 1278-1293.

* Ghiringhelli, G.L., Masarati, P., Mantegazza, P., Nixon, M.W., 1999. Multi-body analysis of a tiltrotor configuration. Nonlinear Dynam. 19 (4), 333-357.

* Guyon, E., Hulin, J.-P., Petit, L., Mitescu, C., 2001. Physical Hydrodynamics. Oxford University Press, p. 528.

* Heil, M., 2004. An efficient solver for the fully-coupled solution of large-displacement fluid-structure interaction problems. Comput. Methods Appl. Mech. Engrg. 193, 1-23.

* Hou, G., Wang, J., Layton, A., 2012. Numerical methods for fluid-structure interaction - A review. Commun. Comput. Phys. 12 (2), 337-377.

* Hübner, B., Walhorn, E., Dinkler, D., 2004. A monolithic approach to fluid-structure interaction using space-time finite elements. Comput. Methods Appl. Mech. Engrg. 193 (23-26), 2087-2104.

* Jürgens, D., 2009. Survey on Software Engineering for Scientific Applications: Reuseable Software, Grid Computing and Application, Tech. Rep. Université technique Carolo-Wilhelmina de Brunswick.

* Kassiotis, C., 2009. Nonlinear Fluid-Structure Interaction: A Partitioned Approach and its Application Through Component Technology (Ph.D. thesis), Université Paris-Est, École Doctorale MODES.

* Küttler, U., Wall, W., 2008. Fixed-point fluid-structure interaction solvers with dynamic relaxation. Comput. Mech. 43, 61-72.

* Landau, L.D., Lifshitz, E.M., 1987. Fluid Mechanics. Course of Theoretical Physics. Pergamon Press.

* Leonard, A., Roshko, A., 2001. Aspects of flow-induced vibration. J. Fluids Struct. 15, 415-425.

* Leroyer, A., 2004. Étude du couplage écoulement/mouvement pour des corps solides ou à déformation imposée par résolution des équations de Navier-Stokes. Contribution à la modélisation numérique de cavitation (Ph.D. thesis), École Centrale de Nantes.

* Leroyer, A., Barré, S., Kobus, J.M., Visonneau, M., 2008. Experimental and numerical investigations of the flow around an oar blade. J. Mar. Sci. Technol. 13 (1), 1-15.

* Leroyer, A., Visonneau, M., 2005. Numerical methods for RANSE simulations of a self-propelled fish-like body. J. Fluids Struct. 20 (7), 975-991.

* Masarati, P., 1999. Comprehensive Multibody AeroServoElastic Analysis of Integrated Rotorcraft Active Controls (Ph.D. thesis), Dipartimento Di Ingegneria Aerospaziale, Politecnico Di Milano.

* Papadakis, G., 2008. A novel pressure-velocity formulation and solution method for fluid-structure interaction problems. J. Comput. Phys. 227 (6), 3383-3404.

* Queutey, P., Deng, G.B., Wackers, J., Guilmineau, E., Leroyer, A., Visonneau, M., 2012. Sliding grids and adaptive grid refinement for RANS simulation of ship-propeller interaction. Ship Technol. Res. 59 (2), 44-58.

* Queutey, P., Visonneau, M., 2007. An interface capturing method for free-surface hydrodynamic flows. Comput. & Fluids 36 (9), 1481-1510.

* Saksono, P.H., Dettmer, W., Peric, D., 2007. An adaptive remeshing strategy for flows with moving boundaries and fluid-structure interaction. Internat. J. Numer. Methods Engrg. 71, 1009-1050.
Söding, H., 2001. How to integrate free motions of solids in fluids. In: 4th Numerical Towing Tank Symposium, Hamburg, Germany.

Uzunoğlu, B., M., T., G., P.W., 2001. Low-Reynolds-number flow around an oscillating circular cylinder using a cell viscous boundary element method. Int. J. Numer. Methods Eng. 50, 2317–2338.

Vikestad, K., Vandiver, J.K., Larsen, C.M., 2000. Added mass and oscillation frequency for a circular cylinder sujected to vortex-induced vibrations and external disturbance. J. Fluids Struct. 14 (7), 1071–1088.

Wakaba, L., Balachandar, S., 2007. On the added mass force at finite Reynolds and acceleration numbers. Theor. Comput. Fluid Dyn. 21 (2), 147–153.

Wall, W.A., Förster, C., Neumann, M., Ramm, E., 2006. Advances in fluid-structure interaction. In: 17th International Conference on the Application of Computer Science and Mathematics in Architecture and Civil Engineering, Weimar, Allemagne.

Wick, T., 2011. Fluid-structure interactions using different mesh motion techniques. Comput. Struct. 89 (13–14), 1456–1467.

Yvin, C., 2014. Fluid-Structure Interaction for Multi-Body Systems. Applications to Complex Joints, Control laws of Actuators and Flexible Systems in the Maritime Field (Ph.D. thesis), ULAM Ecole Centrale de Nantes.