CFD analysis of the flow around the X-31 aircraft at high angle of attack

Aerospace Science and Technology 20 (2012) 38–51

Contents lists available at SciVerse ScienceDirect

Aerospace Science and Technology www.elsevier.com/locate/aescte

CFD analysis of the flow around the X-31 aircraft at high angle of attack O.J. Boelens 1 Department of Flight Physics and Loads, Aerospace Vehicles Division, National Aerospace Laboratory NLR, P.O. Box 90502, 1006 BM, Amsterdam, The Netherlands

a r t i c l e

i n f o

Article history: Available online 29 March 2012 Keywords: CFD X-31 Vortical flow

a b s t r a c t Accurate and cost-effective Computational Fluid Dynamics (CFD) methods play an increasingly important role, for example in the support of fighter aircraft operations. Prior to deployment of such CFD methods they should be well validated and evaluated against state-of-the-art wind tunnel and/or flight test data. The wind tunnel experiments on the detailed X-31 model performed in the DNW Low-Speed-WindTunnel Braunschweig (DNW-NWB) provide an excellent data set for validation and evaluation purposes. This data set has been investigated in the framework of NATO RTO task group AVT-161 ‘Assessment of Stability and Control Prediction Methods for NATO Air & Sea Vehicles’. The National Aerospace Laboratory NLR participated in this task group using its in-house developed flow simulation system ENFLOW, which includes both grid generation tools and a flow solver. The focus of the present paper is on the question to what extend leading edge details, flap gaps, need to be taken into account for the X-31 wind tunnel model to properly simulate the flow around this configuration. To investigate this question, three leading edge configurations have been considered, i.e. one with all leading edge flap gaps, one with only the longitudinal flap gaps and one with no leading edge flap gaps. Results obtained for selected test conditions measured during test run VN01004 (M = 0.18 and Rem.a.c . = 2.07 × 106 ) of the wind tunnel experiments will be discussed. Properly modeling geometrical details of the wind tunnel model at the leading edge is essential in obtaining the vortical flow phenomena observed in the wind tunnel. Analysis of the pitching moment coefficient demonstrates how in case of not resolving geometrical details a seemingly correct behavior is obtained without, however, resolving the underlying flow physics correctly. © 2012 Elsevier Masson SAS. All rights reserved.

1. Introduction Cost-effective Computational Fluid Dynamics (CFD) methods with sufficient accuracy are complementary to experimental methods and play an increasingly important role in simulating manoeuvre conditions. Examples are conditions that cannot be simulated in a wind tunnel or are too dangerous to be performed in flight tests. Stability and control characteristics can be obtained for such conditions using CFD. Prior to the usage of the computed stability and control characteristics the CFD methods should be well validated and evaluated against state-of-the-art wind tunnel and/or flight test data. The proper selection of the level of geometrical detail to be used in a CFD simulation has a large impact on the costeffectiveness. Including more geometrical detail gives rise to a more complex and thus more expensive grid generation task. Incorporating more geometrical detail will also result in grids with more cells. Since the computing time and hence the computational cost depend directly on the number of grid cells, incorporating more geometrical detail will result in more costly simulations.

1

E-mail address: [email protected]. R&D Engineer, Applied Computational Fluid Dynamics.

1270-9638/$ – see front matter © 2012 Elsevier Masson SAS. All rights reserved. http://dx.doi.org/10.1016/j.ast.2012.03.003

Omitting these geometrical details may however result in inaccurate or even wrong results. It is well known that geometrical details of the wing leading edge geometry, i.e. the leading edge snag, have triggered a bifurcation in the flow in the first version of the updated F/A-18, Refs. [6] and [5]. Failure of predicting this bifurcation phenomenon during the design stage resulted in a severe penalty in performance, and gave rise to substantial redesign costs. Another aspect impacting cost-effectiveness is the proper selection of the physical modeling used. Employing a too simple physical modeling may result in not resolving all physics and hence inaccurate or wrong results. Using a too complex physical model will lead to excessive computational times and thus excessive computational costs. The present paper highlights these aspects for the complex flow over the challenging X-31 wind tunnel model geometry. The wind tunnel experiments on the detailed X-31 model, see Fig. 1, performed in the framework of the DLR project SikMa“Simulation of Complex Maneuvers” [7,14,13] provide an excellent data set for validation and evaluation purposes. During the course of this project several wind tunnel entries have been executed in the DNW Low-Speed-Wind-Tunnel Braunschweig (DNW-NWB, 3.25 m × 2.80 m). These entries were specifically aimed at generating high quality data which could be used in a later stage for

O.J. Boelens / Aerospace Science and Technology 20 (2012) 38–51

39

Nomenclature AVT CFD Cp DLR DNW EARSM FAS FMG k L ref M m.a.c . NATO NLR p

Applied Vehicle Technology (one of the seven panels within RTO) Computational Fluid Dynamics pressure coefficient (= ( p − p ∞ )/(1/2ρ∞ u 2∞ )) Deutsches Zentrum für Luft- und Raumfahrt, German Aerospace Center German–Dutch Wind tunnels Explicit Algebraic Reynolds Stress Model Full Approximation Storage Full Multi-Grid turbulent kinetic energy reference length Mach number mean aerodynamic chord North Atlantic Treaty Organization Nationaal Lucht- en Ruimtevaartlaboratorium, Netherlands National Aerospace Laboratory pressure

p∞ RANS Re RTO TNT u u∞ X-LES x y y+

α β

ρ ρ∞ ω

free-stream pressure Reynolds-Averaged Navier–Stokes Reynolds number Research and Technology Organization – scientific arm of NATO Turbulent Non-Turbulent velocity free-stream velocity Extra-Large Eddy Simulation distance along the model body axis, positive aft distance along the span, positive outward Re-like term for flat plate turbulent boundary layer angle of attack, ◦ side-slip angle, ◦ density free-stream density specific turbulent dissipation rate

Fig. 2. Pressure sensor locations: x = 955 mm (60% chord) and x = 1125 mm (70% chord).

Fig. 1. X-31 model in rear sting configuration in DNW Low-Speed-Wind-Tunnel Braunschweig (DNW-NWB) (photo by courtesy of Deutsches Zentrum für Luft- und Raumfahrt, DLR).

the validation of CFD methods. Both steady-state measurements as well as simulations of complex manoeuvres employing a test rig with 6 degrees of freedom have been performed. This data set has been provided by Deutsches Zentrum für Luftund Raumfahrt DLR to the partners participating in NATO RTO task group AVT-161 ‘Assessment of Stability and Control Prediction Methods for NATO Air & Sea Vehicles’ [8,12]. The objectives of this task group are defined as follows: (i) to assess the state-of-the-art in computational fluid dynamics methods for the prediction of static and dynamic stability and control characteristics of military vehicles in the air and sea domains, and (ii) to identify shortcomings of current methods and identify areas requiring further development. The National Aerospace Laboratory NLR participates in this task group using its in-house developed flow simulation system ENFLOW [2], which includes both grid generation tools and a flow solver. This article concentrates on the question to what extend leading edge details, flap gaps, need to be taken into account for the X-31 wind tunnel model to properly simulate the flow around this

configuration. The layout of the article is as follows. The wind tunnel experiments will be discussed in Section 2. The focus will be on a test run which shows a sudden drop in the pitching moment coefficient between angles of attack α of 17◦ and 23◦ , whereas both the lift and drag coefficient increase monotonously for those angles. The grid generation will be discussed in Section 3. Section 4 will discuss the important features of the flow solver ENSOLV, which is part of the flow simulation system ENFLOW. Section 5 will discuss the results obtained at NLR for the above mentioned test run. A section with conclusions (Section 6) completes the paper. 2. Wind tunnel experiments The model used during the wind tunnel experiments [7,14, 13] is based on the X-31 experimental high angle-of-attack aircraft configuration, see Fig. 1. The model is equipped with control devices driven by remotely controlled internal servo-engines to deflect the canard, the two leading edge flaps and the trailing edge flap. Sensors are installed to measure the aerodynamic forces and moments on the model, as well as sectional surface pressure distribution at x = 955 mm (60% chord length) and x = 1125 mm (70% chord), see Fig. 2. The experiments also included steady-state measurements using Pressure Sensitive Paint (PSP), which provides detailed information on the surface pressure distribution on the whole wing. Measurements were performed for a large range of angles of attack and side-slip angles. During the measurements also the effect of different canard settings as well as leading and trailing edge flap settings has been investigated. Test run VN01004 [7] has been selected from the wind tunnel data set as reference for the sectional surface pressure distributions and the aerodynamic force and moment data. This test run constitutes an α -sweep for angles of attack α ranging from −6◦ to 55◦ .

40

O.J. Boelens / Aerospace Science and Technology 20 (2012) 38–51

Fig. 3. Lift, drag and pitching moment coefficient for test run VN01004 (M = 0.18 and Rem.a.c . = 2.07 × 106 ).

For angles of attack α smaller than 16◦ an interval of 2◦ has been used, while an interval of 1◦ has been used for the remainder of the α -sweep. The wind tunnel velocity was 60 m/s, which corresponds to a Mach number of 0.18. The Reynolds number based on

the mean aerodynamic chord was 2.07 million. During this test run all control surfaces were in their respective neutral position. In addition, test run VN01030 [7] has been used. This test run provides Pressure Sensitive Paint data obtained at the same conditions. Fig. 3 shows the lift, drag and pitching moment coefficient for test run VN01004. The present investigation focuses on the behavior of the pitching moment coefficient for angles of attack α between 17◦ and 23◦ and whether this behavior can be simulated with Computational Fluid Dynamics. Between these angles the pitching moment coefficient shows a clear drop, whereas both the lift coefficient and the drag coefficient still increase monotonously. The sectional surface pressure coefficient measured during the same test run is shown in Fig. 4 for angles of attack α ranging from approximately 10◦ to 22◦ at a 2◦ interval. In addition, Fig. 5 shows the corresponding surface pressure coefficient obtained using Pressure Sensitive Paint during test run VN01030. Both figures show that there are some distinct suction peaks on the wing upper surface, especially for angles of attack α larger than 14◦ . Going from wing root to wing tip these peaks will now be discussed. For more details see also Section 5.2 and Ref. [14]. The first peak (around y = 120 mm at x = 955 mm and y = 140 mm at x = 1125 mm in Fig. 4) is due to a vortex originating at the strake. The strength of this vortex increases monotonously for the range of angles of attack shown. Next, a suction peak is observed

Fig. 4. Sectional surface pressure coefficient for test run VN01004 (M = 0.18 and Rem.a.c . = 2.07 × 106 ).

Fig. 5. Surface pressure coefficient obtained using Pressure Sensitive Paint for test run VN01030 (M = 0.18 and Rem.a.c . = 2.07 × 106 ).

O.J. Boelens / Aerospace Science and Technology 20 (2012) 38–51

41

Fig. 6. Changes in the contribution to the pitching moment coefficient C p |n|(x − xm.r . p . )/( L ref A ref ) for the upper wing obtained from the Pressure Sensitive Paint data for test run VN01030 (M = 0.18 and Rem.a.c . = 2.07 × 106 ). Blue indicates a positive contribution, whereas red indicates a negative contribution. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

in Fig. 4 around y = 180 mm at x = 955 mm and, though weak, at y = 200 mm at x = 1125 mm. This vortex is formed at the inner wing leading edge. Its strength grows rapidly going from an angle of attack α of 16◦ to 22◦ . The third peak (around y = 230 mm at x = 955 mm and y = 280 mm at x = 1125 mm in Fig. 4) is caused by a vortex which is formed at the longitudinal flap gap between the wing and the inner flap (gap1 in Fig. 7). The strength of this vortex increases monotonously with the angle of attack α . Finally, a suction peak is observed at x = 1125 around y = 380 mm in Fig. 4. This peak is due to a vortex having its origin at the longitudinal flap gap between the inner flap and the outer flap (gap2 in Fig. 7). This vortex appears at an angle of attack α between 12◦ and 14◦ . The associated pressure peak remains fairly constant at a relatively high level until an angle of attack α of 18◦ . For angles of attack α higher than 18◦ the pressure peak decreases monotonously. At an angle of attack α of 22◦ the suction peak is not present anymore. To further investigate what causes the drop in pitching moment coefficient, Fig. 6 shows the change in the contribution to the pitching moment coefficient C p |n|(x − xm.r . p . )/( L ref A ref ) for the wing upper surface, where C p equals C p (α2 ) − C p (α1 ), xm.r . p . is the location of the moment reference point, |n| is the area of the faces used in the discretization of the Pressure Sensitive Paint data, and L ref and A ref are the reference length and reference area, respectively. The moment reference point is located at xm.r . p . = 937.66 mm. C p |n|(x − xm.r . p . )/( L ref A ref ) gives an indication of the regions on the wing upper surface contributing to changes in the pitching moment coefficient. In Fig. 6, blue indicates a positive contribution, i.e. the pitching moment coefficient increases, whereas red indicates a negative contribution, i.e. the pitching moment coefficient decreases. Based on this figure and the discussion in the previous paragraph it may be concluded that the drop in the pitching moment coefficient is mainly due to the behavior of the vortex originating at gap2 and the vortex formed at the inner wing leading edge. Large negative contributions to the pitching moment can be observed in Fig. 6 for these vortices for angles of attack α between 16◦ and 22◦ . Note also that the strake vortex seems to have a net

positive contribution to the change in pitching moment going from the angles of attack α of 16◦ to 18◦ and a net negative contribution going from the angles of attack α of 18◦ to 20◦ . 3. Grid generation To investigate to what extend leading edge details, flap gaps, need to be taken into account for the X-31 wind tunnel configuration, three different leading edge configurations are considered, i.e.

• the X-31 wind tunnel configuration with all leading edge flap gaps,

• the X-31 wind tunnel configuration with only the longitudinal leading edge flap gaps, and

• the X-31 wind tunnel configuration with no leading edge flap gaps. These three configurations are shown in Fig. 7. It should be noted that during the selected wind tunnel experiment, and therefore also during the simulations, all control surfaces were in their respective neutral position. Structured (multi-block) grids have been generated using NLR’s Cartesian grid mapping technique [1]. The (semi-automatic) grid generation algorithms have been developed at NLR and are part of NLR’s ENFLOW flow simulation system [2]. The method starts with an abstraction of the geometry description, see Fig. 8(a)). In such an abstraction, the geometry including all details is represented by a set of Cartesian blocks. This abstraction is then projected on the real geometry description; see Fig. 8(b)). Next the Navier– Stokes blocks, i.e. the first layer of blocks around the geometry, and the field blocks in physical space are generated automatically. In the resulting topology, the grid dimensions are set manually and the edges are connected automatically. Finally, the grid quality is improved by applying an elliptic smoothing algorithm and the resolution in the Navier–Stokes blocks is improved such that the desired boundary layer resolution is obtained. For the present grids the spacing of the first grid point normal to the solid wall

42

O.J. Boelens / Aerospace Science and Technology 20 (2012) 38–51

Fig. 7. X-31 leading edge configurations considered; X-31 configuration with all leading edge flap gaps (top left figure), X-31 configuration with only the longitudinal leading edge flap gaps (top right figure) and X-31 configuration with no leading edge flap gaps (bottom figure).

Fig. 8. Abstraction of the surface geometry (a) and projected abstraction (b) for the X-31 configuration (G3 (no leading edge flap gaps)).

is 0.005 mm, resulting in a typical y + -value of approximately one. Away from the wall, the spacing increases by a ratio s2 /s1 of 1.1 to ensure smoothness. More details on this grid generation approach can be found in Ref. [1]. First a grid has been generated around the X-31 configuration with no leading edge flap gaps (G3). The sting present in the wind tunnel has not been modeled. In this grid compound faces have been used for the wing upper and lower surface. These compound faces can consist of any number of subfaces. Thus, it is possible to include or reserve the location of the flap gaps and the grid resolution required for the grids with flap gaps in the grid around the X-31 configuration with no leading edge flap gaps. Therefore, all three configurations can be handled with the same topology, except for the blocks in the flap gaps. Next based on this grid, the grids around the X-31 configuration with only the longitudinal leading edge flap gaps (G2) and the X-31 configuration with all leading edge flap gaps (G1) have been constructed by inserting the

respective blocks. Generating the grids in this manner has as major advantage that the grids are virtually identical, except for the flap gaps. Small differences due to the elliptic smoothing are confined to the region very close to the flap gaps. Using this approach, grids around complex configurations can be generated within very short times. For the X-31 configuration with no leading edge flap gaps the grid was completed in slightly more than one week, whereas including first the longitudinal flap gaps and next the span wise flap gaps only took one more day each. The resulting surface grid and grid on the symmetry plane for the X-31 configuration with all leading edge flap gaps (G1) are shown in Fig. 9. In addition, for all three leading edge configurations a detail of the surface grid around gap1, i.e. the longitudinal flap gap between the wing and the inner flap, is shown in Fig. 10. Details on the grids can be found in Table 1. Finally, the blocks in all three grids have been merged. The resulting grids therefore

O.J. Boelens / Aerospace Science and Technology 20 (2012) 38–51

43

Fig. 9. Surface grid and the grid on the symmetry plane for the X-31 configuration with all leading edge flap gaps (G1).

Fig. 10. Detail of the surface grid around gap1; X-31 configuration with all leading edge flap gaps (G1) (top left figure), X-31 configuration with only the longitudinal leading edge flap gaps (G2) (top right figure) and X-31 configuration with no leading edge flap gaps (G3) (bottom figure). The surface grid in the flap gaps is colored red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

consist of blocks with larger grid dimensions, which will result in larger vector lengths and thus shorter computing times when computing on NLR’s SX-8 vector computer. Details on the merged grid can also be found in Table 1. 4. Flow solver 4.1. General description The flow solver ENSOLV, which is part of NLR’s flow simulation system ENFLOW [2], is capable of solving the Euler and Navier–Stokes equations on multi-block structured grids for arbitrary configurations. The configuration can be either fixed or mov-

ing relative to an inertial reference frame, and can be either rigid or flexible. The flow equations are cast into a full conservation form employing the density ρ , the components of the momentum vector ρ u and the total energy per unit volume ρ E as dependent variables. The equations are non-dimensionalized using the freestream static pressure, the free-stream density, the free-stream temperature and a reference length (for example the mean aerodynamic chord). The equations in full conservation form are discretized in space by a second-order accurate, cell-centred, finite-volume method, using multi-block structured grids, central differences, and matrix artificial diffusion. The artificial diffusion consists of a blending

44

O.J. Boelens / Aerospace Science and Technology 20 (2012) 38–51

Table 1 Details of the grids used.

Number of blocks Number of grid cells Number of grid points Number of blocks after merging

Table 2 Computational details of the simulations on G1 (all leading edge flap gaps) grid. G1 (all leading edge flap gaps)

G2 (longitudinal leading edge flap gaps)

G3 (no leading edge flap gaps)

1317 24,899,072

1307 24,737,792

1299 24,651,776

27,826,365

27,641,739

27,542,211

207

197

189

of second-order and fourth-order differences with a Jameson-type shock sensor for the basic flow equations and a TVD discontinuity sensor for the turbulence model equations. For steady flow simulations, the discretized time-dependent system of equations is integrated toward the steady-state solution using a five-stage explicit Runge–Kutta scheme. Local-time stepping and multi-grid acceleration techniques are applied. For timeaccurate simulations, the flow solver uses the dual-time stepping scheme, where for each time step the time-dependent flow equations are integrated in pseudo-time toward a steady-state solution in a similar way as in the steady flow simulation using the same acceleration techniques. 4.2. Turbulence model Several turbulence models are present in the flow solver ENSOLV, including the Turbulent Non-Turbulent (TNT) k–ω model [9,3,4], the EARSM model [4] and a hybrid RANS-LES model for eXtra-Large Eddy Simulation (X-LES) [11,10]. The TNT k–ω model, which has been used for most simulations in the present study, is a variant of the Wilcox k–ω model. The equations of the model are slightly modified by the introduction of a ‘cross diffusion’ term [9]. This modification has been introduced to resolve the dependency of the free-stream value of ω . The model is also extended with a global correction for vortical flows [3,4]. It is well known that the standard model, as with most other two-equation models, over predicts the eddy viscosity within the vortex core which leads to exaggerated diffusion of vorticity. As a consequence the details of the vortex core are lost and low suction peaks with wide vortex bases are a characteristic of the solution. The enhanced model [3] controls the production of turbulent kinetic energy and hence eddy viscosity through an increase in the production of dissipation (ω) within regions of highly rotational flow. A suitable sensor has been used to distinguish between shear layers and vortex cores. This sensor is the ratio between the magnitude of strain-rate and vorticity tensor. In shear layers, the velocity gradient is dominated by the gradient in the normal direction, which results in a ratio of approximately one, while in vortex cores, where the flow experiences pure rotation, the ratio is much less than one. This approach has proven to be effective in producing surface pressure profiles on simple delta wings in good agreement with those of experimental data [3,4]. In addition, to remove the singular behavior of ω at solid boundaries, the equations of the k–ω model are reformulated such that instead of ω the quantity τ = 1/(ω + ω0 ) is used. Here ω0 is a positive constant (default value ω0 L ref /u ∞ = 20, with u ∞ the freestream velocity and L ref the reference length). Finally, the source terms in the k–ω equations are treated explicitly, while a separate time-step is used for the k–ω equations to enhance the efficiency of the scheme. At the solid boundaries, both k and τ are set to zero. To prevent unphysical high values of k near stagnation points, the production term in the k-equation has been limited to a maximum of 5 times the dissipation term in the k-equation.

4 h-grid level 2 h-grid level h-grid level

Number of grid cells

Number of iterations

Computational wall-clock time (use 1 processor of NLR’s SX8)

389,048 3,112,384 24,899,072

2250 1750 1250

1 h 10 m 5 h 38 m 21 h 41 m

Total

29 h 29 m

At the ‘inflow’ parts of the far-field boundary, the free-stream values of the turbulent variables are computed from the freestream turbulent Reynolds number (0.01 in the present simulations) and the free-stream dimensionless turbulent kinetic energy (k/u 2ω = 10−6 in the present simulations). All simulations were run in fully turbulent mode. 4.3. Boundary conditions At the X-31 geometry, a no-slip viscous flow condition (adiabatic solid wall) has been employed. For the upstream, downstream, top, bottom and side far-field faces, a free-stream boundary condition based on Riemann invariants of the locally linearized one-dimensional Euler equations has been used. Since the flow at these faces is subsonic, the value of the ‘incoming’ Riemann invariants is computed using the free-stream values. The remaining invariants are extrapolated from the computational domain. Finally, a symmetry boundary condition has been used at the symmetry plane. For this boundary condition the grid does not necessarily need to be orthogonal to the symmetry plane. 4.4. Details of simulations All simulations have been performed as steady flow simulations. A Full Multi-Grid (FMG) scheme (grid sequencing) has been used to compute the solution on the three grid levels. The solution on a coarse level has been used as initial solution on the next-finer level. The number of iterations on each grid level is shown in Table 2. The Full Approximation Storage (FAS) multi-grid scheme has been used to compute the solution on a specific grid level. One FAS multi-grid level has been used. The simulations have been performed on two processors of NLR’s NEC SX5/8B vector computer. Three to four orders of convergence have been obtained for the root mean square norms. Computational details of the simulations are shown in Table 2. 5. Results 5.1. Test cases Simulations have been performed on the three grids described above. The test conditions used are summarized in Table 3, see also Ref. [7] and Section 2. These test conditions have been selected from test run VN01004. In addition the test conditions at which the corresponding Pressure Sensitive Paint data was obtained are summarized in this table. The PSP test conditions belong to test run VN01030. Note that these data were obtained at slightly different angles of attack α . As stated above, for both test runs the wind tunnel velocity is 60 m/s, which corresponds to a Mach number of 0.18. The Reynolds number based on the mean aerodynamic chord equals 2.07 million. The TNT k–ω model extended with the global correction for vortical flows has been used as turbulence model in all simulations. For test point MP0009 (α = 10.05◦ ) and MP0016 (α = 20.06◦ ) additional simulations have been performed on grid G1

O.J. Boelens / Aerospace Science and Technology 20 (2012) 38–51

45

Table 3 Test conditions. VN01004

VN01030

Force, moment and pressure sensor data ◦

Test point

α,

MP0003 MP0006 MP0009 MP0010 MP0011 MP0012 MP0014 MP0016 MP0018

−1.98 4.03 10.05 12.05 14.05 16.05 18.05 20.06 22.06

PSP data for complete X-31 configuration

β,

◦

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

M

Rem.a.c .

0.18 0.18 0.18 0.18 0.18 0.18 0.18 0.18 0.18

2.07 × 10 2.07 × 106 2.07 × 106 2.07 × 106 2.07 × 106 2.07 × 106 2.07 × 106 2.07 × 106 2.07 × 106 6

◦

Test point

α,

4 10 16 18 20 22 24 26 28

−2.02 3.98 9.98 11.98 13.98 15.98 17.98 19.98 21.98

β,

◦

0.01 0.01 0.01 0.01 0.01 −0.01 −0.01 0.00 0.00

Fig. 11. Contour plot of the total pressure loss for test point VN01004-MP0009 (α = 10.05◦ , M = 0.18 and Rem.a.c . = 2.07 × 106 ) simulated with TNT k–ω with vortex correction on G1 (all leading edge flap gaps). The geometry is colored by the pressure coefficient C p .

(all leading edge flap gaps) employing the standard TNT k–ω model and the Explicit Algebraic Reynolds Stress Model (EARSM). The results of these simulations will be discussed in the following sections. 5.2. Flow structure Before looking to the results in detail, the large-scale vortical flow structure above the wing is discussed. Figs. 11 and 12 show a contour plot of the total pressure loss for test point MP0009 (α = 10.05◦ ) and MP0016 (α = 20.06◦ ), respectively. These results have been obtained on grid G1 (all leading edge flap gaps) using the TNT k–ω model extended with the global correction for vortical flows. For test point MP0009 (Fig. 11) which has a moderate angle of attack α of 10.05◦ , the following vortices are observed: (i) (ii) (iii) (iv) (v)

a a a a a

strake vortex, canard vortex, tip vortex, fuselage vortex, and vortex originating at the inlet region.

The flow on most part of the wing is attached. Only a (weak) signature of the strake vortex can be found in the wing upper surface pressure distribution. For test points having an angle of attack smaller than approximately 12◦ the vortical flow structure is similar to the one discussed here. For test point MP0016 which has a high angle of attack α of 20.06◦ the vortical flow structure is much more complex, see also Ref. [14]. The following vortices can be observed: (i) (ii) (iii) (iv) (v) (vi) (vii)

a strake vortex, an inner wing vortex originating from the wing leading edge, a vortex originating at gap1, an outer wing vortex originating at gap2, a canard vortex, a fuselage vortex, and a vortex originating at the inlet region.

The wing upper surface pressure distribution is highly influenced by the presence of these vortices, showing regions of strong suction. Furthermore it can be seen that the strake vortex and the inner wing vortex originating from the wing leading edge interact

46

O.J. Boelens / Aerospace Science and Technology 20 (2012) 38–51

Fig. 12. Contour plot of the total pressure loss for test point VN01004-MP0016 (α = 20.06◦ , M = 0.18 and Rem.a.c . = 2.07 × 106 ) simulated with TNT k–ω with vortex correction on G1 (all leading edge flap gaps). The geometry is colored by the pressure coefficient C p .

Fig. 13. Sectional surface pressure coefficient simulated on G1 (all leading edge flap gaps) for test point VN01004-MP0009 (α = 10.05◦ , M = 0.18 and Rem.a.c . = 2.07 × 106 ) at x = 955 mm (left figure) and x = 1125 mm (right figure). TNT k–ω : blue line, TNT k–ω with vortex correction: green line, EARSM: red line. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

resulting downstream in a single vortex. The same holds for the vortex originating at gap1 and the outer wing vortex originating at gap2. It should be noted that for vortical flow structures like the one shown, i.e. a cascade of vortices, the mutual influence of vortices is large. For example, changes in the strake vortex will influence the inner wing vortex originating from the wing leading edge, which will then influence the vortex originating at gap1, and so on. Test points having an angle of attack larger than approximately 16◦ show similar vortical flow structures to the one shown in Fig. 12. However, the strength and location of the vortices may differ. For test points having an angle of attack α between approximately 12◦ and 16◦ the vortical flow structure changes from that observed for MP0009 to that observed for MP0016. 5.3. Effect of turbulence model To assess the effect of different turbulence models, simulations have been performed on grid G1 (all leading edge flap gaps) for

test point MP0009 (α = 10.05◦ ) and MP0016 (α = 20.06◦ ) employing both the standard TNT k–ω model, the TNT k–ω model extended with the global correction for vortical flows and the Explicit Algebraic Reynolds Stress Model (EARSM). The results of these simulations are compared with wind tunnel data in Figs. 13 and 14 for MP0009 and Figs. 15 and 16 for MP0016. Figs. 13 and 15 compare the sectional surface pressure coefficient at x = 955 mm and x = 1125 mm, whereas Figs. 14 and 16 compare the surface pressure coefficient. Note that these surface pressure signatures correspond to the vortical flow structures shown in Figs. 11 and 12. For test point MP0009 (α = 10.05◦ ) the results of all three turbulence models agree well among each other as well as with the wind tunnel data. Small regions of high suction are simulated and measured along the wing leading edge and in the wing tip region. Note that for all three turbulence models a weak signature of the strake vortex is present, whereas this signature is absent in the PSP data.

O.J. Boelens / Aerospace Science and Technology 20 (2012) 38–51

47

Fig. 14. Surface pressure coefficient for test point VN01004-MP0009 (α = 10.05◦ , M = 0.18 and Rem.a.c . = 2.07 × 106 ) simulated on G1 (all leading edge flap gaps); PSP data (top left figure), TNT k–ω (top right figure), TNT k–ω with vortex correction (bottom left figure) and EARSM (bottom right figure).

Fig. 15. Sectional surface pressure coefficient simulated on G1 (all leading edge flap gaps) for test point VN01004-MP0016 (α = 20.06◦ , M = 0.18 and Rem.a.c . = 2.07 × 106 ) at x = 955 mm (left figure) and x = 1125 mm (right figure). TNT k–ω : blue line, TNT k–ω with vortex correction: green line, EARSM: red line. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

For test point MP0016 (α = 20.06◦ ) the results can be separated in two groups, i.e. the standard TNT k–ω model and the Explicit Algebraic Reynolds Stress Model on the one hand and the TNT k–ω model extended with the global correction for vortical flows on the other hand. At x = 955 mm all three models give a fairly similar result, which agrees well with wind tunnel data. At x = 1125 mm the standard TNT k–ω model and the Explicit Algebraic Reynolds Stress Model produce similar results. The results of the TNT k–ω model extended with the global correction for vortical flows, however, show the best agreement with the wind tunnel data at this section. All three turbulence models show that the strake vortex and the inner wing vortex originating from the wing leading edge interact and merge downstream into a single vortex. This behavior is

not seen in the Pressure Sensitive Paint data, where the strake vortex seems to extend as an autonomous vortex all the way down to the trailing edge and the signature of the inner wing vortex originating from the wing leading edge disappears some distance aft of x = 955 mm. Due to this different behavior the suction peak associated with the merged vortex observed in the simulations is more outboard than the peak due to the strake vortex. The vortex originating at gap1 is predicted fairly similar by the standard TNT k–ω model and the Explicit Algebraic Reynolds Stress model. The TNT k–ω model extended with the global correction for vortical flows results in a stronger vortex resulting in an earlier breakdown. At x = 1155 mm this breakdown leads to a very good match with the wind tunnel measurement, whereas both other turbulence models show too high suction peaks.

48

O.J. Boelens / Aerospace Science and Technology 20 (2012) 38–51

Fig. 16. Surface pressure coefficient for test point VN01004-MP0016 (α = 20.06◦ , M = 0.18 and Rem.a.c . = 2.07 × 106 ) simulated on G1 (all leading edge flap gaps); PSP data (top left figure), TNT k–ω (top right figure), TNT k–ω with vortex correction (bottom left figure) and EARSM (bottom right figure).

Fig. 17. Sectional surface pressure coefficient for test point VN01004-MP0016 (α = 20.06◦ , M = 0.18 and Rem.a.c . = 2.07 × 106 ) at x = 955 mm (left figure) and x = 1125 mm (right figure) simulated with TNT k–ω with vortex correction. G1 (all leading edge flap gaps): blue line, G2 (longitudinal leading edge flap gaps): green line, G3 (no leading edge flap gaps): red line. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Finally, all three models poorly resolve the outer wing vortex originating at gap2. As was shown in Section 2, at this angle of attack this vortex is very sensitive to small changes in angle of attack. The vortex disappears rapidly going from an angle of attack α of 18.04◦ to 22.06◦ . Based on these results it was decided to use the TNT k–ω model extended with the global correction for vortical flows for the investigation of the effect of the flap gaps. 5.4. Effect of leading edge flap gaps Simulations have been performed for the test conditions selected from test run VN01004 (see Table 3) on the grids G1 (all

leading edge flap gaps), G2 (longitudinal leading edge flap gaps) and G3 (no leading edge flap gaps). Figs. 17 and 18 show an example of the surface pressure coefficient obtained on these three grids for an angle of attack α of 20.06◦ . Fig. 19 shows a comparison of the simulated and experimental pitching moment coefficients. The corresponding curves for the lift and drag coefficient are shown in Figs. 20 and 21, respectively. The simulated lift, drag and pitching moment coefficients have been averaged over the last 150 iterations. The error bars indicate the minimum and maximum value during the last 150 iterations. Not all details of these results will be discussed here. However, looking at the results of these simulations the following remarks can be made:

O.J. Boelens / Aerospace Science and Technology 20 (2012) 38–51

49

Fig. 18. Surface pressure coefficient for test point VN01004-MP0016 (α = 20.06◦ , M = 0.18 and Rem.a.c . = 2.07 × 106 ) simulated with TNT k–ω with vortex correction; PSP data (top left figure), G1 (all leading edge flap gaps) (top right figure), G2 (longitudinal leading edge flap gaps) (bottom left figure) and G3 (no leading edge flap gaps) (bottom right figure).

Fig. 19. Comparison of simulated and experimental pitching moment coefficient for test run VN01004 (M = 0.18 and Rem.a.c . = 2.07 × 106 ); simulations with TNT k–ω with vortex correction. G1 (all leading edge flap gaps): blue line, G2 (longitudinal leading edge flap gaps): green line, G3 (no leading edge flap gaps): red line. The simulated data is averaged over the last 150 iterations. The error bars indicate the minimum and maximum value during the last 150 iterations. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

1. For angles of attack α up to approximately 12◦ , there is no visible effect of the different leading edge configurations. The vortical flow structure is similar to the one described in Section 5.2 for MP0009 (α = 10.05◦ ). The sectional surface pressure coefficients coincide and agree well with the experimental data. The comparison of the surface pressure coefficient

Fig. 20. Comparison of simulated and experimental lift coefficient for test run VN01004 (M = 0.18 and Rem.a.c . = 2.07 × 106 ); simulations with TNT k–ω with vortex correction. G1 (all leading edge flap gaps): blue line, G2 (longitudinal leading edge flap gaps): green line, G3 (no leading edge flap gaps): red line. The simulated data is averaged over the last 150 iterations. The error bars indicate the minimum and maximum value during the last 150 iterations. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

with Pressure Sensitive Paint data also shows that the differences are small. The lift, drag and pitching moment coefficients agree well with the experimental data. 2. For angles of attack α larger than 12◦ , the effect of the different leading edge configurations becomes apparent. (a) Grid G3 (no leading edge flap gaps) gives rise to a vortical flow system consisting only of: (i) a strake vortex, (ii) an

50

O.J. Boelens / Aerospace Science and Technology 20 (2012) 38–51

pitching moment coefficients obtained on both grids start to divert, with the results obtained using grid G1 (all leading edge flap gaps) showing the best agreement with experiment. Simulations obtained using grid G3 (no leading edge flap gaps) initially result in a higher pitching moment coefficient than the other two grids. On this grid more suction on the forward region of the wing is obtained due to a forward shift of the inner wing vortex originating from the wing leading edge. The pitching moment coefficient also shows a drop, however at the slightly lower angle of attack α of 16.05◦ . For an angle of attack α between 18.04◦ and 20.06◦ the pitching moment curve of grid G3 coincides with that of the other two grids despite the different flow topologies. 6. Conclusions

Fig. 21. Comparison of simulated and experimental drag coefficient for test run VN01004 (M = 0.18 and Rem.a.c . = 2.07 × 106 ); simulations with TNT k–ω with vortex correction. G1 (all leading edge flap gaps): blue line, G2 (longitudinal leading edge flap gaps): green line, G3 (no leading edge flap gaps): red line. The simulated data is averaged over the last 150 iterations. The error bars indicate the minimum and maximum value during the last 150 iterations. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

inner wing vortex originating from the wing leading edge, and (iii) a weak outer wing vortex originating at the kink between the inner wing and the outer wing. Consequently, the surface pressure distribution differs significantly from the experimental surface pressure distribution. It should, however, be noted that the strength and the location of the strake vortex is predicted well both at x = 955 mm and x = 1125 mm. (b) Grid G1 (all leading edge flap gaps) and grid G2 (longitudinal leading edge flap gaps) give rise to similar results. The vortical flow system resembles the one discussed in Section 5.2 for MP0016 (α = 20.06◦ ). Both the strake vortex and the inner wing vortex originating from the wing leading edge are predicted well, see for example the suction peaks at x = 955 mm. The downstream behavior of the strake vortex and the inner wing vortex originating from the wing leading edge has been discussed in Section 5.3. On both grids the vortex originating at gap1 can be seen. The introduction of the span wise leading edge flap gaps results in a weakening of this vortex. The pressure distribution found on grid G1 (all leading edge flap gaps) shows a better agreement with the experimental data, both close to the leading edge as well as down stream. Finally, also the outer wing vortex originating at gap2 can be observed. This vortex, however, disappears much faster than in the wind tunnel data. At an angle of attack α of 16.05◦ the agreement in the surface pressure coefficient is good, especially on grid G1 (all leading edge flap gaps). Whereas in the experiment a signature of this vortex can be found until an angle of attack α of 19.98◦ , the vortex has almost completely resolved in the simulations for an angle of attack α of 18.04◦ . 3. The differences described above are reflected in the curve of the pitching moment coefficient, see Fig. 19. Up to an angle of attack α of 12◦ the solutions on all three grids coincide. For larger angles the pitching moment coefficient on grid G1 (all leading edge flap gaps) and G2 (longitudinal leading edge flap gaps) show a similar behavior. The curves have shifted to a higher value in comparison with the experiment, but exhibit the correct behavior regarding the drop in pitching moment coefficient. Note that for an angle of attack α of 22.06◦ the

The present article focused on the question to what extend leading edge details, flap gaps, need to be taken into account for the X-31 wind tunnel model to properly simulate the flow around this configuration. This investigation has been carried out in the framework of NATO RTO task group AVT-161 ‘Assessment of Stability and Control Prediction Methods for NATO Air & Sea Vehicles’. Three leading edge configurations have been considered, i.e. one with all leading edge flap gaps, one with only the longitudinal flap gaps and one with no leading edge flap gaps. Steady flow simulations have been performed employing the flow solver ENSOLV, which is part of NLR’s flow simulation system ENFLOW. Results obtained for selected test conditions measured during test run VN01004 (M = 0.18 and Rem.a.c . = 2.07 × 106 ) of the wind tunnel experiments performed in the DNW Low-SpeedWind-Tunnel Braunschweig (DNW-NWB) have been discussed. Based on these results the following can be concluded:

• NLR’s Cartesian grid mapping technique allows the generation

•

•

•

•

of high quality structured (multi-block) grids around complex configurations within very short turn-around times. Due to the use of compound faces, i.e. faces consisting of any number of subfaces, all three leading edge configurations could be handled with the same topology except for the blocks in the flap gaps. The resulting grids which are virtually identical, except for the flap gaps, have been generated in less than two weeks. For angles of attack α smaller than 12◦ no effect of the different leading edge configurations is observed. In addition, all three turbulence models used, i.e. the standard TNT k–ω model, the TNT k–ω model extended with the global correction for vortical flows and the Explicit Algebraic Reynolds Stress model (EARSM) result in identical surface pressure coefficients (see results for test point MP0009 (α = 10.05◦ )). For angles of attack α larger than 12◦ differences due to the different leading edge configurations occur. For these angles of attack α , a highly complex flow consisting of several vortices is observed. The longitudinal flap gaps are essential in the formation of these vortices. The span wise flap gaps mainly play a role in weakening the vortices originating at the longitudinal flap gap between the wing and the inner flap (gap1) and the inner flap and the outer flap (gap2). With no leading edge flap gaps modeled only a strake vortex and a vortex originating from the inner wing leading edge are observed. For angles of attack α larger than 12◦ the TNT k–ω model extended with the global correction for vortical flows clearly shows the best agreement with the experimental wind tunnel data (see results for test point MP0016 (α = 20.06◦ )). Therefore, this turbulence model has been used in all simulations. The drop in the pitching moment coefficient which occurs for test run VN01004 for angles of attack α between 17◦ and 23◦ is mainly due to the behavior of the vortex originating at the

O.J. Boelens / Aerospace Science and Technology 20 (2012) 38–51

longitudinal flap gap between the inner flap and the outer flap (gap2) and the vortex formed at the inner wing leading edge. • Simulations for both configurations with the longitudinal leading edge flap gaps (grid G1 and grid G2) show a similar drop in the pitching moment coefficient for angles of attack α between 17◦ and 23◦ . The simulations for the configuration with no leading edge flap gaps show, however, also a drop in the pitching moment at an angle of attack α of 16◦ . The physics underlying this drop are clearly different from the physics observed in the experiment. In this case, the vortex originating at the longitudinal flap gap between the inner flap and the outer flap (gap2) is not present, and the drop in pitching moment coefficient is solely due to a forward shift of the inner wing vortex originating from the wing leading edge. • Modeling the geometrical details at the leading edge of the X-31 wind tunnel model properly is therefore essential in obtaining correct results. The analysis of the pitching moment coefficient demonstrates how in case of not resolving geometrical details a seemingly correct behavior is obtained without, however, resolving the underlying flow physics correctly. The work presented here is work in progress. In addition to expanding the analysis on leading edge details presented in this article, in the future simulations for the X-31 wind tunnel model in simple harmonic pitch, plunge and yaw motion, as well as complex manoeuvres will be performed. Experimental data for these motions has been obtained by DLR during a wind tunnel campaign at the Low-Speed-Wind-Tunnel Braunschweig (DNW-NWB). Acknowledgements The author would like to thank all the members of NATO RTO working group AVT-161 ‘Assessment of Stability and Control Prediction Methods for NATO Air & Sea Vehicles’ and particularly Dr. A. Schütte and Prof. R.M. Cummings for the interesting and fruitful discussion during the course of this project. In addition, the author would like to thank Dr. B.B. Prananta and Mr. K.M.J. de Cock for proofreading the article.

51

This work has been conducted under NLR’s programmatic research funding “Kennis voor Beleid en Toepassing”. References [1] O.J. Boelens, K.J. Badcock, S. Görtz, S. Morton, W. Fritz, S.L. Karman Jr., T. Michal, J.E. Lamar, Description of the F-16XL geometry and computational grids used in CAWAPI, Journal of Aircraft 46 (2) (2009) 355–368, http://dx.doi.org/ 10.2514/34957. [2] J.W. Boerstoel, A. Kassies, J.C. Kok, S.P. Spekreijse, ENFLOW, a full-functionality system of CFD codes for industrial Euler/Navier–Stokes flow computations, NLR TP 96286U, NLR, Amsterdam, 1996. [3] F.J. Brandsma, J.C. Kok, H.S. Dol, A. Elsenaar, Leading edge vortex flow computations and comparison with DNW-HST wind tunnel data, NLR-TP-2001-238, NLR, Amsterdam, 2001. [4] H.S. Dol, J.C. Kok, B. Oskam, Turbulence modeling for leading-edge vortex flows, AIAA-2002-0843, 2002. [5] R.M. Hall, S.H. Woodson, J.R. Chambers, Overview of the abrupt wing stall program, Progress in Aerospace Sciences 40 (2004) 417–452. [6] R.M. Hall, S.H. Woodson, J.R. Chambers, Accomplishments of the abrupt-wingstall program, Journal of Aircraft 42 (3) (May–June 2005) 653–660. [7] U. Henne, G. Höhler, Chr. Klein, M. Rein, W. Sachs, A. Schütte, Stationäre Druckmessungen mittels PSP und Druckanbohrungen, sowie Kraft- und Momentenmessungen am X-31 RC-Modell im DNW-NWB, DLR, IB 224-2005 A 16, 2005. [8] A. Jirasek, Application of Volterra functions to X-31 aircraft model motion, AIAA-2009-3630, 2009. [9] J.C. Kok, Resolving the dependence on free-stream values for the k–ω , turbulence model, NLR-TP-99295, NLR, Amsterdam, 1999. [10] J.C. Kok, B.I. Soemarwoto, H. van der Ven, X-LES simulations using a highorder finite-volume scheme, in: S.H. Peng, W. Haase (Eds.), Advances in Hybrid RANS–LES Modelling, in: Notes on Numerical Fluid Mechanics and Multidisciplinary Design, vol. 97, Springer, 2008, pp. 87–96. [11] J.C. Kok, H.S. Dol, B. Oskam, H. van der Ven, Exra-large eddy simulation of massively separated flows, NLR-TP-2003-200, 2003. [12] R. Nangia, X-31 vector aircraft, low speed stability & control, comparisons of wind tunnel data & theory (focus on linear & panel codes), AIAA-2009-3898, 2009. [13] M. Rein, G. Höhler, A. Schütte, A. Bergmann, T. Löser, Ground-based simulation of complex maneuvers of a delta-wing aircraft, Journal of Aircraft 45 (1) (January–February 2008) 286–291. [14] A. Schütte, G. Einarsson, A. Raichle, B. Schöning, M. Orlt, J. Neumann, J. Arnold, W. Mönnich, T. Forkert, Numerical simulation of maneuvering aircraft by aerodynamic, flight mechanics and structural coupling mechanics coupling, AIAA2007-1070, 2007.