Seyed Amin Ghaffari

Stéphane Viazzo

Kai Schneider

Patrick Bontoux
email: bontoux@l3m.univ-mrs.fr

Simulation of forced deformable bodies interacting with two-dimensional incompressible flows: Application to fish-like swimming

Keywords: Fluid interaction with forced deformable bodies, Compact fourth-order direct Poisson solver, Volume penalization, Fish swimming/turning

We present an efficient algorithm for simulation of deformable bodies interacting with two-dimensional incompressible flows. The temporal and spatial discretizations of the Navier-Stokes equations in vorticity stream-function formulation are based on classical fourth-order Runge-Kutta and compact finite differences, respectively. Using a uniform Cartesian grid we benefit from the advantage of a new fourth-order direct solver for the Poisson equation to ensure the incompressibility constraint down to machine zero. For introducing a deformable body in fluid flow, the volume penalization method is used. A Lagrangian structured grid with prescribed motion covers the deformable body interacting with the surrounding fluid due to the hydrodynamic forces and moment calculated on the Eulerian reference grid. An efficient law for curvature control of an anguilliform fish, swimming to a prescribed goal, is proposed. Validation of the developed method shows the efficiency and expected accuracy of the algorithm for fish-like swimming and also for a variety of fluid/solid interaction problems.

Introduction

The quantification and simulation of the flow around biological swimmers is one of the challenges in fluid mechanics [START_REF] Sotiropoulos | Immersed boundary methods for simulating fluid-structure interaction[END_REF]. At the same time bio-inspired design of swimming robots are in growth [START_REF] Rafei | Multi-variable constrained control approach for a threedimensional Eel-like robot[END_REF]. The costs of experimental studies [START_REF] Belkhiri | Modélisation dynamique de la locomotion compliante : Application au vol battant bio-inspiré de l'insecte[END_REF] lead the researchers to develop for efficient predictive numerical algorithms for the hydrodynamic analyses of fish swimming. Difficulties of numerical simulations of fish-like swimming are due to different reasons. One problem is efficient quantification of the kinematics of different species which seems to be far from the proposed simple laws in different studies. However, the main swimming mechanism in the majority of anguilliform fishes consists of a sinusoidal wave enveloped by a profile, which is created by the backbone of the fish and moves from head to tail. The tail beat creates a reversed Kármán street of vortices and generates thrust, leaving thus a momentumless wake back. Efficient simulation of incompressible flows is also an important problem, where the efficiency of the elliptic solver is crucial. The third bottleneck in numerical simulations of fish-like swimming is the coupling of the fluid solver with deformable, moving and rotating bodies. Anguilliform fishes add a constant curvature to their backbone for turning, i.e., they use their body like a rudder for torque generation. Yeo et al. [START_REF] Yeo | Simulation of fish swimming and manoeuvring by an SVD-GFD method on a hybrid meshfree-Cartesian grid[END_REF] studied numerically the straight swimming/cruising and sharp turning manoeuvres in two-dimensions. A carangiform-like swimmer was shown to execute a sharp turn through an angle of 70 o from straight coasting within a space of about one body length. Bergman and Iollo [START_REF] Bergmann | Modeling and simulation of fish-like swimming[END_REF] performed numerical simulations of fish swimming/rotation, looking for food at a prescribed point. They added a radius to the backbone of the fish to perform a rotation. The considered fish is constructed by a complex valued mapping like the Kutta-Joukowski transform. Here we will present a simple law for turning of a anguilliform fish. Our rotation control law [START_REF] Bontoux | An efficient algorithm for simulation of forced deformable bodies interacting with incompressible flows[END_REF] is similar to that presented in [START_REF] Yeo | Simulation of fish swimming and manoeuvring by an SVD-GFD method on a hybrid meshfree-Cartesian grid[END_REF] and [START_REF] Bergmann | Modeling and simulation of fish-like swimming[END_REF], in which the feedback is based on the angle between the line-of-sight and the direction of advance. But instead of adding a radius to the backbone we envisage to work with curvature which seems to be more efficient. We are using the method proposed by Boyer et al. [START_REF] Boyer | Macro-continuous computed torque algorithm for a three-dimensional Eel-like robot[END_REF] which is based on quaternions for efficient description of the fishes backbone kinematics. We are going to apply the rotation control to two-dimensional swimming. Even if due to the shape and deformation style of the fish-like swimmers the surrounding flow is fully three dimensional, most of the fundamental features of swimming are included in twodimensional analyses. For two dimensional incompressible flows the Navier-Stokes equations can be reformulated in terms of vorticity and stream-function which is more efficient. The choice of a finite difference method (FDM) in this paper is related to the use of an immersed boundary method in which a Cartesian grid can be used. Hence the use of FDM is straightforward. Among finite difference methods high order compact methods [START_REF] Hirsh | High order accurate difference solutions of fluid mechanics problems by a compact differencing technique[END_REF][START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF] are more advantageous in terms of accuracy and reasonable cost. We refer to [START_REF] Abide | A 2D compact fourth-order projection decomposition method[END_REF] and [START_REF] Boersma | A 6th order staggered compact finite difference method for the incompressible Navier-Stokes and scalar transport equations[END_REF] for high-order compact discretizations of the incompressible Navier-Stokes equations in primitive variables and to [START_REF] Bontoux | Analyse et optimisation d'une méthode de haute précision pour la résolution des équations de Navier-Stokes instationnaires[END_REF]- [START_REF] Roux | Optimisation of Hermitian methods for Navier-Stokes equations in the vorticity and stream-function formulation[END_REF] and [START_REF] Spotz | High-order compact scheme for the steady stream-function vorticity equations[END_REF] for the vorticity and stream-function formulation. Solving the incompressible Navier-Stokes equations typically implies an elliptic Poisson equation which is the most time consuming part of the algorithm. Direct methods like diagonalization or iterative methods (e.g. red-black PSOR, multi-grid and Krylov subspace methods) can be used. Using high-order discretizations iterative methods are less attractive because the resulted matrices are less sparse, thus the rate of convergences are slow. However iterative methods can cover all types of boundary conditions, see Spotz and Carey [START_REF] Spotz | High-order compact scheme for the steady stream-function vorticity equations[END_REF] for a fourth-order compact discretization of Poisson equation. On the other hand in direct methods the memory limitation is restrictive for simulations on fine grids. Therefore decoupling of the directions by FFT based methods can be advantageous, however this method implies some limitations in the boundary conditions. We are presenting a new fourth-order solver for the Poisson equation which is a combination of a compact finite difference with a sine FFT. The main advantages of our method are fourth-order accuracy, efficiency, the possibility to parallelize and convergence down to zero machine precision. Other advantages and limitations of the proposed solver are discussed in the paper. A difficulty in numerical simulations of fish swimming is the analysis of fluid/solid interaction, which can be handled by strong or loose coupling according to implicit or explicit time advancement, cf. [START_REF] Sotiropoulos | Immersed boundary methods for simulating fluid-structure interaction[END_REF] for a detail discussion. We are using the volume penalization method which was proposed by Angot et al. [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF] and belongs to the so called immersed boundary methods (IBMs). It consists of modeling the immersed body as a porous media, thus getting ride of the Dirichlet boundary conditions by considering both the fluid and the body as one domain with different permeabilities. So one can consider a rectangular solution domain in which the body is immersed and can even move. The penalization method will lead to first-order accuracy near the body but it is one of the most efficient methods in dealing with deformable, moving and rotating bodies in fluids. An extension to include elasticity of the solid in one direction is done in [START_REF] Engels | Two-dimensional simulation of the fluttering instability using a pseudospectral method with volume penalization[END_REF]. We refer to the review of Mittal and Iaccarino [START_REF] Mittal | Immersed boundary methods[END_REF] for a complete classification and description of immersed boundary methods. In the present work we will focus on some numerical aspects of efficient turning laws a topic which is less studied so far. To this end the method of quaternions is adapted to backbone kinematics description. Starting by the code developed in [START_REF] Sabetghadam | Implementation of vortex stretching into the two-dimensional Navier-Stokes equations via arbitrary external straining[END_REF] we are applying compact finite differences to the vorticity stream-function formulation of the Navier-Stokes equations including the penalization term. An efficient direct method is presented for solving the Poisson equation. Thus different numerical aspects of algorithm like accuracy in space and the error introduced by the penalization method will be examined. The code is developed in FORTRAN and is accessible for all [START_REF]The code is developed in FORTRAN and is accessible for all by sending a mail to: ghaffari@L3m.univ-mrs[END_REF]. The paper is organized as follows. First our methodology in-cluding the governing equations, discretization, kinematics of a fish like-swimming and the algorithm for fluid interaction with forced deformable bodies will be presented. Then a validation of the algorithm will be done and errors convergence will be studied. Results for swimming and rotation control will be also reported. Finally, the results will be discussed and some guides for the future works will be addressed.

Methodology

Governing equations of incompressible flow

The governing equations of incompressible flows are the Navier-Stokes equations. In two-dimensional problems the vorticity and stream-function formulation has the advantage that it not only eliminates the pressure variable entirely, but also ensures a divergence-free velocity field (mass conservation, i.e., ∇ • u = 0), if the Poisson equation (2) properly satisfied, see [START_REF] Bontoux | Analyse et optimisation d'une méthode de haute précision pour la résolution des équations de Navier-Stokes instationnaires[END_REF] and [START_REF] Roux | Optimisation of Hermitian methods for Navier-Stokes equations in the vorticity and stream-function formulation[END_REF]. One encounters two scalar valued quantities, i.e., the vorticity ω and the stream-function ψ, instead of the velocity vector and the pressure field, thus it makes the computations more efficient. With this formulation it is possible to use a collocated grid without adding any explicit numerical dissipation which reduces the arithmetics considerably. Therefore we continue with this formulation, but the concepts can also be extended to primitive variable formulation. By taking the curl of the Navier-Stokes equations, one obtains the vorticity transport equation:

∂ t ω + (u • ∇)ω = ν∇ 2 ω + ∇ × F , x ∈ Ω ∈ R 2 (1)
where ω(x, t) = ∇×u = v x -u y denotes the vorticity component which is normal to the considered two-dimensional plane, Ω is the spatial domain of interest, given as an open subset of R 2 , which can be bounded or unbounded in general, u(x, t) is the velocity field, ν = µ/ρ f > 0 is the kinematic viscosity of the fluid, ρ f is the density and F(x, t) is a source term. For a complete description of a particular problem, the above equations need to be complemented to describe an initial/boundary value problem (IBVP). The equation is parabolic in time and the velocity components are

(u, v) = (∂ y ψ, -∂ x ψ), with ψ satisfying a Poisson equa- tion -∇ 2 ψ = ω (2)
which is an elliptic equation in space. The penalization term for unit mass of the fluid reads,

F = -η -1 χ(u -u P ) (3)
where u P (x, t) is the velocity field of the immersed body. The Navier-Stokes equations are written for unit mass of the fluid, therefore the dimension of the source term F is acceleration, i.e., [LT -2 ]. The penalization parameter η is the permeability (porosity) coefficient of the immersed body with dimension [T ]. The mask (characteristic) function χ is dimensionless and describes the geometry of the immersed body

χ(x, t) = { 1 x ∈ Ω b 0 x ∈ Ω f (4)
where Ω f represents the domain of the flow and Ω b represents the immersed body in the domain of the solution.

The solution domain Ω = Ω f ∪ Ω b is governed by the Navier-Stokes equations in the fluid regions and by Darcy's law in the penalized regions, when η → 0.

Spatial discretization

In the present investigation an explicit second-order and an implicit fourth-order compact finite difference methods [START_REF] Hirsh | High order accurate difference solutions of fluid mechanics problems by a compact differencing technique[END_REF] is used for discretization of the spatial derivatives. The advantage of the compact methods over explicit finite differences is illustrated in terms of the scaled modified wavenumber w = k∆x in Fig. 1. For a given periodic function f (x) = e ikx , x ∈ [0, 2π] with known analytical derivatives f ′ (x) = ike ikx , f ′′ (x) = -k 2 e ikx , a numerical approximation of the derivatives at point x j have the form f ′ (x j ) = ik ′ e ikxj and f ′′ (x j ) = -k ′′2 e ikxj . The difference between exact and numerical approximation of the wavenumber is a measure of the discretization error which is purely dispersive for the first derivative and dissipative for the second derivative, if the considered function and discretization is periodic. For the second derivative the scaled modified wavenumbers are compared in Fig. 1 (b) with analytical values given by Lele [START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF] for some explicit and implicit differentiation methods. A good agreement between the numerical approximations of the scaled modified wavenumbers and the analytical values can be observed. It must be noted that the error in terms of the modified wavenumber is not necessarily sensitive to the formal order of the truncation error obtained by Taylor expansion analysis. The desired characteristics of finite difference schemes are better studied by directly optimizing the scheme in Fourier space rather than looking for the lowest truncation error. For example spectral like fivediagonal finite difference schemes designed by Lele [START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF] or Kim [START_REF] Kim | Optimised boundary compact finite difference schemes for computational aeroacoustics[END_REF] are formally fourth order, see Fig. 1 (a). Given the values of a function f on a uniformly spaced mesh x i = (i -1)h, for (i = 1, ..., N ) with h = L x /(N -1). Following [START_REF] Hirsh | High order accurate difference solutions of fluid mechanics problems by a compact differencing technique[END_REF][START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF] a fourth-order approximation of the first and second derivatives are obtained by the classical Padé schemes:

f ′ i-1 + 4f ′ i + f ′ i+1 = 3(f i+1 -f i-1 )/h ( 5
)
f ′′ i-1 + 10f ′′ i + f ′′ i+1 = 12(f i+1 -2f i + f i-1 )/h 2 (6)
for i = 2, ..., N -1, near the boundaries a third-order forward/backward stencil can be used, see [START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF]. A direct solver (LU-decomposition) will be applied to the tri/pentadiagonal system of linear equations along each line. The computational cost of a tridiagonal implicit method is in general three times of the one explicit method and a pentadiagonal linear system twice that of the tridiagonal.

Time integration

A classical fourth-order Runge-Kutta (RK4) method is used for time integration of the penalized vorticity transport equation [START_REF] Taylor | Stability of a viscous liquid contained between two rotating cylinders[END_REF], by assembling all spatial derivatives in the operator k(ω), one has

ω n+1 = ω n + ∆t 6 (k 1 (ω) + 2k 2 (ω * ) + 2k 3 (ω * ) + k 4 (ω * )) (7)
where

k i (ω) = -∂ y ψ i ∂ x ω +∂ x ψ i ∂ y ω +ν∇ 2 ω +∂ x F y -∂ y F x ( 8
)
for i = 1, 2, 3, 4. In each time step Eq. ( 8) must be evaluated four times where Eq. ( 2) must be solved for updating of stream-function (-∇ 2 ψ i = ω * i ). For details and technical discussions of Runge-Kutta methods see [START_REF] Press | Numerical Recipes[END_REF]. However, ∆t is limited by the CF L (Courant-Friedrichs-Lewy) condition which implies that

U ∆t ∆x ≤ CF L ≈ σ i w ′ max ( 9
)
where U is an advection velocity (or a phase speed). In the presence of the nonlinearity in space more attention must be payed. Viscous terms imply a constraint of the form ν∆t

∆x 2 ≤ V SL ≈ σ r w ′′ max (10)
on the time-step, where σ r = 2.9 and σ i = 2.85 are limits of real and imaginary parts of the stable region on complex plan for RK4 method. w ′ max = 1.74 and w ′′ max = 6 are the maximum values of the scaled modified wave numbers for the first and second derivatives calculated via fourthorder Padé scheme, plotted in Fig. 1. It can be seen that w ′ max ∈ [1, π] and w ′′ max ∈ [4, π 2 ] for different approximations of spatial derivatives. Therefore with the use of a high-order method for spatial discretization smaller timesteps must be used. In the presence of moving bodies the displacement of the moving body must not exceed the grid spacing, i.e., ∆t ≤ ∆x/u B . Moreover, by using the explicit penalization method another constraint, ∆t ≤ η, must be respected. Among the four above-mentioned constraints, the smallest ∆t must be chosen. According to Gazzola et al. [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF] with the use of volume penalization method for simulation of fluid/solid interaction the order of accuracy in time is reported to be first-order independent of the time accuracy of the underlying method for time integration of the Navier-Stokes equations. However the larger stability bound of the RK4 is still attractive to enhance the time step. 1: Plots of the scaled modified wavenumber for the first and second derivative, w ′ (w) and w ′′ (w) respectively, versus the scaled wavenumber w = k∆x for different central finite difference methods presented in Lele [START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF] and Kim [START_REF] Kim | Optimised boundary compact finite difference schemes for computational aeroacoustics[END_REF].

for the pressure correction in projection methods. In the vorticity stream-function formulation equation ( 2) has to be solved with Dirichlet boundary condition for vorticity and stream-function. Free slip (ω = 0) boundary conditions in a close rectangular domain (ψ = 0, all around) is applied in all the test cases studied in the present investigation. In the presence of periodic boundary conditions, FFT based direct solvers can be used to efficiently solve the Poisson equation with high accuracy. Even if the flow is not periodic in all directions, like most of the practical problems, in accordance with the boundary conditions for the elliptic equation (homogeneous Dirichlet/Neumann) sine or cosine FFTs can be used, see [START_REF] Kim | Application of a fractional-step method to incompressible Navier-Stokes equations[END_REF]- [START_REF] Orlandi | Fluid flow phenomena: A numerical toolkit[END_REF] and [START_REF] Laizet | A numerical strategy to combine high-order schemes, complex geometry and parallel computing for high resolution DNS of fractal generated turbulence[END_REF]. We are presenting a new direct fourth-order solver for the Poisson equation ( 2) which is a combination of a compact finite difference with a sine FFT. The advantages of our method are fourth-order accuracy, convergence down to zero machine precision, compact tridiagonal stencil, possibility of extension to three dimensions, reduced arithmetics and memory usage in comparison to iterative methods. Moreover the parallelization is straightforward because of decoupling of the operations in different directions. Nearly linear strong scaling (speed up) and efficiency is reported by Laizet and Lamballais [START_REF] Laizet | A numerical strategy to combine high-order schemes, complex geometry and parallel computing for high resolution DNS of fractal generated turbulence[END_REF] for a similar direct solver. They introduced a dual domain decomposition (or pencil) method, in which information along a line is accessible for a CPU by alternative decomposition of the domain in three directions. The limitation of our method (moreover the boundary conditions) is the use of a uniform grid in the direction in which the FFT is applied. When the solver of the parabolic part consists of finite-differences, it is a custom to use a FDM discretization in one direction with-out loss of accuracy and efficiency via a direct tridiagonal solver. The advantage of this approach is the possibility of applying general boundary conditions in that direction and using a refined mesh. The second-order version of this solver can be found in [START_REF] Orlandi | Fluid flow phenomena: A numerical toolkit[END_REF]. For a compact fourth-order collocated discretization of the Poisson equation -∇ 2 ψ = ω, over N x × N y grid points, by using

∂ 2 ψ ∂x 2 = δ 2 x ψ - ∆x 2 12
∂ 4 ψ ∂x 4 + O(∆x 4 ) ( 11
)
where δ 2

x represents a central second-order estimation of the second derivative, for the x-direction we obtain

(δ 2 x - ∆x 2 12
∂ 4 ∂x 4 + ∂ yy )ψ = -ω (12)
because of the presence of the ∆x 2 factor behind the fourthorder derivative, this term cannot be dropped and must be evaluated with second-order accuracy. Therefore, the hole approximation scheme yields fourth-order accuracy. The fourth-order derivative can be evaluated by using the original Poisson equation -∇ 2 ψ = ω, and successive differentiation with respect to x (i.e.,

∂ xx ∂ xx ψ = -∂ xx ∂ yy ψ-∂ xx ω). Replacing ∂ xx by δ 2 x , we find (δ 2 x + ∆x 2 12 δ 2 x ∂ yy + ∂ yy )ψ = -ω - ∆x 2 12 δ 2 x ω (13)
By applying Fourier transform in y direction over Eq. ( 13) and replacing second derivatives by -k 2 y ψ in Fourier space, we have

(δ 2 x - ∆x 2 12 δ 2 x k ′2 y -k ′2 y ) ψ = -ω - ∆x 2 12 δ 2 x ω (14)
Usually the exact wavenumber is replaced by the modified wavenumber k ′2 y which permits to evaluate the difference between the finite-difference and the spectral approximation of the second derivative [START_REF] Orlandi | Fluid flow phenomena: A numerical toolkit[END_REF]. For a fourth-order explicit finite-difference discretization, analytical relation for the scaled modified wavenumber of second derivative is given by Lele [START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF] as follows

k ′2 y = 1 ∆y 2 [ 8 3
( 1 -cos( k y π N y ) ) - 1 6
(

1 -cos( 2k y π N y ) ) ] (15)
Comparison with numerical values in Fig. 1 (b) confirm that Eq. ( 15) is exact. The final tridiagonal system to be solved in Fourier space for each wavenumber of ψ in the y-direction is 16)

β ψi+1,m -(2β + k ′2 y ) ψi,m + β ψi-1,m = -(ω i+1,m + 10ω i,m + ωi-1,m )/12 (
for i = 2, ..., N x -1, where β = ∆x -2 -k ′2 y /12.
In summary, first a one-dimensional direct-FFT of the forcing function is performed in the y-direction. Then for each line in the x-direction the tri-diagonal system ( 16) must be solved to find the solution ψ in wavenumber space. Next an inverse-FFT of the solution is performed. For the real data with zero value at the boundaries (homogeneous Dirichlet boundary condition, i.e., ψ = ω = 0), the natural Fourier transform to use is the sine transform, see e.g. [START_REF] Press | Numerical Recipes[END_REF]. The direction of FDM and FFT can be changed to consider no-slip boundary condition in the y-direction. In order to take into account inflow/outflow boundary conditions the mean flow must be reduced from the total velocity field u = U -U ∞ in the vorticity transport equation (1) to impose ψ = 0 at the boundaries. This is equivalent to move the grid with U ∞ and writing the Navier-Stokes equations in a moving reference frame for the peturbed velocity field u instead of a Galilean inertial frame [START_REF] Rossinelli | GPU accelerated simulations of bluff body flows using vortex particle methods[END_REF].

Algorithm of the fluid interaction with a deformable body 2.2.1. Kinematics of the fish

The geometrically exact theory of nonlinear beams, is developed by Simo [START_REF] Simo | A finite strain beam formulation. The threedimensional dynamic problem[END_REF] and extended for fish vertebral by Boyer et al. [START_REF] Boyer | Macro-continuous computed torque algorithm for a three-dimensional Eel-like robot[END_REF]. In this theory, the beam is considered as a continuous assembly of rigid sections of infinitesimal thickness, i.e., a one-dimensional Cosserat medium. We are summarizing the exact kinematics of the fish backbone in three dimensions for interested readers and future developments, but all the cases in this paper are limited to two-dimensions. Following [START_REF] Boyer | Macro-continuous computed torque algorithm for a three-dimensional Eel-like robot[END_REF], [START_REF] Rafei | Multi-variable constrained control approach for a threedimensional Eel-like robot[END_REF] and [START_REF] Belkhiri | Modélisation dynamique de la locomotion compliante : Application au vol battant bio-inspiré de l'insecte[END_REF] the kinematics of the backbone for Eel-like fishes can be determined by integration along the arclength ξ ∈ [0, l fish ] starting with the head's orientation, position and velocities as boundary conditions. The variation of the orientation is obtained by

∂Q ∂ξ = 1 2 M (Ω)Q ( 17
)
where Q = (cos ϕ 2 , a x sin ϕ 2 , a y sin ϕ 2 , a z sin ϕ 2 ) T are unit normalized (q 2 0 +q 2 1 +q 2 2 +q 2 3 ) 1/2 = 1 quaternions that represent the head frame's orientation with respect to the inertial frame and M (Ω) is an anti-symmetric tensor

M (Ω) =     0 -ω 1 -ω 2 -ω 3 ω 1 0 ω 3 -ω 2 ω 2 -ω 3 0 ω 1 ω 3 ω 2 -ω 1 0     (18)
where Ω = (ω 1 , ω 2 , ω 3 ) T denotes the mean angular velocity. The geometry R = (x, y, z) T in the Galilean reference frame is given by

∂R ∂ξ = Rot(Q)K (19)
where k 2 and k 3 in K = (k 1 , k 2 , k 3 ) T stand for the Eel's backbone transversal curvature and k 1 represents the rate of rotation (twist) of the section around the backbone with the normal aligned with ξ-direction. The rotation matrix in terms of the quaternions is then given by

Rot = 2   q 2 0 + q 2 1 -1 2
q 1 q 2 -q 0 q 3 q 1 q 3 + q 0 q 2 q 1 q 2 + q 0 q 3 q 2 0 + q 2 2 -1 2 q 2 q 3 -q 0 q 1 q 1 q 3 -q 0 q 2 q 2 q 3 + q 0 q 1 q 2 0 + q 2 3 -

1 2   (20)
The variation of linear V = (v 1 , v 2 , v 3 ) T and angular Ω = (ω 1 , ω 2 , ω 3 ) T velocities in the local frame, i.e., the frame attached to the body are given by ∂ ∂ξ

[ V Ω ] = - [ K ∨ Γ ∨ 0 K ∨ ] [ V Ω ] + [ Γ K ] (21)
where (•) represents the time derivative, (∨) stands for the anti-symmetric matrix constructed from a given vector, e.g.,

Γ ∨ =   0 -γ 3 γ 2 γ 3 0 -γ 1 -γ 2 γ 1 0   (22)
where Γ = (γ 1 , γ 2 , γ 3 ) T represents the local transversal shearing of the sections whose first component is the stretching rate along the ξ-direction. The accelerations can also be deduced from the time derivative of Eq. ( 21). For more details we refere to [START_REF] Boyer | Macro-continuous computed torque algorithm for a three-dimensional Eel-like robot[END_REF], [START_REF] Rafei | Multi-variable constrained control approach for a threedimensional Eel-like robot[END_REF] and [START_REF] Belkhiri | Modélisation dynamique de la locomotion compliante : Application au vol battant bio-inspiré de l'insecte[END_REF]. To find the velocities in the frame attached to the head from the velocities V G in the Galilean reference frame and inverse,

(v 1 , v 2 , v 3 ) T = Rot T (v x , v y , v z ) T (23)
can be used. By considering N (1, ..., N points ) discrete points on the Eel's backbone, equations ( 17), ( 19) and [START_REF] Boyer | Macro-continuous computed torque algorithm for a three-dimensional Eel-like robot[END_REF] altogether must be integrated in space by a proper numerical method (N eq = 13 in 3D). We are using a fourthorder Runge-Kutta method for integration and comparisons with a first-order Euler method shows that RK4 is more precise especially when the number of the points along the Eel's backbone is less than N points = 30.

Lagrangian structured grid

The first choice to start the parameterization of the swimmer's body is a symmetric shape. The geometry of a two-dimensional swimmer can be characterized by the half width w(ξ) of the body along its arclength (midline) ξ ∈ [0, l fish ]. Following the work of Kern and Koumoutsakos [START_REF] Kern | Simulations of optimized anguilliform swimming[END_REF] and Carling et al. [START_REF] Carling | Self-propelled anguilliform swimming: simultaneous solution of the two-dimensional Navier-Stokes equations and Newton's laws of motion[END_REF], the half width w(ξ) is defined as

w(ξ) =      √ 2w h ξ -ξ 2 0 ≤ ξ < s b w h -(w h -w t )( ξ-s b st-s b ) 2 s b ≤ ξ < s t w t l fish -ξ l fish -st s t ≤ ξ ≤ l fish (24)
where l fish is the body length, w h = s b = 0.04l fish , s t = 0.95l fish and w t = 0.01l fish . In the mid part of the fish a linear function can be used as in Gazzola et al. [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF]. A structured grid formed by normal to backbone lines with thickness given by ( 24) covers the body. The velocity components of each point on the Lagrangian grid V shape with (I, J) indexes are given by

V shape (I, J) = V BN (I) + r(I, J)Ω BN (I) ( 25
)
where V BN and Ω BN are the linear and angular velocities of the backbone respectively, given by Eq. ( 21). The radius (-w < |r| < w) is aligned with the transversal lines of the structured grid normal to the backbone. Fig. 3 shows an example of the Lagrangian grid covering the fish after deformation in which the corresponding velocities of each point are illustrated. The information of the Lagrangian structured grid covering the deformable body must be transfered to the Eulerian-Cartesian grid by interpolation to find χ(i, j) and u p (i, j). To determine χ(i, j) on the Eulerian grid whose first point (x, y) (1,1) = (0, 0) is located at the origin, the coordinates of each point on the Lagrangian grid X shape (I, J) is divided by ∆x and ∆y.

After applying a correction to the integer part of the results they give the indexes (i, j) of the mask function χ on the Eulerian grid for which χ = 1 is assigned. Following [START_REF] Forestier | Etude par méthode spectral de sillages tridimensionnels en fluide stratifié[END_REF] and [START_REF] Minguez | Large eddy simulation of turbulent flows over complex geometries: Application to a simplified car model[END_REF] the mask is mollified by the Shuman [START_REF] Shuman | Numerical methods in weather prediction: II Smoothing and filtering[END_REF] filter χi,j = (2χ i,j + χ i+1,j + χ i-1,j + χ i,j+1 + χ i,j-1 )/6 [START_REF] Pasquetti | A pseudo-penalization method for high Reynolds number unsteady flows[END_REF] which is equivalent to raised cosine filter in Fourier space, we refer to [START_REF] Pasquetti | A pseudo-penalization method for high Reynolds number unsteady flows[END_REF] for more details. An example of the transfered geometry to the Eulerian grid, after smoothing by Eq. ( 26) is illustrated in Fig. 4 which represents the mask function. The spacing of the grid points on the Lagrangian grid must be fine enough in comparison to ∆X and ∆Y to accurately represent the deformation of the body. We are using two-dimensional linear interpolation

f (x, y) = axy + bx + cy + d (27)
which leads to a 4 × 4 linear system. To determine the coefficients the system is solved by a direct method, i.e., Gauss-Jordan elimination from [START_REF] Press | Numerical Recipes[END_REF]. For all points in the interior of the fish we have χ(i, j) = 1 on the Eulerian grid.

For each point of the Eulerian grid in which χ = 1 the four nearest points of the Lagrangian grid are used to solve the system formed by [START_REF] Minguez | Large eddy simulation of turbulent flows over complex geometries: Application to a simplified car model[END_REF]. For some points due to mollifying by Eq. ( 26) we have 0 < χ < 1, therefore the interpolation automatically becomes an extrapolation. Some points are completely outside of the original Lagrangian shape. An example of the interpolated velocity components on the Eulerian grid is illustrated in Fig. 4. Note that the resulting velocity field u p is not divergence-free, we refer to Gazzola et al. [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF] for a complete theoretical and numerical discussion about this subject. In the present investigation we are not considering this issue.

Hydrodynamic coefficients evaluation

With the use of the penalization method the hydrodynamic forces and the moments acting on the body, which are usually evaluated via surface integrals of the stress tensor σ(u, p) = µ(∇u + (∇u) T )/2 -p I, can computed directly by integrating the penalized velocity over the considered volume (surface in two-dimensions). Thus the forces in [Newton] are given by

F ⋆ = ∂Ωs σ • n dl = lim η→0 ρ f η ∫ Ωs χ(u -u B ) ds + ρ f S pen Ẍcg (28)
for the unit mass (m = ρ f S pen ) of the fluid. By definition F = F ⋆ /m, we have

F ≈ 1 ηS pen ∫ Ωs χ(u -u B ) ds + Ẍcg ( 29
)
The moment in two-dimensions is evaluated by

M Xcg = ∂Ωs (x -X cg ) × σ • n dl = lim η→0 ρ f η ∫ Ωs χ(x -x cg ) × (u -u B ) ds + ρ f ρ b J cg θcg (30) in [N • m],
where J cg = ∫ r 2 dm is the polar moment of inertia taken around center of the mass, n is the unit outward vector to ∂Ω s , X cg is the location of the center of gravity of the immersed body, θ is the angle of rotation with respect to the center of gravity. The dots denote derivatives with respect to time and S pen is the surface of the penalized area.

Denoising of the hydrodynamic coefficients

In dealing with fluid/solid interaction problems, the oscillation of the hydrodynamic forces and moments during successive iterations calculated from Eqs. ( 29) and [START_REF] Laizet | A numerical strategy to combine high-order schemes, complex geometry and parallel computing for high resolution DNS of fractal generated turbulence[END_REF] causes some trouble in correctly predicting the accelerations. The hydrodynamic forces and moments acting on the body are used to calculate the linear and angular accelerations which in turn have impact on the predicted velocity vector and the trajectory of the solid. The oscillations are due to the nature of the penalization method see the discussion in [START_REF] Kolomenskiy | A Fourier spectral method for the Navier-Stokes equations with volume penalization for moving solid obstacles[END_REF], insufficient resolution, the approximative nature of Eqs. ( 29) and [START_REF] Laizet | A numerical strategy to combine high-order schemes, complex geometry and parallel computing for high resolution DNS of fractal generated turbulence[END_REF]. The oscillations are like a noise and lead to an invalid result and even divergence of the simulations. An efficient method to eliminate them is to apply a low-pass filter like exponential smoothing usually used in denoising of time series. This filter is used in [START_REF] Kern | Simulations of optimized anguilliform swimming[END_REF] to denoise the hydrodynamic forces and moments. Simple exponential smoothing does not perform well when there is a trend in the data. In such situations, several methods were devised like second-order (double) exponential smoothing [START_REF] Holt | Forecasting trends and seasonal by exponentially weighted averages[END_REF]

F n = αF n + (1 -α)( F n-1 + b n-1 ) , n = 3, 4, . . . (31) b n = β( F n -F n-1 ) + (1 -β)b n-1 , (α, β) ∈ [0, 1] (32)
where F 1 = F 1 , for n = 2 one can use Eqs ( 31) and ( 32) with α = β = 1. Then α = 1 -(1 -δ) 2 and β = δ 2 /α can be used in which δ is a small band. According to our experience δ = 10 -3 performs well for denoising of the hydrodynamic forces, moments and centroid change of the deformable bodies. However δ = 10 -3 have a strong damping effect, bigger values, i.e., δ = 5 × 10 -3 has less damping effect but there is a risk of divergence in the simulations. A sensitivity analysis must be done for each test case (see also the discussion of the results in Section 3.4).

Body dynamics

The dynamics of an arbitrary solid or deformable body moving in a viscous incompressible fluid is governed by Newton's second law

Σ(F H + F G ) = m Ẍ ( 33
)
where the applied forces can be split into two components; the hydrodynamic forces F H and the forces due to gravity F G . Newton's law can be integrated directly to give the position of the center of gravity as a function of time.

Holding F constant over the discrete physical time step (t n , t n+1 ) yields

∆X cg = 1 2
F n m ∆t 2 + V n ∆t ( 34
)
and V n+1 = V n + Ẍ∆t. The rotational motion is described by Euler's equation

ΣM = d dt (J cg θ) ( 35
)
where M is the applied moment around the reference point. Time integration of Eq. ( 35) regardless of changes in the moment of inertia and M will give the new angle of the body with respect to a given reference ∆θ = 1 2 θn ∆t 2 + θn ∆t [START_REF] Morales | Simulation of confined magnetohydrodynamic flows using a pseudospectral method with volume penalization[END_REF] where θ = M/J and θn+1 = θn + θ∆t (the dotes denote derivation with respect to time). Eqs ( 34) and ( 36) describe a motion with three degrees of freedom (3DOF) for the considered body. In these equations second-order terms can be eliminated as done in [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF] but we are keeping these terms. Eqs. ( 29) and ( 30) provide the fluid forces and moment necessary to integrate the system of ODEs formed by Eqs. [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF] and Eq. [START_REF] Boersma | A 6th order staggered compact finite difference method for the incompressible Navier-Stokes and scalar transport equations[END_REF]. Denoising of forces and moment is done according to [START_REF] Yeo | Simulation of fish swimming and manoeuvring by an SVD-GFD method on a hybrid meshfree-Cartesian grid[END_REF]. Appropriate initial conditions are necessary. In the present computations we use a first-order scheme for time integration of the dynamics equations which seems to be adequate because of the error introduced by penalization method which is also first-order. The same time integration method is also used in [START_REF] Kolomenskiy | A Fourier spectral method for the Navier-Stokes equations with volume penalization for moving solid obstacles[END_REF] and [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF] for the dynamics of the body where the penalization is used. A summary of the algorithm for the fluid/structure interaction is given in Algorithm 1. The flowchart is illustrated in Fig. 5.

Algorithm 1 Fluid/structure interaction 1. Start from an initial condition 2. Body kinematics (a) (Just for the Eel) Create Eel's backbone by integrating Eqs ( 17), ( 19) and ( 21

Validation

The spatial convergence of the flow solver including the penalization term is verified using a Taylor-Couette flow for which an analytical solution is available. Then the solver for solid dynamics is verified via falling cylinder and ellipse in a quiescent fluid. Finally a test case of fish swimming in forward gait is compared with the results of Gazzola et al. [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF].

Spatial convergence for a Taylor-Couette flow

For a rigorous study of the error due to the penalization term added to the Navier-Stokes equation in vorticity and stream-function formulation an exact solution is necessary. The Taylor-Couette configuration is a good choice, first and foremost, because of known Dirichlet boundary conditions everywhere, and secondly, because of the presence of curved walls contrary to other existing analytical solutions usually available for Cartesian domains which thus coincide with the underlying Cartesian grid used to discretize the governing equations. Although the solver is adapted to a Cartesian domain the mask function which represents the penalized area for the Taylor-Couette flow is curved (see Fig. 6) as it is the case for flow around an ellipse or complex geometries which will consider in following. An explicit second-order finite difference method and an implicit fourth-order compact finite difference are used for discretization of the governing equations including the curl of the penalization term ∇ × F. Taylor-Couette flow

T a = R 1 (Ω 2 -Ω 1 ) 2 (R 2 -R 1 ) 3 ν -2 ≈ 1.
[1] consists of a viscous fluid confined between two concentric cylinders with radii (R1, R2) in rotation with different angular velocities (Ω 1 , Ω 2 ). For Taylor numbers below the critical value T a c ≈ 1708, the flow is steady and purely azimuthal, i.e., u z = u r = 0. This state is known as circular Taylor-Couette flow and for which there is an analytical solution [START_REF] Monin | Statistical Fluid Mechanics[END_REF]. The solution is given in cylindrical coordinates, the azimuthal velocity is u θ (r) = Ar + B/r where (r, θ)

∈ [R 1 , R 2 ] × [0, 2π], A = (Ω 2 R 2 2 -Ω 1 R 2 1 )/(R 2 2 -R 2 1 ) and B = R 2 1 R 2 2 (Ω 1 -Ω 2 )/(R 2 2 -R 2 1
) are known. The vorticity between the two cylinders is constant (ω z = 2A) and the stream-function is given by ψ(r) = -Ar 2 /2 -B ln(r) + c 0 where c 0 must be determined with respect to an arbitrary reference point. To use the volume penalization method, the velocity components must be enforced also in the solid regions from known angular velocities (i.e., Ω 1 and Ω 2 ), u θ (r) = rΩ where (r, θ)

∈ [0, R 1 ] ∪ [R 2 , R max ] × [0, 2π].
The vorticity inside the rotating regions is constant and is equal to twice of the domain's angular velocity (ω z = 2Ω) and the stream-function is given by, ψ(r) = Ω/2r 2 + c, where c must be determined for each domain in accordance with c 0 . A unit square domain is considered as the solution domain, the time-step of the RK4 method is calculated by the constraints presented in the Section 2.1.2 and the kinematic viscosity is fixed to ν = 0.01. The radii are chosen R 1 = 0.2 and R 2 = 0.4, respectively. At t = 0 the fluid domain is at rest and the inner-cylinder is set into movement with a fixed angular velocity (Ω 1 = 0.2) while the angular velocity of the outer cylinder is kept equal to zero (Ω 2 = 0). The Taylor number for this configuration (T a = 0.64) is below the critical value, thus the flow is purely azimuthal. The L 1 -error ∥u exact -u num η ∥ for u which is the x-component of the considered velocity field, is calculated for different penalization parameters η and resolutions (N in x and y directions). The simulations are carried out until a steady state is reached, so that the error is independent of the time discretization. The simulations are stopped when the time t end = 10 is reached. Original and mollified mask function at the midline y = 0.5 is illustrated in Fig. 7 (a), comparison of the computed vorticity ω, stream-function ψ and the u velocity component with exact solutions for N = 128 grid points in each direction is plotted in Fig. 7 (a)-(b). The convergence of the L 1 -error of u versus the grid resolution is shown in Fig. 7 (c) for different penalization parameters. Suppose u num η denotes the numerical solution of the penalized equations, for quantifying the numerical error of u num η compared to u exact (the solution to the original Navier-Stokes problem), the error can be estimated by

∥u exact -u num η ∥ ≤ ∥u exact -u η ∥ O( √ η) + ∥u η -u num η ∥ O(∆x p ) (37)
The first term at the right-hand side is the error due to the penalization term and the second term representing the discretization error (p being the formal order of accuracy of the numerical method used to discretize the equation). Where ∥ • ∥ is an appropriate norm. A compromise between two errors is to chose ∆x ≈ √ η, which will lead the to a first-order convergence for the error ∥u exact -u num η ∥ ≤ O(∆x) [START_REF] Nguyen Van Yen | Approximation of the Laplace and Stokes operators with Dirichlet boundary conditions through volume penalization: a spectral viewpoint[END_REF]. The convergence of the L 1error of u versus different penalization parameters is shown in Fig. 7 (d) for different grid resolutions, where the order √ η can be observed. For these calculations the expected formal accuracy is p = 2 and thus the convergence is between first and second order in space as a function of the resolution N , confirming the theoretical analysis of Carbou et al. [START_REF] Carbou | Boundary layer for a penalization method for viscous incompressible flow[END_REF] and the numerical results of Morales et al. [START_REF] Morales | Simulation of confined magnetohydrodynamic flows using a pseudospectral method with volume penalization[END_REF]. We also observe a saturation of the convergence error for large N , due to dominance of the penalization error. An optimal resolution can be found for each η. Decreasing η will lead to an accuracy enhancement in general but for an explicit integration time step ∆t = O(η) is limited by η as discussed in Section 2.1.2. A fine grid also needs a small η as can be seen in the Fig. 7.

Two-way fluid/solid interaction

In this section we attempt to perform a simulation of a two-dimensional cylinder falling in a quiescent fluid to validate the two-way fluid/solid interaction. We compare our results with those of Gazzola et al. [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF] and Namkoong et al. [START_REF] Namkoong | Numerical analysis of twodimensional motion of a freely falling circular cylinder in an infinite fluid[END_REF] which have the same physical parameters. A rigid 2D cylinder of diameter D = 0.005 m with ρ b = 1.01ρ f , is released from rest in a fluid with density ρ f = 996 kg/m 3 and kinematic viscosity ν = 8 × 10 7 m 2 /s and accelerates due to gravity (g = -9.81 m/s 2 ) until it reaches its asymptotic terminal velocity. The domain size is set to (x, y) ∈ [0 , 0.04 m] × [0 , 0.32 m]. The resolution in our simulations is set to 512 × 4096, the penalization parameter η = 10 -3 , the time step ∆t = 10 -3 and filter parameter for denoising of the hydrodynamic coefficients δ = 0.001. In the simulations of Gazzola et al. [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF] the resolution is 1024 × 8192, the penalization parameter η = 10 -4 and the Lagrangian CF L = 0.01. The snapshots of the vorticity isolines generated by the falling cylinder in a quiescent and slightly perturbed fluid are illustrated in Fig. 8 and Fig. 9, respectively. A qualitative agreement with the simulations of Gazzola et al. [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF] can be observed. Comparison of the vorticity structures at t = 13 between the simulation with slightly perturbed initial condition (u ′ ≈ 0.001 × random number × u ∞ ) represented in Fig. 9 (h) and that of quiescent initial condition represented in Fig. 8 (d) shows that presence of a perturbation in the initial condition can trigger the transition in the early stage of falling, i.e., t ≈ 3. Particularly this is important to obtain comparable results with other simulations with different numerical methods where the added numerical dissipation is not necessarily the same. Without adding any initial perturbation the transition can be triggered (e.g., at t ≈ 10) by the numerical errors which are performing like a perturbation (see Fig. 8). This is not controlled, it depends on grid resolution and the numerical implementation and explains the delayed streamwise velocity overshoot and the different transient flow fields. Fig. 10 shows the time evolution of the streamwise velocity obtained with the present method and those of Namkoong et al. [START_REF] Namkoong | Numerical analysis of twodimensional motion of a freely falling circular cylinder in an infinite fluid[END_REF] and Gazzola et al. [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF]. As can be seen the streamwise velocity shows the same dynamics as the reference simulations. In particular the streamwise velocity obtained by simulation with perturbed initial condition overshoots above the terminal velocity and then slows down when the vortices start shedding. The terminal streamwise velocity of the simulation with perturbed initial condition is u streamwise = 0.023 m/s resulting a Reynolds Re ≈ 144 and that of unperturbed initial condition is u streamwise = 0.025 m/s resulting a Reynolds Re ≈ 156. In the former an overshoot can be observed in the streamwise velocity while in the later the overshoot takes place in a bigger time interval or it is entirely eliminated. The terminal velocity differs 8% from the reference terminal velocity in the case of perturbed initial condition and coincide in the case of unperturbed initial condition. The differences are due to different Poisson solvers which is unbounded in the simulation of Gazzola et al. [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF], the boundary conditions which is free-slip and no-penetration in our simulations, different penalization parameters and resolutions. In the author's viewpoint the take-home message here is that the near one relative fluid/solid density leads to a small buoyancy where invalid approximation of the hydrodynamic coefficients specially in the early stages of the falling yields the simulation to a failure. To cope with this challenge the process of filtering of the hydrodynamic coefficients with a proper parameter is devised in the proposed algorithm. For bigger density rations ρ b > 1.5ρ f the difficulty is not serious and can be handled even without denoising of the hydrodynamic coefficients.

Validation of the solid dynamics via a falling ellipse

The validation of the body dynamics interacting with a fluid is done by a falling ellipse due to gravity in a fluid in rest. Different behaviors like steady falling, fluttering, tumbling and chaotic motion can be observed by varying the ellipse's geometry, density ratio and the viscosity of the fluid. These parameters can be summarized in dimensionless moment of inertia

J * cg = 2J cg /(a 4 ρ f ) = π(a 2 + b 2 )(b/2a 3 )(ρ b /ρ f )
and the Reynolds number Re = u t L/ν, where u t is the sedimentation average velocity estimated by

u t = √ 4bg(ρ b /ρ f -1)
see [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF]. The domain of the solution is (x, y) ∈ [0, 5L] × [0, 20L] where L = 2a = 1 and H = 2b = 0. and minor diameters of the ellipse respectively. The resolution of the grid is 512×2048, the polar moment of inertia around the center of gravity J cg = 0.25πab(a 2 + b 2 )ρ b , the initial position (x 0 , y 0 ) = (0.5L x , L y -3a) and the initial angle of the major diameter with respect to the horizon θ 0 = π/4. The density ratio is set to ρ b = 1.538ρ f , the filter parameter for denoising of the hydrodynamic coefficients δ = 0.001, the kinematic viscosity ν = 0.01 m 2 /s, the gravity in y-direction g = -9.81 m/s 2 and the penalization parameter η = 10 -3 . The trajectory and the isolines of the vorticity generated by the falling ellipse in the fluttering regime are illustrated in Fig. 11 at different instants from t = 0.2 up to t = 25. A qualitative agreement with the simulations of Kolomenskiy et al. [START_REF] Kolomenskiy | A Fourier spectral method for the Navier-Stokes equations with volume penalization for moving solid obstacles[END_REF] and Gazzola et al. [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF] can be observed. The corresponding forces and velocity components are plotted in Fig. 12. All of the hydrodynamic coefficients show an oscillatory behavior with a principal frequency f 1 ≈ 0.24. However in the side force a harmonic frequency with f 2 = 2f 1 ≈ 0.48 can be seen which is due to the shedding of the vortices. The chosen reference point for the simulation of the falling ellipse is the center of gravity (cg) for calculation of the polar moment of inertia, rotation angle and the moment. This choice is advantageous for simplification of the Euler Eq [START_REF] Boersma | A 6th order staggered compact finite difference method for the incompressible Navier-Stokes and scalar transport equations[END_REF] to not include the torque due to buoyancy. For the simulations of the fish swimming the buoyancy is set to zero, i.e. ρ b = ρ f . Thus the reference point can coincide with the head which is more suitable for creating the fish starting by the situation of the head.

Fish in forward gait

Anguilliform swimming presented in Gazzola et al. [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF] is considered for further validation of the proposed algorithm. A periodic swimming law is defined by fitting the backbone of the fish to a given curve y(x, t) while keeping the backbone length l fish fixed. Let ξ be the arclength over the curvilinear coordinate of the deformed backbone (0 ≤ ξ ≤ l fish ). For points being uniformly distributed with ∆ξ = l fish /(N -1), over the backbone, y is given by y(x, t) = a(x) sin(2π(x/λ + f t)) [START_REF] Belkhiri | Modélisation dynamique de la locomotion compliante : Application au vol battant bio-inspiré de l'insecte[END_REF] where λ is the wavelength, f represents the frequency of the backbone and the envelope a(x) is given by

a(x) = a 0 + a 1 x + a 2 x 2 (39)
where x is defined by inverting the arclength integral, i.e., ∆x = ∆ξ/ √ 1 + (∂y/∂x) 2 . The wavelength of the fish is defined in accordance with the geometry of the backbone in the Cartesian coordinate. The pointwise curvature of the backbone is needed to be able to use the geometrically exact theory of nonlinear beams. One must switch from the Cartesian system to the curvature, thus second derivative of Eq. ( 38) will lead to

k(ξ, t) = (2a 2 -(2π/λ) 2 a(ξ)) sin(2π(ξ/λ + f t)) + (4π(a 1 + 2a 2 ξ)/λ) cos(2π(ξ/λ + f t)) (40)
where a(ξ) = a 0 + a 1 ξ + a 2 ξ 2 . The parameters used by Kern and Koumoutsakos [START_REF] Kern | Simulations of optimized anguilliform swimming[END_REF] and Gazzola et al. [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF] for the kinematics of the fish are as follows; λ = 1, f = 1, a 2 = 0, a 1 = 0.125/(1 + c), a 0 = 0.125c/(1 + c) and c = 0.03125. The buoyancy is equal to zero, i.e., ρ b = ρ f . The viscosity of the fluid is set to ν = 1.4 × 10 -4 resulting in a Reynolds number approximately Re ≈ 3800, with an asymptotic mean velocity U forward ≈ 0.52. The simulations of Gazzola et al. [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF] are carried out on a rectangular domain (x, y) ∈ [0, 8l fish ] × [0, 4l fish ] with resolution of 4096 × 2048 and penalization parameter η = 10 -4 . We are performing our simulations on a rectangular domain (x, y) ∈ [0, 10l fish ] × [0, 5l fish ] by imposing a penalization parameter inside the body equal to η = 10 -3 with resolution of 2048 × 1024 and 1024 × 512 and ∆t = 10 -3 . The centroid of the fish is initially positioned at x cg = 0.9L x and y cg = 0.5L y . The forward velocities of the center of the mass computed with different methods/parameters are compared in Fig. 14. We impose two degrees of freedom fixing the angular velocity equal to zero. The simulations start with the body and fluid at rest. The motion of the fish is initialized by gradually increasing the amplitude of the backbone through a sinusoidal function (see Fig. 13) from zero to its designated value during the first period T in the reference simulations, i.e., [START_REF] Kern | Simulations of optimized anguilliform swimming[END_REF] and [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF]. Here we are not considering this and start by a sudden movement given by Eq. [START_REF] Belkhiri | Modélisation dynamique de la locomotion compliante : Application au vol battant bio-inspiré de l'insecte[END_REF]. This is the reason why a deviation from the reference solution can be seen in the first period. This deviation will continue systematically until the asymptotic velocity is reached at t = 7.

The reference simulation of Kern and Koumoutsakos [START_REF] Kern | Simulations of optimized anguilliform swimming[END_REF] is based on a body fitted finite volume method which is first-order in time and second-order in space. The Navier-Stokes equations were solved using the commercial package STAR-CD which computes flows on arbitrary Lagrangian-Eulerian grids. The solution of Newtons equations of motion and the resulting grid deformation and motion are implemented in user defined subroutines linked to STAR-CD. The implemented explicit coupling procedure is a staggered integration algorithm proposed by Farhat and Lesoinne [START_REF] Farhat | Two efficient staggered procedures for the serial and parallel solution of three-dimensional nonlinear transient aeroelastic problems[END_REF]. The simulation of Gazzola et al. [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF] is done by a code based on a vortex particle method coupled with Brinkman penalization which handles arbitrarily deforming bodies and specially the corresponding divergent velocity field inside the body. A projection method is used by Gazzola et al. [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF], the resulting Poisson equations for rotational (solenoidal) and potential (divergent) components of the velocity fields are solved by an unbounded domain, FFT based solver in Cartesian grids. A second and fourth order finite difference discretization is used for all other spatial derivatives. The time step is adapted by a Lagrangian CFL condition. The difference on the final forward velocity of the fish reported by Gazzola et al. [START_REF] Gazzola | Simulations of single and multiple swimmers with nondivergence free deforming geometries[END_REF] by taking into account the divergence of the velocity field inside the fish due to deformation of the fish is visible in Fig. 14. Even though the average divergence over the fish volume is zero (i.e. the volume is conserved) locally inside the fish the velocity field is not divergence free. We are not dealing with this issue in this paper because for the moment we are not interested in swimmers with high thickness. In our simulations a grid independent solution is obtained with 2048 × 1024 grid points. The difference of two simulations with 2048 × 1024 and 1024 × 512 grid points can be seen in Fig. 14. Filtering of the hydrodynamic coefficients is necessary to prevent the simulation from divergence and non-physical results. We are using a second-order exponential filtering (31) instead of the firstorder filtering used in [START_REF] Kern | Simulations of optimized anguilliform swimming[END_REF]. This process is like adding a damper to the system therefore a proper value for δ must be chosen via numerical tests for obtaining physical results. We propose values in the range of δ ∈ [0.01, 0.001] for fluid/solid interaction problems, however this can also depend on the manner of non-dimensionalization of the forces. In Fig. 14 the effect of filtering with two filter parameter, i.e., δ = 0.001 and δ = 0.05, can be seen. Filtering with a small filter parameter δ = 0.001 will be more stable but instead will lead to smaller values in the terminal velocity and also a smaller amplitude in its oscillations. In Fig. 15 two snapshots of vorticity isolines at t = 1 and t = 9 with the aforementioned parameters are given.

Application and results

We attempt to propose an efficient law for rotation of an anguilliform swimmer. Fish maneuvering for track- ing a fixed goal is done by adding a constant curvature k offset (θ des , t) all along the fish's backbone ξ ∈ [0, l fish ], to the primary propulsion mode, i.e., k 3 = k(ξ, t) + k offset . However the change of the added curvature k offset via Eq. ( 42) must be gradually, i.e., O(∆t) to perform a physically reasonable rotation. For the fish in forward gait k offset is set equal to zero. For applying a rotation, first a desired curvature k des must be evaluated by the following relation, k des (θ des ) = { -sgn(θ des ) k max |θ des | ≥ θ limit -sgn(θ des ) k max ( θ des θ limit ) 2 else (41) where sgn represents the sign function, i.e., sgn(θ des ) = θ des /|θ des |. For a schematic representation of θ des see Fig. 16. A desired angle θ des should be calculated in each time step according to the position and direction of the head by considering the target. Then by using Eq. ( 41) a desired curvature k des must be found. After that k offset will be evaluated with the following relation,

t 0 0.5 1 -6 -6 -4 -4 -2 -2 0
k n+1 offset (k des ) = { k n offset + ∆k k < k desired k n offset -∆k else (42)
where ∆k = ∆t π/T , k offset will be added to the backbone curvature for performing a rotation. We are using k max = π which is equivalent to turning with a curvature adapted to a semicircle. As in Bergmann and Iollo [START_REF] Bergmann | Modeling and simulation of fish-like swimming[END_REF] we are using θ limit = π/4 . The time derivative of the curvature dk/dt is needed in Eq. ( 21) for velocity calculation and can be calculated numerically. To show the performance of the proposed method a test case of food finding with the above rotation law is performed. The domain size is (x, y) ∈ [0, 5l fish ] × [0, 5l fish ], the resolution is set to 1024 × 1024, the penalization parameter η = 10 initial position of the head (x 0 , y 0 ) = (0.1L x , 0.5L y ) and initial angle of the head θ 0 = 0. Fig. 17 shows snapshots of vorticity isolines obtained during a simulation of swimming fish for finding a food which is located at (x f , y f ) = (0.9L x , 0.5L y ). At t = 0 the fish and the surrounding flow are in rest. After reaching the vicinity (r food = 0.5l fish ) of the food the curvature of the backbone, given by Eq. [START_REF] Sotiropoulos | Immersed boundary methods for simulating fluid-structure interaction[END_REF], will tend to zero by multiplying it with the following function,

C(t) = t f -t t f -t i + 1 2π sin(2π t -t i t f -t i ) , t ∈ [t i , t f ] (43)
which is the mirror of the function presented in Fig. 13, with t i = t reached , t f = t reached +T for gradually decreasing the curvature of the backbone during one period. With the proposed law for rotation which is added to the primary propulsion mode, the fish can perform a 180 o change in swimming direction within an area of about its length.

Conclusion

In this paper an efficient algorithm for simulation of deformable bodies interacting with two-dimensional incompressible flows is proposed. By using a uniform Cartesian grid a new fourth-order direct solver for the solution of the Poisson equation is presented. For introducing a deformable body in fluid flow, the volume penalization method is applied to the solution of the Navier-Stokes equations as a forcing term. Even if the penalization method is first-order in space, an important advantage of this method is that the evaluation of the hydrodynamic coefficients is straightforward. Proper denoising of the hydrodynamic coefficients is vital for dealing with fluid/solid interaction problems via the volume penalization method. An efficient law for curvature control of an anguilliform swimmer looking for food is proposed which is based on geometrically exact theory of nonlinear beams. Validation of the developed method shows the efficiency and expected accuracy of the algorithm for fish-like swimming control and also for a variety of fluid/solid interaction problems. Some perspectives for future developments are to develop high-order immersed boundary methods, adding a multiresolution analyses to the algorithm for grid adaptation, enhancement of rotation law, parallelization and extension to three dimensions. The FORTRAN code is open source and is also accessible upon request [START_REF]The code is developed in FORTRAN and is accessible for all by sending a mail to: ghaffari@L3m.univ-mrs[END_REF].

2. 1 . 3 .

 13 Fourth-order fast Poisson solver In solving the incompressible Navier-Stokes equations an elliptic Poisson equation is frequently encountered which is the most time consuming part of the algorithm. The common case is the pressure Poisson equation normally used with homogeneous Neumann boundary conditions,

Figure 4 :

 4 Figure 4: Interpolated smoothed mask function χ (colored isolines) and the velocity components.

4 .

 4 ) (b) (Just for the Eel) Cover the shape by a Lagrangian structured grid & calculate velocities of each point with Eq. (25) (c) Compute the mask χ(i, j) and smooth it by Eq. (26) (d) Compute the moment of inertia J around the reference point (e) Compute the velocity components of the body u p (i, j), v p (i, j) on the Eulerian grid (Lagrange → Euler) 3. Time integration of flow field via RK4 (a) ω * 0 = ω n For i=1,2,3 (b) Compute k i (ω * i-1 ) from Eq. (8) (c) ω * i = ω n + α i ∆t k i (d) Solve Eq. (2); -∇ 2 ψ i = ω * i for updating (u, v) End For (e) Compute k 4 (ω * 3 ) from Eq. (8) (f) Update vorticity from Eq. (7) (g) Solve Eq. (2); -∇ 2 ψ = ω n+1 Solve for the body dynamics (a) Compute the hydrodynamic coefficients of the body from Eqs. (29) and (30) (b) Denoise the coefficients by Eq. (31) (c) Compute the displacements from Eq. (34) (d) Compute the rotation from Eq. (36) 5. Write necessary data to file 6. If T < T end , Go to step 2 7. End

Figure 5 :

 5 Figure 5: Flowchart of the fluid/solid interaction (FSI) algorithm.

Figure 6 :

 6 Figure 6: Schematic representation of the penalized unit square domain for modeling the Taylor-Couette flow with volume penalization method (χ = 0 represents the fluid domain and χ = 1 the solid domains, respectively). The radius of the inner cylinder is R 1 = 0.15 and that of the outer cylinder is R 2 = 0.4. The angular velocity of the inner cylinder is Ω 1 = 0.2 and that of the outer is equal to zero, ν = 0.01 and Ta = R 1 (Ω 2 -Ω 1 ) 2 (R 2 -R 1 ) 3 ν -2 ≈ 1.

= 13 Figure 8 :

 138 Figure 8: Vorticity isolines (dashed lines are used for negative values) of a falling cylinder in quiescent fluid, g = -9.81 m/s 2 , ρ b /ρ f = 1.01, D = 0.005 m, (x, y) ∈ [0 , 0.04 m] × [0 , 0.32 m] = [0 , 8 D] × [0 , 64 D] and ν = 8 × 10 -7 m 2 /s and Re ≈ 156.

Figure 14 :

 14 Figure 14: Forward velocity U of a 2D anguilliform swimmer's (λ = f = 1). Solid lines indicate the reference simulations performed by Kern and Koumoutsakos (green) [22] and Gazzola et al. (pink and brown) [33]. Dashed lines represent the results with the proposed algorithm.

9 Figure 15 :

 915 Figure 15: Snapshots of vorticity isolines obtained during a simulation where (x, y) ∈ [0 , 10l fish ] × [0 , 5l fish ] by imposing a penalization parameter inside the body equal to η = 10 -3 with resolution of 2048 × 1024, the Reynolds number is approximately Re ≈ 3800.

= 13 Figure 9 :Figure 11 :

 13911 Figure 9: Vorticity isolines (dashed lines are used for negative values) of a falling cylinder in a slightly perturbed fluid, g = -9.81 m/s 2 , ρ b /ρ f = 1.01, D = 0.005 m, (x, y) ∈ [0 , 0.04 m] × [0 , 0.32 m] = [0 , 8 D] × [0 , 64 D], ν = 8 × 10 -7 m 2 /s and Re ≈ 144.

Figure 12 :

 12 Figure 12: Forces (left) and velocity components (right) of a falling ellipse in the fluttering regime, I * = 0.16, ρ b /ρ f = 1.538, a/b = 1/5 and ν = 0.01.

Figure 16 :

 16 Figure16: Schematic representation of desired angle for curvature control in rotation, θ des = θ food -θ Head is the difference of the angles between head's direction and food's angle (-π < θ des < π), picture from Bergmann and Iollo[START_REF] Bergmann | Modeling and simulation of fish-like swimming[END_REF] with a slight modification.

Exact Explicit 2 Explicit 4 Explicit 6 Cmp. 3-D 4 Cmp. 3-D 6 Cmp. 3-D 8 5D Spec. like Lele 5D Spec. like Kim

					10			10
	3		3
	2.5		2.5		8			8
	5 2	D U	1.5 2	w''(w)	4 6			4 6
					2			2
					0			0
					0	0.5	1	1.5	2	2.5	3
								w
		(a) First derivative

Exact w 2 Exp. 2 nd -Analytical Exp. 4 th -Analytical Compact 4 th Anal. Compact 6 th Anal. Explicit 2 nd -order Explicit 4 th -order Compact 4 th -order Compact 6 th -order

(b) Second derivative Figure

  Figure2: Shape of the fish given by Eq. (24) before deformation.

	0	0.2	0.4	0.6	0.8	1
			X
	3.2
	3
					0.25
					0.24
	Y				0.22
	2.8				0.2 0.18
					0.16
					0.14
					0.12
					0.1
					0.08
	2.6				0.06
					0.04
					0.02
					0
		0.6	0.8	1	1.2
			X
	Figure 3: Lagrangian structured grid covering the fish after defor-
	mation and the corresponding velocity vectors of each point, colored by absolute velocity √ u 2 + v 2 .

  Figure 13: Smooth step function Cr(t) = t ′ -sin(2πt ′ )/(2π), t ∈ [t i , t f ] with t ′ = (tt i )/(t f -t i ), t i = 0, t f = 1for gradually starting the motion proposed by Boyer et al.[START_REF] Boyer | Macro-continuous computed torque algorithm for a three-dimensional Eel-like robot[END_REF]. At t = 1 the left-and right-hand limits are equal for the function Cr and its first Cr ′ and second Cr ′′ derivatives.

	6	Cr	6
	4	Cr '	4
		Cr ''
	2		2
			0

Acknowledgment

We give our thanks to professor Boyer and colleagues for useful discussions and sending their code describing Eel's backbone kinematics in forward gait, to Professor Gazzola and colleagues for sending the needed data for comparisons.

(l) t = 15

Figure 17: Snapshots of vorticity isolines obtained during a simulation of swimming fish for finding a food which is located at (x f , y f ) = (0.9Lx, 0.5Ly). At t = 0 the fish and the surrounding flow are in rest. After reaching the vicinity (r = 0.5l fish ) of the food the curvature of the backbone will tends to zero by Eq. ( 43). The domain of the solution is (x, y) ∈ [0 , 5l fish ] × [0 , 5l fish ], the resolution of the grid 1024 × 1024 and kinematic viscosity equal to ν = 1.4 × 10 -4 .