Gauging Centrifugal Instabilities in Compressible Free-Shear Layers Via Nonlinear Boundary Region Equations

Curved free shear layers emerge in many engineering problems involving complex flow geometries, such as the flow over a backward facing step, flows with wall injection in a boundary layer, the flow inside side-dump combustors, or wakes generated by vertical axis wind turbines, among others. Previous studies involving centrifugal instabilities have mainly focused on wall-flows where Taylor instabilities between two rotating concentric cylinders or G\"{o}rtler vortices in boundary layers, resulting from the imbalance between the centrifugal forces and the radial pressure gradients, are generated. Curved free shear layer flows, however, have not received sufficient attention, especially in the nonlinear regime. The present work investigates the development of centrifugal instabilities evolving in a curved free shear layer flow in the nonlinear compressible regime. The compressible Navier-Stokes equations are reduced to the nonlinear boundary region equations (BRE) in a high Reynolds number asymptotic framework wherein the streamwise wavelengths of the disturbances are assumed to be much larger than the spanwise and wall-normal counterparts. We study the effect of the freestream Mach number $\boldsymbol{M_\infty}$, the shear layer thickness $\boldsymbol{\delta}$, the amplitude of the incoming disturbance $\boldsymbol{A}$, and the relative velocity difference across the shear layer $\boldsymbol{\Delta V}$ on the development of these centrifugal instabilities.


Introduction
The stability of curved shear layer flows depends on the velocity difference across the shear layer and the radius of the curvature.For a free shear layer with no curvature, as known as a plane shear layer, the Kelvin-Helmholtz instability is the dominant turbulent mechanism [8].In this case, the two-dimensional disturbances are more unstable than their three-dimensional counterpart [15], and the Kelvin-Helmholtz instability introduces predominantly spanwise oriented vortices.Rayleigh [12] proved that the presence of an inflection point in the basic velocity profile is a necessary condition for the Kelvin-Helmholtz instability.On the other hand, for a curved mixing layer flow with an inflectional velocity profile, the Kelvin-Helmholtz instability mechanism is still present, along with centrifugal instabilities in the form of streamwise oriented Görtler type vortices.
Görtler vortices are mostly known to appear inside a boundary layer flow along a concave surface due to the imbalance between radial pressure gradients and centrifugal forces (e.g., Gortler [4], Hall [5], Swearingen & Blackwelder [16]).For highly curved walls, for example, vortex formation occurs more rapidly and can significantly alter the mean flow causing the laminar flow to transition into turbulence.Under certain conditions, Görtler vortices can be efficient precursors to transition.The growth rate of these counter-rotating streamwise vortical structures depends on the surface curvature and the receptivity of the boundary layer to freestream disturbances and surface imperfections.This type of instability is predominant in many engineering applications.Thus, understanding the processes leading to their development and predicting their occurrence using efficient and tractable methods could potentially advance the overall understanding of turbulence.
Numerous theoretical and numerical studies covered centrifugal instabilities in incompressible curved free shear layer flows.Plesniak et al. [10,11] conducted extensive experimental measurements investigating curved two-stream mixing layers to show how centrifugal effects yield streamwise vortices.The untripped case within this suite of experiments exhibited organized streamwise vorticity, while the tripped case did not.The researchers explain that this is due to the spatially stationary streamwise vortices that provide extra entrainment to the flow in the tripped case (Bell & Mehta [1] showed this for a plane mixing layer).Hu et al. [6] and Liou [7] focused on the effect of the curvature on the inflectional Rayleigh modes, which they found to be minimal, although the curvature excites an unstable three-dimensional disturbance with the amplitude increasing as the streamwise wavenumber decreases.The analytical and numerical study of Otto et al. [9] showed that the unstable modes largely depend on surface curvature.They also employed numerical simulations to solve the parabolic equations, assuming that the wavenumber and Görtler number are both of order one.They found that as the difference between the freestream speeds increased, the layer became more susceptible to centrifugal instabilities.
In the present work, we analyze the development of centrifugal instabilities in highspeed compressible curved free shear layer flows via an efficient numerical algorithm based on the nonlinear boundary region equations (NBREs) -a parabolized version of the Navier-Stokes equations under the assumption that the streamwise wavenumber associated with the disturbances is much smaller than the cross-flow wavenumbers.
The study considers the effects of a wide range of Mach numbers, the amplitude of the freestream disturbance, A, the shear layer thickness, δ, and the velocity difference across the shear layer, ∆V , on the development and growth of these centrifugal instabilities.The study shows, among other things, that the kinetic energy level of the curved shear layer flow is directly proportional to ∆V and A while it is inversely proportional to δ. Increasing A induces larger instability structures, as expected, which may be beneficial for enhancing mixing.The location of the maximum energy moves farther downstream as the freestream Mach number increases.

Scalings
All dimensional spatial coordinates (x * , y * , z * ) are normalized by the spanwise separation λ * , while the dependent variables by their respective freestream values, except the pressure, which is normalized by the dynamic pressure: where λ * is the spanwise wavelength of the disturbances, (u * , v * , w * ) are the velocity components, ρ * the density, p * is pressure, T * temperature, µ * dynamic viscosity, k * thermal conductivity, and all quantities with ∞ at the subscript represent conditions at infinity.Reynolds number based on the spanwise separation, Mach number and Prandtl number are defined as where µ * ∞ , a * ∞ and k * ∞ are freestream dynamic viscosity, speed of sound and thermal conductivity, respectively, and C p is the specific heat at constant pressure.As for boundary layer flows over curved surfaces, we here define the equivalent global Görtler number as where r * is the radius of the curvature.

Boundary region equations: a parabolized form of the Navier-Stokes equations
If the streamwise wavenumber of the disturbances evolving inside the shear layer are much larger that the wavenumbers corresponding to the crossflow directions, then the Navier-Stokes equations can be transformed into a parabolic set of equations in the framework of high Reynolds number asymptotics.For a full compressible, Newtonian flow, the primitive form of the Navier-Stokes equations with non-dimensional variables are considered here in the form is the substantial derivative.The pressure p, the temperature T and the density ρ of the fluid are combined in the equation of state in non-dimensional form, p = ρ T /γM 2 ∞ , assuming that non-chemically-reacting flows are considered.Other notations include the dynamic viscosity µ, and the free-stream Mach number The dynamic viscosity µ and thermal conductivity k are linked to the temperature using a power law in dimensionless form, where b = 0.76 (Ricco & Wu [13]), C p = γR/(γ − 1), γ = 1.4,and P r = 0.72 for air.We re-scale the streamwise distance and time co-ordinate at which the vortex system forms by the following O(1) variables: x = x/R λ , and the time as t = t/R λ .Note that the distance in the wall-normal and spanwise directions are the same, y = ȳ, z = z.Another thing to mention is that, in this region, the crossflow velocity component is small compared to the streamwise velocity component, and pressure variations are negligible.Appropriate dominant balance considerations suggest that the dependent variables in this region must also re-scale as follows: Working out the order-of-magnitude analysis of the Navier-Stokes equations, we obtain the parabolic set of equations, which we refer to as the nonlinear compressible boundary region equations (NCBRE) where ⃗ V is the velocity vector and ∇ c is the crossflow nabla operator: The effect of the wall curvature is contained in the term involving the global Görtler number G λ in the second momentum equation.
A small artificial disturbance is imposed at the inflow boundary in the form: where A is a small amplitude (in this study A = 0.04).λ * is the spanwise wavenumber (dictating the spanwise separation of the centrifugal instabilities), and σ represents the extent of the disturbance in the y direction.In the present work, λ is kept constant at 0.8 cm.Es-Sahli et al. [2] elaborately studied the effect of λ * on the development of centrifugal instabilities in curved free shear layers.The NCBREs are solved using the algorithm developed in Es-Sahli et al. [3] and previously in Sescu and Thompson [14].We slightly adjusted the algorithm to accommodate the free shear layer setting.

Results
In this section, we present and discuss a set of results from the curved free shear layer numerical simulations.The flow domain is split into 'fast' and 'slow' streams, both having velocities V f and V s , respectively, which we define as; where ∆V is the relative velocity difference (here, we set this difference at four levels, 20%, 30%, 40%, and 50%).The temperature in the fast stream is set to T ∞ , while in the slow stream is set to 0.9T ∞ .In our parametric study, we consider three Mach numbers in the fast stream of 2.0, 4.0, and 6.0, respectively (the lowest velocity in the slow stream will correspond to a Mach number of 1.0).We also consider three values for the shear layer thickness δ at the inflow boundary at 0.2, 0.4, and 0.6, where the velocity variation between the fast and slow streams is modeled via a hyperbolic tangent function 0.5(1−tanh (y − y 0 )/δ), whith y 0 representing the location of the shear layer.A similar function is used to model the variation of the temperature in the shear layer, with the same thickness, although in reality the thickness of the thermal layer may be slightly different.The Reynolds number R λ based on the faster freestream and the spanwise separation of the disturbance and the global Görtler number G λ are kept the same for all cases at 10 6 and 2 × 10 5 , respectively (the kinematic viscosity and the curvature of the wall were varied to achieve constant R λ and G λ for all simulations).The grid is uniform in the spanwise direction taking into account that the flow is periodic in this direction, while in the radial direction the grid is stretched towards the top and bottom far-field boundaries.The marching in the streamwise direction is achieved by means of an explixit Euler method, with equally-spaced discretization.At the inflow boundary, centrifugal instabilities are excited by the non-dimensional artificial disturbance in equation 20 imposed on the base flow at the inflow boundary, with the amplitude set at 0.04.
Figures 1-3 present consecutive contour plots depicting the magnitude of crossflow velocity for different Mach numbers and velocity differences ∆V , with the shear layer thickness δ set to 0.2.The color scheme represents the flow velocity, with white at the top indicating the fast stream and black at the bottom indicating the slow stream.These contour plots illustrate the progression of the centrifugal instabilites in the streamwise direction; we hypothesize that the incoming disturbance is sufficiently high to allow the centrifugal instabilities to take over the Kelvin-Helmholtz type instabilities (this mathematical model is not capable of predicting Kelvin-Helmholtz type instabilities because the problem is steady and the system of equations is parabolic).To characterize the shape of the centrifugal instabilities, we identify two types of structures: "primary" and "secondary."The primary flow structure refers to the mushroom-like formation that evolves as the main instability, as exemplified by the fourth panel of the first row in figure 1.On the other hand, the secondary structures are elongated features that emerge from the edges of the primary flow structure, as seen, for example, in the fourth panel of the second row in figure 2.
By comparing the first to the second row of each figure 1, 2 or 3, we notice that by increasing the velocity difference across the shear layer, ∆V , accelerates the development of the mushroom-like primary structures and makes the secondary structures more prominent.The color transition from black to white seen in the bottom rows of each of these figures suggest that the higher the velocity difference across the shear layer, the more intense the mixing in the shear layer is.Comparing the centrifugal instabilities at the same streamwise coordinate across various Mach numbers, we observe a delay in the growth of the mushroom-like structures as M ∞ increases.We also observe that the flow structures become thinner as the Mach number increases.For instance, in the third panel of the first row in figure 1, for which M ∞ = 2, the mushroom structure has already begun to form, whereas in the corresponding panel of figure 3, where M ∞ = 6, the mushroom shape has not yetr developed.This finding suggests that, for high Mach number free shear layer flows, the same mixing efficiency is achieved further downstream.
To quantify the thermal effects of these instabilities, vertical profiles of temperature disturbance T d (x, y, z) = T (x, y, z) − T m (x, y), where T m (x, y) is the spanwise mean component of temperature T (x, y, z), through the center of the mushroom shape are included in figures 4-6.The profiles are compared to each other for different shear layer thicknesses and different velocity difference levels across the shear layer.They all show that increasing the thickness of the shear layer increases the temperature disturbance, and that this increase is slightly less significant at high Mach numbers.As expected, by increasing the velocity difference across the shear layer increases the amplitude levels of the temperature disturbance at all Mach number.Also, the vertical extent of these disturbances seem to increase with increasing the velocity difference, especially towards the slow stream (for example, in figure 4, the top boundary of the disturbance is roughly in y = 1 for all velocity differences, while the bottom boundary is roughly in y = 0.9 for ∆V = 80% and in y = 1.2 for ∆V = 50%).where u m (x, y), v m (x, y), and w m (x, y) are the spanwise mean components of velocity, and z 1 and z 2 are the coordinates of the boundaries in the spanwise direction (note that the energy here is also scaled by ∆V 2 , which implies that the energy variation is the result of disturbances developing in the shear layer).In figures 7-9, we plot E(x) for different values of ∆V and in figure 10 we plot the same energy for different values of δ.The disturbance energy E(x) seems to be directly proportional to ∆V as it is highest for ∆V = 50% and decreases as ∆V is reduced for all considered cases.Moreover, the streamwise location of the energy saturation (the point at which the energy starts to level off) moves farther downstream as ∆V decreases, especially for higher Mach numbers (see figure 9).On the other hand, as seen in figure 10, E(x) is inversely proportional to δ, which may be because viscous effects are more predominant in thicker shear layers; in this figure, we superposed the smallest (solid line) and the    highest (dashed line) velocity difference levels, indicating that the the trend is the same for both (the other two velocity difference levels fall in between).We must point out that the energy reduction due to the increase in the shear layer thickness is rather significant, perhaps because the variation of the shear layer thickness is substantial (th e largest δ is three times greater that the smallest δ).

Conclusions
In this paper, we investigate the nonlinear development of centrifugal instabilities in a compressible curved free shear layer flow using a numerical solution to the boundary region equations, specifically, a parabolized form of the Navier-Stokes equations.Our focus lies on understanding the characteristics of these centrifugal instabilities, which exhibit similarities with the development of G"ortler vortices in boundary layer flows over concave surfaces.The study encompasses variations in the free stream Mach number (M ∞ ), the relative velocity difference between the two streams of the shear layer (∆V ), the shear layer thickness (δ), and the amplitude of the inflow disturbance (A).Upon closer examination of the kinetic energy plots for the M ∞ = 2 case in figure 7, we observe that E(x) is directly proportional to ∆V and inversely proportional to δ across all the considered Mach numbers.However, the increase in shear layer thickness has an insignificant effect on energy reduction, with a mere 1% drop in E resulting from a 100% increase in δ, as evident in figure 10.Furthermore, increasing the amplitude of the inflow disturbance (A) slightly boosts the kinetic energy, with less than 1% increase in E observed for a 100% increase in A. Interestingly, a larger magnitude of A hampers the influence of the relative velocity difference.The energy curves corresponding to different ∆V values exhibit a considerably reduced gap when comparing the A = 0.02 and A = 0.04 cases.Similar trends were observed in the parametric study of centrifugal instability development for the M ∞ = 4 and M ∞ = 6 cases.
Examining the contour plots of crossflow velocity magnitude in figures 1-3 for ∆V = 30%, we find that increasing the disturbance amplitude leads to significant growth in the mushroom-like structure's amplitude and renders the secondary structures more visible, indicating increased mixing for all Mach numbers under consideration.When comparing different Mach numbers at the same crossflow plane location, we observe a slower development of mushroom-like structures as M ∞ increases, suggesting that achieving the same mixing efficiency would require traveling further downstream for higher Mach numbers.Consequently, numerical simulations for higher Mach number cases would incur higher computational costs due to the need for larger grid sizes to maintain comparable mixing efficiency.

Fig. 4
Fig. 4 Profiles of temperature disturbance T d (y) for the M = 2 case.

Fig. 5
Fig. 5 Profiles of temperature disturbance T d (y) for the M = 4 case.

Fig. 6
Fig. 6 Profiles of temperature disturbance T d (y) for the M = 6 case.

Fig. 7
Fig. 7 Vortex energy distribution of different parametric settings for the M = 2 case.

Fig. 8
Fig. 8 Vortex energy distribution of different parametric settings for the M = 4 case.

Fig. 9
Fig. 9 Vortex energy distribution of different parametric settings for the M = 6 case.