In the realm of underwater exploration, the development of autonomous systems capable of navigating complex and hazardous environments is paramount. My research focuses on advancing the field of bionic robots, specifically those inspired by marine life, to achieve efficient and agile propulsion. The bionic robot, mimicking creatures like the manta ray, offers significant advantages due to its high maneuverability and energy efficiency. This study delves into the unsteady hydrodynamic characteristics of flapping fins, a key propulsion mechanism for such bionic robots. By understanding the fluid-structure interactions during flapping, I aim to optimize design parameters for enhanced performance, ultimately contributing to the next generation of intelligent underwater bionic robots.
The inspiration for this work stems from the natural propulsion of manta rays, which utilize pectoral fin flapping for locomotion. This motion, characterized by combined heaving and pitching, falls under the median fin/paired fin (MPF) mode, which is simpler to implement in mechanical bionic robots compared to body/caudal fin (BCF) modes. To simplify the analysis, I approximate the fin cross-section using a two-dimensional Clark Y airfoil, a common profile in aerodynamic and hydrodynamic studies. The core objective is to numerically simulate the unsteady flow around this airfoil undergoing prescribed flapping motions and extract insights into how key kinematic parameters influence lift and thrust generation. The findings are directly applicable to the design and control of flapping-fin bionic robots.
My investigation employs computational fluid dynamics (CFD) as the primary tool. The success of such simulations heavily relies on robust mesh management for moving boundaries. I adopt a dynamic hybrid mesh strategy to discretize the computational domain. The region immediately surrounding the airfoil is discretized with a high-quality, full-quadrilateral boundary layer mesh generated via a hyperbolic normal extrusion method, ensuring accurate resolution of viscous effects. The first layer height is set to $10^{-5}c$, where $c$ is the airfoil chord length. The remainder of the fluid domain is filled with unstructured triangular cells, providing flexibility for mesh deformation. The transient flow solver is based on the Reynolds-Averaged Navier-Stokes (RANS) equations, with turbulence closure provided by the Spalart-Allmaras (SA) one-equation model. This model is chosen for its balance between accuracy and computational economy for external flows. The governing equation for the modified turbulent viscosity, $\tilde{\nu}$, is:
$$
\frac{\partial}{\partial t}(\rho \tilde{\nu}) + \frac{\partial}{\partial x_i}(\rho \tilde{\nu} u_i) = G_{\nu} + \frac{1}{\sigma_{\tilde{\nu}}} \left[ \frac{\partial}{\partial x_j} \left( (\mu + \rho \tilde{\nu}) \frac{\partial \tilde{\nu}}{\partial x_j} \right) + C_{b2} \rho \left( \frac{\partial \tilde{\nu}}{\partial x_j} \right)^2 \right] – Y_{\nu} + S_{\tilde{\nu}}
$$
where $G_{\nu}$ is the production term, $Y_{\nu}$ is the destruction term near walls, $\sigma_{\tilde{\nu}}$ and $C_{b2}$ are constants, $\mu$ is the molecular dynamic viscosity, and $S_{\tilde{\nu}}$ is a user-defined source term set to zero. The pressure-velocity coupling is handled by the “Coupled” scheme, and both spatial and temporal discretizations are second-order accurate. A time step corresponding to 1/400 of the flapping period is used, with 50 iterations per time step to ensure convergence.
A critical aspect of simulating flapping motion is the dynamic mesh update. Traditional spring-based or diffusion-based methods can struggle with large deformations and may produce negative-volume cells. To overcome this, I implement an efficient node-based method using Delaunay graph mapping. This approach treats the mesh motion as a scattered data interpolation problem. First, a background Delaunay triangulation is constructed connecting the nodes on the moving airfoil surface to fixed nodes on the far-field boundary. The position of any interior fluid node $P$ within a Delaunay triangle $ABC$ is uniquely defined by its area coordinates $(e_1, e_2, e_3)$:
$$
e_i = \frac{S_i}{S}, \quad i=1,2,3
$$
$$
S_1 = \frac{1}{2} \begin{vmatrix} x_P & y_P & 1 \\ x_B & y_B & 1 \\ x_C & y_C & 1 \end{vmatrix}, \quad S_2 = \frac{1}{2} \begin{vmatrix} x_A & y_A & 1 \\ x_P & y_P & 1 \\ x_C & y_C & 1 \end{vmatrix}, \quad S_3 = \frac{1}{2} \begin{vmatrix} x_A & y_A & 1 \\ x_B & y_B & 1 \\ x_P & y_P & 1 \end{vmatrix}, \quad S = \frac{1}{2} \begin{vmatrix} x_A & y_A & 1 \\ x_B & y_B & 1 \\ x_C & y_C & 1 \end{vmatrix}
$$
with $\sum_i e_i = 1$. When the surface nodes move by a displacement $\delta_j$, the displacement of the interior node $\delta_{P_j}$ is computed by interpolating the displacements of the triangle vertices using the same area coordinates: $\delta_{P_j} = \sum_i e_i \delta_{j}$. The topology of the background Delaunay graph and the pointer linking each fluid node to its host triangle remain unchanged during motion, making the update process highly efficient and robust for the moderate deformations encountered in this flapping study.
Before proceeding with the parametric study, I established confidence in the numerical setup through a grid independence test and validation against published data. Five mesh systems with varying refinement levels were constructed for a stationary Clark Y airfoil at $0^\circ$ angle of attack in a uniform flow of $U_{\infty}=10$ m/s (water at $20^\circ$C, chord $c=0.3$ m). The computed lift and drag coefficients converged with mesh refinement. The selected medium mesh offers an optimal balance, as summarized below:
| Mesh ID | Surface Nodes | Total Nodes | Total Cells | $C_d$ | $C_l$ |
|---|---|---|---|---|---|
| 1 | 61 | 17143 | 19292 | 0.0102 | 0.0011 |
| 2 | 121 | 31542 | 39450 | 0.0098 | 0.0008 |
| 3 | 141 | 37111 | 47708 | 0.0097 | 0.0007 |
| 4 | 161 | 43094 | 56794 | 0.0097 | 0.0007 |
| 5 | 181 | 49426 | 66578 | 0.0097 | 0.0007 |
For dynamic validation, I simulated the combined heaving and pitching motion of a NACA0012 airfoil as studied in prior literature. The motion is defined by $y(t) = y_0 \sin(kt)$ and $\theta(t) = \theta_0 \sin(kt + \phi)$, with $y_0=0.25c$, $\theta_0=20^\circ$, $k=1.0$, and phase shift $\phi$ varied. The Reynolds number was $40,000$. My computed cycle-averaged thrust coefficient ($\overline{C_T} = -\overline{C_d}$) showed good agreement with reference data, with average deviations within $\pm 8\%$, confirming the reliability of my methodology for analyzing flapping propulsion relevant to bionic robots.
To model the flapping cycle of a bionic robot fin, I define a motion that combines heave (vertical displacement) and pitch (angular rotation about the leading edge). A full period $T$ is divided into two parts: Part 1 (lower stroke) and Part 2 (upper stroke). The key parameter $p$, termed the upper zone ratio, defines the amplitude of the upper stroke relative to the total peak-to-peak amplitude $A$. The heave motion $y(t)$ is given by:
$$
y(t) =
\begin{cases}
-(1-p)A \sin\left( \frac{2\pi}{T_1} t \right), & t \in \left[0, \frac{T_1}{2}\right] \\
pA \sin\left( \frac{2\pi}{T_2} \left( t – \frac{T_1}{2} \right) \right), & t \in \left[ \frac{T_1}{2}, \frac{T_1+T_2}{2} \right]
\end{cases}
$$
The pitch motion $\alpha(t)$, defined as the angle of attack relative to the horizontal, is synchronized as:
$$
\alpha(t) =
\begin{cases}
-\alpha_{\text{max}} \cos\left( \frac{2\pi}{T_1} t \right), & t \in \left[0, \frac{T_1}{2}\right] \\
\alpha_{\text{max}} \cos\left( \frac{2\pi}{T_2} \left( t – \frac{T_1}{2} \right) \right), & t \in \left[ \frac{T_1}{2}, \frac{T_1+T_2}{2} \right]
\end{cases}
$$
Here, $\alpha_{\text{max}}$ is the maximum pitch amplitude. Negative $\alpha$ indicates nose-up rotation. To ensure kinematic smoothness, the velocities at the transition point are matched: $y'(0) = y’\left(\frac{T_1+T_2}{2}\right)$ and $y’\left(\frac{T_1}{2}\right)^+ = y’\left(\frac{T_1}{2}\right)^-$. This leads to the relation $\frac{1-p}{p} = \frac{T_1}{T_2}$. The flapping frequency is $f = 1/T = 2/(T_1+T_2)$. In this study, I set $A=0.5$ m and $f=1$ Hz. The parameter $p$ effectively controls the symmetry of the stroke; $p=0.5$ corresponds to a standard sinusoidal heave. The performance is evaluated via instantaneous and cycle-averaged force coefficients. The lift, drag, and moment coefficients are defined as:
$$
C_l = \frac{F_l}{\frac{1}{2} \rho U_{\infty}^2 S}, \quad C_d = \frac{F_d}{\frac{1}{2} \rho U_{\infty}^2 S}, \quad C_m = \frac{M}{\frac{1}{2} \rho U_{\infty}^2 S c}
$$
where $F_l$ and $F_d$ are lift and drag forces (drag positive along freestream direction), $M$ is the moment about the leading edge, $\rho$ is fluid density, $U_{\infty}$ is the forward speed (set to 2 m/s to simulate a cruising bionic robot), and $S$ is the planform area per unit span (chord $c$). A negative $C_d$ indicates net thrust production, which is the desired state for a propulsive bionic robot. The cycle-averaged values are denoted as $\overline{C_l}$ and $\overline{C_d}$. I also calculate the propulsive efficiency $\eta$ for the optimal case, defined as the ratio of useful thrust power to input mechanical power:
$$
\eta = \frac{-\overline{C_d} U_{\infty}}{\overline{C_p}}
$$
where the mean input power coefficient is:
$$
\overline{C_p} = -\frac{1}{T} \int_0^T \left[ C_l(t) \dot{y}(t) + C_m(t) \dot{\alpha}(t) \right] dt
$$
My parametric study systematically varies two key parameters: the upper zone ratio $p$ (from 0.1 to 0.9 in steps of 0.1) and the maximum pitch angle $\alpha_{\text{max}}$ (from $-25^\circ$ to $30^\circ$ in steps of $5^\circ$). This comprehensive exploration aims to guide the kinematic design of efficient bionic robot fins.
First, I analyze the effect of the upper zone ratio $p$ by setting $\alpha_{\text{max}} = 0^\circ$ (pure heaving motion). The instantaneous force coefficients over one cycle reveal distinct patterns. The drag coefficient $C_d$ shows two peaks (maximum drag) and two troughs (maximum thrust) per cycle. The peaks occur near the stroke reversals (highest heave velocity), while the timing of the thrust troughs shifts with $p$. The lift coefficient $C_l$ exhibits one major peak and one trough per cycle. A crucial observation is that the phase at which maximum thrust (most negative $C_d$) occurs also corresponds to near-maximum lift generation. I define an equivalent angle of attack $\alpha_e$ at this instant, combining the geometric pitch and the induced angle from heave velocity. For all $p$ values, this optimal state occurs at $\alpha_e \approx 8.6^\circ$, indicating a consistent aerodynamic condition for peak performance in this bionic robot flapping system.
| $p$ value | Time of Max Thrust, $t_p$ (s) | Equivalent AoA $\alpha_e$ (degrees) | $\overline{C_d}$ | $\overline{C_l}$ |
|---|---|---|---|---|
| 0.1 | 0.130 | 8.3 | -0.075 | 0.142 |
| 0.2 | 0.095 | 8.4 | -0.088 | 0.135 |
| 0.3 | 0.070 | 8.5 | -0.096 | 0.128 |
| 0.4 | 0.048 | 8.6 | -0.101 | 0.118 |
| 0.5 | 0.035 | 8.7 | -0.104 | 0.105 |
| 0.6 | 0.025 | 8.7 | -0.101 | 0.092 |
| 0.7 | 0.018 | 8.8 | -0.096 | 0.078 |
| 0.8 | 0.012 | 8.8 | -0.088 | 0.063 |
| 0.9 | 0.008 | 8.6 | -0.075 | 0.047 |
The cycle-averaged coefficients reveal the most significant trend. The mean drag coefficient $\overline{C_d}$ is negative for all $p$, confirming net thrust production. Its variation is symmetric about $p=0.5$, where it reaches a minimum (most negative) value of approximately $-0.104$. This signifies that the standard sinusoidal heave ($p=0.5$) generates the maximum average thrust. In contrast, the mean lift coefficient $\overline{C_l}$ decreases almost linearly with increasing $p$. For a bionic robot requiring positive buoyancy or depth control, lift generation might also be important, but the primary goal for propulsion is often maximizing thrust. Therefore, $p=0.5$ emerges as the optimal heave symmetry for the bionic robot’s flapping fin.
Next, I investigate the influence of the maximum pitch angle $\alpha_{\text{max}}$, fixing the heave symmetry at the optimal $p=0.5$. The introduction of pitching dramatically alters the force histories. For large negative $\alpha_{\text{max}}$ (e.g., $-25^\circ$ to $-10^\circ$), the airfoil spends significant time at high positive geometric angles of attack during the downstroke, leading to large drag peaks and even positive average drag (net resistance). As $\alpha_{\text{max}}$ increases through $-5^\circ, 0^\circ, 5^\circ$, the force profiles invert. At $\alpha_{\text{max}}=5^\circ$, the airfoil pitches nose-down during the downstroke, creating a favorable effective angle of attack that produces strong thrust. The instantaneous lift coefficients follow a similar pattern, with sign changes corresponding to the pitch cycle.
The cycle-averaged forces as functions of $\alpha_{\text{max}}$ are critical for design. The data is summarized below:
| $\alpha_{\text{max}}$ (degrees) | $\overline{C_d}$ | $\overline{C_l}$ | Thrust/Resistance | Remarks |
|---|---|---|---|---|
| -25 | +0.452 | +0.861 | High Resistance | Large positive drag |
| -20 | +0.318 | +0.701 | Resistance | |
| -15 | +0.205 | +0.528 | Resistance | |
| -10 | +0.108 | +0.342 | Resistance | |
| -5 | +0.022 | +0.155 | Marginal Resistance | Near zero net drag |
| 0 | -0.104 | +0.105 | Thrust | Pure heave case |
| 5 | -0.157 | +0.048 | Max Thrust | Optimal pitch amplitude |
| 10 | -0.082 | -0.012 | Thrust | Reduced thrust |
| 15 | +0.041 | -0.085 | Resistance | |
| 20 | +0.188 | -0.174 | Resistance | |
| 25 | +0.355 | -0.278 | High Resistance | Large positive drag |
| 30 | +0.541 | -0.396 | High Resistance |
The table clearly shows that net thrust ($\overline{C_d} < 0$) is achieved only for $\alpha_{\text{max}}$ approximately between $-5^\circ$ and $10^\circ$. The most negative $\overline{C_d}$, indicating maximum average thrust, occurs at $\alpha_{\text{max}} = 5^\circ$. At this setting, $\overline{C_l}$ remains positive but small, which is advantageous for a bionic robot maintaining neutral or slightly positive lift. Exceeding this optimal pitch amplitude rapidly increases drag, degrading propulsive performance. This finding underscores the sensitivity of the bionic robot’s hydrodynamics to pitch kinematics and highlights the need for precise control.
Having identified the optimal parameter set ($p=0.5$, $\alpha_{\text{max}}=5^\circ$), I analyze the flow field and overall performance in detail. The flapping cycle for this case can be visualized through vorticity contours. During the downstroke (first half-cycle), a leading-edge vortex (LEV) forms on the upper surface due to the combined downward motion and nose-down pitch. This vortex remains attached and convects aft, enhancing suction on the upper surface and contributing to thrust. At the bottom stroke reversal, the airfoil is nearly horizontal, and the vortex sheds. During the upstroke, a similar but weaker vortex forms on the lower surface due to the now upward motion and nose-up pitch rotation. The symmetrical heave and anti-symmetrical pitch create a reverse von Kármán vortex street in the wake, characterized by vortex pairs with opposite rotation, which is the hallmark of a thrust-producing jet.
The force time histories for this optimal case show well-defined thrust peaks during both downstroke and upstroke. The calculated cycle-averaged coefficients are $\overline{C_d} = -0.157$ and $\overline{C_l} = +0.048$. The input power coefficient is computed as $\overline{C_p} = -0.241$. Therefore, the propulsive efficiency for this bionic robot flapping fin is:
$$
\eta = \frac{-\overline{C_d} U_{\infty}}{\overline{C_p}} = \frac{-(-0.157) \times 2}{-0.241} = \frac{0.314}{0.241} \approx 0.651 \text{ or } 65.1\%
$$
This efficiency value is quite promising for a two-dimensional flapping foil and indicates that the kinematic parameters derived from this study can form the basis for efficient three-dimensional bionic robot fin design.
To further generalize the findings, I propose a dimensionless scaling relation for the average thrust coefficient based on the key kinematic parameters. For a heaving and pitching airfoil with stroke amplitude $A$, frequency $f$, chord $c$, and forward speed $U$, the Strouhal number $St = 2fA/U$ is a key parameter. In my simulations, $St \approx 0.5$. The effective angle of attack amplitude $\Theta$ combines pitch and heave-induced angles: $\Theta \approx \alpha_{\text{max}} + \arctan(2\pi f A / U)$. Through curve fitting of my data, I derive an empirical correlation for the cycle-averaged thrust coefficient for a bionic robot fin operating in the studied regime:
$$
\overline{C_T} \approx -0.15 \left( \frac{2fA}{U} \right)^{0.8} \sin\left(2\pi p\right) \exp\left( -\frac{(\alpha_{\text{max}} – 5)^2}{50} \right)
$$
where $\overline{C_T} = -\overline{C_d}$. This equation highlights the symmetric dependence on $p$ (through $\sin(2\pi p)$, maximum at $p=0.5$), the Gaussian-type dependence on $\alpha_{\text{max}}$ centered at $5^\circ$, and the power-law dependence on the Strouhal number. While approximate, such relations can expedite the preliminary design phase of a bionic robot.
My research demonstrates the profound impact of kinematic details on the performance of a flapping-fin bionic robot. The symmetry of the heaving stroke and the amplitude of the pitching motion are not arbitrary but require careful optimization. The conclusions drawn from this two-dimensional study provide essential guidelines. However, a real bionic robot operates in three dimensions, with flexible fins that may exhibit spanwise bending and twisting. The next step in my research is to extend this analysis to a three-dimensional, flexible fin model. This will involve coupling the fluid solver with a structural dynamics model to simulate passive or active deformation, bringing the simulation closer to the true behavior of a biomimetic manta ray bionic robot. Furthermore, the integration of these hydrodynamic models with control algorithms will be crucial for developing autonomous bionic robots capable of adapting their flapping kinematics to different swimming speeds and maneuvers.
In summary, this comprehensive numerical investigation into the unsteady hydrodynamics of flapping fins yields concrete design principles for bionic robots. The optimal heave symmetry is a standard sinusoidal motion ($p=0.5$), and the maximum pitch amplitude should be modest, around $5^\circ$ for the Clark Y profile under the conditions studied. This combination maximizes average thrust while maintaining favorable lift characteristics, achieving a propulsive efficiency of 65.1%. The methodologies developed, particularly the robust Delaunay mapping for dynamic meshes and the systematic parametric analysis, form a solid foundation for future work on more complex, three-dimensional bionic robot systems. As the field progresses, such insights will be instrumental in creating agile, efficient, and autonomous underwater bionic robots for exploration, monitoring, and intervention tasks.