Modeling and Simulation of Fish-Like Swimming in a Straight-Line Swimming State Using Immersed Boundary Method

Abstract

A self-propelled swimming fish model is established, which can reflect the interaction between fish movement, internal force generated by muscle contraction, and the external force provided by fluid. Using finite element immersed boundary method combined with traditional feedback force method, the self-propelled swimming fish is numerically simulated. Firstly, a self-induced vibration of a cantilever beam immersed in a fluid is one of the benchmarks of fluid-structure interaction, which is used to verify the validity of the numerical method. Secondly, start and cruise process of a single swimming fish in a straight-line swimming state is simulated and analysis of the flow characteristics and fish body movement features is done. The results reveal that the fish gain energy from flow field by the conversion of “C” type and “S” type of fish body.

1. Introduction

Study on the interaction between fluid and floating or flying creatures is usually applied to the bionics and reveals the mechanism of interaction between biological motion and fluid. Investigation of hydrodynamic characteristics of fish swimming can be traced to the invention of ancient boat and oars. In ancient times, people saying “a scull three pulp,” “namely, the scull efficiency is about three times the oar,” shows the efficiency of fish swinging. So, the modeling and simulation of fish-like swimming are of interest in life sciences as well as in engineering applications. More generally, many other technological applications ranging from the sedimentation of granular flows to the simulation of a maneuvering aircraft can benefit from the efficient modeling of the unsteady interaction between the flow and moving objects.

Among the first mathematical results concerning fish-like swimming were those reported by Lighthill [1, 2]. Since then, many other papers were dedicated to the study of inviscid flows of fish swimming [3]. A more recent fish-like swimming example is presented in [4–7]. But the traditional numerical method for biological fluid-structure interaction (FSI) problems has some limitations in large deformation of biological tissue, complicated structure, and calculation efficiency. How to establish an ideal numerical approach to deal with biological tissue in fluid-structure interaction problems has become the focus of researchers.

Motivated by these facts, in this paper we study fish-like swimming in the two-dimensional low Reynolds number case. Even though the experimental results cited are relevant to three-dimensional flows at a Reynolds number such that the flow is in the turbulent regime [8–11], the two-dimensional laminar case is of interest since numerical results for viscous flows and considering the energy provided by muscle in a fish-like swimming model are scarce. Also, the self-propelled fish swimming includes several fields of coupling, such as the fish body kinematics, the fluid dynamics, and structure dynamics as well as the internal force fields provided by muscle contraction, which is a complex FSI problem. Some researchers define the swimming form of the fish [12], but the freedom of swimming is limited. Others [13, 14] use an inverse approach; namely, they start with the measured kinematics and then estimate the hydrodynamic forces on the body and then the muscle forces required to generate the prescribed kinematics. Animal movements result from a complex balance of many different forces. Muscles produce force to move the body; the body has inertial, elastic, and damping properties that may aid or oppose the muscle force; and the environment produces reaction forces back on the body. So, study on hydrodynamic characteristics of fish swimming needs to establish an ideal fish swimming physical model, and the influence of key factors on hydrodynamic characteristics of fish swimming includes fish body shape, internal structure of fish, the energy produced by muscles, and external fluid. This is the main innovation in this paper due to the fact that the force produced by muscles is considered.

Handling well large deformation structure and/or moving boundaries is a major challenge in simulations of real natural flows, which requires creative approaches for mesh generation and boundary condition implementation. Many methods, including dynamical grid [15, 16], immersed boundary methods [17, 18], level set methods [19, 20], and penalization methods [21, 22] have been developed and successfully applied to specifically address this challenge. In particular, major advances in immersed boundary methods, such as the concept of grid generation and the idea of modeling the presence of a solid body by appropriate volume forces, as well as discretizing with a single, fixed, nonboundary conforming mesh system, which are of interest in this work, have made it possible to efficiently study flow problems that are far more complicated than those that traditional computational fluid dynamics methods could handle in the past [23–26].

In this paper, the fish-like swimming as the research object, using finite element immersed boundary method combined with the feedback force method [27–30], a flexible 2D model of freely self-propelled swimming fish that interacts between internal forces, fish body movement, and fluid, is established, and the fish is divided into three sections, namely, head, tail, and tail fin. This model of fish presented here includes an actuated, elastic body, based on that of a 2D fish swimming in a fluid domain. The motion of the body emerges as a balance between internal muscular force and external fluid forces. The hydrodynamic performance of a single swimming fish and fish schooling is numerically simulated.

2. Continuous Formulation of Immersed Boundary Method

In the IB formulation, the fluid is described by the Navier-Stokes equations in the Eulerian form and the structure is described by a viscous hyperelastic model in the Lagrangian form. A physical region Ω contains a fluid and an immersed structure. When the structure moves in the fluid region, let β identify the reference region of the elastic structure and let ∂β denote the boundary of ∇ _s ·P/(N/m³). The physical region of the structure evolves as D(s,t) at time t and the boundary is ∂β _t ; see Figure 1. The coordinate system is different for the two physical spaces; let variable x = (x₁,x₂,…)∈Ω⊂ℜ ^d denote Cartesian coordinates of the fluid in the Eulerian form and let vector variable s = (s₁,s₂,…)∈β⊂ℜ ^d denote Lagrangian positions of the material particles. Let X = X(s,t)∈Ω denote the Cartesian position of material point s at time t. The X():s→X is a mapping function.

Figure 1:

Lagrange and Euler coordinate systems.

In the immersed boundary method formulation, an external force is added to the N-S equation, and this reflects the influence of solid on fluid. The density of fluid is ρ _f and the elastic structure is ρ _s . Assume for simplicity that ρ _f = ρ _s = ρ, so the Navier-Stokes (N-S) equations in Ω are

$\begin{array}{c}\rho \dot{\mathbf{u}}=\nabla ·\mathbf{\sigma }\left(\mathbf{x},t\right)=\nabla ·{\mathbf{\sigma }}_{\mathbf{f}}\left(\mathbf{x},t\right)+{\mathbf{f}}_{\mathbf{s}}\left(\mathbf{x},t\right),\\ \nabla ·\mathbf{u}=\mathrm{0},\end{array}$

(1)

where f _s (x,t) = ∇·σ _s (x,t) is the Eulerian force density that is the action of the elastic structure on the fluid. With the no-slip boundary conditions of solid, the action of the fluid on elastic structure can be written as

$\frac{\partial \mathbf{X}(\mathbf{s},t)}{\partial t}=\mathbf{u}\left(\mathbf{X}\left(\mathbf{s},t\right),t\right).$

(2)

Equations (1)–(2) are the control equations of immersed boundary method and the key of this method is the information transfer between Lagrange variables and Euler variables. The traditional method is to use the Dirac delta approximation functions as δ(x-X(s,t)). We consider an arbitrary region V_δ⊂Ω and we have

${{\displaystyle \int }}_{V_{\delta }}\delta (\mathbf{x}-\mathbf{X}(\mathbf{s},t))d\mathbf{x}=\{\begin{array}{cc}1\hfill & \text{if}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathbf{X}(\mathbf{s},t)\in V_{\delta }\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}$

(3)

In this paper, the four-Cartesian-point two-dimensional form of δ(x-X(s,t)) is used and can be defined as

$\begin{array}{c}{\delta }_h\left(\mathbf{x}\right)=\frac{\mathrm{1}}{h^{\mathrm{2}}}\varphi \left(\frac{x}{h}\right)\varphi \left(\frac{y}{h}\right),\\ {\varphi }_{\mathrm{4}}\left(r\right)=\{\begin{array}{cc}\frac{\mathrm{1}}{\mathrm{8}\left(\mathrm{3}-\mathrm{2}\left|r\right|+\sqrt{\mathrm{1}+\mathrm{4}\left|r\right|-\mathrm{4}{r}^{\mathrm{2}}}\right)},\hfill & \mathrm{0}\le \left|r\right|<\mathrm{1},\hfill \\ \frac{\mathrm{1}}{\mathrm{8}\left(\mathrm{5}-\mathrm{2}\left|r\right|+\sqrt{-\mathrm{7}+\mathrm{12}\left|r\right|-\mathrm{4}{r}^{\mathrm{2}}}\right)},\hfill & \mathrm{1}\le \left|r\right|<\mathrm{2},\hfill \\ \mathrm{0},\hfill & \mathrm{2}\le \left|r\right|,\hfill \end{array}\end{array}$

(4)

where h = Δx = Δy is grid scale. Through δ(x-X(s,t)) we have

${\mathit{f}}_{\mathit{s}}\left(\mathit{x},t\right)={{\displaystyle \int }}_{\beta }\mathit{F}\left(\mathit{s},t\right)\delta \left(\mathit{x}-\mathit{X}\left(\mathit{s},t\right)\right)d\mathit{s},$

(5)

$\mathit{u}(\mathit{X}(\mathit{s},t),t)={{\displaystyle \int }}_{\Omega }\mathit{u}\left(\mathit{x},t\right)\delta \left(\mathit{x}-\mathit{X}\left(\mathit{s},t\right)\right)d\mathit{x}.$

(6)

In (5), F(s,t) is the Lagrangian elastic force density, and the method to solve it will be discussed in the next section.

2.1. Solving F $(\text{s},t\mathbf{)}$

To solve the solid force density, it is needed to establish solid constitutive laws. According to different materials of the structure, two kinds of mechanical models (hyperelastic constitutive model and rigid model) are used in this paper.

2.1.1. Elastic Constitutive Model

The elastic structure is solved by finite element method in Lagrangian coordinate system, and it is convenient to express the stress in the reference (Lagrange) variables. Here, we use the first Piola-Kirchhoff elastic stress tensor field P(s,t). We consider an arbitrary smooth portion V in a reference configuration, and the portion is V _t = X(V,t) at time t. We have

$\begin{array}{c}{{\displaystyle \int }}_{\partial V_t}\mathit{\sigma }\mathit{n}da={{\displaystyle \int }}_{\partial V}\mathit{P}\mathit{N}dA,\\ \mathit{P}(\mathit{s},t)=\mathit{J}\mathit{\sigma }(\mathit{X}(\mathit{s},t),t){\mathit{D}}^{-\mathit{T}}(\mathit{s},t),\end{array}$

(7)

where N is the outer unit normal to the region V, and n is the outer unit normal to the region V _t . σ is the Cauchy stress tensor, and P is the first Piola-Kirchhoff elastic stress tensor. D(s,t) = ∂X(s,t)/∂s = ∇ _s X(s,t) is the deformation gradient associated with the mapping X(s,t), and D- ^T (s,t) is the inverse transpose of D(s,t). J = det⁡D is the Jacobian determinant and J≥0 at any time, so we have da = J F -^TdA.

For any virtual displacement or Lagrangian test function V(s,t), we have that

${{\displaystyle \int }}_{\beta }\mathit{F}(\mathit{s},t)\mathit{V}(\mathit{s},t)d\mathit{s}=-{{\displaystyle \int }}_{\beta }\mathit{P}(\mathit{s},t):{\nabla }_{\mathit{S}}\mathit{V}(\mathit{s},t)d\mathit{s}={{\displaystyle \int }}_{\beta }(\nabla ·\mathit{P}(\mathit{s},t))·\mathit{V}(\mathit{s},t)d\mathit{s}-{{\displaystyle \int }}_{\partial \beta }\left(\mathit{P}\left(\mathit{s},t\right)\mathit{N}\left(s\right)\right)·\mathit{V}\left(\mathit{s},t\right)dA,$

(8)

where F(s,t) is divided into two parts including the effects of both the internal force density ∇ _s ·P and the transmission force density PN(s).

The constitutive law of the elastic structure is defined by hyperelastic models, and the typical model is the St. Venant-Kirchhoff elasticity model (STVK) and neo-Hookean model (INH).

The strain-energy functional W _s ( s ,t) of STVK model is

$W_s(\mathit{s},t)=\frac{{\lambda }_s}{2}tr{\left(\mathit{E}\right)}^2+{\mu }_{s}tr\left({\mathit{E}}^2\right),$

(9)

where E = (C-I)/2 is strain tensor, and C = D ^T D is right Cauchy-Green deformation tensor. The first P-K stress tensor can be written as

$\mathit{P}\left(\mathit{s},t\right)=\mathit{D}·\frac{\partial W_s\left(\mathit{s},t\right)}{\partial \mathit{E}}=\mathit{D}\left({\lambda }_{s}tr\left(\mathit{E}\right)\mathit{I}+2{\mu }_s\mathit{E}\right).$

(10)

The strain-energy functional W_s(s,t) of INH model can be written as

$W_s\left(\mathit{s},t\right)=\frac{{\mu }_s}{2}\left({\mathbf{I}}_{\mathbf{1}}^{\mathbf{C}}-\mathrm{3}\right)-{\mu }_{\text{s}}\left(\ln J\right).$

(11)

The first P-K stress tensor in this case becomes

$\mathbf{P}\left(\mathbf{s},t\right)={\mu }_s\left(\mathbf{D}-{\mathbf{D}}^{-T}\right).$

(12)

On the other hand, the entire physical system including fluid and elastic structure was enforced incompressibility because of the principle of IBM. Elastic material parameters ${\lambda }_s={\nu }_s{E}_s/(1+{\nu }_s)(1-2{\nu }_s)$ and ${\mu }_s=E_s/2(1+{\nu }_s)$ , where ν _s is Poisson ratio and E _s is elastic modulus.

2.1.2. Rigid Model

For the rigid structure, this paper uses the feedback force immersed boundary method, and the equation can be written as

$\mathbf{F}(\mathbf{s},t)={\kappa }_s\left(\mathit{s}\right)\left(\mathit{s}-\mathit{X}\left(\mathit{s},t\right)\right),$

(13)

where κ _s (s) is a penalty coefficient, when κ _s (s)→∞, X(s,t)→s, and the structure can be approximately regarded as a rigid body.

3. Discretization of the Immersed Boundary Formulation

3.1. Finite Difference Discretization of N-S Equations

Using a Cartesian staggered grid, the N-S equation is discretized by finite difference method in space. We assume that Ω is a 1 × 1 square region and is discretized on a regular N × N grid with grid spacing Δx = Δy = h = 1/N, as shown in Figure 2. The pressure p is approximated at the black points of the grid cells, such as (i – 1,j), (i,j), and (i + 1,j), and the components of Eulerian velocity field u = (u₁,u₂) and body force f = (f₁,f₂) are approximated at the white points of the grid cells, such as (i – (1/2),j), (i + (1/2),j), and (i,j + (1/2)). The divergence operator ∇·u, gradient operator ∇p, and Laplace operator ∇²u of N-S equation use second-order accurate finite difference approximations. ∇·u is approximated at the black points of the grid cells. ∇p and ∇²u are approximated at the white points of the grid cells. Nonlinear advection terms for interior nodes are discretized as (assume ϕ is an arbitrary fluid variable)

$\varphi \frac{\partial \varphi }{\partial x}=\frac{{\varphi }_{i,j}}{2\Delta x}({\varphi }_{i+1,j}-{\varphi }_{i-1,j})-\frac{{\varphi }^{+}CF}{\Delta x}({\varphi }_{i+1,j}-3{\varphi }_{i,j}+3{\varphi }_{i-1,j}-{\varphi }_{i-2,j})-\frac{{\varphi }^{-}CF}{\Delta x}\left({\varphi }_{i+2,j}-3{\varphi }_{i+1,j}+3{\varphi }_{i,j}-{\varphi }_{i-1,j}\right),$

(14)

where ${\varphi }^{+}=(\mathrm{1}/\mathrm{2})({\varphi }_{i,j}+\left|{\varphi }_{i,j}\right|)$ , ${\varphi }^{-}=(\mathrm{1}/\mathrm{2})({\varphi }_{i,j}-\left|{\varphi }_{i,j}\right|)$ , and CF = 0.125, ϕ is an arbitrary fluid variable. For boundary nodes, a first-order upwind scheme is used.

Figure 2:

Computational region and grid for FSI-C numerical example.

3.2. Finite Element Discretization of Structure

The finite elements are generated for the structure β. β is composed of elements β _h ^e , and ${\beta }_h={\sum }_{k=\mathrm{1}}^{M_e}{\beta }_h^e$ is the total elements of the mesh. The total number of the elements is M _e . We denote by {s _l }_{l = 1} ^M the nodes of the mesh, and the total number of the nodes is M. The physical positions of the nodes of the mesh are {X _l (t)}_{l = 1} ^M , and the Lagrangian basis functions of the elements of the mesh are {ϕ _l (s)}_{l = 1} ^M . We define an approximation to X(s,t) by

${\mathbf{X}}_h\left(\mathbf{s},t\right)={\displaystyle \sum }_{l=\mathrm{1}}^M{\mathbf{X}}_l\left(t\right){\varphi }_l\left(\mathbf{s}\right),$

(15)

while, in the nodes of the mesh, X _h (s _l ,t) = X _l (t). Thus, using the Lagrangian basis functions {ϕ _l (s)}_{l = 1} ^M , the deformation gradient can be approximated as

${\mathbf{D}}_h\left(\mathbf{s},t\right)=\frac{\partial }{\partial \mathbf{s}}{\mathbf{X}}_h\left(\mathbf{s},t\right)={\displaystyle \sum }_{l=\mathrm{1}}^M{\mathbf{X}}_l\left(t\right)\frac{\partial }{\partial \mathbf{s}}{\varphi }_l\left(\mathbf{s}\right).$

(16)

Using the constitutive laws of the structure, the P _h (s,t) can be computed.

The Lagrangian force density F(s,t) of the weak formulation (8) can be computed, where V(s) is approximated by ${\varphi }_m(\mathbf{s})$ , m = 1,2,…, M, and the Lagrangian discretization of the formulation (8) can be written as

${\displaystyle \sum }_{l=1}^M\left({{\displaystyle \int }}_{\beta }^{}{\varphi }_l\left(\mathit{s}\right){\varphi }_m\left(\mathit{s}\right)d\mathit{s}\right){\mathit{F}}_l\left(t\right)=-{{\displaystyle \int }}_{\beta }^{}{\mathit{P}}_h(\mathit{s},t){\nabla }_{\mathit{s}}{\varphi }_m\left(\mathit{s}\right)d\mathit{s}+{{\displaystyle \int }}_{\partial \beta }^{}{\mathit{P}}_h\left(\mathit{s},t\right)N\left(\mathit{s}\right){\varphi }_m\left(\mathit{s}\right)dA\left(\mathit{s}\right).$

(17)

3.3. Discretion of Fluid-Structure Interaction Equations

The equations of fluid and structure interaction are (5) and (6). The fluid-structure interactions (5) and (6) are approximated by a smoothed function δ _h (x-X(s,t)) in the Cartesian staggered grid, and the Eulerian body force f = (f₁,f₂) is approximated via

$\begin{array}{c}{\mathit{f}}_h\left(\mathit{x},t\right)=\left(f_{h1}\left(\mathit{x},t\right),f_{h2}\left(\mathit{x},t\right)\right),\\ {\left(f_{h1}\right)}_{i-\left(1/2\right),j}={\displaystyle \sum }_{l=1}^M{F}_{l1}\left(t\right){\delta }_h\left({\mathit{x}}_{i-\left(1/2\right),j}-{\mathit{X}}_l\left(t\right)\right)\Delta s_1\Delta s_2,\\ {\left(f_{h1}\right)}_{i-\left(1/2\right),j}={\displaystyle \sum }_{l=1}^M{F}_{l1}\left(t\right){\delta }_h\left({\mathit{x}}_{i-\left(1/2\right),j}-{\mathit{X}}_l\left(t\right)\right)\Delta s_1\Delta s_2.\end{array}$

(18)

The shorthand form can be written as $\mathbf{f}=\hspace{0.17em}\stackrel{⃡}{S}(\mathbf{X})\mathbf{F}$ .

The Eulerian velocity field u = (u₁,u₂) is approximated via

$\begin{array}{c}{\mathit{U}}_h\left(\mathit{s},t\right)=\left(U_{h1}\left(\mathit{s},t\right),U_{h2}\left(\mathit{s},t\right)\right),\\ U_{h1}(\mathit{s},t)={\displaystyle \sum }_{0\le i,j<N}{\left(u_1\right)}_{i-\left(1/2\right),j}{\delta }_h\left({\mathit{x}}_{i-\left(1/2\right),j}-\mathit{X}(\mathit{s},t)\right)h^2,\\ U_{h2}(\mathit{s},t)={\displaystyle \sum }_{0\le i,j<N}{\left(u_2\right)}_{i,j-\left(1/2\right)}{\delta }_h\left({\mathit{x}}_{i,j-\left(1/2\right)}-\mathit{X}(\mathit{s},t)\right)h^2.\end{array}$

(19)

The shorthand form can be written as $\partial \mathbf{X}/\partial t=\hspace{0.17em}\stackrel{⃡}{R}(\mathbf{X})\mathbf{u}$ . Then the equation $\partial \mathbf{X}(\mathbf{s},t)/\partial t=\mathbf{u}(\mathbf{X}(\mathbf{s},t),t)$ can be discretized by finite difference method in time and X _h (s,t) can be computed directly.

3.4. Time Discretization

Let Δt denote the time step, and, at time n·Δt = t ⁿ , we approximate the values of X, u, and p by X ⁿ , u ⁿ , and p ⁿ , respectively. In general the full space-time discretization implies four steps.

1. Given the approximation X ⁿ of X at time t ⁿ , we move the immersed structure, getting ${\mathbf{X}}^{n+(\mathrm{1}/\mathrm{2})}$ .

   First, we preliminary compute ${\tilde{\mathbf{X}}}^{n+\mathrm{1}}$ at time t^{n + 1} via

   $\frac{{\tilde{\mathbf{X}}}^{n+\mathrm{1}}-{\mathbf{X}}^n}{\Delta t}=\hspace{0.17em}\stackrel{⃡}{R}\left({\mathbf{X}}^n\right){\mathbf{u}}^n.$

   (20)

   Then we define an approximation to X at time ${\mathbf{X}}^{n+(\mathrm{1}/\mathrm{2})}$ by

   ${\mathbf{X}}^{n+\left(\mathrm{1}/\mathrm{2}\right)}=\frac{{\tilde{\mathbf{X}}}^{n+\mathrm{1}}+{\mathbf{X}}^n}{\mathrm{2}}.$

   (21)

2. Using ${\mathbf{X}}^{n+(\mathrm{1}/\mathrm{2})}$ at time ${\mathbf{X}}^{n+(\mathrm{1}/\mathrm{2})}$ , we calculate the source term ${\mathbf{f}}^{n+(\mathrm{1}/\mathrm{2})}$ :

   ${\mathbf{f}}^{n+\left(\mathrm{1}/\mathrm{2}\right)}=\hspace{0.17em}\stackrel{⃡}{S}\left({\mathbf{X}}^{n+\left(\mathrm{1}/\mathrm{2}\right)}\right)\mathbf{F}\left({\mathbf{X}}^{n+\left(\mathrm{1}/\mathrm{2}\right)}\right).$

   (22)

3. We find the solutions u^{n + 1} and $p^{n+(\mathrm{1}/\mathrm{2})}$ to the N-S equations with the given load:

   $\rho \left(\frac{{\mathbf{u}}^{n+\mathrm{1}}-{\mathbf{u}}^n}{\Delta t}+\frac{\mathrm{3}}{\mathrm{2}}{\mathbf{u}}^n·{\nabla }_h{\mathbf{u}}^n-\frac{\mathrm{1}}{\mathrm{2}}{\mathbf{u}}^{n-\mathrm{1}}·{\nabla }_h{\mathbf{u}}^{n-\mathrm{1}}\right)=-{\nabla }_h{p}^{n+\left(\mathrm{1}/\mathrm{2}\right)}+{\mu }_{\mathbf{f}}{\nabla }_h^{\mathrm{2}}\left(\frac{{\mathbf{u}}^{n+\mathrm{1}}+{\mathbf{u}}^n}{\mathrm{2}}\right)+{\mathbf{f}}^{n+\left(\mathrm{1}/\mathrm{2}\right)},\begin{array}{c}{\nabla }_h{\mathbf{u}}^{n+\mathrm{1}}=\mathrm{0}.\end{array}$

   (23)

4. We correct ${\tilde{\mathbf{X}}}^{n+\mathrm{1}}$ that is calculated in step (1) at time t^{n + 1} via

   $\frac{{\mathbf{X}}^{n+\mathrm{1}}-{\mathbf{X}}^n}{\Delta t}=\hspace{0.17em}\stackrel{⃡}{R}\left({\mathbf{X}}^{n+\left(\mathrm{1}/\mathrm{2}\right)}\right)\frac{{\mathbf{u}}^{n+\mathrm{1}}+{\mathbf{u}}^n}{\mathrm{2}}.$

   (24)

Because time step-lagged values of u and p are used by the time stepping scheme above, we cannot use that scheme for the initial time step. When n = 0, we deal with it as follows.

1. Given the values of X⁰, u⁰, and f⁰ at time t⁰, we preliminarily compute ${\tilde{\mathbf{X}}}^{\mathrm{1}}$ and ${\tilde{\mathbf{u}}}^{\mathrm{1}}$ via

   $\begin{array}{c}\rho \left(\frac{{\tilde{\mathbf{u}}}^{\mathrm{1}}-{\mathbf{u}}^{\mathrm{0}}}{\Delta t}+{\mathbf{u}}^{\mathrm{0}}·{\nabla }_h{\mathbf{u}}^{\mathrm{0}}\right)=-{\nabla }_h{\tilde{p}}^{\mathrm{1}/\mathrm{2}}+{\mu }_{\mathbf{f}}{\nabla }_h^{\mathrm{2}}\left(\frac{{\tilde{\mathbf{u}}}^{\mathrm{1}}+{\mathbf{u}}^{\mathrm{0}}}{\mathrm{2}}\right)+{\mathbf{f}}^{\mathrm{0}},\\ {\nabla }_h{\tilde{\mathbf{u}}}^{\mathrm{1}}=\mathrm{0},\\ {\mathbf{f}}^{\mathrm{0}}=\stackrel{⃡}{S}\left({\mathbf{X}}^{\mathrm{0}}\right)\mathbf{F}\left({\mathbf{X}}^{\mathrm{0}}\right),\\ \frac{{\tilde{\mathbf{X}}}^{\mathrm{1}}-{\mathbf{X}}^{\mathrm{0}}}{\Delta t}=\stackrel{⃡}{R}\left({\mathbf{X}}^{\mathrm{0}}\right){\mathbf{u}}^{\mathrm{0}},\end{array}$

   (25)

   where we let ${\tilde{p}}^{\mathrm{1}/\mathrm{2}}=\mathrm{0}$ at initial time.

2. We compute ${\mathbf{X}}^{\mathrm{1}/\mathrm{2}}$ and ${\mathbf{u}}^{\mathrm{1}/\mathrm{2}}$ at time $t^{\mathrm{1}/\mathrm{2}}$ via

   ${\mathbf{X}}^{\mathrm{1}/\mathrm{2}}=\frac{{\tilde{\mathbf{X}}}^{\mathrm{1}}+{\mathbf{X}}^{\mathrm{0}}}{\mathrm{2}},\hspace{1em}\hspace{1em}{\mathbf{u}}^{\mathrm{1}/\mathrm{2}}=\frac{{\tilde{\mathbf{u}}}^{\mathrm{1}}+{\mathbf{u}}^{\mathrm{0}}}{\mathrm{2}}.$

   (26)

3. We correct ${\tilde{\mathbf{X}}}^{\mathrm{1}}$ , ${\tilde{\mathbf{u}}}^{\mathrm{1}}$ , and ${\tilde{p}}^{\mathrm{1}/\mathrm{2}}$ that are calculated in step (1) via

   $\rho \left(\frac{{\mathbf{u}}^{\mathrm{1}}-{\mathbf{u}}^{\mathrm{0}}}{\Delta t}+{\mathbf{u}}^{\mathrm{1}/\mathrm{2}}·{\nabla }_h{\mathbf{u}}^{\mathrm{1}/\mathrm{2}}\right)=-{\nabla }_h{p}^{\mathrm{1}/\mathrm{2}}+{\mu }_{\mathbf{f}}{\nabla }_h^{\mathrm{2}}\left(\frac{{\mathbf{u}}^{\mathrm{1}}+{\mathbf{u}}^{\mathrm{0}}}{\mathrm{2}}\right)+{\mathbf{f}}^{\mathrm{1}/\mathrm{2}},$

   (27)
   ${\nabla }_h{\mathbf{u}}^{\mathrm{1}}=\mathrm{0},$

   (28)
   ${\mathbf{f}}^{\mathrm{1}/\mathrm{2}}=\stackrel{⃡}{S}\left({\mathbf{X}}^{\mathrm{1}/\mathrm{2}}\right)\mathbf{F}\left({\mathbf{X}}^{\mathrm{1}/\mathrm{2}}\right),$

   (29)
   $\frac{{\mathbf{X}}^{\mathrm{1}}-{\mathbf{X}}^{\mathrm{0}}}{\Delta t}=\stackrel{⃡}{R}\left({\mathbf{X}}^{\mathrm{1}/\mathrm{2}}\right)\frac{{\mathbf{u}}^{\mathrm{1}}+{\mathbf{u}}^{\mathrm{0}}}{\mathrm{2}}.$

   (30)

When the values of X¹, u¹, and $p^{\mathrm{1}/\mathrm{2}}$ are solved, the full space-time discretization can be used.

4. The Verification Example of IBM

4.1. The Establishment of Computational Model

Figure 2(a) is the computational domain of channel flow. A static rigid cylinder is immersed in fluid, a cantilever beam is attached to the cylinder. The self-induced oscillation of a cantilever beam immersed in a fluid (FSI-C) is one of the benchmarks of fluid-structure interaction. The elastic beam has width h = 0.02 m and length l = 17.5h. The computational domain has length L = 125h and width H = 20h. The center of cylinder is C = (10h, 10h), as the coordinate origin is set in left bottom corner, and the radius of cylinder is r = 2.5h. A control point A fixed at the trailing edge of the structure is with initial position (30h, 10h). The quadrilateral finite element mesh of the beam and uniform Cartesian grid of the fluid are used, as shown in Figure 2(b).

The upper and lower boundary and solid surface is no-slip wall boundary condition.

A parabolic velocity profile is prescribed at the left channel inflow:

$\begin{array}{c}{u}_x\left(\mathrm{0},y\right)=\mathrm{1.5}\stackrel{-}{U}\frac{y\left(H-y\right)}{{\left(H/\mathrm{2}\right)}^{\mathrm{2}}},\\ u_y(\mathrm{0},y,t)=\mathrm{0},\end{array}$

(31)

where $\stackrel{-}{U}$ is the mean inflow velocity.

The starting procedure for the nonsteady tests was to use a smooth increase of the velocity profile against time, namely,

$u_x(\mathrm{0},y,t)=\{\begin{array}{cc}\frac{u_x\left(\mathrm{0},y\right)\left(\mathrm{1}-\cos \left(\pi t/\mathrm{2}\right)\right)}{\mathrm{2}},\hfill & t\le \mathrm{2.0}\hspace{0.17em}\text{s}\hfill \\ u_x(\mathrm{0},y).\hfill & t>\mathrm{2.0}\hspace{0.17em}\text{s}.\hfill \end{array}$

(32)

Zero normal traction and zero tangential slip boundary conditions are used along the outflow. The no-slip condition is prescribed for the fluid on the top and bottom wall as well as cylinder surface.

For this benchmark problem, the fluid is assumed as incompressible and Newtonian. The fluid properties are ρ _f = 1.0 × 10³ kg/m³ and ν _f = 1.0 × 10⁻³ m²/s. The mass density of the elastic beam is the same as the fluid, and a STVK model is used to describe the elasticity of the beam with μ_s = 2.0 × 10⁶ kg/(ms²) and Poisson ratio ν_s = 0.4. The mean inflow velocity is $\stackrel{-}{U}=\mathrm{2.0}\hspace{0.17em}\text{m}/\text{s}$ , and the Reynolds number is $\mathrm{Re}=\mathrm{2}r\stackrel{-}{U}/{\nu }_f=\mathrm{200}$ .

4.2. Analysis of Computational Results

We set κ_s = 1.0 × 10¹¹, and the total time of calculation is 20's. This benchmark is usually used to verify the reliability of numerical algorithm of fluid-structure interaction. The results of the benchmark are summarized in Table 1, and other reference values of Table 1 are derived from paper [31, 32]. Comparisons will be done particularly for one full period of the oscillation of the beam with respect to the position of the point A at the end of the beam structure (see Figure 2). In Table 1, v _x (A) and v _y (A) are the displacement in x- and y-direction of the point A in a period of vibration of beam, and the frequencies f_x and f_y are obtained for the displacements v _x (A) and v _y (A), respectively. The drag and lift forces F _D and F _L can be written as

${\left(F_D,F_L\right)}^T={{\displaystyle \int }}_s{\mathbf{\sigma }}_{\mathbf{f}}\mathbf{n}dS,$

(33)

where S denotes the structure being in contact with the fluid and n is the outer unit normal vector to the integration path with respect to the fluid domain. The column “Unknown” refers to the total degree of freedom number (in space), namely, the sum of unknown for all velocity components, pressure, and displacement components.

Table 1:

Comparison of the computational results between other papers and this paper.

		Unknown	Δt/s	v x (A)/(× 10−3 m)	v y (A)/(× 10−3 m)	F D	F L	fx	fy
Reference	1	942258	1.0e–3	−2.91 ± 2.77	1.47 ± 35.26	459.9 ± 27.92	1.84 ± 157.7	11.63	4.98
	2	72696	2.5e–3	−2.64 ± 2.48	1.38 ± 33.25	451.1 ± 24.57	2.04 ± 150.6	10.73	5.38
	3	72696	5.0e–4	−2.84 ± 2.67	1.28 ± 34.61	452.4 ± 26.19	2.36 ± 152.7	10.84	5.42
	4	304128	2.5e–4	−2.88 ± 2.72	1.47 ± 34.99	460.5 ± 27.74	2.50 ± 153.9	10.93	5.46
	5	324480	2.0e–5	−4.54 ± 4.34	1.50 ± 42.50	467.5 ± 39.50	16.2 ± 188.7	10.12	5.05
Present	1	263410	5.0e–5	−2.22 ± 2.29	2.62 ± 32.50	488.4 ± 54.10	−5.8 ± 176.2	11.22	5.61
	2	1053660	5.0e–5	−2.70 ± 2.27	1.37 ± 32.38	457.3 ± 42.65	−4.1 ± 172.5	10.80	5.40
	3	2167500	5.0e–5	−2.76 ± 2.25	1.32 ± 32.20	448.3 ± 37.20	−1.5 ± 169.5	10.66	5.33

Time dependent values are represented by mean value (me), amplitude (am), and frequency. The mean value and the amplitude are computed from the last period of the oscillations by taking the maximum(max) and minimum(min) values; namely, me = (max⁡ + min⁡)/2 and am = (max⁡-min⁡)/2. In Table 1, the computational results of this paper are all within the range of results obtained by other numerical methods, which verify the reliability of this method in this paper.

5. The Establishment of Self-Propelled Swimming Fish Model

In this paper, in order to simulate the interaction between internal forces, fish body movement, and fluid, a volume force is added and acted on the tail of the fish and the force is simplified from internal forces produced by muscle, which is reflected truly by the real physiological condition of fish swimming. The coordinate of point A is (0.2, 0), which is used to output the fish displacement and velocity results of fish-like movement. The fish is considered as isotropic INH hyperelastic material and using the NACA0012 airfoil represents the body of fish, as shown in Figure 3. The total length is 1 m, and the fish body is divided into three sections, namely, head (0–0.5 m), tail (0.5–0.96 m), and tail fin (0.96–1 m). The coordinate of point A is (0.2, 0), which is used to output the fish displacement and velocity results.

Figure 3:

Model of fish.

Fish-like swim only along the -X direction is considered in this paper; then the control of the self-propelled swimming fish model is as follows.

(1) Head Section (L₁). An elastic constraint is added in Y direction using feedback force method and this control formulation is

$\mathit{F}\left({\mathit{s}}_y,t\right)={\kappa }_s\left({\mathit{s}}_y\right)\left({\mathit{s}}_y-\mathit{X}\left({\mathit{s}}_y,t\right)\right),$

(34)

where κ _s ( s _y ) is a coefficient of elasticity. sy is the Y direction coordinate in reference configuration. The elastic modulus of the head section is set to E₁.

(2) Tail Section (L₂). The internal force produced by muscle is simplified by adding a volumetric force on the tail section of fish body. At first, this volumetric force is loaded in + Y direction, and the tail begins to move. When an arbitrary point of the tail swings to a specified location, the direction of the volumetric force reverses, and the tail begins to swing to -Y direction. And then the direction of the volume force changes again when an arbitrary point of the tail swings to a specified location on the other side. Direction of the volumetric force is repeated to reverse, which simulated the periodic swing of the tail. The distribution of the volumetric force density on the tail section is

$F_y=F_{y\_\max }\times \sin \left(\mathrm{2}\pi \times \left({\mathbf{s}}_x-\mathrm{0.5}\right)\right),$

(35)

where sx is the X direction coordinate in reference configuration. The elastic modulus of the tail section is set to E₂.

(3) Caudal Fin Section (L₃). This section of the body is not controlled, and the elastic modulus is E₃.

6. Numerical Simulation of Single Fish Swimming

6.1. The Establishment of Computational Model

The computational domain has length L = 20 m and width H = 3 m, as shown in Figure 4. A control point P fixed at the head of the fish with initial position (1.6 m, 1.5 m) as the coordinate of left bottom corner is (0, 0). The elastic coefficient in Y-direction on head section is κ _s (s _y ) = 10 × 10⁹.

Figure 4:

Computational region for the single swimming fish.

The fluid is static at the initial time, and the fluid density and kinematic viscosity are ρ _f = 1.0 × 10³ kg/m³ and ν _f = 1.0 × 10⁻⁶ m²/s, respectively. The mass density of the fish is the same as the fluid, and the elastic modulus of head and caudal fin section is E₁ = E₃ = 2.4 × 10⁶ Kg/(ms²) and the tail section is E₂ = 2.4 × 10⁶ Kg/(ms²). The maximum volumetric force density is F_{y_max} = 1.5 × 10⁶ N/(m³).

In the whole computational domain of fluid using four layers adaptive Cartesian grid, the minimum grid spacing is 0.0045 m. The triangular finite element mesh for the fish is used, and the mesh number is 2090, as shown in Figure 5.

Figure 5:

The grid for the swimming fish.

6.2. Analysis of Computational Results

Figure 6 shows fluid motion around the fish during steady swimming. A strong vortex sheds along the tail of caudal fin each time and forms the anti-Kármán vortex street when the tail swings back to the symmetry axis of fish body. Multiple weak secondary vortices are also produced from the head of the fish and move along body and shed finally at the tail of caudal fin.

Figure 6:

The anti-Kármán vortex street to the swimming fish.

Figure 7 shows the swimming speed from t = 0 s to a steady-state speed; in a whole, the form of swimming speed is similar to the results of paper [33]. It is seen that the velocity is gradually increased from t = 0 s to 3.0's, and then the swimming speed is up to a stable status. Also, we can see that the speed is fluctuant and increases, which is mainly affected by the periodical swing of the tail; simultaneously, there is small fluctuation in a period, which is attributed to the shedding of the weak secondary vortices along the fish body.

Figure 7:

The swimming speed against time.

Figure 8 shows the flow distribution and fish shape during start-up process of swimming. The fish body has become “S” shape affected by the fluid dynamic force when the tail is moving opposed to the axis of body, and the body has become “C” shape when tail is moving back to the axis of body. The numerical simulation result of fish-like swimming is similar to the result of experiment [34]. There is a negative acceleration in 0–0.07's, as shown in Figure 7(b). After 0.07s, the fish began to accelerate motion until 0.22's. It can be seen that the tail moves up away to the axis of body from 0.07's, and the tail fin is passively sloping down and the fish body gradually becomes “S” shape until 0.14's. To “S” shape of upmovement, the pressure acted on the two sides of tail fin will contribute to accelerate the motion of fish, which is shown in Figure 9. Subsequently, the fish tail moves down back to the axis of body, and the tail fin is passively sloping up and the fish body gradually becomes “C” shape. From 0.22's to 0.32's, the velocity is fluctuating, the tail is moving down, and fish body is from “C” shape to “S” shape. When fish body of “S” shape continues to move down away to the axis of body, the accelerated motion of fish is also achieved from 0.32's to 0.4's. From 0.4's to 0.5's, the velocity is also fluctuated, and the fish tail is moving up, and fish body is from “C” shape to “S” shape.

Figure 8:

Velocity vector and the pressure contour in start-up process.

Figure 9:

Velocity vector and pressure contour in cruise process.

Figure 9(a) shows velocity vector distribution and fish shape during a steady state of swimming. It is seen that the fish body continues its conversion between the “C” shape and “S” shape when swimming, which is similar to the body move during the start-up process, but the swimming velocity of fish is periodic and consists of two deceleration and two accelerate stages, which is shown in Figure 7(c). From Figure 7(c) and corresponding to Figure 9, it is concluded that the swimming is deceleration from 4.45's to 4.575's in whole when the tail moves up away to the axis of body and is accelerated motion from 4.575's to 4.645's when the tail moves down back to the axis of body. Subsequently, the swimming is deceleration from 4.645's to 4.77's when the tail moves down away to the axis of body and is accelerated motion from 4.77's to 4.84's when the tail moves up back to the axis of body, although there is some small velocity fluctuation in each stage.

Figure 9(b) shows pressure distribution and fish shape during a steady state of swimming. It is shown that the fish body is continuous conversion between “C” shape and “S” shape when the fish is self-propelled swimming. A clockwise swirl vortex is shed from the tail fin when the tail moves up and an anticlockwise swirl vortex comes into being after the tail fin when the tail moves down. The pressure difference on fish tail and tail fin due to conversion of “C” and “S” shape of fish body causes the fish body to continue to gain the energy from the fluid, which achieved the self-propelled swimming forward.

7. Conclusion

In this paper, the finite element immersed boundary method combined with traditional feedback force method is solved numerically with time and space discretization of computational fluid dynamics, and a self-induced vibration of a cantilever beam immersed in a fluid is one of the benchmarks of fluid-structure interaction, which is used to verify the validity of the numerical method.

This paper has established a flexible model of freely self-propelled swimming fish. The fish is considered isotropic neo-Hookean hyperelastic material and using the NACA0012 airfoil represents the body of fish. The advantage of this model is that it could reflect the interaction between fish movement, internal force and the external fluid. Using IBM, the self-propelled swimming of “C” type fish is numerically simulated. The flow field characteristics and fish body movement characteristics in start and cruise process of a single swimming fish is analyzed. The results reveal the fluid dynamics mechanism of the conversion of “C” type and “S” type of fish.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.