## Abstract

The instability of supersonic compression ramp flow is investigated. It is assumed that the Reynolds number is large and that the governing equations are the unsteady triple-deck equations. The mean flow is first calculated by solving the steady equations for various scaled ramp angles *α*, and the numerical results suggest that there is no singularity for increasing ramp angles. The stability of the flow is investigated using two approaches, first by solving the linearized unsteady equations and looking for global modes proportional to *e*^{λt}. In the second approach, the linearized unsteady equations are solved numerically with various initial conditions. Whereas no globally unsteady modes could be found for the range of ramp angles studied, the numerical simulations show the formation of wavepacket type disturbances which grow and convect and reach large amplitudes. However, the numerical results show large variations with grid size even on very fine grids.

## 1. Introduction

The study of interactions of shock wave with boundary layers has a long history and despite many experimental, theoretical and numerical investigations, many fundamental questions remain unanswered in our opinion. One such open question is the precise manner in which instability forms in the separation region in the boundary layer produced by the impinging shock. This instability quickly leads to the breakdown of the laminar flow. A good account of the early experimental findings is given by Chapman *et al*. [1], who found that when the flow was transitional between separation and reattachment, there was a very steep rise in the static wall pressure measurements very close to the transition location. This was very different from what they termed ‘the laminar case’ when the transition was downstream of the reattachment point, and in which the pressure rise was not as marked as in the transitional case. Interaction between a shock wave and a turbulent boundary layer has also received a great deal of attention, see for example the review by Dolling [2]. In [3], it is suggested that the characteristic low-frequency oscillations found in turbulent shock-wave boundary layer interactions is not caused by upstream turbulence, but by the instability of the separation bubble.

The main objective of this work is to investigate the stability of the supersonic flow past a compression ramp in which the governing equations are taken to be the triple-deck equations. Smith [4] in his seminal paper showed how triple-deck theory could be used to describe Tollmien–Schlichting instability in subsonic boundary layers. In supersonic flow, the same linearized triple-deck equations, but with a different interaction law, do not give rise to unstable Tollmien–Schlichting modes.

For supersonic compression ramp flow, numerical solution of the steady equations shows that a recirculation region develops in the vicinity of the corner. For very large ramp angles, secondary separation is also observed [5]. In previous work, Cassel *et al.* [6], when solving the unsteady triple-deck equations for the same problem, noted that an instability forms above a critical ramp angle *α*>3.9 with a wavepacket developing close to the corner, growing in time, but remaining stationary near the corner. It was suggested that this instability is associated with the formation of inflection points of the velocity profile in the recirculation region. On the other hand, a related study by Fletcher *et al.* [7] suggested that the compression ramp flow is convectively unstable for angles *α* in the range 3.2≤*α*<3.7. For larger angles, Fletcher *et al*. [7] suggested that an absolute instability may be responsible for the abrupt local breakdown of the boundary layer. In both these studies, the starting point was the integration of the full unsteady triple-deck equations from a given initial condition until either the steady solution was reached, or an instability was encountered.

One common difficulty which has been encountered by many researchers studying this problem is that the results obtained, and hence the conclusions generated, change from study to study (e.g. [6,7]). This may be because adequately resolved grids have not been used. In addition, the work of Tutty & Cowley [8] would suggest that there may be long wavelength Rayleigh instabilities in the interactive boundary layer formulations which would also complicate grid refinement studies in reaching firm conclusions.

In this work, one major difference from previous studies is that we first compute the basic flow using the steady-state triple-deck equations and then investigate the stability of this basic flow independently. A direct hybrid finite difference and spectral method similar to that described in [5] is used for the computation of the basic flow, and the technique is robust and can compute the flow accurately for large ramp angles. Two distinct approaches are used to investigate the instability of the computed flow. Firstly, a global linear stability analysis, taking disturbances proportional to e^{−λt}, is performed. With this approach, a generalized eigenvalue problem needs to be solved and our results indicate that no unstable eigenvalues λ exist for a wide range of ramp angles, including those giving rise to secondary separation. Secondly, numerical simulations have been carried out with the linearized unsteady equations (linearized about the steady state) using forced disturbances. The simulations suggest that instabilities are present but difficulties were encountered in obtaining grid independent results.

This paper is organized as follows. In §2, we give a brief mathematical description of the problem and discuss the numerical techniques used. In §3, we give the main results of our work, and finally in §4 we give our conclusions.

## 2. Governing equations and numerical techniques

We consider the flow of an ideal compressible polytropic gas past a compression ramp. A supersonic flow approaches the corner with freestream speed , density and pressure . The Reynolds number *Re* and Mach number are defined by
where is the viscosity in the freestream, and *L* is some characteristic length scale (for example, the length of the plate before the start of the corner). We will consider the limit when *Re*≫1 and the Mach number is *O*(1).

A shock wave is generated in the vicinity of the corner and gives rise to an interaction with the boundary layer. A detailed analysis of the interaction based on triple-deck theory is omitted but may be found in many places (e.g. [9]). The interaction is described by the unsteady triple-deck equations given by
2.1
2.2
2.3
and the boundary conditions
2.4
Here, *u*(*x*,*y*,*t*) and *v*(*x*,*y*,*t*) are the scaled non-dimensional velocity components in the lower deck, (*x*,*y*) are the corresponding coordinates, *t* time, *p*(*x*,*t*) the pressure, *A*(*x*,*t*) the displacement function, *f*(*x*) the scaled wall shape, and we have used the Prandtl transposition theorem.

The wall shape used in our work was a smooth corner given by
and *r* is a rounding parameter taken be 0.5.

### (a) Discretization of the equations: basic flow

To calculate the basic flow with ∂*u*/∂*t*=0 in (2.2), equations (2.1)–(2.4) are discretized on a truncated domain using a hybrid finite-difference spectral collocation method as in [5]. Two methods were used for computing the basic flow. The first is as described in [5] and referred to as Method B in that paper. In the second approach, we worked with primitive variables rather than with streamfunction vorticity variables. In the *x*-direction, second-order finite differences are used. In addition, because of the sharp gradients which form in the separated flow regions, we used a clustering of points in regions of high shear stress given by the transformation
where *a* and *b* are two parameters which help control the location of the clustering and density of points in the region. Typically, we used *b*=5, and *a* was chosen to coincide with the position of minimal skin friction from a previous solution. A uniform finite-difference grid in the *X* direction is used so that

In the *y*-direction, the domain was truncated to , and we used Chebychev collocation with the collocation points defined by
Various values were tried and for the results reported the value of was used. Derivatives in *y* may then be calculated by using the Chebychev differentiation matrices [10]. It is noted here that additional care needs to be exercised in computing the higher derivative matrices accurately as, for example, simply squaring the first derivative matrix to obtain the second derivative matrix leads to loss of precision which can be significant for the unsteady problem.

The continuity equation was handled by using
which in the regions of unseparated flow may be discretized as
2.5
where we have defined *u*_{q,j}=*u*(*x*_{q},*y*_{j}) and *v*_{q,j}=*v*(*x*_{q},*y*_{j}). In (2.5), the *M*_{q,j} define the entries of the matrix **M** which is a discrete representation of the integral operator in the continuity equation, and *β*(*x*)=d*X*/d*x*.

In the regions of secondary separation and reattachment, the differencing was adjusted to include a combination of forward and central differences as
with *η*=0 in regions of secondary separation and reattachment.

In a similar way, the *x* derivatives were discretized by using a three-point backward difference in regions of separated flow so that ∂*u*/∂*x* is approximated by
In the separated flow region, this was replaced by a three-point forward difference
The discretization leads to a set of nonlinear discrete equations. Newton linearization was used to solve the nonlinear system and the resulting linear system is of block pentadiagonal form
with suitable adjustments to the matrices at the endpoints *q*=1,2,*m*−1,*m*. Here, *Φ*_{q} is defined as the vector of unknowns at each *x*=*x*_{q} location and
Further details of the matrices and discretization are given in [11]. The linear system was solved directly.

For small values of *α*, the solution can be obtained from an arbitrary initial guess. For larger values of *α*, a previous solution is used to provide the initial guess.

### (b) Global stability techniques

Having obtained a basic flow (which we denote by a suffix B), the instability of the flow was studied using two different methods. In the first approach, we write the perturbations to the flow as
2.6
with similar expressions for the other quantities. Expression (2.6) were substituted into (2.1)–(2.4) which were then linearized for small *ϵ* leading to the unsteady linear perturbation equations
2.7
2.8
2.9
and homogeneous boundary conditions.

The linear system (2.7)–(2.9) was discretized in exactly the same way as for the basic flow. This leads to the solution of the generalized eigenvalue problem
2.10
where **J** is Jacobian matrix of the linear system in the computation of the basic flow, and **H** is a singular diagonal matrix. The generalized eigenvalue problem was solved using either Matlab, or with routines from the ARPACK library.

### (c) Unsteady simulations

Perturbations to the basic flow were also studied by conducting full simulations of the perturbed equations. The perturbations may be written as 2.11 with similar expressions for the other quantities. The linearized equations are similar to (2.7)–(2.9) but with the momentum equation replaced by The same spatial discretization was used. For the time differencing, we used a fully implicit scheme with a second-order backward difference for the time derivative term.

Forced perturbations were introduced into the flow by replacing one of the no-slip conditions by
where *δ*,*μ*,*ξ* and *t*_{d} are parameters affecting the amplitude, location, spread and duration of the disturbance, respectively. Typical values used in the simulations were *δ*=0.001, *μ*=50 000, *t*_{d}=0.015 and *ξ*=0.5.

## 3. Results

### (a) Basic flow

The techniques described in the previous section are fairly robust and allowed the basic flow to be calculated to high angles where secondary separation is present. The use of primitive variables and the handling of the continuity equations in the manner described above is novel. The steady flow problem has been solved by many people but there are relatively few results for large ramp angles. Smith & Khorrami [12] suggest that the solutions cannot be obtained beyond some critical ramp angle because of a singularity which develops as this critical angle is approached. Our results compare well, see also [11], with those of [5] for angles up to *α*=5, although there are slight differences for *α*=7.5, which may be due to the differences in the grid sizes used. For instance, for *α*=4.5, we find separation at *x*_{sep}≈−13.67 and the downstream reattachment at *x*_{reat}≈=9.245. This is consistent with the data presented in [5], whereas in Smith & Khorrami [12] the corresponding values are *x*_{sep}≈−5 and *x*_{reat}≈4. However, for *α*=7.5, the results in [5] give the separation point to be *x*_{sep}≈−52, with reattachment at *x*_{reat}≈18 and the minimal skin friction at *x*_{ms}≈12. Comparison of these values with figure 1 shows that our data disagree with these values. We calculated *x*_{sep}≈−39.34,*x*_{reat}≈12.67 and *x*_{ms}≈8.45. In Korolev *et al.* [5], there are no results shown between *α*=5 and *α*=7.5. The good agreement with [5], at least for angles up to *α*=5, using a different method, would tend to confirm the suggestion in [5] that the [12] results are incorrect. The study in [5] also suggests that there is no singularity for some critical ramp angle, and in fact the solutions can be continued indefinitely, a conclusion that is also consistent with our findings.

Plots of the basic pressure and skin friction are shown in figure 1. These show the pressure developing a long plateau region and the presence of secondary separation with a large drop in skin friction downstream of the corner.

### (b) Global stability results

The mean flow computed earlier was next used in the investigations of the global stability of the flow. The matrices generated are huge and sparse of size (*N*+3)*m*×(*N*+3)*m* with *N*=80, *m*=3001 typically. Thus, not all the spectrum can be evaluated and only the dominant and least stable eigenvalues are computed. Various checks were made for grid independence of the results. Results for the dominant eigenvalue are shown in table 1. Plots of the real and imaginary parts of the spectrum for four different angles are shown in figures 2 and 3. These plots show that there is one dominant eigenvalue close to the origin on the real axis but with real part still positive. Calculations were performed for angles up to *α*=7.75 but no unstable eigenvalues were found! This result is somewhat surprising and unexpected and would suggest that for this range of angles, there is no global linear instability, with the eigenfunction bounded in the spatial direction, with a e^{−λt} type of dependence. Of course, this does not preclude other types of instability modes or global instability at much larger ramp angles. The numerical results of Cassel *et al*. [6] and Fletcher *et al*. [7] suggest that there are instabilities present for all the angles studied in figures 2 and 3.

### (c) Linear stability computations using forced perturbations

A number of computations were carried out using forced perturbations. In figure 4, we show plots of the perturbation skin friction for an angle of *α*=3.6. The results show a number of interesting features. First, convectively unstable disturbances are present. An initial wavepacket forms, grows and is convected downstream. However, secondary wavepackets also form at later times. These enter the separated flow region and increase in amplitude near the point of minimal skin friction. For larger ramp angles, even with infinitesimal small forcing, the wavepackets can reach very large amplitudes.

Various grid size checks were carried out. In figure 5, we show the skin friction for various values of *N* the number of Chebychev points, keeping *m* the number of points in the streamwise direction fixed. These show grid convergence, provided a fine grid is used for small times. In figure 5*b*, it is seen that with a coarse grid, the disturbances reach large amplitudes much quicker, but taking a finer grid suppresses this instability until later times. For large times, say *t*=0.4 in figure 5*d*, there are variations in the region where a secondary wavepacket develops. This is shown more clearly in 6*a*, but it can also be noted that the primary wavepacket features are fully resolved. The same conclusion is reached by keeping a fine grid in *y* with *N*=160 and varying the number of streamwise points (figure 6*b*). The secondary wavepacket features are completely unresolved, even though very fine grids have been used. It is found that a continuous stream of wavepackets form for *α*>3.2.

In order to assess whether the grid-size dependencies observed were due to the numerical methods adopted, some computations were performed in which the no-slip condition was replaced by one appropriate for a moving wall with
replaced by
Apart from the very slight change in boundary conditions, the operator is essentially identical to the non-moving case. The basic flow skin friction is shown in figure 7 for *α*=3.6 for a moving and a non-moving wall. The slight rise in the skin friction for a non-moving wall when *x*<0 disappears and the drop in skin friction after the corner is more pronounced in the moving wall case. Global stability computations for the moving wall case (not presented here) also showed no unstable eigenmodes for the same range of angles. In figure 8, we show the evolution of the disturbance skin friction (mean flow subtracted). As time increases, the disturbance dies down as suggested by the global stability analysis. The same is true for larger angles (figure 8*b*). Grid size checks as shown in figure 9 show fully resolved calculations for the moving wall case. The instabilities found earlier have disappeared.

Finally, additional simulations were carried out for the non-moving wall case in which the linearized equations were rewritten to solve for the total flow and with zero forcing, taking the basic flow as the initial condition. For values of *α*>3.2, we noted that the time signals for the disturbances were similar to the secondary wavepackets seen before for the forced simulations. As there was no forcing, these disturbances arise from perturbations owing to round-off errors.

## 4. Conclusion

A novel method using primitive variables has been successfully used to compute accurate solutions of the steady triple-deck equations for the supersonic compression ramp. Results have been obtained for large angles and are consistent with those presented in [5] for angles up to *α*=5, but there are slight differences for *α*=7.5.

Our unsteady results have shown that despite using a number of different robust numerical techniques, the precise mechanisms causing the breakdown of the unsteady supersonic compression ramp flow, governed by the triple-deck equations, still remain unclear. On the one hand, global stability analysis of the flow suggests that there are no unstable eigenfunctions of the form for a wide range of ramp angles (0<*α*<7.8) studied here. However, simulations of forced (and unforced) disturbances with the linearized disturbance equations do show evidence of some instabilities present for angles *α*>3.2. Difficulties have been experienced in obtaining grid independent and self-consistent results for a number of angles, including *α*=3.6, as disturbances are excited by round-off error. With coarse grids, the disturbances reach large amplitudes much earlier.

The numerical problem as posed does not indicate why there should be so much difficulty in computing the solution. Mathematically, the problem seems fairly innocuous but numerically it appears to be particularly challenging.

Tutty & Cowley [8] suggest that Rayleigh's inflection point criterion is a necessary condition for instability of the triple-deck equations for several interaction laws including the one studied here. This would explain some of the things we have found. However, as also noted in [6], inflection points of the basic profile only occur for larger angles. The profile with *α*=3.6 is free of inflection points and thus instabilities observed in computations here with *α*=3.6 cannot be directly linked with the Tutty & Cowley [8] long-wave Rayleigh instability.

For a moving wall, the evolution of disturbances is consistent with that predicted by the global stability analysis, and disturbances die down eventually.

Having conducted many simulations, it is unclear to us whether this is a physical or numerical instability that we have observed. Our observations are qualitatively similar to [7] but not [6]. The same equations with a slight change in boundary conditions for the moving wall problem show no instabilities. Simulations for other flows governed by the triple-deck equations, such as jet flows, also show similar unexplained behaviour in regions where there are no globally unstable modes [11]. Future work may be usefully directed at examining the non-normality of the operators involved and calculating the pseudo-spectrum to see whether this is able to account for results presented here.

## Acknowledgements

We are grateful to Dr A. Hazel and Prof. K. Cassel for helpful discussions on this work.

## Footnotes

One contribution of 15 to a Theme Issue ‘Stability, separation and close body interactions’.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.