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

Nikolay Mihaylovitch Evstigneev and

Nikolai Alexandrovitch Magnitskii

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/67918

#### Abstract

[27] Bhushan S, Walters DK. A dynamic hybrid RANS/LES modeling framework. Physics of

[28] Lilly DK. A proposed modification of the Germano subgrid-scale closure method. Phys-

[29] Spalart P, Allmaras S. A One-Equation Turbulence Model for Aerodynamic Flows. AIAA

[30] Menter F. Review of the shear-stress transport turbulence model experience from an industrial perspective. International Journal of Computational Fluid Dynamics. 2009;23

[31] Boris JP, Grinstein FF, Oran ES, Kolbe RL. New insights into large eddy simulation. Fluid

[32] Patankar SV, Spalding DB. A calculation procedure for heat, mass, and momentum transfer in three-dimensional parabolic flows. International Journal of Heat Mass Trans-

[33] Barth TJ, Jespersen DC. The Design and Application of Upwind Schemes on Unstruc-

[34] Rhie CM. and Chow WL. Numerical study of the turbulent flow past an airfoil with

[35] Jasak H, Weller HG, Gosman AD. High resolution NVD differencing scheme for arbitrarily unstructured meshes. International Journal for Numerical Methods in Fluids.

[36] Adedoyin AA, Walters DK, Bhushan S. Assessment of modeling and discretization error in finite-volume large eddy simulations. ASME Paper No. IMECE2006-14918. 2006 [37] Breuer M, Jovicic N, Mazaev K. Comparison of DES, RANS and LES for the separated flow around a flat plate at high incidence. International Journal for Numerical Methods in

Fluids. 2012;24:No. 015103

ics of Fluids. 1992;4:pp. 633–635

Paper No. AIAA-92-0439. 1992

fer. 1972;15:pp. 1787–1806

1999;31:pp. 431–449

Fluids. 2003;41:pp. 357–388

Dynamics Research. 1992;10(4–6):pp. 199–228

28 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

tured Seshes. AIAA Paper No. AIAA-89-0366. 1989

trailing edge separation. AIAA Journal. 1983;21:pp. 1525–1532

(4):pp. 305–316

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.

Keywords: laminar-turbulent transition, bifurcations in fluid dynamics, compressible flow bifurcations, Kelvin-Helmholtz instability, Rayleigh-Taylor instability, dynamical systems, direct numerical simulation, eigenvalue solver

© 2017 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## 1. Introduction and definitions

#### 1.1. Introduction

In the chapter, we continue the work started in 2006 with the numerical analysis of laminarturbulent 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., [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].

1. Introduction and definitions

30 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

review of the progress in this field.

1.2. Problem overview

In the chapter, we continue the work started in 2006 with the numerical analysis of laminarturbulent 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

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., [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

1.1. Introduction

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 : u1, ρ<sup>1</sup> and for y < 0 : u2, ρ2. Gravity vector is g ¼ ð0; � g; 0Þ T, directed downward. The analysis gives the following stability condition:

$$k^2 \rho\_1 \rho\_2 (u\_1 - u\_2)^2 < k \lg(\rho\_2^2 - \rho\_1^2). \tag{1}$$

where k is a wavenumber. It is clear, that if u<sup>1</sup> 6¼ u<sup>2</sup> then there exists such k ¼ kcrit, that the flow is unstable and KH instability always occurs (including cases with ρ<sup>1</sup> ¼ ρ<sup>2</sup> and g ¼ 0). If u<sup>1</sup> ¼ u2, then RT instability occurs if ρ<sup>1</sup> > ρ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 initialcondition 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:

$$\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} \tag{2}$$

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(\mu\_1 - \mu\_2)^2} < \frac{1}{4'} \tag{3}$$

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 ¼ γ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 ¼ u<sup>1</sup> � u2. 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″ <sup>ð</sup>ysÞ ¼ 0. If the flow is unstable, then <sup>u</sup>″ ðu � uð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]:

$$
\Theta(k) \sim e^{a(\rho,k,\mathbf{g})t}, a(\rho,k,\mathbf{g}) = \sqrt{A(\rho)\mathbf{g}k} \tag{4}
$$

with

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

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.

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

> Ri <sup>¼</sup> gHðρ<sup>2</sup> � <sup>ρ</sup>1<sup>Þ</sup> ρ0ðu<sup>1</sup> � u2Þ

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

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 ¼ u<sup>1</sup> � u2. 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

<sup>ð</sup>ysÞ ¼ 0. If the flow is unstable, then <sup>u</sup>″

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

The inviscid incompressible RT instability gives exponential growth of the amplitude (θðkÞ) of

<sup>A</sup> <sup>¼</sup> <sup>ρ</sup><sup>2</sup> � <sup>ρ</sup><sup>1</sup> ρ<sup>2</sup> þ ρ<sup>1</sup>

, αðρ, k, gÞ ¼

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi AðρÞgk

q

θðkÞ � e

αðρ, k, gÞt

, ð3Þ

ðu � uðysÞÞ < 0 for some

, ð4Þ

, ð5Þ

linearly in the mixing height H), resulting in the following necessary condition:

32 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

viscosity.

there exists a point ys, that u″

are presented in [33, 41, 42].

infinitesimal perturbations [29]:

with

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 θð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 <sup>e</sup><sup>t</sup> ffiffiffi jkj p . In addition, solutions that grow arbitrarily quickly in the Sobolev space W<sup>k</sup> 2, 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 laminarturbulent 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 \mathbf{u}(\mathbf{x},t)}{\partial t} = \mathbf{F}(\mathbf{u}, \mathbf{x}, \mathbf{R}, t). \tag{6}$$

Let x∈ M ∈ X, t is time, M is a phase space, and ϕ<sup>t</sup> ðxÞ : M · R ! M 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 B∈ M, invariant with respect to ϕ<sup>t</sup> is an attracting set if there exists a neighborhood U of the set B such that for almost all x ∈ U, ϕ<sup>t</sup> ð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:

$$\frac{\partial \mathbf{y}(\mathbf{x},t)}{\partial t} = \mathbf{L}(\mathbf{u}\_{0\prime} \ \mathbf{x}, t \ \ \mathbf{R}\_0) \mathbf{y}(\mathbf{x}, t),$$

where L : M ! M is a linear operator and

$$\mathbf{y}(\mathbf{x},t) = \mathbf{u}\_0(\mathbf{x},t) - \mathbf{u}(\mathbf{x},t).$$

If the spectrum of operator L lies in the left complex half-plane when ∀j∈ Rp, Rj ≤R<sup>0</sup>;j, and if ∃j, Rj > R<sup>0</sup>;<sup>j</sup> 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: ∀j∈ Rp, Rj,n ! Rj,�, 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

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.

Let us introduce some definitions and notations that are required for the analysis of laminarturbulent 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

is a vector of parameters (e.g., similarity complexes in case of fluid dynamics) of size Rp. Let

Definition 1. Compact set B∈ M, invariant with respect to ϕ<sup>t</sup> is an attracting set if there exists a

Definition 3. Regular attractors of the system (6) are stable singular points, stable invariant limited

Definition 6. Bifurcation points of the system (6) are those and only those points in which there is no

<sup>∂</sup><sup>t</sup> <sup>¼</sup> <sup>L</sup>ðu0, <sup>x</sup>, t, <sup>R</sup>0Þyðx, tÞ,

yðx, tÞ ¼ u0ðx, tÞ � uðx, tÞ:

If the spectrum of operator L lies in the left complex half-plane when ∀j∈ Rp, Rj ≤R<sup>0</sup>;j, and if ∃j, Rj > R<sup>0</sup>;<sup>j</sup> with one or few its eigenvalues that are moving in the right half-plane, then the

Definition 7. The infinite sequence of bifurcations: ∀j∈ Rp, Rj,n ! Rj,�, n ! ∞ is known as the

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

<sup>∂</sup><sup>t</sup> <sup>¼</sup> <sup>F</sup>ðu, <sup>x</sup>, <sup>R</sup>, tÞ: <sup>ð</sup>6<sup>Þ</sup>

ðxÞ ! B, t ! ∞.

ðxÞ : M · R ! M is a phase flow on M and R

∂uðx, tÞ

34 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

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

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

Linearize the system (4) on its regular attractor u0ðx, tÞ or the parameter value R0:

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

continuous dependence of the phase portrait of the system by its parameters.

∂yðx, tÞ

Let x∈ M ∈ X, t is time, M is a phase space, and ϕ<sup>t</sup>

cycles, and stable invariant tori of any dimension.

where L : M ! M is a linear operator and

system (6) has a bifurcation at R0.

cascade of bifurcations.

γðtÞ∈ M be a phase trajectory of a solution of Eq. (6).

neighborhood U of the set B such that for almost all x ∈ U, ϕ<sup>t</sup>

1.3. Definitions

general):

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

$$\begin{split} \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_{j}} \left[ \rho u\_{i} \right] &= 0 \\ \frac{\partial}{\partial t} (\rho u\_{i}) + \frac{\partial}{\partial \mathbf{x}\_{j}} \left[ \rho u\_{i} u\_{j} + p \delta\_{ij} - \tau\_{ji} \right] &= \rho \mathbf{g}\_{i'} \quad \text{ (i, j, k)} = 1, 2, 3; \\ \frac{\partial}{\partial t} (\rho E) + \frac{\partial}{\partial \mathbf{x}\_{j}} \left[ \rho u\_{j} E + u\_{j} p - u\_{i} \tau\_{ij} \right] &= \rho u\_{k} \mathbf{g}\_{k'} \\ \rho E &= \frac{1}{2} \rho \mathbf{u}^{2} + \rho \mathbf{e}, \\ p &= (\boldsymbol{\gamma} - 1)(E - 1/2 \rho \mathbf{u}^{2}). \end{split} \tag{7}$$

Here, we are considering some bounded domain Ω⊂R<sup>3</sup> with boundary ∂Ω; scalar functions f are defined as <sup>f</sup> : <sup>Ω</sup> · <sup>½</sup>0;<sup>T</sup> � ! <sup>R</sup>; vector-functions <sup>f</sup> are defined as <sup>f</sup> : <sup>Ω</sup> · <sup>½</sup>0;<sup>T</sup> � ! <sup>R</sup>3, where <sup>T</sup> 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:

$$
\pi\_{\vec{\eta}} = 2\mu \mathbf{S}\_{\vec{\eta}\lambda} \tag{8}
$$

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

$$S\_{i\bar{j}} = \frac{1}{2} \left( \frac{\partial u\_i}{\partial \mathbf{x}\_{\bar{j}}} + \frac{\partial u\_{\bar{j}}}{\partial \mathbf{x}\_{i}} \right) - \frac{1}{3} \delta\_{i\bar{j}} \frac{\partial u\_k}{\partial \mathbf{x}\_k} \,. \tag{9}$$

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\_{\mathbb{C}} ||\mathbf{u}\_{\mathbb{C}}|| \rho\_{\mathbb{C}}}{\mu} \,\,\,\,\tag{10}$$

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

$$
\pi\_{ij} = \frac{2}{R} \mathcal{S}\_{ij}.\tag{11}
$$

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 <sup>δ</sup>. For all results presented here, we use <sup>δ</sup> <sup>¼</sup> <sup>5</sup> � <sup>10</sup>�<sup>7</sup> .

#### 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Þ ¼ k � m ¼ 80.

## 3. Verification

where LC is the characteristic length scale of the problem, jjuCjj is the norm of the characteristic velocity vector, ρ<sup>C</sup> is a characteristic density. With this in mind, we can rewrite Eq. (8) as

<sup>τ</sup>ij <sup>¼</sup> <sup>2</sup>

We will use Eq. (8) if we use the term "viscosity coefficient" and we use Eq. (11) if we use the

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

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

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

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

.

<sup>R</sup> Sij: <sup>ð</sup>11<sup>Þ</sup>

follows:

term "Reynolds number".

2.2. Numerical methods

2.2.1. Main flow solver

2.2.2. Eigenvalue solver

is dimðKÞ ¼ k � m ¼ 80.

ness <sup>δ</sup>. For all results presented here, we use <sup>δ</sup> <sup>¼</sup> <sup>5</sup> � <sup>10</sup>�<sup>7</sup>

36 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

other bifurcations of the main solution cannot be derived.

"burglization" numerical experiments, provided in Ref. [58].

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.

$$\begin{aligned} \rho\_0 &= 1.24, \\ \iota(\mathbf{x}, \mathbf{y})\_0 = \iota(\mathbf{y})\_0 &= \tanh\left(\frac{\mathbf{y} - \mathrm{L}/2}{\delta}\right). \end{aligned} \tag{12}$$

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.


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

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

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:

$$\rho(\mathbf{x},0) = \begin{cases} \rho\_1 - \rho\_m e^{\frac{y-1/4}{L}}; 1/4 > y \ge 0, \\\rho\_2 + \rho\_m e^{\frac{y-1/4}{L}}; 1/2 > y \ge 1/4, \\\rho\_2 + \rho\_m e^{-\frac{y/4-y}{L}}; 3/4 > y \ge 1/2, \\\rho\_1 - \rho\_m e^{-\frac{y-1/4}{L}}; 1 > y \ge 3/4, \\\rho\_1 - \rho\_m e^{\frac{y-1/4}{L}}; 1/4 > y \ge 0, \\\end{cases} \tag{13}$$

$$u\_x(\mathbf{x},0) = \begin{cases} \mathcal{U}\_1 - \mathcal{U}\_m e^{\frac{y-1/4}{L}}; 1/4 > y \ge 0, \\\mathcal{U}\_2 + \mathcal{U}\_m e^{\frac{y-1/4}{L}}; 1/2 > y \ge 1/4, \\\mathcal{U}\_2 + \mathcal{U}\_m e^{\frac{-3/4}{L}}; 3/4 > y \ge 1/2, \\\mathcal{U}\_1 - \mathcal{U}\_m e^{\frac{y-3/4}{L}}; 1 > y \ge 3/4 \end{cases}$$

and

$$
\mu\_y(\mathbf{x},0) = 0.01 \sin\left(4\pi\mathbf{x}\right), \newline p(\mathbf{x}, y, 0) = 2.5. \tag{14}
$$

Here Um ¼ ðU<sup>1</sup> � U2Þ=2, ρ<sup>m</sup> ¼ ðρ<sup>1</sup> � ρ2Þ=2, ρ<sup>1</sup> ¼ 1, ρ<sup>2</sup> ¼ 2, U<sup>1</sup> ¼ 1=2, U<sup>2</sup> ¼ �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Þ ¼ max<sup>x</sup><sup>∈</sup> <sup>Ω</sup> <sup>ρ</sup>u<sup>2</sup> <sup>y</sup>=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 α<sup>y</sup> 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 α<sup>y</sup> and maximum y-direction kinetic energy ky. Ref1 is [60] and Ref2 is data from spectral code.

## 3.2. RT instability verification

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 Ω

ρ<sup>1</sup> � ρme

<sup>ρ</sup><sup>2</sup> <sup>þ</sup> <sup>ρ</sup>me�3=4�<sup>y</sup>

U<sup>1</sup> � Ume

<sup>U</sup><sup>2</sup> <sup>þ</sup> Ume�3=4�<sup>y</sup>

Here Um ¼ ðU<sup>1</sup> � U2Þ=2, ρ<sup>m</sup> ¼ ðρ<sup>1</sup> � ρ2Þ=2, ρ<sup>1</sup> ¼ 1, ρ<sup>2</sup> ¼ 2, U<sup>1</sup> ¼ 1=2, U<sup>2</sup> ¼ �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 perturba-

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 α<sup>y</sup> 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 α<sup>y</sup> and maximum y-direction kinetic energy ky. Ref1 is [60] and Ref2 is data from spectral code.

<sup>U</sup><sup>1</sup> � Ume�y�3=<sup>4</sup>

U<sup>2</sup> þ Ume

<sup>ρ</sup><sup>1</sup> � <sup>ρ</sup>me�y�3=<sup>4</sup>

ρ<sup>2</sup> þ ρme

y�1=4

y�1=4

�yþ1=4

�yþ1=4

<sup>L</sup> ; 1=4 > y ≥ 0;

<sup>L</sup> ; 1=2 > y ≥ 1=4;

<sup>L</sup> ; 3=4 > y ≥ 1=2;

<sup>L</sup> ; 1 > y ≥ 3=4;

<sup>L</sup> ; 1=4 > y ≥ 0;

ð13Þ

<sup>L</sup> ; 1=2 > y ≥ 1=4;

<sup>L</sup> ; 3=4 > y ≥ 1=2;

<sup>L</sup> ; 1 > y ≥ 3=4

uyðx; 0Þ ¼ 0:01 sin ð4πxÞ, pðx, y; 0Þ ¼ 2:5: ð14Þ

<sup>y</sup>=2, and we monitor the growth rate of the amplitude for the

with periodic boundary conditions in all directions. Initial conditions are set as:

8 >>>>>>>><

>>>>>>>>:

8 >>>>>>><

>>>>>>>:

ρðx; 0Þ ¼

38 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

uxðx; 0Þ ¼

and

tions, as kyðtÞ ¼ max<sup>x</sup><sup>∈</sup> <sup>Ω</sup> <sup>ρ</sup>u<sup>2</sup>

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 θð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~ ffiffiffiffiffiffiffiffiffiffiffiffiffi 2πg=k p t, 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 <sup>¼</sup> <sup>α</sup>gAgt2, where <sup>α</sup> is a constant. We consider the domain <sup>Ω</sup> ¼ ½<sup>3</sup> · <sup>1</sup> · <sup>1</sup>� 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.

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

Figure 3. Density for RT instability test case with A ¼ 1=3 and harmonic initial disturbances: τLef t ¼ 0:3 s, τRight ¼ 1:9 s.

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

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.

## 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{aligned} \frac{\partial \mathbf{U}}{\partial t} + \frac{\partial \mathbf{F}}{\partial \mathbf{x}} + \frac{\partial \mathbf{G}}{\partial y} &= \mathbf{0}, \\ \mathbf{U} &= \left[\rho; \rho u\_{\mathbf{x}}; \rho u\_{\mathbf{y}}; \rho \boldsymbol{E}\right]^T, \\ \mathbf{F} &= \left[\rho u\_{\mathbf{x}}; \rho u\_{\mathbf{x}}^2 + p; \rho u\_{\mathbf{x}} u\_{\mathbf{y}}; (\rho \boldsymbol{E} + p) u\_{\mathbf{x}}\right]^T, \\ \mathbf{G} &= \left[\rho u\_{\mathbf{y}}; \rho u\_{\mathbf{y}} u\_{\mathbf{x}}; \rho u\_{\mathbf{y}}^2 + p; (\rho \boldsymbol{E} + p) u\_{\mathbf{y}}\right]^T. \end{aligned} \tag{15}$$

The system can be linearized:

$$\begin{cases} \frac{\partial \mathbf{U}}{\partial t} + \mathbf{A} \frac{\partial \mathbf{U}}{\partial x} + \mathbf{B} \frac{\partial \mathbf{U}}{\partial y} = \mathbf{0}, \\ \mathbf{A} = \partial \mathbf{F} / \partial \mathbf{U}; \mathbf{B} = \partial \mathbf{G} / \partial \mathbf{U}, \end{cases} \tag{16}$$

with eigenvalues being:

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

where c ¼ ffiffiffiffi γp ρ q 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Þ <sup>T</sup> and <sup>n</sup>¼ ð1; <sup>0</sup><sup>Þ</sup> <sup>T</sup>. 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{aligned} \mathbf{U}\_t + \mathbf{A} \mathbf{U}\_x &= 0, \\ \mathbf{x} \in \mathbb{R}; t \in \mathbb{R}^+, \\ \mathbf{U} &= \begin{pmatrix} u\_1, u\_2 \end{pmatrix}^T, \\ u\_1(\mathbf{x}, 0) = 1, u\_2(\mathbf{x}, 0) = \begin{cases} e^{i\mathbf{k}\mathbf{x}}; \mathbf{x} \ge 0 \\ 1; \mathbf{x} < 0. \end{cases} \end{aligned} \tag{17}$$

where k is a wavenumber.

Lemma 1: Let all eigenvalues λ<sup>j</sup> ∈ R of A in Eq. (17) satisfy λ<sup>j</sup> ≥ 0. Then, u<sup>1</sup> and u<sup>2</sup> are constant for x < 0;∀t ∈ R<sup>þ</sup>.

Proof: Let linear operator <sup>A</sup> <sup>¼</sup> a b c d � �; a, b, c, d<sup>∈</sup> <sup>R</sup>. Eigenvalues and eigenvectors are found as:

$$\begin{aligned} \lambda\_{1,2} &= \begin{pmatrix} 1/2d + 1/2a + 1/2\sqrt{(d-a)^2 + 4cb} \\ 1/2d + 1/2a - 1/2\sqrt{(d-a)^2 + 4cb} \end{pmatrix}, \\ H\_{1,2} &= \begin{pmatrix} 2b \\ d - a + \sqrt{(d-a)^2 + 4cb} \\ 1 & d - a - \sqrt{(d-a)^2 + 4cb} \\ 1 & 1 \end{pmatrix} \end{aligned} \tag{18}$$

Let us assume that <sup>λ</sup><sup>1</sup> <sup>&</sup>gt; <sup>0</sup>;λ<sup>2</sup> <sup>≥</sup> 0. We apply characteristic variables <sup>W</sup><sup>0</sup> <sup>¼</sup> <sup>H</sup>�<sup>1</sup> U0:

$$
\begin{pmatrix} W\_{0,1} \\ W\_{0,2} \end{pmatrix} = \begin{pmatrix} \begin{cases} \frac{1}{4} (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^{ik})}{b\sqrt{(d - a)^2 + 4cb}}; x \ge 0. \\\\ \begin{cases} \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^{ik})}{b\sqrt{(d - a)^2 + 4cb}}; x \ge 0. \end{cases} \end{pmatrix} \tag{19}$$

Matrix <sup>A</sup> is diagonalized as <sup>Λ</sup> <sup>¼</sup> <sup>H</sup>�<sup>1</sup> AH. Then, we can rewrite (17) as:

$$H^{-1}\mathbf{U}\_t + H^{-1}A\mathbf{U}\_x = \mathbf{0};$$

$$\mathbf{A} = \operatorname{diag}(\lambda\_i);$$

$$\mathbf{W}\_t + \mathbf{A}\mathbf{W}\_x = \mathbf{0},$$

and the last vector equation can be separated:

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

Then characteristic solution is:

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.

40 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

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

<sup>F</sup> ¼ ½ρux; <sup>ρ</sup>u<sup>2</sup>

∂U ∂t þ

<sup>G</sup> ¼ ½ρuy; <sup>ρ</sup>uyux; <sup>ρ</sup>u<sup>2</sup>

∂U <sup>∂</sup><sup>t</sup> <sup>þ</sup> <sup>A</sup> <sup>∂</sup><sup>U</sup> ∂x

consider a simple initial value problem for system of advection equations:

∂F ∂x þ ∂G <sup>∂</sup><sup>y</sup> <sup>¼</sup> <sup>0</sup>,

U ¼ ½ρ; ρux; ρuy; ρE�

T,

<sup>y</sup> þ p;ðρE þ pÞuy�

λ<sup>1</sup> ¼ u � n � c, λ<sup>2</sup> ¼ u � n, λ<sup>3</sup> ¼ u � n, λ<sup>4</sup> ¼ u � n þ c, ð36Þ

is the speed of sound, n is a unit directional vector. All eigenvalues for all x, y

T,

ð15Þ

ð16Þ

ð17Þ

T:

<sup>T</sup>. In this case, the system of eigenvectors

<sup>x</sup> þ p; ρuxuy;ðρE þ pÞux�

<sup>þ</sup> <sup>B</sup> <sup>∂</sup><sup>U</sup> <sup>∂</sup><sup>y</sup> <sup>¼</sup> <sup>0</sup>,

A ¼ ∂F=∂U; B ¼ ∂G=∂U,

are real since system (16) is hyperbolic. We are considering a special case with the base flow

becomes one dimensional, and we can reduce Eq. (16) to one-dimensional case. Let us first

U<sup>t</sup> þ AU<sup>x</sup> ¼ 0; x ∈ R; t∈ R<sup>þ</sup>, U ¼ ðu1, u2Þ

<sup>u</sup>1ðx; <sup>0</sup>Þ ¼ <sup>1</sup>;u2ðx; <sup>0</sup>Þ ¼ <sup>e</sup>ikx; x <sup>≥</sup> <sup>0</sup>

Lemma 1: Let all eigenvalues λ<sup>j</sup> ∈ R of A in Eq. (17) satisfy λ<sup>j</sup> ≥ 0. Then, u<sup>1</sup> and u<sup>2</sup> are constant for

T,

�

1; x < 0:

,

<sup>T</sup> and <sup>n</sup>¼ ð1; <sup>0</sup><sup>Þ</sup>

4. Bifurcation analysis

The system can be linearized:

with eigenvalues being:

ffiffiffiffi γp ρ q

where k is a wavenumber.

x < 0;∀t ∈ R<sup>þ</sup>.

being parallel to x axis, i.e., u ¼ ðu0; 0Þ

where c ¼

4.1. Phase space for supersonic flow 4.1.1. Inviscid flow phase space properties

$$\mathcal{W}\_{i} = \mathcal{W}\_{0,i}(\mathfrak{x} - \lambda\_{i}t), i = 1, 2.$$

#### 42 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

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

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

$$\begin{cases} A - \because \left( x - \frac{d + a - \sqrt{\left(d - a\right)^2 + 4cb}}{2}t \right) < 0. \\ B - \therefore \left( x - \frac{d + a - \sqrt{\left(d - a\right)^2 + 4cb}}{2}t \right) \ge 0. \end{cases}$$

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

Here:

$$B\pm = \frac{1}{4} \frac{\left(d - a \pm \sqrt{\left(d - a\right)^2 + 4cb}\right) \left(-d + a \pm \sqrt{\left(d - a\right)^2 + 4cb} + 2b \cdot e\right)^{ik} \left(e^{\frac{k + a \pm \sqrt{\left(d - a\right)^2 + 4cb}}{2}}\right)}{b\sqrt{\left(d - a\right)^2 + 4cb}},\tag{20}$$

Now we turn to conservative variables U ¼ HW:

$$
\psi \begin{pmatrix} u\_1 \\ u\_2 \end{pmatrix} = \left( \frac{2b}{d - a + \sqrt{\left(d - a\right)^2 + 4cb}} W\_1 + \frac{2b}{d - a - \sqrt{\left(d - a\right)^2 + 4cb}} W\_2 \right). \tag{21}
$$

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

$$
\begin{pmatrix} u\_1 \\ u\_2 \end{pmatrix} = \begin{pmatrix} \eta\_1 \begin{Bmatrix} W\_1 \\ \mathcal{W}\_1(e^{\vec{\text{rk}}(\mathbf{x}-\lambda\_1 t)}) \end{Bmatrix} + \eta\_2 \begin{Bmatrix} W\_2 \\ \mathcal{W}\_2(e^{\vec{\text{rk}}(\mathbf{x}-\lambda\_2 t)}) \end{Bmatrix} \\ \begin{Bmatrix} W\_1 \\ \mathcal{W}\_1(e^{\vec{\text{rk}}(\mathbf{x}-\lambda\_1 t)}) \end{Bmatrix} + \begin{Bmatrix} W\_2 \\ \mathcal{W}\_2(e^{\vec{\text{rk}}(\mathbf{x}-\lambda\_2 t)}) \end{Bmatrix} \end{pmatrix} \tag{22}
$$

where <sup>η</sup>1, <sup>η</sup><sup>2</sup> are values in Eq. (21) that are independent of <sup>x</sup> and <sup>t</sup>: <sup>η</sup><sup>1</sup> <sup>¼</sup> <sup>2</sup><sup>b</sup> d�aþ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðd�aÞ 2 þ4cb' p <sup>η</sup><sup>2</sup> <sup>¼</sup> <sup>2</sup><sup>b</sup> d�a� ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðd�aÞ 2 <sup>þ</sup>4cb <sup>p</sup> . Let us consider two possible cases:

#### Numerical Analysis of Laminar‐Turbulent Bifurcation Scenarios in Kelvin‐Helmholtz… http://dx.doi.org/10.5772/67918 43

## 1. λ<sup>1</sup> > 0;λ<sup>2</sup> ¼ 0

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

0

{

BBBBBBBBBBBBBBBBBBBB@

{

d � a � þ

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðd � aÞ

� � q

Now we turn to conservative variables U ¼ HW:

u1 u2 � �

¼

η1

0

BB@

p . Let us consider two possible cases:

d � a þ

<sup>2</sup> <sup>þ</sup> <sup>4</sup>cb

<sup>A</sup>� ¼ <sup>1</sup> 4

d � a �

u1 u2 � �

¼

So, the solution can be reformed as:

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðd�aÞ 2 þ4cb

0 B@

Here:

B� ¼ 1 4

<sup>η</sup><sup>2</sup> <sup>¼</sup> <sup>2</sup><sup>b</sup> d�a�

<sup>A</sup> <sup>þ</sup> ; <sup>x</sup> � <sup>d</sup> <sup>þ</sup> <sup>a</sup> <sup>þ</sup>

<sup>B</sup> <sup>þ</sup> ; <sup>x</sup> � <sup>d</sup> <sup>þ</sup> <sup>a</sup> <sup>þ</sup>

<sup>A</sup> � ; <sup>x</sup> � <sup>d</sup> <sup>þ</sup> <sup>a</sup> �

<sup>B</sup> � ; <sup>x</sup> � <sup>d</sup> <sup>þ</sup> <sup>a</sup> �

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðd � aÞ

<sup>2</sup> <sup>þ</sup> <sup>4</sup>cb

b

�d þ a �

b

B@

2b

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðd � aÞ

<sup>2</sup> <sup>þ</sup> <sup>4</sup>cb <sup>q</sup> <sup>W</sup><sup>1</sup> <sup>þ</sup>

> W<sup>1</sup> <sup>W</sup>1ðeikðx�λ1t<sup>Þ</sup>

where <sup>η</sup>1, <sup>η</sup><sup>2</sup> are values in Eq. (21) that are independent of <sup>x</sup> and <sup>t</sup>: <sup>η</sup><sup>1</sup> <sup>¼</sup> <sup>2</sup><sup>b</sup>

� W<sup>1</sup> <sup>W</sup>1ðeikðx�λ1t<sup>Þ</sup>

� �

0 @

42 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

0 @

0 @

> 0 @

� � q

W<sup>1</sup> W<sup>2</sup> � �

q

q

q

q

¼

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðd � aÞ

<sup>2</sup> <sup>t</sup>

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðd � aÞ

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðd � aÞ

<sup>2</sup> <sup>t</sup>

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðd � aÞ

�d þ a �

<sup>2</sup> <sup>þ</sup> <sup>4</sup>cb

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðd � aÞ

<sup>2</sup> <sup>þ</sup> <sup>4</sup>cb

d � a �

� W<sup>2</sup> <sup>W</sup>2ðeikðx�λ2t<sup>Þ</sup>

<sup>2</sup> <sup>þ</sup> <sup>4</sup>cb

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðd � aÞ

q

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðd � aÞ

W<sup>1</sup> þ W<sup>2</sup>

<sup>Þ</sup> <sup>þ</sup> <sup>η</sup><sup>2</sup>

Þ þ

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

<sup>2</sup> <sup>t</sup>

<sup>2</sup> <sup>t</sup>

<sup>2</sup> <sup>þ</sup> <sup>4</sup>cb

<sup>2</sup> <sup>þ</sup> <sup>4</sup>cb

<sup>2</sup> <sup>þ</sup> <sup>4</sup>cb

<sup>2</sup> <sup>þ</sup> <sup>4</sup>cb

q

1 A < 0: 1

CCCCCCCCCCCCCCCCCCCCA

1 A ≥ 0:

1 A < 0:

1 A ≥ 0:

� �

<sup>q</sup> ,

0 � �

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðd � aÞ

þ 2b � e

<sup>q</sup> ,

2b

W<sup>2</sup> <sup>W</sup>2ðeikðx�λ2t<sup>Þ</sup>

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðd � aÞ

Þ

1

Þ

<sup>2</sup> <sup>þ</sup> <sup>4</sup>cb <sup>q</sup> <sup>W</sup><sup>2</sup>

<sup>2</sup> <sup>þ</sup> <sup>4</sup>cb

þ 2b

ik x�dþa�

ffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>ð</sup>d�aÞ2þ4cb p

1

CCA, <sup>ð</sup>22<sup>Þ</sup>

d�aþ

p

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðd�aÞ 2 þ4cb'

CA: <sup>ð</sup>21<sup>Þ</sup>

<sup>2</sup> t

1 CA

ð20Þ

$$\begin{aligned} \text{1.} x \ge \lambda\_1. \Rightarrow \begin{pmatrix} u\_1 \\ u\_2 \end{pmatrix} = \begin{pmatrix} \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{pmatrix}. \\\ \text{2.} x > 0, x < \lambda\_1. \Rightarrow \begin{pmatrix} u\_1 \\ u\_2 \end{pmatrix} = \begin{pmatrix} \eta\_1 W\_1 + \eta\_2 W\_2(e^{ikx}) \\ W\_1 + W\_2(e^{ikx}) \end{pmatrix}. \\\ \text{3.} x \le 0. \Rightarrow \begin{pmatrix} u\_1 \\ u\_2 \end{pmatrix} = \begin{pmatrix} \eta\_1 W\_1 + \eta\_2 W\_2 \\ W\_1 + W\_2 \end{pmatrix}. \end{aligned} \tag{23}$$

2. λ<sup>1</sup> > 0;λ<sup>2</sup> > 0;λ<sup>1</sup> ≥ λ<sup>2</sup>

$$\begin{aligned} 1. & \ge \lambda\_1. \Rightarrow \begin{pmatrix} u\_1\\ u\_2 \end{pmatrix} = \begin{pmatrix} \eta\_1 W\_1(e^{ik(x-\lambda\_1t)}) + \eta\_2 W\_2(e^{ik(x-\lambda\_2t)})\\ W\_1(e^{ik(x-\lambda\_1t)}) + W\_2(e^{ik(x-\lambda\_2t)}) \end{pmatrix}. \\\ 2. & \lambda\_1 \ge x \ge \lambda\_2. \Rightarrow \begin{pmatrix} u\_1\\ u\_2 \end{pmatrix} = \begin{pmatrix} \eta\_1 W\_1 + \eta\_2 W\_2(e^{ik(x-\lambda\_2t)})\\ W\_1 + W\_2(e^{ik(x-\lambda\_2t)}) \end{pmatrix}. \\\ 3. & \ge 0, \mathbf{x} < \lambda\_2. \Rightarrow \begin{pmatrix} u\_1\\ u\_2 \end{pmatrix} = \begin{pmatrix} \eta\_1 W\_1 + \eta\_2 W\_2\\ W\_1 + W\_2 \end{pmatrix}. \\\ 4. & \ge 0, \Rightarrow \begin{pmatrix} u\_1\\ u\_2 \end{pmatrix} = \begin{pmatrix} \eta\_1 W\_1 + \eta\_2 W\_2\\ W\_1 + W\_2 \end{pmatrix}. \end{aligned} \tag{24}$$

After some cumbersome but simple manipulations one can show that η1W<sup>1</sup> þ η2W<sup>2</sup> ¼ 1 and <sup>W</sup><sup>1</sup> <sup>þ</sup> <sup>W</sup><sup>2</sup> <sup>¼</sup> 1. <sup>∀</sup>{a, <sup>b</sup>, <sup>c</sup>, <sup>d</sup>} <sup>∈</sup> <sup>R</sup>. □

Lemma 2: Let there exist at least one eigenvalues λ<sup>0</sup> of A from Eq. (17), that λ<sup>0</sup> < 0. Then values of u<sup>1</sup> or u<sup>2</sup> for x < 0 are not constant for t∈ R<sup>þ</sup>.

Proof: For the proof, we are considering one example and then use the proof of Lemma 1. Let <sup>λ</sup><sup>1</sup> <sup>&</sup>gt; <sup>0</sup>;λ<sup>2</sup> <sup>&</sup>lt; 0 by the Lemma formulation. Let <sup>A</sup> <sup>¼</sup> �1 1 0 1 � �. Then, <sup>λ</sup> <sup>¼</sup> <sup>1</sup> �1 � � and H ¼ 1 2 1 1 0 0 @ 1 A. Using characteristic variables:

$$\begin{pmatrix} W\_{0,1} \\ W\_{0,2} \end{pmatrix} = \begin{pmatrix} 1; \ge < 0 \\ e^{ik\boldsymbol{\varepsilon}}; \boldsymbol{\varepsilon} \ge 0 \\\\ \begin{pmatrix} 1 \\ \boldsymbol{2} \end{pmatrix} \boldsymbol{\varepsilon} \boldsymbol{x} < 0 \\\\ 1 - \frac{1}{2} e^{ik\boldsymbol{\varepsilon}}; \boldsymbol{x} \ge 0 \end{pmatrix} \tag{25}$$

Diagonalizing operator <sup>A</sup>: <sup>Λ</sup> <sup>¼</sup> <sup>H</sup>�<sup>1</sup> <sup>A</sup><sup>H</sup> <sup>¼</sup> 0 1 <sup>1</sup> � <sup>1</sup> 2 ! � �1 1 0 1 ! � � � 1 2 1 1 0 0 @ 1 A ¼ diagðλiÞ

and arriving at two independent equations. The solution in characteristic variables is:

$$\begin{pmatrix} W\_1 \\ W\_2 \end{pmatrix} = \begin{pmatrix} 1; (\mathbf{x} - t) < 0 \\ e^{i\mathbf{k}(\mathbf{x} - t)}; (\mathbf{x} - t) \ge 0 \\\\ \begin{pmatrix} 1 \\ 2 \end{pmatrix}; (\mathbf{x} + t) < 0 \\ 1 - \frac{1}{2} e^{i\mathbf{k}(\mathbf{x} + t)}; (\mathbf{x} + t) \ge 0 \end{pmatrix} . \tag{26}$$

Returning to the original variables:

$$\begin{pmatrix} u\_1 \\ u\_2 \end{pmatrix} = \begin{pmatrix} \begin{cases} & \frac{1}{2}; (\mathbf{x} - t) < 0 \\\\ \frac{1}{2} e^{ik(\mathbf{x} - t)}; (\mathbf{x} - t) \ge 0 \end{cases} + \begin{cases} & \frac{1}{2}; (\mathbf{x} + t) < 0 \\\\ 1 - \frac{1}{2} e^{ik(\mathbf{x} + t)}; (\mathbf{x} + t) \ge 0 \end{cases} \\\\ & \begin{cases} & 1; (\mathbf{x} - t) < 0 \\\\ e^{ik(\mathbf{x} + t)}; (\mathbf{x} - t) \ge 0 \end{cases} \end{pmatrix}. \tag{27}$$

Let us consider the case only for x < 0.

$$\begin{aligned} \text{1.} & \succeq (\mathbf{x} + t) ; \mathbf{x} < 0. \Rightarrow \begin{pmatrix} u\_1 \\ u\_2 \end{pmatrix} = \begin{pmatrix} \mathbf{\hat{3}} - \frac{1}{2} e^{i\mathbf{k} (\mathbf{x} + t)} \\ 1 \end{pmatrix}. \\ & \mathbf{2.} < (\mathbf{x} + t) ; \mathbf{x} < 0. \Rightarrow \begin{pmatrix} u\_1 \\ u\_2 \end{pmatrix} = \begin{pmatrix} 1 \\ 1 \end{pmatrix}. \end{aligned} \tag{28}$$

We now turn to random real variables in matrix A that satisfy the condition ∃λ<sup>i</sup> < 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 A ∈ R<sup>3</sup> · <sup>3</sup> .

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<sup>0</sup> 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 U<sup>R</sup> 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Þ, ∀x ≤ x<sup>0</sup> is unperturbed. Then, Uðx, tÞ, ∀x ≤ x<sup>0</sup> remains unperturbed at ∀t ∈ R<sup>þ</sup>.

Proof: In accordance with the condition of the proposition, the local solution of Eq. (15) at x<sup>0</sup> is a shock wave solution that is located in the first characteristic field, i.e., λ<sup>1</sup> ¼ u � ct, since pL < pR. The second characteristic field is a contact discontinuity with λ<sup>1</sup> ¼ u with the second contact discontinuity being parallel to the flow, since n¼ ð1; 0Þ <sup>T</sup>. The third characteristic field contains an expansion wave, λ<sup>3</sup> ¼ u þ ct. Then, λ<sup>1</sup> ¼ 0, since S ¼ 0, and so the local Jacobi matrix eigenvalues are ordered as:λ<sup>1</sup> ¼ 0 ≤ λ<sup>2</sup> ≤ λ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

Diagonalizing operator <sup>A</sup>: <sup>Λ</sup> <sup>¼</sup> <sup>H</sup>�<sup>1</sup>

Returning to the original variables:

u1 u2 � �

Let us consider the case only for x < 0.

¼

� W<sup>1</sup> W<sup>2</sup> � ¼

44 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

1 2

1:x ≥ ðx þ tÞ; x < 0: )

2:x < ðx þ tÞ; x < 0: )

1 2 e ikðx�tÞ

8 >><

0

BBBBBBBBBBBB@

>>:

<sup>A</sup><sup>H</sup> <sup>¼</sup> 0 1

eikðx�t<sup>Þ</sup>

1 2

<sup>1</sup> � <sup>1</sup> 2 e ikðxþtÞ

and arriving at two independent equations. The solution in characteristic variables is:

(

;ðx � tÞ < 0

;ðx � tÞ ≥ 0

(

þ

eikðxþt<sup>Þ</sup>

u1 u2 !

We now turn to random real variables in matrix A that satisfy the condition ∃λ<sup>i</sup> < 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. □

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<sup>0</sup> 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 U<sup>R</sup> 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:

Lemma 3: Assume that matrix A is diagonalizable. Then Lemmas 1 and 2 are valid for A ∈ R<sup>3</sup> · <sup>3</sup>

8 >><

>>:

1;ðx � tÞ < 0

¼

u1 u2 !

0

BBBBBBB@

<sup>1</sup> � <sup>1</sup> 2

( <sup>1</sup>;ð<sup>x</sup> � <sup>t</sup><sup>Þ</sup> <sup>&</sup>lt; <sup>0</sup>

;ðx � tÞ ≥ 0

;ðx þ tÞ ≥ 0

1 2

<sup>1</sup> � <sup>1</sup> 2 e ikðxþtÞ

;ðx � tÞ ≥ 0

3 2 � 1 2 e ikðxþtÞ

0 @

1

<sup>¼</sup> <sup>1</sup> 1 ! :

;ðx þ tÞ < 0

� �1 1 0 1

1

CCCCCCCA

;ðx þ tÞ < 0

;ðx þ tÞ ≥ 0

1 A: �

1

: ð26Þ

1

CCCCCCCCCCCCA

: ð27Þ

ð28Þ

.

A ¼ diagðλiÞ

0 @

! � �

!

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, M<sup>0</sup> <sup>¼</sup> 2. Wall boundary conditions are set for all ∂Ω<sup>W</sup> 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 ∂Ω<sup>I</sup> on the left plane yz for x ¼ 0 (ρ<sup>0</sup> ¼ 1:24;u¼ ð1; 0; 0Þ T, M<sup>∞</sup> <sup>¼</sup> 2) and outflow boundary <sup>∂</sup>Ω<sup>O</sup> 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 10 � 106 time points after the flow is quasi-stationary. The viscosity coefficient is set to 1 � <sup>10</sup>�<sup>5</sup> . 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.

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

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

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Þ ¼ ∂Ω<sup>0</sup> is the inflow boundary condition, ∂Ωð8;y, zÞ ¼ ∂Ω<sup>1</sup> is the outflow boundary condition, ∂Ωðx, y, 0Þ ¼ ∂Ωðx, y, 1Þ ¼ ∂Ω<sup>2</sup> is the periodic boundary condition in z-direction, and ∂Ωðx, 0, zÞ ¼ ∂Ωðx, 1, zÞ ¼ ∂Ω<sup>3</sup> is the nonreflecting (symmetry) boundary condition. We use simple ghost cell boundary conditions for ∂Ω<sup>2</sup> 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 ∂Ω<sup>0</sup> and are given below:

$$\begin{aligned} \rho|\_{\partial\Omega\_0} &= 1.0, \\ u\_x|\_{\partial\Omega\_0} = 3/2 + 1/2 \tanh\left(\frac{y - 1/2}{\delta} + \mathcal{C}\_z \cos\left(6\pi z\right)\right), \\ u\_y|\_{\partial\Omega\_0} &= 1 \cdot 10^{-5} \sin\left(6\pi y\right), \\ u\_z|\_{\partial\Omega\_0} &= 0, \\ p|\_{\partial\Omega\_0} &= 22.14. \end{aligned} \tag{29}$$

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 · <sup>250</sup> · 250 elements totaling 62:<sup>5</sup> � 106 grid points. The sponge zone for <sup>∂</sup>Ω<sup>1</sup> 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

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

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.

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

46 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

First, we analyze linear stability of the main flow for the described problem setup. For R ≤ 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 ðþ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 τ<sup>a</sup> � Lx=juxj 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 <sup>τ</sup><sup>d</sup> � <sup>L</sup><sup>2</sup> <sup>x</sup>=<sup>μ</sup> <sup>¼</sup> RL<sup>2</sup> <sup>x</sup>. The timescales are equal if <sup>R</sup> � <sup>1</sup> Lxjuxj , 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

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

48 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

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

initial instability.

limited cycle solution, R ¼ 400.

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 τ<sup>a</sup> 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 ∂Ω<sup>0</sup> and initial values are set as follows:

$$\begin{split} |\rho|\_{\partial\Omega\_{0}} &= 1.1 + \frac{1}{10} \text{erf}\left(\frac{y - 1/2}{\delta} + \mathbb{C}\_{z} \cos\left(6\pi z\right)\right), \\ |u\_{x}|\_{\partial\Omega\_{0}} &= 2.0 + \text{erf}\left(\frac{y - 1/2}{\delta} + \mathbb{C}\_{z} \cos\left(6\pi z\right)\right), \\ |p|\_{\partial\Omega\_{0}} &= 22.14 + \rho g y + (\rho\_{\text{max}} - \rho)g. \end{split} \tag{30}$$

The Richardson number (Eq. (3)) is set depending on the problem in range of 1=5 � 1=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.

Numerical Analysis of Laminar‐Turbulent Bifurcation Scenarios in Kelvin‐Helmholtz… http://dx.doi.org/10.5772/67918 51

Figure 11. Projection of limited cycles of period three into υz;υ<sup>y</sup> phase subspace for 2D perturbation problem for R ¼ 2417;Ri ¼ 1=6;Lx ¼ 6.

#### 4.3.1. 3D perturbations

4.2.3. 3D perturbations

values are set as follows:

 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7

u

Feigenbaum-Sharkovskiy cascade.

ρj

uxj


v

<sup>∂</sup>Ω<sup>0</sup> ¼ 1:1 þ

pj

1

<sup>∂</sup>Ω<sup>0</sup> <sup>¼</sup> <sup>2</sup>:<sup>0</sup> <sup>þ</sup> erf <sup>y</sup> � <sup>1</sup>=<sup>2</sup>

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

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

Figure 10. 2D invariant torus (R ¼ 760;Lx ¼ 8) and 3D invariant torus (R ¼ 950;Lx ¼ 8) first Poincaré sections.

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 ∂Ω<sup>0</sup> and initial

<sup>10</sup> erf <sup>y</sup> � <sup>1</sup>=<sup>2</sup>

<sup>∂</sup>Ω<sup>0</sup> ¼ 22:14 þ ρgy þ ðρmax � ρÞg:

The Richardson number (Eq. (3)) is set depending on the problem in range of 1=5 � 1=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

<sup>δ</sup> <sup>þ</sup> Cz cos <sup>ð</sup>6πz<sup>Þ</sup> 

<sup>δ</sup> <sup>þ</sup> Cz cos <sup>ð</sup>6πz<sup>Þ</sup>  ,

ð30Þ

,

and three-dimensional invariant tori for 100 < R < 950; see Figure 10.

50 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

"point\_4.dat" u 1:2

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.

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

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

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.

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.

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.

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.

## 5. Discussion and conclusion

#### 5.1. Discussion

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

Figure 15. Projection of the invariant three-dimensional torus and first Poincaré section for R ¼ 516 to the three-dimen-

Figure 14. Projection of the limited cycle R ¼ 10:5 and invariant torus (with Poincaré section) R ¼ 510 to three-dimen-

sional phase subspace.

sional phase subspace and second section for R ¼ 518.

Figure 13. Coupled problem quasi-stationary solution. Density sections in yz planes with x ¼ 3 and x ¼ 5.

52 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

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 \mathbb{C} \to \dots \to \mathbb{C} \\ \mathfrak{S} \to \dots \to \mathbb{C}11 \to \dots \to \mathbb{C} \tag{31}$$

KH instability with 3D perturbations:

$$P \to \mathbb{C} \to T\mathbb{2} \to T\mathbb{3} \to \mathbb{C} \text{haos.}$$

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

$$P \to \mathbb{C} \to \mathbb{C} \\ 2 \to \mathbb{C} \\ 4 \to \mathbb{C} \\ n \to \dots \to \mathbb{C} \\ 5 \to \mathbb{C} \\ 3 \to \text{Chos}. \tag{32}$$

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

P ! C ! nC ! … ! C ! T2 ! T3 ! … ! nT3 ! … ! T4 ! Chaos: ð33Þ

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:<sup>6</sup> � 109 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.

## Author details

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

54 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

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

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

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

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

P ! C ! T2 ! T3 ! Chaos:

P ! C ! … ! C5 ! … ! C11 ! … ! C ð31Þ

high dimensional attractor. That is clearly the case for coupled problem; see Figure 12.

that were traced by other authors [48–50].

RT instabilities in the setup with the coaxial flow.

is stable.

5.2. Conclusion

bifurcation scenarios.

KH instability with 2D perturbations:

KH instability with 3D perturbations:

Nikolay Mihaylovitch Evstigneev and Nikolai Alexandrovitch Magnitskii\*

\*Address all correspondence to: nikmagn@gmail.com

Federal Research Center "Informatics and Control", Institute for System Analysis of Russian Academy of Science, Russia

## References


[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

[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/

[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] 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),

[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] Manneville P. (2010) Rayleigh-Benard convection: thirty years of experimental, theoreti-

[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

[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] Getling A.V. (1998) Rayleigh-Benard convection: structures and dynamics. World Scien-

[14] Ashwin P., Podvigina O. (2003) Hopf bifurcation with cubic symmetry and instability of

[15] Pfister G., Schmidt H., Cliffet K.A., Mullin T. (1988) Bifurcation phenomena in Taylor-

[16] Mamun C.K., Tuckerman L.S. (1994) Asymmetry and Hopf bifurcation in spherical

[17] Lopez J.M., Marques F. (2005) Finite aspect ratio Taylor–Couette flow: Shil'nikov dynam-

[18] Ma Tian, Wang Shouhong (2006) Stability and bifurcation of the Taylor problem. Arch.

[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 Visu-

[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] Rani H.P., Sheu Tony W.H. (2006) Nonlinear dynamics in a backward-facing step flow.

cal, and modeling work. Springer Tracts Modern Phys. 207, 41–65

convection. Int. J. Bifurc. Chaos. 22(5), 1230018

ABC flow. Proc. R. Soc Lond. A. 459(2035), 1801–1827

Couette Flow. Phys. Fluids. 7(1). doi:10.1063/1.868730

Couette flow in a very short annulus. J. Fluid Mech. 191, 1–18

tific, Singapore. ISBN 978-981-02-2657-2

ics of 2-tori. Phys. D. 211, 168–191

Ration. Mech. Anal. 181, 149–176

alization, Minsk, Belarus

Phys. Fluids. 18, 084101

[6] Sivashinsky G. (1985) Weak turbulence in periodic flows. Phys. D. 17(2), 243–255

56 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

BF01017632

1292–1305


[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)

[40] Fjørtoft T. (1950) Application of integral theorems in deriving criteria of instability of

[41] Mack C. (2009) Global stability of compressible flow about a swept parabolic body. PhD

[42] Henze O., Lemke M., Sesterhenn J. (2015) A parallel and matrix free framework for global

[43] Bernstein I.B., Book D.L. (1983) Effect of compressibility on the Rayleigh-Taylor instabil-

[44] Lezzi A.M., Prosperetti A. (1989) Rayleigh-Taylor instability for adiabatically stratified

[45] Turner L. (2002) Rayleigh-Taylor instabilities and gravity waves in compressible fluids.

[46] Livescu D. (2004) Compressibility effects on the Rayleigh-Taylor instability growth

[47] Guo Yan, Tice Ian (2009) Compressible, inviscid Rayleigh-Taylor instability. Indiana Univ.

[48] Poehlmann A., Richter R., Rehberg I. (2013) Unravelling the Rayleigh-Taylor instability

[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,

[50] Hunter J.K., Thoo J.B. (2011) On the weakly nonlinear Kelvin–Helmholtz instability of

[51] Magnitskii N.A., Sidorov S.V. (2006) New methods for chaotic dynamics. Singapore:

[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

[53] Toro E.F. (1999) Riemann solvers and numerical methods for fluid dynamics. Springer-

[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

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

tangential discontinuities in MHD. J. Hyperbol. Differ. Eqn. 8(4), 691–726

applications. Intech, pp. 133–174. http://dx.doi.org/10.5772/48811

Los Alamos National Laboratory Report LA-UR-02-6439, Los Alamos, USA.

laminar flow and for the baroclinic circular vortex. Geophys. Publ. 17, 1–52

stability analysis of compressible flows. arXiv:1502.03701

58 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

between immiscible fluids. Phys. Fluids. 16(1), 118–127

by stabilization. J. Fluid Mech. 732, R3. doi:10.1017/jfm.2013.424

Math. J. 60(2). doi:10.1512/iumj.2011.60.4193

thesis, Ecole Polytechnique

ity. Phys. Fluids 26, 453

94–121

World Scientific

Verlag, Berlin-Heidelberg.

(5), 721–753. (in Russian)

(in Russian)

fluids. Phys. Fluids A. 1, 1784


## **Interface Instability and Turbulent Mixing**

Jingsong Bai and Tao Wang

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/67750

## Abstract

Richtmyer-Meshkov instability and turbulent mixing are fundamental problems of multi-materials interface dynamics, which mainly focuses on the growth of perturbation on the interface and mixing of different materials. It is very important in many applications such as inertial confinement fusion, high-speed combustion, supernova, etc. In this chapter, we will gain advances in understanding this problem by numerical investigations, including the numerical method and program we used, the verification and validation of numerical method and program, the growth laws and mechanics of turbulent mixing, the effects of initial conditions, the dynamic behavior, and some new phenomenon for Richtmyer-Meshkov instability and turbulent mixing.

Keywords: Richtmyer-Meshkov instability, turbulent mixing, interface dynamics, dynamic behavior

## 1. Introduction

When a shock wave accelerates a perturbed interface separating two different fluids, the Richtmyer-Meshkov (RM) instability will occur. This phenomenon was theoretically analyzed by Richtmyer [1] for the first time in 1960, which was confirmed in experiment by Meshkov [2] in 1969. The main mechanism of the Richtmyer-Meshkov instability (RMI) is the baroclinic vorticity deposition at the interface due to the misalignment of the pressure gradient across the shock and the local density gradient at the interface (∇ρ � ∇p 6¼ 0). At the beginning, the perturbations grow linearly. When entering the nonlinear stage, the perturbations develop into complex structures formed by bubbles (the portion of the light fluid penetrating into the heavy fluid) and spikes shaped with "mushrooms" (the portion of the heavy fluid penetrating into the light fluid), see Figure 1. Afterward, the mushroom structures are eroded and broke up which results in the turbulent mixing eventually. When

Figure 1. Development of the Richtmyer-Meshkov perturbed interface between air and SF6 gases [3].

the incident shock wave impacts the initial perturbed interface, it bifurcates into a transmitted shock wave and a reflected wave. Depending on the material properties of the fluid on both sides of the interface, the reflected wave can be either a shock or a rarefaction wave. The criterion is that when the incident shock wave interacts with the interface from light fluid to heavy fluid, the reflected wave is a shock wave; otherwise, it is a rarefaction wave. If the transmitted shock wave meets a rigid wall and is reflected back to collide with the interface once again, this process is called reshock which can advance the transition to turbulent mixing.

The Richtmyer-Meshkov instability and induced turbulent mixing are very important in a variety of man-made applications and natural phenomena such as inertial confinement fusion (ICF) [4], deflagration-to-detonation transition (DDT) [5], high-speed combustion [6], and astrophysics (i.e., supernova explosions) [7]. For ICF, the ablative shell that encapsulates the deuterium-tritium fuel becomes RM unstable as it is accelerated inward by the ablation of its outer surface by laser or secondary X-ray radiation. The degree of compression achievable in laser fusion experiments is ultimately limited by Richtmyer-Meshkov and Rayleigh-Taylor instabilities. Thus, these instabilities represent the most significant barriers to attaining positive-net-yield fusion reactions in laser fusion facilities. For the supersonic combustion, the Richtmyer-Meshkov instability caused by the interaction of a shock wave with a flame front can greatly promote the mixing of fuel and oxidant and enhance the burning rate. For the supernova explosions, the Richtmyer-Meshkov instability was believed to occur when the outward propagating shock wave generated by the collapsing core of a dying star passes over the helium-hydrogen interface. Observations of the optical output of the supernova 1987A suggest that the outer regions of the supernova were much more uniformly mixed than expected, and indicating significant Richtmyer-Meshkov mixing had occurred [8]. The Richtmyer-Meshkov instability and turbulent mixing are so important and have gained significant attention. However, the turbulent mixing is a complicated three-dimensional unstable flow, which spans a wide range of time-space scales.

## 2. Numerical method and program

By applying the piecewise parabolic method (PPM), volume of fluid (VOF), parallel technique, and so on, we have developed a series of Euler compressible multi-fluid dynamic programs with three orders precision, such as MFPPP [9], MVPPM [10], and multi-viscous-flow and turbulence (MVFT) [11, 12]. MFPPM program is not considering the fluid viscosity, which only solves the Euler equations. MVPPM program is considering the molecular viscosity but it is not changed with temperature. Multi-viscous-flow and turbulence (MVFT) is a large-eddy simulation (LES) program that has four choices of subgrid-scale (SGS) stress models including the Smagorinsky model [13], Vreman model (VM) [14], dynamic Smagorinsky model (DSM) [15], and stretched-vortex model (SVM) [16].

#### 2.1. Governing equations

the incident shock wave impacts the initial perturbed interface, it bifurcates into a transmitted shock wave and a reflected wave. Depending on the material properties of the fluid on both sides of the interface, the reflected wave can be either a shock or a rarefaction wave. The criterion is that when the incident shock wave interacts with the interface from light fluid to heavy fluid, the reflected wave is a shock wave; otherwise, it is a rarefaction wave. If the transmitted shock wave meets a rigid wall and is reflected back to collide with the interface once again, this process is called reshock which can advance the transition to

Figure 1. Development of the Richtmyer-Meshkov perturbed interface between air and SF6 gases [3].

62 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

The Richtmyer-Meshkov instability and induced turbulent mixing are very important in a variety of man-made applications and natural phenomena such as inertial confinement fusion (ICF) [4], deflagration-to-detonation transition (DDT) [5], high-speed combustion [6], and astrophysics (i.e., supernova explosions) [7]. For ICF, the ablative shell that encapsulates the deuterium-tritium fuel becomes RM unstable as it is accelerated inward by the ablation of its outer surface by laser or secondary X-ray radiation. The degree of compression achievable in laser fusion experiments is ultimately limited by Richtmyer-Meshkov and Rayleigh-Taylor instabilities. Thus, these instabilities represent the most significant barriers to attaining positive-net-yield fusion reactions in laser fusion facilities. For the supersonic combustion, the Richtmyer-Meshkov instability caused by the interaction of a shock wave with a flame front can greatly promote the mixing of fuel and oxidant and enhance the burning rate. For the supernova explosions, the Richtmyer-Meshkov instability was believed to occur when the outward propagating shock wave generated by the collapsing core of a dying star passes over the helium-hydrogen interface. Observations of the optical output of the supernova 1987A suggest that the outer regions of the supernova were much more uniformly mixed than expected, and indicating significant Richtmyer-Meshkov mixing had occurred [8]. The Richtmyer-Meshkov instability and turbulent mixing are so important and have gained significant attention. However, the turbulent mixing is a complicated three-dimensional unstable flow, which spans a

turbulent mixing.

wide range of time-space scales.

The governing equations of large-eddy simulation are the Favre-filtered compressible multiviscous-flow Navier-Stokes (NS) equations and are written as in tensor form:

$$\begin{aligned} \frac{\partial \overline{\rho}}{\partial t} + \frac{\partial \overline{\rho} \widetilde{u}\_j}{\partial x\_j} &= 0 \\ \frac{\partial \overline{\rho} \widetilde{u} \widetilde{u}}{\partial t} + \frac{\partial \left( \overline{\rho} \widetilde{u}\_i \widetilde{u}\_j + \overline{p} \delta\_{ij} \right)}{\partial x\_j} &= \frac{\partial \sigma\_{\widetilde{y}}}{\partial x\_j} - \frac{\partial \tau\_{\widetilde{y}}}{\partial x\_j} \\ \frac{\partial \overline{\rho} \overline{E}}{\partial t} + \frac{\partial \left( \overline{\rho} \widetilde{u}\_i \overline{E} + \overline{p} \widetilde{u}\_j \right)}{\partial x\_j} &= \frac{\partial \left( \sigma\_{\widetilde{y}} - \tau\_{\widetilde{u}j} \right) \widetilde{u}\_i}{\partial x\_j} - \frac{\partial \left( q\_j^{\prime} + \mathbf{Q}\_j^{\prime} \right)}{\partial x\_j} \\ \frac{\partial \overline{Y}^{\prime}(s)}{\partial t} + \widetilde{u}\_j \frac{\partial \overline{Y}^{\prime}(s)}{\partial x\_j} &= \frac{\partial}{\partial x\_j} \left( \widetilde{D} \frac{\bar{\mathbf{Y}}^{\prime}(s)}{\partial x\_j} \right) - \frac{\partial \mathbf{Q}\_j^{\prime}}{\partial x\_j} \qquad s = 1, 2, \cdots, N - 1 \end{aligned} \tag{1}$$

here i and j represent the directions of x, y, and z, and the same subscripts denote the tensor summation; ρ, u~kðk ¼ i, jÞ, p, and E are the resolved-scale density, velocity, pressure, and total energy per unit mass; <sup>N</sup> is the number of species; <sup>Y</sup><sup>~</sup> <sup>ð</sup>s<sup>Þ</sup> is the volume fraction of the <sup>s</sup>th fluid and satisfies <sup>X</sup><sup>N</sup> 1 <sup>Y</sup><sup>~</sup> <sup>ð</sup>s<sup>Þ</sup> <sup>¼</sup> 1; <sup>D</sup><sup>~</sup> is the diffusion coefficient set to <sup>D</sup><sup>~</sup> <sup>¼</sup> <sup>ν</sup>=Sc, where <sup>ν</sup> is the kinematic viscosity of the fluid, Sc is the Schmidt number. σij is the deviatoric Newtonian stress tensor, i.e.

$$
\sigma\_{i\bar{j}} = \mu\_l \left( \frac{\partial \tilde{u}\_i}{\partial \mathbf{x}\_{\bar{j}}} + \frac{\partial \tilde{u}\_{\bar{j}}}{\partial \mathbf{x}\_i} - \frac{2}{3} \delta\_{\bar{i}\bar{j}} \left( \frac{\partial \tilde{u}\_k}{\partial \mathbf{x}\_k} \right) \right) \tag{2}
$$

where μ<sup>l</sup> is the dynamic viscosity. ql <sup>j</sup> ¼ �λl∂T<sup>~</sup> <sup>=</sup>∂xj is the resolved heat transport flux per unit time and space, λ<sup>l</sup> ¼ μlcp/Pr<sup>l</sup> is the resolved heat conduction coefficient, cp is the constant pressure specific heat, Pr<sup>l</sup> is the Prandtl number, and T~ is the fluid temperature. The equation of state (EOS) has choices of the ideal gas state form, stiffen gas state form, reduction form of Gruneisen EOS for condensed matter.

For a multi-material mixture case, the average quantities and physical properties of a mixed phase are supposed as a volume-weighted sum over the set of components ρ � �ð Þ<sup>s</sup> , <sup>ρ</sup>u~<sup>i</sup> � �ð Þ<sup>s</sup> , and <sup>ρ</sup>E<sup>~</sup> � �ð Þ<sup>s</sup> used at the initial time, and so on for each, respectively,

$$\overline{\rho} = \sum\_{s=1}^{N} \check{Y}^{(s)} \overline{\rho}^{(s)} \tag{3}$$

$$
\overline{\rho}\check{u}\_i = \sum\_{s=1}^N \check{Y}^{(s)}\left(\overline{\rho}\check{u}\_i\right)^{(s)}\tag{4}
$$

$$\overline{\rho}\tilde{E} = \sum\_{s=1}^{N} \tilde{Y}^{(s)} \left[ \left( \overline{\rho}\overline{e} \right)^{(s)} + \frac{1}{2} \left( \overline{\rho}\tilde{u}\_i^2 \right)^{(s)} \right] \tag{5}$$

$$\frac{1}{\gamma - 1} = \sum\_{s=1}^{N} \frac{\tilde{Y}^{(s)}}{\gamma^{(s)} - 1} \tag{6}$$

$$\mu = \sum\_{s=1}^{N} \mu^{(s)} \,\tilde{Y}^{(s)} \tag{7}$$

$$\overline{D} = \sum\_{s=1}^{N} \overline{D}^{(s)} \tilde{Y}^{(s)} \tag{8}$$

For large-eddy simulation, in the filtering operation, the unresolved-scale motions identified as the "subgrid-scale" are filtered, but the effects of subgrid-scale motions on the resolved-scale motions are retained in the governing equations in the form of SGS turbulence transport fluxes, which must be modeled to complete the closure of LES equations. The SGS stress tensor, the heat, and scalar transport flux are defined as

$$
\pi\_{i\circ} = \overline{\rho} \left( \widehat{u\_i u\_{\circ}} - \tilde{u}\_i \tilde{u}\_{\circ} \right) \tag{9}
$$

$$\mathbf{Q}\_{j}^{T} = \overline{\rho} \left( \overline{c\_{p}} \overline{T} \overline{u}\_{j} - \tilde{c}\_{p} \tilde{T} \tilde{u}\_{j} \right) \tag{10}$$

$$\mathbf{Q}\_{\dagger}^{\chi} = \left( \overleftarrow{\mathbf{\tilde{Y}u\_{\dagger}}} - \mathbf{\tilde{Y}u\_{\dagger}} \right) \tag{11}$$

#### 2.2. Subgrid-scale stress models for LES

#### 2.2.1. Smagorinsky model

The SGS turbulence behavior is assumed to be analogous to the molecular dissipative mechanism, so that the molecular viscosity and diffusion models can be used to simulate the SGS fluxes, and the SGS stress tensor, the heat, and the scalar transport flux are calculated as follows [13]

#### Interface Instability and Turbulent Mixing http://dx.doi.org/10.5772/67750 65

$$\pi\_{\vec{\eta}} = -\mu\_t \left( \frac{\partial \tilde{u}\_i}{\partial \mathbf{x}\_j} + \frac{\partial \tilde{u}\_j}{\partial \mathbf{x}\_i} - \frac{2}{3} \delta\_{\vec{\eta}} \left( \frac{\partial \tilde{u}\_k}{\partial \mathbf{x}\_k} \right) \right) \tag{12}$$

$$Q\_j^T = -\lambda\_t \frac{\partial \tilde{T}}{\partial \mathbf{x}\_j} = -\frac{\mu\_t c\_\mathcal{P}}{\mathbf{Pr}\_t} \frac{\partial \tilde{T}}{\partial \mathbf{x}\_j} \tag{13}$$

$$Q\_j^Y = -D\_t \frac{\partial \tilde{Y}^{(\mathbf{s})}}{\partial \mathbf{x}\_j} = -\frac{\mu\_t}{\overline{\rho} \mathbf{S} c\_t} \frac{\partial \tilde{Y}^{(\mathbf{s})}}{\partial \mathbf{x}\_j} \tag{14}$$

where the SGS turbulent viscosity μ<sup>t</sup> is calculated by the Smagorinsky eddy viscosity model,

$$
\mu\_t = 2C\rho\Delta^2|\overline{S}|\tag{15}
$$

where the dimensionless coefficient <sup>C</sup> <sup>¼</sup> <sup>C</sup><sup>2</sup> <sup>S</sup>, for the isotropic turbulence, the model constant is CS <sup>¼</sup> <sup>0</sup>:17, <sup>Δ</sup> is the grid-filter width, and <sup>j</sup> <sup>S</sup> j ¼ <sup>2</sup>SijSij � �<sup>1</sup>=<sup>2</sup> is the magnitude of the resolved strain-rate tensor,

$$\overline{S}\_{i\dot{\jmath}} = \frac{1}{2} \left( \frac{\partial \tilde{u}\_i}{\partial \mathbf{x}\_{\dot{\jmath}}} + \frac{\partial \tilde{u}\_{\dot{\jmath}}}{\partial \mathbf{x}\_i} \right) \tag{16}$$

#### 2.2.2. Vreman model

The Vreman SGS model is also an eddy viscosity model and is as follows [14]:

$$
\mu\_t = c\overline{\rho} \sqrt{\frac{B\_{\overline{\rho}}}{\alpha\_{\overline{\imath}\overline{\alpha}\_{\overline{\imath}}}}} \tag{17}
$$

with

For a multi-material mixture case, the average quantities and physical properties of a mixed

� �ð Þ<sup>s</sup> , <sup>ρ</sup>u~<sup>i</sup>

<sup>ρ</sup>ð Þ<sup>s</sup> <sup>ð</sup>3<sup>Þ</sup>

� �ð Þ<sup>s</sup> <sup>ð</sup>4<sup>Þ</sup>

<sup>μ</sup>ð Þ<sup>s</sup> <sup>Y</sup><sup>~</sup> ð Þ<sup>s</sup> <sup>ð</sup>7<sup>Þ</sup>

<sup>D</sup>ð Þ<sup>s</sup> <sup>Y</sup><sup>~</sup> ð Þ<sup>s</sup> <sup>ð</sup>8<sup>Þ</sup>

� � <sup>ð</sup>9<sup>Þ</sup>

� �ð Þ<sup>s</sup> ,

ð5Þ

ð6Þ

ð10Þ

ð11Þ

phase are supposed as a volume-weighted sum over the set of components ρ

used at the initial time, and so on for each, respectively,

<sup>ρ</sup> <sup>¼</sup> <sup>X</sup> N

<sup>ρ</sup>u~<sup>i</sup> <sup>¼</sup> <sup>X</sup> N

<sup>Y</sup><sup>~</sup> ð Þ<sup>s</sup>

1 <sup>γ</sup> � <sup>1</sup> <sup>¼</sup> <sup>X</sup>

> <sup>μ</sup> <sup>¼</sup> <sup>X</sup> N

<sup>D</sup> <sup>¼</sup> <sup>X</sup> N

s¼1

s¼1

For large-eddy simulation, in the filtering operation, the unresolved-scale motions identified as the "subgrid-scale" are filtered, but the effects of subgrid-scale motions on the resolved-scale motions are retained in the governing equations in the form of SGS turbulence transport fluxes, which must be modeled to complete the closure of LES equations. The SGS stress

<sup>τ</sup>ij <sup>¼</sup> <sup>ρ</sup> <sup>u</sup>giuj � <sup>u</sup>~iu~<sup>j</sup>

<sup>g</sup>Tuj � <sup>~</sup>cpT~u~<sup>j</sup> � �

<sup>g</sup><sup>j</sup> � <sup>Y</sup><sup>~</sup> <sup>u</sup>~<sup>j</sup> � �

The SGS turbulence behavior is assumed to be analogous to the molecular dissipative mechanism, so that the molecular viscosity and diffusion models can be used to simulate the SGS fluxes, and the SGS stress tensor, the heat, and the scalar transport flux are calculated as

QT <sup>j</sup> ¼ ρ cp

> Q<sup>Y</sup> <sup>j</sup> ¼ Yu

<sup>ρ</sup>E<sup>~</sup> <sup>¼</sup> <sup>X</sup> N

64 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

tensor, the heat, and scalar transport flux are defined as

2.2. Subgrid-scale stress models for LES

2.2.1. Smagorinsky model

follows [13]

s¼1

s¼1

s¼1

<sup>Y</sup><sup>~</sup> ð Þ<sup>s</sup>

<sup>Y</sup><sup>~</sup> ð Þ<sup>s</sup>

<sup>ρ</sup><sup>e</sup> � �ð Þ<sup>s</sup> <sup>þ</sup>

N

s¼1

ρu~<sup>i</sup>

1 <sup>2</sup> <sup>ρ</sup>u~<sup>2</sup> i

� �ð Þ<sup>s</sup> " #

<sup>Y</sup><sup>~</sup> ð Þ<sup>s</sup> <sup>γ</sup>ð Þ<sup>s</sup> � <sup>1</sup>

and <sup>ρ</sup>E<sup>~</sup> � �ð Þ<sup>s</sup>

$$\begin{aligned} \alpha\_{\tilde{i}\tilde{j}} &= \partial\_i \tilde{u}\_{\tilde{j}} = \frac{\partial \tilde{u}\_{\tilde{j}}}{\partial x\_i} \\ \beta\_{\tilde{i}\tilde{j}} &= \Delta\_m^2 \alpha\_{mi} \alpha\_{mj} \\ B\_{\tilde{\rho}} &= \beta\_{11} \beta\_{22} - \beta\_{12}^2 + \beta\_{11} \beta\_{33} - \beta\_{13}^2 + \beta\_{22} \beta\_{33} - \beta\_{23}^2 \end{aligned} \tag{18}$$

The model constant c is related to the Smagorinsky SGS model constant CS by c ≈ 2:5C<sup>2</sup> <sup>S</sup>. The symbol α represents the 3 � 3 matrix of derivatives of the filtered velocity u~. If αijαij equals zero, μ<sup>t</sup> is consistently defined to be zero. In fact, B<sup>β</sup> is an invariant of the matrix β, while αijαij is an invariant of α<sup>T</sup> <sup>α</sup>. If the filter width is the same in each direction, then <sup>Δ</sup><sup>i</sup> <sup>¼</sup> <sup>Δ</sup> and <sup>β</sup> <sup>¼</sup> <sup>Δ</sup><sup>2</sup> αT α.

#### 2.2.3. Dynamic Smagorinsky model

In general, the turbulent kinetic energy is transferred from large scales to small scales in turbulent flows, which is called forward scatter of energy, and then it is dissipated by the viscous action. However, the reversed energy flow from small scales to large scales (the process called backscatter) may also occur. Although the backscatter is just a small range of local phenomenon, it has been shown to be of importance, especially in the transition regime. The physical mechanism of most of the widely used SGS models is the forward scatter of energy; in other words, the SGS models are absolutely dissipative, such as the eddy viscosity models (Smagorinsky and Vreman SGS models). In order to account for the backscatter of energy, several modifications to the eddy viscosity models have been proposed. An improvement is to calculate the eddy viscosity coefficient dynamically which is a function of space and time and can be negative in some regions [15].

The dynamic Smagorinsky model for compressible turbulence is as follows, the model coefficient is

$$\mathbb{C}\_D = \frac{\langle L\_{\hat{\boldsymbol{\eta}}} \boldsymbol{M}\_{\hat{\boldsymbol{\eta}}} \rangle}{\langle \boldsymbol{M}\_{\hat{\boldsymbol{\eta}}} \boldsymbol{M}\_{\hat{\boldsymbol{\eta}}} \rangle} \tag{19}$$

where Lij is the Leonard stress, 〈 〉 indicates the statistics averaging.

$$L\_{ij} = \widehat{\overline{\rho}\widehat{u\_i}\tilde{u\_j}} - \frac{1}{\overline{\overline{\rho}}} \widehat{\left(\widehat{\overline{\rho}\tilde{u\_i}}\widehat{\overline{\rho}\tilde{u\_j}}\right)}\tag{20}$$

$$\rho$$

$$M\_{\vec{\eta}} = -\left[2\hat{\vec{\Delta}}^2 \hat{\vec{\rho}} \left| \hat{\vec{S}} \right| \left( \hat{\vec{S}}\_{\vec{\eta}} - \frac{\delta\_{\vec{\eta}}}{3} \hat{\vec{S}}\_{kk} \right) - 2\overline{\Delta}^2 \overline{\rho} \left| \overline{\vec{S}} \left( \overline{\vec{S}\_{\vec{\eta}} - \frac{\delta\_{\vec{\eta}}}{3} \tilde{\mathbf{S}}\_{kk} \right) \right. \tag{21}$$

the overbars of "¯" and "^" denote the grid filter and test filter, respectively.

#### 2.2.4. Stretched-vortex model

The stretched-vortex model is based on an explicit structural modeling of small-scale dynamics. It can simulate the multiscale compressible turbulence and allows the anisotropy of the flow to be extended to the dissipation scale. In the stretched-vortex model, the flow within a computational grid cell is assumed to result from an ensemble of straight, nearly axisymmetric vortices aligned with the local resolved scale strain or vorticity. The resulting SGS flux terms are [16]

$$
\pi\_{i\dot{j}} = \overline{\rho} \tilde{k} \left( \delta\_{i\dot{j}} - e\_i^\nu e\_{\dot{j}}^\nu \right) \tag{22}
$$

$$Q\_i^T = -\overline{\rho} \frac{\Delta\_c}{2} \tilde{k}^j \mathbb{A} \Big(\delta\_{\vec{\eta}} - e\_i^\nu e\_{\vec{\eta}}^\nu \Big) \frac{\partial \left(\tilde{c}\_p \tilde{T}\right)}{\partial \mathbf{x}\_{\vec{\eta}}} \tag{23}$$

$$Q\_i^Y = -\frac{\Delta\_c}{2} \tilde{k} \mathbb{W} \Big( \delta\_{i\bar{\jmath}} - e\_i^\nu e\_{\bar{\jmath}}^\nu \Big) \frac{\partial \tilde{Y}}{\partial \mathbf{x}\_{\bar{\jmath}}} \tag{24}$$

where <sup>~</sup><sup>k</sup> <sup>¼</sup> ð∞ kc E kð Þdk is the subgrid turbulent kinetic energy, e <sup>ν</sup> is the unit vector aligned with the subgrid vortex axis, kc <sup>¼</sup> <sup>π</sup>=Δ<sup>c</sup> is the cutoff wave number, E kð Þ¼ <sup>K</sup>0ε<sup>2</sup>=<sup>3</sup>k�<sup>5</sup>=<sup>3</sup> exp �2<sup>k</sup> 2 <sup>ν</sup>=ð Þ <sup>3</sup>j~a<sup>j</sup> � � represents the energy spectrum of subgrid vortices, and K<sup>0</sup> is the Kolmogorov prefactor, ε is the local cell-averaged dissipation, <sup>~</sup><sup>a</sup> <sup>¼</sup> <sup>S</sup><sup>~</sup> ijeν i eν <sup>j</sup> is the axial strain along the subgrid vortex axis, <sup>S</sup><sup>~</sup> ij denotes the resolved rate-of-strain tensor.

#### 2.3. Numerical algorithm

For the convenience, in the following sections, the overbar and overtilde of variables are omitted. Operator splitting technique is used to solve the physical problems, described by Eq. (1), into three sub-processes, i.e., the computation of inviscid flux, viscous flux, and heat flux. Eq. (1) can be split into two equations as follows:

$$\begin{aligned} \frac{\partial \rho}{\partial t} + \frac{\partial \rho u\_j}{\partial x\_j} &= 0 \\ \frac{\partial \rho u\_i}{\partial t} + \frac{\partial (\rho u\_i u\_j + p \delta\_{ij})}{\partial x\_j} &= 0 \\ \frac{\partial \rho E}{\partial t} + \frac{\partial (\rho u\_i E + p u\_j)}{\partial x\_j} &= 0 \\ \frac{\partial \rho Y^{(s)}}{\partial t} + u\_j \frac{\partial Y^{(s)}}{\partial x\_j} &= 0 \quad s = 1, 2, \dots, N - 1 \end{aligned} \tag{25}$$

and

local phenomenon, it has been shown to be of importance, especially in the transition regime. The physical mechanism of most of the widely used SGS models is the forward scatter of energy; in other words, the SGS models are absolutely dissipative, such as the eddy viscosity models (Smagorinsky and Vreman SGS models). In order to account for the backscatter of energy, several modifications to the eddy viscosity models have been proposed. An improvement is to calculate the eddy viscosity coefficient dynamically which is a function of space and

The dynamic Smagorinsky model for compressible turbulence is as follows, the model coeffi-

CD <sup>¼</sup> 〈LijMij〉

� <sup>2</sup>Δ<sup>2</sup> ρ

� � <sup>2</sup>

The stretched-vortex model is based on an explicit structural modeling of small-scale dynamics. It can simulate the multiscale compressible turbulence and allows the anisotropy of the flow to be extended to the dissipation scale. In the stretched-vortex model, the flow within a computational grid cell is assumed to result from an ensemble of straight, nearly axisymmetric vortices aligned

S � � � � �

� � � � �

where Lij is the Leonard stress, 〈 〉 indicates the statistics averaging.

66 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

� � � � �

Lij <sup>¼</sup> b ρu~<sup>i</sup> <sup>u</sup>~<sup>j</sup> � <sup>1</sup> <sup>ρ</sup>^ b ρu~<sup>i</sup> b ρu~<sup>j</sup> � �

^ Sij � <sup>δ</sup>ij 3 ^ Skk � �

the overbars of "¯" and "^" denote the grid filter and test filter, respectively.

with the local resolved scale strain or vorticity. The resulting SGS flux terms are [16]

Δc 2

subgrid vortex axis, kc <sup>¼</sup> <sup>π</sup>=Δ<sup>c</sup> is the cutoff wave number, E kð Þ¼ <sup>K</sup>0ε<sup>2</sup>=<sup>3</sup>k�<sup>5</sup>=<sup>3</sup> exp �2<sup>k</sup>

<sup>τ</sup>ij <sup>¼</sup> <sup>ρ</sup>~<sup>k</sup> <sup>δ</sup>ij � <sup>e</sup>

<sup>~</sup>k<sup>½</sup> <sup>δ</sup>ij � <sup>e</sup>

<sup>~</sup>k<sup>½</sup> <sup>δ</sup>ij � <sup>e</sup>

represents the energy spectrum of subgrid vortices, and K<sup>0</sup> is the Kolmogorov prefactor, ε is the

ν i e ν j � �

ν i e ν j � � <sup>∂</sup> <sup>~</sup>cpT<sup>~</sup> � �

ν i e ν j � � ∂Y~

∂xj

∂xj

〈MijMij〉 <sup>ð</sup>19<sup>Þ</sup>

3

<sup>ν</sup> is the unit vector aligned with the

2 <sup>ν</sup>=ð Þ <sup>3</sup>j~a<sup>j</sup> � �

5 ð21Þ

b

Sij � <sup>δ</sup>ij <sup>3</sup> Skk ð20Þ

ð22Þ

ð23Þ

ð24Þ

time and can be negative in some regions [15].

Mij ¼ � <sup>2</sup>^

2.2.4. Stretched-vortex model

where <sup>~</sup><sup>k</sup> <sup>¼</sup>

ð∞ kc Δ 2 ρ^ ^ S � � � � �

> QT <sup>i</sup> ¼ �ρ

> > Q<sup>Y</sup> <sup>i</sup> ¼ � <sup>Δ</sup><sup>c</sup> 2

E kð Þdk is the subgrid turbulent kinetic energy, e

4

cient is

$$\begin{aligned} \frac{\partial \rho}{\partial t} &= 0\\ \frac{\partial \rho u\_i}{\partial t} &= \frac{\partial \sigma\_{ij}}{\partial x\_j} - \frac{\partial \tau\_{ij}}{\partial x\_j} \\ \frac{\partial \rho E}{\partial t} &= \frac{\partial (\sigma\_{ij} - \tau\_{ij})u\_i}{\partial x\_j} - \frac{\partial \left(q\_j^l + Q\_j^T\right)}{\partial x\_j} \\ \frac{\partial \chi^{(s)}}{\partial t} &= \frac{\partial}{\partial x\_j} \left(D \frac{\partial \chi^{(s)}}{\partial x\_j}\right) - \frac{\partial Q\_j^T}{\partial x\_j} \quad s = 1, 2, \dots, N - 1 \end{aligned} \tag{26}$$

For the inviscid flux, the three-dimensional problem can be simplified into one-dimensional (1D) problem in three directions of x, y, and z by the dimension splitting technique,

$$\begin{aligned} \frac{\partial \rho}{\partial t} + \frac{\partial (\rho u)}{\partial x} &= 0 \\ \frac{\partial (\rho u)}{\partial t} + \frac{\partial (\rho u^2 + p)}{\partial x} &= 0 \\ \frac{\partial (\rho E)}{\partial t} + \frac{\partial (\rho u E + p u)}{\partial x} &= 0 \\ \frac{\partial \rho Y(s)}{\partial t} + \mu \frac{\partial Y(s)}{\partial x} &= 0 \qquad s = 1, 2, \dots, N - 1 \end{aligned} \tag{27}$$

And then the one-dimensional Eq. (27) in each direction is resolved by the two-step Lagrange-Remapping algorithm. Also, one time step is divided into four substeps: (i) the piecewise parabolic interpolating of physical quantities, (ii) solving the Riemann problems approximately, (iii) marching of the Lagrange equations, and (iv) remapping the physical quantities back to the stationary Euler meshes. The orders of accuracy of the spatial and temporal schemes are the third and second, respectively, for smooth flows.

#### 2.3.1. Lagrange step of finite volume method

The Lagrange equations in 1D for multi-fluid can be written as

$$\begin{cases} \frac{\partial \tau}{\partial t} - \frac{\partial (r^a u)}{\partial m} = 0 \\ \frac{\partial u}{\partial t} + r^a \frac{\partial p}{\partial m} = 0 \\ \frac{\partial E}{\partial t} + \frac{\partial (r^a u p)}{\partial m} = 0 \\ \frac{\partial Y^{(s)}}{\partial t} = 0 \qquad s = 1, 2, \dots, N - 1 \end{cases} \tag{28}$$

where τ is the specific volume, α is 0, 1, and 2 corresponding to the plane, axial symmetry, and spherical symmetry problems, r is the Lagrange spatial coordinates, m is the Lagrange mass coordinates m ¼ ðr r0 ρr <sup>α</sup>dr, <sup>u</sup> <sup>¼</sup> <sup>∂</sup><sup>r</sup> ∂t , mj�1/2 and mjþ1/2 are the mass coordinates of both sides of the <sup>j</sup>th grid, <sup>Δ</sup><sup>m</sup> <sup>¼</sup> mjþ1/2 � mj�1/2. The mass average of physical quantities <sup>τ</sup>, u, E, <sup>Y</sup>ð Þ<sup>s</sup> � �<sup>n</sup> j for the jth grid can be defined as

$$
\begin{pmatrix} \tau \\ u \\ E \\ Y(s) \end{pmatrix}\_j^n = \frac{1}{\Delta m\_j} \int\_{m\_{j-1/2}}^{m\_{j+1/2}} \begin{pmatrix} \tau(m\_\prime \ t\_n) \\ u(m\_\prime \ t\_n) \\ E(m\_\prime \ t\_n) \\ Y^{(s)}(m\_\prime \ t\_n) \end{pmatrix} dm \tag{29}
$$

Because of the calculations are referred to as scalars, here the averaged physical quantities in a cell are written as a uniformity Q<sup>n</sup> <sup>j</sup> at time t. The piecewise parabolic function Q<sup>n</sup> <sup>j</sup> ð Þ m in cells are constructed to compute the time average of physical quantity Q on both sides of the grid edge mjþ1/2,Q<sup>~</sup> <sup>j</sup>þ1=2,L and <sup>Q</sup><sup>~</sup> <sup>j</sup>þ1=2,R. Then the Riemann problem at the grid edge mjþ1/2 is solved by using double shock wave approximation. After Lagrange marching step, we can obtain distribution of the physical quantities τ<sup>n</sup>þ<sup>1</sup> j n o, unþ<sup>1</sup> j n o, Enþ<sup>1</sup> j n o, and the position of new grids r nþ1 j�1=2 n o at time tnþ1, the pressure pnþ<sup>1</sup> j n o can be computed by equation of state, the <sup>Y</sup>ð Þ<sup>s</sup> � �<sup>n</sup>þ<sup>1</sup> j � � will not change. The marching formula is as follows

$$\begin{cases} r\_{j+1/2}^{n+1} = r\_{j+1/2}^{n} + \tilde{u}\_{j+1/2} \Delta t \\\\ \tau\_{j}^{n+1} = \frac{\left(r\_{j+1/2}^{n+1}\right)^{\alpha+1} - \left(r\_{j-1/2}^{n+1}\right)^{\alpha+1}}{(\alpha+1)\Delta m\_{j}} \\\\ u\_{j}^{n+1} = u\_{j}^{n} - \frac{\Lambda}{2} \frac{\Delta t}{\Delta m\_{j}} \left(\sigma\_{j+1/2} + \sigma\_{j-1/2}\right) \left(\tilde{p}\_{j+1/2} - \tilde{p}\_{j-1/2}\right) \\\\ E\_{j}^{n+1} = E\_{j}^{n} - \frac{\Lambda t}{\Delta m\_{j}} \left(\sigma\_{j+1/2}\tilde{u}\_{j+1/2}\tilde{p}\_{j+1/2} - \sigma\_{j-1/2}\tilde{u}\_{j-1/2}\tilde{p}\_{j-1/2}\right) \\\\ \sigma\_{j+1/2} = \frac{\left(r\_{j+1/2}^{n+1}\right)^{\alpha+1} - \left(r\_{j+1/2}^{n}\right)^{\alpha+1}}{(\alpha+1)\left(r\_{j+1/2}^{n+1} - r\_{j+1/2}^{n}\right)} \end{cases} \tag{31}$$

#### 2.3.2. Remap step

ð28Þ

j

<sup>j</sup> ð Þ m in cells are

dm ð29Þ

, and the position of new grids

can be computed by equation of state, the

for the

And then the one-dimensional Eq. (27) in each direction is resolved by the two-step Lagrange-Remapping algorithm. Also, one time step is divided into four substeps: (i) the piecewise parabolic interpolating of physical quantities, (ii) solving the Riemann problems approximately, (iii) marching of the Lagrange equations, and (iv) remapping the physical quantities back to the stationary Euler meshes. The orders of accuracy of the spatial and temporal

<sup>∂</sup><sup>t</sup> <sup>¼</sup> <sup>0</sup> <sup>s</sup> <sup>¼</sup> <sup>1</sup>, <sup>2</sup>, …, N � <sup>1</sup>

, mj�1/2 and mjþ1/2 are the mass coordinates of both sides of the

1

CCA

where τ is the specific volume, α is 0, 1, and 2 corresponding to the plane, axial symmetry, and spherical symmetry problems, r is the Lagrange spatial coordinates, m is the Lagrange mass

<sup>ð</sup> mjþ1=<sup>2</sup>

0

BB@

τð Þ m, tn u m, t ð Þ<sup>n</sup> E m, t ð Þ<sup>n</sup> Y<sup>ð</sup>s<sup>Þ</sup>

<sup>j</sup> at time t. The piecewise parabolic function Q<sup>n</sup>

, Enþ<sup>1</sup> j n o

ð Þ m, tn

mj�1=<sup>2</sup>

Because of the calculations are referred to as scalars, here the averaged physical quantities in a

constructed to compute the time average of physical quantity Q on both sides of the grid edge mjþ1/2,Q<sup>~</sup> <sup>j</sup>þ1=2,L and <sup>Q</sup><sup>~</sup> <sup>j</sup>þ1=2,R. Then the Riemann problem at the grid edge mjþ1/2 is solved by using double shock wave approximation. After Lagrange marching step, we can obtain distri-

> , unþ<sup>1</sup> j n o

<sup>j</sup>th grid, <sup>Δ</sup><sup>m</sup> <sup>¼</sup> mjþ1/2 � mj�1/2. The mass average of physical quantities <sup>τ</sup>, u, E, <sup>Y</sup>ð Þ<sup>s</sup> � �<sup>n</sup>

schemes are the third and second, respectively, for smooth flows.

68 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

The Lagrange equations in 1D for multi-fluid can be written as

8

>>>>>>>>>>>><

>>>>>>>>>>>>:

<sup>α</sup>dr, <sup>u</sup> <sup>¼</sup> <sup>∂</sup><sup>r</sup> ∂t

> τ u E <sup>Y</sup>ð Þ<sup>s</sup>

1

n

<sup>¼</sup> <sup>1</sup> Δmj

CCA

j

j n o

> j n o

will not change. The marching formula is as follows

0

BB@

at time tnþ1, the pressure pnþ<sup>1</sup>

∂τ ∂t � <sup>∂</sup> <sup>r</sup> <sup>a</sup> ð Þ <sup>u</sup> <sup>∂</sup><sup>m</sup> <sup>¼</sup> <sup>0</sup>

∂u ∂t þ r <sup>a</sup> ∂p <sup>∂</sup><sup>m</sup> <sup>¼</sup> <sup>0</sup>

∂E ∂t þ ∂ r <sup>a</sup> ð Þ up <sup>∂</sup><sup>m</sup> <sup>¼</sup> <sup>0</sup>

<sup>∂</sup>Yð Þ<sup>s</sup>

2.3.1. Lagrange step of finite volume method

coordinates m ¼

r nþ1 j�1=2 n o

<sup>Y</sup>ð Þ<sup>s</sup> � �<sup>n</sup>þ<sup>1</sup> j � �

jth grid can be defined as

ðr r0 ρr

cell are written as a uniformity Q<sup>n</sup>

bution of the physical quantities τ<sup>n</sup>þ<sup>1</sup>

In Lagrange step, the computational cells distort to follow the material motion, and in Remap step, the averaged physical quantities at the distorted Lagrange cells are remapped back to the stationary Euler meshes. The piecewise parabolic interpolation and integral methods in this step are the same as the ones in Lagrange step.

After Remap step, we define the physical quantities in Euler mesh as <sup>ρ</sup>Euler � �<sup>n</sup>þ<sup>1</sup> <sup>j</sup> , uð Þ Euler <sup>n</sup>þ<sup>1</sup> <sup>j</sup> , ð Þ <sup>E</sup>Euler <sup>n</sup>þ<sup>1</sup> <sup>j</sup> , Yð Þ<sup>s</sup> � � Euler � �<sup>n</sup>þ<sup>1</sup> j , the volume of fluid through the grid boundary j þ 1=2 as ΔV� <sup>j</sup>þ1=<sup>2</sup>, and the average of physical quantities as ρ� , u� , p� , <sup>Y</sup>ð Þ<sup>s</sup> � � � �� jþ1=2 . The density after Lagrange step is ρ<sup>n</sup>þ<sup>1</sup> <sup>j</sup> <sup>¼</sup> <sup>1</sup>=τ<sup>n</sup>þ<sup>1</sup> <sup>j</sup> . In Remap step

$$
\Delta m\_{\rangle}^{n+1} = \rho\_{\rangle}^{n+1} V\_{\rangle}^{n+1} \tag{32}
$$

$$
\Delta m\_{\rangle}^{n} = \left(\rho\_{\text{Euler}}\right)\_{\rangle}^{n} V\_{\text{j}}^{n} \tag{33}
$$

$$
\Delta m\_{\dot{\jmath}+1/2}^{\*} = \rho\_{\dot{\jmath}+1/2}^{\*} \Delta V\_{\dot{\jmath}+1/2}^{\*}\tag{34}
$$

$$
\Delta E\_{j+1/2}^\* = \frac{1}{2} \left( u\_{j+1/2}^\* \right)^2 \Delta m\_{j+1/2}^\* + \Delta e\_{j+1/2}^\* \tag{35}
$$

$$\left(\rho\_{\text{Euler}}\right)\_{j}^{n+1} = \left[\Delta m\_{j}^{n+1} - \left(\rho\_{j+1/2}^{\*}\Delta V\_{j+1/2}^{\*} - \rho\_{j-1/2}^{\*}\Delta V\_{j-1/2}^{\*}\right)\right] \Big/ V\_{j}^{n} \tag{36}$$

$$\left(\left(\mu\_{\text{Euler}}\right)\_{j}^{n+1} = \left[\mathfrak{u}\_{j}^{n+1}\Delta m\_{j}^{n+1} - \left(\mathfrak{u}\_{j+1/2}^{\*}\Delta m\_{j+1/2}^{\*} - \mathfrak{u}\_{j-1/2}^{\*}\Delta m\_{j-1/2}^{\*}\right)\right] \Big/\Delta m\_{j}^{n} \tag{37}$$

$$\left( \left( Y^{\{\mathbf{s}\}} \right)\_{\text{Euler}} \right)\_{j}^{n+1} = \left[ \left( Y^{\{\mathbf{s}\}} \right)\_{j}^{n+1} V\_{j}^{n+1} - \left( \left( Y^{\{\mathbf{s}\}} \right)\_{j+1/2}^{\*} \Delta V\_{j+1/2}^{\*} - \left( Y^{\{\mathbf{s}\}} \right)\_{j-1/2}^{\*} \Delta V\_{j-1/2}^{\*} \right) \right] \Big/ V\_{j}^{n} \tag{38}$$

$$\left(\left(E\_{\text{Euler}}\right)\_{j}^{n+1} = \left[E\_{j}^{n+1}\Delta m\_{j}^{n+1} - \left(\Delta E\_{j+1/2}^{\*} - \Delta E\_{j-1/2}^{\*}\right)\right] / \Delta m\_{j}^{n} \tag{39}$$

where Δe� <sup>j</sup>þ1=<sup>2</sup> is the advection of specific energy through the grid boundary <sup>j</sup> <sup>þ</sup> <sup>1</sup>=2.

#### 2.3.3. Viscous flux, heat flux, and scalar flux

The viscous flux, heat flux, and scalar flux of Eq. (26) are calculated based on the computed inviscid flux by using second-order spatial center difference and two-step Runge-Kutta time marching. The first equation of Eq. (26) can be neglected, and which is written in conserved form

$$\frac{\partial \overline{\mathbf{U}}}{\partial t} + \frac{\partial \overline{\mathbf{F}}}{\partial x} + \frac{\partial \overline{\mathbf{G}}}{\partial y} + \frac{\partial \overline{\mathbf{H}}}{\partial z} = \mathbf{0} \tag{40}$$

and

$$\begin{aligned} \mathbf{U} &= \left(\rho u, \rho v, \rho w, \rho \mathbf{E}\right)^{T} \\ \mathbf{\overline{F}} &= \left(-\sigma\_{xx}, -\sigma\_{yx}, -\sigma\_{zx}, -u\sigma\_{xx} - v\sigma\_{yx} - w\sigma\_{zx} + q\_{x}\right)^{T} \\ \mathbf{\overline{G}} &= \left(-\sigma\_{xy}, -\sigma\_{yy}, -\sigma\_{zy}, -u\sigma\_{xy} - v\sigma\_{yy} - w\sigma\_{zy} + q\_{y}\right)^{T} \\ \mathbf{\overline{H}} &= \left(-\sigma\_{xz}, -\sigma\_{yz}, -\sigma\_{zz}, -u\sigma\_{xz} - v\sigma\_{yz} - w\sigma\_{zz} + q\_{z}\right)^{T} \end{aligned} \tag{41}$$

In the Cartesian coordinate system, the spatial derivative terms of Eq. (40) can be dispersed as follows:

$$\mathbf{L}\_{h}\left(\overline{\mathbf{U}}\_{i,j,k}^{n}\right) = \frac{\mathbf{F}\_{i-1/2,j,k}^{n} - \mathbf{F}\_{i+1/2,j,k}^{n}}{\Delta x} + \frac{\overline{\mathbf{G}}\_{i,j-1/2,k}^{n} - \overline{\mathbf{G}}\_{i,j+1/2,k}^{n}}{\Delta y} + \frac{\mathbf{H}\_{i,j,k-1/2}^{n} - \mathbf{H}\_{i,j,k+1/2}^{n}}{\Delta z} \tag{42}$$

Then, Eq. (40) is solved by two Runge-Kutta time marching,

$$\begin{cases} \overline{\mathbf{U}}\_{i,j,k}^{(l)} = \overline{\mathbf{U}}\_{i,j,k}^{\mathrm{E}} + \Delta t \mathbf{L}\_{h} \left( \mathbf{U}\_{i,j,k}^{n} \right) \\ \mathbf{U}\_{i,j,k}^{n+1} = \frac{1}{2} \left[ \overline{\mathbf{U}}\_{i,j,k}^{\mathrm{E}} + \overline{\mathbf{U}}\_{i,j,k}^{(l)} + \Delta t \mathbf{L}\_{h} \left( \overline{\mathbf{U}}\_{i,j,k}^{(l)} \right) \right] \end{cases} \tag{43}$$

The detailed descriptions of numerical algorithm are referred to Ref. [17].

#### 3. Verification and validation

In this section, the validity and reliability of our compressible multi-fluid dynamic programs are to be verified and validated by comparisons with analytical solutions and hydrodynamic interface instability experiments in a shock tube.

Figure 2. Comparison of postshock density (a), pressure (b), and velocity (c) in copper at t ¼ 1.0 μs when a copper pellet collides with a copper target at three different velocities.

## 3.1. Riemann problem of condensed matter

EEuler � �<sup>n</sup>þ<sup>1</sup>

2.3.3. Viscous flux, heat flux, and scalar flux

where Δe�

form

and

follows:

<sup>L</sup><sup>h</sup> <sup>U</sup><sup>n</sup> i,j, k � � ¼ Fn

3. Verification and validation

interface instability experiments in a shock tube.

<sup>i</sup>�1=2,j, <sup>k</sup> � <sup>F</sup><sup>n</sup>

Δx

Then, Eq. (40) is solved by two Runge-Kutta time marching,

U<sup>ð</sup>l<sup>Þ</sup> i,j, <sup>k</sup> <sup>¼</sup> <sup>U</sup><sup>E</sup>

8 ><

>:

U<sup>n</sup>þ<sup>1</sup> i,j, <sup>k</sup> <sup>¼</sup> <sup>1</sup>

iþ1=2,j, k

þ Gn

<sup>2</sup> <sup>U</sup><sup>E</sup>

The detailed descriptions of numerical algorithm are referred to Ref. [17].

i,j, <sup>k</sup> <sup>þ</sup> <sup>Δ</sup>tL<sup>h</sup> <sup>U</sup><sup>n</sup>

i,j, <sup>k</sup> <sup>þ</sup> <sup>U</sup><sup>ð</sup>l<sup>Þ</sup>

In this section, the validity and reliability of our compressible multi-fluid dynamic programs are to be verified and validated by comparisons with analytical solutions and hydrodynamic

<sup>j</sup> <sup>¼</sup> <sup>E</sup><sup>n</sup>þ<sup>1</sup>

70 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

<sup>j</sup> Δm<sup>n</sup>þ<sup>1</sup>

∂U ∂t þ

<sup>U</sup> <sup>¼</sup> <sup>ρ</sup>u,ρv,ρw,ρ<sup>E</sup> � �<sup>T</sup>

∂F ∂x þ ∂G ∂y þ ∂H

F ¼ �σxx, � σyx, � σzx, � uσxx � vσyx � wσzx þ qx � �<sup>T</sup>

G ¼ �σxy, � σyy, � σzy, � uσxy � vσyy � wσzy þ qy � �<sup>T</sup>

H ¼ �σxz, � σyz, � σzz, � uσxz � vσyz � wσzz þ qz � �<sup>T</sup>

In the Cartesian coordinate system, the spatial derivative terms of Eq. (40) can be dispersed as

i,j�1=2, <sup>k</sup> � <sup>G</sup><sup>n</sup>

Δy

i,j, k � �

h i � �

i,j, <sup>k</sup> <sup>þ</sup> <sup>Δ</sup>tL<sup>h</sup> <sup>U</sup><sup>ð</sup>l<sup>Þ</sup>

i,jþ1=2, k

þ Hn

i,j, k

<sup>j</sup> � ΔE�

<sup>j</sup>þ1=<sup>2</sup> is the advection of specific energy through the grid boundary <sup>j</sup> <sup>þ</sup> <sup>1</sup>=2.

The viscous flux, heat flux, and scalar flux of Eq. (26) are calculated based on the computed inviscid flux by using second-order spatial center difference and two-step Runge-Kutta time marching. The first equation of Eq. (26) can be neglected, and which is written in conserved

h i � �

<sup>j</sup>þ1=<sup>2</sup> � <sup>Δ</sup>E�

j�1=2

=Δm<sup>n</sup>

<sup>∂</sup><sup>z</sup> <sup>¼</sup> <sup>0</sup> <sup>ð</sup>40<sup>Þ</sup>

i,j, <sup>k</sup>�1=<sup>2</sup> � <sup>H</sup><sup>n</sup>

i,j, kþ1=2 <sup>Δ</sup><sup>z</sup> <sup>ð</sup>42<sup>Þ</sup>

<sup>j</sup> ð39Þ

ð41Þ

ð43Þ

A copper pellet collides with a copper target with three velocities of 2, 4, and 8 km/s. Using the reduction form of Gruneisen EOS for copper in simulations, the one-dimensional numerical results of postshock density (a), pressure (b), and velocity (c) compare with theoretical solutions, as shown in Figure 2, the solid lines corresponding to numerical results and the dot lines corresponding to the theoretical solutions.

## 3.2. Riemann problem of one-dimensional gas/liquid

At initial time, the region [0, 1.0 cm] is filled with gas with high pressure 1.0 � 108 Pa and density 1.29 g/cm3 . The region [1.0 cm, 5.0 cm] is filled with water with pressure 1.0 � 105 Pa and density 1.0 g/cm3 . The gas and water are all described by Stiffen gas EOS. The left and right boundaries are flow. Figure 3 shows the distributions of the density (a), pressure (b), and velocity (c) at 20 μs. There are a forward shock wave and a backward rarefaction wave after interaction. The pressure and velocity around the interface are well continuous and have no nonphysical oscillation.

## 3.3. Single-mode Richtmyer-Meshkov instability

The two- and three-dimensional single-mode Richtmyer-Meshkov instabilities are numerical simulated by MVPPM program, which also compare with the theoretical model [18]. The initial small perturbation is a sinusoidal one with wavelength 60 mm (global wavelength for a three-dimensional case) and amplitude 1.0 mm. The incident air shock wave with Mach number 1.2 impacts the air/SF6 single-mode interface. Figure 4 shows the comparisons of amplitude of single-mode perturbed interface with the linear impulsive and nonlinear models, the left one corresponding to the two-dimensional (2D) results, and the right one corresponding to the three-dimensional results.

#### 3.4. Simulations of shock tube experiments

The compressible multi-fluid dynamic programs are used to simulated several shock tube experiments of interface instability for validation. These experiments include planar interface and gas cylinder shock tube experiment and planar and cylindrical jelly experiments.

Figure 5 shows the comparison of 2D calculated width of turbulent mixing zone (TMZ) [19] of Leinov's shock tube experiment with reshock with experiment [20] in which the Mach number

Figure 3. Distributions of the density (a), pressure (b), and velocity (c) at 20 μs for one-dimensional Riemann problem of gas/liquid interface.

Figure 4. Comparisons of amplitude of single-mode perturbed interface with the linear impulsive model and nonlinear models, left: two-dimension and right: three-dimension.

Figure 5. TMZ width versus time (t ¼ 0 denotes the reshock arrival to the interface).

of incident air shock wave is 1.2, which impacts the air/SF6 interface, the grid size is 50 μm. Figure 6 shows the 3D calculated TMZ width [21] of Poggi's multi-mode shock tube experiment [22] with reshock in which the Mach number of incident SF6 shock wave is 1.453, which impacts the SF6/air interface. Figures 7 [10] and 8 [23] show the calculated and experimental interface images of AWE's SF6 half-height and double-bump shock tube experiment [24, 25], the Mach number of incident air shock wave is 1.26, respectively. Figure 9 [26] shows the SF6 gas cylinder evolution at different times under the initial air shock wave with the Mach number 1.2, Figure 10 shows the width and height of gas cylinder at different times for experiment and numerical simulations.

Figure 6. Time history of TMZ width.

Figure 3. Distributions of the density (a), pressure (b), and velocity (c) at 20 μs for one-dimensional Riemann problem of

72 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

Figure 4. Comparisons of amplitude of single-mode perturbed interface with the linear impulsive model and nonlinear

gas/liquid interface.

models, left: two-dimension and right: three-dimension.

Figure 5. TMZ width versus time (t ¼ 0 denotes the reshock arrival to the interface).

Figure 7. Two dimension calculated from MVFT and experimental half-height interface at different times.

Figure 8. Three dimension calculated from MVFT and experimental double-bump interface at different times, left column corresponding to experimental images, right column corresponding to 3D calculated images, and middle column corresponding to span-wise average images of 3D calculated results.

Figure 9. Evolutions of SF6 gas cylinder, upper row corresponding to experiment, lower row corresponding to numerical simulation by MVFT program, the time sequences are 0, 50, 190, 330, 470, 610, and 750 μs severally.

For the problems of interface instability with a high density ratio, our hydrodynamic programs are also applicable. For the LLNL's jelly ring experiments [27], the jelly ring only has a sinusoidal periodic initial perturbation at the outer interface and is driven by expansion of the explosion products of a gaseous mixture of C2H2 and O2. The jelly is mainly made of water. The thickness of ring is 15 mm, the mode number of perturbation is 6 and 36, respectively, and the initial amplitude is 1.0 mm. Figure 11 shows the evolutions of the jelly ring at different times from simulated from LLNL' CALE program (a) and our MVPPM program (b), and the experimental image (c) at 600 μs [28]. Figure 12 shows the time histories of radius (a) and the amplitude (b) of the outer and inner interfaces of the jelly ring simulated from CALE and

Figure 10. Width and height of gas cylinder at different times for experimental and numerical simulations.

Figure 11. Evolutions of jelly ring at different simulated from CALE (a: left group) and MVPPM (b: middle group) programs, and the experimental image at 600 μs (c).

Figure 12. Time histories of radius (a) and the amplitude (b) of the outer and inner interfaces of jelly ring simulated from CALE and MVPPM programs.

For the problems of interface instability with a high density ratio, our hydrodynamic programs are also applicable. For the LLNL's jelly ring experiments [27], the jelly ring only has a sinusoidal periodic initial perturbation at the outer interface and is driven by expansion of the explosion products of a gaseous mixture of C2H2 and O2. The jelly is mainly made of water. The thickness of ring is 15 mm, the mode number of perturbation is 6 and 36, respectively, and the initial amplitude is 1.0 mm. Figure 11 shows the evolutions of the jelly ring at different times from simulated from LLNL' CALE program (a) and our MVPPM program (b), and the experimental image (c) at 600 μs [28]. Figure 12 shows the time histories of radius (a) and the amplitude (b) of the outer and inner interfaces of the jelly ring simulated from CALE and

Figure 9. Evolutions of SF6 gas cylinder, upper row corresponding to experiment, lower row corresponding to numerical

simulation by MVFT program, the time sequences are 0, 50, 190, 330, 470, 610, and 750 μs severally.

Figure 8. Three dimension calculated from MVFT and experimental double-bump interface at different times, left column corresponding to experimental images, right column corresponding to 3D calculated images, and middle column

corresponding to span-wise average images of 3D calculated results.

74 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

MVPPM programs [28]. Figure 13 shows the images of jelly ring at 600 μs for mode number 36 from experiment (a), CALE (b), and MVPPM (c) programs [28].

From the above comparisons between our numerical simulations and theoretical model, experiments for the problems of hydrodynamic interface instability from the density ratio low

Figure 13. Images of jelly ring at 600 μs for mode number 36 from experiment (a), CALE program (b), and MVPPM program (c).

to the high density ratio, the numerical results agree well with theory and experiments, so the validity and reliability of our compressible multi-fluid dynamic programs have been verified and validated.

## 4. Interface instability and turbulent mixing

As we know, the physical mechanism for the occurrence of Richtmyer-Meshkov instability is the baroclinic vorticity deposition at the interface resulting from the misalignment of the pressure gradient across the shock front and the local density gradient across the interface. The evolution equation of vorticity is as follows:

$$\frac{d\boldsymbol{\omega}}{dt} = \frac{\nabla\rho \times \nabla p}{\rho^2} + \boldsymbol{\omega} \cdot \nabla \mathbf{u} - \boldsymbol{\omega} \nabla \cdot \mathbf{u} \tag{44}$$

where ω ¼ ∇ � u is the vorticity and viscous terms are neglected. The first term on the right side of Eq. (44) is the baroclinic vorticity production term. The second term on the right side of Eq. (44) is the vortex-stretching term, which is zero in the two-dimensional case, as the vorticity and velocity fields are orthogonal. This term enhances dissipation, resulting in more diffuse and smaller scale structures in the turbulent mixing zone. The third term on the right side is the compression term and does not contribute to the vorticity evolution significantly. The baroclinic vorticity production is much larger when the shock wave impacts the interface and pass through it and constitutes the principal mechanism of the Richtmyer-Meshkov instability.

In this section, we will introduce the growth laws of the Richtmyer-Meshkov instability and the induced turbulent mixing and its dynamic behavior by numerical simulations.

#### 4.1. Growth laws of Richtmyer-Meshkov instability and turbulent mixing

Figure 14 shows the amplitude (a) and growth rate (b) of two-dimensional single-mode Richtmyer-Meshkov instability without reshock for initial perturbation amplitude 1 and 12 mm

Figure 14. Amplitude (a) and growth rate (b) of two-dimensional single-mode Richtmyer-Meshkov instability without reshock.

Figure 15. Schematic of shock tube and computational model.

to the high density ratio, the numerical results agree well with theory and experiments, so the validity and reliability of our compressible multi-fluid dynamic programs have been verified

Figure 13. Images of jelly ring at 600 μs for mode number 36 from experiment (a), CALE program (b), and MVPPM

As we know, the physical mechanism for the occurrence of Richtmyer-Meshkov instability is the baroclinic vorticity deposition at the interface resulting from the misalignment of the pressure gradient across the shock front and the local density gradient across the interface.

where ω ¼ ∇ � u is the vorticity and viscous terms are neglected. The first term on the right side of Eq. (44) is the baroclinic vorticity production term. The second term on the right side of Eq. (44) is the vortex-stretching term, which is zero in the two-dimensional case, as the vorticity and velocity fields are orthogonal. This term enhances dissipation, resulting in more diffuse and smaller scale structures in the turbulent mixing zone. The third term on the right side is the compression term and does not contribute to the vorticity evolution significantly. The baroclinic vorticity production is much larger when the shock wave impacts the interface and pass through it and constitutes the principal mechanism of the Richtmyer-Meshkov

In this section, we will introduce the growth laws of the Richtmyer-Meshkov instability and

Figure 14 shows the amplitude (a) and growth rate (b) of two-dimensional single-mode Richtmyer-Meshkov instability without reshock for initial perturbation amplitude 1 and 12 mm

the induced turbulent mixing and its dynamic behavior by numerical simulations.

4.1. Growth laws of Richtmyer-Meshkov instability and turbulent mixing

<sup>ρ</sup><sup>2</sup> <sup>þ</sup> <sup>ω</sup> � <sup>∇</sup><sup>u</sup> � <sup>ω</sup><sup>∇</sup> � <sup>u</sup> <sup>ð</sup>44<sup>Þ</sup>

4. Interface instability and turbulent mixing

dω

76 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

dt <sup>¼</sup> <sup>∇</sup><sup>ρ</sup> � <sup>∇</sup><sup>p</sup>

The evolution equation of vorticity is as follows:

and validated.

program (c).

instability.

and wavelength 60 mm. After initial shock, the perturbation enters the nonlinear stage quickly. The growth rate increases fast and reaches the highest peak, then will reduce owing to the reduction of the effect of compressibility and the dominant role of flow nonlinearity. At the late times, the amplitude is increasing linearly, and the growth rate remains a constant. For the larger initial amplitude, the perturbation enters the later linear growth stage earlier.

For the multi-mode Richtmyer-Meshkov instability and the induced turbulent mixing with reshock, the initial air shock wave with Mach number 1.2 impacts the air/SF6 interface. The multi-mode initial perturbation is composed of eight dominant mode wavelengths of 0.8, 1.0, 1.25, 1.6, 2.0, 2.5, 3.2, and 4.0 mm superimposed with a very small random disturbance. The shock tube experiment can be referred to Ref. [29]. Figure 15 shows the schematic of shock tube and computational model. The transmitted wave rebounds between the interface and the end-wall of shock tube and produces multiple shock-interface interactions. At about 1.7 ms,

Figure 16. Time history of the TMZ width.

the transmitted shock wave reflected back off the end-wall impacts the interface. Figure 16 shows the time history of TMZ width, the black line denotes the numerically calculated width and the red line is obtained from fitting of the numerical results. As can be seen, after the initial shock, the TMZ width starts to grow as a power law ~t <sup>θ</sup> with the value θ to be determined as 0.352. After the reshock, more energy is deposited onto the interface to promote the development of TMZ, and the TMZ width evolves in time as a negative exponential law ~�e�<sup>t</sup>=<sup>t</sup> � where the value of t \* is 0.519. Then after the following interaction of the reflected rarefaction wave with the interface, the TMZ also evolves as a negative exponential law but with a different factor t \* <sup>¼</sup> 0.875. Under the subsequent impingements with lower and lower intensity, the TMZ width, after a slight reduction caused by the reflected compression wave, evolves in an approximate linear fashion with a velocity of 2.05 m/s. Figure 17 shows the instantaneous images of the TMZ visualized by the volume fraction isosurface YSF6 ¼ 0.1, 0.3, 0.5, 0.7, and 0.9 from blue to orange at different times, the TMZ exhibits a very complex spatial structure.

#### 4.2. Evaluation of different subgrid-scale stress models

In large-eddy simulation, the effect of small scales on large-scale motions is represented by the SGS stress model. Most of the commonly used SGS models assume that the main function of subgrid scales is to remove energy from the large scales and dissipate it through the action of the viscous forces. But, as we know, in fact the energy is also transferred from the small scales to the large scales (backscatter) in a small and local range. The SGS turbulent dissipation, which is the work of SGS stress, represents the energy transfer between resolved and subgrid scales,

$$
\varepsilon\_{\rm scs} = \tau\_{\rm j\overline{}} \overline{\mathbf{S}}\_{\rm ij} \tag{45}
$$

If it is negative, the subgrid scales remove energy from the resolved scales (forward scatter); if it is positive, they release energy to the resolved scales (backscatter). It is easy to see that the eddy viscosity models such as Smagorinsky model, Vreman model, etc. are absolutely dissipative.

Figure 17. Instantaneous images of TMZ visualized by the volume fraction isosurface YSF6 ¼ 0.1, 0.3, 0.5, 0.7, and 0.9 from blue to orange at different times.

the transmitted shock wave reflected back off the end-wall impacts the interface. Figure 16 shows the time history of TMZ width, the black line denotes the numerically calculated width and the red line is obtained from fitting of the numerical results. As can be seen, after the initial

0.352. After the reshock, more energy is deposited onto the interface to promote the development of TMZ, and the TMZ width evolves in time as a negative exponential law ~�e�<sup>t</sup>=<sup>t</sup>

wave with the interface, the TMZ also evolves as a negative exponential law but with a

the TMZ width, after a slight reduction caused by the reflected compression wave, evolves in an approximate linear fashion with a velocity of 2.05 m/s. Figure 17 shows the instantaneous images of the TMZ visualized by the volume fraction isosurface YSF6 ¼ 0.1, 0.3, 0.5, 0.7, and 0.9 from blue to orange at different times, the TMZ exhibits a very complex spatial structure.

In large-eddy simulation, the effect of small scales on large-scale motions is represented by the SGS stress model. Most of the commonly used SGS models assume that the main function of subgrid scales is to remove energy from the large scales and dissipate it through the action of the viscous forces. But, as we know, in fact the energy is also transferred from the small scales to the large scales (backscatter) in a small and local range. The SGS turbulent dissipation, which is the work of SGS stress, represents the energy transfer between resolved and subgrid scales,

If it is negative, the subgrid scales remove energy from the resolved scales (forward scatter); if it is positive, they release energy to the resolved scales (backscatter). It is easy to see that the eddy viscosity models such as Smagorinsky model, Vreman model, etc. are absolutely

\* is 0.519. Then after the following interaction of the reflected rarefaction

\* <sup>¼</sup> 0.875. Under the subsequent impingements with lower and lower intensity,

<sup>θ</sup> with the value θ to be determined as

εSGS ¼ τijSij ð45Þ

�

shock, the TMZ width starts to grow as a power law ~t

78 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

4.2. Evaluation of different subgrid-scale stress models

where the value of t

Figure 16. Time history of the TMZ width.

different factor t

dissipative.

Figure 18. Distribution of calculated SGS turbulent dissipation in a streamwise direction for Smagorinsky and Vreman models and for AWE's SF6 half-height shock tube experiment [24] at two times. Red line corresponding to the Vreman model, green line corresponding to the Smagorinsky model.

Figure 18 shows the distribution of calculated SGS turbulent dissipation in streamwise direction for Smagorinsky and Vreman models and for AWE's SF6 half-height shock tube experiment [24] at two times [11]. The SGS dissipation of the Smagorinsky model is much greater than the Vreman model over 1.5 times; therefore, the dissipation is too great for the Smagorinsky SGS model. The SGS turbulent dissipations of the Vreman model (VM), the dynamic Smagorinsky model (DSM), and the stretched-vortex model (SVM) based on a planar Richtmyer-Meshkov instability with incident Mach number 1.2 are shown in Figure 19 [30]. Before the interfacial flow has developed to be turbulent completely, the dynamic and stretched-vortex models have all predicted the energy backscatter, but the energy backscatter predicted by the dynamic model is

Figure 19. Distribution of the SGS turbulent dissipation in a streamwise direction for the Vreman, dynamic Smagorinsky and stretched-vortex models at four times.

larger and its range is wider. After reshock, the turbulent fluctuations are stronger extremely the turbulent dissipation also increases. The dynamic model's dissipation is the highest, then followed by the Vreman model, and the stretched-vortex model's dissipation is the lowest. At the late time, the SGS dissipation of the dynamic and Vreman models is the same basically, and the dynamic model is also to be absolutely dissipative, yet the stretched-vortex model is still able to predict the local phenomenon of energy backscatter in a small range. So, the dynamic model is poor in representing the energy backscatter. The Vreman and stretched-vortex models are all robust, but the former is absolutely dissipative.

#### 4.3. Effects of the initial conditions on the growth of RMI and the turbulent mixing

#### 4.3.1. Effects of the initial conditions on the growth of single-mode RMI

First, we consider the effects of initial conditions of perturbation amplitude and wavelength on the growth of single-mode Richtmyer-Meshkov instability without reshock [18]. The incident shock wave with Mach number 1.2 impacts the air/SF6 interface. The initial perturbation amplitude and wavelength are listed in Table 1. The calculations are carried out in two dimensions, and the RMI does not develop into turbulent mixing completely. Figure 20 shows the effects of the initial perturbation amplitude on the growth of single-mode RMI for the fixed initial perturbation wavelength λ<sup>0</sup> ¼ 60 mm. The perturbation amplitude (a) and growth rate (b) increase gradually with the increasing of initial amplitude. And when the initial amplitude is much larger, the growth of amplitude enters a linear stage earlier at the late times, and the growth rate remains a constant. The growth rate increases fast and reaches the highest peak at the early times. After the peak, the effect of compressibility is reduced and the flow nonlinearity starts to play a dominant role and causes the growth rate to decay with time. Figure 21 shows the effects of the initial perturbation wavelength on the growth of single-mode RMI for the fixed initial


Table 1. Initial conditions of perturbation for single-mode RMI.

larger and its range is wider. After reshock, the turbulent fluctuations are stronger extremely the turbulent dissipation also increases. The dynamic model's dissipation is the highest, then followed by the Vreman model, and the stretched-vortex model's dissipation is the lowest. At the late time, the SGS dissipation of the dynamic and Vreman models is the same basically, and the dynamic model is also to be absolutely dissipative, yet the stretched-vortex model is still able to predict the local phenomenon of energy backscatter in a small range. So, the dynamic model is poor in representing the energy backscatter. The Vreman and stretched-vortex models are all

Figure 19. Distribution of the SGS turbulent dissipation in a streamwise direction for the Vreman, dynamic Smagorinsky

4.3. Effects of the initial conditions on the growth of RMI and the turbulent mixing

First, we consider the effects of initial conditions of perturbation amplitude and wavelength on the growth of single-mode Richtmyer-Meshkov instability without reshock [18]. The incident shock wave with Mach number 1.2 impacts the air/SF6 interface. The initial perturbation amplitude and wavelength are listed in Table 1. The calculations are carried out in two dimensions, and the RMI does not develop into turbulent mixing completely. Figure 20 shows the effects of the initial perturbation amplitude on the growth of single-mode RMI for the fixed initial perturbation wavelength λ<sup>0</sup> ¼ 60 mm. The perturbation amplitude (a) and growth rate (b) increase gradually with the increasing of initial amplitude. And when the initial amplitude is much larger, the growth of amplitude enters a linear stage earlier at the late times, and the growth rate remains a constant. The growth rate increases fast and reaches the highest peak at the early times. After the peak, the effect of compressibility is reduced and the flow nonlinearity starts to play a dominant role and causes the growth rate to decay with time. Figure 21 shows the effects of the initial perturbation wavelength on the growth of single-mode RMI for the fixed initial

4.3.1. Effects of the initial conditions on the growth of single-mode RMI

80 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

robust, but the former is absolutely dissipative.

and stretched-vortex models at four times.

Figure 20. Effects of the initial perturbation amplitude on the growth of single-mode RMI for the fixed initial perturbation wavelength λ<sup>0</sup> ¼ 60 mm.

Figure 21. Effects of the initial perturbation wavelength on the growth of single-mode RMI for the fixed initial perturbation amplitude a<sup>0</sup> ¼ 3 mm.

perturbation amplitude a<sup>0</sup> ¼ 3 mm. The perturbation amplitude (a) and growth rate (b) reduce gradually at the early times and increase gradually at the late times with the increasing of initial wavelength. And when the initial perturbation strength (ration of initial amplitude to wavelength) is much smaller, the perturbation growth is mainly dependent on the initial perturbation amplitude and slightly dependent on the initial perturbation wavelength at the late times of RMI. The same conclusion can be obtained from three-dimensional calculations.

## 4.3.2. Effects of the initial conditions on the growth of multi-mode RMI and the induced turbulent mixing

For the effects of three-dimensional initial multi-mode conditions, the case of Richtmyer-Meshkov instability is same as the above Section 4.1. Table 2 lists the initial condition of perturbations, the perturbation strength (PS) is defined as the ratio of initial amplitude to wavelength. Turbulent kinetic energy K, dissipation rate ε, and enstrophy Ω are defined as follows [31]:

$$K = \frac{\langle \rho u''\_i u''\_i \rangle}{2 \langle \rho \rangle} + \frac{\langle \tau\_{ii} \rangle}{2 \langle \rho \rangle} \tag{46}$$

$$\varepsilon = \frac{\langle \sigma\_{\vec{\eta}} \eth \mathsf{u} \nu'' , \langle \mathsf{Ox}\_{\vec{\eta}} \rangle }{\langle \rho \rangle} - \frac{\langle \pi\_{\vec{\eta}} \eth \mathsf{u} \nu'' , \langle \mathsf{Ox}\_{\vec{\eta}} \rangle}{\langle \rho \rangle} \tag{47}$$

$$
\Omega = \frac{1}{2} \|\boldsymbol{\omega}\|^2 \tag{48}
$$

where u<sup>00</sup> <sup>i</sup> is the velocity fluctuation and 〈 〉 denotes the transverse plane-average. For largeeddy simulations, the turbulent kinetic energy and dissipation rate both include two parts named as the resolved-scales (the first term on the right side of Eqs. (46) and (47)) and the subgrid-scales (the second term on the right side of Eqs. (46) and (47)).

Figure 22 shows the growth history of TMZ width. Figures 23 and 24 show the time histories of the perk values of turbulent kinetic energy and enstrophy, respectively. For the larger initial perturbation strength, the TMZ grows faster, the turbulent kinetic energy is also larger or the turbulence strength is also stronger, the deposited vorticity is larger too. The development of turbulent mixing has a strong dependence on the initial conditions between the initial shock and the impingement of the first reflected rarefaction wave, after that the evolution of the turbulent mixing has lost the memory of the initial conditions.


Table 2. Model parameters.

Figure 22. Growth history of the TMZ width for different models.

perturbation amplitude a<sup>0</sup> ¼ 3 mm. The perturbation amplitude (a) and growth rate (b) reduce gradually at the early times and increase gradually at the late times with the increasing of initial wavelength. And when the initial perturbation strength (ration of initial amplitude to wavelength) is much smaller, the perturbation growth is mainly dependent on the initial perturbation amplitude and slightly dependent on the initial perturbation wavelength at the late times of

For the effects of three-dimensional initial multi-mode conditions, the case of Richtmyer-Meshkov instability is same as the above Section 4.1. Table 2 lists the initial condition of perturbations, the perturbation strength (PS) is defined as the ratio of initial amplitude to wavelength. Turbulent kinetic energy K, dissipation rate ε, and enstrophy Ω are defined as

> iu<sup>00</sup> i〉 <sup>2</sup>〈ρ〉 <sup>þ</sup> 〈τii〉

> > <sup>i</sup>=∂xj〉

<sup>i</sup>=∂xj〉 〈ρ〉 � 〈τij∂u<sup>00</sup>

<sup>Ω</sup> <sup>¼</sup> <sup>1</sup> 2 j ω j

eddy simulations, the turbulent kinetic energy and dissipation rate both include two parts named as the resolved-scales (the first term on the right side of Eqs. (46) and (47)) and the

Figure 22 shows the growth history of TMZ width. Figures 23 and 24 show the time histories of the perk values of turbulent kinetic energy and enstrophy, respectively. For the larger initial perturbation strength, the TMZ grows faster, the turbulent kinetic energy is also larger or the turbulence strength is also stronger, the deposited vorticity is larger too. The development of turbulent mixing has a strong dependence on the initial conditions between the initial shock and the impingement of the first reflected rarefaction wave, after that the evolution of the

η<sup>0</sup> (mm) 0.07 0.14 0.28 0.56 1.12 PS 0.035 .0.7 0.14 0.28 0.56

<sup>i</sup> is the velocity fluctuation and 〈 〉 denotes the transverse plane-average. For large-

CASE1 CASE2 CASE3 CASE4 CASE5

<sup>2</sup>〈ρ〉 <sup>ð</sup>46<sup>Þ</sup>

〈ρ〉 <sup>ð</sup>47<sup>Þ</sup>

<sup>2</sup> <sup>ð</sup>48<sup>Þ</sup>

RMI. The same conclusion can be obtained from three-dimensional calculations.

<sup>K</sup> <sup>¼</sup> 〈ρu<sup>00</sup>

<sup>ε</sup> <sup>¼</sup> 〈σij∂u<sup>00</sup>

subgrid-scales (the second term on the right side of Eqs. (46) and (47)).

turbulent mixing has lost the memory of the initial conditions.

4.3.2. Effects of the initial conditions on the growth of multi-mode RMI and the

82 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

induced turbulent mixing

follows [31]:

where u<sup>00</sup>

Table 2. Model parameters.

Figure 23. History of turbulent kinetic energy for different models.

Figure 24. History of enstrophy for different models.

#### 4.4. Dynamic characters of multi-mode RM instability and induced turbulent mixing

Figures 25–27 show the spatial profiles of the turbulent kinetic energy, dissipation rate and enstrophy along the motion direction of shock wave at different times before and after reshock individually. They all have a spatial profile similar to Gaussian distribution. The strongest turbulent intensity is located in the center of TMZ. Figures 28–30 show the temporal evolutions of the peak values of the turbulent kinetic energy, dissipation rate, and enstrophy in the spatial profiles, along with their fitted results, respectively. The turbulent kinetic energy, dissipation rate, and enstrophy decay gradually because of dissipation and diffusion. After the initial shock and before the reshock, the turbulent kinetic energy and enstrophy decay with time as a power law ~t <sup>θ</sup> except the dissipation rate which decays with time as an exponential law ~e<sup>t</sup>=<sup>t</sup> . One reason is that the TMZ is not fully developed turbulence before reshock and the other reason is that the flow in TMZ is stronger anisotropic between the transverse direction and the axis direction. After reshock and the first reflected rarefaction wave, they all decay with time as the negative exponential law with the closed decay factors at the same

Figure 25. Spatial profiles of the turbulent kinetic energy at different times. (a) Before reshock and (b) after reshock.

Figure 26. Spatial profiles of the turbulent kinetic energy dissipation rate at different times. (a) Before reshock and (b) after reshock.

Figure 27. Spatial profiles of the enstrophy at different times. (a) Before reshock and (b) after reshock.

Figure 28. Time history of the peak value of turbulent kinetic energy.

4.4. Dynamic characters of multi-mode RM instability and induced turbulent mixing

84 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

time as a power law ~t

law ~e<sup>t</sup>=<sup>t</sup>

after reshock.

Figures 25–27 show the spatial profiles of the turbulent kinetic energy, dissipation rate and enstrophy along the motion direction of shock wave at different times before and after reshock individually. They all have a spatial profile similar to Gaussian distribution. The strongest turbulent intensity is located in the center of TMZ. Figures 28–30 show the temporal evolutions of the peak values of the turbulent kinetic energy, dissipation rate, and enstrophy in the spatial profiles, along with their fitted results, respectively. The turbulent kinetic energy, dissipation rate, and enstrophy decay gradually because of dissipation and diffusion. After the initial shock and before the reshock, the turbulent kinetic energy and enstrophy decay with

<sup>θ</sup> except the dissipation rate which decays with time as an exponential

. One reason is that the TMZ is not fully developed turbulence before reshock and

the other reason is that the flow in TMZ is stronger anisotropic between the transverse direction and the axis direction. After reshock and the first reflected rarefaction wave, they all decay with time as the negative exponential law with the closed decay factors at the same

Figure 25. Spatial profiles of the turbulent kinetic energy at different times. (a) Before reshock and (b) after reshock.

Figure 26. Spatial profiles of the turbulent kinetic energy dissipation rate at different times. (a) Before reshock and (b)

stage, and be similar to the growth of TMZ width. And then, they all decay asymptotically due to no remarkable energy deposition. Therefore, the turbulent mixing behaves in a statistical self-similar pattern. Figure 31 shows the one-dimensional global energy spectra of three velocity components on a log-log scale at three times. The energy spectra of two transverse components of velocity are too close, and there is a difference between transverse and axis components. The turbulent mixing flow is continuous anisotropic yet the anisotropy weakens gradually. That is to say, the development of the turbulent mixing presents a trend of isotropy.

#### 4.5. Numerical study of the elliptic gas cylinder instability

Our group first performed the experimental and numerical investigations of the elliptic gas cylinder instability. As we know that the initial density distribution of the gas cylinder is hard

Figure 29. Time history of the peak value of turbulent kinetic energy dissipation rate.

Figure 30. Time history of the peak value of enstrophy.

to determine in the experiment, which can only give the one-dimensional radial concentration distribution for circular gas cylinder as an approximate Gaussian function. Drawing to the experience on the one-dimensional diffusive interfacial transition layer with finite thickness for circular gas cylinder, we constructed a two-dimensional diffusive interfacial transition layer with finite thickness for elliptic gas cylinder through numerical simulation [32],

$$\rho(\alpha, \beta) = \chi\_0 \rho\_{SF\_6} e^{-\left[\left(a - a\_{\rm min}\right)^2 + \left(\beta - \beta\_{\rm min}\right)^2\right]/\delta^2} \tag{49}$$

$$\frac{(\mathbf{x} - \mathbf{x}\_0)}{\alpha^2} + \frac{(\mathbf{y} - \mathbf{y}\_0)}{\beta^2} = 1\tag{50}$$

Figure 31. One-dimensional global energy spectra of three components of velocity on a log-log scale.

where <sup>x</sup><sup>0</sup> <sup>¼</sup> <sup>y</sup><sup>0</sup> <sup>¼</sup> 0, <sup>α</sup> <sup>∈</sup> [αmin, <sup>α</sup>max], <sup>β</sup> <sup>∈</sup> [βmin, <sup>β</sup>max], <sup>α</sup>min <sup>¼</sup> <sup>β</sup>min <sup>¼</sup> 1.0 � <sup>10</sup>�<sup>5</sup> mm, <sup>α</sup>max <sup>¼</sup> 6.30 mm, βmax ¼ 2.30 mm, and δ ¼ 6.16 mm. χ<sup>0</sup> ¼ 0.71 is the concentration of the elliptic center. The density distribution is plotted in Figure 32. The Mach number of air shock wave is 1.25. The simulations (see Figure 33) reproduce the elliptic gas cylinder instability experiment very well, they achieve to a good agreement qualitatively and quantitatively, some salient features of the vortex pairs are obtained clearly.

Figure 34 shows the vorticity at the center of the core and the distance between the two vortex cores of these two simulations. When the incident shock accelerates the elliptic gas cylinder along the major axis, the absolute vorticity in the vortex core |ωcore| is larger; but the distance between two vortex cores is larger when the incident shock accelerates the cylinder along the minor axis. So, the effect of convergence is stronger when the incident shock accelerates the elliptic gas cylinder along the major axis, for which the gas jet appears at the symmetry center.

#### 4.6. Foundation of new RM instability and mechanism

to determine in the experiment, which can only give the one-dimensional radial concentration distribution for circular gas cylinder as an approximate Gaussian function. Drawing to the experience on the one-dimensional diffusive interfacial transition layer with finite thickness for circular gas cylinder, we constructed a two-dimensional diffusive interfacial transition layer

> e �½ðα�αminÞ 2 þðβ�βminÞ 2 �=δ<sup>2</sup>

<sup>α</sup><sup>2</sup> <sup>þ</sup> <sup>ð</sup><sup>y</sup> � <sup>y</sup>0<sup>Þ</sup>

ð49Þ

<sup>β</sup><sup>2</sup> <sup>¼</sup> <sup>1</sup> <sup>ð</sup>50<sup>Þ</sup>

with finite thickness for elliptic gas cylinder through numerical simulation [32],

ðx � x0Þ

ρðα, βÞ ¼ χ0ρSF<sup>6</sup>

Figure 30. Time history of the peak value of enstrophy.

Figure 29. Time history of the peak value of turbulent kinetic energy dissipation rate.

86 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

At present, the investigation for the Richtmyer-Meshkov instability is performed in a uniform flow field. We first study the Richtmyer-Meshkov instability and turbulent mixing in a

Figure 32. Density images and distributions of the (a) circular and (b) elliptic SF6 gas cylinder initially constructed.

Figure 33. Experimental evolution images and numerical simulation results by MVFT at t ¼ 200, 300, 400, 500, and 600 µs, the experiments corresponding (a), (c), and (e), and the simulations corresponding (b), (d), and (f).

nonuniform flow field by experiment and numerical simulations and find some new phenomena. Figure 35 shows the test section schematic of shock tube. The air shock wave with Mach number 1.27 accelerates the two-mode sinusoidal air/SF6 interface (amplitude A<sup>01</sup> ¼ 5 m, A<sup>02</sup> ¼ 7.5 mm). The experimental Schlieren images (gray images in Figure 38) show the transmitted shock wave and the interface all incline which is different from the familiar RMI. We think this

Figure 34. Vorticity at the center of the cores and the distance between the two vortex cores from simulations.

Figure 35. Initial structure diagram in the shock tube.

nonuniform flow field by experiment and numerical simulations and find some new phenomena. Figure 35 shows the test section schematic of shock tube. The air shock wave with Mach number 1.27 accelerates the two-mode sinusoidal air/SF6 interface (amplitude A<sup>01</sup> ¼ 5 m, A<sup>02</sup> ¼ 7.5 mm). The experimental Schlieren images (gray images in Figure 38) show the transmitted shock wave and the interface all incline which is different from the familiar RMI. We think this

Figure 33. Experimental evolution images and numerical simulation results by MVFT at t ¼ 200, 300, 400, 500, and 600 µs,

the experiments corresponding (a), (c), and (e), and the simulations corresponding (b), (d), and (f).

Figure 32. Density images and distributions of the (a) circular and (b) elliptic SF6 gas cylinder initially constructed.

88 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

Figure 36. Density profiles of nonuniform flows with Gaussian function and uniform flows in a vertical direction.

may be owing to the nonuniformity of flow field. Also, drawing to the experience on the diffusive interfacial transition layer of circular gas cylinder, we construct a nonuniform flow field with a Gaussian distribution of density along the direction perpendicular to the shock motion direction (as shown in Figure 36). The numerical result (shown in Figure 37) confirms this idea [12, 33, 34].

Figure 37. Simulated density distribution in a nonuniform flow field by using Gaussian function.

Figure 38. Experimental schlieren photography images and numerical simulation results by MVFT2D at a certain time series (the sizes of the pictures are ones of the test window [0.038, 0.25 m] [0.0, 0.2 m]).

Figure 38 shows the evolution of the interface and the propagation of the transmitted shock wave, the calculated results agree well with the experiment. Due to the nonuniform flow field of SF6 gas, the propagating velocity of transmitted shock wave in the upper part of shock tube is faster than in the bottom of shock tube, and it forms an oblique shock wave front and the interface. Figure 39 shows the shock front locations between experiment and numerical simulations, Figure 40 shows the perturbation amplitude of experiment, numerical simulation, and theories, they are in good agreement. Figure 41 shows the calculated perturbation amplitude of RM instability for different modes in the initial uniform and nonuniform flow fields. As we can see that, at late times, the growth of small perturbation in a low-density zone catches up and exceeds the large perturbation in a high-density zone, which is opposite to the case of uniform flow field. Figure 42 shows the perturbation amplitudes of four different initial amplitudes 2.5, 5.0, 7.5, and 10 mm in low- and high-density zones of nonuniform flow fields. It shows that the effect of initial amplitude on the growth of RM instability in the nonuniform flow field is different from the case of uniform one, which can be explained by the baroclinic vorticity shown in Figures 43 and 44, the baroclinic vorticity produced in the low-density zone is larger than that in the high-density zone.

Figure 39. Shock front locations of the experiment and calculated results on the three test lines.

Figure 40. Perturbation amplitudes of the experiment, simulation, and comparison with theories.

Figure 38 shows the evolution of the interface and the propagation of the transmitted shock wave, the calculated results agree well with the experiment. Due to the nonuniform flow field of SF6 gas, the propagating velocity of transmitted shock wave in the upper part of shock tube is faster than in the bottom of shock tube, and it forms an oblique shock wave front and the interface. Figure 39 shows the shock front locations between experiment and numerical simulations, Figure 40 shows the perturbation amplitude of experiment, numerical simulation, and theories, they are in good agreement. Figure 41 shows the calculated perturbation amplitude of RM instability for different modes in the initial uniform and nonuniform flow fields. As we can see that, at late times, the growth of small perturbation in a low-density zone catches up and exceeds the large perturbation in a high-density zone, which is opposite to the case of uniform flow field. Figure 42 shows the perturbation amplitudes of four different initial amplitudes 2.5, 5.0, 7.5, and 10 mm in low- and high-density zones of nonuniform flow fields. It shows that the effect of initial amplitude on the growth of RM instability in the nonuniform flow field is different from the case of uniform one, which can be explained by the baroclinic vorticity shown in Figures 43 and 44, the baroclinic vorticity produced in the low-density zone is larger than

Figure 38. Experimental schlieren photography images and numerical simulation results by MVFT2D at a certain time

Figure 37. Simulated density distribution in a nonuniform flow field by using Gaussian function.

90 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

series (the sizes of the pictures are ones of the test window [0.038, 0.25 m] [0.0, 0.2 m]).

that in the high-density zone.

Figure 41. Calculated perturbation amplitude of RM instability for different modes in the initial uniform and nonuniform flow fields.

Figure 42. Perturbation amplitudes of four different initial amplitude 2.5, 5.0, 7.5 and 10 mm in the low and high density zones of nonuniform flow field.

Figure 43. Average vorticity when initial amplitude is 5.0 mm which is in the low and high density zones of nonuniform flow field at 1 ms.

Figure 44. Baroclinic vorticity in the uniform and nonuniform flow fields with the initial amplitude group (A<sup>01</sup> ¼ 2.5 mm, A<sup>02</sup> ¼ 7.5 mm) at 1.0 ms.

Figure 42. Perturbation amplitudes of four different initial amplitude 2.5, 5.0, 7.5 and 10 mm in the low and high density

92 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

Figure 43. Average vorticity when initial amplitude is 5.0 mm which is in the low and high density zones of nonuniform

Figure 44. Baroclinic vorticity in the uniform and nonuniform flow fields with the initial amplitude group (A<sup>01</sup> ¼ 2.5 mm,

zones of nonuniform flow field.

A<sup>02</sup> ¼ 7.5 mm) at 1.0 ms.

flow field at 1 ms.

Figure 45. Density contour of the numerical simulated by MVFT at a time series. The left column with a uniform initial condition, the middle column with a δ<sup>1</sup> nonuniform Gaussian function, and the right column with a δ<sup>2</sup> nonuniform Gaussian function. The small arrow denotes the direction of propagation of the shock wave fronts before reshock the interface.

Figure 46. TMZ width of RM instability in the initial uniform and nonuniform flow fields.

For the above investigations of RM instability in the nonuniform flow field, they are all without reshock. In the following section, we will study the effect of nonuniformity of flow field and the reshock on the RM instability. Here, we consider two kinds of nonuniformity: the nonuniform coefficient δ<sup>1</sup> ¼ 0.6162 m and δ<sup>2</sup> ¼ 0.4961 m. Figure 45 shows the density contour of the numerical simulated by MVFT at a time series, the left column with a uniform initial conditions, the middle column with a δ<sup>1</sup> nonuniform Gaussian function, and the right column with a δ<sup>2</sup> nonuniform Gaussian function. There is a significant difference between the uniform and nonuniform flows before reshock, but the difference decreases in evidence after reshock. Figure 46 shows the TMZ width of RM instability in the initial uniform and nonuniform flow fields. It points out that the growth of the TMZ width for the initial nonuniform flow field is greater than that for the uniform flow field, and the lesser the nonuniform coefficient, the higher the growth rate of TMZ width, but the difference for the three different flow configurations diminishes after reshock.

## 5. Prospect

Now, we investigated the Richtmyer-Meshkov instability and the induced turbulent mixing in fluid flow by using large-eddy simulation, but the turbulent mixing is a complicated threedimensional problem with multiple time-space scales, and the more engineering applications involve the materials mixing with strength. So, the future work will be carried out in the following aspects: (a) the direct numerical simulations with high precision and high resolution on the platform of supercomputer, (b) the interface instability and turbulent materials mixing with strength, this may involve the material properties such as the deformation, fracture, melt, phase transition, material microstructures, etc., and it may make this problem to be more difficult.

## Acknowledgements

Professor Jingsong Bai and Associate professor Tao Wang from the Institute of Fluid Physics, China Academy of Engineering Physics, People's Republic of China and some other from the same institute should be thanked for the contributions and helps to these works, such as Prof. Ping Li, Dr. Liyong Zou, Kun Liu, Mr. Jinhong Liu, Min Zhong and Bing Wang, and Jiaxin Xiao.

The authors would like to thank the support from "Science Challenge Project" (no. TZ2016001), the National Natural Science Foundation of China (nos. 11372294, 11202195, 11072228, and 11532012), the Science Foundation of China Academy of Engineering Physics (nos. 2011B020 2005, 2011A0201002, and 2008B0202011), and the Foundation of National Key Laboratory of Shock Wave and Detonation Physics (no. 9140C670301150C67290).

## Author details

For the above investigations of RM instability in the nonuniform flow field, they are all without reshock. In the following section, we will study the effect of nonuniformity of flow field and the reshock on the RM instability. Here, we consider two kinds of nonuniformity: the nonuniform coefficient δ<sup>1</sup> ¼ 0.6162 m and δ<sup>2</sup> ¼ 0.4961 m. Figure 45 shows the density contour of the numerical simulated by MVFT at a time series, the left column with a uniform initial conditions, the middle column with a δ<sup>1</sup> nonuniform Gaussian function, and the right column with a δ<sup>2</sup> nonuniform Gaussian function. There is a significant difference between the uniform and nonuniform flows before reshock, but the difference decreases in evidence after reshock. Figure 46 shows the TMZ width of RM instability in the initial uniform and nonuniform flow fields. It points out that the growth of the TMZ width for the initial nonuniform flow field is greater than that for the uniform flow field, and the lesser the nonuniform coefficient, the higher the growth rate of TMZ width, but the difference for the three different flow configura-

94 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

Now, we investigated the Richtmyer-Meshkov instability and the induced turbulent mixing in fluid flow by using large-eddy simulation, but the turbulent mixing is a complicated threedimensional problem with multiple time-space scales, and the more engineering applications involve the materials mixing with strength. So, the future work will be carried out in the following aspects: (a) the direct numerical simulations with high precision and high resolution on the platform of supercomputer, (b) the interface instability and turbulent materials mixing with strength, this may involve the material properties such as the deformation, fracture, melt, phase transition, material microstructures, etc., and it may make this problem to be more

Professor Jingsong Bai and Associate professor Tao Wang from the Institute of Fluid Physics, China Academy of Engineering Physics, People's Republic of China and some other from the same institute should be thanked for the contributions and helps to these works, such as Prof. Ping Li, Dr. Liyong Zou, Kun Liu, Mr. Jinhong Liu, Min Zhong and Bing Wang, and Jiaxin Xiao. The authors would like to thank the support from "Science Challenge Project" (no. TZ2016001), the National Natural Science Foundation of China (nos. 11372294, 11202195, 11072228, and 11532012), the Science Foundation of China Academy of Engineering Physics (nos. 2011B020 2005, 2011A0201002, and 2008B0202011), and the Foundation of National Key Laboratory of

Shock Wave and Detonation Physics (no. 9140C670301150C67290).

tions diminishes after reshock.

5. Prospect

difficult.

Acknowledgements

Jingsong Bai1,2 and Tao Wang<sup>1</sup> \*

\*Address all correspondence to: wtao\_mg@163.com

1 Institute of Fluids Physics, China Academy of Engineering Physics, Mianyang, Sichuan, People's Republic of China

2 National Key Laboratory of Shock Wave and Detonation Physics (LSD), Institute of Fluid Physics, China Academy of Engineering Physics, Mianyang, Sichuan, People's Republic of China

## References


[24] Bates KR, Nikiforakis N, Holder D. Richtmyer-Meshkov instability induced by the interaction of a shock wave with a rectangular block of SF6. Physics of Fluids. 2007;19:036101

[10] Bai JS, Wang T, Li P, Zou LY, Liu CL. Numerical simulation of the hydrodynamic instability experiments and flow mixing. Science in China Series G: Physics, Mechanics &

[11] Wang T, Bai JS, Li P, Liu K. Large-eddy simulations of the Richtmyer-Meshkov instability of rectangular interface accelerated by shock waves. Science in China Series G: Physics,

[12] Bai JS, Liu JH, Wang T, Zou LY, Li P, Tan DW. Investigation of the Richtmyer-Meshkov instability with double perturbation interface in nonuniform flows. Physical Review E, Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics. 2010;81(2):

[13] Smagorinsky J. General circulation experiments with the primitive equations. Monthly

[14] Vreman AW. An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic

[15] Moin P, Squires K, Cabot W, Lee S. A dynamic subgrid-scale model for compressible turbulence and scalar transport. Physics of Fluids A (Fluid Dynamics). 1991;3(11):

[16] Hill DJ, Pantano C, Pullin DI. Large-eddy simulation and multiscale modelling of a Richtmyer-Meshkov instability with reshock. Journal of Fluid Mechanics. 2006;557:29–61

[17] Bai JS. High resolution numerical methods and adaptive mesh refinement algorithm for compressible multi-fluid dynamics [PhD thesis]. China Academy of Engineering Physics;

[18] Wang T, Bai JS, Li P, et al. Two and three dimensional numerical investigations of the single-mode Richtmyer-Meshkov instability. Chinese Journal of High Pressure Physics.

[19] Wang T, Bai JS, Li P, Zhong M. The numerical study of shock-induced hydrodynamic

[20] Malamud G, Elabz Y, Leinov E, et al. Two- and three-dimensional numerical simulations of the interaction of a Richtmyer-Meshkov instability induced turbulent mixing zone

[21] Bai JS, Wang T, Liu K, et al. Large-Eddy simulation of the Three-dimensional experiment on Richtmyer-Meshkov instability induced turbulence. International Journal of Astron-

[22] Poggi F, Thorembey MH, Rodriguez G. Velocity measurements in turbulent gaseous mixtures induced by Richtmyer–Meshkov instability. Physics of Fluids. 1998;10:2698–2712 [23] Wang T, Li P, Bai JS, Tao G, Wang B, Zou LY. Large-eddy simulation of the Richtmyer-

Meshkov instability. Canadian Journal of Physics. 2015;93(10):1124–1130

with several re-shocks. In: 10th IWPCTM; 17–21 July 2006, Paris, France. 2006.

theory and applications. Physics of Fluids. 2004;16(10):3670–3681

instability and mixing. Chinese Physics B. 2009;18(3):1127–1135

Astronomy. 2009;52(12):2027–2040

Weather Review. 1963;91:99–164

056302

2746–2757

2003. (in Chinese)

2013;27(2):277–286

omy and Astrophysics. 2012;2:28–36

Mechanics & Astronomy. 2010;53(5):905–914

96 Turbulence Modelling Approaches - Current State, Development Prospects, Applications

