**2.1. Some theorems and assertions**

Here we are talking only about three dimensional evolutionary incompressible Navier-Stokes equations:

$$\begin{aligned} \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}, \end{aligned} \tag{1}$$

*ki* = (*k*0, *k*1, *k*2, ..., *kk*), *i* ∈ **N** (3)

FSM Scenarios of Laminar-Turbulent Transition in Incompressible Fluids 253

∇2*V*(*x*, *<sup>t</sup>*)

 2 *x*∈Ω *dx*

, (6)

, (9)

(4)

(7)

(8)

, (5)

Hence the minimum space dimension that is required to hold all possible flow scales in

 *kk k*0 <sup>3</sup>

The value of *Y* gives us number of degrees of freedom for the given turbulent flow. We can

1 *t t* <sup>0</sup> sup 

**Theorem 01.** *[16] the Hausdorff-Besicovitch dimension of any attractor A in the Navier-Stokes*

*kk k*0

∇2**V**(*x*, *<sup>t</sup>*)

∇2**V**(*x*, *<sup>t</sup>*)

*ν*

 2

 |∇2**V**(*x*, *<sup>t</sup>*)<sup>|</sup> *ν*

> 

3/2

dim *A* ≤ *C* ·

 

1/4 =

where *L* is a macroscopic problem scale. Since all functions in (9) have finite limits then the dimension of an attractor is always limited from above and its dimension is linked with

**Assertion 01.** *If a well-posed initial-boundary value problem for Navier-Stokes equations in bounded domain* <sup>Ω</sup> <sup>∈</sup> **<sup>R</sup>**<sup>3</sup> *with local Lipschitz-continuous boundary <sup>∂</sup>*<sup>Ω</sup> <sup>∈</sup> **<sup>R</sup>**<sup>2</sup> *has a suitable numerical approximation with degrees of freedom Y*<sup>∗</sup> ≥ *Y that convergence to the given problem on every time step of approximation then the approximated numerical attractor A*∗ *converges to the real attractor A.* This assertion can be easily verified using estimates for finite-difference and finite-element approximations of Roger Temam [20]. This assertion answers the first question. Now we

*Y* =

use the value for *�* from [16] where it is shown that *�* is limited by the supremum:

sup **V**(*x*,0)⊂*A*

*where C is an arbitrary constant.* One can derive a more explicit estimate. Since:

*�* ≤ *�* = *ν*

 *� ν*3

dim *<sup>A</sup>* <sup>≤</sup> *<sup>L</sup>*<sup>3</sup> ·

bounded domain **<sup>V</sup>** <sup>∈</sup> **<sup>R</sup>**<sup>3</sup> is defined as:

*�* = *ν* · lim

and we can find the following:

so taking (4) into account:

and

*<sup>t</sup>*→<sup>∞</sup> sup

*equations for a well-posed boundary value problem is limited by:*

*kk* ≤ *kk* =

number of degrees of freedom (4). So we have the following:

and suppose that laminar-turbulent transition is well described by these equations. Here *ρ*<sup>0</sup> is a constant fluid density, *ν* is the kinematic fluid viscosity, **V** is the velocity vector-function and *P* is the scalar pressure function. Main questions of numerical solution trustworthiness of Navier-Stokes equations and its attraction to real system (1) attractor are:

1. Does the change of infinite dimensional system for finite-dimensional alters the attractor and solution?


Some results of attractor approximation and numerical attractors for nonlinear PDEs were obtained in papers [13, 14], but none of them covers such a complicated topic as existence of and numerical trajectories attraction for Navier-Stokes attractors. However in papers [15–19] it is shown that for a well-posed initial-boundary problem there exists an attractor and its dimension and volume is limited from above. This allows us to pose some assumptions on numerical methods that must be used to get a numerical system attractor and that it is close to the attractor of Navier-Stokes equations. To do so we will be following closely to work of R.Temam [16, 20].

Existence of a global attractor for 3D Navier-Stokes equations is proven [15, 16]. Let that attractor *A* exist for a well-posed initial-boundary problem in domain Ω ∈ **R** with local Lipschitz-contonous boundary *<sup>∂</sup>*Ω*<sup>i</sup>* and velocity vector-function **<sup>V</sup>** <sup>∈</sup> **<sup>R</sup>**3. Let *<sup>k</sup>*<sup>0</sup> be the macroscopic wave number that corresponds with the characteristic length scale *L* and *kk* is the Kolmogorov wave number defined as:

$$k\_k = \left(\frac{\epsilon}{\nu^3}\right)^{1/4} \tag{2}$$

*ν* - kinematic fluid viscosity and *�* - rate of turbulent energy dissipation. whole spectra of wave numbers can be described as:

$$k\_i = (k\_0, k\_1, k\_2, \dots, k\_k), i \in \mathbb{N} \tag{3}$$

Hence the minimum space dimension that is required to hold all possible flow scales in bounded domain **<sup>V</sup>** <sup>∈</sup> **<sup>R</sup>**<sup>3</sup> is defined as:

$$Y = \left(\frac{k\_k}{k\_0}\right)^3 \tag{4}$$

The value of *Y* gives us number of degrees of freedom for the given turbulent flow. We can use the value for *�* from [16] where it is shown that *�* is limited by the supremum:

$$\epsilon = \nu \cdot \lim\_{t \to \infty} \sup \left\{ \sup\_{\mathbf{V}(\mathbf{x}, \mathbf{0}) \subset A} \frac{1}{t} \int\_0^t \sup \left| \nabla^2 V(\mathbf{x}, t) \right|^2\_{\mathbf{x} \in \Omega} d\mathbf{x} \right\},\tag{5}$$

and we can find the following:

**Theorem 01.** *[16] the Hausdorff-Besicovitch dimension of any attractor A in the Navier-Stokes equations for a well-posed boundary value problem is limited by:*

$$\dim A \le \mathbb{C} \cdot \frac{k\_k}{k\_0} \,\tag{6}$$

*where C is an arbitrary constant.* One can derive a more explicit estimate. Since:

$$
\epsilon \le \overline{\epsilon} = \nu \left| \nabla^2 \mathbf{V}(\mathbf{x}, t) \right|^2 \tag{7}
$$

and

2 Will-be-set-by-IN-TECH

consider two well studied problems for incompressible Navier-Stokes equations, namely flow over a backward facing step and Rayleigh-Benard convection in cubic cavity. Numerical solutions of these problems for transitional regimes indicated existence of complicated scenarios formed by theory FSM. Thus, it seems reasonable, that there is no unified laminar-turbulent transition scenario, it can be a cascade of stable limit cycles or stable two dimensional or many dimensional tori, but all these scenarios lay in the frameworks of the

Here we are talking only about three dimensional evolutionary incompressible Navier-Stokes

∇ · **V** = 0,

and suppose that laminar-turbulent transition is well described by these equations. Here *ρ*<sup>0</sup> is a constant fluid density, *ν* is the kinematic fluid viscosity, **V** is the velocity vector-function and *P* is the scalar pressure function. Main questions of numerical solution trustworthiness of

1. Does the change of infinite dimensional system for finite-dimensional alters the attractor

Some results of attractor approximation and numerical attractors for nonlinear PDEs were obtained in papers [13, 14], but none of them covers such a complicated topic as existence of and numerical trajectories attraction for Navier-Stokes attractors. However in papers [15–19] it is shown that for a well-posed initial-boundary problem there exists an attractor and its dimension and volume is limited from above. This allows us to pose some assumptions on numerical methods that must be used to get a numerical system attractor and that it is close to the attractor of Navier-Stokes equations. To do so we will be following closely to work of

Existence of a global attractor for 3D Navier-Stokes equations is proven [15, 16]. Let that attractor *A* exist for a well-posed initial-boundary problem in domain Ω ∈ **R** with local Lipschitz-contonous boundary *<sup>∂</sup>*Ω*<sup>i</sup>* and velocity vector-function **<sup>V</sup>** <sup>∈</sup> **<sup>R</sup>**3. Let *<sup>k</sup>*<sup>0</sup> be the macroscopic wave number that corresponds with the characteristic length scale *L* and *kk* is

> *� ν*3 1/4

*ν* - kinematic fluid viscosity and *�* - rate of turbulent energy dissipation. whole spectra of

*kk* =

2. Do numerical approximation errors have crucial affect on the attractor trajectory? 3. Does a numerical solution convergence to a real solution and how close are they?

<sup>0</sup> <sup>∇</sup>*<sup>P</sup>* <sup>=</sup> *<sup>ν</sup>*∇2**V**, (1)

(2)

**2. Construction of high order numerical method for Navier-Stokes**

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

Navier-Stokes equations and its attraction to real system (1) attractor are:

FSM-theory.

**equations**

equations:

and solution?

R.Temam [16, 20].

the Kolmogorov wave number defined as:

wave numbers can be described as:

**2.1. Some theorems and assertions**

*∂***V**

$$k\_k \le \overline{k\_k} = \left(\frac{\overline{\varepsilon}}{\nu^3}\right)^{1/4} = \sqrt{\left(\frac{|\nabla^2 \mathbf{V}(\mathbf{x}, t)|}{\nu}\right)}\tag{8}$$

so taking (4) into account:

$$\dim A \le L^3 \cdot \left(\frac{|\nabla^2 \mathbf{V}(\mathbf{x}, t)|}{\nu}\right)^{3/2},\tag{9}$$

where *L* is a macroscopic problem scale. Since all functions in (9) have finite limits then the dimension of an attractor is always limited from above and its dimension is linked with number of degrees of freedom (4). So we have the following:

**Assertion 01.** *If a well-posed initial-boundary value problem for Navier-Stokes equations in bounded domain* <sup>Ω</sup> <sup>∈</sup> **<sup>R</sup>**<sup>3</sup> *with local Lipschitz-continuous boundary <sup>∂</sup>*<sup>Ω</sup> <sup>∈</sup> **<sup>R</sup>**<sup>2</sup> *has a suitable numerical approximation with degrees of freedom Y*<sup>∗</sup> ≥ *Y that convergence to the given problem on every time step of approximation then the approximated numerical attractor A*∗ *converges to the real attractor A.*

This assertion can be easily verified using estimates for finite-difference and finite-element approximations of Roger Temam [20]. This assertion answers the first question. Now we consider the convergence of the numerical method with certain properties for an arbitrary given well posed initial-boundary problem to the real solution of the problem. We are using the following theorem by Roger Temam [20], page 281:

**Theorem 02.** *Let space dimension is 3 and we have some numerical approximation for Navier-Stokes equations with some conditions on which its stable and discrete elements are h for space and k for time. Then there exists the sequence for h and k* → 0*, that:*

$$
\begin{array}{c}
\mathbf{V}\_{h} \to \mathbf{V} \\
 P\_{h} \to P
\end{array};
\tag{10}
$$

One immediately derives:

*lk* =

Hence *Y* for a developed turbulent regime can be given as:

 *ν*3*L* **V**<sup>3</sup>

1/4 = *ν*3*L*<sup>4</sup> **V**3*L*<sup>3</sup>

*Y* =

is a phase space trajectories attraction manifold we have the following

*Navier-Stokes equations up to the precision order of the numerical method.*

*Sh* · *<sup>∂</sup>***<sup>V</sup>**

**2.2. Numerical method description and consistent mesh adaptation**

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> (**<sup>V</sup>** · ∇) **<sup>V</sup>** <sup>+</sup> *Eu* · ∇*<sup>P</sup>* <sup>=</sup> <sup>1</sup>

**<sup>V</sup>** <sup>=</sup> *<sup>f</sup>*(*<sup>x</sup>*), on *<sup>∂</sup>*Ω<sup>0</sup> <sup>×</sup> (0, *<sup>t</sup>*), **<sup>V</sup>** <sup>=</sup> 0, on *<sup>∂</sup>*Ω<sup>1</sup> <sup>×</sup> (0, *<sup>t</sup>*), *<sup>∂</sup>***<sup>V</sup>**

 *L <sup>L</sup>* · *<sup>R</sup>*3/4 <sup>3</sup>

that corresponds with the maximum number of degrees of freedom for Direct Numerical Simulation. For smaller values of Reynolds (R) number where the flow is in transitional regime the value for (14) is not perfectly true. For the purpose of finding *Y* in this case we are using a spatially constructed numerical procedure which is given below. For now we assume that appropriate values of *h* and *k* for the given *R* are defined and conditions 1 and 2 are true for some selected numerical solution procedure. Taking into account that the attractor

**Assertion 02.** *For a well-posed initial boundary value problem the numerical solution method with correct values of h and k and true Conditions 1 and 2 approximates an attractor of 3D evolutionary*

So we can construct the numerical method with the above mentioned properties and it can be used to analyze laminar-turbulent transition as a nonlinear dynamic system. However we should point out that bifurcation parameters do depend on the numerical method. So the values of these parameters can vary from one method to another for a given bifurcation. We are considering advection initial boundary problems with no external force acting for dimensionless Navier-Stokes equations. The general problem can be described as: Let Ω be a bounded domain with local Lipschiz-continuous boundaries *∂*Ω*i*. One must find velocity vector-function **<sup>V</sup>** : <sup>Ω</sup> <sup>×</sup> [0, *<sup>t</sup>*] <sup>→</sup> **<sup>R</sup>**<sup>3</sup> and scalar pressure function *<sup>P</sup>* : <sup>Ω</sup> <sup>×</sup> [0,*t*] <sup>→</sup> **<sup>R</sup>** such as:

∇ · **V** = 0 in *Q*;

**V**(*x*, 0) = **V**0(*x*, *y*, *z*), in Ω with ∇ · **V**0(*x*, *y*, *z*) = 0.

Similarity criteria are constructed using characteristic macroscopic scales: *Sh* = *f L*/*V* - Strouhal number; *Eu* = 2*P*/(*ρV*2) - Euler number; *R* = *VL*/*ν* - Reynolds number; *L* - macroscopic characteristic scale, *f* - frequency, *V* - characteristic velocity, *P*<sup>0</sup> - reference pressure, *ρ* = *const* - fluid density, *ν* - fluid kinematic viscosity. Various boundary conditions are given as: *∂*Ω<sup>0</sup> - inflow condition, *∂*Ω1- solid wall condition, *∂*Ω<sup>2</sup> - outflow condition and the problem is initialized with initial conditions. Most important similarity criterion is Reynolds number for laminar-turbulent transition. The rest criteria can be used to scale some real problems using *π* -theorem [21] and are not used since we are interested only in nonlinear

*<sup>R</sup>* <sup>∇</sup>2**<sup>V</sup>** in *<sup>Q</sup>* <sup>=</sup> <sup>Ω</sup> <sup>×</sup> (0, *<sup>t</sup>*);

*<sup>∂</sup><sup>n</sup>* = 0, on *∂*Ω<sup>2</sup> × (0, *t*);

(15)

1/4

<sup>=</sup> *<sup>L</sup>* · *<sup>R</sup>*<sup>−</sup>3/4, (13)

FSM Scenarios of Laminar-Turbulent Transition in Incompressible Fluids 255

= *R*9/4, (14)

*strong in L*2(*Q*)*, week in L*∞(*Q*, 0, <sup>Ω</sup>)*; Q* <sup>=</sup> <sup>Ω</sup> <sup>×</sup> [0, *<sup>t</sup>*]*.*

Here: **V**,*P* - velocity vector-function and pressure scalar-function that correspond to a solution for an initial-boundary value problem, **V***h*, *Ph* - velocity vector-function and pressure scalar-function on sequence *h*. *C* is an arbitrary constant. Here we consider one of possible solutions in attractor *A*, since uniqueness of a solution for a problem is not proven in 3D case. Prove of this theorem is given in [20], p.282, we are using only some conditions for the theorem validity:

**Condition 1.** *If a function x* �→ **<sup>V</sup>**(*x*); *<sup>x</sup>* <sup>⊂</sup> **<sup>R</sup>**<sup>3</sup> *has the property* ∇ · **<sup>V</sup>** <sup>=</sup> <sup>0</sup> *, then a discrete function h* �→ **V***<sup>h</sup> must maintain the property* ∇*<sup>h</sup>* · **V***<sup>h</sup>* = 0*.*

**Condition 2.** *For any sequence* **V***<sup>h</sup> the condition:*

$$\sup \left[ \sum\_{h} |\mathbf{V}\_{h+1} - \mathbf{V}\_{h}|^{k+1} - \sum\_{h} |\mathbf{V}\_{h+1} - \mathbf{V}\_{h}|^{k} \right] \le 0,\tag{11}$$

*has to hold.*

Appropriate approximations and numerical procedures where used in [16, 20] . A screw-symmetric numerical operator ∇*<sup>h</sup>* was used for (**V** · ∇) **V** , linear high order operators where used for linear diffusion parts, etc. Time integration was conducted by implicit Crank-Nicolson method for linear part of Navier-Stokes system and explicit second order for nonlinear part. Pressure correction was used so that Condition 1 is true on each time step. However the theorem is true only for *h* and *k* → 0 [19], i.e. when the numerical system phase space dimension goes to ∞ . It is obvious that *h* and *k* are finite for a real numerical method that is applied on computers. Since the attractor dimension is limited from above (9) and number of degrees of freedom is very large but finite (4), one can stipulate that *h* and *k* can be finite if condition (4) and hence (9) are satisfied. It is hard to give a precise estimate but one can approximately evaluate those values by using Kolmogorov turbulent theory for invariant scales.

The least motion scale of turbulence can be defined as *lk* = *k*−<sup>1</sup> *<sup>k</sup>* , where *kk* is given by (2). Since is bounded (7), the dissipation rate for a developed turbulent flow is given approximately as:

$$
\epsilon \propto \frac{\mathbf{V}^3}{L} \tag{12}
$$

One immediately derives:

4 Will-be-set-by-IN-TECH

consider the convergence of the numerical method with certain properties for an arbitrary given well posed initial-boundary problem to the real solution of the problem. We are using

**Theorem 02.** *Let space dimension is 3 and we have some numerical approximation for Navier-Stokes equations with some conditions on which its stable and discrete elements are h for space and k for time.*

> **V***<sup>h</sup>* → **V** *Ph* → *P*

Here: **V**,*P* - velocity vector-function and pressure scalar-function that correspond to a solution for an initial-boundary value problem, **V***h*, *Ph* - velocity vector-function and pressure scalar-function on sequence *h*. *C* is an arbitrary constant. Here we consider one of possible solutions in attractor *A*, since uniqueness of a solution for a problem is not proven in 3D case. Prove of this theorem is given in [20], p.282, we are using only some conditions for the theorem

**Condition 1.** *If a function x* �→ **<sup>V</sup>**(*x*); *<sup>x</sup>* <sup>⊂</sup> **<sup>R</sup>**<sup>3</sup> *has the property* ∇ · **<sup>V</sup>** <sup>=</sup> <sup>0</sup> *, then a discrete function*

*<sup>k</sup>*+<sup>1</sup> <sup>−</sup> ∑ *h*

Appropriate approximations and numerical procedures where used in [16, 20] . A screw-symmetric numerical operator ∇*<sup>h</sup>* was used for (**V** · ∇) **V** , linear high order operators where used for linear diffusion parts, etc. Time integration was conducted by implicit Crank-Nicolson method for linear part of Navier-Stokes system and explicit second order for nonlinear part. Pressure correction was used so that Condition 1 is true on each time step. However the theorem is true only for *h* and *k* → 0 [19], i.e. when the numerical system phase space dimension goes to ∞ . It is obvious that *h* and *k* are finite for a real numerical method that is applied on computers. Since the attractor dimension is limited from above (9) and number of degrees of freedom is very large but finite (4), one can stipulate that *h* and *k* can be finite if condition (4) and hence (9) are satisfied. It is hard to give a precise estimate but one can approximately evaluate those values by using Kolmogorov turbulent theory for invariant

is bounded (7), the dissipation rate for a developed turbulent flow is given approximately

<sup>∝</sup> **<sup>V</sup>**<sup>3</sup>


*k* 

≤ 0, (11)

*<sup>k</sup>* , where *kk* is given by (2). Since

*<sup>L</sup>* , (12)

; (10)

the following theorem by Roger Temam [20], page 281:

*Then there exists the sequence for h and k* → 0*, that:*

*strong in L*2(*Q*)*, week in L*∞(*Q*, 0, <sup>Ω</sup>)*; Q* <sup>=</sup> <sup>Ω</sup> <sup>×</sup> [0, *<sup>t</sup>*]*.*

*h* �→ **V***<sup>h</sup> must maintain the property* ∇*<sup>h</sup>* · **V***<sup>h</sup>* = 0*.* **Condition 2.** *For any sequence* **V***<sup>h</sup> the condition:*

> sup ∑ *h*


The least motion scale of turbulence can be defined as *lk* = *k*−<sup>1</sup>

validity:

*has to hold.*

scales.

as:

$$d\_k = \left(\frac{\upsilon^3 L}{\mathbf{V}^3}\right)^{1/4} = \left(\frac{\upsilon^3 L^4}{\mathbf{V}^3 L^3}\right)^{1/4} = L \cdot \mathbb{R}^{-3/4} \text{.}\tag{13}$$

Hence *Y* for a developed turbulent regime can be given as:

$$Y = \left(\frac{L}{L \cdot R^{3/4}}\right)^3 = R^{9/4},\tag{14}$$

that corresponds with the maximum number of degrees of freedom for Direct Numerical Simulation. For smaller values of Reynolds (R) number where the flow is in transitional regime the value for (14) is not perfectly true. For the purpose of finding *Y* in this case we are using a spatially constructed numerical procedure which is given below. For now we assume that appropriate values of *h* and *k* for the given *R* are defined and conditions 1 and 2 are true for some selected numerical solution procedure. Taking into account that the attractor is a phase space trajectories attraction manifold we have the following

**Assertion 02.** *For a well-posed initial boundary value problem the numerical solution method with correct values of h and k and true Conditions 1 and 2 approximates an attractor of 3D evolutionary Navier-Stokes equations up to the precision order of the numerical method.*

### **2.2. Numerical method description and consistent mesh adaptation**

So we can construct the numerical method with the above mentioned properties and it can be used to analyze laminar-turbulent transition as a nonlinear dynamic system. However we should point out that bifurcation parameters do depend on the numerical method. So the values of these parameters can vary from one method to another for a given bifurcation. We are considering advection initial boundary problems with no external force acting for dimensionless Navier-Stokes equations. The general problem can be described as: Let Ω be a bounded domain with local Lipschiz-continuous boundaries *∂*Ω*i*. One must find velocity vector-function **<sup>V</sup>** : <sup>Ω</sup> <sup>×</sup> [0, *<sup>t</sup>*] <sup>→</sup> **<sup>R</sup>**<sup>3</sup> and scalar pressure function *<sup>P</sup>* : <sup>Ω</sup> <sup>×</sup> [0,*t*] <sup>→</sup> **<sup>R</sup>** such as:

$$\begin{aligned} Sh \cdot \frac{\partial \mathbf{V}}{\partial t} + (\mathbf{V} \cdot \nabla) \, \mathbf{V} + Eu \cdot \nabla P &= \frac{1}{\mathbb{R}} \nabla^2 \mathbf{V} \text{ in } Q = \Omega \times (0, t); \\ \nabla \cdot \mathbf{V} &= 0 \text{ in } Q; \\ \mathbf{V} = f(\vec{x}), \text{ on } \partial \Omega\_0 \times (0, t), \mathbf{V} = 0, \text{ on } \partial \Omega\_1 \times (0, t), \frac{\partial \mathbf{V}}{\partial \vec{n}} = 0, \text{ on } \partial \Omega\_2 \times (0, t); \\ \mathbf{V}(\vec{x}, 0) &= \mathbf{V}\_0(x, y, z), \text{ in } \Omega \text{ with } \nabla \cdot \mathbf{V}\_0(x, y, z) = 0. \end{aligned} \tag{15}$$

Similarity criteria are constructed using characteristic macroscopic scales: *Sh* = *f L*/*V* - Strouhal number; *Eu* = 2*P*/(*ρV*2) - Euler number; *R* = *VL*/*ν* - Reynolds number; *L* - macroscopic characteristic scale, *f* - frequency, *V* - characteristic velocity, *P*<sup>0</sup> - reference pressure, *ρ* = *const* - fluid density, *ν* - fluid kinematic viscosity. Various boundary conditions are given as: *∂*Ω<sup>0</sup> - inflow condition, *∂*Ω1- solid wall condition, *∂*Ω<sup>2</sup> - outflow condition and the problem is initialized with initial conditions. Most important similarity criterion is Reynolds number for laminar-turbulent transition. The rest criteria can be used to scale some real problems using *π* -theorem [21] and are not used since we are interested only in nonlinear dynamics of equations with no regard to a real specific problem. In order to determine necessary number of discrete elements one should take into account two factors; the first is to use the upper bound (14) as a start and derive *�* from the averaged equations and the second is to use modified wave number analysis to calculate how many elements are needed to represent certain wave number and thus necessary number of degrease of freedom. Using this method, see [22] one can numerically determine anisotropic element density in Ω and optimize this number. This approach also satisfies all conditions of theorems and assertions described above. Since the rate of energy dissipation is a function of Reynolds number and particular initial-boundary conditions we are going to compute required *Y* using Reynolds averaging. Let us rewrite momentum equation in (15) in coordinates as:

$$\frac{\partial V\_{l}}{\partial t} + \sum\_{j} \left( V\_{j} \frac{\partial V\_{l}}{\partial \mathbf{x}\_{j}} \right) + \frac{1}{\rho} \frac{\partial P}{\partial \mathbf{x}\_{l}} = \nu \frac{\partial^{2} V\_{l}}{\partial \mathbf{x}\_{l}^{2}}; \text{i.j.} \tag{16}$$

and introduce averaging:

$$\mathbf{V} = \overline{\mathbf{V}} + \mathbf{V}',\tag{17}$$

where summation is performed in accordance with (18). The value *τij*�

*∂*(*V*� *i* )2 *∂xj*

averaging rules (18) as:

index *i*:

trustworthy.

*∂k <sup>∂</sup><sup>t</sup>* <sup>+</sup> ∑ *j*

1 2

*∂k ∂xj*

*Vj* = − ∑ *j*

and perturbation rate of dissipation is described by:

*�* = ∑ *i*,*j*

� *∂V*� *i ∂xj τ*� *ij*� <sup>=</sup> *<sup>ν</sup>* <sup>2</sup> ∑ *i*,*j*

the following correlation to define the Kolmogorov wave number (2):

⎛ ⎝ *R*2 <sup>2</sup> ∑ *i*,*j*

*kk* =

where perturbation kinetic energy is written as:

� *∂ ∂xj* � *ρ*−1*P*� *V*� *<sup>j</sup>* + ∑ *i* � *V*�<sup>2</sup> *<sup>i</sup> V*� *j* � ��

*<sup>k</sup>* <sup>=</sup> <sup>1</sup> <sup>2</sup> ∑ *i* (*V*�

⎛ ⎝ (*V*� *i* )2 *<sup>∂</sup><sup>t</sup>* <sup>+</sup> ∑ *j*

infinitesimal order to the other parts of (22) and then the latter can be simplified by using

*Vj* ⎞ ⎠ + 1 *ρ ∂P*� *∂x V*� *<sup>i</sup>* + 1 <sup>2</sup> ∑ *j*

+ ∑ *j*

One can get the equation for perturbation energy balance applying a summation of (23) by

� *<sup>∂</sup>V*� *i ∂xj τ*� *ij* + *V*� *<sup>j</sup> V*� *i ∂V*� *i ∂xj*

> � *∂V*� *i ∂xj* + *∂V*� *j ∂xi*

> > �<sup>2</sup> ⎞ ⎠

1/4

The latter expression (26) is the exact value of (7) for a given Reynolds number. We are using

� *<sup>∂</sup>V*� *i ∂xj* + *∂V*� *j ∂xi*

where the sum expression is numerically calculated by a test simulation with the maximum possible number of discrete elements, defined by (14). After that we are constructing the isolines of *NY*, calculated by (4), where the dissipation wave number is applied through (27) and *N* is defined by modified wave number analysis, described bellow. Then the mesh that we are using is adopted to satisfy all calculated values of *Y*. Only after all these procedures we can say that dynamic nonlinear analysis results we obtained are true and can be considered

In order to solve an initial-value problem for Navier-stokes numerically we are introducing the following semi discrete scheme based on the fractional step method using high order TVD

*Vi*

+

*∂V*�<sup>2</sup> *<sup>i</sup> V*� *j ∂xj*

FSM Scenarios of Laminar-Turbulent Transition in Incompressible Fluids 257

� = 0.

+ ∑ *i*,*j*

�<sup>2</sup>

� *Tij ∂V*� *i ∂xj*

*<sup>i</sup>* )2, (25)

�

. (26)

, (27)

− *�*, (24)

� is the second

(23)

thus:

$$
\overline{V\_{\dot{l}} + V\_{\dot{j}}} = \overline{V\_{\dot{l}}} + \overline{V\_{\dot{j}}};\\\overline{\partial\_a V} = \partial\_a \overline{V};\\\overline{\overline{V\_{\dot{l}}}} \overline{V\_{\dot{j}}} = \overline{V\_{\dot{l}}} \overline{V\_{\dot{j}}};\\\alpha = \{t; \mathbf{x}\}.\tag{18}$$

Here *V* - averaged and *V*' instantaneous functions of velocity vector-function. By applying (17) to (16) one gets:

$$\frac{\partial \overline{V\_i}}{\partial t} + \sum\_{j} \left( \overline{V\_j} \frac{\partial \overline{V\_i}}{\partial \mathbf{x}\_j} + \overline{V\_j'} \frac{\overline{\partial V\_i'}}{\partial \mathbf{x}\_j} \right) + \frac{1}{\rho} \frac{\partial \overline{P}}{\partial \mathbf{x}\_i} = \sum\_{j} \frac{\partial \overline{\tau\_{ij}}}{\partial \mathbf{x}\_j} \tag{19}$$

where second rank tensor *τij* = *ν <sup>∂</sup>Vi ∂xj* + *<sup>∂</sup>Vj ∂xi* corresponds to the Newtonian fluid. Multiplying (16) on *Vi* and applying (17) one gets:

$$
\overline{\frac{\partial V\_i}{\partial t}V\_i} + \sum\_j \left( \overline{V\_j \frac{\partial V\_i}{\partial \mathbf{x}\_j}V\_i} \right) + \frac{1}{\rho} \frac{\overline{\frac{\partial P}{\partial P}}V\_i} = \sum\_j \overline{\frac{\partial \tau\_{ij}}{\partial \mathbf{x}\_j}V\_i}. \tag{20}
$$

Multiplying (19) on *Vi*:

$$
\frac{\overline{\partial V\_i}}{\partial t} \overline{V\_i} + \sum\_j \left( \overline{V\_i V\_j} \frac{\partial \overline{V\_i}}{\partial \mathbf{x}\_j} \right) + \frac{1}{\rho} \frac{\partial \overline{F}}{\partial \mathbf{x}\_i} \overline{V\_i} = \sum\_j \left[ \left( \frac{\partial \overline{\mathbf{r}\_{ij}}}{\partial \mathbf{x}\_j} + \frac{\partial T\_{ij}}{\partial \mathbf{x}\_j} \right) \overline{V\_i} \right], \tag{21}
$$

one gets the stress equation, where *Tij* = −*Vi* � *Vj* � are the components of virtual (Reynolds) stress tensor. Subtraction of (21) from (20) reads:

$$
\overline{\frac{\partial V\_i'}{\partial t}V\_i'} + \sum\_j \left( \overline{V\_j \frac{\partial V\_i}{\partial \mathbf{x}\_j}} \overline{V\_i} - \overline{V\_j} \frac{\partial \overline{V\_i}}{\partial \mathbf{x}\_j} \overline{V\_i} \right) + \frac{1}{\rho} \frac{\partial \overline{P'}}{\partial \mathbf{x}\_i} \overline{V\_i'} = \sum\_j \left( \frac{\partial \overline{\tau\_{ij}'}}{\partial \mathbf{x}\_j} \overline{V\_i'} - \frac{\partial T\_{ij}}{\partial \mathbf{x}\_j} \overline{V\_i} \right), \tag{22}
$$

where summation is performed in accordance with (18). The value *τij*� *Vi* � is the second infinitesimal order to the other parts of (22) and then the latter can be simplified by using averaging rules (18) as:

$$\begin{split} \frac{1}{2} \left( \frac{\overline{(V\_{i}')^2}}{\partial t} + \sum\_{j} \frac{\partial \overline{(V\_{i}')^2}}{\partial \mathbf{x}\_{j}} \overline{V\_{j}} \right) + \frac{1}{\rho} \frac{\overline{\partial P'}}{\partial \mathbf{x}} \overline{V\_{i}'} + \frac{1}{2} \sum\_{j} \frac{\partial \overline{V\_{i}'^2} \overline{V\_{j}'}}{\partial \mathbf{x}\_{j}} + \\ + \sum\_{j} \left( \overline{\frac{\partial V\_{i}'}{\partial \mathbf{x}\_{j}} \mathbf{r}\_{ij}'} + \overline{V\_{j}'} \overline{V\_{i}'} \frac{\partial \overline{V\_{i}'}}{\partial \mathbf{x}\_{j}} \right) = 0. \end{split} \tag{23}$$

One can get the equation for perturbation energy balance applying a summation of (23) by index *i*:

$$\frac{\partial \mathbf{k}}{\partial t} + \sum\_{j} \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_{j}} \overline{V\_{j}} = -\sum\_{j} \left( \frac{\partial}{\partial \mathbf{x}\_{j}} \left( \rho^{-1} \overline{P'V\_{j}'} + \sum\_{i} \left[ \overline{V\_{i}'^{2}V\_{j}'} \right] \right) \right) + \sum\_{i,j} \left( \overline{T\_{ij}} \frac{\partial \overline{V\_{i}'}}{\partial \mathbf{x}\_{j}} \right) - \varepsilon\_{\prime} \tag{24}$$

where perturbation kinetic energy is written as:

6 Will-be-set-by-IN-TECH

dynamics of equations with no regard to a real specific problem. In order to determine necessary number of discrete elements one should take into account two factors; the first is to use the upper bound (14) as a start and derive *�* from the averaged equations and the second is to use modified wave number analysis to calculate how many elements are needed to represent certain wave number and thus necessary number of degrease of freedom. Using this method, see [22] one can numerically determine anisotropic element density in Ω and optimize this number. This approach also satisfies all conditions of theorems and assertions described above. Since the rate of energy dissipation is a function of Reynolds number and particular initial-boundary conditions we are going to compute required *Y* using Reynolds

**V** = **V** + **V**�

Here *V* - averaged and *V*' instantaneous functions of velocity vector-function. By applying

+ *<sup>∂</sup>Vj ∂xi* 

 + 1 *ρ ∂P ∂xi*

*Vi* = ∑ *j*

� *Vj*

+ *Vj* � *<sup>∂</sup>Vi* � *∂xj*

 *<sup>∂</sup>Vi ∂xj*

*Vi* + *Vj* = *Vi* + *Vj*; *∂αV* = *∂αV*; *ViVj* = *ViVj*; *α* = {*t*; *x*}. (18)

= ∑ *j*

> *∂τij ∂xj*

*∂τij*� *∂xj Vi*

*Vi* = ∑ *j*

> *∂τij ∂xj* + *∂Tij ∂xj*

*∂τij ∂xj*

; *i*, *j* = {1, 2, 3} (16)

, (17)

corresponds to the Newtonian fluid.

 *Vi* 

� are the components of virtual (Reynolds)

� <sup>−</sup> *<sup>∂</sup>Tij ∂xj Vi* 

, (19)

*Vi*. (20)

, (21)

, (22)

averaging. Let us rewrite momentum equation in (15) in coordinates as:

 *Vj ∂Vi ∂xj* + <sup>1</sup> *ρ ∂P <sup>∂</sup>xi* <sup>=</sup> *<sup>ν</sup> <sup>∂</sup>*<sup>2</sup>*Vi ∂x*<sup>2</sup> *i*

 *Vj ∂Vi ∂xj*

*∂Vi <sup>∂</sup><sup>t</sup>* + ∑*<sup>j</sup>*

*∂Vi <sup>∂</sup><sup>t</sup>* <sup>+</sup> ∑ *j*

Multiplying (16) on *Vi* and applying (17) one gets:

*∂Vi ∂t*

*Vi* + ∑ *j*

one gets the stress equation, where *Tij* = −*Vi*

 *Vj ∂Vi ∂xj*

stress tensor. Subtraction of (21) from (20) reads:

*Vi* + ∑ *j*

 *ViVj ∂Vi ∂xj*

 *Vj ∂Vi ∂xj Vi* + 1 *ρ ∂P ∂xi*

> + 1 *ρ ∂P ∂xi*

*∂Vi ∂xj Vi* + 1 *ρ ∂P*� *∂xi Vi* � = ∑ *j*

*Vi* − *Vj*

where second rank tensor *τij* = *ν*

*∂Vi ∂t*

and introduce averaging:

(17) to (16) one gets:

Multiplying (19) on *Vi*:

*∂Vi* � *∂t Vi* � + ∑ *j*

thus:

$$k = \frac{1}{2} \sum\_{i} \overline{(V\_i')^2} \,\tag{25}$$

and perturbation rate of dissipation is described by:

$$\varepsilon = \sum\_{i,j} \left( \frac{\overline{\partial V\_i'}}{\partial \mathbf{x}\_j} \mathbf{r}\_{ij}' \right) = \frac{\nu}{2} \sum\_{i,j} \left( \frac{\partial V\_i'}{\partial \mathbf{x}\_j} + \frac{\partial V\_j'}{\partial \mathbf{x}\_i} \right)^2. \tag{26}$$

The latter expression (26) is the exact value of (7) for a given Reynolds number. We are using the following correlation to define the Kolmogorov wave number (2):

$$k\_k = \left(\frac{R^2}{2} \sum\_{i,j} \overline{\left(\frac{\partial V\_i'}{\partial x\_j} + \frac{\partial V\_j'}{\partial x\_i}\right)^2} \right)^{1/4} \text{.} \tag{27}$$

where the sum expression is numerically calculated by a test simulation with the maximum possible number of discrete elements, defined by (14). After that we are constructing the isolines of *NY*, calculated by (4), where the dissipation wave number is applied through (27) and *N* is defined by modified wave number analysis, described bellow. Then the mesh that we are using is adopted to satisfy all calculated values of *Y*. Only after all these procedures we can say that dynamic nonlinear analysis results we obtained are true and can be considered trustworthy.

In order to solve an initial-value problem for Navier-stokes numerically we are introducing the following semi discrete scheme based on the fractional step method using high order TVD Runge-Kutta forth order method with each step of RK as:

$$\begin{aligned} 1. \mathbf{V}' - \mathbf{V}^n &= -\Delta t (\mathbf{V}^n \cdot \nabla) \mathbf{V}^n; \\ 2. \mathbf{V}'^ - \mathbf{V}' &= \Delta t \Theta \nabla^2 \mathbf{V}'^ ('') \end{aligned}$$

$$\begin{aligned} \text{while } \nabla \cdot \mathbf{V}^{\beta+1} &\neq 0; \begin{cases} 3. \nabla^2 P = -\nabla \cdot \mathbf{V}^{\beta} / \Delta t \\ 4. \mathbf{V}^{\beta+1} = \mathbf{V}'' - \Delta t \nabla P \\ 5. \mathbf{V}^{n+1} = \mathbf{V}^{\beta+1} \end{cases} \end{aligned} \tag{28}$$

instead of integration

one immediately reads:

*Wi d dt***V***<sup>i</sup>* <sup>+</sup> <sup>∑</sup>*<sup>B</sup>*

*u*(*x*) = cos(*kx*) becomes:

write:

 *Si*

*j*=1 **V***f j* **V***f i n*�*j* 

*<sup>δ</sup>x*(cos(*kx*)) = <sup>1</sup>

*k*�

Δ*x* = −*i*

*N*−1 ∑ *j*=0

*cj* exp(

*i*2*πjk*

*<sup>N</sup>* ) = <sup>−</sup>*<sup>i</sup>*

where *s* is half length of the stencil in discrete space and *am* - are the interpolation coefficients of the considered numerical scheme. One can immediately see from (32) that for symmetric schemes only real part presents in the modified wave number since *am* = −*a*−*<sup>m</sup>* and *a*<sup>0</sup> = 0. That explains why symmetric approximation cannot be used for approximation of advection operator in (30). Numerical wave number analysis for WENO scheme was first made in [26] ], but here we use analytical analysis. By applying WENO weights, since WENO scheme is fifth order everywhere (even for discontinuous functions, see [24, 27] one can get interpolation

2Δ*x*

∇ · (**VV**) in (29) and (30) since ∇ · **V** = 0 on every timestep.

*f*(*U*) ·�*ndS* ≈

∑*B j*=1 **V***f <sup>j</sup>* · *n*�*<sup>j</sup>* 

*B* ∑ *j*=1

<sup>Δ</sup>*Sj* <sup>+</sup> /∇*hP*/ <sup>−</sup> <sup>Θ</sup> <sup>∑</sup>*<sup>B</sup>*

Here / / operator is treated by finite element method. Then discrete system (30) is applied in (28) by described numerical methods. Please note that (**V** · ∇)**V** = ∇ · (**VV**) + **V**∇ · **V** =

In order to complete the analysis of numerical scheme and answer all questions positively one must perform modified wave number analysis for used discrete schemes. We are considering first order PDE like on step 1, in (28). Let the differential operator *∂x* be approximated by central differences as *<sup>∂</sup><sup>x</sup>* <sup>=</sup> *<sup>δ</sup><sup>x</sup>* <sup>+</sup> <sup>O</sup>(Δ*x*2) using *<sup>N</sup>* discrete segments each with the length <sup>Δ</sup>*<sup>x</sup>* and let the function be *u*(*x*) = *ck* · exp(*ikx*). One immediately gets analytical solution as *∂xu*(*x*) = *ikck* exp(*ikx*) = *iku*(*x*). For the given central differences approximation the real part

cos(*k*(*<sup>x</sup>* <sup>+</sup> <sup>Δ</sup>*x*)) <sup>−</sup> <sup>1</sup>

Here we can see that the difference of analytical and numerical solutions is in wave number *<sup>k</sup>* that changes to *<sup>k</sup>*� <sup>=</sup> sin(*kx*)/Δ*<sup>x</sup>* and if *<sup>k</sup>*Δ*<sup>x</sup>* << 1 then *<sup>k</sup>*� <sup>=</sup> *<sup>k</sup>* <sup>−</sup> *<sup>k</sup>*3Δ*x*<sup>2</sup> <sup>+</sup> ... . For small *<sup>k</sup>* the result of approximate solution (31) is close to the analytical. But when wave number *k* increases i.e. discretization scale *λ* decreases (*λ* = 2*π*/*k* and if Δ*x* = *λ*/*N* then *k*Δ*x* = 2*π*/*N*), the error grows. So we call *k*� a modified wave number by the numerical scheme. Thus by applying this analysis one can get the minimal undisturbed scale representation of the given numerical discrete scheme. This allows us getting minimal number of finite volumes necessary for the given *k* say calculated *kk*. In general for the given complex function one can

2Δ*x*

= −

*m*=*s* ∑ *m*=−*s*

*f*(*Uj*) ·�*n*Δ*Sj*

*<sup>j</sup>*=<sup>1</sup> [∇*h***V**]

FSM Scenarios of Laminar-Turbulent Transition in Incompressible Fluids 259

*f*

cos(*k*(*x* − Δ*x*)) =

sin(*kx*)

*am* exp(*imk*Δ*x*), (32)

sin(*kx*) Δ*x*

*<sup>j</sup>* ·�*n*Δ*Sj* = 0.

(30)

(31)

Δ*Sj* = 0,

Here **V***<sup>n</sup>* and **V***n*+<sup>1</sup> are the previous timestep and next timestep values of the velocity vector function; *P* is the pressure, Δ*t* - timestep for the given Runge-Kutta stage and Θ is the diffusion parameter and for a forced advection problems equals *R*−1. All other superscripts on **V** are intermediate values of velocity vector function inside the stage. More details on this numerical procedure can be found in [9, 23, 24]. Each step here is just shortly described. On step 1 in (28) advection equations are solved with the condition (11) satisfied. We are using fifth order WENO-type scheme. Several problems were solved for pure advection equation and Burgers equation (in 1D and 2D) before the final variant of WENO scheme was selected. Since the timestep is limited by the accuracy requirements the numerical scheme was explicit. On the second step of (28) diffusion equations are solved by using large stencil approximation with the 6-th order of accuracy. It is possible to apply implicit method here (in (28) in brackets) for natural advection problem described bellow, since values of Θ are of unity magnitude. But for forced advection where Θ = *R*−<sup>1</sup> (for *R* > 100) it is possible to apply explicit time method. Academician Belotserkovsky O.M. [25] suggested a physical interpretation of this fractional step method. Step 1 and 2 are calculating not solenoidal vector field that breaks the mass conservation equation. But if we apply the operator (∇×) on step 1 and 2 for both (28) and Navier-Stokes equations (15), we get the curl transport equations, since ∇×∇*P* = 0. Curl properties are correctly simulated though ∇ · **V** = 0 is not satisfied even for steps 1 and 2. The latter corrected by applying the pressure correction on the third step, where the Poisson equation is solved for pressure scalar function until the solenoidal criteria is met up to the machine accuracy. After that the velocity field is corrected and he velocity vector function is solenoidal on step 4 up to machine accuracy, so condition 1 for correct approximation is satisfied.

Spatial discretization is a combined finite volume for **V** and finite element for P discretization. Since we are focusing on fundamental problems (with simple geometry) the discrete elements are rectangular cuboids, thus it allows us adjusting its dimensions in accordance with the calculated values of *Y*. Finite elements are the same cuboids but variables are stored on vertexes rather than centers of mass as for **V**. So the described system of equations is rewritten for arbitrary convex element *i* with volume *Wi* and Ω = *<sup>N</sup> <sup>i</sup>*=<sup>1</sup> *Wi*; *Wi* ∩ *Wj* = for ∀*i* �= *j* as:

$$\begin{aligned} \label{eq:SDAC-1} \frac{\partial \iint\_{\mathcal{W}\_{l}} \nabla \cdot \mathbf{V} d\mathcal{W} = \frac{1}{\mathcal{W}\_{l}} \oint\_{S\_{l}} \mathbf{V}^{f} \cdot \vec{\text{n}} d\mathcal{S} = 0, \\ \frac{\partial}{\partial t} \iiint\_{\mathcal{W}\_{l}} \mathbf{V} d\mathcal{W} + \oint\_{\mathcal{S}\_{l}} [\mathbf{V}^{f} \mathbf{V}^{f} \cdot \vec{\text{n}}] d\mathcal{S} + \oint\_{\mathcal{S}\_{l}} /\mathbf{P}^{f} / \vec{\text{n}} d\mathcal{S} - \Theta \oint\_{\mathcal{S}\_{l}} [\nabla \mathbf{V}]^{f} d\mathbf{s} = 0, \end{aligned} \tag{29}$$

here *S* - is one of element side square, subscript *f* refers to the face value and �*n* is the unit vector on the side, total number of sides is B. Using discrete elements and applying summation instead of integration

8 Will-be-set-by-IN-TECH

1.**V**� <sup>−</sup> **<sup>V</sup>***<sup>n</sup>* <sup>=</sup> <sup>−</sup>Δ*t*(**V***<sup>n</sup>* · ∇)**V***n*;

5.**V***n*+<sup>1</sup> = **V***β*+<sup>1</sup>

Here **V***<sup>n</sup>* and **V***n*+<sup>1</sup> are the previous timestep and next timestep values of the velocity vector function; *P* is the pressure, Δ*t* - timestep for the given Runge-Kutta stage and Θ is the diffusion parameter and for a forced advection problems equals *R*−1. All other superscripts on **V** are intermediate values of velocity vector function inside the stage. More details on this numerical procedure can be found in [9, 23, 24]. Each step here is just shortly described. On step 1 in (28) advection equations are solved with the condition (11) satisfied. We are using fifth order WENO-type scheme. Several problems were solved for pure advection equation and Burgers equation (in 1D and 2D) before the final variant of WENO scheme was selected. Since the timestep is limited by the accuracy requirements the numerical scheme was explicit. On the second step of (28) diffusion equations are solved by using large stencil approximation with the 6-th order of accuracy. It is possible to apply implicit method here (in (28) in brackets) for natural advection problem described bellow, since values of Θ are of unity magnitude. But for forced advection where Θ = *R*−<sup>1</sup> (for *R* > 100) it is possible to apply explicit time method. Academician Belotserkovsky O.M. [25] suggested a physical interpretation of this fractional step method. Step 1 and 2 are calculating not solenoidal vector field that breaks the mass conservation equation. But if we apply the operator (∇×) on step 1 and 2 for both (28) and Navier-Stokes equations (15), we get the curl transport equations, since ∇×∇*P* = 0. Curl properties are correctly simulated though ∇ · **V** = 0 is not satisfied even for steps 1 and 2. The latter corrected by applying the pressure correction on the third step, where the Poisson equation is solved for pressure scalar function until the solenoidal criteria is met up to the machine accuracy. After that the velocity field is corrected and he velocity vector function is solenoidal on step 4 up to machine accuracy, so condition 1 for correct approximation is

Spatial discretization is a combined finite volume for **V** and finite element for P discretization. Since we are focusing on fundamental problems (with simple geometry) the discrete elements are rectangular cuboids, thus it allows us adjusting its dimensions in accordance with the calculated values of *Y*. Finite elements are the same cuboids but variables are stored on vertexes rather than centers of mass as for **V**. So the described system of equations is rewritten

> *Wi Si*

here *S* - is one of element side square, subscript *f* refers to the face value and �*n* is the unit vector on the side, total number of sides is B. Using discrete elements and applying summation

*Si*

**<sup>V</sup>***<sup>f</sup>* ·�*ndS* <sup>=</sup> 0,

/*P<sup>f</sup>* /�*ndS* <sup>−</sup> <sup>Θ</sup>

*Si*

∇ · **<sup>V</sup>***dW* <sup>=</sup> <sup>1</sup>

[**V***<sup>f</sup>* **<sup>V</sup>***<sup>f</sup>* ·�*n*]*dS* <sup>+</sup>

*<sup>i</sup>*=<sup>1</sup> *Wi*; *Wi* ∩ *Wj* = for ∀*i* �= *j* as:

[∇**V**] *<sup>f</sup> ds* <sup>=</sup> 0,

(29)

for arbitrary convex element *i* with volume *Wi* and Ω = *<sup>N</sup>*

 *Wi*

*Si*

**V***dW* +

(��)

(28)

3.∇2*<sup>P</sup>* <sup>=</sup> −∇ · **<sup>V</sup>***β*/Δ*<sup>t</sup>* 4.**V***β*+<sup>1</sup> <sup>=</sup> **<sup>V</sup>**�� <sup>−</sup> <sup>Δ</sup>*t*∇*<sup>P</sup>*

2.**V**�� <sup>−</sup> **<sup>V</sup>**� <sup>=</sup> <sup>Δ</sup>*t*Θ∇2**V**�

Runge-Kutta forth order method with each step of RK as:

satisfied.

*∂ ∂t Wi* while ∇ · **<sup>V</sup>***β*+<sup>1</sup> �<sup>=</sup> 0;

$$\oint\_{S\_i} f(\mathcal{U}) \cdot \vec{n} dS \approx \sum\_{j=1}^{B} f(\mathcal{U}\_j) \cdot \vec{n} \Delta S\_j$$

one immediately reads:

$$\begin{aligned} \boldsymbol{\Sigma}\_{j=1}^{B} \left[ \mathbf{V}\_{j}^{f} \cdot \vec{n}\_{j}^{\*} \right] \boldsymbol{\Delta S}\_{j} &= \boldsymbol{0}, \\ \boldsymbol{\Delta W}\_{i}^{\frac{d}{d\mathbf{f}}} \mathbf{V}\_{i} + \boldsymbol{\Sigma}\_{j=1}^{B} \left[ \mathbf{V}\_{j}^{f} \mathbf{V}\_{i}^{f} \vec{n}\_{j}^{\*} \right] \boldsymbol{\Delta S}\_{j} + \boldsymbol{\slash \nabla\_{h} \mathbf{P} \boldsymbol{\prime}} \boldsymbol{\prime} - \boldsymbol{\Theta} \boldsymbol{\Sigma}\_{j=1}^{B} \left[ \boldsymbol{\nabla}\_{h} \mathbf{V} \right]\_{j}^{f} \cdot \vec{n} \boldsymbol{\Delta S}\_{j} &= \boldsymbol{0}. \end{aligned} \tag{30}$$

Here / / operator is treated by finite element method. Then discrete system (30) is applied in (28) by described numerical methods. Please note that (**V** · ∇)**V** = ∇ · (**VV**) + **V**∇ · **V** = ∇ · (**VV**) in (29) and (30) since ∇ · **V** = 0 on every timestep.

In order to complete the analysis of numerical scheme and answer all questions positively one must perform modified wave number analysis for used discrete schemes. We are considering first order PDE like on step 1, in (28). Let the differential operator *∂x* be approximated by central differences as *<sup>∂</sup><sup>x</sup>* <sup>=</sup> *<sup>δ</sup><sup>x</sup>* <sup>+</sup> <sup>O</sup>(Δ*x*2) using *<sup>N</sup>* discrete segments each with the length <sup>Δ</sup>*<sup>x</sup>* and let the function be *u*(*x*) = *ck* · exp(*ikx*). One immediately gets analytical solution as *∂xu*(*x*) = *ikck* exp(*ikx*) = *iku*(*x*). For the given central differences approximation the real part *u*(*x*) = cos(*kx*) becomes:

$$\begin{split} \delta\_{\mathbf{x}}(\cos(k\mathbf{x})) &= \frac{1}{2\Delta\mathbf{x}} \cos(k(\mathbf{x} + \Delta\mathbf{x})) - \frac{1}{2\Delta\mathbf{x}} \cos(k(\mathbf{x} - \Delta\mathbf{x})) = \\ &= -\left(\frac{\sin(k\mathbf{x})}{\Delta\mathbf{x}}\right) \sin(k\mathbf{x}) \end{split} \tag{31}$$

Here we can see that the difference of analytical and numerical solutions is in wave number *<sup>k</sup>* that changes to *<sup>k</sup>*� <sup>=</sup> sin(*kx*)/Δ*<sup>x</sup>* and if *<sup>k</sup>*Δ*<sup>x</sup>* << 1 then *<sup>k</sup>*� <sup>=</sup> *<sup>k</sup>* <sup>−</sup> *<sup>k</sup>*3Δ*x*<sup>2</sup> <sup>+</sup> ... . For small *<sup>k</sup>* the result of approximate solution (31) is close to the analytical. But when wave number *k* increases i.e. discretization scale *λ* decreases (*λ* = 2*π*/*k* and if Δ*x* = *λ*/*N* then *k*Δ*x* = 2*π*/*N*), the error grows. So we call *k*� a modified wave number by the numerical scheme. Thus by applying this analysis one can get the minimal undisturbed scale representation of the given numerical discrete scheme. This allows us getting minimal number of finite volumes necessary for the given *k* say calculated *kk*. In general for the given complex function one can write:

$$k' \Delta \mathbf{x} = -i \sum\_{j=0}^{N-1} c\_j \exp(\frac{i2\pi jk}{N}) = -i \sum\_{m=-s}^{m=s} a\_m \exp(imk \Delta \mathbf{x})\_{\prime} \tag{32}$$

where *s* is half length of the stencil in discrete space and *am* - are the interpolation coefficients of the considered numerical scheme. One can immediately see from (32) that for symmetric schemes only real part presents in the modified wave number since *am* = −*a*−*<sup>m</sup>* and *a*<sup>0</sup> = 0. That explains why symmetric approximation cannot be used for approximation of advection operator in (30). Numerical wave number analysis for WENO scheme was first made in [26] ], but here we use analytical analysis. By applying WENO weights, since WENO scheme is fifth order everywhere (even for discontinuous functions, see [24, 27] one can get interpolation

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

coefficients:

$$a\_m^{WENO} = (-1/30; 1/4; 11/5; -1/2; -1/20; 0); \\ s = 3; \mathcal{O}(\Delta x^5). \tag{33}$$

 *<sup>V</sup>n*+<sup>1</sup> <sup>≤</sup> max

�*Vn*� ; *<sup>V</sup>n*−<sup>1</sup> ; ..;

*<sup>β</sup><sup>i</sup>* → ∞.

results about nonlinear dynamics in the problem.

benchmark results but bifurcation sequences as well.

with time step limited by:

if <sup>∀</sup>*β<sup>i</sup>* <sup>=</sup> 0, then *<sup>α</sup><sup>i</sup>*

transition.

Reynolds number.

*<sup>V</sup>n*−*m*+<sup>1</sup>

Δ*t* ≤ min *i*

 

> *αi βi*

The same idea is given in [30]. Proving the lemma is done by considering the combinations of Euler's steps. This lemma gives us constants *α<sup>i</sup>* and *β<sup>i</sup>* that were used to construct the forth order TVD Runge-Kutta method that fulfills the condition (11) and has a maximum stable timestep possible. One should point out that the standard RK4 method with constants (1/6; 1/3; 1/3; 1/6) is not TVD and can't be used for time integration. So it is shown that the presented numerical method is guaranteed to satisfy all given theorems, conditions and assertions and, hence, can be used to describe nonlinear dynamics of laminar-turbulent

The whole work described here took three and half years and 80% of the time was used on calculation. New CUDA technology [31] was applied lately and now all described numerical

**3. Laminar-turbulent transition for the flow over a backward facing step**

One of the best studied problems is the problem of the flow over a backward facing step. It was simulated by many different authors and has lots of benchmark results and even few

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

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

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

methods are calculated using NVIDIA GPUs that greatly accelerate the research.

Δ*t*(*Euler*); (36)

FSM Scenarios of Laminar-Turbulent Transition in Incompressible Fluids 261

By inserting (33) into (32) one can get modified wave number (real and imaginary parts) and derive minimal necessary numbers of elements to represent a given wave number, see fig.1. So one can see that for the 2-nd order central differences scheme at least 8 elements are needed

**Figure 1.** Modified wave number for different schemes to *k*Δ*x* ∈ [0, *π*](left) and minimal number of elements *N* for the correct given wave number representation (right). CD - central differences.

to describe a given wave number correctly, and 2 elements for the 6-th order scheme as well as for the WENO scheme. So now we know that the values that are found in accordance with (27) must be multiplied by two. Applying the same method for other parts of discretization one can point out that the given estimate for two elements is enough and here is omitted for the sake of brevity.

Time integration also requires some care since the condition (11) is in spatio-temporal condition and if CFD literature is known as TVD (TVB for equals sign) condition for nonlinear PDEs [28, 29]. Since WENO approximation is TVB, we are going to give a TVD condition (11) for time integration:

$$\frac{\partial V\_{\text{i}}}{\partial t} = L(V\_{\text{i}}) \,\tag{34}$$

where *L* is the given complex nonlinear TVB differential operator consisting of all steps (28). Simple TVD method for (34) is the Euler's method *Vn*+<sup>1</sup> = *V<sup>n</sup>* + Δ*tL*(*Vn*+*<sup>s</sup>* ). If (34) is the explicit method (i.e. *s* = 0) then the stability criterion is a Courant number Δ*t* ≤ *C*. One can give the following

**Lemma 01.** *If a direct Euler's method* (34) *is applied for a TVB or TVD spatial approximation operator for* Δ*t* ≤ *C and by applying m-stage Runge-Kutta scheme*

$$\begin{aligned} V^{n+1} = \sum\_{i=1}^{m} \left( \mathfrak{a}\_{i} V^{n+1-i} + \Delta t \beta\_{i} L(V^{n+1-i}) \right); \\ \sum\_{i=1}^{m} \mathfrak{a}\_{i} = 1, \forall \mathfrak{a}\_{i\prime} \beta\_{i} \ge 0, \text{on } L(V), \end{aligned} \tag{35}$$

*then the solution is stable in any norm*

$$\left\| \left| V^{n+1} \right| \right\| \le \max\left\{ \left\| \left| V^n \right| \right\| ; \left\| \left| V^{n-1} \right| \right\| ; \dots \left\| \left| V^{n-m+1} \right| \right\| \right\}$$

with time step limited by:

10 Will-be-set-by-IN-TECH

By inserting (33) into (32) one can get modified wave number (real and imaginary parts) and derive minimal necessary numbers of elements to represent a given wave number, see fig.1. So one can see that for the 2-nd order central differences scheme at least 8 elements are needed

**Figure 1.** Modified wave number for different schemes to *k*Δ*x* ∈ [0, *π*](left) and minimal number of elements *N* for the correct given wave number representation (right). CD - central differences.

to describe a given wave number correctly, and 2 elements for the 6-th order scheme as well as for the WENO scheme. So now we know that the values that are found in accordance with (27) must be multiplied by two. Applying the same method for other parts of discretization one can point out that the given estimate for two elements is enough and here is omitted for

Time integration also requires some care since the condition (11) is in spatio-temporal condition and if CFD literature is known as TVD (TVB for equals sign) condition for nonlinear PDEs [28, 29]. Since WENO approximation is TVB, we are going to give a TVD condition (11)

where *L* is the given complex nonlinear TVB differential operator consisting of all steps (28).

explicit method (i.e. *s* = 0) then the stability criterion is a Courant number Δ*t* ≤ *C*. One can

**Lemma 01.** *If a direct Euler's method* (34) *is applied for a TVB or TVD spatial approximation operator*

*αiV<sup>n</sup>*+1−*<sup>i</sup>* + Δ*tβiL*(*Vn*+1−*<sup>i</sup>*

*<sup>∂</sup><sup>t</sup>* <sup>=</sup> *<sup>L</sup>*(*Vi*), (34)

) ;

*<sup>i</sup>*=<sup>1</sup> *<sup>α</sup><sup>i</sup>* <sup>=</sup> 1, <sup>∀</sup>*αi*, *<sup>β</sup><sup>i</sup>* <sup>≥</sup> 0, on *<sup>L</sup>*(*V*), (35)

). If (34) is the

*∂Vi*

Simple TVD method for (34) is the Euler's method *Vn*+<sup>1</sup> = *V<sup>n</sup>* + Δ*tL*(*Vn*+*<sup>s</sup>*

*for* Δ*t* ≤ *C and by applying m-stage Runge-Kutta scheme*

*then the solution is stable in any norm*

*Vn*+<sup>1</sup> = ∑*<sup>m</sup>*

∑*<sup>m</sup>*

*i*=1 

*<sup>m</sup>* <sup>=</sup> (−1/30; 1/4; 11/5; <sup>−</sup>1/2; <sup>−</sup>1/20; 0);*<sup>s</sup>* <sup>=</sup> 3; <sup>O</sup>(Δ*x*5). (33)

coefficients:

the sake of brevity.

for time integration:

give the following

*aWENO*

$$
\Delta t \le \min\_i \frac{\alpha\_i}{\beta\_i} \Delta t (Euler); \tag{36}
$$

if <sup>∀</sup>*β<sup>i</sup>* <sup>=</sup> 0, then *<sup>α</sup><sup>i</sup> <sup>β</sup><sup>i</sup>* → ∞.

The same idea is given in [30]. Proving the lemma is done by considering the combinations of Euler's steps. This lemma gives us constants *α<sup>i</sup>* and *β<sup>i</sup>* that were used to construct the forth order TVD Runge-Kutta method that fulfills the condition (11) and has a maximum stable timestep possible. One should point out that the standard RK4 method with constants (1/6; 1/3; 1/3; 1/6) is not TVD and can't be used for time integration. So it is shown that the presented numerical method is guaranteed to satisfy all given theorems, conditions and assertions and, hence, can be used to describe nonlinear dynamics of laminar-turbulent transition.

The whole work described here took three and half years and 80% of the time was used on calculation. New CUDA technology [31] was applied lately and now all described numerical methods are calculated using NVIDIA GPUs that greatly accelerate the research.
