Computational study of a complex three-dimensional shock boundary-layer interaction

Shock boundary–layer interactions occur in many high-speed aerodynamic flows and they can have a notable impact on design considerations due to the aerodynamic and heat transfer effects. Consequently there is a notable interest in understanding the ability of computational tools to calculate the complex flow fields that can arise in a range of engineering applications. Three-dimensional complex shock boundary layer interaction studies are expensive in both time and computational resources. Although recent studies have begun to focus on the use of more complex computational methods such as large eddy simulations, the aim of this research is to assess the ability of steady Reynolds averaged Navier Stokes turbulence models to simulate the interaction of a planar shock impinging on a cylindrical body under supersonic conditions and to determine if these models have a role to play in engineering design applications. The performance of both eddy viscosity and Reynolds stress models are evaluated relative to an established experimental test case. The impact of Reynolds number and impinging shock strength are also considered. Of the eddy viscosity models it was shown that the Spalart-Allmaras model is unsuitable for this complex interaction and that the k- and Reynolds stress methods both gave notably better agreement with the measured surface static pressures. Overall it was considered that the Reynolds stress method was the best model as it also provided better agreement with the measured surface flow topology. It was concluded that, although a steady Reynolds averaged Navier Stokes approach has known limitations for this type of complex interaction, within an engineering context it can also provide useful results when applied appropriately.


Introduction
The topic of shock boundary layer interactions (SBLI) is a broad subject which has been the focus of extensive research over many years.It is a complex, multifaceted topic, which is not yet fully understood and still causes significant difficulties in the design of a wide range of highspeed vehicles and aerospace components (Allen, Heaslet, & Nitzberg, 1947;Borovoy et al., 2013;Brosh & Kussoy, 1983;Chaplin et al., 2011;Délery, 1999;Donaldson, 1944;Holden, Wadhams, MacLean, & Mundy, 2010).Initial investigations of nominally two-dimensional SBLIs were previously undertaken for aerofoils under supersonic conditions.For an aerofoil under high transonic conditions, a normal shock wave can form and interact with the boundary layer, typically on the suction side (Babinsky & Harvey, 2011).In this interaction, compression waves form near the wall as the boundary layer is unable to support the rise in static pressure across the shock.As the shock strength increases, along with the concomitant adverse pressure gradient, boundary layer separation can occur.Further instances of nominally two-dimensional SBLIs have been investigated extensively, most notably *Corresponding author.Email: d.g.macmanus@cranfield.ac.uk for compression ramps and oblique shocks impinging on flat plates where it has been noted that the interaction depends on the nominal pressure rise across the shock as well as the state and characteristics of the approaching boundary layer (Arnal & Délery, 2004;Délery, Marvin, & Reshotko, 1986;Délery, & Coet, 1990;Délery, 1996).For example, an increase in the shock strength, by changing the ramp angle, can lead to separation of the boundary layer if the pressure rise across the shock is high enough.Separation results in an influence from the SBLI which extends further upstream than with an attached SBLI and an additional shock system of separation and reattachment shocks can arise.Similar effects arise with the impingement of a two-dimensional oblique shock onto a flat plate where the induced increase in the boundary layer thickness generates an additional shock wave ahead of the primary reflected shock which in turn crosses the impinging shock wave (Arnal & Délery, 2004).If the shock is of significant strength, an increase in the thickness of the boundary layer occurs and a separation bubble can also arise.Thus, in nominally two-dimensional SBLIs, shock strength and angle, dependant on Mach number and geometry, along with the Reynolds number, state and shape factor of the approaching boundary layer are key parameters as to the nature of the shock interaction.
Nominally three-dimensional SBLIs are typically more complex and occur in many high-speed aerodynamic flows (Brosh & Kussoy, 1983;Pamadi et al., 2005).Examples include planar shocks impinging onto axisymmetric bodies (Brosh & Kussoy, 1983) and conical shock waves impinging onto flat plates (Gai & Teh, 2000).For these configurations, along the central symmetry plane there are broad similarities with the nominally two-dimensional SBLIs such as the local boundary layer thickening which can lead to separation.Critical factors which affect the likelihood and magnitude of separation are shock strength and shock angle which are typically controlled by the distance of the shock generator from the boundary layer, approaching Mach number and the characteristics of the approaching boundary layer.For these configurations the increase in surface static pressure from the SBLI typically arises several boundary layer thicknesses upstream relative to the nominal ideal shock impingement point and locally strong cross-flows can occur due to the SBLI.Distinct regions arise within the flow field which depend on the initial SBLI and the subsequent evolution of the boundary layer, shock waves, expansion systems, crossflow and indeed the main flow (Chaplin et al., 2011).Adverse pressure gradients at the initial interaction point can be affected by lateral, or azimuthal, flow migration away from the position of peak shock intensity which results in potent cross-flows (Brosh & Kussoy, 1983).In the case of a planar oblique shock impinging onto a cylindrical body (Brosh & Kussoy, 1983), the coalescing of the diffracted shocks on the leeward side resulted in an additional separation region and thereby deflected the flow in a similar manner to an obstacle in the boundary layer (Peake & Tobak, 1982;Sedney & Kitchens, 1975).Further detailed studies of three-dimensional SBLIs are reported by other researchers (Arnal & Délery, 2004;Babinsky & Harvey, 2011;Peake & Tobak, 1982;Saric et al., 1996).

Computational studies of SBLIs
SBLIs comprise a very wide range of configurations and flow regimes and the aerodynamic characteristics of the interaction is highly dependent on the defining conditions.The interactions can encompass aspects such as transition, separation, highly skewed boundary layers, reflected shocks and expansion fans, vortical flows and strongly three-dimensional flow gradients (Babinsky & Harvey, 2011).Furthermore, it is considered that all SBLIs are fundamentally unsteady to some degree (Oliver, Lillard, Schwing, Blaisdell, & Lyrintzis, 2007).Overall these features present a very challenging problem for computational fluid dynamics (CFD).Nevertheless, a broad variety of SBLI problems have been tackled using a range of computational methods with varying degrees of success (Bhagwandin & DeSpirito, 2011;Dolling, 2001;Oliver et al., 2007).These methods have included various forms of Reynolds Averaged Navier Stokes (RANS) approaches (DeBonis et al., 2010;Thivet, 2002;Vallet, 2007), unsteady RANS (URANS) (Barakos, Doerffer, Hirsch, Dussauge, & Babinsky, 2010;Hirsch, 2010a;Sandham, 2010), Detached Eddy Simulations (DES) (Garnier, 2009;Shams & Comte, 2010), Large Eddy Simulations (LES) (DeBonis et al., 2010;Eagle, Driscoll, & Benek, 2012;Hadjadj, 2012) and Direct Numerical Simulations (DNS) (Adams, 2000;Tokura & Maekwa, 2011).Clearly these methods provide different levels of modeling fidelity and sophistication along with concomitant resource demands.Oliver et al. (2007) examined the ability of some steady RANS methods, Spalart-Allmaras (SA), Menter's Shear Stress Transport (SST) and Olsey and Coakley's Lag model, to calculate SBLIs occurring on flat plates and a compression ramp at Mach numbers ranging between 2.87 and 6.0.Simulations of relatively benign SBLIs produce results that are in reasonable agreement with experimental results.However, as shock strength increases, the SBLI becomes more difficult to determine computationally and the difference in accuracy of various RANS models becomes more apparent.It was found that the calculation of too large a separation bubble as well as a region of upstream influence which arises in two-dimensional calculations can be improved with a three-dimensional computational model as this includes relief effects as well as the impact of the wind tunnel walls which arises in many experimental configurations.For these cases, RANS models generally calculate the inviscid flow field as well as the overall effects such as separation and the region of upstream influence to a level deemed adequate for most engineering design work.However, Oliver et al. (2007) showed that the thermal environment and interactions proved more difficult to compute.In addition, the detailed flow visualization and skin-friction-related distributions were adversely affected even if surface pressures correlated closely with experimental data.It was noted that the lowfrequency shock unsteadiness, a result of large streamwise turbulent fluctuations observed in the experiment, was not captured by RANS models due to the averaging of turbulent fluctuations, and thus, there was no physical source for such effects in RANS models.This problem was partially addressed by reducing the grid aspect ratio near the separation region and of the three models tested, the simple formulation of the one-equation SA model providing the lowest fidelity, with the lag model producing results in-between the SA and SST models (Oliver et al., 2007).
Subsequent studies on the ability of steady RANS methods to calculate the SBLI flow field showed that the accuracy depends greatly on the turbulence model (DeBonis et al., 2010).It was reported that for weak SBLIs of planar shocks impinging on a flat plate, differences in the Downloaded by [Cranfield University] at 05:31 01 April 2016 calculated streamwise velocity distributions were within 0.5% of each RANS turbulence model considered for the Menter SST, k-ω variants and Spalart-Allmaras models.However, as the shock intensity increased, the error in all the solutions also increased, although the relative error between each turbulence model was relatively consistent.Grid construction, either structured or unstructured, showed no discernible difference on the calculated flow fields.The most notable difference between models was within the separation region of the SBLI although due to the range of metrics used in the assessment it was difficult to determine a clearly superior method.LES simulations were found to produce similar error levels as RANS methods but the prediction of normal stresses was superior using LES (DeBonis et al., 2010).Investigations of SBLIs with separated flow fields using Goldberg's One-equation Rt model, SA, realizable k-l, Goldberg's realizable q-l, SST, Goldberg's k-ε-Rt and a seven equation second moment closure turbulence model in a thrust vector nozzle (Tian & Lu, 2013) showed that a k-ε model provided the most accurate results when calculating the separation point and shock wave position.
The Reynolds Stress Model (RSM) is a formulation of RANS model that solves the transport equations for Reynolds stresses directly, including an additional equation for dissipation rate, to close the RANS equations.This produces a seven-equation model that has greater computational demands yet is capable of providing higher accuracy for boundary layer flows when compared to the industry standard Eddy-Viscosity (EV) RANS models; SA, k-ε, kω and SST formulations.A study by Vallet (2007) using an RSM model with a 2D compression ramp reported improvements over k-ε models and was able to calculate mean-flow distributions for both separated and attached flows.The inability to correctly simulate experimentally observed low-frequency shock oscillations, also notably absent in computational simulations by Oliver et al. (2007), was again apparent in the foot of the shock-wave Reynolds stress profiles (Vallet, 2007).Further studies, such as those by Mendonca and Sharif (2010) found that surface roughness should not be ignored in SBLI flows due to the upstream influence and movement of the shock as surface roughness is increased.The RSM model was shown to provide the most accurate results although it was also shown that the k-ω model could provide useful results at a decreased computational cost.Comparison of RANS models for SBLIs generated with 2D oblique shocks on a flat plate (Bhagwandin & DeSpirito, 2011), showed that the RSM model is suitable for SBLI flows but that it requires a higher near-wall resolution grid than eddy-viscosity models such as SA and SST methods.This, coupled with the RSM model being a seven-equation model and hence requiring an increased amount of resources relative to the EV models, may be less desirable where rapid solutions are needed.Furthermore, numerical stability issues may arise which suggests that hybrid RANS/LES models may be useful to resolve some of the inherent unsteady effects within SBLIs.DeBonis et al. (2010) suggests that the capability of RSMs to determine the appropriate normal stress distributions warrants further investigation.
Recent case studies continue to prove useful in providing experimental data of complex SBLIs with computational comparisons.Holden et al. (2010) provided experimental results for flat plate and compression ramp SBLI datasets at Mach 4 to 11. Significant differences were reported between SA and SST models, with the former failing to predict separation for compression ramp flow and the latter over-predicting the size of the separated region.Subsequent "blind" code validation studies have provided direction for empirical modifications to RANS codes, such as those suggested by Wilcox (2006).SBLIs generated by crossing shocks from fins for Mach numbers 5 to 8 have been compared experimentally and computationally by Borovoy et al. (2013).The q-ω RANS turbulence model was found to agree well with the experimental results.It was seen that thermocouple sensors were more accurate than optical measurements and, thus, the CFD was in better agreement with the thermocouple results.
Unsteady RANS (URANS) methods have been investigated to address the low-frequency oscillations and, hence, to capture the flow unsteadiness which has been observed experimentally.An experimental oblique shock reflection study at Mach 2.0 (Sandham, 2010) was simulated computationally but failed to capture significant flow unsteadiness and the expected low-frequency oscillations.The study also demonstrated the superior ability of an LES model to capture these effects.Of further note for the RANS methods, is the effect that inlet turbulent intensity has on the size of the separation bubble.A reduction of inlet turbulence intensity from 3% to 0.1% resulted in a separation bubble with twice the length and in better agreement with experimental data (Sandham, 2010).In addition, there was better agreement with the measurements when modeling the system as 3D instead of 2D due to the inclusion of corner flows and the associated aerodynamic effects.This influence of the 3D and corner flows was also reported by Hirsch (2010b) with an oblique shock at Mach numbers 1.7, 2.0 and 2.25.URANS models were found to produce steady flow fields and did not capture high-frequency oscillations at the foot of the shock.Additionally, URANS methods were either unable to predict the natural shock motion or underestimated the effect.Distributions of velocity profiles in the region downstream of the interaction, as well as in the separated regions, were deemed to require improvement, although the upstream velocity profiles were calculated more accurately.It was reported that, to improve the accuracy of SBLIs, improvements need to be made to URANS turbulence models to address 3D separations and corner vortices (Hirsch, 2010b).Barakos et al. (2010) suggests that hybrid DES/LES models and data be used Downloaded by [Cranfield University] at 05:31 01 April 2016 to propose URANS improvements through challenging the assumption of local equilibrium in URANS and the development of new anisotropic stress tensor methods.
Although LES methods are more advanced than RANS computations, studies on a 3D SBLI of an oblique shock impinging on a flat plate at Mach 2.75 concluded that there was little advantage of LES over RANS methods for this particular configuration (Eagle et al., 2012).It was reported that this was partially affected by the three-dimensional nature of the experimental data and the use of two-dimensional computational simulations.Twodimensional studies of oblique shocks interacting with flat plates at Mach 2.75 (Eagle et al., 2012) note the more mature status and understanding of RANS methods as a driver for the similar accuracy and uncertainty in results when comparing RANS with LES.It was further observed that LES methods provide superior predictions of normal stresses when compared to RANS methods.Through experimental observation and computational simulation of a three-dimensional oblique shock interacting with a flat plate at Mach 2.28 with a shock incidence angle (β) of 32.4°, Hadjadj (2012) found excellent agreement between LES predictions and experimental results.The aforementioned low-frequency shock oscillations were observed and helped give credence to the hypothesis that the shock wave and separation bubble act as a coupled system.
Detached Eddy Simulation, originally proposed as a bridge between RANS and LES methods, has the ability to predict separated and complex turbulent flows.Shams and Comte (2010) studied SBLIs in nozzle flows and found very good agreement with experimental results.It was determined that the DES method calculated the unsteady low-frequency shock oscillations and is, thus, a potentially good method of predicting SBLIs.A SBLI accounting for the whole wind tunnel span at Mach 2.3 by Garnier (Garnier, 2009) using DES saw an improvement over 2D LES methods.Corner separations in 3D simulations were, again, noted to affect the results.Low-frequency oscillations, as with Shams and Comte (2010), were present and the results were in good agreement with the experimental data.
Direct Numerical Simulation has the potential to improve the detailed simulations of SBLIs.The computational and economical costs are, however, currently prohibitive for many engineering applications.Attempts have been made, such as by Adams (Adams, 2000), in which a compression ramp (β = 18°) at Mach 3 was directly simulated with a relatively low Reynolds number (Re θ = 1685), to enable calculations with the then current computational resources available.The simulation domain was found to be of insufficiently large enough scale to capture any large-scale shock motion and, therefore, the agreement between experimental and computational data was limited.An oblique SBLI at Mach 2.0 and Re δ * θ = 1000 was predicted using DNS by Tokura and Maekwa (2011).Good agreement between previous studies was found, although comparisons were limited due to the low Reynolds number used in the study.The calculations provide access to flow structures such as lambda-vortices and broadband spectra seen experimentally but not present in most other CFD methods.In general, DNS still requires computational resources beyond current capabilities to be fully realized as an alternative to lower fidelity engineering models.
Overall, the range of CFD methods that have been applied to SBLI flows, indicate that the complexity of the flow field is a major challenge even for the most advanced tools.Although steady RANS calculations using eddy-viscosity-based turbulence models have well-known simplifying assumptions, in some cases they are able to calculate successfully some of the primary flow features.However, there is evidence that key aspects of the interaction such as separation size, reattachment locations and details of the flow topology are not fully resolved.In spite of this, some of the more recent studies indicate that there are only modest benefits of going to the much more demanding methods such as LES or DES.There are noted benefits of adopting a DNS approach, as seen in Tokura and Maekwa (2011), but these methods are particularly resource intensive.Within this context, and from the point of view of pragmatic engineering design, there is still a significant interest in understanding the capabilities of robust CFD tools that provide a measured balance between computational cost and solution fidelity.

Three-dimensional SBLI in multibody configurations
Although the previous experimental and computational work encompasses a wide range of SBLI configurations much of the work has focused on planar shocks impinging on flat plates, glancing and crossing shock interactions, as well as local normal and compression ramp configurations.Relatively little work has been done on the interaction of an impinging oblique shock with a nonplanar body.This is a configuration that is of particular interest within the context of multibody configurations, stores and submunition separation (Chaplin et al., 2011), sabot discard, and two-stage to orbit configurations (Pamadi et al., 2005).
For multibody configurations in close proximity under high-speed flow conditions, aerodynamic interference arises and this can significantly affect the force and moment characteristics (Chaplin, MacManus, & Birch, 2010;Hung, 1985;Wilcox, 1995).The complex flowfield is primarily dominated by the shock and expansion waves, which originate from one body and impinge upon the adjacent body.The interference aerodynamics is further complicated by multiple shock reflections, shock diffraction as well as shock interactions with the viscous body vortex Downloaded by [Cranfield University] at 05:31 01 April 2016 and boundary-layer flows.The induced changes in static pressure and flow angularity across the impinging disturbances modify both the local and overall aerodynamics of a slender body in comparison with the isolated body case.Very limited information is available in the open literature on the effects of mutual interference between slender bodies at high speed.One previous investigation showed that a planar shock impinging on a cone-cylinder body at zero incidence induced changes in normal force and pitching moment coefficient of approximately 0.02 and 0.2 respectively (Wilcox, 1995).These changes were found to increase by up to an order of magnitude when the receiver body was placed at an incidence of σ = 15°.Interference effects of this order would modify the trajectory of the slender body, and would become significant if the slender body pitches (or translates) toward the generator resulting in a collision.
Previous experimental work showed that for multiple bodies in close proximity, the interaction of a conical oblique shock with an ogive-cylinder body at Mach 2.43 could result in a strong SBLI with local flow separations and notable changes in the pressure distributions around the body (Chaplin, 2010;Chaplin et al., 2010).The consequence of this was a significant change in the forces and moments on the body which leads to an adverse change in the body trajectory and ultimately in the two bodies colliding (Chaplin et al., 2010).Additional computational work by Chaplin also showed that steady RANS calculations using the k-w SST turbulence model showed very good agreement with the measured changes in forces and moments (Chaplin et al., 2011).Furthermore, comparisons between the calculated surface static pressure distributions and the measured data using pressure sensitive paints showed that these RANS calculations were able to capture the main aspects of this complex interaction.
In addition to the overall interference loads, it is important to understand the detailed underlying flow physics.Shock-body interactions have been studied previously for a number of pertinent configurations (Chaplin et al., 2011;Derunov, Zheltovodov, & Maksimov, 2008;Fedorov, Malmuth, & Soudakov, 2007;Malmuth & Shalaev, 2004;Volkov & Derunov, 2006).Of particular interest is the work of Brosh, Kussoy, and Hung (1985) and Hung (1985) who investigated a wedge-generated shock passing over a cylinder at M ∞ = 2.85.This particular configuration is the focus of the current research.Brosh et al. (1985) performed an experimental investigation using a configuration which comprised a prismatic wedge shock generator and a cylindrical body onto which the planar oblique shock impinged.
The measurements showed that the impinging shock footprint, in terms of local pressure rise, decreased as the shock diffracted around the body.In addition, the induced pressure rise on the windward reduces quickly due to the impact of expansion waves from the generator forebody.These expansion waves do not diffract to the same extent as the impinging shock and thus the leeward pressure rise associated with the diffracted shock is maintained along the body.Consequently, the difference between the strength and extent of the windward and leeward interactions significantly affects the local normal force distribution over the body.Finally, the windward pressure rise also resulted in a local boundary-layer separation and a double-reflected shock structure around the leeward separation bubble.Both studies found that due to the induced circumferential pressure gradient, a strong azimuthal crossflow occurred which resulted in a local separation on the farside of the receiver body.A similar effect was also noted by Morkovin, Migotsky, Bailey, and Phonney (1952).
This configuration of an oblique shock impinging onto a cylindrical body is of particular interest from the point of view of the specific applications as well as from the understanding of the detailed flow physics.Due to some of the specific flow physics aspects which are different from other SBLI configurations, it is also an interesting case from the point of view of understanding the capability of computational tools to calculate the pertinent flow features and characteristics.As discussed above, CFD has demonstrated mixed capabilities in this regard for other SBLI arrangements such as compression ramps, normal SBLIs, and glancing shocks.The work of Chaplin et al. (2010Chaplin et al. ( , 2011) ) showed that the overall effect of a multibody shock interaction on the forces and moments of an axisymmetric body could be determined using RANS calculations and a k-ω SST turbulence model.
The aim of this work is to examine the performance of different RANS turbulence models for the case of an impinging planar shock onto a cylindrical body.The effect of both eddy-viscosity and Reynolds Stress turbulence models is assessed and the CFD data is compared with the experimental dataset of static pressure distributions, total pressure traverses and oil flow visualizations.Following the assessment of the CFD methods, the impact of shock strength and angle as well as the flow Reynolds number is also examined.Complex 3D SBLIs are underrepresented in the literature in comparison to 2D studies, mainly due to difficulties in accurately representing many of the flow phenomena that occur beyond planar shocks impinging on flat surfaces.Additionally, current trends are in the direction of refining LES and DES methods for such aerodynamic flows which better account for some of the complex flowfield interactions.These, however, come at an increased computational cost in comparison to steady RANS methods.Consequently, it is of interest to evaluate steady RANS methods for these 3D interactions and to determine their appropriateness within the context of an engineering design application.Downloaded by [Cranfield University] at 05:31 01 April 2016

Test case
The configuration of interest in this work was the subject of a previous experimental investigation as reported by Brosh and Kussoy (1983).The configurations comprised a circular cylinder (d = 50.8mm) of 1.016 m length (L), mounted centrally in a 0.254 m wide and 0.381 m high test section of a wind tunnel.The shock generator was a 2D wedge which was positioned above the cylinder at a distance of h/d = 2.97 from the top surface (Figure 1).
The datum test flow conditions were M = 2.88, α = 16°, h = 150.9mm, and a Reynolds number (Re L ) of 18.2 × 10 6 (Table 1) based on the free-stream velocity and length of the cylindrical body.This is the configuration for which the majority of the quantitative experimental data is available (Brosh & Kussoy, 1983).The effect of the flow Reynolds number was also examined across the range of Re L = 7.28 × 10 6 to 58.3 × 10 6 , although for these configurations the free-stream Mach number also changed by modest amounts across the range of M = 2.8 to 2.95, respectively (Table 1).The influence of the shock generator angle (α) was further experimentally examined for an increased alpha of 19°, at the same Mach number (2.88) and Reynolds number (Re L = 18.2 × 10 6 ) as the datum configuration and for a lower Reynolds number (Re L = 7.28 × 10 6 ) with Mach number of 2.80 and α of 13°.The total temperature was kept relatively constant across the configurations and ranged between 101.4 K and 108.2 K.
Surface pressure taps were positioned on the cylinder at 20 mm axial intervals ( X/d = 0.39), with additional azimuthal taps at φ = 90°, 180°and 270°at 0.2 m axial spacing ( X/d = 3.9), to verify flow symmetry.Timeaveraged static and total pressure profiles were measured in the windward (φ = 0°) and leeward (φ = 180°) planes over a vertical distance from the surface of Y/d = 0.394 with a spatial resolution of Y/d = 9.8 × 10 −4 .It was assumed that total temperature was constant throughout the boundary layer in accordance with Kussoy, Horstman, and Acharya (1978).The reported measurement uncertainties in the windward and leeward planes were ± 10% for static pressure, ± 1% for Pitot pressure, ± 6% for static temperature, 12% for density, ± 3% for velocity and ± 1.97 × 10 −3 Y/d for the vertical probe movement in the Y-axis.

Computational model and method
The calculations were performed using Fluent V14.0 (SAS IP Inc., 2011a) using a steady, implicit density-based solver.The flow gradients were calculated using a least squares cell-based method with second order up-winding for the flow and turbulence properties.The convective fluxes are determined using the Roe flux difference splitting scheme.For the turbulence models four eddy-viscosity models were selected; a one equation Spalart-Allmaras model, two equation k-ε realizable and k-ω SST models and a four equation k-ω SST model with a γ -θ transition model (tSST).A seven equation RSM model using the linear-pressure strain sub-model was also evaluated.For the fluid, specific heat capacity (Cp) was set via a threecoefficient polynomial temperature equation, Hilsenralh et al. (1955), which results in γ varying with temperature.Thermal conductivity was calculated by the use of kinetictheory using the molecular and material properties of the fluid to account for the user defined Cp value.The three  coefficient method of Sutherland's law, based on a reference viscosity, temperature and effective temperature, was used to compute viscosity as a function of temperature.For all simulations, the double-precision solver was employed.

CFD domain and boundary conditions
The notation denotes 'X' as the axial axis, 'Y' as the vertical axis and 'φ' as the azimuthal direction around the cylinder and is shown in Figure 2. The computational domain is also highlighted in Figure 1.It comprises half the experimental domain, with the assumption of an X-Y symmetry plane intersecting the cylinder at φ = 0°and φ = 180°to reduce the computational requirements.No-slip wall boundary conditions are defined for the wind tunnel walls, the wedge and the cylinder surfaces.A symmetry plane was used for the 'top' surface of the computation domain which was limited to the same vertical position as the top of, and prior to, the wedge to prevent artificial boundary layer growth which could affect the shock wave characteristics.Experimental surface pressure data was limited axially, from X/d = 1.97 to X/d = 16.34.The computational domain originates as shown in Figure 2 and extends axially X/d = 16.34, with the origin X/d = 7.83 upstream of the wedge leading edge.At the inlet plane the total temperature, total pressure and Mach number are as specified in Table 1.As no information was provided on the turbulence characteristics in the wind-tunnel experiments, the values for turbulence intensity, I, and length scale, l, were determined from approximate empirical correlations based on Re D H and the estimated boundary layer thickness ahead of the shock impingement point, δ 99 (Equations 1, 2) (SAS IP Inc., 2011b).The turbulence intensity, I, was 1.98% and the length scale, l/d, was 0.096.

CFD grid
A structured multiblock approach was used with an Ogrid blocking scheme around the cylinder with a radial expansion ratio of approximately 1.2 to grow the mesh and to provide between 40 and 90 cells across the boundary layer depending on the overall mesh resolution.The axial distribution of the mesh was clustered using an exponential contraction in cell size up to the axial position of the wedge, aft of which uniform cell sizes were used.Three grids of increasing resolution were generated; coarse, medium and fine with a refinement ratio of 1.5 in all three dimensions in accordance with the methodology outlined by Roache (1998).This yielded meshes of approximately 9 × 10 5 (900 k), 3 x10 6 (3M) and 10 × 10 6 (10M) cells.To comply with turbulence models requirements, a y + less than 1 was maintained for all grids and test conditions.
The effect of the spatial discretisation was assessed using the generalized Richardson Extrapolation method (Celik et al., 2008).Grid Convergence Indices (GCIs) were assessed using nodal points of free-stream Mach number and total pressure ratio, P T /P T∞ , at X/d = 9.84 on the windward surface, φ = 0°, of the cylinder.A factor of safety of 2 was used with the data showing the GCI to Downloaded by [Cranfield University] at 05:31 01 April 2016 be less than 1% for all cases and thus the grid was determined as mesh insensitive according to the generalized Richardson Extrapolation.The characteristics of the generated impinging shock was assessed using a position in between the wedge and the cylinder (Y/d = 1.56 and X/d = from 6.89 to 10.83).The static pressure and static temperature distributions in this region were evaluated relative to the theoretical values and the effect of the grid resolution on these parameters was assessed.It was determined that the static pressure and temperature ratios across the shock agreed to within ± 1% of the theoretical value for the fine grid configuration.

Iterative convergence
Iterative convergence was assessed by considering the overall residual parameters throughout the flow domain in addition to integrated flow parameters.Overall, residuals were typically reduced to values below 10 −4 for the continuity, momentum, energy and turbulence parameters.Due to the inherent unsteady nature of SBLIs and the assumption of 'steady' turbulence models, a simulation in which residual values below 10 −4 were maintained for 300 consecutive iterations was considered to have converged.
The convergence of flow characteristics at key locations in the domain was monitored using a variety of metrics.The windward static pressure at locations ahead of the shock impingement position (P1) and in the centre of the interaction (P2), were assessed.Within the accepted convergence tolerance, the reported values of static pressure at these locations varied by typically 0.2% and − 1.6% for P1 and P2 respectively between initial stabilization at approximately 7000 iterations and the final converged value.
Integrated drag and moment coefficients for the full cylinder, C d and C m respectively, were also considered as part of the convergence assessment, and the parameters typically converged to stable values after 5000 iterations where the variation in C d and C m was 0.0005 and 0.0, respectively.A similar characteristic was also observed for y + when averaged over the cylinder the variation with further iterations was in the order of 3.8%.The total mass averaged vorticity magnitude at the outlet plane was also considered and it was found to stabilize to within 1.5%.
For the cases with different Reynolds numbers (Re L , 7.28 × 10 −6 and 58.3 × 10 −6 ), the same domain as the datum configuration, with the 3M cell mesh, was used.For these cases, the y + values over the whole cylinder were assessed and also found to be less than 1.For the configuration where the shock generator wedge angle increased to α = 19°, the geometry was slightly modified to accommodate this change.However, the main parameters of the 3 million cell mesh were preserved to ensure the same number of cells across the boundary layer, overall mesh resolution and axial spacing in the main regions of interest.
Similarly the cylinder y + was preserved at values less than 1.

The datum configuration
The datum configuration was for a Mach number of 2.88, a Re L of 18.2 × 10 6 with a 16 degree wedge shock generator.The various turbulence models provide a range of results for the cylinder boundary layer at the position ahead of the SBLI.The boundary layer thickness (δ 99 ) and compressible displacement thicknesses (δ*) is assessed at X/d = 9.84 at φ = 0°where the boundary layer measurements showed a δ 99 of 12.2 mm and a δ* of 3.5 mm, with uncertainties of ± 3% and ± 12.4%.The SST and tSST models calculated a δ 99 of 12.7 mm and 12.5 mm, respectively, with the latter being within the 3% experimental uncertainty range.The k-ε and RSM models calculated a slightly thicker boundary layer with δ 99 of + 16.9% and + 7.1%, respectively, relative to the measurements.However, the calculated displacement thicknesses showed a different sensitivity where the k-ε and RSM models provided δ* within the uncertainty level of -0.04 mm (-11.4%) and -0.04 mm (-12.3%),respectively.Conversely, the k-ω SST and γ -θ transition k-ω SST models under-predicted δ* by approximately 30%.
The measurements of the surface static pressure ahead of the interaction indicate that there was a variation in the working section flow field conditions which is greater than the quoted uncertainty of the data acquisition system.Consequently, there are variations in the static pressure distributions on the cylinder ahead of the SBLI region, and therefore to enable a clear comparison of the pressure rise and SBLI extent for the CFD and experimental data sets, the cylinder surface static pressure has been nondimensionalized as shown in Equation ( 3): where P ref is the average surface static pressure in the region ahead of the shock across the area between X/d = 2 and 5.5.
The surface static pressure distributions on the windward (φ = 0°) and leeward (φ = 180°) show notable differences in the characteristics of EV and RSM models (Figure 3).On the windward side, the pressure ratio, P*/P T∞ , shows little difference between EV and RSM models and both are broadly in good agreement with the experimental results.Almost all of the data is within experimental uncertainty bounds although there are some key differences (Figure 3).All models calculated a peak pressure ratio which is larger than the nominal experimental data with the SA, SST, tSST and k-ε models over-predicting the windward peak P*/P T∞ by approximately 8.7%, 11.7%, 9.7%, and 9.6%, respectively.All There are small differences between the CFD models with the SA model exhibiting the rise further downstream at X/d = 9.89 and the RSM model further upstream at X/d = 9.61.In the expansion region downstream of the local maximum there are no notable characteristic differences in the CFD results where all of the P*/P T∞ distributions are slightly higher than the measurements and generally at the upper limit of the measurement uncertainty.The experimental data indicates a small region of local recompression at X/d = 11.5 which, although it is within the measurement uncertainty, is not reflected in any of the computations.
On the leeward side (φ = 180°), however, there are more notable differences between the CFD models as well as between the computational results and the measurements (Figure 3).Although all of the models broadly capture the locus of the initial pressure rise, there are large differences in the values of the local maximum in P*/P T∞ .In particular, peak values from the SA, SST and tSST models are 51%, 37% and 81% higher than the measurements, respectively.Most of the pressure ratio distribution from the k-ε calculations were within the measurement uncertainty and the peak, P*/P T∞ was approximately 12% greater than the measurements.Finally, the best results were obtained using the RSM model where the pressure rise and peak value is within -8% of the measured value.
On the leeward side (Figure 3), the measurements also clearly exhibit a double peak in P*/P T∞ due to a recompression following the initial post-shock expansion with local maxima arising at X/d = 11.8 and 13.2.The SA, k-ε and tSST models do not capture this recompression.The SST and RSM simulations calculate this feature with the RSM getting the correct location, X/d = 12.0 and 13.1, and the SST indicating that the local maxima are slightly further aft at X/d = 12.1 and 13.6.The experiments indicate that the peak re-compression has a local relative magnitude of P*/P T∞ = 0.05.Although the SST and RSM models broadly capture the re-compression, the relative magnitude is around P*/P T∞ = 0.065 and 0.04, respectively, or approximately + 30% and -20% of the experimental values.Overall, based on the distributions of the surface P*/P T∞ ratio, the k-ε and RSM models provide the best agreement with the measurements with perhaps slightly better characteristics for the RSM where the additional modeling fidelity provides some very minor advantages in comparison with the more standard EV models.
Figure 4 illustrates the SBLI system in the windward side symmetry plane for the k-ε calculation using the 3M cell grid.The expected main flow features are evident such as the upstream separation shock, a following expansion Downloaded by [Cranfield University] at 05:31 01 April 2016 region, the primary reflected shock and a local region of boundary layer separation.The other models also resolved these main flow features.For this configuration, there is a boundary layer separation which, from the experimental flow visualizations is estimated to extend from X/d = 10.1 to 10.7 along the cylinder symmetry plane.Both the EV and RSM models successfully calculate the occurrence of the separation region although with varying degrees of success.The SA model calculated a separated boundary layer region with an axial extent of 0.40 X/d which is approximately 31% shorter than that measured (0.58 X/d).
For the other EV models, the k-ε, SST and tSST results showed separation lengths of X/d = 0.50, 0.55 and 0.60 which are approximately − 12%, − 4% and + 5% different from the measurements, respectively.The RSM model, however, provides the best calculation with a separation length of 0.58 X/d which is only -0.9% different from the measurements.This is a remarkably good level of agreement between the CFD and the experimental data given the complexity of the flowfield and is in disagreement with the over prediction of separation bubble size seen by Oliver et al. (2007).For the EV models, however, there is no clear relationship between the fidelity of the peak pressure ratio rise on the windward and leeward side and the accuracy of the calculation of the separation region length.
This SBLI is fundamentally three dimensional and the interaction is affected by the diffraction, and ultimate crossing, of the primary impinging shock around each side of the cylinder.These aspects result in a complex distribution of the surface static pressure on the cylinder which is characterized by a predominantly quasi 2D pressure rise on the cylinder windward side, the sweep of the diffracted impinging shock and the notable subsequent local pressure rises on the leeward side due to the crossing of the respective diffracted shocks.As also shown in the experimental results (Figure 5), the diffraction of the shock around the leeward side results in a region where the surface static pressure is affected by the impinging shock wave and its interaction with the curved cylinder surface.On the windward side the SBLI and pressure rise is affected by the local thickening of the boundary layer as well as the reflected shock.With increasing azimuthal angle the local pressure rise decreases as the shock reflection is reduced although this is a notable increase in the pressure rise again on the leeward side.
Where the shock interacts on the leeward side, it influences the upstream flow conditions (Figure 5), beginning experimentally at X/d = 11.0.By comparing this point of initial X/d Leeward Upstream Influence (LUI) in relation to the X/d position aft of the wedge, the CFD models can be compared.The X/d (LUI) position of the RSM model was effectively congruent with experimental data and within measurement uncertainty margins.The SA, k-ε, SST, and tSST EV models were + 12.8%, 1.3%, 3.1%, and 4.7% respectively.The effect of the transition model was modest.The SA model again illustrated its inaccuracy for such SBLIs, whilst the tSST model was slightly worse than the SST model.
The distribution of p/P T∞ over the cylinder highlights some of the aspects of the flow characteristics as the SBLI develops between the initial windward interaction (φ = 0°) and the eventual behavior on the leeward side.The measured distribution of p/P T∞ shows some of the main features although it is noteworthy that the resolution of the experimental data is φ = 10°steps between φ = −10°to 190°(Figure 5, Figure 6).The spatial distribution along the cylinder surface (φ-X/d) highlights the main region of the strong SBLI pressure rise which is broadly concentrated Downloaded by [Cranfield University] at 05:31 01 April 2016  in the region of X/d = 10 to 12 and φ from 0°to ± 90°.The distributions of p/P T∞ shows the reduction in p/P T∞ in the region aft of the main interaction as well as a region of relatively modest static pressure rise which borders the main shock interaction at around X/d = 11 to 12 and φ = 90°to 135°.On the leeward side and the aft region (X/d = 13 to 14), p/P T∞ increases again.This is where the diffracted shock has propagated around the cylinder body, crosses the centerline and augments the static pressure due to the crossing with the diffracted shock from the opposite side.All of these main features are calculated by both the SST and the RSM models (Figure 5, Figure 6) and overall there is good agreement in the topology of the surface static pressure distributions.In the leeward region there is a modest difference in the p/P T∞ distribution between the SST and RSM models where the SST shows regions with stronger rise and pressure gradients, particularly at φ = 140°and φ = 180°(Figure 5).The RSM model calculates more modest pressure gradients in this region and is in better agreement with the measurements (Figure 5, Figure 6).
As seen in the experiment, the computational results also show a strong lateral pressure gradient near the primary interaction at approximately X/d = 10 to 11 (Figure 5).This generates a strong crossflow and causes the flow to exhibit the wake type characteristics with separation (Brosh & Kussoy, 1983).The measured location point of minimum pressure occurs at φ = 110°and X/d = 11.4 (Figure 5).All of the eddy-viscosity models overpredict this azimuthal position of the minimum pressure point with the SA model having the largest difference φ = 140°( + 27%) although the other models also had notable differences, with k-ε, SST and tSST positions of φ = 128°( + 16%), φ = 132°( + 20%), and φ = 135°( + 23%), respectively.The RSM model provided the best results with the determined azimuthal position as φ = 105°( − 5%).The prediction of the lateral cross-flows was also in better agreement on the subsequent leeward side pressure field, with these distributions and pressure gradients also affecting the local surface flow topology.
As part of the experimental studies, the surface flow topology was identified using a Mylar oil visualization method.This visualization highlighted the primary shock induced separation, which sweeps azimuthally around the cylinder, as well as the secondary separation line, notable crossflow regions and partial evidence of a tertiary separation (Figure 7).
The surface flow topology based on the skin friction distributions for the EV and RSM models are presented in Figure 7.In the experiment the flow visualization indicates a strong cross flow region associated with the initial interaction with the impinging shock which is also observed in all four eddy-viscosity models and best represented by the k-ε model.All EV simulations, the first separation line, S 1 , and node of reattachment, NR 1 , seen in the experiment are represented.The SA model, Figure 7(a), captures the fewest aspects of the surface flow features with a secondary separation line, S 2 , as the only other main surface feature present.The k-ε model, Figure 7(b), calculates a surface flow topology which is very similar to that of the SA model, and also fails to simulate the tertiary separation line.The SST, Figure 7(c), and tSST, Figure 7(d), models show similar primary features as with the SA and k-ε models, but also include a tertiary separation line, S 3 , and are more like the experimental data.This tertiary separation is less clear in the SST model although a secondary attachment node, NR 2 , is present in both models and more clearly resolved in the SST model.The RSM model, Figure 7(e), shows the main features which are an improvement over the kε model.The SST model indicates a tertiary reattachment line (R 3 ) in the region between X/d = 12 to 13.5 although the transition model (Figure 7(d)), shows a slightly different structure without a clearly defined R 3 .There is no clear evidence of such a feature in the flow visualization and overall the indications are that the RSM topology is more closely representative of the experimental data.
Overall it was found that for the RANS method, the turbulence model had a notable impact on the results.The turbulence model affected both the flow topology as well as the quantitative aspects of the SBLI.Of the EV models, the k-ω variants best represented the topology while the k-ε models gave the best quantitative agreement.The SA model performed the worst and demonstrated its unsuitability for this type of SBLI.This conclusion has also been reached in the literature by Oliver et al. (2007), Bhagwandin and DeSpirito (2011) and Tian and Lu (2013).The better agreement for the k-ε and k-ω formulations with experimental data, such as static surface pressures and separation position, have been seen by Oliver et al. (2007) and Bhagwandin and DeSpirito (2011) when comparing two and three-dimensional SBLIs on flat plates and compression ramps.Three-dimensional work by DeBonis et al. (2010) of a shock impinging on a flat plate found that kω models produced the lowest error when compared with SA, k-ε, k-ω and SST methods.The current work supports these findings within the context of shock interactions on three-dimensional geometries in addition to highlighting improvements that are possible using an RSM model when compared with the EV k-ε model.Vallet (2007) reported similar findings but for a simpler compression ramp configuration with varying incidence angles.

Effect of Reynolds number
The datum configuration was for a Reynolds number of 18.2 × 10 6 with a freestream Mach number of 2.88.In the experiment, the effect of Reynolds number was considered for values of 7.28 × 10 6 and 58.3 × 10 6 although in the experiment these also required a slight change in the freestream Mach number of 2.80 and 2. (Table 1).Although this therefore intertwines the effects of Reynolds number, shock strength and shock angle it is still of interest to investigate the changes in the SBLI and the ability of the CFD models to calculate these changes.As these configurations were not the main focus of the previous experimental study, there is a reduced amount of experimental measurements available for these configurations with surface pressure data and flow visualizations for the higher Reynolds number (58.3 × 10 6 ) case.
For the configuration with the Reynolds number increased from the datum value of 18.2 × 10 6 to 58.3 × 10 6 , the Mach number was also increased from 2.88 to 2.95 and therefore for a constant wedge angle of 16°there is a theoretical reduction in the shock angle from 34.2°to 33.6°and the static pressure ratio across the shock changes from 2.89 to 2.94.Overall these changes are relatively minor compared with the datum configuration (Figure 3) and both k-ε and RSM models produce results that are in close agreement with each other although neither captures the sudden drop in P*/P T∞ at X/d = 11.5 (Figure 8).On the leeward side, although there is some variation, the k-ε shows better performance in terms of the peak pressure ratio rise (Figure 8   effect of reducing the Reynolds number from 18.2 × 10 6 to 7.28 × 10 6 was partially investigated in the experimental work and computations of this lower Reynolds number case were also performed using the k-ε and RSM models (Figure 9).The calculations broadly show that there was a slight increase in the strength of the interaction with the peak pressure rise on the windward side from the RSM calculations increasing from around 0.14 to 0.16 when Re L increased from 7.28 × 10 6 to 18.2 × 10 6 .The point of the initial pressure rise is slightly further upstream at the lower Downloaded by [Cranfield University] at 05:31 01 April 2016  Re L , but this is within the uncertainty of the arrangement and the slight changes in the flow conditions.On the windward side, the computations similarly show that the RSM model results in a lower peak pressure rise in comparison with the k-ε results.Concomitant with the slightly higher P*/P T∞ on the windward side, both computations similarly show a slightly higher pressure rise on the leeward side.
The measured distributions of static pressure around the cylinder are shown in Figure 10 and Figure 11 for the high Re L of 58.3 × 10 6 and compared with the results from the k-ε and the RMS models.The point of minimum pressure for the high Reynolds case occurs experimentally at φ = 120°and X/d = 11.5.The EV and RSM models predict the minimum pressure position at X/d = 11.8, φ = 130°and X/d = 11.6,φ = 112°, respectively (Figure 10 and Figure 11).In the low Re L case (Re L = 7.28 × 10 6 ), the position of minimum pressure is found to be X/d = 11.4 and 11.2, φ = 120°and 104°for EV and RSM models respectively.Whilst the axial positions are relatively close, the azimuthal positions vary more notably.In the datum simulations (Figure 5, Figure 6), the RSM model was found to calculate the axial position closer to the experimental results and this is borne out again for the high Reynolds case.
The impact of the change in Re L across the range from 7.28 × 10 6 to 58.3 × 10 6 , for both the experimental flow visualizations as well as the k-ε and RSM computations is shown in Figure 12.The differences in the flow topologies are relatively minor and there is no notable change in either flow visualization of computational results.Nevertheless one of the notable features is in the separation bubble region at φ = 0°where the experiments indicate that the separation bubble is larger as Re L is reduced.This characteristic is also observed in both the k-ε and RSM results (Figure 12(a,b,d,e)).For the RSM model there is a slight change in the topology at the highest Re L (in Figure 12(e)), where there is evidence of a tertiary separation line starting at X/d = 13 and φ = 160°.Overall the numerical methods show similar levels of agreement with the experimental data at the increased Reynolds number configuration and the RSM provides better agreement in the surface flow topology than the k-ε simulations.

Effect of shock strength
The results presented in section 5.1 above, indicate that for the datum configuration (M = 2.88, α = 16°, Re L = 18.2 × 10 6 ) that the k-ε and RSM methods provided the best agreement with the surface static pressures along the windward (φ = 0°) and leeward (φ = 180°) centerlines.In addition, the RSM model shows generally better agreement with the wider surface static pressure distributions and the overall flow topology.Although the experiment considered the impact of the shock strength on the interaction, by performing flow visualizations for a configuration where the wedge angle, α, was increased to 19°, the centerline pressure measurements were not acquired.
Nevertheless, given the confidence in the application of the k-ε and RSM methods to this problem, it is of interest to assess the impact of the increase in the shock strength and change in incidence angle, on the simulated flow fields.The k-ε and RSM methods both show similar changes in the centerline P*/P T∞ (Figure 13 and Figure 14) where the increase in shock strength from 0.14 to 0.20 for the RSM cases is broadly in agreement with the relative theoretical increase of 40% for an inviscid interaction.Of course, the overall levels of P*/ P T∞ are lower than the inviscid value by approximately 30% based on the measurements for the datum configuration.Although the relative changes in P*/ P T∞ on the windward side are broadly as expected, on the leeward side the impact of the increased shock strength is attenuated through the interaction so that the increase in this region is only approximately 13% for both the k-ε and RSM methods.
The increase in the impinging shock strength and change in shock angle has a small impact on the surface flow topology as illustrated in the flow visualizations (Figure 15) but was more difficult to resolve computationally.This has been noted by Oliver et al. (2007) and DeBonis et al. (2010) with Vallet (2007) observing significant failures in the k-ε model's SBLI prediction ability at higher ramp incidence angles in comparison to RSM methods.The flow topology effects are likely due to the increased size of the separation bubble: when α is increased from 16°to 19°, there is a slight increase in the size of the primary separation at φ = 0°and there are minor changes to the loci of the secondary and tertiary separation lines, S 2 and S 3 , respectively.The k-ε model shows some of these characteristic changes, but the structure associated with the primary reattachment line (R 1 ) and the absence of the tertiary separation (S 3 ) indicates that the agreement with the experimental data is perhaps worse than the datum configuration.As with the datum case, the RSM model shows better agreement with the oil visualization also for this case of α = 19°.

Conclusion
A shock-boundary layer interaction test case based on the complex configuration of planar oblique wave impinging onto the curved surface of a prismatic cylinder was computationally investigated.The investigation is based around an established test case under a nominal freestream Mach number of 2.88 and the impact of Reynolds number and shock strength and angle were also considered.The study assessed the performance of a range of computational models and evaluated the results relative to quantitative experimental data as well as flow visualizations.The research considered a range of eddy viscosity RANS models including the Spalart-Allmaras, k-ε Realisable, kω Shear Stress Transport and the k-ω SST model with a γ -θ transition model.In addition a Reynolds Stress Model was also evaluated.A computational approach was developed to establish mesh independence as well as appropriate iterative convergence.Downloaded by [Cranfield University] at 05:31 01 April 2016 The main flow features of such a SBLI were observed; separation bubble, leeward upstream influence, lateral cross-flow, as well as primary and secondary flow separations.Of the EV models tested it was shown that the SA model is unsuitable for complex 3D SBLIs as it did not capture many important flow features.The surface flow topologies were better represented by k-ω variants (SST and the γ -θ SST model) whilst the magnitude was better represented by the k-ε two-equation model.
The RSM seven-equation model improved on the EV models in terms of the flow features such as the separation bubble size, and azimuthal pressure distribution.It replicated the leeward separation lines closer to the experimental results in comparison with the EV models but had difficulty, however, in predicting the leeward focal point.In terms of the magnitude of the local peak pressure values, both the k-ε EV and the RSM models provided the best results, with the RSM providing a better distribution of the flow features.
The sensitivity to changes in Reynolds number, shock strength and incidence were assessed and the RSM computations broadly agreed with the oil flow visualizations.The main aspects of the interaction were found to be relatively insensitive to Reynolds number with a slightly stronger interaction observed at the lowest Reynolds number where the static pressure rise was calculated to be approximately 10% greater.The effect of increasing shock strength showed a slightly larger primary separation region and at the primary impingement location, the static pressure rise increased as expected relative to the changes in the impinging shock.However, on the leeward side, the computations indicate that the interaction attenuates the influence of the impinging shock and the relative pressure rise is notably smaller than that on the windward side.Within the context of multibody shock interactions, this is a notable observation as it indicates that the impact on windward and leeward sides are not simply related and therefore it is expected that the overall resulting body forces and moments will similarly be sensitive to the shock strength.
All of the computational models were tested using a steady-state solver and, given that it is known that such SBLIs typically have low-frequency shock unsteadiness, the steady RANS method may have limitations in predicting this class of SBLIs.However, in spite of this, the computational methods investigated here show that steady RANS methods, specifically k-ε and RSM models, are capable of calculating the overall flow features and quantitative changes associated with an oblique shock impinging on an axisymmetric body.Clearly this conclusion comes with some caveats due to the known assumptions of the steady RANS method, however, given the relative computational cost in comparison with LES or DNS calculations, it indicates that within an engineering context it may be able to play a useful role.The S-A and k-ω variants were found to be unsatisfactory.
Beyond the aforementioned limitations of using steady RANS methods, improvements to the application of EV and RSM models to complex 3D SBLIs could be made by using unsteady-RANS to account for the low-frequency shock unsteadiness experimentally observed and inherent to this type of flow.Furthermore, existing literature suggests that LES and DES models will come to dominate the prediction of this class of flows, however, comparison of results between the different models could yield improvements to RANS methods and thus enable further improvements to this class of model which is typically computationally cheaper and continue their use in applied engineering contexts.

Figure 1 .
Figure 1.Schematic layout of experimental test setup with computational domain shown.

Figure 2 .
Figure 2. Model domain and boundary conditions.

Table 1 .
Brosh and Kussoy (1983)geometric and nominal flow conditions of tested computational models.