Numerical study of the propulsive performance of two-dimensional pitching foils at very high frequencies: scaling laws and turbulence e ﬀ ects

Purpose – This paper aims to analyze the propulsive performance of small-amplitude pitching foils at very high frequencies with double objectives: to ﬁ nd out scaling laws for the time-averaged thrust and propulsive ef ﬁ ciency at very high frequencies; and to characterize the Strouhal number above which the effect of turbulenceon the meanvalues cannot be neglected. Design/methodology/approach – The thrust force and propulsive ef ﬁ ciency of a pitching NACA0012 foil at high reduced frequencies ( k ) and a Reynolds number Re = 16 000 are analyzed using accurate numerical simulations, both assuming laminar ﬂ ow and using a transition turbulence model. The time-averaged results are validated with available experimental data for k up to about 12 (Strouhal number, St, up to 0.6). This study also compares the present numerical results with the predictions of theoretical models and existing numerical results. For a foil pitching about its quarter chord with amplitude a 0 = 8 o , the reduced frequency is varied here up to k = 30 (St up to 2), much higher than in any previous numerical or experimental work. Findings – For this pitch amplitude, turbulence effects are found negligible for St . 0.8, and affecting less than 10% to the time-averaged thrust coef ﬁ cient C T for larger St Linear potential theory fails for very large k , even for the small pitch amplitude considered, particularly for the power coef ﬁ cient, and therefore for the propulsive ef ﬁ ciency. It is found that C T (cid:1) St 2 for large St, in agreement with recent models, and the propulsive ef ﬁ ciency decays as 1/k, in disagreement with the linear potential theory. Originality/value – Pitching foils are as ef ﬁ cient and energy harvesting devices. Their performance at very high reduced frequencies has not been suf ﬁ ciently analyzed before. The provide accurate numerical simulations to discern when turbulence is relevant for the computation of the time-averaged thrust and ef ﬁ ciency and how their scaling with the reduced frequency is affected in relation to the laminar- ﬂ ow predictions. This is relevant because some small-amplitude theoretical models predict high propulsive ef ﬁ ciency of pitching foils at very high computations Picasso frequencies over certain ranges of the structural parameters, and only very accurate numerical simulations may decide on these predictions. the of ( (cid:7) is in this laminar and turbulent numerical simulations at a low-to-moderate Reynolds number ( Re 16000), and ﬁ xed ( and axis ( a – 1/2). of for small-amplitude pitching motions relation . for time-averaged coef ﬁ cient propulsive ef ﬁ ciency.


Introduction
The study of the aerodynamic performance of pitching foils is of great physical interest to understand the propulsion mechanisms of many aquatic and aerial animals and also of technological interest to design flapping propellers for bio-inspired aquatic and aerial selfpropelled robots, among some other applications such as wind turbines and rotorcrafts. Consequently, it has been the subject of numerous investigations by many researchers (see, e.g. the recent review by Wu et al., 2020). Some of these works have been dedicated to the interesting study of transient behaviour and dynamic stall associated to highamplitude pitching motion (Lee and Gerontakos, 2004;Lind and Jones, 2016;Geng et al., 2018). Here we focus on the no less interesting study of modeling the time-averaged thrust and propulsive efficiency of small-amplitude pitching motion, for which there also exist a large number of recent experimental and numerical studies (among them, Mackowski and Williamson, 2015;Floryan et al., 2017;Mackowski and Williamson, 2017;Moored and Quinn, 2018;Senturk and Smits, 2019;Ayancik et al., 2019;Thakor et al., 2020;Mahbub and Muhammad, 2020;Rezaei and Taha, 2021). In these previous works, however, the range of high reduced frequencies has not been sufficiently explored owing to the experimental and numerical difficulties. Recently, Fernandez-Feria and Sanmiguel-Rojas (2020) presented results from laminar-flow numerical simulations to find out how the pitch axis location affects to the aerodynamic performance of a pitching foil, for a wide range of pitching frequencies and small amplitude.
The characterization of the propulsive performance at high frequencies is of great interest for several reasons. First, because this high-frequency limit could decide about the validity of existing simple models, either theoretical or based on experimental and/or numerical data, which mainly differ at high reduced frequencies (Fernandez-Feria and Sanmiguel-Rojas, 2019). Second, because high-fidelity numerical simulations for smallamplitude in this limit could elucidate the reduced frequency (or the Strouhal number) above which turbulence effects may play a relevant role in the computation of the propulsive performance of two-dimensional pitching foils at the low-to-moderate Reynolds numbers of interest in bio-inspired propulsion. And third, and more relevant, because recent theoretical results have shown that high propulsive efficiency may be reached for quite large values of the reduced frequency when the foil pitches about some particular axis (Fernandez-Feria, 2017), so it would be of interest to dispose of reliable numerical methods and results to characterize the propulsive performance of pitching foils at very high reduced frequencies. This is the main aim of the present study, where such accurate numerical method is provided, and where the mean thrust and propulsive efficiency of pitching foils are characterized at very high reduced frequencies, much higher than in any previous published work to our knowledge, for a Reynolds number Re = 16000. We have selected this value because it is in the range of interest in bio-inspired propulsion (Wu, 2011), and because, according to the recent laminar flow numerical simulations by Senturk and Smits (2019), the average thrust and power are relatively insensitive to Reynolds numbers greater than about 16 000. In addition, these authors claim that the boundary layer on the pitching foil remains laminar for this, and even higher, Reynolds number. However, we found here that Scaling laws and turbulence effects turbulence in the boundary layer may become relevant to characterize the propulsion performance as the Strouhal number is increased for this Reynolds number, even for smallamplitude pitching motion when the flow never detaches the boundary layer. The problem and the numerical method are described in Sections 2 and 3, respectively. The numerical results are validated against available experimental and numerical data up to moderately high reduced frequencies in Section 4, together with a discussion on the accuracy of the time-averaged results as the frequency increases. Results in a wide range of reduced frequencies are presented and discussed in Section 5, including a discussion of the role played by turbulence in the computation of the thrust force and propulsive efficiency. Finally, a summary of the main results and the most relevant conclusions are given in the Section 6.

Formulation of the problem
We consider a two-dimensional foil of chord-length c immersed in a uniform current of speed U of an incompressible fluid of density r and viscosity m pitching about a point P located at a distance ca/2 from the mid-chord ( Figure 1), with kinematics: where v is the angular frequency and a 0 the pitch amplitude. The problem is governed by the following nondimensional parameters: a 0 , the pivot point location, a, the Reynolds number, Re = rUc/m , and the reduced frequency: (2) where f is the pitching frequency in Hz. We shall also use, instead of k, the Strouhal number based on the trailing-edge amplitude A = (1-a) csin a 0 : valid for a pitch axis located upstream of the mid chord (-1 # a # 0). In most of the numerical simulations reported in the present work we shall use a NACA(National Advisory Committee for Aeronautics) 0012 airfoil pitching at a quarter-chord length (a = -1/2) for Re = 16000 and a = 8°(except for some comparisons with previous experimental data, for which other values of Re and a 0 are used), focusing on very large values of k. As usual, the thrust and power coefficients are defined as: where T and P = -M(da/dt) are the thrust and input power per unit span, respectively, with M the pitching moment per unit span. The corresponding time-averaged quantities are as follows: with t = 1/f the cycle's period. Time-averaged coefficients over a number of cycles will be used in the numerical computations. The propulsive efficiency is defined as: Some results will be given in terms of the thrust and power scaled with the maximum oscillating speed fA instead of the stream speed U: These coefficients are related to equations (5)-(6) as:

Numerical method
The numerical method and grids employed in this work are similar to that used in Fernandez-Feria and Sanmiguel-Rojas (2020), with the main difference that here we also use turbulence models in addition to laminar flow. It is shown that these grids allow to reach higher reduced frequencies with precision and to capture correctly the turbulence scales at the boundary layers. For this reason, only the distinctive features of the numerical method are described below with some detail. Turbulence is modeled with the four equations transition g -Re u SST model, specially developed for transitional flows (Langtry et al., 2006;Langtry and Menter, 2009;Malan et al., 2009). It is a combination of the k -v SST model coupled with intermittency g and transition onset Reynolds number, Re u , which is the critical Reynolds number where the intermittency starts. The first two equations of the model are similar to those of the standard k -v SST, plus an equation for intermittency and a fourth one for transition momentum thickness Reynolds number. This method has been shown to capture accurately the turbulent shear stresses in unsteady flows at moderate Reynolds numbers in the rage considered in the present work (Gorji et al., 2014). More recently Tank et al. (2017) found that the transition SST model predicts correctly the lift and drag for a stationary NACA0012 Scaling laws and turbulence effects airfoil at a Reynolds number 5 Â 10 4 for angles of attack below 8°. Wellenberg et al. (2019) found a good agreement for the lift and drag values of a NACA0015 airfoil pitching at low frequencies for angles of attack lower than 8°and a Reynolds number 1.5 Â 10 6 . Thus, the transition turbulence model is used in the present study since it is capable of giving highly accurate predictions of the onset and the amount of flow separation under adverse pressure gradients by the inclusion of transport effects into the formulation of the eddy-viscosity, and it has the best precision-to-computational cost ratio. To speed up the convergence for the pitching airfoil, the simulations are initialized with an unsteady solution with the airfoil at rest for both the laminar flow and the SST turbulence model. We use a circular computational domain centered at the pivot point ( Figure 2). The radius of the domain is 25c, which is sufficiently large to impose the inlet/outlet boundary conditions without perturbing the flow close to the airfoil (Wellenberg et al., 2019). The airfoil oscillates without deforming the mesh by rotating the mesh as a rigid body inside a circular interface with a radius of 23c. A similar dynamic mesh technique with Ansys-Fluent for a problem of pitching and plunging plates in tandem was used by Ortega-Casanova and Fernandez-Feria (2019). However, for studies where the shape of the solid is changing, morphing methods combined with the dynamic mesh option of Ansys-Fluent work better (Abdessemed et al., 2018).
More details about the numerical implementation in Ansys-Fluent and about the mesh are given in Fernandez-Feria and Sanmiguel-Rojas (2020). The main features of the three meshes used in the sensitivity analysis, and some characteristics values of the convergence of the turbulent model used here, are given in Table 1. This table also includes the results for the time-averaged thrust coefficient computed with both laminar flow and the g -Re u SST turbulence model for St = 2, Re = 16000, a = -1/2 and a 0 = 8°, corresponding to a reduced frequency k^30, which is the most adverse case analyzed in the present study. Although the differences in C T between meshes #2 and #3 are below 2%, only mesh #3 guarantees that the y þ is below unity. This is a necessary condition to reproduce with accuracy the velocity profile in the viscous sub-layer in order to capture correctly the adverse pressure gradients. Therefore, all calculations reported in the following sections were performed with mesh #3 for both laminar and turbulent flows.

Comparison with previous numerical, experimental and theoretical results for low St
Once the mesh has been selected to obtain numerically accurate results, we first validate them by comparing with experimental results of C T by Mackowski and Williamson (2015). Particularly for a NACA0012 airfoil pitching about its quarter chord at Re = 12000 for a 0 = 2°and for reduced frequencies up to k = 12, corresponding to St^0.2 ( Figure 3). As the maximum value of the Strouhal number is quite small, only numerical results from laminar flow simulations with mesh #3 are plotted (the values of C T obtained with the transition g -Re u SST turbulence model practically coincide; see §5 below). In Figure 3 we also plot results from two approaches of the linear potential theory (Garrick, 1936;Fernandez-Feria, 2016), corrected with the quasi-static drag C D0 = 0.0373.

caling laws and turbulence effects
It is observed that the time-averaged results from the present numerical simulations and the experiments by Mackowski and Williamson practically coincide up to k % 6, and also with the results from the linear potential theory of Fernandez-Feria, while Garrick's results slightly overestimate both as k increases. For k Z 6 the numerical results are slightly larger than the experimental time-averaged results. This is partly due to the time-averaging of the experimental data, which may lose accuracy as the frequency increases. (Except otherwise specified, here we report results from simulations with 25 to 50 periods, depending on St, and computing the time-averaged quantities over the last ten cycles.) Actually, this was corroborated in Fernandez-Feria and Sanmiguel-Rojas (2020) by comparing numerical and experimental results for the thrust amplitude C max T À C min T =2, which does not need timeaveraging of the experimental (or numerical) data. The agreement shown in Figure 3 of that reference for Reynolds number 16 600 and three different values of the pitching amplitude is very good even for the highest values of the reduced frequency considered for each a 0 . In addition to this averaging effect, it must be noted that three-dimensional effects in the experimental flow, that are not solved in the present 2D simulations, become more relevant as the reduced frequency increases (Moriche et al., 2016).
We further compare the present numerical results with very recent numerical simulations by Senturk and Smits (2019) for a NACA0012 airfoil pitching about a = -0.5. In particular, we select time-averaged results for Re = 16000 and a pitching amplitude a 0 = 8° (Figure 4). (As shown by Senturk and Smits (2019), the results for Re = 16000 and Re = 32000 are practically identical.) Although the agreement between both sets of numerical results is excellent at low frequencies, up to St % 0.4, which corresponds to k % 6, the results start disagreeing for St Z 0.4, especially for C T . The main source for this disagreement is most likely to be the computation of the time-averaged quantities, which as St increases has to be performed after an increasing number of cycles has been computed, so that the temporal signal has evolved sufficiently towards a periodic state ( Figure 5). The time-averaged results of Senturk and Smits (2019) were performed over the final four cycles in computations of just ten cycles. For comparison sake, the present time-averaged results in Figure 4 are also performed over the last four cycles, but in computations of both 10 and 25 cycles, with the former ones approaching the results of Senturk and Smits. As aforementioned, in all the computations reported below the timeaveraged quantities are performed over the last ten cycles, in computations of 25-50 cycles, depending on St. This has a high computational cost: in parallel computing with 16 processors intel E5-2670 at 2.6 GHz, 64 GB of RAM memory and InfiniBand interconnection of 54 Gbits per second, the laminar flow cases take around three days and about twice for turbulent flow cases.
Also plotted in Figure 4 are the results from the linear potential theory of Theodorsen (1935) and Garrick (1936) for C P , and the approaches of Garrick (1936) and Fernandez-Feria (2016) for C T , with their corresponding propulsive efficiency h , both sharing the same C P . In both cases we have subtracted an offset drag C D0 to C T , computed by extrapolating to k ! 0 the numerical results for low k. In particular, we have selected the computations for k = 0.25, 0.5 and 0.75 [not shown in Figure 4(a)] to obtain C D0 À C T k ¼ 0 ð Þ ' 0:06, in close agreement with the value obtained by Senturk and Smits (2019) for these values of Re and a 0 (C D0 ¼ 0:0185 þ 4:6= ffiffiffiffiffi Re p ' 0:055). It should be noted that this quasi-static drag C D0 has nothing to do with the drag coefficient of a stationary foil for the same Reynolds number (i.e. the Blasius-type drag for high Re), which has also been computed to be 0.031. Now, for this higher pitch amplitude a 0 = 8°, the agreement with the linear potential theory is not so good as in Figure 3 for a 0 = 2°. But the most striking feature is the disagreement in C P , with the potential theory greatly underestimating the numerical results HFF 32,5 for moderate and large frequencies (though the pitch amplitude a 0 = 8°is no so large). As a consequence, Garrick's propulsive efficiency always tends to 1/2 as St increases, greatly overestimating the numerical results, in spite of the fact that for these relatively large frequencies and higher amplitudes C T by Garrick (1936) performs better than the impulsebased approach by Fernandez-Feria (2016). caling laws and turbulence effects

Results for high frequencies and discussion
Numerical results up to very high values of the reduced frequency are presented in this section, computed with both the laminar and turbulence models. The same case just considered, i.e. a 0 = 8°, Re = 16000 and a = -0.5, is selected, but now varying St up to 2, which according to equation (3) corresponds to k^30 for these values of a 0 and a. Figure 6 shows the numerical results for the time-averaged coefficients scaled with St 2 , as defined in equation (8), and the propulsive efficiency. Before discussing the different correlations also shown in this figure, it must be observed that the effect of turbulence in the time-averaged coefficients is noticeable only above a certain Strouhal number, being this effect more relevant in the thrust than in the power coefficient. The origin of these differences in the mean values can be appreciated in Figures 7 and 8, where the temporal evolution of C T and the vorticity field are plotted, respectively, for the cases St = 0.2 and St = 2. For St = 0.2, the instantaneous C T (t) computed with the laminar flow, and the corresponding vorticity field at a given instant of time, are practically indistinguishable from those computed with the turbulent model. For St = 2, however, the vorticity fields are quite different in the wake (Figure 8). This difference has a relatively small effect in the thrust coefficient [ Figure 7(b)], and even less in its mean value [ Figure 6(a)], but it is significantly larger than for St = 0.2.
The effect of turbulence on the forces that the fluid exerts on the airfoil is better appreciated in Figure 9, where the pressure and skin friction distributions on the airfoil obtained with the laminar and turbulent codes are compared for the same cases considered in Figure 8. For St = 0.2, the laminar and turbulent distributions are practically identical, differing for St = 2, especially for the pressure coefficient, which has a greater influence on the mean lift (not considered here) than on the time-averaged thrust. It is interesting to note that the generation of turbulence strongly depends on the Strouhal number St, in spite of the fact that the Reynolds number based on the chord-length and the upstream velocity is the same in all the simulations. For this value Re = 16000 (and even for Re = 32000), Senturk and Smits (2019) claim that the boundary layer remains laminar, but it is clear in Figure 10 that turbulence generation in the boundary layer depends on St for a fixed value of Re. In this Figure 10 we plot the turbulence kinetic energy for St = 0.2 (Figure 10(a)) and for St = 2 (Figures 10(b)-(c)) for the same instant of time plotted in Figures  8 and 9. In the first case the turbulence intensity is very weak, being only appreciable far downstream in the wake. Consequently, it has no practical effect on the aerodynamic force and moment on the foil because only the vorticity field in, or close to, the boundary layer, together with the vorticity shed at the trailing edge in its close vicinity, affect the forces and moment on a flapping foil (Martín-Alc antara et al., 2015). Note that very weak turbulence kinetic energy is also observed in the downstream part of the boundary layers in Figure 10(a), denoting a slight flow separation, but it is so weak that it has no effect at all on the thrust plotted in Figure 7(a). However, for St = 2 (Figures 10(b)-(c)), the turbulence intensity is much higher in the wake and, what is more relevant, in the boundary layer. Thus, turbulence generation starts almost from the leading edge, undoubtedly due to the high frequency oscillations of the leading edge pivoting about one quarter chord downstream, and affecting to the instantaneous thrust plotted in Figure 7(b). Nonetheless, the influence of this greater turbulence intensity in the time-averaged thrust and power remains relatively small.
To better characterize the Strouhal number above which turbulence affects the computed time-averaged values, Figure 11 plots the relative error in the values ofĈ T resulting from both models, Scaling laws and turbulence effects where the superscript T is for the turbulence model and the superscript L for laminar flow. This figure shows that the effect of turbulence on the time-averaged thrust triggers when St Z 0.8, but the difference between both formulations remains less than 10% forĈ T , and much smaller forĈ P .
We turn now to the correlations shown in Figure 6. Results from linear potential theory are not included in Figure 6 because, as discussed at the end of the previous section, Theodorsen's power coefficient becomes quite inaccurate as the frequency increases, which affects to the corresponding propulsive efficiency. Instead, we plot the correlations by Senturk and Smits (2019) forĈ T andĈ P , which these authors obtain from their laminar flow numerical simulations for this case and increasing Reynolds numbers, but considering only St # 0.6. These correlations, which are based on the previous scaling by Floryan et al. (2017), ð Þ , and the corresponding h ¼Ĉ T =Ĉ P , with b ¼ 2:9 À 24= ffiffiffiffiffi Re p ; C D0 ¼ 0:0185 þ 4:6= ffiffiffiffiffi Re p and e ¼ 7:1 À 24= ffiffiffiffiffi Re p . These correlations are specific for the pitching amplitude a 0 = 8°and the pitching axis at a = -0.5, working What seems clear from the present high frequency numerical results is thatĈ T tends to a constant value, e say, as St (or k for a given a 0 ) increases, which must depend on the pitch axis location a, and to a lesser degree on Re for high values of Re, but which is almost Relative difference of C T computed with the laminar flow and the turbulent model for the cases plotted in Figure 6 HFF 32,5 independent of the pitch amplitude a 0 since C T scales as St 2 $a 2 0 . For the present case with a = -1/2 (and Re = 16000), Figure 6(a) shows that e(a = -0.5) % 2, instead of the value b /2 % 1.35 of the correlation by Senturk and Smits (2019) for this Re. On the other hand,Ĉ P is linear with St for large k, but now the proportionality constant depends on a and a 0 , and to a lesser degree on Re (for large Re). Taking into account the scaling by Senturk and Smits (2019) and the definition (3) of St, it results that the quantity C * P Ĉ P 1 À a ð Þsina 0 =St is a function of only the pitch axis location a for large St, g(a) say, if Re is high enough. This quantity is plotted in Figure 6(d), from which g(a = -0.5) % 3.3. Figure 6(d) also shows that the correlation by Senturk and Smits (2019)  Finally, the propulsive efficiency for large frequencies (and large Re) can be approximated, from the present numerical simulations by: This correlation is also plotted in Figure 6(c) as a function of St for a = -1/2.

Summary and conclusions
The propulsive performance of a NACA0012 pitching foil in the limit of high reduced frequencies (k ) 1) is investigated in this paper from laminar and turbulent numerical simulations at a low-to-moderate Reynolds number (Re = 16000), and fixed pitch amplitude (a 0 = 8°) and pitch axis location (a = -1/2). The following summary of the main contributions of the present paper are thus exclusively for small-amplitude pitching motions and in relation to time-averaged quantities. The first conclusion is that turbulence intensity becomes appreciable, even in the boundary layer on the oscillating foil, as St increases for this Reynolds number. However, the effect of turbulence on the computed time-averaged thrust coefficient C T is less than 10% for the highest k considered, corresponding to a Strouhal number St = 2, being inappreciable for St . 0.8. Its effect is even less for the computed time-averaged power coefficient and propulsive efficiency. In relation to the scaling of the average values, it is found that C T = 2St 2 ð Þ ! e a ð Þ for k !1, or large St, in qualitative agreement with the model by Senturk and Smits (2019), based on a previous scaling by Floryan et al. (2017) for small St and fixed pitch axis location. It is found that C P 1 À a ð Þsina 0 = 2St 3 ð Þ ! g a ð Þ as k, or St, increases, also in accordance with the model by Senturk and Smits (2019), and in clear disagreement with the linear potential theory, that predicts C P $St 2 for large k. This last behavior is however correct for small St. Finally, the propulsive efficiency decays as h(a)/k for large k. The characterization of the functions of the pitching axis location e(a), g(a) and h(a) is left for future work.