- Email: [email protected]
Detached-eddy simulation of flow past a pitching NACA 0015 airfoil with pulsed actuation
Accepted Manuscript Detached-eddy simulation of flow past a pitching NACA 0015 airfoil with pulsed actuation
Liang Wang, Song Fu
PII: DOI: Reference:
S1270-9638(16)30881-1 http://dx.doi.org/10.1016/j.ast.2017.06.002 AESCTE 4058
To appear in:
Aerospace Science and Technology
Received date: Revised date: Accepted date:
19 October 2016 2 June 2017 3 June 2017
Please cite this article in press as: L. Wang, S. Fu, Detached-eddy simulation of flow past a pitching NACA 0015 airfoil with pulsed actuation, Aerosp. Sci. Technol. (2017), http://dx.doi.org/10.1016/j.ast.2017.06.002
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Detached-Eddy Simulation of flow past a pitching NACA 0015 airfoil with pulsed actuation Liang Wang1, [, and Song Fu1,] Tsinghua University, Beijing, China, 100084
1
This study conducted Detached-Eddy Simulations (DES) of actively controlled flow over a 0.35 m chord NACA0015 airfoil at an incidence angle of 11o and a chord Reynolds number of 1 million. The the uncontrolled, natural flow were simulated by using DES-type methods, with a comparison of different underlying RANS models and subgrid-scale stress models in LES mode. Results from these computations were compared with experimental observations, enabling their reliable assessment through the detailed investigation of the Reynolds stresses as well as the separation and reattachment. It was found that among all five DES-type methods, only the Spalart-Allmaras-based Improved Delayed DES (IDDES-SA) captured the separation point as measured in the experiments. The classical vortex-shedding and the shear-layer flapping modes for airfoil flows with shallow separation were also extracted from the IDDES-SA results by using Dynamic Mode Decomposition. With better understanding of unsteady flow features, effective control practices were illustrated. The pulsed blowing actuation was adopted according to the experimental settings. The flow-field around jet orifices was resolved, instead of the use of actuation boundary conditions. It was found that for all actuation frequencies investigated in this article, the saturation of lift-to-drag ratio was yielded already at intermediate momentum coefficients, two-fifths of the actuation amplitude used in the experiment; the optimal actuation frequency was 0.6 among all cases, consistent with the dominant frequency of the baseline flow. The suppression of separation bubble led to a maximum lift-to-drag ratio enhancement of 194%, as well as a sharp decrease of all Reynolds stress components. To
[ Assistant professor, School of Aerospace Engineering. ] Professor, School of Aerospace Engineering. Email: [email protected] 1
further reduce the strong unsteadiness introduced by excitation, a phase shift among jet orifices was applied and evaluated. The pulsed blowing with the jet orifices in an anti-phased manner was found to be very effective to damp the fluctuations of the lift around its mean value, though it led to a small decrease in mean lift.
Nomenclature Į
= angle of attack, deg
Ahole
= area of jet orifice, 7.85E-7 m2
Awing
= area of wing surface per orifice, 5.25E-3 m2
c
= chord length, m
CDES
= constant taken as 0.628 for the SA model and 0.550 for the SST model
CL
= lift coefficient
Cp
= pressure coefficient
Cȝ
= jet momentum coefficient
ǻt
= time step, s
DES
= Detached-Eddy Simulation
DDES
= Delayed DES
IDDES
= Improved DDES
fd
= blending function
fex
= frequency of periodic excitation, 1 / s
F*
= fex × c sinĮ /U, non-dimensional excitation frequency
HRLM
= Hybrid RANS/LES Method
k
= turbulent kinetic energy, m2 / s2
Rec
= Reynolds number based on c
R.M.S.
= Root Mean Square
SA
= Spalart-Allmaras
SMG
= Smagorinsky
SST
= Shear Stress Transport
St
= f c sinĮ/ U, Strouhal number
U
= freestream velocity, m / s
ua , ujet
= amplitude of the excitation velocity, excitation velocity, m / s
WALE
= Wall-Adapting Local Eddy-viscosity
xsep
= separation location, m
ȝ
= dynamic viscosity, ȝ = ȡȞ
Ȟ
= kinetic viscosity, m2 / s
2
Ȟt
= eddy viscosity, m2 / s
Ȧ
= specific turbulence dissipation rate, 1 / s
IJ
= t (U / c), time in convective units
Subscripts
= free-stream value
w
= wall
I. Introduction The determination of flow separation regions is a key task in the design process of wings and turbines, which beyond the calibration regime of the conventional RANS (Reynolds-Averaged Navier-Stokes) models that are not always consistent in their prediction capability. Though LES (Large-Eddy Simulation) is a viable alternative to RANS for accurate computations, the near-wall resolution required for LES makes it expensive. This is the rationale behind the continual development and application of HRLM (Hybrid RANS/LES Method), which combines advantages from both sides. In order to exhibit both LES and RANS behavior, a simple modification of the Spalart-Allmaras model1 was proposed by Spalart2 such that the computation behaves as LES in the unsteady "detached" regions and as RANS in the "attached" near-wall regions of the flow. This scheme, known as DES (Detached-Eddy Simulation), has been widely applied for a plethora of external flows since its inception. Following this development, two improved DES variants, DDES (Delayed DES)3 and IDDES (Improved DDES)4, have been further proposed by the coworkers of Spalart. This approach was also applied by Travin et al.6 and further by Gritskevich et al.22 to the RANS SST (Shear Stress Transport) model5. However, in DES-type methods, the turbulent boundary layers upstream of a free shear layer are typically in RANS mode and therefore do not contain any resolved turbulence. As a consequence, the development of turbulence in the shear layer will start from the intrinsic instability of the shear layer and will not be forced by turbulence coming from the boundary layers. This resembles the case for laminar boundary-layer flows upstream, where there is no sufficient turbulence due to slow shear-layer development, known as the "Grey-Area" problem. To address this problem, Mockett et al.7 used the WALE SGS (subgrid-scale stress) model8 to replace the conventional Smagorinsky model, resulting in a
3
dramatic acceleration of RANS/LES transition. Note that this approach retained generality, as demonstrated for several free-shear flow cases including a plane shear layer, a backwards-facing step and a round jet. Recently, Siauw et al.9, 10 experimentally studied the actively controlled flow over a 0.35 m chord NACA0015 airfoil at an incidence angle of 11o and a chord Reynolds number of 1 million. With the use of a pulsed-jet actuation, they successfully eliminated the trailing edge separation and therefore improve the aerodynamic performance of the airfoil. It was then selected as a standard test case by the EU-China project MARS (Manipulation of Reynolds Stress for Separation Control and Drag Reduction). The researchers within this project tried to numerically validate the experimental findings using HRLM. However, their results agree poorly with the experimental data for the uncontrolled case where the flow separation is mild11. Besides, in the experiment the pulsating frequency is of 1 Hz while the maximum jet velocity is of 200 m/s. The very low actuation frequency makes the pulsed jet similar to continuous jet blowing while the actuation intensity is comparably too high for real applications. The primary objective of this work is thus to numerically investigate the effects of the flow-control parameters, such as intensity and frequency, on the mildly separated flow on a NACA0015 airfoil, by using DES-type methods. The secondary objective is to assess the DES-type methods used for the baseline, natural flow case, with a comparison of different choices of underlying RANS model as well as SGS model in LES mode. This paper is organized as follows: Section 2 presents the mathematical formulations including DES-type method (DDES and IDDES variants) methods as well as alternative SGS models; Section 3 gives numerical details, especially for mesh generation and the modeling of excitation mechanisms; the numerical results of both baseline and controlled flows, including flow quantities and the subsequent spectra analyses, are presented in Section 4; the conclusions of this study are summarized in Section 5.
4
II. Mathematical Formulation In this section, one-equation SA (Spalart-Allmaras) model1 and k–Ȧ SST model5 are presented first to start the introduction of the DES type methods used here. Then, we express the DES method with DDES and IDDES variants. Finally, the alternative SGS (Smagorinsky and WALE) models in LES Mode for DDES are given. The details of the model formulations can be found in the literatures. A. The RANS Spalart-Allmaras Model and Menter’s SST Model The SA model and the SST model are the work-horses turbulence models of the aerospace industry. Both of them are eddy viscosity-based models. The SA model is based on a well-calibrated transport equation for a working quantity, as expressed as DUQ~ Dt
~ 1 ° w Cb1 US Q~ ® V ° wx j ¯
2 ~ 2 ª § wQ~ · ½° wQ~ º ~ ¸ ¾ C w Uf w §¨ Q ·¸ « P P » Cb 2 U ¨¨ ¸ wx j »¼ ©d ¹ «¬ © wx j ¹ °¿
(1)
~
where S is a function of the vorticity, Cb1, Cb2 and Cw the model constants and fw the empirical functions. The turbulent Reynolds stress tensor is then modelled by the eddy viscosity Pt with Pt
(2)
UQ f v
1
The SST model requires the solution of transport equations for k and Ȧ: wU k wk º w ª « U u j k P V k Pt » W S E *U kZ wt wx j »¼ ij ij wx j «¬
(3)
and w UZ wZ º w ª « Uu jZ P V Z Pt » wx j ¬« wt wx j »¼
PZ E kZ 2 2 1 F1
UV Z 2 wk wZ Z wx j wx j
(4)
where S is the mean strain rate, and F1 is a blending function. The eddy-viscosity ȝt is calculated from Pt
§ U k U a1k · min ¨ , ¸ © Z :F2 ¹
5
(5)
where a1=0.31 and F2 is also a blending function. The SA and SST model coefficients and empirical functions can be found in Ref.1 and Ref.5, respectively. B. IDDES Method It is designed in DES that the entire boundary layer is simulated by using RANS. The k equation in RANS, Eq. (3), can be rewritten as w U k w Uu j k wt wx j
w wx j
ª§ Pt «¨ P Vk ¬«©
· wk º U k 1.5 W S » ¸ ij ij LRANS ¹ wx j ¼»
(6)
where LRANS is the turbulence length scale for RANS. It equals to k0.5/Z for the SST model. IDDES4 is developed from DDES3 and it combines DDES and the wall-modeled LES (WMLES), with LRANS in Eq. (6) substituted by a length scale LIDDES
fd 1 fe LRANS 1 fd LLES
(7)
Here, the length scale for LES, LLES, is further written as the grid scale ' multiplied by a model constant CDES. fd and fe are both blending functions. It is seen from Eq. (7) that once fe equals to zero, LIDDES can be expressed as
LIDDES LDDES
fd LRANS 1 fd LLES
(8)
and IDDES reverts to DDES. Ref. 4 and Ref. 22 present the details of all blending functions. C. Alternative Model Formulations in LES Mode for DDES For DES-type models in LES mode, the following form can be derived under the assumption of local equilibrium (i.e. equality of the dissipation and generation terms of the underlying transport equations):
X SGS
Cm '
2
Dm (u )
(9)
where ǻ is the sub-grid characteristic length scale, and Dm (u) stands for a differential operator, acting on the resolved velocity field. For the Smagorinsky model
6
Dm
Ds
2Sij Sij , Cm Cs | 0.18ˈSij
1 wui wu j ( ) 2 wx j wxi
(10)
There are mainly two drawbacks in this model: one is that the differential operator Dm (u) doesn't vanish in near-wall regions, which gives a non-zero value of Ȟt as soon as a velocity gradient exists; from the formula, we can see it is only related to the strain rate of the turbulent structure but not to their rotation rate8. To solve these two problems, Mockett et al.7 adapted the WALE model8 that can be expressed as Dm
Sijd
Dw
( S ijd S ijd )3/2
(11)
( Sij Sij )5/2 ( S ijd S ijd )5/4
1 2 1 ( gij g 2ji ) G ij g kk2 2 3
, gij2
gik g kj
, gij
wui wx j
(12)
To modify DDES to behave like the WALE model in LES mode, we leave the length scale substitution, Eq. (8), unmodi¿ed and introduce an additional function to substitute the corresponding term for Dm (u) in the LES mode region only. The velocity gradient invariant in the underlying RANS model, S*W-DDES, is substituted by SW* DDES * S RANS
* * S RANS f d pos LRANS LLES S RANS BW SW*
2:ij :ij
, pos x Bw
0 , if x d 0 ® ¯1 , if x ! 0
CW2 / CS2
(13) (14) (15)
here Cw and Cs are model constants. Readers may refer to Ref.8 for all the coefficients. III. Numerical methods and grid generation A. Flow solver An implicit pressure-based solver with a fully conservative approximation of the governing equations is employed in the flow simulation. The code is based on curvilinear coordinates and uses cell-centered collocated storage arrangement on semi-block structured grids for all quantities. The one-equation SA
7
model1 and k–Ȧ SST model5 are used as the underlying models for the DDES3 and IDDES4presented here. The CDES parameter is set to 0.628 for the SA model and to 0.550 for the SST model, as calibrated against the decay of isotropic turbulence12. In order to best address the inconsistency in the demands posed by RANS and LES on the numerical scheme for the convective fluxes, the blending function of Travin et al.6 is employed to assure an appropriate switch between a higher-order TVD scheme and central differencing. The diffusive fluxes are approximated using a second order central scheme and for time discretization second order backward differencing is applied. The continuity equation is conserved by the SIMPLE algorithm whereby the decoupling of pressure and velocity is prevented through a modified Rhie & Chow interpolation13. The code is parallelized via domain decomposition and the data interchange between processors is realized through the standardized MPI-library. RANS solutions obtained with the SA or SST models are used as initial condition for the DDES and the IDDES. B. Flow configuration The numerical test model (see Fig. 1a) represents a NACA 0015 airfoil at a chord Reynolds number of about 1 million in response to the deployment and removal of actuation orifices positioned at 30% of the chord length c ( = 0.35 m) of the airfoil. The influence of the incidence angle on the separation was carefully analyzed in the experiments9, 10. It was found that at 11o of incidence the separation was two-dimensional with minimal vibration. The chosen flow condition is such that jet deployment corresponds to complete flow attachment over the airfoil; jet removal will cause flow separation up to an extent defined by the uncontrolled separated state. For the pulsed-blowing flow control, an array of 44 holes is deployed at 0.3c down-stream of the leading edge of the airfoil. This array of holes occupies the central one-third spanwise portion of the airfoil with the holes spaced apart in the span with an interval of 0.043 c. The jets are realized by blowing air into tubes connecting to the holes distributed on the airfoil surface. The tubes have a diameter of 0.001 m, and the jet Àow from the tubes has a pitch angle of 30o and a skew angle of 60o, as shown in Fig. 1b.
8
The Àow area at the jet exit on the airfoil surface is an ellipse with a major axis length of 0.002 m and a minor axis length of 0.001 m.
(a) (b) Figure 1. Schematics of the NACA0015 airfoil (a); of the pulsed-jet actuator position and orientation (b).
C. Computational Mesh In the previous work14, two 2-D grids of 0.10 million (see Fig. 2b) and 0.14 million nodes, respectively, were generated to test the grid independence. The 3-D grid was based on an expansion of the 2-D mesh slice into the spanwise direction. 101 layers of the 2-D mesh were combined. As a result of refining the grid, the turbulent viscosity is slightly lowered in the free shear layer. This indicates a more LES-like simulation than with baseline grid. However, it already provides a good representation of the turbulent structures on the baseline grid. It is also proven that the use of periodic boundaries is valid in case of an adequate spanwise domain size. The computational domain adopted is of size LC×LN×LZ = 17.1c×7.4c×0.18c and the total number of grid is 13.5 million nodes. The spanwise length is designed for four holes to be uniformly distributed. According to Spalart16, both the focus region (above the airfoil while downstream the separation point at x/c = 0.7) and the departure region (wake), shown in Fig. 2b in red, need to be resolved sufficiently to be in LES mode. ǻx+, ǻy+ and ǻz+ are no larger than 33.7 in the focus region. The non-dimensional wall distance of the first cell center remains below y+ =1 on the entire surface.
9
(a) block topology
(b) meshes around the airfoil
(c) zoom-in at trailing edge
(d) jet orifice configuration
(e) x-y plane mesh nearby jet orifice Figure 2. Computational mesh.
D. Numerical setup and parameters
RANS simulations are performed prior to the unsteady simulation using the same turbulence model as used for the DESs. All the DESs have been run for more than 20 convective time units, IJ = t U / c, before the evaluation of averages and the storage of the surface data was initiated. The sampling times are at least 30 convective time units for the DESs.
10
The prescribed non-dimensional time step size of ǻt = 2.5 ×10í2 c sinĮ /U guarantees CFLmax < 1.0 in the separation region17 as well as providing accurate resolution up to St = 2. The spectra will be displayed up to 1/3 of the Nyquist frequency, which is St = 20. Roughly 2 orders of magnitude residual reduction are reached within 10 outer iterations per time step. At the wind tunnel entry all flow quantities including the velocity components and turbulent properties are prescribed. The level of turbulence at the inflow is set to Tu = (2k/3) 0.5 / U = 0.5% and the turbulent viscosity ȝt / ȝ = 0.1. Previous investigations by the present authors23, 24 led to the choice of these values. At the outflow a convective boundary condition is used that allows unsteady flow structures to be transported outside the domain. The complete airfoil surface is modeled as a non-slip boundary condition. The wind tunnel walls are approximated as inviscid walls. Periodic boundary conditions are set in the spanwise boundaries. The resolved actuation orifices (see Fig. 2d) are considered in the present paper. The actuation velocity inlet is placed at its bottom. The time-dependent velocity for the pulsed actuation (see Fig. 3) is modeled as follows:
W W DC T 1 1 1 1 u jet ( x, t ) ua ( x) [ tanh( p )] [ tanh( p )] (20) 2 2 ts 99 2 2 ts 99 with
t pU f c sin D , t s 99 f sc , DC ,T , W p * F Uf c sin D (21) where ua is the amplitude velocity of the time-dependent excitation, F* the non-dimensional actuation T
frequency, T the cycle duration of the pulse, ts99 the time the pulsed jet achieving its maximum value, sc the smoothing parameter of the flank of the pulse, DC the duty cycle of the pulse, IJp the dimensionless pulse time given in convective units, and tp the pulse time. tp is reset after each cycle duration so that a periodic behavior is generated. Note that the jet exits are set as solid wall surfaces in the baseline computations.
11
Figure 3. Modeling of the time-dependent part of the pulsed actuation.
IV. Results and discussions In the following, it first presents the results of the baseline, natural flow, with the assessment of five different DES-type models through the detailed investigation of the Reynolds stresses as well as the separation and reattachment. Second, the results of the controlled simulations are given. These will be divided into results of the investigation of the excitation intensity as well as frequency, and results of a phase shift among jet orifices. A. Baseline Flow The various results reported herein correspond to the five different models employed during the study. Results of the DDES and IDDES computations are labeled DDES and IDDES, respectively. Then, results labeled SA and SST respectively correspond to those predictions in which the one-equation SA model and k–Ȧ SST model are used as the underlying RANS models. Finally, the results labeled by SMG have been obtained using the standard SGS model (Smagorinsky model) while by WALE using WALE SGS model in LES mode. For instance, DDES-SA-WALE represents the DDES computations with the SA RANS model as well as the WALE SGS model.
12
(a) DDES-SA-SMG
(b) DDES-SA-WALE
(c) DDES-SST-SMG
(d) IDDES-SST
(e) IDDES-SA
Figure 4. Instantaneous flow using the Ȝ2 vortex core criterion at Ȝ2 = - 500 s-2, shaded by dimensionless velocity magnitude.
Fig. 4 compares the instantaneous turbulence structures (iso-surface of the second invariant of velocity gradient). It is seen that the grid resolution remains sufficiently fine enough to resolve the intrinsic small-scale turbulent structures in the shear layer for all five models. The calculated separation locations from the five models are given in Table 1, which are de¿ned using the zero-Cf criterion in the distribution of the streamwise component of skin friction, Cfx, as shown in Fig. 5. The measured separation location (estimated by oil-Àow visualization), xsep, at 68% of the chord length is fairly well predicted
by
both
IDDES-SA
and
DDES-SA-WALE,
whereas
both
DDES-SA-SMG
and
DDES-SST-SMG models predict a delayed separation compared to experimental results, and therefore there is little turbulence generation and flow structures shown in Fig. 4a and Fig. 4c. On the contrary,
13
IDDES-SST gives an over-prediction in the bubble length (see Fig. 5) as well as in the turbulence structures (see Fig. 4d).
Figure 5. Time- and spanwise-averaged skin-friction distribution along the airfoil. Table 1. Comparisons of separation location between computations and measurements.
Methods x
sep
DDES-SA-SMG
DDES-SA-WALE
DDES-SST-SMG
IDDES-SST
IDDES-SA
Measurements
0.849c
0.676c
0.886c
0.573c
0.649c
0.687c
Figure 6. Time- and spanwise-averaged surface pressure distribution along the airfoil.
14
The pressure distribution contributes the most to the lift enhancement. The time-averaged surface pressure plot Cp in Fig. 6, including zoomed views near the trailing edge, indicates that IDDES-SA gives accurate prediction in the trailing-edge separated-flow region that is simulated with LES. The overshoot of both DDES-SA-SMG and DDES-SST-SMG in Cp at trailing edge can be attributed to their over-prediction in xsep and vice-versa for the IDDES-SST results. Moreover, the initial dip in all DES results in the leading-edge section can be attributed to the tripping applied at 0.4% c from the leading edge in the experiment.
Figure 7. Time- and spanwise-averaged streamwise velocity profiles at x/c= 0.84 (left) and 0.97 (right).
Figure 8. Time- and spanwise-averaged streamwise normal stress < u'u'> profiles at x/c= 0.84 (left) and 0.97 (right).
15
Fig. 7 and Fig. 8 compare the time-averaged streamwise velocity profiles and streamwise normal stress, respectively, along the wall-normal direction at x/c = 0.84 and 0.97 with the available experimental results. Note that the two profiles are located across the separation bubble. It is observed that the measured profiles are in errors, where the velocity does not show negative sign for the reverse flow in the separation bubble. Due to capturing incorrect separation location, IDDES-SST predicts too an
intensive reverse flow in the separation bubble, while both DDES-SA-SMG and DDES-SST-SMG show no negative sign in velocity profiles and pretty weak turbulence intensity. IDDES-SA and DDES-SA-WALE give similar velocity profiles, whereas their discrepancy is huge in Reynolds stresses. Although there are some gaps between the results by IDDES-SA and measurements, the consistent distribution trends is obtained. On the contrary, turbulence intensity predicted by DDES-SA-WALE is close to zero.
Figure 9. Comparison of streamwise velocity (left) and shear stress (right) profiles at x/c=1.98 downstream the trailing edge.
Fig. 9 displays the time-averaged streamwise velocity and shear stress profiles one chord length downstream the trailing edge. Only IDDES-SA predicts agreeable velocity and Reynolds stress distributions. Considering the larger separation provided by IDDES-SST, it over-predicts the wake width compared with the experiments. It is vice-versa for both DDES-SA-SMG and DDES-SST-SMG results. However, DDES-SA-WALE also gives an under-prediction in wake width, though it predicts the accurate separation location upstream. Note that there exist slight “wriggles” in both IDDES computations, which can be attributed to the higher grid-density requirement of the IDDES approach.
16
(a)
(b)
Figure 10. Sampling time series of the global lift coefficient for (a) IDDES-SA and (b) DDES-SA-WALE. IJ = t (U/ c), time in convective units.
(a)
(b)
Figure 11. Spectra of the global lift coefficient for (a) IDDES-SA and (b) DDES-SA-WALE.
To further compare the results of IDDES-SA and of DDES-SA-WALE, spectra analyses on the unsteady flow quantities are carried out by using Discrete Fourier Transform (DFT). Fig. 11 presents the power spectra density of the global lift coefficient as a function of the Strouhal number, St = f csinĮ /U, based on the projected airfoil height csinĮ. It should be noted that the smallest frequency (or the St number in Fig. 11) is related to the total time of the sampling (see Fig. 10), and the largest frequency is determined by the sampling frequency, which equals to 1/ǻt = 40 U / (csinĮ). It is observed from Fig. 10 that the R.M.S. value of global lift coefficient sharply decreases from 0.042 for IDDES-SA to 0.0063 for DDES-SA-WALE results, and from Fig. 11 that the peak value in the DDES-SA-WALE spectra is extraordinarily intensive compared to that in IDDES-SA results. The high-frequency fluctuations of
17
DDES-SA-WALE results (see Fig. 11b) may be regarded as the spurious numerical noises that cannot be dissipated with the low level of eddy viscosity in the separation bubble (see Fig. 8b). Moreover, recall that DDES-SA-WALE predicts accurate separation location however near-zero Reynolds stress in the separation bubble (see Fig. 8). We can conclude that the separated flow simulated by DDES-SA-WALE is actually laminar. It is noted that this approach is designed to counteract excess eddy viscosity in free shear layer, which has addressed the "Grey Area" problem in cases that the free shear layer is away from the wall. However, it also eliminates the eddy viscosity at the near-wall region in this shallow-separation case, resulting in no turbulence appearances downstream the separation point.
(a) (b) Figure 12. SA-IDDES computations: Dynamic Mode Decomposition, Koopman spectrum, for the uncontrolled (a) and the controlled (b) cases.
(a) (b) Figure 13. SA-IDDES computations: iso-contours of spanwise-averaged pressure for the regular mode (a) and the shear-layer flapping mode (b).
Therefore, the SA-IDDES results have been selected in the subsequent analysis. The Dynamic Mode Decomposition (DMD) is applied to link a spatial structure (a dynamic mode) to a given frequency and to
18
compare with the spectral analysis performed above. The DMD results are based on a sequence of 600 spanwise-averaged pressure snapshots in the focus region (see Fig. 2b), with a constant sampling period ǻt = 0.075 c sinĮ /U. The non-dimensional frequency resolution of DMD is thus around St of 0.012. Fig. 12a depicts the eigenspectrum of the Koopman operator. As can be seen, the dynamics are affected primarily by the characteristic frequencies St = 0.60 and St =0.034. The former is characteristic of the vortex shedding, also called the regular mode18, while the latter can be recognized as the shear-layer flapping mode associated with Strouhal numbers on the order of 0.0219. Consequently, it is useful to observe the modes patterns in the three directions of space, and to spatially correlate them with the spectral analysis. In Fig. 13, a three-dimensional view with iso-surface of pressure normalized between Ѹ 0.001 and +0.001 and the corresponding top view in the spanwise direction in order to show the shape of both the vortex shedding and flapping modes. Overall, it can be concluded that the SA-IDDES results for the baseline setup are reliable enough to allow further numerical investigation with pulsed-blowing actuation for flow control. B. Controlled Flow On the basis of the uncontrolled flow, the pulsed blowing flow-control mechanism is applied. All controlled computations use the baseline flow solutions as initial flow conditions. The blowout angle, duty cycle and geometrical parameters are chosen according to the measurement9. The other actuation parameters, such as the jet momentum coefficient, Cȝ = 2 (Ahole /Awing)×(ua /U)2, and the reduced excitation frequency, F*= fex × c sinĮ /U, are varied. Table 2 provides a summary of the parameter variations. Previous investigations by the present authors15 led to the choice of these values. It is expected that with an optimal F*, at the two lower Cȝ values, the lift can be augmented and an increase between the two values can be observed. At the higher VR values, at some point the saturation in the gain in lift-to-drag ratio is expected. Moreover, the present computation can well reproduce the transient process in the controlled-flow experiment9 from original separated flow to a statistically attached flow, as shown in Fig.
19
14, where Cȝ = 0.67% and F* = 0.0017. It is seen that the jets are switched on at t = 0, and at about t = 0.1 s, the trailing-edge separation is eliminated with a reduced drag coefficient. Table 2. Overview of the investigated parameters for pulsed-blowing actuation.
F* 0.3, 0.6, 0.9, 1.2
Cȝ , % 0.0067, 0.0268, 0.107, 0.24, 0.67
Phase shift, deg I. 0, 180, 0, 180; II. 0, 90, 180, 270; III. 0, 270, 90, 180; IV. 0, 180, 90, 270;
Figure 14. Transient distribution of drag coefficient after the jets are switched on att = 0.
An overview of the lift-to-drag ratios for all actuated cases in comparison to the baseline flow case is given in Fig. 15. It is seen that for all F* values, saturation is yielded already at intermediate momentum coefficients of Cȝ = 0.107%; and the optimal F* is of 0.6 among all cases. Therefore, the results of the case with Cȝ = 0.107% and F* = 0.6, are selected to study the excitation mechanism in the following.
Figure 15. Mean lift-to-drag ratio for all actuated cases.
20
It is seen from Fig. 16 that three-dimensional vortex structures are formed at the orifices and the mild trailing-edge separation is dramatically suppressed. The time-averaged mean flow patterns are shown in Fig. 17, respectively for the baseline flow and for the controlled flow. Recall that with the jets switched off, it is verified that the separation occurs at xsep = 0.7 c by the IDDES-SA computation.
(a) 1/2T (b) T Figure 16. Instantaneous snapshot of the controlled flow.
Figure 17. Time-averaged mean flow patterns on the airfoil surface with (right) and without (left) flow control.
The mean surface pressure distribution is displayed in Fig. 18 in comparison with the measured data for both the baseline and the excited flows. Good agreement is obtained for both cases. In response to the pulsed blowing, higher Cp is predicted on the pressure side of the jet-controlled airfoil. Because of the increase of Cp on the pressure side and its decrease on the suction side in the excited Àow, it contributes to a larger lift coefficient than the baseline Àow.
21
Figure 18. Time- and spanwise-averaged surface pressure distribution with and without flow control.
Figure 19 compares the velocity profiles of the baseline and the controlled Àows at x/c = 0.971 and x/c = 1.98 with the available experimental results. The calculations provide reasonably good comparisons with the measurements for both the baseline and controlled Àows. The controlled Àow has augmented velocity pro¿les, to resist the adverse pressure gradient at the tailing edge. Besides, the wake becomes narrow in the controlled Àow case. For the turbulence statistics at the streamwise locations upstream the trailing edge, it is seen from Fig. 20 that although there are some gaps between calculations and measurements, the consistent distribution trends. For the controlled flows, the peak value of the Reynolds stresses decreases dramatically and their locations shifts closer to the wall due to the flow attachment.
Figure 19. Time-averaged streamwise velocity profiles (left) near the tailing edge, at x/c = 0.846, and (right) in the wake, at x/c = 1.98.
22
Figure 20. Time-averaged streamwise normal stress profiles at x/c = 0.84 (left) and 0.97 (right).
Due to the pulsed blowing actuation at a duty cycle, the jet orifices are all inactive at the same time period (see Fig. 16). Therefore, we conducted a phase shift among jet orifices, which has only a local effect on the evolution of vortices induced by jets, as shown in Fig. 21. A list of the phases of each orifice is given in Table 2. The results show the most effective way is that the orifices work anti-phased (case I): while one orifice is blowing out, the adjacent ones are inactive and vice-versa. The results of case I with F* = 0.6 and Cȝ = 0.107% are thus selected to study the excitation mechanism in the following.
Figure 21. The evolution of vortices induced by the anti-phased jets.
23
Figure 22 compares the dynamic behaviors of the anti-phased and the normal cases in the time-dependent progression of the lift. For clarity, the time signal of the pulse is only depicted for the un-shifted case. It is seen that using the pulsed blowing actuation in combination with a phase shift does not improve the lift. However, the R.M.S. value decreases from 0.048 to 0.021 for the phase-shifted case. The reduction of the fluctuation range by 56% could be more valuable than a slightly higher lift of ¨ CL = 4.2%. The signi¿cant changes in dynamics become obvious by means of frequency spectra shown in Fig. 23. It is seen that the frequency spectrum of the controlled flow (without phase shift) also exhibits a peak at St of 0.6, but due to the periodic perturbation, the peak is amplified compared to the uncontrolled flow (see Fig. 11a); for the phase-shifted case, the ¿rst dominant amplitude at the Strouhal number of the excitation is 1 order lower, and the next higher harmonics are all damped. It is noted that the outcome of a phase
shift also depends on the spacing between two neighboring jet orifices, which will be studied by the present authors to fully eliminate the harmonics.
Figure 22. Influence of a phase shift of 180o on the lift coefficient.
24
(a) (b) Figure 23. Spectra of the global lift coefficient for the controlled flow without (a) and with (b) phase shift.
Again, DMD is used in order to extract the main physical modes of the controlled flow. The eigen-spectrum of the Koopman operator is given in Fig. 12b. Recall that for the uncontrolled flow case (see Fig. 12a), the vortex shedding mode corresponds the most energy-containing mode at St of 0.60, while the shear-layer flapping mode is the second in energy rank at St of 0.034. It is seen that the vortex shedding mode is dramatically depressed with the scaled energy from 0.051 to 0.00076, while the scaled energy of the flapping mode decreases from 0.053 to 0.035. Figure 24 depicts the deconstructions of both modes.
Figure 24. Iso-contours of spanwise-averaged pressure for the regular mode (a) and the shear-layer flapping mode (b) of the controlled flow.
V. Conclusion A comprehensive numerical study of active flow control on a near-stall NACA0015 airfoil is conducted in this work. Using an in-house finite-volume based flow solver, detached-eddy simulations are carried out to evaluate the effectiveness of pulsed actuation mechanism to control the flow, delay separation, and increase the lift-to-drag ratio.
25
First of all, for the natural flow with mild separation, it is necessary to make an assessment of the five different DES-type models through the detailed investigation of the Reynolds stresses as well as the separation and reattachment. It is found that with the same mesh resolution, only the Spalart-Allmaras-based IDDES predicts flow separation and turbulent structures as demonstrated experimentally. The comparison between IDDES-SA and IDDES-SST results indicates that the underlying turbulence model is the determining factor for the prediction of the separation and the shear layer, because SST model provides less eddy viscosity in the attached flows, in consistent with the conclusion made by Mockett for several flow cases7. Both DDES methods with Smagorinsky SGS model (DDES-SA-SMG and DDES-SST-SMG) encounter the "Grey Area" problem that results in a delayed separation prediction, which was also reported by Durrani and Qin for the case of an A-airfoil with shallow flow separation20. Moreover, the application of WALE SGS model in LES mode for DDES proves to be inappropriate for such a mildly separated flow case, in which the free shear layer is close to the wall. Consequently, spectra analyses based on the IDDES-SA results are carried out by using Discrete Fourier Transform. Subsequently, Dynamic Mode Decomposition is applied to identify spatial structures to the different Strouhal frequencies. The classical vortex-shedding mode at St of 0.6, and the shear-layer flapping mode at St of 0.034, are both observed in the baseline flows, demonstrating again the correctness of the IDDES-SA computations. With better understanding of unsteady flow features, effective control practices are illustrated. It is observed from the IDDES-SA results that the elimination of the trailing-edge separation has been achieved by momentum injection into the boundary layer through the near-wall vortex motion induced by the jet flow. The suppression of separation bubble leads to a maximum lift-to-drag ratio enhancement of 194%, as well as a sharp decrease of all Reynolds stress components. The parameters are varied such that the non-dimensional excitation intensity (Cȝ) ranges from 0.0067% to 0.67% and the reduced frequency (F*) ranges from 0.3 to 1.2. It is found that for all F* values, the saturation of lift-to-drag ratio is yielded already at intermediate momentum coefficients of Cȝ = 0.107%, two-fifths of the actuation amplitude used 26
in the experiment10; the optimal F* is 0.6 among all cases, consistent with the dominant frequency of the baseline flow. This indicates that the trailing-edge separation is prevented not only by the increase of the circulation around the entire airfoil through the pulsed jet-injection, but also by the excitation of frequency-specific longitudinal vortices that form and keep the flowclose tothe surface due to increased mixing. A phase shift among jet orifices is investigated as a second parameter. The pulsed blowing with the jet orifices in an anti-phased manner is found to be very effective to damp the fluctuations of the lift around its mean value, though it leads to a small decrease in mean lift. Frequency spectra reveal that a significant reduction of the first dominant amplitude takes place. It suggests that for an effective use of pulsed actuation on such a near-stall airfoil, not only is an optimal choice of the ‘classical’ parameters (intensity and frequency) necessary but also the unsteady parameter has to be taken into account. Further investigations such as a variation of the duty cycle will be conducted. In all, the IDDES-SA computations have reproduced all major experimentally observed features of the present airfoil flow and its response to the perturbation provoked by the pulsed jet injected. The perturbation causes specific changes to the evolution of the vortex structures in the separated shear layer that requires well resolved. Therefore, RANS simulations cannot reproduce such a strong sensitivity to the actuation at certain frequencies21. Considering its relatively low computational cost, it is of great potential for AFC by using a well-calibrated HRLM in combination with optimization algorithms. Acknowledgments This work was supported by the National Key Basic Research Program of China (2014CB744801), the EU-China Joint Projects MARS (266326) and DRAGY (690623), the NSFC Grants 11572177, 51376106 & 11272183, and the Tsinghua University initiative Scientific Research Program (2014z21020). The authors also acknowledge technical discussion with Prof. Qibing Li at Tsinghua University. References 1.
Spalart P R, Allmaras S. A one-equation turbulence model for aerodynamic flows. Proc. 30th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, USA, 1992.
27
2.
Spalart P R, Jou W H, Strelets M, Allmaras, S R. Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach. Advances in DNS/LES, 1997, 1: 4-8.
3.
Spalart P R, Deck S, Shur M L, Squires, K D, Strelets M K, Travin A. A new version of detached-eddy simulation, resistant to ambiguous grid densities. Theoretical and computational fluid dynamics, 2006, 20(3): 181-195.
4.
Shur M L, Spalart P R, Strelets M K, Travin A K. A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capabilities. International Journal of Heat and Fluid Flow, 2008, 29(6): 1638-1649.
5.
Menter F R. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA journal, 1994, 32(8): 1598-1605.
6.
Travin A, Shur M, Strelets M M, Spalart P R. Physical and numerical upgrades in the detached-eddy simulation of complex turbulent flows. Advances in LES of complex flows. Springer Netherlands, 2002: 239-254.
7.
Mockett C, Fuchs M, Garbaruk A, et al. Two non-zonal approaches to accelerate RANS to LES transition of free shear layers in DES. Progress in Hybrid RANS-LES Modelling. Springer International Publishing, 2015: 187-201.
8.
Nicoud F, Ducros F. Subgrid-scale stress modelling based on the square of the velocity gradient tensor. Flow, Turbulence and Combustion, 1999, 62(3): 183-200.
9.
Siauw W L, Bonnet J P, Tensi J, Noack B R, Cattafesta L. Transient dynamics of the flow around a NACA 0015 airfoil using fluidic vortex generators. International Journal of Heat and Fluid Flow, 2010, 31(3): 450-459.
10. Siauw W L. Transient process of separation and attachment over a NACA 0015 airfoil controlled by fluidic
vortex generators. Ph.D. thesis, Poitiers, 2008. 11. Wang W, Siouris S, Qin N. Hybrid RANS/LES for active flow control. Aircraft Engineering & Aerospace
Technology, 2014, 86(3):179-187. 12. Mockett C. A comprehensive study of detached-eddy simulation. Ph.D. thesis, Berlin University of
Technoledge, 2009. 13. Knacke T. Potential effects of Rhie & Chow type interpolations in airframe noise simulations. Accurate and
efficient aeroacoustic prediction approaches for airframe noise, VKI lecture notes. Von Kármán Institute for Fluid Dynamics, Rhode-Saint-Genese, Belgium, 2013. 14. Wang L, Li L Y, Fu S. Detached-Eddy Simulation of flow past a pitching NACA 0015 airfoil with pulsed
actuation. In: 54th AIAA Aerospace Sciences Meeting, 2016: 1832. 15. Wang L, Li L Y, Fu S. Numerical investigation of active flow control on a pitching NACA 0015 airfoil using
Detached-Eddy Simulation. In: 52nd AIAA Aerospace Sciences Meeting, 2014:0765. 16. Spalart P R, Streett C. Young Person’s Guide to Detached-Eddy Simulation Grids. NASA CR 2001-211032,
2001. 17. Shur M L, Spalart P R, Strelets M K. LES-based evaluation of a microjet noise reduction concept in static and
flight conditions. Journal of Sound and Vibration, 2011, 330(17): 4083-4097. 18. Sigurdson L W. The structure and control of a turbulent reattaching flow. Journal of Fluid Mechanics, 1995,
298: 139-165.
28
19. Kiya M S, Shimizu M, Mochizuki O. Sinusoidal forcing of a turbulent separation bubble. Journal of Fluid
Mechanics,1997, 342:119-139. 20. Durrani N, Qin N. Behavior of Detached-Eddy Simulations for mild airfoil trailing-edge separation. Journal of
Aircraft, 2011,48:193-202. 21. Jakarliü, S., Jester-Zürker, R., Tropea, C., 2001. Proc. 9th ERCOFTAC/IAHR/COST Workshop on refined
turbulence modelling. 22. Gritskevich M S, Garbaruk A V, Schütze J, Menter, F R. Development of DDES and IDDES Formulations for
the k-Ȧ Shear Stress Transport Model. Flow, Turbulence and Combustion, 2012, 88(3): 431 431-449. 23. Wang L, Mockett C, Knacke T and Thiele F. Detached-Eddy Simulation of Landing-Gear Noise. In 19th
AIAA/CEAS Aeroacoustics Conference, 2013, Berlin, Germany, AIAA Paper-2013-2069. 24. Hu, R., Wang, L., and Fu, S., “Investigation of the Coherent Structures in Flow behind a Backward-facing
Step,” International Journal of Numerical Methods for Heat and Fluid Flow, Vol. 26, No. 3, 2016, pp. 1050-1068.
29
Liang Wang, Song Fu
PII: DOI: Reference:
S1270-9638(16)30881-1 http://dx.doi.org/10.1016/j.ast.2017.06.002 AESCTE 4058
To appear in:
Aerospace Science and Technology
Received date: Revised date: Accepted date:
19 October 2016 2 June 2017 3 June 2017
Please cite this article in press as: L. Wang, S. Fu, Detached-eddy simulation of flow past a pitching NACA 0015 airfoil with pulsed actuation, Aerosp. Sci. Technol. (2017), http://dx.doi.org/10.1016/j.ast.2017.06.002
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Detached-Eddy Simulation of flow past a pitching NACA 0015 airfoil with pulsed actuation Liang Wang1, [, and Song Fu1,] Tsinghua University, Beijing, China, 100084
1
This study conducted Detached-Eddy Simulations (DES) of actively controlled flow over a 0.35 m chord NACA0015 airfoil at an incidence angle of 11o and a chord Reynolds number of 1 million. The the uncontrolled, natural flow were simulated by using DES-type methods, with a comparison of different underlying RANS models and subgrid-scale stress models in LES mode. Results from these computations were compared with experimental observations, enabling their reliable assessment through the detailed investigation of the Reynolds stresses as well as the separation and reattachment. It was found that among all five DES-type methods, only the Spalart-Allmaras-based Improved Delayed DES (IDDES-SA) captured the separation point as measured in the experiments. The classical vortex-shedding and the shear-layer flapping modes for airfoil flows with shallow separation were also extracted from the IDDES-SA results by using Dynamic Mode Decomposition. With better understanding of unsteady flow features, effective control practices were illustrated. The pulsed blowing actuation was adopted according to the experimental settings. The flow-field around jet orifices was resolved, instead of the use of actuation boundary conditions. It was found that for all actuation frequencies investigated in this article, the saturation of lift-to-drag ratio was yielded already at intermediate momentum coefficients, two-fifths of the actuation amplitude used in the experiment; the optimal actuation frequency was 0.6 among all cases, consistent with the dominant frequency of the baseline flow. The suppression of separation bubble led to a maximum lift-to-drag ratio enhancement of 194%, as well as a sharp decrease of all Reynolds stress components. To
[ Assistant professor, School of Aerospace Engineering. ] Professor, School of Aerospace Engineering. Email: [email protected] 1
further reduce the strong unsteadiness introduced by excitation, a phase shift among jet orifices was applied and evaluated. The pulsed blowing with the jet orifices in an anti-phased manner was found to be very effective to damp the fluctuations of the lift around its mean value, though it led to a small decrease in mean lift.
Nomenclature Į
= angle of attack, deg
Ahole
= area of jet orifice, 7.85E-7 m2
Awing
= area of wing surface per orifice, 5.25E-3 m2
c
= chord length, m
CDES
= constant taken as 0.628 for the SA model and 0.550 for the SST model
CL
= lift coefficient
Cp
= pressure coefficient
Cȝ
= jet momentum coefficient
ǻt
= time step, s
DES
= Detached-Eddy Simulation
DDES
= Delayed DES
IDDES
= Improved DDES
fd
= blending function
fex
= frequency of periodic excitation, 1 / s
F*
= fex × c sinĮ /U, non-dimensional excitation frequency
HRLM
= Hybrid RANS/LES Method
k
= turbulent kinetic energy, m2 / s2
Rec
= Reynolds number based on c
R.M.S.
= Root Mean Square
SA
= Spalart-Allmaras
SMG
= Smagorinsky
SST
= Shear Stress Transport
St
= f c sinĮ/ U, Strouhal number
U
= freestream velocity, m / s
ua , ujet
= amplitude of the excitation velocity, excitation velocity, m / s
WALE
= Wall-Adapting Local Eddy-viscosity
xsep
= separation location, m
ȝ
= dynamic viscosity, ȝ = ȡȞ
Ȟ
= kinetic viscosity, m2 / s
2
Ȟt
= eddy viscosity, m2 / s
Ȧ
= specific turbulence dissipation rate, 1 / s
IJ
= t (U / c), time in convective units
Subscripts
= free-stream value
w
= wall
I. Introduction The determination of flow separation regions is a key task in the design process of wings and turbines, which beyond the calibration regime of the conventional RANS (Reynolds-Averaged Navier-Stokes) models that are not always consistent in their prediction capability. Though LES (Large-Eddy Simulation) is a viable alternative to RANS for accurate computations, the near-wall resolution required for LES makes it expensive. This is the rationale behind the continual development and application of HRLM (Hybrid RANS/LES Method), which combines advantages from both sides. In order to exhibit both LES and RANS behavior, a simple modification of the Spalart-Allmaras model1 was proposed by Spalart2 such that the computation behaves as LES in the unsteady "detached" regions and as RANS in the "attached" near-wall regions of the flow. This scheme, known as DES (Detached-Eddy Simulation), has been widely applied for a plethora of external flows since its inception. Following this development, two improved DES variants, DDES (Delayed DES)3 and IDDES (Improved DDES)4, have been further proposed by the coworkers of Spalart. This approach was also applied by Travin et al.6 and further by Gritskevich et al.22 to the RANS SST (Shear Stress Transport) model5. However, in DES-type methods, the turbulent boundary layers upstream of a free shear layer are typically in RANS mode and therefore do not contain any resolved turbulence. As a consequence, the development of turbulence in the shear layer will start from the intrinsic instability of the shear layer and will not be forced by turbulence coming from the boundary layers. This resembles the case for laminar boundary-layer flows upstream, where there is no sufficient turbulence due to slow shear-layer development, known as the "Grey-Area" problem. To address this problem, Mockett et al.7 used the WALE SGS (subgrid-scale stress) model8 to replace the conventional Smagorinsky model, resulting in a
3
dramatic acceleration of RANS/LES transition. Note that this approach retained generality, as demonstrated for several free-shear flow cases including a plane shear layer, a backwards-facing step and a round jet. Recently, Siauw et al.9, 10 experimentally studied the actively controlled flow over a 0.35 m chord NACA0015 airfoil at an incidence angle of 11o and a chord Reynolds number of 1 million. With the use of a pulsed-jet actuation, they successfully eliminated the trailing edge separation and therefore improve the aerodynamic performance of the airfoil. It was then selected as a standard test case by the EU-China project MARS (Manipulation of Reynolds Stress for Separation Control and Drag Reduction). The researchers within this project tried to numerically validate the experimental findings using HRLM. However, their results agree poorly with the experimental data for the uncontrolled case where the flow separation is mild11. Besides, in the experiment the pulsating frequency is of 1 Hz while the maximum jet velocity is of 200 m/s. The very low actuation frequency makes the pulsed jet similar to continuous jet blowing while the actuation intensity is comparably too high for real applications. The primary objective of this work is thus to numerically investigate the effects of the flow-control parameters, such as intensity and frequency, on the mildly separated flow on a NACA0015 airfoil, by using DES-type methods. The secondary objective is to assess the DES-type methods used for the baseline, natural flow case, with a comparison of different choices of underlying RANS model as well as SGS model in LES mode. This paper is organized as follows: Section 2 presents the mathematical formulations including DES-type method (DDES and IDDES variants) methods as well as alternative SGS models; Section 3 gives numerical details, especially for mesh generation and the modeling of excitation mechanisms; the numerical results of both baseline and controlled flows, including flow quantities and the subsequent spectra analyses, are presented in Section 4; the conclusions of this study are summarized in Section 5.
4
II. Mathematical Formulation In this section, one-equation SA (Spalart-Allmaras) model1 and k–Ȧ SST model5 are presented first to start the introduction of the DES type methods used here. Then, we express the DES method with DDES and IDDES variants. Finally, the alternative SGS (Smagorinsky and WALE) models in LES Mode for DDES are given. The details of the model formulations can be found in the literatures. A. The RANS Spalart-Allmaras Model and Menter’s SST Model The SA model and the SST model are the work-horses turbulence models of the aerospace industry. Both of them are eddy viscosity-based models. The SA model is based on a well-calibrated transport equation for a working quantity, as expressed as DUQ~ Dt
~ 1 ° w Cb1 US Q~ ® V ° wx j ¯
2 ~ 2 ª § wQ~ · ½° wQ~ º ~ ¸ ¾ C w Uf w §¨ Q ·¸ « P P » Cb 2 U ¨¨ ¸ wx j »¼ ©d ¹ «¬ © wx j ¹ °¿
(1)
~
where S is a function of the vorticity, Cb1, Cb2 and Cw the model constants and fw the empirical functions. The turbulent Reynolds stress tensor is then modelled by the eddy viscosity Pt with Pt
(2)
UQ f v
1
The SST model requires the solution of transport equations for k and Ȧ: wU k wk º w ª « U u j k P V k Pt » W S E *U kZ wt wx j »¼ ij ij wx j «¬
(3)
and w UZ wZ º w ª « Uu jZ P V Z Pt » wx j ¬« wt wx j »¼
PZ E kZ 2 2 1 F1
UV Z 2 wk wZ Z wx j wx j
(4)
where S is the mean strain rate, and F1 is a blending function. The eddy-viscosity ȝt is calculated from Pt
§ U k U a1k · min ¨ , ¸ © Z :F2 ¹
5
(5)
where a1=0.31 and F2 is also a blending function. The SA and SST model coefficients and empirical functions can be found in Ref.1 and Ref.5, respectively. B. IDDES Method It is designed in DES that the entire boundary layer is simulated by using RANS. The k equation in RANS, Eq. (3), can be rewritten as w U k w Uu j k wt wx j
w wx j
ª§ Pt «¨ P Vk ¬«©
· wk º U k 1.5 W S » ¸ ij ij LRANS ¹ wx j ¼»
(6)
where LRANS is the turbulence length scale for RANS. It equals to k0.5/Z for the SST model. IDDES4 is developed from DDES3 and it combines DDES and the wall-modeled LES (WMLES), with LRANS in Eq. (6) substituted by a length scale LIDDES
fd 1 fe LRANS 1 fd LLES
(7)
Here, the length scale for LES, LLES, is further written as the grid scale ' multiplied by a model constant CDES. fd and fe are both blending functions. It is seen from Eq. (7) that once fe equals to zero, LIDDES can be expressed as
LIDDES LDDES
fd LRANS 1 fd LLES
(8)
and IDDES reverts to DDES. Ref. 4 and Ref. 22 present the details of all blending functions. C. Alternative Model Formulations in LES Mode for DDES For DES-type models in LES mode, the following form can be derived under the assumption of local equilibrium (i.e. equality of the dissipation and generation terms of the underlying transport equations):
X SGS
Cm '
2
Dm (u )
(9)
where ǻ is the sub-grid characteristic length scale, and Dm (u) stands for a differential operator, acting on the resolved velocity field. For the Smagorinsky model
6
Dm
Ds
2Sij Sij , Cm Cs | 0.18ˈSij
1 wui wu j ( ) 2 wx j wxi
(10)
There are mainly two drawbacks in this model: one is that the differential operator Dm (u) doesn't vanish in near-wall regions, which gives a non-zero value of Ȟt as soon as a velocity gradient exists; from the formula, we can see it is only related to the strain rate of the turbulent structure but not to their rotation rate8. To solve these two problems, Mockett et al.7 adapted the WALE model8 that can be expressed as Dm
Sijd
Dw
( S ijd S ijd )3/2
(11)
( Sij Sij )5/2 ( S ijd S ijd )5/4
1 2 1 ( gij g 2ji ) G ij g kk2 2 3
, gij2
gik g kj
, gij
wui wx j
(12)
To modify DDES to behave like the WALE model in LES mode, we leave the length scale substitution, Eq. (8), unmodi¿ed and introduce an additional function to substitute the corresponding term for Dm (u) in the LES mode region only. The velocity gradient invariant in the underlying RANS model, S*W-DDES, is substituted by SW* DDES * S RANS
* * S RANS f d pos LRANS LLES S RANS BW SW*
2:ij :ij
, pos x Bw
0 , if x d 0 ® ¯1 , if x ! 0
CW2 / CS2
(13) (14) (15)
here Cw and Cs are model constants. Readers may refer to Ref.8 for all the coefficients. III. Numerical methods and grid generation A. Flow solver An implicit pressure-based solver with a fully conservative approximation of the governing equations is employed in the flow simulation. The code is based on curvilinear coordinates and uses cell-centered collocated storage arrangement on semi-block structured grids for all quantities. The one-equation SA
7
model1 and k–Ȧ SST model5 are used as the underlying models for the DDES3 and IDDES4presented here. The CDES parameter is set to 0.628 for the SA model and to 0.550 for the SST model, as calibrated against the decay of isotropic turbulence12. In order to best address the inconsistency in the demands posed by RANS and LES on the numerical scheme for the convective fluxes, the blending function of Travin et al.6 is employed to assure an appropriate switch between a higher-order TVD scheme and central differencing. The diffusive fluxes are approximated using a second order central scheme and for time discretization second order backward differencing is applied. The continuity equation is conserved by the SIMPLE algorithm whereby the decoupling of pressure and velocity is prevented through a modified Rhie & Chow interpolation13. The code is parallelized via domain decomposition and the data interchange between processors is realized through the standardized MPI-library. RANS solutions obtained with the SA or SST models are used as initial condition for the DDES and the IDDES. B. Flow configuration The numerical test model (see Fig. 1a) represents a NACA 0015 airfoil at a chord Reynolds number of about 1 million in response to the deployment and removal of actuation orifices positioned at 30% of the chord length c ( = 0.35 m) of the airfoil. The influence of the incidence angle on the separation was carefully analyzed in the experiments9, 10. It was found that at 11o of incidence the separation was two-dimensional with minimal vibration. The chosen flow condition is such that jet deployment corresponds to complete flow attachment over the airfoil; jet removal will cause flow separation up to an extent defined by the uncontrolled separated state. For the pulsed-blowing flow control, an array of 44 holes is deployed at 0.3c down-stream of the leading edge of the airfoil. This array of holes occupies the central one-third spanwise portion of the airfoil with the holes spaced apart in the span with an interval of 0.043 c. The jets are realized by blowing air into tubes connecting to the holes distributed on the airfoil surface. The tubes have a diameter of 0.001 m, and the jet Àow from the tubes has a pitch angle of 30o and a skew angle of 60o, as shown in Fig. 1b.
8
The Àow area at the jet exit on the airfoil surface is an ellipse with a major axis length of 0.002 m and a minor axis length of 0.001 m.
(a) (b) Figure 1. Schematics of the NACA0015 airfoil (a); of the pulsed-jet actuator position and orientation (b).
C. Computational Mesh In the previous work14, two 2-D grids of 0.10 million (see Fig. 2b) and 0.14 million nodes, respectively, were generated to test the grid independence. The 3-D grid was based on an expansion of the 2-D mesh slice into the spanwise direction. 101 layers of the 2-D mesh were combined. As a result of refining the grid, the turbulent viscosity is slightly lowered in the free shear layer. This indicates a more LES-like simulation than with baseline grid. However, it already provides a good representation of the turbulent structures on the baseline grid. It is also proven that the use of periodic boundaries is valid in case of an adequate spanwise domain size. The computational domain adopted is of size LC×LN×LZ = 17.1c×7.4c×0.18c and the total number of grid is 13.5 million nodes. The spanwise length is designed for four holes to be uniformly distributed. According to Spalart16, both the focus region (above the airfoil while downstream the separation point at x/c = 0.7) and the departure region (wake), shown in Fig. 2b in red, need to be resolved sufficiently to be in LES mode. ǻx+, ǻy+ and ǻz+ are no larger than 33.7 in the focus region. The non-dimensional wall distance of the first cell center remains below y+ =1 on the entire surface.
9
(a) block topology
(b) meshes around the airfoil
(c) zoom-in at trailing edge
(d) jet orifice configuration
(e) x-y plane mesh nearby jet orifice Figure 2. Computational mesh.
D. Numerical setup and parameters
RANS simulations are performed prior to the unsteady simulation using the same turbulence model as used for the DESs. All the DESs have been run for more than 20 convective time units, IJ = t U / c, before the evaluation of averages and the storage of the surface data was initiated. The sampling times are at least 30 convective time units for the DESs.
10
The prescribed non-dimensional time step size of ǻt = 2.5 ×10í2 c sinĮ /U guarantees CFLmax < 1.0 in the separation region17 as well as providing accurate resolution up to St = 2. The spectra will be displayed up to 1/3 of the Nyquist frequency, which is St = 20. Roughly 2 orders of magnitude residual reduction are reached within 10 outer iterations per time step. At the wind tunnel entry all flow quantities including the velocity components and turbulent properties are prescribed. The level of turbulence at the inflow is set to Tu = (2k/3) 0.5 / U = 0.5% and the turbulent viscosity ȝt / ȝ = 0.1. Previous investigations by the present authors23, 24 led to the choice of these values. At the outflow a convective boundary condition is used that allows unsteady flow structures to be transported outside the domain. The complete airfoil surface is modeled as a non-slip boundary condition. The wind tunnel walls are approximated as inviscid walls. Periodic boundary conditions are set in the spanwise boundaries. The resolved actuation orifices (see Fig. 2d) are considered in the present paper. The actuation velocity inlet is placed at its bottom. The time-dependent velocity for the pulsed actuation (see Fig. 3) is modeled as follows:
W W DC T 1 1 1 1 u jet ( x, t ) ua ( x) [ tanh( p )] [ tanh( p )] (20) 2 2 ts 99 2 2 ts 99 with
t pU f c sin D , t s 99 f sc , DC ,T , W p * F Uf c sin D (21) where ua is the amplitude velocity of the time-dependent excitation, F* the non-dimensional actuation T
frequency, T the cycle duration of the pulse, ts99 the time the pulsed jet achieving its maximum value, sc the smoothing parameter of the flank of the pulse, DC the duty cycle of the pulse, IJp the dimensionless pulse time given in convective units, and tp the pulse time. tp is reset after each cycle duration so that a periodic behavior is generated. Note that the jet exits are set as solid wall surfaces in the baseline computations.
11
Figure 3. Modeling of the time-dependent part of the pulsed actuation.
IV. Results and discussions In the following, it first presents the results of the baseline, natural flow, with the assessment of five different DES-type models through the detailed investigation of the Reynolds stresses as well as the separation and reattachment. Second, the results of the controlled simulations are given. These will be divided into results of the investigation of the excitation intensity as well as frequency, and results of a phase shift among jet orifices. A. Baseline Flow The various results reported herein correspond to the five different models employed during the study. Results of the DDES and IDDES computations are labeled DDES and IDDES, respectively. Then, results labeled SA and SST respectively correspond to those predictions in which the one-equation SA model and k–Ȧ SST model are used as the underlying RANS models. Finally, the results labeled by SMG have been obtained using the standard SGS model (Smagorinsky model) while by WALE using WALE SGS model in LES mode. For instance, DDES-SA-WALE represents the DDES computations with the SA RANS model as well as the WALE SGS model.
12
(a) DDES-SA-SMG
(b) DDES-SA-WALE
(c) DDES-SST-SMG
(d) IDDES-SST
(e) IDDES-SA
Figure 4. Instantaneous flow using the Ȝ2 vortex core criterion at Ȝ2 = - 500 s-2, shaded by dimensionless velocity magnitude.
Fig. 4 compares the instantaneous turbulence structures (iso-surface of the second invariant of velocity gradient). It is seen that the grid resolution remains sufficiently fine enough to resolve the intrinsic small-scale turbulent structures in the shear layer for all five models. The calculated separation locations from the five models are given in Table 1, which are de¿ned using the zero-Cf criterion in the distribution of the streamwise component of skin friction, Cfx, as shown in Fig. 5. The measured separation location (estimated by oil-Àow visualization), xsep, at 68% of the chord length is fairly well predicted
by
both
IDDES-SA
and
DDES-SA-WALE,
whereas
both
DDES-SA-SMG
and
DDES-SST-SMG models predict a delayed separation compared to experimental results, and therefore there is little turbulence generation and flow structures shown in Fig. 4a and Fig. 4c. On the contrary,
13
IDDES-SST gives an over-prediction in the bubble length (see Fig. 5) as well as in the turbulence structures (see Fig. 4d).
Figure 5. Time- and spanwise-averaged skin-friction distribution along the airfoil. Table 1. Comparisons of separation location between computations and measurements.
Methods x
sep
DDES-SA-SMG
DDES-SA-WALE
DDES-SST-SMG
IDDES-SST
IDDES-SA
Measurements
0.849c
0.676c
0.886c
0.573c
0.649c
0.687c
Figure 6. Time- and spanwise-averaged surface pressure distribution along the airfoil.
14
The pressure distribution contributes the most to the lift enhancement. The time-averaged surface pressure plot Cp in Fig. 6, including zoomed views near the trailing edge, indicates that IDDES-SA gives accurate prediction in the trailing-edge separated-flow region that is simulated with LES. The overshoot of both DDES-SA-SMG and DDES-SST-SMG in Cp at trailing edge can be attributed to their over-prediction in xsep and vice-versa for the IDDES-SST results. Moreover, the initial dip in all DES results in the leading-edge section can be attributed to the tripping applied at 0.4% c from the leading edge in the experiment.
Figure 7. Time- and spanwise-averaged streamwise velocity profiles at x/c= 0.84 (left) and 0.97 (right).
Figure 8. Time- and spanwise-averaged streamwise normal stress < u'u'> profiles at x/c= 0.84 (left) and 0.97 (right).
15
Fig. 7 and Fig. 8 compare the time-averaged streamwise velocity profiles and streamwise normal stress
intensive reverse flow in the separation bubble, while both DDES-SA-SMG and DDES-SST-SMG show no negative sign in velocity profiles and pretty weak turbulence intensity. IDDES-SA and DDES-SA-WALE give similar velocity profiles, whereas their discrepancy is huge in Reynolds stresses. Although there are some gaps between the results by IDDES-SA and measurements, the consistent distribution trends is obtained. On the contrary, turbulence intensity predicted by DDES-SA-WALE is close to zero.
Figure 9. Comparison of streamwise velocity (left) and shear stress
Fig. 9 displays the time-averaged streamwise velocity and shear stress
16
(a)
(b)
Figure 10. Sampling time series of the global lift coefficient for (a) IDDES-SA and (b) DDES-SA-WALE. IJ = t (U/ c), time in convective units.
(a)
(b)
Figure 11. Spectra of the global lift coefficient for (a) IDDES-SA and (b) DDES-SA-WALE.
To further compare the results of IDDES-SA and of DDES-SA-WALE, spectra analyses on the unsteady flow quantities are carried out by using Discrete Fourier Transform (DFT). Fig. 11 presents the power spectra density of the global lift coefficient as a function of the Strouhal number, St = f csinĮ /U, based on the projected airfoil height csinĮ. It should be noted that the smallest frequency (or the St number in Fig. 11) is related to the total time of the sampling (see Fig. 10), and the largest frequency is determined by the sampling frequency, which equals to 1/ǻt = 40 U / (csinĮ). It is observed from Fig. 10 that the R.M.S. value of global lift coefficient sharply decreases from 0.042 for IDDES-SA to 0.0063 for DDES-SA-WALE results, and from Fig. 11 that the peak value in the DDES-SA-WALE spectra is extraordinarily intensive compared to that in IDDES-SA results. The high-frequency fluctuations of
17
DDES-SA-WALE results (see Fig. 11b) may be regarded as the spurious numerical noises that cannot be dissipated with the low level of eddy viscosity in the separation bubble (see Fig. 8b). Moreover, recall that DDES-SA-WALE predicts accurate separation location however near-zero Reynolds stress in the separation bubble (see Fig. 8). We can conclude that the separated flow simulated by DDES-SA-WALE is actually laminar. It is noted that this approach is designed to counteract excess eddy viscosity in free shear layer, which has addressed the "Grey Area" problem in cases that the free shear layer is away from the wall. However, it also eliminates the eddy viscosity at the near-wall region in this shallow-separation case, resulting in no turbulence appearances downstream the separation point.
(a) (b) Figure 12. SA-IDDES computations: Dynamic Mode Decomposition, Koopman spectrum, for the uncontrolled (a) and the controlled (b) cases.
(a) (b) Figure 13. SA-IDDES computations: iso-contours of spanwise-averaged pressure for the regular mode (a) and the shear-layer flapping mode (b).
Therefore, the SA-IDDES results have been selected in the subsequent analysis. The Dynamic Mode Decomposition (DMD) is applied to link a spatial structure (a dynamic mode) to a given frequency and to
18
compare with the spectral analysis performed above. The DMD results are based on a sequence of 600 spanwise-averaged pressure snapshots in the focus region (see Fig. 2b), with a constant sampling period ǻt = 0.075 c sinĮ /U. The non-dimensional frequency resolution of DMD is thus around St of 0.012. Fig. 12a depicts the eigenspectrum of the Koopman operator. As can be seen, the dynamics are affected primarily by the characteristic frequencies St = 0.60 and St =0.034. The former is characteristic of the vortex shedding, also called the regular mode18, while the latter can be recognized as the shear-layer flapping mode associated with Strouhal numbers on the order of 0.0219. Consequently, it is useful to observe the modes patterns in the three directions of space, and to spatially correlate them with the spectral analysis. In Fig. 13, a three-dimensional view with iso-surface of pressure normalized between Ѹ 0.001 and +0.001 and the corresponding top view in the spanwise direction in order to show the shape of both the vortex shedding and flapping modes. Overall, it can be concluded that the SA-IDDES results for the baseline setup are reliable enough to allow further numerical investigation with pulsed-blowing actuation for flow control. B. Controlled Flow On the basis of the uncontrolled flow, the pulsed blowing flow-control mechanism is applied. All controlled computations use the baseline flow solutions as initial flow conditions. The blowout angle, duty cycle and geometrical parameters are chosen according to the measurement9. The other actuation parameters, such as the jet momentum coefficient, Cȝ = 2 (Ahole /Awing)×(ua /U)2, and the reduced excitation frequency, F*= fex × c sinĮ /U, are varied. Table 2 provides a summary of the parameter variations. Previous investigations by the present authors15 led to the choice of these values. It is expected that with an optimal F*, at the two lower Cȝ values, the lift can be augmented and an increase between the two values can be observed. At the higher VR values, at some point the saturation in the gain in lift-to-drag ratio is expected. Moreover, the present computation can well reproduce the transient process in the controlled-flow experiment9 from original separated flow to a statistically attached flow, as shown in Fig.
19
14, where Cȝ = 0.67% and F* = 0.0017. It is seen that the jets are switched on at t = 0, and at about t = 0.1 s, the trailing-edge separation is eliminated with a reduced drag coefficient. Table 2. Overview of the investigated parameters for pulsed-blowing actuation.
F* 0.3, 0.6, 0.9, 1.2
Cȝ , % 0.0067, 0.0268, 0.107, 0.24, 0.67
Phase shift, deg I. 0, 180, 0, 180; II. 0, 90, 180, 270; III. 0, 270, 90, 180; IV. 0, 180, 90, 270;
Figure 14. Transient distribution of drag coefficient after the jets are switched on att = 0.
An overview of the lift-to-drag ratios for all actuated cases in comparison to the baseline flow case is given in Fig. 15. It is seen that for all F* values, saturation is yielded already at intermediate momentum coefficients of Cȝ = 0.107%; and the optimal F* is of 0.6 among all cases. Therefore, the results of the case with Cȝ = 0.107% and F* = 0.6, are selected to study the excitation mechanism in the following.
Figure 15. Mean lift-to-drag ratio for all actuated cases.
20
It is seen from Fig. 16 that three-dimensional vortex structures are formed at the orifices and the mild trailing-edge separation is dramatically suppressed. The time-averaged mean flow patterns are shown in Fig. 17, respectively for the baseline flow and for the controlled flow. Recall that with the jets switched off, it is verified that the separation occurs at xsep = 0.7 c by the IDDES-SA computation.
(a) 1/2T (b) T Figure 16. Instantaneous snapshot of the controlled flow.
Figure 17. Time-averaged mean flow patterns on the airfoil surface with (right) and without (left) flow control.
The mean surface pressure distribution is displayed in Fig. 18 in comparison with the measured data for both the baseline and the excited flows. Good agreement is obtained for both cases. In response to the pulsed blowing, higher Cp is predicted on the pressure side of the jet-controlled airfoil. Because of the increase of Cp on the pressure side and its decrease on the suction side in the excited Àow, it contributes to a larger lift coefficient than the baseline Àow.
21
Figure 18. Time- and spanwise-averaged surface pressure distribution with and without flow control.
Figure 19 compares the velocity profiles of the baseline and the controlled Àows at x/c = 0.971 and x/c = 1.98 with the available experimental results. The calculations provide reasonably good comparisons with the measurements for both the baseline and controlled Àows. The controlled Àow has augmented velocity pro¿les, to resist the adverse pressure gradient at the tailing edge. Besides, the wake becomes narrow in the controlled Àow case. For the turbulence statistics at the streamwise locations upstream the trailing edge, it is seen from Fig. 20 that although there are some gaps between calculations and measurements, the consistent distribution trends. For the controlled flows, the peak value of the Reynolds stresses decreases dramatically and their locations shifts closer to the wall due to the flow attachment.
Figure 19. Time-averaged streamwise velocity profiles (left) near the tailing edge, at x/c = 0.846, and (right) in the wake, at x/c = 1.98.
22
Figure 20. Time-averaged streamwise normal stress
Due to the pulsed blowing actuation at a duty cycle, the jet orifices are all inactive at the same time period (see Fig. 16). Therefore, we conducted a phase shift among jet orifices, which has only a local effect on the evolution of vortices induced by jets, as shown in Fig. 21. A list of the phases of each orifice is given in Table 2. The results show the most effective way is that the orifices work anti-phased (case I): while one orifice is blowing out, the adjacent ones are inactive and vice-versa. The results of case I with F* = 0.6 and Cȝ = 0.107% are thus selected to study the excitation mechanism in the following.
Figure 21. The evolution of vortices induced by the anti-phased jets.
23
Figure 22 compares the dynamic behaviors of the anti-phased and the normal cases in the time-dependent progression of the lift. For clarity, the time signal of the pulse is only depicted for the un-shifted case. It is seen that using the pulsed blowing actuation in combination with a phase shift does not improve the lift. However, the R.M.S. value decreases from 0.048 to 0.021 for the phase-shifted case. The reduction of the fluctuation range by 56% could be more valuable than a slightly higher lift of ¨ CL = 4.2%. The signi¿cant changes in dynamics become obvious by means of frequency spectra shown in Fig. 23. It is seen that the frequency spectrum of the controlled flow (without phase shift) also exhibits a peak at St of 0.6, but due to the periodic perturbation, the peak is amplified compared to the uncontrolled flow (see Fig. 11a); for the phase-shifted case, the ¿rst dominant amplitude at the Strouhal number of the excitation is 1 order lower, and the next higher harmonics are all damped. It is noted that the outcome of a phase
shift also depends on the spacing between two neighboring jet orifices, which will be studied by the present authors to fully eliminate the harmonics.
Figure 22. Influence of a phase shift of 180o on the lift coefficient.
24
(a) (b) Figure 23. Spectra of the global lift coefficient for the controlled flow without (a) and with (b) phase shift.
Again, DMD is used in order to extract the main physical modes of the controlled flow. The eigen-spectrum of the Koopman operator is given in Fig. 12b. Recall that for the uncontrolled flow case (see Fig. 12a), the vortex shedding mode corresponds the most energy-containing mode at St of 0.60, while the shear-layer flapping mode is the second in energy rank at St of 0.034. It is seen that the vortex shedding mode is dramatically depressed with the scaled energy from 0.051 to 0.00076, while the scaled energy of the flapping mode decreases from 0.053 to 0.035. Figure 24 depicts the deconstructions of both modes.
Figure 24. Iso-contours of spanwise-averaged pressure for the regular mode (a) and the shear-layer flapping mode (b) of the controlled flow.
V. Conclusion A comprehensive numerical study of active flow control on a near-stall NACA0015 airfoil is conducted in this work. Using an in-house finite-volume based flow solver, detached-eddy simulations are carried out to evaluate the effectiveness of pulsed actuation mechanism to control the flow, delay separation, and increase the lift-to-drag ratio.
25
First of all, for the natural flow with mild separation, it is necessary to make an assessment of the five different DES-type models through the detailed investigation of the Reynolds stresses as well as the separation and reattachment. It is found that with the same mesh resolution, only the Spalart-Allmaras-based IDDES predicts flow separation and turbulent structures as demonstrated experimentally. The comparison between IDDES-SA and IDDES-SST results indicates that the underlying turbulence model is the determining factor for the prediction of the separation and the shear layer, because SST model provides less eddy viscosity in the attached flows, in consistent with the conclusion made by Mockett for several flow cases7. Both DDES methods with Smagorinsky SGS model (DDES-SA-SMG and DDES-SST-SMG) encounter the "Grey Area" problem that results in a delayed separation prediction, which was also reported by Durrani and Qin for the case of an A-airfoil with shallow flow separation20. Moreover, the application of WALE SGS model in LES mode for DDES proves to be inappropriate for such a mildly separated flow case, in which the free shear layer is close to the wall. Consequently, spectra analyses based on the IDDES-SA results are carried out by using Discrete Fourier Transform. Subsequently, Dynamic Mode Decomposition is applied to identify spatial structures to the different Strouhal frequencies. The classical vortex-shedding mode at St of 0.6, and the shear-layer flapping mode at St of 0.034, are both observed in the baseline flows, demonstrating again the correctness of the IDDES-SA computations. With better understanding of unsteady flow features, effective control practices are illustrated. It is observed from the IDDES-SA results that the elimination of the trailing-edge separation has been achieved by momentum injection into the boundary layer through the near-wall vortex motion induced by the jet flow. The suppression of separation bubble leads to a maximum lift-to-drag ratio enhancement of 194%, as well as a sharp decrease of all Reynolds stress components. The parameters are varied such that the non-dimensional excitation intensity (Cȝ) ranges from 0.0067% to 0.67% and the reduced frequency (F*) ranges from 0.3 to 1.2. It is found that for all F* values, the saturation of lift-to-drag ratio is yielded already at intermediate momentum coefficients of Cȝ = 0.107%, two-fifths of the actuation amplitude used 26
in the experiment10; the optimal F* is 0.6 among all cases, consistent with the dominant frequency of the baseline flow. This indicates that the trailing-edge separation is prevented not only by the increase of the circulation around the entire airfoil through the pulsed jet-injection, but also by the excitation of frequency-specific longitudinal vortices that form and keep the flowclose tothe surface due to increased mixing. A phase shift among jet orifices is investigated as a second parameter. The pulsed blowing with the jet orifices in an anti-phased manner is found to be very effective to damp the fluctuations of the lift around its mean value, though it leads to a small decrease in mean lift. Frequency spectra reveal that a significant reduction of the first dominant amplitude takes place. It suggests that for an effective use of pulsed actuation on such a near-stall airfoil, not only is an optimal choice of the ‘classical’ parameters (intensity and frequency) necessary but also the unsteady parameter has to be taken into account. Further investigations such as a variation of the duty cycle will be conducted. In all, the IDDES-SA computations have reproduced all major experimentally observed features of the present airfoil flow and its response to the perturbation provoked by the pulsed jet injected. The perturbation causes specific changes to the evolution of the vortex structures in the separated shear layer that requires well resolved. Therefore, RANS simulations cannot reproduce such a strong sensitivity to the actuation at certain frequencies21. Considering its relatively low computational cost, it is of great potential for AFC by using a well-calibrated HRLM in combination with optimization algorithms. Acknowledgments This work was supported by the National Key Basic Research Program of China (2014CB744801), the EU-China Joint Projects MARS (266326) and DRAGY (690623), the NSFC Grants 11572177, 51376106 & 11272183, and the Tsinghua University initiative Scientific Research Program (2014z21020). The authors also acknowledge technical discussion with Prof. Qibing Li at Tsinghua University. References 1.
Spalart P R, Allmaras S. A one-equation turbulence model for aerodynamic flows. Proc. 30th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, USA, 1992.
27
2.
Spalart P R, Jou W H, Strelets M, Allmaras, S R. Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach. Advances in DNS/LES, 1997, 1: 4-8.
3.
Spalart P R, Deck S, Shur M L, Squires, K D, Strelets M K, Travin A. A new version of detached-eddy simulation, resistant to ambiguous grid densities. Theoretical and computational fluid dynamics, 2006, 20(3): 181-195.
4.
Shur M L, Spalart P R, Strelets M K, Travin A K. A hybrid RANS-LES approach with delayed-DES and wall-modelled LES capabilities. International Journal of Heat and Fluid Flow, 2008, 29(6): 1638-1649.
5.
Menter F R. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA journal, 1994, 32(8): 1598-1605.
6.
Travin A, Shur M, Strelets M M, Spalart P R. Physical and numerical upgrades in the detached-eddy simulation of complex turbulent flows. Advances in LES of complex flows. Springer Netherlands, 2002: 239-254.
7.
Mockett C, Fuchs M, Garbaruk A, et al. Two non-zonal approaches to accelerate RANS to LES transition of free shear layers in DES. Progress in Hybrid RANS-LES Modelling. Springer International Publishing, 2015: 187-201.
8.
Nicoud F, Ducros F. Subgrid-scale stress modelling based on the square of the velocity gradient tensor. Flow, Turbulence and Combustion, 1999, 62(3): 183-200.
9.
Siauw W L, Bonnet J P, Tensi J, Noack B R, Cattafesta L. Transient dynamics of the flow around a NACA 0015 airfoil using fluidic vortex generators. International Journal of Heat and Fluid Flow, 2010, 31(3): 450-459.
10. Siauw W L. Transient process of separation and attachment over a NACA 0015 airfoil controlled by fluidic
vortex generators. Ph.D. thesis, Poitiers, 2008. 11. Wang W, Siouris S, Qin N. Hybrid RANS/LES for active flow control. Aircraft Engineering & Aerospace
Technology, 2014, 86(3):179-187. 12. Mockett C. A comprehensive study of detached-eddy simulation. Ph.D. thesis, Berlin University of
Technoledge, 2009. 13. Knacke T. Potential effects of Rhie & Chow type interpolations in airframe noise simulations. Accurate and
efficient aeroacoustic prediction approaches for airframe noise, VKI lecture notes. Von Kármán Institute for Fluid Dynamics, Rhode-Saint-Genese, Belgium, 2013. 14. Wang L, Li L Y, Fu S. Detached-Eddy Simulation of flow past a pitching NACA 0015 airfoil with pulsed
actuation. In: 54th AIAA Aerospace Sciences Meeting, 2016: 1832. 15. Wang L, Li L Y, Fu S. Numerical investigation of active flow control on a pitching NACA 0015 airfoil using
Detached-Eddy Simulation. In: 52nd AIAA Aerospace Sciences Meeting, 2014:0765. 16. Spalart P R, Streett C. Young Person’s Guide to Detached-Eddy Simulation Grids. NASA CR 2001-211032,
2001. 17. Shur M L, Spalart P R, Strelets M K. LES-based evaluation of a microjet noise reduction concept in static and
flight conditions. Journal of Sound and Vibration, 2011, 330(17): 4083-4097. 18. Sigurdson L W. The structure and control of a turbulent reattaching flow. Journal of Fluid Mechanics, 1995,
298: 139-165.
28
19. Kiya M S, Shimizu M, Mochizuki O. Sinusoidal forcing of a turbulent separation bubble. Journal of Fluid
Mechanics,1997, 342:119-139. 20. Durrani N, Qin N. Behavior of Detached-Eddy Simulations for mild airfoil trailing-edge separation. Journal of
Aircraft, 2011,48:193-202. 21. Jakarliü, S., Jester-Zürker, R., Tropea, C., 2001. Proc. 9th ERCOFTAC/IAHR/COST Workshop on refined
turbulence modelling. 22. Gritskevich M S, Garbaruk A V, Schütze J, Menter, F R. Development of DDES and IDDES Formulations for
the k-Ȧ Shear Stress Transport Model. Flow, Turbulence and Combustion, 2012, 88(3): 431 431-449. 23. Wang L, Mockett C, Knacke T and Thiele F. Detached-Eddy Simulation of Landing-Gear Noise. In 19th
AIAA/CEAS Aeroacoustics Conference, 2013, Berlin, Germany, AIAA Paper-2013-2069. 24. Hu, R., Wang, L., and Fu, S., “Investigation of the Coherent Structures in Flow behind a Backward-facing
Step,” International Journal of Numerical Methods for Heat and Fluid Flow, Vol. 26, No. 3, 2016, pp. 1050-1068.
29