12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022

HYDRODYNAMICS OF A COW-NOSED RAY (RHINOPTERA JAVANICA) IN
FORWARD SWIMMING AT A HIGH REYNOLDS NUMBER

Dong Zhang
Department of Engineering Mechanics
Tsinghua University
Beijing 100084, China
email: dong zhangyx@126.com

Wei-Xi Huang
Department of Engineering Mechanics
Tsinghua University
Beijing 100084, China
email: hwx@tsinghua.edu.cn

ABSTRACT(2016) simulated the swimming of stingrays at two different
swimming speed. The fast-swimming stingray was found 12% more efficient than the slow-swimming one. It was the firsttime to show the existence of LEV on stingrays and the LEV
was related to the thrust and efficiency enhancement. Zhang
et al. (2022a) performed a numerical study on the hydrodynamic enhancement mechanism in batoid fish swimming.
LEV, TEV and tip vortex (TV) were observed to construct themain vortex structures. By using derivative-moment transformation (DMT)- based force analysis, at low speed, TV contributes most to the thrust, whereas the LEV and TEV are unfavourable to the thrust generation. On the contrary, at a fast
speed, the role of LEV changes to a favourable effect on thethrust. The effect of wavenumber was shown to influence theeffective angle of attack of the fins. A smaller wavenumber
leads to a larger effective angle of attack, and thus the strongest
vortices. However, wavenumber of 0.4 has an optimal balance
between the thrust enhancement and weakness effect of different local structures. Also, works focusing on more topics
of batoid fish can be found, such as the kinematics (Fish et al.,2016), flexibility (Menzer et al., 2022), the role of cephalic fins
(Zhao et al., 2021) of manta rays and hydrodynamic benefitsof protruding eyes and mouth on stingrays (Mao et al., 2021).
Numerical simulation is an important method in investigating the hydrodynamics of aquatic animals. However, thehigh Reynolds number (Re) associated with the complex bodymorphology and swimming kinematics are the major challenges for numerical simulation because it requires very finegrids to resolve the boundary layer and the flow. Most of theabove studies used low Reynolds numbers to present the realistic swimming flow. However, batoid fish swims at higher
Reynolds numbers. For example, cow-nosed rays swim at
Re = 5 × 104 - 3 × 105 (Zhang et al., 2022a). Manta rays swim
at Re = 6.85 × 105 − 7.71 × 106 (Fish et al., 2016). These motions have located in the turbulent regime for the flow over aflat plate. Moreover, we have concluded the Reynolds number of recent works about more species and models in figure 1(Zhang et al., 2022b) and found that the high-Re investigation
is rare, especially the numerical work. Zhang et al. (2018a,b)
conducted some high-Re simulations, but 2D models were em-

Numerical simulation is an important method to investigate the hydrodynamics of aquatic animals. However, the highReynolds number (Re) simulation is a major challenge for suchmoving boundary problems and most work is limited the Re
below 104 . In the present study, we investigated the hydrodynamics and flow field of a forward swimming cow-nosed ray(Rhinoptera javanica) at Re = 1.48 × 105 . The computationalmodel was constructed based on the biological data. A hybrid
immersed boundary/wall-model approach for large eddy simulation was employed to simulate the high-Reynolds number
turbulent flow around the cow-nosed ray. The hydrodynamic
characteristics and vortex evolution of the swimmer were investigated. The effects of Re were also analyzed.

INTRODUCTION
Among aquatic animals, batoid fish are known for theiroutstanding swimming abilities, such as high-speed swimming, high-efficiency cruising, and maneuverability. Batoid
fish are special aquatic animals that have a dorsoventrally flattened body and expanded pectoral fins. They adopted pectoralfin-based locomotion as their sole propulsion strategy. The
hydrodynamic mechanisms behind their high-speed and highefficiency propulsion are of particular interest to the design ofhigh-performance underwater vehicles.
Batoid fish is known to have various locomotion modes
from undulatory to oscillatory (Rosenberger, 2001). One parameter to distinguish the two modes is the wavenumber travelling on the fins. Undulatory motion is defined to have morethan one wave present on the fins at a time, whereas oscillatory motion are thought to have less than half a wave. The
typical species of the former include skates and most stingrays.
Manta, cownose, eagle and bat rays employ the latter mode.
Dewey et al. (2012) studied the hydrodynamics of thebatoid-inspired oscillating propulsion using a mechanical flexible fin at Re = 11600. DPIV was used to display the wake andthe interconnected vortex rings was observed in a 2D form,i.e., leading edge vortices (LEVs) and trailing edge vortices(TEVs) become entangled with one another. Bottom Ii et al.

1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022

x(x f , y f ,t) = x f
y(x f , y f ,t) = y f (1 − (1 − k)| f (x f ,t)|y f /SL)cosθ (x f , y f ,t)z(x f , y f ,t) = z f + y f (1 − (1 − k)| f (x f ,t)|y f /SL)sinθ (x f , y f ,t)(1)

θ (x f , y f ,t) = θmax y f /SL · f (x f ,t)f (x f ,t) = sin(ωt − 2πW x f /BL)

(2)

where (x f , y f ) defines the neutral plane of the pectoral fins ina flat position (figure 2a). At time t, the shape is describedby the deformed coordinates (x, y, z). W = BL/λ is the nondimensional wavenumber, which is used to present the chordwise deformation. λ is the wavelength of the chordwise travelling wave. ω is the angular frequency. In order to fit the spanwise deformation, k and θmax are used to control the parameters, whose meanings can be explained by figure 2b. It is ası withsumed that the spanwise curve can be treated as an arc oT
a radius of l. θ is the deviation angle, defined by the straightline oT joining the fin root and tip. The parameter k means theı , i.e., k = 2lsinθmax /(l · 2θmax ),ratio of oT to its arc length oT
when the fin tip reaches the highest position. θmax is obtainedby substituting y f = SL into the expression of z(x f , y f ,t) inequation (1), by providing the amplitude A. In our simulations, for the amplitude of A =0.6BL, a value of θmax = 0.41
and k = 0.972 can be determined for the model.
In order to simulate the high-Reynolds number turbulent
flow over the complex moving boundaries, a hybrid immersedboundary (IB)/wall-model approach for large eddy simulation
is employed. Details about this approach can be found in ourprevious work of Ma et al. (2019), so only brief descriptions
are introduced here. The flow solver is governed the by filtered incompressible Navier-Stokes equations with an additional momentum force fi .

Figure 1. The Reynolds numbers used in recent works ofswimming and flying; the circle represents high-fidelity numerical work (labelled by N ) and the fork represents experimental work (labelled by E). Adapted with permission from
Zhang et al. (2022b).

Figure 2. (a) The cow-nosed ray model in a flat position andthe coordinate system for the kinematic model; (b) The spanwise deformation of the pectoral fin in the front view at different snapshots.

∂ ũi ∂ ũi ũ j∂ p̃ ∂ τ̃i j
1 ∂ ũi
+
=
−
+
+ fi
∂t
∂xj
∂ xi ∂ x j
Re ∂ x j ∂ x j

ployed in their work. To our knowledge, most numerical workof aquatic animals were limited the Re below 104 . It is difficult to obtain the flow field information, which is significantly
important in understanding the flow mechanisms of aquatic
swimming, at such high Re. In the present work, we conducted a simulation of the flow of a cow-nosed ray swimmingat Re = 1.48 × 105 .

∂ ũi
=0
∂ xi

(3)

where ũi is the filtered velocity components, p̃ is the filteredpressure. The large-scale structures of the turbulent flow field
are directly resolved by the grid as the SGS structures andmodelled using the Smagorinsky model. The additional force
in the governing equations is calculated by a penalty method(Huang et al., 2011) on the Lagrangian IB points.
To implement the wall model, the thin boundary layer
equation was solved on an embedded mesh refined along thewall-normal direction and a dynamic matching procedure wasadopted. This solver has been successfully applied in simulations of flapping propulsion, and validations can be also foundin the previous studies (Ma et al., 2019, 2021).

COMPUTATIONAL MODELA three-dimensional (3D) model of the cow-nosed ray
was constructed based on several section profiles, which canbe seen from the previous study (Huang et al., 2020). By therecordings of free-swimming cow-nosed ray specimens at an
aquarium using a high-resolution digital camera, fin kinematic
data were analysed. Some key points, such as the swimmers’
head, tail base and pectoral fin tip were recorded. It was observed that the main deformation utilized by cow-nosed rays
can be decomposed into spanwise deformation and chordwise
deformation. We used a rotation around the longitudinal axis
(x-axis) of the body and a chord-wise travelling wave to mathematical describe the above deformations, respectively. The
mathematical approximation is introduced as follows (Zhang
et al., 2022a),

RESULTS
We conducted the simulations at the Reynolds number fora realistic cow-nosed ray , Re = 1.48 × 105 . In order to detect some differences at such a high Reynolds number, anothersimulations were also performed at Re = 1.48 × 104 for comparison. The frequency f and tip-to-tip amplitude A are set at2

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022
0.7 Hz and 0.6BL, respectively, resulting in a Strouhal numberof 0.39. It is well known that many swimming animals are observed to cruise in a narrow range of Strouhal numbers (0.2 ∼0.4) (Triantafyllou et al., 1991). Such parameters are also in
the range of those observed in most batoid fishes (Rosenberger,
2001; Zhang et al., 2022a).

Hydrodynamic force
Figure 3a presents the instantaneous thrust which was
normalized by 0.5ρu2 BL2 . Around the mid-upstroke (t/T =0.3) or mid-downstroke (t/T = 0.8), it produces a thrust peak,whereas it produces a drag peak as the fin reverses its stroke.
By decomposing the total force into pressure (figure 3b) andviscous components (figure 3c), we found the main sourceof the thrust comes from the pressure force. Compared withthe hydrodynamic performance at a lower Re (1.48 × 104 ),the thrust improves significantly with the increase of Re. The
thrust enhancement comes from both the pressure and viscous
force. From figure 3b the pressure force peak is enhanced by
83 %, and the viscous force decreases by 60 %. Those are
due to that the high Re condition changes the flow around theswimmer, which is analysed in the following section. Notably,
although the force value varies along with Re, the temporalvariation tendency changes little. This coincides with previous study (Bozkurttas et al., 2009) as the Reynolds number
increases from 540 to 6300, for a bluegill sunfish. Our result
indicates that this tendency remains the same for the flow at
the Reynolds number of the order of 105 .

Vortex dynamics
Figure 4 displays the vortex evolution of the swimming
cow-nosed ray at Re = 1.48 × 105 , visualized by the isosurface of Q = 500 and colored by streamwise velocity. Figures
4a - 4c denote the downstroke, whereas figures 4d - 4f denote the upstroke. At t/T = 0.17, the boundary layer separates at the leading edge and formed a leading-edge vortex(LEV). The LEV structure breaks up into several co-rotating
vortical structures. Similar phenomenon was also detected inrevolving insect wings (Chen et al., 2017). However, there
form more vortices in the present work. At t/T = 0.33, as the
LEVs convect downstream, the vortex cores are less diffused,especially at the tip region. At t/T = 0.5, there are no obvious coherent LEV structures, associating with the shed of tip
vortices (TVs). Small-scale turbulent structures are likely to
be dominant at the posterior of the cow-nosed ray and it forms
hairpin vortices in the wake. At the dorsal of the cow-nosed
ray, it also generates the hairpin vortices as seen in the enlarged view of figure 5. Two hairpin vortices packets, one containing five hairpins, another containing two hairpins, can beclearly observed. Such structures are usually detected in wall
turbulence (Adrian, 2007). The existence of the hairpin on
aquatic animals indicates that the flow is turbulent at the realistic Reynolds numbers. During the upstroke (t/T = 0.67 − 1),
the broken LEVs keep convecting to the wake. Meanwhile, therepeated process takes place at the lower side of the pectoral
fins.
For comparison, figure 6 gives the vortex structures at
Re = 1.48 × 104 . The LEVs are also observed to break into
several structures. However, larger-scale coherent structures
dominate the flow, even in the wake when the structures areshed from the surface. Notably, at the dorsal of the cow-nosed
ray, no hairpin vortices are formed.
Figure 7 displays the pressure distribution at the two
Reynolds numbers during the downstroke. The attachment of

Figure 3. Time histories of the (a) total thrust force, (b) pressure force component and (c) friction force component of aswimming cow-nosed ray at Re = 1.48 × 104 and 1.48 × 105 .

LEV produces a low-pressure region near the leading edge.
We can observe the moving of such low-pressure region withthe convection of LEVs. The multiple LEVs create a long narrow band, which can be seen in figures 7a and 7d. At the dorsalof the surface, the existence of hairpin vortex produces someshort low-pressure bands at Re = 1.48 × 105 , whereas suchband does not occur at Re = 1.48 × 104 . Another difference isthat it produces a stronger pressure fluctuation, which gives theevidence that the flow at such Re is turbulent. To describe thepressure difference due to the Reynolds number, the spanwise
vorticity contours and pressure distribution along the chord at
the plane of z/SL = 0.31, 0.62 are shown in figure 8. As shown3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022

Figure 4. The instantaneous vortical structures of the swimming cow-nosed ray at Re = 1.48 × 105 , visualized by the isosurface of
Q = 500 and colored by streamwise velocity. ur e f is the reference velocity at incoming flow.

in figure 8a, the attachment of LEV is truly observed to generate a pressure valley, corresponding the low-pressure regionin figure 7 and the valley value is smaller for the high-Re flow.
From the figure 8b, we can see that the pressure difference between the leading edge and trailing edge is larger for the highRe case, creating a higher thrust force. Moreover, the stronger
fluctuation shown in figure 8b provides visible evidence of thetransition from laminar to turbulence.
From the shear distribution (figure 9) of the surface, it apparently shows that the high shear region locates at the leadingedge of the model at Re = 1.48 × 105 . At Re = 1.48 × 104 ,almost half of the whole body are occupied by the high shear
region and the shear peak is higher than that at high Reynolds
number, which gives an explanation to its high viscous force.

Figure 5. The enlarged view of vortical structures to show
the multiple LEVs and hairpin vortices at Re = 1.48 × 105 .

CONCLUSIONSIn the present study, we investigated the hydrodynamics and flow field of a forward swimming cow-nosed ray(Rhinoptera javanica) at Re, as high as 1.48 × 105 , which isin the Re range of real swimming. We constructed a computational model based on biological data. A hybrid immersed boundary/wall-model approach for large eddy simulation was employed to simulate the high-Reynolds number
turbulent flow. We traced the evolution of the LEV and foundthe hairpin vortices at the dorsal and wake of cow-nosed rays.
The LEV to hairpin vortices transformation suggested a laminar to turbulent flow transition, whereas it does not take place
at 1.48 × 104 . Such transition leads to a change in the pressure
and shear force distributions. As a result, the pressure component of the thrust force improves and the viscous component
decreases; thus, The thrust force is enhanced as the Reynolds
number increases. As stated in the introduction, the Reynolds
number is limited to 104 for most previous studies. The present
work gave an evidence that the force at low-Re condition canreflect the real hydrodynamics, at least below 105 . However,when the attention is focused at the characteristics related tofluctuations, such as acoustics, high Re simulation is needed.
Our work gives an insight into the understanding of the real-

Figure 6. The vortical structures of the swimming cow-nosed
ray at Re = 1.48×104 , visualized by the isosurface of Q = 300
and colored by the streamwise velocity.

412th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022

Figure 7. The surface pressure distributions at (a),(d) t/T = 0.17, (b),(e) t/T = 0.33, (c),( f ) t/T = 0.5, at (a − c) Re = 1.48 × 104
and (d − f ) Re = 1.48 × 105 .

Figure 8. The spanwise vorticity contours and pressure distribution along the chord at the plane of z/SL = 0.31, 0.62, at (a) Re =1.48 × 104 and (b) Re = 1.48 × 105 .

Figure 9. The surface shear distributions at (a),(d) t/T = 0.17, (b),(e) t/T = 0.33, (c),( f ) t/T = 0.5, at (a − c) Re = 1.48 × 104 and(d − f ) Re = 1.48 × 105 .

512th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan, July 19–22, 2022
Ma, M., Huang, W. X. & Xu, C. X. 2019 A dynamic wallmodel for large eddy simulation of turbulent flow over complex/moving boundaries based on the immersed boundary
method. Physics of Fluids 31 (11), 115101.
Ma, M., Huang, W. X., Xu, C. X. & Cui, G. X. 2021 A
hybrid immersed boundary/wall-model approach for largeeddy simulation of high-reynolds-number turbulent flows.
International Journal of Heat and Fluid Flow 88, 108769.
Mao, Qian, Zhao, Jiazhen, Liu, Yingzheng & Sung, Hyung Jin
2021 Hydrodynamic benefits of protruding eyes and mouth
in a self-propelled flexible stingray. Physics of Fluids 33 (8),081915.
Menzer, Alec, Gong, Yuchen, Fish, Frank E & Dong, Haibo
2022 Bio-inspired propulsion: Towards understanding therole of pectoral fin kinematics in manta-like swimming.
Biomimetics 7 (2), 45.
Rosenberger, Lisa J 2001 Pectoral fin locomotion in batoid
fishes: undulation versus oscillation. Journal of Experimental Biology 204 (2), 379–394.
Triantafyllou, MS, Triantafyllou, GS & Gopalkrishnan, R
1991 Wake mechanics for thrust generation in oscillating
foils. Physics of Fluids A: Fluid Dynamics 3 (12), 2835–
2837.
Zhang, D., Huang, Q., Pan, G., Yang, L. & Huang, W. X.
2022a Vortex dynamics and hydrodynamic performance enhancement mechanism in batoid fish oscillatory swimming.
Journal of Fluid Mechanics 930, A28.
Zhang, Dong, Pan, Guang, Chao, Liming & Yan, Guoxin
2018a Mechanisms influencing the efficiency of aquatic locomotion. Modern Physics Letters B 32 (25), 1850299.
Zhang, Dong, Pan, Guang, Chao, Liming & Zhang, Ya 2018b
Effects of reynolds number and thickness on an undulatory
self-propelled foil. Physics of Fluids 30 (7), 071902.
Zhang, Dong, Zhang, Jun-Duo & Huang, Wei-Xi 2022b Physical models and vortex dynamics of swimming and flying:a review. Acta Mechanica pp. 1–40.
Zhao, Jiazhen, Mao, Qian, Pan, Guang, Huang, QiaoGao &Sung, Hyung Jin 2021 Hydrodynamic benefit of cephalicfins in a self-propelled flexible manta ray. Physics of Fluids
33 (8), 081906.

istic turbulent flow of aquatic swimming and may inspire thedesign of biomimetic propulsors in the real environment.

ACKNOWLEDGEMENTS
The authors acknowledge funding support from the National Natural Science Foundation of China (Grant Numbers
12102227 and 11772172) and China Postdoctoral Science
Foundation (Grant Number 2022M711754).

REFERENCES
Adrian, Ronald J 2007 Hairpin vortex organization in wall turbulence. Physics of fluids 19 (4), 041301.
Bottom Ii, RG, Borazjani, I, Blevins, EL & Lauder, GV 2016
Hydrodynamics of swimming in stingrays: numerical simulations and the role of the leading-edge vortex. Journal of
Fluid Mechanics 788, 407–443.
Bozkurttas, M, Mittal, R, Dong, H, Lauder, GV & Madden,
P 2009 Low-dimensional models and performance scalingof a highly deformable fish pectoral fin. Journal of Fluid
Mechanics 631, 311–342.
Chen, Di, Kolomenskiy, Dmitry, Nakata, Toshiyuki & Liu,
Hao 2017 Forewings match the formation of leading-edge
vortices and dominate aerodynamic force production in revolving insect wings. Bioinspiration & biomimetics 13 (1),016009.
Dewey, Peter A, Carriou, Antoine & Smits, Alexander J 2012
On the relationship between efficiency and wake structureof a batoid-inspired oscillating fin. Journal of fluid mechanics 691, 245–266.
Fish, F. E., Schreiber, C. M., Moored, K. W., Liu, G., Dong,
H. & Bart-Smith, H. 2016 Hydrodynamic performance ofaquatic flapping: efficiency of underwater flight in the
manta. Aerospace 3 (3), 20.
Huang, Q., Zhang, D. & Pan, G. 2020 Computationalmodel construction and analysis of the hydrodynamics ofa rhinoptera javanica. IEEE Access 8, 30410–30420.
Huang, W. X., Chang, C. B. & Sung, H. J. 2011 An improved
penalty immersed boundary method for fluid–flexible bodyinteraction. Journal of Computational Physics 230 (12),5061–5079.

6