Computational Analysis of Propulsion Performance of Modified Pitching Motion Airfoils in Laminar Flow

The thrust generation performance of airfoils with modified pitching motion was investigated by computational fluid dynamics (CFD) modeling two-dimensional laminar flow at Reynolds number of 10. The effect of shift distance of the pitch axis outside the chord line (R), reduced frequency (k), pitching amplitude (θ), pitching profile, and airfoil shape (airfoil thickness and camber) on the thrust generated and efficiency were studied. The results reveal that the increase in R and k leads to an enhancement in thrust generation and a decrease in propulsive efficiency. Besides, there exists an optimal range of θ for the maximum thrust and the increasing θ induces a rapid decrease in propulsive efficiency. Six adjustable parameters (K) were employed to realize various nonsinusoidal pitching profiles. An increase in K results in more thrust generated at the cost of decreased propulsive efficiency. The investigation of the airfoil shape effect reveals that there exists an optimal range of airfoil thickness for the best propulsion performance and that the vortex structure is strongly influenced by the airfoil thickness, while varying the camber or camber location of airfoil sections offers no benefit in thrust generation over symmetric airfoil sections.


Introduction
The oscillation of an airfoil is a common form of propulsion for many animals, such as flying animals (birds, insects) and swimming species (fish, aquatic mammals).Because of its significance in studying the fundamental mechanism of flight and design of microaerial vehicles (MAVs), flapping airfoils aerodynamics has drawn much attention in recent years.Much work has been done on oscillating airfoils to examine the effects of various elements on thrust generation and propulsive efficiency.Typical flapping motions such as plunging, pitching, and combined plunging and pitching motions are defined by Rozhdestvensky and Ryzhov [1].Detailed studies related to flapping airfoils with pure pitching motion can be found in Garrick [2], Koochesfahani [3], Sarkar and Venkatraman [4], and Xiao and Liao [5].The general conclusion is that for the pure pitching motion, positive thrust force can be generated only at a relatively high reduced frequency ( > 8), and the thrust generated by pure pitching motion is too low to be useful and practical.
Flapping airfoils with pure plunging motion have been extensively investigated both numerically and experimentally by Lewin and Haj-Hariri [6], Young and Lai [7], Heathcote et al. [8], and Ashraf et al. [9].The thrust generation of an airfoil undergoing pure plunging motion is limited because of the leading edge vortex at high reduced frequency [6].It seems that a combined plunging and pitching motion is more attractive and much research has been conducted on such motion [9][10][11][12][13][14][15][16][17].The general conclusion from these investigations is that the thrust and propulsive efficiency generated by the oscillating airfoil strongly depend on the kinematic parameters such as flapping frequency, amplitude, and phase difference between pure plunging and pitching motion.A peak propulsion efficiency is observed at a phase angle shift of 90 ∘ (pitching motion leading plunging motion by 90 ∘ ) and the propulsion efficiency drops rapidly as soon as the dynamic stall occurs.
According to the studies discussed above, the application of pure pitching and plunging motions is limited due to their low thrust generation from a practical point of view, while the combined plunging and pitching motion may be limited because of its mechanical complexity.Go and Hao [18] proposed a modified pitching motion (Figure 1) and experimentally studied its propulsion characteristics in a water tunnel, and the modified pitching motion was also used by Platzer et al. [17] for MAV design.The difference between this modified pitching motion and the traditional one is the location of the pitch axis.For the traditional study of pitching airfoil aerodynamics, the pitch axis position is restricted within the chord line, while for the modified pitching motion the pitch axis lies outside the airfoil chord line and locates in front of the leading edge of the airfoil, which results in a combined plunge, pitch, and back-and-forth motion.Compared with the combined plunging and pitching motion, the modified pitching motion has an advantage in terms of mechanical simplicity and 3D motion combination effect.However, the experimental study by Go and Hao [18] is very limited because only the motion parameters of reduced frequency and pitch axis are considered and these parameters are chosen in a very small range.Besides, no numerical simulations have been conducted on the modified pitching airfoil yet.
The objective of this work is to systematically evaluate and quantify the thrust and efficiency of airfoils undergoing the proposed modified pitching motion.This modified pitching motion has the potential to increase the thrust generation through the change of motion parameters (the shift of the pitch axis by a distance  outside of the chord line, reduced frequency , pitching amplitude , and pitching profile) and airfoil shape (airfoil thickness and camber).Besides, the present study also provides useful data to understand the functioning of Platzer's MAV [17].

Computational Details
2.1.Flow Solver.The CFD package CFX 13.0 with an unsteady incompressible and viscous flow solver was applied to simulate the unsteady flow field around the airfoil section.The CFX code with second-order upwind spatial discretization and a second-order-accurate backward implicit scheme for transient terms was used for all the simulations.The mass flow was evaluated such that a pressure-velocity coupling was achieved by the Rhie and Chow [19] algorithm.In the present study, all the cases of pitching airfoils were investigated at the Reynolds number of 10 4 and the flow field was assumed to be laminar.The same flow solver has also been presented and employed in our previous work (Lu et al. [20]).

Dynamic Mesh Strategy.
The mesh was composed of two subdomains with a sliding grid interface, and both of the two zones were meshed by a structured grid.A no-slip boundary condition was applied along the airfoil surface.The inner zone had a diameter of 12 and the far-field boundary was set to at least 35 from the airfoil such that its effect on the flow around the moving airfoil was negligible [21].The dynamic mesh technique was employed for inner zone to realize the pitching motion whereas the outer zone was fixed, as employed by Wang et al. [22] and Ashraf et al. [9].The dynamic mesh motion was controlled by a CFX expression language (CEL) subroutine developed and attached to the CFX solver.With CEL employed, the specification of the motion of nodes on boundary or subdomain regions of the mesh can be realized.The same mesh motion strategy has also been used in our previous work (Lu et al. [20]).Yang et al. [23] reported that the O-type mesh has higher computational efficiency than the C-type mesh but yields almost the same aerodynamic forces as those resulting from the latter.Thus an O topology mesh was used in the inner zone.The height of the first row of cells was set at a distance to the wall of 10 −5 corresponding to  + less than one.

Motion Description.
The pure sinusoidal pitching motion of airfoil section is given by where   is the mean angle of attack and  is the pitching frequency.Sarkar and Venkatraman [4] pointed out that the pitching amplitude, mean angle of attack, reduced frequency, and pitch axis position restricted within the chord line have significant effect on thrust generation.They reveal that at a given oscillation frequency the average thrust decreases for higher values of mean angle of attack of the airfoil.As a result,   = 0 ∘ is chosen in this paper.In the present study, the effects of pitching profile were studied with various nonsinusoidal pitching motions, which are governed by equations expressed as follows: ( By tuning an adjustable parameter , one can gradually change the designed pitching profile from a sawtooth wave ( = −1) to a square wave ( → ∞).As a result, in effect  can be viewed as a measure of how rapidly the pitching airfoil reverses direction.The sample waveform shapes for  = −1, −0.97, −0.85, 0, 1, 2, and 10 at  = 10 ∘ are shown in Figure 2.

Model Validation.
At first we performed a rigorous validation study to assure satisfactory independency of the force predictions with respect to both mesh and time discretizations for the present study.Three sets of mesh sizes were considered: a coarse mesh (48,000 cells and 280 nodes on the airfoil), a medium mesh (96,000 cells and 400 nodes on the airfoil), and a fine mesh (192,000 cells and 500 nodes on the airfoil).Similar observations were also made for the temporal resolution analysis.Three levels of time steps per cycle were also considered (using the medium mesh): 1200, 2400, and 4800.
Validation tests were performed for two different typical operating points in the parametric space of pitch axis distance , reduced frequency , and pitching amplitude .Both of these two cases involve dynamic stall and leading edge vortex shedding (LEVS), but at different , , and .Results for mesh and time refinements are shown in Table 1 with time averaged thrust coefficient, power coefficient, and efficiency considered.Figure 3 shows the time history of   at  = 0.5,  = 20 ∘ , and  = 3.As can be seen, for both cases, 2400 time steps per cycle yield time-accurate predictions for both mean and instantaneous values.It is also shown that the mediummesh resolution provides satisfactory accuracy in space.For both test cases, differences between medium and fine mesh results are negligible, with variations of less than 2% on cycle averaged values, and within 4% on instantaneous peak force coefficients.
The further validation against the existing literature has been presented in our previous work (Lu et al. [20]), in which the pure plunging motion of a NACA0012 foil at  0 = 0.175 and different reduced frequencies ( = 0-8,  = 2/ ∞ ) at Re = 2 × 10 4 is simulated.The computational results are in close agreement with the Navier-Stokes predictions of Young and Lai [16] and experimental results by Heathcote et al. [8].

Motion Parameters.
The motion parameters have a strong effect on flow dynamics through their direct effect on the airfoil inertia.To systematically evaluate and quantify how motion parameters affect the thrust performance of an airfoil undergoing the proposed modified pitching motion, a wide range of  (from 0.5 to 2.5),  (from 2 to 7), and  (from 5 ∘ to 30 ∘ ) were considered.Besides, the effect of various nonsinusoidal motions on the propulsion performance of the pitching airfoil was also investigated.In this section all the computations were performed on an oscillating NACA0012 airfoil.

Effect of Reduced Frequency and Pitch Axis Distance.
First we consider varying the reduced frequency  and pitch axis distance  for a given . Figure 4   motion ( = −0.33)and modified pitching motion are shown in Figure 4.For traditional pitching motion the flow changes from drag producing to thrust producing at about  = 5, and even in such condition the thrust generated is rather insignificant.However, the proposed modified pitching motion has great significance in thrust generation enhancement; that is, for a given (, ) the increasing  induces much more thrust generated.As expected  mean noticeably increases with  and  mean appears to be a quadratic function of , and this agrees with the study by Go and Hao [18].
Sarkar and Venkatraman [4] reported that the thrust generation by an oscillating airfoil mainly depends on the qualitative wake structure behind the trailing edge of the airfoil.When a pitching airfoil operates in a thrust producing regime, the wake structure behind the airfoil is opposite in sense to the Karman vortex pattern.The wake structure called "jet" gives rise to a net thrust on the airfoil.At a fixed , an increase in  and  results in larger plunge amplitude and maximum pitch rate, leading to stronger reverse Von Karman vortex and averaged jet in the wake (not present here for conciseness).Thus, the increasing  and  lead to thrust generation enhancement.Meanwhile, the increasing  and  create a larger maximum effective angle of attack and bring a negative consequence; that is, the larger  and higher  cause more severe leading edge vortex shedding, which tends to reduce propulsive efficiency (Figure 4) [16].
Figure 5 shows the time histories of   and   for one flapping cycle (with / = 0 corresponding to  = 0 ∘ and the airfoil pitching downwards) at  = 10 ∘ considering the effects of  and .The comparison between the curves for the cases of  = 0.5, , and 1.5 at  = 3 reveals that the instantaneous thrust coefficient at larger  is systematically larger than that at smaller .Besides, the increasing  induces a significant increase in maximum thrust coefficient and input power coefficient, and this can be attributed to the higher maximum pitching rate induced by the larger .Similar effect of the increasing  at fixed  on   and   can also be found in Figure 5.For all the cases shown in Figure 5, the maximum   occurs as the airfoil pitches towards the mean position and the minimum   occurs as the airfoil pitches away from the mean position.This can be explained by the airfoil motion and the corresponding flow physics around the airfoil.The formation of reverse Von Karman vortex in the wake, which causes the thrust generation, begins as the airfoil pitches towards the mean position (/ = 0.25-0.5 and 0.75-1), and the process is almost complete when the airfoil reaches the mean position.As the airfoil pitches away from the mean position, the reverse Von Karman vortex has been shed into the wake and moves downstream.Besides, the occurrence of leading edge vortex shedding contributes to drag as it acts on a section of the airfoil surface that is facing downstream and thus causes a backward suction, leading to low level of   .As a result, most of the thrust is generated as the airfoil pitches towards the mean position.

Effect of Pitching
Amplitude.The primary objective of this section is to study the effect of pitching amplitude  on the modified pitching airfoil propulsion performance.While at fixed  and , changing  affects crucial characteristics such as the maximum pitch rate and effective angle of attack profile experienced by the airfoil.Figure 6 shows the variation of  mean and  with  for different  at  = 3.The time histories of   and   for one flapping cycle at  =  and  = 3 for different  are shown in Figure 7. Figure 6 clearly shows that when a pitching airfoil operates in a thrust producing regime, there exists an optimal range of pitching amplitude for maximum thrust.At relatively small , an increase in  leads to more thrust generated, and the effect of increasing  is similar to that of increasing  and .However, as  reaches some degree the increase in  brings no improvement in thrust generation while significant increase in power consumption.The vorticity fields for  = 5 ∘ , 10 ∘ , and 15 ∘ at  =  and  = 3 as the airfoil pitching upward and passing though the mean position shown in Figure 8 reveal that the larger  results in more vigorous reverse Von Karman vortex street, and we believe that this is the reason that an increase in  can enhance the thrust performance.However, the increasing  induces larger scales of leading edge vortices, which lead to greater fluctuation of thrust coefficient and contribute to drag as it acts on a section of the airfoil surface that is facing downstream and causes a backward suction.Besides, when the airfoil is pitching up or pitching down after passing through the mean position, the larger  leads to higher angle of attack and pressure at the pressure surface, resulting in more drag generated.Therefore, drag producing regimes are observed between / = 0-0.25 and 0.5-0.75 at  = 25 ∘ (Figure 7), which cause the decrease in  mean .However, the increasing  always induces a significant increase in power consumption at each , resulting in a rapid decrease in .

Effect of Nonsinusoidal Motion.
This section aims to study the effect of nonsinusoidal pitching motion of a   NACA0012 airfoil on its propulsion performance.Computations were conducted with five values of  = −0.97,−0.85, 0, 1, and 2 at  =  and  = 10 ∘ .The effects of nonsinusoidal pitching motion on  mean and  are presented in Figure 9. Figure 10 shows the vorticity field during the upstroke for different  at / = 0.75 for the case of  = ,  = 10 ∘ , and  = 3.The variations of  mean and  mean show a similar trend; that is, both  mean and  mean increase with .The propulsion efficiency firstly presents a flat distribution at relatively low  and then steadily decreases with the increasing .It is evident that the nonsinusoidal motion noticeably influences the pitching airfoil propulsion performance.Both  mean and  mean increase with  at fixed  and , while it is the other way round in terms of .
According to (2), increasing  leads to a higher maximum pitch rate, which causes stronger reverse Von Karman vortex (Figure 10) and averaged jet in the wake.As a result, the increasing  leads to more thrust generated.Meanwhile, the increasing  also results in a larger maximum effective angle of attack and causes more severe leading edge vortex shedding (Figure 10), which tends to reduce propulsive efficiency [16].It is also observed that the growth of  mean with is not very noticeable at  < 0 at a fixed .However, while  > 0, a faster rate of  mean growth with  is found.This is because the maximum pitch rate grows more and more rapidly as  increases.
The effect of nonsinusoidal pitching motion on flow structure lies not only in the strength of leading edge vortex and the reverse Von Karman vortex street in the wake, but also in the formation time of the vortices.From Figure 10 we can find that the increasing  induces an earlier occurrence of LEV, meaning that the LEV occurs at a lower /, while the trend for trailing edge vortex (TEV) is the other way round.Besides, as  increases, a notable change in the motion of TEV is noticed.The vertical distance between the upper row of counterclockwise rotating vortices and the lower row of clockwise rotating vortices decreases with the increasing .
To examine these results in detail, the instantaneous   and   in one oscillation period at  = 10 ∘ and  =  for various nonsinusoidal motions are presented in Figure 11.As expected the pitching profile noticeably affects the variation of   and   .The increasing  induces a great increase in maximum thrust coefficient and input power coefficient, and this can be attributed to the higher maximum pitch rate induced by the larger .In contrast its effect on the minimum thrust coefficient and input power coefficient is quite limited.Thus  mean and  mean increase with at fixed  and .The larger  results in the slower pitching reversal; that is, the larger  induces the airfoil to oscillate at relatively high pitching amplitude for a longer time, during which the angular velocity and thrust coefficient remain at a low level (Figure 11).

Effect of Airfoil Thickness.
To gain some insight into the effect of airfoil thickness, simulations with five symmetric airfoils (NACA0006, NACA0012, NACA0018, NACA0024, and NACA0030) were carried out at  = 10 ∘ and  = .
The time averaged thrust coefficient  mean and propulsion efficiency  for all the tested airfoil thicknesses are presented in Figure 12.As can be seen, for each  there exists an optimal range of airfoil thickness for the best propulsion performance, and the maximum thrust and efficiency are produced by the NACA0018 section.Besides, the effect of airfoil thickness on the thrust generated is more pronounced at higher .Again, we find that for each airfoil  mean increases with , and  tends to decrease with the increasing .
Figure 13 shows time histories of   and   for one flapping cycle at  = ,  = 10 ∘ , and  = 3 for airfoils with different thicknesses.We can find that the NACA0018 section generates the gentlest change in thrust, and a further increase in thickness results in a systematical decrease in instantaneous thrust coefficient.Besides, the increasing airfoil thickness induces a significant lag in the maximum thrust coefficient, meaning that the maximum thrust coefficient is obtained at a higher /.Despite the apparent differences in the variation of thrust coefficient, the resulting power coefficient varies little between the five sections, as seen in Figure 13.The greatest difference in the input power coefficient between airfoils with different thicknesses occurs at about / = 0.5 and 1, when the airfoil passes through the mean position and gets the maximum pitch rate.
In order to understand the physics behind the change in thrust performance due to increased airfoil thickness, the vorticity fields at various instances during the upstroke (/ = 0.25-0.75)are displayed in Figure 14.Previous researchers have reported that the thrust generation by an oscillating airfoil is mainly dependent on the qualitative wake structure downstream of the trailing edge.A time averaged thrust force is generated as a reverse Von Karman vortex street occurs in the wake, as shown in Figure 14.The formation of the leading edge vortex (LEV) varies with the thickness of the airfoil section.The thin airfoil has an effectively sharper leading edge compared to the rounder leading edge of the thicker airfoils.Thus the separated flow area on the airfoil surface for NACA0006 is much larger compared with that for the NACA0018 and NACA0030 airfoils at the same instant.For NACA0006 airfoil, the large scale leading edge separation dominates pressure distribution of airfoil surface and results in the great fluctuation of thrust coefficient (Figure 13).For NACA0012 and NACA0018 airfoils, the frontal area available for the suction pressure to react upon is greater than that of NACA0006 airfoil and thus more thrust is generated by thicker airfoils, as reported by Lentink and Gerritsma [24].However, the vorticity field for NACA0030 at / = 0.4-0.6 shown in Figure 14 reveals that with a further increase in the thickness of the section, the LEV starts to form only beyond the maximum thickness point, which leads to lower thrust coefficient throughout the cycle (Figure 13).As observed in Figure 13, the power coefficient   is not much altered due to varying airfoil section thickness; thus, an optimally thick airfoil section would provide more thrust and propulsion efficiency.

Effect of Airfoil
Camber.The objective of this section is to study the camber effect on the pitching airfoil propulsion performance.Three cambered sections (NACA2612, NACA4312, and NACA4612) were chosen for the consideration of both camber and camber location.The  mean and  for the cambered sections are compared with those of the NACA0012 (i.e., same thickness of 0.12 but with no camber) in Figure 15 with sinusoidal pitching motion used at  =  and  = 10 ∘ .As can be seen, differences between the  mean for different airfoils at each  are quite small, as well as .Thus, in contrast to the effect of motion parameters and thickness, varying the camber or camber location of airfoil sections offers little or no benefit in thrust generation over symmetric airfoil sections.The time histories of   and   in one oscillation period at  = 10 ∘ and  =  for each airfoil are presented in Figure 16 and the corresponding pressure coefficient distributions of the airfoils at / = 0.2 and 0.8 are shown in Figure 17.It is observed that the airfoil camber noticeably affects the variation of   and   .The only obvious difference in the results of a pitching cambered airfoil to that of a symmetric airfoil is that the cambered airfoil produces an asymmetric thrust coefficient during a pitching cycle.This is because the increase in thickness of upper surface and decrease in thickness on lower surface of the airfoil caused by the camber affect the pressure distribution of airfoil surface (Figure 17).Thus more thrust is generated in the upstroke of the pitching motion and approximately the same amount less is generated in the downstroke.As a result, the overall time averaged thrust remains almost the same.

Summary and Conclusions
The investigation of the thrust generation potential of airfoils with modified pitching motion has been carried out using 2D laminar flow calculations for various NACA airfoils at Re = 10 4 .The effects of motion parameters (the shift of the pitch axis by a distance  outside of the chord line, reduced frequency , pitching amplitude , and pitching profile) and airfoil shape (airfoil thickness and camber) were studied.The thrust generated and efficiency of airfoils undergoing the proposed modified pitching motion were systematically evaluated and quantified.The results reveal that the increase in  and  leads to more thrust generated.However, the trend for propulsive efficiency  is the other way round.It is also observed that when a pitching airfoil operates in a thrust producing regime, there exists an optimal range of pitching amplitude for maximum thrust.However, the increasing  always induces a rapid decrease in .To realize various nonsinusoidal pitching motions of airfoil, six adjustable parameters  = −0.97,−0.85, 0, 1, 1.5, and 2 were employed, among which  = 0 corresponds to a pure sinusoidal motion.An increase in  results in a better thrust generation performance at the cost of a decrease in , especially for  > 0.
Then the effect of airfoil shape on modified pitching airfoil propulsion performance was investigated with varying airfoil thickness and camber considered.The thickness study was performed on 2D NACA symmetric airfoils with 6%, 12%, 18%, 24%, and 30% thick sections.It is found that for each  there exists an optimal range of airfoil thickness for the best propulsion performance, and the maximum thrust and efficiency are produced by the NACA0018 section.Besides, the vortex structure is strongly influenced by the airfoil thickness.The study of camber effect indicates that varying the camber or camber location of airfoil sections offers no benefit in thrust generation over symmetric airfoil sections.Plunge position of airfoil :

2 MathematicalFigure 1 :
Figure 1: (a) Traditional pitching motion; (b) modified pitching motion with pitch axis at distance  to airfoil leading edge.

Figure 2 :
Figure 2: Variation of instantaneous () profile in one period for different  at  = 10 ∘ .

Figure 4 :
Figure 4: Variation of mean thrust coefficient and propulsion efficiency with  at  = 10 ∘ .

Figure 5 :
Figure 5: Variation of thrust and power coefficients with time at  = 3 and  = 10 ∘ for different .

Figure 13 :Figure 14 :
Figure 13: Variation of thrust and power coefficients with time at  = ,  = 10 ∘ , and  = 3 for airfoils with different thicknesses.

Figure 15 :Figure 16 :
Figure 15: Variation of mean thrust coefficient and propulsion efficiency with reduced frequency at  =  and  = 10 ∘ for airfoils with different camber.