Changes in rays’ swimming stability due to the phase difference between left and right pectoral fin movements
Analytical targetsTwo species of undulation motion rays with different pectoral fin shapes bred in KAIYUKAN were analyzed: sharpnose stingray Dasyatis acutirostra and pitted stingray Dasyatis matsubarai (Fig. 1a,b). Blender 2.7925 was used to construct stingray models from pictures26,27 as accurately as possible; Blender is a free and open-source 3D creation suite used to make realistic characters for movies, etc. Detailed information on how to construct models using Blender is provided in our previous paper28. To focus on the effects of pectoral fin movements, we did not consider the body’s shape as in the previous studies12,29. The height and disk width (WD) of all models were set to 0.01 m and 0.44 m, respectively, considering the previous studies30,31. The disk length of each model was determined from WD, referring to the aspect ratio of the rays’ photographs26,27; the disk length (LD) of D. acutirostra and D. matsubarai are 0.348 m and 0.344 m, respectively.Figure 1Analytical targets and description of motion. (a) Analytical model of D. acutirostra. (b) Analytical model of D. matsubarai. (c) Description of motion, (d) the relationship between any two points on the surface before and after the deformation.Full size imageMotionThe motion was given to satisfy the following equations:$$z = left{ {begin{array}{*{20}l} {A;sin left( {omega left( {t – kTleft( {frac{{angl{text{e}}left( {x_{i} ,y_{i} } right) – 10^{{text{o}}} }}{{All;angl{text{e}}}} – theta } right)} right)h_{1} h_{2} } right.} hfill & {left( {10^{{text{o}}} le angl{text{e}}left( {x_{i} ,y_{i} } right) le 170^{{text{o}}} } right)} hfill \ {A;sin left( {omega left( {t – kTleft( {frac{{350^{{text{o}}} – left( {angl{text{e}}left( {x_{i} ,y_{i} } right) – 10^{{text{o}}} } right)}}{{All;angl{text{e}}}}} right)} right)h_{1} h_{2} } right.} hfill & {left( {190^{{text{o}}} le angl{text{e}}left( {x_{i} ,y_{i} } right) le 350^{{text{o}}} } right)} hfill \ end{array} } right.$$
(1)
$$begin{array}{c}{h}_{1}=a{r}_{i}^{3}+b{r}_{i}^{2}+c{r}_{i}end{array}$$
(2)
$$h_{2} = left{ {begin{array}{*{20}l} {dleft( {angleleft( {x_{i} ,y_{i} } right) – 10^{ circ } } right)^{2} + eleft( {angleleft( {x_{i} ,y_{i} } right) – 10^{ circ } } right)} hfill & {left( {10^{ circ } le angleleft( {x_{i} ,y_{i} } right) le 170^{ circ } } right)} hfill \ {dleft( {350^{ circ } – left( {angleleft( {x_{i} ,y_{i} } right) – 10^{ circ } } right)} right)^{2} + eleft( {350^{ circ } – left( {angleleft( {x_{i} ,y_{i} } right) – 10^{ circ } } right)} right)} hfill & {left( {190^{ circ } le angleleft( {x_{i} ,y_{i} } right) le 350^{ circ } } right)} hfill \ end{array} } right.$$
(3)
$$begin{array}{c}{left({r}_{i}-{r}_{i-1}right)}^{2}+{left({z}_{i}-{z}_{i-1}right)}^{2}={left({r}_{i}^{mathrm{^{prime}}}-{r}_{i-1}^{mathrm{^{prime}}}right)}^{2}+{left({z}_{i}^{mathrm{^{prime}}}-{z}_{i-1}^{mathrm{^{prime}}}right)}^{2}end{array}$$
(4)
$$begin{array}{c}angleleft({x}_{i},{y}_{i}right)= angleleft({x}_{i}^{^{prime}},{y}_{i}^{^{prime}}right).end{array}$$
(5)
Equation (1) represents the amount of movement of the model surface in the z-axis direction, where (A) is the amplitude of the pectoral fin tip, (omega) is the angular velocity, (t) is time, (k) is the wavenumber, (T) is the period, angle(({x}_{i},{y}_{i})) is the angle made by the line connecting the center of rotation and any point (({x}_{i},{y}_{i})) on the model surface with the x-axis, Allangle is the range where the motion is given (160°), and (theta) is the phase difference between the movements of the right and left pectoral fins (Fig. 1c). ({h}_{1}) is the weighting from the center to the radial direction: it is necessary to set the amplitude at the ray’s center to zero and increase the amplitude toward the pectoral fin tip (a = 119.786, b = -7.957, c = 0.498). ({h}_{2}) is the weighting in the circumferential direction: it is necessary to increase the amplitude from the anterior to the tip of the pectoral fin and decrease the amplitude from the tip of the pectoral fin to the posterior (d = − 1.563 × 10–4, e = 0.025). Equation (4) is the condition in which the distance between two neighboring points in the same radial direction is equal before and after the movement (Fig. 1d). (r) is the distance between the center of rotation and any point (({x}_{i},{y}_{i})), defined as (sqrt{{x}_{i}^{2}+{y}_{i}^{2}}). Equation (5) defines angle(({x}_{i},{y}_{i})) as being constant before and after the move (Fig. 1c). The variables after the move are marked with ‘. Variables used in the analysis are A = 0.089 m, T = 0.499, k = 1.270, and ω = 12.599 rad/s. Videos of the created motion from the front and the side are shown in the “Supplement” (Supplement Movies 3, 4).Analytical conditionsAnalysis cases were conducted with eight conditions: two types of pectoral fin shape (Fig. 1a,b) and four types of phase difference (0 (T), 0.25 (T), 0.5 (T), and 0.75 (T)). These conditions were set for investigating the effects of phase differences between left and right pectoral fin movements on swimming and how these effects vary with pectoral fin shape.Numerical methodsA CFD simulation of the ray models in the water flow was performed using OPENFOAM, an open-source finite volume method CFD toolbox32, to calculate the forces acting on the rays in each axial direction and the moment around each axis. The governing equations were the continuity equation and the three-dimensional incompressible Reynolds-averaged Navier–Stokes equation, expressed by:$$begin{array}{*{20}c} {nabla cdot u = 0} \ end{array}$$
(6)
$$begin{array}{*{20}c} {frac{partial u}{{partial {text{t}}}} + nabla cdot left( {uu} right) = – nabla p + nabla cdot left( {vnabla u} right) + nabla cdot left[ {nu left{ {left( {nabla u} right)^{T} – frac{1}{3}nabla cdot uI} right}} right], } \ end{array}$$
(7)
where (u) is the velocity vector, t is the time, p is the static pressure divided by the reference density, (nu) is the kinematic viscosity, and I is the unit tensor. The Reynolds number was defined regarding the previous studies3 as:$$begin{array}{c}{R}_{e}=frac{U{L}_{D}}{nu },end{array}$$
(8)
where (U) (/ms) is the given flow speed3 (1.5 × LD/ms), ({L}_{D}) (m) is the length of the ray models, and (nu) is the kinematic viscosity of water at 20 °C (1.0 × 10–6 m2/s). The Reynolds number in this study is 1.8 × 105; considering this, we used the k–ω shear stress turbulence model33,34. The k–ω shear stress turbulence model is a type of Reynolds-averaged Navier–Stokes equation (RANS) turbulence model that is widely used to calculate for the fish swimming flow35,36,37. The overset grid method was used in this study; it is a generic implementation of overset meshes. For both static and dynamic cases, cell-to-cell mapping between multiple, disconnected mesh regions is employed to generate a composite domain38,39. This method permits complex mesh motions and interactions without the penalties associated with deforming meshes. The process is described in detail by Noack40. The calculation volume was 5.4 WD in length, 5.4 WD in height, and 5.4 WD in width (Fig. 2a,b). A hexahedral volume mesh was created using the snappyHexMesh of OPENFOAM. The fluid region was divided into two parts: the overset region and the background region (Fig. 2a,b). The overset region moves and transforms to match the motion of the ray and was made with fine meshes around the analysis target and coarse meshes in the outlying areas; a one-layer boundary layer mesh was created around the analysis target. The overset region shape is an ellipsoid (Fig. 2a,b). The minimum mesh volume is 7.3 × 10–10 (m3), and the maximum mesh volume is 2.6 × 10–2 (m3). The total number of meshes was 9.0 × 105. At the outlet boundary, the average static relative pressure was set to 0 Pa. The surfaces of the fish model were formed into non-slip surfaces.Figure 2Meshes for CFD simulation and differences in force between different meshes. (a) Meshes at the coronal plane of the whole fluid region. (b) Frontal cross-section of the fluid region at the green line in (a). The red region is the overset region. (c,d) Comparison of the instantaneous drag coefficient and the moment coefficient around the y-axis of D. matsubarai between the coarse, fine, and dense mesh.Full size imageThe drag coefficient ({C}_{D}left(tright)), the lateral force coefficient ({C}_{l}left(tright)), the lift coefficient ({C}_{L}left(tright)), the moment coefficient around the x-axis ({C}_{mx}left(tright)), the moment coefficient around the y-axis ({C}_{my}left(tright)) and the moment coefficient around the z-axis ({C}_{mz}left(tright)) were calculated as:$$begin{array}{c}{C}_{D}left(tright)=frac{Dleft(tright)}{frac{1}{2}rho {U}^{2}{L}_{D}{W}_{D}}end{array}$$
(9)
$$begin{array}{c}{C}_{l}left(tright)=frac{lleft(tright)}{frac{1}{2}rho {U}^{2}{L}_{D}{W}_{D}}end{array}$$
(10)
$$begin{array}{c}{C}_{L}left(tright)=frac{Lleft(tright)}{frac{1}{2}rho {U}^{2}{L}_{D}{W}_{D}}end{array}$$
(11)
$$begin{array}{c}{C}_{mx}left(tright)=frac{{M}_{psi }left(tright)}{frac{1}{2}rho {U}^{2}{L}_{D}^{2}{W}_{D}}end{array}$$
(12)
$$begin{array}{c}{C}_{my}left(tright)=frac{{M}_{phi }left(tright)}{frac{1}{2}rho {U}^{2}{L}_{D}^{2}{W}_{D}}end{array}$$
(13)
$$begin{array}{c}{c}_{mz}left(tright)=frac{{M}_{theta }left(tright)}{frac{1}{2}rho {U}^{2}{L}_{D}^{2}{W}_{D}},end{array}$$
(14)
where (Dleft(tright)) is the calculated drag, (lleft(tright)) is the calculated lateral force, (Lleft(tright)) is the calculated lift, ({M}_{psi }left(tright)) is the calculated moment around the x-axis, ({M}_{phi }left(tright)) is the calculated moment around the y-axis, ({M}_{theta }left(tright)) is the calculated moment around the z-axis, and (rho) (kg/m3) is the density of water at 20 °C (998 kg/m3). As shown in a previous study41. the propulsive efficiency (eta) is defined as the ratio of output power ({P}_{o}) to input power ({P}_{e}) which can be written as:$$begin{array}{c}{P}_{o}left(tright)=frac{1}{T}{int }_{0}^{T}Dleft(tright)Udtend{array}$$
(15)
$$begin{array}{c}{P}_{e}left(tright)=frac{1}{T}{int }_{0}^{T}left[Dleft(tright)dot{x}left(tright)+lleft(tright)dot{y}left(tright)+Lleft(tright)dot{z}left(tright)right]dtend{array}$$
(16)
$$begin{array}{c}eta =frac{{P}_{o}}{{P}_{e}}.end{array}$$
(17)
An in-house program calculated the forces acting on rays in each axial direction and the moment around each axis. The numerical method’s validity and reliability were verified by comparing previous experimental and numerical analytical studies of heaving and pitching on airfoil naca001341. A high degree of similarity to previous studies was confirmed; the mean difference in the propulsive efficiency from the previous study of analysis was 5%, and the difference from the previous study of the experiment was 9%. Detailed information such as mesh, length, and velocity, of this analysis method’s verification is provided in the “Supplement”.A grid sensitivity study was conducted using three meshes: coarse, fine, and dense. The coarse mesh has 8.1 × 105 elements, the fine mesh has 9.0 × 105 elements, and the dense mesh has 9.9 × 105 elements. The analysis was conducted using a condition with no phase difference of D. matsubarai. As shown in Fig. 2c,d, the drag coefficient and the moment coefficient around the y-axis are almost the same when the mesh is fine and when the mesh is dense. The mean drag and propulsive efficiency error of fine and coarse meshes are 2.7% and 3.5%, respectively. The fine mesh was used in all simulation cases considering accuracy. We used the same meshes for all cases. More
