Reynolds-Stress Budgets in an Impinging Shock Wave/Boundary-Layer Interaction

ρ density τ Reynolds-stress tensor ij θ boundary-layer momentum thickness ζ second viscosity B law of the wall intercept constant C specific heat at constant pressure p e internal energy e total energy I , I , I digital filter length scale in x , y , and z directions, respectively x y z k turbulent kinetic energy L length scale M Mach number N , N digital filter size in y and z directions y z p pressure P r Prandtl number Re Reynolds number s span S mean strain-rate tensor ij s instantaneous strain-rate tensor ij T temperature t time t , t viscous stress tensor ij u , v , w non-dimensional velocity in x , y , and 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, or 3 imp impingement location int interaction inv inviscid ref , ∞ reference value same as freestream value sep separation vis viscous w wall Conventions ¯ Reynolds averaged ˜ Favre averaged Superscripts ′ Reynolds decomposition ′′ Favre decomposition ∗ non-dimensional value T transpose tot total tr trace I. Introduction Literature is replete with experimental, numerical, and modeling studies of shock wave/boundary-layer 1–3 interactions (SBLIs). This interest is warranted because SBLIs are ubiquitous in transonic, supersonic, and hypersonic speed regimes and thus relevant to commercial, military, and space vehicles of the past, present, and future. SBLIs are known to create excessive unsteady aerothermal loads that can compromise structural integrity, cause component failure, and result in loss of control. In particular and of relevance to the 2 of 25 American Institute of Aeronautics and Astronautics present work is the adverse pressure gradient in the propulsion flowpath, which can cause flow separation, distortion, and loss in engine efficiency. Because of this, SBLIs are one of the most actively researched phenomena in high speed flows, yet far from being completely understood, and hence numerical simulation and modeling of SBLIs are a continuing challenge.

The unsteady nature of SBLIs and the resulting three-dimensional (3D) separation, often including corner separation (i.e., propulsion flowpaths are often rectangular), make their prediction extremely difficult.

Investigating the source of low-frequency unsteadiness of the reflected (or separation) shock in an SBLI, which is typically one or two orders of magnitude lower than the frequencies within the incoming turbulent 1, 4, 5 boundary layer, is important for flow control strategies. However, understanding the locally anisotropic and inhomogeneous turbulence field is key to accurate prediction of separation via numerical simulations and modeling of SBLIs.

While the corner influence is not inherent to the incident/reflected shock configuration, it is present in practical situations where sidewalls exist, e.g. wind tunnels, inlets, isolators, etc. Also, it is important to address the fact that the majority of experimental configurations involving compression ramp and inci- dent/reflected shock exhibit some level of corner influence due to corner separation, which renders portions of the flow near corners 3D, while the centerline remains two-dimensional (2D). Due to such a prevalence and coupling, corner influence has become an area of active research as shown by the works of Babinsky et 6 7, 8 9 al., Benek et al., and Eagle at al., which includes experimental as well as numerical investigations.

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, where an incident oblique shock impinges on a turbulent boundary layer. Depending on the strength of the incident shock, the flow may not separate, incipiently separate, or completely separate. In figure 1, the flow is shown to separate. The streamwise length of the separation bubble is L . The interaction length obtained by extending the inviscid portions of the incident and sep reflected shocks to the wall is L . The development of the separation results in the forward migration of int the reflected shock, of course this is in a steady sense. In reality, the separation bubble has been shown to shrink and grow in size, which causes the reflected shock to oscillate aft (small separation) and forward (large separation) with a characteristic low frequency. Numerous experimental and computational works have shown this to be between Strouhal number 0 . 02 − 0 . 05.

Figure 1. Two-dimensional anatomy of an impinging shock wave/boundary-layer interaction.

3 of 25 American Institute of Aeronautics and Astronautics 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.

This is mainly due to their ease of use, and also because of 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. However, scale-resolving methods offer significant gains in accuracy.

A list of scale-resolving approaches is compiled 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 Pad´ e-type filter to dissipate energy in the high spatial wavenumber 12, 13 range. 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 in 2011 by the American Institute of Aeronautics and Astronautics (AIAA) with an intent to share, assess, and determine the most promising SBLI prediction methods. 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 and 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 relative accuracy of a method 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 into 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. 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 a nonlinear algebraic equation or by using a Reynolds-stress transport model. However, among these approaches, none has emerged as clearly superior.

In the authors’ previous work, terms in the exact equation of turbulent kinetic energy were studied for a developing turbulent boundary layer and an SBLI. However, due to increased interest in Reynolds- stress models, here the focus is shifted on the exact equation of Reynolds stress transport. Understanding of the various terms in the exact equation 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 turbulent boundary 18 19 20 layers, channel flows, and square ducts. The present work is the first comprehensive examination of the Reynolds-stress transport terms within an SBLI flowfield. Our objective is twofold: 4 of 25 American Institute of Aeronautics and Astronautics

• Validate the simulation with an SBLI experimental database to ensure consistency of the computed quantities.

• Study the budgets of Reynolds stress and other relevant turbulence quantities with an intent to inform future model development.

II. Governing Equations In this section the Reynolds stress transport equation will be discussed along with the equations pertinent to the ILES formulation and the digital filter approach, which is used to obtain turbulent fluctuations at the inflow.

A. Reynolds-stress Budget The Reynolds stress is defined as ′′ ′′ ¯ ρτ = − ρu u (1) ij i j The Reynolds stress transport 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.

∂ ∂ ′′ ′′ ′′ ′′ ( ρu u ) + ( ρu u ˜ u ) = P + D − ¯ ρ + Π + M (2) k ij ij ij ij ij i j i j ∂t ∂x k Each term on the right-hand side of equation 2 is defined as ∂ ˜ u ∂ ˜ u j i ′′ ′′ ′′ ′′ P = − ρu u − ρu u Production (3) ij i k j k ∂x ∂x k k ν T P D = D + D + D Diffusion (4) ij ij ij ij [ ] ∂ ν ′′ ′′ D = u t + u t Molecular Diffusion (5) kj ki ij i j ∂x k [ ] ∂ T ′′ ′′ ′′ D = − ρu u u Turbulent Diffusion (6) ij i j k ∂x k [ ] ∂ P ′′ ′′ ′ ′ D = − p u δ + p u δ Pressure Diffusion (7) jk ik ij i j ∂x k ( ) ′′ ′′ ∂u ∂u j i ′ Π = p + Pressure Strain (8) ij ∂x ∂x j i ′′ ′′ ∂u ∂u j i ¯ ρ = t + t Dissipation (9) ij ki kj ∂x ∂x k k ( ) ( ) ¯ ¯ ∂ t ∂ ¯ p ∂ t ∂ ¯ p kj ki ′′ ′′ M = u − + u − Mass Flux (10) ij i j ∂x ∂x ∂x ∂x k j k i 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 + ζ δ (11) ij ij ij ∂x k In equation 11, ζ is obtained by relating it to μ . Such an assumption is valid for monatomic gases and widely used in computational fluid dynamics.

ζ = − μ (12) 5 of 25 American Institute of Aeronautics and Astronautics B. Compressible Navier-Stokes Equations The compressible Navier-Stokes equations in non-dimensional form are given by ∂ρ ∂ + ( ρu ) = 0 (13) i ∂t ∂x i ∂ ∂ 1 ( ρu ) + ( ρu u + pδ − t ) = 0 (14) i i j ij ij ∂t ∂x Re j [ ] ∂ ∂ 1 1 ( ρe ) + u ( ρe + p ) − ( u t ) + q = 0 (15) 0 i 0 i ij i ∂t ∂x Re ( γ − 1) P rM Re i The total energy and heat flux are defined as 1 ∂T e = e + u u q = − μ (16) 0 i i i 2 ∂x i The non-dimensionalization is performed using the following definitions where * (asterisk) represents non-dimensional quantities. Except equations 17a and 17b, the * has been dropped for simplicity and all quantities are non-dimensional in this paper, unless stated otherwise.

x u tV ρ i i ref ∗ ∗ ∗ ∗ x = u = t = ρ = (17a) i i L V L ρ ref ref p T μ e ∗ ∗ ∗ ∗ p = T = μ = e = (17b) 2 2 ρ V T μ V ref ref ref ref ref The reference conditions are the upstream 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 p ref ref using Sutherland’s law and a perfect gas is assumed.

ρ V L μ C V ref ref ref p ref √ Re = P r = M = (18) γ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 (19) ∂t ∂x ∂y ∂z ∂x ∂y ∂z T where U is the solution vector defined as U = [ ρ, ρu, ρv, ρw, ρe ] . The inviscid and viscous flux vectors are defined as     ρu 0         ρu + p t xx         F = F = t (20)  ρuv    inv vis xy   Re       t ρuw xz     μ ∂T ( ut + vt + wt ) + ( ρe + p ) u xx xy xz 2 ( γ − 1) P rM ∂x The above equations are transformed into curvilinear coordinates, and using the strong conservation form 21, 22 the following is obtained 6 of 25 American Institute of Aeronautics and Astronautics ˆ ˆ ˆ ˆ ˆ ˆ ˆ ∂ U ∂ F ∂ G ∂ H ∂ F ∂ G ∂ H inv inv inv vis vis vis ˆ + + + = + + + S (21) ∂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 ) (22) inv x inv y inv z inv J ˆ F = ( ξ F + ξ G + ξ H ) (23) vis x vis y vis z vis J ˆ ˆ ˆ ˆ and similarly G , H , G , and H .

inv inv vis vis C. Unsteady Inflow Boundary Method In a previous study, budget of turbulent kinetic energy was investigated by Vyas et al. using a similar approach to that taken here. However, in that work, the generation of the supersonic turbulent boundary 23–25 layer was achieved by using the counterflow force model. The approach, although functional, required a large streamwise domain to facilitate a transition from the specified laminar boundary-layer inflow profile to a turbulent boundary layer. Thus, an alternate approach where such laminar-to-turbulent transition can be avoided by specifying an unsteady boundary-layer profile at the inflow is highly desired. This will also allow for smaller mesh sizes or possibly transferring some of the mesh points into the area of interest, i.e., the shock wave/boundary-layer interaction region.

1. The Digital Filter Approach The digital filter approach was originally proposed by Klein et al., where the filtering operation was performed in 3D space. Later, Veloudis et al. investigated specification of varying filter coefficients in the wall-normal direction to allow for varying length scales. Xie and Castro simplified and sped up the approach by performing the filtering operation on the 2D inflow plane and correlating the calculated 2D field with the one at the previous timestep to account for the length scale in streamwise direction using Taylor’s 29 30 hypothesis. The approach was further improved by Touber and Sandham and the current work closely follows their implementation of the digital filter.

Assume a set of p random numbers, { r } , such that the set has zero mean and unit variance.

k 1 ≤ k ≤ p p p ∑ ∑ r = r /p = 0 r r = r /p = 1 (24) k k k k k k =1 k =1 A filter operator can be defined, where N is a positive integer related to the filter length scale and { b } , a set of real numbers, are the filter coefficients.

j − N ≤ j ≤ N N ∑ v ≡ F ( r ) = b r (25) k N k j k + j j = − N The filter operator, equation 25, is linear and non-recursive. Thus the averaging and filtering operations commute, and knowing the properties of { r } in equation 24, the following is true.

k N ∑ r = 0 r r = b b (26) k k k + q j j − q j = − N + q 28 30 Following the approach of Xie and Castro and Touber and Sandham, the correlation function was chosen to be exponential as opposed to Gaussian, assumed by Klein et al. in the original work. Xie and Castro argued that the choice of exponential form will produce an energy-decay rate in the inertial subrange which corresponds to slope of − 2, rather than − 5 / 3. However, the behavior of large-scale structures is realistically modelled. Such an assumption is valid since the digital filter is operational at the inflow plane 7 of 25 American Institute of Aeronautics and Astronautics and the domain has enough length upstream of the region of interest to allow the flow to recover the correct velocity spectra in the inertial subrange. The two-point correlation at the y − z inflow plane is defined as: ( ) πr R ( r, α, t ) = exp − (27) 2 I In equation 27, r is the separation variable, α = y = z = 0 , t = 0 is the time, and I is the length scale defined such that I = n ∆ α . So equation 28 follows such that q ∆ α represents separation distance and q α I α is the interval step.

( ) v v π | q | k k + q R ( q ∆ α ) = = exp − (28) vv v v 2 n k k Using equations 26 and 28, filter coefficients can be computed by solving: N ∑ b b j j − q ( ) π | q | j = − N + q = exp − (29) N ∑ 2 n b j j = − N Xie and Castro obtained a solution for filter coefficient b as: k ( ) ˜ b π | k | k ˜ b = , where b ' exp − (30) k k   (1 / 2) n N ∑ ˜   b j j = − N In the above, N is defined as N ≥ 2 n . Thus, the above one-dimensional approach can be extended to α α calculate a 2D filter coefficient at the y − z inflow plane as: b = b b (31) jk j k Thus, an approximation for a 2D filter coefficient becomes: [ ( )] | j | | k | exp − π + n n y z b = (32) jk   1 / 2 [ ( )] N N y  z  ∑ ∑ | j | | k | exp − 2 π +  n n  y z k = − N j = − N z y In general, the following steps are taken within the digital filter algorithm: 1. Generate a set of p random numbers with zero mean and unit variance corresponding to the 2D inflow mesh. Here, the Box-Muller theorem is used to combine two independent and uniformly distributed sets of numbers a and b in the range (0 , 1) into c and d , which are independent and normally distributed, and also satisfy zero mean and unit variance conditions. They can be calculated as below, however, either one can be used.

√ √ c = − 2 ln( a ) cos(2 πb ) d = − 2 ln( a ) sin(2 πb ) (33) 2. Select the relevant length scales, I , I , and I . Since the current implementation is 2D, I and I are x y z y z converted into an equivalent number of grid points using the grid spacing, i.e., n = I / ∆ α .

I α α 8 of 25 American Institute of Aeronautics and Astronautics

3. Now the filter coefficients can be calculated using equation 32. Filter bounds are calculated using N = 2 n .

α I α 4. The filter coefficients calculated in equation 32 can now be used to filter the zero mean and unit variance normally distributed random numbers calculated in equation 33. Thus, we impose relevant y − z length scales on the field of these normally distributed random numbers using equation 25.

5. The length scale in the x direction is imposed by correlating the old time step ( t − ∆ t ) to the current time step ( t ) as demonstrated by Xie and Castro to make filtering a 2D process in contrast to the 3D filtering originally proposed by Klein et al. In equation 34, ∆ t is the time step and t is the L Lagrangian time scale, such that t = I / ¯ u .

L x √ ( ) ( ) π ∆ t π ∆ t t − ∆ t t t v = v exp − + v 1 − exp − (34) k k k 2 t t L L 6. The final step involves the transformation originally proposed by Lund et al. to obtain a time- dependent inflow velocity field. Components of a are defined by the Reynolds-stress tensor, τ , ij ij obtained from a previous RANS turbulent flat plate calculation performed on the same mesh as the LES simulation, at identical Mach and Reynolds numbers, using the k −  turbulence model as described in Gerolymos.

t u = u + a v (35) i i ij j   √ τ 0 0 √     τ /a τ − a 0 a = (36) 21 11 22 ij 21   √ 2 2 τ /a ( τ − a a ) /a τ − a − a 31 11 32 21 31 22 33 31 32 7. Lastly, the calculation of thermodynamic fluctuations was addressed by Touber and Sandham by invoking the Strong Reynolds Analogy (SRA). However, they noted that the validity of such an as- sumption is debatable, since it holds true in a weak sense as shown by Guarini et al. in a DNS simulation of Mach 2.5 supersonic turbulent boundary layer. Since the goal here is to provide an approximate first guess at the inflow plane, this approach was determined to be sufficient and is shown below: ′ ′ 2 T ¯ u u 2 2 = − ( γ − 1) M where M = (37) ¯ u T γRT ′ 8. Invoking SRA also means that the pressure fluctuations in the boundary layer are negligible, i.e., p = 0.

Thus, the following relation for the density fluctuations results: ′ ′ ρ T = − (38) ¯ ρ T This concludes the procedure that produced the time-dependent fluctuating inflow plane. The filter parameters are provided in table 1. The separation unit spacing (∆ α ) in the y direction was obtained by picking an average spacing in the log-law region of the boundary layer, whereas grid spacing was picked in the z direction.

III. Simulation Methodology A. Numerical Schemes Spatial derivatives are calculated using a sixth-order compact spectral-like finite-difference scheme of Lele and later adapted in the FDL3DI code by Visbal and Gaitonde. An eighth-order low-pass Pad´ e-type non-dispersive spatial filter is applied after each sub-iteration to maintain stability and dissipate energy 9 of 25 American Institute of Aeronautics and Astronautics Table 1. Digital filter parameters Parameter x y z a Length Scale, I 0.5 0.5 0.5 α b ∆ α - 0.01 0.0195 Equivalent Gridpoints, M = 2 I /∆ α - 50 25 α α a Non-dimensionalized by δ b α = x, y, z 12, 13 in the high spatial wavenumber range where the turbulent energy spectrum is poorly resolved. The filter coefficient, α , which determines sharpness of the filter cutoff was set to 0 . 45 based on Bisek and f 37 38 Garmann et al. The implicit time integration was performed with the second-order Beam-Warming scheme using two Newton-like sub-iterations and approximate factorization. The non-dimensional time step for this problem was 0.002.

The use of an explicit sub-grid scale model is completely avoided in the present simulations since the intent is to study the budget of the exact equation of the Reynolds-stress transport. Such an approach was taken by Morgan et al., where flow conditions similar to the present study were used for Re = 4800 θ and comparisons were made with the DNS of Pirozzoli and Bernardini which showed that at low to moderate Reynolds number, omitting an explicit sub-grid scale model is suitable. Kawai et al. showed that inclusion of an explicit sub-grid scale model in addition to the low-pass filtering, as done in the present study, introduced excessive numerical dissipation of the resolved turbulence, which may underpredict turbulence statistics sought here. A comparison of high-fidelity implicit and sub-grid scale model LES for airfoils at low Reynolds number performed by Garmann et al. showed no benefit in using an explicit sub-grid model.

Thus, high-fidelity ILES, i.e., combination of high-order spatial scheme and high-order filter, is an attractive approach for this work as it is paired with well-resolved meshes in the near-wall region.

B. Boundary Conditions and Mesh Parameters The flow conditions represent the experiments performed at the Institut Universitaire des Syst` emes Ther- miques Industriel (IUSTI) in Marseille, France. 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 several orders of magnitude larger than the fine mesh used here would be re- quired. Thus, the simulated Reynolds number was varied, one matching that of the experiment and two others obtained by reducing the experimental Reynolds number by half successively. Reducing the Reynolds number allows the present work to be practical and achievable by ILES. Such an approach was also taken 42 40 43 by Mullenix and Gaitonde, Pirozzoli and Bernardini, and Visbal et al. In order to have consistent flow relations, the pressure was reduced by a factor of half and one-quarter while maintaining the temperature, when the Reynolds number was reduced. This assures that the Mach number and the velocity scales remain the same between the experiment and the simulation.

At the inflow of the computational domain, the digital filter was used to generate turbulent fluctuations and provide a turbulent boundary layer. The interior of the flowfield was initialized with the same RANS solution as the one used to impose the inflow boundary condition using the digital filter. The wall was treated as no-slip and adiabatic. 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-Hugoniot relations were used to impose the oblique shock generated by a wedge in the experiment. Such a simulated shock was achieved by maintaining the farfield boundary condition for x < x at the pre-shock conditions identified shock in table 2, while x > x were set to the post-shock conditions obtained via oblique shock relations at shock Mach 2.29 and 8.0-degree angle of deflection.

In figure 2, the mesh schematic is presented where the coordinates are non-dimensionalized by the length scale δ . For the purpose of clarity, the coarse mesh is presented and every fourth point is shown. The mesh can be divided into two distinct sections in the streamwise direction: 1) the constant x -spacing section is where the SBLI occurs and 2) the coarsening to the outflow section. Both 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 10 of 25 American Institute of Aeronautics and Astronautics

Table 2. Flow conditions upstream of the interac- tion region at ( x − x ) /δ = − 5 . 7 imp 99 Property Experiment Simulation (nominal) M 2.29 2.29 ∞ U , m/s 545.0 545.0 ∞ P , Pa 50663 12600, 25300, 50663 ∞ T , K 300 300 δ , mm 10 . 0 11 . 0 θ , mm 0 . 87 0 . 87 ∗ δ , mm 3 . 0 3 . 1 Re 53420 14500, 29000, 58000 δ Re 4640 1150, 2300, 4600 θ Figure 2. Mesh topology for SBLI simulation. Coordinates non-dimensionalized by the experimental δ .

th th Every 8 point in the x direction and 16 point in the y direction is shown for clarity.

Table 3. Mesh parameters Coarse Medium Fine Domain size a x × y × z 30 × 4 × 5 30 × 4 × 5 30 × 4 × 5 Computational points N × N × N 1025 × 257 × 257 1537 × 257 × 257 2049 × 257 × 257 x y z 6 6 6 N 67 . 7 × 10 101 . 5 × 10 135 . 35 × 10 total Constant region N 945 1417 1889 x ∆ x 0.0212 0.0141 0.0105 At ( x − x ) /δ = − 5 . 7 imp 99 + ∆ x 15 10 7 N 187 187 187 y,bl − 3 − 3 − 3 ∆ y 1 × 10 1 × 10 1 × 10 w + ∆ y 0.68 0.68 0.68 ∆ z 0.0195 0.0195 0.0195 + ∆ z 13 13 13 a x, y, z are non-dimensionalized by δ parameters for coarse, medium, and fine meshes.

11 of 25 American Institute of Aeronautics and Astronautics IV. Results A. A Comparison with the Experiment Normalized velocity, shown for various mesh levels in figure 3(a), is in agreement with the viscous sublayer.

Coarse, medium, and fine meshes at Re = 4600 lined on top of each other and showed a slope that matched θ the expected logarithmic law (log law) of the wall ( κ = 0 . 41 and B = 5 . 5), however an offset was present.

Garnier et al. also reported such a behavior where LES simulations were performed at Re = 4600 and θ the offset was attributed to the underprediction of the friction velocity, which showed small improvement with mesh refinement. Thus, this may indicate a lack of grid resolution in the log-law region. This is further supported by the coarse mesh solution at Re = 2300, which is in agreement with the log law, DNS of θ 40 45 Pirozzoli and Bernardini, and LES of Eitel-Amor at Re = 4400. A comparison of skin friction profiles θ for the three mesh levels (figure 3(b)) show good agreement in the region upstream of the SBLI and within the SBLI, but some variation was present in the reattachment region. Thus, the current mesh resolution, although not as refined as the DNS presented in figure 3, is deemed suitable for the present simulation at experimental Reynolds number.

(a) van Driest transformed velocity profiles (b) Skin friction Figure 3. Upstream boundary layer at ( x − x ) /δ = − 5 . 7 .

imp 99 A comparison of interaction lengths, presented in table 4, show that the ILES and DNS simulations consistently underpredict the interaction region, including the present work. All simulations were done with periodic boundary conditions in the spanwise direction using an appropriate spanwise domain to avoid artificial amplification of the separation bubble. Thus, the rather good agreement obtained by Touber and 30 46 Sandham can be attributed to the choice of the subgrid-scale model (Mixed-Time-scale). A separate RANS study was performed to better understand the interaction behavior relative to the choice of domain in the spanwise direction, i.e., periodic or half-span with one sidewall. Results, although omitted here for brevity, showed that the periodic simulations underpredicted the interaction length. The tunnel where the experiment was performed had an aspect ratio of 1 . 39 and the inverse viscous aspect ratio, δ /s , of − 2 6 5 . 88 × 10 , thus characterized as quasi-two-dimensional, but it is suggested by the authors of the present work that the experimental data was still influenced by sidewalls. Such a relation was hypothesized by Babinsky et al. between the separation length and inverse viscous aspect ratio. At large and moderate inverse viscous aspect ratios, the corner separation and sidewall effects become critical in prediction of the interaction region at the tunnel centerline. Agostini et al. noted that the interaction region can be underpredicted by up to 20%. In the LES simulation of Larchevˆ eque et al., it was concluded that the ◦ ◦ ◦ interaction region produced by a 9 . 5 deflection case matched experiments with 8 and 8 . 8 deflections when performed by two different groups of researchers. Thus, only at sufficiently large span can such an experiment become two-dimensional and then a choice of periodic boundary condition would be accurate. So it can be 12 of 25 American Institute of Aeronautics and Astronautics argued that the present ILES calculation simulates SBLI in such a large span wind tunnel and captures key physics void of corner separation and sidewalls. Hence making it more attractive for studying turbulence modeling.

Table 4. Comparison of previous and present simulations with the experimental data.

a Authors Nominal Re Interaction Length Simulation θ Experiment 4600 4.60 - Touber and Sandham 4600 4.80 Periodic/LES-SGS Agostini et al. 4600 3.45 Periodic/LES-SGS Morgan et al. 4800 2.90 Periodic/ILES Morgan et al. 2300 2.90 Periodic/ILES Pirozzoli and Bernardini 2300 2.89 Periodic/DNS Present work 4600 2.94 Periodic/ILES Present work 2300 3.06 Periodic/ILES a Non-dimensionalized by δ , i.e., L /δ 99 int 99 A qualitative comparison with the experiment is presented in figure 4, where the coordinates have been scaled by the respective interaction lengths. Key features like the undisturbed boundary-layer, turning of the flow over the separation bubble, and the boundary-layer recovery can be observed in figures 4(a) and 4(b). Locations of high stress also coincide, however the magnitudes differ due to the aforementioned underprediction of the interaction length. Undisturbed boundary-layer stress profiles, figure 5, show that the streamwise stress is in good agreement with the data and the near-wall peak is captured, but disagreement exists in the outer part of the boundary layer. The discrepancy begins at y/δ = 0 . 3, which happens to be the boundary-layer displacement thickness. An opposite trend was observed with the wall-normal stress, i.e., better agreement was obtained in the outer part of the boundary layer. Touber and Sandham reported similar behavior in their work, even though they used a zonal digital filter approach in the wall-normal direction. A separate study, results not included here, examined the impact of digital filter length scales (cf. table 1) and possible improvements in generation of fluctuations at the inflow plane. Among items under further consideration are zonal implementation to impose the relevant length scales in the wall-normal direction and in situ calculation of relevant length scales in the streamwise direction.

B. Reynolds-stress Budgets 1. Validation of Budgets The budget of the Reynolds-stress transport and the imbalance in the budget was calculated at the undis- turbed upstream location ( x − x ) /δ = − 5 . 7 for Re = 2300 and 4600. The imbalance highlights the imp 99 θ relative magnitude of each term and should be zero if production, transfer, and dissipation mechanisms that govern the energy cascade are appropriately captured. A non-zero value is an indication of the error in the solution. To calculate the imbalance, unsteady term was assumed to be negligible in equation 2 and the convective term was subtracted from the right-hand side to yield equation 39.

∂ ′′ ′′ Imbalance = P + D − ¯ ρ + Π + M − ( ρu u ˜ u ) = 0 (39) ij ij ij ij ij k i j ∂x k The turbulent kinetic energy (TKE) budget was calculated by taking one-half of the trace of the Reynolds- stress budget. The present TKE budget is similar to the Mach 2.5 DNS of Guarini et al. and the Mach 4 DNS of Sinha et al. Both examined a turbulent boundary layer and are thus relevant to the current work. A comparison of TKE budget with the finely-resolved incompressible LES of Eitel-Amor (figure 6) at Re = 4400 shows good agreement, especially in the viscous sublayer and log-law regions. Differences in the θ + peak of production and troughs of molecular and turbulent diffusion within the buffer region (10 ≤ y ≤ 30) can be attributed to the vast difference in mesh resolution. The terms related to compressibility, i.e., pressure strain and turbulent mass flux were small and thus omitted from figure 6. The imbalance reached a peak in the buffer region with a value of 0 . 056, approximately four times larger than the LES of Eitel-Amor.

13 of 25 American Institute of Aeronautics and Astronautics ′′ ′′ k = tr ( ρu u ) (40) i j Budgets are plotted for the normal stresses and the primary shear stress in figure 7. To lend confidence to the present calculations, the results are plotted alongside a turbulent boundary layer spectral DNS of ¨ Schlatter and Orl¨ u. The present results also look similar to the turbulent channel flow DNS of Mansour 19 18 20 et al., Spalart’s DNS of a turbulent boundary layer, and turbulent square duct DNS of Huser et al.

It is noted however that these works examining the budget relied on the incompressible form of equation 2.

Moreover, the pressure-velocity correlation was expressed by combining the pressure diffusion (equation 7) and pressure strain (equation 8) terms as shown in equation 41, so some discrepancy in figure 7 is due to ¨ the fact that the available Schlatter and Orl¨ u data is of the total velocity-pressure correlation.

′ ′ ∂p ∂p tot P ′′ ′′ Π = D + Π = − u − u (41) ij ij ij i j ∂x ∂x j i Since the following comparison is presented in the undisturbed part of the upstream boundary layer at ′′ ′′ ′′ ′′ ′′ ′′ ( x − x ) /δ = − 5 . 7, ρu w and ρv w were negligible as expected. The ρu u budget shows molecular imp 99 diffusion in balance with dissipation at the wall, however, away from the wall and in the buffer region, peak production is balanced by dissipation, pressure strain, and diffusion terms. In the log-law region, molecular and turbulent diffusion become less dominant, thus production is balanced by dissipation and pressure strain ′′ ′′ ′′ ′′ terms. In the ρv v and ρw w budgets, figures 7(b) and 7(c), the pressure strain term assumed a dominant role and governed the transfer of energy to these normal stresses. Huser et al. concluded that near the wall pressure strain and pressure diffusion are of opposite value and thus cancel, so the pressure-velocity should be modeled using equation 41, however that is only true for off-diagonal components of these two terms. While off-diagonal terms do not vanish at the wall, the diagonal pressure strain naturally goes to zero and pressure diffusion is negligible at the wall as shown by the present calculation. Knowing this behavior Mansour et al. suggested that a possible model for pressure-velocity correlation should be based on the trace of equation 41. However, the role of the pressure strain term as an energy transfer mechanism is the reason that the term (equation 8) is directly modeled in recent contemporary Reynolds-stress turbulence models.

′′ ′′ ′′ ′′ For both ρv v and ρw w , the pressure strain term is balanced by dissipation in the region away from ′′ ′′ the wall. In the case of ρv v the peak in the pressure strain occured at the beginning of the log-law ′′ ′′ region, indicating wall-normal ejection of energy packets, while the peak in ρw w occurs at the end of the ′′ ′′ viscous sublayer. Finally, the budget of ρu v , examined in figure 7(d), shows that in the viscous sublayer region a positive dissipation is balanced by a negative molecular diffusion and away from the wall a negative production is balanced by positive turbulent diffusion, pressure diffusion, and pressure strain. This may seem counterintutive but because shear stress can either be positive or negative, a negative value of larger magnitude means higher shear stress and vice versa.

′′ ′′ The effect of the Reynolds number scaling is investigated by comparing ρu u at three simulated Reynolds numbers (figure 8) for the coarse mesh. The lowest peak imbalance in the sum of the budget for the smallest Reynolds number was indicative of the fact that the coarse mesh was sufficiently resolved. But for the higher Reynolds numbers, small increases in the peak imbalance were observed, which corresponded to the smaller dissipation term. Thus, suggesting that the simulations at higher Reynolds numbers were approaching the limit where ILES would still be valid and may potentially benefit from an increase in the mesh resolution.

The largest Reynolds number had a peak in the budget imbalance of 0 . 09 for the coarse mesh. The peak in the imbalance dropped to 0 . 08 for the fine mesh, implying that the meshes used in the present work are sufficient for ILES and significantly large numbers of grid points would be necessary to drive the imbalance to zero. This was further supported by the similarity in wall properties between the present fine mesh and 53 53 Poggie’s Grid 2 and Grid 3. Poggie performed a detailed mesh resolution study for a turbulent boundary layer at similar flow conditions, where Re ranged from approximately 1300 − 2000 and mesh sizes ranged θ 7 10 from 1 . 1 × 10 to 3 . 3 × 10 . Although budgets were not computed in his work, the Reynolds stress profiles showed no change beyond Grid 3. Thus, from this point onward, the fine mesh results will be presented for the budgets.

14 of 25 American Institute of Aeronautics and Astronautics 2. Budgets in Shock wave/Boundary-layer Interaction The budgets of Reynolds stress within the SBLI region are presented next. Figures 9 to 11 show compar- isons at the pre-reflected-shock, post-reflected-shock, separation bubble, and post-impinging-shock locations.

Since the friction velocity becomes small in the SBLI region, normalizing the budgets by ρ u /ν becomes w w τ ineffective. In the absence of a new means of normalization, this section omitted the aforementioned nor- malization in favor of the non-dimensionalization employed in equations (17a) and (17b), i.e., freestream conditions ρ U /ν . The ordinate axis still utilizes the inner scaling since the friction velocity only has a ∞ ∞ ∞ + first-order effect on y .

Budgets in figure 7 are replotted in figure 9 with the freestream normalization (left panel) and compared ′′ ′′ with the post-reflected-shock budgets (right panel). The ρu u component in figures 9(a) and 9(b) showed a threefold increase in the peak production, which was balanced by the turbulent diffusion and an amplification in the pressure strain term. The increased dissipation at the wall was balanced by the molecular diffusion and pressure strain terms. Previously inactive terms like turbulent mass flux and convection also showed a significant increase. This resulted in a net gain in the budget imbalance for the component, however this was expected. The unsteady term in equation 2 was assumed to be zero in order to calculate the imbalance, while such an approach was valid for the undisturbed boundary layer, the assumption wouldn’t be accurate in the presence of characteristically low-frequency oscillations of the reflected shock/separation region. The abrupt increase in the production leads to a lag in the energy transfer mechanisms, followed by the increased anisotropy and the budget imbalance. In other words, the budget imbalance could be thought of as the unsteady contribution, which becomes significant in the interaction region.

′′ ′′ ′′ ′′ ρv v and ρw w , are shown in figures 9(c) to 9(f). Notice that the The secondary normal stresses, post-reflected-shock budgets have ordinate axes an order of magnitude larger than the pre-reflected-shock budgets. The proportional behavior of pressure diffusion and pressure strain terms remained the same for ′′ ′′ ρv v in the near wall region, while away from the wall, convection, turbulent diffusion, and production ′′ ′′ terms came to prominence. For ρw w , the near-wall dissipation was balanced by molecular diffusion and pressure strain terms, a relationship that remained proportional. Here, however, only the convection term became important away from the wall. It is important to point out that the role pressure strain and pressure diffusion terms play, as agents of energy transfer, remained unchanged and became amplified in the presence + of the shock. The artifact present at y ∼ 250 was due to the profiles crossing the reflected shock, thus ′′ ′′ ′′ ′′ significant for ρv v and absent for ρw w due to the 2D nature of the shock.

Figure 10 shows budgets in the separation bubble (left panel) and the post-impinging-shock recovery ′′ ′′ region (right panel). Examining ρu u in figures 10(a) and 10(b), it was clear that in the near-wall region molecular diffusion and pressure strain terms balanced dissipation, which established the key role pressure strain played in the interaction region and it’s residual importance in the post-impinging-shock recovery + region. The peak production at y ∼ 50 corresponds to the high stress region (figure 4(c)) above the separation bubble. This high rate of production was balanced by pressure strain, turbulent transport, and + convection terms. Further away from the wall at y ∼ 100, another area of activity corresponds to the impinging shock, where pressure diffusion and pressure strain terms become prominent and balanced by + production and turbulent mass flux terms. Lastly, the fluctuation at y ∼ 250 was due to the reflected shock. Similarly in the recovery region a two-peak behavior was observed for the production term. The + first peak at y ∼ 10 was representative of a typical turbulent boundary layer, but lower in magnitude since + the boundary layer was not fully recovered. While the second peak at y ∼ 300 was due to the high stress region that continued past the interaction region. The secondary normal stresses showed a similar behavior as ( x − x ) /δ = − 2 . 0, however the terms in the recovery region had lower magnitudes as expected.

imp 99 ′′ ′′ Finally, the primary shear stress, ρu v , is presented in figure 11 for all four locations. Notice that the ordinate axis is an order of magnitude higher for budgets in the interaction region. It was also clipped to highlight terms with smaller magnitudes, however the near-wall peak values of pressure strain and pressure diffusion terms are stated on the plots. Both terms maintain proportionality at the wall and cancel each other, but a large increase in the magnitude was observed. Away from the wall and in the undisturbed boundary layer, pressure strain and pressure diffusion terms jointly balanced the negative production term, whereas in the interaction and recovery regions a new layer of high pressure strain and pressure diffusion developed, albeit with reversed signs compared to the wall. This behavior was most prominent in the post- reflected-shock region and diminished at aft locations. Previously benign terms like turbulent transport and turbulent mass flux also played a role to balance the negative production term, especially in the region confined between the separation bubble and the crossing of the impinging and reflected shocks.

15 of 25 American Institute of Aeronautics and Astronautics

V. Conclusions

Implicit large-eddy simulations (ILES) of a two-dimensional shock wave/boundary-layer interaction (SBLI) were performed. Comparisons with the experiment showed that the interaction region was underpre- dicted. This was attributed to the role sidewalls and corner influences played on the tunnel centerline, despite the fact that the experiment was claimed to be void of three-dimensional effects. This conclusion was sup- ported by contemporary studies of other authors and could be validated by performing a three-dimensional, full-span, two-sidewall simulation, currently under consideration. Thus, the present simulation, which em- ployed periodic boundary conditions in the spanwise direction was truly two-dimensional and valuable for turbulence model development and validation.

Budgets were computed for the Reynolds-stress transport equation and validated against previous highly- resolved LES and DNS simulations of a turbulent boundary layer, channel flow, and a square duct. Grid resolution in the y and z directions was deemed sufficient for ILES, thus grid resolution in only the x direction was studied. It was shown that the coarse mesh was adequate for the current flow conditions. However, an investigation to assess the ability of the coarse mesh to resolve the scale of turbulence by reducing the Reynolds number of the simulation by half and one-quarter showed that the budget imbalance dropped as the Reynolds number was reduced. This indicated that the imbalance in the budget shown for the undisturbed boundary layer, at the largest Reynolds number, was due to the mesh not supporting the smallest of the scales to dissipate the energy. Interestingly, this also meant that the approach of using a high-order low-pass filter in lieu of a subgrid-scale model was valid as it added minimal numerical dissipation, but this was not quantified here.

Budgets within the SBLI showed an advantage in studying the pressure strain and pressure diffusion terms separately. They cancelled each other at the wall but their behavior was different away from the wall.

While the pressure diffusion term remained small but non-zero for the normal stress budgets, it assumed an important role in the primary shear stress budgets and in one instance it was the same order of magnitude as the pressure strain term near the wall. Thus, Reynolds stress turbulence models could show improvement by accounting for the pressure diffusion term rather than the current practice of some model developers to assume it negligible or include it into the model for turbulent diffusion. Moreover, the pressure strain term continued to play a key role in the transfer of energy from the primary normal stress to the secondary normal stresses, which makes it an important term to model, as is the case with Reynolds stress turbulence models in use today.

Acknowledgments

The authors would like to thank NASA’s Transformative Aeronautics Concepts Program and Trans- formational Tools and Technology Project for its generous support. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.

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. 0496 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.

Babinsky, H., Oorebeek, J., and Cottingham, T. G., “Corner Effects in Reflecting Oblique Shock-wave/Boundary-layer Interactions,” No. 0859 in 51st Aerospace Sciences Meeting, American Institute of Aeronautics and Astronautics, January 2013.

Benek, J. A., Casimir J. Suchyta, I., and Babinsky, H., “The Effect of Wind Tunnel Size on Incipient Shock Bound- ary Layer Interaction Experiments,” No. 0862 in 51st Aerospace Sciences Meeting, American Institute of Aeronautics and Astronautics, January 2013.

Benek, J. A., Casimir J. Suchyta, I., and Babinsky, H., “The Effect of Wind Tunnel Size and Shock Strength on 16 of 25 American Institute of Aeronautics and Astronautics Incident Shock Boundary Layer Interaction Experiments,” No. 3336 in 44th Fluid Dynamics Conference, American Institute of Aeronautics and Astronautics, June 2014.

Eagle, W. E., Driscoll, J. F., and Benek, J. A., “Experimental Investigation of Corner Flows in Rectangular Supersonic Inlets with 3D Shock-Boundary Layer Effects,” No. 0857 in 49th Aerospace Sciences Meeting, American Institute of Aeronautics and Astronautics, January 2011.

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.

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.

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.

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.

Vyas, M. A., Waindim, M., and Gaitonde, D. V., “Budget of Turbulent Kinetic Energy in a Shock Wave/Boundary-layer Interaction,” No. 3187 in 46th Fluid Dynamics Conference, American Institute of Aeronautics and Astronautics, June 2016.

Spalart, P. R., “Direct Simulation of a Turbulent Boundary Layer up to R = 1410,” Journal of Fluid Mechanics , θ Vol. 187, May 1988, pp. 61–98.

Mansour, N. N., Kim, J., and Moin, P., “Reynolds Stress and Dissipation Rate Budgets in a Turbulent Channel Flow,” Journal of Fluid Mechanics , Vol. 194, June 1988, pp. 15–44.

Huser, A., Biringen, S., and Hatay, F. F., “Direct Numerical Simulation of Turbulent Flow in a Square Duct: Reynolds Stress Budgets,” No. 0187 in 32nd Aerospace Sciences Meeting, American Institute of Aeronautics and Astronautics, January 1994.

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.

Klein, M., Sadiki, A., and Janicka, J., “A Digital Filter Based Generation of Inflow Data for Spatially Developing Direct Numerical or Large Eddy Simulations,” Journal of Computational Physics , Vol. 186, February 2003, pp. 652–665.

Veloudis, I., Yang, Z., McGuirk, J. J., Page, G. J., and Spencer, A., “Novel Implementation and Assessment of a Digital Filter Based Approach for the Generation of LES Inlet Conditions,” Flow, Turbulence and Combustion , Vol. 79, March 2007, pp. 1–24.

Xie, Z.-T. and Castro, I. P., “Efficient Generation of Inflow Conditions for Large Eddy Simulation of Street-scale Flows,” Flow, Turbulence and Combustion , Vol. 81, April 2008, pp. 449–470.

Pope, S. B., Turbulent Flows , Cambridge University Press, 8th ed., 2011.

Touber, E. and Sandham, N. D., “Oblique Shock Impinging on a Turbulent Boundary Layer: Low-Frequency Mecha- nisms,” No. 4170 in 38th Fluid Dynamics Conference, American Institute of Aeronautics and Astronautics, June 2008.

Lund, T. S., Wu, X., and Squires, K. D., “Generation of Turbulent Inflow Data for Spatially-developing Boundary-layer Simulations,” Journal of Computational Physics , Vol. 140, 1998, pp. 233–258.

Gerolymos, G. A., “Implicit Multiple-grid Solution of the Compressible Navier-Stokes Equations Using k −  Turbulence Closure,” AIAA Journal , Vol. 28, No. 10, October 1990, pp. 1707–1717.

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.

Lele, S. K., “Compact Finite Difference Scheme with Spectral-like Resolution,” Journal of Computational Physics , Vol. 103, August 1992, pp. 16–42.

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.

Bisek, N. J., “High-order Implicit Large-eddy Simulations of a Supersonic Corner Flow over a Compression Ramp,” No.

3335 in 44th Fluid Dynamics Conference, American Institute of Aeronautics and Astronautics, June 2014.

Garmann, D. J., Visbal, M. R., and Orkwis, P. D., “Comparative Study of Implicit and Subgrid-scale Model Large-eddy Simulation Techniques for Low-Reynolds Number Airfoil Applications,” International Journal for Numerical Methods in Fluids , Vol. 71, August 2013, pp. 1546–1565.

Beam, R. M. and Warming, R. F., “An Implicit Factored Scheme for the Compressible Navier-Stokes Equations,” AIAA Journal , Vol. 16, No. 4, April 1978, pp. 393–402.

17 of 25 American Institute of Aeronautics and Astronautics Morgan, B., Duraisamy, K., Nguyen, N., Kawai, S., and Lele, S. K., “Flow Physics and RANS Modelling of Oblique Shock/Turbulent Boundary Layer Interaction,” Journal of Fluid Mechanics , Vol. 729, July 2013, pp. 231–284.

Pirozzoli, S. and Bernardini, M., “Direct Numerical Simulation Database for Impinging Shock Wave/Turbulent Boundary- layer Interaction,” AIAA Journal , Vol. 49, No. 6, June 2011, pp. 1307–1312.

Kawai, S., Shankar, S. K., and Lele, S. K., “Assessment of Localized Artificial Diffusivity Scheme for Large-eddy Simu- lation of Compressible Turbulent Flows,” Journal of Computational Physics , Vol. 229, November 2010, pp. 1739–1762.

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.

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.

Garnier, E., Sagaut, P., and Deville, M., “Large Eddy Simulation of Shock/Boundary-layer Interaction,” AIAA Journal , Vol. 40, No. 10, October 2002, pp. 1935–1944.

¨ Eitel-Amor, G., Orl¨ u, R., and Schlatter, P., “Simulation and Validation of a Spatially Evolving Turbulent Boundary Layer up to Re = 8300,” International Journal of Heat and Fluid Flow , Vol. 47, 2014, pp. 57–69.

θ Inagaki, M., Kondoh, T., and Nagano, Y., “A Mixed-Time-Scale SGS Model With Fixed Model Parameters for Practical LES,” Journal of Fluids Engineering , Vol. 127, January 2005, pp. 1–13.

Agostini, L., Larchevˆ eque, L., Dupont, P., Debi` eve, J.-F., and Dussauge, J.-P., “Zones of Influence and Shock Motion in a Shock/Boundary-layer Interaction,” AIAA Journal , Vol. 50, No. 6, June 2012, pp. 1377–1387.

Larchevˆ eque, L., Dupont, P., de Martel, E., Garnier, E., and Debieve, J.-F., “Experimental and Numerical Study of Unsteadiness in Boundary Layer/ Shock Wave Interaction,” Turbulence and Interactions , , No. 110, 2010, pp. 263–269.

Touber, E. and Sandham, N. D., “Large-eddy Simulations of Low-frequency Unsteadiness in a Turbulent Shock-induced Separation Bubble,” Theoretical and Computational Fluid Dynamics , Vol. 23, May 2009, pp. 79–107.

Sinha, K., Mart´ ın, M. P., and Candler, G. V., “Assessment of the k −  Turbulence Model for Compressible Flows using Direct Simulation Data,” No. 0730 in 39th Aerospace Sciences Meeting, American Institute of Aeronautics and Astronautics, January 2001.

¨ 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.

C´ ecora, R.-D., Radespiel, R., Eisfeld, B., and Probst, A., “Differential Reynolds-stress Modeling for Aeronautics,” AIAA Journal , Vol. 53, No. 3, March 2015, pp. 739–755.

Poggie, J., Bisek, N. J., and Gosse, R., “Resolution Effects in Compressible, Turbulent Boundary Layer Simulations,” Computer and Fluids , Vol. 120, July 2015, pp. 57–69.

18 of 25 American Institute of Aeronautics and Astronautics (a) u (b) v √ ′ ′ (c) u u √ ′ ′ (d) v v ′ ′ (e) − u v Figure 4. Mean velocity, normal stress, and shear stress comparisons at Re = 4600 . Simulation is shown with θ labeled lines and overlaid on experimental data contours. Both share the same legend.

19 of 25 American Institute of Aeronautics and Astronautics √ √ ′ 2 ′ 2 ′ ′ (a) u and v (b) − u v Figure 5. Stress profiles comparison at ( x − x ) /δ = − 5 . 7 for Re = 4600 . Markers without lines represent imp 99 θ experimental data.

Figure 6. Budget of the turbulent kinetic energy at ( x − x ) /δ = − 5 . 7 . Present results (solid lines) at imp 99 Re = 4600 compared with LES of Eitel-Amor (dashed lines) at Re = 4400 .  : Production, 4 : Molecular θ θ diffusion, O : Turbulent diffusion, ♦ : dissipation, and Σ : Sum. Budgets have been normalized by ρ u /ν .

w w τ 20 of 25 American Institute of Aeronautics and Astronautics ′′ ′′ ′′ ′′ (a) ρu u (b) ρv v ′′ ′′ ′′ ′′ (c) ρw w (d) ρu v Figure 7. Reynolds-stress budget at ( x − x ) /δ = − 5 . 7 . Present results (solid lines) at Re = 4600 compared imp 99 θ ¨ with DNS of Schlatter and Orl¨ u (dashed lines) at Re = 4000 .  : Production, 4 : Molecular diffusion, O : θ Turbulent diffusion, . : Pressure diffusion, / : Pressure strain, ♦ : dissipation, © : Turbulent mass flux, X: Convection, and Σ : Sum. Terms smaller than three orders of magnitude are omitted from the plots for clarity.

Budgets have been normalized by ρ u /ν .

w w τ 21 of 25 American Institute of Aeronautics and Astronautics ′′ ′′ ′′ ′′ (a) ρu u at Re = 1150 (b) ρu u at Re = 2300 θ θ ′′ ′′ (c) ρu u at Re = 4600 θ ′′ ′′ Figure 8. Reynolds-stress budget at ( x − x ) /δ = − 5 . 7 for the coarse mesh. ρu u for present results imp 99 compared at Re = 1150 , 2300 , and 4600 .  : Production, 4 : Molecular diffusion, O : Turbulent diffusion, . : θ Pressure diffusion, / : Pressure strain, ♦ : dissipation, © : Turbulent mass flux, X: Convection, and Σ : Sum.

Terms smaller than three orders of magnitude are omitted from the plots for clarity. Budgets have been normalized by ρ u /ν .

w w τ 22 of 25 American Institute of Aeronautics and Astronautics ′′ ′′ ′′ ′′ (a) ρu u at ( x − x ) /δ = − 5 . 7 (b) ρu u at ( x − x ) /δ = − 2 . 0 99 99 imp imp ′′ ′′ ′′ ′′ (c) ρv v at ( x − x ) /δ = − 5 . 7 (d) ρv v at ( x − x ) /δ = − 2 . 0 imp 99 imp 99 ′′ ′′ ′′ ′′ (e) ρw w at ( x − x ) /δ = − 5 . 7 (f) ρw w at ( x − x ) /δ = − 2 . 0 imp 99 imp 99 Figure 9. Reynolds normal stress budget for pre-reflected-shock (left panel) and post-reflected-shock (right panel). Coarse mesh at Re = 4600 .  : Production, 4 : Molecular diffusion, O : Turbulent diffusion, . : Pressure θ diffusion, / : Pressure strain, ♦ : dissipation, © : Turbulent mass flux, X: Convection, and Σ : Sum. Terms smaller than three orders of magnitude are omitted from the plots for clarity. Budgets are non-dimensionalized by the freestream ( ρ U /ν ), so inner scaling is not applied.

∞ ∞ ∞ 23 of 25 American Institute of Aeronautics and Astronautics ′′ ′′ ′′ ′′ (a) ρu u at ( x − x ) /δ = − 0 . 90 (b) ρu u at ( x − x ) /δ = 1 . 44 99 99 imp imp ′′ ′′ ′′ ′′ (c) ρv v at ( x − x ) /δ = − 0 . 90 (d) ρv v at ( x − x ) /δ = 1 . 44 imp 99 imp 99 ′′ ′′ ′′ ′′ (e) ρw w at ( x − x ) /δ = − 0 . 90 (f) ρw w at ( x − x ) /δ = 1 . 44 imp 99 imp 99 Figure 10. Reynolds normal stress budget for separation (left panel) and post-impinging-shock (right panel).

Coarse mesh at Re = 4600 .  : Production, 4 : Molecular diffusion, O : Turbulent diffusion, . : Pressure diffusion, θ / : Pressure strain, ♦ : dissipation, © : Turbulent mass flux, X: Convection, and Σ : Sum. Terms smaller than three orders of magnitude are omitted from the plots for clarity. Budgets are non-dimensionalized by the freestream ( ρ U /ν ), so inner scaling is not applied.

∞ ∞ ∞ 24 of 25 American Institute of Aeronautics and Astronautics ′′ ′′ ′′ ′′ (a) ρu v at ( x − x ) /δ = − 5 . 7 (b) ρu v at ( x − x ) /δ = − 2 . 0 imp 99 imp 99 ′′ ′′ ′′ ′′ (c) ρu v at ( x − x ) /δ = − 0 . 9 (d) ρu v at ( x − x ) /δ = 1 . 44 imp 99 imp 99 Figure 11. Reynolds shear stress budget for all locations in figures 9 and 10. Coarse mesh at Re = 4600 .

θ  : Production, 4 : Molecular diffusion, O : Turbulent diffusion, . : Pressure diffusion, / : Pressure strain, ♦ : dissipation, © : Turbulent mass flux, X: Convection, and Σ : Sum. Terms smaller than three orders of magnitude are omitted from the plots for clarity. Budgets are non-dimensionalized by the freestream ( ρ U /ν ), so inner ∞ ∞ ∞ scaling is not applied.

25 of 25 American Institute of Aeronautics and Astronautics