Michel Bergmann
email: michel.bergmann@inria.fr

Numerical modeling of a self propelled dolphin jump out of water

A computational model is developed to investigate the jump of a self propelled dolphin out of water. This model relies on the Navier-Stokes equations where a fictitious domain approach with the volume penalization method is used for the fluid structure coupling, and the continuous surface force approach is used to model the water-air interface, the later being tracked in a level-set framework. The dolphin geometry is based on freely available data from the literature. While the body deformation is imposed, the leading linear and angular displacements are computed from the Newton's laws. Numerical simulations show that it is necessary to generate large propulsives forces to allow the jump out of water. When the dolphin is out of water, its trajectory follows a pure ballistic one.

Introduction

The study of underwater swimming animals is interesting to understand the possible effects related to the evolution, in particular to make link between the geometry of the body, the unsteady nature of the flow, and the unsteady swimming deformations imposed by muscles. The later can be optimized by using the energy enclosed in incoming vortices, in active (Gopalkrishnan et al. 1994, Liao et al. 2003, Fish & Lauder 2006, Triantafyllou et al. 2002) or passive (Beal et al. 2006) ways, or using vortices possibly generated by neighbor fishes (Bergmann & Iollo 2011, Maertens et al. 2017, Li et al. 2022, Gao & Triantafyllou 2018). Vortex induced control by the caudal tails of vortices generated along the fish body (Zhu et al. 2002) is also another mechanism to improve the swimming efficiency. The drag reduction and thrust generation are linked to the generation of an inverse Von Kármán street, as observed for simple manipulated circular cylinder wake flow (Bergmann et al. 2006).

A detailed comprehension of all these mechanisms may help to design new concept of Autonomous Underwater Vehicles (AUVs). The swimming of a fish has been considered, following the pioneer but controversial work of (Gray 1936), later on called the Gray's paradox, as being a very efficient mean of locomotion. A review from a fluid mechanics point of view of the Gray's paradox can be found in (Bale et al. 2014). Also, this type of motion allows a maneuverability never reached by the previous generation of AUVs, of torpedo type [START_REF] Bozkurttas | Mining Smartness from Nature (CIMTEC[END_REF].

As already mentioned, dolphin's kick swimming is an efficient way to move in water. Human swimmers, especially during competitions, try to mimic the dolphin's kick during the fully immersed stage (von Loebbecke, Mittal, Mark & Hahn 2009, von Loebbecke, Mittal, Fish & Mark 2009b). A comparison of human and cetaceans is given in (von Loebbecke, Mittal, Fish & Mark 2009a).

The swimming of an animal in water is a complex and beautiful phenomenon that is easy to observe in nature, but difficult to model numerically. The complexity comes from the reciprocal coupling between the fluid and the animal's body. This phenomenon is even more complex when the animal swims close to the surface, because it is necessary to take into account the deformation of the water surface and the forces related to the surface tension. The complexity is greatest when the animal crosses the surface, typically when it jumps out of the water. High-performance propulsion together with trajectory control to optimize the position at the early jump are required to jump out of water. The jumping dynamics of a variety of aquatic animals are described in (Chang et al. 2019), and of focus on the study of the archer fish jump can be found in (Mendelson & Techet 2020). Examples of jump spinning (twist) maneuvering can be found in (Fish et al. 2006).

The goal of this paper is thus to study the dolphin jump in a numerical framework. Only few numerical simulations of dolphin like swimming have been reported in the literature, most of them have been cited above. One reason is that it is difficult to define a realistic dolphin geometry. In this study we will use the freely available geometry of the Lagenorhynchus obliquidens dolphin proposed in (Tanaka et al. 2019). The other reason is, as previously mentioned, the complexity of the numerical modeling. Indeed, it is necessary to model the unsteady fluid structure interactions with large interface deformations and possible large density ratios (body and water interface deformations). In this study we consider a fictitious domain approach where the computational mesh, i.e. grid nodes, does not necessarily follow the interfaces. The dolphin body is considered with the Volume Penalization Method (Angot et al. 1999) and the fluid-structure interface is followed in a Lagrangian way, and the water-air interface is tracked in a Eulerian way using a Level set function (Osher & Sethian 1988[START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF]. Even if it is not a dominant force at the dolphin's scale, we consider the surface tension using the Continuous Surface Force Method (CSF) originally introduced in (Brackbill et al. 1992). Other fictitious domain approaches can be used to compute large density ratios flows, an as for instance the one recently proposed in (Nangia, Griffith, Patankar & Bhalla 2019, Nangia, Patankar & Bhalla 2019). For simplicity reasons, no skin effects are considered in this study, even if they can play an important role in passive drag reduction (Pavlov 2006).

The paper is organized as follows. After having introduced the modeling and numerical methods in §2, including the description of the flow configuration, the governing equations and the associated numerical schemes, we will present the numerical results for dolphin jumps out of water in §3. Some conclusions and perspectives are finally given in §4.

Modeling and numerical methods

Flow configuration

A sketch of the flow configuration is given in figure 1. The whole domain is

Ω = Ω + f ∪ Ω - f ∪ Ω s ,
where Ω s is the domain for the structure, i.e. the dolphin, Ω + f is the domain filled with water, and Ω - f is the domain filled with air. The two fluid domains are separated by an interface Γ f , and the structure and the fluid domains are separated by an interface Γ s . The intersection of Γ f and Γ s is noted Σ t = Γ f ∩ Γ s . For three-dimensional problems, Σ t could be a line, usually called the triple line. The external boundary of the domain is noted Γ e . In what follows, the domain Ω is a three-dimensional Cartesian box, and the external boundary is Γ e = Γ top ∪Γ bottom ∪Γ lef t ∪Γ right ∪Γ f ront ∪Γ back , where each boundary is defined as indicated by the subscripts. The density of the dolphin is ρ s , the density and dynamic viscosity of water are ρ + and µ + , and the density and dynamic viscosity of air are ρ -and µ -.

The different domains and interfaces introduced before also depend on time, and are mathematically defined by two level-set (Osher & Sethian 1988[START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF]) scalar functions ψ f and ψ s . Here, the structure is arbitrarily defined by ψ s > 0, the interface Γ s by ψ s = 0 and the fluids (water and air) by ψ s < 0. Similarly, the water is arbitrarily defined by ψ f > 0, the interface Γ f by ψ f = 0, and the air by ψ f < 0.

Governing equations

Without loss of generality, the density ρ ∈ R and dynamic viscosity µ ∈ R of the fluids can be written:

ρ ≡ ρ(ψ f ) = ρ -+ H(ψ f )(ρ + -ρ -), (1)
µ ≡ µ(ψ f ) = µ -+ H(ψ f )(µ + -µ -). ( 2
)
where H is the Heaviside function, i.e. H(x < 0) = 0, H(x > 0) = 1, and we can define H(x = 0) = 0.5. The velocity fields is u =∈ R 3 , p ∈ R is the pressure field and g ∈ R 3 is the gravity acceleration vector. The viscous part of the stress tensor is

D(u) = ∇u + ∇ T u 2 .
In what follows, we note u = (u, v, w), where u, v and w denote the velocity components in the x, y and z directions respectively. The incompressible Navier-Stokes equations for both fluids in domain

Ω + f ∪ Ω - f are: ∂u ∂t + (u • ∇)u = - 1 ρ ∇p + 1 ρ ∇ • 2µD(u) + g in Ω + f ∪ Ω - f , (3)
∇ • u = 0 in Ω + f ∪ Ω - f -, (4)
with initial conditions u(x, t = 0) = u 0 and p(x, t = 0) = p 0 , boundary conditions on the external boundary for velocity u(x, t) and pressure p(x, t) for x ∈ Γ e , boundary conditions on the structure boundary Γ s , and conditions through the bi-fluid interface Γ f . On the fluid-structure interface Γ s , we have:

u(x, t) = u(x, t) on Γ s . (5)
where the velocity u(x, t) will be described later on. Two jump conditions have also to be satisfied throught the bi-fluid interface Γ f . The first kinematic condition is:

[u] = 0 throught Γ f . (6)
The second dynamic jump condition traduces the equilibrium between pressure forces, the friction forces and the surface tension:

[-pI + 2µD(u)] • n = σκn throught Γ f , (7)
where σ est the value of the surface tension, κ is the curvature of the interface Γ f and n is the unit normal to interface Γ f pointed to the air.

Fictitious domain approach

The numerical resolution of the Navier-Stokes equations ( 3) and (4) requires a discretization on a mesh. The difficulty is that the interfaces Γ f and Γ s , and thus the domains Ω + f , Ω - f and Ω s , are time dependent. Different approaches to taken into account interface conditions (5), ( 6) and ( 7) can be envisioned. This first class of methods are based on body-fitted grids. In these methods, degrees of freedom are put on the interfaces, and it is thus possible to impose directly the interface conditions (5), ( 6) and (7). Mesh deformation can be handle with arbitrary Lagrangian-Eulerian (ALE) method. These methods are accurate, but require mesh adaptation and a mesh partitioning for parallel computations.

A second class of methods, adopted in this study, is based on fictitious domain approaches. In these approaches, the interfaces and associated domains, do not covered the same mesh nodes at each time step: the interface can cross a fixed mesh, and no interface markers are thus required. Simple meshes like Cartesian ones can be used. The drawback is that the accuracy at interfaces can be degraded, and an extra work has to be performed to model the interfaces.

The goal of fictitious approaches is to add extra force-like terms in the momentum equations (3) to take into account the conditions (5), ( 6) and ( 7). The condition on the fluid-structure interface is modeled with an extra term s and the two conditions on the bi-fluid interface are modeled with the extra term f .

The system (3), ( 4), ( 5), ( 6) and ( 7) is thus recasted in a system written in the whole domain Ω:

∂u ∂t + (u • ∇)u = - 1 ρ ∇p + 1 ρ ∇ • 2µD(u) + g + s + f in Ω, (8)
∇ • u = 0 in Ω, (9)
The methods used to model the terms s and f have been inspired by the Immersed Boundary Method (IBM) introduced by Peskin (Peskin 1972). In this continuous framework, an infinite force applied on a Dirac support function is regularized around the interface. This is mainly done for numerical reasons that will be explained later on.

In this study, the fluid structure interface is modeled using the Volume Penalization introduced in (Angot et al. 1999) §2.3.1, and the bi-fluide interface with surface tension is modeled using the Continuous Surface Force (CSF) method introduced in (Brackbill et al. 1992) §2.3.2.

2.3.1.

Modeling of the conditions on the fluid-structure interface Γ s Among the most popular fictitious approches are the Immersed Boundary Methods (IBM) originally introduced in (Peskin 1972) and later on used in several studies (Mittal & Iaccarino 2005, Mittal et al. 2008[START_REF] Shirgaonkar | [END_REF]) and the Volume Penalization (VP) Method introduced in (Angot et al. 1999) and used for instance in biolocomotion problems (Bergmann & Iollo 2011, Bergmann, Iollo & Mittal 2014, Bergmann & Iollo 2016). Another approche combining IBM and VP has been developed in (Bergmann, Hovnanian & Iollo 2014). Here, we used the VP method with

s = χ K ( u -u), (10)
where K 1 is the penalization parameter, χ = H(ψ s ) is the characteristic function, and u is the velocity of the body (on the boundary and inside the body). In domain Ω ± we have ψ s < 0 and thus χ = 0, and the Navier-Stokes equations are recovered. In domain Ω s , we have ψ s > 0, and thus χ = 1, and the dominant term in (3) is χ K ( uu) leading to u → K→0 u in Ω s . Here, we chose K = 10 -8 . The computation of u will be detailed in §2.5.

2.3.2.

Modeling of the conditions on the bi-fluid interface Γ f The CSF (Continuum Surface Force) has been originally introduced in (Brackbill et al. 1992) and developed in (Sussman et al. 1994). The CSF is to model the surface tension term as an extra volume force f in the momentum equation. The jump

[-pI + 2µD(u)] • n = σκn throught Γ f , (11)
is considered in an integral form and we have:

f = 1 ρ σκδn. ( 12
)
where σ is the surface tension, κ the curvature of the interface, n the normal to the interface pointing to the air, and δ the Dirac distribution.

The normal n and the curvature κ are:

n = ∇ψ f ∇ψ f , (13)
κ = ∇ • ∇ψ f ∇ψ f . ( 14
)
The problem is relaxed by regularization of the Dirac distribution δ in a narrow band with width

H (ψ f ) =        0 si ψ f < -, 1 2 1 + ψ f + 1 π sin( πψ f ) si |ψ f | ≤ , 1 si ψ f > . (15)
The regularized Dirac distribution is:

δ (ψ f ) = dH (ψ f ) dψ f =    0 si |ψ f | > , 1 2 1 + cos( πψ f ) si |ψ f | ≤ . ( 16
)
The extra term due to the surface tension is thus finally

f = 1 ρ σκδ n. (17)
Finally, due to the continuous framework introduced beforehand, the jump [u] = 0 throught Γ f in equation ( 6) is naturally satisfied.

Interface tracking

The model will be closed with the computations of functions ψ s and ψ f . The definition of these functions is not unique, and due to its hyper-regularity we chose ψ f to be the signed distance function, i.e. ∇ψ f (x) = 1 ∀x. Let d(x) ≥ 0 be the minimal distance from a given point x to the interface Γ f , we have:

     ψ f (x) = 0 if x ∈ Γ f , ψ f (x) = +d if x ∈ Ω + f , ψ f (x) = -d if x ∈ Ω - f . (18)
The same holds for ψ s . The signed distance allows to simplify the computation of the normal and the curvature of the interface, and determine in a more precise way the distance ε defining the narrow band for CSF regularization. These interfaces, represented by functions ψ f and ψ s , can be transported in Eulerian or Lagrangian ways. Due to possible large deformations and break downs, the bi-fluid interface Γ f will be defined with a level set function that is able to track such interface with topology changes. The level set ψ f is thus transported by:

∂ψ f ∂t + w • ∇ψ f = 0 in Ω ± f , (19)
where w can be any desired velocity fields satisfying w(x, t) = u(x, t) ∀x ∈ Γ f . For simplicity reasons we chose w(x, t) = u(x, t) ∀x ∈ Ω, but other choices like the extension velocities can be used (Adalsteinsson & Sethian 1999).

The signed distance behavior of the function ψ f can be lost integrating ( 19), and it is thus necessary time to time to reinitialize ψ f to recover the signed distance function ∇ψ f (x) = 1. One of the most popular algorithm is the first order fast marching (Sethian 1996[START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF], extended later to higher orders (Ahmed et al. 2011). Another class of method was introduced by Sussman (Sussman et al. 1994). The goal is to solve a partial differential equation to recover the property ∇ψ f (x) = 1. This relaxation method is:

∂ψ f ∂τ + sign(ψ 0 f )(|∇ψ f | -1) = 0, (20)
ψ f (x, τ = 0) = ψ 0 f = ψ f (x, t), (21)
for a temporal horizon 0

≤ τ ≤ τ F such that |∇ψ f (x, τ F )| = 1. The sign function is sign(ψ f ) =      1, if ψ f > 0, 0, if ψ f = 0, -1, if ψ f < 0, (22) or a regularized version of it.
In what follows, the dolphin surface Γ s will be approximated by a mesh (i.e. with markers), and it is thus more convenient to use a Lagrangian transport for any mesh point

x h : dx h dt = u h , (23)
where the definition and computation of u will be defined in §2.5.

The signed distance function is recovered computing the minimal distance to the interface:

ψ s (x) = min y∈Γs x -y 2 S(x). ( 24
)
where S(x) denote the sign applied on a point x, with S(x) > 0 inside the body, and S(x) < 0 outside. This sign function can be computed with simple geometric arguments, from the outward normal to the body. The triple line Σ t , defined as being the intersection between Γ s and Γ f , has to be modeled. Indeed, without special treatment, due to non-slip boundary condition on Γ s , the line would never move on Γ s . In this study we choose the macroscopic model introduced by Cox (Cox 1986). A local force proportional to the deviation of the equilibrium angle (Patankar 2003) is added to allow a displacement of the triple line.

Finally, since the flows under consideration are turbulent, a turbulence model introduced in (Vreman 2004) is added.

Dolphin model

The geometry of the Dolphin is built from the stl file given in the supplementary material of (Tanaka et al. 2019). Indeed, the authors have performed a very productive work given lot of material to work with. The geometry is reproduced in Figure 2. The two first figures are the exact geometry given by (Tanaka et al. 2019), and the two last figures are the geometry a little bit remodeled to fit our numerical framework with regular surface mesh.

The body velocity is

u(x, t) = u(x, t) + u θ (x, t) + u(x, t) ∀x ∈ Ω s , (25)
where u(x, t) is the linear velocity, u θ (x, t) is the angular velocity and u(x, t) is the deformation velocity. While the deformation u(x, t) has to be imposed by swimmer muscles, the linear and angular velocities are the results of the loads generated by the fluid on the body, and are computed from the Newton's laws.

In order to avoid adding extra forces and torques due to deformation, a Procrustes analysis is performed. The deformation must not introduce any linear and angular displacements. We thus compute the linear and angular displacements induced by the imposed deformation, and subtract them to the final admissible deformations. This is actually done in Figure 4.

The deformation velocity u(x, t) can be easily computed following surface markers in a Lagrangian way, after having performed the procrustes analysis. To recover u(x, t) values in Ω s , interpolation is performed from values at the boundary Γ s .

The swimming law is defined by a deformation of the swimmer midline, also called the backbone. As we will see later on, the dolphin jump can be decomposed in three steps.

• The first one is an underwater generation of a propulsive force using a periodic swimming with a large amplitude and frequency. • The second one corresponds to the jump out of the water, i.e. when the dolphin is in the air: the dolphin just imposes (active or passive) a constant curvature along the midline.

• The third and last step is another periodic swimming law with a small amplitude when the dolphin is back in water.

We consider a backbone deformation in the plan (0, x, y) of the third figure in figure 2, where point O is at the front head of the dolphin and x positive to the right. We consider moreover that the midline for the steady body is 0 ≤ x ≤ , y = 0, where is the length of the swimmer.

Most of the fishes impose the following periodic swimming law (Barrett et al. 1999):

y(x, t) = a(x) sin(k x -ω t), (26)
with a constant phase velocity c p = ω/k, where k = 2π/λ is the wavenumber, corresponding to wavelength λ, ω is the circular frequency of oscillations, and where a(x) is the envelop defined by a(x) = A/2(c 0 + c 1 (x -1) + c 2 (x 2 -1)). Other kinds of swimming laws can be found in (Smits 2019) and (Han et al. 2020). The maximal tail amplitude A/2 is an important parameter for the swimming (Lighthill 1970). Here, we chose c 0 = 1, c 1 = -0.825 and c 2 = 1.625 for a unit length fish with 0 ≤ x ≤ 1. These parameters are used to mimic a thunniform like swimming. Several midline deformations are presented in Figure 3 for one stroke, with A = 0.33. Instead imposing directly the swimming law (26), it is more simple to impose a curvature κ(s) along the midline with curvilinear abscissa s. The curvature corresponding to the swimming law y(s, t) = a(s) sin(k s -ω t) for the first and third stage of the jump is:

κ(s, t) = ∂ 2 y ∂s 2 1 + ∂y ∂s 2 3 2 . ( 27
)
The general curvature imposed on the midline is:

κ(s, t) = p(t)κ(s, t) + q(t)κ, ( 28
)
where κ is the curvature imposed for the jump, i.e. when the dolphin is in the air, and 0 ≤ p(t) ≤ 1 and 0 ≤ p(t) ≤ 1 are two functions, possibly piecewise defined, allowing to go smoothly from step 1 to step 2 and from step 2 to step 3. This swimming law is apply to the dolphin midline (the backbone). All body surface mesh nodes are displaced with respect to Euler Bernoulli beam deformation, i.e. each orthogonal section to the backbone for the undeformed shape remains orthogonal to the backbone during the deformation. Eight snapshots of the deformed body shape corresponding to (28) over one periodic stroke (p = 1 and q = 0) are presented in Figure 4. Note that a procrustes analysis is performed. Finally, the linear and angular motions are obtained from the Newton's law

         m du dt = F ext , dJω dt = M ext , (29)
where m and J are the mass and inertia matrix of the dolphin, and u and ω denote the linear and angular velocities.

The forces and the torques are computed by

F ext = - Γs T(u, p) n dx + g, (30)
M ext = - Γs r ∧ T(u, p) n dx, ( 31
)
where T(u, p) = -pI + µ(∇u + ∇u T ) is the stress tensor, n is the unit outward vector to Γ s , and r = xx G with x G the dolphin center of mass.

The rotation velocity is given by u θ = ω ∧ r.

Note that other methods can be used to computed the forces and torques in fictitious domain approach when no surface markers are available. For instance, the surface integrals in ( 30) and ( 31) can be transformed in a control volume integrals (Nangia et al. 2017).

Numerical resolution

The Navier-Stokes equations ( 8) and ( 9), and associated equations like transport and reinitialization of the level set function, are discretized in space on an uniform Cartesian mesh using finite differences method. The temporal discretization of the Navier-Stokes equations is based on a fractional step method introduced in (Chorin 1968) and (Temam 1969). The numerical method used in this paper follows (Bergmann & Iollo 2011, Bergmann, Hovnanian & Iollo 2014, Bergmann, Iollo & Mittal 2014, Bergmann & Iollo 2016) and we recall thus only the main steps.

All the differential operators are discretized in space using centered second order schemes, except the convective terms where up-winding is necessary for stability reasons. The convective terms in the momentum equations are discretized with a classical third order finite differences scheme, while the convective term in the level set transport equation ( 19) is discretized with a WENO5 scheme (Liu et al. 1994, Jiang & Shu 1996). A total variation diminishing third order Runge-Kutta (RK3-TVD) scheme is used to discretized equation ( 19) in time. The reinitialization equation ( 20) is solved using the the subcell fix method proposed in [START_REF] Russo | [END_REF] with extensions developed in (Luddens et al. 2015). It is not necessary to perform a reinitialization step at each time step, and we use the criterion introduced in (Luddens et al. 2015) to decide when it is necessary to perform a reinitialization of the level set function.

The integrals over Γ s involved to evaluate the forces and torques are discretized using the body the surface mesh (that is a tessellation of the surface Γ s ). Since the dolphin's surface mesh nodes do not coincide with the fluid mesh, the stress tensor is interpolated from fluid mesh to the body surface mesh.

Numerical validations

The numerical approach presented above has already been validated onto several benchmark problems in our previous publications. Several validations in three dimensions complement numerous validations already performed for two dimensional problems. Our previous validations consider mono-fluid problems, i.e. without bi-fluid interface. A first validation is the sedimentation of a sphere (Bergmann, Iollo & Mittal 2014) for several sphere diameters and fluid viscosities. For bio-inspiration problems involving large deformations, another validation has been performed in (Bergmann & Iollo 2016) onto an eel swimming originally proposed by (Kern & Koumoutsakos 2006) and studied numerically thereafter in (Bhalla et al. 2013) and (Bhalla et al. 2014). For both validation problems, the numerical results are in good agreements with reference ones.

In what follow, a three-dimensional bi-fluid validation is performed for the entry of a sphere in water. This problem has been studied experimentally (Aristoff et al. 2010) and (Truscott et al. 2014), and numerically with the commercial software ABAQUS in (Ahmadzadeh et al. 2014).

A sphere with diameter D = 0.0254 m impacts the water surface at V 0 = 2.17 m/s. The density of the sphere is ρ s = 1.14ρ + . The densities and dynamic viscosities of the fluids are ρ + = 1000 kg/m 3 and µ + = 10 -3 P a.s for the water, and ρ -= 1.2 kg/m 3 and 1.87 • 10 -5 P a.s for the air. The surface tension is σ = 0.0728 N/m.

The computational domain under consideration is a Cartesian box with x ∈ [-6, 6] cm, y ∈ [-20, 10] cm and z ∈ [-6, 6] cm. At rest, the air-water interface Γ f is located at y I = 0 cm, and the initial position of the center of mass of the sphere is

x G = (x G , y G , z G ) = (0, D/2, 0) m.
The computational domain is discretized by 240 × 600 × 240 ≈ 35 • 10 6 cells, and the step size is h = ∆x = ∆y = ∆z = 0.05 cm. The time step is defined with a CFL condition equal to 0.3.

Non-slip boundary conditions are applied on all external boundaries. The flow is initially at rest, i.e. u(0, x) = 0.

Figure 5 presents the image sequences showing the water-entry cavity formed by the sphere until t = 65 ms where the Pinch-Off, that is characterized by a cavity closure event, occurs. Results are in good agreement with experimental results (Aristoff et al. 2010) where the Pinch-Off is observed at t = 63.5 ms, and with the numerical results (Ahmadzadeh et al. 2014) where the Pinch-Off is observed at t = 64.9 ms.

A time history comparison of experimental (Aristoff et al. 2010), numerical ABAQUS simulation (Ahmadzadeh et al. 2014) and present numerical results of sphere center depth are shown in Figure 6. The present numerical results are in good agreements with reference results, especially with the experimental ones.

Numerical results

Numerical configuration and parameters

The densities and dynamic viscosities of the fluids are ρ + = 1000 kg/m 3 and µ + = 10 -3 P a.s for the water, and ρ -= 1.2 kg/m 3 and 1.87 • 10 -5 P a.s for the air. The surface tension is σ = 0.0728 N/m.

The dolphin geometry is extracted from the supplementary material from (Tanaka et al. 2019) and the geometry is reproduced in Figure 2. Here, the length of the dolphin is imposed to be = 3 m, leading to total mass m = 316 kg with ρ S = 1000 kg/m 3 . The computational domain under consideration is x ∈ [-12, 4] m, y ∈ [-3, 9] m and z ∈ [-2, 2] m. At rest, the air-water interface Γ f is located at y I = 4.5 m, and the center of mass of the dolphin is x G = (x G , y G , z G ) = (0, 0, 0) m. The initial orientation of the dolphin is θ = (θ x , θ y , θ z ) = (0, 0, θ 0 ). While the initial angular velocity is ω = 0, the initial dolphin linear velocity is u = (u G , v G , w G ) = (0, v 0 , 0), where v 0 can be a positive constant to decrease the distance to reach a sufficient velocity for a real jump, and thus to decrease the size of the computational domain.

The velocity boundary condition on the external boundaries Γ ext are the following. On the side boundaries (left, right, front, back), periodic conditions are imposed, u f ront = u back and u lef t = u right . No slip boundary conditions are imposed on the bottom u bottom = (0, 0, 0), and streamline type of boundary conditions are imposed on the top, i.e. v top = 0,

∂u top ∂y = ∂w top ∂y = 0.
The dolphin surface is discretized with a regular mesh composed by 200 × 200 = 40, 000 cells, see two bottom figures in 2. Even if large skewness of mesh near fins and tails are observed, the mesh is a tessellation of the body surface on which integrals ( 30) and ( 31) can be discretized. This mesh is similar to those used in (Bergmann & Iollo 2016, Bergmann, Iollo & Mittal 2014) for carangiform and thunniform swimmers.

The computational domain is discretized by 640 × 480 × 160 ≈ 50 • 10 6 cells, and the step size is h = ∆x = ∆y = ∆z = 0.025 m. The time step is defined with a CFL condition equal to 0.3.

Due to the large dolphin deformations, boundary layers are detached and it is thus not necessary to use a costly boundary layer mesh. The turbulence model introduced in (Vreman 2004) is used. Other turbulence model may be used to study the laminarturbulence transition [START_REF] Riedeberger | High Performance Computing in Science and Engineering '11[END_REF].

Artificial waves at the air-water interface are generated by a Kelvin-Helmholtz instability using initial conditions u(x, 0) = 0 m/s in the water (x ∈ Ω + f ) and u(x, 0) = 0 m/s in the air (x ∈ Ω - f ). From the initial water surface position at y I = 4.5 m, we define two other locations: the first one is over water, y + I = y I + 0.6 m, and the second one is under water, y - I = y I -0.8 m. We denote by T = 0.2 s a transition period over that p goes linearly from 0 to 1 or from 1 to 0.

The dolphin body deformation is defined according to the following swimming law.

(i) In order to produce enough propulsive effect to get out of the water starting only 4.5 m under the surface, we imposed a sinusoidal law with large amplitude and frequency, i.e. A = 1.5 m and f = ω/(2π) = 4 Hz , with initial linear transition with p(0) = 0 and p(T ) = 1. This swimming law is performed until the dolphin center of mass y G reaches y + I at time t + I , and we thus have p(t) = 1 ∀t ∈ [T, t + I ]. The periodic swimming law is then linearly damped to reach p(t + I + T ) = 0. (ii) When the periodic swimming is damped from t = t + I , a mean curvature κ = 1/R with R = 2 m is progressively added to the whole midline. The smoothing function q depends on the vertical position of the center of mass such that maximal curvature is obtained when the dolphin is in the air. When y G = y I , at the end of the jump, the function q is damped and the curvature tends to zero.

(iii) When the dolphin re-enter into the water due to the gravity, more precisely when y G < y - I at t = t - I , another periodic swimming with A = 1 m and f = 2 Hz is activated with a linear transition with p(t - I ) = 0 and p(t) = 1 for t ≥ t - I + T .

Mathematically, the function p depends on some temporal intervals defined by:

p(t) =                  t/T if 0 ≤ t < T, 1 if T ≤ t < t + I , (t + I + T -t)/T if t + I ≤ t < t + I + T, 0 if t + I + T ≤ t < t - I , (t -t - I )/T if t - I ≤ t < t - I + T, 1 if t ≥ t - I + T, (32)
and the function q is linked to the position y G of the dolphin center of mass for t ≤ t - I :

q(t) = 1 max 1, 1 (10 -4 + minmod{0.7(t -t + I ), 1} (y G (t) -y - I )) 2 , ( 33
)
where

minmod{a, b} =      a if |a| < |b|, ab > 0 b if |b| < |a|, ab > 0 0 otherwise, (34)
This handcrafted law has been developed from several numerical experiments to mimic as best as possible the dolphin jump.

Numerical results

A numerical simulation with θ 0 = 60 • and V = 4 m/s is performed using 256 processors. The initial air velocity is u(x, 0) = -5 m/s for (x) ∈ Ω - f . The temporal evolution of functions p and q are given in figure 7(b), where the different time t - I , t I , t + I and the other ones necessary to defined functions p and q are plotted using vertical lines. These lines depend on the vertical position y G of the dolphin center of mass following (32) and (33), see figure 7(a). The periods with p(t) ≥ 0 are highlighted in blue, and the period where q(t) ≥ 0 is highlighted in red. Blank zones thus denote undeformed body shape.

Figure 8 shows the temporal evolution of the linear and angular position of the dolphin, the linear velocity, and the forces and torques exerted by the fluid on the body. The vertical velocity of the dolphin goes from v G = 4 m/s at t = 0 s to v G = 6.4 m/s at t = t + I , see figure 8(b). Note that we have chosen an initial velocity v G > 0 to make use of a smaller computational domain. After this initial acceleration, the obtained velocity allows the dolphin to go out of water (figure 8(a)). When the dolphin is the air (in the red zone in figure 8), the dolphin vertical velocity v G and the angle θ z start to decrease due to gravity effect and the curvature. After re-entry into water, a new sinusoidal swimming law is performed and the the dolphin swims in the water with a position θ z ≈ -10 • along which a slight acceleration is observed.

Eight snapshots {S i } i=1,..8 uniformly extracted during the actual numerical simulation are presented in figures 9, 10 and 11. These snapshots correspond to {t i } i=1,..8 with t i = 0.45 (i -1) + 0.6 s and are plotted in figure 8 using vertical black lines. While figures 9 and 10 present the dolphin jump visualisation of the water surface from different views, figure 11 presents snapshots of the Q-crietrion.

The two first snapshots S 1 and S 2 correspond to the first periodic swimming law for propulsion. Snapshots S 3 and S 4 correspond to a period in the air where maximal curvature κ is imposed. The fifth snapshot, S 5 is extracted for almost undeformed body shape when re-entering in water. This deformation corresponds to q κ = 0.007 m -1 and p = 0. Finally, snapshots S 6 , S 7 and S 8 correspond to the second periodic swimming law when the dolphin is back in water. In water, the vortices generated by the self propelled dolphin are similar to whose observed for most of fish like swimmers (Zhu et al. 2002, Triantafyllou et al. 2000, Liao et al. 2003, Bergmann & Iollo 2016). Other kind of vortices are generated in the shear layer at the water-air interface. Indeed, the initial water and air velocity are respectively u = 0 m/s and u = -5 m/s.

Trajectories of the dolphin center of mass corresponding to different initial velocity v 0 and angle θ 0 are presented in figure 12. As expected, the highest jumps correspond to highest initial velocities. The correlation between initial velocity and the jump can not be easily determined. Indeed, the periodic swimming law is not stopped at the (35) where (x 0 , y 0 ) is initial projectile position, α is the initial trajectory orientation at (x 0 , y 0 ), v 0 is the initial velocity along the α direction at (x 0 , y 0 ), g = 9.81 m/s 2 is the gravity acceleration and t is the time. Figure 13 shows the real dolphin trajectory and several ballistic trajectories corresponding to different initial positions (x 0 , y 0 ). For y 0 = y I at t = 0.944 s, i.e. when the dolphin center of mass reach the water surface, large deviations between the dolphin and ballistic trajectories are observed. This is due to the fact that the dolphin is still accelerating (see figure 8). For y 0 = y + I = y I + 0.6 m at t = 1.117 s, the periodic swimming law is progressively stopped. A part of the dolphin is still in water, the acceleration is positive, and the deviation between the real trajectory and the ballistic trajectory is lower. Finally, for y 0 = y I + 0.8 m at t = 1.23 s the acceleration tends slowly to zero and the ballistic trajectory almost perfectly matches the real dolphin trajectory, until the dolphin goes back to water at y G = y I . This is an obvious observation since the forces exerted by the air on the dolphin (F X ≈ 6 N , see figure 8) are quite low comparing to the inertial effect of a 318 kg dolphin with velocity around 6 m/s.

Conclusions

In this paper we have developed a numerical model of the dolphin jump. This model deals with complex physical phenomena such as fluid-structure interactions and liquid interface involving large deformations with possible topological changes (bubbles and droplets). The numerical model relies on fictitious domain approaches where the mesh used to solve the fluid does not follow the interfaces. The CSF method is used to model the water-air interface, and the Volume Penalization method is used to take into account the dolphin's body. The numerical results are qualitatively in good agreements with the dolphin jumps observed in nature. The dolphin has to generate large forces to reach a velocity threshold necessary for the jump. The jump is similar to a ballistic projectile where the resistive forces are neglected. This hypothesis is justified since the computed resistive forces are only around 6 Newtons for a 320 kg dolphin with a velocity around 6 m/s. We are now working on the optimization of the initial periodic swimming law to generate maximal acceleration to produce the highest jump.

Figure 1 .

 1 Figure 1. Two-dimensional sketch of the general three-dimensional flow configuration.

Figure 2 .

 2 Figure 2. Geometry of the Dolphin built with the stl file from supplementary material of Tanaka et al. (2019). The two bottom figures show the mesh used in this study.

Figure 3 .

 3 Figure 3. Periodic swimming law: backbone deformation with A = 0.33.

Figure 4 .

 4 Figure 4. Temporal evolution of the dolphin shape over one stroke period τ . Should be read from left to right, top to bottom. t = 0, T = τ /8, T = 2τ /8, T = 3τ /8, T = 4τ /8, T = 5τ /8, T = 6τ /8, T = 7τ /8

Figure 5 .

 5 Figure 5. Image sequences showing the water-entry cavity formed by the sphere. Times since the sphere center passed the free surface (t = 0) are shown until the Pinch-Off at t = 65 ms.

Figure 6 .

 6 Figure 6. Temporal evolution of the vertical position y G of the center of mass of the sphere.

Figure 7 .

 7 Figure 7. Temporal evolution of functions p and q and the vertical position of the dolphin center of mass y G .

Figure 8 .

 8 Figure 8. Temporal evolutions of dolphin linear and angular position and velocity, and aerodynamic forces and torques.

  (a) t 1 = 0.60 s. (b) t 2 = 1.05 s. (c) t 3 = 1.50 s. (d) t 4 = 1.95 s.(e) t 5 = 2.40 s.(f) t 6 = 2.85 s.(g) t 7 = 3.30 s.(h) t 8 = 3.75 s.

Figure 9 .

 9 Figure 9. Snapshots of the dolphin with water surface for θ 0 = 60 • and v 0 = 4 m/s. Top view.

Figure 10 .

 10 Figure 10. Snapshots of the dolphin with water surface for θ 0 = 60 • and v 0 = 4 m/s. Front view.

Figure 11 .

 11 Figure 11. Snapshots of the dolphin with water surface for θ 0 = 60 • and v 0 = 4 m/s. Front view.

Figure 12 .

 12 Figure 12. Comparison of the dolphin jump with different initial conditions v 0 and θ 0 . Jump is from right to left.

Figure 13 .

 13 Figure 13. Comparison of the dolphin jump (v 0 = 4 m/s and θ 0 = 70 • ) with ballistic trajectories. Jump is from right to left.

Acknowledgements

Computer time for this study was provided by the computing facilities MCIA (Mésocentre de Calcul Intensif Aquitain) of the Université de Bordeaux and of the Université de Pau et des Pays de l'Adour.