## 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

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

where

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

where

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

The inviscid incompressible RT instability gives exponential growth of the amplitude (

with

where *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

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

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

Definition 1. *Compact set* *, invariant with respect to* *is an attracting set if there exists a neighborhood* *of the set* *such that for almost all*

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

where

If the spectrum of operator

Definition 7. *The infinite sequence of bifurcations:* *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]:

(7) |

Here, we are considering some bounded domain

where

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

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

### 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

## 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

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

The solution is performed on **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.

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 |

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

(13) |

and

Here

For our simulation, we use a 3D version of the problem and set up grid **Figure 2**. One can see a very reasonable agreement with reference data. Small difference in the evolution of growth rate

### 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

We follow the problem setup in [67] and investigate the mixing length **Figure 3**.

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

## 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:

(15) |

The system can be linearized:

with eigenvalues being:

where

where

Lemma 1: *Let all eigenvalues* *of* *in* Eq. (17) *satisfy* *. Then,* *and* *are constant for*

Proof: Let linear operator

(18) |

Let us assume that

(19) |

Matrix

and the last vector equation can be separated:

Then characteristic solution is:

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

Here:

Now we turn to conservative variables

It is clear from the obtained solution that eigenvectors in Eq. (21) are independent of

where

(24) |

After some cumbersome but simple manipulations one can show that

Lemma 2: *Let there exist at least one eigenvalues* *of* *from* Eq. (17), *that* *. Then values of* *or* *for* *are not constant for*

Proof: For the proof, we are considering one example and then use the proof of Lemma 1. Let

Diagonalizing operator

Returning to the original variables:

(27) |

Let us consider the case only for

We now turn to random real variables in matrix

Lemma 3: *Assume that matrix* *is diagonalizable. Then Lemmas 1 and 2 are valid for*

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

Proposition 1: *Let there be such solution to the system* (15)*, that at some line* *there is a stationary shock wave with* *and* *is unperturbed. Then,* *remains unperturbed at*

Proof: In accordance with the condition of the proposition, the local solution of Eq. (15) at

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

We are performing the simulation for quasi‐stationary solution on the grid of **Figure 5**.

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

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

The perturbations are set depending on the problem setup, either 2D (in this case

#### 4.2.1. Eigenvalue analysis

First, we analyze linear stability of the main flow for the described problem setup. For **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 **Figure 9**. It was found that the length of the domain in the

Let us consider time scales in the problem. The maximum physical process time **Figure 7**) demonstrates about

#### 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

#### 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 **Figure 10**.

### 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

(30) |

The Richardson number (Eq. (3)) is set depending on the problem in range of **Figure 11**. Interesting to note that this cycle was found in the reverse cascade due to the limited length,

#### 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

The starting point is **Figure 14**.

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 **Figure 15**. This torus becomes singular (by period doubling bifurcations on one of the frequencies). This is shown in **Figure 16** for

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

## 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

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:

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