A Comparative Numerical Study on the Performances and Vortical Patterns of Two Bioinspired Oscillatory Mechanisms: Undulating and Pure Heaving

The hydrodynamics and energetics of bioinspired oscillating mechanisms have received significant attentions by engineers and biologists to develop the underwater and air vehicles. Undulating and pure heaving (or plunging) motions are two significant mechanisms which are utilized in nature to provide propulsive, maneuvering, and stabilization forces. This study aims to elucidate and compare the propulsive vortical signature and performance of these two important natural mechanisms through a systematic numerical study. Navier-Stokes equations are solved, by a pressure-based finite volume method solver, in an arbitrary Lagrangian-Eulerian (ALE) framework domain containing a 2D NACA0012 foil moving with prescribed kinematics. Some of the important findings are (1) the thrust production of the heaving foil begins at lower St and has a greater growing slope with respect to the St; (2) the undulating mechanism has some limitations to produce high thrust forces; (3) the undulating foil shows a lower power consumption and higher efficiency; (4) changing the Reynolds number (Re) in a constant St affects the performance of the oscillations; and (5) there is a distinguishable appearance of leading edge vortices in the wake of the heaving foil without observable ones in the wake of the undulating foil, especially at higher St.


Introduction
The hydrodynamics and energetics of bioinspired oscillating mechanisms have received significant attentions by engineers and biologists over the past decades. Regarding the development of biomimetic robots, autonomous underwater vehicles (AUVs), and Micro Air Vehicles (MAVs), efficient control surfaces and propulsion mechanisms are required to produce propulsive forces and/or moments [1]. The oscillating mechanisms which are used by animals in nature can provide forces and moments efficiently [2][3][4]. This is due to the evolution process of these animals and their adaptation with operating medium [5]. Heaving (or plunging) and undulation are two significant mechanisms which are commonly used by body and/or fins of aquatic animals and wings of aerobic animals for propulsion, maneuvering, and stabilization purposes [2]. For example, in some fish species (some kinds of carangiform and Thunniform swimmers, e.g., some of the dolphins, tunas, or sharks), the caudal fin has pure heaving motions, while in some other species (anguilliform and subcarangiform swimmers like spinny dogfish shark) it has undulating motions [5,6].
In propulsion, flying animals generally use oscillatory mechanisms to generate both thrust and lift, while many swimming species use oscillatory mechanisms primarily to generate thrust [7]. Animals moving in water have usually either longitudinally slender (thin) or vertically (rarely laterally) flat body and caudal fins, which gives them a higher flexibility for propulsion, a lower drag to move against, and, hence, less need for exertion [2]. In many species, caudal fin provides whole or a bulk of the propulsive, maneuvering, and stabilizing forces. The caudal fin of these animals has a foillike geometry in 2D cross sections. An oscillating (rigid heaving or flexible undulating) caudal fin produces thrust under certain conditions [8]. The thrust is produced through the formation of an average velocity profile behind the oscillating fin which is in the form of a jet [9]. Because of the influence of the periodic excitation from the fin, a downstream moving staggered array of vortices, closely resembling a Karman vortex street but with reverse rotational direction, is formed in the wake [9]. The physics of thrust generation with oscillating motions was first presented by Durand [10] and after that observed in a number of experimental studies, for example, [11][12][13]. Rosen [14], Videler et al. [15], Lauder [5,16], Drucker and Lauder [17], and Tytell and Lauder [18] are also some of the researchers who observed experimentally similar patterns in the wake of different fish species. However, numerical studies have provided more details on wake structures and effective parameters of undulatory motion of flexible body or fins [19][20][21][22] as well as vortical pattern in oscillating heaving foils [7,9,[23][24][25].
Freymuth [11], Koochesfahani [12], Jones et al. [25], and Zhang et al. [26] showed through experiments and simulations that the wakes of oscillating airfoils can be characterized as drag-producing, neutral, or thrust producing mechanisms depending on the heaving/undulating frequency and amplitude. The general conclusion from the existing investigations is that the thrust force and propulsive efficiency strongly depend on the kinematic parameters of oscillations. Furthermore, the wake structure of the oscillatory mechanisms is closely related to the nondimensional parameter Strouhal number (St = 2 / ), where is the free-stream velocity, is the amplitude of oscillations, and is the oscillating frequency.
In this paper a NACA0012 foil, which oscillates with heaving and undulating mechanisms, is studied in three low, medium, and high Reynolds numbers, respectively, Re = 4000, Re = 40000, and Re = 400000 (Re = /], where is the characteristic length equal to the chord length of the foil and ] is the kinematic viscosity of the fluid). The oscillating frequency is varied from 0.001 to 0.05 for Re = 4000, from 0.01 to 0.5 for Re = 40000, and from 0.1 to 5 for Re = 400000, which are equivalent to the St in the range of [0.05-2.5] with maximum amplitude excursion of 0.1 . In present study, simulations are carried out for a typical foil geometry (NACA0012) which is commonly used in many applications as control surfaces or propulsors. In addition, a number of tests are accomplished to validate and evaluate the numerical solver in comparison with experimental data.
Several experimental and numerical works have separately studied the wake structure, energetics, and flow characteristics of undulating and heaving foils including the lift, drag, and propulsion efficiency [12,13,[26][27][28][29][30][31][32][33][34][35]. However, this study attempts to present a comparison for the energetics, performance, and wake structures of the heaving and undulating mechanisms in an extended range of Strouhal numbers. Such comparisons could provide a wealth of data and consequences about the advantages and disadvantages of these mechanisms and also are helpful to select and design the appropriate mechanism for the specific applications. Furthermore, remarkable findings are brought out from this comparison which are presenting in Section 4. The range of Strouhal numbers reported here, that is, 0.05 < St < 2.5, for the undulating motion is wider than the well-known range in literature, that is, 0.2 < St < 0.4 [36]. This is in order for more observations and answering the following questions: why do the animals use the undulatory propulsion mechanisms in the reported range? And why do not they use lower or higher ranges of Strouhal numbers? In addition, the presented qualitative and quantitative data are useful in the designing and developing of bioinspired oscillatory mechanisms. Current research also brings out a comparison of the especial deformability of the undulating fins over the rigidity of the heaving fins which mostly exist in aquatic animals. Owing to the difficulty in testing the swimming parameters of live fish, experimental studies on the oscillating caudal fins are rare and hence such systematic numerical studies which can simulate the real conditions are essential.
The rest of this paper is organized as follows. The governing equations as well as methodology and implementation of numerical approach, including the enforcement of boundary conditions and the mesh motion and correction procedure, are described in the next section. In addition, the validation tests and calculation of parameters are presented. Then, the results including thrust, consumed power, efficiency, and wake structure of the heaving and undulating mechanisms are presented and discussed in detail. Finally, the major conclusions are presented.

Fluid-Solid Interaction
Problem. The oscillations of bodies such as heaving and undulation of foils involve moving and deforming boundaries, respectively. The arbitrary Lagrangian-Eulerian (ALE) approach [37] is employed to encounter the fluid and moving body interactions. From the ALE viewpoint, the nodes of the computational domain mesh may be moved with the continuum, same as normal Lagrangian manner, be held fixed in Eulerian fashion, or be moved in some arbitrarily specified way to give a continuous rezoning capability [38]. In addition, the simulation of conservation equations on the moving meshes requires the geometric conservation law (GCL) to be satisfied [39].

Governing Equations of Fluid.
The ALE formulations of the unsteady incompressible viscous flow, governed by the integral forms of Navier-Stokes equations, are as follows.
Conservation of mass: Conservation of momentum: where, in the above equations, is the fluid density, ⃗ V is the Cartesian velocity vector of the fluid, ( ) is an arbitrary  volume whose boundary ( ) = ( ) moves with the mesh velocity ⃗ V , ⃗ is the specific body force vector, and denotes the Cauchy stress tensor which can be calculated from where is the Kroenecker delta function, and are dynamic stresses (shear and normal stresses) which for an incompressible, isotropic, Newtonian fluid are expressed as where is the dynamic viscosity of the fluid. The two important nondimensional numbers in this problem are Reynolds number (Re), which is the ratio of inertial over viscous forces and characterizes steady inline motion of heaving and undulating foils, and Strouhal (St) number which is the ratio of unsteady to inertial forces and describes the structure of wakes (or vortical patterns) behind the heaving and undulating foils which oscillate with frequency (or period ) [40,41].

Simulation Algorithm.
The overall procedure for carrying out the computational simulations of the fluid-solid interaction problem of oscillating foil can be divided into the following major steps (Figure 1(a)): (i) Construction of the geometry and generating the mesh for foil and computational domain.
(ii) Defining the geometry and kinematics into the solver.
(iii) Defining the appropriate initial and boundary conditions.
(iv) Performing the unsteady CFD computations and obtaining required hydrodynamic forces.
(v) Updating the motion or deformation of the foil's geometry and moving or modifying the mesh with dynamic mesh resolution techniques.
(vi) Continuing the unsteady solution until steady periodical drag, lift, and moment values are observed or more periods of oscillations.

Solver.
The implicit method is conducted for the temporal discretization of the equations. As well, the second-order upwind scheme is utilized for the convective fluxes which 4 Applied Bionics and Biomechanics require the determination of the gradient of scalar properties [42]. Also, the least squares scheme is used in the solver for the gradient evaluation associated with the Total Variation Diminishing (TVD) limiter [43,44] to prevent numerical oscillations. In addition, a second-order accurate, centraldifference approximation is applied to the diffusion terms (dynamic stresses) in momentum equations. The collocated grid approach is employed in which all of the flow variables are stored at the same locations and thus, the number of coefficients that must be computed and stored is minimized. The pressure-velocity coupling is achieved by utilizing the Pressure-Implicit with Splitting of Operators (PISO) scheme [45] which, by using a relationship between velocity and pressure corrections, enforces the mass conservation and obtains the pressure field. Furthermore, a segregated pressure-based solver is employed to sequentially solve the governing equations. Figure 1(b) illustrates the flowchart of the fluid solver which is outlined as follows: (i) Starting the calculations at the new time step using the latest solution of velocity components and pressure field as starting estimates for current time.
(ii) Sequentially solving the momentum equations by utilizing the values of the pressure and face mass fluxes from previous time step.
(iii) Solving the pressure correction equation using the recently obtained velocity field and the mass flux.
(iv) Correction of face mass fluxes, pressure, and the velocity field using the pressure correction obtained from previous step.
(v) Solving the second pressure correction equation and correcting both velocities and pressure again.
(vi) Checking for the convergence of the equations; if the convergence was not reached returning to step 2 and repeating; otherwise, proceeding to next time step.
The convergence criterion is conducted so that the Root Mean Squares (RMS) of the pressure and the velocities are less than 10 −5 in the computational domain nodes.

Mesh Construction.
Combined structured and unstructured meshes of 20436 nodes and 30535 elements are constructed around the NACA0012 foil ( Figure 2) and also for the validation purposes, two other meshes with the same geometrical constructions are constructed by approximately half and twice the number of grid nodes of the mentioned mesh. The physical phenomena of fluid mechanics such as boundary layer and velocity and pressure gradients are essential factors to the construction of the mesh. The laminar boundary layer thickness is estimated with Blasius equation [46]. In order to decrease the computational costs of dynamic mesh refinements (as are described in coming sections), the solution domain is composed of two parts: a nondeforming region which remains unchanged during the deformations (Figure 2(b)) and a deforming region (circle) which is altered by the flexible motions of the foil (Figure 2(a)).

Boundary Conditions.
A limited solution domain is created around the body, whose boundaries are far enough to not affect the results at far-field boundaries, particularly the velocity field and the wake structure. The foil is located at 6 from the inlet in the axial direction ( is the chord length equal to 1 m), 8 from each transverse far-field, and 20 from the outlet. Simulations are carried out in Re = 4000, 40000, and 400000 for St in the range of [0.05, 2.5] to study the variations of instantaneous and time averaged forces, consumed power, efficiency, and flow patterns versus St. Inlet velocity, outlet pressure, and far-field free-stream velocity are prescribed. The inner oscillating zone and the outer steady zone in the domain mesh are created with a common overlapping interface describing an interior zone.

Mesh Motion and Correction
Techniques. The foil deforms flexibly by undulation motions and makes the surrounding flow strongly unsteady. In addition, the computational grids are rigidly moved with the heaving motions. To simulate this unsteady flow, dynamic grid techniques should be employed along with the corresponding unsteady flow solver. The boundary-conforming methods (dynamic grid methods) are still more popular in CFD community including the simulations of the external biofluid dynamics [26]. The mostly used approach in this category is "spring-analogymethod" (SAM) [47], which is an algorithm to move the grid nodes surrounding the deforming solid walls supposing the cell edges as springs. For large deformations, other mesh refinement methods such as adaptive local refinement [48,49] should be used accompanied with SAM to resolve the mesh distortion problems. On the other hand, rigid-body motions use the body fitted mesh approach [50] to impose the translational motions of heaving in such a way that whole domain translates with the translations of the rigid boundary.
where is the moved -coordinate of the point ( , ) because of heaving with amplitude . In all of the simulations, the normalized maximum heaving amplitude is 0.1, in accordance with the maximum excursion of the undulating motions.
The formulations of the undulatory motions for the chord line or the backbone of the virtual fish are as (6) (Figure 3(b)). Consider with the coefficients 0 = 0.02, 1 = −0.08, and 2 = 0.16 to match the experimental curve of a typical carangiform swimmer [21,51]. In all of the simulations, the normalized wavelength / is 1.

Wake Visualization.
The most important characteristics of the flow are location of the vortex cores and distribution and breakdown of the vortices generated from the foil's surface. Two methods are employed to visualize the vortical patterns of the moving/deforming foils: vorticity criterion and criterion [52]. Vorticity criterion uses the values of vorticity computed by The criterion is based on the value of the second invariant of the velocity gradient tensor ∇ ⃗ V [52] as follows: where Ω and are the antisymmetric and symmetric components of ∇ ⃗ V, respectively, which are defined as follows: Physically, can be a balance between the strain rates ( = √ ) and the rotation (Ω = √Ω Ω ). Thus, positive values of indicate regions where the strength of rotation (rotation rate) dominates the strength of strain (strain rate) and hence, -isosurfaces can denote the vortex envelopes [52].

Forces, Consumed Power, and Efficiency Calculations.
Momentum transfer of the foil to the surrounding water and vice versa is via drag, lift, and thrust [40] which are produced from the pressure and velocity distributions (provided from the fluid solver) on foil's body due to the flapping and undulating motions. The foil has unsteady or time dependent periodic lateral movements and can also have unsteady or steady inline motions. Inline unsteady motions are known as accelerating or braking. For the steady propulsion, the forces and moments acting on the foil are balanced. The pressure forces are calculated by integration of the pressure distribution (which is normal to the surface) over the foil's surface. For an incompressible flow over an impenetrable surface, the normal stress force is zero and integration of only the viscous shear stresses over surface results in the viscous forces. The integrations are as follows: where V and are the viscous and pressure forces, respectively, Γ is foil's surface, ( ̸ = ) are tangential components of the stress tensor on the surface, is the unit normal vector on the surface element ( ), is the pressure value on the foil. Thrust and drag forces are the sum of all forces in the direction or the counter direction of motion, respectively.
The propulsion performance of an oscillating airfoil is represented with three important parameters: mean input power ( ), mean inline force ( ), and propulsion efficiency ( ). The mean input power is expressed as where ℎ( )/ and ( ) are the transvers displacement rate and the force component in the direction of -axis, respectively. / and ( ) are rates of changes in the pitching angle and the value of pitching moment in plane, respectively. Finally, the overall propulsion efficiency is calculated by

Validation, Grid Dependency Study, and Time Step
To verify the accuracy, convergence, and stability of the solver, a series of validation tests are performed along with the grid dependency tests. Also, the value of time step for each simulation is selected in such a way that the solver can capture the smallest grid size with respect to the inlet and moving grid velocities.

Validation of Solver.
A high-aspect-ratio foil, illustrated at Figure 13, with chord length , moves at constant forward speed , performing a heaving motion, ℎ( ), of amplitude ℎ 0 and frequency , and a pitching motion, ( ), of amplitude 0 and frequency . The pitching motion has a phase lead with respect to the heaving motion, which is denoted by . The one-third-chord point is the pivot point. The true angle of attack profile can be calculated mathematically as follows: The foil oscillates in a flapping motion around one-third chord with the heaving amplitude of 0.75 ( = 0.1 m). Two cases are used for the validation purposes. The maximum pitch angle for = 90 ∘ can be determined approximately as follows [30]: In  Figures 17 and 18 summarize the results for the grid dependence study. Here the results are presented for the case of Re = 40000 using three grid sizes of 10124, 20436, and 40732 nodes. All tests of Cases 1 and 2 have repeated to evaluate the dependence of results on domain resolution. Based on this grid dependence study, it is concluded that the medium grid size of 20436 is sufficient for carrying out the simulations which are insensitive to grid size.

Results and Discussions
The simulations are carried out for the cases of heaving and undulatory oscillations at Re = 4000, 40000, and 400000, a maximum trailing edge peak to peak amplitude of 0.

Instantaneous Inline and Transverse Forces and Pitching
Moment. The instantaneous inline and transverse hydrodynamic forces and pitching moment, which are shown in Figure 4, are presented to bring out qualitative insights of the differences in time varying results at different Strouhal numbers. The foils are assumed to be restrained by rigid bases, where the net hydrodynamic forces and the pitching moment are absorbed by these hypothetical bases. Therefore, the inline forces, shown in Figure 4(a), are net forces that would be available to accelerate the foil either forward or backward, depending on the sign of the mean time value, when the hypothetical bases are removed. The forward or positive acceleration occurs when the propulsive or thrust force exceeds the resistive or drag force. In contrast, when the drag force dominates the thrust force, the negative or backward acceleration occurs.
From Figure 4, it can be found that, to reach a certain value of inline force (thrust), much energy is lost due to oscillations of the traverse and opposing inline forces and moment in each cycle. Since all the existing thrust producing mechanisms suffer from the energy losses, the more efficient mechanism among them should be selected. The energy losses mostly originate from the vortical structures of oscillating mechanisms. The vortical structure is created due to interactions of the foil body and its surrounding fluid. Hence, the differences between the oscillation mechanisms would originate from the differences of their wake structures which inherently are due to their different kinematics and rigidity/flexibility. The approximate ratios of the peak amplitudes of oscillations of the instantaneous inline and transverse forces and the moment of the heaving mechanism over the undulating mechanism are, respectively, 0.

Undulation-Re = 4000
Undulation-Re = 40000 Undulation-Re = 400000 Heaving-Re = 4000 Heaving-Re = 400000 Heaving-Re = 40000 as shown in Figure 5(a), the curves corresponding to the undulating foil descend with lower slopes than the heaving foil. The differences of descending slopes show that, in given Reynolds and Strouhal numbers, heaving motions produce higher thrust. In particular, at higher St, the produced thrust by heaving foil is much more than the undulating foil in a given Re. It means the undulating mechanism has some limitations to produce high thrust forces. However, considering solely the amount of thrust production is a single attitude to the performance. Thus, for a comprehensive evaluation of the performance of these mechanisms, other parameters such as power consumption and efficiency should also be considered.
Comparing the curves shown in Figure 5, increasing the Reynolds number of flow from 40000 to 400000 causes the frequency of undulating and heaving motions to be increased  in order to make the Strouhal number remain constant. Hence, as it is seen, the time averaged inline forces are increased due to increasing the Re. In contrast, by decreasing the Re and keeping the St constant, the produced inline forces are also decreased. In addition, the change of Re causes the movement of the steady state Strouhal numbers St * and St * , forward or backward, respectively, by decreasing or increasing Re.

Consumed Power.
The consumed power ( Figure 6) also descends monotonously with increasing the St, due to the direction of the -axis. However, the meaning of single sign of the power is that, to oscillate a foil, always some power must be consumed. As Figure 6 indicates, the power curve of the undulating foil descends with a lower slope than the heaving foil. In other words, in given St and Re, especially at higher St, the heaving mechanism needs much more power than the undulating mechanism. However, at low St (approximately St ≤ 0.3), nearly same powers are required to oscillate the foils with both undulating and heaving kinematics. It is also confirmed by looking at the similar vortical patterns shown in Figures 8 and 11. Furthermore, by comparison of the rates of changes in the power curves with inline force curves, in a given Re, it could be found that the difference between the descending slopes of power curves grows more intensely. The significances of power consumption and energy costs are highlighted in the autonomous underwater vehicles, because of limited availabilities to energy in their operational domains and also their limiting design conditions such as weight. In addition, increasing the Re from 4000 to 40000 and keeping the St constant causes increase of the power consumption of the oscillating foils with great differences in higher Re. After the peak efficiency point, at higher St, the value of efficiency decreases monotonously from the peak propulsive efficiency. The ascending and descending behaviors of the efficiency curves, from lower to higher Strouhal numbers, can be attributed to the vortical patterns of the oscillation mechanisms. However, the focus on these features requires scrutinized investigations of the vortical patterns and also the evolution of the vortices from forming to decaying stages at different Strouhal numbers.

Efficiency.
A comparison of the efficiency curves of heaving and undulatory mechanisms in a given Re brings out the fact that, at low Strouhal numbers, where the efficiency is less than the peak efficiency point of the undulating foil, the value of efficiency for heaving oscillations is greater than the efficiency of the undulating oscillations. After the peak efficiency point of the heaving foil, at higher Strouhal numbers, the efficiency of undulatory oscillations reaches higher values than heaving oscillations. The difference between the efficiency values of heaving and undulating foils is greater at regions near the peak efficiency point of undulating foil, while this difference decreases by receding from that point. Mostly higher efficiency behavior of undulation mechanism, comparing with the heaving mechanism, shows a good potential of applying this mechanism to design of the underwater vehicles and the stabilizing control surfaces.
By comparing the efficiency curves in different Re, it is obvious that changes in Re move the peak efficiency points and alters the efficiency values. For undulating motion, increasing the Re from 4000 to 400000 increases . But this behavior does not occur for heaving motion, since the maximum is for Re = 40000 and the minimum is for Re = 4000.

Wake
Structures. The wake of oscillating foils has already been studied experimentally using Particle Image Velocimetry (PIV) and numerically using the CFD simulations. These studies have showed the generation of vortices in the downstream wake and few works have considered comparatively the vortical patterns of different oscillation mechanisms. This paper does not aim to deeply investigate the wake structures of undulating and heaving mechanisms, but some important comparative and qualitative features at low and high Strouhal numbers are presented to recognize how significant discrepancies between the performance and the effects of these oscillation mechanisms might originate from their vortical patterns. However, to discover the wake phenomenon and its nature and also to distinguish in detail the qualitative and quantitative discrepancies of wake of these mechanisms, more deep researches are required to focus on the wake structure, the vortex-body or vortex-vortex interactions, and the evolution of the vortices from formation to decay. From low St vorticity contours (Figures 8 and 10) and values (Figures 12(a) and 12(c)), some similarities especially at trailing edge region can be observed in the vortical patterns of heaving and undulating mechanisms. This is attributed to the weakness of the vortices. In contrast, there are many discrepancies between their vortical patterns at high St (Figures 9, 11, 12(b), and 12(d)). Furthermore, Figures 8  and 10 show that, in the vicinity of trailing edge, distinct vortices cannot be observed at low St. On the other hand, these distinct vortices at low St are formed in a specific time instance and also at especial distance from the trailing edge. However, the delay on the formation of distinct vortices in the heaving foil is greater than the undulating foil. It is probably due to the stronger vortices in the undulating foil's wake (as denoted by maximum and minimum values on the top of each plot). The generated jets behind the foils and hence the produced thrust from the pairs of counter rotating vortices are also weak at low St. Furthermore, the vortices disappear in a short length behind the trailing edge because of the effects of the viscosity of medium.
The wake pattern of both oscillatory mechanisms depends on the operating Strouhal number. At higher Strouhal numbers, stronger vortices are produced and they are shed into the downstream at the tip of the trailing edge. The produced vortices in the leading and trailing edges of the heaving foil surface are both strong at higher Strouhal numbers. From Figures 9 and 11, it could be inferred that the produced jets behind the foils from the pairs of counter rotating vortices are strong due to the strong vortices. The produced vortices at higher St disappear in further distances behind the trailing edge in comparison to the lower St. The most significant difference between the vortical patterns of the heaving and undulating foils is the production of leading edge vortices by heaving oscillations. This feature can considerably affect the power consumption and the efficiency of the oscillation mechanisms (Figures 6 and 7). It is worth to note that, in some applications, the vortical signature and trace of propulsion or stabilizing mechanisms are particularly significant. As mentioned, the positive values of indicate regions where the strength of rotation (rotation rate) dominates the strength of strain (strain rate) and hence, -isosurfaces could denote the vortex envelopes. For example, the criterions shown in Figure 12 confirm our findings about the vortical patterns at low and high Strouhal numbers. In addition, better quantitative comparisons could be achieved from the criterion. From the values indicated at top of Figure 12, the ratios of maximum values of undulating over heaving oscillations are approximately 8.79 and 1.14, at low and high St, respectively. In addition, the ratios of low over high St maximum values of for undulating and heaving oscillations are approximately 33.1 and 254.47, at low and high St, respectively. It might be suggested that the lower power consumption and higher efficiency behaviors of undulating oscillations at higher Strouhal numbers are directly related to its vortical structure and the strengths of the formed vortices as well as their arrangements. Furthermore, as mentioned by vorticity criterion, the pattern of undulating foil shows more regular pattern. The primary reason is the appearance of the leading edge vortices for the heaving foil and their interactions with each other and/or with downstream vortices.
Studies of the unsteady hydrodynamics of oscillating motions of flexible and rigid fins have suggested that a rich set of phenomena exist, depending on the nondimensional frequency of oscillations, the wavelength of the excitation, and the aspect ratio of the fin. In some cases, [54], Zhang et al. [26], and also this study, simple wake structures have been observed that bear a strong resemblance to the structure of coflowing jets and wakes. In other cases, such as the work of Moored et al. [55], Dewey et al. [56], and Borazjani and Sotiropoulos [57] bifurcating wakes are seen, and both cases appear to correspond to a peak in efficiency. Moored et al. [55] developed an alternative framework, based on hydrodynamic wake. They found that local optima in propulsive efficiency occur when the driving frequency of a flapping fin matches the hydrodynamic resonant frequency of the jet profile and there can be multiple wake resonant frequencies and modes corresponding to multiple peaks in efficiency resonance theory [55]. However, in our cases, due to the geometry and kinematic parameters, the single row of vortices [57] called 2S wake pattern [55] with a peak efficiency point for each mechanism is observed.

Conclusions
In the present study, the simulations of two important natural propulsive mechanisms (pure heaving and undulating) were carried out to elucidate and compare the features of thrust generation, power consumption, and efficiency as well as vortical patterns, by systematically varying the Strouhal and Reynolds numbers. This study helps to get insights into the behaviors of these mechanisms in producing the propulsive and stabilizing forces for engineering applications. The simulations were carried out in Re = 4000, 40000, and 400000, with the maximum trailing edge peak to peak amplitude of The time histories of the instantaneous inline and transverse hydrodynamic forces and moment, the time averaged inline force and consumed power, the efficiency, and the contours of vorticity and criterion were presented and discussed.
In general, the time dependent forces and moments show two peaks corresponding to the forward and backward strokes of the foil tip, in each cycle. The upstroke and downstroke peaks of the heaving and undulating mechanisms are similar at low St, but the difference between them considerably increases with increasing St and a greater value for heaving motions. Strouhal numbers of oscillations determine Comparative observational features of the vortical patterns of both undulating and heaving foils show that the wake pattern of each motion depends on the Strouhal number. At low Strouhal numbers, the downstream wakes in the closely near regions of the foils' trailing edges show more similarities to each other. Moreover, at low St, the distinct vortices are formed with a delay. At higher Strouhal numbers, stronger vortices are produced and they are shed into the downstream at the tip of the trailing edge for both undulating and heaving foils. The formation of the vortices is directly related to the fluctuations of the pressure coefficients. The jumps in pressure coefficients represent high pressure gradient regions in which the vortices shed from the surface of the foil. The most significant difference between the pressure coefficient graphs and vortical patterns of the heaving and undulating foils is appearance of leading edge and surrounding vortices in the wake of the heaving foil which could considerably affect the power consumption and efficiency as well as the trace of oscillations. In addition, the criterions confirm our findings about the vortical patterns.