Numerical Analysis of Laminar‐Turbulent Bifurcation Scenarios in Kelvin‐Helmholtz and Rayleigh‐Taylor Instabilities for Compressible Flow

Abstract

In the chapter, we are focused on laminar-turbulent transition in compressible flows triggered by Kelvin-Helmholtz (KH) and Rayleigh-Taylor (RT) instabilities. Compressible flow equations in conservation form are considered. We bring forth the characteristic feature of supersonic flow from the dynamical system point of view. Namely, we show analytically and confirm numerically that the phase space is separated into independent subspaces by the systems of stationary shock waves. Floquet theory analysis is applied to the linearized problem using matrix-free implicitly restarted Arnoldi method. All numerical methods are designed for CPU and multiGPU architecture using MPI across GPUs. Some benchmark data and features of development are presented. We show that KH for symmetric 2D perturbations undergoes cycle bifurcation scenarios with many chaotic cycle threads, each thread being a Feigenbaum-Sharkovskiy-Magnitskii (FShM) cascade. With the break of the symmetry, a 3D instability develops rapidly, and the bifurcations includes Landau-Hopf scenario with computationally stable 4D torus. For each torus, there exist threads of cycles that can develop chaotic regimes, so the flow is more complicated and difficult to study. Thus, we present laminar-turbulent development of compressible RT and KH instabilities as the bifurcations scenarios.

1. Introduction and definitions

1.1. Introduction

In the chapter, we continue the work started in 2006 with the numerical analysis of laminar‐turbulent transition using dynamical system approach. The results of the previous work are partially collected in [1]. One must point out that the bifurcation analysis of laminar‐turbulent transition becomes a very interesting subject of research that opens a new way of looking at the problem of turbulence. This can be traced by the growing amount of publications in this field [2]. From mechanical point of view, one must study not only secondary flow (that bifurcates from the primary flow) but all other secondary flows and ways of their formation until they lead to the “turbulent” flow. At this point, the definition of “turbulence” must be introduced. It usually changes from paper to paper since there is no strict definition. For example, statistical definition is usually given in engineering papers and attractor‐related definition is given in mathematical papers. We stick to the latter definition and introduce definitions after a short review of the progress in this field.

1.2. Problem overview

One of first ideas of bifurcation formation of complex flow structures originates from the 2D periodic forced flow problem by A.N. Kolmogorov [3], formulated in 1958 in the “Selected topics in analysis” seminar at Lomonosov Moscow State University. In 1960, the problem of the stability for the primary flow was solved by Meshalkin and Sinai [4] using chain fractions. It was stated that if the aspect ratio of the domain $\alpha =Y/X<1$, then the flow is unstable. It was later confirmed [5, 6] that the flow is indeed unstable and forms oscillatory solutions. Analysis of laminar‐turbulent transition was conducted later in Refs. [7, 8]. It was found that the system undergoes supercritical pitchfork bifurcation with further formation of limited cycle and torus through Hopf bifurcation. The flow complication continuous on invariant cycles (or tori if the symmetry is not exploited) cascades with Feigenbaum‐Sharkovkskiy sequence up to the torus period three with sequential emergence of the chaos. Another classical problem is the Rayleigh‐Benard convection problem. The linear stability was analyzed in 1916 by Lord Rayleigh himself [9]. The emergence of secondary stationary and oscillatory flows was later confirmed by many authors, e.g., [10–12]. Last two references contain bifurcation scenarios for 2D and 3D rectangular domains for laminar‐turbulent transition (“frozen turbulence” in 2D case since stretching of vortex tubes is locked). See detailed monograph [13] for more information. Bifurcations, as laminar‐turbulent transition mechanism, are also considered for ABC flow [14]. The appearance of Landau‐Hopf bifurcations up to the stable invariant torus of co‐dimension three with Feigenbaum bifurcations on two‐dimensional invariant tori was observed. It is known that the formation of the turbulence occurs through bifurcations in Cuette‐Taylor flow [15–18] and experimental bifurcation mode analysis in Ref. [19]. It is shown that the flow undergoes either subcritical or supercritical pitchfork bifurcation with further flow development through series of limited cycles or tori up to co‐dimension two. Analytical research of the Cuette‐Taylor problem was conducted in Ref. [18]. Nonlinear theory of secondary flow is constructed and its stability is analyzed. Bifurcation analysis of the flow over a backward‐facing step problem is presented in Refs. [20–22]. One shows the appearance of stability loss due to the Taylor‐Gortler instability (corresponds to the Hopf bifurcation) with further formation of the stable invariant three‐dimensional torus. We were able to trace the formation of 3D torus of period two before chaos in [22]. It was also found that the bifurcation scenario of limited circles and limited two‐dimensional and three‐dimensional tori is an initial stage of transition to chaos in 2D Poiseuille flow [23, 24].

In this chapter, we concentrate on compressible flow features and instabilities triggered by Kelvin‐Helmholtz (KH) and Rayleigh‐Taylor (RT) mechanisms. Rayleigh‐Taylor instability is the instability of an interface when a lighter fluid (gas) supports a heavier one in a gravitational field or some external potential. Direction of the potential gradient vector is from the heavier fluid to the lighter one. It can also occur when a lighter fluid pushes a heavier one. It was first studied experimentally by Lord Rayleigh [25] and later, theoretically, by Taylor [26]. Kelvin‐Helmholtz instability is the instability of a shear layer (a tangential discontinuity for inviscid fluid) that can occur when there is velocity shear in a single continuous fluid or if there is a velocity difference across the interface between two fluids. It was initially studied by Helmholtz [27] and Lord Kelvin [28]. Those two instabilities are often considered together. Indeed, RT instability causes movement of adjusted fluids in different directions with the appearance of the shear layer that is subject to KH instability.

The problem of RT and KH instabilities was first attacked by linear stability analysis. We give brief information about it, following [29–31] and a highly recommended book [32]. The first stability conditions for the RT and KH instability were derived from the problem of interface stability of inviscid potential flow of two fluids. The setup is as follows: interface is at $y=0$, flow properties are for $y>0:u_1,{\rho }_1$ and for $y<0:u_2,{\rho }_2$. Gravity vector is $g=(\mathrm{0,}-g{\mathrm{,0})}^T$, directed downward. The analysis gives the following stability condition:

$k^2{\rho }_1{\rho }_2{(u_1-u_2)}^2<kg({\rho }_2^2-{\rho }_1^2),$E1

where $k$ is a wavenumber. It is clear, that if $u_1\ne u_2$ then there exists such $k=k_{crit}$, that the flow is unstable and KH instability always occurs (including cases with ${\rho }_1={\rho }_2$ and $g=0$). If $u_1=u_2$, then RT instability occurs if ${\rho }_1>{\rho }_2$, e.g., the lighter fluid supports the heavier one. This result is obtained from very simple problem formulation that is merely adequate. In most cases of linear analysis, Squire’s theorem [32–34] is applied in order to drive the problem of linear analysis of a flat‐parallel 3D flow to a 2D flow analysis. The first improvement of the result (Eq. (1)) for KH is the consideration of boundary condition problem instead of initial‐condition one (e.g., flow of two co‐directional fluids with the shear). In this case, one can derive the critical wavelength ($\lambda =2\pi /k$, ${\lambda }_{crit}=2\pi /k_{crit}$) for instability [35] to be:

${\lambda }_{crit}=\frac{\pi {\rho }_0{(u_1-u_2)}^2}{g({\rho }_2-{\rho }_1)}{\left[1+\frac{u_1+u_2}{\sqrt{2(u_1^2+u_2^2)}}\right]}^{-1}$E2

It is clear that in this case, the critical wavelength will be longer than for the initial value problem case. The necessary instability condition for the stratified shear flow instability is derived from the energy balance equation and consideration of continuous shear velocity and density profiles [32]. Then, from the stream‐function formulation, one derives Taylor‐Goldstein equation, which governs the vertical structure of the wave perturbation as a function of density gradient. The condition for instability is balanced with the velocity vertical gradient (both of them varying linearly in the mixing height $H$), resulting in the following necessary condition:

$Ri=\frac{gH({\rho }_2-{\rho }_1)}{{\rho }_0{(u_1-u_2)}^2}<\frac{1}{4},$E3

where $Ri$ is called Richardson number, [36]. Extension of the analysis is dedicated to the viscous KH or RT instabilities. In this case the mode analysis shows [37–39], that the neutral curve for KH instability may have shorter regions of stability, depending on the values of viscosity and wavenumber. It is interesting to note that viscosity may play a destabilizing role[38], where air‐water stratified system is considered. This analysis usually done, again, by introducing Squire’s reduction, showing that the most dangerous three‐dimensional disturbance is flat (two‐dimensional). Inclusion of the surface tension additionally stabilizes the problem [38] in the region of small wavenumbers; however, this subject is beyond the scope of the chapter. In the later stages of instability, viscosity stabilizes the flow. The RT instability remains unstable under inverse condition (Eq. (1)) despite the action of viscosity.

The RT and KH instabilities in compressible fluid are more complex to study. One of the first results of KH instability is presented in [33] for ideal gas with basic thermodynamic state characterized by the constant speed of sound $c=\gamma P/\rho$, where $P$ is pressure, $\rho$ is density, and $\gamma$ is the adiabatic index (ratio of specific heats). It is noted that there are three types of disturbances associated with the eigenvalues of the linearized operator: subsonic ($M<1$), sonic ($M=1$), and supersonic ($M>1$), where Mach number $M=u/c$, $u$ is the velocity of the basic parallel flow, i.e., $u=u_1-u_2$. Supersonic and sonic disturbances are responsible for another type of instability called Richtmyer‐Meshkov instability that is not considered here. So, considering subsonic disturbances, it was noticed that for a wavenumber $k$, the increase of $M$, complex $kM$ increases, and it leads to the decrease of such wavemodes that are linearly unstable. It is shown that the Fjortoft’s theorem [40] (let there exists a point $y_s$, that $u^{″}(y_s)=0$. If the flow is unstable, then $u^{″}(u-u(y_s))<0$ for some point $y\in [y_{min},y_{max}]$) also presents sufficient condition for stability to subsonic disturbances. Results of eigenmodes and growth rates for KH instabilities for compressible flow are presented in [33, 41, 42].

The inviscid incompressible RT instability gives exponential growth of the amplitude ($\theta (k)$) of infinitesimal perturbations [29]:

$\theta (k)\sim e^{\alpha (\rho ,k,g)t},\alpha (\rho ,k,g)=\sqrt{A(\rho )gk},$E4

with

$A=\frac{{\rho }_2-{\rho }_1}{{\rho }_2+{\rho }_1},$E5

where $\alpha$ is the temporal growth rate, $k$ is the spatial wavenumber, and A is called the Atwood number. The growth rate is an important property of such instabilities. It can be estimated as the magnitude of the least stable eigenvalue. The growth rate of a linearized system can be detected experimentally or numerically while the secondary solution remains stable. With further complication of the flow, i.e., development of solution bifurcating branches, the value of growth rate changes. In accordance with Eq. (4), the growth of modes is more intensive for high $k$. Inclusion of viscosity $\mu$ in the model results in the formation of maximum growth rate for specific $k(\mu )$ with $\theta (k)\to 0$ with $k\to \infty$. So the growth rate and the mixing layer length become functions of time that can be compared with experimental or numerical results. The mixing length is a function of initial wavenumber disturbance, Atwood number, gravity constant, and diffusion coefficient.

Compressibility effects of the RT instability include both stabilizing and destabilizing effects compared to the incompressible flow [31, 43–46]. For infinite domains, the growth rate $\theta (k)$ obtained for the compressible case is bounded by the growth rates obtained for the corresponding incompressible flows. This result is not affected by the presence of surface tension or viscosity. As the speed of sound is increased, the rate of growth increases toward the value obtained for the corresponding constant density incompressible flow, while as $\gamma$ increases, $\theta (k)$ decreases toward the value obtained for the corresponding incompressible flow with exponentially varying density. Another effect of compressibility is the decrease of the average local Atwood number for onset of instability. It is interesting to note [47] that the growth rate for compressible RT instability of two immiscible inviscid barotropic fluids was found to be $e^{t\sqrt{\left|k\right|}}$. In addition, solutions that grow arbitrarily quickly in the Sobolev space $W_2^k$, which leads to an ill‐posedness result for the linearized problem were presented in [47]. Here, it is shown that the nonlinear problem does not admit reasonable estimates of solutions for small time in terms of the initial data. The compressible inviscid RT instability must be regularized, e.g., with viscosity.

There are also some papers on the bifurcation analysis of RT and KH instabilities. Very good experimental results for RT instability are presented in [48], where stabilization is performed using zimuthally rotating magnetic field for ferrofluids. A subcritical pitchfork bifurcation was found during the stability loss, where this bifurcation depended on the magnetic field strength as a parameter. Furthermore, it is shown on phase space portrait (Figure 7 in [48]) that further increase of complexity goes through Andronov‐Hopf bifurcation. Another paper containing some information on bifurcation for KH instability is [49], where authors obtained reliable results (DNS and eigenvalue analysis) confirming the emergence of self‐oscillatory solution. Hence, Andronov‐Hopf bifurcation for KH instability is confirmed. Another paper on KH bifurcation analysis in MHD flow is presented in Ref. [50], where Andronov‐Hopf bifurcation was found.

Finally, we can conclude the following. Initial stability loss is well studied for RT and KH instabilities. One can use such data as stability criteria, growth rates, and eigenmodes to verify the linear analysis in numerical methods. There is evidence of bifurcation route to chaos in RT and KH instabilities, but this analysis is very complex and only initial stages (pitchfork and Andronov‐Hopf bifurcations) are observed. The subject of this chapter is the analysis of compressible flow from the dynamical systems point of view and the attempt to bring together data on the bifurcation analysis of laminar‐turbulent transition of KH and RT instabilities.

1.3. Definitions

Let us introduce some definitions and notations that are required for the analysis of laminar‐turbulent transition from the nonlinear dynamic system point of view. For some of these definitions, we are using [51, 52]. We are considering a system (infinite‐dimensional, in general):

$\frac{\partial u(x,t)}{\partial t}=F(u,x,R,t)\mathrm{.}$E6

Let $x\in \mathbb{M}\in X$, $t$ is time, $\mathbb{M}$ is a phase space, and ${\varphi }^t(x):\mathbb{M}\times \mathbb{R}\to \mathbb{M}$ is a phase flow on $\mathbb{M}$ and $R$ is a vector of parameters (e.g., similarity complexes in case of fluid dynamics) of size $Rp$. Let $\gamma (t)\in \mathbb{M}$ be a phase trajectory of a solution of Eq. (6).

Definition 1. Compact set $B\in \mathbb{M}$, invariant with respect to ${\varphi }^t$ is an attracting set if there exists a neighborhood $U$ of the set $B$ such that for almost all $x\in U$, ${\varphi }^t(x)\to B,t\to \infty$.

Definition 2. Attractor of the system (6) is an indecomposable attracting set, [51, 52].

Definition 3. Regular attractors of the system (6) are stable singular points, stable invariant limited cycles, and stable invariant tori of any dimension.

Definition 4. Singular attractor of Eq. (6) is an attractor that is not regular.

Definition 5. System (6) has a chaotic solution if and only if it has a singular attractor.

Definition 6. Bifurcation points of the system (6) are those and only those points in which there is no continuous dependence of the phase portrait of the system by its parameters.

Linearize the system (4) on its regular attractor $u_0(x,t)$ or the parameter value $R_0$:

$\frac{\partial y(x,t)}{\partial t}=\text{L}(u_0,x,t,R_0)y(x,t),$E34

where $\text{L}:\mathbb{M}\to \mathbb{M}$ is a linear operator and

$y(x,t)=u_0(x,t)-u(x,t)\mathrm{.}$E35

If the spectrum of operator $\text{L}$ lies in the left complex half‐plane when $\forall j\in Rp,R_j\le R_{\mathrm{0,}j},$ and if $\exists j,R_j>R_{\mathrm{0,}j}$ with one or few its eigenvalues that are moving in the right half‐plane, then the system (6) has a bifurcation at $R_0$.

Definition 7. The infinite sequence of bifurcations: $\forall j\in Rp,R_{j,n}\to R_{j,*},n\to \infty$ is known as the cascade of bifurcations.

Any singular attractor is a limit of a cascade of bifurcations of regular attractors. For example, the simplest singular Feigenbaum attractor is a nonperiodic trajectory, which is the limit of the cascade of period doubling bifurcations of stable cycles. So, any cascade of bifurcations is known as a scenario of transition to chaos.

Definition 8. Let the trajectories of the Navier‐Stokes dynamical system in the phase space be located at least on one singular attractor (there can be more attracting states since the system exhibits multistability phenomenon), then such solution of the system is called turbulent.

2. Governing equations and numerical methods

2.1. Governing equations

We are considering conservation laws for viscous gas dynamics. This system can be written in the following differential form [53]:

$\begin{array}{c}\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x_j}\left[\rho u_j\right]=0\\ \frac{\partial }{\partial t}\left(\rho u_i\right)+\frac{\partial }{\partial x_j}\left[\rho u_i{u}_j+p{\delta }_{ij}-{\tau }_{ji}\right]=\rho g_i,\{i,j,k\}=\mathrm{1,2,3};\\ \frac{\partial }{\partial t}\left(\rho E\right)+\frac{\partial }{\partial x_j}\left[\rho u_{j}E+u_{j}p-u_i{\tau }_{ij}\right]=\rho u_k{g}_k,\\ \rho E=\frac{1}{2}\rho u^2+\rho e,\\ p=(\gamma -1)(E-1/2\rho u^2)\mathrm{.}\end{array}$E7

Here, we are considering some bounded domain $\Omega \subset {\mathbb{R}}^3$ with boundary $\partial \Omega$; scalar functions $f$ are defined as $f:\Omega \times [\mathrm{0,}\mathcal{T}]\to \mathbb{R}$; vector‐functions $f$ are defined as $f:\Omega \times [\mathrm{0,}\mathcal{T}]\to {\mathbb{R}}^3$, where $\mathcal{T}$ is some defined finite time. Then, $\rho E$ is a scalar function of the full gas energy; $e$ is a scalar function of the internal gas energy; $\gamma$ is the adiabatic index (value 1.4 is used for calculations); $p$ is a scalar pressure function; $u$ is a gas velocity vector function; $\rho$ is a scalar function of gas density; $g$ is an external forcing. We assume the Einstein summation rule. We also assume that the gas is Newtonian with the following viscous tensor:

${\tau }_{ij}=2\mu S_{ij},$E8

where $\mu$ is a dynamic gas viscosity and strain rate tensor is defined as:

$S_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)-\frac{1}{3}{\delta }_{ij}\frac{\partial u_k}{\partial x_k}\mathrm{.}$E9

Please note that we ignore secondary viscosity coefficient due to complexity. We use the integral form of Eq. (7) for problems that may have discontinuous solutions. We use the main bifurcation parameter that is considered in the paper—the Reynolds number. It is found as:

$R=\frac{L_C\left|\right|u_C\left|\right|{\rho }_C}{\mu },$E10

where $L_C$ is the characteristic length scale of the problem, $\left|\right|u_C\left|\right|$ is the norm of the characteristic velocity vector, ${\rho }_C$ is a characteristic density. With this in mind, we can rewrite Eq. (8) as follows:

${\tau }_{ij}=\frac{2}{R}{S}_{ij}\mathrm{.}$E11

We will use Eq. (8) if we use the term “viscosity coefficient” and we use Eq. (11) if we use the term “Reynolds number”.

This method of analysis is based on the construction of phase space and Poincaré sections. The method is analogous to [51, 1]. In order to construct multidimensional Poincaré sections, we take data from different points of the physical space and use slices with the given gap thickness $\delta$. For all results presented here, we use $\delta =5\cdot {10}^{-7}$.

2.2. Numerical methods

In this section, we give a brief description of numerical methods that are used to perform simulation of problems and bifurcation analysis. We are forced to relay to the numerical methods since one can basically derive only primary or secondary solutions analytically. All other bifurcations of the main solution cannot be derived.

2.2.1. Main flow solver

The flow solver is constructed using finite volume‐finite element method. The inviscid part is solved using modified WENO schemes of order 7 or 9; for more information, see [54]. The Riemann solver is a Godunov exact solver, see [55]. For viscous part of equations, we use finite element nodal method that is described for Poisson equation in Ref. [12]. Time integration is explicit (WENO)‐implicit (FEM) method of the third order. The method is implemented on multiGPU architecture. The efficiency of the method and specific features of implementation are outlined in Ref. [56]. The inflow and outflow boundary conditions are set using characteristic reconstruction with damping zones; for more information see Ref. [57]. All discretization of the problems is done using results of modified linear analysis for WENO schemes and “burglization” numerical experiments, provided in Ref. [58].

2.2.2. Eigenvalue solver

In order to study the linear stability of the problem, we are using implicitly restarted Arnoldi (IRA) method designed to be used on multiGPU architecture. We use matrix‐free method to call the main solver and analytical linearization of governing equations in order to construct Jacobi matrix‐vector operator. For troublesome problems, we apply either exponential or Cayley transformations. See [12, 59] for more details on implementation features and application examples. In all calculations, unless specified otherwise, we use $k=10$ target eigenvalues with additional Krylov vectors totaling $m=8$, i.e., Krylov subspace dimension is $\text{dim}(\mathcal{K})=k\cdot m=80$.

3. Verification

Verification of the main numerical scheme is conducted in [55, 58], and IRA eigensolver is verified in Ref. [59]. In this section, we present verification results for KH and RT instabilities.

3.1. KH instability verification

The primary problem that can be used as a benchmark is the numerical calculation of the growth rate. It can be done using IRA eigenvalue solver for the periodic KH problem in X‐direction and nonreflecting boundary condition in Y‐direction with $\Omega =[L\times H]$ and initial mixing layer thickness $\delta$. This problem originated from [33] and was numerically calculated in Refs. [41, 42]. The domain length in X‐direction should fit at least one entire wavelength of the desired mode and was chosen, analogously to [42], to be exactly one wavelength of the least stable mode scaled by mixing layer thickness $\delta$. Assuming a wavenumber $k$ of the least stable mode, the length of the domain in X‐direction would result as $L=2\pi \delta /k$. The domain length in Y‐direction is chosen as $H=8L$.

$\begin{array}{c}{\rho }_0=\mathrm{1.24,}\\ u{(x,y)}_0=u{(y)}_0=\text{tanh}\left(\frac{y-L/2}{\delta }\right)\mathrm{.}\end{array}$E12

Pressure is chosen as to satisfy Mach number that will vary depending on the particular problem.

The solution is performed on $32\times 256$ grid with $k=8$ and $m=5$ IRA parameters, that results in the $dim(\mathcal{K})=40$. The results are presented in Table 1 for different Mach numbers and modulus of the least stable eigenvector is presented in Figure 1. We can see good convergence properties of the growth rate for all least stable modes. The eigenvector clearly indicates the formation of two symmetric billows.

Figure 1.

Modulus of the least stable normalized eigenvector of total energy for KH 2D test problem. The problem is rotated 90° CW.

M	Computed α	Analytical α [33]	α [41]	α [42]
0.1	0.187312	0.187	0.188	0.189
0.5	0.1411	0.141	0.1411	0.1411
0.9	0.054743	0.055	0.054723	0.0547

Table 1.

Results of the validation for growth rates α as functions of Mach number.

Another test case considered is with regard to [60], where a 2D periodic KH problem is simulated with very high accuracy. For this purpose, benchmark code in [60] was used with $4096\times 4096$ grid points. The problem is set as KH instability of shear flow in rectangular unit domain $\Omega$ with periodic boundary conditions in all directions. Initial conditions are set as:

$\begin{array}{l}\rho (x\mathrm{,0})=\{\begin{array}{c}{\rho }_1-{\rho }_m{e}^{\frac{y-1/4}{L}};1/4>y\ge \mathrm{0,}\\ {\rho }_2+{\rho }_m{e}^{\frac{-y+1/4}{L}};1/2>y\ge 1/\mathrm{4,}\\ {\rho }_2+{\rho }_m{e}^{-\frac{3/4-y}{L}};3/4>y\ge 1/\mathrm{2,}\\ {\rho }_1-{\rho }_m{e}^{-\frac{y-3/4}{L}};1>y\ge 3/\mathrm{4,}\end{array}\\ u_x(x\mathrm{,0})=\{\begin{array}{c}{U}_1-U_m{e}^{\frac{y-1/4}{L}};1/4>y\ge \mathrm{0,}\\ U_2+U_m{e}^{\frac{-y+1/4}{L}};1/2>y\ge 1/\mathrm{4,}\\ U_2+U_m{e}^{-\frac{3/4-y}{L}};3/4>y\ge 1/\mathrm{2,}\\ U_1-U_m{e}^{-\frac{y-3/4}{L}};1>y\ge 3/4\end{array}\end{array}$E13

and

$\begin{array}{c}{u}_y(x\mathrm{,0})=0.01\sin (4\pi x),p(x,y\mathrm{,0})=\mathrm{2.5.}\end{array}$E14

Here $U_m=(U_1-U_2)/2$, ${\rho }_m=({\rho }_1-{\rho }_2)/2$, ${\rho }_1=1$, ${\rho }_2=2$, $U_1=1/2$, $U_2=-1/2$, $L=0.025$.

For our simulation, we use a 3D version of the problem and set up grid $256\times 256\times 256$ such that $u_z(x\mathrm{,0})=0$ with periodic boundary conditions in $z$‐direction. The problem is executed until $\mathcal{T}=1.5$ is reached. We also compare these results with spectral code that is used for incompressible Navier‐Stokes simulation. We monitor the maximum y‐direction kinetic energy of perturbations, as $k_y(t)={{\displaystyle \text{max}}}_{x\in \Omega }\rho u_y^2/2$, and we monitor the growth rate of the amplitude for the unstable $u_y$ mode using an FFT. The results are presented in Figure 2. One can see a very reasonable agreement with reference data. Small difference in the evolution of growth rate ${\alpha }_y$ in the initial stage can be explained by the usage of a 3D code instead of 2D as in reference data.

Figure 2.

Growth rate ${\alpha }_y$ and maximum y‐direction kinetic energy $k_y$. Ref1 is [60] and Ref2 is data from spectral code.

3.2. RT instability verification

Verification data on RT instability are freely available in the literature, e.g., [61–67] and many others. One of the benchmarks is the turbulent mixing layer length. It is known that the RT process is basically divided into three stages: (1) linear stage related to the linear stability analysis where initial perturbations $\theta (k)$ are growing exponentially as in Eq. (4); (2) nonlinear stage where the formation of spikes and bubbles occurs. In this stage bubbles are traveling toward heaver fluid and spikes, toward the lighter one. Growth rate of bubbles is proportional to ${\theta }_B$~$\sqrt{2\pi g/k}t$, and growth rate of spikes is proportional to ${\theta }_J$~$g{t}^2$; (3) stage of dissipation of regular structures. In this stage, KH instability is developing on the contacts moving with shear velocity and destabilization of lighter bubbles in heavier fluid due to the secondary RT instability. The onset of turbulence depends on the initial perturbation wavenumber and viscosity coefficient.

We follow the problem setup in [67] and investigate the mixing length $L_m$ as a function of time with $L_m=\alpha gAg{t}^2$, where $\alpha$ is a constant. We consider the domain $\Omega =[3\times 1\times 1]$ with gravity direction in $-x$, and we use $900\times 300\times 300$ number of elements. We test two setups with random and sinusoidal initial perturbations. The latter results are presented in Figure 3.

Figure 3.

Density for RT instability test case with $A=1/3$ and harmonic initial disturbances: ${\tau }_{Left}=0.3\text{s}$, ${\tau }_{Right}=1.9\text{s}$.

In Figure 4, we can see the mixing length for random perturbations and for sinusoidal perturbations with comparison to the literature on different length of waves. One can see that the results are close to the ones in the literature. We identified the constants $\alpha$ for sinusoidal perturbations: ${\alpha }_{short}=0.0271$ for short wavelength and ${\alpha }_{short}=0.057$ for long wavelength. These results are very close to those provided in the literature.

Figure 4.

Mixing length of RT instability for random (left) and sinusoidal (right) initial perturbations. Ref 1,2—[67], Ref 3—[61, 63], Ref4—[62].

4. Bifurcation analysis

4.1. Phase space for supersonic flow

4.1.1. Inviscid flow phase space properties

Let us consider the following vector form of inviscid system (7) for 2D case:

$\begin{array}{c}\frac{\partial U}{\partial t}+\frac{\partial F}{\partial x}+\frac{\partial G}{\partial y}=0,\\ U=[\rho ;\rho u_x;\rho u_y;\rho E{]}^T,\\ F=[\rho u_x;\rho u_x^2+p;\rho u_x{u}_y;(\rho E+p)u_x{]}^T,\\ G=[\rho u_y;\rho u_y{u}_x;\rho u_y^2+p;(\rho E+p)u_y{]}^T\mathrm{.}\end{array}$E15

The system can be linearized:

$\begin{array}{c}\frac{\partial U}{\partial t}+\text{A}\frac{\partial U}{\partial x}+\text{B}\frac{\partial U}{\partial y}=0,\\ \text{A}=\partial F/\partial U;\text{B}=\partial G/\partial U,\end{array}$E16

with eigenvalues being:

${\lambda }_1=u\cdot n-c,{\lambda }_2=u\cdot n,{\lambda }_3=u\cdot n,{\lambda }_4=u\cdot n+c,$E36

where $c=\sqrt{\frac{\gamma p}{\rho }}$ is the speed of sound, $n$ is a unit directional vector. All eigenvalues for all $x,y$ are real since system (16) is hyperbolic. We are considering a special case with the base flow being parallel to $x$ axis, i.e., $u=(u_0{\mathrm{,0})}^T$ and $n{=(\mathrm{1,0})}^T$. In this case, the system of eigenvectors becomes one dimensional, and we can reduce Eq. (16) to one‐dimensional case. Let us first consider a simple initial value problem for system of advection equations:

$\begin{array}{c}{U}_t+\text{A}{U}_x=\mathrm{0,}\\ x\in \mathbb{R};t\in {\mathbb{R}}^{+},\\ U=(u_1,u_2{)}^T,\\ u_1(x\mathrm{,0})=\mathrm{1,}{u}_2(x\mathrm{,0})=\{\begin{array}{c}{e}^{ikx};x\ge 0\\ 1;x<0.\end{array},\end{array}$E17

where $k$ is a wavenumber.

Lemma 1: Let all eigenvalues ${\lambda }_j\in \mathbb{R}$ of $\mathtt{A}$ in Eq. (17) satisfy ${\lambda }_j\ge 0$. Then, $u_1$ and $u_2$ are constant for $x<\mathrm{0,}\forall t\in {\mathbb{R}}^{+}$.

Proof: Let linear operator $\mathtt{A}=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)$; $a,b,c,d\in \mathbb{R}$. Eigenvalues and eigenvectors are found as:

$\begin{array}{c}{\lambda }_{\mathrm{1,2}}=\left(\begin{array}{c}1/2d+1/2a+1/2\sqrt{{(d-a)}^2+4cb}\\ 1/2d+1/2a-1/2\sqrt{{(d-a)}^2+4cb}\end{array}\right),\\ H_{\mathrm{1,2}}=\left(\begin{array}{cc}\frac{2b}{d-a+\sqrt{{(d-a)}^2+4cb}};& \frac{2b}{d-a-\sqrt{{(d-a)}^2+4cb}}\\ 1& 1\end{array}\right)\end{array}$E18

Let us assume that ${\lambda }_1>\mathrm{0,}{\lambda }_2\ge 0$. We apply characteristic variables $W_0=H^{-1}{U}_0$:

$\left(\begin{array}{c}{W}_{\mathrm{0,1}}\\ W_{\mathrm{0,2}}\end{array}\right)=\left(\begin{array}{c}\{\begin{array}{c}\frac{1}{4}\frac{(d-a+\sqrt{{(d-a)}^2+4cb})(-d+a+\sqrt{{(d-a)}^2+4cb}+2b)}{b\sqrt{{(d-a)}^2+4cb}};x<0.\\ \frac{1}{4}\frac{(d-a+\sqrt{{(d-a)}^2+4cb})(-d+a+\sqrt{{(d-a)}^2+4cb}+2b\cdot e^{ikx})}{b\sqrt{{(d-a)}^2+4cb}};x\ge 0.\end{array}\\ \{\begin{array}{c}\frac{1}{4}\frac{(d-a-\sqrt{{(d-a)}^2+4cb})(-d+a-\sqrt{{(d-a)}^2+4cb}+2b)}{b\sqrt{{(d-a)}^2+4cb}};x<0.\\ \frac{1}{4}\frac{(d-a-\sqrt{{(d-a)}^2+4cb})(-d+a-\sqrt{{(d-a)}^2+4cb}+2b\cdot e^{ikx})}{b\sqrt{{(d-a)}^2+4cb}};x\ge 0.\end{array}\end{array}\right)$E19

Matrix $\mathtt{A}$ is diagonalized as $\Lambda =H^{-1}AH$. Then, we can rewrite (17) as:

$\begin{array}{c}{H}^{-1}{U}_t+H^{-1}A{U}_x=0;\\ \mathbf{\Lambda }=diag({\lambda }_i);\\ {\mathbf{W}}_t+\mathbf{\Lambda }{\mathbf{W}}_x=0,\end{array}$

and the last vector equation can be separated:

$${(W_i)}_t+{\lambda }_i{(W_i)}_x=0;i=\mathrm{1,2.}$$

Then characteristic solution is:

$$W_i=W_{\mathrm{0,}i}(x-{\lambda }_{i}t),i=\mathrm{1,2.}$$

with the application of initial data (Eq. (19)):

$\begin{array}{c}\left(\begin{array}{c}{W}_1\\ W_2\end{array}\right)=\\ \left(\begin{array}{c}\{\begin{array}{c}A+;(x-\frac{d+a+\sqrt{{(d-a)}^2+4cb}}{2}t)<0.\\ B+;(x-\frac{d+a+\sqrt{{(d-a)}^2+4cb}}{2}t)\ge 0.\end{array}\\ \{\begin{array}{c}A-;(x-\frac{d+a-\sqrt{{(d-a)}^2+4cb}}{2}t)<0.\\ B-;(x-\frac{d+a-\sqrt{{(d-a)}^2+4cb}}{2}t)\ge 0.\end{array}\end{array}\right)\end{array}$

$$A\pm =\frac{1}{4}\frac{(d-a\pm +\sqrt{{(d-a)}^2+4cb})(-d+a\pm \sqrt{{(d-a)}^2+4cb}+2b)}{b\sqrt{{(d-a)}^2+4cb}},$$

Here:

$B\pm =\frac{1}{4}\frac{(d-a\pm \sqrt{{(d-a)}^2+4cb})(-d+a\pm \sqrt{{(d-a)}^2+4cb}+2b\cdot e^{ik(x-\frac{d+a\pm \sqrt{{(d-a)}^2+4cb}}{2}t)})}{b\sqrt{{(d-a)}^2+4cb}},$E20

Now we turn to conservative variables $U=HW$:

$\left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{c}\frac{2b}{d-a+\sqrt{{(d-a)}^2+4cb}}{W}_1+\frac{2b}{d-a-\sqrt{{(d-a)}^2+4cb}}{W}_2\\ W_1+W_2\end{array}\right)\mathrm{.}$E21

It is clear from the obtained solution that eigenvectors in Eq. (21) are independent of $x-{\lambda }_{i}t$. So, the solution can be reformed as:

$\left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{c}{\eta }_1\{\begin{array}{c}{W}_1\\ W_1(e^{ik(x-{\lambda }_{1}t)})\end{array}+{\eta }_2\{\begin{array}{c}{W}_2\\ W_2(e^{ik(x-{\lambda }_{2}t)})\end{array}\\ \{\begin{array}{c}{W}_1\\ W_1(e^{ik(x-{\lambda }_{1}t)})\end{array}+\{\begin{array}{c}{W}_2\\ W_2(e^{ik(x-{\lambda }_{2}t)})\end{array}\end{array}\right),$E22

where ${\eta }_1,{\eta }_2$ are values in Eq. (21) that are independent of $x$ and $t$: ${\eta }_1=\frac{2b}{d-a+\sqrt{{(d-a)}^2+4cb’}}$ ${\eta }_2=\frac{2b}{d-a-\sqrt{{(d-a)}^2+4cb}}$. Let us consider two possible cases:

1. ${\lambda }_1>\mathrm{0,}{\lambda }_2=0$

   $\begin{array}{c}1.x\ge {\lambda }_1\mathrm{.}\Rightarrow \left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{c}{\eta }_1{W}_1(e^{ik(x-{\lambda }_{1}t)})+{\eta }_2{W}_2(e^{ikx})\\ W_1(e^{ik(x-{\lambda }_{1}t)})+W_2(e^{ikx})\end{array}\right)\mathrm{.}\\ 2.x>\mathrm{0,}x<{\lambda }_1\mathrm{.}\Rightarrow \left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{c}{\eta }_1{W}_1+{\eta }_2{W}_2(e^{ikx})\\ W_1+W_2(e^{ikx})\end{array}\right)\mathrm{.}\\ 3.x\le 0.\Rightarrow \left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{c}{\eta }_1{W}_1+{\eta }_2{W}_2\\ W_1+W_2\end{array}\right)\mathrm{.}\end{array}$E23

2. ${\lambda }_1>\mathrm{0,}{\lambda }_2>\mathrm{0,}{\lambda }_1\ge {\lambda }_2$

$\begin{array}{c}1.x\ge {\lambda }_1\mathrm{.}\Rightarrow \left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{c}{\eta }_1{W}_1(e^{ik(x-{\lambda }_{1}t)})+{\eta }_2{W}_2(e^{ik(x-{\lambda }_{2}t)})\\ W_1(e^{ik(x-{\lambda }_{1}t)})+W_2(e^{ik(x-{\lambda }_{2}t)})\end{array}\right)\mathrm{.}\\ 2.{\lambda }_1\ge x\ge {\lambda }_2\mathrm{.}\Rightarrow \left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{c}{\eta }_1{W}_1+{\eta }_2{W}_2(e^{ik(x-{\lambda }_{2}t)})\\ W_1+W_2(e^{ik(x-{\lambda }_{2}t)})\end{array}\right)\mathrm{.}\\ 3.x\ge \mathrm{0,}x<{\lambda }_2\mathrm{.}\Rightarrow \left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{c}{\eta }_1{W}_1+{\eta }_2{W}_2\\ W_1+W_2\end{array}\right)\mathrm{.}\\ 4.x<0.\Rightarrow \left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{c}{\eta }_1{W}_1+{\eta }_2{W}_2\\ W_1+W_2\end{array}\right)\mathrm{.}\end{array}$E24

After some cumbersome but simple manipulations one can show that ${\eta }_1{W}_1+{\eta }_2{W}_2=1$ and $W_1+W_2=1$. $\forall \{a,b,c,d\}\in \mathbb{R}$. □

Lemma 2: Let there exist at least one eigenvalues ${\lambda }_0$ of $\mathtt{A}$ from Eq. (17), that ${\lambda }_0<0$. Then values of $u_1$ or $u_2$ for $x<0$ are not constant for $t\in {\mathbb{R}}^{+}$.

Proof: For the proof, we are considering one example and then use the proof of Lemma 1. Let ${\lambda }_1>\mathrm{0,}{\lambda }_2<0$ by the Lemma formulation. Let $\mathtt{A}=\left(\begin{array}{cc}-1& 1\\ 0& 1\end{array}\right)$. Then, $\lambda =\left(\begin{array}{c}1\\ -1\end{array}\right)$ and $H=\left(\begin{array}{cc}\frac{1}{2}& 1\\ 1& 0\end{array}\right)$. Using characteristic variables:

$\left(\begin{array}{c}{W}_{\mathrm{0,1}}\\ W_{\mathrm{0,2}}\end{array}\right)=(\begin{array}{c}\{\begin{array}{c}1;x<0\\ e^{ikx};x\ge 0\end{array}\\ \{\begin{array}{c}\frac{1}{2};x<0\\ 1-\frac{1}{2}{e}^{ikx};x\ge 0\end{array}\end{array})$E25

Diagonalizing operator $\mathtt{A}$: $\mathrm{\Lambda }=H^{-1}\text{A}H=\left(\left(\begin{array}{cc}0& 1\\ 1& -\frac{1}{2}\end{array}\right)\cdot \left(\begin{array}{cc}-1& 1\\ 0& 1\end{array}\right)\right)\cdot \left(\begin{array}{cc}\frac{1}{2}& 1\\ 1& 0\end{array}\right)=diag({\lambda }_i)$ and arriving at two independent equations. The solution in characteristic variables is:

$\left(\begin{array}{c}{W}_1\\ W_2\end{array}\right)=(\begin{array}{c}\{\begin{array}{c}1;(x-t)<0\\ e^{ik(x-t)};(x-t)\ge 0\end{array}\\ \{\begin{array}{c}\frac{1}{2};(x+t)<0\\ 1-\frac{1}{2}{e}^{ik(x+t)};(x+t)\ge 0\end{array}\end{array}).$E26

Returning to the original variables:

$\left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{c}\{\begin{array}{c}\frac{1}{2};(x-t)<0\\ \frac{1}{2}{e}^{ik(x-t)};(x-t)\ge 0\end{array}+\{\begin{array}{c}\frac{1}{2};(x+t)<0\\ 1-\frac{1}{2}{e}^{ik(x+t)};(x+t)\ge 0\end{array}\\ \{\begin{array}{c}1;(x-t)<0\\ e^{ik(x+t)};(x-t)\ge 0\end{array}\end{array}\right).$E27

Let us consider the case only for $x<0$.

$\begin{array}{c}1.x\ge (x+t);x<0.\Rightarrow \left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{c}\frac{3}{2}-\frac{1}{2}{e}^{ik(x+t)}\\ 1\end{array}\right)\mathrm{.}\\ 2.x<(x+t);x<0.\Rightarrow \left(\begin{array}{c}{u}_1\\ u_2\end{array}\right)=\left(\begin{array}{c}1\\ 1\end{array}\right)\mathrm{.}\end{array}$E28

We now turn to random real variables in matrix $\mathtt{A}$ that satisfy the condition $\exists {\lambda }_i<0$. Solutions (27) and (22) only depend on matrix variables and are identical in the other. So, this particular example is valid for any values in matrix $\mathtt{A}$ and only depends on signs of eigenvalues of $\mathtt{A}$. □

Lemma 3: Assume that matrix $\mathtt{A}$ is diagonalizable. Then Lemmas 1 and 2 are valid for $\mathtt{A}\in {\mathbb{R}}^{3\times 3}$.

Proof: The proof is straightforward, since the operator only depends on eigenvalues. Then follow the proofs for Lemmas 1 and 2. □

Now we turn our attention to Eq. (15). We assume one space dimension of the system and consider a Riemann problem for some line $x=x_0$ and consider variables on the left ($L$) and on the right ($R$) from this line. We assume that $U_L$ is an unperturbed constant vector and that $U_R$ contains perturbations. We also assume that $p_L<p_R$ and that there is a stationary shock wave with shock speed $S=0$ at line $x=x_0$. In this case, we can conclude the following:

Proposition 1: Let there be such solution to the system (15), that at some line $x=x_0$, there is a stationary shock wave with $p_L<p_R$ and $U(x\mathrm{,0}),\forall x\le x_0$ is unperturbed. Then, $U(x,t),\forall x\le x_0$ remains unperturbed at $\forall t\in {\mathbb{R}}^{+}$.

Proof: In accordance with the condition of the proposition, the local solution of Eq. (15) at $x_0$ is a shock wave solution that is located in the first characteristic field, i.e., ${\lambda }_1=u-ct$, since $p_L<p_R$. The second characteristic field is a contact discontinuity with ${\lambda }_1=u$ with the second contact discontinuity being parallel to the flow, since $n{=(\mathrm{1,0})}^T$. The third characteristic field contains an expansion wave, ${\lambda }_3=u+ct$. Then, ${\lambda }_1=0$, since $S=0$, and so the local Jacobi matrix eigenvalues are ordered as:${\lambda }_1=0\le {\lambda }_2\le {\lambda }_3$.

Since the system is hyperbolic and the local Jacobi matrix has linear independent eigenvalues, we can transform, locally, system (15) to the characteristic variables with Jacobi matrix diagonalization. We can then apply Lemmas 3 and 1 to complete the proof.

It is more difficult to prove the same proposition for a general case; however, the same conclusion is valid for 3D problem in case of a flat plane shock wave solution. It concludes that the phase space of the system for hyperbolic equations is separated into independent or slightly dependent regions by shock waves. If the shock wave is stationary, then no perturbations are transformed through it. We will demonstrate this with the numerical solution of the problem with the formation of a stationary shock wave.

4.1.2. Numerical simulation

We wish to analyze properties of the phase space of Eq. (7) in case of appearance of strong shocks and confirm results of the previous subsection. More results can be found in Ref. [55]. We are considering a rectangular domain $\Omega$, sized $8\times 3\times 3$. A small rectangular body $B$ is located in the center of the domain. This results in the formation of the front stationary detached shock wave in front of the body in case of the supersonic flow. Gas is homogeneous in $\Omega \backslash B$ at initial conditions with $\rho =\mathrm{1.24,}u{=(\mathrm{1,0,0})}^T,M_0=2$. Wall boundary conditions are set for all $\partial {\Omega }_W$ in accordance with considering equations, except for inflow and outflow boundaries. No‐slip conditions are set for viscous gas flow, and slip conditions are set for Eq. (7) with $\mu =0$. We impose two unphysical boundary conditions: inflow boundary $\partial {\Omega }_I$ on the left plane $yz$ for $x=0$ (${\rho }_0=\mathrm{1.24,}u{=(\mathrm{1,0,0})}^T,M_{\infty }=2$) and outflow boundary $\partial {\Omega }_O$ on the right plane $yz$ for $x=8$. These boundary conditions are set with the procedure described in Section 2.2.

We are performing the simulation for quasi‐stationary solution on the grid of $500\times 200\times 200$ elements. Phase space portraits are constructed using $10\cdot {10}^6$ time points after the flow is quasi‐stationary. The viscosity coefficient is set to $1\cdot {10}^{-5}$. We do not perform convergence analysis for viscous flow since this solution is used to show the validity of the phase space separation. Result of the 2D middle section of the numerical Schlieren is presented in Figure 5.

Figure 5.

Schlieren picture of the supersonic Mach 2 flow over a body. Points indicate phase space probe points.

One can see point numbers in the section that demonstrate different phase space portraits from these points. (Figure 6).

Figure 6.

Projections of phase space portraits in different points of the physical domain. Top left to right: point 1, point 2, bottom left to right: point 6, point 4.

We can monitor no perturbations in the inflow part with the formation of the fixed point in the phase space. This confirms the results of the previous subsection. At point 2, near the front shock wave, one can see the formation of the limited cycle of period 11. Perturbations from the body are substantially weakened resulting in a low dimensional process in this point. Please note that the area in which point 2 is located is isolated by shock wave configuration. The low dimensional attractor in this area is formed by the oscillations of the attached shock wave and contact discontinuity. At point 6, we can see an almost symmetrical two‐dimensional invariant torus that corresponds to the quasi‐periodic solution. This is the type of von‐Karman sheet solution generated by the flow over the body in subsonic region of the flow. In this area, we expected location of Andronov‐Hopf and secondary Hopf bifurcations. Another separated area with point 4 was tested, where we can monitor chaotic solution. The dimension of the chaos is greater than two, since we are able to monitor chaotic solution in the second Poincaré section. Other points gave qualitatively identical results.

4.2. Kelvin‐Helmholtz problem bifurcation scenario

In this subsection, we are considering the bifurcation scenario of the KH problem. The setup of the problem is the following. The domain is $\Omega =[8\times 1\times 1]$, the flow direction is from left to right. The following boundary conditions are imposed: $\partial \Omega (0,y,z)=\partial {\Omega }_0$ is the inflow boundary condition, $\partial \Omega (\mathrm{8,}y,z)=\partial {\Omega }_1$ is the outflow boundary condition, $\partial \Omega (x,y,0)=\partial \Omega (x,y,1)=\partial {\Omega }_2$ is the periodic boundary condition in $z$‐direction, and $\partial \Omega (x,0,z)=\partial \Omega (x,1,z)=\partial {\Omega }_3$ is the non‐reflecting (symmetry) boundary condition. We use simple ghost cell boundary conditions for $\partial {\Omega }_2$ and $\partial {\Omega }_3$. For $\partial {\Omega }_1$, we use appropriately chosen absorbing boundary conditions with characteristic far‐field treatment and we use characteristic boundary conditions for $\partial {\Omega }_0$. The initial conditions are identical to the boundary conditions on $\partial {\Omega }_0$ and are given below:

$\begin{array}{c}\rho {|}_{\partial {\Omega }_0}=\mathrm{1.0,}\\ u_x{|}_{\partial {\Omega }_0}=3/2+1/2\text{tanh}\left(\frac{y-1/2}{\delta }+C_z\cos (6\pi z)\right),\\ u_y{|}_{\partial {\Omega }_0}=1\cdot {10}^{-5}\sin (6\pi y),\\ u_z{|}_{\partial {\Omega }_0}=\mathrm{0,}\\ p{|}_{\partial {\Omega }_0}=\mathrm{22.14.}\end{array}$E29

The perturbations are set depending on the problem setup, either 2D (in this case $C_z=0$) or 3D (in this case $C_z=1$). The gravitational force is absent. The perturbations are not set in the classical sense, since those must be set along the $x$‐direction. However, we noticed that this kind of initial perturbations is sufficient to trigger instability. For computation, we are using a grid of $1000\times 250\times 250$ elements totaling $62.5\cdot {10}^6$ grid points. The sponge zone for $\partial {\Omega }_1$ is chosen to be 20 elements in $x$‐direction. We use the Reynold number ($R$) as bifurcation parameter with further details given below.

4.2.1. Eigenvalue analysis

First, we analyze linear stability of the main flow for the described problem setup. For $R\le 157$ the flow remained stable. The diffusion is sufficient enough to suppress instability. This solution corresponds to the stable fixed point in the phase space. With the increase of $R>158$, the flow loses stability with the formation of the initial billows. We are able to show the leading unstable eigenvector for this case; see Figure 7. As one can see, the most dangerous eigenvector is a 2D symmetrical one. This corresponds with the conclusions from Orr‐Sommerfeld equations for flat‐plane viscous flow stability. However, one can see, that the third most dangerous eigenvector has already a 3D structure. The formation of the initial billows corresponds to the Andronov‐Hopf bifurcation (two complex conjugate eigenvalues are crossing the imaginary axis). The leading eigenvalues for this problem are presented in Figure 8. One can see the monodromy matrix eigenvalues as well. One eigenvalue lies on a unit circle at the point $(+\mathrm{1,0}i)$ on the complex plane. This eigenvalue corresponds to the zero Lyapunov exponent of the dynamical system that represents the cycle direction. This limited cycle is shown in Figure 9. It was found that the length of the domain in the $x$‐direction ($L_x$) is another parameter. With the increase of the domain up to $16$, the first critical Reynolds number gradually changes to $R_c=25$. It can be explained by the wavelength in the $x$‐direction. It is interesting to note that using $L_x=4$, the inviscid flow is numerically stable (assuming constant perturbation amplitude in Eq. (29)). This leads to the conclusion that the diffusion may pay a positive role in initial instability.

Figure 7.

Modulus of the first and the third most dangerous eigenvectors of $u_y$ for the KH coaxial shear flow problem.

Figure 8.

Largest real eigenvalues for the KH problem (left) and monodromy matrix largest magnitude eigenvalues for limited cycle solution, $R=400$.

Figure 9.

Projection of the limited cycle for $R=1000$ and cycle period 11 for $R=1140$.

Let us consider time scales in the problem. The maximum physical process time $\tau$ can be roughly estimated as advection time scale ${\tau }_a~L_x/\left|u_x\right|$ for advection‐dominated flow. Providing that, this time is smaller for $L_x=4$, initial perturbations for inviscid flow are dominated by advection timescale rather than diffusion timescale ${\tau }_d~L_x^2/\mu =\mathbf{R}{\mathbf{L}}_x^2$. The timescales are equal if $R~\frac{1}{L_x\left|u_x\right|}$, so with the increase of $L_x$, one gets a smaller value of $R$ to have diffusion play a substantial role in the process. For a fixed value of $R$, the diffusion is decreasing with the decrease of $L_x$. The mechanism works as follows. Initially, the diffusion is dominated by setting a small value of $R$ with the obvious stabilizing effect to all wavenumbers. One can see that the first unstable eigenvector (Figure 7) demonstrates about ${\lambda }_{critical}$~$12\pi$. This wavelength is dominated by diffusion for small $R$. The destabilizing role of diffusion starts playing a role with the increase of $R$, resulting in sequential bifurcations. The further route to chaos depends on the initial perturbations and bifurcation parameters.

4.2.2. 2D perturbations

With the further increase of the Reynolds number, we can monitor the formation of singular limited cycles. These cycles are low dimensional attractors that are formed by a continuous set of period doubling bifurcations. All these cycles lie in the Feigenbaum cascade. We are able to monitor the formation of cycles until cycle period 11, in Figure 9 at around $R=1140$. This cycle originates from the Sharkovskii order and is one of the key cycles in Feigenbaum‐Sharkovskiy‐Magnitskii (FShM) cascade [51]. Unfortunately, we are unable to trace the cycle period three which is the last cycle in Sharkovskii direct order. It is known that chaos emerges after this cycle. With further increase of $R$ value, we witnessed the reversed cascade up to the solution of a limited cycle. All these cycles are corresponding to the 2D instabilities of the primary flow resulting in low dimensional chaos and correspond to the oscillation of billows in physical space in $xy$ plane. It is noted that the increase of domain length in $x$‐direction, or decrease of base velocity $u_x$ results in the change of the bifurcation scenario. This again can be explained by the domination of ${\tau }_a$ with the increase of $R$. We confirmed this proposition by increasing the $L_x$ up to $16$ for the fixed value of $R=2500$ thus giving the diffusion scale more time to develop. While cycle was formed for $L_x=\mathrm{8,}R=2500$, the chaotic solution emerges for $L_x=\mathrm{16,}R=2500$, so the value of $L_x$ becomes another bifurcation parameter. Moreover, in this case a two dimensionality of billows is destroyed with the formation of chaotic flow. We are unable to trace any regular attractors after that. This can be the result of disturbing diffusion action on the flow through screw symmetric viscous tensor components.

4.2.3. 3D perturbations

We found that the initial 2D instability is reformed for the 3D perturbations. This can be explained since the flow is no longer flat parallel. The change in the inflow conditions resulted in the formation of chaotic solution right after the flow is developing. It was difficult to trace any bifurcations further with only one regular structure found being a cycle, two‐dimensional and three‐dimensional invariant tori for $100<R<950$; see Figure 10.

Figure 10.

2D invariant torus ($R=\mathrm{760,}{L}_x=8$) and 3D invariant torus ($R=\mathrm{950,}{L}_x=8$) first Poincaré sections.

4.3. Coupled (Kelvin‐Helmholtz and Rayleigh‐Taylor) problem bifurcation scenario

We are considering the problem of stratified gas flow with velocity shear layer. This is a coupled KH and RT problem that can be characterized by the density difference and application of the gravitational force. We set up the problem analogously to the previous one but we have one more parameter that is the Richardson number. The boundary on $\partial {\Omega }_0$ and initial values are set as follows:

$\begin{array}{c}\rho {|}_{\partial {\Omega }_0}=1.1+\frac{1}{10}\text{erf}\left(\frac{y-1/2}{\delta }+C_z\cos (6\pi z)\right),\\ u_x{|}_{\partial {\Omega }_0}=2.0+\text{erf}\left(\frac{y-1/2}{\delta }+C_z\cos (6\pi z)\right),\\ p{|}_{\partial {\Omega }_0}=22.14+\rho gy+({\rho }_{max}-\rho )g\mathrm{.}\end{array}$E30

The Richardson number (Eq. (3)) is set depending on the problem in range of $1/5-1/8$, where $H=4\delta$ is the length of the mixing layer and ${\rho }_{max}=1.2$. We also consider two options—$C_z=0$ or $C_z=1$. The results of the 2D perturbations are analogous to the previous subsection with the formation of limited cycles period five and three presented in Figure 11. Interesting to note that this cycle was found in the reverse cascade due to the limited length, $L_x=6$. The increase of $L_x$ up to $10$ leads to the formation of chaotic solution. It was noted that the process was more related to the shear flow of gases with nonzero density ratio. The development of RT instability was not detected; thus, these series can be qualified as pure KH instability problem. The initial stage of instability is formed by low dimensional attractors with confirmed inverse Feigenbaum‐Sharkovskiy cascade.

Figure 11.

Projection of limited cycles of period three into ${\upsilon }_z\mathrm{,}{\upsilon }_y$ phase subspace for 2D perturbation problem for $R=\mathrm{2417,}Ri=1/\mathrm{6,}{L}_x=6$.

4.3.1. 3D perturbations

Results for 3D perturbations are presented in Figures 12 and 13. One can clearly see the formation of initial billows of KH instability, followed by RT instability. The latter can be clearly noticed in later stages of the flow; see Figure 13. With $L_x$ acting as bifurcation parameter in mind, we fix the length and increase value of $R$ until the solution is a regular attractor. Please mind that the initial flow is linearly unstable, and small values of $R$ may give a destabilizing effect.

Figure 12.

Coupled problem quasi‐stationary solution. Density isosurfaces.

Figure 13.

Coupled problem quasi‐stationary solution. Density sections in $yz$ planes with $x=3$ and $x=5$.

The starting point is $R=1$ for which the solution is stationary and totally diffusion dominated. This solution corresponds to the stable point in the physical space. At around $R=10.5$, the first bifurcation of the stationary flow is detected with the formation of the limited cycle; see Figure 14.

Figure 14.

Projection of the limited cycle $R=10.5$ and invariant torus (with Poincaré section) $R=510$ to three‐dimensional phase subspace.

After that the influence of viscosity compressed the region of parameters forming chaotic attractors around the cycle with small amplitude. These attractors are, supposedly, low dimensional singular cycles. But the small amplitude of the solution over the cycle made it impossible to analyze. The next attractor that we are able to detect is the limited torus, presented in Figure 14 with the Poincaré section indicated. Close resemblance to the cycle may indicate that the attractor was formed as the result of multistability from the cycle by Hopf bifurcation. This indicates the existence of two irrational frequencies in the system. Further increase of the Reynolds number up to $R=516$ resulted in the other Hopf bifurcation with the formation of the three‐dimensional invariant torus, presented in Figure 15. This torus becomes singular (by period doubling bifurcations on one of the frequencies). This is shown in Figure 16 for $R=518$. However, this cascade of period doubling bifurcations is reversed to the original 3D torus. The next bifurcation that could be traced at $R=520.5$ is another Hopf bifurcation leading to the formation of the four‐dimensional invariant torus.

Figure 15.

Projection of the invariant three‐dimensional torus and first Poincaré section for $R=516$ to the three‐dimensional phase subspace and second section for $R=518$.

Figure 16.

Projection of the invariant four‐dimensional torus into three‐dimensional phase subspace (top left) and sequential first, second, and third sections in the phase space for $R=520.5$ (left to right, up to down). Only parts of sections are shown.

Further increase of the Reynolds number leads to the chaotic solution that corresponds to the dense field of points in phase subspace projections up to about $R=2100$. With further increase of $R$, we observed the formation of inverse cascades. For verification purposes, the domain length was increased again up to $L_x=16$ for fixed $R=2500$. This resulted in the chaotic solution.

5. Discussion and conclusion

5.1. Discussion

This research for the bifurcation scenarios in RT and KH instabilities started in 2013 and continues still. From this research, we obtained information that raises some questions that must be addressed. The first question is the stability of the KH and RT problems for low Reynolds number. It is known [38] that viscosity may result in instability. However, we found two ways of this instability to occur. Starting with very low Reynolds number results in linear stability of the flow. When linear stability is violated, we can witness a rapidly growing instability that results in highly suppressed cascade of bifurcations. We are still unable to deduce this cascade. However, from a certain value of $R$, the solution becomes regular with the formation of regular attractor (either invariant cycle or torus). This indicates that the cascade of bifurcations that is induced by diffusion is reversed at some point resulting in the final return to the attractor with the maximum basin of attraction. Those are the bifurcations that were traced by other authors [48–50].

Another topic of discussion is the Landau‐Hopf scenario that was found in coupled problem. This means that the system triggers secondary instabilities that have irrational frequencies. It is known, say from [68], that the existence of secondary oscillations on KH instability is triggered by density stratification. This results in the formation of instabilities of outer filament or inner vortex type. These instabilities might have irrational frequencies resulting in the formation of high dimensional attractor. That is clearly the case for coupled problem; see Figure 12.

Another difficulty is related to the limited length of the domain during presented bifurcation analysis. We show in Section 4.1, that the information for subsonic flow is translated back to the upwind direction. This results in the global phase space formation. On the other hand, the length of the domain becomes a bifurcation parameter since physical process evolves with the background flow and is cut at the point of flow leaving the domain. The situation for the initial valued problem is different since the process will eventually relax (e.g., densities are swapped for RT instability), and phase space solution will be represented by a fixed stable point. It means that the process consists of direct and inverse cascades of bifurcations starting and finishing at the point. In most cases, the initial point is linearly unstable and the final point is stable.

5.2. Conclusion

In this paper, we show some results regarding RT and KH instabilities and phase space properties for supersonic solution with shock waves. We show that the phase space is separated for inviscid flow and give proof of this fact for 1D gas flow. This fact is numerically demonstrated, and some results are obtained that confirm the validity of this proposition for multidimensional case. We also considered linear stability and bifurcation analysis for KH and RT instabilities in the setup with the coaxial flow.

We introduce the following notations of regular attractors to write down bifurcation scenarios: $P$ is a point, $Cn$ is a limited cycle of period $n$, $nTM$ is an invariant torus of dimension $M$ and period $n$ (on any frequency). For fixed values of domain length, we obtained the following bifurcation scenarios.

KH instability with 2D perturbations:

$P\to C\to \mathrm{\dots }\to C5\to \mathrm{\dots }\to C11\to \mathrm{\dots }\to C$E31

KH instability with 3D perturbations:

$P\to C\to T2\to T3\to Chaos\mathrm{.}$E39

Coupled problem (KH and RT instability) with 2D perturbations:

$P\to C\to C2\to C4\to Cn\to \mathrm{..}\to C5\to C3\to Chaos\mathrm{.}$E32

Coupled problem (KH and RT instability) with 3D perturbations:

$P\to C\to nC\to \mathrm{\dots }\to C\to T2\to T3\to \mathrm{\dots }\to nT3\to \mathrm{\dots }\to T4\to Chaos\mathrm{.}$E33

The following scenarios show the following. First, the existence of FShM scenario in KH instability is confirmed by Eqs. (31) and (32). In Eq. (31), we notice the direct and inverse cascades with the Sharkovskiy sequence are not fully developed. In Eq (32), we are able to show the full FShM cascade up to limited cycle of period three. However, these low dimensional attractors are only the initial trigger mechanism of laminar‐turbulent transition in such complicated problems. It can be clearly seen in Eq. (33) where a Landau‐Hopf scenario is developed up to a regular attractor represented by the 4D invariant torus. The existence of computationally stable 4D invariant torus is a remarkable fact. It took $2.6\cdot {10}^9$ time samples to analyze and about 3.5 months to calculate. We can notice the appearance of singular 3D tori in the cascade. Most likely, some part of incomplete FShM scenario with direct and inverse cascades is involved in these singular tori. It is interesting to note that we found no resonance tori during the bifurcation analysis that are typical for other problems.