Three-dimensional wake topology and propulsive performance of low-aspect-ratio pitching-rolling plates

Flapping motion is widely utilized by many types of biological propulsors, including fish pectoral fins^1–3 and insect/bird wings.^4,5 Such motion generally consists of an oscillatory rolling motion about a fixed joint with a simultaneous change in the geometric pitch angle with respect to the spanwise axis. In previous studies, this pitching-rolling motion was simplified as a two-dimensional pitching and/or heaving problem by assuming that the propulsor’s aspect ratio (AR) was sufficiently large. The associated wake structures and hydrodynamic performance have been experimentally^6–8 and computationally^9–13 studied in detail. An increase in the Strouhal number (St) has been found to induce a transition of the wake pattern from a von Kármán vortex street to an inverse von Kármán vortex street that characterizes propulsive wakes. At higher St numbers, the shed vortex pairs propagate at an angle to the streamwise axis, and the inverse von Kármán vortex pattern gives rise to an asymmetric wake.^14–18 However, for finite-aspect-ratio plates undergoing identical motions, the asymmetric wakes are absent due to the suppression of vortex coupling by the three dimensionality and the symmetric circulation distribution along the interconnected vortex loops.¹⁹ 

Numerous studies have been conducted to understand the wake structures and the propulsive performance of finite-aspect-ratio plates undergoing pitching and/or heaving motions. In particular, flow visualization was used by von Ellenrieder et al.²⁰ to investigate the wake structure produced by low-aspect-ratio plates (AR = 3) under pitching-heaving motion at Re = 164 and Strouhal numbers from 0.2 to 0.4. The wake model in their work indicated that a pair of merged vortex rings shed each half cycle, forming a zig-zag chain. This observation was then numerically verified by Blondeaux et al.²¹ using a low Strouhal number (St = 0.175) flapping plate. However, at a higher Strouhal number (St = 0.35), they found that the flow wakes consisted of two sets of vortical structures. A similar wake pattern transition was reported by Buchholz and Smits²² in a study on low-aspect-ratio pitching panels (AR = 0.54). Using dye visualization and digital particle image velocimetry (DPIV), a chain of vortical structures was observed at St = 0.23 and Re = 640. With increasing Strouhal number, the wake began to bifurcate and formed two branches of horseshoe-like vortices. Further exploration²³ has revealed that although distinct wake structures were observed with varying Strouhal number, the key features of the wake pattern are highly robust, even for Reynolds numbers on the order of 10⁴. Moreover, Dong and his colleagues^24,25 numerically studied the wake topologies of pitching-heaving ellipsoidal foils over a much broader range of parameters, including varying aspect ratios from 0.64 to 5.09, Strouhal numbers from 0.4 to 1.2, and Reynolds numbers in the range of 100 ≤ Re ≤ 400. Similar bifurcated wake structures have been observed and found to be quite robust for all low aspect ratios.

The aforementioned 2D and 3D pitching-heaving studies have given insights into the hydrodynamic mechanism in fish caudal fins and marine mammal flukes. In particular, when a rigid foil undergoes a pitching-heaving motion, symmetric spanwise flow is generally expected due to the unchanged instantaneous effective angles of attack along the foil span. However, when a foil performs a combined pitching and rolling motion about a fixed hinge, the effective angle of attack with respect to the incoming flow varies along the foil span. As a result, an asymmetric leading-edge vortex with an increase in strength from the foil root to its tip is expected to be formed during the pitching-rolling motion. However, the fundamental flow physics and wake structures of low-aspect-ratio pitching-rolling foils are not fully investigated. Only a few studies can be found on bio-inspired pitching-rolling propulsors, and most of them have focused on the propulsive performance. For instance, Techet²⁶ experimentally studied a robotic pitching-rolling foil with an NACA 0012 cross section and found that the high thrust and high propulsive efficiency of the foil could only be achieved when the foil flapped with a smaller maximum angle of attack (α_max = 30°) and high Strouhal number (St = 0.6). The high rolling amplitude did not appear to improve the thrust but could reduce the power losses. In another experimental study, Bandyopadhyay et al.²⁷ explored low-aspect-ratio flapping fins (AR = 3) at chord Reynolds numbers between 3.6 × 10³ and 1.5 × 10⁵ in a tow tank. The authors demonstrated that low-aspect-ratio pitching-rolling fins can only produce thrust within a bounded Strouhal number (0.2 < St < 1.5) and that their optimal efficiency can be achieved by tailoring the Strouhal number and pitch amplitude. In addition to these experimental studies, the wake structures of bio-inspired propulsion have only been discussed in relation to flow simulations of fish pectoral fins (Dong et al.²⁸ and Liu et al.²⁹) and rigid foil with an aspect ratio of 2.55 (Bozkuttas³⁰). These works indicated the existence of distinct flow features differing from those of the aforementioned pitching and/or heaving plates. In particular, the strong spanwise flow caused by the rolling motion resulted in the formation of stronger leading-edge and tip vortices, and the interactions among those vortices formed substantially more complicated vortex structures that convected into the wake. Their results also qualitatively demonstrated that pitching-rolling plates could serve as a better canonical model for investigating the hydrodynamics of bio-inspired flapping propulsion. Additional quantitative measurements and analysis are needed to understand the vortex dynamics and hydrodynamic performance of pitching-rolling propulsors.

To this end, we extend our previous studies on pitching-heaving plates^24,25 to investigate the wake topology and propulsive performance of low-aspect-ratio plates undergoing a pitching-rolling motion over a range of Strouhal numbers, Reynolds numbers, aspect ratios, and phase differences between pitching and rolling motions. The three-dimensional incompressible Navier-Stokes equations are solved numerically; therefore, all unsteady viscous effects and spanwise effects are included. The current paper begins by describing the simulation methodology and grid independence study. This is followed by a detailed discussion of the wake topology and propulsive performance. Finally, the conclusions are given in Sec. IV.

In the current study, flow fields were generated by a direct numerical simulation of the three-dimensional unsteady, viscous incompressible Navier-Stokes equations, as written in the following:

where u_i (i = 1, 2, 3) are the velocity components in the x-, y-, and z-directions, respectively; p is the pressure; and Re is the Reynolds number.

The equations are discretized using a second-order central difference scheme on a non-uniform Cartesian mesh, where the velocity and pressure are collocated at the cell centers. The unsteady equations are solved using a fractional step method, which provides second-order accuracy in time. An Adams-Bashforth scheme and an implicit Crank-Nicolson scheme are employed for the discretization of the convective terms and diffusion terms, respectively. Boundary conditions on immersed bodies are imposed through a “ghost-cell” procedure, which can handle bodies with finite or zero thicknesses. Thus, flow simulations with complex boundaries are conducted on stationary non-body-conformal Cartesian grids. This arrangement eliminates the need for the complicated re-meshing algorithms that are usually employed by conventional Lagrangian body-conformal methods. Additional details and flow validations for this code can be found in Ref. 31. The current flow solver has also been used to model canonical rotational plates,^32,33 pitching-heaving plates,^24,25 pitching-rolling plates,³⁴ ray-inspired underwater vehicles,²⁹ the flapping wing of a cicada,^4,35 and the pectoral fin of a bluegill sunfish.^2,28

The current study employs elliptic plates with a finite aspect ratio (AR) and a plate thickness of 3% of the chord length (a similar plate geometry was also used in Ref. 25). The elliptic plate is defined by two major diameters: the chord length c and the span length S. The AR of the plate is defined as the square of the span (S) divided by the planform area of the plate (A_plan) and equals 4S/πc. The surface of the plate is represented by a fine unstructured mesh with triangular elements.

The three-dimensional plate undergoes combined rolling and pitching motions in the yz-plane, as shown in Fig. 1. In Fig. 1(a1), the plate is rolling upward from the lowest position (“i”) to the highest position (“iii”) with an instantaneous change in pitching angle. The maximum pitching angle is reached at the mid-upstroke position (“ii”). To better illustrate the pitching-rolling motion, a schematic of the corresponding mid-chord and mid-span line kinematics are also provided in Figs. 1(a2) and 1(a3), respectively. Similar to the upstroke motion, Figs. 1(b1)-1(b3) present the downward rolling motion of the plate and the minimum pitching angle is obtained at the mid-downstroke position (“iv”).

The rolling motion of the plate with respect to the x-axis is given by

where ϕ(t) is the instantaneous rolling position at time t, A_ϕ is the rolling amplitude, and f is the flapping frequency. The negative sign in front of A_ϕ means that the rolling motion of the plate starts from the lowest position of a stroke.

The pitching motion with respect to spanwise axis is defined as

where θ(t) is the instantaneous pitching position, A_θ is the pitching amplitude, ψ is the phase difference angle between pitching and rolling motion (ψ = 90° for the baseline case), and θ_bias is the static pitching bias angle and is set to 0° in the current study for generating zero mean lift.

Hence, the Reynolds number is defined as Re = U_∞c/ν, based on the incoming flow velocity (U_∞) and the chord length (c). The Strouhal number is defined as St = 2A_ϕR_avgf/U_∞, where A_ϕ is the rolling amplitude and R_avg is the average rotational radius (R_avg = S/2) in the spanwise direction. This definition is consistent with previous experimental studies.^27,36

The kinematic parameters adopted in this work are, for instance, a reasonable representation of the rotation of a bluegill sunfish pectoral fin.^28,37 Table I provides a concise summary of all the parameters involved and their ranges. The time history of the plate rotational angle is shown in Fig. 2.

The grids employed in the current study are designed to provide high resolution in the region around the flapping plates as well as in the wake region, as shown in Fig. 3. At the left-hand boundary, we provide a constant inflow velocity boundary condition. The right-hand boundary is the outflow boundary, which is provided with a zero streamwise gradient boundary condition for the velocity, allowing the vortices to convect out of this boundary without significant reflections. The zero-stress boundary condition is provided at all lateral boundaries. A homogeneous Neumann boundary condition is used for the pressure at all boundaries.

To quantify the hydrodynamic performance, the forces and power coefficients are defined as

where T, L, Z, and P are the thrust, lift, spanwise force, and hydrodynamic power, respectively, and A_plan, the planform area of the plate, is equal to πcS/4. The instantaneous hydrodynamic power (P) is defined as the surface integration of the inner product between the pressure and the velocity in each discretized element.

Moreover, the propulsive efficiency is defined as

where $T̄$ is the cycle-averaged thrust and $P̄$ is the cycle-averaged hydrodynamic power, in which only the positive power is considered.

The computational domain size of all simulations was chosen as 30c × 25c × 25c. This choice was based on our experience with the simulation of such flows and test simulations on a number of different domain sizes. To maintain consistent grid resolutions in both the y- and z-directions, the nominal grid size employed in the current simulations ranged from 235 × 233 × 113 for the smaller aspect-ratio plate (AR = 1.27) to 235 × 273 × 145 for the larger aspect-ratio plate (AR = 2.55). Fig. 3 shows a typical grid used in the current study for the small-aspect-ratio case. As observed in this figure, non-uniform grids with two layers of refined meshes were used in the current study. Very dense meshes were used in a cuboidal region around the plate in all three directions, with the smallest resolution being Δx = 0.03c. Outside of this region, a layer with slightly coarser meshes (Δx = 0.05c) was arranged in all three directions to resolve the complex wake structures behind the plates. Beyond this layer, the grid is rapidly stretched in the y-direction and the z-direction. In the x-direction, the stretching is rapid upstream of the plate, where we do not expect any streamwise gradients. In the wake region, however, the stretching factor was maintained at 3% to ensure a relatively high streamwise resolution. This arrangement can limit the numerical dispersion that is associated with the use of central difference schemes on highly stretched meshes.³⁸ 

Comprehensive studies have been conducted to assess the effect of the spatial and temporal grid resolutions on the salient features of the computed flow and to demonstrate that the chosen grids produce accurate results. Spatial grid refinement studies were conducted by simultaneously doubling the grid in all three dimensions in the refined zones. The overall grid size for this refined grid was 283 × 273 × 169 (approximately 13.1 × 10⁶ grid points). Temporal independence studies were conducted by halving the Δt of the original case while using the same nominal grids. All these simulations were conducted on a Xeon E5-2609 CPU with 128 GB of core memory, and each nominal grid simulation required anywhere from 96 to 120 CPU hours.

Table II compares the key hydrodynamic quantities for the spatial and temporal grid independence studies conducted for the AR = 1.27, St = 0.6, Re = 200, and ψ = 90° case. The table shows that the spatial and temporal grid refinements lead to a difference in the mean thrust and efficiency of less than 1% and a variation in the root mean square (r.m.s.) lift and spanwise force of approximately 2% at most. This clearly demonstrates that the hydrodynamic forces computed in the current study were grid independent. In addition to the above comparison of the hydrodynamic forces and efficiency, it is useful to assess the effect of the grid on the flow development in the wake. In Fig. 4, we compare the wake profiles for these three different grids. Figs. 4(a) and 4(b) show the mean streamwise (⁠$u¯1−U∞$⁠) and transverse (⁠$u¯2$⁠) velocity profiles at the mid-span, and Fig. 4(c) shows the spanwise mean velocity profiles at the center of the xz-plane. In all of these plots, we note that the differences between the profiles computed using the three different grids are negligible. Also plotted in Fig. 4(d) are profiles of the fluctuating kinetic energy, defined as $(u1′2¯+u2′2¯+u3′2¯)/2$⁠, in the mid-span plane. A comparison of these profiles is an even more aggressive test of grid dependence because the fluctuations inherently contain more information from the smaller spatial and temporal scales in the flow. The difference between the three sets of profiles is small, which clearly establishes the fidelity and accuracy of the current simulations.

To better evaluate the effects of different grids on the instantaneous vortical structure, we also quantified the size of the vortex loop by measuring the slides of the 3-D wake topology in both near and far wakes at t/T = 4.25, as shown in Fig. 5. Specifically, the distances between each pair of vortices were measured. These distances represent the diameters of the vortex loops formed in the downstream of the plate. Fig. 5(a) shows that the distances between the inner vortex pair and the outer vortex pair which are denoted as d₁ and d₂, respectively, in which the center of each vortex is chosen as the point with the highest value of the streamwise vorticity (ω_x). Similarly, in Fig. 5(b), the distances between another two vortex pairs shed during the previous half stroke are denoted as d₃ and d₄, respectively. Table III compares these distance quantities for the spatial and temporal grid independence studies conducted for the case with AR = 1.27, St = 0.6, Re = 200, and ψ = 90°. The table shows that the spatial and temporal grid refinements lead to a difference in the near-wake vortex distances of less than 2% and a variation in the far-wake vortex distances of at most approximately 4%. This clearly indicates that the instantaneous vortical structures computed using the current setup were also grid independent.

In this section, we first present a detailed discussion of the wake topology as well as the force production of a circular plate undergoing a pitching-rolling motion. Following this discussion, the effect of the Strouhal number, Reynolds number, aspect ratio, and phase difference angle on the vortex dynamics and associated propulsive performance are discussed for the range of parameters listed in Table I.

The wake topology and force production for the flapping circular plate (AR = 1.27) are examined at St = 0.6 and Re = 200 with ψ = 90° in this section. Fig. 6 shows the vortex shedding process near the plate trailing edge at nine phases. The corresponding plate kinematics can be found in Fig. 2. The shell and core of the vortex structures are visualized using the Q-criterion³⁹ with two iso-surface values, namely, Q = 1 (in grey) and 6 (in color), respectively. The vortex cores are colored according to ω_x to indicate the vorticity convection direction. The senses of the vorticity in each section are indicated by arrows according to the right-hand rule.

As the plate rolls downward from the middle point to the lowest point (Figs. 6(a)-6(c)), a “C”-shaped vortex loop is found to be formed by the root vortex (V₁), the tip vortex (V₂), and the trailing-edge vortex (TEV). Note that the rolling motion creates different shear rates between the plate root and tip and affects the development of the vortex loop. The strength asymmetry between V₁ and V₂ can be easily observed in Fig. 6(c) when the plate reaches to the lowest point. When the tip of the plate starts rolling upward from the lowest position to the middle position of the cycle, the leading-edge of the plate is also pitching up (as shown in Figs. 6(d)-6(f)). During this period, newly developed vortices, resulting from the trailing-edge shear layer (V₃) and the leading-edge vortex (V₄), merge together and form a new inner “C”-shaped vortex loop with an opposite sense to the outer loop (V₁ and V₂). It is worth noting that the newly formed inner “C”-shaped vortex loop is a completely different structure from the contrails generated by the previous vortex loop. Additional visualization of the formation and evolution of the wake topology can be found in the supplementary video (see Fig. 6).

Note that the “double-C”-shaped vortex structures generated around the plate trailing edge are distinct from the observations in previous studies on purely pitching,^22,23,40,41 purely heaving,¹⁹ and pitching-heaving^21,24,25 panels/foils. In those previous studies, only single interconnected vortex loops were observed around the trailing edge during each half cycle, and the loops exhibited a symmetry about the spanwise central plane.^19,25 In the current study, however, the rolling motion results in a faster development of the leading-edge vortex (V₄). A strong leading-edge vortex is formed within only one quarter cycle and is then stretched and tilted along the tip connecting to the trailing-edge shear layer (V₃). This flow feature distinguishes the current wake structures from those in the previous pitching and/or heaving studies.

Along with the plate continuously moving upward, both “C”-shaped vortex loops started to shed outward from the plate’s trailing edge and connected with each other, as shown in Figs. 6(g)-6(i). As a result, a “double-C”-shaped vortex ring was formed and convected downstream. At the same time, a new outer “C”-shaped vortex loop starts to form. During the second half cycle, another “double-C”-shaped vortex ring (as shown in Fig. 6(a), which was formed a cycle earlier) with the opposite sense formed in the same manner. A pair of “double-C”-shaped vortex rings were then produced from each flapping cycle and resulted in a bifurcated wake pattern in the downstream (Fig. 7(a)). The side view (Fig. 7(b)) of the bifurcating wake pattern is similar to that observed in previous studies on pitching and/or heaving plates.^19,23,24 However, according to the top view (Fig. 7(c)), two key features differentiate the current wake structures from those of previous studies. First, in previous pitching and/or heaving studies,^19,23,24 the vorticity generated by the streamwise edges dominates the evolution of the wake. The stronger tip vortices induce a compression on the associated spanwise vortices. As the vortex rings convect downstream, a spanwise narrowing of the wake can be observed from the top view (i.e., Fig. 16 in Ref. 23 and Fig. 9 in Ref. 24). However, a similar narrowing process was not found in the downstream wake of the present work. Instead, the size of the vortex rings gradually increased as they convected downstream, as can be observed in Fig. 7(c) by comparing the size of the shed vortex rings R₁, R₃, R₅, and R₇ in the formation sequence. This difference can be explained by the increased spanwise flow caused by the rolling motion. For the pitching-rolling plate, the intensification of the leading-edge vortex and the diminishment of the root vortex attenuate the compression conditions found in the pitching-heaving plates²⁴ and result in the elongation of the vortex rings. Second, as the vortex rings convect downstream, a slight wake deflection in the spanwise direction can be observed from the top view of the wake topology (Fig. 7(c)). Furthermore, note that there are two sets of vortex “contrails” (see, for instance, the rings in Fig. 7(b)) at its upstream and downstream ends that extend toward its two adjacent rings. As the vortices convect downstream, these vortex contrails become weaker and gradually annihilate. Simultaneously, the “double-C”-shaped vortex rings gradually evolve into single-loop vortex rings (see, for instance, ring R₃ in Fig. 7).

To better understand the structure of the “double-C”-shaped vortex ring and the interaction between two adjacent vortex rings, several yz-plane slices were cut into the wake when the plate was at the middle of its upstroke (t/T = 4.25), and the streamwise vorticity contours (ω_x) are plotted in Fig. 8. Red and blue indicate the counter-clockwise and clockwise rotation directions, respectively.

In the 2-D contours of Fig. 8, the contained vortices of the “double-C”-shaped rings R₈ and R₇ are labeled as V₁ ∼ V₄ and V₅ ∼ V₈, respectively. Figs. 8(a)-8(d) show the detailed structures of the inner and outer “C” loops of R₈, in which Figs. 8(a)-8(c) are at the upstream of the trailing edge and Fig. 8(d) is slightly behind the plate trailing edge. Figs. 8(e)-8(h) illustrate the vortical structure of the “double-C”-shaped ring R₇, which was shed from the previous half cycle with the opposite sense of R₈. Figs. 8(e) and 8(f) show the interaction between R₇ and R₈. It can be observed that the outer “C” loop (V₁ and V₂) of R₈ is connected with the inner “C” loop (V₇ and V₈) of R₇ by a pair of contrails in the upstream. The interaction between R₇ and R₆ can also be observed in Figs. 8(g) and 8(h), where another pair of contrails with alternating signs is in the downstream between the outer “C” loop of R₇ and the inner “C” loop of R₆. The wake pattern begins to repeat in the far wake (Fig. 8(i)). Furthermore, as shown in Figs. 8(j)-8(l), due to the viscous dissipation effect at such a low Reynolds number, the “double-C” vortical structures gradually evolved into single-loop rings.

Another distinct feature that can be observed in Fig. 8 is that the strengths of the vortex cores for both the inner and outer “C” loops are non-equivalent. This is because the rolling motion of the plate created a difference in the shear rate between the plate root and tip. As the plate rolled either upward or downward, the formed vortex loops were strengthened near the plate tip and weakened about the root, thus leading to an asymmetric ring shape. Consequently, the downstream wake pattern is found to be slightly defected from the midline in the span towards the plate tip, as shown in Fig. 7(c).

Fig. 9 shows the temporal variation in the force coefficients over two cycles when the forces reached a periodic stage. To provide a better understanding of the vortex forming/shedding process, the superimposed pitching and rolling motion curves are also included in Fig. 9. The cycle-averaged value of each force component is indicated by the dashed-dotted line. In Fig. 9(a), the thrust peaks occurred twice during each cycle at the instant when the plate was near the center of its trajectory. A slight drag was produced when the plate started to reverse its rolling direction. The maximum and mean thrust coefficients for this case were found to be 3.98 and 1.46, respectively. Figs. 9(b) and 9(c) show the other two force components produced by the plate. Both forces have equivalent positive and negative variations over a cycle. The calculated mean lift and spanwise forces are 0.0 and 0.1, respectively. Note that the peak values of the lift (⁠$CLmax=5.35$⁠) and spanwise force (⁠$CZmax=3.19$⁠) components are in a similar range to the maximum peak thrust. This is different from the findings of previous pitching-heaving studies,^24,25,42 in which the peak thrust is significantly smaller than the peak lift. In addition, two peaks of the thrust and lift production can be observed during each half stroke. The first peak has a smaller magnitude than the second peak. This phenomenon could be related to the double-loop vortex shedding process close to the plate trailing edge. It is worth noting that the double peak in the force production observed in the pitching-rolling plates is in agreement with previous studies on bio-inspired flapping wings/fins, including experimental measurements of insect wings^43,44 and numerical simulations of insect model wings^45,46 and fish pectoral fins.^28,47

In this section, we examined the effect of the Strouhal number on the wake topology and propulsive performance of the circular plate (AR = 1.27) at Re = 200 and ψ = 90°. Fig. 10 presents both side and top views of the wake topology for St = 0.4 and 0.8, respectively. This plot can be examined in conjunction with the same plots for the St = 0.6 case in Fig. 7.

As shown in Figs. 10(a) and 10(b), the lower Strouhal number case (St = 0.4) shows a decreased vorticity strength. In this case, the “double-C”-shaped vortex rings rapidly evolved to single vortex rings. In contrast, the higher Strouhal number case (St = 0.8) presents a delay process for this wake evolution, as shown in Figs. 10(c) and 10(d). In addition, the concomitant increase in vorticity strength led to the enhancement of the mutual induction between two adjacent vortex rings. As a consequence, the oblique angle (α) between two sets of vortex rings increased with increasing Strouhal number, which can be observed by comparing Figs. 10(a) and 10(c) with Fig. 7(b). This change in oblique angle can affect the efficiency of momentum transport in the downstream, implying the existence of an optimal Strouhal number in terms of propulsive efficiency. To further quantify the inclination angle of the vortex formation, both the wake oblique angle (α) and the vortex ring orientation angle (β) with respect to the wake centerline were measured in the near wake. As tabulated in Table IV, the oblique angle of the bifurcated wake increases monotonically with the Strouhal number. On the other hand, the angle of the vortex ring orientation initially follows an increasing trend, reaches a maximum value at St = 1.0, and then decreases abruptly. As the shed vortex ring convecting to downstream, the strength of the vortex ring is gradually decreasing due to the viscous dissipation. As a result, both the oblique angle and vortex ring orientation angle are decreased with the downstream distance. This phenomenon is especially clear for the high Strouhal number cases, i.e., at St = 0.8 Fig. 10(c).

Another salient feature that must be noted here is the wake deflection along the mid-span, which can be easily observed in Figs. 10(b), 10(d), and 7(c). For lower Strouhal number cases (i.e., St = 0.4 and 0.6), the wake always deflects towards to the plate tip direction. For the higher Strouhal number case (St = 0.8), however, the wake initially deflects to the tip direction in the near wake, and then gradually deflects back in the far wake. A higher self-induced vortex ring velocity for the higher Strouhal number cases might be responsible for this wake deflection phenomenon. In addition, the viscous effects in the far wake could also play a role in wake evolution. Furthermore, as the vortex rings convect downstream, the elongation of the vortex rings is proportional to the increment of the Strouhal number.

Fig. 11 shows the variation in the mean thrust coefficient and propulsive efficiency as a function of Strouhal number for the baseline case. First, it can be observed from this figure that for all cases, the thrust monotonically increases with the Strouhal number (0.4 ≤ St ≤ 1.2). Such behavior has been well documented for both two-dimensional^6,7 and low-aspect-ratio^23,24 pitching and/or heaving foils. In addition, the propulsive efficiency exhibits an initial rapid increase and then gradually decreases with increasing Strouhal number. The maximum value of the efficiency is reached at St = 0.6, which is higher than the range usually considered optimal for swimming and flying animals.^48,49 This is probably due to the particular choice of parameters, such as the maximum pitching angle and Reynolds number, in the current study. The relatively low Reynolds number in the current study leads to a proportionately large shear drag that has to be overcome for the plate to produce a net thrust. This tends to push the optimal Strouhal numbers to higher values. As shown in Ref. 50, even at high Reynolds number, depending on the maximum angle of attack, the optimal Strouhal numbers can vary from 0.3 to 0.6 for a rigid flapping foil.

In this section, two additional simulations, one at Re = 100 and the other at 400, have been conducted for the baseline case to assess the effect of the Reynolds number on the wake topology and instantaneous force generation. Note that the higher Reynolds number case was conducted on the finer spatial meshes.

Fig. 12 shows a perspective view of the vortex topology for these cases, and this can be compared with the corresponding plot in Fig. 6(f). For the lower Reynolds number of 100, the formation of the double-loop vortex around the plate trailing edge remained evident, although the loops rapidly dissipated after they were shed into the wake. This behavior is consistent with the increased viscous effect in this case. The higher Reynolds number case, on the other hand, shows many of the features observed for the case at Re = 200 in both the near and the far flow field. In particular, the wake was also found to exhibit “double-C”-shaped rings propagating downstream. This indicates that the basic vortex dynamics of these low-aspect-ratio flapping foils are insensitive to changes in the Reynolds number between 100 and 400.

Fig. 13 shows the time-varying force coefficients of the plate flapping at an aspect ratio of 1.27 for Reynolds numbers ranging between 100 and 400 with St = 0.6. The corresponding Re = 200 case is re-plotted for comparison. As shown in Fig. 13(a), the thrust coefficient increases with the Reynolds number. Specially, the peak thrust coefficient for the Re = 100 case is approximately 3.45, whereas that for the Re = 400 case is approximately 4.35, which amounts to an increase in magnitude of over 26%. The mean thrust coefficients for the Re = 100 and Re = 400 case are 1.17 and 1.65, respectively. For the lift coefficient (Fig. 13(b)), the double peaks appear for all Reynolds numbers but with slight differences in magnitude. In Fig. 13(c), the spanwise force exhibits a smaller negative peak with increasing Reynolds number. The mean values of both the lift and spanwise force coefficients are approximately zero.

In this section, we examine the effect of the plate aspect ratio (AR) on the wake topology and instantaneous force production. Fig. 14 shows the vortical structures for the elliptical plate with AR = 1.91 and 2.55 at St = 0.6, Re = 200, and ψ = 90°. This plot can be compared with the corresponding plot in Fig. 6(f). To illustrate the effect of plate aspect-ratio on the inner-loop development in Fig. 14, the 2D streamwise vorticity contours (ω_x) and spanwise vorticity contours (ω_z) are shown in Figs. 15(a)-15(c) and Figs. 15(d)-15(f), respectively. It shows that the size of the inner-loop in streamwise direction (which consists of V₃ and V₄) increases along with the increment of plate aspect-ratios (Figs. 15(a)-15(c)), however, the spanwise vorticity magnitude of the inner loop decreases (Figs. 15(d)-15(f)).

Fig. 16 shows the time history of the force coefficients for the case of different aspect ratios at Re = 200 and St = 0.6. As the aspect ratio increases, the thrust peak increases and the magnitudes of the other two transverse forces decrease. The mean thrust coefficients for the AR = 1.91 and AR = 2.55 case are 1.62 and 1.77, respectively. This is a clear indication of what we expect to be a relatively high propulsive efficiency for a larger aspect-ratio plate.

In Subsections III A–III D, we discussed the wake topology and force generation of pitching-rolling plates with phase difference angle (ψ) 90° between pitching and rolling motion at various Strouhal numbers, Reynolds numbers, and aspect-ratios. In the current section, we examine the effect of the phase difference angle on the vortical structures and hydrodynamic performance for an AR = 1.27 circular plate at St = 0.6 and Re = 200.

Fig. 17 shows both side and top views of the wake topology for ψ = 100°, 110°, and 120°, respectively. This plot can be examined in conjunction with the same view for the ψ = 90° cases in Figs. 7 and 10. Due to the self-induced vortex ring velocity, the wake topology of ψ = 90° cases shown in Subsections III A and III B (i.e., Figs. 10(b), 10(d), and 7(c)) presents a slight wake deflection in the spanwise direction towards the plate tip. However, as the phase difference angle (ψ) increases from 90° to 120°, the wake gradually deflects from the tip direction to the root direction as illustrated in the right column of Fig. 17. This wake deflection angle (γ) can reach up to −7° in ψ = 120° case. Here, the positive and negative signs of γ represent the wake deflection toward the plate tip and toward the plate root, respectively. Another salient feature of the wake topology is the wake oblique angle (α) increasing with the increment of phase difference angle. The vorticity iso-surface shows that slight changes in the phase difference between pitching and rolling motion can enhance the vorticity strength of the shed vortex ring and thus lead to an increased mutual induction between two adjacent vortex rings. The obtained wake oblique angles and wake deflection angles with respect to the phase difference angles are listed in Table V.

The enhancement of the vorticity strength due to the increased phase difference is also evident in the force generation. Fig. 18 shows the instantaneous force coefficients for ψ ranging between 100° and 120° together with the baseline case (ψ = 90°). This plot shows that the peak values of all three force components are significantly increased when the phase difference angle is varied. It should be noted that the production of large force peaks can be extremely useful for performing rapid maneuvers. Thus, the phase difference angle can be treated as a potential control parameter for bio-robotic design. In addition to the instantaneous force history, the cycle-averaged performance of different phase angles is evaluated and shown in Fig. 19. The figure shows that a slight change in the phase difference angle can increase both the cycle-averaged thrust and efficiency compared to the baseline case (ψ = 90°). In particular, the interpolation of these data indicates that the optimal propulsion efficiency can be obtained for a phase difference angle of approximately 105° (pitching leads the rolling motion), which can lead to the enhancement of the cycle-averaged thrust and efficiency of up to 23% and 15%, respectively. It is worth noting that current optimal phase difference angle (around 105°) is different from that found in the pitching-heaving foils.⁵¹ The previous numerical investigation of pitching-heaving foils found that the foils perform most efficiently when the phase shift angle was set up between 70° and 90°. This different phase angle preference is because of the unique wake structures produced by pitching-rolling plates comparing to that formed by the pitching-heaving plate. As shown in Fig. 17, the wake deflection angles (γ) are affected by the change of the phase difference angle. When the phase difference angle is above 90°, the wake convects to the downstream with the minimum deflection. Fig. 19 also shows that the plate reaches to the maximum efficiency between 100° and 110°.

Simulations of flow past thin pitching-rolling plates have been conducted using an immersed boundary solver. Despite the simplified kinematics adopted to represent natural flapping propulsors, the current study is expected to provide some general insights into the hydrodynamics of low-aspect-ratio propulsors.

The simulations have shown that, in the near wake close to the plate trailing edge, a distinct double-loop vortical structure with opposite vortex sense can be observed. As the vortex loops shed into the wake, these two loops become connected with each other and form an interconnected “double-C”-shaped ring during each half flapping cycle. The periodically shed vortex rings eventually formed a bifurcating wake structure with a slight deflection in the spanwise direction to the plate tip. In addition, the contained vortex ring size gradually increases as it convects downstream. Persistence of the aforementioned wake formation is shown for a range of Strouhal numbers (0.4 ∼ 1.2), Reynolds numbers (100 ∼ 400), and plate aspect ratios (1.27 ∼ 2.55). Further analysis has shown that the formation and shedding process of this unique vortical structure resulted in double peaks in the force history. In addition, the cycle-averaged thrust was also found to increase monotonically with the Strouhal number, and the maximum propulsive efficiency can be achieved at a Strouhal number of approximately 0.6. Changing both the Reynolds number and plate aspect ratio affected the peak values of the force production but had little influence on the mean force magnitude. Our further investigation of the effects of the phase difference angle between pitching and rolling motions indicated that the propulsive performance can be further improved if an appropriate phase difference angle is selected.

This research is supported by AFOSR Grant No. FA9550-12-1-0071, NSF Grant No. CEBT-1313217, and ONR Grant No. MURI N00014-14-1-0533.