Budget of Turbulent Kinetic Energy in a Shock Wave Boundary-Layer Interaction

ρ density ρ counterflow force model charge density c τ , τ Reynolds-stress tensor ij θ boundary-layer momentum thickness θ frequency of applied voltage f ϑ counterflow force effective duty cycle ζ second viscosity C specific heat at constant pressure p D counterflow force model coefficient—ratio of electrical to inertial forces c e total energy e electronic charge c k turbulent kinetic energy L length scale M Mach number p pressure P r Prandtl number Re Reynolds number S mean strain-rate tensor ij s instantaneous strain-rate tensor ij T temperature t time t viscous stress tensor ij u , v , w non-dimensional velocity in x , y , z directions, respectively V velocity magnitude x, y, z non-dimensional coordinates Subscripts 0 total condition i , j , k index notation, equal to 1, 2, 3 inv inviscid ref , ∞ reference value same as freestream value vis viscous w wall Conventions ¯ Reynolds averaged ˜ Favre averaged Superscripts ′ Reynolds decomposition ′′ Favre decomposition ∗ non-dimensional value + non-dimensionalized by inner scales tr transpose I. Introduction Shock wave/boundary-layer interactions (SBLIs) have a rich history of research and advancement over 1–3 the course of the last 60 years in the areas of experimentation, modeling, and simulation. SBLIs are ubiquitous in transonic, supersonic, and hypersonic flight, and thus relevant to the majority of commercial, military, and space vehicles of the past, present, and future. They can be found on the external body, like the vehicle nose, wings, and/or tail, and internal to the propulsion flowpath, which includes mixed compression inlets, diffusers, and isolators. Aerothermal loads created by SBLIs can compromise structural integrity, cause component failure, and result in loss of control. Adverse pressure gradients in the propulsion flowpath can cause flow distortion and lead to loss in engine efficiency. Because of these characteristics, study of SBLIs is one of the most active areas of research in high speed flows. But despite the collection of work pursued in 2 of 20 American Institute of Aeronautics and Astronautics industry, academia, and government laboratories, a complete understanding of SBLIs remains illusive, and hence their prediction via modeling and simulation is difficult.

The challenges in modeling and simulation are multifold—SBLIs are unsteady in nature, the resulting separation is highly three-dimensional (3D), and corners often separate (e.g., propulsion flowpaths are often rectangular). A combination of these lead to a locally inhomogeneous flowfield. A prominent debate in SBLI research involves the source of low-frequency unsteadiness of the reflected (or separation) shock, which is typically one or two orders of magnitude lower than the frequencies within the incoming turbulent boundary 1, 4, 5 layer. Although much experimental and computational work has been done to study the source of this unsteadiness no clear consensus has been reached. Questions remain: Is the source of low-frequency 6 7, 8 unsteadiness located upstream or downstream of the interaction? Or is it a coupled system influenced by sources present upstream, downstream, and within the separation bubble? An elementary understanding of separation and corner influence exists, however modeling improvements from such understanding has yet to be realized. It is this uncertainty in the knowledge of SBLI physics which make them extremely difficult to model and simulate. An example of a canonical two-dimensional (2D) SBLI flowfield is shown in figure 1. Fluctuations in the incoming turbulent boundary layer can be classified as upstream sources, while the unsteadiness of the near-wall/viscous portion of the reflected shock and separation bubble can be classified as downstream sources.

Figure 1. Two-dimensional anatomy of a shock wave/boundary-layer interaction In most practical applications, Reynolds-averaged Navier-Stokes (RANS) computational fluid dynamics (CFD) solvers coupled with turbulence models are used to calculate flowfields where SBLIs are present.

Although as the name suggests, Reynolds averaging or taking a mean of the Navier-Stokes (NS) equations renders them futile, by definition, if the interest lies within exploring the low-frequency unsteadiness driving the reflected shock. Thus, RANS CFD is impervious to unsteadiness by its fundamental design, which could be a factor affecting SBLI predictions. Nonetheless, RANS CFD has been widely used in industry and academia to obtain SBLI predictions due in part to its ease of use, and also to the challenges presented by scale-resolving methods like hybrid RANS/large-eddy simulations (LES), LES, and direct numerical simulations (DNS) in the form of available computer resources and lengthy simulation times.

Yet another approach would be unsteady RANS (URANS), where a global timestep is used to march the solution forward. If it is the low-frequency unsteadiness of the reflected shock that we wish to study, URANS can possibly capture such phenomenon as the reflected shock frequency is on the order of hundreds 1, 7 of Hertz. The incoming turbulent boundary layer as well as the separation bubble exhibit frequencies in 1, 9 the range of 10-40 kHz. Since such a separation of scale exists between the two, in theory, URANS can capture the low-frequency unsteadiness. However, URANS would still require the use of turbulence models, which account for large variations in predictions. Thus, URANS would not be fundamentally superior to RANS.

3 of 20 American Institute of Aeronautics and Astronautics

A description of scale-resolving approaches is provided by Georgiadis et al. in a paper that provides a summary of current practices in LES. In LES, large-scale structures are resolved and a sub-grid scale model is employed to model the scales which cannot be resolved by the mesh. A subset of LES is implicit LES (ILES). Like LES, ILES calculates the large-scale turbulent structures, but it does not explicitly model the smallest scales, instead it uses a high-order low-pass filter to dissipate energy in the high spatial wavenumber 12, 13 range where the turbulent energy spectrum is poorly resolved. So, the use of an explicit sub-grid scale model is completely avoided. Thus, ILES is an attractive approach for this work as it is not as expensive as DNS, but does provide a seamless changeover to DNS as the mesh resolution is refined.

A workshop was organized by the American Institute of Aeronautics and Astronautics (AIAA) with an intent to share, assess, and determine the most promising SBLI prediction methods in 2011. The results obtained by the participants were compiled by DeBonis et al. in which comparisons with experimental data and error metrics were presented. While the majority of solutions were obtained with RANS methods, some were obtained with hybrid RANS/LES, LES, and DNS. It was concluded that the turbulence model played a significant role in variations among different RANS solutions . In general, the Menter Shear Stress 14 15 Transport (SST) and Wilcox models, which are k − ω based two-equation models, and the one-equation Spalart-Allmaras (SA) model, compare well with each other, however the error in all solutions increased as the adverse pressure gradient becomes stronger and the size of the separation increased . It was also found that the scale-resolving methods provided some of the best and the worst solutions, clearly indicating that scale-resolving methods are feasible and accurate but that more development is needed.

Perhaps the most interesting revelation Debonis et al. presented was the fact that the accuracy of a solution was not consistent for different variables of interest within the same solution, i.e. high prediction accuracy in ¯ u velocity did not guarantee the same accuracy for ¯ v velocity. This, combined with the other observations above, shows shortfalls of the current one- and two-equation turbulence models. Common among most turbulence models in use today is the Boussinesq eddy-viscosity approximation, which establishes a linear relation between the Reynolds-stress tensor, τ , and the mean strain-rate tensor, S . Such a ij ij relationship does not exist in shock-separated flows such as SBLI, flows where rapid changes in mean strain- rate occur, and where secondary flows are present—all examples of flows which are inhomogeneous and anisotropic.

Some efforts have been made to incorporate the above effects in the existing turbulence models by modifying the turbulent kinetic energy, k , and dissipation rate,  , equations to account for inhomogeneity and anisotropy. Hamlington and Dahm addressed this by replacing the mean strain-rate tensor in a standard two-equation approach with a new mean strain-rate tensor that accounts for flow history, thus allowing the Reynolds-stress tensor to adjust over a finite lag while retaining the simplicity of a two-equation formulation.

Sinha et al. included additional terms representing shock unsteadiness in the k and  equations. A linear analysis was used to model the unsteadiness in the k −  model by realizing that a positive fluctuation in the streamwise velocity leads to the reflected/separation shock motion downstream, while a negative fluctuation in the streamwise velocity causes the reflected/separation shock to move upstream. The improved turbulence model predicted k over a Mach number range of 1.29-3.0 reasonably well and showed improvement over a ˜ ′′ ′′ realizable k −  model with a realizability constraint of 0 ≤ u u ≤ 2 k . Numerous other improvements have been suggested to the one- and two-equation turbulence model formulations, though discussing them is beyond the scope of this paper. Another approach is to completely bypass the linear relationship between the Reynolds-stress tensor and the mean strain-rate tensor in one- and two-equation formulations for nonlinear constitutive relations. In some of the more advanced techniques, a direct prescription of the Reynolds-stress tensor is sought using nonlinear algebraic equation or by using Reynolds-stress transport model. However, among these approches, none has emerged as a clear winner yet.

It is the authors’ belief that a physical understanding of the various terms in the exact equation of turbulent kinetic energy transport and their interactions with each other would shed light on the fundamental mechanisms present in SBLI. This knowledge may be used to improve the current turbulence models and/or propose new models. In the past, such efforts involved DNS studies of the turbulent kinetic energy budget on 18 19 flat plate boundary layers or channel flows with no direct comparison with SBLI flowfields. In the present work, we take a step back from modeling and use ILES to study the budget of turbulent kinetic energy and other relevant turbulence quantities with an objective of informing future turbulence model development for better predictions of SBLIs.

4 of 20 American Institute of Aeronautics and Astronautics II. Governing Equations In this section the turbulent kinetic energy transport equation will be discussed along with the equations pertinent to the ILES formulation and the counterflow force model, which is used to obtain a turbulent boundary layer.

II.A. Turbulent Kinetic Energy The turbulent kinetic energy is defined as ˜ ′′ ′′ k = u u (1) i i The transport of turbulent kinetic energy, a standard part of two-equation turbulence models, is given by equation 2. The first term on the left-hand side represents the unsteady term, while the second term represents the convection—together, the left-hand side is the substantial derivative. The budget terms are on the right-hand side.

∂ (¯ ρk ) ∂ (¯ ρ ˜ u k ) j + = P + T + D − ¯ ρ + D + Π + M (2) ν p ∂t ∂x j Each term on the right-hand side of equation 2 is defined as ∂ ˜ u i ˜ ′′ ′′ P = − ¯ ρ u u Production (3) i j ∂x j ( ) ′′ ′′ ′′ T = − ρu u u Turbulent Transport (4) j i i ,j ( ) ′′ D = t u Molecular Diffusion (5) ν ij i ,j ′′ ¯ ρ = t u Dissipation (6) ij i,j ( ) ′′ ′ D = − p u Pressure Diffusion (7) p i ,i ′′ ′ Π = p u Pressure Dilatation (8) i,i ′′ M = u (¯ t − ¯ p ) Mass Flux (9) ij,j ,i i The production term represents the rate of transfer of kinetic energy from the mean flow to the turbu- lence or alternatively, generation of turbulent kinetic energy, which results from the mean velocity gradients in the flowfield. The propagation of the turbulent kinetic energy by the means of triple correlation of velocity fluctuations is given by the turbulent transport term. The molecular diffusion term, also referred to as viscous diffusion, is responsible for molecular transport of the turbulent kinetic energy. The conversion of the turbulent kinetic energy to thermal internal energy is represented by the dissipation term. The pressure diffusion term is a form of transport occurring due to diffusion mechanism as a result of pressure and velocity-gradient interaction. Finally, the compressibility terms, pressure dilatation and mass flux , appears in the compressible form of the turbulent kinetic energy transport equation.

Here t is the viscous stress tensor based on the instantaneous strain-rate tensor s . And δ is the ij ij ij Kronecker delta.

∂u k t = 2 μs + ζ δ (10) ij ij ij ∂x k In equation 10, ζ is obtained by relating it to μ . Such an assumption is valid for monatomic gases and widely used in computational fluid dynamics.

ζ = − μ (11) 5 of 20 American Institute of Aeronautics and Astronautics II.B. Compressible Navier-Stokes Equations The compressible Navier-Stokes equations in non-dimensional form are given by ∂ρ + ∇ · ( ρ V ) = 0 (12) ∂t ∂ ( ρ V ) 1 + ∇ · ( ρ VV + p I ) − ∇ · τ = D ρ E (13) c c ∂t Re [ ] ∂ ( ρe ) 1 1 + ∇ · ( ρe + p ) V − ( V · τ ) − Q = D ρ V · E (14) 0 c c ∂t Re ( γ − 1) P rM Re The velocity, conduction heat transfer, and electric field vectors are defined as       u Q E x x       V = Q = E = (15)  v   Q   E  y y w Q E z z The right-hand side of the equations 13 and 14 represent source terms for the dielectric barrier discharge actuator to be discussed in section II.C. The non-dimensionalization is performed using the following defini- tions where * represents non-dimensional quantities. Except equations 16a and 16b, the * has been dropped for simplicity and all quantities are non-dimensional in this paper, unless stated otherwise.

tV ρ p T ref ∗ ∗ ∗ ∗ t = ρ = p = T = (16a) L ρ ρ V T ref ref ref ref μ u v w ∗ ∗ ∗ ∗ μ = u = v = w = (16b) μ V V V ref ref ref ref The reference conditions are the freestream conditions and the length scale, L , is the same as the boundary-layer thickness, both are discussed later in section III.B. The non-dimensional parameters Reynolds number, Prandtl number, and Mach number are defined below. The specific heat at constant pressure is C and κ is the thermal conductivity constant. The molecular viscosity, μ is calculated using the p ref ref Sutherland’s law and a perfect gas is assumed.

ρ V L μ C V ref ref ref p ref Re = P r = M = √ (17) γp ref μ κ ref ref ρ ref The above Navier-Stokes equations can be expressed in flux-vector form as ∂ U ∂ F ∂ G ∂ H ∂ F ∂ G ∂ H inv inv inv vis vis vis + + + = + + + S (18) ∂t ∂x ∂y ∂z ∂x ∂y ∂z tr where U is the solution vector defined as U = [ ρ, ρu, ρv, ρw, ρe ] , the inviscid and viscous flux vectors are defined as     ρu         τ ρu + p xx         F =   F =  τ  (19) inv ρuv vis xy     Re     τ ρuw xz     1 ∂T ( uτ + vτ + wτ ) + ( ρe + p ) u xx xy xz 2 ( γ − 1) P rM ∂x 6 of 20 American Institute of Aeronautics and Astronautics The above equations are transformed into curvilinear coordinates and using the strong conservation form 20, 21 the following is obtained ˆ ˆ ˆ ˆ ˆ ˆ ˆ ∂ U ∂ F ∂ G ∂ H ∂ F ∂ G ∂ H inv inv inv vis vis vis ˆ + + + = + + + S (20) ∂t ∂ξ ∂η ∂ζ ∂ξ ∂η ∂ζ ˆ ˆ such that U = U / J and S = S / J where J represents the Jacobian of the transformation. The transformed flux vectors are defined as ˆ F = ( ξ F + η G + ζ H ) (21) inv x inv y inv z inv J ˆ F = ( ξ F + η G + ζ H ) (22) vis x vis y vis z vis J ˆ ˆ ˆ ˆ and similarly G , H , G , and H .

inv inv vis vis II.C. Counterflow Force Model The counterflow force model, which models the effect of dielectric barrier discharge (DBD) actuator, is used to trip the laminar boundary layer to turbulent. The model was developed by Shyy et al. and later 23 24, 25 used by Gaitonde et al. in a subsonic flow control application. Mullenix et al. and Waindim et 26, 27 al. used it to generate supersonic turbulent boundary layers on a flat plate for use as the inflow to shock wave/boundary-layer interaction problems.

The source terms, on the right-hand side of the Navier-Stokes equations, represent the effect of a DBD actuator which is composed of a submerged electrode excited with 5kV of AC at 5 kHz as described by Shyy et al. The model development was rooted in experimental work and thus the model is empirical in nature. The original intent of the model was to simulate flow control in numerical calculations, however it is employed as a tripping mechanism to obtain a turbulent boundary layer in the present work.

One of the key parameters employed in the above model is the counterflow force coefficient, defined as the ratio of electrical and inertial forces.

ρ e EL c c ref D = (23) c ρ U ref ref Here ρ is the discharge density and e is the electronic charge. The counterflow force is related to the c c coefficient as F = D θ ϑαρ E δ (24) e c f c cr Here θ is the frequency of the applied voltage, ϑ is the duty cycle to recover an effective mean force, α f is a factor of collision efficiency, and δ is a function of value one where the local electric field exceeds the cr breakdown electric field and zero elsewhere. The magnitude of the electric field vector varies linearly as | E | = E + k x + k y (25) 0 1 2 E = V /d , where V is the applied voltage and d is the streamwise distance between two electrodes. Thus, 0 e e the formulation allows for a peak value in the near-wall region and diminishing to zero at a specified distance away from the wall. The constants k and k can be calculated using a and b , the extents of the force field 1 2 in y and x direction respectively. The various parameters discussed above are listed in table 1 and remain unchanged from the work of Mullenix et al.

1 − E 1 − E 0 0 k = k = (26) 1 2 b a 7 of 20 American Institute of Aeronautics and Astronautics

Table 1. Counterflow force model parameters Parameter Value D 25.0 c θ 3000.0 f − 6 ϑ 67 × 10 α 1.0 ρ 1.0 c E 0.3125 a 0.03125 b 0.125 III. Simulation Methodology III.A. Numerical Schemes A fifth-order bandwidth- and order-optimized weighted essentially non-oscillatory (WENO) scheme was used with Roe’s formulation to compute the inviscid fluxes. The viscous fluxes were calculated using a sixth- order compact spectral-like scheme. The implicit time integration was performed with the second-order Beam-Warming scheme using two sub-iterations and approximate factorization. The non-dimensional time step for this problem was 0.002.

III.B. Boundary Conditions and Mesh Parameters The flow conditions represent the experiments performed at the Ohio State University’s Gas Dynamics and Turbulence Lab. Table 2 shows a comparison of experimental and simulated flow conditions. It is noted that to resolve the scales of turbulence corresponding to the experimental Reynolds number, a mesh various orders of magnitude larger than the fine mesh used here would be required. Thus, the simulated Reynolds number was lowered to be practical and achievable by ILES. Such an approach was also taken by Mullenix 25 31 32 and Gaitonde, Touber and Sandham, and Visbal et al. In order to have consistent flow relations, the pressure was reduced by a factor of 10 while maintaining the temperature. This assures that the Mach number and the velocity scales remain the same between the experiment and the simulation. This was done so that the frequencies associated with the velocity scales remain comparable between the experiment and the simulation. It will be discussed later that the reduction of Reynolds number by a factor of 10 might not be sufficient and further scaling is necessary for the present flow conditions and meshes.

Table 2. Flow conditions Property Experiment Simulation M 2.33 2.33 ∞ U , m/s 556.0 556.0 ∞ P , Pa 23,511.0 2351.1 ∞ T , K 295.6 295.6 T , K 269.7 269.7 w − 3 − 3 δ , m 5 . 3 × 10 5 . 3 × 10 Re 175,202.0 17,520.2 δ At the inflow, a laminar boundary layer profile, calculated using a compressible Blasius solution was imposed. The interior of the flowfield was also initialized with the same Blasius solution. The wall was treated as no-slip and the ratio of expected adiabatic wall temperature to freestream static temperature was 1.95. The outflow and farfield boundaries were obtained by extrapolating values from the interior. A periodic boundary condition was applied in the z direction. The Rankine-Hugonoit relations were used to impose the oblique shock generated by a wedge in the experiment. Such a simulated shock was achieved by 8 of 20 American Institute of Aeronautics and Astronautics manintaining the farfield boundary condition for x < x at the pre-shock conditions identified in table shock 2, while x > x were set to conditions obtained via oblique shock relations at Mach 2.33 and angle of shock attack set to 9 degrees. The shock location, x , was set to 23.19 to obtain the inviscid shock impingement shock location of x ∼ 62 . 0.

Figure 2. Mesh topology for SBLI simulation Table 3. Mesh parameters Coarse Medium Fine Domain size a x × y × z 95 × 25 × 5 95 × 25 × 5 95 × 25 × 5 Computational points N × N × N 917 × 167 × 133 1101 × 201 × 167 1301 × 251 × 201 x y z 7 7 7 N 2 . 025 × 10 3 . 696 × 10 6 . 564 × 10 total Counterflow force region N 100 100 100 x x 2.6 2.6 2.6 center x 1.0 1.0 1.0 length ∆ x 0.010 0.010 0.010 Constant region N 666 833 1000 x ∆ x 0.0997 0.0846 0.0665 + ∆ x 23.51 20.55 16.35 x =65 N 101 125 151 y,bl − 3 − 3 − 3 ∆ y 1 . 96 × 10 1 . 96 × 10 1 . 96 × 10 w + ∆ y 0.4621 0.4764 0.4822 x =65 ∆ z 0.0376 0.0299 0.0250 + ∆ z 8.86 7.26 6.15 x =65 a x, y, z are non-dimensionalized by δ In figure 2, the mesh schematic is presented where the coordinates are non-dimentionalized by the length scale L . For the purpose of clarity, every fourth point is shown. The mesh can be divided into five distinct sections: 1) the inflow section where the mesh is refined until it reaches the counterflow force section, 2) the counterflow force section where a constant axial spacing is maintained with the actuator located at x = 2 . 6, 3) a coarsening section, 4) the constant spacing section is where the SBLI occurs, and 5) the coarsening to outflow section. All sections have the same hyperbolic tangent grid clustering in the y direction and a constant spacing in the z direction. Table 3 shows a list of parameters for coarse, medium, and fine 9 of 20 American Institute of Aeronautics and Astronautics

meshes. Although all three meshes were investigated, only the fine mesh results are presented here. The mesh sensitivity study of Mullenix and Gaitonde was the basis for this decision.

IV. Results In this section, the validation of the upstream boundary layer is presented first and followed by detailed turbulent kinetic energy budgets of the SBLI region.

IV.A. Validation of upstream boundary layer In order to obtain a qualitative description of the developing turbulent boundary layer, the Q-criterion, as described by Jeong and Hussain, was computed and the iso-surfaces of Q = 0 . 5 are presented in figure 3 colored by ¯ u velocity. The effect of the counterflow force actuator can be seen at the beginning of the computational domain followed by the laminar-to-turbulent transition phase, which results in coherent structures. Among the coherent structures are hairpin vortices, figure 4, peaking from the smaller structures in the boundary layer. The change in the boundary layer is noticeable beyond the shock wave/boundary-layer interaction. Towards the end of the computational domain the structures disappear due to grid coarsening.

In a mesh refinement study, Mullenix et al. has shown a similar behavior at Q = 1 . 0 for a variety of meshes, where it was also noted that refining of the mesh led to a higher concentration of resolved coherent structures and a quicker transition process.

Figure 3. Iso-surface of Q-criterion at Q = 0.5 colored by ¯ u velocity The remainder of the results that follow are time averaged for six flow-through times. The time-averaged results were compared three flow-through times apart to ascertain statistical convergence of mean quantities and agrees with observations of Bisek, where a similar approach was used to determine convergence. The skin friction through the laminar-to-turbulent boundary layer transition is presented in figure 5(a) along with 35, 36 the theoretical laminar and turbulent flat plate skin friction profiles. It is clear that the centerline skin 10 of 20 American Institute of Aeronautics and Astronautics Figure 4. Hairpin vortices, Iso-surface of Q-criterion at Q = 0.5 colored by ¯ u velocity friction coefficient oscillates about the theoretical turbulent curve beyond x = 60, however the span-averaged skin friction was marginally lower. It is noted that in the work of Bisek a similar comparison was presented for three mesh levels, coarse, medium, and fine at approximately the same flow conditions. As the mesh was refined, a quicker transition to turbulent skin friction and a better match with the theoretical curve was obtained. Since the fine mesh results presented here are obtained with half the coarse mesh resolution used by Bisek, it is concluded that a finer mesh would provide an improvement over the current prediction, but at a significant additional computational cost. A comparison of the centerline and span-averaged velocity profiles with the inner layer and the logarithmic law of the wall is provided in figure 5(b). The velocity was transformed to the incompressible form using the van Driest transformation. Both the centerline and span-averaged profiles at Re = 3270 show the expected behavior of the turbulent boundary layer and agree θ ¨ with the highly-resolved (3 billion point) incompressible DNS simulation of Schlatter and Orl¨ u. The value Re = 3270 represents x = 56 . 0, the station upstream of the area where the SBLI, to be discussed in section θ IV.B, occurs. Thus, this location was an important landmark to validate the developed boundary layer before the SBLI region.

¨ Next, the Reynolds shear stress is compared with the incompressible DNS of Schlatter and Orl¨ u in figure 6 for Re = 3270 and 3630. Both DNS profiles show self-similarity and are consistent with the θ expected behavior of Reynolds shear stress for a fully developed equilibrium boundary layer. The predicted profiles in the present simulation for both Re show a trend where Reynolds shear stress is underpredicted θ when normalized by friction velocity ( u ) and overpredicted when the normalized Reynolds shear stress is τ scaled by local density (¯ ρ/ ¯ ρ ). Also, the Reynolds shear stress at Re = 3270 peaks above 1 . 0 in the log w θ layer. This is a concern for two potential reasons: 1) it is a sign of inadequate resolution at the boundary layer edge and 2) the boundary layer has not reached an equilibrium state at that station. This was also 39 40 observed by Gatski and Erlebacher and Maeder et al., although their work was inconclusive on the 18 41 cause of this behavior. However, Guarini et al. and Pirozzoli et al. later showed with their highly resolved Mach 2.25 and 2.5 DNS simulations, respectively, that the compressible flow simulations should follow established incompressible flow data as proposed by Morkovin, albeit density scaling was necessary as shown by Guarini et al. It is important to note that the predicted Reynolds shear stress at Re = 3630 θ matches better with the DNS suggesting that the boundary layer is approaching an equilibrium state albeit + + at a longer axial distance. Moreover, it was found that in the viscous sublayer − uv = 0 . 0008 y for both Re discussed above, which was also shown by Patel et al., where the authors made a thorough comparison θ of various turbulence models with experimental data.

11 of 20 American Institute of Aeronautics and Astronautics (a) Skin friction coefficient (b) van Driest transformed velocity profiles Figure 5. Centerline and span-averaged characteristics of the upstream boundary layer Figure 6. Span-averaged Reynolds shear stress at Re = 3270 (red) and 3630 (blue) θ The turbulent kinetic energy and Reynolds normal stresses show a similar behavior as the Reynolds ′ ′ ′ ′ ¨ shear stress. Compared to the incompressible DNS of Schlatter and Orl¨ u, figure 7(a), v v and w w are + + ′ ′ ′ ′ ′ ′ underpredicted in the range of 10 . y . 100, while u u was overpredicted. Beyond y ∼ 100, v v and w w ′ ′ are closer to the DNS, yet underpredicted and u u , contrary to the buffer layer region, was underpredicted.

41 18 As discussed before, Pirozzoli et al. and Guarini et al. have shown that scaling the Reynolds stress by local density improves the comparison of compressible predictions with incompressible DNS. However, ′ ′ ′ ′ as shown in figure 7(b), improvement was only realized in the v v and w w components for the buffer ′ ′ layer and beginning of the log layer regions. Density scaling does however show that u u was consistently overpredicted rather than what was observed for the unscaled comparison in figure 7(a).

Finally, the budget of turbulent kinetic energy is presented in figure 8 to validate the profiles forward of the SBLI region. The budget terms are normalized by ¯ ρ u /ν . The production, dissipation, molecular w w τ ¨ diffusion, and turbulent transport are compared with Schlatter and Orl¨ u and presented at a matching Re .

θ Although, the present simulation follows the trend of the incompressible DNS, the magnitude of the terms 12 of 20 American Institute of Aeronautics and Astronautics 2 2 (a) Normalized by u (b) Normalized by u and scaled by ¯ ρ/ ¯ ρ w τ τ Figure 7. Span-averaged turbulent kinetic energy and Reynolds normal stress in the budget are larger. This can be attributed to the lack of mesh resolution at the simulated Reynolds number (cf. table 2), suggesting that a finer mesh or lowering the simulated Reynolds number may improve comparison with the DNS.

Figure 8. Span-averaged budget of the turbulent kinetic energy at Re = 3270 θ IV.B. Budgets for Shock Wave/Boundary-layer Interaction The analysis of the shock wave/boundary-layer interaction (SBLI) will focus on budgets of the turbulent kinetic energy at key locations in the SBLI region. Comparisons of the mean flow and separation charac- teristics with the experimental measurements are omitted here and readers are directed to the work of Mullenix and Gaitonde, where successes and shortfalls of the predictions using the present methodology were discussed.

Figure 9 shows contours of ¯ u in the vicinity of the SBLI at the spanwise centerline and identifies key locations, listed in table 4, which will be used during the discussion of the results. The dashed lines mark 13 of 20 American Institute of Aeronautics and Astronautics theoretical extension, to the wall, of the mean impinging and the reflected shocks. Under the confluence of these shocks, a separation bubble is evident in the interaction region. A comparison of the centerline and span-averaged skin friction, presented in figure 10, show that between x = 58 . 2 and x = 60 . 9, the skin friction becomes negative—an indication of flow separation.

Figure 9. Contours of ¯ u at the domain centerline Figure 10. Skin friction in the SBLI region Table 4 shows a list of key locations for which the budget will be presented and the corresponding Re θ where applicable. The positions are relative to the mean impinging and reflected shocks and defined in figure 9.

The budgets are presented to indicate the behavior of the terms as the flow approaches the interaction region, within the interaction, and after. In all plots, production, dissipation, molecular diffusion, and turbulent transport terms are shown. As before, the budget terms are normalized by ¯ ρ u /ν . The pressure w w τ diffusion, pressure dilatation, and mass flux terms were small such that they would appear as zero on each plot, thus they are omitted.

14 of 20 American Institute of Aeronautics and Astronautics Table 4. Key locations Station Location x Re θ 1 Incoming flat plate boundary layer 56.0 3270 2 Upstream of the reflected shock 56.6 3280 3 Downstream of the reflected shock 57.5 3300 4 Separation bubble 59.5 - 5 Downstream of the impinging shock 63.4 - 6 Recovered flat plate boundary layer 70.0 - (a) Upstream, x = 56 . 6 and Re = 3280 (b) Downstream, x = 57 . 5 and Re = 3300 θ θ Figure 11. Span-averaged budget upstream and downstream of the mean reflected shock (a) Separation bubble, x = 59 . 5 (b) Relaxation region, x = 63 . 4 Figure 12. Span-averaged budget upstream and downstream of the mean impinging shock The effect of the interaction region was observed as early as station x = 56 . 6 (figure 11(a)), where a marginal increase in the budget terms can be observed (cf. with figure 8). The increase in magnitude for all four quantities was evident at station x = 57 . 5, which was just aft of the extension, to the wall, of the 15 of 20 American Institute of Aeronautics and Astronautics mean reflected shock position. Particularly, the production and dissipation terms increased by an order of magnitude. The turbulent transport term (cf. figures 11(a) and 11(b)) shows that the trough in the + buffer layer ( y = 14 . 5) has moved closer to the wall and a new peak has developed at the beginning of the log layer. This fundamentally alters the near wall dynamics responsible for transfer and propagation of the turbulent kinetic energy. Interestingly, at this station, peaks and troughs for production, molecular diffusion, and dissipation terms remain at approximately the same location in the wall normal direction as the station forward of the reflected shock.

Figure 12 shows two locations, station x = 59 . 5 within the separation bubble and x = 63 . 4 past the theoretical extension, to the wall, of the mean impinging shock. Note that the abscissa has been extended + to y = 600 to accommodate the active region of the boundary layer. In the separation bubble, figure 12(a), the peak in the production term moved away from the wall, out of the buffer layer, and into the + beginning of the log layer. The shape of the peak broadened and spanned up to y = 200, well into the log + layer. The turbulent transport term also shifted away from the wall with a peak at y = 170, showing a coupling between the production and the turbulent transport terms. Also the peak in the sublayer has grown in magnitude for the turbulent transport term, a sign of activity in the near wall region. The molecular diffusion and dissipation terms showed no change in dynamics, but did show a change in magnitude. In the relaxation region, figure 12(b), a lower peak in the production term, reminiscent of the boundary layer + before the reflected shock has developed at y = 6 . 5, an indication of the flow returning back to the + undisturbed boundary layer state. Also, a larger peak at y = 150 was also observed which corresponds to the high production region evident in figure 15. The turbulent transport term behavior also showed near wall behavior similar to the undisturbed boundary layer, a peak in the sublayer and a trough in the buffer + layer. A trough near y = 150 seems to coincide with a peak in the production at the same location.

Further downstream, the boundary layer begins to approach an equilibrium form and the budgets, pre- sented in figure 13, showed identical near wall behavior to that of station x = 56 . 0, but with marginally higher magnitudes. Although not shown here, the production and turbulent transport terms remained active 2 + 3 in the log and outer layer regions (10 < y < 10 ) with secondary peaks and troughs.

Figure 13. Span-averaged budget of the turbulent kinetic energy at x = 70 . 0 The ratio of production-to-dissipation of the turbulent kinetic energy, P / , can be often found in the 43, 44 literature pertaining to turbulence modeling. In wall-bounded flows, it can be used to confirm a local equilibrium in the log layer of a turbulent boundary layer. In figure 14, P / is plotted for the incompressible ¨ DNS of Schlatter and Orl¨ u along with various stations of the present simulation. The DNS data shows that in the log layer the ratio approach the value of unity as expected, an indication of local equilibrium.

A comparison of the incoming turbulent boundary layer at station x = 56 . 0 ( Re = 3270), upstream of the θ SBLI region, shows that the present simulation does not attain a complete equilibrium form in the log layer.

Nonetheless, the profiles at two upstream stations ( x = 56 . 0 and x = 56 . 6) and one downstream station 16 of 20 American Institute of Aeronautics and Astronautics ( x = 70 . 0) matched reasonably well with the DNS in the viscous sublayer and exhibits self-similarity. In the SBLI region, stations x = 57 . 5 and x = 59 . 5, a dramatic change in the ratio is evident and the local flowfield is far from equilibrium as expected.

Figure 14. The production-to-dissipation ratio, P / Figure 15. Iso-surface of the production term at 0 . 5 , 2 . 0 , and 4 . 0 and colored by dissipation 17 of 20 American Institute of Aeronautics and Astronautics

Iso-surfaces of the production term are presented in figure 15 at three levels, 0.5, 2.0, and 4.0. They are colored by dissipation to show regions where dissipation may balance production of the turbulent kinetic energy for each iso-surface. The values are normalized by ¯ ρ u /ν . The amplification of the production term w w τ as flow approaches the SBLI is clear. The large production rate around the mean impinging and reflected shocks was local and the production drops quickly away from the shocks. The region between the two shocks show large magnitude of production, which is not balanced by dissipation. The production in the relaxation region continues at elevated levels and subsides gradually downstream.

V. Conclusion

The budget for turbulent kinetic energy has been presented forward, through, and aft of an SBLI.

Validation of the ILES approach by means of comparison with the DNS showed that the fine mesh considered in the current simulations was inadequate to accurately resolve the turbulent scales in the buffer layer and log layer even at a reduced Reynolds number. However, the current simulation predicted matching trends with that of the DNS when scaled by local density.

Comparison with the DNS also suggested that the boundary layer may not have reached a complete equilibrium forward of the SBLI, thus a longer axial domain for the boundary layer development, aft of the counterflow force actuator, is necessary. Alternatively, newer methods of specifying unsteady inflow conditions are being investigated.

The budgets presented for the SBLI provide insight into the dynamic nature of the interaction region, i.e., the quickly changing nature of the production and the turbulent transport terms and more subdued nature of the dissipation and the viscous diffusion terms. In particular, the turbulent transport term showed shifting of peaks and troughs in the wall normal direction. A new peak also developed immediately aft of the reflected shock suggesting a jump-start in the transport process to compensate for the increase in the production term. Large magnitude of the budget terms persisted in the interaction region, but began to subside in the relaxation region. At the farthest downstream location the budget looked near identical to that of the incoming boundary layer upsteam of the SBLI. Overall, considering the difference in mesh resolution between the present ILES versus the DNS, the ILES predictions captured the trends reasonably well.

Acknowledgments

The authors would like to thank our colleague, Dr. Dennis Yoder, for countless discussions regarding turbulence modeling, and NASA’s Transformative Aeronautics Concepts Program and Transformational Tools and Technology Project for its generous support.

References

Dolling, D. S., “Fifty Years of Shock-Wave/Boundary-Layer Interaction Research: What Next?” AIAA Journal , Vol. 39, No. 8, August 2001, pp. 1517–1531.

Zheltovodov, A. A., “Some Advances in Research of Shock Wave/Turbulent Boundary-layer Interactions,” No. 496 in 44th Aerospace Sciences Meeting, American Institute of Aeronautics and Astronautics, January 2006.

Gaitonde, D. V., “Progress in Shock Wave/Boundary Layer Interactions,” No. 2607 in 43rd Fluid Dynamics Conference, American Institute of Aeronautics and Astronautics, June 2013.

Souverein, L. J., Dupont, P., Debi` eve, J.-F., Dussauge, J.-P., van Oudheusden, B. W., and Scarano, F., “Effect of Interaction Strength on Unsteadiness in Turbulent Shock Wave Induced Separations,” AIAA Journal , Vol. 48, No. 7, July 2010, pp. 1480–1493.

Dussauge, J.-P., Dupont, P., and Debi` eve, J.-F., “Unsteadiness in Shock Wave Boundary Layer Interactions with Sepa- rations,” Aerospace Science and Technology , Vol. 10, No. 2, December 2006, pp. 85–91.

Ganapathisubramani, B., Clemens, N. T., and Dolling, D. S., “Effect of Upstream Boundary Layer on the Unsteadiness of Shock-induced Separation,” Journal of Fluid Mechanics , Vol. 585, April 2007, pp. 369–394.

Dupont, P., Haddad, C., and Debi` eve, J.-F., “Space and Time Organization in a Shock-induced Separated Boundary Layer,” Journal of Fluid Mechanics , Vol. 559, December 2006, pp. 255–277.

Li, J., Priebe, S., Grube, N., and Mart´ ın, M. P., “Conditional Analysis of the Unsteadiness in Shock Wave/Turbulent Boundary-layer Interactions,” No. 0436 in 52nd Aerospace Sciences Meeting, American Institute of Aeronautics and Astronau- tics, January 2014.

Clemens, N. T. and Narayanaswamy, V., “Low-frequency Unsteadiness of Shock Wave/Turbulent Boundary Layer Inter- actions,” The Annual Review of Fluid Mechanics , Vol. 46, January 2014, pp. 469–492.

18 of 20 American Institute of Aeronautics and Astronautics DeBonis, J. R., Oberkampf, W. L., Wolf, R. T., Orkwis, P. D., Turner, M. G., Babinsky, H., and Benek, J. A., “Assessment of Computational Fluid Dyanmics and Experimental Data for Shock Boundary-layer Interactions,” AIAA Journal , Vol. 50, No. 4, April 2012, pp. 891–903.

Georgiadis, N. J., Rizzetta, D. P., and Fureby, C., “Large-eddy Simulations: Current Capabilities, Recommended Prac- tices, and Future Research,” AIAA Journal , Vol. 48, No. 8, August 2010, pp. 1772–1784.

Gaitonde, D. V. and Visbal, M. R., “Pad´ e-type Higher-order Boundary Filters for the Navier-Stokes Equations,” AIAA Journal , Vol. 38, No. 11, November 2000, pp. 2103–2112.

Visbal, M. R., Morgan, P. E., and Rizzetta, D. P., “An Implicit LES Approach Based on High-order Compact Differencing and Filterning Schemes,” No. 4098 in 16th Computational Fluid Dynamics Conference, American Institute of Aeronautics and Astronautics, June 2003.

Menter, F. R., “Zonal Two Equation k − ω Turbulence Models for Aerodynamic Flows,” No. 2906 in 24th Fluid Dynamics Conference, American Institute of Aeronautics and Astronautics, July 1993.

Wilcox, D. C., Turbulence Modeling for CFD , DCW Industries, Inc., 3rd ed., July 2010.

Hamlington, P. E. and Dahm, W. J. A., “Reynolds Stress Closure for Nonequilibrium Effects in Turbulent Flows,” Physics of Fluids , Vol. 20, No. 11, November 2008.

Sinha, K., Mahesh, K., and Candler, G. V., “Modeling Shock Unsteadiness in Shock/Turbulence Interaction,” Physics of Fluids , Vol. 15, No. 8, August 2003, pp. 2290–2297.

Guarini, S. E., Moser, R. D., Shariff, K., and Wray, A., “Direct Numerical Simulation of a Supersonic Turbulent Boundary Layer at Mach 2.5,” Journal of Fluid Mechanics , Vol. 414, January 2000, pp. 1–33.

Huang, P. G., Coleman, G. N., and Bradshaw, P., “Compressible Turbulent Channel Flows: DNS Results and Modelling,” Journal of Fluid Mechanics , Vol. 305, August 1995, pp. 185–218.

Gaitonde, D. V. and Visbal, M. R., “High-order Schemes for Navier-Stokes Equations: Algorithm and Implementation into FDL3DI,” Technical Report 3060, Air Force Research Laboratory, August 1998.

Anderson, D. A., Tannehill, J. C., and Pletcher, R. H., Computational Fluid Mechanics and Heat Transfer , Computational Methods in Mechanics and Thermal Sciences, Hemisphere Publishing Corporation, 1984.

Shyy, W., Jayaraman, B., and Andersson, A., “Modeling of Glow Discharge-induced Fluid Dynamics,” Journal of Applied Physics , Vol. 92, No. 11, December 2002, pp. 6434–6443.

Gaitonde, D. V., Visbal, M. R., and Roy, S., “Control of Flow Past a Wing Section with Plasma-based Body Forces,” No. 5302 in 36th Plasmadynamics and Lasers Conference, American Institute of Aeronautics and Astronautics, June 2005.

Mullenix, N. J., Gaitonde, D. V., and Visbal, M. R., “Spatially Developing Supersonic Turbulent Boundary Layer with a Body-Force-Based Method,” AIAA Journal , Vol. 51, No. 8, August 2013, pp. 1805–1819.

Mullenix, N. J. and Gaitonde, D. V., “Analysis of Unsteady Behavior in Shock/Turbulent Boundary Layer Interac- tions with Large-Eddy Simulations,” No. 0404 in 51st Aerospace Sciences Meeting, American Institute of Aeronautics and Astronautics, January 2013.

Waindim, M., Yentsch, R. J., and Gaitonde, D. V., “A Body-force Based Method to Generate Supersonic Equilibrium Turbulent Boundary Layer Profiles,” No. 0940 in 52nd Aerospace Sciences Meeting, American Institute of Aeronautics and Astronautics, January 2014.

Waindim, M. and Gaitonde, D. V., “Results and Analysis of Implicit Large Eddy Simulations of Equilibrium Spatially Developing Turbulent Boundary Layers at Multiple Mach Numbers,” No. FEDSM2014-21391 in Fluids Engineering Division Summer Meeting, American Society of Mechanical Engineers, August 2014.

Mullenix, N. J. and Gaitonde, D. V., “A Bandwidth and Order Optimized WENO Interpolation Scheme for Compressible Turbulent Flows,” No. 0366 in 49th Aerospace Sciences Meeting, American Institute of Aeronautics and Astronautics, January 2011.

Visbal, M. R. and Gaitonde, D. V., “On the Use of High-order Finite-difference Schemes on Curvilinear and Deforming Meshes,” Journal of Computational Physics , Vol. 181, May 2002, pp. 155–185.

Webb, N., Clifford, C., and Samimy, M., “Preliminary Results on Shock Wave/Boundary-layer Interaction Control Using Localized Arc Filament Plasma Actuators,” No. 3426 in 41st Fluid Dynamics Conference, American Institute of Aeronautics and Astronautics, June 2011.

Touber, E. and Sandham, N. D., “Comparison of three Large-eddy Simulations of Shock-induced Turbulent Separtion Bubbles,” Shock Waves , Vol. 19, August 2009, pp. 469–478.

Visbal, M. R., Rizzetta, D. P., and Mathew, J., “Large-eddy Simulations of Flow Past a 3-D Bump,” No. 0917 in 45th Aerospace Sciences Meeting, American Institute of Aeronautics and Astronautics, January 2007.

Jeong, J. and Hussain, F., “On the Identification of a Vortex,” Journal of Fluid Mechanics , Vol. 285, July 1995, pp. 69–94.

Bisek, N. J., “High-order Implicit Large-eddy Simulations of a Supersonic Corner Flow,” No. 0588 in 52nd Aerospace Sciences Meeting, American Institute of Aeronautics and Astronautics, January 2014.

White, F. M. and Christoph, G. H., “A Simple Theory of the Two-dimensional Compressible Turbulent Boundary Layer,” Journal of Basic Engineering , Vol. 94, No. 3, September 1972, pp. 636–642.

White, F. M., Viscous Fluid Flows , Series in Mechanical Engineering, McGraw-Hill, 2nd ed., 1991.

van Driest, E. R., “On Turbulent Flow Near a Wall,” Journal of Aeronautical Sciences , Vol. 23, No. 11, January 1956, pp. 1007–1011.

¨ Schlatter, P. and Orl¨ u, R., “Assessment of Direct Numerical Simulation Data of Turbulent Boundary Layers,” Journal of Fluid Mechanics , Vol. 659, July 2010, pp. 116–126.

Gatski, T. B. and Erlebacher, G., “Numerical Simulation of a Spatially Evolving Supersonic Turbulent Boundary Layer,” Technical Memorandum 2002-211934, NASA Langley Research Center, September 2002.

Maeder, T., Adams, N. A., and Kleiser, L., “Direct Simulation of Turbulent Supersonic Boundary Layers by an Extended Temporal Approach,” Journal of Fluid Mechanics , Vol. 429, August 2001, pp. 187–216.

19 of 20 American Institute of Aeronautics and Astronautics Pirozzoli, S., Grasso, F., and Gatski, T. B., “Direct Numerical Simulation and Analysis of a Spatially Evolving Supersonic Turbulent Boundary Layer at M=2.25,” Physics of Fluids , Vol. 16, No. 3, March 2004, pp. 530–545.

Morkovin, M. V., “Effects of Compressibility on Turbulent Flows,” The Mechanics of Turbulence , edited by A. Favre, Gordon and Breach Science Publishers, New York, 1964, pp. 367–380.

Patel, V. C., Rodi, W., and Scheuerer, G., “Turbulence Models for Near-wall and Low Reynolds Number Flows: A Review,” AIAA Journal , Vol. 23, No. 9, January 1985, pp. 1308–1319.

Mansour, N. N., Kim, J., and Moin, P., “Near-wall k −  Turbulence Modeling,” AIAA Journal , Vol. 27, No. 8, 1989, pp. 1068–1073.

20 of 20 American Institute of Aeronautics and Astronautics