Abstract

This paper describes a detached-eddy simulation (DES) for the flow over a wall-mounted hump. The Reynolds number based on the hump chord is Re𝑐=9.36Γ—105 with an in-let Mach number of 0.1. Solutions of the three-dimensional Reynolds-averaged Navier-Stokes (RANS) procedure are obtained using the Wilcox π‘˜βˆ’πœ” equations. The DES results are obtained using the model presented by Bush and Mani and are compared with RANS solutions and experimental data from NASA's 2004 Computational Fluid Dynamics Validation on Synthetic Jets and Turbulent Separation Control Workshop. The DES procedure exhibited a three-dimensional flow structure in the wake, with a 13.65% shorter mean separation region compared to RANS and a mean reattachment length that is in good agreement with experimental measurements. DES predictions of the pressure coefficient in the separation region also exhibit good agreement with experiment and are more accurate than RANS predictions.

1. Introduction

Simulation of β€œsteady’’ and β€œunsteady’’ flow of aerodynamic bodies has matured a great deal over the past decade. Aerodynamic performance and flow structures can be predicted with acceptable accuracy except in the complex flow regions of mixing and at off-design conditions near stall. In these regions, flow structures with multiple eddies that mingle and mix are often not predicted well due to inadequate computational grid density and a breakdown of turbulence models.

The focus of the current effort has been to validate the Navier-Stokes code and to further develop the detached-eddy viscosity model using the wall-mounted hump grid. This case proved difficult for most of the participants of the NASA workshop, especially those using two or even three-dimensional RANS techniques [1–3]. Most, if not all, participants using Reynolds-averaged Navier-Stokes predicted too large a separation region, which clearly showed the limitations presented with simply using a RANS procedure. Efforts by Morgan et al. [1] to surpass the limitations of RANS models with high-order spacial and temporal methods applied to the Navier-Stokes and π‘˜βˆ’πœ– equations showed little improvement, especially in the separation region of the flow. Similar results were obtained by Balakumar [3] using a fifth-order accurate weighted essentially nonoscillatory (WENO) scheme for space discretization and a third-order, total variation diminishing (TVD) Runge-Kutta scheme for time integration. Israel et al. [4] evaluated a direct numerical simulation (DNS) as well as a flow simulation methodology (FSM) where the RANS closure equations are combined with a contribution function that is used to determine the magnitude of the turbulent stress term. They observed good agreement between FSM, DNS, and the experimental data and suggest that differences between the numerical predictions and experiment can be attributed to discrepancies in the numerical representation of the experiment rather than deficiencies in the numerical methodology. Miller and Seitz [5] compared a SAS model proposed by Menter, Kuntz, and Bender [6] as well as DES and found that both the methods produced results in good agreement with experiment and prove to be more appropriate for industrial applications since they are less computationally intensive than DNS and LES simulation.

2. Governing Equations

The unsteady, Favre-averaged governing flow-field equations for an ideal, compressible gas in the right-handed, Cartesian coordinate system using primary variables are used. The three-dimensional continuity, momentum, and energy equations can be written in conservative form as follows: πœ•πœŒ+πœ•ξ€·πœ•π‘‘πœŒπ‘’π‘—ξ€Έπœ•π‘₯𝑗=0,πœ•πœŒπ‘’π‘–+πœ•ξ€·πœ•π‘‘πœŒπ‘’π‘—π‘’π‘–ξ€Έπœ•π‘₯𝑗=βˆ’πœ•π‘ƒπœ•π‘₯𝑖+πœ•Μ‚πœπ‘–π‘—πœ•π‘—,πœ•πΈ+πœ•ξ€·πœ•π‘‘πœŒπ‘’π‘—π»ξ€Έπœ•π‘₯𝑗=πœ•πœ•π‘₯π‘—ξ‚Έπ‘’π‘–Μ‚πœπ‘–π‘—+ξ‚΅πœ‡+πœ‡Pr𝑇Prπ‘‘ξ‚Άπœ•β„Žπœ•π‘₯𝑗,Μ‚πœπ‘–π‘—=πœπ‘–π‘—+πœπ‘…π‘–π‘—=ξ€·πœ‡+πœ‡π‘‡ξ€Έπ‘†π‘–π‘—,𝑆𝑖𝑗=12ξ‚΅πœ•π‘’π‘–πœ•π‘₯𝑗+πœ•π‘’π‘—πœ•π‘₯π‘–ξ‚Άβˆ’23πœ•π‘’π‘˜πœ•π‘₯π‘˜π›Ώπ‘–π‘—.(1)

Since we are dealing with compressible flows, we require an equation of state to relate the energy with pressure and enthalpy 1𝐸=πœŒπ‘’+2𝑉2ξ‚„,1𝐻=β„Ž+2𝑉2=π›Ύπ‘ƒπ›Ύβˆ’1𝜌+12𝑉2=𝐸+π‘ƒπœŒ.(2)

Additional governing equations as developed by Wilcox [7, 8] are used for the transport of turbulent kinetic energy and turbulence dissipation rate in regions of the flow where the computational grid or global time-step size cannot resolve the turbulent eddies. The equations for the turbulent mixing energy and specific dissipation rate can be seen below and are Wilcox 2006 [7] model πœ‡π‘‡=πœŒπ‘˜πœ”,(3)πœ•πœŒπ‘˜+πœ•π‘‘πœ•πœŒπ‘’π‘—π‘˜πœ•π‘₯𝑗=πœπ‘…π‘–π‘—πœ•π‘’π‘–πœ•π‘₯π‘—βˆ’π›½π‘˜πœ•πœŒπœ”π‘˜+πœ•π‘₯π‘—ξ‚Έξ€·πœ‡+πœŽβˆ—πœ‡π‘‡ξ€Έπœ•π‘˜πœ•π‘₯𝑗,(4)πœ•πœŒπœ”+πœ•π‘‘πœ•πœŒπ‘’π‘—πœ”πœ•π‘₯𝑗=ξ‚€π›Ύπœ”π‘˜ξ‚πœπ‘…π‘–π‘—πœ•π‘’π‘–πœ•π‘₯π‘—βˆ’π›½πœ”πœŒπœ”2+πœ•πœ•π‘₯π‘—ξ‚Έξ€·πœ‡+πœŽπœ‡π‘‡ξ€Έπœ•πœ”πœ•π‘₯𝑗,(5) where π›½π‘˜=9100,π›½πœ”=3,𝜎40βˆ—=121,𝜎=2.(6)

In regions of the flow where the larger-scale eddies can be resolved with the computational grid, techniques borrowed from large-eddy simulation are used to represent the viscous shear and turbulent viscosity. The large-eddy subgrid model described by Smagorinsky [9] is modified according to the detached-eddy considerations described by Strelets [10] and Bush and Mani [11] πœ‡π‘‡=πœŒπ‘™π‘™π‘’βˆšπ‘˜,(7) where 𝑙𝑙𝑒 is an eddy length scale proportional to the grid/time-step filter width, Ξ”: π‘™π‘™π‘’ξƒ©βˆš=minπ‘˜πœ”,π›½π‘˜πΆdesΞ”ξƒͺ.(8) In addition, the dissipation term, π›½π‘˜πœŒπ‘˜πœ”, of the turbulent kinetic energy transport equation (4) is limited by the eddy length scale, 𝑙𝑙𝑒, according to π›½π‘˜πœŒπ‘˜πœ”β‡’π›½π‘˜ξƒ©βˆšπœŒπ‘˜βˆ—maxπœ”,π‘˜π›½π‘˜πΆdesΞ”ξƒͺ,(9) where 𝐢des is a proportionality coefficient ξ‚€βˆšΞ”=max𝑑π‘₯,𝑑𝑦,𝑑𝑧,π‘’βˆ—π‘‘π‘‘,ξ‚π‘˜βˆ—π‘‘π‘‘.(10)

Equation (10) determines the smallest eddies that can be resolved. The first three terms determine what the minimum eddy size the grid can support. The two time terms are the convection velocity and subgrid scale turbulence, respectively. They are included to ensure that the time step is small enough to resolve the unsteady effects. Equation (7) can be made to resemble (3) using the following: πœ”π΅ξƒ©βˆš=maxπœ”,π‘˜π›½π‘˜πΆdesΞ”ξƒͺ,πœ‡π‘‡=πœŒπ‘˜πœ”π΅.(11)

3. Numerical Techniques

The conservation of mass, momentum, and energy equations are solved using a Lax-Wendroff control-volume, time-marching scheme as developed by Ni [12], Dannenhoffer [13], and Davis et al. [14, 15]. Numerical solutions of unsteady flows can be performed with either the explicit [12] or a dual time-step procedure [16]. These techniques are second-order accurate in time and space. A multiple-grid convergence acceleration scheme [12] is used for steady, Reynolds-averaged Navier-Stokes solutions and the inner convergence loop of unsteady simulations using the dual time-step scheme. The approach is called MBFLO and has two-dimensional [17], axi-symmetric [18] (with and without swirl), and three-dimensional [19] versions. The three-dimensional procedure and results for a DES using the Bush and Mani algorithm for flow over a wall-mounted hump are described here.

The combined second- and fourth-difference dissipation model of Jameson [20] is used in the current procedure for both the mean flow and turbulence equations. The fourth-difference dissipation is scaled by the inverse of the absolute value of the mean strain rate squared. This function decays the numerical dissipation in all viscous flow regions, including boundary layers, wakes, large eddies, and secondary flows.

Parallelization is performed using the Message Passing Interface (MPI) library [21]. Figure 1 and Table 1 show the typical speed-up and associated efficiencies as functions of the number of processors. The configuration used to generate this data was similar to that shown below in the results section where the computational grid consisted of 1.80 million grid points. The data was generated on a Linux cluster consisting of 3.6GHz Intel Xeon processors. Figure 1 shows that a speed-up factor of 18.06 is realized with 20 processes yielding a 90% parallel efficiency. In Table 1, we can see that efficiencies of 90% and higher can be obtained if no less than 100,000 grid points per process are used.

4. Results

The Navier-Stokes codes has been verified with analytical data for a series of standard test cases such as steady inviscid flow over circular bumps, turbulent flow over airfoils, and laminar and turbulent flow over a flat plate, as well as axi-symmetric flows [17–19]. The focus of the current investigation is to demonstrate and validate the DES model for a separated flow and to determine what advantage may exist, in terms of accuracy, for DES compared to three-dimensional URANS and its applicability in the design process.

The simulation of a turbulent flat plate, constructed in such a way as to develop the inflow boundary conditions for the wall-mounted hump, was initially performed. The length of the flat plate was determined through preliminary CFD tests, so the predicted boundary layer thickness was essentially that of the experiment. The flow conditions were also set to match the wall-mounted hump with Ma=0.1 and Re𝑐=9.36Γ—105.

For the flat plate mesh, the computational domain extended from βˆ’2.80≀π‘₯/𝑐≀4.25 in the streamwise direction and from 0.0≀𝑦/𝑐≀0.909 in the normal direction. There were 121 points used in the streamwise direction with a maximum stretching ratio of Ξ”π‘₯=1.2. The wall spacing normal to the surface was set to 3.6Γ—10βˆ’6 and corresponds to Δ𝑦+=0.25. A maximum stretching ratio of 1.2 was achieved with 177 points clustered near the flat plate surface.

The simulation was performed by initializing the inlet with a uniform flow field where the nondimensional velocity components were set to 𝑒=1,𝑣=0, and 𝑀=0 and nondimensional density to 𝜌=1. In the area over the flat plate, the velocity profile was initialized using the log-law profile, and density and energy were recomputed to keep the pressure at the freestream. During the computation, the inflow boundary condition held total pressure and total temperature constant, while at the exit, static pressure was held constant. Along the flat plate surface from the inlet to the leading edge of the flat plate, an inviscid wall boundary condition was imposed with an adiabatic no-slip wall over the entire flat plate surface. The upper wall of the flow domain was also set to an inviscid wall.

To verify the predicted solution, the velocity profile was plotted in terms of the inner layer velocity and length variables 𝑒+ and 𝑦+, respectively, and compared to the empirical viscous sublayer, log-law, and 1/7th power law relations as shown in Figure 2(a). Figure 2(b) shows the predicted velocity and experimentally measured profiles at π‘₯/𝑐=βˆ’2.8, which corresponds to the inlet of the wall-mounted hump.

Once the inlet conditions were generated using the flat plate, the wall-mounted hump test case for turbulent flow was simulated and compared to experimental data from the NASA workshop [22] that focused on synthetic jets and turbulent separation control. The hump geometry was constructed to simulate a 20% thick Glauert-Goldschmied airfoil with a chord length of 0.472m (1.38ft.), a maximum height of 0.054m (0.176ft.), and a span of 0.584m (1.92ft.). The experimental data was obtained for Ma∞=0.10 and Re𝑐=9.36Γ—105 based on the chord. The test case that was investigated was without flow control, where there was no blowing or suction in the slot. The numerical results obtained used a three-dimensional structured grid smoothed over the slot and with the top wall shape adjusted to approximately account for side plate blockage, as recommended by the workshop [22]. The grid shown in Figure 3 extended from π‘₯/𝑐=βˆ’9.20 to 3.5 in the streamwise direction, with the hump located from π‘₯/𝑐=0.0 to 1.0 and from 𝑦/𝑐=0.0 to 0.9 in the normal direction. The mesh contained 155,937 grid points with 881 in the streamwise direction and 177 normal to the wall. An initial RANS simulation was performed using three planes in the spanwise direction for a grid that consisted of 467,811 points. This allowed a nominally two-dimensional flow simulation to be performed and compared to the experiment. For true three-dimensional simulations, the original 881Γ—177 grid described was used and extruded in the spanwise direction 0.15 chord lengths and can be seen in Figure 4. The spanwise direction was meshed using 49 points with uniform spacing and corresponds to roughly two boundary layer thicknesses based on the inlet.

With the aid of the turbulent flat plate simulation, the inflow velocity and density profiles for the wall-mounted hump case were generated and used to initialize the flow domain over the wall-mounted hump and helped decrease the computational time required to obtain a solution. The top of the flow domain was set with an inviscid wall boundary condition with an adiabatic no-slip wall along the south boundary starting at the inlet and leading over the hump section. Static pressure was once again held constant at the exit of the flow domain. For the nominally two-dimensional case using the steady RANS solver, inviscid walls were imposed for the spanwise boundaries, whereas periodic boundary conditions were used for the three-dimensional unsteady simulation. The turbulent freestream intensity, Tu, and the dissipation length scale, 𝑙dis, were set to 1.00% and 0.06, respectively. The DES coefficient used was that suggested by the original authors of a value of 0.6 [11].

The temporal periodicity and unsteady behavior was initially studied for the three-dimensional DES case. The information obtained was then used to determine the number of time-steps necessary to resolve a minimum of 10 periodic cycles for the time-averaged DES at the given global time-step. Figures 5 and 6(a) show the signal history for the instantaneous π‘₯-, 𝑦-, and, 𝑧-surface forces as well as the Power Spectral Density (PSD) as a function of the number of time-steps. The PSD of surface forces was plotted here as a function of the number of time-steps to more readily see periodicity in the flow. In Figure 6(a), we see that the peak spectral signal is repeated approximately every 400 time-steps, which corresponds to a frequency of 420Hz. The corresponding Strouhal number based on the hump height and freestream speed is approximately 0.65. Once the period was determined, this case was run for approximately 10 periodic cycles and time-averaged.

Figure 7 shows the time-averaged pressure coefficient distribution resulting from the RANS and time-averaged DES simulations. The numerical results predicted using the RANS procedure produced higher pressure levels in the separation region, 0.65≀π‘₯/𝑐≀1.3, when compared to the experimental data. This is consistent when compared to other RANS simulations [1–3] in the NASA workshop. Researchers who used shear-stress transport (SST) RANS models showed similar predictions for the pressure leading up to the hump as well as higher pressures in the separation region. The baseline RANS cases run by Morgan, Rizzetta, and Visbal [1] also showed similar pressure predictions for the flow leading up to the hump and in the separation region. RANS simulations run by Krishnan, Squires, and Forsythe [2] also showed similar predictions as well as those made by Balakumar [3] and Ε ariΔ‡ et al. [23]. The current DES case showed similar comparison with the experimental data in the acceleration region of the flow, 0.00≀π‘₯/𝑐≀0.65, with significant improvement in the separation region located between 0.65≀π‘₯/𝑐≀1.3.

A comparison of the time-averaged skin-friction coefficient, 𝐢𝑓, is shown in Figure 8. Once again, there is good agreement with the steady RANS prediction upstream of the hump to π‘₯/𝑐=0 and over the hump to the beginning of the separation region π‘₯/𝑐=0.65. The DES case slightly underpredicts the skin friction in the acceleration region of the flow with a max value of approximately 0.005. Similar results were shown by Morgan et al. [24] using an implicit large-eddy simulation. Both simulations accurately predicts the onset of the separation region at approximately π‘₯/𝑐=0.65. In the separation region, we see the predicted reattachment point for the RANS case is at π‘₯/𝑐=1.2 and the time-averaged DES at π‘₯/𝑐=1.105. The DES matches the experimental data in the separation region much better than the RANS procedure. Table 2 shows an overview of the separation and reattachment locations for the RANS and DES procedures compared to the experimental data. Here, we can see that both procedures accurately predict the onset of separation with the RANS procedure overpredicting the size of the separation bubble. The DES, however, properly predicts the reattachment point and the separation bubble size.

Figure 9 shows a stream trace velocity comparison. Here we can more easily see the flow reattachment locations and the effect the DES has on the flow. In Figure 9(b), we see that the streamline plot clearly shows that the RANS solution has a longer separation bubble than that observed experimentally in Figure 9(a). The time-averaged DES streamlines, Figure 9(c), show significant improvement for the predicted separation bubble length, with a 13.65% shorter mean separation region compared to RANS and a mean reattachment length that is in good agreement with experimental measurements.

Time-averaged velocity profiles at π‘₯/𝑐=1.0,1.1,1.2, and 1.3, corresponding to locations within the separation region and slightly downstream of it, were also compared with experimental data. The peak reverse flow velocity predicted at π‘₯/𝑐=1.0 in Figure 10(a) is close to the experimentally measured velocity and slightly lags in the region away from the wall. Referring to the velocity profiles at π‘₯/𝑐=1.1 and 1.2 in Figures 10(b) and 10(c), it can be seen that the experimentally measured flow reattachment point is at π‘₯/π‘β‰ˆ1.1, while MBFLO RANS and DES procedures predicted the reattachment point at π‘₯/π‘β‰ˆ1.2 and 1.105, respectively. In the near wall region of the flow, the DES procedure more accurately compares with the experimental velocity profiles at all four locations.

The DES and RANS profiles of the 𝑣-velocity component at π‘₯/𝑐=1.0 and 1.1 in Figure 11 show good qualitative comparison at π‘₯/𝑐=1.0 and 1.1 and match the experimentally measured values very well at π‘₯/𝑐=1.2 and 1.3. It should be noted that the 𝑣-component of the velocity is a magnitude smaller than the 𝑒-component and small changes can more readily be seen.

Figure 12 shows contour plots comparing the spanwise spatially averaged vorticity and for the RANS and instantaneous DES. The instantaneous vorticity contours in Figure 12(b) show a large range of resolved eddies consistent with large-eddy simulation treatment.

Instantaneous spanwise vorticity contour slices at 𝑦=3.6576, 4.2672, and 5.4864cm. normal to the wall are shown in Figure 13. Here, we can more easily see the three-dimensionality that has formed in the wake of the wall-mounted hump. Figure 14 presents instantaneous vorticity isosurfaces for the DES prediction. Here we can clearly see the separated shear layer downstream of the hump.

The detached-eddy simulation procedure, described by Bush and Mani [11], has been shown to be consistently more accurate than standard RANS. It predicts well the experimentally measured flow quantities such as pressure coefficient, surface skin friction, reattachment length, and mean velocity profiles. The RANS procedure predicted a delayed reattachment point, which indicates reduced turbulent mixing inside the separation region. Attempts at using higher-order numerical techniques, applied to RANS procedures [1, 3], have shown similar results. It should be noted that the DES procedure used a single constant model coefficient of 0.6; although additional computations using various coefficients would be essential to further examine the DES procedure, it was beyond the scope of the current investigation.

5. Conclusion

A general Reynolds-averaged/detached-eddy simulation procedure was applied to the prediction of the flow over a wall-mounted hump. The initial results using the RANS and DES procedures compared well with experiment as well as other participants of NASA’s 2004 Computational Fluid Dynamics Validation Workshop. Like other participants using RANS models, the onset of separation was accurately predicted and the reattachment point was overpredicted. The RANS procedure also overpredicted the mean pressure, skin friction, and velocity profiles in the separation zone. The DES procedure using the Bush and Mani model showed much better results. The three-dimensional structures resolved in the wake of the DES improved the local flow physics in the separation region and the predictions of the mean pressure distribution, skin friction, and streamwise velocity.

Nomenclature

E:Total energy
π‘’βˆΆInternal energy
𝐻∢Total enthalpy
β„ŽβˆΆStatic enthalpy
π‘˜βˆΆTurbulent kinetic energy
π‘βˆΆPressure
Pr∢ Prandtl number
Pr𝑑:Turbulent Prandtl number
𝑆𝑖𝑗:Mean strain-rate tensor
𝑒𝑖:Velocity component
π‘‰βˆΆVelocity magnitude
π›½π‘˜,πœ”: Constants defined by π‘˜βˆ’πœ” equations
πœ‡βˆΆCoefficient of viscosity
πœ‡π‘‡:Turbulent coefficient of viscosity
πœ”βˆΆTurbulent dissipation frequency
πœ”π΅: DES filtered turbulent dissipation frequency
Μ‚πœπ‘–π‘—:Total shear stress tensor
πœπ‘–π‘—:Laminar shear stress tensor
πœπ‘…π‘–π‘—: Reynolds shear stress tensor
𝜌∢Density
Tu∢ Freestream turbulence level as %
𝑙𝑙𝑒:Eddy length scale
Ξ”βˆΆ DES filter width
𝜎,πœŽβˆ—: Constants defined by π‘˜βˆ’πœ” equations
𝐢des: DES coefficients.

Acknowledgments

The authors would like to thank Dr. John Clark and the managers of the turbine branch at the Wright-Patterson Air Force Research Laboratory in Dayton, Ohio, for their support of this effort under contract 09-S590-0009-20-C1. Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-AC04-94AL85000.