HYDRODYNAMIC ANALYSIS OF BIONIC CHIMERICAL WING PLANFORMS INSPIRED BY MANTA RAY EIDONOMY

In this paper, inspired by the external morphology of a manta ray (Mobula alfredi), four chimerical wing planforms are designed to assess its gliding performance. The planforms possess an arbitrary combination of extra hydrodynamic features like tubercles at the leading edge (L.E.) and trailing edge (T.E.) inspired by humpback whale's flippers and flukes, respectively, as longitudinal ridges inspired by whale shark's economy. In addition, another planform is designed to investigate the possible effects of manta ray's injuries (geometric deficiency) generated by predator's attacks or boat strikes on its locomotion (gliding) performance. In this regard, turbulent flow physics involved in the problem is numerically simulated at different angles of attack (AoA). High Reynolds number, 10, corresponding to the swimming of a juvenile manta ray at an average speed equals one m/s. The results show that the manta ray-inspired planform with L.E. undulations exhibits a superior performance at high AoAs than its other counterpart variants. In addition, the results demonstrate that injuries on the manta ray's body can noticeably modify hydrodynamics and, as a result corresponding hydrodynamical forces and moments acting on the swimming animal in the gliding phase.


INTRODUCTION
The external morphology of aquatic creatures plays a crucial role in the hydrodynamics of locomotion [1]. As an example, the presence of tubercles on the leading edge of the humpback whale's flippers leads to the formation of streamwise vortices, which ultimately delays the separation onset and development [2] [3]. In general, for all marine creatures, their external geometry along with the presence of special geometrical features are majorly responsible for some favorable hydrodynamic characteristics in swimming like denticles on the sharkskin [4], ridges on the leatherback turtle's carapace [5], ridges on the whale shark's body [6] and ventral pleats on the belly of a humpback whale [7], to name a few.
The Manta ray, Mobula alfredi, is one of the most charismatic swimming animals in oceanic islands and reefs (Fig.1). These species belong to the Mobulidae family and can be typically found in Indo-west Pacific oceans with a confirmed range of distribution in Indonesia coasts and the red sea, extended to an expected range in the 'Persian Gulf' at the south of Iran [8]. Their mature size can reach 2.7-3.5 meters on average, depending on the male/female categories. Its life span is about 40 years, with a weight reaching 0.7 tons [8]. In contrast to fishes, including sharks, these filter-feeding swimmers have developed a disproportionately large brain to their body weight similar to mammals, which gives them a higher level of Indonesian Journal of Engineering and Science, Vol. 2, No. 3, 2021 Taheri https://doi.org/10.51630/ijes.v2i3. 25 12 capabilities in functionality and behavior [8]. In general, the swimming of a manta ray is performed via a combination of two modes [9]: flapping of its left and right pectoral fins [10] [11] and gliding mode provided by its large and extended planform or sometimes just as a pure gliding mode [12] [13]. As aforementioned, the external geometry of a manta ray, including wing planform shape and the shape of hydrofoil ribs, directly affects the swimming performance of the animal in both flapping and gliding modes. For instance, Luo et al. has successfully used an optimization technique to optimize the shape of the hydrofoil sections for a given manta ray-inspired planform in the pure gliding mode [9]. Here, we implement some biomimetic changes in the overall economy (planform shape) of a manta ray-inspired wing to grasp the underlying physics of the problem. In this regard, non-uniform bionic L.E. and T.E. undulations are applied to the manta-ray inspired wing planform to investigate hydrodynamic consequences towards an optimum design strategy. Furthermore, hydrodynamic effects of typical injuries on the manta ray's gliding performance are also assessed by considering an injured body shape of a Maldives best-known manta ray, so-called 'Babaganoush,' dated May 2019 [14]. In the following sections, details are presented.

CHIMERICAL MANTA RAY BODY MODELS
To construct the chimerical wings based on the mana ray's body planform in the present paper, a set of hydrofoils (i.e., eleven cross-sections of the Mobula alfredi body model) at selected planes perpendicular to the lateral axis are numerically generated by high resolution in MATLAB and imported into the SolidWorks CAD environment [17]. In this regard, an optimized hydrofoil section shape introduced by Luo et al. is adopted here [9]. Then, a total number of 2 guide-curves representing L.E. and T.E. of the manta ray planform (in the midplane defined by y=0) are fed to the SolidWorks CAD environment. Finally, the final 3D model variants of the manta ray planform models are constructed by applying a 'Lofting process,' as shown in Fig. 2 and Fig. 3. The latter process is performed by stitching successive cross-section curves of the geometry guided by the guide curves. As shown in Fig. 2 and Fig. 3, The wing planforms possess an arbitrary combination of some selected hydrodynamic features like tubercles at L.E. and T.E. inspired by humpback whale's flippers [3] and flukes, respectively, as well as longitudinal ridges inspired by whale shark's economy [6].
Considering aforementioned hydrodynamic specific features, a total number of 4 bioinspired wings without a dihedral angle is designed corresponding to planforms exhibiting: smooth LE-smooth TE (i.e., real Mobula alfredi planform), wavy LE-wavy TE, wavy LEsmooth TE, smooth LE-wavy TE, as summarized in Table 1 (illustrated in Fig. 3). With the aid 13 of these hydro-planforms, hydrodynamic effects of these prominent geometric features such as the evolution of vortical structures, flow separation evolutions, and post-stall behavior of these 'swimming wings' can be assessed, as discussed in the present paper. In general, geometric deficiencies can be present in the manta ray's external morphology (economy) due to injuries or bites on T.E. of their pectoral fins caused by predator's attacks (e.g. sharks/mammals), boat strikes while feeding or swimming close to the oceanic water surface. As explained in detail by Stevens et al., manta rays exhibit a noticeable capability to quickly heal these body wounds and regenerate large portions of the missing flesh [8]. Many examples of manta rays with these types of injuries have been observed by Stevens et al. over the years in the Maldives and Mozambique manta ray populations [8]. As mentioned earlier, one of the famous examples with a 'geometric deficiency' at T.E. is a reef manta ray (Mobula alfredi), namely 'Babaganoush.' The manta ray has been injured by a speedboat strike and has been observed for the first time in the Maldives by 'The Manta Trust' research team having large and deep injuries at T.E. on November 2018 [14]. However, the manta was incredibly healed in the next show-up in May 2019, as observed by' The Manta Trust.' In the present paper, the body planform of 'Babaganoush' in May 2019 is considered to assess hydrodynamic consequences of these kinds of 'geometric deficiencies' (Fig. 3.c).

NUMERICAL METHODOLOGY
To pursue the goals of the present study, a numerical computation campaign consisting of 90 individual flow simulations is planned, including simulations of the designed planforms (Table 1) and validation test cases (section 3.3). For manta ray simulations, the inflow velocity is set based on a prescribed Reynolds number of 10 6 . In general, manta rays are among relatively slow swimmers in oceans and reefs. For instance, an oceanic manta ray (Mobula birostris) has been recorded with a swimming speed of 0.46-2.51 m/s while turning [18] and about 0.25-0.47 m/s in the 'foraging phase' [19]. On the other hand, giant mantas (Mobula tarapacana) can dive with up to a six m/s descending speed to the depth of 2000 m with a low temperature [20] by keeping the brain warm enough compared to its surrounding tissues provided by their particular 'network of blood vessels [8]. It was also reported that a reef manta ray (Mobula alfredi) could reach a 432 m depth [21]. More recently, reef manta rays have reached a deeper depth of 672 m during nighttime [22].
On the other hand, water stream currents in oceans and reefs are another essential factor that should be considered. In general, current ocean speed is fastest near the ocean surface with about 2.5 m/s, and gulf current speed is about 1.8 m/s [23]. As an example, merging of the Pacific and Indian ocean waters generate one of the strongest reef currents on the planet around Indonesian Journal of Engineering and Science, Vol. 2, No. 3, 2021 Taheri https://doi.org/10.51630/ijes.v2i3. 25 15 the Indonesian islands with the speed of about four m/s, which makes it difficult to swim and hover (keeping a position unchanged in 3D space) for marine creatures like mantas and also human divers. As a result, oceanic and reef mantas can typically experience a wide range of AoA ( ) generated by a vector summation of their swimming speed and natural stream current speeds. It is worth mentioning that manta rays are also facing a broad range of AoAs in performing a vast range of maneuvers and feeding strategies such as 'cyclone' and 'somersault' feedings, to name a few.
Here to perform the planned simulations at different AoAs ( ), ranging from 0 to 40 degrees, including deep-stall region, the model is placed at a fixed position in space, and the relevant AoAs are implemented via setting of the freestream blowing angle for the streamwise x-axis (Fig. 4). For the simulations with a given AoA ( ), the inflow velocity field is applied by the following formula: The coordinate system adopted for the upcoming manta ray gliding flow simulations

Computational Domain and Grid Generation
Computational grids are generated in the SolidWorks meshing tool environment [17] [24]. It takes advantage of an adaptive local grid generation on a global Cartesian mesh to capture fine details of the bionic wing's geometries and its corresponding boundary layers [24]. After applying grid convergence tests, a well-converged grid with about the ultimate 2 million elements is utilized for all upcoming flow simulations over all variants of the manta ray geometry. As an example, Fig. 5 shows the grid generated around the bionic planform with wavy LE-wavy TE. It exhibits a well-resolved body-fitted mesh clustering to the wall, capable of capturing the geometry's delicate features.
As one can see in the figure, to minimize the boundary effects, the computational domain is extended about three times of the body length in all directions, except in the downstream, which is extended four times the body length to capture wakes and tip vortices more clearly. Fig. 6 depicts top-views of the computational grids generated around all variants of the bionic planforms. As shown in Fig. 5 and Fig. 6, applying three levels of clustering and an adaptive meshing technique, a smooth transition between a relatively coarse mesh in the outer-flow region and the fine mesh in the near-body zone is achieved.

Turbulent Flow Simulations and Settings
Incompressible Navier-Stokes equations govern Single-phase turbulent flows over juvenile and mature manta rays. Here, the equations are solved 6 10 Re  using Reynold Averaged Navier-Stokes (RANS) technique [24]. In this regard, turbulence is treated with Lam-Bremhorst low-Reynolds (LB-LRN, hereafter) by SolidWorks Flow Simulation (SFS) solver [24] [25]. The name of 'low-Reynolds number' shows the capability of the method to treat both low-speed boundary layer (B.L.) and high-speed regions in a single framework without wall-modeling. This is achieved in the SFS solver by applying at least ten nodes in the B.L. and considering wall distances in the turbulence model calculations [24]. In general, Navier-Stokes equations can be expressed in the tensor notation, as below [24]: where i S and *  is defined as: t  Furthermore, k is turbulent eddy viscosity and kinetic energy, respectively, while ij  denotes the Kronecker delta. In LB-LRN model, k and  are coupled as [24,25]: where the parameters are adjusted, as below: To calculate t  , the following formula is adopted [26]: In the above expressions, where the parameters are calculated as: 2 ,.
where y is the wall distance. Finally turbulence time and length scales can be computed as the following [26]: In SFS, fluid flow equations are solved using the finite-volume method by applying operator-splitting, conjugate gradient, and multigrid and SIMPLE computing techniques [17] [24]. Here for turbulent flow simulations over the bionic planforms, inlet velocity is adjusted with equation (1), and other boundaries are set as 'outflows.' Inflow turbulence intensity and length are set as 0.1% and 1.210 -4 m, respectively. Forces and moments, including: {Fx, Fy, Fz, Mx, My, and Mz} acting on the bionic planforms, are continuously monitored to ensure the satisfaction of all convergence criteria. SFS solver automatically applies the measures based on the selected global and local goal functions (e.g., defined as velocity, pressure, forces, and moments) to achieve the lowest residual errors [17].

Validation Test Cases
To investigate the performance of the proposed numerical strategy, two preliminary cases are considered first before jumping to the turbulent flow simulations over the manta rayinspired planforms, including a low-aspect-ratio wing ( 1 AR = ) and an initial manta planform.

Low-Aspect Ratio Wings
Low-aspect ratio wings possess different characteristics compared to the high-aspect-ratio ones with long spans. By decreasing the aspect ratio (defined as 2 / AR b c = , where b is the wingspan and c is the averaged chord of the wing), trailing/tip vortices become more influential over a large portion of the wing and modify the aero/hydrodynamic performance of the wing. As a result, by decreasing AR , the lift slope decreases [27]. Here, turbulent flows over a low-aspect-ratio wing ( 1 AR = ) with NACA 0012 airfoil sections are simulated at a high Reynolds number 5 Re 1.5 10  (Fig. 7). In this regard, a well-converged grid with about 1 million elements is utilized. The computational domain is extended to 4 and 10 times the wingspan in x − and x + directions, respectively. In addition, the domain is extended 3 times of the wing span in the lateral directions to minimize the boundary effects. Fig. 7. Formation of the tip vortices on a low-aspect-ratio wing with NACA 0012 airfoil section, visualized by tracer particle dynamics and vorticity iso-contours at 12  = Fig. 7 shows the formation of the wingtip vortices over the low-aspect-ratio wing at AoA equals to12 . In this regard, a tracer particle study has been performed. To do so, the tracer particles, about 300 water spherical particles with a diameter of 0.0001 m, are continually released from the wing surface and convected downstream by the background flow field

19
(having the same velocity as the local flow field). For the calculations, ideal reflection has been applied for fluid particle-solid interactions, e.g., in the separation zones on the wing. As shown in Fig. 7, tip vortices extend far downstream with a spiral motion visualized by tracer particles colored by axial velocity. Intersections of the counter-rotating tip vortices are also visible on a plane positioned /4 xb = downstream of the wing visualized by iso-contours of vorticity, 100 x  = 1/s. Fig. 8 depicts the lift coefficient curve versus AoA for the adopted wing compared to the experimental curve measured by Chen et al. [28] and the theoretical curve by Lowry and Polhamus [29], which exhibits a perfect agreement.

Initial Manta-Ray Planform (Comparative Test Case)
In this section, turbulent flows over a manta ray-inspired 'test case' planform with an optimized airfoil cross-section proposed by Luo et al. [9] are simulated. In this regard, the test case geometry is constructed by the lofting process via providing L.E. and T.E. guide curves and a total number of eleven airfoil cross-sections between left and right wingtips similar to the strategy aforementioned in section (2) for the designed bionic chimerical planforms. Here, to perform numerical simulations of the turbulent flows, a well-converged grid with a maximum of 1.5 million elements is utilized. Fig. 9 (top) shows grid generation around the test case geometry. As one can see in the figure, all detail features of the geometry are captured by the generated body-fitted grid via three levels of refinements.  20 Similar to the adopted setting explained in sub-section (3.3.1), a tracer particle study was performed to capture tip vortex topology. Fig. 9 (Bottom) shows the topology of the counterrotating tip/trailing vortices 10  = . As one can see in the figure, the platform generates a relatively large vorticity region ( 1 x  = 1/s) downstream, induced by tip/trailing vortices 10  = . Fig. 10 depicts lift coefficient versus AoA obtained by the present study, which exhibits a good agreement over a full range of AoA, 0 10   , with the Luo et al. results [9].

RESULTS AND DISCUSSIONS
In this section, complete results of the turbulent flow simulations for different variants of the designed chimerical manta ray-inspired planforms (Table 1

Chimerical Planforms (Comparative Study)
As discussed in detail in section (2), in the present study, the main layout of the platforms has been designed based on a reef manta ray (Mobula alfredi) eidonomy. The final geometries are constructed by adding some extra topological features like L.E. and T.E. undulations inspired by humpback whales [2] [3] and longitudinal ridges inspired by whale sharks [6]. In this regard, geometrical features of the related species, for example, T.E. undulations of the humpback whale's fluke, are digitalized, normalized, and finally rescaled and implemented at T.E. of the designed manta ray-inspired planform (Fig. 2). A similar procedure is also applied for L.E. undulations inspired by humpback whale's flippers (Fig. 2). The effects of AoA variations on the evolution of tip/trailing vortices for these planforms are shown in Fig. 11 via tracer particle dynamics. As one can see in the figure, by increasing AoA, tip vortices are getting more distorted and tilted upward. Furthermore, Iso-contours of vorticity defined by the axial vorticity, 1 x  = 1/s, at a plane positioned at /2 xb = also show that by increasing AoA, cross-sections of tip/trailing vortices are getting larger and intensified.
For example, Fig. 12 shows the formation and evolution of the vortical structures by increasing AoA for the bio-inspired wing planforms with smooth LE-smooth TE (I) compared to wavy LE-wavy TE (II). The vortical structures are captured here by λ2 -criterion (λ2=-1). As one can see in the figure, a pair of wall-bounded horseshoe vortices with a relatively large scale (comparable to the planform chord) is generated over the wing at α = 20°. By increasing AoA, the vortex pairs become more significant and finally merge at α = 36° and produce a large dominant vortical structure. As it is also evident in the figure, the presence of L.E. and T.E. undulations modify formation, pattern, and evolution of the minor and significant vortical structures at all AoAs. Taheri https://doi.org/10.51630/ijes.v2i3. 25 22 Fig. 13 shows mean separation zones at different AoAs for all bio-inspired planforms introduced in Table 1. As shown in Fig. 12, L.E. and T.E. undulations with their particular peak and trough topology (here based on the humpback whale's flipper and fluke, respectively) considerably affect fluid dynamic characteristics over the wing in the gliding phase. As a result, the topology of separation zones is modified, as shown in Fig. 13. It is worth mentioning that an intensive back-flow is induced and embraced by the dominant horseshoe vortex pairs, which leads to the formation of two significant disconnects (not merged) separation zones on the left and right sides of the bio-inspired swimming wings. For example, in the planform with smooth LE-smooth TE at α = 24°, two disconnect separation zones are generated at the centers of the horseshoe vortex pairs, visible in Fig. 12 Fig. 13 23 By comparing the pattern and topology of the separation zones generated on the bioinspired planforms at all AoAs in Fig.13, one can conclude that a minor separation zone is generated in the planform with wavy LE-smooth TE Fig. 14

An Injured Manta Ray (Babaganoush on May 2019)
In section (2), Babaganoush, a well-known injured reef manta ray, was introduced. In this section, hydrodynamic consequences of its injury are assessed by CFD simulations. Fig. 15-a showed the injured body of Babaganoush in May 2019. As one can see in the figure, the aft Taheri https://doi.org/10.51630/ijes.v2i3. 25 24 part of the body near the centerline was lost in the boat strike accident in 2018. In the present calculations, a 3D model of the body is constructed by digitalizing the aft body pattern and exporting it to the SolidWorks CAD environment. Fig. 15-b shows an isometric view of the 3D model adopted for the upcoming simulations. Fig. 15-c also shows mesh generation around the body. Fine details of the 'geometric deficiency' are captured by appropriate three levels of refinements as before, which optimizes the accuracy of the computations and computational cost. Fig. 16 depicts streamlines over the body of Babaganoush at α = 12°. As one can see in the figure, a recirculation zone exhibiting chaotic dynamics is generated right after the injured part of the body, similar to the objects with blunt aft-body experiencing dominant pressure drag, in addition to the separation zones forming on the wing planform at this AoA. This feature is translated to higher drag coefficients almost at all AoAs than the healthy Babaganoush (base), as shown in Fig. 17. Lift coefficient variations in AoA for the injured and healthy Babaganoush are also depicted in Fig. 17. As shown in Fig. 17, lift coefficients of the injured body are modified and exhibit lower/higher values over specific intervals of AoA. Fig. 18 shows intersections of trailing/tip vortices on a plane positioned at x/b=2 visualized by iso-contours of vorticity, ωx=±1 1/s, at different AoAs from the back-side of the manta. At α = 24°, dissymmetry of the counterrotating vortices is more pronounced. Increasing AoA cross-sections of tip and trailing vortices merger are getting larger and intensified (Fig. 18). The final important topic is the nonequilibrium state, corresponding to an unbalance generation of forces and moments by the swimming animal. An injured manta, such as Babaganoush, experiences unbalance forces and moments exerted on the body, induced by the 'geometric deficiency.' This makes it more difficult for the animal to maintain a so-called 'trimmed swimming,' similar to the concept of a 'trimmed flight' for an airplane. For example, Fig. 19 shows the absolute values of rolling moments versus AoA exerted on the injured Babaganoush body compared to a healthy one. As one can see, the injured animal experiences higher levels of rolling moments at all AoAs. Therefore, the injured manta should perform asymmetric flapping and hold asymmetric curvatures to the left and right pectoral fins, e.g., dihedral angle, in the gliding and maneuvering phases to maintain the 'trimmed' swimming state.

CONCLUSION
The present paper assessed the effects of L.E. and T.E. geometric undulations on a manta ray-inspired planform. In this regard, four chimerical wings without applying dihedral angles were designed with optional combinations of tubercles at L.E. and T.E. extracted from humpback whale's flipper and fluke, respectively. The results showed that the designed wing planform with wavy LE-smooth TE exhibits a superior separation control and hydrodynamic performance (lift generation) in the gliding phase, especially at the post-stall region. The results also confirmed that the relative position of peaks and troughs in the L.E. and T.E. curves should be designed and optimized in a single framework (not individually) to achieve maximum performance. It is worth mentioning that the main benefits of implementing undulations at T.E. are achieved in the 'flapping' mode rather than the 'gliding' mode, as observed in the oscillatory motion of humpback whale's flukes in the vertical direction (up/down oscillatory motions). Numerical turbulent flow simulations of an injured reef manta ray, Babaganoush, showed that the 'geometric deficiency' leads to a higher drag coefficient and brings a higher level of Indonesian Journal of Engineering and Science, Vol. 2, No. 3, 2021 Taheri https://doi.org/10.51630/ijes.v2i3. 25 27 disequilibrium of forces and moments, making it more difficult for the animal to maintain the so-called 'trimmed swimming.'

ACKNOWLEDGEMENT
The author would like to sincerely acknowledge every effort by institutes, organizations, and individuals to protect 'Manta Rays' worldwide.