# Numerical methods for RANSE simulations of a self-propelled fish-like body

Alban Leroyer, Michel Visonneau

► To cite this version:

Alban Leroyer, Michel Visonneau. Numerical methods for RANSE simulations of a self-propelled fish-like body. Journal of Fluids and Structures, 2005, 20 (7), pp.975-991. 10.1016/j.jfluidstructs.2005.05.007. hal-00699458

HAL Id: hal-00699458

https://hal.science/hal-00699458

Submitted on 1 Mar 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Numerical methods for RANSE simulations of a
self-propelled fish-like body

A. Leroyer, M. Visonneau

Laboratoire de Mécanique des Fluides, CNRS-UMR 6598, Ecole Centrale de Nantes, B.P. 92101,
44321 Nantes Cedex 3, France

This paper deals with numerical simulations on a self-propelled fish-like body. Firstly, the methods used to achieve such simulations based on the solution of the Reynolds-Averaged-Navier-Stokes Equations are described: resolution of the Newton's law (parameterization, use of a description of the rotation based on a quaternion), regridding strategies linked to body motion and deformation, fluid-motion coupling methods. In the second part, realistic applications are shown. The interaction between the fluid and the body is fully solved since the evolution of the body's position is not imposed but computed. Only the imposed deformation of the shape and its interaction with the fluid produce forces on the body and, consequently, motion.

Keywords: Biomimesis; Biohydrodynamics; Navier-Stokes solver; Newton's laws; Fluid-structure interaction

1. Introduction

During the last two decades, great advances have been achieved in turbulence modelling, mesh generation and linear
system resolution. Therefore, flow solvers have become more general, robust and accurate owing to the large amount of
research work and to the growth of storage capacity and computers power. Now, the maturity of algorithms makes
possible integration and coupling of new physical features to deal with more and more realistic applications: multi-
phase flows, cavitation, fluid-structure interaction, etc. Flows around moving bodies is one of these interdisciplinary
fields of interest (Triantafyllou et al., 2000; Liu and Kawachi, 1999; Chen and Doi, 2002). Recently, this feature has
been implemented into the ISIS flow solver, developed by the CFD Department of the Fluid Mechanics Laboratory.
Bodies can be rigid or flexible with a prescribed law of deformation. Whereas the former can lead to application related
to slamming and seakeeping, these latter configurations can be very useful to study underwater locomotion in order to
understand the remarkable solutions developed by living species over million of years. This paper is devoted to one of
them: the self-propelled motion of fish-like bodies.

E-mail address: Michel.Visonneau@ec-nantes.fr
## 2. Flow solver

The ISIS flow solver, developed by EMN (Equipe Modélisation Numérique, i.e., the CFD Department of the Fluid Mechanics Laboratory), uses the incompressible Unsteady Reynolds-Averaged Navier-Stokes equations (URANS). The solver is based on the finite volume method to build the spatial discretization of the transport equations. The face-based method is generalized to 2-D, rotationally symmetric, or 3-D unstructured meshes for which nonoverlapping control volumes are bounded by an arbitrary number of constitutive faces. The velocity field is obtained from the momentum conservation equations and the pressure field is extracted from the mass conservation constraint, or continuity equation, transformed into a pressure equation. In the case of turbulent flows, additional transport equations for modelled variables are solved in a form similar to the momentum equations and they can be discretized and solved using the same principles.

### 2.1. Conservation equations

For an incompressible flow of viscous fluid under isothermal conditions, mass, momentum and volume fraction conservation equations can be written as (using the generalized form of Gauss' theorem):

$$ \frac{\partial}{\partial t} \int_V \rho dV + \int_S \rho (\vec{U} - \vec{U}_d) \cdot \vec{n} dS = 0, \qquad (1a) $$

$$ \frac{\partial}{\partial t} \int_V \rho U_i dV + \int_S \rho U_i (\vec{U} - \vec{U}_d) \cdot \vec{n} dS = \int_S (\tau_{ij} I_j - p I_i) \cdot \vec{n} dS + \int_V \rho g_i dV, \qquad (1b) $$

where $V$ is the domain of interest, or control volume, bounded by the closed surface $S$ moving at the velocity $\vec{U}_d$ with a unit normal vector $\vec{n}$ directed outward. $\vec{U}$ and $p$ represent respectively the velocity and pressure fields. Density is defined by $\rho$, $\tau_{ij}$ and $g_i$ are the components of the viscous stress tensor and the gravity, whereas $I_j$ is a vector whose components vanished, except for the component $j$ which is equal to unity.

### 2.2. Turbulence modelling

Several turbulence closures are included in the flow solver, ranging from linear eddy-viscosity-based models to full second-order closures (Duvigneau et al., 2003). For these studies, the one-equation model of turbulence (Spalart and Allmaras, 1992) and different other closures based on $k-\omega$ formulation have been used.

### 2.3. Numerical framework

#### 2.3.1. Spatial discretization

All the flow variables are stored at geometric centers of the arbitrary-shaped cells. Surface and volume integrals are evaluated according to second-order accurate approximations by using the values of integrand that prevail at the center of the face $f$, or cell $C$, and neighbor cells. The various fluxes appearing in the discretized equations are built using centered and/or upwind schemes. For example, the convective fluxes are obtained by two kinds of upwind schemes. A first scheme available in the flow solver, (HD) for hybrid differencing, is a combination of upwind (UD) and centered (CD) schemes. Contrary to a practical approach (Demirdić and Muzaferija, 1995; Ferziger and Perić, 2002) where CD/UD blending is fixed with a global blending factor for all faces of the mesh, the HD scheme results from a local blending factor based on the signed Peclet number at the face. An other upwind scheme which is implemented in ISIS, is the Gamma Differencing Scheme (GDS) (Jasak, 1996). Through a normalized variable diagram (NVD) analysis (Leonard, 1988), this scheme enforces local monotonicity and convection boundedness criterium (CBC) (Gaskell and Lau, 1988). A pressure equation is obtained in the spirit of the Rhie and Chow procedure (Rhie and Chow, 1983) and momentum and pressure equations are solved in an segregated way like in the well-known SIMPLE coupling procedure.

#### 2.3.2. Temporal discretization

The temporal discretization is based on a three-step upwind discretization ensuring a second-order accuracy. The convection and diffusion are treated implicitly.
### 3. Resolution of the Newton's laws

The aim is to calculate the temporal evolution of the kinematic characteristics of a rigid or flexible body submitted to forces acting on it (hydrodynamic forces, gravity, etc.).

The first step consists in the definition of an initial inertial frame of reference $\mathcal{R}_0 = (O_0, \beta_0)$, linked to the physical space, assimilated to a Galilean referential. The origin $O_0$ of this referential is fixed to the mass center of the body at the initial time. The desired shape of the body is then imposed on $\mathcal{R}_0$. A rigid motion moving $\mathcal{R}_0$ to a new referential $\mathcal{R}_1$ is then carried out to place the body in space (see Fig. 1). Characteristics of the body motion are given by this transformation $\mathcal{R}_0 \to \mathcal{R}_1$ (position and orientation).

Considering the material system $S$ defined by the domain $\mathcal{D}_t$, the Newton's second law leads to the two Eqs. (2) and (3), where $G(t)$ is the center of gravity of the system, $\vec{R}$ the total force acting on $S$ and $\vec{M}_{O_1}$ the torque.

Conservation of the linear momentum:

$$ \overline{\mathcal{A}(S/\mathcal{R}_0)} = \left. \frac{d}{dt} \right|_{\mathcal{R}_0} \iiint_{\mathcal{D}_t} \rho \overline{V(M, S/\mathcal{R}_0)} dv = \overline{R}. \qquad (2) $$

Conservation of the angular momentum at $O_1$:

$$ \delta(O_1, S/\mathcal{R}_0) = \left. \frac{d}{dt} \right|_{\mathcal{R}_0} \iiint_{\mathcal{D}_t} \overline{O_1 M} \times \rho \overline{V(M, S/\mathcal{R}_0)} dv + \overline{V(O_1/\mathcal{R}_0)} \times M_S \overline{V(G(t)/\mathcal{R}_0)} = \overline{M}_{O_1}. \qquad (3) $$

In the framework of this general 3-D resolution, classical description of the body orientation by three successive rotations (like Euler Angles), written as $(\psi, \theta, \phi)$, is unsuitable because of singular configurations for which the triplet $(\psi, \theta, \phi)$ is not unique. The use of a quaternion eliminates these problems. This technique has already been used previously by McDonald and Whitfield (1997). Quaternion space $\mathbb{H}$ is a 4-D algebra spanned by a real axis e and three orthogonal imaginary axis denoted by $i, j, k$. It can be roughly considered as a spatial extension of the plane representation of complex numbers. One can show that for any rotation around the axis directed by the unit vector $\vec{u}$ with an angle $\theta$ is associated only one quaternion $q = \cos(\theta/2)e + u \sin(\theta/2)$ on the $\mathbb{H}$ basis $(e, i, j, k)$. This quaternion (and its time derivative) can be related with the instantaneous rotation vector $\vec{\Omega}$. This relation makes possible to calculate $q$ knowing $\vec{\Omega}$. Therefore, using Eqs. (2) and (3) on the studied body, the deformation-motion decomposition leads to the following coupled system:

$$ \begin{bmatrix} 1 & 0 & 0 & 0 & z_G & -y_G \\ 0 & 1 & 0 & -z_G & 0 & x_G \\ 0 & 0 & 1 & y_G & -x_G & 0 \\ 0 & -z_G & y_G & \frac{A1}{M_S} & \frac{-F1}{M_S} & \frac{-E1}{M_S} \\ z_G & 0 & -x_G & \frac{-F1}{M_S} & \frac{B1}{M_S} & \frac{-D1}{M_S} \\ -y_G & x_G & 0 & \frac{-E1}{M_S} & \frac{-D1}{M_S} & \frac{C1}{M_S} \end{bmatrix} \frac{d}{dt} \begin{bmatrix} \dot{x} \\ \dot{y} \\ \dot{z} \\ \alpha \\ \beta \\ \gamma \end{bmatrix} = \frac{1}{M_S} \begin{bmatrix} \vec{R} \cdot \vec{x}_0 - \overrightarrow{S_{\text{Res}} \cdot \vec{x}_0} \\ \vec{R} \cdot \vec{y}_0 - \overrightarrow{S_{\text{Res}} \cdot \vec{y}_0} \\ \vec{R} \cdot \vec{z}_0 - \overrightarrow{S_{\text{Res}} \cdot \vec{z}_0} \\ \overrightarrow{M_{O_1} \cdot \vec{x}_0} - \overrightarrow{S_I \cdot \vec{x}_0} - \overrightarrow{S_{\text{Mom}} \cdot \vec{x}_0} \\ \overrightarrow{M_{O_1} \cdot \vec{y}_0} - \overrightarrow{S_I \cdot \vec{y}_0} - \overrightarrow{S_{\text{Mom}} \cdot \vec{y}_0} \\ \overrightarrow{M_{O_1} \cdot \vec{z}_0} - \overrightarrow{S_I \cdot \vec{z}_0} - \overrightarrow{S_{\text{Mom}} \cdot \vec{z}_0} \end{bmatrix}, \qquad (4) $$

$$ \frac{d}{dt} \begin{bmatrix} x \\ y \\ z \end{bmatrix} = \begin{bmatrix} \dot{x} \\ \dot{y} \\ \dot{z} \end{bmatrix}, \qquad (5) $$

$$ \frac{d}{dt} \begin{bmatrix} q_0 \\ q_1 \\ q_2 \\ q_3 \end{bmatrix} = \frac{1}{2} \begin{bmatrix} 0 & -\alpha & -\beta & -\gamma \\ \alpha & 0 & -\gamma & \beta \\ \beta & \gamma & 0 & -\alpha \\ \gamma & -\beta & \alpha & 0 \end{bmatrix} \begin{bmatrix} q_0 \\ q_1 \\ q_2 \\ q_3 \end{bmatrix}, \qquad (6) $$

$$ q_0^2 + q_1^2 + q_2^2 + q_3^2 = 1, \qquad (7) $$

where $M_S$ means the mass of the body, $(A1, B1, C1, D1, E1, F1)$ are the elements of the inertia matrix of the body (see (4)), $(x, y, z)$ the coordinates of the vector $\overrightarrow{O_0O_1} = \vec{X}$ on the basis $\beta_0$, $(\dot{x}, \dot{y}, \dot{z})$ the coordinates of the vector $\overrightarrow{V(O_1/\mathcal{R}_0)}$ on
<figure><img src="image_2.png" /><figcaption>Fig. 1. Deformation-motion splitting.</figcaption></figure>

the basis $\beta_0$, $(\alpha, \beta, \gamma)$ the coordinates of the vector $\overrightarrow{\Omega_0^I}$ on the basis $\beta_0$, $(q_0, q_1, q_2, q_3)$ the coordinates of $\mathbf{q}$ on the basis associated to $\beta_0$, $\{\overrightarrow{R}, \overrightarrow{M_{O_1}}\}$ the screw of external forces acting on the body, $\overrightarrow{S_I}$ the inertial source term, $(\overrightarrow{S_{Res}}, \overrightarrow{S_{Mom}})$ the source term induced by the possible deformation of the body.

In the foregoing, the expression for $\overrightarrow{S_I}$ is given by

$$ 
\begin{aligned}
\overrightarrow{S_I} &= \iiint_{\mathcal{G}_I} \rho(\overrightarrow{O_1M} \cdot \overrightarrow{\Omega_0^I})(\overrightarrow{O_1M} \times \overrightarrow{\Omega_0^I}) dv, \\
\overrightarrow{S_I} \cdot \overrightarrow{\chi_0} &= \alpha\gamma F1 - \alpha\beta E1 + (\gamma^2 - \beta^2) D1 + \beta\gamma (C1 - B1), \\
\overrightarrow{S_I} \cdot \overrightarrow{\chi_0} &= \alpha\beta D1 - \beta\gamma F1 + (\alpha^2 - \gamma^2) E1 + \alpha\gamma (A1 - C1), \\
\overrightarrow{S_I} \cdot \overrightarrow{\zeta_0} &= \beta\gamma E1 - \alpha\gamma D1 + (\beta^2 - \alpha^2) F1 + \alpha\beta (B1 - A1);
\end{aligned}
\quad (8) $$

and the expressions $\overrightarrow{S_{Res}}$ and $\overrightarrow{S_{Mom}}$ are defined by

$$ \overrightarrow{\frac{S_{Res}}{M_S}} = \overrightarrow{\Gamma(G, S/\mathcal{R}_1)} + 2 \overrightarrow{\Omega_0^I} \times \overrightarrow{V(G, S/\mathcal{R}_1)} + \overrightarrow{\Omega_0^I} \times (\overrightarrow{\Omega_0^I} \times \overrightarrow{O_1G}); \quad (9) $$

$$ \overrightarrow{S_{Mom}} = \iiint_{\mathcal{G}_I} \rho \overrightarrow{O_1M} \times \overrightarrow{\Gamma(M, S/\mathcal{R}_1)} dv + 2 \iiint_{\mathcal{G}_I} \rho(\overrightarrow{V(M, S/\mathcal{R}_1)} \cdot \overrightarrow{O_1M}) \overrightarrow{\Omega_0^I} dv \\ - 2 \iiint_{\mathcal{G}_I} \rho(\overrightarrow{O_1M} \cdot \overrightarrow{\Omega_0^I}) \overrightarrow{V(M, S/\mathcal{R}_1)} dv. \quad (10) $$

Eq. (7) is redundant because a quaternion describing rotation has always a unit length. However, the numerical resolution of Eq. (6) does not maintain exactly this length and a renormalization is then necessary. An internal mesh of flexible bodies is needed to compute additional dynamic terms induced by the imposed deformation.

## 4. Regridding strategies

To implement bodies motion in a flow solver, the mesh must be adapted to the new position of bodies in time. In order to keep an appropriate grid, three complementary methods have been integrated:

(i) spring analogy regridding procedure,

(ii) rigid transformation of the mesh,

(iii) analytical weighted regridding approach.

These different regridding strategies can be mixed which provides a flexibility very useful when one wants to combine deformation and rigid body motion modes.
## 4.1. Spring analogy technique

With this method, the mesh is represented by a mechanical structure composed by fictitious lineal and torsional springs (Farhat et al., 1998). The new mesh is obtained by resolving a pseudo-structural system after imposing the new position of nodes belonging to bodies and boundaries. This technique is very convenient because of its capacity to deal with any deformation. The quality of the final mesh remains quite good even for relatively large displacements. Furthermore, the connectivity is conserved. However, the computational cost is quite expensive because of the resolution of three additional elliptic equations.

## 4.2. Rigid transformation

Here, all the nodes are moved with the solid motion of the body, coming from the resolution of the Newton's laws or from an imposed law. This approach is operational only for a simulation with one body in an infinite fluid domain, but tolerates movement of arbitrary amplitudes.

## 4.3. Analytical weighted regridding

The aim of the technique is to set a direct method without solving any system in the case of peculiar motion (rigid or deriving from beam kinematics).

It is a kind of extrapolation of the rigid transformation. But, in this case, the displacement of each node is weighted by a factor $k_p$ varying between one and zero according to the distance of the body (one for the nodes belonging to the body, zero for the nodes of other boundaries). In the case of multi-body simulations, one coefficient is attached to each node for each motion. The values of this factor are calculated at the beginning of the simulation by solving a Laplacian operator.

As far as rigid bodies are concerned, the weighting is directly applied to the parameters of the rigid transformation $\mathcal{R}_0 \to \mathcal{R}_1$. Concerning beam-like bodies, the procedure is slightly more complex: a “virtual” beam crossing the whole domain and overlapping the neutral line of the studied beam-like body is defined. In the undeformed state, each node is associated to the cross-section of this virtual beam embedded in it. With the Kirchoff assumptions—the cross-sections remain plane, unstretched and normal to the neutral line—the new mesh is obtained by weighting with the factor $k_p$ the rotation of each cross section (see Figs. 2 and 3). In fact, this new position called “virtual position” and written $\vec{X}_v$ is relaxed with $\vec{X}_r$ which represents the initial position $\vec{X}_{ref}$ possibly turned by an adapted rotation in order to limit the deformation. Therefore, the final position of the node ($\vec{X}_{def}$) is given by

$$ \vec{X}_{def} = k_X \vec{X}_v + (1 - k_X) \vec{X}_r \quad \text{with} \quad k_X = \frac{\sin\left(\frac{\pi}{2}(2k_p - 1)\right) + 1}{2}. \qquad (11) $$

An analytical smoothing is also carried out to improve the quality of this body-fitted mesh. An extension of this approach for shell-like bodies is also feasible but has not yet been implemented.

## 4.4. Choice of the method

For simulations dealing with an isolated body in an infinite domain, the displacement of the mesh by a rigid transformation seems to be the most appropriate method. This technique is not expensive in terms of CPU time and

<figure><img src="image_3.png" /><figcaption>Fig. 2. Definition of the “virtual” beam.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Analytical weighting for beam-like body.</figcaption></figure>

enables motion of unlimited amplitudes in the physical space. When, for instance, the interaction between the motion of the ship and the free-surface flow computed with free-surface capturing strategies is considered, it is mandatory to keep as much as possible unchanged the position of the refined grid domain describing the interface between air and water, to avoid any numerical smearing of this interface. Therefore, in that case, a regridding strategy based on a spring analogy or solid-motion weighting can be used for the three degrees of rotation (yaw, roll and pitch) and one degree of vertical translation (heave) in the vicinity of the body and a solid body motion should be preferred to displace the body according to the two remaining horizontal degrees of freedom. Most of the time, a solid-motion weighted strategy can be used when the amplitudes of these degrees of freedom are moderate, and this approach should be preferred because the cost in terms of CPU time of a spring analogy regridding technique which requires the solution of three elliptical equations is very high.

In the case of a flexible body, the use of spring analogy technique is compulsory if the imposed deformation is not based on beam kinematics. This is not yet so damaging since the imposed deformation needs to be applied only at each time step. As a matter of fact, only the position in space is modified during the nonlinear iterations. The latter depending on a rigid motion, the two other methods which have insignificant CPU time can be used to update the body position in the nonlinear loop. For beam-like body, the analytical weighted technique can favourably replace the spring analogy concerning the deformation part.

# 5. Fluid-motion coupling

In the case of prescribed motions, the body is displaced at every time step. After recomputing the new mesh, the flow is solved. Therefore, there is no real coupling when motion is imposed, due to the lack of fluid feedback on the body position. The problem is more complex when motion is solved with Newton's law. Actually, body kinematics are linked to the flow at the same time step by forces acting on them. In the first method, called weakly coupled method, hydrodynamic forces and moments provided by the RANS equations are calculated only at the end of every time step. They are used to obtain the new position of the body. With this approach, instabilities may appear as soon as the body density is close to the density of the surrounding fluid, situations which are frequent in hydrodynamics for floating bodies, for instance. This justifies a stronger coupling.

## 5.1. Nonlinear coupling

In order to stabilize coupling, the estimated flow during the nonlinear iteration is used to calculate forces and moments acting on bodies and then new body position. This process is iterated until nonlinear RANSE and Newton's laws residuals are reduced by an appropriate factor (here 10³). This is described by the algorithm shown in Fig. 4.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Nonlinear coupling.</figcaption></figure>

Even if this procedure is definitely more stable for the test-cases involving solid bodies, divergent oscillations may remain for flexible bodies. To achieve a stable coupling, a specific treatment of the equations of motion needs to be applied.

## 5.2. Stabilization of the nonlinear coupling

As shown by Söding (2001), the instability comes from the dependence of fluid forces on body acceleration. This force is usually split in a first part which does not depend on the acceleration and a second part linearly linked with it which is called the added mass term:

$$f = \bar{f} - a\gamma \quad \text{with } \gamma: \text{body acceleration.}$$

The estimation of the added mass coefficients would require the solution of one Poisson equation for each degree of freedom at each nonlinear iteration, which would lead to a dramatic increase of the computational cost. Since such an accurate evaluation is not compulsory, this physical property is simply used to modify the resolution of the Newton's laws equations in order to produce a stable coupling.
Look at the 1-D problem. The equation of the motion is reduced to

$$m\gamma = f.$$

The two sides of this equation are increased of an added mass term in which the right-hand side is evaluated at the previous nonlinear iteration. The equation to solve is then

$$(m + a)\gamma^{\text{nlit}} = f + a\gamma^{\text{nlit} - 1}.$$

At convergence ($\gamma^{\text{nlit}} = \gamma^{\text{nlit} - 1}$), the initial equation is recovered. The term $a$ needs not be accurately assessed, since it does not affect the converged solution. It must only be chosen just high enough to ensure stability.

# 6. Applications

The methods described before have been used to simulate flows around an elongated eel-like form. These works are performed in the framework of the CNRS Project named ROBEA. The aim of this interdisciplinary project is to design an efficient flexible underwater robot based on fish-like locomotion. The actual objectives of our flow computations are mainly to build a simplified model linking deformation and motion, which will be used for the real-time command of each vertebra of the robot.

## 6.1. Geometry and mesh

Since the project is only in its first stage, the shape and the length are not entirely designed. Therefore, the actual simulations are performed on a temporary 1 m long geometry (see Fig. 5) whereas the final prototype will be about 2 m long with a slightly different shape. The motion is constrained to remain 2-D with 3 d.o.f.: two translations ($T_x, T_y$) and one rotation ($R_z$). The computational domain is then limited to the upper half-space with a symmetry condition on the plane $z = 0$.

Two meshes have been tested. The first one, named M65e3 contains 65000 cells and is fitted to a wall-function approach for turbulence models (see Fig. 6). The second one, named M12e4, contains 120000 cells which makes possible to use a priori more accurate near wall low-Reynolds number turbulence closures.

## 6.2. Description of the prescribed law of deformation

These analytical laws have to impose the temporal evolution of the fish-body shape in the referential $\mathcal{R}_0$. To simplify the definition, the position and the orientation of the body head remain fixed. This referential is then displaced to fulfill the Newton's law.

### 6.2.1. The "Forward" law

This law aims at representing the periodic evolution of the eel shape during straightforward locomotion (see Fig. 7). In this configuration, the curvature named $c_F$ is specified by Eqs. (12) and (13):

$$c_F(s, t) = A(s) \sin \left[ 2\pi \left( \frac{s}{\lambda} - \frac{t}{T} \right) \right], \qquad (12)$$

$$A(s) = A_2 s^2 + A_1 s + A_0. \qquad (13)$$

<figure><img src="image_6.png" /><figcaption>Fig. 5. Visualization of the eel-body shape.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 6. Visualization of the fluid domain (M65e3 mesh).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Visualization of the “Forward” law.</figcaption></figure>

The main purpose of this law is to reach a periodic state. Therefore, to avoid too long a transient state, bodies can be launched with an initial speed. At the beginning of the simulation, the law of deformation is progressively prescribed by weighting the curvature by a relaxation coefficient evolving smoothly from 0 to 1 in $n_r$ periods (here $n_r = 1$).

The characteristics of the “Forward” law are given by (i) the wavelength $\lambda$, (ii) the period $T$, and (iii) the coefficients of the polynomial $A0, A1$ and $A2$.

### 6.2.2. The “Turn” law

For this maneuvering, the curvature denoted by $c_T$ is imposed:

$$c_T(x, t) = \alpha(t) c_F(x, t) + \frac{1 - \alpha(t)}{R_T} \qquad (14)$$

The turn maneuvering has been divided into three phases, as detailed in the Table 1. During the turn phase, a constant radius of curvature $R_T$ is simply imposed ($\alpha(t) = 0$). Before and after the turn maneuvering, the “Forward” law is prescribed. Then $\alpha(t) = 1$. Fig. 8 shows the evolution of the body shape during the P1 and P3 transient phase of the turn.
<table><caption>Table 1
Subdivision of the “Turn” law</caption><thead><tr><th>Phase</th><th>P1: Start of the turn</th><th>P2: Turn</th><th>P3: End of the turn</th></tr></thead><tbody><tr><th>Timing</th><td>α<sub>T1</sub> T</td><td>α<sub>T2</sub> T</td><td>α<sub>T3</sub> T</td></tr><tr><th>α</th><td>1 → 0</td><td>0 — 0</td><td>0 → 1</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Visualization of P1 and P3 phases of the “Turn2” law.</figcaption></figure>

<table><caption>Table 2
Values of the parameters</caption><thead><tr><th>Law</th><th>λ</th><th>T</th><th>A2</th><th>A1</th><th>A0</th><th>R<sub>T</sub></th><th>α<sub>T1</sub></th><th>α<sub>T2</sub></th><th>α<sub>T3</sub></th></tr></thead><tbody><tr><th>Forward</th><td>1</td><td>1</td><td>2</td><td>0.5</td><td>1</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><th>Turn1</th><td>1</td><td>1</td><td>2</td><td>0.5</td><td>1</td><td>0.5</td><td>0.5</td><td>0.5</td><td>0.5</td></tr><tr><th>Turn2</th><td>1</td><td>1</td><td>2</td><td>0.5</td><td>1</td><td>0.5</td><td>0.5</td><td>1</td><td>0.5</td></tr></tbody></table>

## 6.3. Simulations and results

As suggested by the geometry and the beam-like kinematics imposed to the body, analytical weighting regridding procedure described before is used to build a body-fitted mesh at each nonlinear iteration. The displacement in space of the referential $\mathcal{R}_1$ is carried out by a rigid transformation. A nonlinear coupling with an added-mass term is required to ensure a stable coupling. One “Forward” law and two “Turn” laws have been performed. All the parameters are given in Table 2.

### 6.3.1. Straightforward self-propelled locomotion

Here, the “Forward” law is imposed. Several turbulence closures have been tested on the two meshes M65e3 and M12e4. Furthermore, a laminar simulation has been performed on the M12e4 mesh to evaluate the modifications implied by turbulence. Results given without any comment concern the mesh M12e4 with the $k-\omega$ Wilcox model.

At the beginning of the simulations, the fish is swimming towards $-\mathcal{X}$. An initial speed of $-0.6 \text{ m s}^{-1}$ is imposed to the body. The evolution of the solved motion is described in Fig. 9. As it can be expected, the longitudinal speed along the mean forward axis tends towards a constant mean speed (Fig. 10). One can notice that the transient state induces a small rotation of the body which is conserved after reaching the periodic state. In the actual configuration, this angle is about $1^\circ$ (Fig. 11).

As it can be seen in Fig. 12, the global pressure force resulting both from the deformation of the body and its motion acts favourably for the forward motion, whereas the viscous force tends to slow down the eel. The global external force
<figure><img src="image_10.png" /><figcaption>Fig. 9. Representation of the solved motion.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Evolution of the head longitudinal speed.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Evolution of the head angle.</figcaption></figure>

(i.e. the sum of the two previous ones) tends towards a mean zero value when the body has reached a periodic state, which is not surprising since the mean speed of the eel becomes constant.

Figs. 13 and 14 show the transversal velocity with respect to longitudinal velocity for the head and the tail of the eel. As expected, the amplitude of transversal velocity is larger near the tail. In both cases, one can observe a limit cycle corresponding to the periodic regime.

The value of the mean longitudinal speed is given for some test-cases in Table 3. Influence of the turbulence modelling and, to a lesser extent, of the mesh can be globally analysed. Firstly, the laminar flow leads to a higher mean speed than for all turbulent flows (about 25%). Differences between all turbulence modellings are quite small, the maximum discrepancy for the mean speed being about 5%. The higher differences are observed for the mesh M12e4 adapted to near-wall low-Reynolds number turbulence closures (see Fig. 15), whereas all the tested wall-function turbulence modelling produce quite similar results. For each tested turbulence closure, the near-wall low-Reynolds number simulation (mesh M12e4) always provides a higher mean speed than the wall-function-based closure (mesh M65e3).

The nondimensional Reynolds number and Strouhal number are also given in Table 3 for the periodic state. $Re = UL/v$ is based on the length L of the body, the mean speed U and the kinematic viscosity v. $St = fA/U$ is calculated using f, the frequency of the imposed deformation and A, the total transversal motion of the tail in the
<figure><img src="image_13.png" /><figcaption>Fig. 12. Evolution of forces acting on the body.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Transversal speed versus longitudinal speed of head.</figcaption></figure>

head-fixed reference frame. It should be noticed that, if *A* is calculated using the total transversal motion of the tail with respect to the longitudinal axis, the Strouhal numbers are slightly lower.

### 6.3.2. *Turn maneuvers*

These simulations aim at showing the possibilities of reproducing not only straightforward swimming but also fast-maneuvering. Here are represented two turn maneuvers with different P2 times.

As it can be seen in Fig. 16, the longer is the phase P2, the wider is the rotation of the turn. The speed modulus of the head is shown in Fig. 17. It can be noticed that after the turn, the body reaches again its periodic state in its new direction very quickly. The curvature radius observed during the turn is about that of the body.
<figure><img src="image_15.png" /><figcaption>Fig. 14. Transversal speed versus longitudinal speed of tail.</figcaption></figure>

<table><caption>Table 3
Characteristics of the periodic state</caption><thead><tr><th>Mesh</th><th>Turbulence closure</th><th>Mean longitudinal speed</th><th>Re</th><th>St</th></tr></thead><tbody><tr><td>M12e4</td><td>Laminar</td><td>0.777 m s<sup>-1</sup></td><td>7.77 &times; 10<sup>5</sup></td><td>0.335</td></tr><tr><td>M12e4</td><td>Spalart-Allmaras</td><td>0.620 m s<sup>-1</sup></td><td>6.20 &times; 10<sup>5</sup></td><td>0.420</td></tr><tr><td>M12e4</td><td>Near wall k-&omega;(SST) (Menter)</td><td>0.649 m s<sup>-1</sup></td><td>6.49 &times; 10<sup>5</sup></td><td>0.400</td></tr><tr><td>M12e4</td><td>Near wall k-&omega; (Wilcox)</td><td>0.631 m s<sup>-1</sup></td><td>6.31 &times; 10<sup>5</sup></td><td>0.412</td></tr><tr><td>M12e4</td><td>Near wall k-&omega;(BSL) (Menter)</td><td>0.637 m s<sup>-1</sup></td><td>6.37 &times; 10<sup>5</sup></td><td>0.408</td></tr><tr><td>M65e3</td><td>Wall function k-&omega;(SST) (Menter)</td><td>0.637 m s<sup>-1</sup></td><td>6.37 &times; 10<sup>5</sup></td><td>0.408</td></tr><tr><td>M65e3</td><td>Wall function k-&omega; (Wilcox)</td><td>0.624 m s<sup>-1</sup></td><td>6.24 &times; 10<sup>5</sup></td><td>0.417</td></tr><tr><td>M65e3</td><td>Wall function k-&omega;(BSL) (Menter)</td><td>0.630 m s<sup>-1</sup></td><td>6.30 &times; 10<sup>5</sup></td><td>0.413</td></tr></tbody></table>

At the end of the turn, one can observe a slight way back of the direction of the body which is due to the transition between the turn and the straightforward swimming (Fig. 16). A modification of the law can be considered to avoid this problem, although one must recall that a real eel probably does not impose instantaneously a uniform curvature to its whole body to turn. A more detailed study of this point will be performed in the future to achieve a more realistic maneuver. Fig. 18 depicts the positions of the body and the pressure field in the plane of symmetry $Z = 0$ at different times: $t = 1s$ corresponds to the beginning of the P1 phase. Fig. 19 shows the pressure field on the body shape at the middle of the P2 phase.

## 7. Conclusion

Among numerous species which can be observed in Nature, fishes have always fascinated human beings with their incredible ease of coupling instinctively Newton's laws and Navier-Stokes equations in order to swim. The efficiency
<figure><img src="image_16.png" /><figcaption>Fig. 15. Evolution of the longitudinal head speed with near-wall low-Reynolds number turbulence closures.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16. Motion of the head for the different laws.</figcaption></figure>

and the great maneuverability of their locomotion function developed during their evolution provide a tremendous application field for fluid engineering illustrated by the growing number of related studies.

With the numerical methods described in this paper, it is now possible to simulate flows around “virtual fishes”. By imposing a realistic temporal evolution of the shape of the body, the response of the fluid induced by the deformation creates force and then motion. As shown by the applications which are illustrated in this paper, straightforward swimming but also fast maneuvers like turns can now be simulated. Apart from the shape of the body, the quality of the biomimesis only depends on agreement of the imposed laws with the real studied species. These laws have been extracted from various studies related to kinematics analysis (Gray, 1932; D’Août and Aerts, 1999; Wardle et al., 1995;
<figure><img src="image_18.png" /><figcaption>Fig. 17. Head speed modulus for the different laws.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 18. Horizontal view of the pressure in the plane of symmetry $Z = 0$.</figcaption></figure>

Videler, 1993). Simulations provide quantities which can hardly be measured on living species, like drag and efficiency, and then the optimization of the parameters involved in the law of deformation for several kinds of maneuverings is a natural subject of investigation for a research team already involved in shape optimization and optimal flow control (Duvigneau and Visonneau, 2004a, b). These future studies may contribute to a better understanding of the physical mechanisms involved in the art of swimming. To explore, analyse and appropriate the various solutions given by Nature represent one of the biomimesis purposes. This study is a first modest step in that direction even if the understanding and the modeling of all the sophisticated solutions developed by a real fish to take into account its environment will not be reached in a near future.
<figure><img src="image_20.png" /><figcaption>Fig. 19. Pressure field on the body at $t = 2.5$ s.</figcaption></figure>

## Acknowledgements

The authors gratefully acknowledge the scientific committee of CINES (project dmn2050) and IDRIS (project 1308) for the attribution of CPU time. This study is partly supported by the CNRS Project Robea.

## References

Chen, Z., Doi, Y., July 2002. Numerical study on propulsion by undulating motion in laminar-turbulent flow. 24th Symposium on Naval Hydrodynamics, Fukuoka, Japan.

D'Août, K., Aerts, P., 1999. A kinematic comparison of forward and backward swimming in the eel anguilla anguilla. Journal of Experimental Biology 202, 1511–1521.

Demirdiç, I., Muzaferija, S., 1995. Numerical method for coupled fluid flow, heat transfer and stress analysis using unstructured moving meshes with cells of arbitrary topology. Computer Methods in Applied Mechanics and Engineering 125, 235–255.

Duvigneau, R., Visonneau, M., 2004a. Hybrid genetic algorithms and artificial neural networks for complex design optimization in cfd. International Journal for Numerical Methods in Fluids 44 (11), 1257–1270.

Duvigneau, R., Visonneau, M., 2004b. Simulation and optimization of aerodynamic stall control using a synthetic jet. Second AIAA Conference on Flow Control, Portland, USA.

Duvigneau, R., Visonneau, M., Deng, G., 2003. On the role played by turbulence closures for hull shape optimization at model and full scale. Journal of Marine Science and Technology 153 (8).

Farhat, C., Degand, C., Koobus, B., Lesoinne, M., 1998. Torsional springs for two-dimensional dynamic unstructured fluid meshes. Computer Methods in Applied Mechanics and Engineering 163, 231–245.

Ferziger, J., Perić, M., 2002. Computational Methods for Fluid Dynamics, third ed. Springer, Berlin.

Gaskell, P., Lau, A., 1988. Curvature compensated convective transport: SMART a new boundedness preserving transport algorithm. International Journal for Numerical Methods in Fluids 8, 617–641.

Gray, J., 1932. Study in animal locomotion I—the movement of the fish with special reference to the eel. Journal of Experimental Biology 10, 88–104.

Jasak, H., 1996. Error analysis and estimation in the finite volume method with application to fluid flows. Ph.D. Thesis, Imperial College, University of London.

Leonard, B., 1988. Simple high-accuracy resolution program for convective modelling of discontinuities. International Journal for Numerical Methods in Fluids 8, 1291–1318.
Liu, H., Kawachi, K., 1999. A numerical study of undulatory swimming. Journal of Computational Physics 155, 223-247.

McDonald, H., Whitfield, D., 1997. Self-propelled maneuvering underwater vehicles. 21st Symposium on Naval Hydrodynamics, Washington, USA, pp. 478-489.

Rhie, C., Chow, W., 1983. A numerical study of the turbulent flow past an isolated airfoil with trailing edge separation. AIAA Journal 17, 1525-1532.

Söding, H., 2001. How to integrate free motions of solids in fluids. Fourth Numerical Towing Tank Symposium, Hamburg, Germany.

Spalart, P., Allmaras, S., 1992. A one-equation turbulence model for aerodynamic flows. AIAA 30th Aerospace Sciences Meeting. AIAA Paper 92-0439, Reno, NV, USA.

Triantafyllou, M.S., Triantafyllou, G.S., Yue, D.K.P., 2000. Hydrodynamics of fishlike swimming. Annual Review of Fluid Mechanics 32, 33-53.

Videler, J., 1993. Fish Swimming. Fish and Fisheries Series 10. Chapman & Hall, London.

Wardle, C., Videler, J., Altringham, J., 1995. Tuning in to fish swimming waves: body from, swimming mode and muscle function. Journal of Experimental Biology 198, 1629-1636.