Numerical Simulation of Supersonic Jet Impingement on Inclined Plate

Supersonic jet impinging on an inclined plate is explored numerically by solving Reynolds Averaged Navier Stokes (RANS) equations using a commercial solver. The numerical uncertainty is quantified by grid convergence index parameter analysis. Simulations capture crisply all the essential features of the flow field including jet boundary, jet shock, upper tail shock, front and behind plate shocks. Good agreement between computational and experimental results for different plate angles forms the basis of further analysis. Four different two equations turbulence models predict similar jet structures. Effect of plate angles, pressure ratio, nozzle-plate distance on impinged jet structures are evaluated from numerical results. It has been observed that with increase in the plate angle, the maximum pressure zone in the plate changed from crescent shape to round shape. Well resolved RANS simulations are capable of capturing finer details of supersonic impinging jets on inclined plate.

Study of supersonic jets impinging on an inclined flat plate is important for scientific investigation as well as practical applications. The problem of jet impingement in inclined plate for aerospace applications appears in the design of jet deflector, multistage rocket separation at higher attitude, rocket test-stand environment, plume ducting system of canisterised missiles, space module attitude-control thruster operation, among others. Such flow field contains many complex fluid dynamics phenomena like shock/shock interactions and shock/boundary layer interactions. The jet structures change greatly depending on the inclined angle of the plate, jet pressure ratios and the nozzle-plate distance. It has been observed that under certain flow conditions, localized pressure peaks damaged the impingement plate. Inclined jet impingement exhibits more complex features than the perpendicular jet impingement1,2. The maximum pressure on the inclined plate can be several times larger than that on the perpendicular plate because of the complex shock-shock interactions present. Moreover, the stagnation bubble may disappear when the plate angle decreases. Indeed, it is a great challenge to understand the physics of these flows and correctly predict the heat and pressure loads on the impinging plate.

Most of the former studies on supersonic jet impingement deal with experimental investigation of perpendicular impingement on a flat plate1-4. These studies, explained the existence of stagnation bubbles that appears in the vicinity of the plate surface, affecting the pressure peaks on the plate surface and the stability of the jet. One of the earlier important work on supersonic jet impingement on inclined plate is the experimental study by Lamont and Hunt 5. From the pressure measurement on the plate surface and shadowgraph visualization it has become clear that the plate inclination has a strong influence on the pressure distribution. Kim and Chang6 and Wu et al 7,8 simulated numerically the experimental study of Lamout and Hunt by solving Euler and Navier Stokes equations respectively and obtained qualitative agreement with experimental results. These numerical studies could provide explanations of the complex shock interactions. Nakai et al 9,10 conducted an important experimental investigation of supersonic jet impingement on inclined plate with different inclination angles, pressure ratios and nozzle plate distances. They have measured detailed pressure distribution in the plate surface with Pressure Sensitive Paints (PSP) and used schlieren methods to visualize the flow field. Based on the schlieren images, different shock wave structures were classified into three major types. According to this classification, the Type I flowfield is mostly observed when the plate is nearly perpendicular to the jet, the nozzle-plate distance is large, and the pressure ratio is low. Three shock waves (upper tail shock in the upstream area, lower tail shock in the downstream area and bowl shaped plate shock in the middle area) appear over the plate in addition to the jet shock at the jet boundary. As the plate angle decreases (θ<500) the wall jet expands along the plate surface and generates an intermediate tail shock wave emanating from the triple point generated at the end of the plate shock wave and Type II flowfield occurs. A convex plate shock region is formed where the plate shock and barrel shock interact in the upstream region and a “stagnation bubble” region appears. As the plate angle decreases further (θ<300), the intermediate tail shock wave merges with the barrel shock wave and give rise to Type III flow field. Furthermore, Nakai et al 9 showed that the properties of the flow fileds, such as the number of the pressure peaks and the existence of circulation area etc., can be forecast roughly by the inclined angle of the plate, and the nozzle-plate distance regardless of the pressure ratio.

Because of high computational requirements, The CFD studies of jet impingements were limited to perpendicular impingement with low pressure ratios. McIlroy and Fujii11 have numerically simulated the experimental condition of Nakai et al 9,10 for pressure ratios of 4.5, 7.4 and 10.1 and presented a detailed flow structure of this complex problem. Three dimensional Navier Stokes equations with Baldwin Lomax turbulence model12 were solved by discretizing the inviscid terms with a second order accurate Simple High-resolution Upwind Scheme (SHUS)13 and viscous terms with a second order central differencing. It was observed from their studies that stagnation point of the main stream, strong shock waves in the upstream area, reattachment of the detached flow and interaction between the intermediate tail shock and the boundary layer are some of the important factors for the localized pressure peaks on the plate surface. The flow fields were classified into six types by the presence of these pressure peaks. Yoshinori et al14 also simulated the experimental conditions of Nakai et. al.9,10 by employing both Reynolds Averaged Navier Stokes (RANS) simulation and Implicit Large Eddy Simulation (ILES). The RANS simulations for this case are similar to that employed by McIlroy and Fujii 11; whereas higher (7th) order accurate shock capturing Weighted Compact Nonliner Scheme (WCNS)15 (a combination of WENO and Compact schemes) for inviscid term and 6th order central differencing for viscous term are employed for ILES studies without any explicit LES sub-grid scale model. Although, a great detail of the flow structures are explained, the comparison of both RANS and ILES with experimental results are only qualitative.It is clear from the above discussion that adequate studies are not presented in the literature to understand the flow feature of impingement of supersonic jet in an inclined plate and to demonstrate the quantitative prediction capability of the CFD codes for this complex fluid dynamic problem.

In the present work, the experimental conditions of Nakai et. al.9,10 are simulated numerically using commercial CFD software16 for different pressure ratios, inclination angles and the nozzle exit to impingement plate distance. The computed flow profiles are compared with the experimental value to demonstrate the capability of CFD codes with industry standard physical models.

Schematic of experimental condition for which the simulations are carried out is shown in Fig 1. The experiments are carried out in a small induction-type wind tunnel. Supersonic cold jet from a large low-pressure chamber passed through Mach 2.2 conical nozzle and impinged on an inclined plate placed in the test section. The total chamber pressure of the flow is kept constant at 100 kPa and the test section is connected to a vacuum chamber with pressure level of 10 kPa. Ambient pressure near the flat plate is varied to obtain the necessary pressure ratio of the nozzle exit and ambient

Figure 1.(a) Schematic picture of nozzle with flat plate (b) Covergent divergent nozzle. Schematic of experimental condition9,10.

pressure. The throat diameter (Dt) and the nozzle exit diameter (Dj) are 5 mm and 7.08 mm respectively. The location and the angle of the plate in the test section are controlled electrically to give the position L/Dj from 1.0 to 4.5 and the inclination angle (θ) in the range of 30° to 60°. (θ=90° corresponds to the case of the flat plate being perpendicular to the jet flow). A Thin Layer Chromatography (TLC) plate is glued over the inclined plate to measure the pressure and temperature along the plate using pressure and temperature sensitive paints. Flow visualization with the Schlieren system is also carried out in the experiment. Further details of the experiment are available in Nakai9-10.

3.1 Methodology:

Commercial CFD software, CFX16 is used for the simulation. It solves 3-D Reynolds Averaged Navier Stokes (RANS) equation on structured grid based on finite volume approach. It also solves one of the following turbulence models viz. k-$ϵ$ , k-$\omega$ or SST turbulence model, etc along with the RANS equations. The software has four major modules

a) ICEM-CFD, imports CAD geometry or creates geometry and generates unstructured volume meshing based on the user input

b) preprocessor - sets up the boundary condition and initial field condition

c) solver manager - solves the flow field based on the grid and the boundary condition and

d) postprocessor - visualizes and extracts the results

3.1.1 Governing Equations

The appropriate system of equations governing the turbulent flow of a compressible gas may be written as:

Continuity equation:

$\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial {x}_{k}}\left(\rho \text{\hspace{0.17em}}{u}_{k}\right)=0$ k = 1,2,3

Momentum equation:

$\frac{\partial }{\partial t}\left(\rho \text{\hspace{0.17em}}{u}_{i}\right)+\frac{\partial }{\partial {x}_{k}}\left(\rho \text{\hspace{0.17em}}\text{\hspace{0.17em}}{u}_{i}\text{\hspace{0.17em}}{u}_{k}\right)+\frac{\partial P}{\partial {x}_{i}}=\frac{\partial \left({\tau }_{ik}\right)}{\partial {x}_{k}}$, i,k = 1,2,3

Energy equation:

$\frac{\partial }{\partial \text{\hspace{0.17em}}t}\left(\rho \text{\hspace{0.17em}}\text{\hspace{0.17em}}E\right)+\frac{\partial }{\partial \text{\hspace{0.17em}}{x}_{k}}\left(\rho \text{\hspace{0.17em}}\text{\hspace{0.17em}}{u}_{k}H\right)=-\frac{\partial }{\partial \text{\hspace{0.17em}}{x}_{k}}\left({u}_{j}\text{\hspace{0.17em}}{\tau }_{jk}\right)+\frac{\partial {q}_{k}}{\partial \text{\hspace{0.17em}}{x}_{k}}$, j,k = 1,2,3

where,$\rho$ , ui, p, E and H are the density, velocity components, pressure and total energy and total enthalpy respectively. In eddy viscosity models, the stress tensor is expressed as a function of turbulent viscosity ( ${\mu }_{t}$ ).Turbulent kinetic energy (k), Turbulent dissipation rate($ϵ$ ) and Specific dissipation rate($\omega$ ) are defined as:

Turbulent kinetic energy k,

$k=\overline{{{u}^{\prime }}_{i}{{u}^{\prime }}_{i}}/2$

Turbulent dissipation rate $\epsilon$ ,
$\epsilon \equiv \nu \overline{\frac{\partial {{u}^{\prime }}_{i}}{\partial {x}_{j}}\left(\frac{\partial {{u}^{\prime }}_{i}}{\partial {x}_{j}}+\frac{\partial {{u}^{\prime }}_{j}}{\partial {x}_{i}}\right)}$

Specific dissipation rate $\omega$ ,
$\omega \approx \frac{\epsilon }{k}$

Turbulent kinetic energy (k) equation:

$\frac{\partial }{\partial t}\left(\rho \text{\hspace{0.17em}}\text{​}k\right)+\frac{\partial }{\partial {x}_{k}}\left(\rho \text{\hspace{0.17em}}{u}_{k}k\right)=\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\partial }{\partial {x}_{k}}\left(\left(\frac{{\mu }_{l}}{Pr}+\frac{{\mu }_{t}}{{\sigma }_{K}}\right)\text{\hspace{0.17em}}\frac{\partial k}{\partial {x}_{k}}\right)+{S}_{k}$

Rate of dissipation of turbulent kinetic energy ($ϵ$ ) equation:

$\frac{\partial }{\partial t}\left(\rho \text{\hspace{0.17em}}\text{​}\epsilon \right)+\frac{\partial }{\partial {x}_{k}}\left(\rho \text{\hspace{0.17em}}{u}_{k}\epsilon \right)=\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\partial }{\partial {x}_{k}}\left(\left(\frac{{\mu }_{l}}{Pr}+\frac{{\mu }_{t}}{{\sigma }_{\epsilon }}\right)\text{\hspace{0.17em}}\frac{\partial \epsilon }{\partial {x}_{k}}\right)+{S}_{\epsilon }$

and $\mu$ = ${\mu }_{l}$ + ${\mu }_{t}$ is the total viscosity; ${\mu }_{l}$ , ${\mu }_{t}$ being the laminar and turbulent viscosity and Pr is the Prandtl number. The source terms Sk and Se of the K and e equation are defined as

${S}_{K}={\tau }_{ik}\frac{\partial {u}_{i}}{\partial {x}_{k}}-\rho \text{\hspace{0.17em}}\epsilon$ and    ${S}_{\epsilon }=\text{\hspace{0.17em}}\text{\hspace{0.17em}}C{}_{\epsilon 1}\tau {}_{ik}\frac{\partial {u}_{i}}{\partial {x}_{k}}-{C}_{\epsilon 2}\frac{\rho \text{\hspace{0.17em}}{\epsilon }^{2}}{K}$

where turbulent shear stress is defined as

${\tau }_{ik}={\mu }_{t}\left(\frac{\partial {u}_{i}}{\partial {x}_{k}}+\frac{\partial {u}_{k}}{\partial {x}_{i}}\right)$

Laminar viscosity (μl) is calculated from Sutherland law as

Where T is the temperature and μref, Tref and S are known values.

The turbulent viscosity μt is calculated as

${\mu }_{t}={c}_{\mu }\frac{\rho \text{\hspace{0.17em}}{k}^{2}}{\epsilon }$

The coefficients involved in the calculation of μ are taken as
cμ =0.09,         C$ϵ$ 1=1.44,                    C$ϵ$2 = 1.92
${\sigma }_{k}$ =1.0,          ${\sigma }_{3}$ =1.3,                       ${\sigma }_{c}$ = 0.9

The heat flux qk is calculated as  ${\theta }_{\kappa }=\text{​}\text{\hspace{0.17em}}-\lambda \text{\hspace{0.17em}}\frac{\partial Τ}{\partial {\xi }_{\kappa }}$ , $\lambda$ is the thermal conductivity

3.1.2 k-ω Turbulence model

The turbulent viscosity is calculated as function of k and ω17.

${\mu }_{t}=f\left(\frac{\rho k}{\omega }\right)$
Turbulent kinetic energy (k) equation:

$\frac{\partial }{\partial t}\left(\rho k\right)+\frac{\partial }{\partial {x}_{i}}\left(\rho k{u}_{i}\right)=\frac{\partial }{\partial {x}_{j}}\left({\Gamma }_{k}\frac{\partial k}{\partial {x}_{j}}\right)+{G}_{k}-{Y}_{k}$
Specific dissipation rate (ω) equation:

$\frac{\partial }{\partial t}\left(\rho \omega \right)+\frac{\partial }{\partial {x}_{i}}\left(\rho \omega {u}_{i}\right)=\frac{\partial }{\partial {x}_{j}}\left({\Gamma }_{\omega }\frac{\partial \omega }{\partial {x}_{j}}\right)+{G}_{\omega }-{Y}_{\omega }$

Where Gk is turbulence production due to viscous force and similar to Sk, ,${\Gamma }_{k}=\mu +\frac{{\mu }_{t}}{{\sigma }_{k}}$ ,${G}_{w}=\alpha \frac{\omega }{k}{G}_{K}$ ,    and ${\Gamma }_{w}=\mu +\frac{{\mu }_{t}}{{\sigma }_{w}}$ of the k and ω equations respectively. Where ${\beta }^{1}=0.09$ , $\alpha =5/9,$ ,$\beta =0.075$ and ${\sigma }_{k}={\sigma }_{w}=2$  are the model constants.
The turbulent viscosity ${\mu }_{t}$ is calculated as ${\mu }_{t}={c}_{\mu }\frac{\rho \text{\hspace{0.17em}}{k}^{}}{\omega }$

3.1.3 SST Turbulence Model

To retain the robust and accurate formulation of Wilcox’s k- ω model in the near wall region, and to take advantage of the freestream independence of the k- ε model in the outer part of the boundary layer, Menter 18 blended both the models through a switching function. k- ε model was transformed into Wilcox’s k- ω formulation and was multiplied by (1-F1) and added to original k- ω model multiplied by F1 . The blending function F1 will be one in the near wall region and zero away from the surface. In the second step, the definition of eddy viscosity was modified in the following way to account for the transport of the principal turbulent shear stress

( $\tau =\overline{-\rho {u}^{\prime }{v}^{\prime }}$ )
${\nu }_{t}=\frac{{a}_{1}k}{\mathrm{max}\left({a}_{1}\omega ;\Omega {F}_{2}\right)}$

3.1.4 Renormalized group (RNG) k-ε model:

An improved method for rapidly strained flows based on rigorous statistical technique 19 is also used in the present calculations. The k and ε equations are given below,

where Gk is the generation of turbulence kinetic energy due to the mean velocity gradients, calculated as ${G}_{k}=-\rho \overline{{{u}^{\prime }}_{i}{{u}^{\prime }}_{j}}\frac{\partial {u}_{j}}{\partial {x}_{i}}$     and YM represents compressibility effects given by,

.

The turbulent Mach number Mt is given by ${\text{M}}_{t}=\sqrt{k/{a}^{2}}$  ,
The model constants are taken as,

The additional term in ε equation Rε is given as
, ${\eta }_{0}=4.38$   and $\beta =0.012$ .

3.2 Discretisation of Governing Equations

The solver utilizes a finite volume approach, in which the conservation equations in differential form are integrated over a control volume described around a node, to obtain an integral equation. The pressure integral terms in the momentum integral equation and the spatial derivative terms in the integral equations are evaluated using finite element approach. An element is described with eight neighboring nodes. The advective term is evaluated using upwind differencing with physical advection correction. The set of discretised equations form a set of algebraic equations:

$\stackrel{}{A}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\stackrel{\to }{x}\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}b$  where $\stackrel{\to }{x}$  is the solution vector. The solver uses an iterative procedure to update an approximated ${x}_{n}$  (solution of x at n­th time level) by solving for an approximate correction x'  from the equation $A\text{\hspace{0.17em}}\text{\hspace{0.17em}}\stackrel{\to }{{x}^{\prime }}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\stackrel{\to }{R}$ , where $\stackrel{\to }{R}\text{\hspace{0.17em}}=\stackrel{\to }{b}\text{\hspace{0.17em}}-\text{\hspace{0.17em}}\text{\hspace{0.17em}}A\text{\hspace{0.17em}}{\stackrel{\to }{x}}_{n}$  is the residual at nth time level. The equation $A\text{\hspace{0.17em}}\text{\hspace{0.17em}}\stackrel{\to }{{x}^{\prime }}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\stackrel{\to }{R}$  is solved approximately using an approach called Incomplete Lower Upper factorization method. An algebraic multigrid method is implemented to reduce low frequency errors in the solution of the algebraic equations. Maximum residual (= ${\phi }_{j}^{n+1}-f\left({\phi }_{j}^{n+1},\text{\hspace{0.17em}}{\phi }_{j}^{n}\right)$ ) < 10-4 is taken as convergence criteria.

3.3 Computational Grid and inflow parameters.

The schematic of computational domain is shown in Fig 2.  It consists of CD nozzle and an inclined plate. To capture all the essential features of the flow field, the computational domain is extented upto 150Dj, 75Dj, 100Dj and 75Dj (Dj is the nozzle exit diameter) in the longitudinal direction, upper, lower and lateral directions boundary respectively.  In the simulation, X-axis is taken along the longitudinal direction, whereas Y and Z axis are taken along the height and width as shown in the

Figure 2.Computational domain showing computational boundary and co-ordinate directions.

Fig. 2. The origin is taken at the center of the nozzle exit plane. A multiblock structured grid consisting of 5.3 million grids is generated in the computational domain. The grid distribution in the X-Z and Y-Z plane is shown in Fig. 3. A very fine grid (y+ ~ 5) is employed near the plate surface and nozzle exit plane. Number of numerical simulations are carried out by considering different plate angles (θ), pressure ratios (PR = pj/p where pj and p are nozzle exit and freestream pressure respectively) and the distance between the nozzle exit and plate (L/Dj). The inflow parameters and simulation matrix are shown in Table 1. In the nozzle inlet plane, stagnation temperature (298 K) and stagnation pressures (100 kPa) are specified and the atmospheric conditions are specified at far field, outflow and the rest of the inflow domain. The simulations are carried out for cold flow condition.  In the simulation, the plate angle (θ), pressure ratio (pj / p) and flat plate location (L/Dj) are varied in the range 30° - 60°, 1.2 - 10.0 and 1.0 – 3.0 respectively.

Figure 3.Computational grid in the jet vane flow field (a) X-Y plane (b) Y-Z plane at nozzle exit.

Table 1. Inflow parameters and simulation matrix.

Taking the advantage of the symmetry, only one half of the computational domain is simulated. The grid independence of the results is demonstrated by comparing the radial surface pressure distribution on the plate  between coarse grid (2.6 Million) and fine grid (5.2 Million) for plate inclination (θ) 30°,  Pressure ratio (pj /p) of 7 and plate location  (L/Dj) of 2.0 in Fig. 4. The distance is nondimensionalized by nozzle exit radius (Rn), whereas the pressure is nondimensionalized by the nozzle chamber pressure (Pc).A very good comparison of surface pressure between the two grids is obtained. The quantification of error due to grid in the form of Grid Convergence Index (GCI) is also presented in the figure.

Figure 4.Comparison of surface pressure with different grids.

The main source of numerical error in CFD for steady state boundary value problem is iterative convergence or grid convergence error.  This error can be estimated by carrying out the simulation in two different grids namely; the coarse grid and the fine grid.  The simplest of such estimate is given by the relative difference  $\epsilon$ = (f2-f1)/f120, where f represent any quantity of interest and the indices 2 and 1 refer to the fine and coarse grid solution respectively (In the present calculation, surface pressure is taken as the parameter of interest). Roache21 has proposed a grid-convergence index (GCI) as an error based on uncertainty estimate of the numerical solution as

$GCI={F}_{s}\frac{1\epsilon 1}{{\left({h}_{2}/{h}_{1}\right)}^{p}-1}$

where h is the order of grid spacing, p is the order of accuracy of numerical scheme and Fs is a factor of safety.  Roache 22 has suggested Fs = 3 for minimal of two grid calculations.  For the present calculation p is equal 2 with h2/h1 equal to 2, GCI is order of ε.  Maximum errors between two grid are within 5 %.  This analysis indicates that the grid is adequate to capture most of features of the flow and the solution in grid independent.  The comparison with the experimental data revealed that the simulation captured the plume impingement region in the plate very crisply. In the downstream region, the numerical results show sharp shock reflections compared to the experimental data. It is not clear whether the measurements points are fine enough to capture these shock reflections accurately. Mach number distribution and pressure distribution around the plate are shown in Figs. 5 and 6 to depict the complex three dimensional structure of the flow field.  Mach number distribution is presented in the symmetry plane in Figure 5. Pressure Isosurface near nozzle exit with streamlines on flat plate colors with pressure ratio is shown in Fig 6(a) while a composite picture of pressure ratio contours on the inclined plate and density contour on the plane perpendicular to the plate is presented in Fig. 6(b). The jet boundary, jet shock, upper tail shock, front and behind plate shocks are captured crisply.

Figure 5.Mach number distribution at symmetry plane for θ = 30° and pressure ratio ( pj /p∞) = 7.0.

Figure 6.Three dimensional structure of the flow field. (a) Pressure isosurface near nozzle exit with streamlines on flat plate colors with pressure ratio (b) composite picture of pressure ratio contours on the inclined plate and density contour on the plane perpendicular to the plate.

The numerical and experimental Schlieren for different plate inclinations are presented in Fig. 7. All the finer features of the under expanded jet (pj /p = 7.0) structure as seen in the experimental results  are  captured by the numerical Schlieren. The expansion waves generated from the under expanded nozzle get reflected from the jet boundary and transform into compression waves which coalesce to form barrel shape jet shock. Mach disc is formed very close to the plate. With the increase in the plate angle, the flow type changed from Type

Figure 7.Comparison of schlieren picture (Density gradient) for pj /p = 7 and L/Dj = 2.0 for different plate inclinations (θ = 30°, 40° 50° and 60°).

III to Type I as explained by Nakai9,10 et al. The crescent shape pressure peak for lower inclination angle got changed into round shaped pressure peak.  The maximum error of location of the Mach disc between the computation and experiment is less than 10 %. Jet structures obtained from the RANS simulations are comparable with Large Eddy Simulation results11 available in the literature.

The pressure contour on the plate for different plate angles are shown in Fig. 8. The location and pattern of impingements are changing marginally with the change of plate angles. Computed surface pressures for different plate inclinations (θ = 30°, 40°, 50°, and 60°) are compared with the experimental data in Fig. 9. The computed pressure peaks and features within the first shock cell match very well with the experimental result. The shock reflections downstream of the mach discs present in the numerical results are not observed in the experiment. To resolve the difference, a new simulation has been carried out for θ = 40° (Fig. 9(b)) with very fine grid of 7.62 million size compared to the 4.07 million grid considered earlier. The new grid is very fine in the downstream of the impingement point.

Figure 8.Pressure distributions on the plate with different plate angles.

Figure 9.Surface pressure comparison for pj /p = 7 and L/Dj = 2.0 with different plate angles (a) θ =30°, (b) θ = 40°, (c) θ = 50°, (d) θ = 60°.

The comparison of surface pressure with two different grids is shown in Fig. 9(b). It could be observed that the result did not change with this fine grid. It may be noted that Grid Convergence Index (GCI) distribution shown in Fig. 4 is very small which indicate the adequacy of the grid distribution for this problem.  The cause of the difference between experiment and numerical results particularly in the downstream of the impingement point could not be explained. Four different turbulence models namely; k-ε, k-ω, RNG k-ε and SST were studied to assess the predictive capabilities of the models to simulate the plume impingement flow field in inclined plate. The computed surface pressures with different turbulence models for θ = 30° are shown in Fig. 10. Performances of all the turbulence models are similar. Parametric simulations are also carried out to study the effect of the nozzle exit to free stream pressure ratio (pj /p) on the plate surface pressure. The comparison of the surface pressure for three different pressure ratios 4, 7, 10 for θ = 45° are shown in Fig 11 The structure of the jet remains almost the same. With more expansion for higher ratios, the peak pressure reduces. Also, the pressure levels downstream of the Mach disc are lower for higher pressure ratios. Simulations are also carried out to with different

Figure 10.Surface pressure for pj /p= 7, L/Dj = 2.0, and θ = 30° with different turbulence models.

Figure 11. Surface pressure for θ = 45°, L/Dj = 2.0 with different pressure ratios.

positioning of the plate at L/Dj = 1.0, 2.0, and 3.0 to find out the effect of the position of the plate on the surface pressure. The pressure field for plate inclination (θ) 45° and pressure ratio (PR) 8.4 for different plate locations are shown in Fig 12. As expected, with the increase of distance from nozzle, the surface pressure in the impingement zone reduces significantly. The variation of computed surface pressure along the plate length is presented in Fig 13. Peak pressure reduces to four times for L/Dj = 2.0 compared to that of 1.0. Also, the sizes of the shock cells change for different plate distances. The peak surface pressure for different plate locations are presented in Fig 14. The reduction in the peak pressure for the near location of plates is more compared to the distant location as found from the changing slope of the curve.

Numerical simulations are carried out to simulate supersonic jet impingement on an inclined plate.

Figure 12.Pressure field for (a) L/Dj = 1, (b) L/Dj = 2 and (c) L/Dj = 3 for θ = 45°, pj /p= 8.4.

Figure 13.Comparisons of Surface pressure for different locations of plate (L/Dj) from nozzle exit.

Figure 14.Variation of peak pressure ratio with different L/Dj from nozzle exit.

Three dimensional RANS equations are solved along with two equation turbulence model using commercial CFD code. Grid independence of the solution is demonstrated and the numerical errors are quantified through Grid Convergence Index parameter. Simulations capture all the features of the flow fields including jet boundary, jet shock, upper tail shock, front and behind plate shocks. Jet structures obtained from the RANS equations are comparable with Large Eddy simulation results available in the literature. Computed pressure on the plate and numerical Schlieren matches very well with experimental results for different plate angles. The simulation predicts smaller shock cells downstream of the mach disc not captured by the experiment.  All the four different two equations turbulence models employed in the study predict the same jet structure. The crescent shape maximum pressure zone at the plate for lower inclination angle changes into round shape for higher inclination angles. With increase in the pressure ratios and nozzle-plate distance, the peak pressure and shock reflections reduces significantly. It has been demonstrated that well resolved RANS simulations are adequate to capture the finer details of impinging jet in an inclined plate.

1. Ginzburg, I.P.; Semilentenco, B.G.; Terpigorer, V.S. & Uskov, V.N. Some singularities of supersonic underexpanded jet interaction with a plane obstacle. J. Eng. Phys.,1970,19(3), 1081-1084. [Full text via CrossRef]

2. Kalghatgi, G.T. & Hunt, B.L. Experiments on the impingement of a supersonic jet on a flat plate, Aeronaut. Q., 1976, 27(2), 169-185.

3. Gummer, J.H. & Hunt, B.L. The impingement of a uniform, axisymmetric, supersonic jet on a perpendicular flat plate.  Aeronaut. Q., 1971, 12(4), 403-420.

4. Pamadi, B.L. On the impingement of supersonic jet on a normal flat surface. Aeronaut. Q., 1982, 33(3), 199-218.

5. Lamont, P.J. & Hunt, B.L. The impingement of underexpanded, axisymmetric jets on perpendicular and inclined flat plates. J. Fluid Mech., 1980, 100(3), 471-511. [Full text via CrossRef]

6. Kim, K.H., & Chang, K.S., Three-dimensional dtructure of a dupersonic jet impinging on an inclined plate. J. Space. Rockets, 1994, 31(5), 778–782. [Full text via CrossRef]

7. Wu, J.; Tang, L.; Tong, X.; Luke, E. & Cinnella, P. Numerical simulations of jet impingement and jet separation. Dept. of Aerospace Engineering, Mississippi State Univ., Rept.MSSU-ASE-00-2, Mississippi State, Aug. 2000.

8. Wu, J.; Tang, L.; Luke, E. A.; Tong, X.-L. & Cinnella, P. Comprehensive numerical study of jet-flow impingement over flat plates. J. Space. Rockets, 2002, 39(3), 357-366. [Full text via CrossRef]

9. Nakai, Y.; Fujimatsu, N. & Fujii, K. Flow classification of the under-expanded super sonic jet impinging on a flat plate. AIAA Paper 2003-3467. doi: 10.2514/6.2003-3467

10. Nakai, Y.; Fujimatsu, N. & Fujii, K. Experimental study of underexpanded supersonic jet impingement on an inclined flat plate.  AIAA Journal, 2006, 44(11), 2691-2699. [Full text via CrossRef]

11. McIlroy, K. & Fujii, K. Computational analysis of supersonic underexpanded jets impinging on an inclined flat plate.  AIAA Paper 2007-3859. [Full text via CrossRef]

12. Baldwin, B.S. & Lomax, H. Thin layer approximation and algebraic model for separated turbulent flows.  AIAA Paper 78-257. [Full text via CrossRef]

13. Shima, E. & Jounouchi, T. Role of CFD in aeronautical engineering (No. 14) -AUSM type upwind schemes. In the Proceedings of 14th NAL Symposium on Aircraft Computational Aerodynamics, NAL SP-34, 1996. pp.7-12.

14. Yoshinori, G.; Nonomura, T.& Kansai, M. Detailed analysis of flat plate pressure peaks created by supersonic jet impingements. AIAA Paper 2009-1289. [Full text via CrossRef]

15. Deng, X.G. & Zhang, H. Developing high-order weighted compact nonlinear schemes. J. Comput. Phys., 2000, 165(1), 22-44. [Full text via CrossRef]

16. ANSYS CFX, Release 11.0: Installation and Overview, January, 2007.

17. Wilcox, D.C. Multiscale model for turbulent flows.  AIAA Journal, Year 26(11), 1311-1320. [Full text via CrossRef]

18. Yakhot, V., & Orszag, S.A. Renormalization Group Analysis of turbulence. J. Sci. Comput., 1986, 1(1), pp.3-57. [Full text via CrossRef]

19. Menter, F.R. Performance of popular turbulence models for attached and separated adverse pressure gradient flows. AIAA Journal, Year 30(8), 2066-2072. [Full text via CrossRef]

20. Cutler, A.D.; Danehy, P.M.; O’Byrne, S.; Rodrigues C.G. & Drummond, J.P. Supersonic combustion experiment for CFD model development and validation. AIAA Paper 2004-266. [Full text via CrossRef]

21. Roache, P.J. Error bars for CFD. AIAA Paper No. 2003-0408. [Full text via CrossRef]

22. Roache, P.J. Verification and validation in computational Science and EngineeringHermon  Publishers, New Mexico, 1998.

 Mr Malsur Dharavath obtained his ME (Aerospace Engineering) from Indian Institute of Science (IISc), Bengaluru. He is working in the Defence Research and Development Laboratory (DRDL), Hyderabad. He has about 6 journal papers and 7 conference papers to his credit. His research ares includes: Scramjet propulsion system, combustion modeling, free and confined supersonic jets applicable to missile propulsion. Dr Debasis Chakraborty obtained his PhD (Aerospace Engineering) from Indian Institute of Science (IISc), Bengaluru. Presently, he is working as Technology Director, Computational Dynamics (CFD), DRDL, Hyderabad. He was the guest editor of special issue of Defence Science Journal on CFD. His research areas includes: CFD, aerodynamics, high speed combustion and propulsion. He has about 48 journal papers and 60 conference papers to his credit.