- Email: [email protected]
Dynamic hybridization of MILES and RANS for predicting airfoil stall
Available online at www.sciencedirect.com
Computers & Fluids 37 (2008) 161–169 www.elsevier.com/locate/compfluid
Dynamic hybridization of MILES and RANS for predicting airfoil stall Ichiro Nakamori *, Toshiaki Ikohagi
1
Institute of Fluid Science, Tohoku University, 2-1-1, Katahira, Aoba-ku, Sendai, Miyagi 980-8577, Japan Received 29 September 2005; received in revised form 11 September 2006; accepted 30 June 2007 Available online 30 August 2007
Abstract Hybridization comprised of an algebraic turbulence model based on the Reynolds average Navier–Stokes (RANS) equations with a monotonically integrated large eddy simulation (MILES) is proposed to simulate transient fluid motion during separation and vortex shedding at high Reynolds numbers. The proposed hybridization utilizes the Baldwin–Lomax model with the Degani–Schiff modification as the RANS model in the near-wall region and a MILES far from the wall. Although the hybridizationis assumed to be a MILES with wall modeling, the transition line between the RANS and the MILES modes is determined by the turbulent intensity that is dominated by the large eddies in the grid-scale. This hybrid model is applied to the flows past three different types of airfoils (NACA633-018, NACA631-012 and NACA64A-006) near stall, at a chord Reynolds number of Re = 5.8 · 106. These airfoils are classified as trailingedge-stall, leading-edge-stall and thin-airfoil-stall airfoils, respectively. The computed results are compared with wind tunnel experiments. The hybrid model successfully demonstrates accurate stall angle and surface pressure distribution predictions near the stall for each type of airfoil. The airfoil simulation results confirmed that the hybrid model provides a better prediction than the RANS model for unsteady turbulent flows with separation and vortex shedding simulations. 2007 Elsevier Ltd. All rights reserved.
1. Introduction Recently, CFD has become the preferred method for predicting aerodynamic performance during aircraft design. In particular, aerodynamic analysis at stall with massive separation is especially important when evaluating an aircraft’s ability to take off and land. However, current Reynolds average Navier–Stokes (RANS) models cannot accurately predict turbulent flow with vortex shedding or massive separation, because it is difficult to develop a RANS model that is suitable for a turbulent flow with a high degree of unsteadiness. Large eddy simulations (LES), however, could prove useful in predicting unsteady aerodynamic characteristics, because they are independent of the Reynolds average, thereby allowing for the resolution of highly unsteady fluid motion. Unfortunately, *
1
Corresponding author. Tel.: +81 22 217 5228; fax: +81 22 217 5229. E-mail address: [email protected] (I. Nakamori). Tel.: +81 22 217 5228; fax: +81 22 217 5229.
0045-7930/$ - see front matter 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2007.06.004
supercomputers are not currently able to resolve all eddies, because eddies near walls decrease in size with an increase in the Reynolds number. This problem could be avoided by employing a wall function for a specific turbulent flow. However, the wall function generally cannot be employed when the boundary layer includes the laminar to turbulence transition point. Another feasible option for bypassing the enormous CPU and storage requirement and resolving the wall turbulence is hybridization. This involves conducting a RANS near the wall and an LES for the detached eddies, as conducted in detached eddy simulation (DES) [4,5,14]. In DES, the eddy viscosity value switches between the Spalart–Allmaras (SA) model [15] and the LES mode using one eddy viscosity equation. The designated point is determined by the grid-scale. Because the DES switches from RANS to LES at a designated point away from the wall, the glitch point in the velocity profile occurs near the switching point. An example of this characteristic is demonstrated in Nikitin et al. [12]. In addition, the subgrid-scale (SGS) model in LES
162
I. Nakamori, T. Ikohagi / Computers & Fluids 37 (2008) 161–169
is effective for prediction accuracy analysis only when sufficient grid resolution is used in all directions: the streamwise direction, the spanwise direction, and the normal direction to the wall. However, the SGS model is inappropriate for practical grid resolution, such that the aspect ratio of the streamwise grid size and the grid size in the normal direction to the wall (or the aspect ratio of the spanwise grid size and the grid size in the normal direction to the wall) is extremely high near the wall. In the present study, in order to reduce the glitch point in the velocity profile, the switching point at which the turbulent intensity of the resolved detached eddy reaches a certain fixed level is noted. Further, a monotonically integrated large eddy simulation (MILES) [2] is applied as an LES after the switching point far from the wall. This MILES uses inherent numerical dissipation to damp out the SGS turbulence. The MILES leads to the DNS when a finer grid system is used. To save CPU time, the Baldwin–Lomax [1] algebraic model with the modification by Degani and Schiff [7] (BLDS model) is employed as a RANS model in conjunction with the MILES. The concept of the Degani–Schiff modification coincides with the directivity of the present hybrid model from the point of view that the excessive eddy viscosity of the RANS mode is decreased in the separated vortical flow. The objective of this study is to obtain accurate prediction of stall with massive separation by using less computational time in comparison with large eddy simulations.The hybrid model is then applied to flows past various types of airfoils near stall, in order to determine its applicability for predicting aerodynamic mean values and unsteady characteristics. The experimental aerodynamic characteristics of airfoils have previously been described by McCullough and Gault [11]. Each airfoil is categorized by its stalling characteristics as a trailing-edge-stall airfoil, a leading-edge-stall airfoil or a thin-airfoil-stall airfoil. Because it is challenging to predict the necessary angle of attack to produce a stall by CFD, the improvement of prediction accuracy of the stall angle of attack is very significant.
2. Numerical method 2.1. Governing equations and numerical scheme The Navier–Stokes equations are solved in an inertial system. The primitive variables in the solution vector for the RANS mode use the Reynolds average and those in the solution vector for the LES mode use the spatial average. The curvilinear coordinates are used to transform the equations: oQ oE oF oG 1 oE v oF v oG v þ þ þ ¼ þ þ : ot on og of Re on og of
ð1Þ
The transformed Q and flux vectors are given by 1 Q ¼ ðq; qu; qv; qw; eÞ J ^ b þ ny F b þ nz G nx E E¼ ; J ^ b þ gy F b þ gz G gx E ; F¼ J ^ b þ fy F b þ fz G fx E G¼ ; J where 2
qu
3
^ b þ ny F b þ nz G nx E ; J ^ b þ gy F b þ gz G gx E Fv ¼ ; J ^ b þ fy F b þ fz G fx E Gv ¼ ; J Ev ¼
2
qv
3
7 7 6 2 6 6 qu þ p 7 6 quv 7 7 7 6 6 2 7 7 6 E¼6 6 quv 7; F ¼ 6 qv þ p 7; 7 7 6 6 4 quw 5 4 qvw 5 ðe þ pÞu ðe þ pÞv 3 2 0 7 6 sxx 7 6 7 6 7; 6 sxy Ev ¼ 6 7 7 6 sxz 5 4 usxx þ vsxy þ wsxz 3 2 0 7 6 sxy 7 6 7 6 7; s Fv ¼ 6 yy 7 6 7 6 syz 5 4
2
ð2Þ
qw
3
7 6 6 quw 7 7 6 7 G ¼6 6 qvw 7; 7 6 4 qw2 þ p 5 ðe þ pÞw
usxy þ vsyy þ wsyz 3 0 7 6 sxz 7 6 7 6 7: 6 syz Gv ¼ 6 7 7 6 szz 5 4 usxz þ vsyz þ wszz 2
ð3Þ In these equations, J is the Jacobian of the coordinate transformation, q is the density, u, v, and w are the fluid velocity components, e is the total energy per unit volume, p is the pressure, and qi are the heat fluxes. The viscous stress tensors sij can be expressed as oui ouj 2 ouk þ dij sij ¼ ðl þ lt Þ ; ð4Þ oxj oxi 3 oxk where l is the molecular viscosity coefficient and lt is the eddy viscosity. The heat fluxes are calculated using 1 l lt oc2 þ ; ð5Þ qi ¼ c 1 Pr Prt oxi where Pr is the Prandtl number (Pr = 0.72), Prt is the turbulent Prandtl number (Prt = 0.9), c is the specific ratio (c = 1.4), and c is the speed of sound. The convective fluxes in Eq. (3) are discretized using Roe’s finite difference splitting [13] and the fifth-order compact nonlinear scheme [8].
I. Nakamori, T. Ikohagi / Computers & Fluids 37 (2008) 161–169
The viscous fluxes in Eq. (3) are discretized using the second-order accurate central differences. Approximate factorization implicit time integration is used to produce a CFL number that is greater than one and to simulate transient motion in flow-time. In order to maintain second-order accuracy in time, a second-order backward scheme for time integration is combined with the Newton-iterative method and alternative direction implicit factorization. 2.2. Turbulence model The proposed hybrid method blends the RANS model and the MILES. However, the SGS model is not employed in the MILES. Therefore, a comparably simple RANS model, such as the Baldwin–Lomax zero-equation model, should be chosen. The eddy viscosity, lt, used in the RANS region is evaluated using the BL model, which is given by ( lt;inner ¼ ql2 jxj; if d < d c ; ð6Þ lt ¼ lt;outer ¼ 0:0168qC cp F wake F Kleb ; if d > d c : Here, d is the normal distance from the wall, and dc is the location at which the inner value is equal to the outer value. jxj is the magnitude of the vorticity vector. In Prandtl’s study, the mixing length l is defined as l ¼ jd½1 expðd þ =26Þ;
ð7Þ +
where j is the von Karman constant and d is the dimensionless distance from the wall. The function Fwake determines the outer length scale that is associated with the magnitude of the vorticity vector, jxj, and the difference between the maximum and minimum absolute value of velocity udif. Fwake is given by C wk d max u2dif F wake ¼ min y max F max ; ; ð8Þ F max where udif ¼ F max
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2 þ v2 þ w2
max
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2 þ v2 þ w2
min
less than 14, then the value of the turbulent eddy viscosity is reset to lt = 0 along the g line. The present study applies the Degani and Schiff (DS) [7] modification to manage separated flow to some extent. Even though the DS modification was originally designed for crossflow separation, it can be used to improve the BL model’s ability to determine the outer length scale in a massive separated flow around an airfoil, as described later. The present paper defines a well-defined peak in F(d) as a peak at which F(d) drops to less than 90% of its local maximum value. In order to avoid an incorrect peak value for F(d), the first four points from the wall are skipped while searching from the wall outward. In addition, if the inner velocity product ug=2 Æ ug=g1 is less than 0, the code stops searching because an adverse flow will occur at g = g1. Although the DS modification improves the prediction accuracy near the separation, it is still possible to overestimate the eddy viscosity after the separation point. After the DS modification has been instituted, the eddy viscosity is decreased by the present hybrid model to some extent when the computational cell deals with large eddies. The proposed hybrid method detects the turbulent motion in and around the eddy and determines the switch point at which RANS switches to MILES. MILES is then applied outward from this transition point. The new model utilizes the same model constant used in the dynamic Smagorinsky model [9,10] as a threshold in defining the conjugation line between RANS and MILES. The dynamic model constant is given as ðCDÞ2 ¼
1 M ij Lij ; 2 M ij M ij
ð11Þ
where 1 oui ouj S ij ¼ þ ; 2 oxj oxi
jSj ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffi 2S ij S ij ;
M ij ¼ a2 jSjS ij jSjS ij
;
and
1 ¼ maxðF ðdÞÞ j
1 Lij ¼ ui uj ui uj dij ðuk uk uk uk Þ: 3
and F ðdÞ ¼ ljxj:
163
ð9Þ
Here, dmax is the value of d that produces a local maximum value for F(d). The function FKleb represents the influence of the intermittence of the edge of the turbulent boundary layer, which is given by the empirical formula " 6 #1 C Kleb d : ð10Þ F Kleb ¼ 1 þ 5:5 d max The BL model utilizes the following empirical coefficients: j = 0.4, A = 26, Cwk = 0.25, Ccp = 1.6, and CKleb = 0.3. The treatment by the original BL model of the transition between the laminar and turbulent boundary layer is now addressed. If the maximum dimensionless turbulent eddy viscosity (lt/l1)max along the g line normal to the wall is
ð12Þ
In these equations, D is the filter width, which is defined as 1
D ¼ ð1=J Þ3 , and a is set to 2. The over-bar, denotes the test filter width. For example, u is computed as D2 o2 u o2 u o2 u u ¼ u þ þ : ð13Þ 24 ox2 oy 2 oz2 It is known that the magnitude of the model constant C is on the order of 0.1 in the turbulent boundary layer if the computational cell is able to resolve the large-scale turbulence structure near the wall [9]. Otherwise, coefficient C is approximately 0. In the present study, the eddy viscosity value is defined as follows: RANS ; if d < d 1 ; lt ð14Þ lt ¼ 0; if d > d 1 ;
I. Nakamori, T. Ikohagi / Computers & Fluids 37 (2008) 161–169
Coefficient of dynamic model, C
164
C = C1
μt of RANS
Eddy viscosity, μ t
d = d1
μt of hybrid model
0 0
Distance from the wall, d
Fig. 1. Explanatory diagram of the hybridization.
The simulation was first performed on the flow around the NACA633-018 airfoil, which is classified as a trailingedge-stall airfoil. During the experiment [11], the flow stalled as the turbulent separation point near the trailing edge moved upstream. Fig. 2 shows a comparison of the surface pressure distributions between the hybrid model and the BLDS model at an angle of 12.58. The results for the hybrid model with a threshold C1 of 0.13 agrees well with the experimental results, but a threshold C1 of 0.10 produces a slight difference in the surface pressure distribution on the suction side. When the threshold C1 is set to 0, which converts the turbulence model into a MILES, the prediction results differ significantly from the experimental results. Therefore, the threshold C1 is set to 0.13 in the present study. The two-dimensional computations using the hybrid method and the BLDS model coincide well with the experimental data. However, when the angle of attack a = 16.53, as in Fig. 3, the two-dimensional computations fail to accurately predict the surface pressure distributions.
where lRANS is the eddy viscosity value calculated by the t BLDS model, and lSGS is set to 0 because MILES is emt ployed without an explicit SGS model. d1 is the distance at which coefficient C becomes greater than the threshold C1, while traveling away from the wall. An explanatory diagram of the hybridization is shown in Fig. 1, and the threshold C1 has been empirically set to 0.13.
The objective of this experiment was to validate the hybrid model’s ability to accurately predict the aerodynamics near a stall. This assessment was performed by comparing the computational and experimental results. The experiments were previously conducted by McCullough and Gault [11] using the NACA test wind tunnel, and the angle of attack, and lift and drag coefficients, which were adjusted for the effects of the wind tunnel, were detailed in the report. The present study employed the modified angle of attack in order to simulate the flows around the airfoils. Three different types of airfoils were used: the NACA633-018 airfoil (trailing-edge stall), the NACA631012 airfoil (leading-edge stall), and the NACA64A-006 airfoil (thin-airfoil stall). The tests were conducted at a free-stream Mach number of M1 = 0.17 and a Reynolds number of Re = 5.8 · 106. The C-type topology grid system was used. This grid system contains 113 points on the suction side, 73 points on the pressure side, 35 points along the cross-section (wake region), 85 points in the normal direction to the wall, and 40 points in the spanwise direction. The spanwise dimension Lz/c of the computational domain was set at Lz/c = 1.0. These dimensions are sufficient to eliminate contamination due to the use of the periodic boundary condition in the spanwise direction [3]. The time step Dt was fixed at 5T/1000, and the time-average values were computed over a range of 200 times the flow-time T, which is defined as T = c/U1 using the code length c and the free-stream velocity U1.
-8
Experiment Hybrid(3D, C =0.13) Hybrid(3D, C =0.10) Hybrid(3D, C =0.0) BLDS(3D)
-7 -6 Cp
-5 -4 -3 -2 -1 0 1 0
0.2
0.4
0.6
0.8
1.0
x/c -9 -8 Experiment Hybrid(2D) BLDS(2D)
-7 -6 -5 Cp
3. Results and discussion
-9
-4 -3 -2 -1 0 1 0
0.2
0.4
0.6
0.8
1.0
x/c Fig. 2. Comparison of surface pressure distributions on the NACA633018 airfoil near stall. Numerical results with the hybrid model and the B–L model with the D–S modification at an angle of attack of 12.58. (a) Three-dimensional model. (b) Two-dimensional model.
I. Nakamori, T. Ikohagi / Computers & Fluids 37 (2008) 161–169
For comparison, the simulation results obtained using the SA model and the DES are shown in Fig. 3. All of the two-dimensional simulations conducted using the present hybrid model, the BLDS model, and the SA model fail to predict the overall surface pressure distributions. In addition, the three-dimensional simulations using the BLDS model and the DES give a smaller plateau region due to the shorter separation region, as compared with the experimental data. As shown in Fig. 6(b), no distinguishable three-dimensionality appears in the separated region using the DES. Thus, the DES does not take the effect of the three-dimensional flow into account, which leads to a degree of error in the surface pressure distribution in the separated region. If a finer grid system is used in the DES, then the DES will work as an LES in the region nearer to the wall boundary, and the flow prediction will be accurate. The eddy viscosity of the RANS in the separation region inhibits the growth of the three-dimensional
165
detached eddies, and the eddy viscosity is decreased by the two-step process of the Degani–Schiff modification and the hybrid model. As a result, the present hybrid model gives an almost perfect surface pressure prediction, although it is performed on the same grid resolution. From the point of view of the mesh dependency, the present hybrid model has an advantage over the DES and the BLDS model. The difference in the pressure distribution between the hybrid model (3D) and the BLDS model (3D) is due to a difference in the starting point of the separation on the suction side. Fig. 4 shows a comparison of the surface friction distributions on the suction side of the NACA633-018 airfoil at angles of attack of a = 12.58 and a = 16.53. Both distributions indicate the effect of the transition to turbulence at x/c = 0.018, with the friction coefficients becoming local maxima. The hybrid model, however, predicts the starting point of separation farther upstream than that predicted by the BLDS model at an angle of attack of a = 16.53. These
-9 -8
20
Experiment Hybrid(3D) BLDS(3D) DES(3D)
-7 -6 Cp
-5
Hybrid(3D) BLDS(3D)
15
-4 10
cf Re
0.5
-3 -2 -1
5
0 1
0
0.0
0.2
0.4
0.6
0.8
1.0 0.0
x/c
0.2
0.4
0.6
0.8
x/c
-9
20
-8 Experiment Hybrid(2D) BLDS(2D) SA(2D)
-7 -6 -4
0.5
Cp
-5
Hybrid(3D) BLDS(3D)
15
10
cf Re
-3 -2
5
-1 0
0
1 0.0
0.2
0.4
0.6
0.8
1.0
x/c Fig. 3. Comparison of surface pressure distributions on the NACA633018 airfoil near stall. Numerical results with the hybrid model and the B–L model with the D–S modification at an angle of attack of 16.53. (a) Three-dimensional model. (b) Two-dimensional model.
0.0
0.2
0.4
0.6
0.8
x/c Fig. 4. Comparison of surface friction distributions on the suction side of the NACA633-018 airfoil after stall. Numerical results with the hybrid model and the B–L model with the D–S modification. (a) a = 12.58. (b) a = 16.53.
166
I. Nakamori, T. Ikohagi / Computers & Fluids 37 (2008) 161–169
locations are x/c = 0.61 and x/c = 0.66, respectively. These two values differ due to the hybrid model’s diminution of the eddy viscosity to near zero away from the wall in the separation region. Fig. 5 shows a comparison of the timeaverage eddy viscosity distributions on the suction side with respect to the distance from the wall at x/c = 0.2, x/ c = 0.4 and x/c = 0.5 at the angles of attack a = 12.58 and a = 16.53. The hybrid and BLDS models have identical eddy viscosity profiles at x/c = 0.2, where the flow should be handled in the RANS mode, because no separation appears near the leading edge in this case. In addition, Fig. 5(a) shows a gradual diminution of the eddy viscosity downstream. However, the hybridization decreases the eddy viscosity downstream, as evidenced in Fig. 5(b), in which the angle of attack is a = 16.53. This means that an important feature of the hybridization is that the present model automatically controls the switching point of LES/ RANS according to separation. Fig. 6 shows the visualized flow field around the NACA633-018 airfoil at an angle of attack a = 16.53.
The contours are colored according to the pressure distribution, and the isosurfaces of the spanwise velocity component w = ±0.02 are interposed. When the angle of attack is a = 16.53, the beginning of the separation is x/c = 0.61, as shown in Fig. 4(b), but the present model predicts a disturbed flow structure upstream of the separation point, and the resolved large eddy contributes to the hybridization’s diminution of the eddy viscosity. The two-dimensional computations in Figs. 2 and 3 are similar because no three-dimensional resolvable eddy is solved. When the proposed hybrid model is unable to resolve the detached eddy, the turbulence model automatically returns to the RANS mode, and the two-dimensional computations do not produce any three-dimensional eddies. Therefore, the two-dimensional computations of the hybrid method and the BLDS model produce nearly identical pressure distributions, as shown in Figs. 2 and 3. Next, the grid dependency examination in the spanwise direction is performed at an angle of attack of a = 16.53, and the resulting surface pressure distributions are shown
400 x/c = 0.2 Hybrid BLDS
200
<μt / μ∞>
<μt / μ∞>
400
0 0.02 0.04 d/c
0.06
0.00
400
0.02 0.04 d/c
0.06
400 x/c = 0.4 Hybrid BLDS
200
<μt / μ∞>
<μt / μ∞>
200
0 0.00
0
x/c = 0.4 Hybrid BLDS
200
0 0.00
0.02 0.04 d/c
0.06
0.00
400
0.02 0.04 d/c
0.06
400 x/c = 0.5 Hybrid BLDS
200
<μt / μ∞>
<μt / μ∞>
x/c = 0.2 Hybrid BLDS
0
x/c = 0.5 Hybrid BLDS
200
0 0.00
0.02 0.04 d/c
0.06
0.00
0.02 0.04 d/c
0.06
Fig. 5. Comparison of eddy viscosity distributions along with the distance from the wall on the suction side of the NACA633-018 airfoil. (a) a = 12.58. (b) a = 16.53.
I. Nakamori, T. Ikohagi / Computers & Fluids 37 (2008) 161–169
167
Fig. 6. Comparison of surface friction distributions on the suction side of the NACA633-018 airfoil at an angle of attack of a = 16.53. (a) Hybrid model. (b) DES.
-9 -8 Experiment Nz=40 20
-7 -6 -5 Cp
in Fig. 7. The computations with the hybrid model are performed in the spanwise domain of Lz/c = 1.0. The numbers of cells in the spanwise direction are Nz = 20 and 40, respectively. The results for these two cases are nearly identical. This indicates a converged solution with respect to the spanwise resolution. Finally, Fig. 8 shows the lift coefficient distributions predicted by the hybrid model, Cl, with respect to the angle of attack, a, using 256 · 84 · 40 grids for the three-dimensional simulation. The three-dimensional simulation generates maximum lift near a = 14.59, which is the experimental stall angle. In order to verify the stall prediction for the NACA631-012 (leading-edge stall) and the NACA64A-006 (thin-airfoil stall) airfoils, additional similar three-dimensional simulations were performed on the 256 · 84 · 40 grids. Fig. 9 shows the simulation curves and the experimental results for the lift coefficients. The stall predictions were nearly perfect for both cases. The mechanism of thin-airfoil stall is different from that of other types of stall. For the NACA64A-006 air-
-4 -3 -2 -1 0 1 0.0
0.2
0.4
0.6
0.8
1.0
x/c Fig. 7. Comparison of surface pressure distributions on the NACA633018 airfoil at an angle of attack 16.53 using the hybrid model on 20 and 40 grids in the spanwise direction.
168
I. Nakamori, T. Ikohagi / Computers & Fluids 37 (2008) 161–169
1.0
tion of the long bubbles over the suction side at the larger angle of attack. At that time, the lift coefficient curve gives the extremal value near the angle of a = 9 in Fig. 9. These results demonstrate that the proposed hybrid model is able to predict different types of stalls around the airfoils.
0.8
4. Summary
1.6 1.4
Cl
1.2
0.6
Lift Experiment Hybrid(3D)
0.4 0.2 0.0 0
2
4
6
8
α
10 12 14 16 18 20
Cl
Fig. 8. Comparison of computed and experimental lift coefficient for the NACA633-018 airfoil.
1.8 1.6 1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0
Experiment Hybrid
A new hybrid model composed of the RANS model with the MILES (fifth-order compact scheme) was proposed, and its feasibility was assessed via the turbulent flow around the airfoil and stall simulations. Unlike the existing hybrid method, the proposed method does not require the conjugation point at which the RANS mode is switched to the LES mode to be fixed. This feature is important when computing the flow for different types of stall, because the laminar/turbulent separation appears when a stall is present. The simulations demonstrate that the hybrid model can successfully predict the trailing-edge stall and the surface pressure distributions when compared to the RANS model. In addition, for leading-edge and thin-airfoil stalls, the proposed method generates an adequate angle of stall prediction. This means that the present hybrid model is able to predict different types of stalls around the airfoils. Acknowledgement
0
2
4
6
8
10 12 14 16 18 20 22
The computation was performed using the SX-5 supercomputer of the Institute of Fluid Science, Tohoku University.
α
References 1.0 0.8
Cl
0.6 0.4
Experiment Hybrid
0.2 0.0 0
2
4
6
8
10
12
α
Fig. 9. Comparison of computed and experimental lift coefficients for the NACA631-012 airfoil and the NACA64A-006 airfoil.
foil, the short bubbles appear near the leading edge at angle of a = 5.5. The short bubbles promote the forma-
[1] Baldwin BS, Lomax H. Thin layer approximation and algebraic model for separated turbulent flows. AIAA Paper, 78-0257, 1978. [2] Boris JP, Grinstein FF, Oran ES, Kolbe RL. New insights into large eddy simulation. Fluid Dyn Res 1992;10:199–228. [3] Chapman DR. Computational aerodynamics development and outlook. AIAA J 1979;17:1293–313. [4] Constantinescu G, Chapelet MC, Squires KD. Turbulence modeling applied to flow over a sphere. AIAA J 2003;41(9):1733–43. [5] Constantinescu G, Squires KD. Numerical investigations of flow over a sphere in the subcritical and supercritical regimes. Phys Fluids 2004;16(5):1449–66. [7] Degani D, Schiff LB. Computation of turbulent supersonic flows around pointed bodies having crossflow separation. J Comput Phys 1986;66:173–96. [8] Deng X, Zhang H. Developing high-order weighted compact nonlinear schemes. J Comput Phys 2000;165:22–44. [9] Germano M, Piomelli U, Moin P, Cabot WH. Dynamic subgrid-scale eddy viscosity model. Phys Fluids A 1991;3(7):1760–5. [10] Lilly DK. A proposed modification of the germano subgrid-scale closure method. Phys Fluids A 1992;4:633–5. [11] McCullough GB, Gault DE. Examples of three representative types of airfoil section stall at low speeds. NACA-TN-2502, 1951. [12] Nikitin NV, Nicoud F, Wasistho B, Squires KD, Spalart PR. An approach to wall modeling in large-eddy simulations. Phys Fluids 2001;12:1629–32.
I. Nakamori, T. Ikohagi / Computers & Fluids 37 (2008) 161–169 [13] Roe PL. Approximate Riemann solvers, parameter vectors and difference schemes. J Comput Phys 1981;43:357–72. [14] Shur ML, Spalart PR, Strelets MK, Travin AK. Detached-eddy simulation of an airfoil at high angle of attack. In: Rodi W, Laurence D, editors. Proceedings of the fourth international symposium on
169
engineering turbulencce modelling and measurements. Amsterdam: Elsevier; 1999. p. 669–78. [15] Spalart PR, Allmaras SR. A one-equation turbulence model for aerodynamic flows. La Recherche Aerospatiale, 1, 1994. p. 5–21; also AIAA Paper 92-0439, 1992.
Computers & Fluids 37 (2008) 161–169 www.elsevier.com/locate/compfluid
Dynamic hybridization of MILES and RANS for predicting airfoil stall Ichiro Nakamori *, Toshiaki Ikohagi
1
Institute of Fluid Science, Tohoku University, 2-1-1, Katahira, Aoba-ku, Sendai, Miyagi 980-8577, Japan Received 29 September 2005; received in revised form 11 September 2006; accepted 30 June 2007 Available online 30 August 2007
Abstract Hybridization comprised of an algebraic turbulence model based on the Reynolds average Navier–Stokes (RANS) equations with a monotonically integrated large eddy simulation (MILES) is proposed to simulate transient fluid motion during separation and vortex shedding at high Reynolds numbers. The proposed hybridization utilizes the Baldwin–Lomax model with the Degani–Schiff modification as the RANS model in the near-wall region and a MILES far from the wall. Although the hybridizationis assumed to be a MILES with wall modeling, the transition line between the RANS and the MILES modes is determined by the turbulent intensity that is dominated by the large eddies in the grid-scale. This hybrid model is applied to the flows past three different types of airfoils (NACA633-018, NACA631-012 and NACA64A-006) near stall, at a chord Reynolds number of Re = 5.8 · 106. These airfoils are classified as trailingedge-stall, leading-edge-stall and thin-airfoil-stall airfoils, respectively. The computed results are compared with wind tunnel experiments. The hybrid model successfully demonstrates accurate stall angle and surface pressure distribution predictions near the stall for each type of airfoil. The airfoil simulation results confirmed that the hybrid model provides a better prediction than the RANS model for unsteady turbulent flows with separation and vortex shedding simulations. 2007 Elsevier Ltd. All rights reserved.
1. Introduction Recently, CFD has become the preferred method for predicting aerodynamic performance during aircraft design. In particular, aerodynamic analysis at stall with massive separation is especially important when evaluating an aircraft’s ability to take off and land. However, current Reynolds average Navier–Stokes (RANS) models cannot accurately predict turbulent flow with vortex shedding or massive separation, because it is difficult to develop a RANS model that is suitable for a turbulent flow with a high degree of unsteadiness. Large eddy simulations (LES), however, could prove useful in predicting unsteady aerodynamic characteristics, because they are independent of the Reynolds average, thereby allowing for the resolution of highly unsteady fluid motion. Unfortunately, *
1
Corresponding author. Tel.: +81 22 217 5228; fax: +81 22 217 5229. E-mail address: [email protected] (I. Nakamori). Tel.: +81 22 217 5228; fax: +81 22 217 5229.
0045-7930/$ - see front matter 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2007.06.004
supercomputers are not currently able to resolve all eddies, because eddies near walls decrease in size with an increase in the Reynolds number. This problem could be avoided by employing a wall function for a specific turbulent flow. However, the wall function generally cannot be employed when the boundary layer includes the laminar to turbulence transition point. Another feasible option for bypassing the enormous CPU and storage requirement and resolving the wall turbulence is hybridization. This involves conducting a RANS near the wall and an LES for the detached eddies, as conducted in detached eddy simulation (DES) [4,5,14]. In DES, the eddy viscosity value switches between the Spalart–Allmaras (SA) model [15] and the LES mode using one eddy viscosity equation. The designated point is determined by the grid-scale. Because the DES switches from RANS to LES at a designated point away from the wall, the glitch point in the velocity profile occurs near the switching point. An example of this characteristic is demonstrated in Nikitin et al. [12]. In addition, the subgrid-scale (SGS) model in LES
162
I. Nakamori, T. Ikohagi / Computers & Fluids 37 (2008) 161–169
is effective for prediction accuracy analysis only when sufficient grid resolution is used in all directions: the streamwise direction, the spanwise direction, and the normal direction to the wall. However, the SGS model is inappropriate for practical grid resolution, such that the aspect ratio of the streamwise grid size and the grid size in the normal direction to the wall (or the aspect ratio of the spanwise grid size and the grid size in the normal direction to the wall) is extremely high near the wall. In the present study, in order to reduce the glitch point in the velocity profile, the switching point at which the turbulent intensity of the resolved detached eddy reaches a certain fixed level is noted. Further, a monotonically integrated large eddy simulation (MILES) [2] is applied as an LES after the switching point far from the wall. This MILES uses inherent numerical dissipation to damp out the SGS turbulence. The MILES leads to the DNS when a finer grid system is used. To save CPU time, the Baldwin–Lomax [1] algebraic model with the modification by Degani and Schiff [7] (BLDS model) is employed as a RANS model in conjunction with the MILES. The concept of the Degani–Schiff modification coincides with the directivity of the present hybrid model from the point of view that the excessive eddy viscosity of the RANS mode is decreased in the separated vortical flow. The objective of this study is to obtain accurate prediction of stall with massive separation by using less computational time in comparison with large eddy simulations.The hybrid model is then applied to flows past various types of airfoils near stall, in order to determine its applicability for predicting aerodynamic mean values and unsteady characteristics. The experimental aerodynamic characteristics of airfoils have previously been described by McCullough and Gault [11]. Each airfoil is categorized by its stalling characteristics as a trailing-edge-stall airfoil, a leading-edge-stall airfoil or a thin-airfoil-stall airfoil. Because it is challenging to predict the necessary angle of attack to produce a stall by CFD, the improvement of prediction accuracy of the stall angle of attack is very significant.
2. Numerical method 2.1. Governing equations and numerical scheme The Navier–Stokes equations are solved in an inertial system. The primitive variables in the solution vector for the RANS mode use the Reynolds average and those in the solution vector for the LES mode use the spatial average. The curvilinear coordinates are used to transform the equations: oQ oE oF oG 1 oE v oF v oG v þ þ þ ¼ þ þ : ot on og of Re on og of
ð1Þ
The transformed Q and flux vectors are given by 1 Q ¼ ðq; qu; qv; qw; eÞ J ^ b þ ny F b þ nz G nx E E¼ ; J ^ b þ gy F b þ gz G gx E ; F¼ J ^ b þ fy F b þ fz G fx E G¼ ; J where 2
qu
3
^ b þ ny F b þ nz G nx E ; J ^ b þ gy F b þ gz G gx E Fv ¼ ; J ^ b þ fy F b þ fz G fx E Gv ¼ ; J Ev ¼
2
qv
3
7 7 6 2 6 6 qu þ p 7 6 quv 7 7 7 6 6 2 7 7 6 E¼6 6 quv 7; F ¼ 6 qv þ p 7; 7 7 6 6 4 quw 5 4 qvw 5 ðe þ pÞu ðe þ pÞv 3 2 0 7 6 sxx 7 6 7 6 7; 6 sxy Ev ¼ 6 7 7 6 sxz 5 4 usxx þ vsxy þ wsxz 3 2 0 7 6 sxy 7 6 7 6 7; s Fv ¼ 6 yy 7 6 7 6 syz 5 4
2
ð2Þ
qw
3
7 6 6 quw 7 7 6 7 G ¼6 6 qvw 7; 7 6 4 qw2 þ p 5 ðe þ pÞw
usxy þ vsyy þ wsyz 3 0 7 6 sxz 7 6 7 6 7: 6 syz Gv ¼ 6 7 7 6 szz 5 4 usxz þ vsyz þ wszz 2
ð3Þ In these equations, J is the Jacobian of the coordinate transformation, q is the density, u, v, and w are the fluid velocity components, e is the total energy per unit volume, p is the pressure, and qi are the heat fluxes. The viscous stress tensors sij can be expressed as oui ouj 2 ouk þ dij sij ¼ ðl þ lt Þ ; ð4Þ oxj oxi 3 oxk where l is the molecular viscosity coefficient and lt is the eddy viscosity. The heat fluxes are calculated using 1 l lt oc2 þ ; ð5Þ qi ¼ c 1 Pr Prt oxi where Pr is the Prandtl number (Pr = 0.72), Prt is the turbulent Prandtl number (Prt = 0.9), c is the specific ratio (c = 1.4), and c is the speed of sound. The convective fluxes in Eq. (3) are discretized using Roe’s finite difference splitting [13] and the fifth-order compact nonlinear scheme [8].
I. Nakamori, T. Ikohagi / Computers & Fluids 37 (2008) 161–169
The viscous fluxes in Eq. (3) are discretized using the second-order accurate central differences. Approximate factorization implicit time integration is used to produce a CFL number that is greater than one and to simulate transient motion in flow-time. In order to maintain second-order accuracy in time, a second-order backward scheme for time integration is combined with the Newton-iterative method and alternative direction implicit factorization. 2.2. Turbulence model The proposed hybrid method blends the RANS model and the MILES. However, the SGS model is not employed in the MILES. Therefore, a comparably simple RANS model, such as the Baldwin–Lomax zero-equation model, should be chosen. The eddy viscosity, lt, used in the RANS region is evaluated using the BL model, which is given by ( lt;inner ¼ ql2 jxj; if d < d c ; ð6Þ lt ¼ lt;outer ¼ 0:0168qC cp F wake F Kleb ; if d > d c : Here, d is the normal distance from the wall, and dc is the location at which the inner value is equal to the outer value. jxj is the magnitude of the vorticity vector. In Prandtl’s study, the mixing length l is defined as l ¼ jd½1 expðd þ =26Þ;
ð7Þ +
where j is the von Karman constant and d is the dimensionless distance from the wall. The function Fwake determines the outer length scale that is associated with the magnitude of the vorticity vector, jxj, and the difference between the maximum and minimum absolute value of velocity udif. Fwake is given by C wk d max u2dif F wake ¼ min y max F max ; ; ð8Þ F max where udif ¼ F max
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2 þ v2 þ w2
max
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2 þ v2 þ w2
min
less than 14, then the value of the turbulent eddy viscosity is reset to lt = 0 along the g line. The present study applies the Degani and Schiff (DS) [7] modification to manage separated flow to some extent. Even though the DS modification was originally designed for crossflow separation, it can be used to improve the BL model’s ability to determine the outer length scale in a massive separated flow around an airfoil, as described later. The present paper defines a well-defined peak in F(d) as a peak at which F(d) drops to less than 90% of its local maximum value. In order to avoid an incorrect peak value for F(d), the first four points from the wall are skipped while searching from the wall outward. In addition, if the inner velocity product ug=2 Æ ug=g1 is less than 0, the code stops searching because an adverse flow will occur at g = g1. Although the DS modification improves the prediction accuracy near the separation, it is still possible to overestimate the eddy viscosity after the separation point. After the DS modification has been instituted, the eddy viscosity is decreased by the present hybrid model to some extent when the computational cell deals with large eddies. The proposed hybrid method detects the turbulent motion in and around the eddy and determines the switch point at which RANS switches to MILES. MILES is then applied outward from this transition point. The new model utilizes the same model constant used in the dynamic Smagorinsky model [9,10] as a threshold in defining the conjugation line between RANS and MILES. The dynamic model constant is given as ðCDÞ2 ¼
1 M ij Lij ; 2 M ij M ij
ð11Þ
where 1 oui ouj S ij ¼ þ ; 2 oxj oxi
jSj ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffi 2S ij S ij ;
M ij ¼ a2 jSjS ij jSjS ij
;
and
1 ¼ maxðF ðdÞÞ j
1 Lij ¼ ui uj ui uj dij ðuk uk uk uk Þ: 3
and F ðdÞ ¼ ljxj:
163
ð9Þ
Here, dmax is the value of d that produces a local maximum value for F(d). The function FKleb represents the influence of the intermittence of the edge of the turbulent boundary layer, which is given by the empirical formula " 6 #1 C Kleb d : ð10Þ F Kleb ¼ 1 þ 5:5 d max The BL model utilizes the following empirical coefficients: j = 0.4, A = 26, Cwk = 0.25, Ccp = 1.6, and CKleb = 0.3. The treatment by the original BL model of the transition between the laminar and turbulent boundary layer is now addressed. If the maximum dimensionless turbulent eddy viscosity (lt/l1)max along the g line normal to the wall is
ð12Þ
In these equations, D is the filter width, which is defined as 1
D ¼ ð1=J Þ3 , and a is set to 2. The over-bar, denotes the test filter width. For example, u is computed as D2 o2 u o2 u o2 u u ¼ u þ þ : ð13Þ 24 ox2 oy 2 oz2 It is known that the magnitude of the model constant C is on the order of 0.1 in the turbulent boundary layer if the computational cell is able to resolve the large-scale turbulence structure near the wall [9]. Otherwise, coefficient C is approximately 0. In the present study, the eddy viscosity value is defined as follows: RANS ; if d < d 1 ; lt ð14Þ lt ¼ 0; if d > d 1 ;
I. Nakamori, T. Ikohagi / Computers & Fluids 37 (2008) 161–169
Coefficient of dynamic model, C
164
C = C1
μt of RANS
Eddy viscosity, μ t
d = d1
μt of hybrid model
0 0
Distance from the wall, d
Fig. 1. Explanatory diagram of the hybridization.
The simulation was first performed on the flow around the NACA633-018 airfoil, which is classified as a trailingedge-stall airfoil. During the experiment [11], the flow stalled as the turbulent separation point near the trailing edge moved upstream. Fig. 2 shows a comparison of the surface pressure distributions between the hybrid model and the BLDS model at an angle of 12.58. The results for the hybrid model with a threshold C1 of 0.13 agrees well with the experimental results, but a threshold C1 of 0.10 produces a slight difference in the surface pressure distribution on the suction side. When the threshold C1 is set to 0, which converts the turbulence model into a MILES, the prediction results differ significantly from the experimental results. Therefore, the threshold C1 is set to 0.13 in the present study. The two-dimensional computations using the hybrid method and the BLDS model coincide well with the experimental data. However, when the angle of attack a = 16.53, as in Fig. 3, the two-dimensional computations fail to accurately predict the surface pressure distributions.
where lRANS is the eddy viscosity value calculated by the t BLDS model, and lSGS is set to 0 because MILES is emt ployed without an explicit SGS model. d1 is the distance at which coefficient C becomes greater than the threshold C1, while traveling away from the wall. An explanatory diagram of the hybridization is shown in Fig. 1, and the threshold C1 has been empirically set to 0.13.
The objective of this experiment was to validate the hybrid model’s ability to accurately predict the aerodynamics near a stall. This assessment was performed by comparing the computational and experimental results. The experiments were previously conducted by McCullough and Gault [11] using the NACA test wind tunnel, and the angle of attack, and lift and drag coefficients, which were adjusted for the effects of the wind tunnel, were detailed in the report. The present study employed the modified angle of attack in order to simulate the flows around the airfoils. Three different types of airfoils were used: the NACA633-018 airfoil (trailing-edge stall), the NACA631012 airfoil (leading-edge stall), and the NACA64A-006 airfoil (thin-airfoil stall). The tests were conducted at a free-stream Mach number of M1 = 0.17 and a Reynolds number of Re = 5.8 · 106. The C-type topology grid system was used. This grid system contains 113 points on the suction side, 73 points on the pressure side, 35 points along the cross-section (wake region), 85 points in the normal direction to the wall, and 40 points in the spanwise direction. The spanwise dimension Lz/c of the computational domain was set at Lz/c = 1.0. These dimensions are sufficient to eliminate contamination due to the use of the periodic boundary condition in the spanwise direction [3]. The time step Dt was fixed at 5T/1000, and the time-average values were computed over a range of 200 times the flow-time T, which is defined as T = c/U1 using the code length c and the free-stream velocity U1.
-8
Experiment Hybrid(3D, C =0.13) Hybrid(3D, C =0.10) Hybrid(3D, C =0.0) BLDS(3D)
-7 -6 Cp
-5 -4 -3 -2 -1 0 1 0
0.2
0.4
0.6
0.8
1.0
x/c -9 -8 Experiment Hybrid(2D) BLDS(2D)
-7 -6 -5 Cp
3. Results and discussion
-9
-4 -3 -2 -1 0 1 0
0.2
0.4
0.6
0.8
1.0
x/c Fig. 2. Comparison of surface pressure distributions on the NACA633018 airfoil near stall. Numerical results with the hybrid model and the B–L model with the D–S modification at an angle of attack of 12.58. (a) Three-dimensional model. (b) Two-dimensional model.
I. Nakamori, T. Ikohagi / Computers & Fluids 37 (2008) 161–169
For comparison, the simulation results obtained using the SA model and the DES are shown in Fig. 3. All of the two-dimensional simulations conducted using the present hybrid model, the BLDS model, and the SA model fail to predict the overall surface pressure distributions. In addition, the three-dimensional simulations using the BLDS model and the DES give a smaller plateau region due to the shorter separation region, as compared with the experimental data. As shown in Fig. 6(b), no distinguishable three-dimensionality appears in the separated region using the DES. Thus, the DES does not take the effect of the three-dimensional flow into account, which leads to a degree of error in the surface pressure distribution in the separated region. If a finer grid system is used in the DES, then the DES will work as an LES in the region nearer to the wall boundary, and the flow prediction will be accurate. The eddy viscosity of the RANS in the separation region inhibits the growth of the three-dimensional
165
detached eddies, and the eddy viscosity is decreased by the two-step process of the Degani–Schiff modification and the hybrid model. As a result, the present hybrid model gives an almost perfect surface pressure prediction, although it is performed on the same grid resolution. From the point of view of the mesh dependency, the present hybrid model has an advantage over the DES and the BLDS model. The difference in the pressure distribution between the hybrid model (3D) and the BLDS model (3D) is due to a difference in the starting point of the separation on the suction side. Fig. 4 shows a comparison of the surface friction distributions on the suction side of the NACA633-018 airfoil at angles of attack of a = 12.58 and a = 16.53. Both distributions indicate the effect of the transition to turbulence at x/c = 0.018, with the friction coefficients becoming local maxima. The hybrid model, however, predicts the starting point of separation farther upstream than that predicted by the BLDS model at an angle of attack of a = 16.53. These
-9 -8
20
Experiment Hybrid(3D) BLDS(3D) DES(3D)
-7 -6 Cp
-5
Hybrid(3D) BLDS(3D)
15
-4 10
cf Re
0.5
-3 -2 -1
5
0 1
0
0.0
0.2
0.4
0.6
0.8
1.0 0.0
x/c
0.2
0.4
0.6
0.8
x/c
-9
20
-8 Experiment Hybrid(2D) BLDS(2D) SA(2D)
-7 -6 -4
0.5
Cp
-5
Hybrid(3D) BLDS(3D)
15
10
cf Re
-3 -2
5
-1 0
0
1 0.0
0.2
0.4
0.6
0.8
1.0
x/c Fig. 3. Comparison of surface pressure distributions on the NACA633018 airfoil near stall. Numerical results with the hybrid model and the B–L model with the D–S modification at an angle of attack of 16.53. (a) Three-dimensional model. (b) Two-dimensional model.
0.0
0.2
0.4
0.6
0.8
x/c Fig. 4. Comparison of surface friction distributions on the suction side of the NACA633-018 airfoil after stall. Numerical results with the hybrid model and the B–L model with the D–S modification. (a) a = 12.58. (b) a = 16.53.
166
I. Nakamori, T. Ikohagi / Computers & Fluids 37 (2008) 161–169
locations are x/c = 0.61 and x/c = 0.66, respectively. These two values differ due to the hybrid model’s diminution of the eddy viscosity to near zero away from the wall in the separation region. Fig. 5 shows a comparison of the timeaverage eddy viscosity distributions on the suction side with respect to the distance from the wall at x/c = 0.2, x/ c = 0.4 and x/c = 0.5 at the angles of attack a = 12.58 and a = 16.53. The hybrid and BLDS models have identical eddy viscosity profiles at x/c = 0.2, where the flow should be handled in the RANS mode, because no separation appears near the leading edge in this case. In addition, Fig. 5(a) shows a gradual diminution of the eddy viscosity downstream. However, the hybridization decreases the eddy viscosity downstream, as evidenced in Fig. 5(b), in which the angle of attack is a = 16.53. This means that an important feature of the hybridization is that the present model automatically controls the switching point of LES/ RANS according to separation. Fig. 6 shows the visualized flow field around the NACA633-018 airfoil at an angle of attack a = 16.53.
The contours are colored according to the pressure distribution, and the isosurfaces of the spanwise velocity component w = ±0.02 are interposed. When the angle of attack is a = 16.53, the beginning of the separation is x/c = 0.61, as shown in Fig. 4(b), but the present model predicts a disturbed flow structure upstream of the separation point, and the resolved large eddy contributes to the hybridization’s diminution of the eddy viscosity. The two-dimensional computations in Figs. 2 and 3 are similar because no three-dimensional resolvable eddy is solved. When the proposed hybrid model is unable to resolve the detached eddy, the turbulence model automatically returns to the RANS mode, and the two-dimensional computations do not produce any three-dimensional eddies. Therefore, the two-dimensional computations of the hybrid method and the BLDS model produce nearly identical pressure distributions, as shown in Figs. 2 and 3. Next, the grid dependency examination in the spanwise direction is performed at an angle of attack of a = 16.53, and the resulting surface pressure distributions are shown
400 x/c = 0.2 Hybrid BLDS
200
<μt / μ∞>
<μt / μ∞>
400
0 0.02 0.04 d/c
0.06
0.00
400
0.02 0.04 d/c
0.06
400 x/c = 0.4 Hybrid BLDS
200
<μt / μ∞>
<μt / μ∞>
200
0 0.00
0
x/c = 0.4 Hybrid BLDS
200
0 0.00
0.02 0.04 d/c
0.06
0.00
400
0.02 0.04 d/c
0.06
400 x/c = 0.5 Hybrid BLDS
200
<μt / μ∞>
<μt / μ∞>
x/c = 0.2 Hybrid BLDS
0
x/c = 0.5 Hybrid BLDS
200
0 0.00
0.02 0.04 d/c
0.06
0.00
0.02 0.04 d/c
0.06
Fig. 5. Comparison of eddy viscosity distributions along with the distance from the wall on the suction side of the NACA633-018 airfoil. (a) a = 12.58. (b) a = 16.53.
I. Nakamori, T. Ikohagi / Computers & Fluids 37 (2008) 161–169
167
Fig. 6. Comparison of surface friction distributions on the suction side of the NACA633-018 airfoil at an angle of attack of a = 16.53. (a) Hybrid model. (b) DES.
-9 -8 Experiment Nz=40 20
-7 -6 -5 Cp
in Fig. 7. The computations with the hybrid model are performed in the spanwise domain of Lz/c = 1.0. The numbers of cells in the spanwise direction are Nz = 20 and 40, respectively. The results for these two cases are nearly identical. This indicates a converged solution with respect to the spanwise resolution. Finally, Fig. 8 shows the lift coefficient distributions predicted by the hybrid model, Cl, with respect to the angle of attack, a, using 256 · 84 · 40 grids for the three-dimensional simulation. The three-dimensional simulation generates maximum lift near a = 14.59, which is the experimental stall angle. In order to verify the stall prediction for the NACA631-012 (leading-edge stall) and the NACA64A-006 (thin-airfoil stall) airfoils, additional similar three-dimensional simulations were performed on the 256 · 84 · 40 grids. Fig. 9 shows the simulation curves and the experimental results for the lift coefficients. The stall predictions were nearly perfect for both cases. The mechanism of thin-airfoil stall is different from that of other types of stall. For the NACA64A-006 air-
-4 -3 -2 -1 0 1 0.0
0.2
0.4
0.6
0.8
1.0
x/c Fig. 7. Comparison of surface pressure distributions on the NACA633018 airfoil at an angle of attack 16.53 using the hybrid model on 20 and 40 grids in the spanwise direction.
168
I. Nakamori, T. Ikohagi / Computers & Fluids 37 (2008) 161–169
1.0
tion of the long bubbles over the suction side at the larger angle of attack. At that time, the lift coefficient curve gives the extremal value near the angle of a = 9 in Fig. 9. These results demonstrate that the proposed hybrid model is able to predict different types of stalls around the airfoils.
0.8
4. Summary
1.6 1.4
Cl
1.2
0.6
Lift Experiment Hybrid(3D)
0.4 0.2 0.0 0
2
4
6
8
α
10 12 14 16 18 20
Cl
Fig. 8. Comparison of computed and experimental lift coefficient for the NACA633-018 airfoil.
1.8 1.6 1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0
Experiment Hybrid
A new hybrid model composed of the RANS model with the MILES (fifth-order compact scheme) was proposed, and its feasibility was assessed via the turbulent flow around the airfoil and stall simulations. Unlike the existing hybrid method, the proposed method does not require the conjugation point at which the RANS mode is switched to the LES mode to be fixed. This feature is important when computing the flow for different types of stall, because the laminar/turbulent separation appears when a stall is present. The simulations demonstrate that the hybrid model can successfully predict the trailing-edge stall and the surface pressure distributions when compared to the RANS model. In addition, for leading-edge and thin-airfoil stalls, the proposed method generates an adequate angle of stall prediction. This means that the present hybrid model is able to predict different types of stalls around the airfoils. Acknowledgement
0
2
4
6
8
10 12 14 16 18 20 22
The computation was performed using the SX-5 supercomputer of the Institute of Fluid Science, Tohoku University.
α
References 1.0 0.8
Cl
0.6 0.4
Experiment Hybrid
0.2 0.0 0
2
4
6
8
10
12
α
Fig. 9. Comparison of computed and experimental lift coefficients for the NACA631-012 airfoil and the NACA64A-006 airfoil.
foil, the short bubbles appear near the leading edge at angle of a = 5.5. The short bubbles promote the forma-
[1] Baldwin BS, Lomax H. Thin layer approximation and algebraic model for separated turbulent flows. AIAA Paper, 78-0257, 1978. [2] Boris JP, Grinstein FF, Oran ES, Kolbe RL. New insights into large eddy simulation. Fluid Dyn Res 1992;10:199–228. [3] Chapman DR. Computational aerodynamics development and outlook. AIAA J 1979;17:1293–313. [4] Constantinescu G, Chapelet MC, Squires KD. Turbulence modeling applied to flow over a sphere. AIAA J 2003;41(9):1733–43. [5] Constantinescu G, Squires KD. Numerical investigations of flow over a sphere in the subcritical and supercritical regimes. Phys Fluids 2004;16(5):1449–66. [7] Degani D, Schiff LB. Computation of turbulent supersonic flows around pointed bodies having crossflow separation. J Comput Phys 1986;66:173–96. [8] Deng X, Zhang H. Developing high-order weighted compact nonlinear schemes. J Comput Phys 2000;165:22–44. [9] Germano M, Piomelli U, Moin P, Cabot WH. Dynamic subgrid-scale eddy viscosity model. Phys Fluids A 1991;3(7):1760–5. [10] Lilly DK. A proposed modification of the germano subgrid-scale closure method. Phys Fluids A 1992;4:633–5. [11] McCullough GB, Gault DE. Examples of three representative types of airfoil section stall at low speeds. NACA-TN-2502, 1951. [12] Nikitin NV, Nicoud F, Wasistho B, Squires KD, Spalart PR. An approach to wall modeling in large-eddy simulations. Phys Fluids 2001;12:1629–32.
I. Nakamori, T. Ikohagi / Computers & Fluids 37 (2008) 161–169 [13] Roe PL. Approximate Riemann solvers, parameter vectors and difference schemes. J Comput Phys 1981;43:357–72. [14] Shur ML, Spalart PR, Strelets MK, Travin AK. Detached-eddy simulation of an airfoil at high angle of attack. In: Rodi W, Laurence D, editors. Proceedings of the fourth international symposium on
169
engineering turbulencce modelling and measurements. Amsterdam: Elsevier; 1999. p. 669–78. [15] Spalart PR, Allmaras SR. A one-equation turbulence model for aerodynamic flows. La Recherche Aerospatiale, 1, 1994. p. 5–21; also AIAA Paper 92-0439, 1992.