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

HAL
open science

Numerical simulation of a non linear coupled
fluid-structure problem by explicit finite element-finite
volume coupling

Jean-François Sigrist, Christian Lainé, Bernard Peseux

► To cite this version:

Jean-François Sigrist, Christian Lainé, Bernard Peseux. Numerical simulation of a non linear coupled fluid-structure problem by explicit finite element-finite volume coupling. 2005 ASME Pressure Vessel and Piping Conference, Jul 2005, Denver, United States. 10.1115/PVP2005-71227. hal-01008557

HAL Id: hal-01008557

https://hal.science/hal-01008557

Submitted on 6 Nov 2018

**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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
NUMERICAL SIMULATION OF A NON LINEAR COUPLED FLUID-STRUCTURE PROBLEM
BY EXPLICIT FINITE ELEMENT-FINITE VOLUME COUPLING

Jean François SIGRIST

Service Scientifique et Technique
DCN Propulsion

44620 LA MONTAGNE, France

jean-francois sigrist@dcn.fr

Christian LAINE

Service Scientifique et Technique
DCN Propulsion

44620 LA MONTAGNE, France

Bernard PESEUX

G&M - Pole Dynamique et Couplages

Ecole Centrale de Nantes - 1, rue de La Noë

44321 Nantes Cedex 3, France

ABSTRACT

The present paper deals with the numerical simulation of a coupled non linear
fluid-structure problem by explicit coupling between a finite element structure
code and a finite volume fluid code. This numerical study is carried out in
order to develop robust and general coupling with FE and CFD commercial
code for industrial applications.

A geometrically simple non linear coupled problem is presented in order to
validate the numerical approach. The structure non linear problem is solved
with a finite element technique, using a iterative implicit algorithm for time
integration. The fluid problem is solved using standard numerical techniques
(finite volume approach, implicit splitting operator scheme). The whole
coupled problem is solved with a commercial CFD code: a dedicated FE
structure code is developed in the CFD code together with coupling (in time,
in space) procedures.

The proposed method is validated in the case of a incompressible inviscid
fluid, for which the coupled problem is solved with an analytical solution. The
present study gives a reference test case for a full scale fluid-structure model.
Industrial applications can now be considered by coupling commercial FE and
FV codes with general coupling code.

INTRODUCTION

The present paper deals with the numerical simulation of
coupled fluid-structure problem, in which the coupling process is
bases on mechanical exchanges between each sub-problem. In the
present study other fluid-structure coupling phenomena, such as
thermal effect, are not taken into account. At the fluid structure
interface, fluid and structure are coupled via fluid forces and structure
displacement, as sketched by Fig. 1.

Over the past years many numerical methods have been
developed in order to solve linear and non linear coupled fluid-

structure problems, using finite element and boundary element
techniques [4] [15] [16].

<figure><img src="image_3.png" /><figcaption>Figure 1. Mechanical Fluid/structure coupling [20]</figcaption></figure>

Nowadays, the numerical approach tends to propose general
coupling algorithms with specific solvers for each sub-problems, with
various coupling strategies [5] [8] [17] [20], depending on the
physical coupling phenomenon.

As the numerical design naval propulsion structures needs
taking fluid/structure interaction into account, DCN Propulsion
launched, with Ecole Centrale Nantes G&M, a R&D study in order to
apply numerical methods on industrial projects. From the academic
point of view, the coupling process seems not to pose any particular
dificulties. From the industrial point of view, an investigation of the
various aspects of the numerical fluid-structure coupling on a simple
case is necessary before considering a full-scale coupling with
existing commercial codes for industrial structure design.

The fluid-structure coupled problem studied in the present
paper is described by Fig. 2.<fnref n="1" /> The numerical resolution of the coupled

<footnote n="1">¹ It is an extension in the 3D case of the non linear 2D problem studied in a previous publication [21].</footnote>
problem is performed using an explicit coupling techniques with a CFD fluid code and FE structure code, as represented by Fig. 3.

<figure><img src="image_4.png" /><figcaption>Figure 2. 1D structure-3D fluid non linear coupled problem</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3. Fluid and structure code coupling: general principle</figcaption></figure>

## 1. STRUCTURE PROBLEM

The structure problem is described by the following beam non linear equation of motion, taking into account geometrical non linearity [12] [26]. The coupled reaction/bending equation of motion are formulated in the relative frame as:

$$ \rho_s S \frac{\partial^2 u}{\partial t^2} + \frac{\partial}{\partial z} \left\{ ES \left[ \frac{\partial u}{\partial z} + \frac{1}{2} \left( \frac{\partial v}{\partial z} \right)^2 \right] \right\} = 0 \quad (1) $$

$$ \rho_s S \frac{\partial^2 v}{\partial t^2} + \frac{\partial^2}{\partial z^2} \left( EI \frac{\partial^2 v}{\partial z^2} \right) - \frac{\partial}{\partial z} \left\{ ES \left[ \frac{\partial u}{\partial z} + \frac{1}{2} \left( \frac{\partial v}{\partial z} \right)^2 \right] \frac{\partial v}{\partial z} \right\} = - \rho_s S \gamma \quad (2) $$

The boundary conditions are:

$$ u_{z=0} = 0 \qquad v_{z=0} = 0 \qquad \frac{\partial v}{\partial z_{z=0}} = 0 \quad (3) $$

for the beam clamped end, and:

$$ \left( \frac{\partial u}{\partial z} + \frac{1}{2} \left( \frac{\partial v}{\partial z} \right)^2 \right)_{z=L} = 0 \qquad \frac{\partial^2 v}{\partial z^2} \Big|_{z=L} = 0 \qquad \frac{\partial^3 v}{\partial z^3} \Big|_{z=L} = 0 \quad (4) $$

for the beam free end.

### 1.1. Space discretisation: finite element method

The numerical resolution of the structure problem is performed using the finite element method [1]. The variational formulation of the non linear problem is obtained by multiplying Eq. (1) and (2) by virtual displacements $\delta u$ and $\delta v$, integrating the equations over $[0, L]$, using part integration, leading to:

$$ \int_{0}^{L} \rho_s S \frac{\partial^2 u}{\partial t^2} \delta u + \int_{0}^{L} ES \frac{\partial u}{\partial z} \frac{\partial \delta u}{\partial z} + \int_{0}^{L} \frac{ES}{2} \left( \frac{\partial v}{\partial z} \right)^2 \frac{\partial \delta u}{\partial z} = 0 \quad (5) $$

$$ \int_{0}^{L} \rho_s S \frac{\partial^2 v}{\partial t^2} \delta v + \int_{0}^{L} EI \frac{\partial^2 v}{\partial z^2} \frac{\partial \delta v}{\partial z} + \int_{0}^{L} ES \left[ \frac{\partial u}{\partial z} + \frac{1}{2} \left( \frac{\partial v}{\partial z} \right)^2 \right] \frac{\partial v}{\partial z} \frac{\partial \delta v}{\partial z} = 0 \quad (6) $$

$$ \int_{0}^{L} ES \left[ \frac{\partial u}{\partial z} + \frac{1}{2} \left( \frac{\partial v}{\partial z} \right)^2 \right] \frac{\partial v}{\partial z} \frac{\partial \delta v}{\partial z} = - \int_{0}^{L} \rho_s S \gamma \delta v \quad (7) $$

The space discretisation uses 2-nodes/3-degrees of freedom beam finite element, as represented by Fig. 4. The unknown displacement are approximated with linear shape functions for $u$, and cubic shape functions for $v$, that is writing $u = \langle N_u \rangle \{U\}$ and $v = \langle N_v \rangle \{V\}$ with:

$$ N_u = \langle 1 - s, s \rangle \quad (8) $$

$$ N_v = \langle 1 - 3s^2 + 2s^3, I(s - 2s^2 + s^3), 3s^2 - 2s^3, I(-s^2 + s^3) \rangle \quad (9) $$

<figure><img src="image_6.png" /><figcaption>Figure 4. Beam finite element with two nodes, three degrees of freedom per node</figcaption></figure>
Using a GALERIN method (i.e. writing an identical approximation for unknown displacement u, v and virtual displacement with the same shape functions), the variational principle given by Eqs. (5) and (6) is written in the following discrete form:

$$ \begin{bmatrix} M_u & 0 \\ 0 & M_v \end{bmatrix} \begin{Bmatrix} \ddot{U} \\ \ddot{V} \end{Bmatrix} + \begin{bmatrix} K_u & 0 \\ 0 & K_v \end{bmatrix} \begin{Bmatrix} U \\ V \end{Bmatrix} + \begin{bmatrix} R_u(V) \\ R_v(U) \end{bmatrix} = \begin{bmatrix} 0 \\ \Phi_v \end{bmatrix} \quad (10) $$

In Eq. (10), the mass matrices are given by $[M_s] = \int_0^L \rho_s S \{N_s\} \{N_s\}$.

and $[M_v] = \int_0^L \rho_s S \{N_v\} \{N_v\}$, the stiffness matrices are expressed as

$[K_s] = \int_0^L ES \left\{ \left( \frac{\partial N_u}{\partial z} \right) \left\{ \frac{\partial N_u}{\partial z} \right\} \right\}$ and $[K_v] = \int_0^L EI \left\{ \left( \frac{\partial^2 N_v}{\partial z^2} \right) \left\{ \frac{\partial^2 N_v}{\partial z^2} \right\} \right\}$.

and the loading vector is written $\{\Phi_s\} = -\int_0^L \rho_s S \gamma \{N_s\}$. With the shape functions defined by Eqs. (8) and (9), these matrices can be analytically calculated. The non linear terms are:

$$ \{R_s(V)\} = \int_0^L \frac{ES}{2} \left\{ \left( \frac{\partial N_u}{\partial z} \right) \langle V \rangle \left\{ \left( \frac{\partial N_v}{\partial z} \right) \left\{ \frac{\partial N_u}{\partial z} \right\} \right\} \langle V \rangle \right\} \quad (11) $$

$$ \begin{aligned} \{R_s(U, V)\} &= \int_0^L ES \left\{ \left( \frac{\partial N_v}{\partial z} \right) \left\{ \left\{ \left( \frac{\partial N_u}{\partial z} \right) \{U\} \right\} + \frac{1}{2} \langle V \rangle \left\{ \left( \frac{\partial N_u}{\partial z} \right) \left\{ \left( \frac{\partial N_v}{\partial z} \right) \langle V \rangle \right\} \left\{ \left( \frac{\partial N_v}{\partial z} \right) \langle V \rangle \right\} \right\} \right\} \end{aligned} \quad (12) $$

## 1.2. Time integration: non linear implicit scheme

The dynamic non linear system given by Eq. (10) can be formulated as follows:

$$ M \ddot{X} + K(X) = \Phi - M \Gamma \quad (13) $$

with $K(X)$ denoting the non linear stiffness, $\Phi$ external forces$^2$, and $-M\Gamma$ being the inertial term$^3$. If fluid forces are taken into account, $\Phi$ is obtained from the discrete form of the following virtual work $\delta W = \int_{\Gamma(t)} \delta u \cdot \bar{\sigma} \cdot \bar{n} d\Gamma$ where $\bar{\sigma}$ is the fluid stress tensor, $\bar{n}$ is the outer normal at the fluid-structure interface $\Gamma$.

The numerical time integration of Eq. (12) can be performed with classical techniques, using explicit or implicit techniques.

In explicit techniques [7], a three points numerical scheme is used to evaluate the acceleration at time $t_n$. For instance, the centered finite difference scheme gives the approximation:

$$ \ddot{X}_n = \frac{X_{n+1} - 2X_n + X_{n-1}}{\delta t^2} \quad (14) $$

Using the expression of Eq. (14) in the general equation of motion leads to:

$$ X_{n+1} = \delta t^2 M^{-1} \times [\phi_n - M\Gamma_n - K_n] + 2X_n - X_{n-1} \quad (15) $$

In Eq. (15), all terms in the right side are known, allowing a straightforward calculation of $X_{n+1}$. The explicit technique can be coupled with the finite volume approach of the fluid problem in a staggered resolution of the coupled problem. Though the implementation of the coupling is rather easygoing in the explicit approach, it suffers from stability limitations. This leads to the use of an implicit technique.

In implicit technique [23], a two points numerical scheme is used; the NEWMARK scheme gives for example:

$$ \ddot{X}_{n+1} = \frac{X_{n+1} - X_n}{\beta \delta t^2} - \frac{\dot{X}_n}{\beta \delta t} - \left(\frac{1}{2\beta} - 1\right) \ddot{X}_n \quad (16) $$

The substitution of Eq. (16) in the equation of motion leads to an implicit equation in terms of the unknown displacement $X_{n+1}$. An iterative procedure based on fixed point algorithm is performed to obtain $X_{n+1}$. The algorithm starts with a linear calculation (with the linear part of $K(X)$), which gives an estimation of the structure displacement $X_{n+1}^0$. Iterations are performed with the following approximation:

$$ K(X_{n+1}^0) = K(X_{n+1}^i) + \left( \frac{\partial K}{\partial X_{n+1,j}} \right) \cdot (X_{n+1}^i - X_{n+1}^j) + O(\delta X_{n+1}^i)^2 \quad (17) $$

which gives:

$$ X_{n+1}^{i+1} = \tilde{K}_{n+1,j}^{-1} \times [\Phi_{n+1} + K(X_{n+1}^i) + \left( \frac{\partial K}{\partial X_{n+1,j}} \right) \cdot X_{n+1}^i + \frac{M}{\beta \delta t^2} X_n + \frac{M}{\beta \delta t} \dot{X}_n + M \left( \frac{1}{2\beta-1} \right) \ddot{X}_n] \quad (18) $$

where $\tilde{K}_{n+1,j} = \frac{M}{\beta \delta t^2} + \left( \frac{\partial K}{\partial X_{n+1,j}} \right)$. The iteration process is stopped

when a convergence criterion is satisfied; a simple criterion can be written as:

$$ \frac{\| X_{n+1}^{i+1} - X_{n+1}^{i} \|}{\| X_{n+1}^{i} \|} \le \varepsilon \quad (19) $$

<footnote>$^2$ In the present paper, the external forces are the fluid forces, which will be computed from a CFD code using a finite volume approach, see §.2.</footnote>

<footnote>$^3$ The proposed $-M\Gamma$ notation in Eq. (12) is an alternative way of writing the dynamic loading $\{\Phi_s\}$ appearing in Eq. (10).</footnote>
where $\epsilon \ll 1$. This approach is equivalent to a NEWTON procedure, and converges under known assumptions. The non linear scheme requires the calculation of the tangent stiffness matrix $\frac{\partial K}{\partial X}$. In the present case, the matrix is expressed as:

$$ \frac{\partial K}{\partial X} = \begin{bmatrix} K_u & \frac{\partial R_v}{\partial V} \\ \frac{\partial R_v}{\partial U} & K_v + \frac{\partial R_v}{\partial V} \end{bmatrix} \qquad (20) $$

with:

$$ \frac{\partial R_v}{\partial V} = \int_0^L ES \left\{ \left\langle \frac{\partial N_v}{\partial z} \right\rangle \left\langle V \right\rangle \left\{ \left\langle \frac{\partial N_v}{\partial z} \right\rangle \left\langle \frac{\partial N_v}{\partial z} \right\rangle \right\} \right\rangle \quad (21) $$

$$ \frac{\partial R_v}{\partial U} = \int_0^L ES \left\{ \left\langle \frac{\partial N_v}{\partial z} \right\rangle \left\langle \frac{\partial N_v}{\partial z} \right\rangle \left\{ \left\langle \frac{\partial N_v}{\partial z} \right\rangle V \right\} \right\} \quad (22) $$

and:

$$ \frac{\partial R_v}{\partial V} = \int_0^L ES \left\{ \left\langle \frac{\partial N_v}{\partial z} \right\rangle \left\langle V \right\rangle \left\{ \left\langle \frac{\partial N_v}{\partial z} \right\rangle \left\langle \frac{\partial N_v}{\partial z} \right\rangle \left\langle \frac{\partial N_v}{\partial z} \right\rangle \right\} \right\rangle + \int_0^L ES \left\{ \left\langle \frac{\partial N_v}{\partial z} \right\rangle \left\langle \frac{\partial N_v}{\partial z} \right\rangle \left\langle U \right\rangle \left\langle \frac{\partial N_v}{\partial z} \right\rangle \right\rangle + \int_0^L ES \left\{ \left\langle \frac{\partial N_v}{\partial z} \right\rangle \left\langle V \right\rangle \left\{ \left\langle \frac{\partial N_v}{\partial z} \right\rangle \left\langle \frac{\partial N_v}{\partial z} \right\rangle \left\langle V \right\rangle \left\langle \frac{\partial N_v}{\partial z} \right\rangle \right\rangle \right\} \quad (23) $$

The non-linear terms of these matrices are computed using a GAUSS-LEGENDRE procedure [1], that is using a numerical integration scheme based on the following approximation:

$$ \int_{1}^{+1} \psi(\zeta) d\zeta \approx \sum_{i=1}^{L} w_i \psi(\zeta_i) \quad (24) $$

with $\zeta_i$ Gauss integration point and $w_i$ the associated integration weight. The numerical scheme (24) using L points is exact for the polynomial functions of degree $2L-1$: the numerical integration procedure used in the present case will be based on a 5 points scheme. Furthermore, the iterative scheme given by Eq. (18) needs inverting matrix $\frac{M}{\beta\delta t^2} + \left(\frac{\partial K}{\partial X}\right)_{a+1,j}$ at each time step, which is done with a LU factorization.

The above finite element procedure is implemented as FORTRAN subroutines, which will be integrated as a simple structural code in the commercial CFD code used for the numerical resolution of the fluid problem. Coupling procedures presented in §.3 will also be implemented as FORTRAN subroutine of the CFD code.

## 2. FLUID PROBLEM

The fluid problem is described by the general NAVIER-STOKES equation. The conservation equation are integrated over a moving control volume $\Omega(t)$ of boundary $\tilde{\Omega}(t)$, using the LIEBNITZ rule

and the GAUSS theorem. This leads to the global conservation equation over an arbitrary control volume [25]. in ALE formulation [19]:

* Mass conservation

$$ \frac{d}{dt} \left( \int_{\Omega(t)} \rho d\Omega \right) + \int_{\Omega(t)} \rho (v_j - v_j^*) dn_j = \int_{\Omega(t)} S_{\rho} d\Omega \quad (25) $$

* Momentum conservation

$$ \frac{d}{dt} \left( \int_{\Omega(t)} \rho v_i d\Omega \right) + \int_{\Omega(t)} \rho v_i (v_j - v_j^*) dn_j = - \int_{\Omega(t)} p d\nu n_j + \int_{\Omega(t)} \mu \left( \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i} \right) dn_j + \int_{\Omega(t)} S_{\rho} d\Omega \quad (26) $$

In Eqs. (26) and (26), $S_\rho$ and $S_{v_j}$ stand respectively for mass and momentum sources: the fluid problem can then be formulated in the moving frame. The fluid unknown are the pressure and velocity fields $p, v_j$, but as the conservation equation are integrated over a moving control volume, Eqs. (25) and (26) show another unknown fluid that is the grid velocity $v_j^*$: another equation has to be used to close the fluid problem. This supplementary equation is the space law conservation [6] which is written:

$$ \frac{d}{dt} \left( \int_{\Omega(t)} d\Omega \right) = \int_{\Omega(t)} v_j^* dn_j \quad (27) $$

The fluid problem is then fully characterized by Eqs. (25) to (27)$^4$. Additional boundary conditions are to be taken into account such as imposed pressure or velocity. The coupling with the structure problem is expressed as $\bar{v}|_r = \frac{\partial \bar{u}}{\partial t}|_r$ where $\bar{u}$ is the structure displacement and $\Gamma$ is the fluid-structure interface.

### 2.1. Space discretisation: finite volume method

The general form of an integrated conservation equation for a fluid unknown $\phi$ is the following one:

$$ \frac{d}{dt} \left( \int_{\Omega(t)} \rho \phi d\Omega \right) + \int_{\Omega(t)} \rho (v_j - v_j^*) \phi dn_j = \int_{\Omega(t)} S_{\rho} d\Omega + \int_{\Omega(t)} \Gamma \frac{\partial \phi}{\partial x_j} dn_j \quad (28) $$

Equation (28) is discretised using a finite volume method. The fluid domain is divided in elementary fluid cells (Fig. 5 gives a typical representation of a computational grid for a 2D cartesian problem), for which a integrated conservation equation is written. A moving mesh technique is used, the elementary grid velocity is deduced from known node displacement in order to satisfy Eq. (27).

<footnote>$^4$ In the case of a compressible fluid, the fluid state law $\rho(p)$ allows a closure of the problem.</footnote>
[6]. The grid velocity $v_i^*$ appearing in the right side of Eq. (28) becomes then a source term on the left hand side. Integrated conservation equation are written for each cell: P is the cell center, WW, W, SW, S, SE, E, EE, NE, N are the centers of neighboring cells.

<figure><img src="image_7.png" /><figcaption>Figure 5. 2D fluid computational grid</figcaption></figure>

The time dependant in Eq. (28) are approximated by an EULER scheme, i.e. writing:

$$ \frac{d}{dt} \left( \int_{\Omega_P(t)} \rho \phi d\Omega \right) \approx \frac{\rho \Omega}{\delta t} (\phi_P^{n+1} - \phi_P^{n-1}) \quad (29) $$

The convective and diffusive terms are calculated using finite difference approximation, leading to the global expression [9]:

$$ \int_{\Omega_P(t)} \rho v_j \phi - \Gamma \frac{\partial \phi}{\partial x_j} dn_j \approx \sum_M C_M^{n-1} \phi_M^{n-1} - D_M^{n-1} \phi_M^{n-1} \quad (30) $$

where M stands for the neighboring points of cell $\Omega_P$.

Taking into account Eqs. (29) and (30) leads to the following algebraic non linear system $[A(\phi)]\{\phi\} = \{b\}$, that is:

$$ A_p \phi_p^{n+1} = \sum_M A_p \phi_p^{n+1} + b_p \phi_p^n \quad (31) $$

## 2.2. Time integration: PISO algorithm

From discretisation practices presented above, it is possible to produce the following algebraic system resulting from the finite volume discretisation of the momentum equation, Eq. (32), and PÖSSON equations, Eq. (33).

$$ A_p v_{i,P}^{n+1} = H(v_i^{n+1}) + D(p^{n+1}) + S_v^{n-1} \quad (32) $$

$$ A_p p_p^{n+1} = G(p^{n+1}) + S_p(v_i^{n+1}) \quad (33) $$

Equations (32) and (33) couple the pressure and velocity unknowns. An implicit iterative scheme is used to calculate the new

values of the velocity and pressure field [13]. The algorithm starts with an estimation of the velocity field, by solving the following non linear system, obtained from Eq. (33) by using the pressure field computed at the previous time step:

$$ A_p v_{i,P}^{n+1,(0)} = H(v_i^{n+1,(0)}) + D(p^{n+1}) + S_v^{n+1} \quad (34) $$

Equations (32) and (33) are then used to produce successive correction, by using the iterative relations:

$$ \left\{ \begin{aligned} A_p v_{i,P}^{n+1,(q)} &= H(v_i^{n+1,(q)}) + D(p^{n+1,(q)}) + S_v^{n+1} \\ A_p p_p^{n+1,(q)} &= G(p^{n+1,(q)}) + S_p(v_i^{n+1,(q)}) \end{aligned} \right. \quad (35) $$

The iterative scheme is stopped when a convergence criterion is satisfied or when a maximum number of inner iterations is reached.

The fluid problem is solved with the numerical principles expose before, using the commercial CFD code Star-CD [28]. An elementary validation of the moving mesh technique and fluid force calculation is performed in the elementary case of an oscillating cylinder in a confined annular space. The FRITZ model [10] gives the added mass coefficient for small amplitude motions in the case of a perfect fluid. Figure 6 gives a comparison of the analytical and computed added mass coefficient for several confinement ratios $\alpha$.

<figure><img src="image_8.png" /><figcaption>Figure 6. Elementary validation of the moving mesh and pressure force calculation with the Star-CD CFD code</figcaption></figure>

## 3. COUPLED PROBLEM

The fluid and structure problem are solved with the numerical principles presented in the previous paragraphs. The numerical methods for fluid are available as such in the CFD code STAR-CD, which will be used for the present study; the structural code will be implemented in the STAR-CD code as FOTRAN subroutines, together with coupling procedures.

The numerical solving strategy of the STAR-CD code does not allow numerical exchanges within a time iteration with PISO

<footnote n="5">⁵ A numerical study currently in progress focuses on the calculation of added mass and added damping coefficient in the case of a viscous fluid. Preliminary results show a good agreement with the analytical model from CHEN and comparison with other numerical results [3] will be carried out.</footnote>
algorithm: coupling with a structure code will thus be possible only with an explicit technique as detailed below.

## 3.1. Coupling in time

As mentioned before, the time coupling strategy will be based on a staggered, explicit coupling procedure, which is represented by Fig. 7.

At time step $t_n$, structure displacement $u_n$ is known. From previous fluid geometry $\Omega_n$, it is possible to calculate the new fluid geometry $\Omega_n$ and solve the space law conservation equation to obtain the grid velocity field $v_n$ (step [1] in the coupling procedure). The calculations proceed with the numerical resolution of the fluid problem on fluid domain $\Omega_n$, requiring inner iterations for the PISO algorithm (step [2]). Fluid forces on the fluid/structure $\phi_n$ interface are then deduced and transferred on the structure problem (step [3]). The advance in time for the structural problem is then carried out, with inner iterations for the implicit scheme (step [4]).

<figure><img src="image_9.png" /><figcaption>Figure 7. FE structure code and FV fluid code explicit time coupling</figcaption></figure>

From the structure point of view, the implicit approach requires the knowledge of external forces at time step $t_{n+1}$, whereas the explicit coupling procedure gives fluid forces from $t_n$ to $t_{n+1}$; a linear extrapolation is used to a better prediction of the fluid forces; the external forces used in the structure advance are $\tilde{\phi}_n = 2\phi_n - \phi_{n-1}$.

For small time step size, or small mesh displacement, the computed fluid forces show high frequency numerical oscillations. This can induce numerical instabilities for the structure part. This problem is addressed by introducing a filtering process within the coupled procedure<fnref n="6" />. A first order filter is used. It is described by its transfer function:

$$H(p) = \frac{1}{\tau \cdot p + 1} \qquad (36)$$

where $p$ is the LAPLACE operator, and $\omega_e = 2\pi / \tau$ is the filtering pulsation. Using a bilinear transform of the LAPLACE variable $p$, that is $p = \frac{2}{\delta t} \cdot \frac{z-1}{z+1}$, and using the z-transform [24], one deduces the following discrete relation:

$$\phi_n^* = \frac{\delta t}{2\tau + \delta t} \cdot (\phi_n + \phi_{n-1} - \phi_{n-1}^*) + \frac{2\tau}{2\tau + \delta t} \cdot \phi_{n-1}^* \qquad (37)$$

where $\phi$ and $\phi^*$ are the non-filtered and filtered discrete values, $\delta t$ is the time step size<fnref n="7" />.

<figure><img src="image_10.png" /><figcaption>Figure 8. Discrete numeric filter</figcaption></figure>

## 3.2. Coupling in space

The space coupling procedure handles numerical exchanges between the fluid and structure problems. Structure displacement are transferred from the structure mesh to the fluid mesh, and fluid forces are projected from fluid boundaries to structure mesh, as sketched by Fig. 9.

<figure><img src="image_11.png" /><figcaption>Figure 9. FE structure code and FV fluid code space coupling</figcaption></figure>

As fluid and structure meshes can be very different (a few finite element are enough to describe the structure dynamic, an important

<footnote n="6">⁶ From a numerical point of view, the existence of numerical oscillations in the computed fluid forces can also induce numerical structure instabilities, even in implicit coupling techniques. Many authors use in this case a blending factor approach. The computed fluid forces are corrected with previous</footnote>

<footnote>forces, for instance as $\tilde{\phi}_n = \theta\phi_n + (1 - \theta)\phi_{n-1}$ with $0 < \theta < 1$ [17]. This correction is equivalent to the filtering/extapolation procedure used in the present study.</footnote>

<footnote n="7">⁷ In the 2D case, this numerical treatment proved to be efficient and precise, even when fluid forces are dominant over structure inner stiffness forces [21]. The major drawback of this technique lies on its empirical approach: several calculation with various filtering frequencies have to be performed in order to demonstrate that numerical results are not affected (up to a certain point) by this numerical approach.</footnote>
number of fluid finite volume are required to solve the fluid problem, with local refinement to describe boundary layer), the space coupling procedure uses interpolation techniques in the displacement force exchanges between the two sub-problems [2].

Since the fluid problem is re-meshed at each time step to take into account structure deformation, the fluid finite volume can undergo great shape deformation. Many re-meshing coupling techniques can be developed to preserve the fluid mesh quality [22]. As the fluid geometry is rather simple in the present case, it is possible to develop a re-meshing procedure that produces little deformation of the fluid cells, based on a purely geometrical approach of the problem.

Figure 10 shows a 2D half slice of the fluid mesh, in the initial state and a deformed state; fluid cells length ratios are preserved even for large amplitude inner cylinder motion<fnref n="8" />.

<figure><img src="image_12.png" /><figcaption>Figure 10. Initial and deformed mesh</figcaption></figure>

## 4. VALIDATION OF THE COUPLING TECHNIQUE

The numerical coupling procedure is validated in the case of a perfect incompressible fluid. Under these assumptions, a linearisation of the NAVIER-STOKES equations leads to the following formulation of the fluid problem, in the non deforming relative frame [11]:

$$ 
\begin{align} 
\Delta p &= 0 \nonumber \\ 
\frac{\partial p}{\partial z(z_0)} &= \frac{\partial p}{\partial z(z_h)} = 0 \tag{38} \\ 
\frac{\partial p}{\partial r(r_h)} &= -\rho_F \gamma \cos \theta \nonumber 
\end{align} 
$$

Coupling conditions with the structure problem are written as:

$$ 
\begin{align} 
\varphi_h(z) &= -\int_0^{\pi} p(R, \theta, z) \cos \theta R d\theta \nonumber \\ 
\frac{\partial p}{\partial \bar{c}^r(r_h)} &= -\rho_F \left( \frac{\partial^2 v}{\partial \bar{c}^2} + \gamma \right) \cos \theta \tag{39} 
\end{align} 
$$

The fluid problem has an analytical solution that can be written in terms of pressure, from which fluid force on the elastic beam is deduced as written as:

$$ 
\begin{align} 
\varphi_h(z) &= -\pi R \sum_{i=0}^{\infty} P_i(R) \cos(q_i z) \times \nonumber \\ 
&\quad \int_{0}^{\pi} \frac{\partial^2 v}{\partial t^2}(\zeta, t) \cos(q_i \zeta) d\zeta + \rho_F \pi R^2 \gamma \tag{40} 
\end{align} 
$$

This leads to the calculation of an added mass matrix $M_d$ in the left hand side of structure equation of motion and a displaced fluid force vector $M_d \Gamma$ on the right hand side:

$$ (M + M_h) \ddot{Y} + K(X) = -(M - M_d) \Gamma \quad (41) $$

The structure problem taking into account fluid effect is then solved with numerical techniques presented in §.1. The geometrical and physical characteristics of the coupled problem are $R = 0.1$ m, $R^* = 0.2$ m, $L = 1$ m, $e = 0.02$ m, $\rho_S = 8000$ kg/m³, $E = 6.04 \times 10^8$ Pa. The whole system is subjected to a sine wave acceleration:

$$ \gamma(t) = \begin{cases} \gamma \times \sin\left(2\pi \frac{t}{\tau}\right) & \text{if } 0 \le t \le \tau \\ 0 & \text{if } \tau < t \end{cases} \quad (42) $$

Figure 11 shows the analytical solution and CFD calculation with explicit FE/FV coupling in terms of beam free end displacement. The parameters of the imposed acceleration is are $\gamma = 50$ m/s², $\tau = 100$ ms. The numerical results show a good agreement between the two methods, which validates the explicit coupling procedure.

<figure><img src="image_13.png" /><figcaption>Figure 11. Non linear beam response: free end displacement $v(L,t)$ computed with the analytical solution and CFD calculation with non linear FE/FV coupling technique</figcaption></figure>

Figure 12 gives a comparison of two coupled calculation with the previous acceleration parameters: the elastic beam is coupled with a light and heavy fluid. The numerical simulations illustrates the added mass and displaced mass effect: the period of oscillation is increased, the amplitude oscillation decreased in the coupling with the heavy fluid, compared to the coupling with the light fluid.

<footnote n="8">⁸ For very large displacement, i.e. when the inner cylinder almost comes into contact with the outer one, the proposed re-meshing technique will fail (other ones too !...) to preserve mesh quality. In this case, a numerical technique based on birth/death fluid cells approach would be useful. Such a method is available in the STAR-CD code [28] and is currently investigated in the present coupled problem.</footnote>
<figure><img src="image_14.png" /><figcaption>Figure 12. Non linear beam response free end displacement v(L,t) computed with CFD calculation. Elastic beam coupled with water ($\rho=1000 \text{ kg/m}^3$) and with air ($\rho=1 \text{ kg/m}^3$)</figcaption></figure>

## CONCLUSION

A general explicit coupling between a structure finite element code and a fluid finite volume code has been presented and applied to the numerical simulation of a simple non linear fluid-structure coupled problem.

Various aspect of the coupling procedure have been detailed: numerical methods for fluid and structure sub-problems. coupling practices, including explicit staggered time coupling, fluid re-meshing strategy, data exchanges with numerical space interpolation and time filtering techniques.

The procedure is validated in the simple case of a perfect inviscid fluid with a numerical solution of the coupled problem. The numerical fluid/structure simulation with a more general fluid model can be considered. Current numerical studies lying on the exposed principles are performed for various cases (viscous fluid, with two fluid phases, etc...).

The present study gives a reference test case for a full scale fluid-structure model. Coupling commercial FE code (such as ANSYS) and FV code (such as STAR-CD) with general coupling code (such as MpCCI) for industrial applications can now be considered.

The industrial application under consideration are related to coupled problems where the fluid-structure system, initially at rest, is excited by imposed displacement or acceleration. Fluid-structure dynamic are governed by added mass and damping effect, for which a numerical explicit resolution gives good results, without too a heavy modeling and computational cost. For coupled problems in which fluid flow induces continuous energy exchanges from fluid to structure leading to structure instabilities, explicit techniques might not be effective. Further studies will tend to investigate implicit coupling methods with appropriate numerical tools and numerical strategies and perform industrial calculation in such cases.

## REFERENCES

[1] K.J. BATHE. *Finite Element Procedures in Engineering Analysis*. Prentice-Hall, Englewood Cliffs, 1982.

[2] A. BECKERT. Coupling Fluid (CFD) and Structural (FE) models Using Finite Interpolation Elements. *Aerospace Science and technology*, 4 (1), 13-22, 2000.

[3] Z. BENDJEDOU, E. LONGATE, A. ADOBES, M. SOULI. Steam Generator Tube Vibrations Experimental determination versus ALE Computation of Fluidelastic Forces. Pressure Vessel and Piping, Cleveland, 20-24 July 2003.

[4] F. BEROT, B. PESEUX. Vibro-Acoustic Behavior of Submersed Cylindrical Shells: Analytical Formulation and Numerical model. *Journal of Fluids & Structures*, 12, 959-1003, 1998.

[5] F. CASADEI, J.P. HALLEUX, A. SALA, F. CHILLE. Transient Fluid Structure Interaction Algorithm for Large Industrial Applications. Computer Methods in Applied Mechanical Engineering. 190, 3081-3110, 2001.

[6] I. DEMIRDZIC, M. PERIC. Space Conservation law in Finite Volume Calculations of Fluid Flow. *International Journal for Numerical Methods in Fluids*, 8, 1037-1050, 1988.

[7] M.A. DOKAINISH, K. SUBBARAJ. A Survey of Direct Time-Integration Methods in Computational Dynamics-I: Explicit Methods. *Computer & Structures*, 32, 1371-1386, 1989.

[8] C. FARHAT, M. LESOINNE, N. MAMAN. Mixed Explicit/Implicit Time Integration of Coupled Aeroelastic Problems. *International Journal of Numerical Methods in Fluids*, 21, 807-835, 1995.

[9] J.H. FERZIGER, M. PERIC. *Computational Methods for Fluid Dynamics*. Springer-Verlag, 1999.

[10] R.J. FRITZ. The Effect of Liquids on the Dynamic Motion of Immersed Solids. *Journal of Engineering for Industry*, 167-173, 1972.

[11] R.J. GIBERT. Vibration des structures. Interaction avec les fluides. Sources d'excitation aléatoires. Collection de la Direction des Etudes et Recherches d'Electricité de France, vol. 69, Eyrolles, 1986.

[12] J. HINOY. YOSHIMURA, N. ANANTHANARAYANA. Vibration Analysis of Non-Linear Beams Subjected to a Moving Load Using the Finite Element Method. *Journal of Sound and Vibration*, 100, 477-491, 1985.

[13] R.I. ISSA. Solution of The Implicitly Discretised Fluid Flow Equation by Operator Splitting. *Journal of Computational Physics*, 62, 40-65, 1985.

[14] E. DE LANGRE. *Fluides et Solides*. Les éditions de l'Ecole Polytechnique, 2002.

[15] J. MAKERLE. Fluid-Structure Interaction Problems, Finite Element Approach and Boundary Elements Approaches. A Bibliography. *Finite Elements in Analysis and Design*, 31, 231-240, 1999.

[16] J.J.P. MORAND, R. OHAYON. *Fluid Structure Interaction*. Wilcy & Sons, 1995.

[17] J.M.T. PENROSE, C.J. STAPLES. Implicit Fluid-Structure Coupling for Simulation of Cardiovascular Problems. *International Journal for Numerical Methods in Fluids*, 1, 1-13, 2001.

[19] J. SARRATE, A. HUERTA, J. DONEA. Arbitrary Lagrangian Eulerian Formulation for Fluid Multi-Rigid Bodies Interaction Problems. *Computation Mechanics, New Trend and Applications*, 1998.

[20] M. SCHÄFER. Coupled Fluid-Solid Problems Survey on Numerical Approaches and Applications. Pressure Vessel and Piping. Cleveland, 20-24 July 2003.

[21] J.F. SIGRIST et al. Numerical Simulation of Fluid Structure Problems by Coupling Fluid Finite Volume and Structure Finite Element or Modal Approach. Flow Induced Vibrations, Paris, 5-9 July 2004.

[22] R.M.C. SO, Y. LRU, Y.G. LAI. Mesh Shape Preservation for Flow-Induced Vibrations Problems. International Mechanical Engineering Congress and Exposition. New orleans, 17-22 November 2002.

[23] K. SUBBARAJ, M.A. DOKAINISH. A Survey of Direct Time-Integration Methods in Computational Dynamics-II: Implicit Methods. *Computer & Structures*, 32, 1387-1401, 1989.

[24] Y. THOMAS. *Signaux et systèmes linéaires*. Masson, 1992.

[25] Z.U.A. WARSI. Conservation form of the Navier Stokes Equations in General Non Steady Coordinates. *American Institute of Aeronautics and Astronautics Journal*, 19, 240-242, 1980.

[26] K. WASHIZU. Variational Methods in Elasticity and Plasticity. Pergamon, 1982.

[27] T. WINTERGERSTE. Fluid-Structure Interaction with STAR-CD and PERMAS using MpCCI. *BENCHmark* (NAFEMS Publication), April 2002.

[28] STAR-CD Methodology-User Guide. CD Adapco Group, 2003.