The Interaction Mechanisms of Swimming Biomimetic Fish Aligned in Parallel Using the Immersed Boundary Method

Abstract

In natural environments, fish almost always swim in groups. Investigating the coupled mechanism of biomimetic fish exhibiting autonomous swimming capabilities advances our understanding of fish schooling phenomena and simultaneously aids in refining the structural and formation configurations of underwater robotic vehicles. This work innovatively develops an algorithm based on the Direct-Forcing Immersed Boundary Method (DF-IBM) and implements it in an efficient, modular software program written in C++. The program accelerates the calculation process by using a multigrid method. Validation against a benchmark case of flow around a cylinder, with comparison to data from the existing literature, verifies the program’s precision with discrepancies of less than 3.6%. Based on this algorithm, the paper analyzes the incompressible viscous flow during the movement of parallel-aligned biomimetic fish. It uncovers the interaction between the fish’s motion and the surrounding flow field and also reveals the hydrodynamic mechanisms of the group motion of the parallel-aligned biomimetic fish. The flow field under varying spacing and phases between the parallel-aligned biomimetic fish proves that the interaction between the flow fields induced by the two fish bodies becomes increasingly significant when decreasing the lateral spacing from $1.4L$ to $0.6L$. Notably, an initial lateral convergence of the fish bodies is observed, followed by a sideways swimming pattern at a particular pitch angle, accompanied by a decrement in their forward swimming velocity as they approach each other. Additionally, this study compares flow field alterations in parallel-aligned biomimetic fish with identical lateral spacing but opposing flapping phases. The findings indicate that, irrespective of the phase, the fish exhibit an initial convergence followed by a sideways motion at a specific pitch angle. However, due to disparities in the tail’s flow field, a larger pitch angle is generated when the fish swim in unison. All the findings above will provide a solid theoretical foundation for the design and optimization of underwater robotic vehicles.

1. Introduction

Applying the swimming motion of fish to biomimetic research not only uncovers the mechanical and morphological mechanisms of biological movement but also offers important insights for biomimetic applications. Moreover, research on the swimming mechanisms of biomimetic robot fish can enhance the maneuverability, propulsion, and adaptability of underwater robots in complex environments. They provide critical technical support for the economic, scientific, and sustainable development and utilization of marine resources. Therefore, numerous quantitative analyses of the flow characteristics around fish bodies have been conducted to figure out the fundamental mechanisms of fish swimming.

Experimental methods [1,2], especially the ones using digital particle image velocimetry (DPIV), are widely applied in measuring the flow field. DPIV has been utilized to study and analyze swimming zebrafish, visualizing different vortex structures induced by three different swimming states, including the C-start rapid escape, steady forward swimming, and a transition from forward swimming to gliding [3]. They also demonstrated the kinematic and dynamic characteristics of the fish’s midline in these three states. Tytell [4] employed DPIV and high-speed cameras to study the kinematics and dynamics under Anguilliform mode during forward swimming. It is found that in a state of stable swimming, the tail flow consisted of lateral jets with minimal net axial momentum, reflecting a balance between thrust and drag. During acceleration, the jets reoriented downstream, increasing the axial momentum of the fluid. The swimming of needlefish under the Anguilliform mode was also experimentally investigated by using high-speed cameras to analyze the kinematic characteristics of the body trunk, dorsal fin, caudal fin, pectoral fins, and anal fin at different swimming speeds [5]. However, experimental methods involve significant investments, are time-consuming, and have relatively large experimental errors.

The Immersed Boundary Method (IBM), initially proposed by Peskin [6], is a classic non-body-fitted mesh method. As illustrated in Figure 1, in the IBM mesh approach, the fluid domain and solid domain are discretized using a fixed background Eulerian mesh and a moving Lagrangian mesh, respectively. The Eulerian mesh, applied throughout the problem domain, does not require a consideration of the solid boundaries, thereby simplifying the pre-processing mesh generation process and eliminating the need for regenerating the fluid mesh. Numerical interactions between the fluid and solid are accomplished through interpolation using a Delta function. This method does not require the generation of complex body-fitted meshes, and it simplifies the treatment of interfaces while offering high computational efficiency. Therefore, it has extensive applications in biofluid dynamics.

Exploring swimming mechanisms from the perspective of the hydromechanics of external fluid has also drawn significant attention. An advanced moving IBM tailored for unstructured grids, including polyhedral and triangular, was reported. The inherent adaptability of unstructured grids facilitates the utilization of finer grid resolutions and the meshing of intricate rigid geometries [7]. Borazjani et al. [8,9,10] combined the basic principles of IBM with finite difference methods to comparatively analyze the differences in tail vortices between Anguilliform and Carangiform modes and studied the swimming performance of fish under these two modes. The IBM method was also used to investigate the hydrodynamic characteristics of a self-propelled flexible plate and a biomimetic manta ray, revealing the propulsion mechanisms and movement patterns that achieve efficient swimming speeds [11,12]. An integrated IBM-based parallel computing for hydrodynamics with structural dynamics simulations was applied to study vortex-induced oscillations of a 3D slender flexible body under an inclined flow [13]. It can be noticed that the spacing and phase differences play an important role in the energy harvesting of dual-flapping wings of different shapes. Therefore, Ma et al. [14] employed numerical simulations to examine the effects of those two parameters, revealing significant changes in energy absorption efficiency compared to a single wing, whether in tandem or staggered arrangements. With the rapid development of computational technology, a numerical framework combining deep learning with the Immersed Boundary-Lattice Boltzmann Method was developed to study the adaptive behavior of self-propelled fish in complex environments, simulating fish in incompressible viscous flows and training models to optimize task performance [15]. In addition, a sharp-interface IBM was also developed to simulate interactions between fluid flow and deformable moving bodies, with results showing that the algorithm performed well in simulating flow movement around deformable and moving objects [16]. Bontoux et al. [17] proposed an efficient algorithm for simulating the interaction between deformable bodies and 2D, incompressible flows, and based on this, they introduced a curvature calculation method for efficiently controlling eel-like swimming toward specified targets. The validation results indicate that their algorithm offers high computational efficiency and accuracy for fish-like swimming control and various fluid–structure interactions. Despite the development of numerous numerical methods aimed at studying the swimming of biomimetic fish, accurately capturing the moving boundary of the fish body remains a significant challenge.

Compared to traditional methods for computing fish swimming, the Direct-Forcing Immersed Boundary Method (DF-IBM) applied in this work offers several advantages. Specifically, DF-IBM does not require the generation of body-fitted grids, thereby avoiding the complexities associated with meshing intricate geometries and the computational challenges that come with it. This makes DF-IBM particularly well suited for handling complex shapes. Additionally, the incorporation of multigrid acceleration techniques enhances computational efficiency. These features collectively make DF-IBM a powerful tool for simulating fluid–structure interactions in biological systems. To uncover the mechanisms behind autonomous swimming in fish schools, this paper develops an algorithm based on DF-IBM to examine the hydrodynamic characteristics of parallel-aligned biomimetic fish swimming autonomously in still water. The influence of lateral spacing and phase differences on swimming performance metrics such as forward speed, drag, and the structures of tail vortices is systematically analyzed.

2. Numerical Methodology

The water surrounding swimming fish can be considered an incompressible fluid. Therefore, 2D dimensionless incompressible continuity and Navier–Stokes (N–S) equations are required to simulate the flow field around the swimming fish, which are given as follows:

$$\nabla ·\mathit{U}=0$$

(1)

$${\displaystyle \frac{\partial \mathit{U}}{\partial \mathrm{t}}}+\left(\mathit{U}\cdot \nabla \right)\mathit{U}=-\nabla p+{\displaystyle \frac{1}{\mathrm{Re}}}{\nabla }^2\mathit{U}+f$$

(2)

where $\mathit{U}\left(u,v\right)$ represents the velocity in the flow field; $p$ represents the pressure; $\mathrm{Re}$ represents the Reynolds number, which characterizes the relative magnitude of inertial forces to viscous forces; and $f$ represents the volumetric force indicating the external force applied to the fluid, specifically the interaction force between the object’s boundary and the flow field.

The moving boundary of the fish body is processed using the DF-IBM proposed by Fadlun et al. [18], which integrates the finite difference method. In this approach, the effect of the boundary on the flow field is also replaced by the volumetric force source term in the N–S equation. The discretized N-S equation is expressed as follows:

$${\displaystyle \frac{u^{n+1}-u^n}{\Delta t}}=rh{s}^{n+{\displaystyle \frac{1}{2}}}+f^{n+{\displaystyle \frac{1}{2}}}$$

(3)

where $rhs$ represents the sum of the convective term, diffusive term, and pressure gradient term, and $f$ is the volumetric force source term. During computation, since the boundary points applied in the IBM often do not perfectly align with the grid points in the computational domain, this paper uses an interpolation function as applied in previous work [18]. In order to enhance computational efficiency, the Poisson equation is solved using a multigrid method, and the convective term is discretized using a high-order weighted essentially non-oscillatory (WENO) scheme [18].

Time discretization is carried out using the step projection method [19]. At any given time step $n$, assuming the velocity field $u^n$ and pressure field $p^n$ are known, the first intermediate velocity $u^{\ast }$ is calculated by

$${\displaystyle \frac{u^{\ast }-u^n}{\Delta t}}={\displaystyle \frac{3}{2}}{A}^n-{\displaystyle \frac{1}{2}}{A}^{n-1}$$

(4)

where $A^n=-\nabla \cdot \left(uu\right)+{\displaystyle \frac{{\nabla }^{2}u}{\mathrm{Re}}}$. Considering the influence of the pressure term, we have

$${\displaystyle \frac{u^{\ast \ast }-u^{\ast }}{\Delta t}}=-\nabla p^{n+1}$$

(5)

Finally, solving the pressure Poisson equation:

$${\nabla }^2{p}^{n+1}={\displaystyle \frac{1}{\Delta t}}\nabla \cdot u^{\ast }$$

(6)

The force source term is then computed to obtain the velocity at the next time step $u^{n+1}$.

Since 94% of the total computational time is allocated to resolving the Poisson equation [20], this study employs the multigrid method to accelerate the computational process inherent to the Poisson equation. The fundamental principle behind the multigrid method is to use coarser grids to more efficiently eliminate low-frequency errors inherent in finer grids, thereby accelerating the convergence rate.

3. Validation of Numerical Method

3.1. Cylinder Wake

In this study, the accuracy of the proposed DF-IBM is validated using the 2D flow around a cylinder, which is a commonly used benchmark to verify the accuracy of numerical methods for simulating the autonomous swimming of biomimetic fish. Numerical simulations are performed at a Reynolds number of $\mathrm{Re}=100$, and the results indicate that vortices form behind the cylinder, shedding periodically to create a Kármán vortex street. During this process, the drag and lift coefficients, $C_d$ and $C_l$, of the 2D cylinder, also undergo periodic variations, which are simplified as follows:

$$C_d={\displaystyle \frac{2F_D}{\rho U_{\infty }^{2}D }}$$

(7)

$$C_l={\displaystyle \frac{2F_L}{\rho U_{\infty }^{2}D }}$$

(8)

where $\rho$ is the fluid density, set to 1; $U_{\infty }$ is the velocity of the incoming flow; $D$ is the cylinder diameter; and $F_D \mathrm{and} F_L$ are the drag and lift forces acting on the cylinder, respectively. In DF-IBM, these coefficients can be calculated using the following equations:

$$F_D=-\underset{0}{\overset{m}{{\displaystyle \int }}}{f}_{x}dx$$

(9)

$$F_L=-\underset{0}{\overset{m}{{\displaystyle \int }}}{f}_{y}dy$$

(10)

where $m$ represents the boundary length, while $f_x$ and $f_y$ represent the Lagrangian forces at boundary points in the $x$ and $y$ directions, respectively. For the explicit DF-IBM, these forces are the sums of force at all boundary points during the iteration process, and Ω represents the computational domain encompassing the entire boundary.

The specific conditions for the cylinder flow simulation are given as follows: the Neumann boundary condition is applied at the outlet, the slip boundary condition is applied at the upper and lower boundaries, and the no-slip boundary condition is applied on the cylinder surface. With the cylinder diameter being set as the characteristic length, the diameter is thus non-dimensionalized as $D$ = 1. The computational domain is set to 12$D$ × 24$D$, with a time step of dt = 0.0005, and the cylinder boundary is discretized into 320 points. Figure 2a,b represent vorticity maps of the cylinder wake achieved in this study and in reference [21] at Re = 100, and Figure 2c is the corresponding streamline obtained in this study, where consistent periodic Kármán vortex streets are observed behind the cylinder.

As shown in Figure 3, once the calculation stabilizes, the drag and lift coefficients exhibit periodic fluctuations. These oscillations are caused by the periodic shedding of vortices, leading to periodic fluctuation in the coefficients.

Table 1. A comparison of the drag and lift coefficients at a Reynolds number of 100.

	$\overline{{\mathit{C}}_{\mathit{d}}}$	$\overline{{\mathit{C}}_{\mathit{d}}}$ Error	${\mathit{C}}_{\mathit{l}\mathbf{m}\mathbf{a}\mathbf{x}}$
This study	1.39	/	0.33
Lai and Peskin [22]	1.44	3.6%	0.33
Su [23]	1.40	0.7%	/

compares the average drag coefficient $\overline{C_d}$ and the maximum lift coefficient $C_{l\max }$ obtained in this study with the corresponding results from Lai and Peskin [22] and Su [23]. The maximum relative error is only 3.6%, demonstrating that the results from this study align well with the reference data, confirming the program’s effectiveness in simulating stationary cases.

3.2. Autonomous Swimming of Biomimetic Fish

Before applying the IBM to investigate the swimming dynamics of tandem bio-inspired fish, it is essential to validate the accuracy of the proposed method in handling moving boundaries. Cai et al. [24,25] presented a novel implicit IBM for simulating incompressible viscous flow over complex stationary or moving solid boundaries. The method enforces the no-slip constraint at solid boundaries by treating the boundary force as a Lagrange multiplier, avoiding excessive constraints on the time step. It then introduces an additional moving force equation to determine the boundary force, making it highly suitable for moving boundary problems. Since this study focuses on tandem swimming of bio-inspired eels, a geometry model of a single bio-inspired eel is used for the swimming validation.

The numerical simulations of the fish body in this study are based on the profile of an eel-like fish shown in Figure 4 [26]. The fish body has a length of $L=1$ and a thickness of $D=0.04$, and boundary conditions are applied to ensure that no elastic deformation occurs during oscillation. The fish body is divided into three regions along the horizontal axis: the head, the trunk, and the tail. The shape of the fish body in a horizontal cross-section is described by the following function:

$$w\left(s\right)=\left\{\begin{array}{l}\sqrt{2w_{h}s-s^2} , 0\le s<s_b\\ w_h-\left(w_h-w_t\right){\left({\displaystyle \frac{s-s_b}{s_m-s_b}}\right)}^2 , s_b\le s\le s_m\\ w_t{\displaystyle \frac{L-s}{L-s_m}} , s_m\le s\le L\end{array}\right.$$

(11)

where $w_h=s_b=0.04L$, $s_m=0.95L$, $w_t=0.01L$, and $s$ is the arc length of the spinal centerline.

The movement mode used in this study is based on the flapping pattern [26,27,28,29,30,31]. The dimensionless equation, with the body length as the characteristic length, is expressed through the oscillation equation of the body centerline as follows:

$$y_b\left(s,t\right)=A_m{\displaystyle \frac{\left(\frac{s}{L}+b\right)}{1+b}}\sin \left(2\pi \left({\displaystyle \frac{s}{L}}-{\displaystyle \frac{t}{T}}\right)\right)$$

(12)

where $y_b\left(s,t\right)$ represents the lateral displacement of the centerline, $A_m$ is the amplitude coefficient, $b$ is the amplitude parameter, and $T$ is the oscillation period of the fish body. In this study, $A_m=0.125$ and $b=0.03125$. Figure 5a shows the profile variation in an eel fish during one oscillation cycle. Figure 5b shows the contour lines of the fish body’s oscillation over one cycle. The angle $\theta$ is defined as the angle between the axis of the fish body coordinate system and the global ground coordinate system, representing the pitch angle of the fish. As shown in Figure 5c, two coordinate systems are established in this work to display the autonomous propulsion of the biomimetic fish, which are a global ground coordinate system ($oxy$) and a body-fixed coordinate system ($o_b{x}_b{y}_b$) on the fish. The origin of the body-fixed coordinate system is located at the center of mass, with $x_b$ representing the line from the fish’s head to its center of mass, and $O_b{y}_b$ situated perpendicular to $O_b{x}_b$.

Prior to investigating the swimming patterns of parallel biomimetic fish, it is imperative to validate the simulation of a single biomimetic eel’s locomotion. For this purpose, the kinematic parameters are set as follows: the viscosity coefficient $\mu =1.4\times {10}^{-4}$, the body oscillation frequency $f=1$, and the tail amplitude $A_m=0.125$. The computational domain, as shown in Figure 6, measures $4L\times 16L$, with a grid resolution of 512 × 2048 cells. Symmetry boundary conditions are applied around the entire flow field, allowing fluid to pass freely through these boundaries. The body boundary is discretized into 402 points, and the fish starts its motion from the right side of the flow field, moving toward the left. The time step for flow field calculations is set to ${\displaystyle \frac{T}{1000}}$, where $T$ represents the oscillation period of the fish body. This numerical case aims to explore the force states experienced by the single fish during autonomous forward swimming and the instantaneous changes in the surrounding fluid dynamics.

In the process of autonomous swimming, the undulatory motion of the fish body serves as an intrinsic energy source that drives interactions with the surrounding flow field. The deformations of the fish body generate forces on the fluid, to which the flow field responds with opposing forces, facilitating both the translational and rotational movements of the fish within the plane. Beginning from its initial configuration, the fish initiates undulatory motion, and upon reaching a stable cruising state, the evolution of the flow field is systematically analyzed.

Figure 7 and

Table 2. Comparison of steady cruising speeds.

	$\mathbf{Forward} \mathbf{Speed} {\mathit{U}}_{\Vert }$	$\mathbf{Lateral} \mathbf{Speed} {\mathit{U}}_{\perp }$	$\mathbf{Error} \mathbf{of} {\mathit{U}}_{\Vert }$	$\mathbf{Error} \mathbf{of} {\mathit{U}}_{\perp }$
Study Result	0.541	0.042	/	/
Kern [27]	0.540	0.040	0.18%	5.00%
Gazzola [28]	0.550	0.040	1.64%	5.00%

show the comparison of the self-propelled swimming speeds of the fish body over time with the results achieved in previous reports [27,28]. It can be observed that the horizontal velocity $U_{\Vert }$ of the fish’s center of mass in the x-direction reaches a steady cruising speed of 0.541 after approximately 12 oscillation cycles. The translational velocity $U_{\perp }$ in the y-direction begins to exhibit periodic variations after about one cycle, oscillating around the zero point with a stable amplitude of approximately 0.042. The discrepancies between the results obtained from this study and those reported in previous works, which are no greater than 5%, robustly confirm the accuracy and viability of the model proposed in this work. In addition, the numerical results also indicate that the autonomous forward swimming process of the fish is not along a straight line but rather follows a wave-like path. It is evident that after achieving a steady cruising speed, the oscillation frequency of the $u_c$ curve is twice that of the $y_c$ curve changes. Based on the final stable swimming speed of the fish body, the calculated Reynolds number is Re = 3864, which is also consistent with the results reported by Kern [27] and Gazzola [28].

Figure 8 shows the reverse Kármán vortex street structure in the wake of a swimming fish obtained in this work. During the undulatory motion, the fish body exerts forces on the surrounding fluid. Under the interaction of fluid–structure coupling, the fluid reciprocates by applying a forward force, or thrust, to the fish body, ultimately resulting in forward motion. This simulation result realistically replicates the swimming process of a single fish in nature. As seen from this figure, during autonomous forward swimming, the downward movement of the fish tail generates a counterclockwise positive vortex, while the upward movement creates a clockwise negative vortex, forming a classic reverse Kármán vortex street structure. This type of reverse Kármán vortex street in the wake has been extensively observed in the swimming processes of fish and is considered a factor associated with thrust generation. After reaching a steady cruising speed, the biomimetic fish exhibits interactions between two vortices of the same period but in opposite directions.

4. Results and Discussion

4.1. Initial and Boundary Conditions

In nature, fish almost always swim in groups, and within a certain range, the vortex fields generated by the tails of multiple fish can interact with each other, thereby affecting the swimming performance of the fish. The size of the computational domain, displayed in Figure 9, is set as $4\times 16,$ with the corresponding grids consisting of $512\times 2048$. The body boundary is discretized into 402 points. Symmetry boundary conditions are applied to the outer boundaries of the flow field, allowing fluid to pass freely through the boundaries. The fish swims from the right side to the left side of the flow field, with a time step of the computation set to $T/1000$.

4.2. Parametric Investigations

This section primarily examines the lift and drag coefficients, as well as the pressure and vortex changes around two parallel-aligned fish, during their autonomous swimming with varying lateral spacings and different oscillation phases, given specific motion parameters.

4.2.1. Effect of Fish Spacings

This study firstly examines two parallel-aligned fish flapping in phase, as shown in Figure 10, under the following conditions: viscosity coefficient $\mu =1.4\times {10}^{-4}$, oscillation frequency $f=1$, and amplitude coefficient $A_m=0.125$.

Table 3. Computational parameters for parallel-aligned biomimetic fish at different spacings.

NO.	1	2	3	4	5
Lateral spacing $d$	0.6$L$	0.8$L$	1.0$L$	1.2$L$	1.4$L$

presents the computational parameters for different lateral spacings $d$.

During autonomous swimming with different lateral spacings, flow fields near the tails of Fish 1# and Fish 2# interact, affecting the hydrodynamic parameters of the fish bodies. Figure 11 shows the variation in the center of mass of the biomimetic fish at different lateral spacings. Overall, as the lateral spacing decreases, the flow field disturbance between the fish intensifies. At a spacing of $d=1.4L$, the centers of mass of the two fish align horizontally along the forward direction, indicating minimal interaction between their flow fields. When the spacing is reduced to $d=0.6L$, the distance between the centers of mass initially decreases and then increases. At this smaller spacing, the flow velocity between the two fish increases, causing them to move closer due to pressure. As the pressure between them builds, the coupled effects between the induced flow field and the movements of fish drive the fish apart.

As the fish bodies move forward, it is necessary to compare their resultant velocity $V_{mag}$. The resultant velocity is oriented along the direction of the fish’s forward movement and is the vector sum of the horizontal and vertical velocities. Figure 12 shows the time variation in the resultant velocity at different lateral spacings. It can be observed that when the lateral spacing is $d=0.6L$, the speed of the fish decreases as they approach each other. This reduction in speed is due to the merging of the shedding vortices behind the tails of Fish 1# and Fish 2#, which leads to a decrease in thrust.

Figure 13 presents the contour plots of pressure and vorticity around the fish bodies at $t=16T$ for different lateral spacings. Overall, as the lateral spacing between the two fish increases, the interference of flow fields between Fish 1# and Fish 2# decreases. From the pressure contour plots, it is evident that the larger spacing gives rise to less interaction between the two fish, a trend that is even more apparent in the vorticity plots on the right. When the lateral spacing is $d=0.6L$, the vortices behind the tails of Fish 1# and Fish 2# merge into larger vortices that flow backward, and two opposing pitch angles of fish are induced under the combined effects of the flow field and the autonomous swimming. Based on Figure 11, it can be seen that as the spacing gradually increases, the pitch angles of the two fish gradually decrease. When the lateral spacing is $d=1.4L$, the two biomimetic fish swim to the left along the horizontal direction with minimal mutual influence.

The energy consumed by the moving structures is another important parameter, particularly for the application of underwater biomimetic robotic fish. According to the powerful calculation strategy proposed by Albuquerque et al. [32], the energy consumption can be attributed to three aspects, including the pitching motion $E_p$, the translation motion along the x-axis $E_x$, and the translation motion along the y-axis $E_y$. In this work, the energy consumptions at spacings of $d=0.6L$ and $d=1.4L$ at $t=16T$ are estimated in

Table 4. Energy consumptions at different spacings.

	${\mathit{E}}_{\mathit{p}}$	${\mathit{E}}_{\mathit{x}}$	${\mathit{E}}_{\mathit{y}}$	$\mathit{E}$
$d=0.6L$	0.1302	0.685	0.321	1.136
$d=1.4L$	0.011	0.742	0.046	0.799

. It has already been found that the interaction between the two fish bodies at $d=0.6L$ is more significant than that at $d=1.4L$; thus, the pitch angles of the two fish are relatively larger at $d=0.6L$. Therefore, the corresponding energy consumption induced by the pitching motion at $d=0.6L$ (0.1302) is significantly higher than that at $d=1.4L$ (0.011). This table also shows that, at any given spacing, the energy consumption induced by the pitching motion is much lower than that induced by the translation motion along the x-axis and somewhat lower than that induced by the translation motion along the y-axis. Generally, the overall energy consumption at $d=0.6L$ (1.136) is much higher than that at $d=1.4L$ (0.799). The analysis of the data in Table 4 reveals that the energy consumption of underwater biomimetic robotic fish can be effectively decreased by increasing the spacing between two motion robotic fish.

4.2.2. Effect of Oscillation Phases

Figure 14 shows the in-phase and out-of-phase autonomous swimming of two parallel-aligned biomimetic fish at a given lateral spacing of $d=0.6L$ between them. The changes in speed, pitch angle, and vorticity during the autonomous forward swimming of the two fish are calculated to investigate the coupling effects of the flow fields around the two fish under different oscillation phases.

The variation in the centroid positions of the parallel-aligned biomimetic fish during in-phase and out-of-phase autonomous swimming is illustrated in Figure 15. The results show that, for both in-phase and out-of-phase oscillations at a lateral spacing of $d=0.6L$, the lateral distance between the centroids initially decreases and then increases. This behavior occurs because, during forward swimming, the fluid velocity between the two fish is higher than that outside of them. According to Bernoulli’s principle, the lower pressure in the middle causes the fish to move closer together. As the distance $d$ decreases to a critical point, the pressure between the two fish bodies increases, resulting in opposite pitch angles, as shown in Figure 16, that cause the fish to move apart again due to the interaction forces.

To further investigate the interaction of flow fields at different oscillation phases between parallel-aligned biomimetic fish, Figure 17 illustrates the vorticity fields around the fish at various time steps from $t=2T$ to $t=14T$ during in-phase oscillation. Overall, in the autonomous forward swimming of the parallel in-phase oscillating fish, the lateral spacing $d$ between Fish 1# and Fish 2# initially decreases and then increases, consistent with the trend observed in Figure 15. From $t=2T$ to $t=8T$, the lateral spacing between the biomimetic fish gradually decreases, with the vorticity generated by their tails increasingly interacting, causing the tail vorticity to come closer together. At $t=10T$, vortices with the same direction generated by both fish begin to merge, and this merging creates a reactive force that acts on the fish. At $t=12T$, a certain angle is developed between the two biomimetic fish, with the vortices in the wake merging effectively. At $t=14T$, the vorticity behind both fish is fully merged and then moves backward. The merged vorticity is significantly stronger than the initial vorticity induced by a single fish, and the two parallel-aligned biomimetic fish swim sideways at opposite pitch angles.

In comparison to the vorticity produced during in-phase oscillation, the vorticity produced during out-of-phase oscillation is more symmetrical, as shown in Figure 18. Overall, from $t=2T$ to $t=14T$, the vorticity field produced during out-of-phase oscillation is symmetric about the central axis between the two parallel-aligned fish. The trend in lateral spacing between the fish is similar to that during in-phase oscillation, decreasing initially and then increasing. However, unlike in-phase oscillation, the vorticity fields behind the tails of the two fish during out-of-phase oscillation do not merge. From $t=6T$ to $t=10T$, the vorticity fields behind Fish 1# and Fish 2# are of opposite signs and disperse backward after colliding. Due to the increased pressure field between the fish and the influence of the vorticity field behind, two opposite pitch angles of fish are generated. As shown in Figure 16, at $t=14T$, the pitch angles induced by the two fish in phase are greater than those induced by the two fish out of phase.

5. Conclusions

In order to accelerate the calculation process of fish swimming underwater, this work innovatively develops an algorithm based on the DF-IBM and implements it in an efficient, modular software program written in C++. The autonomous swimming of two parallel-aligned biomimetic fish with different spacings and fixed distances at varying phases is thoroughly studied. The main findings are summarized as follows:

(1) Smaller lateral spacing between the two fish increases the mutual interaction between the flow fields and the fish. When $d=0.6L$, the fish initially exhibit a tendency to move toward each other, followed by swimming laterally away at a certain pitch angle. Additionally, as the fish approach closer proximity, a reduction in their forward swimming velocity is noted.

(2) The flow field variations in parallel-aligned biomimetic fish with the same lateral spacing at $d=0.6L$ under different phases are also investigated. The results indicate that, regardless of their swimming phases, the fish initially exhibit a behavior of approaching each other, followed by swimming laterally away at a specific pitch angle. However, it is noteworthy that due to structural variations in the tail flow field, the pitch angle becomes more pronounced during in-phase swimming.

Generally, this study advances our understanding of the swimming mechanisms of underwater bionic fish schools, which is pivotal for elucidating collective behavior within aquatic ecosystems. It serves as a catalyst for the development of sophisticated underwater robotics tailored for applications such as environmental monitoring and search-and-rescue operations. Furthermore, the insights gained from this study are instrumental in formulating energy-efficient propulsion strategies, ultimately paving the way for innovative designs in marine vehicles. Ultimately, it is recommended that future endeavors concentrate on the intricate swimming dynamics within a fish school comprising more than three individuals.