- Email: [email protected]
Transonic wing stall of a blended flying wing common research model based on DDES method
CJA 729 12 November 2016 Chinese Journal of Aeronautics, (2016), xxx(xx): xxx–xxx
No. of Pages 12
1
Chinese Society of Aeronautics and Astronautics & Beihang University
Chinese Journal of Aeronautics [email protected] www.sciencedirect.com
4
Transonic wing stall of a blended flying wing common research model based on DDES method
5
Tao Yang a,b,*, Li Yonghong b, Zhang Zhao b, Zhao Zhongliang b, Liu Zhiyong b
3
6 7
8
a b
State Key Laboratory of Aerodynamics, China Aerodynamics Research and Development Center, Mianyang 621000, China High Speed Aerodynamic Institute, China Aerodynamics Research and Development Center, Mianyang 621000, China
Received 29 September 2015; revised 24 December 2015; accepted 31 January 2016
9
11 12
KEYWORDS
13
Delayed detached eddy simulation; Flying wing; Vortex lift; Vortex breakdown; Wing stall
14 15 16 17 18
19
Abstract Numerical simulation of wing stall of a blended flying wing configuration at transonic speed was conducted using both delayed detached eddy simulation (DDES) and unsteady Reynolds-averaged Navier-Stokes (URANS) equations methods based on the shear stress transport (SST) turbulence model for a free-stream Mach number 0.9 and a Reynolds number 9.6 106. A joint time step/grid density study is performed based on power spectrum density (PSD) analysis of the frequency content of forces or moments, and medium mesh and the normalized time scale 0.010 were suggested for this simulation. The simulation results show that the DDES methods perform more precisely than the URANS method and the aerodynamic coefficient results from DDES method compare very well with the experiment data. The angle of attack of nonlinear vortex lift and abrupt wing stall of DDES results compare well with the experimental data. The flow structure of the DDES computation shows that the wing stall is caused mainly by the leeward vortex breakdown which occurred at x/xcr = 0.6 at angle of attack of 14°. The DDES methods show advantage in the simulation problem with separation flow. The computed result shows that a shock/vortex interaction is responsible for the wing stall caused by the vortex breakdown. The balance of the vortex strength and axial flow, and the shock strength, is examined to provide an explanation of the sensitivity of the breakdown location. Wing body thickness has a great influence on shock and shock/ vortex interactions, which can make a significant difference to the vortex breakdown behavior and stall characteristic of the blended flying wing configuration. Ó 2016 Production and hosting by Elsevier Ltd. on behalf of Chinese Society of Aeronautics and Astronautics. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/ licenses/by-nc-nd/4.0/).
* Corresponding author at: State Key Laboratory of Aerodynamics, China Aerodynamics Research and Development Center, Mianyang 621000, China. E-mail address: [email protected] (Y. Tao). Peer review under responsibility of Editorial Committee of CJA.
Production and hosting by Elsevier
1. Introduction
20
The development of studying and understanding the phenomenon of vortex lift has greatly impacted the aerodynamic configuration of the modern fighters, and the aerodynamic configuration using delta wing with large sweep angle has become one of the research topics.1–3 This kind of aerodynamic layout has the evident advantages of maneu-
21
http://dx.doi.org/10.1016/j.cja.2016.10.018 1000-9361 Ó 2016 Production and hosting by Elsevier Ltd. on behalf of Chinese Society of Aeronautics and Astronautics. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
22 23 24 25 26
CJA 729 12 November 2016
No. of Pages 12
2 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85
verability, agility and stall incidence increase.4,5 The occurrence of shocks on delta wings introduces complex shock/ vortex interactions, particularly at moderate to high angles of incidence. These interactions can make a significant difference to the vortex breakdown behavior. With the increase in angle of attack, flight vehicles would suffer abrupt wing stall due to the condition change such as adverse pressure gradient. The reason for abrupt wing stall is nonlinear lift loss caused by the leeward vortex breakdown, and the mechanism and impact factors are very complicated.5–7 The strengthening of the shock which stands off the sting as the incidence is increased can lead to a shock/vortex interaction triggering breakdown. The location of breakdown can shift upstream by as much as 30% of the chord at a single 1° incidence interval2,3 due to this interaction. Much literature has reported studies of this problem via theoretical analysis, experiments and numerical computations.8 Because the wing stall and vortex breakdown phenomenon involve the computation of flow separation, the accuracy of unsteady Reynolds-averaged Navier-Stokes (URANS) equations method is insufficient and thus causes the deviation between the calculated stall angle and experimental data.9,10 As a hybrid method applied to simulating this kind of separation flow, delayed detached eddy simulation (DDES) takes the advantages of higher precision and acceptable computation consumption. It has been applied to many problems, such as shock vibration of transonic wing, cavity flow, and wing separation flow.11–15 Aerodynamic efficiency defined as MaCL/CD has been steadily increasing from 13 in 1960 to 16 in the mid-1990s16, where Ma is stream velocity Mach number, CL is lift coefficient and CD is drag coefficient. Higher aerodynamic efficiency will benefit the flight range and fuel economy, and thus aircraft designers pursue higher cruising speed as well as lower drag and higher lift to drag ratios. As an example the cruising speed increases from Ma = 0.78 for A320 to about Ma = 0.85 for A380. Blended-wing-body (BWB) configuration takes an advantage of lift to drag ratios compared with traditional configurations. The developing BWB aircraft of Boeing Company aims at a design cruising Ma = 0.85. The low aspect ratio blended flying wing configuration has more potential on increasing cruising speed. The most typical example is X47b which is defined as a cruising speed at high subsonic speed and some literature explicitly refers to it as about Ma = 0.9. The measured data from 2.4 m 2.4 m transonic wind tunnel at transonic conditions 17 (Ma = 0.8 and 0.9) showed a sudden jump of the breakdown location toward the wing apex of a blended flying wing common research model when a critical angle of incidence was reached. All the computation study was conducted at Ma = 0.9 because the flow condition at this Mach number was more similar to the design cruising speed and the shock wave/ vortex interaction phenomenon was more typical. DDES model was employed to study the transonic stall problem, and by comparing with URANS and experimental data, the accuracy and applicability were validated. The flow structures and the role of vortex breakdown at stall angles of attack were briefly analyzed, too.
Y. Tao et al. 2. Mathematical formulation and numerical method
86
2.1. Governing equations and numerical procedure
87
The governing equations, which are the three-dimensional time-dependent compressible Navier-Stokes equations, can be written in terms of generalized coordinates as
88
^ @ E^ @ F^ @ G^ @ E^V @ F^V @ G^V @Q þ þ þ þ þ ¼ @n @g @f @t @n @g @f
ð1Þ
b is the vector of conserved variables, density, momenwhere Q b are inviscid flux b F, b G tum and total energy per unit volume, E, b b b terms and, E V , F V , G V are viscous flux terms. A general, 3D transformation between the Cartesian variables (x, y, z) and the generalized coordinated (n, g, f) is implied. The equation of state for an ideal gas is used and the molecular viscosity is assumed to obey Sutherland’s law. To nondimensionalize the equations, we use the free-stream variables including the density, temperature, speed of sound, and root chord of the blended wing as characteristic scales. The initial and boundary conditions are presented as follows. The initial condition is set as the free-stream quantities. The far field boundary conditions are treated by local onedimensional Riemann-invariants. No-slip and adiabatic conditions are applied on the aircraft surface. The governing Eq. (1) is discretized using finite volume method and the time term is discretized by an implicit approximate-factorization method with sub-iterations to ensure the second-order accuracy Simon et al.18. In this paper, the convective terms were adopted by second-order central difference scheme and a fourth-order artificial viscous term was introduced to restrain the numerical oscillation Lu et al.19 and Wang et al.20. To capture the discontinuity caused by shockwave, a second-order upwind scheme with the Roe_s flux-difference splitting is introduced into the inviscid flux. The artificial dissipation is also turned off in the region where the upwind scheme works. A binary sensor function Ui+1/2 at cell face i + 1/2 is used for the detection of shockwaves. Ui+1/2 is determined by the pressure and density curvature criteria proposed by Hill et al.21
89 90 91
93
94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124
iþ1
ap 2 ¼ maxðaip ; aiþ1 p Þ aiþ1=2 ¼ maxðaiq ; aiþ1 q q Þ piþ1 2pi þ pi1 ¼ piþ1 þ 2pi þ pi1 q 2qi þ qi1 aiq ¼ iþ1 q þ 2q þ q
ð2Þ 126 127
aip
iþ1
i
ð3Þ
i1
where aip and aiq represent the pressure p and density r relative curvatures at cell center, respectively. Uiþ1=2 is 1, when apiþ1=2 > c1 and aiþ1=2 > c2; but zero, otherwise. The 3D verq sion of this detection is used in the simulations. Similar to the treatment Hill et al.21, the values of c1 and c2 that proved to give the best results are chosen as 0.01. Based on this detection, the Roe’s 2nd-order upwind flux only operates at the cells in the vicinity of shock waves.
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
129 130 131 132 133 134 135 136 137
CJA 729 12 November 2016
No. of Pages 12
Transonic wing stall of a blended flying wing common research model CDES ¼ CDES1 F1 þ CDES2 ð1 F1 Þ
138
2.2. Turbulence modeling
139
DDES simulation is implemented in the present work for turbulence closure.22–24 The k-x based shear stress transport (SST) model was designed to give highly accurate predictions of the onset and the amount of flow separation under adverse pressure gradients by the inclusion of transport effects into the formulation of the eddy-viscosity. We can refer to the Ref.23 for details on the constants and the quantities involved. The original DES method depends on the mesh scale seriously. If the mesh is too intensive within the boundary layer, the RANS will prematurely switch to large eddy simulation (LES) inducing Reynolds modeled stress depletion (MSD) and grid induced separation (GIS).25 Focusing on this problem, Menter improved the SST-DES model to SST-DDES method.23 The governing equations of the SST-DDES model read as23 pffiffiffiffiffi @ðqkÞ þ r ðqUkÞ ¼ r ½ðl þ rk lt Þrk þ Pk q k3 =lDDES @t ð4Þ
140 141 142 143 144 145 146 147 148 149 150 151 152 153 154
156 157
159
@ðqxÞ þ r ðqUxÞ ¼ r ½ðl þ rx lt Þrx @t rk rx q þ a Pk bqx2 þ 2ð1 F1 Þqrx2 x lt
ð5Þ
160 162 163 164 165 166 167 168 169 171 172
174
a1 k lt ¼ q maxða1 x; F2 SÞ
177 178 180 181
183 184 185 186 188
ð6Þ
where k is the turbulent kinetic energy, U the velocity, x the specific dissipation rate, m the kinematic viscosity, t the time, l viscosity coefficient, lt the turbulence viscosity coefficient, Pk the turbulence production term, lDDES the DDES length scale, and S the magnitude of the strain rate tensor. F1 and F2 denote the SST blending functions which read as follows F1 ¼ tanhðarg41 Þ
ð7Þ
! ! pffiffiffi k 500m 4qrx2 k ; arg1 ¼ min max ; Cl xdw d2w x d2w CDkx
ð8Þ
175
CDkx
rk rx ; 1010 ¼ max 2qrx2 x
F2 ¼ tanh arg22 ! pffiffiffi 2 k 500m arg2 ¼ max ; Cl xdw d2w x
ð9Þ ð10Þ ð11Þ
Here dw is the distance to the nearest wall. The production term reads as follows: Pk ¼ minðlt S ; 10 Cl qkxÞ 2
ð12Þ
The DDES length scale reads as follows:
189 190 192
lDDES ¼ lRANS fd maxð0; lRANS lLES Þ
ð13Þ
193 195
lLES ¼ CDES hmax
ð14Þ
196
198
lRANS
pffiffiffi k ¼ Cl x
3
ð15Þ
ð16Þ
Here hmax is the maximum edge length of the cell. Finally, the empiric blending function fd is computed with the use of the following relations: fd ¼ 1 tanh½ðCD1 rd Þ
CD2
mt þ m qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rd ¼ j2 d2w 0:5ðS2 þ X2 Þ
ð17Þ
199 201 202 203 204 205 207 208
ð18Þ
Here X is the magnitude of vorticity tensor. The model constants read as follows: Cl = 0.09, k = 0.41, a1 = 0.31, CDES1 = 0.78, CDES2 = 0.61, CD1 = 20, CD2 = 3. All the constants are computed by a blend from the corresponding constants of the k-x model via a = a1F1 + a2 (1 F1). a1 = 5/9, a2 = 0.44, b1 = 0.075, b2 = 0.0828, rk1 = 0.85, rk2 = 1, rx1 = 0.5, rx2 = 0.856.
210 211 212 213 214 215 216 217
3. Model and simulation cases
218
3.1. Model
219
In order to meet the needs of future aircraft aerodynamic test and research, the relevant domestic institutions independently have designed a blended flying wing common research model, as a low-aspect ratio flying wing shape of general research platform.17 The basic geometry parameters of the low aspect ratio flying wing model are shown in Fig 1: the leading edge sweepback angle is 65°, the rear sweepback angle is ±47°, the length of the whole model cr is 15.32 m, the averaged aerodynamic chord length cref is 9.56 m, and the distance between the moment reference point and the leading edge is 6.9 m. To compare with the experiments, the model used for simulation and wind tunnel test scaled to 1/19.
220
3.2. Mesh accuracy and time step
232
Structured mesh is applied and far-field boundary condition is set 15 times of the root chord length to the model surface and non-reflection condition. Unsteady simulation was employed; Mach number is 0.9, and Reynolds number based on the averaged aerodynamic chord length is 9.6 106. Each case is calculated for 10,000 steps. The first 3000 steps are used to gain stable flow field while 3000–10,000 steps are the effective data for analysis. Before calculating physical problem by unsteady computation, at first the density of the mesh used should be considered if it is enough for the spatial flow field structures and whether the time step is proper for revealing unsteady fluctuation information. In this paper, the method proposed in Ref.22 is used to check the mesh independence and the time step. Power spectrum density (PSD) analysis is applied to analyzing the aerodynamic force and the moment time history. For a given amount of the mesh, the physical time step should be short enough to keep the frequency of fluctuations unchanged. The independence of mesh dense and time step is studied at the angle of attack of 14°. As shown in Fig. 2, three mesh types of different densities are adopted, namely coarse mesh (5 million), medium mesh (10 million), fine mesh (20 million), and the normalized time scale is taken as t* = 0.020, t* = 0.010 and t* = 0.005, where t ¼ DtUc 1 and St ¼ Ufc1 . Dt is real time
233
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
221 222 223 224 225 226 227 228 229 230 231
234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256
CJA 729 12 November 2016
No. of Pages 12
4
Y. Tao et al.
Fig. 1
257 258 259 260 261 262
scale is less than 0.010, the main frequency of the pitching moment almost remains invariable at St = 1.85. Thus t* = 0.010 is proper for the time normalization for medium mesh. The grid refinement factor is defined as22 1=d N1 r¼ ð19Þ N2
263
where N1 and N2 are the number of grid points in the finer and coarser of the two grids being compared respectively and d is the spatial dimension (three in our case). The grid convergence index (GCI) is defined as Fs f f1 GCI ¼ p 2 ð20Þ r f1
271
where Fs is a safety factor of 3.0, p is the order of accuracy for the spatial scheme employed, and f1, f2 are the finer and coarser grid global solution quantities (CL, CD, Cm) considered. The GCI values for the force and moment coefficients are provided in Table 1 for simulations at a = 14° and t* = 0.01. The GCI values are relatively low when comparing the three grids, and slightly decrease when comparing with the medium and fine grids. Through the GCI values we can see that the grid volume has little effect on the aerodynamic coefficients of the model, and the medium grids can afford the simulation.
278
4. Results and analysis
289
4.1. Comparison of simulation results of DDES and URANS
290
In order to compare DDES and URANS on the simulation ability of the transonic aerodynamic features of the low aspect ratio flying wing configuration, the SST turbulent modeling is applied to both DDES and URANS. Based on the results of simulation and experiment data17 as shown in Fig. 4, the lift increases linearly when the angle of attack is less than 4°, while evident vortex lift appears at the angle of attack higher than 4°. At the angle of attack smaller than 12°, both methods can
291
264 265 266 267 268
270
272 273 274 275 277
279 280 281 282 283 284 285 286 287 288
Basic geometry parameters of model.
step length, U1 is stream velocity, c is root chord length, and f is frequency of fore and moment. It takes 100 physical time steps for the free-stream flow to pass the root chord length when t* = 0.01. For medium mesh, the PSD of pitching moment time history of different time steps is analyzed as shown in Fig 3. When normalized time
Fig. 2
Computation grid.
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
292 293 294 295 296 297 298
CJA 729 12 November 2016
No. of Pages 12
Transonic wing stall of a blended flying wing common research model
Fig. 3 Variation of pitching moment power spectrum density with time step for medium density mesh. 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333
reveal the vortex lift well. However, when the angle of attack is higher than 14°, the lift drops down abruptly, namely abrupt wing stall, and after that it goes up gradually. This phenomenon can be revealed from the curves of drag and pitching moment, which is the abrupt drop of drag and the increase in nose-up pitching moment. It can be found from the results of the two methods that DDES can capture the abrupt wing stall and gain the aerodynamic force which matches the experimental data well. Although URANS reveals the abrupt wing stall too, the calculated stall incidence is 2° larger than the experiment results. Thus it can be gained that DDES is better than URANS to reveal the abrupt wing stall and has the advantage of simulating this kind of separation flow. The flow field structures at stall incidence of a = 14° are analyzed using Q criterion (namely the second invariant tensor of velocity gradient) to describe the vortex system structure. The instantaneous Q criterion isosurface calculated from these two methods is shown in Fig. 5. It can be seen from Fig. 5(a) that in the flow field calculated by URANS method, two main vortexes in the whole lee side did not break down until it was dissipated at the rear of the model. This corresponds to the aerodynamic force of non-abrupt wing stall calculated by URANS. Fig. 5(b) represents the vortex structures of the flow field calculated by DDES method. It can be found that the two main vortexes break down at the location around x/xcr = 0.6 on the lee side, which matches the location where the lift coefficient drops down abruptly, where x is streamwise location and xcr is chord length. The pressure coefficient is larger downstream of the vortex breakdown location than that upstream. Flow field simulated by DDES matches the experiment results qualitatively, and this method can explain the abrupt wing stall phenomenon reasonably. Thus the results analyzed below are gained by DDES method. The calculated streamline on the surface is compared with the oil visualization result, as shown in Fig. 6. Three typical
Table 1
5
states are selected and analyzed, which are a = 10°, 14°, 16° corresponding to the states before, at and after abrupt wing stall respectively. The black line in the oil visualization represents the separation line, while the red one means the reattachment line. Based on Fig. 6, the calculated results in most regions match the experimental data well. In different states, the streamlines on the surface change with the variation of the vortex spatial structures. When a = 10°, two evident reattachment lines and a separation line appear on the surface alternately, and there exists a separation line even at the leading edge. It can be revealed from the vortex counter that strong leading edge vortex appears on the lee side and induces the second vortex, which increases the suction and vortex lift. When a = 14°, an evident reattachment line and a separation line appear on the upper surface, and the leading edge. Evident stacking appears in the middle and tail section of the separation line, and based on the vortex counter, the leading edge vortex has broken up at this location where there is 13% body length ahead compared with that of a = 10°. Due to the breakdown of the leading edge vortex, the vortex lift of this region decreases dramatically which results in the lift coefficient decrease and positive pitching moment. When a = 16°, an evident reattachment line and a separation line appear on the upper surface and a separation line exists near the leading edge. The flow feature is similar to that of a = 14°, but the vortex broken location is ahead within small region. Compared with the a = 14° case, due to the development of leading edge vortex, the increase in vortex lift is more than the decrease caused by the vortex broken ahead, and thus the vortex lift and the lift coefficient of the whole vehicle go up with the increase in angle of attack.
CL
CD
Cm
Coarse Medium Fine
0.582 0.641 0.643
0.1392 0.1581 0.1578
0.1233 0.0928 0.0923
335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365
4.2. Vortex force component and role of vortex breakdown
366
The variation of the lift characteristic of the low aspect ratio flying wing configuration studied in this paper is directly related to the formation, development and breakdown of the leading edge vortex on the lee side. To view the variation of the vortex lift in detail, the leading edge suction analogy method proposed by Polhamus26 is applied to decomposing the lift into potential lift and vortex lift as follows:
367
CL ¼ CLp þ CLv
ð21Þ
where CLp indicates potential lift coefficient and CLv the vortex lift coefficient. Because the lift at zero-incidence is not zero due to the curve of the airfoil cross section of the studied flying wing configuration, the lift decomposition equation of Polhamus is revised in this paper that a constant CL0 which represents the zero-incidence lift coefficient is introduced into the equation to make it more applicable to the problem dis-
Force and moment coefficients for the three grids used in convergence study at a = 14°.
Grid
334
GCI CL
CD
Cm
0.496204 0.015272
0.664591 0.00929
1.34035 0.21876
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
368 369 370 371 372 373 374 376 377 378 379 380 381 382 383
CJA 729 12 November 2016
No. of Pages 12
6
Y. Tao et al.
Fig. 4
Fig. 5
Comparison of aerodynamic coefficients.
Instantaneous Q criterion isosurface (a = 14°, Q = 10).
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
CJA 729 12 November 2016
No. of Pages 12
Transonic wing stall of a blended flying wing common research model
Fig. 6 384 385 386 388 389 390 391 392 393 394 395 397 398 399 400
Comparison of calculated streamline on surface and oil visualization results.
cussed herein. Thus the final equation for potential lift coefficient is CLp ¼ Kp sin a cos a þ CL0 2
ð22Þ
where Kp is the lift-curve slope given by small-angle-of-attack potential-flow lifting surface theory, and accounts for the true boundary condition, and cos2 a arises from the assumption of a Kutta-type flow condition at the leading edge. While the equation for the vortex lift is used by the result of typical leading edge suction analogy method. CLv ¼ Kv cos a sin2 a
7
ð23Þ
where Kv sin2 a represents the potential-flow leading-edge suction, therefore the vortex normal force, and cos a indicates the component in the lift direction.
The flying wing model used here has a low aspect ratio equal to 1.54, and the leading edge sweepback angle is 65°, the rear sweepback angle is ±47°, =65°, aspect ratio 1.54. For those parameters, Kp = 1.9 and Kv = 2.7 can be obtained with Polhamus approach26 based on the Variations figure of Kp and Kv for delta wings provided by Ref.26. Because the Polhamus approach is suitable for incompressible flow mainly, a correction should be made for the compressible effect. According to Glauert rule, compressibility has no effect on the aerodynamic force of delta wing, when the Mach number is lower enough. C0L
CL ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 Ma21
401 402 403 404 405 406 407 408 409 410 411 412
ð24Þ
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
414
CJA 729 12 November 2016
No. of Pages 12
8 415 416 417
419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434
Y. Tao et al.
Take the assumption that CLp and CLv obey Glauert rule too. C0Lp
CLp ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; 1 Ma21
C0Lv
CLv ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 Ma21
ð25Þ
So Kp ¼ 4:3589 and Kv ¼ 6:1942 at Ma1 ¼ 0:9. Fig. 7 shows the comparison of lift coefficients predicted by Polhamus approach and Glauert rule and obtained from experiments. The predicted lift is much higher than the experiment results. So for transonic flow, this does not hold true in the case of such high number flow. Because the traditional Polhamus approach and Glauert rule have low precision to decompose vortex force component at transonic flow, the least square method (LSM) was employed to calibrate the Kp and Kv in this study. Lift coefficient curve shows that the slope is basically linear at angle of attack ranging from 2° to 4°. In this range, the least square method is used to fit the Eq. (22) and then Kp and CL0 were calibrated as at Ma1 ¼ 0:9. Then the potential lift coefficient curve CLp can be achieved as shown in Fig 7.
Fig. 7
Vortex lift coefficients for blending wing configuration.
The true vortex lift coefficient curve named as CL CLp curve can be obtained by Coefficient CL minus the potential lift coefficient CLp . In this range from 4° to 12°, the least square method is used to fit the Eq. (23) and then Kv was calibrated. In this case Kv equals to 3.149; the detail results are shown in Fig. 7. From the fitting results of the vortex lift coefficient curve, the fitting curve does not match the original curve very well but it still represents its main characteristics. Since the lift coefficients of CFD match well with the experimental result at low angle of attack, CLpDDES was replaced by CLpLSM . The vortex lift coefficients of DDES method were computed by CLDDES CLpLSM as shown in Fig. 7. The vortex lift coefficient calculated by the method above is shown in Fig. 7 where it can be found that the vortex lift remains zero when the angle of attack is less than 4° meaning the leading edge vortex has not formed yet; however, the vortex lift increases dramatically when the angle of attack is more than 4°. With the increase in angle of attack, the vortex lift coefficient suffers abrupt drop down when the angle of attack is more than 12°, which causes the abrupt decrease in the lift coefficient of flying wing configuration model. To explain the relation between the abrupt decrease and the vortex breakdown and the mechanism of vortex breakdown, the detail of the flow field is investigated. Fig. 8 shows the pressure distribution p/pref through vortex cores at a = 12, 14°, which indicates that shocks/vortex interactions do exist in the flow. These shocks/vortex interactions influence the flow structures obviously especially at larger angles of incidence. At a = 12°, two adverse pressure gradients at location of x=xcr ¼ 0:33 and x=xcr ¼ 0:77 affect the vortex strength to a different extent. At the first location, the adverse pressure gradient is low ((@p=@xaxial ¼ 0:43), which has a small influence on the vortex stability. However, at the second location where a strong normal shock exists, the pressure through vortex cores has a sudden rise and the adverse pressure gradient is high ((@p=@xaxial ¼ 4:17), which indicate that strong shocks/vortex interactions exist in these areas. At a = 14°, two adverse pressure gradients also exist in the flow field. Compared to that of a = 12°, the location of the second adverse pressure gradient moves forward to x=xcr ¼ 0:58, which is closer to the head of the configuration and the (@p=@xaxial increases to 7.26. Due to the high adverse pressure gradient, the main leading-edge vortex breakdown is homologous with the Q criterion isosurface (Fig. 9).
Fig. 8
Pressure distribution through vortex cores.
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479
CJA 729 12 November 2016
No. of Pages 12
Transonic wing stall of a blended flying wing common research model 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504
From the CFD results of flow structure, there are shocks/ vortex interactions at both lower incident and higher incidence, with a stronger interaction for higher incidence, and the interaction causes a considerable weakening of the vortex core, which may result in vortex breakdown. There seems to be a criterion for the vortex to feel the effect of the shock and remain coherent or breakdown. Ref.25 demonstrated an importance criterion parameter for vortex breakdown caused by shocks/vortex interaction, which is a combination of tangential velocity Uh and axial velocity Uaxial of the vortex core. The swirl ratio or the Rossby number is a nondimensional parameter defined as the ratio of axial velocity and tangential velocity of vortex core, as defined by Ro ¼ Uaxial =Uh . It is used as a criterion to measure the vortex intensity and the susceptibility of the vortex to shock-induced breakdown. The maximum swirl velocity of the vortex and the maximum axial velocity of the vortex cores are used for the calculation of Rossby number in this investigation. When the leading-edge vortex passes through a normal shock, the swirl velocity changes little, while the axial velocity will decrease distinctly. As a result, the Rossby number will decrease and the stability of the vortex decrease makes the vortex easier to breakdown. Spall and robinson found that if Ro > 1.4, the vortex is stable, and if 0.9 < Ro < 1.4, the vortex is unstable, while if Ro < 0.9, the vortex breaks down
Fig. 9
Instantaneous Q criterion isosurface (Q = 10).
9
when Rossby number was used as the criterion parameter for leading-edge vortex breakdown studies of delta wings. Using the criterion, the Rossby numbers at a = 12° and 14° were calculated which is shown in Fig. 10. From the result, it is clear that shock has a great influence on the vortex stabilities. For a = 12°, at the locations of x=xcr ¼ 0:33 and x=xcr ¼ 0:5, the Rossby numbers without exception decrease obviously. As the shocks and shocks/vortex interactions are weak at these locations, the vortex keeps stable even though the vortex strength decreases. When the vortex passes through the second adverse pressure gradient, shocks/vortex interactions become stronger and the vortex becomes unstable. Subsequently, the Rossby number decreases gradually to 0.9 at the location of x=xcr ¼ 0:93, which means the vortex is broken down here. At a = 14°, shocks/vortex interactions are stronger compared to those of a = 12°, Rossby number reduces rapidly at the location of x=xcr ¼ 0:57, and the vortex becomes unstable and then breaks down. Fig. 11 shows the shock structures in the flow field at a = 12° and 14°; due to the high curve slope at the leading edge, the flow accelerates to a local peak value and makes a weak shock produced at the location of x=xcr ¼ 0:25 at a = 12°, and a strong normal shock can be observed around the location of x=xcr ¼ 0:75. Near the central span-wise section of the configuration, two cross flow shocks appear as the interactions of the main leading-edge vortex and second vortex at this area. As the angle of incidence increases to 14°, the weak shock almost disappears at the location of x=xcr = 0.25 compared to a = 12°, while two new normal shocks introduce at the location around xxcr ¼ 0:50 0:55 and the vortex breaks down as it passes through the normal shocks and the cross flow shocks destructs too at this angle of attack. Fig. 12 shows pressure coefficient distribution at the symmetry plane on the wing for both angles of incidence. At a = 12°, around the location of x=xcr ¼ 0:75, a normal shock is there. As angles of incidence increase to 14°, this normal shock moves to x=xcr ¼ 0:55 which makes the vortex breaks down after the vortex passes through this shock. At the rear of the configuration, weak shocks can be observed for both angles of incidence due to the high curve slope. Through comparison, we found that this model had a less breakdown resistant to angles of incidence compared with
Fig. 10
Rossby number distribution against root chord.
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546
CJA 729 12 November 2016
No. of Pages 12
10
Y. Tao et al.
Fig. 12 Pressure coefficient distribution at symmetry plane on the wing for both angles of incidence.
Fig. 11
547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569
Shock structures in flow field.
the other 65° swept delta wing configurations proposed in Ref.3,6 In transonic flow field, the wing-body thickness ratio has great influence on the shock behavior on delta wings and the occurrence of shocks introduces complex shock/vortex interactions which are responsible for vortex breakdown. In addition, almost all the published investigation models with flying-wing configurations have a relatively small wing-body thickness ratio, and the flow differences around flying-wing configuration with a moderate and a small wing-body thickness ratio are not entirely described and understood. Based on the previous study of the 65° swept flying-wing configuration with a moderate 0.16 wing-body thickness ratio (marked as Body-orig), through reducing the wing-body thickness (marked as Body-thin) and keeping the round leading edge and outer wing geometry identical, we investigated wingbody thickness effects on the aerodynamic and vortex flow characteristics of the flying-wing configuration at Ma = 0.9. Fig. 13 shows comparisons of typical cross-sectional shapes between the original wing-body thickness ratio model and the reduced wing-body thickness ratio model. Fig. 14 shows comparisons of lift and pitching moment coefficients versus incidence relative to the ‘original’ and ‘reduced’ models. For a < 12°, effect of inner wing thickness
Fig. 13
Typical cross-sectional shapes of two configurations.
on lift coefficient is very weak. However, as angle of incidence increases, the lift coefficient of ‘original’ model becomes relatively flat up to 15° incidence, whereas the lift coefficient of ‘thin’ model has a linear increase up to 20°. A similar phenomenon occurs in the pitching moment. Abrupt nose up in pitching moment of ‘reduced’ model is at 20°, while the ‘original’ one at 12°. That is to say ‘reduced’ model can delay stall and vortex breakdown by an angle of 8° compared with the original model at Mach number 0.9. From the pressure distribution through vortex cores at a = 14° in Fig. 15, it is clear that accelerating pressure exists up to x=xcr ¼ 0:65 of the ‘reduced’ model, and the location of abrupt pressure rise is far from the head of the configuration compared to the ‘original’ model, which explains why the ‘reduced’ model has a larger stall angle of incidence. From this investigation we can see that wing-body thickness has great influence on the introduction of shock and shock/vortex interactions which can make a significant difference to the vortex breakdown behavior and aerodynamic characteristic of the flying-wing configuration, which should be considered in the aerodynamic design process.
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590
CJA 729 12 November 2016
No. of Pages 12
Transonic wing stall of a blended flying wing common research model
11
5. Conclusions
591
Numerical simulation of abrupt stall phenomena of blended flying wing configuration at transonic speed was conducted using both DDES and URANS methods at transonic speed. Based on the analysis of the aerodynamic coefficients, flow structures and streamlines, the following conclusions can be drawn:
592
(1) The sudden motion in breakdown location observed in mean lift, drag and pitching moment coefficients of experiments is due to a shock/vortex interaction. (2) The DDES methods perform more precisely than the URANS method (the onset angle of the breakdown movement was predicted about 2° greater than the measurements) and the aerodynamic coefficient results from DDES method were compared very well with the experiment results and the shock strength or axial flow in the vortex was well predicted by the DDES methods. (3) Wing-body thickness has great influence on shock and shock/vortex interactions, which can make a significant difference to the vortex breakdown behavior and stall characteristic of the blended flying wing configuration.
598
593 594 595 596 597
599 600 601 602 603 604 605 606 607 608 609 610 611 612
Fig. 14 Comparisons of lift and pitching moment coefficients vs incidence.
Fig. 15 cores.
Comparison of pressure distribution through vortex
Acknowledgment
613
This work was funded by the National Natural Science Foundation of China (No. 11372337).
614
References
616
1. Schiavetta LA, Boelens OJ, Fritz W. Analysis of transonic flow on a slender delta wing using CFD. Reston: AIAA; 2006, Report No.: AIAA-2006-3171. 2. Elsennnr A, Hoeijmakers HWM. An experimental study of the flow over a sharp-edged delta wing at subsonic and transonic speeds. Paris: AGARD; 1991, Report No.: AGARD-CP-494. 3. Chu J, Luckring JM. Experimental surface pressure data obtained on a 65° delta wing across Reynolds number and Mach number ranges: Volume 1-sharp leading edge.. Washington (DC): NASA; 1996, Report No.: NASA TM 4645. 4. Hill DJ, Pantano C, Pullin DI. Large-eddy simulation and multiscale modeling of a Richtmyer-Meshkov instability with reshock. J Fluid Mech 2006;557:29–61. 5. Fang BR. Aero-configuration design of aero plan. Beijing: Aviation Industry Press; 1997 [in Chinese]. 6. Jobe CE. Vortex breakdown location over 65° delta wings empiricism and experiment. Aeronaut J 2004;108:475–82. 7. Kalkhoran IM, Smart MK. Aspects of shock wave-induced vortex breakdown. Progress Aerospace Sci 2000;36(1):63–95. 8. Fritz W, Cummings RM. What was learned from the numerical simulations for the VFE-2. Reston: AIAA; 2008, Report No.: AIAA-2008-0399. 9. Schiavetta L, Boelens O, Crippa S. Sting effects on vortex breakdown in transonic delta wing flows. Reston: AIAA; 2008, Report No.: AIAA-2008-0395. 10. Schiavetta LA, Boelens OJ, Crippa S. Shock effects on delta wing vortex breakdown. J Aircraft 2009;46(3):903–14. 11. Spalart PR, Jou WH, Strelets M, Allmaras SR. Comments on the feasibility of les for wings, and on a hybrid rans/les approach, advances in DNS/LESProceedings of the 1st AFOSR Int. Conf. on DNS/LES. Columbus (OH): Greyden Press; 1997.
617
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
615
618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647
CJA 729 12 November 2016
12 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672
12. Viswanathan A, Squires KD, Forsythe JR. Detached-eddy simulation around a forebody at high angle of attack. Reston: AIAA; 2003, Report No.: AIAA-2003-0263. 13. Lv ZY, Zhu LG. Study on forms of vortex breakdown over delta wing. Chin J Aeronaut 2004;17(1):13–6. 14. Forsythe JR, Squires KD, Wurtzler KE, Spalart PR. Detachededdy simulation of the F-15E at high alpha. Reston: AIAA; 2002, Report No.: AIAA-2002-0591. 15. Morton SA, Steenman MB, Cummings RM, Forsythe JR. DES grid resolution issues for vertical flows on a delta wing and an F18C. Reston: AIAA; 2003, Report No.: AIAA-2003-1103. 16. Denning RM, Allen JE, Armstrong FW. The broad delta airliner. Aeronaut J 2003;107(1075):547–58. 17. Su JC, Huang Y, Li YH. Research on flow characteristics of lowaspect-ratio flying wing at transonic speed. Acta Aerodynamic Sinica 2015;33(3):307–12 [in Chinese]. 18. Simon F, Deck S, Guillen P, Sagaut P, Merlen A. Numerical simulation of the compressible mixing layer past an axisymmetric trailing edge. J Fluid Mech 2007;591:215–53. 19. Lu XY, Wang SW, Sung HG. Large-eddy simulations of turbulent swirling flows injected into a bump chamber. J Fluid Mech 2005;527:171–95. 20. Wang SW, Yang V, Hsiao G, Hsieh SY. Large-eddy simulations of gas-turbine swirl injector flow dynamics. J Fluid Mech 2007;583:99–122.
No. of Pages 12
Y. Tao et al. 21. Hill DJ, Pantano C, Pullin DI. Large-eddy simulation and multi scale modeling of a Richtmyer-Meshkov instability with reshock. J Fluid Mech 2006;557:29–61. 22. Menter FR, Kuntz M, Langtry R. Ten years of experience with the SST turbulence model. In: Proceedings of the 4th international symposium on turbulence heat and mass transfer; 2003. p. 625_32. 23. Menter FR. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J 1994;32(8):1598–605. 24. Jian Liu, Haisheng Sun, Zhitao Liu, Zhixiang Xiao. Numerical investigation of unsteady vortex breakdown past 80°/65°doubledelta wing. Chin J Aeronaut 2014;27(3):521–30. 25. Spalart PR. Detached-eddy simulation. Annu Rev Fluid Mech 2009;41:181–202. 26. Polhamus EC. Application of the leading edge suction analogy of vertex lift to the drag due to lift of sharp edge delta wings. Washington (DC): NASA; 1968, Report No.: TN-D-4739. Tao Yang is an associate research fellow in China Aerodynamics Research and Development Center. His main research interests are vortex aerodynamics, unsteady flows and their control. Zhao Zhongliang is a research fellow in China Aerodynamics Research and Development Center. His main research interests are unsteady flows, virtual flight testing and flight control.
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697
No. of Pages 12
1
Chinese Society of Aeronautics and Astronautics & Beihang University
Chinese Journal of Aeronautics [email protected] www.sciencedirect.com
4
Transonic wing stall of a blended flying wing common research model based on DDES method
5
Tao Yang a,b,*, Li Yonghong b, Zhang Zhao b, Zhao Zhongliang b, Liu Zhiyong b
3
6 7
8
a b
State Key Laboratory of Aerodynamics, China Aerodynamics Research and Development Center, Mianyang 621000, China High Speed Aerodynamic Institute, China Aerodynamics Research and Development Center, Mianyang 621000, China
Received 29 September 2015; revised 24 December 2015; accepted 31 January 2016
9
11 12
KEYWORDS
13
Delayed detached eddy simulation; Flying wing; Vortex lift; Vortex breakdown; Wing stall
14 15 16 17 18
19
Abstract Numerical simulation of wing stall of a blended flying wing configuration at transonic speed was conducted using both delayed detached eddy simulation (DDES) and unsteady Reynolds-averaged Navier-Stokes (URANS) equations methods based on the shear stress transport (SST) turbulence model for a free-stream Mach number 0.9 and a Reynolds number 9.6 106. A joint time step/grid density study is performed based on power spectrum density (PSD) analysis of the frequency content of forces or moments, and medium mesh and the normalized time scale 0.010 were suggested for this simulation. The simulation results show that the DDES methods perform more precisely than the URANS method and the aerodynamic coefficient results from DDES method compare very well with the experiment data. The angle of attack of nonlinear vortex lift and abrupt wing stall of DDES results compare well with the experimental data. The flow structure of the DDES computation shows that the wing stall is caused mainly by the leeward vortex breakdown which occurred at x/xcr = 0.6 at angle of attack of 14°. The DDES methods show advantage in the simulation problem with separation flow. The computed result shows that a shock/vortex interaction is responsible for the wing stall caused by the vortex breakdown. The balance of the vortex strength and axial flow, and the shock strength, is examined to provide an explanation of the sensitivity of the breakdown location. Wing body thickness has a great influence on shock and shock/ vortex interactions, which can make a significant difference to the vortex breakdown behavior and stall characteristic of the blended flying wing configuration. Ó 2016 Production and hosting by Elsevier Ltd. on behalf of Chinese Society of Aeronautics and Astronautics. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/ licenses/by-nc-nd/4.0/).
* Corresponding author at: State Key Laboratory of Aerodynamics, China Aerodynamics Research and Development Center, Mianyang 621000, China. E-mail address: [email protected] (Y. Tao). Peer review under responsibility of Editorial Committee of CJA.
Production and hosting by Elsevier
1. Introduction
20
The development of studying and understanding the phenomenon of vortex lift has greatly impacted the aerodynamic configuration of the modern fighters, and the aerodynamic configuration using delta wing with large sweep angle has become one of the research topics.1–3 This kind of aerodynamic layout has the evident advantages of maneu-
21
http://dx.doi.org/10.1016/j.cja.2016.10.018 1000-9361 Ó 2016 Production and hosting by Elsevier Ltd. on behalf of Chinese Society of Aeronautics and Astronautics. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
22 23 24 25 26
CJA 729 12 November 2016
No. of Pages 12
2 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85
verability, agility and stall incidence increase.4,5 The occurrence of shocks on delta wings introduces complex shock/ vortex interactions, particularly at moderate to high angles of incidence. These interactions can make a significant difference to the vortex breakdown behavior. With the increase in angle of attack, flight vehicles would suffer abrupt wing stall due to the condition change such as adverse pressure gradient. The reason for abrupt wing stall is nonlinear lift loss caused by the leeward vortex breakdown, and the mechanism and impact factors are very complicated.5–7 The strengthening of the shock which stands off the sting as the incidence is increased can lead to a shock/vortex interaction triggering breakdown. The location of breakdown can shift upstream by as much as 30% of the chord at a single 1° incidence interval2,3 due to this interaction. Much literature has reported studies of this problem via theoretical analysis, experiments and numerical computations.8 Because the wing stall and vortex breakdown phenomenon involve the computation of flow separation, the accuracy of unsteady Reynolds-averaged Navier-Stokes (URANS) equations method is insufficient and thus causes the deviation between the calculated stall angle and experimental data.9,10 As a hybrid method applied to simulating this kind of separation flow, delayed detached eddy simulation (DDES) takes the advantages of higher precision and acceptable computation consumption. It has been applied to many problems, such as shock vibration of transonic wing, cavity flow, and wing separation flow.11–15 Aerodynamic efficiency defined as MaCL/CD has been steadily increasing from 13 in 1960 to 16 in the mid-1990s16, where Ma is stream velocity Mach number, CL is lift coefficient and CD is drag coefficient. Higher aerodynamic efficiency will benefit the flight range and fuel economy, and thus aircraft designers pursue higher cruising speed as well as lower drag and higher lift to drag ratios. As an example the cruising speed increases from Ma = 0.78 for A320 to about Ma = 0.85 for A380. Blended-wing-body (BWB) configuration takes an advantage of lift to drag ratios compared with traditional configurations. The developing BWB aircraft of Boeing Company aims at a design cruising Ma = 0.85. The low aspect ratio blended flying wing configuration has more potential on increasing cruising speed. The most typical example is X47b which is defined as a cruising speed at high subsonic speed and some literature explicitly refers to it as about Ma = 0.9. The measured data from 2.4 m 2.4 m transonic wind tunnel at transonic conditions 17 (Ma = 0.8 and 0.9) showed a sudden jump of the breakdown location toward the wing apex of a blended flying wing common research model when a critical angle of incidence was reached. All the computation study was conducted at Ma = 0.9 because the flow condition at this Mach number was more similar to the design cruising speed and the shock wave/ vortex interaction phenomenon was more typical. DDES model was employed to study the transonic stall problem, and by comparing with URANS and experimental data, the accuracy and applicability were validated. The flow structures and the role of vortex breakdown at stall angles of attack were briefly analyzed, too.
Y. Tao et al. 2. Mathematical formulation and numerical method
86
2.1. Governing equations and numerical procedure
87
The governing equations, which are the three-dimensional time-dependent compressible Navier-Stokes equations, can be written in terms of generalized coordinates as
88
^ @ E^ @ F^ @ G^ @ E^V @ F^V @ G^V @Q þ þ þ þ þ ¼ @n @g @f @t @n @g @f
ð1Þ
b is the vector of conserved variables, density, momenwhere Q b are inviscid flux b F, b G tum and total energy per unit volume, E, b b b terms and, E V , F V , G V are viscous flux terms. A general, 3D transformation between the Cartesian variables (x, y, z) and the generalized coordinated (n, g, f) is implied. The equation of state for an ideal gas is used and the molecular viscosity is assumed to obey Sutherland’s law. To nondimensionalize the equations, we use the free-stream variables including the density, temperature, speed of sound, and root chord of the blended wing as characteristic scales. The initial and boundary conditions are presented as follows. The initial condition is set as the free-stream quantities. The far field boundary conditions are treated by local onedimensional Riemann-invariants. No-slip and adiabatic conditions are applied on the aircraft surface. The governing Eq. (1) is discretized using finite volume method and the time term is discretized by an implicit approximate-factorization method with sub-iterations to ensure the second-order accuracy Simon et al.18. In this paper, the convective terms were adopted by second-order central difference scheme and a fourth-order artificial viscous term was introduced to restrain the numerical oscillation Lu et al.19 and Wang et al.20. To capture the discontinuity caused by shockwave, a second-order upwind scheme with the Roe_s flux-difference splitting is introduced into the inviscid flux. The artificial dissipation is also turned off in the region where the upwind scheme works. A binary sensor function Ui+1/2 at cell face i + 1/2 is used for the detection of shockwaves. Ui+1/2 is determined by the pressure and density curvature criteria proposed by Hill et al.21
89 90 91
93
94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124
iþ1
ap 2 ¼ maxðaip ; aiþ1 p Þ aiþ1=2 ¼ maxðaiq ; aiþ1 q q Þ piþ1 2pi þ pi1 ¼ piþ1 þ 2pi þ pi1 q 2qi þ qi1 aiq ¼ iþ1 q þ 2q þ q
ð2Þ 126 127
aip
iþ1
i
ð3Þ
i1
where aip and aiq represent the pressure p and density r relative curvatures at cell center, respectively. Uiþ1=2 is 1, when apiþ1=2 > c1 and aiþ1=2 > c2; but zero, otherwise. The 3D verq sion of this detection is used in the simulations. Similar to the treatment Hill et al.21, the values of c1 and c2 that proved to give the best results are chosen as 0.01. Based on this detection, the Roe’s 2nd-order upwind flux only operates at the cells in the vicinity of shock waves.
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
129 130 131 132 133 134 135 136 137
CJA 729 12 November 2016
No. of Pages 12
Transonic wing stall of a blended flying wing common research model CDES ¼ CDES1 F1 þ CDES2 ð1 F1 Þ
138
2.2. Turbulence modeling
139
DDES simulation is implemented in the present work for turbulence closure.22–24 The k-x based shear stress transport (SST) model was designed to give highly accurate predictions of the onset and the amount of flow separation under adverse pressure gradients by the inclusion of transport effects into the formulation of the eddy-viscosity. We can refer to the Ref.23 for details on the constants and the quantities involved. The original DES method depends on the mesh scale seriously. If the mesh is too intensive within the boundary layer, the RANS will prematurely switch to large eddy simulation (LES) inducing Reynolds modeled stress depletion (MSD) and grid induced separation (GIS).25 Focusing on this problem, Menter improved the SST-DES model to SST-DDES method.23 The governing equations of the SST-DDES model read as23 pffiffiffiffiffi @ðqkÞ þ r ðqUkÞ ¼ r ½ðl þ rk lt Þrk þ Pk q k3 =lDDES @t ð4Þ
140 141 142 143 144 145 146 147 148 149 150 151 152 153 154
156 157
159
@ðqxÞ þ r ðqUxÞ ¼ r ½ðl þ rx lt Þrx @t rk rx q þ a Pk bqx2 þ 2ð1 F1 Þqrx2 x lt
ð5Þ
160 162 163 164 165 166 167 168 169 171 172
174
a1 k lt ¼ q maxða1 x; F2 SÞ
177 178 180 181
183 184 185 186 188
ð6Þ
where k is the turbulent kinetic energy, U the velocity, x the specific dissipation rate, m the kinematic viscosity, t the time, l viscosity coefficient, lt the turbulence viscosity coefficient, Pk the turbulence production term, lDDES the DDES length scale, and S the magnitude of the strain rate tensor. F1 and F2 denote the SST blending functions which read as follows F1 ¼ tanhðarg41 Þ
ð7Þ
! ! pffiffiffi k 500m 4qrx2 k ; arg1 ¼ min max ; Cl xdw d2w x d2w CDkx
ð8Þ
175
CDkx
rk rx ; 1010 ¼ max 2qrx2 x
F2 ¼ tanh arg22 ! pffiffiffi 2 k 500m arg2 ¼ max ; Cl xdw d2w x
ð9Þ ð10Þ ð11Þ
Here dw is the distance to the nearest wall. The production term reads as follows: Pk ¼ minðlt S ; 10 Cl qkxÞ 2
ð12Þ
The DDES length scale reads as follows:
189 190 192
lDDES ¼ lRANS fd maxð0; lRANS lLES Þ
ð13Þ
193 195
lLES ¼ CDES hmax
ð14Þ
196
198
lRANS
pffiffiffi k ¼ Cl x
3
ð15Þ
ð16Þ
Here hmax is the maximum edge length of the cell. Finally, the empiric blending function fd is computed with the use of the following relations: fd ¼ 1 tanh½ðCD1 rd Þ
CD2
mt þ m qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rd ¼ j2 d2w 0:5ðS2 þ X2 Þ
ð17Þ
199 201 202 203 204 205 207 208
ð18Þ
Here X is the magnitude of vorticity tensor. The model constants read as follows: Cl = 0.09, k = 0.41, a1 = 0.31, CDES1 = 0.78, CDES2 = 0.61, CD1 = 20, CD2 = 3. All the constants are computed by a blend from the corresponding constants of the k-x model via a = a1F1 + a2 (1 F1). a1 = 5/9, a2 = 0.44, b1 = 0.075, b2 = 0.0828, rk1 = 0.85, rk2 = 1, rx1 = 0.5, rx2 = 0.856.
210 211 212 213 214 215 216 217
3. Model and simulation cases
218
3.1. Model
219
In order to meet the needs of future aircraft aerodynamic test and research, the relevant domestic institutions independently have designed a blended flying wing common research model, as a low-aspect ratio flying wing shape of general research platform.17 The basic geometry parameters of the low aspect ratio flying wing model are shown in Fig 1: the leading edge sweepback angle is 65°, the rear sweepback angle is ±47°, the length of the whole model cr is 15.32 m, the averaged aerodynamic chord length cref is 9.56 m, and the distance between the moment reference point and the leading edge is 6.9 m. To compare with the experiments, the model used for simulation and wind tunnel test scaled to 1/19.
220
3.2. Mesh accuracy and time step
232
Structured mesh is applied and far-field boundary condition is set 15 times of the root chord length to the model surface and non-reflection condition. Unsteady simulation was employed; Mach number is 0.9, and Reynolds number based on the averaged aerodynamic chord length is 9.6 106. Each case is calculated for 10,000 steps. The first 3000 steps are used to gain stable flow field while 3000–10,000 steps are the effective data for analysis. Before calculating physical problem by unsteady computation, at first the density of the mesh used should be considered if it is enough for the spatial flow field structures and whether the time step is proper for revealing unsteady fluctuation information. In this paper, the method proposed in Ref.22 is used to check the mesh independence and the time step. Power spectrum density (PSD) analysis is applied to analyzing the aerodynamic force and the moment time history. For a given amount of the mesh, the physical time step should be short enough to keep the frequency of fluctuations unchanged. The independence of mesh dense and time step is studied at the angle of attack of 14°. As shown in Fig. 2, three mesh types of different densities are adopted, namely coarse mesh (5 million), medium mesh (10 million), fine mesh (20 million), and the normalized time scale is taken as t* = 0.020, t* = 0.010 and t* = 0.005, where t ¼ DtUc 1 and St ¼ Ufc1 . Dt is real time
233
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
221 222 223 224 225 226 227 228 229 230 231
234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256
CJA 729 12 November 2016
No. of Pages 12
4
Y. Tao et al.
Fig. 1
257 258 259 260 261 262
scale is less than 0.010, the main frequency of the pitching moment almost remains invariable at St = 1.85. Thus t* = 0.010 is proper for the time normalization for medium mesh. The grid refinement factor is defined as22 1=d N1 r¼ ð19Þ N2
263
where N1 and N2 are the number of grid points in the finer and coarser of the two grids being compared respectively and d is the spatial dimension (three in our case). The grid convergence index (GCI) is defined as Fs f f1 GCI ¼ p 2 ð20Þ r f1
271
where Fs is a safety factor of 3.0, p is the order of accuracy for the spatial scheme employed, and f1, f2 are the finer and coarser grid global solution quantities (CL, CD, Cm) considered. The GCI values for the force and moment coefficients are provided in Table 1 for simulations at a = 14° and t* = 0.01. The GCI values are relatively low when comparing the three grids, and slightly decrease when comparing with the medium and fine grids. Through the GCI values we can see that the grid volume has little effect on the aerodynamic coefficients of the model, and the medium grids can afford the simulation.
278
4. Results and analysis
289
4.1. Comparison of simulation results of DDES and URANS
290
In order to compare DDES and URANS on the simulation ability of the transonic aerodynamic features of the low aspect ratio flying wing configuration, the SST turbulent modeling is applied to both DDES and URANS. Based on the results of simulation and experiment data17 as shown in Fig. 4, the lift increases linearly when the angle of attack is less than 4°, while evident vortex lift appears at the angle of attack higher than 4°. At the angle of attack smaller than 12°, both methods can
291
264 265 266 267 268
270
272 273 274 275 277
279 280 281 282 283 284 285 286 287 288
Basic geometry parameters of model.
step length, U1 is stream velocity, c is root chord length, and f is frequency of fore and moment. It takes 100 physical time steps for the free-stream flow to pass the root chord length when t* = 0.01. For medium mesh, the PSD of pitching moment time history of different time steps is analyzed as shown in Fig 3. When normalized time
Fig. 2
Computation grid.
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
292 293 294 295 296 297 298
CJA 729 12 November 2016
No. of Pages 12
Transonic wing stall of a blended flying wing common research model
Fig. 3 Variation of pitching moment power spectrum density with time step for medium density mesh. 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333
reveal the vortex lift well. However, when the angle of attack is higher than 14°, the lift drops down abruptly, namely abrupt wing stall, and after that it goes up gradually. This phenomenon can be revealed from the curves of drag and pitching moment, which is the abrupt drop of drag and the increase in nose-up pitching moment. It can be found from the results of the two methods that DDES can capture the abrupt wing stall and gain the aerodynamic force which matches the experimental data well. Although URANS reveals the abrupt wing stall too, the calculated stall incidence is 2° larger than the experiment results. Thus it can be gained that DDES is better than URANS to reveal the abrupt wing stall and has the advantage of simulating this kind of separation flow. The flow field structures at stall incidence of a = 14° are analyzed using Q criterion (namely the second invariant tensor of velocity gradient) to describe the vortex system structure. The instantaneous Q criterion isosurface calculated from these two methods is shown in Fig. 5. It can be seen from Fig. 5(a) that in the flow field calculated by URANS method, two main vortexes in the whole lee side did not break down until it was dissipated at the rear of the model. This corresponds to the aerodynamic force of non-abrupt wing stall calculated by URANS. Fig. 5(b) represents the vortex structures of the flow field calculated by DDES method. It can be found that the two main vortexes break down at the location around x/xcr = 0.6 on the lee side, which matches the location where the lift coefficient drops down abruptly, where x is streamwise location and xcr is chord length. The pressure coefficient is larger downstream of the vortex breakdown location than that upstream. Flow field simulated by DDES matches the experiment results qualitatively, and this method can explain the abrupt wing stall phenomenon reasonably. Thus the results analyzed below are gained by DDES method. The calculated streamline on the surface is compared with the oil visualization result, as shown in Fig. 6. Three typical
Table 1
5
states are selected and analyzed, which are a = 10°, 14°, 16° corresponding to the states before, at and after abrupt wing stall respectively. The black line in the oil visualization represents the separation line, while the red one means the reattachment line. Based on Fig. 6, the calculated results in most regions match the experimental data well. In different states, the streamlines on the surface change with the variation of the vortex spatial structures. When a = 10°, two evident reattachment lines and a separation line appear on the surface alternately, and there exists a separation line even at the leading edge. It can be revealed from the vortex counter that strong leading edge vortex appears on the lee side and induces the second vortex, which increases the suction and vortex lift. When a = 14°, an evident reattachment line and a separation line appear on the upper surface, and the leading edge. Evident stacking appears in the middle and tail section of the separation line, and based on the vortex counter, the leading edge vortex has broken up at this location where there is 13% body length ahead compared with that of a = 10°. Due to the breakdown of the leading edge vortex, the vortex lift of this region decreases dramatically which results in the lift coefficient decrease and positive pitching moment. When a = 16°, an evident reattachment line and a separation line appear on the upper surface and a separation line exists near the leading edge. The flow feature is similar to that of a = 14°, but the vortex broken location is ahead within small region. Compared with the a = 14° case, due to the development of leading edge vortex, the increase in vortex lift is more than the decrease caused by the vortex broken ahead, and thus the vortex lift and the lift coefficient of the whole vehicle go up with the increase in angle of attack.
CL
CD
Cm
Coarse Medium Fine
0.582 0.641 0.643
0.1392 0.1581 0.1578
0.1233 0.0928 0.0923
335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365
4.2. Vortex force component and role of vortex breakdown
366
The variation of the lift characteristic of the low aspect ratio flying wing configuration studied in this paper is directly related to the formation, development and breakdown of the leading edge vortex on the lee side. To view the variation of the vortex lift in detail, the leading edge suction analogy method proposed by Polhamus26 is applied to decomposing the lift into potential lift and vortex lift as follows:
367
CL ¼ CLp þ CLv
ð21Þ
where CLp indicates potential lift coefficient and CLv the vortex lift coefficient. Because the lift at zero-incidence is not zero due to the curve of the airfoil cross section of the studied flying wing configuration, the lift decomposition equation of Polhamus is revised in this paper that a constant CL0 which represents the zero-incidence lift coefficient is introduced into the equation to make it more applicable to the problem dis-
Force and moment coefficients for the three grids used in convergence study at a = 14°.
Grid
334
GCI CL
CD
Cm
0.496204 0.015272
0.664591 0.00929
1.34035 0.21876
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
368 369 370 371 372 373 374 376 377 378 379 380 381 382 383
CJA 729 12 November 2016
No. of Pages 12
6
Y. Tao et al.
Fig. 4
Fig. 5
Comparison of aerodynamic coefficients.
Instantaneous Q criterion isosurface (a = 14°, Q = 10).
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
CJA 729 12 November 2016
No. of Pages 12
Transonic wing stall of a blended flying wing common research model
Fig. 6 384 385 386 388 389 390 391 392 393 394 395 397 398 399 400
Comparison of calculated streamline on surface and oil visualization results.
cussed herein. Thus the final equation for potential lift coefficient is CLp ¼ Kp sin a cos a þ CL0 2
ð22Þ
where Kp is the lift-curve slope given by small-angle-of-attack potential-flow lifting surface theory, and accounts for the true boundary condition, and cos2 a arises from the assumption of a Kutta-type flow condition at the leading edge. While the equation for the vortex lift is used by the result of typical leading edge suction analogy method. CLv ¼ Kv cos a sin2 a
7
ð23Þ
where Kv sin2 a represents the potential-flow leading-edge suction, therefore the vortex normal force, and cos a indicates the component in the lift direction.
The flying wing model used here has a low aspect ratio equal to 1.54, and the leading edge sweepback angle is 65°, the rear sweepback angle is ±47°, =65°, aspect ratio 1.54. For those parameters, Kp = 1.9 and Kv = 2.7 can be obtained with Polhamus approach26 based on the Variations figure of Kp and Kv for delta wings provided by Ref.26. Because the Polhamus approach is suitable for incompressible flow mainly, a correction should be made for the compressible effect. According to Glauert rule, compressibility has no effect on the aerodynamic force of delta wing, when the Mach number is lower enough. C0L
CL ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 Ma21
401 402 403 404 405 406 407 408 409 410 411 412
ð24Þ
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
414
CJA 729 12 November 2016
No. of Pages 12
8 415 416 417
419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434
Y. Tao et al.
Take the assumption that CLp and CLv obey Glauert rule too. C0Lp
CLp ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; 1 Ma21
C0Lv
CLv ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 Ma21
ð25Þ
So Kp ¼ 4:3589 and Kv ¼ 6:1942 at Ma1 ¼ 0:9. Fig. 7 shows the comparison of lift coefficients predicted by Polhamus approach and Glauert rule and obtained from experiments. The predicted lift is much higher than the experiment results. So for transonic flow, this does not hold true in the case of such high number flow. Because the traditional Polhamus approach and Glauert rule have low precision to decompose vortex force component at transonic flow, the least square method (LSM) was employed to calibrate the Kp and Kv in this study. Lift coefficient curve shows that the slope is basically linear at angle of attack ranging from 2° to 4°. In this range, the least square method is used to fit the Eq. (22) and then Kp and CL0 were calibrated as at Ma1 ¼ 0:9. Then the potential lift coefficient curve CLp can be achieved as shown in Fig 7.
Fig. 7
Vortex lift coefficients for blending wing configuration.
The true vortex lift coefficient curve named as CL CLp curve can be obtained by Coefficient CL minus the potential lift coefficient CLp . In this range from 4° to 12°, the least square method is used to fit the Eq. (23) and then Kv was calibrated. In this case Kv equals to 3.149; the detail results are shown in Fig. 7. From the fitting results of the vortex lift coefficient curve, the fitting curve does not match the original curve very well but it still represents its main characteristics. Since the lift coefficients of CFD match well with the experimental result at low angle of attack, CLpDDES was replaced by CLpLSM . The vortex lift coefficients of DDES method were computed by CLDDES CLpLSM as shown in Fig. 7. The vortex lift coefficient calculated by the method above is shown in Fig. 7 where it can be found that the vortex lift remains zero when the angle of attack is less than 4° meaning the leading edge vortex has not formed yet; however, the vortex lift increases dramatically when the angle of attack is more than 4°. With the increase in angle of attack, the vortex lift coefficient suffers abrupt drop down when the angle of attack is more than 12°, which causes the abrupt decrease in the lift coefficient of flying wing configuration model. To explain the relation between the abrupt decrease and the vortex breakdown and the mechanism of vortex breakdown, the detail of the flow field is investigated. Fig. 8 shows the pressure distribution p/pref through vortex cores at a = 12, 14°, which indicates that shocks/vortex interactions do exist in the flow. These shocks/vortex interactions influence the flow structures obviously especially at larger angles of incidence. At a = 12°, two adverse pressure gradients at location of x=xcr ¼ 0:33 and x=xcr ¼ 0:77 affect the vortex strength to a different extent. At the first location, the adverse pressure gradient is low ((@p=@xaxial ¼ 0:43), which has a small influence on the vortex stability. However, at the second location where a strong normal shock exists, the pressure through vortex cores has a sudden rise and the adverse pressure gradient is high ((@p=@xaxial ¼ 4:17), which indicate that strong shocks/vortex interactions exist in these areas. At a = 14°, two adverse pressure gradients also exist in the flow field. Compared to that of a = 12°, the location of the second adverse pressure gradient moves forward to x=xcr ¼ 0:58, which is closer to the head of the configuration and the (@p=@xaxial increases to 7.26. Due to the high adverse pressure gradient, the main leading-edge vortex breakdown is homologous with the Q criterion isosurface (Fig. 9).
Fig. 8
Pressure distribution through vortex cores.
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479
CJA 729 12 November 2016
No. of Pages 12
Transonic wing stall of a blended flying wing common research model 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504
From the CFD results of flow structure, there are shocks/ vortex interactions at both lower incident and higher incidence, with a stronger interaction for higher incidence, and the interaction causes a considerable weakening of the vortex core, which may result in vortex breakdown. There seems to be a criterion for the vortex to feel the effect of the shock and remain coherent or breakdown. Ref.25 demonstrated an importance criterion parameter for vortex breakdown caused by shocks/vortex interaction, which is a combination of tangential velocity Uh and axial velocity Uaxial of the vortex core. The swirl ratio or the Rossby number is a nondimensional parameter defined as the ratio of axial velocity and tangential velocity of vortex core, as defined by Ro ¼ Uaxial =Uh . It is used as a criterion to measure the vortex intensity and the susceptibility of the vortex to shock-induced breakdown. The maximum swirl velocity of the vortex and the maximum axial velocity of the vortex cores are used for the calculation of Rossby number in this investigation. When the leading-edge vortex passes through a normal shock, the swirl velocity changes little, while the axial velocity will decrease distinctly. As a result, the Rossby number will decrease and the stability of the vortex decrease makes the vortex easier to breakdown. Spall and robinson found that if Ro > 1.4, the vortex is stable, and if 0.9 < Ro < 1.4, the vortex is unstable, while if Ro < 0.9, the vortex breaks down
Fig. 9
Instantaneous Q criterion isosurface (Q = 10).
9
when Rossby number was used as the criterion parameter for leading-edge vortex breakdown studies of delta wings. Using the criterion, the Rossby numbers at a = 12° and 14° were calculated which is shown in Fig. 10. From the result, it is clear that shock has a great influence on the vortex stabilities. For a = 12°, at the locations of x=xcr ¼ 0:33 and x=xcr ¼ 0:5, the Rossby numbers without exception decrease obviously. As the shocks and shocks/vortex interactions are weak at these locations, the vortex keeps stable even though the vortex strength decreases. When the vortex passes through the second adverse pressure gradient, shocks/vortex interactions become stronger and the vortex becomes unstable. Subsequently, the Rossby number decreases gradually to 0.9 at the location of x=xcr ¼ 0:93, which means the vortex is broken down here. At a = 14°, shocks/vortex interactions are stronger compared to those of a = 12°, Rossby number reduces rapidly at the location of x=xcr ¼ 0:57, and the vortex becomes unstable and then breaks down. Fig. 11 shows the shock structures in the flow field at a = 12° and 14°; due to the high curve slope at the leading edge, the flow accelerates to a local peak value and makes a weak shock produced at the location of x=xcr ¼ 0:25 at a = 12°, and a strong normal shock can be observed around the location of x=xcr ¼ 0:75. Near the central span-wise section of the configuration, two cross flow shocks appear as the interactions of the main leading-edge vortex and second vortex at this area. As the angle of incidence increases to 14°, the weak shock almost disappears at the location of x=xcr = 0.25 compared to a = 12°, while two new normal shocks introduce at the location around xxcr ¼ 0:50 0:55 and the vortex breaks down as it passes through the normal shocks and the cross flow shocks destructs too at this angle of attack. Fig. 12 shows pressure coefficient distribution at the symmetry plane on the wing for both angles of incidence. At a = 12°, around the location of x=xcr ¼ 0:75, a normal shock is there. As angles of incidence increase to 14°, this normal shock moves to x=xcr ¼ 0:55 which makes the vortex breaks down after the vortex passes through this shock. At the rear of the configuration, weak shocks can be observed for both angles of incidence due to the high curve slope. Through comparison, we found that this model had a less breakdown resistant to angles of incidence compared with
Fig. 10
Rossby number distribution against root chord.
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546
CJA 729 12 November 2016
No. of Pages 12
10
Y. Tao et al.
Fig. 12 Pressure coefficient distribution at symmetry plane on the wing for both angles of incidence.
Fig. 11
547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569
Shock structures in flow field.
the other 65° swept delta wing configurations proposed in Ref.3,6 In transonic flow field, the wing-body thickness ratio has great influence on the shock behavior on delta wings and the occurrence of shocks introduces complex shock/vortex interactions which are responsible for vortex breakdown. In addition, almost all the published investigation models with flying-wing configurations have a relatively small wing-body thickness ratio, and the flow differences around flying-wing configuration with a moderate and a small wing-body thickness ratio are not entirely described and understood. Based on the previous study of the 65° swept flying-wing configuration with a moderate 0.16 wing-body thickness ratio (marked as Body-orig), through reducing the wing-body thickness (marked as Body-thin) and keeping the round leading edge and outer wing geometry identical, we investigated wingbody thickness effects on the aerodynamic and vortex flow characteristics of the flying-wing configuration at Ma = 0.9. Fig. 13 shows comparisons of typical cross-sectional shapes between the original wing-body thickness ratio model and the reduced wing-body thickness ratio model. Fig. 14 shows comparisons of lift and pitching moment coefficients versus incidence relative to the ‘original’ and ‘reduced’ models. For a < 12°, effect of inner wing thickness
Fig. 13
Typical cross-sectional shapes of two configurations.
on lift coefficient is very weak. However, as angle of incidence increases, the lift coefficient of ‘original’ model becomes relatively flat up to 15° incidence, whereas the lift coefficient of ‘thin’ model has a linear increase up to 20°. A similar phenomenon occurs in the pitching moment. Abrupt nose up in pitching moment of ‘reduced’ model is at 20°, while the ‘original’ one at 12°. That is to say ‘reduced’ model can delay stall and vortex breakdown by an angle of 8° compared with the original model at Mach number 0.9. From the pressure distribution through vortex cores at a = 14° in Fig. 15, it is clear that accelerating pressure exists up to x=xcr ¼ 0:65 of the ‘reduced’ model, and the location of abrupt pressure rise is far from the head of the configuration compared to the ‘original’ model, which explains why the ‘reduced’ model has a larger stall angle of incidence. From this investigation we can see that wing-body thickness has great influence on the introduction of shock and shock/vortex interactions which can make a significant difference to the vortex breakdown behavior and aerodynamic characteristic of the flying-wing configuration, which should be considered in the aerodynamic design process.
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590
CJA 729 12 November 2016
No. of Pages 12
Transonic wing stall of a blended flying wing common research model
11
5. Conclusions
591
Numerical simulation of abrupt stall phenomena of blended flying wing configuration at transonic speed was conducted using both DDES and URANS methods at transonic speed. Based on the analysis of the aerodynamic coefficients, flow structures and streamlines, the following conclusions can be drawn:
592
(1) The sudden motion in breakdown location observed in mean lift, drag and pitching moment coefficients of experiments is due to a shock/vortex interaction. (2) The DDES methods perform more precisely than the URANS method (the onset angle of the breakdown movement was predicted about 2° greater than the measurements) and the aerodynamic coefficient results from DDES method were compared very well with the experiment results and the shock strength or axial flow in the vortex was well predicted by the DDES methods. (3) Wing-body thickness has great influence on shock and shock/vortex interactions, which can make a significant difference to the vortex breakdown behavior and stall characteristic of the blended flying wing configuration.
598
593 594 595 596 597
599 600 601 602 603 604 605 606 607 608 609 610 611 612
Fig. 14 Comparisons of lift and pitching moment coefficients vs incidence.
Fig. 15 cores.
Comparison of pressure distribution through vortex
Acknowledgment
613
This work was funded by the National Natural Science Foundation of China (No. 11372337).
614
References
616
1. Schiavetta LA, Boelens OJ, Fritz W. Analysis of transonic flow on a slender delta wing using CFD. Reston: AIAA; 2006, Report No.: AIAA-2006-3171. 2. Elsennnr A, Hoeijmakers HWM. An experimental study of the flow over a sharp-edged delta wing at subsonic and transonic speeds. Paris: AGARD; 1991, Report No.: AGARD-CP-494. 3. Chu J, Luckring JM. Experimental surface pressure data obtained on a 65° delta wing across Reynolds number and Mach number ranges: Volume 1-sharp leading edge.. Washington (DC): NASA; 1996, Report No.: NASA TM 4645. 4. Hill DJ, Pantano C, Pullin DI. Large-eddy simulation and multiscale modeling of a Richtmyer-Meshkov instability with reshock. J Fluid Mech 2006;557:29–61. 5. Fang BR. Aero-configuration design of aero plan. Beijing: Aviation Industry Press; 1997 [in Chinese]. 6. Jobe CE. Vortex breakdown location over 65° delta wings empiricism and experiment. Aeronaut J 2004;108:475–82. 7. Kalkhoran IM, Smart MK. Aspects of shock wave-induced vortex breakdown. Progress Aerospace Sci 2000;36(1):63–95. 8. Fritz W, Cummings RM. What was learned from the numerical simulations for the VFE-2. Reston: AIAA; 2008, Report No.: AIAA-2008-0399. 9. Schiavetta L, Boelens O, Crippa S. Sting effects on vortex breakdown in transonic delta wing flows. Reston: AIAA; 2008, Report No.: AIAA-2008-0395. 10. Schiavetta LA, Boelens OJ, Crippa S. Shock effects on delta wing vortex breakdown. J Aircraft 2009;46(3):903–14. 11. Spalart PR, Jou WH, Strelets M, Allmaras SR. Comments on the feasibility of les for wings, and on a hybrid rans/les approach, advances in DNS/LESProceedings of the 1st AFOSR Int. Conf. on DNS/LES. Columbus (OH): Greyden Press; 1997.
617
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
615
618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647
CJA 729 12 November 2016
12 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672
12. Viswanathan A, Squires KD, Forsythe JR. Detached-eddy simulation around a forebody at high angle of attack. Reston: AIAA; 2003, Report No.: AIAA-2003-0263. 13. Lv ZY, Zhu LG. Study on forms of vortex breakdown over delta wing. Chin J Aeronaut 2004;17(1):13–6. 14. Forsythe JR, Squires KD, Wurtzler KE, Spalart PR. Detachededdy simulation of the F-15E at high alpha. Reston: AIAA; 2002, Report No.: AIAA-2002-0591. 15. Morton SA, Steenman MB, Cummings RM, Forsythe JR. DES grid resolution issues for vertical flows on a delta wing and an F18C. Reston: AIAA; 2003, Report No.: AIAA-2003-1103. 16. Denning RM, Allen JE, Armstrong FW. The broad delta airliner. Aeronaut J 2003;107(1075):547–58. 17. Su JC, Huang Y, Li YH. Research on flow characteristics of lowaspect-ratio flying wing at transonic speed. Acta Aerodynamic Sinica 2015;33(3):307–12 [in Chinese]. 18. Simon F, Deck S, Guillen P, Sagaut P, Merlen A. Numerical simulation of the compressible mixing layer past an axisymmetric trailing edge. J Fluid Mech 2007;591:215–53. 19. Lu XY, Wang SW, Sung HG. Large-eddy simulations of turbulent swirling flows injected into a bump chamber. J Fluid Mech 2005;527:171–95. 20. Wang SW, Yang V, Hsiao G, Hsieh SY. Large-eddy simulations of gas-turbine swirl injector flow dynamics. J Fluid Mech 2007;583:99–122.
No. of Pages 12
Y. Tao et al. 21. Hill DJ, Pantano C, Pullin DI. Large-eddy simulation and multi scale modeling of a Richtmyer-Meshkov instability with reshock. J Fluid Mech 2006;557:29–61. 22. Menter FR, Kuntz M, Langtry R. Ten years of experience with the SST turbulence model. In: Proceedings of the 4th international symposium on turbulence heat and mass transfer; 2003. p. 625_32. 23. Menter FR. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J 1994;32(8):1598–605. 24. Jian Liu, Haisheng Sun, Zhitao Liu, Zhixiang Xiao. Numerical investigation of unsteady vortex breakdown past 80°/65°doubledelta wing. Chin J Aeronaut 2014;27(3):521–30. 25. Spalart PR. Detached-eddy simulation. Annu Rev Fluid Mech 2009;41:181–202. 26. Polhamus EC. Application of the leading edge suction analogy of vertex lift to the drag due to lift of sharp edge delta wings. Washington (DC): NASA; 1968, Report No.: TN-D-4739. Tao Yang is an associate research fellow in China Aerodynamics Research and Development Center. His main research interests are vortex aerodynamics, unsteady flows and their control. Zhao Zhongliang is a research fellow in China Aerodynamics Research and Development Center. His main research interests are unsteady flows, virtual flight testing and flight control.
Please cite this article in press as: Tao Y et al. Transonic wing stall of a blended flying wing common research model based on DDES method, Chin J Aeronaut (2016), http://dx.doi.org/10.1016/j.cja.2016.10.018
673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697