Open access peer-reviewed chapter

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

Written By

Nikolay Mihaylovitch Evstigneev and Nikolai Alexandrovitch Magnitskii

Submitted: 09 November 2016 Reviewed: 16 February 2017 Published: 26 July 2017

DOI: 10.5772/67918

From the Edited Volume

Turbulence Modelling Approaches - Current State, Development Prospects, Applications

Edited by Konstantin Volkov

Chapter metrics overview

2,023 Chapter Downloads

View Full Metrics


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.


  • laminar‐turbulent transition
  • bifurcations in fluid dynamics
  • compressible flow bifurcations
  • Kelvin‐Helmholtz instability
  • Rayleigh‐Taylor instability
  • dynamical systems
  • direct numerical simulation
  • eigenvalue solver

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 α=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., [1012]. 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 [1518] 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. [2022]. 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 [2931] 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:u1,ρ1 and for y<0:u2,ρ2. Gravity vector is g=(0,g,0)T, directed downward. The analysis gives the following stability condition:


where k is a wavenumber. It is clear, that if u1u2 then there exists such k=kcrit, that the flow is unstable and KH instability always occurs (including cases with ρ1=ρ2 and g=0). If u1=u2, then RT instability occurs if ρ1>ρ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 [3234] 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 (λ=2π/k, λcrit=2π/kcrit) for instability [35] to be:


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:


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 [3739], 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=γP/ρ, where P is pressure, ρ is density, and γ 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=u1u2. 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 ys, that u(ys)=0. If the flow is unstable, then u(uu(ys))<0 for some point y[ymin,ymax]) 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 (θ(k)) of infinitesimal perturbations [29]:




where α 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 μ in the model results in the formation of maximum growth rate for specific k(μ) with θ(k)0 with k. 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, 4346]. For infinite domains, the growth rate θ(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 γ increases, θ(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 et|k|. In addition, solutions that grow arbitrarily quickly in the Sobolev space W2k, 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):


Let xMX, t is time, M is a phase space, and ϕt(x):M×RM is a phase flow on M and R is a vector of parameters (e.g., similarity complexes in case of fluid dynamics) of size Rp. Let γ(t)M be a phase trajectory of a solution of Eq. (6).

Definition 1. Compact set BM, invariant with respect to ϕt is an attracting set if there exists a neighborhood U of the set B such that for almost all xU, ϕt(x)B,t.

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 u0(x,t) or the parameter value R0:


where L:MM is a linear operator and


If the spectrum of operator L lies in the left complex half‐plane when jRp,RjR0,j, and if j,Rj>R0,j with one or few its eigenvalues that are moving in the right half‐plane, then the system (6) has a bifurcation at R0.

Definition 7. The infinite sequence of bifurcations: jRp,Rj,nRj,*,n 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]:


Here, we are considering some bounded domain ΩR3 with boundary Ω; scalar functions f are defined as f:Ω×[0,T]R; vector‐functions f are defined as f:Ω×[0,T]R3, where T is some defined finite time. Then, ρE is a scalar function of the full gas energy; e is a scalar function of the internal gas energy; γ is the adiabatic index (value 1.4 is used for calculations); p is a scalar pressure function; u is a gas velocity vector function; ρ 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:


where μ is a dynamic gas viscosity and strain rate tensor is defined as:


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:


where LC is the characteristic length scale of the problem, ||uC|| is the norm of the characteristic velocity vector, ρC is a characteristic density. With this in mind, we can rewrite Eq. (8) as follows:


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 δ. For all results presented here, we use δ=5107.

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 dim(K)=km=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 Ω=[L×H] and initial mixing layer thickness δ. 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 δ. Assuming a wavenumber k of the least stable mode, the length of the domain in X‐direction would result as L=2πδ/k. The domain length in Y‐direction is chosen as H=8L.


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

The solution is performed on 32×256 grid with k=8 and m=5 IRA parameters, that results in the dim(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.

MComputed αAnalytical α [33]α [41]α [42]

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×4096 grid points. The problem is set as KH instability of shear flow in rectangular unit domain Ω with periodic boundary conditions in all directions. Initial conditions are set as:




Here Um=(U1U2)/2, ρm=(ρ1ρ2)/2, ρ1=1, ρ2=2, U1=1/2, U2=1/2, L=0.025.

For our simulation, we use a 3D version of the problem and set up grid 256×256×256 such that uz(x,0)=0 with periodic boundary conditions in z‐direction. The problem is executed until 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 ky(t)=maxxΩρuy2/2, and we monitor the growth rate of the amplitude for the unstable uy 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 α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 αy and maximum y‐direction kinetic energy ky. 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., [6167] 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 θ(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 θB~2πg/kt, and growth rate of spikes is proportional to θJ~gt2; (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 Lm as a function of time with Lm=αgAgt2, where α is a constant. We consider the domain Ω=[3×1×1] with gravity direction in x, and we use 900×300×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: τLeft=0.3s, τRight=1.9s.

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 α for sinusoidal perturbations: αshort=0.0271 for short wavelength and α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:


The system can be linearized:


with eigenvalues being:


where c=γpρ 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=(u0,0)T and n=(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:


where k is a wavenumber.

Lemma 1: Let all eigenvalues λjR of A in Eq. (17) satisfy λj0. Then, u1 and u2 are constant for x<0,tR+.

Proof: Let linear operator A=(abcd); a,b,c,dR. Eigenvalues and eigenvectors are found as:


Let us assume that λ1>0,λ20. We apply characteristic variables W0=H1U0:


Matrix A is diagonalized as Λ=H1AH. Then, we can rewrite (17) as:


and the last vector equation can be separated:


Then characteristic solution is:


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




Now we turn to conservative variables U=HW:


It is clear from the obtained solution that eigenvectors in Eq. (21) are independent of xλit. So, the solution can be reformed as:


where η1,η2 are values in Eq. (21) that are independent of x and t: η1=2bda+(da)2+4cb η2=2bda(da)2+4cb. Let us consider two possible cases:

  1. λ1>0,λ2=0


  2. λ1>0,λ2>0,λ1λ2


After some cumbersome but simple manipulations one can show that η1W1+η2W2=1 and W1+W2=1. {a,b,c,d}R. □

Lemma 2: Let there exist at least one eigenvalues λ0 of A from Eq. (17), that λ0<0. Then values of u1 or u2 for x<0 are not constant for tR+.

Proof: For the proof, we are considering one example and then use the proof of Lemma 1. Let λ1>0,λ2<0 by the Lemma formulation. Let A=(1101). Then, λ=(11) and H=(12110). Using characteristic variables:


Diagonalizing operator A: Λ=H1AH=((01112)(1101))(12110)=diag(λi) and arriving at two independent equations. The solution in characteristic variables is:


Returning to the original variables:

(u1u2)=({ 12;(xt)<012eik(xt);(xt)0+{ 12;(x+t)<0112eik(x+t);(x+t)0 { 1;(xt)<0eik(x+t);(xt)0 ).E27

Let us consider the case only for x<0.


We now turn to random real variables in matrix A that satisfy the condition λ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 A and only depends on signs of eigenvalues of A. □

Lemma 3: Assume that matrix A is diagonalizable. Then Lemmas 1 and 2 are valid for AR3×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=x0 and consider variables on the left (L) and on the right (R) from this line. We assume that UL is an unperturbed constant vector and that UR contains perturbations. We also assume that pL<pR and that there is a stationary shock wave with shock speed S=0 at line x=x0. In this case, we can conclude the following:

Proposition 1: Let there be such solution to the system (15), that at some line x=x0, there is a stationary shock wave with pL<pR and U(x,0),xx0 is unperturbed. Then, U(x,t),xx0 remains unperturbed at tR+.

Proof: In accordance with the condition of the proposition, the local solution of Eq. (15) at x0 is a shock wave solution that is located in the first characteristic field, i.e., λ1=uct, since pL<pR. The second characteristic field is a contact discontinuity with λ1=u with the second contact discontinuity being parallel to the flow, since n=(1,0)T. The third characteristic field contains an expansion wave, λ3=u+ct. Then, λ1=0, since S=0, and so the local Jacobi matrix eigenvalues are ordered as:λ1=0λ2λ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 Ω, sized 8×3×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 Ω\B at initial conditions with ρ=1.24,u=(1,0,0)T,M0=2. Wall boundary conditions are set for all Ω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 μ=0. We impose two unphysical boundary conditions: inflow boundary ΩI on the left plane yz for x=0 (ρ0=1.24,u=(1,0,0)T,M=2) and outflow boundary Ω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×200×200 elements. Phase space portraits are constructed using 10106 time points after the flow is quasi‐stationary. The viscosity coefficient is set to 1105. 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 Ω=[8×1×1], the flow direction is from left to right. The following boundary conditions are imposed: Ω(0,y,z)=Ω0 is the inflow boundary condition, Ω(8,y,z)=Ω1 is the outflow boundary condition, Ω(x,y,0)=Ω(x,y,1)=Ω2 is the periodic boundary condition in z‐direction, and Ω(x,0,z)=Ω(x,1,z)=Ω3 is the non‐reflecting (symmetry) boundary condition. We use simple ghost cell boundary conditions for Ω2 and Ω3. For Ω1, we use appropriately chosen absorbing boundary conditions with characteristic far‐field treatment and we use characteristic boundary conditions for Ω0. The initial conditions are identical to the boundary conditions on Ω0 and are given below:


The perturbations are set depending on the problem setup, either 2D (in this case Cz=0) or 3D (in this case Cz=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×250×250 elements totaling 62.5106 grid points. The sponge zone for Ω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 R157 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 (+1,0i) 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 (Lx) is another parameter. With the increase of the domain up to 16, the first critical Reynolds number gradually changes to Rc=25. It can be explained by the wavelength in the x‐direction. It is interesting to note that using Lx=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 uy 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 τ can be roughly estimated as advection time scale τa~Lx/|ux| for advection‐dominated flow. Providing that, this time is smaller for Lx=4, initial perturbations for inviscid flow are dominated by advection timescale rather than diffusion timescale τd~Lx2/μ=RLx2. The timescales are equal if R~1Lx|ux|, so with the increase of Lx, 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 Lx. 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 λcritical~12π. 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 ux results in the change of the bifurcation scenario. This again can be explained by the domination of τa with the increase of R. We confirmed this proposition by increasing the Lx up to 16 for the fixed value of R=2500 thus giving the diffusion scale more time to develop. While cycle was formed for Lx=8,R=2500, the chaotic solution emerges for Lx=16,R=2500, so the value of Lx 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=760,Lx=8) and 3D invariant torus (R=950,Lx=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 Ω0 and initial values are set as follows:


The Richardson number (Eq. (3)) is set depending on the problem in range of 1/51/8, where H=4δ is the length of the mixing layer and ρmax=1.2. We also consider two options—Cz=0 or Cz=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, Lx=6. The increase of Lx 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 υz,υy phase subspace for 2D perturbation problem for R=2417,Ri=1/6,Lx=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 Lx 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 Lx=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 [4850].

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:


KH instability with 3D perturbations:


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


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


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.6109 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.


  1. 1. Evstigneev N.M., Magnitskii N.A. (2012) Chapter 10: FSM scenarios of laminar‐turbulent transition in incompressible fluids. In Nonlinearity, bifurcation and chaos – theory and applications. Intech, Rijeka, pp. 251–280.
  2. 2. Charles Li  Y. (2007) Chaos in partial differential equations, Navier‐Stokes equations and turbulence. ICCM. III, 110–122
  3. 3. Arnold V.I., Meshalkin L.D. (1958–1959) The seminar of A.N. Kolmogorov on selected topics in analysis. Usp. Mat. Nauk. 15, 247–250
  4. 4. Meshalkin L.D., Sinai I.G. (1961) Investigation of the stability of a stationary solution of a system of equations for the plane movement of an incompressible viscous liquid. J. Appl. Math. Mech. 25(4), 1700–1705
  5. 5. Yudovich V.I. (1973) Natural oscillations arising from loss of stability in parallel flows of a viscous liquid under long‐wavelength periodic disturbances. Fluid Dyn. 8, 26. doi:10.1007/BF01017632
  6. 6. Sivashinsky G. (1985) Weak turbulence in periodic flows. Phys. D. 17(2), 243–255
  7. 7. Lucas D., Kerswell R. (2014) Spatiotemporal dynamics in two‐dimensional Kolmogorov flow over large domains. J. Fluid Mech. 750, 518–554. doi:10.1017/jfm.2014.270
  8. 8. Evstigneev N.M., Magnitskii N.A., Silaev D.A. (2015) Qualitative analysis of dynamics in Kolmogorov’s problem on a flow of a viscous incompressible fluid. Differ. Eqn. 51(10), 1292–1305
  9. 9. Lord Rayleigh O.M. F.R.S. (1916) On convection currents in a horizontal layer of fluid, when the higher temperature is on the under side. Philos. Mag. Ser. 6(32), 192, 529–546
  10. 10. Manneville P. (2010) Rayleigh‐Benard convection: thirty years of experimental, theoretical, and modeling work. Springer Tracts Modern Phys. 207, 41–65
  11. 11. Paul Supriyo, Verma Mahendra K., Wahi Pankaj, Reddy Sandeep K., Kumar Krishna (2012) Bifurcation analysis of the flow patterns in two‐dimensional Rayleigh‐Benard convection. Int. J. Bifurc. Chaos. 22(5), 1230018
  12. 12. Evstigneev N.M. (2016) Laminar‐turbulent bifurcation scenario in 3D Rayleigh‐Benard convection problem. Open J. Fluid Dyn. 6, 496–539. doi:10.4236/ojfd.2016.64035
  13. 13. Getling A.V. (1998) Rayleigh‐Benard convection: structures and dynamics. World Scientific, Singapore. ISBN 978‐981‐02‐2657‐2
  14. 14. Ashwin P., Podvigina O. (2003) Hopf bifurcation with cubic symmetry and instability of ABC flow. Proc. R. Soc Lond. A. 459(2035), 1801–1827
  15. 15. Pfister G., Schmidt H., Cliffet K.A., Mullin T. (1988) Bifurcation phenomena in Taylor‐Couette flow in a very short annulus. J. Fluid Mech. 191, 1–18
  16. 16. Mamun C.K., Tuckerman L.S. (1994) Asymmetry and Hopf bifurcation in spherical Couette Flow. Phys. Fluids. 7(1). doi:10.1063/1.868730
  17. 17. Lopez J.M., Marques F. (2005) Finite aspect ratio Taylor–Couette flow: Shil’nikov dynamics of 2‐tori. Phys. D. 211, 168–191
  18. 18. Ma  Tian, Wang Shouhong (2006) Stability and bifurcation of the Taylor problem. Arch. Ration. Mech. Anal. 181, 149–176
  19. 19. Furukawa H., Horikoshi H., Ohazama N., Watanabe T. (2012, June 25–28) PIV analysis of mode bifurcation in Taylor vortex flow. In 15th International Symposium on Flow Visualization, Minsk, Belarus
  20. 20. Barkley D., Gomes M.G.M., Henderson R.D. (2002) Three‐dimensional instability in flow over a backward‐facing step. J. Fluid Mech. 473, 167–190. doi:10.1017/S002211200200232X
  21. 21. Rani H.P., Sheu Tony W.H. (2006) Nonlinear dynamics in a backward‐facing step flow. Phys. Fluids. 18, 084101
  22. 22. Evstigneev N.M., Magnitskii N.A., Sidorov S.V. (2009) On the nature of turbulence in a problem on the motion of a fluid behind a ledge. Differ. Eqn. 45(1), 68–72
  23. 23. Soibelman I., Meiron D.I. (1991) Finite‐amplitude bifurcations in plane Poiseuille flow: two‐dimensional Hopf bifurcation. J. Fluid Mech. 229, 389–416
  24. 24. Casas P.S., Jorba A. (2012) Hopf bifurcations to quasi‐periodic solutions for the two‐dimensional plane Poiseuille flow. Commun. Nonlinear Sci. Numer. Simul. 17(7), 2864
  25. 25. Lord Rayleigh O.M. F.R.S. (1883) Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density. Proc. Lond. Math. Soc. XIV, 70–177
  26. 26. Taylor G.I. (1950) The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. Proc. R. Soc. Lond. 201, 192–196
  27. 27. von Helmholtz Hermann (1868) Uber discontinuierliche Flussigkeits‐Bewegungen. Monatsberichte der Koniglichen Preussische Akademie der Wissenschaften zu Berlin. 23, 215–228
  28. 28. Lord Kelvin (William Thomson). (1871) Hydrokinetic solutions and observations. Philos Mag. 42, 362–377
  29. 29. Drazin P.G. (2002) Introduction to hydrodynamic stability. Cambridge University Press, Cambridge, United Kingdom
  30. 30. Kull H.J. (1991) Theory of the Rayleigh‐Taylor instability. Phys. Rep. 206(5), 197–325. North‐Holland
  31. 31. Book D.L. (1984) Rayleigh‐Taylor instability in compressible media. NRL memorandum report no. 5373. Washington, DC: Naval Research Laboratory
  32. 32. Schmid P.J., Henningson D.S. (2001) Stability and transition in shear flows. Springer, Springer-Verlag New York, lnc.
  33. 33. Blumen W. (1970) Shear layer instability of an inviscid compressible fluid. J. Fluid Mech. 40(Part 4), 769–781
  34. 34. Hesla Todd I., Pranckh Ferdinand R., Preziosi Luigi (1986) Squire’s theorem for two stratified fluids. Phys. Fluids 29(9), 2808–2811. doi:10.1063/1.865478
  35. 35. Cushman‐Roisin B. (2005) Kelvin‐Helmholtz instability as a boundary‐value problem. Environ Fluid Mech. 5, 507–525
  36. 36. Richardson L.F. (1920) The supply of energy from and to atmospheric eddies. Proc. R. Soc. London. A97, 354–373
  37. 37. Kumar K., Tuckerman L.S. (1994) Parametric instability of the interface between two fluids. J. Fluid Mech. 279, 49–68
  38. 38. Funada T., Joseph D.D. (2001) Viscous potential flow analysis of Kelvin‐Helmholtz instability in a channel. J. Fluid Mech. 445, 263–283
  39. 39. Yoshikawa H., Wesfreid J.E. (2011) Oscillatory Kelvin‐Helmholtz instability. Part 1. A viscous theory. J. Fluid Mech. 675, 223–248
  40. 40. Fjørtoft T. (1950) Application of integral theorems in deriving criteria of instability of laminar flow and for the baroclinic circular vortex. Geophys. Publ. 17, 1–52
  41. 41. Mack C. (2009) Global stability of compressible flow about a swept parabolic body. PhD thesis, Ecole Polytechnique
  42. 42. Henze O., Lemke M., Sesterhenn J. (2015) A parallel and matrix free framework for global stability analysis of compressible flows. arXiv:1502.03701
  43. 43. Bernstein I.B., Book D.L. (1983) Effect of compressibility on the Rayleigh‐Taylor instability. Phys. Fluids 26, 453
  44. 44. Lezzi A.M., Prosperetti A. (1989) Rayleigh‐Taylor instability for adiabatically stratified fluids. Phys. Fluids A. 1, 1784
  45. 45. Turner L. (2002) Rayleigh‐Taylor instabilities and gravity waves in compressible fluids. Los Alamos National Laboratory Report LA‐UR‐02‐6439, Los Alamos, USA.
  46. 46. Livescu D. (2004) Compressibility effects on the Rayleigh‐Taylor instability growth between immiscible fluids. Phys. Fluids. 16(1), 118–127
  47. 47. Guo Yan, Tice Ian (2009) Compressible, inviscid Rayleigh‐Taylor instability. Indiana Univ. Math. J. 60(2). doi:10.1512/iumj.2011.60.4193
  48. 48. Poehlmann A., Richter R., Rehberg I. (2013) Unravelling the Rayleigh‐Taylor instability by stabilization. J. Fluid Mech. 732, R3. doi:10.1017/jfm.2013.424
  49. 49. Ilak M., Schlatter P., Bagheri S., Henningson D.S. (2012) Bifurcation and stability analysis of a jet in cross‐flow: onset of global instability at a low velocity ratio. J. Fluid Mech. 686, 94–121
  50. 50. Hunter J.K., Thoo J.B. (2011) On the weakly nonlinear Kelvin–Helmholtz instability of tangential discontinuities in MHD. J. Hyperbol. Differ. Eqn. 8(4), 691–726
  51. 51. Magnitskii N.A., Sidorov S.V. (2006) New methods for chaotic dynamics. Singapore: World Scientific
  52. 52. Magnitskii N.A. (2012) Chapter 6: Universality of transition to chaos in all kinds of nonlinear differential equations. In Nonlinearity, bifurcation and chaos – theory and applications. Intech, pp. 133–174.
  53. 53. Toro E.F. (1999) Riemann solvers and numerical methods for fluid dynamics. Springer‐Verlag, Berlin-Heidelberg.
  54. 54. Evstigneev N.M. (2016) On the construction and properties of WENO schemes order five, seven, nine, eleven and thirteen. Part 1. Construction and stability. Comput. Res. Model. 8(5), 721–753. (in Russian)
  55. 55. Evstigneev N.M., Magnirskii N.A. (2012) On phase space peculiarities of gas dynamics equations for a supersonic initial‐boundary value problem. Proc. ISA RAS. 62(4), 85–102. (in Russian)
  56. 56. Evstigneev N.M., Ryabkov O.I. (2016) Application of multi GPU+CPU architecture for the direct numerical simulation of laminar‐turbulent transition. Numer. Methods Program. 17, 55–64. (in Russian)
  57. 57. Evstigneev N.M., Magnitskii N.A. (2014) On the initial stage of Kelvin‐Helmholtz instability for laminar‐turbulent transition in viscous gas flow. Proc. ISA RAS. 64(3), 41–52. (in Russian)
  58. 58. Evstigneev N.M. (2016) On the construction and properties of WENO schemes order five, seven, nine, eleven and thirteen. Part 2. Comput. Res. Model. 8(6), 885–910. (in Russian)
  59. 59. Evstigneev N.M. (2017) Implementation of implicitly restarted Arnoldi method on multiGPU architecture with application to fluid dynamics problems. Commun. Comput. Inform. Sci. (accepted for publication)
  60. 60. McNally C.P., Lyra W., Passy J.C. (2012) A well‐posed Kelvin‐Helmholtz instability test and comparison. Astrophys. J. Suppl. Ser. 201(2), 1–17
  61. 61. Youngs  D.L. (1984) Numerical simulation of turbulent mixing by Rayleigh‐Taylor instability. Phys. D. 12, 32–44
  62. 62. Dimonte G, Youngs D.L., et al. (2004) A comparative study of the turbulent Rayleigh‐Taylor instability using high‐resolution three‐dimensional numerical simulations: the alpha‐group collaboration. Phys. Fluids. 16, 1668–1693
  63. 63. Kucherenko, Yu. A. et al. (2003) Experimental investigation into the self‐similar mode of mixing of different density gases in the Earth’s gravitational field. Laser Part. Beams. 21, 385–388
  64. 64. Cabot W.H., Cook A.W. (2006) Reynolds number effects on Rayleigh‐Taylor instability with possible implications for type Ia supernovae. Nat. Phys. 2(8), 562–568
  65. 65. Belotserkovskii O.M., Oparin A.M. (2000) Numerical experiment on turbulence: from order to chaos. Moscow: Nauka. (in Russian)
  66. 66. Sin’kova O.G., et al. (1997) Three‐dimensional DNS of gravitational turublence mixing. In 6th International Workshop on the Physics of Compressible Turbulent Mixing, Marseille, France, pp. 470–479
  67. 67. Yanilkin Y.V., Statsenko V.P., Kozlov V.I. (2009) Mathematical simulation of turbulent mixing in compressible fluids. Nauka, Sarov, Russia.
  68. 68. Lecoanet D., McCourt M., Quataert E., Burns K.J., Vasil G.M., Oishi J.S., Brown B.P., Stone J.M., O’Leary R.M. (2015) A validated nonlinear Kelvin‐Helmholtz benchmark for numerical hydrodynamics. Mon. Not. R. Astron. Soc. 455(4), 4274–4288

Written By

Nikolay Mihaylovitch Evstigneev and Nikolai Alexandrovitch Magnitskii

Submitted: 09 November 2016 Reviewed: 16 February 2017 Published: 26 July 2017