**3.1. Initial-boundary value problem, mesh adaptation and benchmark verification**

The geometry and boundary conditions are taken from [32] with small adaptation, since laminar-turbulent transition is investigated as a dynamic system so we compare not only benchmark results but bifurcation sequences as well.

Domain Ω with local Lipschiz-continuous boundaries *∂*Ω*<sup>i</sup>* is represented by a rectangular channel divided by two unequal parts in Cartesian coordinates. The length of the domain spans in X axis direction. The first part of the channel has length *L*<sup>1</sup> = 2.0, the second part of the channel length is *L*<sup>2</sup> = 10.0, so the whole length of the domain is *L* = *L*<sup>1</sup> + *L*<sup>2</sup> = 12.0. Height spans in Y axis and the second part of the domain is higher than the first one by the size of the step *h* = 0.6 with the first part height of *H* = 0.9. So the second part of the domain height is *H* + *h* = 1.5. Width of the domain *W* = 3.5 spans in Z axis direction and is the same for both part. The geometry is almost identical to [32] with *L*<sup>1</sup> = 10.0. The formed step causes the flow to create recirculation zones inside Ω and transit to turbulence with the growth of Reynolds number.

There are many different boundary conditions available in papers for this particular problem [33–35, 38, 39]. All boundaries in step Y axis direction *∂*Ω<sup>1</sup> ∈ min*y*(Ω) are given by solid walls with no slip condition. The upper boundary *∂*Ω<sup>2</sup> ∈ max*y*(Ω) is given as either wall or symmetry. Boundaries *∂*Ω<sup>3</sup> ∈ max*z*(Ω) and *∂*Ω<sup>4</sup> ∈ min*z*(Ω) vary from one work to another. The boundary *∂*Ω<sup>3</sup> is given as a solid wall but *∂*Ω<sup>4</sup> as symmetry plane in [33]. Boundaries in Z direction in some other papers, i.e. [35, 40], are given by symmetry planes thus modeling a semi-2D problem and decreasing number of mesh elements due to the lack of side boundary layers. We choose boundary conditions like in [32] where boundaries in Z direction are no slip walls. Boundary *∂*Ω<sup>5</sup> ∈ min*x*(Ω) is the inflow boundary where Poiseuille laminar solution for R=200 in the channel with the same width and height. Boundary *∂*Ω<sup>6</sup> ∈ max*x*(Ω) is the outflow boundary. We introduce smooth outflow boundary conditions [41] with the outflow buffer so that unphysical outflow conditions are neglected. Pressure is set on the outflow boundary as a reference value. Initial conditions in whole domain Ω for velocity vector-function and scalar pressure function obtained by laminar stationary solution of this problem with *R* = 200 so that ∇ · **V** = 0 at t=0. Reynolds number is defined as:

$$R = V^0 L / \upsilon\_\prime L = 2WH / (\mathcal{W} + H). \tag{37}$$

**Figure 2.** Calculated values of *lk* for *R* = 1000. Central section and literal 3D sections are shown.

functions is considered as well as friction wall coefficients.

with high precision using shareware GetData Graph Digitizer 2.24.

presented data, especially with DNS data (rectangles) from [45].

boundaries in Z direction.

with reference data.

While solving various modifications of the problem we changed height of the step in accordance with the geometry presented for particular benchmark tests and here all detail data of these geometries is omitted due to brevity. Integral and statistical data of **V** and *P*

FSM Scenarios of Laminar-Turbulent Transition in Incompressible Fluids 263

First we overseen qualitative comparison of results for R=1500 with side walls with DNS [35], [42], [37] and physical model [43] results. Velocity isolines sections are presented in fig.3. One can clearly see the 3D structure of the flow plus the recirculation vertexes that appear near all walls. This agrees well with the presented papers. To compare qualitatively periodic boundary results we used papers [44] and [35]. The flow is much more simple for periodic

For the quantitative analysis we considered some papers: [44] for *R* = {100; 389; 1000}, [34] for *R* = 150 − 800, [35] for *R* = 1500, [45] in *Cf* comparison for *R* = 5500 (using LES model) and [38, 46] for virtual stress tensor comparison for *R* = 1800 (using LES model). All data agreed well for benchmarks with maximum difference of 7%. Data was extracted from papers

Detail results are omitted for the sake of brevity, but some results are presented here. Reattachment vortex length ration to the height of the step is presented in fig.4 and compared with results from various sources, i.e. [34]. One can clearly see that the results are well agreed

More quantitative results are presented in fig. 5, where a mean wall friction coefficient is presented. The dynamic eddy viscosity model was used to solve the flow for *R* = 5500, and was compared with other LES and DNS results. One can see that the result agree well with

Here *V*<sup>0</sup> is the inflow velocity normalized by the flow rate, like in [32]. Due to the fact that DNS is very computational demanding we are limiting the maximum R by analyzing literature results. It is shown in [32] that a "chaotic attractor" that corresponds to noise frequencies band appear at *R* > 1500 so we set maximum *R* = 1500. For the purpose of minimizing computational efforts we are considering three meshes adopted for *R* = 700,*R* = 1000 and *R* = 1500. The rest Reynolds number meshes are lying inside the segment. Using an upper bound estimate (14) and taking modified wave number analysis in account one gets total mesh numbers as: 5 040 813(R=700), 11 246 827(R=1000) and 24 004 984 (R=1500). We set *L* = 1 <sup>3</sup> 7 from the given geometry and since we know *R*, we can evaluate minimal Kolmogorov linear scale of one element *lk*:

$$l\_k = \{0.024495281; 0.018745885; 0.013830488\}\,\tag{38}$$

and, hence, in Cartesian coordinates we have *M* = {491; 641; 869}, *N* = {62; 81; 109}, *K* = {144; 188; 254}, elements for the upper bound estimate. Time step is selected 0.005 for all calculations [**?** ]. We performed these three control calculations on fine grids and calculated isosurfaces of minimal scales *lk* by (27). Example of these isosurfaces for *R* = 1000 is shown in fig.2.

It is clearly visible that we must use fines mesh with element size of 0.01577 of boundary layers and recirculation zone, but the rest of the flow has s much lighter demand for small scales. So final meshes where chosen to be:

$$\begin{aligned} G1 &= \{245X54X83\}, R = 700; \Delta\mathbf{x}\_{\mathfrak{d}} = \{0.0677; 0.01957; 0.03964\} \approx 103^{\mathfrak{d}};\\ G2 &= \{339X72X19\}, R = 1000; \Delta\mathbf{x}\_{\mathfrak{d}} = \{0.04696; 0.01577; 0.028548\} \approx 143^{\mathfrak{d}};\\ G3 &= \{560X100X177\}, R = 1500; \Delta\mathbf{x}\_{\mathfrak{d}} = \{0.04696; 0.01577; 0.028548\} \approx 215^{\mathfrak{d}};\end{aligned} \tag{39}$$

where cubic values are shown for required memory estimate. After that all other calculations for the given problem are performed on these (39) three grids that satisfy necessary conditions of theorems. In order to perform benchmark tests boundaries *∂*Ω3,4 where changed accordingly. Values of are the same *lk* for inside flow region and near periodic boundaries so our grid resolution (39) is overdensed for test calculations with periodic boundaries. So new mesh was generated in the same manner through calculation of *lk* but it contained less number of *lk* elements, so benchmark calculations took less time to solve.

12 Will-be-set-by-IN-TECH

Z direction in some other papers, i.e. [35, 40], are given by symmetry planes thus modeling a semi-2D problem and decreasing number of mesh elements due to the lack of side boundary layers. We choose boundary conditions like in [32] where boundaries in Z direction are no slip walls. Boundary *∂*Ω<sup>5</sup> ∈ min*x*(Ω) is the inflow boundary where Poiseuille laminar solution for R=200 in the channel with the same width and height. Boundary *∂*Ω<sup>6</sup> ∈ max*x*(Ω) is the outflow boundary. We introduce smooth outflow boundary conditions [41] with the outflow buffer so that unphysical outflow conditions are neglected. Pressure is set on the outflow boundary as a reference value. Initial conditions in whole domain Ω for velocity vector-function and scalar pressure function obtained by laminar stationary solution of this

Here *V*<sup>0</sup> is the inflow velocity normalized by the flow rate, like in [32]. Due to the fact that DNS is very computational demanding we are limiting the maximum R by analyzing literature results. It is shown in [32] that a "chaotic attractor" that corresponds to noise frequencies band appear at *R* > 1500 so we set maximum *R* = 1500. For the purpose of minimizing computational efforts we are considering three meshes adopted for *R* = 700,*R* = 1000 and *R* = 1500. The rest Reynolds number meshes are lying inside the segment. Using an upper bound estimate (14) and taking modified wave number analysis in account one gets total mesh numbers as: 5 040 813(R=700), 11 246 827(R=1000) and 24 004 984 (R=1500). We set *L* = 1 <sup>3</sup>

from the given geometry and since we know *R*, we can evaluate minimal Kolmogorov linear

and, hence, in Cartesian coordinates we have *M* = {491; 641; 869}, *N* = {62; 81; 109}, *K* = {144; 188; 254}, elements for the upper bound estimate. Time step is selected 0.005 for all calculations [**?** ]. We performed these three control calculations on fine grids and calculated isosurfaces of minimal scales *lk* by (27). Example of these isosurfaces for *R* = 1000 is shown

It is clearly visible that we must use fines mesh with element size of 0.01577 of boundary layers and recirculation zone, but the rest of the flow has s much lighter demand for small scales. So

where cubic values are shown for required memory estimate. After that all other calculations for the given problem are performed on these (39) three grids that satisfy necessary conditions of theorems. In order to perform benchmark tests boundaries *∂*Ω3,4 where changed accordingly. Values of are the same *lk* for inside flow region and near periodic boundaries so our grid resolution (39) is overdensed for test calculations with periodic boundaries. So new mesh was generated in the same manner through calculation of *lk* but it contained less

number of *lk* elements, so benchmark calculations took less time to solve.

*<sup>G</sup>*<sup>1</sup> <sup>=</sup> {245*X*54*X*83} , *<sup>R</sup>* <sup>=</sup> 700; <sup>Δ</sup>*x<sup>α</sup>* <sup>=</sup> {0.0677; 0.01957; 0.03964} <sup>≈</sup> 1033; *<sup>G</sup>*<sup>2</sup> <sup>=</sup> {339*X*72*X*119} , *<sup>R</sup>* <sup>=</sup> 1000; <sup>Δ</sup>*x<sup>α</sup>* <sup>=</sup> {0.04696; 0.01577; 0.028548} <sup>≈</sup> 1433; *<sup>G</sup>*<sup>3</sup> <sup>=</sup> {560*X*100*X*177} , *<sup>R</sup>* <sup>=</sup> 1500; <sup>Δ</sup>*x<sup>α</sup>* <sup>=</sup> {0.04696; 0.01577; 0.028548} <sup>≈</sup> 2153;

*R* = *V*0*L*/*ν*, *L* = 2*WH*/(*W* + *H*). (37)

*lk* = {0.024495281; 0.018745885; 0.013830488} , (38)

7

(39)

problem with *R* = 200 so that ∇ · **V** = 0 at t=0. Reynolds number is defined as:

scale of one element *lk*:

final meshes where chosen to be:

in fig.2.

**Figure 2.** Calculated values of *lk* for *R* = 1000. Central section and literal 3D sections are shown.

While solving various modifications of the problem we changed height of the step in accordance with the geometry presented for particular benchmark tests and here all detail data of these geometries is omitted due to brevity. Integral and statistical data of **V** and *P* functions is considered as well as friction wall coefficients.

First we overseen qualitative comparison of results for R=1500 with side walls with DNS [35], [42], [37] and physical model [43] results. Velocity isolines sections are presented in fig.3. One can clearly see the 3D structure of the flow plus the recirculation vertexes that appear near all walls. This agrees well with the presented papers. To compare qualitatively periodic boundary results we used papers [44] and [35]. The flow is much more simple for periodic boundaries in Z direction.

For the quantitative analysis we considered some papers: [44] for *R* = {100; 389; 1000}, [34] for *R* = 150 − 800, [35] for *R* = 1500, [45] in *Cf* comparison for *R* = 5500 (using LES model) and [38, 46] for virtual stress tensor comparison for *R* = 1800 (using LES model). All data agreed well for benchmarks with maximum difference of 7%. Data was extracted from papers with high precision using shareware GetData Graph Digitizer 2.24.

Detail results are omitted for the sake of brevity, but some results are presented here. Reattachment vortex length ration to the height of the step is presented in fig.4 and compared with results from various sources, i.e. [34]. One can clearly see that the results are well agreed with reference data.

More quantitative results are presented in fig. 5, where a mean wall friction coefficient is presented. The dynamic eddy viscosity model was used to solve the flow for *R* = 5500, and was compared with other LES and DNS results. One can see that the result agree well with presented data, especially with DNS data (rectangles) from [45].

**Figure 5.** Mean wall friction coefficient *Cf* for *R* = 5500, line - current method, pints - reference data

FSM Scenarios of Laminar-Turbulent Transition in Incompressible Fluids 265

**Figure 6.** Reynolds stress correlations comparison for *R* = 1800 with experimental data at length

are formed by velocity vector components. Under infinite dimensional phase space here we consider finite dimensional phase space produced by the numerical system (30) whose space dimension is greater than the attractor dimension of the system (15) for the given initial-boundary value problem with a finite preset Reynolds number. It is important to notice that there exists a hysteresis of a solution if one approaches a fixed R from different sides, say solution *G* exists for *R*<sup>0</sup> if we travel to it by *R*<sup>1</sup> → *R*0, *R*<sup>1</sup> < *R*<sup>0</sup> but does not exist for *R*<sup>0</sup> if *R*<sup>1</sup> > *R*0, i.e. see [49]. So we are only considering the following cascade of bifurcation

The solution is laminar for *R* from 100 up to 736 but the time of stationary solution formulation increases as *R* grows. For the stationary laminar solution we can monitor a fixed point in infinite dimensional phase space and in all five three dimensional subspaces of velocity vector function. Starting from a system exhibits one frequency mode regime. A limited cycle *C*<sup>1</sup> is formed from every stationary point in phase space and has projections in all subspaces. One

At the point of *R* = 850, the cycle looses stability and forms a stable two dimensional invariant torus *T*<sup>2</sup> = *C*<sup>1</sup> ⊗ *C*<sup>2</sup> as the result of Andronov-Hopf bifurcation. This torus is located in all infinite dimensional phase space and can be found in all subspaces for points *Pi*. Presumably, this is due to the fact of incompressibility what expresses in elliptic operator for pressure

*x*/*h* = {4; 6; 15} in the center (W/2).

parameters *R*<sup>1</sup> < *R*<sup>2</sup> < ... < *Rn* = 1500.

can see a projection of the cycle for the point *P*<sup>1</sup> in fig. 7.

**Figure 3.** Instantaneous velocity vector-function modulus isolines in sections for R=1500.

**Figure 4.** Mean reattachment length to the step height, compared with reference data, line - current method

Comparison for Reynolds stress correlations at *R* = 1800 with data from [38] presented in fig.6. One can see a very good agreement between DNS and experimental data.

More benchmark details are omitted for the sake of brevity, see [23, 47].

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

The idea of analysis was used many times in our research [7, 8, 23, 48] and corresponds with [32]. We saved data for velocity function in format *Vx*, *Vy*, *Vz*. Saving data for all domain was impossible since there's not enough disk space (it requires about 6,43 terabytes for one fixed Reynolds number), so we choose several points and saved data from them. Points in Cartesian coordinates normalized by length in each direction are: *p*1 = {0.1, 0.5, 0.5}; *p*2 = {0.2, 0.5, 0.5}; *p*3 = {0.5, 0.5, 0.5}; *p*4 = {0.7, 0.5, 0.5}; *p*5 = {0.8, 0.1, 0.1}. We are using Reynolds number as the bifurcation parameter for the problem and forming sets of three dimensional phase subspaces of the whole infinite dimensional phase space. The subspaces

14 Will-be-set-by-IN-TECH

**Figure 3.** Instantaneous velocity vector-function modulus isolines in sections for R=1500.

**Figure 4.** Mean reattachment length to the step height, compared with reference data, line - current

fig.6. One can see a very good agreement between DNS and experimental data.

More benchmark details are omitted for the sake of brevity, see [23, 47].

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

Comparison for Reynolds stress correlations at *R* = 1800 with data from [38] presented in

The idea of analysis was used many times in our research [7, 8, 23, 48] and corresponds with [32]. We saved data for velocity function in format *Vx*, *Vy*, *Vz*. Saving data for all domain was impossible since there's not enough disk space (it requires about 6,43 terabytes for one fixed Reynolds number), so we choose several points and saved data from them. Points in Cartesian coordinates normalized by length in each direction are: *p*1 = {0.1, 0.5, 0.5}; *p*2 = {0.2, 0.5, 0.5}; *p*3 = {0.5, 0.5, 0.5}; *p*4 = {0.7, 0.5, 0.5}; *p*5 = {0.8, 0.1, 0.1}. We are using Reynolds number as the bifurcation parameter for the problem and forming sets of three dimensional phase subspaces of the whole infinite dimensional phase space. The subspaces

method

**Figure 5.** Mean wall friction coefficient *Cf* for *R* = 5500, line - current method, pints - reference data

**Figure 6.** Reynolds stress correlations comparison for *R* = 1800 with experimental data at length *x*/*h* = {4; 6; 15} in the center (W/2).

are formed by velocity vector components. Under infinite dimensional phase space here we consider finite dimensional phase space produced by the numerical system (30) whose space dimension is greater than the attractor dimension of the system (15) for the given initial-boundary value problem with a finite preset Reynolds number. It is important to notice that there exists a hysteresis of a solution if one approaches a fixed R from different sides, say solution *G* exists for *R*<sup>0</sup> if we travel to it by *R*<sup>1</sup> → *R*0, *R*<sup>1</sup> < *R*<sup>0</sup> but does not exist for *R*<sup>0</sup> if *R*<sup>1</sup> > *R*0, i.e. see [49]. So we are only considering the following cascade of bifurcation parameters *R*<sup>1</sup> < *R*<sup>2</sup> < ... < *Rn* = 1500.

The solution is laminar for *R* from 100 up to 736 but the time of stationary solution formulation increases as *R* grows. For the stationary laminar solution we can monitor a fixed point in infinite dimensional phase space and in all five three dimensional subspaces of velocity vector function. Starting from a system exhibits one frequency mode regime. A limited cycle *C*<sup>1</sup> is formed from every stationary point in phase space and has projections in all subspaces. One can see a projection of the cycle for the point *P*<sup>1</sup> in fig. 7.

At the point of *R* = 850, the cycle looses stability and forms a stable two dimensional invariant torus *T*<sup>2</sup> = *C*<sup>1</sup> ⊗ *C*<sup>2</sup> as the result of Andronov-Hopf bifurcation. This torus is located in all infinite dimensional phase space and can be found in all subspaces for points *Pi*. Presumably, this is due to the fact of incompressibility what expresses in elliptic operator for pressure

Further increase of *R* leads to more complex topology of cycles that are forming the torus. The

FSM Scenarios of Laminar-Turbulent Transition in Incompressible Fluids 267

process of this complication is shown in fig.9 at point *P*<sup>3</sup> and in fig.10 at point *P*5.

**Figure 9.** Torus complication in subspace projection at point *P*3. *R* = {851; 883}

**Figure 10.** Torus complication in subspace projection at point *P*5, *R* = {851; 882; 883; 885}

projections in different subspaces are shown in fig.9, 10 for *R* ≥ 883.

dimensional torus structure in fig.12 in sections for R=883.

This process continues up to *R* = 883. It was found that at this point the two dimensional torus looses stability in the whole infinite dimensional space and forms a three dimensional invariant torus through the Andronov-Hopf bifurcation. This torus is formed by the topological multiplication of three irrational frequency cycles *T*<sup>3</sup> = *C*<sup>1</sup> ⊗ *C*<sup>2</sup> ⊗ *C*3. Its

Since the space dimension in one point of the phase subspace is the same as the formed torus dimension one should increase the phase subspace dimension. It can be clearly seen on sections fig.11 where a new "coiling" can be monitored around the old cycle in plane section. To do so we are taking data from another point thus forming a four dimensional subspace and then we are performing two sections by two dimensional planes. The first section in *P*<sup>5</sup> (plus data from *P*<sup>2</sup> marked with asterisk) is shown in fig.12 on the left. The torus structure can be seen by performing another section by the plane *W* = 0.005. One can see the three

**Figure 7.** Point *P*1, stable limit cycle, R=740

and hence the acoustic speed is infinite. An example of this torus is shown in fig.8. It is clearly seen that the system is very sensible for the values of bifurcation parameters. It is shown in fig.8, that a change of *R* by 1.0 changes the attractor from *C*<sup>1</sup> to *T*2. Comparing these results with the mentioned work [32] we can say that the formation of the first cycle appears almost at the same Reynolds number, here we have *R* = 737 and in [32] the value is R=735. The formation of two-frequent mode in [32] appears at R=855, but there's no stable phase space trajectory available, since numerical errors started dominating. And further investigation in [32] is performed using frequency analysis that indicated creation of the other independent frequency. The recent report in TsAGI [50] of Sibgatullin I.N. indicated that other initial-boundary problems exhibit the resembling scenario by which an invariant torus is formatted in phase space. But all reports and papers indicate that a chaotic behavior follows the formation of two dimensional invariant tori. We suppose that this is due to the fact that numerical methods, used in the papers and reports, don't meet necessary conditions for theorems and assertions outlined here and so we continued the numerical analysis.

**Figure 8.** *R* = 849 Phase subspace cycle projection and section, point P4(left); *R* = 850 Phase subspace torus projection and section, point P4(right)

Further increase of *R* leads to more complex topology of cycles that are forming the torus. The process of this complication is shown in fig.9 at point *P*<sup>3</sup> and in fig.10 at point *P*5.

**Figure 9.** Torus complication in subspace projection at point *P*3. *R* = {851; 883}

16 Will-be-set-by-IN-TECH

and hence the acoustic speed is infinite. An example of this torus is shown in fig.8. It is clearly seen that the system is very sensible for the values of bifurcation parameters. It is shown in fig.8, that a change of *R* by 1.0 changes the attractor from *C*<sup>1</sup> to *T*2. Comparing these results with the mentioned work [32] we can say that the formation of the first cycle appears almost at the same Reynolds number, here we have *R* = 737 and in [32] the value is R=735. The formation of two-frequent mode in [32] appears at R=855, but there's no stable phase space trajectory available, since numerical errors started dominating. And further investigation in [32] is performed using frequency analysis that indicated creation of the other independent frequency. The recent report in TsAGI [50] of Sibgatullin I.N. indicated that other initial-boundary problems exhibit the resembling scenario by which an invariant torus is formatted in phase space. But all reports and papers indicate that a chaotic behavior follows the formation of two dimensional invariant tori. We suppose that this is due to the fact that numerical methods, used in the papers and reports, don't meet necessary conditions for theorems and assertions outlined here and so we continued the numerical analysis.

**Figure 8.** *R* = 849 Phase subspace cycle projection and section, point P4(left); *R* = 850 Phase subspace

**Figure 7.** Point *P*1, stable limit cycle, R=740

torus projection and section, point P4(right)

**Figure 10.** Torus complication in subspace projection at point *P*5, *R* = {851; 882; 883; 885}

This process continues up to *R* = 883. It was found that at this point the two dimensional torus looses stability in the whole infinite dimensional space and forms a three dimensional invariant torus through the Andronov-Hopf bifurcation. This torus is formed by the topological multiplication of three irrational frequency cycles *T*<sup>3</sup> = *C*<sup>1</sup> ⊗ *C*<sup>2</sup> ⊗ *C*3. Its projections in different subspaces are shown in fig.9, 10 for *R* ≥ 883.

Since the space dimension in one point of the phase subspace is the same as the formed torus dimension one should increase the phase subspace dimension. It can be clearly seen on sections fig.11 where a new "coiling" can be monitored around the old cycle in plane section.

To do so we are taking data from another point thus forming a four dimensional subspace and then we are performing two sections by two dimensional planes. The first section in *P*<sup>5</sup> (plus data from *P*<sup>2</sup> marked with asterisk) is shown in fig.12 on the left. The torus structure can be seen by performing another section by the plane *W* = 0.005. One can see the three dimensional torus structure in fig.12 in sections for R=883.

**Figure 11.** Plane sections of phase subspaces in *P*<sup>4</sup> and *P*<sup>5</sup> for *R* = (882.5 → 883) from left to right

**Figure 13.** Section of the four dimensional phase subspaces by additional planes at point *P*<sup>3</sup> + *P*2(∗)

results show that a double period three dimensional torus is formed in all infinite dimensional phase space for *R* = 883.8. Further investigation is impossible due to the exponential growth of bifurcations. The second section view is already next to incomprehensible for *R* = 890 and consists of a filled black squire. In conclusion we found the following scenario of

It can be seen in (40) that the scenario exhibits initial stage of Landau-Hopf scenario (Hopf bifurcation cascades) and initial stage of FSM scenario, since the Feigenbaum scenario is progressing after the three dimensional torus. It corresponds with the initial stage of FSM scenario. One should point out that a three dimensional torus is stable at least for the time of

One of the best studied and widely analyzed problems of fluid mechanics is the Rayleigh-Benard natural convection problem. The problem has been considered by many and has lots of results in the field of nonlinear bifurcation analysis, analytical, numerical and

**4.1. The Oberbeck-Boussinesq approximation of Rayleigh-Benard convection,**

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

**4. Laminar-turbulent transition for Rayleigh-Benard convection**

*C*<sup>1</sup> → *T*<sup>2</sup> → *T*<sup>3</sup> → *T*3⊗<sup>2</sup> →, (40)

FSM Scenarios of Laminar-Turbulent Transition in Incompressible Fluids 269

(top) and *P*<sup>5</sup> + *P*2(∗) (bottom)

laminar-turbulent transition:

experimental.

numerical simulation for 1, 5 · 107 timesteps at *<sup>R</sup>* <sup>=</sup> 883.

**dimensionless form and benchmarks**

**Figure 12.** Section of the phase subspace in Point *P*<sup>3</sup> (top) and *P*<sup>5</sup> (bottom) plus *P*<sup>2</sup> data (marked asterisk), *R* = 883, on the left. Section of the four dimensional phase subspace by the additional plane, on the right

Three dimensional invariant torus losses its stability and through the doubling period bifurcation forms a double period three dimensional torus for *R* = 883.8 It requires huge amount of data to get results shown in fig.13. Depicted results required about 375Gb of data per one point and could only be completed by numerical methods on GPUs or on massive parallel clusters. All different sections are presented in fig.13 for four dimensional phase subspaces at *P*<sup>3</sup> + *P*<sup>2</sup> and *P*<sup>5</sup> + *P*2.

It can be clearly seen that the torus topology is more complicated, i.e. see fig.13 at point *P*5, and fig.12 - second section, although the bifurcation parameter only changed by 0.091%. These

18 Will-be-set-by-IN-TECH

**Figure 11.** Plane sections of phase subspaces in *P*<sup>4</sup> and *P*<sup>5</sup> for *R* = (882.5 → 883) from left to right

**Figure 12.** Section of the phase subspace in Point *P*<sup>3</sup> (top) and *P*<sup>5</sup> (bottom) plus *P*<sup>2</sup> data (marked asterisk), *R* = 883, on the left. Section of the four dimensional phase subspace by the additional plane, on the right

Three dimensional invariant torus losses its stability and through the doubling period bifurcation forms a double period three dimensional torus for *R* = 883.8 It requires huge amount of data to get results shown in fig.13. Depicted results required about 375Gb of data per one point and could only be completed by numerical methods on GPUs or on massive parallel clusters. All different sections are presented in fig.13 for four dimensional phase

It can be clearly seen that the torus topology is more complicated, i.e. see fig.13 at point *P*5, and fig.12 - second section, although the bifurcation parameter only changed by 0.091%. These

subspaces at *P*<sup>3</sup> + *P*<sup>2</sup> and *P*<sup>5</sup> + *P*2.

**Figure 13.** Section of the four dimensional phase subspaces by additional planes at point *P*<sup>3</sup> + *P*2(∗) (top) and *P*<sup>5</sup> + *P*2(∗) (bottom)

results show that a double period three dimensional torus is formed in all infinite dimensional phase space for *R* = 883.8. Further investigation is impossible due to the exponential growth of bifurcations. The second section view is already next to incomprehensible for *R* = 890 and consists of a filled black squire. In conclusion we found the following scenario of laminar-turbulent transition:

$$T\_1 \to T\_2 \to T\_3 \to T\_{3 \otimes 2} \to \text{ } \tag{40}$$

It can be seen in (40) that the scenario exhibits initial stage of Landau-Hopf scenario (Hopf bifurcation cascades) and initial stage of FSM scenario, since the Feigenbaum scenario is progressing after the three dimensional torus. It corresponds with the initial stage of FSM scenario. One should point out that a three dimensional torus is stable at least for the time of numerical simulation for 1, 5 · 107 timesteps at *<sup>R</sup>* <sup>=</sup> 883.
