Hydrodynamic analysis of fin–fin interactions in two-manta-ray schooling in the vertical plane

This study investigates the interaction of a two-manta-ray school using computational fluid dynamics simulations. The baseline case consists of two in-phase undulating three-dimensional manta models arranged in a stacked configuration. Various vertical stacked and streamwise staggered configurations are studied by altering the locations of the top manta in the upstream and downstream directions. Additionally, phase differences between the two mantas are considered. Simulations are conducted using an in-house developed incompressible flow solver with an immersed boundary method. The results reveal that the follower will significantly benefit from the upstroke vortices (UVs) and downstroke vortices depending on its streamwise separation. We find that placing the top manta 0.5 body length (BL) downstream of the bottom manta optimizes its utilization of UVs from the bottom manta, facilitating the formation of leading-edge vortices (LEVs) on the top manta’s pectoral fins during the downstroke. This LEV strengthening mechanism, in turn, generates a forward suction force on the follower that results in a 72% higher cycle-averaged thrust than a solitary swimmer. This benefit harvested from UVs can be further improved by adjusting the phase of the top follower. By applying a phase difference of π/3 to the top manta, the follower not only benefits from the UVs of the bottom manta but also leverages the auxiliary vortices during the upstroke, leading to stronger tip vortices and a more pronounced forward suction force. The newfound interaction observed in schooling studies offers significant insights that can aid in the development of robot formations inspired by manta rays.


Introduction
Aquatic animal swimming has inspired scientists to design robotic devices such as autonomous underwater vehicles (AUVs), and this expanding interest motivates many to study the mechanism behind thrust and lift generation across a range of swimmers.While there is considerable interest in mimicking the body-caudal-fin (BCF) mode of propulsion used by ray-finned fishes like tuna [1][2][3][4], the desire to replicate the median-paired-fin (MPF) mode of motion seen in batoid fishes, such as manta rays, has also been on the rise [5,6].In BCF swimming, fish utilize traveling-wave undulation of the body midline, adopting the caudal fin (CF) as the main propulsor [4].In contrast, in MPF swimming of batoids, the bird-like oscillatory flapping of the pectoral fins has been identified as the primary propulsion mechanism [4,7].Despite their propulsive kinematics differences, BCF and MPF swimmers have both been observed to swim in groups or schools [8][9][10].Much work has been done to understand the benefits of schooling among BCF swimmers [11][12][13], however, little is known about the benefits of schooling among MPF swimmers.This paper aims to address the intriguing knowledge gap surrounding the unique swimming styles, thrust-producing mechanism, and hydrodynamic interactions among ray-like swimmers in a school, as they may potentially derive distinct benefits from collective swimming.
Extensive research has been conducted to evaluate the hydrodynamic performance of ray-like swimmers, focusing on uncovering the fundamental thrust-producing mechanism employed by batoids.Clark and Smits developed a platform that enabled an oscillating batoid-inspired model fin that recreates some important features of a juvenile manta ray's fin.They operated this model fin at a range of Strouhal numbers while maintaining a low Reynolds number at 680 to preserve the vortex structures visualized through a dye system.They observed a wake pattern near the model fin that bifurcates in the direction of oscillation, while at the fin tip region, the tip vortices (TVs) on each side interlock each other [14].Expanding on this model by adding spanwise curvature, Moored et al introduced a high dorsoventral oscillatory amplitude to a batoid-inspired model fin, realized by a pulley system [15].Visualizing the vortex structure through dye in low Reynolds number flows, they observed a development in leading-edge vortices (LEVs) towards the tip, the formation and function of which have been studied across a range of flapping motions [16,17].
Employing computational fluid dynamics as a tool, researchers have also been able to perform hydrodynamic analysis of ray-like swimmers for models possessing various kinematic and geometric characteristics.A ray-inspired robot was simulated and studied by Liu et al., which adopts ray-inspired flapping motion operating at a constant velocity [18].Bifurcating wake in the vertical direction was reported, and higher bifurcations were observed toward the fin tip region compared to lower-order bifurcations near the basal region.A cownose ray model was implemented by Bianchi et al, with a prescribed oscillatory flapping kinematics to the 3D ray-like fin approximated as a series of NACA foils [19].The reported wake structures resembling a reverse von Karman street, in which the vortices are similar to those reported by Liu et al.Instead of bio-inspired robots, Bottom et al studied a real stingray model with features obtained by laser-scanning a freshwater stingray specimen [20].They performed simulations in two different flow velocities with corresponding kinematics in high Reynolds-number flows.From the results, LEV formation was consistently found throughout low and high swimming speeds, and they reported a higher intensity in LEV circulation strength when swimming fast.To further categorize the connections of vortex structures to hydrodynamic performance, Zhang et al performed a numerical study on a cownose ray model, with both chordwise and spanwise deformations prescribed in the flapping kinematics, immersed in a flow with Re = 2220 [21].By integrating the domain enclosing each vortex structure, they highlighted the relationship between hydrodynamic forces and the vortex structures.Specifically, they reported that the LEV mainly contributes to lift force, whereas TVs are the primary source of thrust production, and trailing edge vortices mainly induce drag.Furthermore, while maintaining a constant tip-to-tip amplitude, they introduced a certain degree of asymmetry in upstroke and downstroke amplitudes similar to the reconstruction works by Fish et al and found an increased circulation strength of LEVs associated with thrust production.Important wake vortices due to LEV development have been visualized and analyzed by Fish et al who investigated a realistic manta ray model by comparing different numerical simulations with high-resolution video footage of steady swimming manta rays in open water [22].Their findings revealed a significant asymmetry in pectoral fin flapping motion, with higher amplitude during the upstroke and a shorter duration during the downstroke.The resulting wake exhibited an asymmetric pattern, with a forward jet observed in the upstroke vortex (UV).In contrast, the vortex rings generated in downstroke produce strong backward jets, indicating the generation of high thrust associated with the production of downstroke vortex (DV).Smaller-scale vortex tubes connected the UVs and DVs, bifurcating vertically in the wake, with minimal wake dissipation observed in the lateral direction.Clear bifurcating wake structures can also be noticed in the vertical plane.
Regarding schooling behavior among aquatic animals, some of the hydrodynamic interactions and their benefits have been well-documented among BCF swimmers, through both experimental and computational approaches.Experimentally Herskin and Steffensen examined a school of juvenile sea bass, Dicentrachus labrax, in comparable sizes swimming in a large flow tank.Comparing to a solitary swimmer, they observed lower oxygen consumption rates and tail beat frequencies in the trailing fish, indicating the presence of an energy-saving mechanism through schooling [11].Similar observations on lower oxygen consumption through schooling have also been reported involving other swimmers, such as horse mackerels and grey mullets [12,23].Li et al designed and employed a pair of bio-mimetic robotic fish and experimentally studied their interaction [24].They found that by adopting tail-beat phase differences on the following fish in response to different streamwise spacing, the trailing fish exhibits similar energysaving behavior.They later conducted experiments with frees-swimming real fish using 32 pairs of goldfish and discovered a similar vortex phase matching trend.Numerical investigations of two-dimensional (2D) fish schools revealed some of the flow mechanisms that enhanced the performance of BCF swimmers.In a dense diamond school, the wall effect and block effect have been identified by Pan and Dong [13], by which the near-body vortex and pressure interactions helped enhance the average performance of the school, especially the performance of the leading and lateral fish.By applying a phase difference to the trailing fish, they also observed an increased propulsive efficiency in the trailing fish, harvesting energy from the vortex flow [25].By increasing the following distance of the trailing fish, some beneficial mechanisms can be observed by which the downstream swimmer captures upstream vortices.By placing a cylinder in a flow tank where a trout swam freely, Liao et al observed that the trout would position itself behind the cylinder and alter its body kinematics to better capture the shed vortices to improve its power consumption [26,27].Similar benefits and associated vortex dynamics have been documented in a computational study conducted by Verma et al [28].In pairs of three-dimensional (3D) 'smart' fish-like swimmers, the follower was found to position themselves in the path of the leader's vortices to save energy by utilizing the vortex-induced forward flow velocities, and a range of optimal following positions were identified in the wake of the leader.Testing different positions of paired 3D swimmers in planar stagger formation, Seo and Mittal reported a body drag reduction in the follower exceeding 10% when separated from the leader by a distance of 0.36 and 0.5 body lengths (BL) in the streamwise direction [29].By adjusting the phases of the follower, 13% and 12% increments in thrust force and propulsive efficiency were reached in their 2-fish system, featuring a LEV enhancement observed on the trailing fish.More types of formations involving different numbers of fish have been investigated in a computational study conducted by Park and Sung, who simulated finlike swimmers in various 2D arrangements, including tandem, diagonal, triangular, and diamond [30].In the tandem and diagonal formations, the follower was able to derive 14% and 6% of power saving compared to a solitary swimmer, respectively, through bodyvortex interaction, by which the follower intercepted the cores of vortices shed by the leader, and benefited by the reduced pressure created by the vortex.
The analyses mentioned above predominantly focused on BCF swimmers, with strong interactions observed in the horizontal plane where the CF flaps laterally.This suggests that the region for hydrodynamic interactions may correspond to the direction of the propulsor's flapping motion.Since manta rays flap their fins in the vertical direction, strong interactions can be expected in the vertical plane.Gao et al investigated a pair of 3D manta ray models gliding in tandem, parallel (side-by-side), and vertical formations at various attack angles [31].They discovered that in a vertical arrangement, drag reduction benefits alternate between the gliding pair: the top model experiences drag reduction when the attack angle is positive, while the bottom model benefits from drag reduction when the attack angle is negative.However, this investigation only considered gliding behavior.It did not consider the vortex dynamics as a result of the unique ray-like kinematics, and the potential of vortex-capturing interactions between two manta rays.Further understanding is required regarding how MPF swimmers can effectively utilize schooling to obtain hydrodynamic advantages.Hence, it is crucial to analyze the flow field and interactions between two manta rays with bio-inspired kinematics to establish a correlation between MPF swimmer kinematics and their schooling arrangements.
This study focuses on a two-manta-ray school, explicitly exploring the spatial positioning and phase difference of the follower relative to the leader, to advance our knowledge of fish schooling in the context of MPF swimmers.The main objective is to investigate the optimal group configuration for two manta rays in the vertical plane and analyze the underlying hydrodynamics of this arrangement.The paper is organized as follows: section 2 introduces the 3D model, school configurations, and the numerical method.Numerical simulation results and our analyses are presented in section 3 revealing the effects of stagger formations and phase differences on performance.Section 4 concludes the study, highlighting its contribution to understanding fish schooling among MPF swimmers.The results of this study provide valuable information about the hydrodynamic interactions that occur in close group configurations of manta rays.This paves the way for future research on larger groups, phase lag, and flexibility.These findings have practical applications in optimizing swarm formations and designing robots that are inspired by batoid fishes.

3D manta ray model
The 3D manta ray model adopted in this study is from our previous study [22,32], in which the model fidelity and complexity are captured accurately using to the video-captured biological motion.The kinematics is shown in figure 1(a) with fin tip trajectory over a cycle, with the downstroke motion depicted in red, and the upstroke in blue.A high asymmetry is observed between the downstroke and the upstroke, resulting in asymmetric force generation and wake features [22,32].The tip-to-tip flapping amplitude is denoted as A. The fin kinematics are broken down into bending (θ B ) and pitching angles (θ P ), expressed as a function of time in equations ( 1 The geometric quantities discussed in this work associated with the model are shown in figure 1(b), where BL is the BL, S is the streamwise separation, and V is the vertical separation, which is kept at 0.5 BL throughout this study.In the present study, the BL is used as the distance unit between the protuberance of the eyes and the end of the pectoral fins, following previous work found in [21,22,32].The schooling cases in this study are constructed by displacing the top manta-ray model (blue) in the swimming direction while the bottom manta (red) remains stationary, ensuring that both swimmers are in the same vertical (y-z) plane.Specifically, when the top manta moves downstream, it becomes the follower of a vertically staggered school, and when it moves upstream, the bottom manta becomes the follower instead.As indicated in figure 1, both manta ray models swim toward the negative z-direction, and the positive z-direction indicates the downstream direction.
A total of 18 cases are designed and described in table 1.The stacked configuration is defined as the baseline case.Vertical stagger cases are then constructed by displacing only the top swimmer in the stacked configuration downstream (positive z-direction) or upstream (negative z-direction).As stated previously, the asymmetricity in the flapping motion introduces asymmetric downstroke and upstroke wakes which inevitably will impact the interactions achieved by the follower.By displacing the top swimmer downstream and upstream, the effect of downstroke and upstroke wake on the follower can be studied.For clarity, in table 1, the cases with the top manta displaced downstream are denoted by 'D' , and the cases with the top manta displaced upstream are denoted by 'U' .The effects of phase differences are investigated by applying phase differences on the top manta in case D2, which we will show to improve hydrodynamic performance the most, and the phase-varying cases are denoted by 'PD' .

Numerical implementation and domain setup
To simulate the fluid flow with complex biological motion, an in-house finite-difference-based sharpinterface immersed-boundary method direct numerical simulation solver is employed in this work, using Cartesian-grid.The governing equation for fluids in this work is the incompressible Navier-Stokes equations, shown in equation ( 3), where u i is the velocity component, p is the pressure scalar, and Re is the Reynolds number.
The purpose of this study is to visualize the hydrodynamic interaction between two manta rays, in which the main vortex structures should be identified.To achieve a balance between attainable computational feasibility and high results fidelity, Re is set as 1200 in this study, where the inertia effect dominates [21,32].Non-dimensional flow parameters Reynolds number, Re and the Strouhal number St, are defined in equation ( 4), where ν is kinematic viscosity, and f is the flapping frequency.To maintain consistency with previous works, the Strouhal number is defined as the effective fin-tip Strouhal number, consistent with definitions given in [18,21,22].Given a fintip amplitude of 1.1 BL, free-stream velocity U ∞ of 1.25 BL/cycle was found to balance the cycle-averaged thrust and drag by Menzer et al [32], leading to a St of 0.88.
The computational domain has dimensions of 21 BL × 20 BL × 22 BL, in which the school swims in a negative z-direction.The spacing is equivalent in three spatial directions.The densest region has a grid size of ∆x = 0.008 BL, with a semi-dense region located immediately downstream of the dense area.The suitable grid spacing is realized by utilizing a topological local mesh refinement (TLMR) technique [34].This technique has been applied in previous studies on bio-inspired flights and animal swimming [32,39,40].Two structured-mesh nested blocks are employed.A larger parent block with semidense grid spacing encloses the far-field wake topology, while a smaller block with the densest grid spacing fits closely to each body to capture the near-body vortices.Figure 2(a) shows these two blocks' relative positions and the boundary conditions of the computational domain.In figure 2(a) the grid density is reduced for ease of view.In the case where the top swimmer changes location, the corresponding mesh block moves with the model.An incoming flow velocity is applied in the z-direction, and a homogeneous Neumann boundary condition is applied for pressure at all boundaries.
Grid and time-step size independent studies were conducted to determine the discretization of space and time that yielded accurate results while maintaining efficient solve times.Figure 2(b) shows the comparison of forward force coefficients (C Z ), defined as the total force in the z-direction (F Z ) nondimensionalized by 1  2 ρU 2 ∞ BL 2 , in three sets of grid sizes.Figure 2(c) shows the comparison of C Z in three timestep resolutions.The peak difference of C Z values between the nominal grid (adopted in this paper) and the fine grid is 1.8%, and the difference between the medium and high resolution for timestep is less than 2%.This demonstrates that the hydrodynamic force in the current study is grid-and timestepindependent.Hydrodynamic performance is evaluated through non-dimensional derived variables: hydrodynamic force and power consumption coefficients, and propulsive efficiency.The forces are calculated in our solver by integrating the pressure and viscous stress on the surface of a swimmer.From the previous finding by Fish et al [22], the pectoral fins, as outlined in figure 3(b), primarily generate thrust force (F T ) in the swimming direction for the manta ray.Conversely, the body predominantly generates drag force (F D ) to balance against F T .The thrust and drag coefficients, C T and C D , are therefore evaluated by normalizing F T and F D as shown in equation ( 5).Hydrodynamic power consumption, p w , is calculated as where σ is the local stress tensor on a body surface element with an area of dS, n is the normal vector of the surface element, and u rel is the velocity vector of the surface element relative to the local flow.Power consumption coefficients and the propulsive efficiency are defined in equation (5). ( In equation ( 5), CT and CPW represent the cycleaveraged thrust production and power consumption.The propulsive efficiency (η) is calculated by measuring the ratio of the pectoral fins' thrust generation to the swimmer's total power.

Results and discussions
This section presents the hydrodynamic performance for individual models in each case described in section 2.1.Section 3.1 presents a solidary swimmer for comparison purposes.Section 3.2 presents the hydrodynamics performance in the stacked baseline case.Sections 3.3 and 3.4 study the highest performance case in upper and lower wake configurations, respectively.Section 3.5 investigates the effects of the phase difference between two manta ray models in the best performance configuration.To ensure the periodicity of forces and interactions, all cases were simulated up to nine cycles.After the initial three cycles, we observed consistent periodicity in the force history.Specifically, when comparing the peak and cycle-averaged thrust, differences of less than 0.3% and 1% were found between the third and fourth cycles, and merely 0.01% and 0.46% differences were found between the fourth and the ninth cycles.This indicates that the force production and vortex interactions have reached a periodic state and are suitable to analyze.

Solitary swimmer
For comparison purposes, a numerical simulation of a solitary manta ray is performed, and hydrodynamic parameters are plotted in figure 3. Consistent with the findings by Fish et al [22], adopting viscous model, solitary manta ray model swimming exhibits dorsoventrally bifurcating wake structures (figure 3(a)), leaving key vortex structures as UVs, DVs and auxiliary vortices (AVs) interconnecting UVs and DVs.Key vortex features include LEVs and TVs, with smaller vortex tubes connecting them to UVs and DVs.The solitary manta ray model is partitioned into fin and body parts, as shown in figure 3(b).The body section is defined when near zero movement (<0.2% flapping motion) is present, which is consistent with the morphology reported by Schaefer and Summers [45] based on batoid skeletal structure.The time history of instantaneous performance coefficients contributed by both parts of the model is shown in figure 3(c

Stack configuration
This section presents the results for stack configuration following the same simulation set-up as in the previous section.In a style similar to figure 3(c), instantaneous thrust coefficients over a motion cycle are plotted in figure 4(a) for each swimmer in  stack configuration in black alongside that of a solitary swimmer in red.Performance differences are observed over a cycle.Averaging the performance differences in time, the top manta exhibits a 9.1% enhancement in the total forward force over a cycle, with a 4.2% higher thrust.The bottom manta displays a 23.1% lowered thrust relative to a solitary swimmer.From figure 4(a), the largest performance difference is observed at t/T = 0.16.The fins of the bottom manta ray exhibit a 25% lower fin thrust compared to a solitary manta ray, and the top swimmer experiences 33% higher fin thrust.To understand the performance difference, figure 4(b) presents the surface pressure (C P ) contours along with the normalized flow pressure iso-surface corresponding to C P = +0.8 , depicted in orange.Notice that the dorsal surfaces of the fin tips of both mantas are surrounded by high-pressure flow, colored orange, correspondingly the fin-tip dorsal surfaces are seen with high surface pressure.Due to the pressure-field enhancement in the narrow gap between the bottom manta's fin tips and the top manta, the bottom manta fin tips show a higher surface pressure than those of the top manta.At the instant depicted in figure 4(b), the mantas are in the early stage of the downstroke, and the fin tips encounter high pressure due to fin twisting.In the stacked configuration, the bottom swimmer's fin tips are located under the fins of the top swimmer, separated only by a narrow gap.The high-pressure region under the top manta merges with that on the bottom manta fin tip, resulting in an enhanced high-pressure region at the front of the fin tip of the bottom swimmer, exerting a backward force and therefore lowering thrust as observed in figure 4(a).On the other hand, the enhanced highpressure region results in a higher surface pressure value on the slightly backward-facing ventral surface of the top manta, which helps it to generate higher thrust.

Effect of leader downstroke wake on the hydrodynamic advantage to the follower
This section examines upper stagger configurations, which are constructed by displacing the top manta ray model in the downstream direction (case D1-D6), such that it operates in the upper wake region of the bottom leader as a follower.Figure 5 T and η * in black and red lines.From figure 5(a), the top following manta exhibits higher thrust production as it moves downstream, with a peak enhancement of 72% observed when it is 0.5 BL downstream.As its position is moved further downstream, the follower's thrust enhancement becomes less significant, and solitary performance is observed on both swimmers when it is around 0.9 BL downstream.When the separation becomes 1 BL, the follower exhibits lower thrust production.The body drag for the follower first decreases as S increases, with a minimum body drag found when it is 0.625 BL downstream and increases with higher separation.Less drastic changes are observed on the bottom leading manta, whereas its thrust production steadily increases as the top following manta moves downstream.The leader in upper stagger cases also shows generally higher body drag compared to a solitary swimmer.Figure 5(b) plots both the normalized school-averaged thrust coefficient and the propulsive efficiency (η * ) in black and red lines, respectively.Their peak values can be identified when the follower is displaced 0.5 BL downstream.The propulsive efficiency exhibits a similar trend as the thrust production.From the observations above, it can be noticed that case D2 is the best performing case.With only a cost of 15% leader thrust reduction, this formation increases the follower thrust production by 72%, while also reducing the follower body drag by 14%, and achieves an improvement in schoolaveraged propulsive efficiency by 28%.
To identify the mechanism responsible for the cycle-averaged performance improvement, instantaneous performance plots are studied to unveil the temporal development of such follower advantage.This is visualized by plotting the thrust and drag coefficients as a function of time over a typical flapping cycle, as shown in figure 6, for both swimmers in case D2, compared to those of a solitary swimmer plotted in red.Instantaneous thrust coefficients are plotted in figure 6(a), and drag coefficients are in figure 6(b).From figure 6(a), while the leader exhibits a similar time history of thrust production as that of a solitary swimmer, the follower displays a distinctive thrust production from t/T = 0.06 to t/T = 0.28.Specifically, it sustains positive thrust production even when both the solitary swimmer and the leader exhibit significant thrust valleys, at t/T = 0.16.At this moment of a cycle, the fins of the solitary swimmer have a low tilting angle, resulting in a late LEV formation which causes reverse thrust production [22].The follower, in the same motion as the solitary swimmer and the leader, does not experience the thrust valley, indicating the presence of beneficial interactions.Figure 6(b) further indicates a lower body drag experienced by the follower throughout the majority of a cycle.As for the leader, close to solitary thrust production is observed, while the body drag is reduced in downstroke and increased in upstroke, resulting in the overall increased body drag shown in figure 5(a).
Three time-instants, t/T = 0.06, t/T = 0.16, and t/T = 0.28 are selected around the first thrust peak to investigate the underlying mechanism for the follower's thrust enhancement in case D2, and the corresponding vortex dynamics are plotted to visualize the flow interaction in figure 7. Vortex topology is visualized by plotting iso-surface with Q-criterion of values at 500 in white as the vortex shell, and at 800 in blue as the vortex core.In figure 7(a), the leader sheds UV at the onset of the downstroke, though no interactions are yet observed, and thus little performance differences are present.In a later time-instant, the follower's pectoral fins intercept the UV from the leader.Clearly shown in figure 7(b), the UV from the bottom leader collides with the top follower's fin, dissipating the UV and forming a LEV on the upper surface of the follower.At the same time, a high thrust production is observed in figure 6(a).The UV-LEV interaction then dissipates at t/T = 0.28, as shown in figure 7(c).As the LEV on top diminishes and the UV is above the fin, the thrust improvement experienced by the follower decreases at this time (figure 6(a)).
To more precisely study the effect of this interaction on body forces, an x-slice parallel to the yzplane is made at the location where the interception of UV takes place at three time instants consistent with figure 7. X-vorticity (ω x ) contours are plotted in figures 8(a1)-(c1).The temporal development of vorticity is as follows: the leader sheds UV  in the early downstroke, and then it is intercepted by the fin of the follower, with a larger LEV core on the follower than on the leader.As it continues to its mid-downstroke in figure 8(c1), LEV on the follower shows similar ω x strength and thickness as those of the leader, indicating the subsiding of the effect of vortexcapturing mechanism.While the follower possesses some auxiliary x-vortices due to UV dissipation at the moment, the thrust improvement is not as significant.The corresponding pressure zone contours are shown in figures 8(a2)-(c2).A uniform highpressure region is present at the front of both manta ray models and concentrated low-pressure regions are observed in the fin tip regions for both swimmers (see figure 8(a2)).The pressure difference can be observed on the follower in figure 8(b2) at the exact location where LEV enhancements are observed.The low-pressure region behind the leader merged with the leading edge of the follower fin, enhancing the negative pressure region in front of the follower.As reported by Peng et al [46], low-pressure regions can be mutually enhancing each other when in close proximity, and hence the hydrodynamic force production of the fish.This mechanism is consistent with the observation in figure 8(b2).Compared to the leadingedge negative-pressure region in the leader, that of the follower shows a larger region with higher amplitude, which is a strong implication for suction force forward.The location of the low-pressure region collides with the LEV observed in figure 8(b1).At the lower surface region of the leading edge of the follower fin, a larger positive pressure region is also observed, enhancing thrust by exerting a forward-facing force on the fin.As the downstroke proceeds, both the leader and the follower have similar pressure contours with alike magnitudes (see figure 8(c2)).This suggests the subsiding thrust enhancement in figure 6 is connected to the dissipation of the UV as shown in figure 7(c).To examine the resultant surface pressure contours, and to ensure this interaction is not limited to the region near the slice-cut, surface pressure (C P ) contours are shown in figures 8(a3)-(c3).Figures 8(a3) and (c3) show no major surface pressure contour differences between the leader and the follower, which are consistent with the similar pressure zone shown in figures 8(a2) and (c2), respectively.A noticeable surface pressure contour difference is shown in figure 8(b3) for the follower.A large region of negative surface pressure is observed on the upper surface of the pectoral fins at the leading edge, which is consistent with the pressure zone slice depicted in figure 8(b2).The follower fin's low-pressure region extends through the mediolateral direction indicating that the low-pressure zone shown in figure 8(b2) is not limited to the x-plane slice shown.Hence, the vortex-interaction-based LEV enhancement shown in figure 8(b1) extends along the fin, likely coinciding with the low surface pressure region in figure 8(b3).
To quantify the strength of UV-LEV interaction shown in figure 8(b1), a series of equally spaced slices are made along the fin leading edge at the time instant t/T = 0.16 as denoted by the group of red lines in figure 9(c) to show the LEV core.This time-instant corresponds with the large thrust enhancement in figure 6(a) and spatially aligns with the low-pressure regions found in figure 8(b3).Figures 9(a) and (b) shows contours of ω x on these slices, as ω x is assumed to be dominant in the LEV generation [18].Vortex boundaries are shown in transparent white with a Qcriterion value at 400.

Effect of leader upstroke wake on the hydrodynamic advantage to the follower
This section examines the lower stagger configurations, which are constructed by displacing the top swimmer upstream from the stacked formation (case U1-U6).Thus, the bottom manta becomes the follower of the top leading manta, and we now study the interactions present for a follower swimming in the downstroke wake of the leader.Figure 10(a) summarizes the normalized cycle averaged thrust and drag coefficients, shown in black and red lines, of each swimmer in all lower stagger cases.The normalized school-averaged thrust coefficient and propulsive efficiency over a motion cycle is plotted in figure 10(b).
In figure 10(a), the thrust coefficients as a function of S for the follower display a parabolic shape that resembles that of the top follower in figure 5(a).As the streamwise separation increases, the following manta gains higher thrust with a maximum of about 1.31 in case U2, as well as lower body drag with a minimum of 0.87 in case U3.Thrust improvement on the follower decreases when the separation distance increases beyond 0.5 BL, and lower than solitary performance is displayed when S is higher than 0.77 BL.The leader's thrust production monotonically increases as S increases although all of them are slightly lower than those of a solitary swimmer except when S = 1 BL.Schooling advantages are obtained when S is less than 0.75 BL with higher school-averaged propulsive efficiency, and the maximum propulsive efficiency and forward force coefficient are found in case U2.
Next, case U2 is studied with instantaneous (a) thrust and (b) body drag coefficients of each swimmer plotted in figure 11 alongside those by a solitary swimmer, shown in black and red lines, respectively.Clear concordance of force production by the schooling and solitary mantas can be observed during the downstroke, and a major performance difference is shown in the upstroke from t/T = 0.66 to t/T = 0.81, when the follower generates higher thrust.Specifically, a local thrust peak is observed at t/T = 0.75, contrary to a solitary swimmer found to be approaching its second thrust peak after a local minimum at the same time.This local thrust peak for the schooling manta at t/T = 0.75 indicates the follower obtains its most significant thrust increase during the upstroke.The thrust enhancement dissipates as they reach the end of a motion cycle when both swimmers exhibit thrust production like the solitary swimmer.As for the body drag coefficient shown in figure 11(b), both swimmers show higher body drag during the downstroke, and in the upstroke, the follower shows lower body drag throughout.The body drag differences between the follower and the solitary swimmer are more significant in upstroke than in downstroke, which results in the lower body drag shown in figure 10(a).
To study the mechanism responsible for thrust enhancements exhibited by the following manta, three time-instants during the upstroke, at t/T = 0.66, t/T = 0.75, and t/T = 0.81, are examined.Vortex dynamics are visualized in figure 12 by plotting the iso-surfaces for Q-criterion values of 200 and 350 for the vortex shell in the core in white and blue, respectively.Shown in figure 12(a), the leader sheds a DV that propagates downstream, and is then captured by the follower's fin in a later time instant as in figure 12(b).At t/T = 0.75, a large LEV is observed on its fin, while the DV from the leader dissipates greatly.Vortex separation is found in the proximal direction, while attachment is found in the distal direction near the tip.This is caused by the location of DV.The strong circulation in DV is intercepted by the follower's pectoral fin, forming strong LEV towards the fin tip, while the downstream DV propagation tears off the proximal LEV.The vortex separation morphs into a vortex tube (VT), peals the tip portion of the LEV off the fin tip, and then develops into the dual VT shown in figure 12(c).
In a style similar to figure 8, an x-plane slice where the DV interception occurs is made in the same three time instants.The x-vorticity ω x contours are plotted through figures 13(a1)-(3).Both swimmers exhibit similar ω x contour in figure 13(a1), with the DV clearly identified for the leader.At a later instant, the DV from the leader collides with the lower surface of the follower's fin.The counter-rotating DV and LEV results in an upward jet, strengthening the LEV on the following manta.The intercepted part of DV maintains high vorticity as shown in figure 13(c1), resulting in the persistence of high-strength LEV.This is consistent with figure 11(a) at t/T = 0.83 when the follower exhibits a second thrust improvement peak.The corresponding fluid pressure contours are shown in figures 13(a2)-(c2).In figure 13(a2), the weak positive-pressure regions and strong negativepressure regions were found respectively on the dorsal and ventral surfaces of both the follower and the leader, exerting a forward-facing force.As the DV is captured by the follower, it enhances the fin LEV as shown in figure 12(b), and correspondingly the negative-pressure region on the leading-edge lower side of the bottom follower was enhanced, as shown in figure 13(b2).In figure 13(c2), the low-pressure region for the follower is more significant than that for the leader, with a contour strength and shape similar to the LEV shell shown in figure 13(c1).Surface pressure contours (C P ) on the lower surface at the corresponding time instants are plotted in figures 13(a3)-(c3).The enhanced LEV observed for the follower at t/T = 0.75 exerts a suction force with enhanced magnitude and extending over more surface area, compared to that experienced by the leader, as shown in figures 13(b3).An enhanced negative surface pressure contour on the follower is also observed in figure 13(c3).The LEV on the follower is enhanced by the DV from the top manta, forming a low-pressure region in the front of its pectoral fins and enhancing its tip low-pressure region.
Next, a series of x-plane slices are made along the mediolateral direction on the fin tip at t/T = 0.75 to quantify the strength of DV-LEV interaction as performed for the UV-LEV interaction in figure 9. Slices are made from mid-span (0.68 BL) to the fin tip (0.96 BL), where major interactions and LEV enhancements can be observed (figures 12 and 13).Figure 14 displays the ω x contour slices in the spanwise direction for the top leading manta and the bottom following manta in figures 14(a1)-( 2), and for the solitary swimmer in figure 14(b) as a comparison.Consistent with the previous observation, the leader exhibits similar x-vorticity along the LEV to a solitary swimmer.Figure 14(c) shows the normalized circulation strengths along the mediolateral direction for a solitary swimmer, and both swimmers in case U2.Along the mediolateral direction on the follower's fin, the vorticity strength first decreases and then increases before a sharp cutoff towards the tip.From figure 14(a2), this trend can be correlated.As the tear-off LEV merges to the tip vortex, the vortex strength reduces as the vortex cross-section shrinks, and the later vorticity strength increase is from the enhanced tip vorticity.The trend at the fin tip region, from 0.8 BL to 0.96 BL, is seen among all swimmers, but the follower exhibits higher fin tip vorticity magnitudes, indicating its gaining an advantage.Sharp cut-off of vorticity strength is present for all swimmers at x = 0.96 BL, where the tip vortex sheds and, therefore, lower ω x strength is present.

Effect of phase difference
It can be concluded from sections 3.3 and 3.4 that among all schooling cases, case D2 is the bestperforming schooling configuration as the UV-LEV interaction is significant.To study the robustness of the thrust enhancement mechanism herein, phase shift is applied to the follower of case D2, and the results are presented in this section.Without any phase shift, in the flapping motion of a manta ray model, the onset of downstrokes indicates the beginning of a cycle motion.A phase lead (ϕ > 0 • ) indicates the flapping motion preceding the baseline cycle motion.In phase-different cases, the phase of the follower in case D2 is adjusted by an increment of π/3 (60 • ), and the normalized cycleaveraged hydrodynamic performance coefficients are shown in figure 15(a).Consistently in the formation of D2, a thrust enhancement is experience by the follower regardless of its phase relative to the leader.A follower thrust enhancement even higher than that of D2 is achieved when the follower exhibits a 60 • phase lead (case PD1), in which case the follower thrust production reaches a peak 2.3 times higher than that of a solitary swimmer, while its body drag is also slightly reduced.The phase difference also brings slight thrust enhancement for the bottom leading manta.and red lines, respectively.Comparing the thrust production plots of the follower of case D2 to that of case PD1, the one from case PD1 exhibits two thrust peaks of similar magnitudes, in both its upstroke and downstroke, instead of just a singular downstroke peak.A downstroke thrust enhancement is observed at t/T = 0.16, the mechanism of which has been discussed in section 3.3.This indicates the robustness of the UV-LEV interaction.Further, the top manta ray gains a second thrust enhancement in upstroke, which is unprecedented in the previous cases.The high magnitude of upstroke thrust production at t/T = 0.65 is shown to be comparable to the downstroke thrust production.The thrust improvement in the upstroke peak in figure 15(c) is then studied.
To identify the mechanism responsible for the enhanced upstroke thrust production in case PD1, the time instant when the follower exhibits the second thrust peak, t/T = 0.65 in figure 15(c), is examined.The followers from both cases D2 and PD1, at the same time-instant of their motion cycle, are displayed in figure 16.Iso-surfaces are used to visualize flow interaction with Q-criterion values of 500 for the vortex shell in white, and 800 for the vortex core in blue, as shown in figure 16(a).It is observed that the follower from case PD1 exhibits a more significant LEV attachment than that of case D2 does, whereas the LEV of the latter appears to be scattered (see figure 16(a2)).Such LEV attachment for the follower of case PD1 indicates a low-pressure region is likely to be present around its fin tip region.Figure 16(b) plots the surface pressure C P contour on the lower surface of the pectoral fins of the follower in both cases.Comparing the surface pressure contour between figures 16(b1) and (b2), a larger negative-pressure region can be observed at the fin tip region, with a higher strength of the negative-pressure zone.This surface pressure difference is consistent with the vortex strengths observed in figures 16(a1) and (a2).This indicates that when the phase shift is applied on the  follower, it can maintain the LEV attachment during upstroke, in comparison to the scattered vortices observed in case D2.
To reveal the mechanism enhancing the LEV attachment on the follower after a phase shift, an xplane slice is made amid the LEV to show the fluid pressure contours in figure 17, for three time-instants, (a) t/T = 0.50, (b) t/T = 0.65, and (c) t/T = 0.82.In figure 17(a), a shared high-pressure region connects the fin tips of the leader and follower at t/T = 0.31.At this moment, the leader has undergone half of its downstroke exhibiting instantaneous peak thrust coefficient, with a high local pitch angle so that its dorsal surface faces downstream, whereas the follower transitions from downstroke to upstroke with fin tips still descending.Figure 17(b) shows a similar fluid pressure contour as a shared negative-pressure region is present between two swimmers.The bottom leader is starting its reversal at the end of the downstroke, as the fin tip sheds a row of AVs that reduce flow pressure in the wake.At that exact moment, the top follower amid its upstroke forms LEV.As the pressure in the vortex core is lower than the surrounding fluid, a local low fluid pressure field is present.The adjacency between pressure fields created by the AVs from the leader and the LEV of the follower results in the enhancement of leading-edge suction, enabling thrust improvement on the follower.In figure 17(c), little interaction is observed as both swimmers are in upstroke, forming a high-pressure region on the upper surface of the fin and a low-pressure region on the lower surface of the fin respectively.With no adjacent pressure region of the same sign, no enhanced pressure region is observed, and therefore little interaction is found at this moment.

Conclusions
In this study, we employed an immersed-boundarymethod flow solver to investigate the hydrodynamic performance of a two-manta-ray school in the vertical plane, focusing on the effects of spatial arrangement and phase difference.Two manta ray models were arranged in a stacked formation, and then the top manta ray was displaced in the streamwise direction to create vertical stagger formations.We also introduced a phase difference to the configuration that exhibited highest thrust enhancement.Through the computational analysis, we observed that the baseline case of the stacked formation had a detrimental effect on the bottom swimmer.During the early stage of the downstroke, the fin tip of the bottom manta was in close proximity to the lower surface of the top manta, resulting in an enhanced positive pressure region between them.This enhanced high-pressure region, while slightly reducing drag on the top manta, simultaneously reduced the thrust on the bottom one.By displacing the top manta by less than one BL either downstream or upstream from the stacked formation, we achieved drag reduction and thrust enhancement in various stagger formations.The vortex-capturing mechanism behind these improvements was examined.The follower always exhibited the largest variation in performance and achieved higher thrust enhancement, regardless of whether it is above or below the leader.
Among the different configurations tested, the most pronounced interactions were observed when the streamwise separation was 0.5 BL (case D2 and case U2), where both followers achieved notable increases in cycle-averaged thrust of 72.0% and 30.6%, respectively.When the top manta ray model moved downstream and assumed the follower position, it experienced improvements during the downstroke by capturing the UV shed by the leader.As confirmed by circulation calculations, this resulted in stronger downstroke fin-root LEVs than those observed in the leader and solitary manta rays.The intensified circulation of the LEV increased the forward suction force on the fins of the top manta, leading to enhanced thrust generation.
Conversely, when the top manta ray model moved upstream, the bottom manta ray became the follower.In this scenario, the follower benefitted during the upstroke by intercepting the DV shed by the leader.This resulted in the formation of enhanced LEV near the fin tip region, as verified by circulation calculations.The presence of this strengthened LEV enabled the follower to generate higher thrust compared to the leader.Comparing case D2 to case U2, it was observed that while the follower could benefit from both UV and DV interception, the advantage obtained from UV interception was more significant than that from DV interception.The interesting observation that the optimal streamwise spacing remains 0.5 BL for both UV and DV interception, is attributed to the advection of the vortices and the positions of the fin tips of the follower relative to the leader.UVs and DVs from the leader bifurcate at similar rates and are carried downstream by the same incoming flow.Thus, a staggered follower can capture the UVs and DVs optimally at similar streamwise spacings from the leader.
To explore the effect of phase differences on thrust enhancement, phase differences were incrementally applied to the optimal thrust enhancement case (case D2) in intervals of π/3.The highest level of thrust enhancement was observed when the follower led the leader in phase by π/3, resulting in a cycle-average thrust production 33% higher than that of the follower in case D2.Further analysis of the follower's thrust production time history revealed a unique upstroke thrust peak previously unseen, in addition to the downstroke enhancement in case D2 (as discussed in section 3.3).The vortex dynamics demonstrated a strengthened upstroke LEV near the fin tip region, attributed to the capture of AVs.This resulted in a better attachment of LEV to the fin tips, accompanied by a large negative surface pressure value on the fins, and therefore an enhanced forward suction on the follower in the ϕ = π/3 case.
These findings provide valuable insights into the hydrodynamic benefits of schooling in batoid fishes, bridging the gap between oscillatory swimming and schooling behavior.Moreover, they unveil the underlying mechanisms responsible for these benefits in manta rays.This increased understanding serves as a foundation for the strategic deployment of multiple bio-inspired manta-like robots in collective operations, opening new possibilities in the field of underwater robotics.

Figure 1 .
Figure 1.(a) Oscillatory kinematics of the manta ray model with superimposed fin tip trajectory of a cycle.(b) Schematics illustrating the definitions of body length (BL) and spatially related parameters.The negative z-direction aligns with the direction of swimming.

Figure 3 .
Figure 3. (a) Side view of the 3D wake structure at t/T = 0.625 for the manta ray model, visualized through iso-surface by Q-criterion of values at 5 for vortex shell in white and 40 for vortex core in blue.(b) Fin and body definition of a manta ray model, and (c) their thrust and drag coefficients in black and red lines, respectively.Only one of the fins is plotted due to symmetry.
), in which the black line indicates the C T and the red line indicates the body C D .The grey region, extending 0.46 of a motion cycle, indicates the period of the downstroke, while the rest indicates the upstroke (figure 1(a)).

Figure 4 .
Figure 4. (a) Instantaneous thrust coefficients of a solitary swimmer in red line and both swimmers in stacked formation in black lines, observed over one cycle.(b) Surface pressure (CP) contours of two manta rays, with orange iso-surface indicating flow regions corresponding to CP of +0.8.

Figure 5 .
Figure 5. (a) Normalized cycle-averaged thrust and drag coefficients for both swimmers in upper stagger cases, and (b) school-averaged thrust coefficient and propulsive efficiency.The dotted lines are for those of a solitary swimmer.
(a) summarizes the cycle-averaged thrust and drag coefficients, C * T and C * D in black and red lines respectively, for each individual swimmer normalized by the those of a solitary swimmer.Figure 5(b) shows the normalized school-averaged thrust coefficient and the propulsive efficiency, C *

Figure 6 .
Figure 6.Instantaneous (a) thrust and (b) drag coefficient plots as a function of time for both swimmers in case D2 (top following manta moved 0.5BL downstream) plotted in solid and dotted lines, and for a solitary swimmer plotted in red.

Figure 9 .
Figure 9. Eight spaced x-slices showing the dominant vortex formations for (a1) the leader and (a2) the follower in case D2, and (b) the solitary swimmer for comparison.The circulation strengths of the vortices on the slices are plotted in (c).

Figure 9 (
c) plots the circulation strength (|Γ x |) nondimensionalized by the product of free-stream velocity U ∞ , and BL, against mediolateral location away from the centerline, nondimensionalized by BL.The strength of the circulation of the LEV at these specific locations are calculated as Γ x = ˜S ω x dS, where ω x is the x-vorticity at the given time instant, and S is the cross-sectional area of the vortex core with |ω x | > 70.Figures 9(a1)-(2) represents the leader and follower in case D2, and figure 9(b) represents the solitary swimmer, respectively.The general trend of larger LEV cores is observed on the follower, while those of the leader and a solitary swimmer are similar.The vortex strength plot in figure 9(c) confirms the observations in figures 9(a) and (b), that the circulation strengths for the solitary swimmer and the leader are close throughout the span.The follower's leading-edge vortex strength is stronger than that of the leader across the measured span length.Specifically, a steady increase up to 0.6 BL away from the centerline is found, consistent with the region of low-surface pressure in figure8(b3).The UV-LEV interaction leads to a stronger vorticity on the follower, and therefore larger LEV-induced negative pressure region thereby enhancing the follower's propulsive performance.

Figure 10 .
Figure 10.(a) Normalized averaged thrust and drag coefficients for both swimmers in lower stagger cases, and (b) school-averaged forward force coefficient and propulsive efficiency.The dotted lines are for those of a solitary swimmer.

Figure 11 .
Figure 11.(a) Thrust and (b) body drag coefficients for each swimmer in case U2 (top manta moved 0.5BL upstream) in black lines.Red lines indicate those for a solitary swimmer.
Figure 15(b) plots the school averaged thrust coefficient and propulsive efficiency, plotted in black and red curves, respectively.Both curves peak at a phase difference of 60 • , which indicates that the thrust enhancement obtained from the phase shift does not require extra power consumption.Figure 15(c) shows the instantaneous thrust coefficients of the followers in cases PD1 and D2, and of a solitary swimmer, represented by blue, black,

Figure 14 .
Figure 14.Eight equally x-slices made on the fin across the mediolateral showing dominant vortex formations on (a1)-(2) the leader and follower in the case U2, and on (b) the solitary swimmer.The vortex circulation strengths on the slices are plotted in (c).

Figure 15 .
Figure 15.(a) Normalized cycle averaged thrust and drag coefficients for both swimmers as a function of phase difference, and (b) instantaneous thrust coefficients by the solitary swimmer (red), and the follower from case PD1 (blue), and case D2 (black).

Figure 17 .
Figure 17.Lateral view on a flow slice showing the temporal development of pressure contours at (a) t/T = 0.31, (b) t/T = 0.5 and (c) t/T = 0.66.

Table 1 .
Description of cases.