**4.1. The Oberbeck-Boussinesq approximation of Rayleigh-Benard convection, dimensionless form and benchmarks**

One of the possible mathematical models for this problem is the Oberbeck-Boussinesq approximation. Here we are closely following [51] and assume, that fluid physical properties ( *ν*,*β* ) are only linear functions of temperature perturbations. The fluid density can be given

#### 20 Will-be-set-by-IN-TECH 270 Nonlinearity, Bifurcation and Chaos – Theory and Applications FSM Scenarios of Laminar-Turbulent Transition in Incompressible Fluids <sup>21</sup>

as a function of temperature perturbation as:

$$
\rho = \rho\_0 \left( 1 - \beta (T - T\_0) \right),
\tag{41}
$$

solution routine is adopted for CPU calculations [53] and Geometric Multigrid method for

In order to satisfy conditions of the theorems we must adopt mesh by introducing Reynolds number through Prandtl and Payleigh numbers. We are using paper results [55, 56] that

for the range of 0.9 <sup>&</sup>lt; *Pr* <sup>&</sup>lt; 2 and 1 · 105 <sup>&</sup>lt; *Ra* <sup>&</sup>lt; <sup>1</sup> · 109. It is also known [57, 58] that transition form "soft" turbulence (where some frequencies can be determined in frequency analysis) to "hard" turbulence (where frequency-amplitude response becomes a coloured noise) occurs at *Racr* <sup>≈</sup> <sup>4</sup> · 107 for *Pr* <sup>&</sup>gt; 0.9. Assuming that it is impossible to make any quantitative analysis after similarity criteria greater than critical values we take relation (44) and find Reynolds

Using (45) and applying it to (27) one can get the following mesh adopted variants for cubic

*G*4 = {250*X*250*X*100} for *Pr* > 4.0; Δ*xmax* = 4.5Δ*xwall*;*r*/*h* = 4; *G*5 = {400*X*400*X*47} for *Pr* > 5.0; Δ*xmax* = 3.5Δ*xwall*;*r*/*h* = 14;

where *r* is the cylinder radius and *h* is the cylinder height. Immersed boundary is used for cylindrical approximation, for more information see [59]. For other geometry domains grid is

Benchmarking the problem requires many different domain configurations and different dimensionless variations of the equations (43). We skip that for the sake of brevity and can recommend book by professor Getling A.V. [51] for more information on dimensionless forms, physical background and analytical analysis. It is known from the linear minimodal approximation of the problem [51] that the flow with wall boundary conditions on Ω<sup>0</sup> is more stable than that with periodic conditions or two dimensional problems. For the infinite *Pr* number the critical value of Rayleigh criterion for the first instability is *Racr* = 1707.762 for the wave number *<sup>k</sup>* <sup>=</sup> 3.117 in case of wall boundary conditions and *Racr* <sup>=</sup> <sup>27</sup>*π*4/4 <sup>≈</sup> 657.511 with *k* = 2.221 for periodic boundary conditions. For benchmark verification several papers

Solution in cylindric domain for *Ra* = 2000 − 31000 and *r*/*h* = 4 with zero initial conditions is shown in fig.14, top and with initial perturbation *Vz* = cos (2*πz*/(2*r*)) /100 – in fig.14 bottom, where one can see the appearance of absolutely different solution for perturbed initial

Some results in rectangular domain where also considered and compared with [61] with very good agreement. More results where compared with [51, 60] for *r*/*h* = 14 for range of

*G*1 = {250*X*250*X*250} for *Pr* = 0.9; Δ*xmax* = 6.1Δ*xwall*; *G*2 = {216*X*216*X*216} for *Pr* > 1.0; Δ*xmax* = 5.2Δ*xwall*; *G*3 = {185*X*185*X*185} for *Pr* > 1.25; Δ*xmax* = 5.0Δ*xwall*;

*Rmax* � *Ra*0.44

specified in the same manner and omitted here for the sake of brevity.

conditions with the prototype function cos(*kx*), see [51].

and cylindrical domains after calibration simulations:

*<sup>R</sup>* � *Ra*0.44*Pr*<sup>−</sup>0.76, (44)

FSM Scenarios of Laminar-Turbulent Transition in Incompressible Fluids 271

*cr* min(*Prcr*)−0.76 = 2856. (45)

(46)

GPU [54].

number as:

where considered [60].

indicate the following relation is true:

where *ρ* - fluid density, *T* - fluid temperature, *β* - fluid thermal expansion coefficient, *ρ*0,*T*<sup>0</sup> - mean values of fluid density and temperature. It is assumed under Oberbeck-Boussinesq approximation that density only changes due to temperature difference and, thus, causing buoyancy, yet fluid is considered incompressible. Temperature emission due to friction is also neglected. Introducing (41) to Navier-Stokes equations (1), assuming temperature passive advection-diffusion and taking gravity vector in Cartesian coordinates as *g* = {0; 0; −1} in account, one gets:

$$\begin{array}{c} \nabla \cdot \mathbf{V} = 0, \\ \frac{\partial \mathbf{V}}{\partial t} + (\mathbf{V} \cdot \nabla) \mathbf{V} + \rho\_0^{-1} \nabla P = \nu \nabla^2 \mathbf{V} + \vec{g}\beta (T - T\_0), \end{array} \tag{42}$$
 
$$\begin{array}{c} \frac{\partial T}{\partial t} + \mathbf{V} \nabla T = \chi \nabla^2 T, \end{array} \tag{42}$$

here *χ* is a fluid thermal conductivity coefficient. There are many scales can be chosen that make (42) dimensionless for the Rayleigh-Benard convection problem. One of the most common ways [51] is to use time scale *τ* as *τ* = *h*2/*ν*, where *h* is the length between two planes with given temperature difference. Another way is to associate *τ* with the momentum transport by viscous terms and then *τ* = *h*2/*ν*. Introducing additional dimensionless similarity criteria one can formulate a Rayleigh-Benard convection problem: find vector-function **<sup>V</sup>** : <sup>Ω</sup> <sup>×</sup> [0, *<sup>t</sup>*] <sup>→</sup> **<sup>R</sup>**3, scalar pressure function *<sup>P</sup>* : <sup>Ω</sup> <sup>×</sup> [0, *<sup>t</sup>*] <sup>→</sup> **<sup>R</sup>** and and scalar temperature function *T* : Ω × [0, *t*] → **R** that satisfy the following initial-boundary value problem:

$$\begin{aligned} \frac{\partial \mathbf{V}}{\partial t} + \left(\mathbf{V} \cdot \nabla\right) \mathbf{V} + \nabla P &= \nabla^2 \mathbf{V} + \mathbf{R} a Pr^{-1} (T - T\_0) \cdot (0; 0; -1)^T \text{ in } Q; \\ \nabla \cdot \mathbf{V} &= 0 \text{ in } Q = \Omega \times (0, t); \\ \frac{\partial T}{\partial t} + \mathbf{V} \nabla T &= Pr^{-1} \nabla^2 T \text{ in } Q; \\ \mathbf{V} &= 0, \partial T / \partial \vec{\mathbf{n}} = 0, \text{ on } \partial \Omega\_0 \times (0, t); \mathbf{V} = 0, T = T\_d \text{ on } \partial \Omega\_1; \\ \mathbf{V}(\vec{x}, 0) &= \mathbf{V}\_0(\vec{x}), \nabla \cdot \mathbf{V}\_0 = 0; T(\vec{x}, 0) = T\_0(\vec{x}) \text{ in } \Omega. \end{aligned} \tag{43}$$

Here: Ω is a bounded domain with local Lipschiz-continuous boundary *∂*Ω*i*; *t* is time; *T* is a fluid temperature; *T*<sup>0</sup> is a reference fluid temperature; *Pr* = *ν*/*χ* is the dimensionless Prandtl number; *Ra* = *gβh*3Δ*T*/(*νχ*) is the Rayleigh number; *h* is height in Ω between *∂*Ω<sup>1</sup> ; Δ*T* is a temperature difference between *∂*Ω1.

There are two main types of boundary conditions. First consider boundary conditions whose plane is parallel to the temperature gradient, i.e. Ω0. One usually chooses either periodic boundary conditions for temperature and velocity or wall boundary conditions with Neumann type for temperature. On other planes, namely Ω1, temperature gradient by Dirichlet boundary conditions is set with wall no slip boundary for velocity.

Numerical solution method differs from (28) only in temperature equation and in diffusion part of Navier-Stokes equations. We applied implicit five diagonal matrix solution method [52] for all diffusion parts of (43) thus accuracy drops from 6-th order to 4-th order in space. It can be shown by the modified wave number analysis that number of elements remains the same just with a little lost of accuracy. In order to solve matrix equation [*A*] [*X*] = [*B*] that arises from implicit method for diffusion parts of equation a five diagonal fast factorization solution routine is adopted for CPU calculations [53] and Geometric Multigrid method for GPU [54].

20 Will-be-set-by-IN-TECH

where *ρ* - fluid density, *T* - fluid temperature, *β* - fluid thermal expansion coefficient, *ρ*0,*T*<sup>0</sup> - mean values of fluid density and temperature. It is assumed under Oberbeck-Boussinesq approximation that density only changes due to temperature difference and, thus, causing buoyancy, yet fluid is considered incompressible. Temperature emission due to friction is also neglected. Introducing (41) to Navier-Stokes equations (1), assuming temperature passive advection-diffusion and taking gravity vector in Cartesian coordinates as *g* = {0; 0; −1} in

∇ · **V** = 0,

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> **<sup>V</sup>**∇*<sup>T</sup>* <sup>=</sup> *<sup>χ</sup>*∇2*T*,

here *χ* is a fluid thermal conductivity coefficient. There are many scales can be chosen that make (42) dimensionless for the Rayleigh-Benard convection problem. One of the most common ways [51] is to use time scale *τ* as *τ* = *h*2/*ν*, where *h* is the length between two planes with given temperature difference. Another way is to associate *τ* with the momentum transport by viscous terms and then *τ* = *h*2/*ν*. Introducing additional dimensionless similarity criteria one can formulate a Rayleigh-Benard convection problem: find vector-function **<sup>V</sup>** : <sup>Ω</sup> <sup>×</sup> [0, *<sup>t</sup>*] <sup>→</sup> **<sup>R</sup>**3, scalar pressure function *<sup>P</sup>* : <sup>Ω</sup> <sup>×</sup> [0, *<sup>t</sup>*] <sup>→</sup> **<sup>R</sup>** and and scalar temperature function *T* : Ω × [0, *t*] → **R** that satisfy the following initial-boundary

> *<sup>∂</sup><sup>t</sup>* <sup>+</sup> (**<sup>V</sup>** · ∇) **<sup>V</sup>** <sup>+</sup> <sup>∇</sup>*<sup>P</sup>* <sup>=</sup> <sup>∇</sup>2**<sup>V</sup>** <sup>+</sup> *RaPr*−1(*<sup>T</sup>* <sup>−</sup> *<sup>T</sup>*0) · (0; 0; <sup>−</sup>1)*<sup>T</sup>* in *<sup>Q</sup>*; ∇ · **V** = 0 in *Q* = Ω × (0, *t*);

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> **<sup>V</sup>**∇*<sup>T</sup>* <sup>=</sup> *Pr*−1∇2*<sup>T</sup>* in *<sup>Q</sup>*; **V** = 0, *∂T*/*∂n* = 0, on *∂*Ω<sup>0</sup> × (0, *t*); **V** = 0, *T* = *T<sup>α</sup>* on *∂*Ω1; **V**(*x*, 0) = **V**0(*x*), ∇ · **V**<sup>0</sup> = 0; *T*(*x*, 0) = *T*0(*x*) in Ω.

Here: Ω is a bounded domain with local Lipschiz-continuous boundary *∂*Ω*i*; *t* is time; *T* is a fluid temperature; *T*<sup>0</sup> is a reference fluid temperature; *Pr* = *ν*/*χ* is the dimensionless Prandtl number; *Ra* = *gβh*3Δ*T*/(*νχ*) is the Rayleigh number; *h* is height in Ω between *∂*Ω<sup>1</sup> ; Δ*T* is a

There are two main types of boundary conditions. First consider boundary conditions whose plane is parallel to the temperature gradient, i.e. Ω0. One usually chooses either periodic boundary conditions for temperature and velocity or wall boundary conditions with Neumann type for temperature. On other planes, namely Ω1, temperature gradient by

Numerical solution method differs from (28) only in temperature equation and in diffusion part of Navier-Stokes equations. We applied implicit five diagonal matrix solution method [52] for all diffusion parts of (43) thus accuracy drops from 6-th order to 4-th order in space. It can be shown by the modified wave number analysis that number of elements remains the same just with a little lost of accuracy. In order to solve matrix equation [*A*] [*X*] = [*B*] that arises from implicit method for diffusion parts of equation a five diagonal fast factorization

Dirichlet boundary conditions is set with wall no slip boundary for velocity.

<sup>0</sup> <sup>∇</sup>*<sup>P</sup>* <sup>=</sup> *<sup>ν</sup>*∇2**<sup>V</sup>** <sup>+</sup>*<sup>g</sup>β*(*<sup>T</sup>* <sup>−</sup> *<sup>T</sup>*0), *<sup>∂</sup><sup>T</sup>*

*ρ* = *ρ*<sup>0</sup> (1 − *β*(*T* − *T*0)), (41)

(42)

(43)

as a function of temperature perturbation as:

*∂***V**

*<sup>∂</sup><sup>t</sup>* + (**<sup>V</sup>** · ∇)**<sup>V</sup>** <sup>+</sup> *<sup>ρ</sup>*−<sup>1</sup>

*∂T*

account, one gets:

value problem:

*∂***V**

temperature difference between *∂*Ω1.

In order to satisfy conditions of the theorems we must adopt mesh by introducing Reynolds number through Prandtl and Payleigh numbers. We are using paper results [55, 56] that indicate the following relation is true:

$$R \simeq R a^{0.44} Pr^{-0.76} \text{ \AA} \tag{44}$$

for the range of 0.9 <sup>&</sup>lt; *Pr* <sup>&</sup>lt; 2 and 1 · 105 <sup>&</sup>lt; *Ra* <sup>&</sup>lt; <sup>1</sup> · 109. It is also known [57, 58] that transition form "soft" turbulence (where some frequencies can be determined in frequency analysis) to "hard" turbulence (where frequency-amplitude response becomes a coloured noise) occurs at *Racr* <sup>≈</sup> <sup>4</sup> · 107 for *Pr* <sup>&</sup>gt; 0.9. Assuming that it is impossible to make any quantitative analysis after similarity criteria greater than critical values we take relation (44) and find Reynolds number as:

$$R\_{\text{max}} \simeq R a\_{cr}^{0.44} \min(Pr\_{cr})^{-0.76} = 2856. \tag{45}$$

Using (45) and applying it to (27) one can get the following mesh adopted variants for cubic and cylindrical domains after calibration simulations:

$$\begin{aligned} \text{G1} &= \{250X250X250\} \text{ for } Pr = 0.9; \Delta x\_{\text{max}} = 6.1 \Delta x\_{\text{wall}};\\ \text{G2} &= \{216X216X216\} \text{ for } Pr > 1.0; \Delta x\_{\text{max}} = 5.2 \Delta x\_{\text{wall}};\\ \text{G3} &= \{185X185X185\} \text{ for } Pr > 1.25; \Delta x\_{\text{max}} = 5.0 \Delta x\_{\text{wall}};\\ \text{G4} &= \{250X250X100\} \text{ for } Pr > 4.0; \Delta x\_{\text{max}} = 4.5 \Delta x\_{\text{wall}}; r/h = 4;\\ \text{G5} &= \{400X400X47\} \text{ for } Pr > 5.0; \Delta x\_{\text{max}} = 3.5 \Delta x\_{\text{wall}}; r/h = 14; \end{aligned}$$

where *r* is the cylinder radius and *h* is the cylinder height. Immersed boundary is used for cylindrical approximation, for more information see [59]. For other geometry domains grid is specified in the same manner and omitted here for the sake of brevity.

Benchmarking the problem requires many different domain configurations and different dimensionless variations of the equations (43). We skip that for the sake of brevity and can recommend book by professor Getling A.V. [51] for more information on dimensionless forms, physical background and analytical analysis. It is known from the linear minimodal approximation of the problem [51] that the flow with wall boundary conditions on Ω<sup>0</sup> is more stable than that with periodic conditions or two dimensional problems. For the infinite *Pr* number the critical value of Rayleigh criterion for the first instability is *Racr* = 1707.762 for the wave number *<sup>k</sup>* <sup>=</sup> 3.117 in case of wall boundary conditions and *Racr* <sup>=</sup> <sup>27</sup>*π*4/4 <sup>≈</sup> 657.511 with *k* = 2.221 for periodic boundary conditions. For benchmark verification several papers where considered [60].

Solution in cylindric domain for *Ra* = 2000 − 31000 and *r*/*h* = 4 with zero initial conditions is shown in fig.14, top and with initial perturbation *Vz* = cos (2*πz*/(2*r*)) /100 – in fig.14 bottom, where one can see the appearance of absolutely different solution for perturbed initial conditions with the prototype function cos(*kx*), see [51].

Some results in rectangular domain where also considered and compared with [61] with very good agreement. More results where compared with [51, 60] for *r*/*h* = 14 for range of

22 Will-be-set-by-IN-TECH 272 Nonlinearity, Bifurcation and Chaos – Theory and Applications FSM Scenarios of Laminar-Turbulent Transition in Incompressible Fluids <sup>23</sup>

**4.2. Nonlinear dynamics of laminar-turbulent transition.**

subspace. This cycle projection for point *p*<sup>1</sup> is shown in fig16.

calculations.

*Ra* <sup>=</sup> 2.67820 · 105.

*4.2.1. First calculation series.*

The analysis is absolutely identical to the one used for the problem of the flow over backward facing step. For this point we select five points with relative Cartesian coordinates *p*<sup>1</sup> = {0.5; 0.5; 0.5}; *p*<sup>2</sup> = {0.4; 0.5; 0.5}; *p*<sup>3</sup> = {0.5; 0.8; 0.5}; *p*<sup>4</sup> = {0.8; 0.5; 0.667}; *p*<sup>5</sup> = {0.16; 0.167; 0.1}. The subspaces of infinite dimensional phase space are constructed by velocity vector functions in various points and their combination. Similarity criteria of *Ra* and *Pr* numbers are used as bifurcation parameters. We took fixed values for Prandtl number and increased *Ra* number for every *Pr* number, thus we considered five various series of

FSM Scenarios of Laminar-Turbulent Transition in Incompressible Fluids 273

Prandtl number is set as 1.866666666. While we change *Ra* from 1 to 14100 we could see the formation of various stationary solutions with various recirculation zones formation in the domain with the complete correspondence to [62]. Laminar solution can be seen for *Ra* < 2.5 · 105 in the whole domain, that corresponded to the fixed stationary point in the infinite dimensional phase space and in all subspaces. However this fixed point jumped from one position to another in the phase space as a function of *Ra* number. This stationary point looses stability for *Ra* <sup>=</sup> 2.5 · 105 and a limited cycle is formed in the phase space and in each phase

**Figure 16.** Phase 3D subspace and its sections for *Pr* = 1.86666. From left to right: limited cycle, point

This cycle looses stability at *Ra* <sup>=</sup> 2.67820 · 105 and forms two dimensional invariant torus through the Andronov-Hopf bifurcation. This torus projection in the three dimensional phase subspace and its plane section are depicted in fig. 16. This torus is lying in the whole infinite

The two-dimensional tori of double and quadruple periods are formed with the further increasing of Rayleigh number (fig. 17 and fig.18 ). Point *p*<sup>5</sup> in fig.18 depicts the torus

Further increasing of Ra number led to chaotic instability. This can be related to the limited accuracy of the numerical method and solution trajectories slipped from one to another in such sensible dynamic system. The cascade of bifurcations agreed well with experimental

*<sup>p</sup>*1, *Ra* <sup>=</sup> 2.5 · 105; two dimensional invariant torus, point *<sup>p</sup>*1; Poincare section of the torus,

dimensional phase space and has projections in all of its subspaces.

projection near wall boundaries for *Ra* <sup>=</sup> 2.9133225210 · 105.

**Figure 14.** Velocity vector function section *Ra* = 2000 − 31000, *r*/*h* = 4. Zero initial conditions on the top and perturbed initial conditions on the bottom.

*Ra* = {2565, 4500, 6500, 16850}. We could clearly see the process of different pattern formation (round, circular, "vetruvian man", hexagon, irregular). For the analysis of nonlinear dynamics we considered a cubic domain absolutely identical to the physical experiment, presented in the book by P.G. Frik [62]. Experiment was made in a copper cubic domain with the height of 40mm. Horizontal boundaries where thermally stabilized and vertical boundaries formed a constant gradient. The experiment was aimed to investigate frequencies and attractors using thermal differential pares. Numerical simulation for *Ra* <sup>=</sup> <sup>2</sup> · 105 and *Pr* <sup>=</sup> 7 is shown in fig.15.

**Figure 15.** Velocity vectors in 3D and plane section (1;0;0), for *Ra* <sup>=</sup> <sup>2</sup> · 105 and *Pr* <sup>=</sup> 7. Every six vector is shown.

## **4.2. Nonlinear dynamics of laminar-turbulent transition.**

The analysis is absolutely identical to the one used for the problem of the flow over backward facing step. For this point we select five points with relative Cartesian coordinates *p*<sup>1</sup> = {0.5; 0.5; 0.5}; *p*<sup>2</sup> = {0.4; 0.5; 0.5}; *p*<sup>3</sup> = {0.5; 0.8; 0.5}; *p*<sup>4</sup> = {0.8; 0.5; 0.667}; *p*<sup>5</sup> = {0.16; 0.167; 0.1}. The subspaces of infinite dimensional phase space are constructed by velocity vector functions in various points and their combination. Similarity criteria of *Ra* and *Pr* numbers are used as bifurcation parameters. We took fixed values for Prandtl number and increased *Ra* number for every *Pr* number, thus we considered five various series of calculations.

### *4.2.1. First calculation series.*

22 Will-be-set-by-IN-TECH

**Figure 14.** Velocity vector function section *Ra* = 2000 − 31000, *r*/*h* = 4. Zero initial conditions on the

*Ra* = {2565, 4500, 6500, 16850}. We could clearly see the process of different pattern formation (round, circular, "vetruvian man", hexagon, irregular). For the analysis of nonlinear dynamics we considered a cubic domain absolutely identical to the physical experiment, presented in the book by P.G. Frik [62]. Experiment was made in a copper cubic domain with the height of 40mm. Horizontal boundaries where thermally stabilized and vertical boundaries formed a constant gradient. The experiment was aimed to investigate frequencies and attractors using thermal differential pares. Numerical simulation for *Ra* <sup>=</sup> <sup>2</sup> · 105 and *Pr* <sup>=</sup> 7 is shown in

**Figure 15.** Velocity vectors in 3D and plane section (1;0;0), for *Ra* <sup>=</sup> <sup>2</sup> · 105 and *Pr* <sup>=</sup> 7. Every six vector

top and perturbed initial conditions on the bottom.

fig.15.

is shown.

Prandtl number is set as 1.866666666. While we change *Ra* from 1 to 14100 we could see the formation of various stationary solutions with various recirculation zones formation in the domain with the complete correspondence to [62]. Laminar solution can be seen for *Ra* < 2.5 · 105 in the whole domain, that corresponded to the fixed stationary point in the infinite dimensional phase space and in all subspaces. However this fixed point jumped from one position to another in the phase space as a function of *Ra* number. This stationary point looses stability for *Ra* <sup>=</sup> 2.5 · 105 and a limited cycle is formed in the phase space and in each phase subspace. This cycle projection for point *p*<sup>1</sup> is shown in fig16.

**Figure 16.** Phase 3D subspace and its sections for *Pr* = 1.86666. From left to right: limited cycle, point *<sup>p</sup>*1, *Ra* <sup>=</sup> 2.5 · 105; two dimensional invariant torus, point *<sup>p</sup>*1; Poincare section of the torus, *Ra* <sup>=</sup> 2.67820 · 105.

This cycle looses stability at *Ra* <sup>=</sup> 2.67820 · 105 and forms two dimensional invariant torus through the Andronov-Hopf bifurcation. This torus projection in the three dimensional phase subspace and its plane section are depicted in fig. 16. This torus is lying in the whole infinite dimensional phase space and has projections in all of its subspaces.

The two-dimensional tori of double and quadruple periods are formed with the further increasing of Rayleigh number (fig. 17 and fig.18 ). Point *p*<sup>5</sup> in fig.18 depicts the torus projection near wall boundaries for *Ra* <sup>=</sup> 2.9133225210 · 105.

Further increasing of Ra number led to chaotic instability. This can be related to the limited accuracy of the numerical method and solution trajectories slipped from one to another in such sensible dynamic system. The cascade of bifurcations agreed well with experimental

**Figure 17.** Phase 3D subspace and its sections for *Pr* = 1.86666. From left to right: double period two-dimensional invariant torus, point *<sup>p</sup>*1, *Ra* <sup>=</sup> 2.767858170 · 105; Poincare section of the double period torus; quadruple period two-dimensional invariant torus, point *<sup>p</sup>*1, *Ra* <sup>=</sup> 2.9133225210 · 105.

**Figure 19.** Phase projection and its sections for *Pr* = 1.6129. From left to right: the two dimensional invariant torus projection on (*Vx*, *Vy*) plane for point *<sup>p</sup>*1, *Ra* <sup>=</sup> 1.366 · 106; Poincare section of the torus at

double period two-dimensional invariant torus that can be seen in fig.19. Further analysis was

FSM Scenarios of Laminar-Turbulent Transition in Incompressible Fluids 275

Prundtl number is set 1.354839. The limit cycle in system phase space is formed from the stationary point at *Ra* <sup>=</sup> 1.286 · 106 and its projection on (*Vx*, *Vy*) plane is shown in fig.20.

**Figure 20.** Phase subspace projection to (*Vx*, *Vy*) plane for *Pr* = 1.354839 at point *p*1. From left to right: stable cycle, *Ra* <sup>=</sup> 1.286 · 106; double period stable cycle, *Ra* <sup>=</sup> 1.296 · 106; quintuple period stable cycle,

Next doubling period bifurcation occurs at *Ra* <sup>=</sup> 1.296 · 106 with the formation of double period limit cycle, see fig.20. This cycle is lying in all infinite dimensional phase space and in all subspaces. With the increasing of *Ra* number the cycle suffers cascades of bifurcations in accordance with the FSM scenario. The quintuple period cycle is formed for *Ra* <sup>=</sup> 1.306 · 106

A triple period cycle can be seen at *Ra* <sup>=</sup> 1.308 · 106 in all subspaces and its projection on (*Vx*, *Vy*) plane is shown in fig.21. It means that the system suffered the whole subharmonic cascade of bifurcations and now there are all other unstable cycles exist in the system in accordance with the theory FSM. In conclusion one can say that there exist multiple scenarios of laminar-turbulent transition in Rayleigh-Benard convection as functions of *Pr*/*Ra* ratio:

2. Landau-Hopf bifurcation scenario with doubling period bifurcation on n-dimensional

1. Landau-Hopf bifurcation scenario that forms sets of many-dimensional tori *Tn*.

point *<sup>p</sup>*1; Poincare section of the double period torus at point *<sup>p</sup>*1, *Ra* <sup>=</sup> 1.36905 · 106.

limited due to the numerical noise.

who's projection on (*Vx*, *Vy*) plane is shown in fig.20.

invariant torus (Landau-Hopf scenario + FSM scenario).

*4.2.3. Third calculation series.*

*Ra* <sup>=</sup> 1.306 · 106.

**Figure 18.** Phase 3D subspace and its sections for *Pr* = 1.86666. From left to right: 1/2 (zoomed) Poincare section of the quadruple period two-dimensional torus, point *<sup>p</sup>*1, *Ra* <sup>=</sup> 2.9133225210 · 105; quadruple period two-dimensional invariant torus, point *<sup>p</sup>*5, *Ra* <sup>=</sup> 2.9133225210 · 105; Poincare section of the quadruple period torus, point *p*1.

results from [62] with only difference in exact values for *Ra* number. In [62] it could be clearly seen the formation of limit cycle and torus (two irrational frequencies), but further investigations led to chaos so it was impossible to tell from the experimental data of further bifurcations for the formed torus.

## *4.2.2. Second calculation series.*

*Pr* number is fixed as 1.61290. The fluid motion is stationary up to *Ra* <sup>=</sup> 1.361 · 106 when a limit cycle is formed from the stationary point in the whole infinite dimension phase space. Further increasing of *Ra* up to 1.365 · 106 led to the doubling period bifurcation with formation of double period limit cycle and, immediately, formation of the two-dimensional invariant torus through the Andronov-Hopf bifurcation.

This torus can be clearly seen for *Ra* <sup>=</sup> 1.366 · 106 on the 3D subspace projection to the (*Vx*, *Vy*) plane in fig.19. Another doubling period bifurcation occurs at *Ra* <sup>=</sup> 1.36905 · 106 forming a

**Figure 19.** Phase projection and its sections for *Pr* = 1.6129. From left to right: the two dimensional invariant torus projection on (*Vx*, *Vy*) plane for point *<sup>p</sup>*1, *Ra* <sup>=</sup> 1.366 · 106; Poincare section of the torus at point *<sup>p</sup>*1; Poincare section of the double period torus at point *<sup>p</sup>*1, *Ra* <sup>=</sup> 1.36905 · 106.

double period two-dimensional invariant torus that can be seen in fig.19. Further analysis was limited due to the numerical noise.

### *4.2.3. Third calculation series.*

24 Will-be-set-by-IN-TECH

**Figure 17.** Phase 3D subspace and its sections for *Pr* = 1.86666. From left to right: double period two-dimensional invariant torus, point *<sup>p</sup>*1, *Ra* <sup>=</sup> 2.767858170 · 105; Poincare section of the double period

**Figure 18.** Phase 3D subspace and its sections for *Pr* = 1.86666. From left to right: 1/2 (zoomed) Poincare section of the quadruple period two-dimensional torus, point *<sup>p</sup>*1, *Ra* <sup>=</sup> 2.9133225210 · 105; quadruple period two-dimensional invariant torus, point *<sup>p</sup>*5, *Ra* <sup>=</sup> 2.9133225210 · 105; Poincare section of

results from [62] with only difference in exact values for *Ra* number. In [62] it could be clearly seen the formation of limit cycle and torus (two irrational frequencies), but further investigations led to chaos so it was impossible to tell from the experimental data of further

*Pr* number is fixed as 1.61290. The fluid motion is stationary up to *Ra* <sup>=</sup> 1.361 · 106 when a limit cycle is formed from the stationary point in the whole infinite dimension phase space. Further increasing of *Ra* up to 1.365 · 106 led to the doubling period bifurcation with formation of double period limit cycle and, immediately, formation of the two-dimensional invariant

This torus can be clearly seen for *Ra* <sup>=</sup> 1.366 · 106 on the 3D subspace projection to the (*Vx*, *Vy*) plane in fig.19. Another doubling period bifurcation occurs at *Ra* <sup>=</sup> 1.36905 · 106 forming a

the quadruple period torus, point *p*1.

bifurcations for the formed torus.

torus through the Andronov-Hopf bifurcation.

*4.2.2. Second calculation series.*

torus; quadruple period two-dimensional invariant torus, point *<sup>p</sup>*1, *Ra* <sup>=</sup> 2.9133225210 · 105.

Prundtl number is set 1.354839. The limit cycle in system phase space is formed from the stationary point at *Ra* <sup>=</sup> 1.286 · 106 and its projection on (*Vx*, *Vy*) plane is shown in fig.20.

**Figure 20.** Phase subspace projection to (*Vx*, *Vy*) plane for *Pr* = 1.354839 at point *p*1. From left to right: stable cycle, *Ra* <sup>=</sup> 1.286 · 106; double period stable cycle, *Ra* <sup>=</sup> 1.296 · 106; quintuple period stable cycle, *Ra* <sup>=</sup> 1.306 · 106.

Next doubling period bifurcation occurs at *Ra* <sup>=</sup> 1.296 · 106 with the formation of double period limit cycle, see fig.20. This cycle is lying in all infinite dimensional phase space and in all subspaces. With the increasing of *Ra* number the cycle suffers cascades of bifurcations in accordance with the FSM scenario. The quintuple period cycle is formed for *Ra* <sup>=</sup> 1.306 · 106 who's projection on (*Vx*, *Vy*) plane is shown in fig.20.

A triple period cycle can be seen at *Ra* <sup>=</sup> 1.308 · 106 in all subspaces and its projection on (*Vx*, *Vy*) plane is shown in fig.21. It means that the system suffered the whole subharmonic cascade of bifurcations and now there are all other unstable cycles exist in the system in accordance with the theory FSM. In conclusion one can say that there exist multiple scenarios of laminar-turbulent transition in Rayleigh-Benard convection as functions of *Pr*/*Ra* ratio:

1. Landau-Hopf bifurcation scenario that forms sets of many-dimensional tori *Tn*.

2. Landau-Hopf bifurcation scenario with doubling period bifurcation on n-dimensional invariant torus (Landau-Hopf scenario + FSM scenario).

**6. References**

416-33.

World Scientific.

, 41, 11, pp. 1550-58.

3, pp. 318-322.

28, pp714-725.


Equations, Academic Press, NY, 55-73

Methods.// Proc. R. Soc. Lond.v.434, pp23-39.

Cambridge: Clay Mathematics Institute.


systems. //Differential Equations, 44, 12, pp. 1618-27.

Commun. Nonlinear Sci. Numer. Simul., 15, p. 2851-2859.

equations.//Mat. Sb., Volume 188, Number 6, Pages 47-56.

[1] Carlson A., Jaffe A. and Wiles A. (Editors). (2006) The Millennium Prize Problems.

FSM Scenarios of Laminar-Turbulent Transition in Incompressible Fluids 277

[2] Sadovnichii V. and Simo S.(Editors). (2002) Modern Problems of Chaos and Nonlinearity.

[3] Magnitskii N.A. (2008) Universal theory of dynamical chaos in nonlinear dissipative systems of differential equations. // Commun. Nonlinear Sci. Numer. Simul., 13, pp.

[4] Magnitskii N.A. (2007) Universal theory of dynamical and spatio-temporal chaos in complex systems. // Dynamics of Complex Systems 1, 1, pp. 18-39. (in Russian) [5] Magnitskii N.A. (2008) New approach to analysis of Hamiltonian and conservative

[6] Magnitskii N.A. (2010) On topological structure of singular attractors of nonlinear systems of differential equations. //Differential Equations, 46, 11, pp.1551-1560. [7] Magnitskii N.A., Sidorov S.V. (2006) New Methods for Chaotic Dynamics. - Singapure:

[8] Magnitskii N.A., Sidorov S.V. (2005) On transition to diffusion chaos through a subharmonic cascade of bifurcations of two-dimensional tori. // Differential Equations

[9] Evstigneev N.M., Magnitskii N.A., Sidorov S.V. (2009) On nature of laminar-turbulent flow in backward facing step problem. //Differential equations, v.45, 1, pp.69-73. [10] Evstigneev N.M., Magnitskii N.A., Sidorov S.V. (2009) On the nature of turbulence in Rayleigh-Benard convection.//Differential Equations V.45, N.6, pp.909-912. [11] Evstigneev N.M, Magnitskii N.A. and Sidorov S.V. (2010) Nonlinear dynamics of laminarturbulent transition in three dimensional Rayleigh-Benard convection.//

[12] Evstigneev N.M., Magnitskii N.A. (2010) On possible scenarios of the transition to turbulence in Rayleigh-Benard convection.// Doklady Akademii Nauk, Vol. 433, No.

[13] Filatov A.N., Ipatova V.M. (1996) On globally stable difference schemes for the barotropic vorticity equation on a sphere // Russian J. Numer. Anal. Math. Modelling.

[14] Ipatova V.M. (1997) Attractors of approximations to non-autonomous evolution

[15] Foias C., Temam R.(1987) The connection between the Navier-Stokes equations, dynamical systems, and turbulence theory, in Directions in Partial Differential

[16] Temam R. (1991) Approximation of Attractors, Large Eddy Simulations and Multiscale

[17] Ladyzhenskaya O.A. (1985) On the finiteness of the dimention of bounded invarian sets for the Navier-Stokes equations and other related dissipative systems.// J.Soviet Math.

[18] Chepyzhov V., Vishik M. (2002) Attractors for Equations of Mathematical Physics.//

Amer. Math. Soc. Colloq. Publ., Vol. 49, Amer. Math. Soc., Providence, RI.

**Figure 21.** Triple period stable cycle at point, *Pr* <sup>=</sup> 1.354839, *Ra* <sup>=</sup> 1.308 · 106; Stream lines in <sup>Ω</sup> for the same *Pr* and *Ra*.

3. FSM scenario with subharmonic cascade of bifurcations of stable limit cycles.
