Effect of active oscillation of local surface on the performance of low Reynolds number airfoil

JID:AESCTE AID:105774 /FLA

[m5G; v1.282; Prn:13/02/2020; 14:11] P.1 (1-12)

Aerospace Science and Technology ••• (••••) ••••••

1

Contents lists available at ScienceDirect

67 68

2 3

Aerospace Science and Technology

4

69 70 71

5

72

6

www.elsevier.com/locate/aescte

7

73

8

74

9

75

10

76

11 12 13

Effect of active oscillation of local surface on the performance of low Reynolds number airfoil

16 17

81

School of Aerospace Engineering, Beijing Institute of Technology, Beijing, 10081, PR China

83

82 84 85

19

a r t i c l e

i n f o

a b s t r a c t

86 87

21 22 23 24 25 26 27 28 29 30 31 32

79

Juanmian Lei ∗ , Jiawei Zhang, Jianping Niu

18 20

78 80

14 15

77

Article history: Received 16 August 2019 Received in revised form 19 December 2019 Accepted 8 February 2020 Available online xxxx Communicated by Qiulin Qu Keywords: Low Reynolds number Laminar separation Flow control Active surface oscillation Lift enhancement Drag reduction

In low Reynolds number conditions, laminar boundary layer separation and the consequent stall are usual, which are harmful to the aerodynamic characteristics of an airfoil. Energy addition to the flow by active oscillation of local surface can effectively restrain the flow separation and improve the performance. In the present paper, flow control based on active oscillation of local surface on an E387 airfoil is studied with numerical methods. The Reynolds number is Re = 30000 and the location of the oscillation region is a part of the upper surface near the leading edge. The respective relationships between the flow control effect and excitation parameters such as oscillation amplitude and excitation frequency are obtained. The flow mechanisms such as the behavior of the separation vortex and the variation of the equivalent shape are analyzed in both time-averaged and unsteady aspects. Furthermore, flow control effect at different angles of attack is also investigated. The results show that laminar separation can be effectively restrained, and the lift enhancement and drag reduction are evident. © 2020 Elsevier Masson SAS. All rights reserved.

88 89 90 91 92 93 94 95 96 97 98

33

99

34

100

35

101

36 37

Introduction

38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58

In aerospace research field, Reynolds number in the range of 10 –106 is usually called low Reynolds number [1]. In this situation, most of the flow area is laminar, and the flow contains relatively small momentum, so it easily separates from the wall if there exists an adverse pressure gradient. Transition to turbulence usually occurs after the laminar boundary layer separation, and a laminar separation bubble (LSB) will form with the reattachment [1–3]. After laminar separation, the decrease of lift and increase of drag mean the occurrence of stall, and the nonlinear and hysteresis behaviors of the vortex shedding will result in complex aerodynamic characteristics of the flight vehicle [4]. Various kinds of flow control method have been put forward to restrain laminar separation and improve the aerodynamic characteristics of flight vehicles operating at low Reynolds number [5,6]. One of them is virtual shaping, such as blowing, suction and synthetic jet/zero net mass-flux jet [7–9]. Another is local active deformation, such as the vibrating ribbon or wall vibration [10,11]. 4

59 60 61 62 63 64 65 66

*

Corresponding author. E-mail address: [email protected] (J. Lei).

https://doi.org/10.1016/j.ast.2020.105774 1270-9638/© 2020 Elsevier Masson SAS. All rights reserved.

The plasma flow control is also a useful means [12]. However, these flow control technologies are mainly designed for conventional flight vehicles and may be not fully suitable at low Reynolds number. Recently, a flow control method based on active oscillation of the flexible skin/membrane gains growing attention. With the active oscillation of local surface, energy is added into the near-wall fluid so the it can remain more attached to the wall. Munday and Jacob [13,14] carried out a series of experimental researches, in which the upper surface of a NACA 4415 airfoil was set to be actively oscillating so that the camber was altered. At different Reynolds numbers and different angles of attack, locations of the separation point and ranges of the separation region were measured based on flow visualization techniques. The results showed a notable reduction in the range of separation region at certain oscillation amplitudes and reduced frequencies. The results were satisfactory but further research was still needed. Katam [15] studied the cases in Ref. [14] with numerical simulation approaches. The upper surface between 20% and 80% chord was set to be the oscillating zone, and the parameter settings of deformation were approximately consistent with the wind tunnel test case. The results were in qualitative agreements with those in Ref. [14], and this work was also one of the pioneers of numerical research in this field. In the research of Kang [16], the computational fluid dynamics and computational structure mechanics were coupled. With numerical simulation methods, the self-excited vibration and the aerodynamic characteristics of a local flexible wing at low Reynolds numbers were studied. The fluid-structure inter-

102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

AID:105774 /FLA

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

[m5G; v1.282; Prn:13/02/2020; 14:11] P.2 (1-12)

J. Lei et al. / Aerospace Science and Technology ••• (••••) ••••••

2

67

Nomenclature A CD CDf C Dp CL Cp Fc Fv Re T V∞ W c cf f s t t y+

amplitude of oscillation, dimensionless drag coefficient, dimensionless friction drag coefficient, dimensionless pressure drag coefficient, dimensionless lift coefficient, dimensionless pressure coefficient, dimensionless convective flux vector viscous flux vector Reynolds number, dimensionless period of the flow field, dimensionless freestream velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m/s conserved variable vector airfoil chord length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m skin friction coefficient, dimensionless frequency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Hz mesh first cell height, dimensionless physical time, dimensionless physical time step size, dimensionless dimensionless first cell height

68 69

Greek

70

α τ τw ρ ω

angle of attack . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . degree pseudo time, dimensionless wall shear stress . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . N/m2 density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . kg/m3 angular frequency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . rad/s

∞ sep reatt ref

80 81 82 83 84 85 86

29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

87 88 89 90 91 92

26 28

75

79

25 27

74

78

midpoint of the oscillation region ratio of a value of the flow control case to its baseline case dimensionless value; time-averaged value

–

73

77

left end point of the oscillation region right end point of the oscillation region predefined or baseline value far-field boundary separation reattachment reference

Superscript * ’

72

76

Subscript 1 2 0

71

action process was analyzed. The results showed that the local flexible wing had a notably larger lift-drag ratio compared with the rigid wing. The research took fluid-structure interaction into consideration and this made it more practical. Furthermore, Kang [17] carried out a series of parametric studies about the flow around a NACA 0012 airfoil with local flexible structure at low Reynolds numbers. The effects of different oscillation frequencies, amplitudes and equilibrium positions on the aerodynamic characteristics and the evolution of the flow field structure were analyzed. Di [18] conducted numerical computations about the flow around a NACA 0012 airfoil with a vibration diaphragm in fully stalled condition, and the results were presented at different angles of attack and Re = 1.2 × 105 . The optimum parameters of the control methods were found to achieve the highest lift-to-drag ratio. Jones [19] used the Direct Numerical Simulations (DNS) method to explain the flow mechanism in the periodic surface morphing of an airfoil, and the flow control could almost eliminate the separation and significantly improve the aerodynamic performance. With DNS method, more flow details could be obtained and the flow mechanism was clearer. Although there is some literature about the flow control technology based on active oscillation of local surface, more comprehensive analyses about the effect of various factors are still needed. In the present research, numerical simulations are conducted to study the effect of active oscillation of local surface on the laminar separation and the performance of a low speed airfoil Eppler E387 at a low Reynolds number. The respective relationships between the flow control effect and different excitation parameters such as the oscillation amplitude and frequency are obtained. The flow mechanisms such as the behavior of the separation vortex and the variation of equivalent shapes are analyzed in both timeaveraged and unsteady aspects. The flow control effect at different angles of attack is also investigated. In this paper, the introduction of the numerical method and the verification are provided in Section 1. Then the definition of the oscillation scheme is described in Section 2. Section 3 consists of the numerical results and the discussions, including the effect of excitation amplitudes, frequencies and the flow control effect at different angles of attack. The flow mechanisms are also analyzed. Some useful conclusions are given in Section 4.

1. Numerical simulation methods

93 94

1.1. Governing equations, discrete scheme and turbulence model

95 96

Direct numerical simulation (DNS), large eddy simulation (LES) and Reynolds averaged N-S equations (RANS) methods can be employed in the numerical simulation of low Reynolds number flow [4,20]. For a complex configuration in engineering, the RANS method is an affordable approach and usually has satisfying accuracy. In this paper, the transient low Reynolds number flow field around an airfoil is studied, and the two-dimensional unsteady RANS (URANS) equation in its integral form is adopted, solved with the unsteady dual-time step method:

∂ ∂τ



∂ Wd + ∂t





Wd −

97 98 99 100 101 102 103 104 105 106 107

(Fc − F v )dS = 0

(1)

In Eq. (1), the vector W contains the conserved variables, Fc and Fv are the convective flux terms and viscous flux terms, respectively. The τ and t are the pseudo time and physical time, respectively. The  denotes the volume of the mesh cell and S is the boundary of the cell. For low Reynolds number flow, the free stream velocity is fairly low, and the flow can be dealt with as incompressible. In the present research, the incompressible N-S equations are solved with the finite volume method, and the pressure-velocity coupling is solved with the PISO (Pressure-Implicit with Splitting of Operators) [21] scheme. The spatial gradient is discretized with the cell-based least square method; the pressure, momentum, energy and time are discretized with second-order upwind scheme; and the turbulence relevant variables are discretized with first-order upwind scheme. The freestream velocity and pressure are given on the inlet and outlet boundaries, separately. No slip and adiabatic condition is defined on the wall boundary. Nonlinear flow phenomena, such as the boundary layer separation and transition to turbulence may occur in the low Reynolds number flow. In the present study, the intermittency transition model based on the SST k-ω model is utilized. This model is a further development on the basis of the γ -Re θ model [22,23]. The former avoids solving the Re θ equation, so it is less expensive.

108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:105774 /FLA

[m5G; v1.282; Prn:13/02/2020; 14:11] P.3 (1-12)

J. Lei et al. / Aerospace Science and Technology ••• (••••) ••••••

3

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12 14 15 16 17 18 19 20

More importantly, this model is Galilean invariant, so it can be applied to surfaces that move relative to the coordinate system. The model is appropriate for the prediction of transition in the wallbounded flows. 1.2. Mesh deformation technique

21 22 23 24 25 26 27 28

In the simulation of the flow around an airfoil with local surface oscillation, the mesh should be deformed according to the shape change of the airfoil. Some widely used mesh deformation methods [24,25] are not so suitable. In consideration of the characteristics of the present problem, a mesh deformation method based on the interpolation of displacement and diffusion is adopted. The expression is Eq. (2):



29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51

y =

(x1 ≤ x ≤ x2 , | y | ≤ y 0 )

f (x, t ) y∞ − y y∞ − y0

f (x, t )

(x1 ≤ x ≤ x2 , y 0 < | y | ≤ y ∞ )

(2)

In Eq. (2),  y is the displacement of the mesh point in ydirection; f (x, t) defines the displacement of the oscillation area as a function of the x-coordinate along the airfoil chord and time t. The x1 and x2 are the x-coordinates of the left and right end points of the oscillation area, respectively. The y ∞ is the y-coordinate of the far-field boundary, and y 0 is a predefined y-coordinate. According to Eq. (2), the y-direction displacement of the wall boundary generates the periodic movement of the oscillation area. For any near-wall mesh point in fluid, whose y-coordinate is between − y 0 and y 0 , the displacement is set to be equal to the corresponding wall boundary mesh point. For the other mesh points, of which the absolute values of the y-coordinate are between y 0 and y ∞ , the displacements are obtained by interpolation with Eq. (2). With this mesh deformation method, the relative location between the near-wall mesh point and the corresponding wall boundary mesh point is almost constant, and the mesh quality can keep well. This straightforward technique requires a very small calculation cost, and the mesh can still be fully recovered after several deformation periods.

52 53

1.3. Numerical methods verification

54 55 56 57 58 59 60 61 62 63 64 65 66

78

Fig. 1. Overall and near-wall mesh around NACA 0012 airfoil.

13

Comparing with using a scale-resolving model (like LES, DES, DNS, etc.), numerical simulation based on the URANS model is not always reliable, especially when predicting the time-dependent three-dimensional vortex shedding flow structure that also includes transition to turbulence. Thus a verification must be done and the URANS 2D results and the scale-resolving model results must match well, including the time averaged and the main transient features. Then the URANS 2D model can be used for parametric investigation. Jones et al. [26] carried out direct numerical simulations about the low Reynolds number flow around a NACA 0012 airfoil. In the present research, this case is simulated with URANS 2D methods

79

first and the results are compared with the DNS results [26], as a part of the validation of the URANS 2D numerical methods adopted in this paper. In the research of Ref. [26], the Reynolds number based on the airfoil chord is Re = 50000, and the angle of attack is α = 5◦ . In this paper, two sets of C-topology structured meshes with different node distributions and cell densities are generated for the test, one with node number 733*154 (G1) and the other with 445*111 (G2). The distance between the outer boundary and the airfoil is about 50 times chord length. The first cell height is s = 1.5 × 10−4 c, and the dimensionless first cell height y + is 1. Sketches of the mesh are shown in Fig. 1. The physical time step size is t = 1.5 × 10−3 c/V ∞ . The timemarching calculation is from t = 0 to t = 36.5c/V ∞ . The residuals decrease by 3 orders of magnitude in each physical time step, and the flow field shows periodic behaviors. Regular and periodic variations can be found in the lift coefficient C L and drag coefficient C D curves, as shown in Fig. 2. After the solution reaches convergence, one period in Fig. 2 contains about 250 time steps. The flow field is time-averaged on 10 unit time (10c/V ∞ ). Comparisons of the results are shown in Table 1 and Fig. 3. As can be seen from Table 1, for the time-averaged lift coefficient, the difference between the results of DNS 2D and DNS 3D is fairly large, which is about 20%; while the URANS result is very close to the DNS 2D, and the difference is about 3%. For the time-averaged drag coefficient, the difference between the results of DNS 2D and DNS 3D is about 14%, and the difference between the URANS result and the DNS 2D is about 10%–14%. The URANS results of the time-averaged separation point are close to the DNS 2D, and the time-averaged reattachment point is close to the DNS 3D. The coarse mesh G2 has larger errors in the solution of the reattachment point and vortex shedding frequency. In Fig. 3(b), the skin friction coefficient c f is defined as:

τw cf = 1 ρ V2 2 ∞ ∞

80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115

(3)

In Eq. (3), τ w is the wall shear stress, ρ∞ and V ∞ are the free stream density and free stream velocity, relatively. The positive direction of τ w is defined as from the leading edge to the trailing edge along airfoil surface. Thus, negative c f indicates local reverse flow. As can be seen from Fig. 3, for wall pressure coefficient and skin friction coefficient, the results of the URANS are fairly consistent with DNS, and the main flow field features and their distributions can be well predicted. In general, the differences in the URANS 2D results between mesh G1 and G2 are slight. According to the above results, the URANS 2D method adopted in this paper is proper for the prediction of the main features in low Reynolds number flow. In the 1990s, Selig et al. [27] conducted a series of tests and researches about 34 types of low speed airfoil in the UIUC lowturbulence subsonic wind tunnel. As a high-performance airfoil de-

116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

AID:105774 /FLA

[m5G; v1.282; Prn:13/02/2020; 14:11] P.4 (1-12)

J. Lei et al. / Aerospace Science and Technology ••• (••••) ••••••

4

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80 81

15 16

Fig. 2. Time histories of aerodynamic coefficients of URANS 2D G1. (a): lift coefficient; (b): drag coefficient.

17 18 19

10−4 c, and the dimension-

Table 1 Comparisons among the DNS and URANS results.

20 21

Case

CL

CD

Separation xsep

Reattachment xreatt

Vortex shedding frequency f

DNS 2D [26] DNS 3D [26] URANS 2D G1 URANS 2D G2

0.499 0.621 0.485 0.483

0.0307 0.0358 0.0264 0.0276

0.151 0.099 0.177 0.161

0.582 0.607 0.605 0.631

3.37

22 23 24 25

− 2.62 1.99

26 27 28 29 30

Table 2 Mesh parameters. Mesh

Total cells (k)

Points along airfoil

Points along wake

Points from wall to far field

Coarse Medium Fine

50 107 159

150 215 260

100 155 195

100 145 175

31 32 33 34 35 36 37 38 39 40 41 42 43

signed by Richard Eppler, E387 is appropriate for thermal duration sailplanes and UAVs, so it is selected for study in this paper. The relative thickness of this airfoil is 9.06% and the camber is 3.90%. The chord of the test model is c = 0.3048 m, and the Reynolds number based on the chord length is Re = 30000. The turbulence intensity is less than 0.1% in the wind tunnel. (1) Mesh independence verification

44 45 46 47

To study the effect of different mesh on the results, three sets of C-topology structured mesh are generated for the test, as are listed in Table 2 and sketched in Fig. 4.

The first cell height is s = 2.6 × less first cell height is y + = 1. On the basis of the test condition, the reference time is defined as t ref = c/V ∞ . All the dimensionless time-relative variables are based on t ref . In the mesh independence verification, the dimensionless time step is t = 0.002. As the angle of attack in a practical cruising flight is usually small, the angle of attack α = 4◦ is selected in the study, and the time history of unsteady lift coefficient is examined. As can be seen from Fig. 5, the results of the medium mesh and fine mesh agree well, and obvious periodic regularities can be noticed. While the result of the coarse mesh is quite different from those of the formers. Taking accuracy and computational cost into account, the medium mesh is selected for the following simulations. (2) Time step independence verification The dimensionless time steps t = 0.015, 0.002 and 0.0003 are selected to study the effect on the results. As can be seen from Fig. 6, the results of time steps t = 0.002 and 0.0003 agree well, while the result of t = 0.015 appears to be quite different. Taking accuracy and computational cost into account, time step t = 0.002 is selected for the following simulation. (3) Verification with the wind tunnel test results With the verified medium mesh and time step t = 0.002, numerical simulations of the transient flow field around an E387 airfoil are carried out at Reynolds number Re = 30000 and several angles of attack. As only the lift coefficient data is provided in

82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113

48

114

49

115

50

116

51

117

52

118

53

119

54

120

55

121

56

122

57

123

58

124

59

125

60

126

61

127

62

128

63

129 130

64 65 66

Fig. 3. Time-averaged aerodynamic coefficients. (a): wall pressure coefficient; (b): skin friction coefficient. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

131 132

JID:AESCTE AID:105774 /FLA

[m5G; v1.282; Prn:13/02/2020; 14:11] P.5 (1-12)

J. Lei et al. / Aerospace Science and Technology ••• (••••) ••••••

5

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

Fig. 4. Overall and near-wall mesh around E387 airfoil.

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27 28

93

Fig. 5. Unsteady lift coefficients of different meshes.

94

29

95

30

96

31

97

32

98

33

99

Fig. 7. Comparison of the C L results.

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107 108

42 43

Fig. 6. Unsteady lift coefficients of different time step sizes.

109 110

44 45 46 47 48 49 50 51

Ref. [27], the time-averaged lift coefficients are compared with the wind tunnel test results, as are shown in Fig. 7. As can be seen from Fig. 7, in the range of the angle of attack studied, the lift coefficient result of the numerical simulation agrees well with that of the experiment. From the verifications shown above, the numerical methods, mesh and time step size are appropriate for the numerical simulations in the present research.

52 53

2. Definition of the oscillation scheme

54 55 56 57 58 59 60 61 62 63 64 65 66

In some previous studies [13,14], the oscillation area is defined as from 20% to 80% chord, which covers the flow separation region, and the flow control is straightforward and effective. However, in engineering practice, some high-altitude long endurance airships or UAVs [28] are usually equipped with solar panels on the upper surface of the wing, so as to acquire continuous power for longer endurance. In consideration of this, also for the sake of structural strength, the oscillation region should be small and keep clear of the main part of the airfoil. In the present work, the oscillation region is defined as the segment from the leading edge to 10% chord on the upper surface of the airfoil. As the oscillation region is near the leading edge, its

111

Fig. 8. Coordinate system and deformation scheme.

112 113

disturbance could naturally travel downstream and affect the development of the flow around the airfoil. As are demonstrated in Fig. 8, the origin O of the coordinate system is fixed at the leading edge of the airfoil, with the x-axis pointing rearward along the chord, the y-axis pointing upward and perpendicular to the x-axis. Eq. (4) is the expression of the deformation:

114 115 116 117 118 119 120 121

 y = f (x, t ) = A x¯ sin(ωt )

(4)

122 123

In Eq. (4),  y is the magnitude of deformation in the direction perpendicular to the chord, and A is the amplitude of the oscillation. The x¯ is the dimensionless x coordinate, relating to the distribution of the deformation along chord. The ω is the angular frequency of oscillation, with the relationship ω = 2π f , and t is time. The distribution of deformation along chord is defined as

 x¯ = 1 −

2(x∗ − x) x2 − x1

124 125 126 127 128 129 130

2

131

(5)

132

JID:AESCTE

AID:105774 /FLA

[m5G; v1.282; Prn:13/02/2020; 14:11] P.6 (1-12)

J. Lei et al. / Aerospace Science and Technology ••• (••••) ••••••

6

1 2 3 4 5 6 7

Table 3 Excitation amplitudes. Number

A (×10−3 c)

Number

A (×10−3 c)

1 2 3 4

0.1 0.2 0.5 1

5 6 7 8

2 5 10 15

8 9 10 11 12 13 14 15 16

In Eq. (5), the x1 , x2 and x∗ are the x coordinates of the left end point, right end point and the midpoint of the oscillation region, separately. In the present study, x1 = 0, x2 = 0.1c, and x∗ = 0.05c. The comparison among the shapes of the rigid airfoil and the deformed airfoils is shown in Fig. 8. The red and blue lines indicate the positive and negative maximum deformation positions, separately. The black line is the baseline shape of the airfoil.

17 18

3. Results and discussions

19 20

3.1. Effect of the excitation amplitude

21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

In the present section, 24 cases are studied, including 3 excitation frequencies and 8 excitation amplitudes. The excitation frequencies are f = 0.5 f 0 , f = f 0 and f = 1.5 f 0 , separately; and f 0 means the vortex shedding frequency of the rigid airfoil. The excitation amplitudes of the conditions are listed in Table 3. After convergence, the time statistics of flow field is made on 10 vortex shedding periods, in order to remove the unsteady effect and obtain the time-averaged flow field characteristics. The time-averaged lift coefficient C L and drag coefficient C D are then divided by those of the rigid airfoil, being denoted as C L and C D , separately. Then the flow control effect can be seen clearly. As can be seen in Fig. 9, comparing with the baseline condition, the lift coefficient has a significant increment by up to 49%, and drag coefficient decreases by up to 39%. The lift-drag ratio increases by 118%. When the excitation amplitude is smaller than

A = 2 × 10−3 c, the results vary significantly with the amplitude. When A = 5 × 10−3 c, the flow control effects are fairly satisfying in the conditions concerned. The C L − t and C D − t curves at different excitation amplitudes when f = 0.5 f 0 , f 0 and 1.5 f 0 are shown in Figs. 10, 11 and 12, respectively. The amplitudes of flow parameter fluctuation become obviously larger with the excitation amplitude. As the periodic deformation excitation is inherently a disturbance, with the increase of the disturbance magnitude, the flow parameters fluctuate more strongly. Only when the energy input from the disturbance is big enough, can effective influence and modulation be acted on the flow field. As a result, the frequency of flow parameter becomes equal to that of the excitation. On different excitation frequencies, when the excitation amplitude is A = 5 × 10−3 c, the flow parameters vary periodically with time, with the frequency being equal to that of the excitation. The lift is larger and drag is smaller. The flow field in condition f = f 0 , A = 5 × 10−3 c is analyzed below in detail. As can be seen from Fig. 13(a), for the rigid airfoil, there exists a big separation zone on the upper surface and no reattachment occurs. Besides, the low-pressure area on the upper surface is small. When f = f 0 , A = 5 × 10−3 c, the separation region significantly decreases, with the separation point moving rearward and reattachment occurring, forming a flat and short separation bubble. Besides, the length of the low-pressure suction region increases and the pressure is obviously lower, so the lift increases. In Fig. 14, the equivalent shape is drawn based on the streamline on the boundary of the time-averaged separation zone. Take the separation bubble and the airfoil as a whole, and the qualitative analysis of the airfoil aerodynamic characteristics can be obtained based on the equivalent shape. The dash-dotted lines in Fig. 14 indicate the camber of the airfoil. For the rigid airfoil, comparing with the original airfoil shape, the equivalent shape has an increased thickness and a much smaller camber. The increase of thickness leads to larger form drag, and the decrease of camber leads to smaller lift. Thus, the rigid airfoil has a bad performance.

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112

47

113

48

114

49

115

50

116

51

117

52

118

53

119

54

120

55

121

56

122

57

123

58

124

59

125

60

126

61

127

62

128

63

129

64

130 131

65 66

Fig. 9. Flow control effects for different excitation amplitudes.

132

JID:AESCTE AID:105774 /FLA

[m5G; v1.282; Prn:13/02/2020; 14:11] P.7 (1-12)

J. Lei et al. / Aerospace Science and Technology ••• (••••) ••••••

7

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80 81

15

Fig. 10. Time histories of C L and C D when f = 0.5 f 0 .

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96 97

31

Fig. 11. Time histories of C L and C D when f = f 0 .

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112 113

47

Fig. 12. Time histories of C L and C D when f = 1.5 f 0 .

48 50 51 52 53 54 55 56 57 58

Fig. 13. Time-averaged pressure coefficient contours and streamlines.

59 60 61 62 63 64 65 66

114 115

49

Fig. 14. Equivalent shapes.

With flow control, the equivalent shape is very close to the original shape of the airfoil; the lift coefficient is satisfying and the low Reynolds number performance is improved. As can be seen in Fig. 15(a), in the condition with flow control, the pressure difference between the upper and lower surfaces increases, leading to the increase of the area enclosed by the C p curve, so lift is enhanced, and the center of pressure is closer to the leading edge. The c f curve in the flow control condition shows a separation point (c f = 0) closer to the trailing edge, and the length of the separation region sufficiently decreases. Besides, the reattachment occurs, which is in agreement with the flow field shown in Fig. 13. The difference in the locations of the transition to turbulence in these two conditions can also be seen in Fig. 15(b). For rigid airfoil, the absolute value of the c f abruptly increases at about 0.9c, (although not very clear), indicating the transition to turbulence. In the flow control condition, this transition occurs

116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

8

AID:105774 /FLA

[m5G; v1.282; Prn:13/02/2020; 14:11] P.8 (1-12)

J. Lei et al. / Aerospace Science and Technology ••• (••••) ••••••

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80 81

15 16

Fig. 15. Time-averaged aerodynamic coefficients.

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96

31

97

32

98

Fig. 16. Unsteady wall pressure coefficient.

99

33 34

100

Table 4 Coefficients of total drag and drag components.

35 36 37

Rigid airfoil f = f 0 , A = 5 × 10−3 c

38 39

101

CD

C Dp

CDf

0.05087 0.03137

0.03677 0.01661

0.01410 0.01476

Fig. 17. Unsteady vorticity magnitude contours and streamlines of the rigid airfoil.

43 44 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

104 105 107

41

45

103

106

40 42

102

at about 0.63c, which is sufficiently earlier than the former. The early transition may be a result of the disturbance produced by the surface oscillation, and it also leads to the expansion of the turbulent flow region. The wall shear stress and skin friction in the turbulent boundary layer are usually larger than those in the laminar boundary layer, as a result of the enhanced entrainment and mixing, so the expansion of the turbulent flow region may result in the increase of friction drag. However, in the situation of the present research, with the existence of the separation bubble, a reversed flow region is generated, and the wall shear stress is in a direction opposite to the free stream, so the friction drag is decreased, instead. In fact, the difference in the friction drag between the rigid airfoil condition and the flow control condition is small, as is shown in Table 4. On the other hand, the pressure drag is mainly affected by the equivalent shape, and significant differences exist between the rigid airfoil and the flow control conditions, so there is an evident difference in the total drag. For the rigid airfoil in Fig. 16(a), the pressure coefficient keeps approximately stable in one period of vortex shedding, and only small variations can be found around the trailing edge. As can be seen from Fig. 17, the flow past the upper surface of the rigid airfoil does not reattach after separation and forms a large sepa-

ration vortex — the primary vortex, of which the length is about 60% chord, and the vortex core moves periodically in a small range near the trailing edge. The primary vortex keeps attached to the airfoil and a periodically shedding secondary vortex is generated near the trailing edge. The primary vortex covers a large part of the upper surface, and its position and vorticity magnitude are relatively stable. As a result, the pressure distribution on much of the upper surface is also relatively stable, and only near the trailing edge the unsteady characteristics is obvious (Fig. 16(a)). The above phenomena agree well with the behavior of the trailing edge laminar separation bubble (TLSB) described in Ref. [29]. As a comparison, the condition with flow control is shown in Fig. 16(b). When t = 0, surface in the oscillating region moves upward and passes the equilibrium position. When t = 0.25T , the displacement reaches maximum and the local camber increases, resulting in a suction peak on the C p curve. When t = 0.5T , surface in the oscillating region returns to the equilibrium position and moves downward, and the suction peak mentioned above decreases. The pressure of the upper airfoil surface increases, and correspondingly the pressure difference between the upper and lower surfaces decreases. When t = 0.75T , the displacement reaches maximum in the negative direction, and the local camber decreases, resulting in local high-pressure. In one period, the vortices indicated by the low-pressure peaks keep moving downstream, and finally shed from the trailing edge.

108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:105774 /FLA

[m5G; v1.282; Prn:13/02/2020; 14:11] P.9 (1-12)

J. Lei et al. / Aerospace Science and Technology ••• (••••) ••••••

9

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91 92

26

Fig. 18. Flow control effects on different excitation frequencies.

27 29 30

93 94

28

95

Table 5 Excitation frequencies.

96

31

Number

f / f0

Number

f / f0

Number

f / f0

97

32

1 2 3 4 5 6

0 0.25 0.50 0.75 1.00 1.25

7 8 9 10 11 12

1.50 1.75 2.00 2.25 2.50 2.75

13 14 15 16 17

3.00 3.25 3.50 3.75 4.00

98

33 34 35 36

99 100 101 102 103

37

104

38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

3.2. Effect of the excitation frequency According to the above results, when A = 5 × 10−3 c, larger lift and smaller drag can be obtained, and the lift-to-drag performance is relatively good. In the present section, the oscillating amplitude is defined as A = 5 × 10−3 c, and the effect of the excitation frequency is studied. The excitation frequencies vary from 0 (rigid airfoil) to 4 f 0 , and 17 conditions are studied, as are listed in Table 5. The time-averaged results are analyzed first. As can be seen from Fig. 18(a), when the excitation frequency f is between 0.25 f 0 and 2 f 0 , the increase of lift is evident comparing with the rigid airfoil, and the maximum increment can be 54%. When f > 2 f 0 , the lift enhancement effect decreases suddenly, and lift finally becomes similar to the value of the rigid airfoil. A tendency of first decreasing and then increasing occurs in the drag characteristics, as is shown in Fig. 18(b). The minimum drag locates around f = f 0 , and the value is 62% of the rigid airfoil. When f > 2 f 0 , the drag reduction effect almost disappears, and the difference between the flow control and the rigid airfoil conditions is slight. As a result, the largest lift-to-drag ratio appears around f = f 0 , and the lift-to-drag ratio is 119% larger than the rigid airfoil, which evidently demonstrates the effect of the present flow control method. When f > f 0 , the lift-to-drag ratio decreases with the excitation frequency, and finally converges to that in the rigid airfoil condition. From the results shown above, when the excitation

Fig. 19. Time-averaged streamlines, pressure coefficient contours and equivalent shapes.

105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121

Fig. 20. Time-averaged pressure coefficients.

122 123

frequency is f = f 0 , the best control effect can be obtained and the increase of C L /C D is the maximum. As can be seen from Fig. 19, with the increase of the excitation frequency, the length and thickness of the time-averaged separation bubbles on the upper surface also increase. As is interesting, when f = 2 f 0 , the separation bubble consists of 4 vortices, including 2 primary vortices and 2 secondary vortices, resulting in a long low-pressure area on the rear of the upper surface (also shown in Fig. 20). Besides, when f = f 0 , the equivalent shape

124 125 126 127 128 129 130 131 132

JID:AESCTE

10

AID:105774 /FLA

[m5G; v1.282; Prn:13/02/2020; 14:11] P.10 (1-12)

J. Lei et al. / Aerospace Science and Technology ••• (••••) ••••••

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79 80

14 15

Fig. 21. Unsteady streamlines and velocity magnitude contours in one excitation period.

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109 110

44 45

Fig. 22. Flow control effects at different angles of attack.

47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

111 112

46

is very close to the original airfoil shape, and the variations of camber and thickness are the slight, leading to fine performance. When f = 2 f 0 or f = 3 f 0 , the equivalent shapes have much difference comparing with the original airfoil, and the cambers become rather flat. The differences in the equivalent shapes can help to explain the dependencies of the C L ’ and C D ’ on the oscillation frequency. Then the unsteady flow field in Fig. 21 is analyzed. For the rigid airfoil, there exists a large separation area on the upper surface, and it consists of a primary vortex and a secondary vortex induced near the trailing edge. The flow velocity in the separation area is very low, and the suction of the upper surface is weak. When f = f 0 , the separation area significantly reduces, and a series of small vortices travel downstream along the upper surface. The vortex-induced flow acceleration (red areas in Fig. 21) leads to lower pressure on the upper surface and produces larger lift. When f = 2 f 0 , as the excitation frequency is higher, more vortices are generated and the interval becomes smaller. Some vortices begin to merge. Although the flow separation is widespread compared

113 114 115 116 117 118 119 120

Fig. 23. Time-averaged streamlines.

121 122

with the f = f 0 case, the vortex-induced flow acceleration area is longer and closer to the trailing edge, so the suction area is larger, leading to even larger lift. When f = 3 f 0 , with the higher excitation frequency, the vortices are more crowded compared with the f = 2 f 0 case and travel away from the wall. Some small vortices merge into a big one and the flow velocity is very low in it. The size and the flow structures of the separation zone is similar to the rigid airfoil case, so the pressure distributions of these two cases are also similar. The time-averaged pressure distributions in Fig. 20 also support the analyses.

123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:105774 /FLA

[m5G; v1.282; Prn:13/02/2020; 14:11] P.11 (1-12)

J. Lei et al. / Aerospace Science and Technology ••• (••••) ••••••

11

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93 94

28 29

Fig. 24. Time-averaged wall pressure coefficients.

31 32 33 34 35 36 37

Table 6 Flow natural frequencies and excitation frequencies. Angle of attack α (◦ )

Natural frequency f 0 (Hz)

f = 0.5 f 0 (Hz)

f = f0 (Hz)

f = 1.5 f 0 (Hz)

f = 2 f0 (Hz)

2 4 6

15.97 13.26 10.97

7.99 6.63 5.49

15.97 13.26 10.97

23.96 19.89 16.46

31.94 26.52 21.94

38 39 40

3.3. Flow control effect at different angles of attack

41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

95 96

30

The above results and analyses are all based on the angle of attack α = 4◦ . In this section, a preliminary study of the applicability of the flow control method at different angles of attack is conducted. As are shown in Table 6, the conditions combining different angles of attack and excitation frequencies are studied. The f 0 is the flow natural frequency of the rigid airfoil, and its value is equal to the reciprocal of vortex shedding 1/ T , obtained from the CFD simulation result at the corresponding angle of attack. For the rigid airfoil, the frequency f 0 usually decreases with the angle of attack. As can be seen in Fig. 22, the increase of lift and decrease of drag are evident at each angle of attack studied. When the excitation frequency f is around f 0 , the increase of lift-drag ratio can be better. The effect of lift-drag ratio enhancement becomes better with the increase of the angle of attack, and the increment of lift-drag ratio can be about 86%, 116% and 134% at the angles of attack α = 2◦ , 4◦ and 6◦ , respectively. The results shown above demonstrate that the flow control method in this paper is effective at different angles of attack studied. As can be seen from Fig. 23, for the rigid airfoil, the separation bubble on the rear of the upper surface becomes larger with the increase of the angle of attack, and the location of the separation point moves upstream. With flow control, however, the separation bubble at a larger angle of attack is smaller than that at a smaller angle of attack. As a result, the lift-drag ratio enhancement

is larger at a higher angle of attack, and the flow control effect is better. As can be seen in Fig. 24, for all the 3 angles of attack, with flow control, the pressure coefficient on the front part of the upper surface decreases, and an overall increase occurs on the lower surface, resulting in the enhancement of lift. Besides, the center of pressure moves forward with flow control. 4. Conclusions

97 98 99 100 101 102 103 104 105 106

In this paper, flow control by means of active oscillation of local surface on an E387 airfoil is studied. The chord Reynolds number is Re = 30000 and the oscillation region is on the upper surface near the leading edge. The flow control effects at different oscillation amplitudes, excitation frequencies and angles of attack are investigated. Some conclusions could be drawn: (1) With flow control, lift enhancement and drag reduction are evident at the angle of attack α = 4◦ . The maximum increment of lift, decrement of drag, and increment of lift-drag ratio can be 49%, 39% and 118%, respectively. When the oscillation amplitude A is less than 2 × 10−3 c, the variation of aerodynamic coefficient with the oscillation amplitude is obvious; and the lift-drag ratio is nearly maximum when A = 5 × 10−3 c. (2) The best flow control effect occurs when the excitation frequency is equal to the frequency of vortex shedding in the separated shear layer, namely the flow natural frequency. When the excitation frequencies are f < 2 f 0 , the lift enhancement and drag reduction effect is satisfactory. When f > 2 f 0 , with the increase of f , the flow control becomes ineffective. (3) With flow control, the traveling vortices and the resulting flow acceleration and low-pressure region are the direct cause for the lift increase. The laminar separation bubble leads to the variation of the camber and thickness of the equivalent shape, effecting the lift and drag. (4) For the rigid airfoil, the frequency of the vortex shedding decrease with the increase of the angle of attack. At different an-

107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

AID:105774 /FLA

12

1 2 3

[m5G; v1.282; Prn:13/02/2020; 14:11] P.12 (1-12)

J. Lei et al. / Aerospace Science and Technology ••• (••••) ••••••

gles of attack, the maximum increase of lift-drag ratio all occurs near f = f 0 , with the maximum increments being 86%, 116% and 134% at α = 2◦ , 4◦ and 6◦ , separately.

4 5

Declaration of competing interest

6 7 8 9

The authors declare that they do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.

10 11

References

12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

[1] P.B.S. Lissaman, Low-Reynolds-number airfoils, Annu. Rev. Fluid Mech. 15 (1) (1983) 223–239. [2] H.P. Horton, Laminar Separation Bubbles in Two and Three Dimensional Incompressible Flow, Queen Mary University of London, 1968. [3] M. Gaster, The structure and behaviour of laminar separation bubbles, 1967, TIL. [4] R. Kojima, T. Nonomura, A. Oyama, et al., Large-eddy simulation of lowReynolds-number flow over thick and thin NACA airfoils, J. Aircr. 50 (1) (2013) 187–196. [5] M. Gad-el-Hak, Flow Control: Passive, Active, and Reactive Flow Management, Cambridge University Press, 2006. [6] R.D. Joslin, D.N. Miller, Fundamentals and Applications of Modern Flow Control, American Institute of Aeronautics and Astronautics, 2009. [7] D. Kamari, M. Tadjfar, A. Madadi, Optimization of SD7003 airfoil performance using TBL and CBL at low Reynolds numbers, Aerosp. Sci. Technol. 79 (2018) 199–211. [8] S. Cao, Y. Li, J. Zhang, et al., Lagrangian analysis of mass transport and its influence on the lift enhancement in a flow over the airfoil with a synthetic jet, Aerosp. Sci. Technol. 86 (2019) 11–20. [9] H. Zhang, S. Chen, Q. Meng, et al., Flow separation control using unsteady pulsed suction through endwall bleeding holes in a highly loaded compressor cascade, Aerosp. Sci. Technol. 72 (2018) 455–464. [10] D. Neuburger, I. Wygnanski, The use of a vibrating ribbon to delay separation on two-dimensional airfoils: some preliminary observations, in: Proc. Workshop II on Unsteady Separated Flow, 1988, pp. 333–341. [11] S.K. Sinha, Flow separation control with microflexural wall vibrations, J. Aircr. 38 (3) (2001) 496–503. [12] X. Meng, H. Hu, X. Yan, et al., Lift improvements using duty-cycled plasma actuation at low Reynolds numbers, Aerosp. Sci. Technol. 72 (2018) 123–133. [13] D. Munday, J. Jacob, Active control of separation on a wing with oscillating camber, J. Aircr. 39 (1) (2002) 187–189.

[14] N.J. Pern, Jacob J. Comparison, Between zero mass flux flow control methods for low speed airfoil separation control, in: 35th AIAA Fluid Dynamics Conference and Exhibit, 2005, p. 4884. [15] V. Katam, R. LeBeau, J. Jacob, Simulation of separation control on a morphing wing with conformal camber, in: 35th AIAA Fluid Dynamics Conference and Exhibit, 2005, p. 4880. [16] W. Kang, J.Z. Zhang, P.H. Feng, Aerodynamic analysis of a localized flexible airfoil at low Reynolds numbers, Commun. Comput. Phys. 11 (4) (2012) 1300–1310. [17] W. Kang, P. Lei, J. Zhang, et al., Effects of local oscillation of airfoil surface on lift enhancement at low Reynolds number, J. Fluids Struct. 57 (2015) 49–65. [18] G. Di, Z. Wu, D. Huang, The research on active flow control method with vibration diaphragm on a NACA0012 airfoil at different stalled angles of attack, Aerosp. Sci. Technol. 69 (2017) 76–86. [19] G. Jones, M. Santer, G. Papadakis, Control of low Reynolds number flow around an airfoil using periodic surface morphing: a numerical study, J. Fluids Struct. 76 (2018) 95–115. [20] A. Gross, H. Fasel, Numerical investigation of separation for airfoils at low Reynolds numbers, in: 40th Fluid Dynamics Conference and Exhibit, 2010, p. 4736. [21] R.I. Issa, Solution of the implicitly discretised fluid flow equations by operatorsplitting, J. Comput. Phys. 62 (1) (1986) 40–65. [22] F.R. Menter, R.B. Langtry, S.R. Likki, et al., A correlation-based transition model using local variables—part I: model formulation, J. Turbomach. 128 (3) (2006) 413–422. [23] R.B. Langtry, F.R. Menter, S.R. Likki, et al., A correlation-based transition model using local variables—part II: test cases and industrial applications, J. Turbomach. 128 (3) (2006) 423–434. [24] F.J. Blom, Considerations on the spring analogy, Int. J. Numer. Methods Fluids 32 (6) (2000) 647–668. [25] J. Niu, J. Lei, T. Lu, Numerical research on the effect of variable droop leadingedge on oscillating NACA 0012 airfoil dynamic stall, Aerosp. Sci. Technol. 72 (2018) 476–485. [26] L.E. Jones, R.D. Sandberg, N.D. Sandham, Direct numerical simulations of forced and unforced separation bubbles on an airfoil at incidence, J. Fluid Mech. 602 (2008) 175–207. [27] M.S. Selig, Summary of low speed airfoil data, SoarTech, 1995. [28] D. Pande, D. Verstraete, Impact of solar cell characteristics and operating conditions on the sizing of a solar powered nonrigid airship, Aerosp. Sci. Technol. 72 (2018) 353–363. [29] P. Bai, F. Li, Q. Liu, et al., Evolution of the non-linear and unsteady low Reynolds number laminar separation bubble around the airfoil with small angle of attack, in: 46th AIAA Fluid Dynamics Conference, 2016, p. 4337.

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112

47

113

48

114

49

115

50

116

51

117

52

118

53

119

54

120

55

121

56

122

57

123

58

124

59

125

60

126

61

127

62

128

63

129

64

130

65

131

66

132