A modified Leishman–Beddoes model for airfoil sections undergoing dynamic stall at low Reynolds numbers

Journal of Fluids and Structures 93 (2020) 102852

Contents lists available at ScienceDirect

Journal of Fluids and Structures journal homepage: www.elsevier.com/locate/jfs

A modified Leishman–Beddoes model for airfoil sections undergoing dynamic stall at low Reynolds numbers ∗

Johan Boutet a , , Grigorios Dimitriadis a , Xavier Amandolese b,c a

Aerospace and Mechanical Engineering Department, University of Liège, Quartier Polytech 1, Allée de la découverte, 9, 4000 Liège, Belgium b LadHyX, CNRS-Ecole Polytechnique, F-91128 Palaiseau, France c Structural Mechanics and Coupled Systems Laboratory, Conservatoire National des Arts et Métiers, 2 Rue Conté, 75003 Paris, France

article

info

Article history: Received 6 February 2019 Received in revised form 21 October 2019 Accepted 18 December 2019 Available online xxxx Keywords: Leishman–Beddoes Dynamic stall Wind tunnel testing

a b s t r a c t A modified Leishman–Beddoes (Leishman and Beddoes, 1989) model is proposed to predict the aerodynamic load responses of airfoils undergoing dynamic stall at low Reynolds, low Mach numbers and low to very high equivalent reduced pitch rates. The modifications include using Wagner function aerodynamics for the attached flow model, defining and using a time-varying reduced pitch rate, equating the time delays applied to the normal force, the separation point location and the vortex shedding and using Sheng and Galbraith’s dynamic stall onset criteria (Sheng et al., 2006). The steady and unsteady force and pitching moment measurements of three airfoils with distinct stall mechanisms (a flat plate, a NACA0012, and a NACA0018) are used to demonstrate the validity of this new model. Dynamic tests consist of oscillating the airfoils in pitch around the quarter chord with a mean angle A0 = 10◦ and with different prescribed reduced frequencies from k = 0.015, to k = 0.16, and amplitudes from, A = 5◦ to 20◦ , at Reynolds numbers of the order of Re ≃ 1.8 × 104 . The resulting equivalent reduced pitch rates range from r ′ = 0.001 to 0.06. The computations of the values of the different model parameters are demonstrated and it is observed that the values of the time delay parameters change continuously and smoothly with time instead of jumping discontinuously at discrete time instances as proposed by Leishman and Beddoes (1989). It is shown that the modified Leishman–Beddoes model is in good agreement with experiments for low to moderate equivalent reduced pitch rates (r ′ ≤ 0.04). It also significantly improves the dynamic load prediction in this range of pitch rates, in comparison to the original Leishman–Beddoes model. Nevertheless, as for the original model, the modified model presented here becomes inaccurate at the highest equivalent reduced pitch rates of r ′ > 0.05. © 2019 Published by Elsevier Ltd.

1. Introduction The current popularity of small unmanned aerial vehicles has created significant interest in steady and unsteady aerodynamics at low Reynolds numbers. In particular, the development of flapping flight micro-drones requires an improved understanding of the mechanism of dynamic stall at such Reynolds numbers. ∗ Corresponding author. E-mail address: [email protected] (J. Boutet). https://doi.org/10.1016/j.jfluidstructs.2019.102852 0889-9746/© 2019 Published by Elsevier Ltd.

2

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

The phenomenon of dynamic stall refers to the unsteady flow behavior experienced by an airfoil for which the effective angle of attack varies in time in the stall region. This is a matter of great concern for wings operating at high angles of attack, helicopter rotor blades and wind turbine blades. Many studies have been devoted to the dynamic stall process (see for example the works of McCroskey et al. (1982, 1976)) and various dynamic stall models have been proposed (a review can be found in Larsen et al. (2007) and Dimitriadis (2017)). Several key points of the dynamic stall mechanism can be highlighted. The first one is a delay of the flow separation process above the static stall angle of attack. Beyond that point, a concentrated vortex may form near the leading edge which produces an overshoot in lift and a nose down pitching moment when the leading edge vortex is swept towards the trailing edge. Furthermore, significant hysteresis in the dynamic forces and moment also occurs. The dynamic stall process and associated fluid loads strongly depend on the airfoil geometry, mean angle of attack, reduced frequency and amplitude of the pitch or pitch-plunge motion. The Reynolds number, Mach number and turbulence characteristics of the incoming flow are also of great importance. The modeling of dynamic stall is challenging because it involves detached and turbulent flow. Typical approaches range from semi-empirical, such as the Leishman–Beddoes (LB) (Leishman and Beddoes, 1989) or the ONERA models, Dat et al. (1979), McAlister et al. (1984) to Computational Fluid Dynamic (CFD) methods. Examples of CFD solutions for 2-D and 3-D dynamic stall can be found in Srinivassan et al. (1993), Ekaterinaris et al. (1995) or Spentzos et al. (2007) works. The disadvantage of such approaches is that they are very computationally expensive and are therefore unsuitable for design purposes. Empirical models are much faster but they must be well calibrated for each individual airfoil, Mach number and Reynolds number. The popular Leishman–Beddoes dynamic stall semi-empirical model was developed initially for helicopter rotor applications and tuned for Mach numbers higher than 0.3. Sheng et al. (2008a,b, 2006) showed that at lower Mach numbers the LB model required modifications in order to correspond better to experimental observations. In particular, the computation of dynamic stall onset and vortex onset must be modified at such Mach numbers. The modifications by Sheng et al. (2008a,b, 2006) concerned low Mach numbers but intermediate or even high Reynolds numbers, corresponding to wind turbine applications. The present paper shows that further modifications are needed for low and very low Reynolds numbers and proposes a new version of the LB model. This new model aims to rethink the steady and unsteady parameters of the LB model and their computation from the corresponding experimental data. In particular, the effects of the frequency and amplitude of vibration (reduced pitch rate) on the LB model’s unsteady parameters is explored. Furthermore, an extension of the analysis by Sheng et al. (2008a,b, 2006) is carried out on the relationship between the vortex onset angle and the reduced pitch rate. The validity of this modified LB model is assessed for three wing sections with distinct static stall mechanisms : a flat plate, a NACA 0012 airfoil, and a NACA 0018 airfoil at Reynolds numbers of the order of 104 . 2. Experimental setup The experiments were performed in a closed-loop low-speed wind tunnel at LadHyX, which has a rectangular testsection 0.26 m wide and 0.24 m high. Tests have been performed at a mean flow velocity U ≃ 7 m/s, for which the non-uniformity in the test section is less than 1% and the turbulence intensity is close to 1.2%. 2.1. Three wing models Three wing sections were chosen for the experiments, a flat plate, a NACA0012 and a NACA0018. For intermediate Reynolds numbers, Re ≃ 106 , those wing sections are supposed to exhibit three distinct stall behaviors (Gault, 1957), a smooth thin-airfoil stall for the flat plate, a sharp leading-edge stall for the NACA0012 and a moderate trailing-edge stall for the NACA0018. At the low Reynolds number range of the present study, Re ≃ 104 , those stall classification could be still relevant, even if both the NACA models are also affected by a laminar separation bubble process at low angles of attack. This will be discussed in Section 3. The flat plate model was made from a carbon fiber plate of span l = 130 mm, chord c = 35 mm and thickness t = 1.5 mm. This model featured a rectangular cross section with a thickness-to-chord ratio of 4.3%. The leading and trailing edges were not rounded off, so that they remained sharp. The NACA0012 and NACA0018 models of span l = 130 mm and chord c = 38 mm were produced by means of 3D printing. Hand sanding was used to smooth the residual roughness caused by the 3D printing process. The three wing models are shown in Fig. 1 and their dimensions are summarized in Table 1. Each model was mounted horizontally in the test section. A small distance from the wall was allowed at one end to limit boundary layer perturbations. An end plate was installed at the other extremity of the model to suppress three-dimensional flow effects. 2.2. Measurement setup The harmonic pitching motion was directly driven by a brushless motor (Maxon flat motor EC 60, 100 W). The motor was controlled by a digital position controller (EPOS2 24/5) using Proportional/Integral/Derivative (PID) control. The unsteady loads acting on the wings were measured by a six-axis force/torque sensor (ATI Nano43) mounted between the motor and the models. A 24-bit data acquisition system furnished by Muller-BBM was used to receive the analog

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

3

Fig. 1. Side view of the tested wings.

Fig. 2. NACA0012 in mounted configuration inside the wind tunnel. Table 1 Wings dimensions and Reynolds numbers. Wing

Chord [mm]

Span [mm]

Reynolds

Flat plate NACA0012 NACA0018

35 38 38

130 130 130

17500 18400 18400

transducer signals and convert them to force and moment using the calibration matrix provided by ATI. Fig. 2 shows the mounting of the motor, force/torque sensor and wing (from left to right respectively) in the wind tunnel. The motion imposed by the motor was sinusoidal of the form

α (t) = A0 + A sin 2π ft where α (t) is the instantaneous pitch angle, A0 the mean pitch value, A the pitch amplitude and f the frequency. A laser displacement sensor (Keyence LB-11W) was used to measure the instantaneous angle of attack simultaneously with the unsteady forces and moments. All data were recorded with a sampling rate of 1024 Hz and the acquisition time was equal to or higher than 100 times the period of the harmonic pitching motion.

4

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852 Table 2 Equivalent reduced pitch rate for the flat plate as a function of the frequency and amplitude of motion in degree. f [hz]

1

2.5

5

7.5

k [–]

0.015

0.038

0.075

0.11

0.15

0.0013 0.0026 0.0039 0.0052

0.0033 0.0066 0.099 0.0133

0.0065 0.0131 0.0196 0.0262

0.096 0.0192 0.0288 0.0384

0.0131 0.0262 0.0393 0.0524

A A A A

= 5◦ = 10◦ = 15◦ = 20◦

10

Table 3 Equivalent reduced pitch rate for the NACA 0012 and 0018 as a function of the frequency and amplitude of motion in degree. f [hz]

1

2.5

5

7.5

k [–]

0.017

0.043

0.086

0.13

0.17

0.0015 0.003 0.045 0.0059

0.0038 0.0075 0.0113 0.015

0.0075 0.015 0.0225 0.03

0.0113 0.0227 0.034 0.0454

0.0148 0.0297 0.0445 0.0593

A A A A

= 5◦ = 10◦ = 15◦ = 20◦

10

2.3. Data processing The loads of interest were the normal force, N, chordwise force, C and pitching moment, M, as well as their coefficient forms N C M cn = , cc = , cm = 1/2ρ U 2 S 1/2ρ U 2 S 1/2ρ U 2 Sc where ρ is the air density, U the wind tunnel airspeed and S = c × l the wing surface. An end plate was installed on the free wingtip to reduce three-dimensional flow effects such that it can be assumed that cn , cc and cm are 2D load coefficients. Direct unsteady load measurements on moving models are subjected to two types of error. The first source is the natural vibration of the wing and motor supports, which in this case were designed to be as light and as stiff as possible. The second source of error involves the combination of static weight and dynamic inertia effects. Dynamic tare subtraction and low-pass filtering were applied to the data during post-processing in order to remove these errors. Dynamic tare was systematically identified at wind-off conditions prior to each dynamic test for a given amplitude and frequency of motion. A brick-wall low-pass filter was used to filter the dynamic loops at frequencies above 40 Hz for the flat plate, NACA0012 and NACA0018. This cut-off frequency eliminated the first resonant frequency of the full setup, while it preserved up to the 7th harmonic of the unsteady aerodynamic load for frequencies up to 5 Hz, i.e. a reduced pulsation close 0.08 in our experiments (see Tables 2 and 3). Nevertheless, some of the harmonics were lost due to the filter at the higher motion frequencies, 7.5 and 10 Hz. The complete unsteady measurement and data reduction procedure is demonstrated for the normal force measurement in Fig. 3 for a pitching motion of amplitude A = 15◦ , frequency f = 2.5 Hz and wind tunnel velocity U ≃ 7 m/s. The reduced frequency is defined as k=

ωb

U where ω = 2π f and b = c /2 is the half-chord. For the example of Fig. 3, k = 0.043. It must be noted that the wind-off load measurements included added mass aerodynamic effects but their amplitude was negligible at such low reduced frequencies. Therefore, it was assumed that structural inertial loads constituted the only measurable load component at U = 0 m/s. The procedure was the following: 1. Simultaneous measurement of the pitching motion, normal force (Fig. 3(b)), chordwise force and pitching moment at wind-off conditions. 2. Low pass filtering of the wind-off normal force (Fig. 3(c)), chordwise force and pitching moment. 3. Calculation of the statistics of the wind-off normal force (Fig. 3(d)), chordwise force and moment over at least 50 consecutive cycles of motion. 4. Simultaneous measurement of the pitching motion (Fig. 3(a)), normal force (Fig. 3(e)), chordwise force and pitching moment at wind-on conditions. 5. Low pass filtering of the wind-on normal force (Fig. 3(f)), chordwise force and pitching moment. 6. Calculation of the statistics of the wind-on normal force (Fig. 3(g)), chordwise force and moment over at least 50 consecutive cycles of motion. 7. Subtraction of the wind-off (dynamic tare) forces and moment from the wind-on loads (Fig. 3(h)).

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

5

Fig. 3. Summary of the unsteady measurement process for a NACA0012 airfoil with a pitching motion of amplitude 15◦ , frequency f = 2.5 Hz, wind velocity U = 7.25 m/s and reduced frequency k = 0.043.

6

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

8. Calculation of the resulting dynamic loops for the normal (Fig. 3(i)), chordwise and moment coefficients. The statistics used to show the results are the median and the interquartile region defined as the area between the first and third quartiles. At each considered time point on the cycles :

• The first quartile is the value below which lie the first lower 25% of all measured points at this specific time instance of the cycle.

• The median is the value below which lie the first lower 50% of all measured points at this specific time instance of the cycle.

• The third quartile is the value below which lie the first lower 75% of all measured points at this specific time instance of the cycle. 3. Static stall experiments A first set of experiments was carried out in order to measure the aerodynamic forces acting on the wings at static conditions. The measurements concerned the normal and chordwise force coefficients and the pitching moment coefficient at angles of attack ranging from 0◦ to 45◦ . The Reynolds number was constant at Re ≃ 1.8 × 104 . For each fixed angle of attack, the average values of the load coefficients were calculated from 10 seconds of measurements acquired at a sampling frequency of 1024 Hz. Fig. 4 plots the variation of the normal force and pitching moment coefficients with angle of attack. Measured data are plotted along with the results obtained using semi-empirical reconstruction models that will be discussed in detail in Section 5.2 and the classical Leishman–Beddoes reconstruction model (Leishman and Beddoes, 1989) with Eqs. (21) and (23). As mentioned earlier, the three airfoil shapes were chosen because they have different static stall mechanisms at intermediate and high Reynolds numbers: smooth thin airfoil stall (Flat plate), sharp leading edge stall (NACA0012) and moderate trailing edge stall (NACA0018). The thin airfoil stall behavior is confirmed for the flat plate which exhibit a short linear region followed by a smooth stall behaviors, i.e. the lift smoothly moves away from the linear evolution, with neither local maximum nor decreasing region up to 40◦ . The NACA 0012 model exhibit a sharper stall behavior for which a local maximum in lift is noticeable close to α = 9◦ . As reported by (Broeren and Bragg, 2001), a leading edge stall mechanism is characterized by an abrupt and sudden drop in normal force coefficient; such a drop cannot be seen in Fig. 4(c) and in the absence of pressure measurements, it is not possible to confirm whether the NACA 0012 experienced a leading edge stall at the low Reynolds number tested here. The NACA 0018 exhibits a smoother stall behavior (compared with the NACA 0012), which also occurs at a rather low angle of attack, close to α = 5.5◦ . Again it is difficult to say whether the stall mechanism is different between the NACA 0018 and NACA 0012, but for both airfoils a laminar separation bubble process occurs at low angle of attack, and it might certainly affects the stall mechanisms. Indeed, one can notice that the normal force for the two NACA airfoils is not linear at low angles of attack. This phenomenon has also been observed by Lutz et al. (2001) for the NACA 0009 airfoil, and is attributed to a laminar separation bubble at low angles of attack. This nonlinearity is ignored in the present work, as will be discussed in Section 5.2.1. The classic reconstruction model proposed by Leishman and Beddoes fails to represent the normal force coefficient at higher angle of attack as can be seen on Figs. 4(a), 4(c) and 4(e) because it does not allow for s < 0.04, as seen in Eq. (21). Furthermore, the static moment predictions of the classical Leishman Beddoes model do not fit the experimental data as well as the cubic spline model proposed in this work. 4. Dynamic stall experiments 4.1. Prescribed pitch motion The frequencies chosen for the pitching oscillations were f = [1, 2.5, 5, 7.5, 10] Hz. The corresponding ranges of reduced frequency are k = [0.015, 0.15] for the flat plate and k = [0.016, 0.16] for the two NACA wings. The mean pitch angle was always set to A0 = 10◦ and the oscillation amplitude to A = [5, 10, 15, 20]◦ . As a consequence, the maximum 2D wind tunnel blockage coefficient, defined as c sin(30◦ )/Hv , where Hv is the height of the test section, was always less than 8.4%. Therefore, no blockage corrections were necessary throughout the present study. Higher frequency and amplitude combinations were not possible due to motor limitations. All tests were carried out at a mean flow velocity of U ≃ 7.5 m/s for the flat plate and U ≃ 7.25 m/s for the two NACA wings. For pitch ramp tests, whereby the pitch angle increases in time at constant velocity, the reduced pitch rate is defined as (Sheng et al., 2008a, 2006, 2008b) b (1) U where α˙ is the constant pitch velocity in radians per second. For sinusoidal motion, the pitch velocity oscillates so that the reduced pitch rate is a function of time, i.e. r = α˙

r(t) =

ωb U

A cos ωt = kA cos ωt

(2)

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

Fig. 4. Normal force cn and pitching moment cm coefficient variation with angle of attack from static tests.

7

8

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

Fig. 5. Aerodynamic loads for a NACA0012 oscillating with amplitude A = 10◦ , frequency f = 5 Hz and mean angle A0 = 10◦ and Reynolds number Re ≃ 1.8 × 104 .

Where k is the reduced frequency and A is the amplitude expressed in radians. Sheng et al. (2006) define an equivalent (or maximum) reduced pitch rate as r ′ = kA

(3)

The reduced frequency and equivalent reduced pitch rate values for the dynamic stall tests carried out in the context of the present work are tabulated in Table 2 for the flat plate and 3 for the two NACA wings. 4.2. Dynamic stall onset As described by McCroskey et al. (1976), dynamic stall involves a stall delay, a complex reattachment process and, in some cases, a Leading Edge Vortex (LEV), which is shed near the leading edge and travels downstream. In that case, the dynamic stall onset angle is defined as the angle at which the vortex is shed. If no LEV is involved, which is the case for the flat plate airfoil in this paper, the dynamic stall onset angle can be defined as the angle at which the lift versus instantaneous pitch angle curve starts to depart from linearity. Dynamic stall onset does not occur at the same instantaneous angle of attack as static stall, it is delayed due to several different phenomena (see for example McCroskey et al. (1976)). It is therefore important to determine the dynamic stall onset angle for all the different experiments presented here. Sheng et al. (2006) stated that the dynamic stall onset involving a Leading Edge Vortex can be identified from experimental plots of the aerodynamic loads against instantaneous angle of attack by looking for one or more of the following:

• A change of slope in the normal force coefficient; • A local maximum in the upstroke section of the chordwise force coefficient; • A sudden drop in the pitching moment coefficient by ∆cm = 0.05. Figs. 5(a)–5(c) show an example of the experimentally measured median and interquartile envelopes for the normal force, chordwise force and pitching moment coefficients respectively, plotted against instantaneous pitch angle. The airfoil is the NACA0012, oscillating with mean angle of A0 = 10◦ , amplitude A = 10◦ and frequency f = 5 Hz. On each of the figures, the dynamic stall onset angle is identified, according to the three criteria given above.

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

9

Fig. 6. Comparison between equivalent reduced pitch rate and reduced pitch rate with angle of dynamic stall onset for a NACA0012.

Sheng et al. (2006) showed that there is a degree of uncertainty in the measurement of the angle of dynamic stall onset. In their case, the dispersion was of the order of around one degree. Figs. 5(a)–5(c) show a greater dispersion for the onset angle (of the order of 3 degrees). Of the three criteria, the change of slope in cn is difficult to locate precisely and the drop ∆cm = 0.05 seems arbitrary. The maximum of the chordwise force coefficient, Fig. 5(b), is the most robust criterion for locating the dynamic stall onset. In fact, the same criterion proved to be the most robust for the NACA 0012 and 0018 airfoils, in all the conditions tested during the course of the present work; it will therefore be used for these two wings throughout the rest of the paper. The criterion used for the flat plate airfoil will be described in the next subsection. Sheng et al. (2008a, 2006) related the dynamic stall angle for airfoils at low Mach numbers to the constant reduced pitch rate for pitch ramp experiments and to the equivalent reduced pitch rate for pitch oscillation tests. They plotted graphs of the variation of the dynamic stall angle, αds , against r ′ (evaluated from Eq. (3)) and derived a bilinear relationship between these two quantities. The alternative proposed in the present work is to relate the dynamic stall angle to the instantaneous reduced pitch rate, as given by Eq. (2). The procedure is the following: 1. For a particular choice of airfoil, amplitude and frequency, plot cc (t) against α (t) and identify the point where cc (t) is a maximum. 2. Evaluate the time at which the maximum value of cc (t) occurs, tds , and the corresponding pitch angle αds = α (tds ), pitch velocity α˙ ds = α˙ (tds ) and reduced pitch rate r(tds ). 3. Plot αds against r(tds ). 4. Choose a test with different amplitude and frequency values and repeat from step 1.. Fig. 6 compares the results obtained from plotting αds against r ′ and r(tds ) for the NACA 0012 airfoil. At low dynamic stall angles the two plots are nearly identical but, for the highest stall angles, there is a difference of up to 0.005 between r ′ and r(tds ). The advantage of plotting αds against r(tds ) rather than r ′ is that the pitch rate is calculated from instantaneous values and therefore applies to any general type of motion, not just sinusoidal oscillations. This characteristic is particularly useful when using αds (r) to model the aeroelastic responses of airfoils at high angles of attack, which are generally not sinusoidal. During such a response, if the flow is initially attached, dynamic stall will occur when α (t) and r(t) are such that α (t) = αds (r(t)) and α (t)α˙ (t) > 0. As mentioned earlier, Sheng et al. (2006) derived a bilinear relationship between dynamic stall angle and reduced pitch rate, as shown in Fig. 7(a). The figure plots the dynamic stall angles obtained by Sheng and Galbraith for the NACA 0012 undergoing pitch ramp up tests at a Reynolds number of 1.5 × 106 , along with the data measured during the course of the present work for the same airfoil, against r. The relationship between dynamic stall onset angle and reduced pitch rate can indeed be represented as bilinear, with the break between the two linear regions occurring around r = 0.01. The dynamic stall angles identified at Re ≃ 1.8 × 104 are all much lower than the ones obtained by Sheng and Galbraith due to the much more laminar nature of the boundary layer at this low Reynolds number. The bilinear model used by Sheng and Galbraith is not the only possible representation of the αds (r) relationship. Fig. 7(b) shows that the same data can also be represented by a continuous function, in this case a square root of the form

αds (r) =

−bds +

√

b2ds − 4ads (cds − r) 2ads

(4)

where ads , bds , cds are coefficients to be evaluated by a curve fit and αds is the angle of dynamic stall onset in degrees. Again, the advantages of using a continuous model for αds (r) become evident when modeling general motions that are not pitch ramps or sinusoidal oscillations, in particular aeroelastic responses. The absence of a discontinuity makes the

10

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

Fig. 7. Comparison of two curve fitting methods for the dynamic stall onset angle applied to NACA0012 at high and low Reynolds.

Fig. 8. Dynamic stall onset angle variation with instantaneous reduced pitch rate for the NACA 0012 and 0018 airfoils. Table 4 Curve fit coefficients for each airfoils. ads bds cds

N0012

N0018

Flat plate

1.0053 × 10−4 0.0017 −0.0157

5.206 × 10−5 0.0013 −0.0095

−2.9713 × 10−5 0.0035 −0.0259

time integration of the equations of motion simpler and less prone to numerical instabilities. Fig. 8 plots the dynamic stall onset angle function αds (r) for the NACA0012 and NACA0018 airfoils, along with the square root curve fits. Table 4 gives the coefficient values ads , bds and cds used for the curve fit for each airfoil. 4.3. Flat plate dynamic stall onset criteria The chordwise force response, cc (t), for the flat plate did not exhibit a clear maximum, so that this criterion could not be used for determining the dynamic stall angle. A different criterion was used, based on observing the variation of the normal force coefficient variation with pitch angle, cn (α ), at both static and dynamic conditions. Fig. 9(a) shows that the static normal force curve can be seen as an asymptote of two straight lines, one with a high slope before stall and one with a low slope after stall. A static stall angle can therefore be defined as the angle at which these two straight lines intersect. Fig. 9(b) demonstrates that similar behavior occurs during pitch oscillations. On the upstroke, cn (α (t)) can be seen as an asymptote of two straight lines. Consequently, the dynamic stall angle is defined as the angle αds at which these two lines intersect. The variation of the stall angle αds with instantaneous reduced pitch rate r for the flat plate is shown in Fig. 10. The function αds (r) is curve-fitted using the square root representation of Eq. (4). The coefficient values are given in Table 4.

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

11

Fig. 9. Static and dynamic stall criterion for the flat plate.

Fig. 10. Stall angle variation with instantaneous reduced pitch rate for the flat plate.

Note that αds (r) curves upwards for the flat plate while it curves downwards for the NACA 0012 and 0018 airfoils (see Fig. 8). 5. Modified Leishman–Beddoes model The Leishman Beddoes model applies to rigid airfoils with pitch and/or plunge degrees of freedom undergoing forced motion, usually sinusoidal or constant pitch-up. It can represent stall delay, trailing edge separation and the shedding of leading edge vortices. Its inputs are the airfoil geometry such as the chord and the position of the pitch axis, flight condition and motion and its outputs are the unsteady aerodynamic loads as a function of time. The modified version of the model proposed here can be visualized as a block diagram (see Fig. 11) containing the following blocks:

• Attached flow module: This module computes the circulatory and impulsive normal loads acting on the airfoil, cnc and cnI respectively. It will take as input :

– The pitch and/or plunge motion of the rigid airfoil section; – The chord; – The position of the pitch axis.

• Detached flow module: This module represents the nonlinear effects of separation and stall delay on the aerodynamic load responses. It will use as input : – The experimental static stall characteristics for pitching moment and normal force. Flow separation is represented by the location of the separation point s and the position of the center of pressure g which are obtained from the static experimental data. Leishman and Beddoes (1989) used a preset function to fit the experimental result. In the current paper, the functions s and g are curve fitted with a series of cubic splines.

12

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

Fig. 11. Modified Leishman Beddoes overview.

The position of the center of pressure is also defined directly as a function of the angle of attack g(α ) instead of a function of the separation point g(s) as previously defined by Leishman and Beddoes (1989). • Dynamic stall module: This module creates a function αds (r(t)) that estimates the angle of attack at which vortex shedding will begin, if there is one. It also calculates a characteristic unsteady time lag Ta applied to the angle of attack and to the separation point U αe − α ′ α˙′ =

(5)

s˙′ =

(6)

b Ta U ss (αe ) − s′ b

Ta

where αe is the effective angle of attack due to the motion, to be defined in Section 5.1.1. Finally, this module also calculates the time taken by the vortex to shed over the airfoil, Tv l . It will take as input : – The airfoil experimental dynamic stall characteristics at different amplitudes and reduced frequencies. The angle of vortex shedding onset αds was defined by Sheng et al. (2006) to replace the Leishman and Beddoes (1989) criterion at low Reynolds numbers Sheng et al. (2006) curve fitted the result with a set of two linear function. In the present model, the function αds is fitted with the inverse of a second order polynomial. • Unsteady flow module: This final module computes the aerodynamic loads acting on the airfoil. 5.1. Attached flow module 5.1.1. Circulatory lift In the classical Leishman–Beddoes model, the circulatory loading associated with the unsteady attached flow is computed by means of indicial compressible aerodynamic response functions. As the flow is assumed incompressible in the present work, the step change in lift coefficient, for an airfoil undergoing a step change in downwash ∆w ≪ U, is expressed using the Wagner function, Φ (t), as follows : cnc (t) = a0 Φ (t)

∆w

(7)

U

where a0 is the lift curve slope of the airfoil, which is usually approximated by 2π for thin airfoils, while Φ (t) is Jones (1938) approximation of the Wagner function

Φ (t) = 1 − Ψ1 e−

ϵ1 U b

t

− Ψ2 e −

ϵ2 U b

t

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

13

Fig. 12. Rigid thin plate airfoil scheme.

with Ψ1 = 0.165, Ψ2 = 0.335, ϵ1 = 0.0455, ϵ2 = 0.3. Duhamel’s principle can be applied to Eq. (7) in order to express a continuous lift response as the time integral of infinitesimal step responses (Fung, 1993) cnc (t) = a0

(

w(0)

Φ (t) +

U

1 ∂w (τ )

t

∫

∂τ

U

0

Φ (t − τ ) dτ

) (8)

where w (t) is the downwash at time t and τ is an integration variable. The troublesome can be removed by applying integration by parts. cnc (t) = a0

(

w(t) U

Φ (0) −

1 ∂ Φ (t − τ )

t

∫ 0

∂τ

U

w (τ ) d τ

∂w (τ ) ∂τ

term inside the integral

) (9)

The downwash w (t) must now be computed as a function of the kinematics of the wing. In this work, the modified Leishman–Beddoes model is presented in a general formulation involving pitch and plunge degrees of freedom. When the model will be compared with experiments in Section 6.2, the plunge will be set to zero. Fig. 12 defines the plunge, h(t), and pitch, α (t), degrees of freedom, the chord, c, and the distance between the pitch center and the half-chord, xe . The downwash becomes

˙ + α˙ (t)d w(t) = U α (t) + h(t)

(10)

where d = (1/2 − a)b and a = xe /b. Note that this definition of w (t) assumes low pitch and plunge velocities and small angles of attack. This definition is retained for higher angles of attack, but the lift is modified by the semi-experimental result from the detached flow module, as explained in Section 5.2 and it is assumed that this modification also represents the effect of high angle of attack. Applying a high angle correction to the added mass terms from Eqs. (17) and (18) has a negligible effect. From the downwash, an effective angle of attack is defined as

αe =

w

U After combining Eqs. (9) and (10), the circulatory normal force coefficient becomes cnc (t) a0

( =

α+

h˙ U

t

∫

Ψ1 ϵ 1

+ 0 t

∫

Ψ2 ϵ 2

+ 0

+ U b U b

e e

α˙ d

)

U −ϵ1 U b

−ϵ2 U b

Φ (0) (

(t −τ )

(t −τ )

(

α (τ ) + α (τ ) +

˙ τ) h( U

˙ τ) h( U

+

α˙ (τ )d

+

α˙ (τ )d

)

U U

)

(11)

dτ dτ

(12)

Following Lee et al. (1999), the continuous unsteady circulatory lift coefficient can be written as follows.

˙ (t) cnc (t) = Cq˙ + Dq + Ez + χ Φ

(13)

z˙ = Wz + Fq

(14)

where

[

q = h(t)

[

α (t)

]T ]T

z = z1 (t) z2 (t) z3 (t) z4 (t) ] a0 [ ˙ (0) U Φ (0) + dΦ ˙ (0) D = Φ U a0 U [ Ψ ϵ 2 ) ( Ψ ϵ2 E = − 1b 1 − 2b 2 Ψ1 ϵ1 1 − ϵ1 db b

)] ( Ψ2 ϵ2 1 − ϵ2 db

14

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

C =

a0 U

χ =− W = −

Φ (0) 1

[

a0 U U b

]

(h(0) + dα (0))

diag (ϵ1 , ϵ2 , ϵ1 , ϵ2 )

[

1 0

F =

d

1 0

]T

0 1

0 1

while the pitching moment coefficient around the pitching axis is obtained from c (t) = cm

1

(

1

2

2

) + a cnc (t)

(15)

The full derivation can be found in Appendix A. There is also an additional moment coefficient associated with the circulatory flow, given by c cm (t) = − 1

1 2

(

πb

)

1

−a

2

U

α˙ = Bm q˙

(16)

As this term is circulatory, it will be affected by flow separation. 5.1.2. Impulsive loads The non-circulatory sectional lift and moment coefficients, also known as the added mass effect, can be computed from the non-circulatory terms derived by Theodorsen (1935) and presented by Fung (1993) cnI = I cm =

π b (¨ U2

πb

πb

h − a b α¨ +

)

[

2U 2

ah¨ −

(

α˙ = An q¨ + Bn q˙ ] ¨ bα¨ = Am q

U) 1

a2 +

8

(17) (18)

As these terms are non-circulatory, they are assumed to be unaffected by flow separation. One should note that the terms associated with the accelerations were removed with the dynamic tare as explained in Section 2.2. However, these terms are negligible at the reduced frequency used in the current paper. 5.2. Detached flow module This module aims to compute the non linear effect of flow separation on the normal force and pitching moment coefficients. To that end, empirical functions s(α ) and g(α ) are developed from steady experimental data. Furthermore, this module represents the stall delay phenomenon by applying a time delay Ta to the instantaneous effective angle of attack. 5.2.1. Computation of s(α ) The Leishman–Beddoes model uses Kirchhoff–Helmholtz theory (Leishman, 2006) to represent the normal force coefficient acting on a static airfoil in an airflow as cn = a0 (α − α0 )

(

1+

√

s(α )

)2 (19)

2

where a0 is the lift curve slope, α is the angle of attack, α0 is the zero lift angle and s is the position of the separation point on the airfoil, as defined in Fig. 13. Eq. (19) is a simple model to compute the nonlinear force and moment behavior due to a progressive trailing edge separation process. According to Leishman and Beddoes [1], it can also be used for other separation process, such as the leading-edge stall, by promoting some separation at the trailing edge to model the nonlinear behavior in the force and moment response. By inverting relation (19), an expression for the separation point s can be deduced from experimental steady normal load data. s(α ) =

( √ 2

cn a0 (α − α0 )

)2 −1

(20)

This expression gives the values of s at all the discrete angles of attack at which cn was measured but an analytic function s(α ) is required by the model. Leishman and Beddoes (1989) chose the relation

{ s=

1 − 0.3 exp {(α − α1 )/S1 } 0.04 + 0.66 exp {(−α + α1 )/S2 }

α ≤ α1 α > α1

(21)

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

15

Fig. 13. Flow separation scheme for a flat plate.

Fig. 14. Comparison between the two NACAs steady normal lift and the chosen lift curve slope at low angle of attack.

where coefficients S1 and S2 define static stall characteristics and α1 defines the break point corresponding to s = 0.7. These coefficients are chosen to fit the data obtained from Eq. (20). The function of Eq. (21) can represent data from a wide set of lift curves but is not general. As an example, Leishman– Beddoes method overestimate the lift at high angle of attack as shown on Fig. 4 because the separation cannot be lower than s = 0.04. Sheng et al. (2006), and Pierce (1996) used lookup tables to compute the separation point from the experimental results. In the current paper, the static separation function s is curve fitted by a set of cubic splines to obtain the function ss (α ), where the subscript s denotes steady conditions. This choice affords much greater flexibility for modeling s(α ) data from any lift curve in comparison to Eq. (21) or the lookup tables proposed by Sheng and Galbraith and Pierce. The computation of the cubic splines is detailed in the Appendix B. The value of the lift curve slope a0 is determined from the static experimental data for each airfoil. For the flat plate, it is simply the slope of the linear region of the normal lift at low angle of attack which yields a0 = 6.18. For the two NACA airfoils, the experimental lift curve slope is not constant at low angles of attack due to the laminar separation bubble phenomenon. In order to define a lift curve slope, a line was traced starting at the origin and being tangent to the experimental lift curve, as shown in Fig. 14. The slope of this line was a0 = 5.2 for the NACA 0012; the same slope was used for the NACA 0018. This choice for a0 also has the advantage that the separation point function s never exceeds the value of s = 1 (separation point located at the trailing edge). Fig. 15 demonstrates the proposed fitting process on the NACA 0012 static lift curve of Fig. 4(c). Eq. (20) is applied to this lift curve in order to obtain the data points (crosses) in Fig. 15(a). These points are then curve-fitted by cubic splines in order to obtain the solid line. The figure also plots the best fit obtained from Eq. (21) (line and circles), demonstrating that the cubic spline approach follows the data much better for all angles α > 5◦ . Function ss (α ) can then be substituted into Eq. (19) to reconstruct the normal force curve. Fig. 15(b) shows that the reconstruction obtained from the cubic spline fit is much more accurate than the one obtained from Eq. (21), especially for high angle of attack α > 25◦ . Furthermore, as cubic splines are analytic, the function ss (α ) is also analytic. It should be noted that both curve fits in Fig. 15 are inaccurate for angles α < 5◦ because the Kirchhoff–Helmholtz model only represents a ‘‘trailing edge’’ stall mechanism. Nevertheless, function s(α ) only aims to represent the separation process, so the impact of the laminar separation bubble near the leading edge has been ignored when computing the curve fit of ss (α ), see Appendix B. 5.2.2. Computation of g(α ) Under attached flow conditions and thin airfoil theory assumptions, it results that the aerodynamic center of a 2D wing section lies on the quarter chord, which means that the pitching moment around c /4 is constant with angle of attack. As the flow separates this is no longer the case and the aerodynamic center starts to move downstream, as shown

16

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

Fig. 15. Curve fitting function s(α ) for the NACA0012 wing using cubic splines.

Fig. 16. Position of the point of application of the normal force under separated flow conditions.

in Fig. 16. Leishman and Beddoes (1989) suggested the following empirical relation to represent the non-dimensional distance between the point of application of the lift and the quarter chord g(α ) =

0 cm − cm1/4 1 /4

cn

(22)

0 where cm1/4 is the measured moment coefficient at the quarter chord, cn is the measured normal force coefficient and cm 1/4 is the zero lift moment coefficient at the quarter chord. As the moment was measured at the quarter chord during the 0 course of the present experiments, the moment measured at low angles of attack was nearly constant and cm1/4 ≃ cm1/4 (0). 0 The values were cm = 0.0174, cm0 1/4 = 0.006 and cm0 1/4 = 0.006365 for the flat plate, the NACA0012 and the NACA0018, 1/4 respectively. As in the case of the separation point position s(α ), Eq. (22) gives g(α ) at the discrete angles of attack at which the normal force and pitching moment coefficients were measured from static tests. In order to obtain an analytical expression for g(α ), Leishman and Beddoes proposed to represent the data of Eq. (22) using the function

g(s) = K0 + K1 (1 − s) + K2 sin(π sm )

(23)

Parameters K0 , K1 , K2 and m are coefficients used to fit the data. As s is a function of α , g(s) in Eq. (23) is essentially g(α ). In the present work, the empirical relation g is computed directly as a function of the angle of attack and cubic splines, defined in Appendix B, are used in order to curve fit gs (α ) from experimental data, as shown in Fig. 17. Again, cubic splines can fit the data with better accuracy than Eq. (23), as cubic splines are more flexible. Furthermore gs (α ) can be easily extended to negative angles of attack, which is not as straightforward for g(s) from Eq. (23). In Fig. 17(a), gs (α ) is equal to zero at angles α < 5◦ and around α = 10◦ , denoting fully attached flow. In the interval 5◦ < α < 10◦ , gs (α ) decreases down to −0.03 and then increases again. This phenomenon is assumed to be due to a laminar separation bubble at this low Reynolds number; it is ignored in the computation of g(α ) because the model does not aim to represent the effects of separation bubbles. Once the function gs (α ) has been evaluated, the pitching moment coefficient around the pitch axis can be computed from cm =

( ( 1

1

2

2

) ) + a − gs (α ) cn + cm0 1/4

(24)

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

17

Fig. 17. Curve fitting function g(α ) for the NACA0012 wing using cubic splines.

Fig. 18 plots the functions ss (α ) and gs (α ) measured during the course of the present experiments for all the wings, along with the square root curve fits. These curve fits were used to calculate the reconstructed load coefficients of Fig. 4. As a conclusion to the detached flow module section, it can be noted that the positions of the separation point, s(α ), and of the center of pressure, g(α ), could also be treated as purely empirical parameters that do not differentiate between the physical phenomena of stall and laminar separation bubbles. In that case, it would be possible to allow for values less than 1 for s and other than 0 for g at low angles of attack. This choice was not made within the context of the present work in order to focus on the stall process. 5.3. Dynamic stall module This module models dynamic stall, stall delay and the effect of the shedding and motion of the leading edge vortex. In particular, the time instance at which the vortex is shed and the time it takes to clear the trailing edge are estimated from the load measurements obtained during the dynamic tests. Furthermore, the aerodynamic load overshoot due to the passage of the vortex is included in the model. Vortex shedding is modeled only for the NACA 0012 and 0018 airfoils. One of the clear signs of the shedding and passage of a leading edge vortex is a sudden increase in the slope of the cn (α (t)) curve on the upstroke, after dynamic stall has occurred. Such behavior was not observed in any of the flat plate load responses and therefore the vortex shedding module was not used for this airfoil. 5.3.1. Dynamic stall and stall delay The functions ss (α ) and gs (α ) only apply to static stall and so do the normal and pitching moment coefficients of Eqs. (19) and (24). During pitching and plunging motion, the positions of the separation point and center of pressure must be calculated as functions of the effective angle of attack of Eq. (11). Furthermore, as mentioned earlier, dynamic stall occurs at instantaneous pitch angles higher than those at which static stall occurs, due to the stall delay effect. The Leishman–Beddoes model represents stall delay by applying a time lag to the instantaneous separation point, such that s˙′ =

U ss (αe ) − s′ b

Ta

where s′ is the lagged separation point and Ta is the time lag. For the standard LB model, the position of the center of pressure g is computed as a function of the lagged separation point gs (s′ ). In the present work, the center of pressure displacement is computed as a function of the angle of attack instead. A lagged equivalent angle α ′ is introduced such that

α˙ ′ =

U αe − α ′ b

Ta

Hence, the position of the center of pressure is gs (α ′ ). Sheng et al. (2008b) defined the time lag Ta as the slope of the second straight line fitting the dynamic stall onset data, as shown in Fig. 19(a). As the present modeling of αds (r) is based on a square root curve fit, the definition of Ta is generalized as Ta (r) =

dαds dr

= √

1 b2ds

− 4ads (cds − r )

2π 180

(25)

18

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

Fig. 18. Determination of separation function s(α ) and force application parameter g(α ) from the static tests.

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

19

Fig. 19. Two methods used to compute time lag Ta for a NACA0012 at Reynolds ≃ 1.8 × 104 .

Fig. 20. Evolution of lag parameter Ta with reduced pitch rate.

with αds expressed in radians, so that Ta varies continuously with r, as it is possible to deduce from Fig. 19(b) for the NACA 0012 wing. Fig. 20 plots the variation of Ta with r for all the airfoils studied in the present work. It can be noted that the non dimensional stall delay Ta decreases with the reduced pitch rate for the two NACAs and increases for the flat plate. 5.3.2. Leading edge vortex A vortex is created at the leading edge and starts to shed over the airfoil when the following conditions are met

αe ≥ αds (r) α˙ e αe ≥ 0

(26) (27)

Hence, the shedding of the vortex occurs during the upstroke and is related to the onset of dynamic stall. Once the vortex starts to shed, a non dimensional time variable tv = (t − t0 )

U

(28) b is defined to keep track of the position of the vortex over the airfoil, where t is the current time and t0 is the last time at which the stall onset criterion was met. The vortex clears the trailing edge when tv = tv l , the total time it takes for the vortex to travel over the entire surface of the airfoil. A non-dimensional version of tv l is given by Tv l =

U b

tv l

5.3.3. Estimation of Tv l The parameter Tv l represents the non dimensional time taken by the vortex to shed over the airfoil. It can be estimated with these different methods :

20

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

Fig. 21. Determination of the time taken by the vortex to travel over the airfoil tv l , for a NACA0012 airfoil and Reynolds number Re ≃ 1.8 × 104 . Table 5 Vortex passage time curve fit coefficients for each airfoil. aT v l bT v l cT v l

Flat plate

N0012

N0018

– – –

0.166 0.0794 −3.1128 × 10−4

0.036 0.07 −1.2074 × 10−4

• Time-resolved particle image velocimetry (PIV) to identify the vortex and its motion; • Unsteady pressure data to identify the vortex pressure effects along the chord; • Infer the position of the vortex from the measured unsteady aerodynamic load responses. In the current paper, only the last method can be used as there is no PIV or pressure data available for the test cases. The time taken in seconds by the vortex to shed over the airfoil, tv l , can be obtained from two time instances:

• The dynamic stall onset time, computed in Section 4.2. • The time when the normal load is maximum, assumed to be the time when the vortex reaches the trailing edge of the airfoil. An example of these two time instances are shown in the normal force plot of Fig. 21(a) and in the α (t) plot of 21(b). Estimating tv l for all the test cases and using equation 5.3.2 to plot the corresponding Tv l against instantaneous reduced pitch frequency leads to Fig. 22(a). It can be seen that Tv l can be curve-fitted as a continuous hyperbolic function of r. Plotting 1/Tv l against r for the NACA0012 (Fig. 22(b)) results in a square-root behavior, so that a reasonable choice for the curve fit is 2aTvl

Tv l =

−bTvl +

√

b2T vl

− 4aTvl (cTvl − r)

Table 5 gives the values of the curve fit coefficients aTvl , bTvl and cTvl for the two NACA airfoils studied in this work. No coefficients are tabulated for the flat plate, since vortex shedding is not modeled for this airfoil. Looking at Fig. 22(a), Tv l increases significantly for very low reduced pitch rates. In their original paper (Leishman and Beddoes, 1989), Leishman and Beddoes simply stated that this parameter should be constant for a given Mach number and did not study its variation with reduced pitch rate. This phenomenon is perhaps linked to a shielding effect of the airfoil as discussed by Choudhry et al. (2014) in the case of ramp tests with the airfoil stopping at the end of the ramp motion. 5.4. Unsteady flow module This module calculates the final aerodynamic loads applied on the wing section. It assembles the loads due to the attached flow, the trailing edge separation and the vortex shedding.

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

21

Fig. 22. Variation of Tv l with r for the two NACA wings (left), 1/Tv l vs r for the NACA 0012 (right).

5.4.1. Loads associated with trailing edge separation As in the classical formulation of the LB model (Leishman and Beddoes, 1989), the normal force due to the unsteady f trailing edge separation process, cn , can be estimated by means of a modified formulation of the Kirchhoff–Helmholtz relation, using the circulatory normal force computed by the attached flow module with the lagged static flow separation location s′ , i.e.

( cnf

=

cnc

√ )2 s′

1+

(29)

2 f

The associated pitching moment coefficient, cm , is then obtained using the attached cnc , an expression for variation of the center of pressure with respect to the delayed angles of attack α ′ and the moment induced by the wake, so that f cm

[( ( =

1

1

2

2

)

)

+ a − g s (α ) ′

cnc

−

1

(

2

1 2

) −a

πb U

α˙

](

1+

√ )2

2

s′

+ cm0 1/4

(30)

5.4.2. Circulation change Leishman (1989) added a vortex-induced normal lift coefficient cnv to represent the effects of the vortex on the total lift. It is described as an accumulation of vorticity that is lost in the wake when vortex shedding occurs. This vortex lift is represented by the following differential equation c˙nv = c˙v −

U cnv

(31)

b Tv

where cv =

⎧ ⎨

cnc

⎩

0

−

f cn

=

cnc

( 1−

(

√

1+ s′ 2

)2 )

tv ≤ 2Tv l

(32)

tv > 2Tv l

As this assumes that the vortex lift cnv is linked to vortex shedding, the position of the vortex has an effect on the v corresponding moment coefficient cm , so that

[

v cm = −0.25 1 − cos

(

π tv

)]

cnv

Tv l

(33)

In the current work, the interpretation of the vortex-induced normal lift coefficient cnv is slightly modified by removing the condition that cv = 0 when tv > 2Tv l . It is not directly linked to vortex shedding effects but mostly represents the delay in the change of circulation around the airfoil for abrupt motion. This lagged change of circulation is estimated from the difference between the attached and separated lift, i.e.

⎛ cv = cnc − cnf = cnc ⎝1 −

(

1+

⎞ √ )2 s′ ⎠

2

(34)

22

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852 Table 6 Aerodynamic load overshoot coefficients for each airfoil. B1 B2

Flat plate

N0012

N0018

– –

0.015 0.075

0.01 0.1

This change of circulation is also allowed to decay exponentially with a non dimensional time constant Tv . By combining the circulation buildup and decay one gets an added normal force that can be computed using Eq. (31). This force represents unsteady effects on the circulation caused by high pitch rates; cnv decreases with pitch rate, becoming equal to zero at static conditions. In the present, work, the decay time of this circulation was chosen equal to the dynamic stall time delay, i.e. Tv = Ta . The normal force cnv is linked to the circulation around the airfoil and applied to the center of pressure of the airfoil. Therefore, the pitching moment coefficient due to the excess accumulation of circulation in the vicinity of the airfoil is given by v cm =

( ( 1

1

2

2

) ) + a − gs (α ′ ) cnv

(35)

5.4.3. Vortex shedding effects The aerodynamic load overshoot caused by the motion of the vortex over the airfoil can be computed following Sheng et al. (2008a) assumption that this overshoot is proportional to the difference between the lagged s(α ′ ) and non-lagged s(αe ) separation points, i.e.

( ) ∆cnv = B1 s′ − ss (αe ) Vx v

(36)

v

∆cm = −B2 ∆cn

(37)

where B1 and B2 are parameters controlling the amplitude of the overshoot and Vx is a modified shape function proposed in this paper.

⏐ ⏐

(

Vx = ⏐⏐sin

)⏐ π tv ⏐⏐ 2Tv l ⏐

(38)

The shape function is periodic in order to represent the possibility of multiple vortex shedding and its effects on the aerodynamic loads. The coefficients B1 and B2 are assumed to be independent of the pitch rate and estimated to best fit the available data; they are tabulated in Table 6 for the two NACA airfoils (the flat plate does not feature vortex shedding). 5.5. Complete equations Finally, the total unsteady aerodynamic loads can be calculated. First the state space variables are computed from the following state space equations z˙ = Wz + Fq V ss (αe ) − s′ s˙′ = b Ta V α − α′ e ′

α˙ = c˙v n

b

Ta V cnv

= c˙v −

b Tv

Then, the aerodynamic loads are defined by the relations cnI = I cm =

π b (¨

h − a b α¨ +

U2

πb

[

2U 2

)

ah¨ −

(

a2 +

πb

α˙

U) ] 1

α¨

8

˙ (t) cnc = Cq˙ + Dq + Ew + r Φ

( cnf

f cm

=

cnc

1+

s′

2

[( ( =

√ )2

1

1

2

2

)

)

+ a − gs (α ) ′

cnc

−

1 2

(

1 2

) −a

πb U

α˙

](

1+

√ )2

2

s′

+ cm0 1/4

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

23

Table 7 Original Leishman–Beddoes parameters for each airfoil. a0 S1 S2

α1◦

K0 K1 K2 0 cm 1/4 cn1 Tp Tf Tv Tv l

v

( ( 1

1

2

2

)

Flat plate

N0012

N0018

6.18 0.0705 0.0936 7.89 −0.1475 −0.0027 0 0.0174

5.2 0.0494 0.0635 8.5 0 −0.1125 0.0013 0.006

5.2 0.0494 0.0423 5.8 0 −0.0719 −0.0356 0.006

∞

0.5687 1.7 8.8 8.8 8.2

0.4585 1.7 8.2 4 12

1.7 8.2 4 12

)

+ a − gs (α ) cnv ( ) ∆cnv = B1 s′ − ss (αe ) Vx ∆cmv = −B2 ∆cnv cm =

′

Finally, the total loads are cn = cnI + cnf + cnv + ∆cnv I cm = cm + cmf + cmv + ∆cmv

v Note that, for the flat plate airfoil, ∆cnv = 0 and ∆cm = 0 since vortex shedding does not occur. However, cnv and cmv are not zero since they are here interpreted as terms linked to circulation change around the airfoil for abrupt motions and not to vortex shedding.

6. Validation of the modified Leishman–Beddoes model 6.1. Original Leishman–Beddoes The results of the modified Leishman–Beddoes model presented in the present work are compared with experimental measurements and predictions from the original Leishman–Beddoes model (Leishman and Beddoes, 1989). The code used to compute the LB model comes from Dimitriadis (2017). The parameters used for the original Leishman–Beddoes model for each airfoil are presented in Table 7. It should be noted that the parameter cn1 was set to ∞ for the flat plate in order to remove any vortex shedding in the Leishman–Beddoes model similarly to what was done for the modified LB model presented in this paper. 6.2. Comparison Given the large number of test cases, only three reduced pitch rates have been chosen for each airfoil in order to compare the predictions obtained from the modified Leishman–Beddoes model presented in this paper (plain line), the original Leishman–Beddoes model (Leishman and Beddoes, 1989) (dot dashed line) and the experimental results (dashed line). Figs. 23–25 respectively show the comparison of the normal force and pitching moment coefficient for the flat plate, the NACA0012 and the NACA0018. Note that the case f = 10 hz, A = 20◦ was not shown for the NACA0018 because the experimental results were unreliable. For the low reduced pitch rate there is good agreement between both models and the experimental results for the flat plate (see Figs. 23(a) and 23(b)). For the NACA0012 and NACA0018 sections (see Figs. 24(a), 24(b) and Figs. 25(a),25(a)), the modified Leishman–Beddoes model is in good agreement with the experiments and improves significantly the prediction of the normal force overshoot, nose down pitching moment and hysteresis effect, in comparison to the original LB Model. Qualitatively, the agreement between both models and the experimental results deteriorates with increasing pitch rate. Nevertheless the modified Leishman–Beddoes model is still acceptable for intermediate pitch rates (see Figs. 23(c), 23(d), 24(c), 24(d), 25(c) and 25(d)). In particular, the dynamic loop shapes for both normal force and pitching moment obtained from the modified model better match the experiments. Nevertheless both models underestimate the normal force

24

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

Fig. 23. Comparison between model and experiment for the flat plate.

overshoot (see Figs. 23(c), 24(c), 25(c)) and for both the NACA sections, the nose down pitching moment is overestimated with the original LB model and underestimated with the modified formulation (see Figs. 23(d), 24(d) and 25(d)). Both model predictions become less accurate at the highest pitch rate. For the flat plate section the modified Leishman– Beddoes model gives better results than the original one for the normal force dynamic loop (Fig. 23(e)), but the amount of hysteresis in the pitching moment is strongly overestimated. For the NACA sections both the models underestimate significantly the normal force overshoot (see Figs. 24(e) and 25(e)), overestimate the pitching moment dynamic loop and,

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

25

Fig. 24. Comparison between model and experiment for the NACA0012.

as for the intermediate pitch rate, the nose down pitching moment is overestimated with the original LB model and underestimated with the modified formulation (see Figs. 24(f) and 25(f)). Concerning the normal force overshoot for medium and high pitch rates, the modified model gives slightly better results than the original one. Nevertheless the model could be improved using B1 and B2 coefficients in Eqs. (36) and (37), dependent on the pitch rate.

26

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

Fig. 25. Comparison between model and experiment for the NACA0018.

Concerning the pitching moment, the main issue for both models is to better capture the hysteretic effect for moderate and high pitch rates. This could be achieved with a better representation of the variation of the center of pressure position during the dynamic stall process. A better prediction of the nose down pitching moment could also be achieved with a more accurate identification of the Tv l nondimensional parameter.

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

27

7. Conclusions This work deals with the phenomenon of dynamic stall at low Reynolds numbers and low to very high equivalent reduced pitch rates, presents experimental measurements and attempts to model the observed aerodynamic load responses using the Leishman–Beddoes model, in both its original form and a modified formulation developed specifically for this Reynolds number range. It is shown that the modified model is more accurate than the original at low and intermediate reduced pitch rates, although both models become increasingly inaccurate at the highest pitch rate values considered here. The main differences between the present modified model and the original Leishman–Beddoes or the one proposed by Sheng et al. (2006) are:

• Sheng and Galbraith assume that the variation of dynamic stall onset angle with reduced pitch rate is bilinear whereas, in the present work, this variation is represented by a continuous quadratic function.

• Consequently, the time delay Ta is a continuous function of reduced pitch rate in the present work, while it is constant in the model by Sheng and Galbraith.

• The time delays Ta and Tv are taken to be equal to each other in the present work, as they vary smoothly in time. In both the original Leishman–Beddoes model and in Sheng and Galbraith’s version Ta and Tv are distinct and either constant or discontinuous functions of time. • The positions of the separation point and center of pressure during static stall are curve fitted by cubic splines in the present work, while they are fitted by discontinuous exponential functions in the original model and a lookup table in the Sheng and Galbraith version. These modifications result in a more robust and more stable dynamic stall model, which represents more accurately the experimental data obtained in the wind tunnel at very low Reynolds numbers. In fact, the same improvements can be applied to higher Reynolds number ranges but such an application is beyond the scope of the present work. The main inaccuracies in the model predictions concern the underestimation of the normal force overshoot due to the leading edge vortex and the overestimation of the thickness of the pitching moment hysteresis loop at the highest pitch rates. The first of these problems can be overcome by tuning the values of the B1 and B2 parameters in Eqs. (36) and (37), such that the leading edge vortex effect can be represented more accurately. However, there is no current physical model that can be used to adjust the values of these two parameters as functions of the pitch rate. In the original and modified Leishman–Beddoes models, B1 and B2 are constant but, clearly, this approach cannot be used at the very wide reduced pitch range considered in the present work. Hence, any further improvements to the modified Leishman–Beddoes model must concentrate on developing a robust way of estimating B1 and B2 as functions of the pitch rate. The overestimation of the thickness of the pitching moment hysteresis loop at the highest pitch rate is independent of the occurrence of stall, i.e. the thickness is overestimated both when the flow is attached and when it is separated, by both versions of the Leishman–Beddoes model. It is believed that this inaccuracy is related to the added mass terms, whose magnitude is overpredicted at the highest pitch rate by both the Wagner and Leishman–Beddoes impulsive load estimators. Improvements to the modified model should then also concentrate on better formulations for the added mass terms at high pitch rates. As a final thought, it should be noted that the highest equivalent reduced pitch rate for which results are presented in this paper is much higher than the pitch rates usually used in the Leishman–Beddoes and Onera model literature. Here, the highest reduced pitch rate is of 0.0529 for the flat plate and 0.0585 for the NACA wings, while the vast majority of comparisons between load predictions and measurements in the literature are carried out at reduced pitch rates that rarely exceed 0.02. Therefore, the present modified model is quite satisfactory for the reduced pitch rate range that is usually employed in the literature. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments The authors would like to acknowledge the financial support of the European Union (ERC Starting Grant NoVib 307265) and the technical support provided by the Laboratoire de’Hydrodynamique de l’Ecole Polytechnique (LadHyX). Appendix A. Circulatory lift computation The circulatory normal force coefficients was expressed as follow in Section 5.1.1 cnc (t) a0

( =

α+

h˙ U

+

α˙ d U

)

Φ (0)

28

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852 t

∫ +

Ψ1 ϵ1

U

Ψ2 ϵ2

U

0 t

∫ + 0

b b

e e

−ϵ1 U b

−ϵ2 U b

(t −τ )

(t −τ )

( (

α (τ ) + α (τ ) +

˙ τ) h( U

˙ τ) h( U

+

α˙ (τ )d

+

α˙ (τ )d

)

U

)

U

dτ dτ

(A.1)

Following Lee et al. (1999), the following changes of variables are performed in order to eliminate the integrals from Eq. (12) t

∫

e

z1 (t) =

−ϵ1 U

(t −τ )

b

h(τ ) dτ

z3 (t , y) =

t

∫

e

−ϵ1 U b

(t −τ )

α (τ ) dτ

−ϵ2 U

(t −τ )

α (τ ) dτ

0

0 t

∫ z2 (t) =

e

−ϵ2 U

(t −τ )

b

h(τ ) dτ

z4 (t , y) =

(A.2)

t

∫

e

0

b

0

where the variables zi (t) are called aerodynamics states. Substituting these definitions in Eq. (12) results in cnc (t) a0

α˙ d

( ) h(0) α (0)d ˙ = α+ + Φ (0) − Φ (t) + U U U U ( ) h αd U U ˙ (0) + − Ψ1 ϵ12 2 z1 − Ψ2 ϵ22 2 z2 +Φ U U b b ( ) ( ) U ϵ1 d U ϵ2 d + Ψ1 ϵ1 1− z 3 + Ψ2 ϵ 2 1− z4 h˙

(

b

)

b

b

(A.3)

b

First order differential equations for the aerodynamic state variables can be obtained by applying Leibnitz’s integral rule, such that z˙1 (t) = h − z˙2 (t) = h −

ϵ1 U b ϵ2 U

z1 (t)

z˙3 (t) = α −

z2 (t)

z˙4 (t) = α −

ϵ1 U b ϵ2 U

z3 (t) (A.4) z4 (t)

b b Finally, the continuous unsteady circulatory lift coefficient can be written as follows

˙ (t) cnc (t) = Cq˙ + Dq + Ez + χ Φ

(A.5)

z˙ = Wz + Fq

(A.6)

where

]T α (t)

[

q = h(t)

]T

[

z = z1 (t) z2 (t) z3 (t) z4 (t) ] a0 [ ˙ (0) U Φ (0) + dΦ ˙ (0) D = Φ U a0 U [ Ψ ϵ 2 ( ) Ψ ϵ2 E = − 1b 1 − 2b 2 Ψ1 ϵ1 1 − ϵ1 db b [ ] a0 C = Φ (0) 1 d U a0 χ = − (h(0) + dα (0)) U U W = − diag (ϵ1 , ϵ2 , ϵ1 , ϵ2 ) b

[ F =

1 0

1 0

0 1

( ) Ψ2 ϵ2 1 − ϵ2 db

]

]T

0 1

while the pitching moment coefficient around the pitching axis is obtained from c cm (t) =

1 2

(

1 2

) + a cnc (t)

(A.7)

Appendix B. Cubic spline curve fit The functions s(α ) and g(α ) are curve fitted with a series of cubic splines defined by a chosen number of N nodes as follows

⎧ ⎨ i1,0 α 3 + j1,0 α 2 + k1,0 α + l1,0 υ (α ) = in,n−1 α 3 + jn,n−1 α 2 + kn,n−1 α + ln,n−1 ⎩ i 3 2 N ,N −1 α + jN ,N −1 α + kN ,N −1 α + lN ,N −1

for α ≤ α0 for αn−1 ≤ α ≤ αn for α ≥ αN

(B.1)

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

29

Fig. B.26. Curve fitting function s(α ) and g(α ) for the NACA0012 wing showing the cubic splines nodes.

dυ (α )

Each node is defined by its position and the slope of the function going through the node (υ (αn ), dα n ). These nodes are fixed in order to get a solution for the parameters in , jn , kn and ln . The positions of the nodes are chosen manually. Two nodes are placed at each extremity of the dataset to avoid unwanted extrapolation. Nodes are placed at angles of attack where the slope changes abruptly in order to prevent exaggerated oscillations of the spline. The slopes at each node are obtained by a least squares fit of the experimental data points (αi , fi ) using the following error criterion e=

E ∑

(υ (xi ) − υi )2

(B.2)

i=1

except for the slopes of the first and last nodes, which are set manually. For the computation of s(α ), the effects of the laminar bubble are ignored by removing the first few data points of s(α ) up to the first local maximum of ss (α ). At this point, we assume that ‘‘trailing edge’’ separation is the dominant effect. The first node is fixed manually at α0 = 0◦ . The airflow is expected to be fully attached at zero angle of attack and s(α ) should decrease as the angle of attack becomes positive or negative. Consequently, the first node is defined as

α0 = 0◦ , s(α0 ) = 1,

ds(α0 ) dα

=0

(B.3)

The complete node locations for s(α ) are shown in Fig. B.26(a) for the NACA0012. For the computation of g(α ), the effects of the laminar bubble are removed by ignoring the same data points as for the computation of s(α ). Again, the first node is set manually. The center of pressure is assumed to be located at the quarter chord at zero angle of attack and should slowly move toward the mid-cord as the angle of attack increases. The first node is then

α0 = 0, g(α0 ) = 0,

dg(α0 ) dα

=0

(B.4)

Fig. B.26(b) plots the node locations for the NACA0012 and for the g(α ) function.

References Broeren, A.P., Bragg, M.B., 2001. Unsteady stalling characteristics of thin airfoils at low Reynolds number. Prog. Astronaut. Aeronaut. 195, 191–213. Choudhry, A., Leknys, R., Arjomandi, M., Kelso, R., 2014. An insight into the dynamic stall lift characteristics. Exp. Therm Fluid Sci. 58, 188–208. Dat, R., Tran, C.T., Petot, D., 1979. Modèle phénoménologique de décrochage dynamique sur profil de pale d’hélicoptère. In: 16e Colloque D’aérodynamique Appliquée, Lille, France. Dimitriadis, G., 2017. Introduction to Nonlinear Aeroelasticity. John Wiley & Sons. Ekaterinaris, J., Srinivasan, G., McCroskey, W., 1995. Present Capabilities of Predicting Two-Dimensional Dynamic Stall. Tech. Rep. 552, AGARD. Fung, Y.C., 1993. An Introduction to the Theory of Aeroelasticity. Dover Publications, Inc., Mineola, New York. Gault, D.E., 1957. A Correlation of Low-Speed, Airfoil-Section Stalling Characteristics with Reynolds Number and Airfoil Geometry. Tech. Rep., NACA TN-3963. Jones, R.T., 1938. Operational Treatment of the Nonuniform-Lift Theory in Airplane Dynamics. Tech. Rep., NACA-TN-667. Larsen, J.W., Nielsen, S.R., Krenk, S., 2007. Dynamic stall model for wind turbine airfoils. J. Fluids Struct. 23 (7), 959–982. Lee, B., Price, S., Wong, Y., 1999. Nonlinear aeroelastic analysis of airfoils: Bifurcation and chaos. Prog. Aerosp. Sci. 35 (3), 205–334.

30

J. Boutet, G. Dimitriadis and X. Amandolese / Journal of Fluids and Structures 93 (2020) 102852

Leishman, J.G., 1989. State-space model for unsteady airfoil behavior and dynamic stall. In: 30th AIAA Structures, Structural Dynamics and Materials Conference, Mobile, AL. Leishman, J.G., 2006. Principles of Helicopter Aerodynamics. Cambridge university press. Leishman, J.G., Beddoes, T., 1989. A semi-empirical model for dynamic stall. J. Am. Helicopter Soc. 34 (3), 3–17. Lutz, T., Würz, W., Wagner, S., 2001. Numerical optimization and wind-tunnel testing of low reynolds number airfoils. In: Mueller, T. (Ed.), Fixed and Flapping Wing Aerodynamics for Micro Air Vehicle Applications, vol. 195. American Institute of Aeronautics and Astronautics, Inc., Reston, VA, pp. 169–188. McAlister, K.W., Lambert, O., Petot, D., 1984. Application of the ONERA Model of Dynamic Stall. Tech. Rep., NASA-TP-2399. McCroskey, W.J., Carr, L.W., McAlister, K.W., 1976. Dynamic stall experiments on oscillating airfoils. AIAA J. 14 (1), 57–63. McCroskey, W.J., McAlister, K.W., Carr, L.W., Pucci, S., 1982. An Experimental Study of Dynamic Stall on Advanced Airfoil Sections. Volume 1. Summary of the Experiment. Tech. Rep., NASA TM-84245. Pierce, K.G., 1996. Wind Turbine Load Prediction Using the Beddoes-Leishman Model for Unsteady Aerodynamics and Dynamic Stall (Master’s thesis). The University of Utah. Sheng, W., Galbraith, R.M., Coton, F., 2006. A new stall-onset criterion for low speed dynamic-stall. J. Solar Energy Eng. 128 (4), 461–471. Sheng, W., Galbraith, R., Coton, F., 2008a. A modified dynamic stall model for low Mach numbers. J. Solar Energy Eng. 130 (3), 031013. Sheng, W., Galbraith, R., Coton, F., 2008b. Prediction of dynamic stall onset for oscillatory low-speed airfoils. J. Fluids Eng. 130 (10), 101204. Spentzos, A., Barakos, G., Badcock, K., Richards, B., Coton, F., Galbraith, R., 2007. Computational fluid dynamics study of three-dimensional dynamic stall of various planform shapes. J. Aircr. 44 (4). Srinivassan, G.R., Ekaterinaris, J.A., McCroskey, W.J., 1993. Dynamic stall of an oscillating wing. Part 1: Evaluation of turbulence models. In: AIAA Applied Aerodynamics Conference, 11th, Monterey, CA, Aug. 9-11. Theodorsen, T., 1935. General Theory of Aerodynamic Instability and the Mechanism of Flutter. Tech. Rep., NACA-TR-496.