Soe, M., Vahala, G., Pavlo, P., Vahala, L., & Chen, H. (1998).Thermal lattice Boltzmann

Shan, X. (1997). Simulation of Rayleigh-Benard convection using a lattice Boltzmann method

Son, G. & Dhir, V.K. (1998). Numerical simulation of a single bubble during particle nucleate

Shan, X. & Chen, H. (1993). Lattice Boltzmann model for simulating flows with multiple

Swift, M.R., Orlandini, E., Osborn,W.R., Yeomeans, J.M. (1996). Lattice Boltzmann simulations of liquid-gas and binary fluid systems, Phys.Rev. E54.pp.5041~5052. Tomiyama, A., Sou, A.; Minagawa, H., & Sakaguchi, T. (1993). Numerical analysis of a single

Vahala, G., Pavlo, P., Vahala, L., & Martys, N. S. (1998). Thermal lattice-Boltzmann models

Vahala, L., Wah, D., Vahala, G., Carter, J., & Pavlo, P. (2000).Thermal lattice Boltzmann simulation for multispecies fluid equilibration[J]. Phys. Rev E, 62: 507-516. Wittke, D.D. & Chao, T.B. (1967). Collapse of vapour bubbles with translatory motion, J.

Yan, Y., & Li, W. (2006). Numerical modelling of a vapors bubble growth in uniformly

Zheng, H.W., Shu, C., & Chew, Y.T.(2006). A Lattice Boltzmann for multiphase flows with

superheated liquid, Int J. of Numerical Methods for Heat & Fluid Flow, Vol.16.

(TLBM) for compressible flows[J]. Int. J. Mod. Phys C, 9: 1247-1261.

phases and components, Phys. Rev. E47(3).pp.1815~1819.

bubble by VOF method, JSME Int J. Series B. Vol.36.No.1.

large density ratio, J. Comput. Phys. Vol.218. pp.353~371.

simulations of variable Prandtl number turbulent flows [J]. Phys. Rev E, 57(4): 4227-

boiling on a horizontal surface, Proceedings of 11th IHTC. Kyongiu, Korea, August

Vol. 41, No.2, pp. 165~189.

J.Comput.Phys.206:16-47.

Transfer 13, 657-666 .

pp.1119~1127.

Applied Physics.Vol.25.No.4.pp.293~500.

Tech..Rep.No.16, Cambridge, MA.

[J]. Phys.Rev.E, 55: 2780-2788

Heat Transfer. Vol.89. pp.17~24.

No.7. pp.764~778.

23-28, Vol.2. pp.533~538.

138.

4237.

boundary conditions of free surface flows with heat and mass transfer, Numerical Heat Transfer (An International Journal of Computation and Methodology), Part B,

rising through a quiescent liquid, Numerical Heat Transfer (An International

simulation of in compressible two-phase flows at high density ratio,

contact angle during nucleate pool boiling, Int, J. Heat and Mass Transfer, 50 : 127-

Basic tendency in computational fluid dynamics (CFD) consists in development of black box software for solving scientific and engineering problems. Numerical methods for solving nonlinear partial differential equations in black box manner should satisfy to the requirements:

a) the least number of the problem-dependent components

b) high computational efficiency

c) high parallelism

d) the least usage of the computer resources.

We continue with the 2D (*N* = 2) Navier–Stokes equations governing flow of a Newtonian, incompressible viscous fluid. Let Ω ∈ **R***N* be a bounded, connected domain with a piecewise smooth boundary *∂*Ω. Given a boundary data, the problem is to find a nondimensional velocity field and nondimensional pressure such that:

a) continuity equation

$$
\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0,
\tag{1}
$$

b) X-momentum

$$
\frac{\partial u}{\partial t} + \frac{\partial (u^2)}{\partial x} + \frac{\partial (vu)}{\partial y} = -\frac{\partial p}{\partial x} + \frac{1}{\text{Re}} \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right),
\tag{2}
$$

c) Y-momentum

$$
\frac{
\partial \boldsymbol{v}
}{
\partial t
} + \frac{
\partial (
\boldsymbol{u}
\boldsymbol{v}
)
}{
\partial \boldsymbol{x}
} + \frac{
\partial (
\boldsymbol{v}^2
)
}{
\partial \boldsymbol{y}
} = -\frac{
\partial \boldsymbol{p}
}{
\partial \boldsymbol{y}
} + \frac{1}{
\mathsf{Re}} \left(
\frac{
\partial^2 \boldsymbol{v}
}{
\partial \boldsymbol{x}^2
} + \frac{
\partial^2 \boldsymbol{v}
}{
\partial \boldsymbol{y}^2
}
\right) \,.
\tag{3}
$$

Reynold number Re is defined as

$$\text{Re} = \frac{\rho u\_s l\_s}{\mu}\_{\prime}$$

where *ρ* and *μ* are density and viscosity, respectively. Choice of the velocity scale *us* and geometric scale *ls* depends on the given problem.

where *β* is an exact solution. Choice of the preconditioner *Q* so

for Solving Navier–Stokes Equations on Structured Grids

guarantees geometric convergence rate of the Uzawa iterations �

intensive study Benzi et al. (2005).

benchmark and applied problems.

simplified Navier–Stokes equations.

a) X-momentum and mass conservation equations

*∂*(*u*2) *∂x* + *∂*(*vu*)

� 1

0

*∂u ∂t* +

⎧ ⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎩

until the convergent solution will be obtained.

b) continuity equation (1).

saddle point problems.

�

�*β* − *β*(*k*+1)

**2. Remarks on solvers for simplied Navier–Stokes equations**

�*I* − *Q*−1*BA*−1*BT*�

� � *qk*+1

Unfortunately the preconditioner *Q* is strongly problem-dependent component of the Uzawa algorithm. Additional problem arises at formulation of the boundary conditions for *Q*. As a rule, the preconditioner has some relaxation parameters and determination of their optimal values is sufficiently difficult problem. Now construction of the preconditioner is subject of

177 Convergence Acceleration of Iterative Algorithms

3) *development of new approaches for convergence acceleration* of iterative algorithms for solving

The main obstacles to be overcome are execution time requirements and the generation of computational grids in complex three-dimensional domains Benzi et al. (2005). Recently convergence acceleration technique based on original pressure decomposition has been proposed for structured grids Martynenko (2009). The technique can be used in black box software. The chapter represents detailed description of the approach and its application for

Limited characteristics of the first computers and absence of efficient numerical methods put difficulties for simulation of fluid flows based on the full Navier–Stokes equations. As a result, computational fluid dynamics started from simulation of the simplest flows described by the

As an example, we consider 2D laminar flow between parallel plates. Figure 1 represents geometry of the problem. Assuming that the pressure is not changed across the flow (*p*�

in case of *L* � 1), full Navier–Stokes equations can be reduced to the simplified form:

*u*(*t*, *x*, *y*) *dy* =

*∂y* = − *∂p*

� 1

0

Since the mass conservation equation follows from the continuity equation (1), system (7) must be solved first. Solution of system (7) gives velocity components *u* and pressure *p*. After that the continuity equation (1) is used for determination of *v*. The computations are repeated

Let us consider solution of system (7) in details. Assume that an uniform computational grid (*h* = *hx* = *hy*) is generated. Linearized finite-differenced equations with block unknowns

*∂x* + 1 Re � *∂*2*u ∂x*2 +

*u*(*t*, 0, *y*) *dy*

*∂*2*u ∂y*2 �

� *q* < 1

� �.

*y* = 0

, (7)

� �*β* − *β*(0)

Equations (1)–(3) can be rewritten in the operator form

$$\begin{cases} \mathcal{N}(\vec{V}) + \nabla P = F \\ \nabla \vec{V} = G \end{cases} \tag{4}$$

where N is nonlinear convection-diffusion operator, *F* and *G* are source terms, *V* and *P* are velocity and pressure, respectively. It is assumed that the operator N accounts boundary conditions. Note that 2D and 3D Navier–Stokes equations can be written as equation (4), where first and second equations abbreviate momentum and continuity equations. Linearized discrete Navier–Stokes equations can be written in the matrix form

$$
\begin{pmatrix} A \ B^T \\ B & 0 \end{pmatrix} \begin{pmatrix} \mathfrak{a} \\ \mathfrak{g} \end{pmatrix} = \begin{pmatrix} f \\ \mathfrak{g} \end{pmatrix} \tag{5}
$$

in which *α* and *β* represent the discrete velocity and discrete pressure, respectively. Here nonsymmetric *A* is a block diagonal matrix corresponding to the linearized discrete convection-diffusion operator N . The rectangular matrix *BT* represents the discrete gradient operator while *B* represents its adjoint, the divergence operator.

Large linear system of saddle point type (5) cannot be solved efficiently by standard methods of computational algebra. Due to their indefiniteness and poor spectral properties, such systems represent a significant challenge for solver developers Benzi et al. (2005).

Preconditioned Uzawa algorithm enjoys considerable popularity in computational fluid dynamics. The iterations for solving the saddle point system (5) are given by

$$\begin{cases} A\mathfrak{a}^{(k+1)} = -\mathsf{B}^{T}\mathfrak{f}^{(k)} + \mathsf{f} \\\\ \mathsf{Q}\mathfrak{f}^{(k+1)} = \mathsf{Q}\mathfrak{f}^{(k)} + \left(\mathsf{B}\mathfrak{a}^{(k+1)} - \mathsf{g}\right)' \end{cases} \tag{6}$$

where the matrix *Q* is some preconditioner.

Preconditioned Uzawa algorithm (6) defines the following way for improvement of the solvers for the Navier–Stokes equations:

1) *development of numerical methods for solving the boundary value problems.*

Uzawa iterations require fast numerical inversion of the matrices *A* and *Q*. Now algebraic and geometric multigrid methods are often used for the given purpose Wesseling (1991). Multigrid methods give algorithms that solve sparse linear system of N unknowns with *O*(N) computational complexity for large classes of problems. Variant of geometric multigrid methods with the problem-independent transfer operators for black box or/and parallel implementation is proposed in Martynenko (2006; 2010).

2) *development of preconditioning.*

Error vector in Uzawa iterations satisfies to the condition

$$||\mathfrak{B} - \mathfrak{B}^{(k+1)}|| \lesssim ||I - \mathbb{Q}^{-1}BA^{-1}B^{T}|| \cdot ||\mathfrak{B} - \mathfrak{B}^{(k)}||\_{\prime}$$

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

� N (*V* ) + ∇*P* = *F*

where N is nonlinear convection-diffusion operator, *F* and *G* are source terms, *V* and *P* are velocity and pressure, respectively. It is assumed that the operator N accounts boundary conditions. Note that 2D and 3D Navier–Stokes equations can be written as equation (4),

> � �*α β* � = � *f g* �

in which *α* and *β* represent the discrete velocity and discrete pressure, respectively. Here nonsymmetric *A* is a block diagonal matrix corresponding to the linearized discrete convection-diffusion operator N . The rectangular matrix *BT* represents the discrete gradient

Large linear system of saddle point type (5) cannot be solved efficiently by standard methods of computational algebra. Due to their indefiniteness and poor spectral properties, such

Preconditioned Uzawa algorithm enjoys considerable popularity in computational fluid

*Qβ*(*k*+1) = *Qβ*(*k*) + (*Bα*(*k*+1) − *g*)

Preconditioned Uzawa algorithm (6) defines the following way for improvement of the

Uzawa iterations require fast numerical inversion of the matrices *A* and *Q*. Now algebraic and geometric multigrid methods are often used for the given purpose Wesseling (1991). Multigrid methods give algorithms that solve sparse linear system of N unknowns with *O*(N) computational complexity for large classes of problems. Variant of geometric multigrid methods with the problem-independent transfer operators for black box or/and parallel

�*I* − *Q*−1*BA*−1*BT*�

� · � �*β* − *β*(*k*)

� �,

where first and second equations abbreviate momentum and continuity equations. Linearized discrete Navier–Stokes equations can be written in the matrix form

> �*A BT B* 0

systems represent a significant challenge for solver developers Benzi et al. (2005).

*Aα*(*k*+1) = −*BTβ*(*k*) + *f*

dynamics. The iterations for solving the saddle point system (5) are given by

1) *development of numerical methods for solving the boundary value problems.*

operator while *B* represents its adjoint, the divergence operator.

⎧ ⎨ ⎩

implementation is proposed in Martynenko (2006; 2010).

Error vector in Uzawa iterations satisfies to the condition

�*β* − *β*(*k*+1)

� � �

�

where the matrix *Q* is some preconditioner.

solvers for the Navier–Stokes equations:

2) *development of preconditioning.*

∇*V* = *G* , (4)

(5)

, (6)

Equations (1)–(3) can be rewritten in the operator form

where *β* is an exact solution. Choice of the preconditioner *Q* so

$$\left\| \left\| I - Q^{-1}BA^{-1}B^{T} \right\| \right\| \leqslant q < 1$$

guarantees geometric convergence rate of the Uzawa iterations

$$||\mathfrak{B} - \mathfrak{B}^{(k+1)}|| \le q^{k+1} ||\mathfrak{B} - \mathfrak{B}^{(0)}||.$$

Unfortunately the preconditioner *Q* is strongly problem-dependent component of the Uzawa algorithm. Additional problem arises at formulation of the boundary conditions for *Q*. As a rule, the preconditioner has some relaxation parameters and determination of their optimal values is sufficiently difficult problem. Now construction of the preconditioner is subject of intensive study Benzi et al. (2005).

3) *development of new approaches for convergence acceleration* of iterative algorithms for solving saddle point problems.

The main obstacles to be overcome are execution time requirements and the generation of computational grids in complex three-dimensional domains Benzi et al. (2005). Recently convergence acceleration technique based on original pressure decomposition has been proposed for structured grids Martynenko (2009). The technique can be used in black box software. The chapter represents detailed description of the approach and its application for benchmark and applied problems.

#### **2. Remarks on solvers for simplied Navier–Stokes equations**

Limited characteristics of the first computers and absence of efficient numerical methods put difficulties for simulation of fluid flows based on the full Navier–Stokes equations. As a result, computational fluid dynamics started from simulation of the simplest flows described by the simplified Navier–Stokes equations.

As an example, we consider 2D laminar flow between parallel plates. Figure 1 represents geometry of the problem. Assuming that the pressure is not changed across the flow (*p*� *y* = 0 in case of *L* � 1), full Navier–Stokes equations can be reduced to the simplified form: a) X-momentum and mass conservation equations

$$\begin{cases} \frac{\partial u}{\partial t} + \frac{\partial (u^2)}{\partial x} + \frac{\partial (vu)}{\partial y} = -\frac{\partial p}{\partial x} + \frac{1}{\text{Re}} \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) \\\\ \int\_0^1 u(t, x, y) \, dy = \int\_0^1 u(t, 0, y) \, dy \end{cases} \tag{7}$$

b) continuity equation (1).

Since the mass conservation equation follows from the continuity equation (1), system (7) must be solved first. Solution of system (7) gives velocity components *u* and pressure *p*. After that the continuity equation (1) is used for determination of *v*. The computations are repeated until the convergent solution will be obtained.

Let us consider solution of system (7) in details. Assume that an uniform computational grid (*h* = *hx* = *hy*) is generated. Linearized finite-differenced equations with block unknowns

is used for computation of pressure by the iterative method

for Solving Navier–Stokes Equations on Structured Grids

(*k*) *i* − *p*

(1)

guess can be given by perturbation of the first one, for example *p*

(*k*) *i* − *p*

*F*(*k*)

(0)

acceleration of the iterative methods intended for full Navier–Stokes equations.

*p*(*t*, *x*, *y*, *z*) = *px*(*t*, *x*) + *py*(*t*, *y*) + *pz*(*t*, *z*)

*pxyz*(*t*, *x*, *y*, *z*) = −*px*(*t*, *x*) − *py*(*t*, *y*) − *pz*

(*n*+1)

(*k*−1) *i*

179 Convergence Acceleration of Iterative Algorithms

where superscript *k* denotes the secant method iterations. Note that the approach requires

problems. Usually it is required several secant iterations to reduce error of the discrete mass

Note that in 2D case the system (9) can be solved by direct methods, i.e. without the secant iterations. However in 3D case the direct methods require unpractical computational efforts

As contrasted to the Uzawa algorithm (6), the method does not require some preconditioner(s), relaxation parameter(s), extra computer memory and has high convergence rate. Unfortunately, basic assumption *p* = *p*(*t*, *x*) does not allow apply the method directly for solving full Navier–Stokes equations (1)–(3). Accounting the attractive properties, the algorithm for solving the simplified Navier–Stokes equations can be used for convergence

Reduction of system (5) to the saddle point system with zero block of the least size is popular approach in CFD. For example, similar reduction based on special unknown ordering is used

In order to apply the abovementioned approach for solving full Navier–Stokes equations, it is necessary artificially extract «one-dimensional parts of pressure» from the pressure field. For the given purpose, let add and subtract items *px*(*t*, *x*), *py*(*t*, *y*) and *pz*(*t*, *z*) depending only on

−*px*(*t*, *x*) − *py*(*t*, *y*) − *pz*

where superscripts *x*, *y* and *z* denote dependence of the functions on the spatial variables. Let

*p*(*t*, *x*, *y*, *z*) = *px*(*t*, *x*) + *py*(*t*, *y*) + *pz*(*t*, *z*) + *pxyz*(*t*, *x*, *y*, *z*). (10)

*F*(*k*)

*i* , *k* = 1, 2, . . . ,

*i* . First starting guess can be obtained by extrapolation. For

*i* = 2*pi*−1 − *pi*−2 and compute *F*(0). Second starting

(1)

*i* , but the secant method is direct solver for linear

(*t*, *z*) + *p*(*t*, *x*, *y*, *z*)

(*t*, *z*) + *p*(*t*, *x*, *y*, *z*).

 ,

*i* = 1.001*p*

(0)

*i* . It gives *F*(1)

.

*i* − *F*(*k*−1) *i*

*p* (*k*+1) *i* = *p*

example, for uniform grid we obtain *p*

Function *F* depends almost linearly on *p*

in Vanka smoother Vanka (1986).

one spatial variable, i.e.

us introduce a new function

Finally the pressure can be represented as

(0) *i* and *p*

conservation equation down to roundoff error.

due to five-diagonal structure of the coefficient matrix.

**3. Principle of formal decomposition of pressure**

+

two starting guesses *p*

(a) Geometry of problem about the flow between parallel plates (b) Block ordering of unknowns

#### Fig. 1. Flow between parallel plates

ordering shown on Figure 1 are written as

$$\begin{cases} a\_j u\_{i j - 1}^{(n+1)} + b\_j u\_{ij}^{(n+1)} + c\_j u\_{ij + 1}^{(n+1)} = p\_i^{(n+1)} + d\_j \\\\ \sum\_{j = 1}^{\mathcal{N}\_y} u\_{ij}^{(n+1)} = \frac{1}{h} G\_0 \end{cases} \tag{8}$$

where

$$G\_0 = \int\limits\_{0}^{1} u(t, 0, y) \, dy$$

is the given inlet mass flow rate and superscript *n* denotes time layer. Missing the superscript (*n* + 1), the system (8) can be rewritten in the matrix form

$$
\begin{pmatrix} b\_1 \ c\_1 & \cdots & 1 \\ a\_2 \ b\_2 \ c\_2 & \cdots & 1 \\ a\_3 \ b\_3 \ c\_3 & \cdots & 1 \\ & a\_4 \ b\_4 & \cdots & 1 \\ & & \cdots & \cdots & \cdots & \cdots \\ 1 & 1 & 1 & 1 & \cdots & 0 \\ \end{pmatrix} \begin{pmatrix} u\_{i1} \\ u\_{i2} \\ u\_{i3} \\ u\_{i4} \\ \vdots \\ u\_i \end{pmatrix} = \begin{pmatrix} d\_1 \\ d\_2 \\ d\_3 \\ \vdots \\ d\_4 \\ h^{-1}G\_0 \end{pmatrix}. \tag{9}
$$

Comparison of systems (9) and (5) shows that solution of the simplified Navier–Stokes equations (7) also is reduced to solution of the saddle point system. The principal difference between systems (9) and (5) consists in size of zero block in the coefficient matrix. Since the zero block in system (9) has the least size 1 × 1 because of the pressure is independent on *y*, efficient iterative algorithms for solution of system (9) have been proposed and developed. The most promising of them is secant method Briley (1974), where error of the mass conservation equation

$$F(p\_i^{(n+1)}) = \sum\_{j=1}^{N\_y} u\_{ij}^{(n+1)} (p\_i^{(n+1)}) - \frac{1}{h} G\_0$$

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

(b) Block ordering of unknowns

(*n*+1) *i* + *dj*

, (8)

. (9)

*h G*0

(a) Geometry of problem about the flow between parallel plates

> *u*(*n*+1) *ij* + *cj*

*G*0 = � 1

*b*1 *c*1 ··· 1 *a*2 *b*2 *c*2 ··· 1 *a*3 *b*3 *c*3 ··· 1 *a*4 *b*4 ··· 1

................. 1111 ··· 0

0

is the given inlet mass flow rate and superscript *n* denotes time layer. Missing the superscript

⎛

*ui*1 *ui*2 *ui*3 *ui*4 ··· *pi*

⎞

⎛

*d*1 *d*2 *d*3 *d*4 ··· *h*−1*G*0 ⎞

⎟⎟⎟⎟⎟⎟⎠

⎜⎜⎜⎜⎜⎜⎝

⎟⎟⎟⎟⎟⎟⎠ =

⎜⎜⎜⎜⎜⎜⎝

Comparison of systems (9) and (5) shows that solution of the simplified Navier–Stokes equations (7) also is reduced to solution of the saddle point system. The principal difference between systems (9) and (5) consists in size of zero block in the coefficient matrix. Since the zero block in system (9) has the least size 1 × 1 because of the pressure is independent on *y*, efficient iterative algorithms for solution of system (9) have been proposed and developed. The most promising of them is secant method Briley (1974), where error of the mass

> *u*(*n*+1) *ij* (*p*

(*n*+1) *i* ) − 1

*h G*0

N*y* ∑ *j*=1 ⎞

⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠

*u*(*n*+1) *ij*+1 = *p*

*u*(*n*+1) *ij* = 1

*u*(*t*, 0, *y*) *dy*

N*y* ∑ *j*=1

Fig. 1. Flow between parallel plates

where

conservation equation

ordering shown on Figure 1 are written as

⎧ ⎪⎪⎪⎪⎨ *aj u*(*n*+1) *ij*−1 + *bj*

(*n* + 1), the system (8) can be rewritten in the matrix form

⎛

⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝

*F*(*p* (*n*+1) *i* ) =

⎪⎪⎪⎪⎩

is used for computation of pressure by the iterative method

$$p\_i^{(k+1)} = p\_i^{(k)} - \frac{p\_i^{(k)} - p\_i^{(k-1)}}{F\_i^{(k)} - F\_i^{(k-1)}} F\_i^{(k)} \, \qquad k = 1, 2, \ldots, n$$

where superscript *k* denotes the secant method iterations. Note that the approach requires two starting guesses *p* (0) *i* and *p* (1) *i* . First starting guess can be obtained by extrapolation. For example, for uniform grid we obtain *p* (0) *i* = 2*pi*−1 − *pi*−2 and compute *F*(0). Second starting guess can be given by perturbation of the first one, for example *p* (1) *i* = 1.001*p* (0) *i* . It gives *F*(1) . Function *F* depends almost linearly on *p* (*n*+1) *i* , but the secant method is direct solver for linear problems. Usually it is required several secant iterations to reduce error of the discrete mass conservation equation down to roundoff error.

Note that in 2D case the system (9) can be solved by direct methods, i.e. without the secant iterations. However in 3D case the direct methods require unpractical computational efforts due to five-diagonal structure of the coefficient matrix.

As contrasted to the Uzawa algorithm (6), the method does not require some preconditioner(s), relaxation parameter(s), extra computer memory and has high convergence rate. Unfortunately, basic assumption *p* = *p*(*t*, *x*) does not allow apply the method directly for solving full Navier–Stokes equations (1)–(3). Accounting the attractive properties, the algorithm for solving the simplified Navier–Stokes equations can be used for convergence acceleration of the iterative methods intended for full Navier–Stokes equations.

Reduction of system (5) to the saddle point system with zero block of the least size is popular approach in CFD. For example, similar reduction based on special unknown ordering is used in Vanka smoother Vanka (1986).

#### **3. Principle of formal decomposition of pressure**

In order to apply the abovementioned approach for solving full Navier–Stokes equations, it is necessary artificially extract «one-dimensional parts of pressure» from the pressure field. For the given purpose, let add and subtract items *px*(*t*, *x*), *py*(*t*, *y*) and *pz*(*t*, *z*) depending only on one spatial variable, i.e.

$$\begin{aligned} p(t, \mathbf{x}, \mathbf{y}, \mathbf{z}) &= p^{\mathbf{x}}(t, \mathbf{x}) + p^{\mathbf{y}}(t, \mathbf{y}) + p^{\mathbf{z}}(t, \mathbf{z}) \\ &+ \left( -p^{\mathbf{x}}(t, \mathbf{x}) - p^{\mathbf{y}}(t, \mathbf{y}) - p^{\mathbf{z}}(t, \mathbf{z}) + p(t, \mathbf{x}, \mathbf{y}, \mathbf{z}) \right), \end{aligned}$$

where superscripts *x*, *y* and *z* denote dependence of the functions on the spatial variables. Let us introduce a new function

$$p^{xyz}(t, \mathbf{x}, y, z) = -p^x(t, \mathbf{x}) - p^y(t, y) - p^z(t, z) + p(t, \mathbf{x}, y, z).$$

Finally the pressure can be represented as

$$p(t, \mathbf{x}, y, z) = p^x(t, \mathbf{x}) + p^y(t, y) + p^z(t, z) + p^{xyz}(t, \mathbf{x}, y, z). \tag{10}$$

Fig. 2. Driven cavity and location of the control volumes *V*1 and *V*2

*ht*

181 Convergence Acceleration of Iterative Algorithms

*ht*

*u*(*n*+1/2) *ij* − *u*(*n*)

*∂x* − *∂*(*vu*) *∂y* + 1 Re *∂*2*u ∂x*2 +

*u*(*t*, *x*, *y*) *dy* = 0 ,

It is easy to see that intermediate velocity field *V*(*n*+1/2) is independent on pressure.

This disadvantage can be compensated partially by the pressure decomposition (10). Application of the decomposition requires two mass conservation equations for 2D problems. Integration of the continuity equation (1) over the control volumes *V*1 and *V*2 shown on

1

0

*ht*

where *ht* is time semispacing, *V*(*n*+1/2) is intermediate velocity field and *n* is a time layer. Stage I consists in solution of the momentum equations without pressure gradients. For

*ij*

Stage II: Δ*p* = ∇*V*(*n*+1/2)

= −(*V*(*n*)

= −∇*p* ,

*ht*

,

*∂*2*u ∂y*2

(*n*)

*ij*

*v*(*t*, *x*, *y*) *dx* = 0 . (13)

∇)*V*(*n*) + Re−1Δ*V*(*n*) ,

= *ψij* , (11)

. (12)

Stage I: *V*(*n*+1/2) − *V*(*n*)

for Solving Navier–Stokes Equations on Structured Grids

Stage III: *V*(*n*+1) − *V*(*n*+1/2)

simplicity X-momentum can be written as

where *ψij* is the given function defined as

Figure 2 gives

*ψij* =

> 1

0

− *∂*(*u*2)

Representation (10) will be called a principle of formal decomposition of pressure. *Basic idea of the method consists in application of the efficient numerical methods developed for the simplified Navier–Stokes equations for determination of part of pressure* (*i.e. for px*(*t*, *x*) + *py*(*t*, *y*) + *pz*(*t*, *z*)). Fast computation of part of pressure results in reduction of total computational efforts needed for full Navier–Stokes equations.

In spite of simplicity of the representation (10), it is necessary to comment the principle of formal decomposition of pressure:

*Remark 1*. All items *px*(*t*, *x*), *py*(*t*, *y*), *pz*(*t*, *z*) and *pxyz*(*t*, *x*, *y*, *z*) have no physical meaning, but physical meaning has their sum. In follows, the items *px*(*t*, *x*), *py*(*t*, *y*) and *pz*(*t*, *z*) will be called as «one-dimensional components of the pressure», and *pxyz*(*t*, *x*, *y*, *z*) as «multidimensional component». The quotes «» will indicate absence of the physical meaning of the «pressure components».

*Remark 2*. In *N*-dimensional case (*N* = 2, 3) pressure is represented as sum of *N* + 1 «components», therefore the method requires *N* extra conditions for determination of the «one-dimensional components». The convergence acceleration technique uses *N* mass conservation equations as a priori information of physical nature.

*Remark 3*. In spite of representation of the pressure as sum of *N* + 1 «components», all momentum equations have only two «pressure» gradients. For example, for X-momentum we obtain

$$\begin{split} \frac{\partial p}{\partial \mathbf{x}} &= \frac{\partial}{\partial \mathbf{x}} \Big( p^{\mathbf{x}}(t, \mathbf{x}) + p^{y}(t, \mathbf{y}) + p^{z}(t, \mathbf{z}) + p^{xyz}(t, \mathbf{x}, \mathbf{y}, \mathbf{z}) \Big) \\ &= \frac{\partial p^{\mathbf{x}}}{\partial \mathbf{x}} + \frac{\partial p^{xyz}}{\partial \mathbf{x}}. \end{split}$$

*Remark 4*. Efficiency of the acceleration technique depends strongly on the flow nature. For directed fluid flows (for example, flows in nozzles, pipes etc.) gradient of one of «one-dimensional component of pressure» *px*(*t*, *x*), *py*(*t*, *y*) or *pz*(*t*, *z*) is dominant. In this case impressive reduction of computational work is expected as compared with traditional algorithms (i.e. *px*(*t*, *x*) = *py*(*t*, *y*) = *pz*(*t*, *z*) = 0). However for rotated flows (for example, flow in a driven cavity) the approach shows the least efficiency.

*Remark 5*. In 3D case the method will be more efficient than in 2D case.

*Remark 6*. Velocity components and corresponding «one-dimensional components» in equation (10) are computed only in coupled manner. Velocity components and «multidimensional component» *pxyz*(*t*, *x*, *y*, *z*) in equation (10) can be computed in decoupled (segregated) or coupled manner.

*Remark 7*. Gradients of the «one-dimensional components» can be obtained in analytical form for explicit schemes. Implicit schemes require formulation of an auxiliary problem for determination of gradients of the «one-dimensional components».

#### **4. Development of explicit schemes**

First, consider modification of the explicit schemes using well-known benchmark problem about rotated flow in a driven cavity (Figure 2). Let a staggered grid with grid spacing *hx* and *hy* has been generated. Classical three-stage splitting scheme is represented as

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

Representation (10) will be called a principle of formal decomposition of pressure. *Basic idea of the method consists in application of the efficient numerical methods developed for the simplified Navier–Stokes equations for determination of part of pressure* (*i.e. for px*(*t*, *x*) + *py*(*t*, *y*) + *pz*(*t*, *z*)). Fast computation of part of pressure results in reduction of total computational efforts needed

In spite of simplicity of the representation (10), it is necessary to comment the principle of

*Remark 1*. All items *px*(*t*, *x*), *py*(*t*, *y*), *pz*(*t*, *z*) and *pxyz*(*t*, *x*, *y*, *z*) have no physical meaning, but physical meaning has their sum. In follows, the items *px*(*t*, *x*), *py*(*t*, *y*) and *pz*(*t*, *z*) will be called as «one-dimensional components of the pressure», and *pxyz*(*t*, *x*, *y*, *z*) as «multidimensional component». The quotes «» will indicate absence of the physical meaning

*Remark 2*. In *N*-dimensional case (*N* = 2, 3) pressure is represented as sum of *N* + 1 «components», therefore the method requires *N* extra conditions for determination of the «one-dimensional components». The convergence acceleration technique uses *N* mass

*Remark 3*. In spite of representation of the pressure as sum of *N* + 1 «components», all momentum equations have only two «pressure» gradients. For example, for X-momentum

*Remark 4*. Efficiency of the acceleration technique depends strongly on the flow nature. For directed fluid flows (for example, flows in nozzles, pipes etc.) gradient of one of «one-dimensional component of pressure» *px*(*t*, *x*), *py*(*t*, *y*) or *pz*(*t*, *z*) is dominant. In this case impressive reduction of computational work is expected as compared with traditional algorithms (i.e. *px*(*t*, *x*) = *py*(*t*, *y*) = *pz*(*t*, *z*) = 0). However for rotated flows (for example,

*Remark 6*. Velocity components and corresponding «one-dimensional components» in equation (10) are computed only in coupled manner. Velocity components and «multidimensional component» *pxyz*(*t*, *x*, *y*, *z*) in equation (10) can be computed in decoupled

*Remark 7*. Gradients of the «one-dimensional components» can be obtained in analytical form for explicit schemes. Implicit schemes require formulation of an auxiliary problem for

First, consider modification of the explicit schemes using well-known benchmark problem about rotated flow in a driven cavity (Figure 2). Let a staggered grid with grid spacing *hx* and

*hy* has been generated. Classical three-stage splitting scheme is represented as

(*t*, *z*) + *pxyz*(*t*, *x*, *y*, *z*)

*px*(*t*, *x*) + *py*(*t*, *y*) + *pz*

conservation equations as a priori information of physical nature.

flow in a driven cavity) the approach shows the least efficiency. *Remark 5*. In 3D case the method will be more efficient than in 2D case.

determination of gradients of the «one-dimensional components».

for full Navier–Stokes equations.

formal decomposition of pressure:

of the «pressure components».

*∂p ∂x* = *∂ ∂x*

(segregated) or coupled manner.

**4. Development of explicit schemes**

= *∂px ∂x* + *∂pxyz ∂x* .

we obtain

Fig. 2. Driven cavity and location of the control volumes *V*1 and *V*2

$$\begin{aligned} \text{Stage I:} \quad & \frac{V^{(n+1/2)} - V^{(n)}}{h\_t} = -(V^{(n)} \nabla)V^{(n)} + \text{Re}^{-1} \Delta V^{(n)},\\ \text{Stage II:} \quad & \Delta p = \frac{\nabla V^{(n+1/2)}}{h\_t},\\ \text{Stage III:} \quad & \frac{V^{(n+1)} - V^{(n+1/2)}}{h\_t} = -\nabla p,\end{aligned}$$

where *ht* is time semispacing, *V*(*n*+1/2) is intermediate velocity field and *n* is a time layer. Stage I consists in solution of the momentum equations without pressure gradients. For simplicity X-momentum can be written as

$$\frac{u\_{ij}^{(n+1/2)} - u\_{ij}^{(n)}}{h\_t} = \psi\_{ij},\tag{11}$$

where *ψij* is the given function defined as

$$\psi\_{ij} = \left( -\frac{\partial(u^2)}{\partial x} - \frac{\partial(vu)}{\partial y} + \frac{1}{\text{Re}} \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) \right)\_{ij}^{(n)}.\tag{12}$$

It is easy to see that intermediate velocity field *V*(*n*+1/2) is independent on pressure. This disadvantage can be compensated partially by the pressure decomposition (10). Application of the decomposition requires two mass conservation equations for 2D problems. Integration of the continuity equation (1) over the control volumes *V*1 and *V*2 shown on Figure 2 gives

$$\int\_{0}^{1} u(t, \mathbf{x}, \mathbf{y}) \, d\mathbf{y} = 0, \qquad \int\_{0}^{1} v(t, \mathbf{x}, \mathbf{y}) \, d\mathbf{x} = 0. \tag{13}$$

Approximation of the mass conservation equations on the staggered grid is given by

$$h\_y \sum\_{j=1}^{N\_y} u\_{ij}^{(m)} = 0 \,,\tag{14}$$

$$h\_{\chi} \sum\_{i=1}^{N\_{\chi}} v\_{ij}^{(m)} = 0 \,, \tag{15}$$

where *m* = *n*, *n* + 1/2, *n* + 1 and N*x* = 1/*hx*, N*y* = 1/*hy*. As contrasted with equation (11) in the classical approach, the velocity component *u* and «one-dimensional component of pressure» *px* should satisfy to the system

$$\begin{cases} \frac{u\_{ij}^{(n+1/2)} - u\_{ij}^{(n)}}{h\_t} = -\left(\frac{\partial p^x}{\partial x}\right)\_i^{(n+1/2)} + \psi\_{ij} \\\ h\_y \sum\_{j=1}^{N\_y} u\_{ij}^{(n+1/2)} = 0 \end{cases} \tag{16}$$

i.e. *u* and *px* are computed in the coupled manner using the discrete mass conservation equation (14).

It is clear that the system (16) can be written in form of (5), where *A* is the diagonal matrix for explicit schemes. This fact allows obtain analytic solution of the saddle point system (16). Multiplication of the first equation in system (16) on *hy* and summation give

$$\frac{1}{h\_{\mathbf{f}}} \left( h\_y \sum\_{j=1}^{\mathcal{N}\_y} u\_{ij}^{(n+1/2)} - h\_y \sum\_{j=1}^{\mathcal{N}\_y} u\_{ij}^{(n)} \right) = -\sum\_{j=1}^{\mathcal{N}\_y} h\_y \left( \frac{\partial p^x}{\partial \mathbf{x}} \right)\_i^{(n+1/2)} + h\_y \sum\_{j=1}^{\mathcal{N}\_y} \psi\_{ij}. \tag{17}$$

Left-hand side of the equation equals zero due to equation (14). Furthermore

$$\sum\_{j=1}^{\mathbf{N}\_y} h\_y \left( \frac{\partial p^{\mathbf{x}}}{\partial \mathbf{x}} \right)\_i^{(n+1/2)} = \left( \frac{\partial p^{\mathbf{x}}}{\partial \mathbf{x}} \right)\_i^{(n+1/2)} \sum\_{j=1}^{\mathbf{N}\_y} h\_y = \left( \frac{\partial p^{\mathbf{x}}}{\partial \mathbf{x}} \right)\_i^{(n+1/2)}$$

because (*px*)� *i* is independent on *j* and N*y* ∑ *j*=1 *hy* = 1 is dimensionless height of the cavity. Equation (17) is reduced to

$$\left(\frac{\partial p^{\mathbf{x}}}{\partial \mathbf{x}}\right)\_{i}^{(n+1/2)} = h\_{\mathbf{y}} \sum\_{j=1}^{\mathbf{N}\_{\mathbf{y}}} \boldsymbol{\psi}\_{ij} \,. \tag{18}$$

,

Fig. 3. Ratio of execution time at flow simulation in the driven cavity

*u*(*n*+1/2) *ij* = *u*(*n*)

(*n*+1/2) *i* = (*px*)

Accounting decomposition (10), other stages of the algorithm are written as

Stage III: *V*(*n*+1) − *V*(*n*+1/2)

For the numerical experiment law of the lid motion is taken as

(*px*)

for Solving Navier–Stokes Equations on Structured Grids

0 is some constant. It is clear that solution of system (19) is

*ij* − *hthy*

183 Convergence Acceleration of Iterative Algorithms

(*n*+1/2) *i*−1 + *hxhy*

Velocity component *u* and «one-dimensional component of pressure» *px* are computed in the coupled manner saving explicit nature of the computation. Other velocity components are

Stage II: Δ*pxy* = ∇*V*(*n*+1/2)

In the stages only «multidimensional component» *pxy* is used for computation of the velocity

Reynolds number Re = 1000 is based on the cavity height and the lid velocity max *U*(*n*) *w* = 1. Staggered uniform grid *hx* = *hy* = *h* = 1/200, *ht* = *h*/5 is used for the flow simulation. Ratio

*U*(*n*) *w* = min *n*

*ht*

100 ; 1 .

N*y* ∑ *j*=1

*ψij* + *htψij* ,

N*y* ∑ *j*=1 *ψij* .

*ht*

= −∇*pxy* .

,

where (*px*)

field.

(*n*+1/2)

computed in the similar way.

Substitution of the equation into system (16) gives a new form of the system

$$\begin{cases} \frac{u\_{ij}^{(n+1/2)} - u\_{ij}^{(n)}}{h\_t} = -h\_y \sum\_{j=1}^{\mathcal{N}\_y} \Psi\_{ij} + \Psi\_{ij} \\\\ \frac{(p^x)\_i^{(n+1/2)} - (p^x)\_{i-1}^{(n+1/2)}}{h\_x} = h\_y \sum\_{j=1}^{\mathcal{N}\_y} \Psi\_{ij} \end{cases} \tag{19}$$

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

*u*(*m*)

where *m* = *n*, *n* + 1/2, *n* + 1 and N*x* = 1/*hx*, N*y* = 1/*hy*. As contrasted with equation (11) in the classical approach, the velocity component *u* and «one-dimensional component of

= −

i.e. *u* and *px* are computed in the coupled manner using the discrete mass conservation

It is clear that the system (16) can be written in form of (5), where *A* is the diagonal matrix for explicit schemes. This fact allows obtain analytic solution of the saddle point system (16).

> N*y* ∑ *j*=1 *hy* � *∂px ∂x*

�(*n*+1/2)

= *hy*

*ij*

(*n*+1/2) *i*−1

N*y* ∑ *j*=1

N*y* ∑ *j*=1

= −*hy*

= *hy*

N*y* ∑ *j*=1

N*y* ∑ *j*=1 *ψij*

*ψij* + *ψij*

*hy* =

*i*

N*y* ∑ *j*=1

�(*n*+1/2)

*i*

Substitution of the equation into system (16) gives a new form of the system

*hx*

*u*(*n*+1/2) *ij* − *u*(*n*)

*ht*

�*∂px ∂x*

�(*n*+1/2)

+ *ψij*

�(*n*+1/2)

+ *hy*

�(*n*+1/2)

*ψij* . (18)

, (19)

*i*

*hy* = 1 is dimensionless height of the cavity.

N*y* ∑ *j*=1

,

*ψij* . (17)

*i*

�*∂px ∂x*

*i*

*ij* = 0 , (14)

*ij* = 0 , (15)

, (16)

Approximation of the mass conservation equations on the staggered grid is given by

*hy* N*y* ∑ *j*=1

*hx* N*x* ∑ *i*=1 *v* (*m*)

*ij*

*u*(*n*+1/2) *ij* = 0

Multiplication of the first equation in system (16) on *hy* and summation give

Left-hand side of the equation equals zero due to equation (14). Furthermore

= � *∂px ∂x*

� *∂px ∂x*

*u*(*n*) *ij*

⎞ ⎠ = −

N*y* ∑ *j*=1

pressure» *px* should satisfy to the system

equation (14).

because (*px*)�

1 *ht* ⎛ ⎝*hy*

> N*y* ∑ *j*=1 *hy* � *∂px ∂x*

Equation (17) is reduced to

N*y* ∑ *j*=1 ⎧ ⎪⎪⎪⎪⎪⎨

*u*(*n*+1/2) *ij* − *u*(*n*)

*ht*

⎪⎪⎪⎪⎪⎩

*u*(*n*+1/2) *ij* − *hy*

�(*n*+1/2)

*i*

*i* is independent on *j* and

⎧ ⎪⎪⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎪⎪⎩

(*px*)

(*n*+1/2) *i* − (*px*)

*hy* N*y* ∑ *j*=1

Fig. 3. Ratio of execution time at flow simulation in the driven cavity where (*px*) (*n*+1/2) 0 is some constant. It is clear that solution of system (19) is

$$\begin{aligned} u\_{ij}^{(n+1/2)} &= u\_{ij}^{(n)} - h\_t h\_y \sum\_{j=1}^{\mathcal{N}\_y} \psi\_{ij} + h\_t \psi\_{ij}, \\\\ (p^x)\_i^{(n+1/2)} &= (p^x)\_{i-1}^{(n+1/2)} + h\_x h\_y \sum\_{j=1}^{\mathcal{N}\_y} \psi\_{ij}. \end{aligned}$$

Velocity component *u* and «one-dimensional component of pressure» *px* are computed in the coupled manner saving explicit nature of the computation. Other velocity components are computed in the similar way.

Accounting decomposition (10), other stages of the algorithm are written as

$$\begin{aligned} \text{Stage II:} & \qquad \Delta p^{xy} = \frac{\nabla V^{(n+1/2)} }{h\_t} \\ \text{Stage III:} & \qquad \frac{V^{(n+1)} - V^{(n+1/2)}}{h\_t} = -\nabla p^{xy} .\end{aligned}$$

,

In the stages only «multidimensional component» *pxy* is used for computation of the velocity field.

For the numerical experiment law of the lid motion is taken as

$$\mathcal{U}\_w^{(n)} = \min\left(\frac{n}{100}; 1\right).$$

Reynolds number Re = 1000 is based on the cavity height and the lid velocity max *U*(*n*) *w* = 1. Staggered uniform grid *hx* = *hy* = *h* = 1/200, *ht* = *h*/5 is used for the flow simulation. Ratio

work for approximated solution of the auxiliary problem is negligible small as compared with the total efforts. Note that the equations of the auxiliary problem are not pressure-linked. To illustrate influence of the auxiliary problem on convergence rate, we use Uzawa algorithm (6) for simulation of stationary flow in the driven cavity starting the iterand zero: *u*(0) = 0, *v*(0) = 0 and *p*(0) = 0. Accounting zero boundary conditions for *v*, first equation of system (6)

185 Convergence Acceleration of Iterative Algorithms

Re � *∂*2*u ∂x*2 +

> 1 Re � *∂*2*u ∂x*2 +

and *v* = 0. Finally both problem (22) and (23) are reduced to systems of linear algebraic

�*b*� < 10−7 .

The computations are performed with Re = 100 on uniform staggered grid 101 × 101 (*hx* =

Figure 4 represents solution of the Navier–Stokes equations in "stream function–vorticity" (+) Ghia et al. (1982), primitive variables formulations (—) and solutions of equations (22) and (23) in the middle section of the cavity (*x* = 0.5) at Re = 100. It is easy to see that use of the mass conservation equations in the auxiliary problem makes it possible to obtain more

�*Ax* − *b*�

accurate approximation to solution of the full Navier–Stokes equations (1)–(3).

Fig. 4. Distribution of the velocity component *u* in the middle section of the cavity

*dx* +

*∂*2*u ∂y*2 �

> *∂*2*u ∂y*2 �

(22)

(23)

*∂*(*u*2) *∂x* = 1

*∂x* = − *dpx*

*u*(*x*, *y*) *dy* = 0

and *v* = 0. In the auxiliary problem the system (20) takes the form

*∂*(*u*2)

equations *Ax* = *b*. For clearness these equations are solved until

⎧ ⎪⎪⎪⎨

for Solving Navier–Stokes Equations on Structured Grids

⎪⎪⎪⎩

� 1 0

is reduced to

*hy* = 1/100).

of the execution time *T*(*n*) *m* /*T*(*n*) *c* is used as a criterion of the convergence acceleration, where *T*(*n*) *m* and *T*(*n*) *c* are execution time for abovementioned and classical approaches, respectively. Figure 3 shows result of the numerical test. Obtained result for *n* = 200

$$\frac{1}{200} \sum\_{n=1}^{200} T\_m^{(n)} / T\_c^{(n)} = 0.81$$

illustrate the least acceleration efficiency arising at simulation of the rotated flows.

#### **5. Development of implicit schemes**

Application of the pressure decomposition (10) for improvement of the implicit schemes requires solution of an auxiliary problem because of the «pressure» gradients can not be determined in explicit form such as equation (18).

#### **5.1 Auxiliary problem**

Auxiliary problem is intended for fast computation of the «one-dimensional components» *px*(*t*, *x*), *py*(*t*, *y*) and *pz*(*t*, *z*) in decomposition (10). It is assumed that the solution of the auxiliary problem will be close to the solution of the Navier–Stokes equations.

Formulation of the auxiliary problem is based on replacement of the continuity equation (1) by the mass conservation equations. For example, for the driven cavity (Figure 2) the auxiliary problem with the mass conservation equations (13) instead of the continuity equation (1) takes the form:

a) X-momentum and mass conservation equations

$$\begin{cases} \frac{\partial u}{\partial t} + \frac{\partial (u^2)}{\partial x} + \frac{\partial (vu)}{\partial y} = -\frac{\partial p^x}{\partial x} - \left[ \frac{\partial p^{xy}}{\partial x} \right] + \frac{1}{\text{Re}} \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) \\\quad \int u(t, x, y) \, dy = 0 \end{cases} \tag{20}$$

b) Y-momentum and mass conservation equations

$$\begin{cases} \frac{\partial v}{\partial t} + \frac{\partial (uv)}{\partial x} + \frac{\partial (v^2)}{\partial y} = -\frac{\partial p^y}{\partial y} - \left[\frac{\partial p^{xy}}{\partial y}\right] + \frac{1}{\text{Re}} \left(\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}\right) \\ \quad \int v(t, x, y) \, d\mathbf{x} = 0 \end{cases} \tag{21}$$

where square brackets mean that the «pressure» gradients (*pxy*)� *x* and (*pxy*)� *y* are fixed (i.e. its values are taken from previous iteration). Braces mean that the momentum and mass conservation equations are solved only in coupled manner.

Since the systems (20) and (21) are similar to the simplified Navier–Stokes equations (7), the systems can be solved by the same numerical methods. Main difference consists in stopping criterion: auxiliary problem can be solved approximately, i.e. it is necessary to perform several iterations of line Seidel method with the secant iterations. As a result, extra computational

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

of the execution time *T*(*n*) *m* /*T*(*n*) *c* is used as a criterion of the convergence acceleration, where *T*(*n*) *m* and *T*(*n*) *c* are execution time for abovementioned and classical approaches, respectively.

Application of the pressure decomposition (10) for improvement of the implicit schemes requires solution of an auxiliary problem because of the «pressure» gradients can not be

Auxiliary problem is intended for fast computation of the «one-dimensional components» *px*(*t*, *x*), *py*(*t*, *y*) and *pz*(*t*, *z*) in decomposition (10). It is assumed that the solution of the

Formulation of the auxiliary problem is based on replacement of the continuity equation (1) by the mass conservation equations. For example, for the driven cavity (Figure 2) the auxiliary problem with the mass conservation equations (13) instead of the continuity equation (1) takes

*∂x* −

*∂y* − � *∂pxy ∂y*

its values are taken from previous iteration). Braces mean that the momentum and mass

Since the systems (20) and (21) are similar to the simplified Navier–Stokes equations (7), the systems can be solved by the same numerical methods. Main difference consists in stopping criterion: auxiliary problem can be solved approximately, i.e. it is necessary to perform several iterations of line Seidel method with the secant iterations. As a result, extra computational

� *∂pxy ∂x*

� + 1 Re

� + 1 Re �*∂*2*u ∂x*2 +

� *∂*2*v ∂x*2 + *∂*2*u ∂y*2 �

*∂*2*v ∂y*2 �

*x* and (*pxy*)�

, (20)

, (21)

*y* are fixed (i.e.

*T*(*n*) *m* /*T*(*n*) *c* = 0.81

Figure 3 shows result of the numerical test. Obtained result for *n* = 200

1 200

**5. Development of implicit schemes**

**5.1 Auxiliary problem**

⎧ ⎪⎪⎪⎨

*∂u ∂t* +

� 1

0

*∂v ∂t* +

� 1 0

⎪⎪⎪⎩

⎧ ⎪⎪⎪⎨

⎪⎪⎪⎩

the form:

determined in explicit form such as equation (18).

a) X-momentum and mass conservation equations

*∂*(*u*2) *∂x* + *∂*(*vu*)

*u*(*t*, *x*, *y*) *dy* = 0

b) Y-momentum and mass conservation equations

*∂*(*uv*) *∂x*

*v*(*t*, *x*, *y*) *dx* = 0

+ *∂*(*v*2)

conservation equations are solved only in coupled manner.

where square brackets mean that the «pressure» gradients (*pxy*)�

200 ∑ *n*=1

illustrate the least acceleration efficiency arising at simulation of the rotated flows.

auxiliary problem will be close to the solution of the Navier–Stokes equations.

*∂y* = − *∂px*

*∂y* = − *∂py*

work for approximated solution of the auxiliary problem is negligible small as compared with the total efforts. Note that the equations of the auxiliary problem are not pressure-linked.

To illustrate influence of the auxiliary problem on convergence rate, we use Uzawa algorithm (6) for simulation of stationary flow in the driven cavity starting the iterand zero: *u*(0) = 0, *v*(0) = 0 and *p*(0) = 0. Accounting zero boundary conditions for *v*, first equation of system (6) is reduced to

$$\frac{\partial(u^2)}{\partial x} = \frac{1}{\text{Re}} \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) \tag{22}$$

and *v* = 0. In the auxiliary problem the system (20) takes the form

$$\begin{cases} \frac{\partial (u^2)}{\partial x} = -\frac{dp^x}{dx} + \frac{1}{\text{Re}} \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) \\\\ \int\_0^1 u(x, y) \, dy = 0 \end{cases} \tag{23}$$

and *v* = 0. Finally both problem (22) and (23) are reduced to systems of linear algebraic equations *Ax* = *b*. For clearness these equations are solved until

$$\frac{||Ax - b||}{||b||} < 10^{-7}\,\text{.}$$

The computations are performed with Re = 100 on uniform staggered grid 101 × 101 (*hx* = *hy* = 1/100).

Figure 4 represents solution of the Navier–Stokes equations in "stream function–vorticity" (+) Ghia et al. (1982), primitive variables formulations (—) and solutions of equations (22) and (23) in the middle section of the cavity (*x* = 0.5) at Re = 100. It is easy to see that use of the mass conservation equations in the auxiliary problem makes it possible to obtain more accurate approximation to solution of the full Navier–Stokes equations (1)–(3).

Fig. 4. Distribution of the velocity component *u* in the middle section of the cavity

(a) Isolines of stream function (b) Isobars

*u*(*t*, *x*, *y*) *dy* =

*y*

Numerical experiments show that execution time can be reduced in ∼ 400 times for the given problem (staggered grid 101 × 1401, unpreconditioned Uzawa algorithm, Re = 800). Figure 6 explains the impressive reduction of the computational efforts. It is easy to see that pressure is changed mainly in *x* direction except small subdomain near attachment point of bottom eddy (i.e. *p*(*x*, *y*) ≈ *px*(*x*)). Since the «one-dimensional component of the pressure» *px*(*x*) is computed in the auxiliary problem, the proposed algorithm is very efficient for solving the

Authors *lB lT wT xTL xTR* Nodes

present 6.1000 5.6300 0.28 4.8400 10.4700 141501

Table 1. Comparison of results of the flow simulation over backward-facing step (Re = 800)

Gartling (1990) 6.1000 5.6300 – 4.8500 10.4800 129681 Gresho et al. (1993) 6.0820 5.6260 – 4.8388 10.4648 245760 Gresho et al. (1993) 6.1000 5.6300 – 4.8600 10.4900 8000 Keskar & Lin (1999) 6.0964 5.6251 – 4.8534 10.4785 3737

Barton (1997) 6.0150 5.6600 – 4.8200 10.4800

0

Redefining velocity components to be zero inside the step, we obtain the following mass

187 Convergence Acceleration of Iterative Algorithms

 H

*u*(*t*, 0, *y*) *dy* ,

*u*(*t*, *L*, *ξ*) − *u*(*t*, 0, *ξ*)

 *dξ* .

0

Fig. 6. Stationary flow over a backward-facing step

for Solving Navier–Stokes Equations on Structured Grids

 H

0

*v*(*t*, *x*, *y*) *dx* = −

conservation equations for the given problem

 L

0

Table 1 represents comparison of obtained results.

problem.

#### **5.2 Main problem**

Accounting the pressure decomposition (10), the momentum equations in the main problem are written as

$$
\frac{\partial u}{\partial t} + \frac{\partial (u^2)}{\partial x} + \frac{\partial (vu)}{\partial y} = -\left[\frac{dp^x}{d\mathbf{x}}\right] - \frac{\partial p^{xy}}{\partial \mathbf{x}} + \frac{1}{\mathbf{Re}} \left(\frac{\partial^2 u}{\partial \mathbf{x}^2} + \frac{\partial^2 u}{\partial y^2}\right), \tag{24}
$$

$$\frac{\partial v}{\partial t} + \frac{\partial (uv)}{\partial x} + \frac{\partial (v^2)}{\partial y} = -\left[\frac{dp^y}{dy}\right] - \frac{\partial p^{xy}}{\partial y} + \frac{1}{\text{Re}} \left(\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}\right),\tag{25}$$

where square brackets mean that the «pressure» gradients (*px*)� and (*py*)� are fixed (i.e. the gradients have been computed in the auxiliary problem (such as equation (20) and (21) for the driven cavity)). Main problem consists of momentum (24), (25) and continuity equations (1). Algorithm for simulation of the flows with given mass flow rate can be represented as:

Stage I: *auxiliary problem*: several iterations of line (2D) or plane (3D) Seidel method with the secant iterations

Stage II: *main problem*: iterations of basic method (SIMPLE, Uzawa or Vanka iterations, etc.)

Stage III: check convergence, continue (go to 1) if necessary

#### **5.3 Flow over a backward-facing step**

The next benchmark problem about backward-facing step flow is used for illustration of the impressive convergence acceleration for the directed fluid flows.

Consider the stationary laminar flow over a backward-facing step, which is another well studied test case. Figure 5 shows the geometry of the flow. The fact that the solution of the incompressible Navier–Stokes equations over a backward-facing step at Re = 800 is steady and stable has been confirmed in a number of recent works.

No-slip boundary conditions are imposed on the step and the upper and lower walls, a parabolic velocity *u* profile is specified at the channel inlet (*v* = 0), and zero natural boundary conditions (*v* = 0 and *u*� *x* = 0) are imposed at the channel outlet. The Reynolds number Re is based on the channel height (H = 1) and the average inlet velocity in the parabolic profile. The channel length is L = 14.

Fig. 5. Geometry of problem about the backward-facing step flow

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

Accounting the pressure decomposition (10), the momentum equations in the main problem

where square brackets mean that the «pressure» gradients (*px*)� and (*py*)� are fixed (i.e. the gradients have been computed in the auxiliary problem (such as equation (20) and (21) for the driven cavity)). Main problem consists of momentum (24), (25) and continuity equations (1). Algorithm for simulation of the flows with given mass flow rate can be represented as:

Stage I: *auxiliary problem*: several iterations of line (2D) or plane (3D) Seidel method with the

Stage II: *main problem*: iterations of basic method (SIMPLE, Uzawa or Vanka iterations, etc.)

The next benchmark problem about backward-facing step flow is used for illustration of the

Consider the stationary laminar flow over a backward-facing step, which is another well studied test case. Figure 5 shows the geometry of the flow. The fact that the solution of the incompressible Navier–Stokes equations over a backward-facing step at Re = 800 is steady

No-slip boundary conditions are imposed on the step and the upper and lower walls, a parabolic velocity *u* profile is specified at the channel inlet (*v* = 0), and zero natural boundary

is based on the channel height (H = 1) and the average inlet velocity in the parabolic profile.

*x* = 0) are imposed at the channel outlet. The Reynolds number Re

 *∂*2*u ∂x*2 +

 *∂*2*v ∂x*2 + *∂*2*u ∂y*2

*∂*2*v ∂y*2  , (24)

, (25)

 *dpx dx* − *∂pxy ∂x* + 1 Re

 *dpy dy* − *∂pxy ∂y* + 1 Re

**5.2 Main problem**

*∂u ∂t* +

*∂v ∂t* +

**5.3 Flow over a backward-facing step**

secant iterations

conditions (*v* = 0 and *u*�

The channel length is L = 14.

*∂*(*u*2) *∂x* + *∂*(*vu*) *∂y* = −

*∂*(*uv*) *∂x*

+ *∂*(*v*2) *∂y* = −

Stage III: check convergence, continue (go to 1) if necessary

impressive convergence acceleration for the directed fluid flows.

and stable has been confirmed in a number of recent works.

Fig. 5. Geometry of problem about the backward-facing step flow

are written as

Fig. 6. Stationary flow over a backward-facing step

Redefining velocity components to be zero inside the step, we obtain the following mass conservation equations for the given problem

$$\int\_{0}^{\mathcal{H}} u(t, x, y) \, dy = \int\_{0}^{\mathcal{H}} u(t, 0, y) \, dy$$

$$\int\_{0}^{\mathcal{L}} v(t, x, y) \, dx = -\int\_{0}^{y} \left( u(t, L, \vec{\xi}) - u(t, 0, \vec{\xi}) \right) d\xi \,.$$

Numerical experiments show that execution time can be reduced in ∼ 400 times for the given problem (staggered grid 101 × 1401, unpreconditioned Uzawa algorithm, Re = 800). Figure 6 explains the impressive reduction of the computational efforts. It is easy to see that pressure is changed mainly in *x* direction except small subdomain near attachment point of bottom eddy (i.e. *p*(*x*, *y*) ≈ *px*(*x*)). Since the «one-dimensional component of the pressure» *px*(*x*) is computed in the auxiliary problem, the proposed algorithm is very efficient for solving the problem.

Table 1 represents comparison of obtained results.

Table 1. Comparison of results of the flow simulation over backward-facing step (Re = 800)

flow in the catalyst. Diffusion-dominant nature of fluid flow in the microcatalyst simplifies the grid generation. Example of the simplest computational grid for this problem is shown on

189 Convergence Acceleration of Iterative Algorithms

Nonuniform staggered grid 385 × 3150 is used for the flow simulation (Re = 350). Figure 9 represents distribution of the stream function and pressure near first column of the needles. Chess order of the needle location results in eddy-free flow inside the microcatalyst. However

> 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35

Recently the numerical methods for fluid flow prediction have been classified into two categories: density-based and pressure-based. For the pressure-based approach, methods are

4.40 4.50 4.60

0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85

(b) Isobars

Figure 8. No-slip conditions are approximated exactly on the needle surfaces.

intensive eddy formation after last column of the needles is observed (Figure 10).

0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85

for Solving Navier–Stokes Equations on Structured Grids

(a) Isolines of stream function

Fig. 9. Flow near first column of the needles

**5.5 Compressible ow in laval micronozzle**

0.00

Fig. 10. Eddy formation after last column of the needles

0.10

0.20

0.30

0.00

0.10

0.20

0.30

Fig. 8. Staggered grid in the microcatalyst

#### **5.4 Flow in microcatalyst**

Proposed approach has been used for simulation of incompressible fluid flows in microcatalyst. The microcatalyst represents 2D channel with iridium-covered needles located in chess order as shown on Figure 7.

Redefining velocity components to be zero inside the needles, there is no remarkable difference in formulation of the auxiliary problem for flow over backward-facing step and for 14 Will-be-set-by-IN-TECH

Proposed approach has been used for simulation of incompressible fluid flows in microcatalyst. The microcatalyst represents 2D channel with iridium-covered needles located

Redefining velocity components to be zero inside the needles, there is no remarkable difference in formulation of the auxiliary problem for flow over backward-facing step and for

Fig. 7. Geometry of the microcatalyst

Fig. 8. Staggered grid in the microcatalyst

in chess order as shown on Figure 7.

**5.4 Flow in microcatalyst**

flow in the catalyst. Diffusion-dominant nature of fluid flow in the microcatalyst simplifies the grid generation. Example of the simplest computational grid for this problem is shown on Figure 8. No-slip conditions are approximated exactly on the needle surfaces.

Nonuniform staggered grid 385 × 3150 is used for the flow simulation (Re = 350). Figure 9 represents distribution of the stream function and pressure near first column of the needles. Chess order of the needle location results in eddy-free flow inside the microcatalyst. However intensive eddy formation after last column of the needles is observed (Figure 10).

#### **5.5 Compressible ow in laval micronozzle**

Recently the numerical methods for fluid flow prediction have been classified into two categories: density-based and pressure-based. For the pressure-based approach, methods are

Fig. 10. Eddy formation after last column of the needles

where

follows

*∂ ∂t* � 1

*py*, *pz* and *pxyz* are fixed.

0

�*x*

0

*U* =

*ρu*2 + *p* − 4

*ρuv* − *y*¯*y* Re *∂v ∂y*¯ − *�*

*ρvu* − *y*¯*y* Re *∂u ∂y*¯ − *�*

> 3 *y*¯*y* Re *∂v ∂y*¯ + 2 3 *�* Re � *∂u*

*ρv*2 + *p* − 4

*J dξ dy* + *�*

*E* =

*F* =

⎛

for Solving Navier–Stokes Equations on Structured Grids

⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝

⎛

⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝

Parameter *�* is the micronozzle width-to-length ratio.

*ρ*(*t*, *ξ*, *y*)

X-momentum should be updated using pressure

⎛

*ρ ρu ρv ρi*

⎞

⎟⎟⎠ , *H* =

191 Convergence Acceleration of Iterative Algorithms

*ρu*

*∂x*¯ + *y*¯*x*

Pe � *∂T*

*ρv*

*ρvi* − *y*¯*y* Pe *∂T ∂y*¯

First mass conservation equation is obtained by integration of the continuity equation as

�*ρu J* �� � � � *x*

In the auxiliary problem for incompressible flows, iterations of line (2D) or plane (3D) Seidel method are stopped then the velocity component satisfies to the mass conservation equation. Computation of compressible flows requires updating of thermophysical properties (density *ρ*, coefficient of viscosity in Re and heat conductivity coefficient in Pe) using updated pressure in the line or plane. In 3D case values of thermophysical properties of the fluid for

*p*(*t*, *x*, *y*, *z*) = *px*(*t*, *x*)+[*py*(*t*, *y*) + *pz*(*t*, *z*) + *pxyz*(*t*, *x*, *y*, *z*)],

temperature *T* and equation of state. Here square brackets mean that the pressure components

Figure 12 represents isobars in the Laval micronozzles. It is easy to see that the isobars are almost vertical lines near throat and in supersonic part of the micronozzle. It means that the pressure is changed mainly along the micronozzle axis. In other words, «one-dimensional component of the pressure» *px* in decomposition (10) is dominant in this problem. For

� 1

0

Re � *∂v*

*∂x*¯ + *y*¯*x*

Re � *∂v*

⎛

⎞

⎟⎟⎠ ,

*∂v ∂y*¯ �

*∂v ∂y*¯ �

> *∂u ∂y*¯ �

*∂x*¯ + *y*¯*x*

⎞

⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ ,

⎞

⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ .

*dy* = 0 .

⎜⎜⎝

*∂u ∂y*¯ � + 2 3 *y*¯*y* Re *∂v ∂y*¯

*∂x*¯ + *y*¯*x*

*∂T ∂y*¯ �

*∂x*¯ + *y*¯*x*

*dy* − *�* � 1

0

� *ρu J* �� � � � 0

⎜⎜⎝

3 *�* Re � *∂u*

*ρui* − *�*

classified into coupled and segregated (decoupled). Density-based algorithms traditionally are used to compute high speed compressible flows. Computational methods for low Mach number compressible flows are an active research field in recent years. The pressure-velocity coupling problem discussed earlier for incompressible flows are also encountered in the methods when used for low-speed applications.

Pressure decomposition (10) shows that there are not pure density-based and segregated solvers because of the velocity components and corresponding «one-dimensional components of pressure» (i.e. (*u*,*px*), (*v*,*py*) and (*w*,*pz*)) always are computed in the coupled manner. «Multidimensional component» *pxyz* in (10) can be computed by coupled or segregated method using density-based or pressure-based approach.

Consider application of the pressure decomposition for simulation of compressible flow in flat Laval micronozzle. Width of subsonic part of the micronozzle is 1 mm. Grid generation is based on mapping of the non-dimensional physical domain with nonuniform grid onto computational domain (unit square) with uniform grid. Direct (ABCD → A¯ B¯C¯ D¯ ) and reverse (ABCD ← A¯ B¯ C¯ D¯ ) mappings are shown on Figure 11, where the function *ϕ*(*x*) describes the micronozzle profile. The mappings can be given by

$$\mathfrak{x} = \mathfrak{x}, \qquad \mathfrak{y} = -\frac{1}{\beta} \ln \left( 1 - (1 - e^{-\beta}) \frac{\mathcal{Y}}{q(\mathfrak{x})} \right) \mathfrak{y}$$

where (*x*, *y*) and (*x*¯, *y*¯) are spatial variables in physical and computational domains, respectively. Parameter *β* > 0 is intended for the grid refinement near solid wall. Jacobian (*J*) of the mapping

$$J = \begin{vmatrix} \bar{\boldsymbol{x}}\_{\mathcal{X}} & \bar{\boldsymbol{x}}\_{\mathcal{Y}} \\ \bar{\boldsymbol{y}}\_{\mathcal{X}} \ \bar{\boldsymbol{y}}\_{\mathcal{Y}} \end{vmatrix} = \frac{1 - e^{-\beta}}{\beta} \frac{e^{\beta \mathcal{G}}}{\boldsymbol{\varrho}(\bar{\boldsymbol{x}})}$$

is non-singular (*J* �= 0). In addition, *J* → 1/*ϕ*(*x*¯) at *β* → 0 for uniform grid in *y* direction.

Fig. 11. Non-dimensional physical and computational domains

Finally, non-dimensional compressible Navier–Stokes equations in the computational domain are written as

$$\frac{\partial}{\partial t}\left(\frac{\mathbf{U}}{J}\right) + \epsilon \frac{\partial}{\partial \vec{x}}\left(\frac{E}{J}\right) + \epsilon \frac{\partial}{\partial \vec{y}}\left(\frac{\vec{y}\_{\vec{x}}E}{J}\right) + \frac{\partial F}{\partial \vec{y}} = \frac{H}{J}\lambda$$

where

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

classified into coupled and segregated (decoupled). Density-based algorithms traditionally are used to compute high speed compressible flows. Computational methods for low Mach number compressible flows are an active research field in recent years. The pressure-velocity coupling problem discussed earlier for incompressible flows are also encountered in the

Pressure decomposition (10) shows that there are not pure density-based and segregated solvers because of the velocity components and corresponding «one-dimensional components of pressure» (i.e. (*u*,*px*), (*v*,*py*) and (*w*,*pz*)) always are computed in the coupled manner. «Multidimensional component» *pxyz* in (10) can be computed by coupled or segregated

Consider application of the pressure decomposition for simulation of compressible flow in flat Laval micronozzle. Width of subsonic part of the micronozzle is 1 mm. Grid generation is based on mapping of the non-dimensional physical domain with nonuniform grid onto computational domain (unit square) with uniform grid. Direct (ABCD → A¯ B¯C¯ D¯ ) and reverse (ABCD ← A¯ B¯ C¯ D¯ ) mappings are shown on Figure 11, where the function *ϕ*(*x*) describes the

> *β* ln

> >

is non-singular (*J* �= 0). In addition, *J* → 1/*ϕ*(*x*¯) at *β* → 0 for uniform grid in *y* direction.

Finally, non-dimensional compressible Navier–Stokes equations in the computational domain

 *y*¯*xE J*

 + *∂F ∂y*¯ = *H J* ,

respectively. Parameter *β* > 0 is intended for the grid refinement near solid wall.

where (*x*, *y*) and (*x*¯, *y*¯) are spatial variables in physical and computational domains,

 = 1 − *e*−*β β*

1 − (1 − *e*

−*β*) *y ϕ*(*x*) ,

*eβy*¯ *ϕ*(*x*¯)

methods when used for low-speed applications.

method using density-based or pressure-based approach.

micronozzle profile. The mappings can be given by

Jacobian (*J*) of the mapping

are written as

*x*¯ = *x* , *y*¯ = − 1

*J* = *x*¯*x x*¯*y y*¯*x y*¯*y*

Fig. 11. Non-dimensional physical and computational domains

*∂ ∂t* *U J* + *� ∂ ∂x*¯ *E J* + *� ∂ ∂y*¯

$$\mathcal{U} = \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ \rho i \end{pmatrix}, \qquad H = \begin{pmatrix} 0 \\ 0 \\ 0 \\ S \end{pmatrix},$$

$$\mathbf{E} = \begin{pmatrix} \rho u \\ \rho u^2 + p - \frac{4}{3} \frac{\varepsilon}{\text{Re}} \left( \frac{\partial u}{\partial \bar{x}} + \bar{y}\_x \frac{\partial u}{\partial \bar{y}} \right) + \frac{2}{3} \frac{\bar{y}\_y}{\text{Re}} \frac{\partial v}{\partial \bar{y}} \\\ \rho uv - \frac{\bar{y}\_y}{\text{Re}} \frac{\partial v}{\partial \bar{y}} - \frac{\varepsilon}{\text{Re}} \left( \frac{\partial v}{\partial \bar{x}} + \bar{y}\_x \frac{\partial v}{\partial \bar{y}} \right) \\\ \rho ui - \frac{\varepsilon}{\text{Re}} \left( \frac{\partial T}{\partial \bar{x}} + \bar{y}\_x \frac{\partial T}{\partial \bar{y}} \right) \end{pmatrix}.$$

$$F = \begin{pmatrix} \rho v \\ \rho vu - \frac{\ddot{y}\_y}{\text{Re}} \frac{\partial u}{\partial \ddot{y}} - \frac{\epsilon}{\text{Re}} \left( \frac{\partial v}{\partial \ddot{x}} + \ddot{y}\_x \frac{\partial v}{\partial \ddot{y}} \right) \\\\ \rho v^2 + p - \frac{4}{3} \frac{\ddot{y}\_y}{\text{Re}} \frac{\partial v}{\partial \ddot{y}} + \frac{2}{3} \frac{\epsilon}{\text{Re}} \left( \frac{\partial u}{\partial \ddot{x}} + \ddot{y}\_x \frac{\partial u}{\partial \ddot{y}} \right) \\\\ \rho vi - \frac{\ddot{y}\_y}{\text{Pe}} \frac{\partial T}{\partial \ddot{y}} \end{pmatrix}.$$

Parameter *�* is the micronozzle width-to-length ratio.

First mass conservation equation is obtained by integration of the continuity equation as follows

$$\frac{\partial}{\partial t} \int\_0^1 \int\_0^x \frac{\rho(t, \xi, y)}{J} \, d\xi \, dy + \epsilon \int\_0^1 \left( \frac{\rho u}{J} \right) \Big|\_{x} \, dy - \epsilon \int\_0^1 \left( \frac{\rho u}{J} \right) \Big|\_{0} \, dy = 0 \dots$$

In the auxiliary problem for incompressible flows, iterations of line (2D) or plane (3D) Seidel method are stopped then the velocity component satisfies to the mass conservation equation. Computation of compressible flows requires updating of thermophysical properties (density *ρ*, coefficient of viscosity in Re and heat conductivity coefficient in Pe) using updated pressure in the line or plane. In 3D case values of thermophysical properties of the fluid for X-momentum should be updated using pressure

$$p(t, \mathbf{x}, y, z) = p^{\mathbf{x}}(t, \mathbf{x}) + [p^y(t, y) + p^z(t, z) + p^{xyz}(t, \mathbf{x}, y, z)]\_{\mathbf{x}}$$

temperature *T* and equation of state. Here square brackets mean that the pressure components *py*, *pz* and *pxyz* are fixed.

Figure 12 represents isobars in the Laval micronozzles. It is easy to see that the isobars are almost vertical lines near throat and in supersonic part of the micronozzle. It means that the pressure is changed mainly along the micronozzle axis. In other words, «one-dimensional component of the pressure» *px* in decomposition (10) is dominant in this problem. For

**X**

for Solving Navier–Stokes Equations on Structured Grids

**X**

(b) Isobars

5.2 5.4 5.6 5.8 6 6.2

**Y**

0.2

0.4

0.6

0.8

1

193 Convergence Acceleration of Iterative Algorithms

(a) Isolines of stream function

 1

*u*(*t*, *xi*, *y*) *dy* = max

*u*(*t*, *xi*, *y*) *dy* = min

*i*

*i* 1

at decelerated plunger motion. As a result, the auxiliary problem ensures expansion of perturbations caused by the plunger on all domain in each iteration. Uniform staggered grid 201× 2401 (*hx* = *hy* = 1/200) is used for simulation the fluid flow around the moving plunger. Reynolds number Re = 200 is based on the plunger velocity and distance between the plates. Stream function isolines and isobars near the moving plunger are shown on Figure 14.

Proposed convergence acceleration technique based on the pressure decomposition (10) should be incorporated with well-known algorithms for solving Navier–Stokes equations. Multigrid methods having (almost) optimal convergence rate for many applications seem to be the most promising solvers for many CFD problems. Our purpose is development of multigrid method with the least number of problem-dependent components for using in black

To overcome problem of robustness, the Navier–Stokes equations (4) should be adapted for the multigrid algorithm Martynenko (2006). Adaptation of the Navier–Stokes equations (so called Σ-modification) consists in representation of the velocity *V* and pressure *P* as sum of

 1

*u*˜(*t*, *xi*, *y*) *dy*

*u*˜(*t*, *xi*, *y*) *dy*

0

0

0

 1

0

**6. Pressure decomposition in geometric multigrid methods**

at accelerated or uniform plunger motion and in form

5.2 5.4 5.6 5.8 6 6.2

Fig. 14. Flow near the moving plunger

**Y**

0

formulated as

box software.

**6.1 Nonlinear multigrid iterations**

0.2

0.4

0.6

0.8

1

Fig. 13. Starting location of the plunger

the given problem, the auxiliary problem makes it possible to compute the most «part of pressure» (i.e. *px*(*t*, *x*) + *py*(*t*, *y*)) and corresponding change of thermophysical properties of the fluid based on simplified (pressure-unlinked) momentum equations in primitive variables formulation and mass conservation equations.

Fig. 12. Isobars in Laval micronozzle (sizes in meters)

#### **5.6 Flows with unspecied mass rate**

Previously fluid flows with the given mass rate have been simulated and analyzed. In some applications mass flow rate cannot be given in advance. The simplest example of such flows is problem about moving plunger. Immovable plunger is located between parallel plates filled by incompressible fluid. Figure 13 shows geometry of the problem and starting location of the plunger. Motion of the plunger causes the fluid flow. It is clear that mass flow rate depends on the plunger speed. Algorithm for simulation of the flows with unspecified mass rate should be modified as:

Stage I: *main problem*: iterations of basic method (SIMPLE, Uzawa or Vanka iterations, etc.)

Stage II: check convergence; continue if necessary

Stage III: *auxiliary problem*: several iterations of line (2D) or plane (3D) Seidel method with the secant iterations; continue (go to 1)

Assume that *u*˜*ij*, *v*˜*ij* and *p*˜ *xy ij* are approximation to the solution of the Navier–Stokes equations obtained after iterations of the basic method. Then the mass conservation equation can be

Fig. 14. Flow near the moving plunger

formulated as

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

the given problem, the auxiliary problem makes it possible to compute the most «part of pressure» (i.e. *px*(*t*, *x*) + *py*(*t*, *y*)) and corresponding change of thermophysical properties of the fluid based on simplified (pressure-unlinked) momentum equations in primitive variables

Previously fluid flows with the given mass rate have been simulated and analyzed. In some applications mass flow rate cannot be given in advance. The simplest example of such flows is problem about moving plunger. Immovable plunger is located between parallel plates filled by incompressible fluid. Figure 13 shows geometry of the problem and starting location of the plunger. Motion of the plunger causes the fluid flow. It is clear that mass flow rate depends on the plunger speed. Algorithm for simulation of the flows with unspecified mass rate should

Stage I: *main problem*: iterations of basic method (SIMPLE, Uzawa or Vanka iterations, etc.)

Stage III: *auxiliary problem*: several iterations of line (2D) or plane (3D) Seidel method with the

obtained after iterations of the basic method. Then the mass conservation equation can be

*ij* are approximation to the solution of the Navier–Stokes equations

Fig. 13. Starting location of the plunger

formulation and mass conservation equations.

Fig. 12. Isobars in Laval micronozzle (sizes in meters)

Stage II: check convergence; continue if necessary

*xy*

secant iterations; continue (go to 1)

Assume that *u*˜*ij*, *v*˜*ij* and *p*˜

**5.6 Flows with unspecied mass rate**

be modified as:

$$\int\_0^1 \mathfrak{u}(t, \mathfrak{x}\_{i\prime} y) \, dy = \max\_i \int\_0^1 \mathfrak{u}(t, \mathfrak{x}\_{i\prime} y) \, dy$$

at accelerated or uniform plunger motion and in form

$$\int\_0^1 \mathfrak{u}(t, \mathfrak{x}\_{i\prime} y) \, dy = \min\_i \int\_0^1 \mathfrak{u}(t, \mathfrak{x}\_{i\prime} y) \, dy$$

at decelerated plunger motion. As a result, the auxiliary problem ensures expansion of perturbations caused by the plunger on all domain in each iteration. Uniform staggered grid 201× 2401 (*hx* = *hy* = 1/200) is used for simulation the fluid flow around the moving plunger. Reynolds number Re = 200 is based on the plunger velocity and distance between the plates. Stream function isolines and isobars near the moving plunger are shown on Figure 14.

#### **6. Pressure decomposition in geometric multigrid methods**

Proposed convergence acceleration technique based on the pressure decomposition (10) should be incorporated with well-known algorithms for solving Navier–Stokes equations. Multigrid methods having (almost) optimal convergence rate for many applications seem to be the most promising solvers for many CFD problems. Our purpose is development of multigrid method with the least number of problem-dependent components for using in black box software.

#### **6.1 Nonlinear multigrid iterations**

To overcome problem of robustness, the Navier–Stokes equations (4) should be adapted for the multigrid algorithm Martynenko (2006). Adaptation of the Navier–Stokes equations (so called Σ-modification) consists in representation of the velocity *V* and pressure *P* as sum of

initial Navier–Stokes equations, i.e.

**6.2 Multigrid structure**

(*N* = 1, 2, 3) coarse grids *G*1

1. all coarse grids *G*1

2. the finest grid *G*0

finest grid *G*0

**6.3 Multigrid cycle**

1.

performed. Finally we obtain *Gl*

*Ru*(*t*, *x*, *y*) = − *∂u*ˆ

for Solving Navier–Stokes Equations on Structured Grids

*Rv*(*t*, *x*, *y*) = − *∂v*ˆ

Assume that a finest staggered grid *G*0

1, *G*1

as than the mesh size on the finest grid. 4. control volume on the coarse grids *G*1

coarse grid generation is given in Martynenko (2006).

1, *G*1

2,..., *G*1

*Ruv*(*t*, *x*, *y*) = − *∂u*ˆ

195 Convergence Acceleration of Iterative Algorithms

*∂x* − *∂*(*v*ˆ*u*ˆ)

*∂x* − *∂*(*v*ˆ2)

Additional convection terms in Σ-modified momentum equations are result of nonlinear nature of the convection-diffusion operator N , i.e. N ∗ �= N . Approximation of the source terms *Ruv*, *Ru* and *Rv* defines the accuracy, monotonicity and conservatism of the numerical solutions. Approximation of other terms in the modified equations defines only multigrid

Recently variant of the geometric multigrid methods with the problem-independent transfer operators (so-called Robust Multigrid Technique) has been proposed and developed Martynenko (2006; 2010). The problem-independent restriction and prolongation operators

3*N* with the following properties:

1, *G*1

2,..., *G*1

*k*, *l* = 0, 1, . . . , *L*+, *k* = 1, 3*Nl* computational grids (so called

3*N* have no common points, i.e. *G*1

3. all grids are similar to each other, but a mesh size on the coarse grids is three times as large

3*N* coarser grids. The coarse grid generation is finished when no further coarsening can be

a multigrid structure), where *L*+ is number of the coarsest level and *N* = 2, 3. Details of the

Multigrid cycles for simulation of flows with given and unspecified mass flow rates are shown on Figure 15. There are two kinds of smoothing on the multigrid structure. One of them (marked as ◦) is intended for solving Σ-modified Navier–Stokes equations with computation of correction of the «multidimensional component of pressure». Another smoothing (marked as ) is intended for computation of correction of the «one-dimensional component of

2,..., *G*1

1, *G*1

*∂t* − *∂*(*u*ˆ2)

*∂t* − *∂*(*u*ˆ*v*ˆ)

convergence rate because *cu* → 0, *cv* → 0 and *cp* → 0 for convergent solution.

are result of the multiple coarse grid corrections on subgrids of the finest grid.

Robust Multigrid Technique is based on representation of the finest grid *G*0

2,..., *G*1

1 is the union of all coarse grids *G*1

The coarse grid generation is further recurrently repeated: each grid *G*1

*∂x* − *∂v*ˆ *∂y* ,

> *∂*2*u*ˆ *∂y*2 ,

*∂*2*v*ˆ *∂y*2 ,

1 has been generated in the domain. Coarsening in

*n* ∩ *G*1

3*N* , i.e. *G*0

3*N* is union of 3*N* control volumes on the

1, *G*1

2,..., *G*1

3*N* gives

1 as union of 3*N*

*m* = ∅ , *n* �= *m*.

1 = 3*N k*=1 *G*1 *k* .

*∂y* − *∂p*ˆ *∂x* + 1 Re *∂*2*u*ˆ *∂x*2 +

*∂y* − *∂p*ˆ *∂y* + 1 Re *∂*2*v*ˆ *∂x*2 +

two functions

$$
\vec{V} = \mathsf{C}\_{\vec{V}} + \vec{\mathcal{V}}\_{\prime} \qquad P = \mathsf{C}\_{P} + \mathsf{P}\_{\prime\prime}
$$

where discrete analogues of the functions *CV* and *CP* will be coarse grid corrections and discrete analogues of the functions *V* ˆ and *P*ˆ will be approximations to the solutions in the following multigrid iterations.

As a result, the Navier–Stokes equations (4) can be rewritten in the Σ-modified form

$$\begin{cases} \mathcal{N}(\mathsf{C}\_{\vec{V}} + \vec{\hat{V}}) + \nabla(\mathsf{C}\_{P} + \hat{P}) = F \\ \qquad \nabla(\mathsf{C}\_{\vec{V}} + \vec{\hat{V}}) = G \end{cases}$$

Since N (*CV* + *V* ˆ) = N ∗(*CV* ) + N (*V* ˆ), where N ∗ �= N for the nonlinear operator N , we obtain

$$\begin{cases} \mathcal{N}^\*(\mathbb{C}\_{\vec{V}}) + \nabla \mathbb{C}\_P = F^\* \\ \qquad \nabla \mathbb{C}\_{\vec{V}} = G^\* \end{cases} \tag{26}$$

where *F*∗ = *F* − N (*V* ˆ) − ∇*P*ˆ and *G*∗ = *G* − ∇*V* ˆ .

It is clear that main difference between Σ-modified and initial forms of the Navier–Stokes equations consists of the nonlinear convection-diffusion operator (N ∗ in equation (26) instead of N in equation (4)) and source terms (*F*∗ and *G*∗ in equation (26) instead of *F* and *G* in equation (4)). Note that Σ-modification does not require some linearization of the Navier–Stokes equations. Therefore modified Navier–Stokes equations with other transport equations can be solved in coupled manner on all coarse grids.

For example, 2D Σ-modified Navier–Stokes equations are written as:

a) Σ-modified continuity equation

$$\frac{\partial \mathcal{C}^{\mu}}{\partial \mathbf{x}} + \frac{\partial \mathcal{C}^{\upsilon}}{\partial y} = \mathcal{R}^{\mu \upsilon}(t, \mathbf{x}, y)\_{\upsilon}$$

b) Σ-modified X-momentum

$$\begin{split} \frac{\partial \mathcal{L}^{\mu}}{\partial t} + \frac{\partial (\mathcal{c}^{\mu})^{2}}{\partial \mathbf{x}} + 2 \frac{\partial (\mathcal{C} \mathcal{c}^{\mu})}{\partial \mathbf{x}} + \frac{\partial (\mathcal{C} \mathcal{c}^{\upsilon})}{\partial \mathbf{y}} + \frac{\partial (\mathcal{C} \mathcal{c}^{\mu})}{\partial \mathbf{y}} + \frac{\partial (\mathcal{C}^{\mu} \mathcal{c}^{\upsilon})}{\partial \mathbf{y}} &= 0 \\ = -\frac{\partial \mathcal{c}^{p}}{\partial \mathbf{x}} + \frac{1}{\mathsf{Re}} \left( \frac{\partial^{2} \mathcal{c}^{\mu}}{\partial \mathbf{x}^{2}} + \frac{\partial^{2} \mathcal{c}^{\mu}}{\partial \mathbf{y}^{2}} \right) + \mathsf{R}^{\mu}(t, \mathbf{x}, \mathbf{y}) \end{split}$$

c) Σ-modified Y-momentum

$$\frac{\partial \mathbf{c}^{v}}{\partial t} + \frac{\partial (\hat{\mathbf{u}} \mathbf{c}^{v})}{\partial \mathbf{x}} + \frac{\partial (\hat{\mathbf{c}} \mathbf{c}^{\mu})}{\partial \mathbf{x}} + \frac{\partial (\mathbf{c}^{\mu} \mathbf{c}^{v})}{\partial \mathbf{x}} + \frac{\partial (\mathbf{c}^{v})^{2}}{\partial y} + 2 \frac{\partial (\hat{\mathbf{c}} \mathbf{c}^{v})}{\partial y} = 0$$

$$= -\frac{\partial \mathbf{c}^{p}}{\partial y} + \frac{1}{\mathbf{Re}} \left( \frac{\partial^{2} \mathbf{c}^{v}}{\partial \mathbf{x}^{2}} + \frac{\partial^{2} \mathbf{c}^{v}}{\partial y^{2}} \right) + \mathbf{R}^{v}(t, \mathbf{x}, y)$$

where discrete analogues of the functions *cu*, *cv* and *cp* will be coarse grid corrections and discrete analogues of the functions *u*ˆ, *v*ˆ and *p*ˆ will be approximations to the solutions in the following multigrid iterations. Source terms in the Σ-modified equations coincide with the initial Navier–Stokes equations, i.e.

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

where discrete analogues of the functions *CV* and *CP* will be coarse grid corrections and

� N ∗(*CV* ) + ∇*CP* = *F*∗

ˆ . It is clear that main difference between Σ-modified and initial forms of the Navier–Stokes equations consists of the nonlinear convection-diffusion operator (N ∗ in equation (26) instead of N in equation (4)) and source terms (*F*∗ and *G*∗ in equation (26) instead of *F* and *G* in equation (4)). Note that Σ-modification does not require some linearization of the Navier–Stokes equations. Therefore modified Navier–Stokes equations with other transport

*∂y* = *Ruv*(*t*, *x*, *y*),

*∂*(*uc*ˆ *v*) *∂y*

*∂*(*cucv*) *∂x*

where discrete analogues of the functions *cu*, *cv* and *cp* will be coarse grid corrections and discrete analogues of the functions *u*ˆ, *v*ˆ and *p*ˆ will be approximations to the solutions in the following multigrid iterations. Source terms in the Σ-modified equations coincide with the

*∂*2*cu ∂y*2

+

�

+

�

*∂*2*cv ∂y*2

*∂*(*vc*ˆ *u*) *∂y*

*∂*(*cv*)2 *∂y*

+

+ 2

+ *Rv*(*t*, *x*, *y*),

+ *Ru*(*t*, *x*, *y*),

*∂*(*cucv*) *∂y* =

*∂*(*vc*ˆ *v*) *∂y* =

As a result, the Navier–Stokes equations (4) can be rewritten in the Σ-modified form

ˆ , *P* = *CP* + *P*ˆ ,

ˆ) + ∇(*CP* + *P*ˆ) = *F*

∇(*CV* + *V*

ˆ) = *G* .

ˆ and *P*ˆ will be approximations to the solutions in the

ˆ), where N ∗ �= N for the nonlinear operator N , we

∇*CV* = *G*∗ , (26)

*V* = *CV* + *V*

N (*CV* + *V*

ˆ) − ∇*P*ˆ and *G*∗ = *G* − ∇*V*

equations can be solved in coupled manner on all coarse grids. For example, 2D Σ-modified Navier–Stokes equations are written as:

> *∂cu ∂x* + *∂cv*

+ 2

+

*∂*(*uc*ˆ *u*) *∂x*

*∂*(*vc*ˆ *u*) *∂x*

+

� *∂*2*cu ∂x*2 +

+

� *∂*2*cv ∂x*2 +

⎧ ⎨ ⎩

ˆ) = N ∗(*CV* ) + N (*V*

two functions

Since N (*CV* + *V*

where *F*∗ = *F* − N (*V*

a) Σ-modified continuity equation

*∂cu ∂t* +

*∂cv ∂t* + *∂*(*cu*)2 *∂x*

= − *∂cp ∂x* + 1 Re

*∂*(*uc*ˆ *v*) *∂x*

= − *∂cp ∂y* + 1 Re

b) Σ-modified X-momentum

c) Σ-modified Y-momentum

obtain

discrete analogues of the functions *V*

following multigrid iterations.

$$R^{\mu\upsilon}(t, \mathbf{x}, y) = -\frac{\partial \hat{\mathfrak{u}}}{\partial \mathbf{x}} - \frac{\partial \hat{\mathfrak{v}}}{\partial y},$$

$$R^{\mu}(t, \mathbf{x}, y) = -\frac{\partial \hat{\mathfrak{u}}}{\partial t} - \frac{\partial (\hat{\mathfrak{u}}^2)}{\partial \mathbf{x}} - \frac{\partial (\vartheta \hat{\mathfrak{u}})}{\partial y} - \frac{\partial \hat{\mathfrak{p}}}{\partial \mathbf{x}} + \frac{1}{\text{Re}} \left( \frac{\partial^2 \hat{\mathfrak{u}}}{\partial \mathbf{x}^2} + \frac{\partial^2 \hat{\mathfrak{u}}}{\partial y^2} \right),$$

$$R^{\upsilon}(t, \mathbf{x}, y) = -\frac{\partial \vartheta}{\partial t} - \frac{\partial (\vartheta \vartheta)}{\partial \mathbf{x}} - \frac{\partial (\vartheta^2)}{\partial y} - \frac{\partial \hat{\mathfrak{p}}}{\partial y} + \frac{1}{\text{Re}} \left( \frac{\partial^2 \vartheta}{\partial \mathbf{x}^2} + \frac{\partial^2 \vartheta}{\partial y^2} \right).$$

Additional convection terms in Σ-modified momentum equations are result of nonlinear nature of the convection-diffusion operator N , i.e. N ∗ �= N . Approximation of the source terms *Ruv*, *Ru* and *Rv* defines the accuracy, monotonicity and conservatism of the numerical solutions. Approximation of other terms in the modified equations defines only multigrid convergence rate because *cu* → 0, *cv* → 0 and *cp* → 0 for convergent solution.

#### **6.2 Multigrid structure**

Recently variant of the geometric multigrid methods with the problem-independent transfer operators (so-called Robust Multigrid Technique) has been proposed and developed Martynenko (2006; 2010). The problem-independent restriction and prolongation operators are result of the multiple coarse grid corrections on subgrids of the finest grid.

Assume that a finest staggered grid *G*0 1 has been generated in the domain. Coarsening in Robust Multigrid Technique is based on representation of the finest grid *G*0 1 as union of 3*N* (*N* = 1, 2, 3) coarse grids *G*1 1, *G*1 2,..., *G*1 3*N* with the following properties:

The coarse grid generation is further recurrently repeated: each grid *G*1 1, *G*1 2,..., *G*1 3*N* gives 3*N* coarser grids. The coarse grid generation is finished when no further coarsening can be performed. Finally we obtain *Gl k*, *l* = 0, 1, . . . , *L*+, *k* = 1, 3*Nl* computational grids (so called a multigrid structure), where *L*+ is number of the coarsest level and *N* = 2, 3. Details of the coarse grid generation is given in Martynenko (2006).

#### **6.3 Multigrid cycle**

Multigrid cycles for simulation of flows with given and unspecified mass flow rates are shown on Figure 15. There are two kinds of smoothing on the multigrid structure. One of them (marked as ◦) is intended for solving Σ-modified Navier–Stokes equations with computation of correction of the «multidimensional component of pressure». Another smoothing (marked as ) is intended for computation of correction of the «one-dimensional component of

**X**

**X**

(b) Isobars (*n* = 50)

**X**

(d) Isobars (*n* = 400)

**X**

(f) Isobars (*n* = 1000)

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

**Y**

0

**Y**

0

**Y**

0

0.2

0.4

0.6

0.8

1

0.2

0.4

0.6

0.8

1

0.2

0.4

0.6

0.8

1

197 Convergence Acceleration of Iterative Algorithms

(a) Isolines of stream function (*n* = 50)

**X**

(c) Isolines of stream function (*n* = 400)

**X**

(e) Isolines of stream function (*n* = 1000)

Fig. 16. Flow picture in the driven cavity (*n* = 50, 400, 1000)

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

for Solving Navier–Stokes Equations on Structured Grids

**Y**

0

**Y**

0

**Y**

0

0.2

0.4

0.6

0.8

1

0.2

0.4

0.6

0.8

1

0.2

0.4

0.6

0.8

1

pressure» on the finest grid. In other words solution of the modified auxiliary problem is additional smoothing on the finest grid.

(a) Flows with unspecified mass flow rate (b) Flows with the given mass flow rate

#### Fig. 15. Multigrid cycles

#### **6.4 Numerical test**

The algorithm is tested by simulation of unsteady flow in a driven cavity at Re = 1000. Law of the lid motion is taken as

$$\mathcal{U}\_w^{(n)} = \min\left(\frac{n}{50}; 1\right),$$

and stopping criterion is posed as

$$\max\_{ij} \left| \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right|\_{ij}^{(n+1)} < 10^{-10}.$$

Six-level multigrid structure with uniform finest staggered grid 501 × 501 (*ht* = *hx* = *hy* = 1/500) is used for the test. Correction of the «multidimensional component of pressure» is computed using Uzawa algorithm with diagonal preconditioning.

Figures 16–18 show evolution of the flow in the cavity. The main vortex is located near the lid after finish of the lid acceleration. Then the vortex moves to upper right corner under influence of the lid motion. After that the vortex moves along diagonal of the cavity to the center (Figure 16). Motion of the main vortex generates two additional vortices (Figure 17). The first vortex is formed in the lower corner, but the second vortex is formed on right vertical wall of the cavity. Location of the corner vortex is stable because of its motion is limited by the cavity walls. However the wall vortex can moving under influence of the main vortex. Corner and wall vortices agglomerate in common vortex. Agglomerated vortex tends to stable corner position, but a new corner vortex is generated in the left lower corner of the cavity. Figure 18 represents close-to-steady flow picture in the cavity.

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

pressure» on the finest grid. In other words solution of the modified auxiliary problem is

The algorithm is tested by simulation of unsteady flow in a driven cavity at Re = 1000. Law

 *n* 50 ; 1 ,

(*n*+1)

< 10−10.

*ij*

Six-level multigrid structure with uniform finest staggered grid 501 × 501 (*ht* = *hx* = *hy* = 1/500) is used for the test. Correction of the «multidimensional component of pressure» is

Figures 16–18 show evolution of the flow in the cavity. The main vortex is located near the lid after finish of the lid acceleration. Then the vortex moves to upper right corner under influence of the lid motion. After that the vortex moves along diagonal of the cavity to the center (Figure 16). Motion of the main vortex generates two additional vortices (Figure 17). The first vortex is formed in the lower corner, but the second vortex is formed on right vertical wall of the cavity. Location of the corner vortex is stable because of its motion is limited by the cavity walls. However the wall vortex can moving under influence of the main vortex. Corner and wall vortices agglomerate in common vortex. Agglomerated vortex tends to stable corner position, but a new corner vortex is generated in the left lower corner of the cavity. Figure 18

*U*(*n*) *w* = min

max *ij*

computed using Uzawa algorithm with diagonal preconditioning.

represents close-to-steady flow picture in the cavity.

 *∂u ∂x* + *∂v ∂y*

(b) Flows with the given mass flow rate

additional smoothing on the finest grid.

(a) Flows with unspecified mass flow

rate

Fig. 15. Multigrid cycles

of the lid motion is taken as

and stopping criterion is posed as

**6.4 Numerical test**

Fig. 16. Flow picture in the driven cavity (*n* = 50, 400, 1000)

**X**

**X**

(b) Isobars (*n* = 5000)

**X**

(d) Isobars (*n* = 10000)

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

**Y**

0

**Y**

0

Work supported by Russian Foundation for the Basic Research (project no. 09-01-00151). I wish to express a great appreciation to professor M.P. Galanin (Keldysh Institute of Applied Mathematics of Russian Academy of Sciences), who have guided and supported the

«Part of pressure» (i.e. sum of the «one-dimensional components» in decomposition (10)) can be computed using the simplified (pressure-unlinked) Navier–Stokes equations in primitive variables formulation and the mass conservation equations. «One-dimensional components of pressure» and corresponding velocity components are computed only in coupled manner. As a result, there are not pure segregated algorithms and pure density-based approach on structured grids. Proposed method does not require preconditioners and relaxation

0.2

0.4

0.6

0.8

1

0.2

0.4

0.6

0.8

1

199 Convergence Acceleration of Iterative Algorithms

(a) Isolines of stream function (*n* = 5000)

**X**

(c) Isolines of stream function (*n* = 10000)

Fig. 18. Flow picture in the driven cavity (*n* = 5000, 10000)

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

for Solving Navier–Stokes Equations on Structured Grids

**Y**

0

**Y**

0

**7. Acknowledgements**

researches.

**8. Conclusion**

0.2

0.4

0.6

0.8

1

0.2

0.4

0.6

0.8

1

Fig. 17. Flow picture in the driven cavity (*n* = 2250, 3000, 3500)

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

**Y**

0

**Y**

0

**Y**

0

0.2

0.4

0.6

0.8

1

0.2

0.4

0.6

0.8

1

0.2

0.4

0.6

0.8

1

**X**

(b) Isobars (*n* = 2250)

**X**

(d) Isobars (*n* = 3000)

**X**

(f) Isobars (*n* = 3500)

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

**X**

(a) Isolines of stream function (*n* = 2250)

**X**

(c) Isolines of stream function (*n* = 3000)

**X**

(e) Isolines of stream function (*n* = 3500)

Fig. 17. Flow picture in the driven cavity (*n* = 2250, 3000, 3500)

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

**Y**

0

**Y**

0

**Y**

0

0.2

0.4

0.6

0.8

1

0.2

0.4

0.6

0.8

1

0.2

0.4

0.6

0.8

1

Fig. 18. Flow picture in the driven cavity (*n* = 5000, 10000)

#### **7. Acknowledgements**

Work supported by Russian Foundation for the Basic Research (project no. 09-01-00151). I wish to express a great appreciation to professor M.P. Galanin (Keldysh Institute of Applied Mathematics of Russian Academy of Sciences), who have guided and supported the researches.

#### **8. Conclusion**

«Part of pressure» (i.e. sum of the «one-dimensional components» in decomposition (10)) can be computed using the simplified (pressure-unlinked) Navier–Stokes equations in primitive variables formulation and the mass conservation equations. «One-dimensional components of pressure» and corresponding velocity components are computed only in coupled manner. As a result, there are not pure segregated algorithms and pure density-based approach on structured grids. Proposed method does not require preconditioners and relaxation

**10**

*Russia*

**Neural Network Modeling of**

**Hydrodynamics Processes**

*Voronezh State Technical University*

Sergey Valyuhov, Alexander Kretinin and Alexander Burakov

Many of the computational methods for equation solving can be considered as methods of weighted residuals (MWR), based on the assumption of analytical idea for basic equation solving. Test function type determines MWR specific variety including collocation methods, least squares (RMS) and Galerkin's method. MWR algorithm realization is basically reduced to nonlinear programming which is solved by minimizing the total equations residual by selecting the parameters of test solution. In this case, accuracy of solving using the MWR is defined by approximating test function properties, along with degree of its conformity with its initial partial differential equations, representing a continuum solution of mathematical

On fig. 1, computing artificial neural network (ANN) is presented in graphic form, illustrating process of intra-network computations. The input signals or the values of input

**1. Introduction**

physics equations.

Fig. 1. Neural network computing structure

parameters. Pressure decomposition is very efficient acceleration technique for simulation of directed fluid flows.

#### **9. References**

## **Neural Network Modeling of Hydrodynamics Processes**

Sergey Valyuhov, Alexander Kretinin and Alexander Burakov *Voronezh State Technical University Russia*

#### **1. Introduction**

26 Will-be-set-by-IN-TECH

200 Hydrodynamics – Optimizing Methods and Tools

parameters. Pressure decomposition is very efficient acceleration technique for simulation

Barton I.E. (1997) The entrance effect of laminar flow over a backward-facing step geometry,

Benzi, M.; Golub, G.H.; Liesen, J. (2006) Numerical solution of saddle point problems, *Acta*

Briley, W.R. (1974) Numerical method for predicting three-dimensional steady viscous flow in

Gartling D. (1990) A test problem for outflow boundary conditions-flow over a backward-facing step, *Int. J. for Num. Meth. in Fluids*, Vol. 11, pp. 953-967. Ghia, U.; Ghia, K.N.; Shin, C.T. (1982) High-Re solutions for incompressible flow using the

Gresho, P.M.; Gartling, D.K.; Torczynski, J.R.; Cliffe, K.A.; Winters, K.H.; Garratt, T.G.; Spence,

Keskar, J.; Lin, D.A. (1999) Computation of laminar backward-facing step flow at Re=800 with

Martynenko, S.I. (2006) Robust Multigrid Technique for black box software, *Comp. Meth. in*

Martynenko, S.I. (2009) A physical approach to development of numerical methods for solving

Martynenko, S.I. (2010) Potentialities of the Robust Multigrid Technique, *Comp. Meth. in Appl.*

Vanka S.P. (1986) Block-implicit multigrid solution of Navier–Stokes equations in primitive

Navier-Stokes equations and a multigrid method, *J. Comp. Phys.*, Vol. 48, pp.387-411.

A.; Goodrich, J.W. (1993) Is a steady viscous incompressible two-dimensional flow over a backward-facing step at Re=800 stable? *Int. J. for Num. Meth. in Fluids*, Vol. 17,

a spectral domain decomposition method, *Int. J. for Num. Meth. in Fluids*, Vol. 29,

Navier-Stokes equations in primitive variables formulation, *Int. J. of Comp. Science and*

*Int. J. for Num. Meth. in Fluids*, Vol. 25, pp. 633-644.

ducts, *J. Comp. Phys.*, Vol. 14, pp.8-28.

*Appl. Math.*, Vol. 6, No. 4, pp.413-435.

variables, *J. Comp. Phys.*, Vol. 65, pp.138-158.

Wesseling, P. (1991) *An Introduction to Multigrid Methods*, Wiley, Chichester.

*Math.*, Vol. 2, No. 4, pp.291-307.

*Math.*, Vol. 10, No. 1, pp.87-94.

of directed fluid flows.

*Numerica*, pp. 1-137.

pp. 501-541.

pp. 411-427.

**9. References**

Many of the computational methods for equation solving can be considered as methods of weighted residuals (MWR), based on the assumption of analytical idea for basic equation solving. Test function type determines MWR specific variety including collocation methods, least squares (RMS) and Galerkin's method. MWR algorithm realization is basically reduced to nonlinear programming which is solved by minimizing the total equations residual by selecting the parameters of test solution. In this case, accuracy of solving using the MWR is defined by approximating test function properties, along with degree of its conformity with its initial partial differential equations, representing a continuum solution of mathematical physics equations.

On fig. 1, computing artificial neural network (ANN) is presented in graphic form, illustrating process of intra-network computations. The input signals or the values of input

Fig. 1. Neural network computing structure

Neural Network Modeling of Hydrodynamics Processes 203

parameters. Thus, ANN training under outlet functional is composed of inner parameters definition using trial solution (1) for meeting the equation (3) goal and its solution is realized

Usually total squared error at net outlets is presented as an objective function at neural net training and an argument is the difference between the resulted '*s*' net outlet and the real value that is known a priori. This approach to neural net utilization is generally applied to the problems of statistical set transformation along with definition of those function values unknown a priori (net outlet) from argument (net inlet). As for simulation issues, they refer to mathematical representation of the laws of physics, along with its modification to be applied practically. It is usually related to necessity for developing a digital description of the process to be modeled. Under such conditions we will have to exclude the a priori known computation result from the objective function and its functional task. Objective function during the known law simulation, therefore, shall only be defined by inlet data and

1 2

Use of neuronet's method of weighted residuals (NMWR) requires having preliminary systematic study for each specific case, aimed at: 1) defining the number of calculation nodes (i.e. the calculation grid size); 2) defining number of neurons within the network, required for obtaining proper approximation power; 3) choosing initial approximations for training neural network test solution ; 4) selecting additional criteria in the goal function for training procedure regularization in order to avoid possible solution non-uniformity; 5) analyzing the possibilities for applying multi-criteria optimization algorithms to search neural network solution parameters (provided that several optimization criteria are

Artificial neural network used for hydrodynamic processes studying is presented by two fundamentally different approaches. The first is the NMWR used for direct differential hydrodynamics equations solution. The NMWR description and its example realization for Navier-Stokes equations solution is presented in papers (Kretinin, 2006; Kretinin et al., 2008). These equations describe the 2D laminar isothermal flow of viscous incompressible liquid. In the paper (Stogney & Kretinin, 2005), the NMWR is used for simulating flows within a channel with permeable wall. Neural network solution results of hydrodynamic equations for the computational zone consisting of two sub-domains are presented below. One is rotating, while another is immobile. In this case, for NMWR algorithm realization

specifying the conjugate conditions at the two sub-domains border is not required.

in pipeline in the event of medium leakage through the wall hole.

**2.1 NMWR application: Preliminary studying**

In the second approach, neural network structures are applied to computational experiment results approximation obtained by using traditional methods of computational hydrodynamics and for obtaining of hydrodynamic processes multifactor approximation models. This approach is illustrated by hydrodynamics processes neural network modeling

There are specific ANN training programs such as STATISTICA NEURAL NETWORKS or NEURAL TOOLBOX in the medium of MATLAB, adjusting the parameters of the network

*s s*

*E yf* **x** . (4)

2

*S*

through the corresponding modification of functional quality equation (2) training.

law simulated:

available).

variables are distributed and "move" along the connections of the corresponding input together with all the neurons of hidden layer. The signals may be amplified or weakened by being multiplied by corresponding coefficient (weight or connection). Signals coming to a certain neuron within the hidden layer are summed up and subjected to nonlinear transformation using so-called activation function. The signals further proceed to network outputs that can be multiple. In this case the signal is also multiplied by a certain weight value, i.e. sum of neuron output weight values within the hidden layer as a result of neural network operation. Artificial neural networks of similar structure are capable for universal approximation, making possible to approximate arbitrary continuous function with any required accuracy.

To analyze ANN approximation capabilities, perceptron with single hidden layer (SLP) was chosen as a basic model performing a nonlinear transformation from input space to output space by using the formula (Bishop, 1995):

$$y(\mathbf{w}, \mathbf{x}) = \sum\_{i=1}^{q} v\_i f\_{\sigma} \left( b\_i + \sum\_{j=1}^{n} w\_{ij} x\_j \right) + b\_0 \,\,\,\tag{1}$$

where *n* **x R** is network input vector, comprised of *j x* values; *q* – the neuron number of the single hidden layer; *s* **w R** – all weights and network thresholds vector; *w – ij* weight entering the model nonlinearly between *j*-m input and *i*-m neuron of the hidden layer; *i v* – output layer neuron weight corresponding to the *i*-neuron of the hidden layer; 0 , *i b b* – thresholds of neurons of the hidden layer and output neuron; *f*σ – activation function (in our case the logistic sigmoid is used). ANN of this structure already has the universal approximation capability, in other words it gives the opportunity to approximate the arbitrary analog function with any given accuracy. The main stage of using ANN for resolving of practical issues is the neural network model training, which is the process of the network weight iterative adjustment on the basis of the learning set (sample) , , , 1,..., *n ii i* **x xR** *y ik* in order to minimize the network error – quality functional

$$J(\mathbf{w}) = \sum\_{i=1}^{k} Q\left(f\_{\varepsilon}(\mathbf{w}, i)\right),\tag{2}$$

where **w** – ANN weight vector; 2 *Q f* ( ,) , *i f i* **w w** – ANN quality criterion as per the *i*training example; *f* **w wx** , , *i y y i i* – *i*-example error. For training purposes the statistically distributed approximation algorithms may be used based on the back error propagation or the numerical methods of the differentiable function optimization.

#### **2. Neuronet's method of weighted residuals for computer simulation of hydrodynamics problems**

Let us consider that a certain equation with exact solution *y*( ) *x*

$$L(\overline{y}) = 0 \tag{3}$$

for non-numeric value *ys* equation (3) presents an arbitrary **xs** within the learning sample. We have *L*(*y*)*=R* with substitution of approximate solution (1) into equation (3), with *R* as equation residual. *R* is continuous function *R=f*(**w,x**), being a function of SLP inner 202 Hydrodynamics – Optimizing Methods and Tools

variables are distributed and "move" along the connections of the corresponding input together with all the neurons of hidden layer. The signals may be amplified or weakened by being multiplied by corresponding coefficient (weight or connection). Signals coming to a certain neuron within the hidden layer are summed up and subjected to nonlinear transformation using so-called activation function. The signals further proceed to network outputs that can be multiple. In this case the signal is also multiplied by a certain weight value, i.e. sum of neuron output weight values within the hidden layer as a result of neural network operation. Artificial neural networks of similar structure are capable for universal approximation, making possible to approximate arbitrary continuous function with any

To analyze ANN approximation capabilities, perceptron with single hidden layer (SLP) was chosen as a basic model performing a nonlinear transformation from input space to output

1 1

where *n* **x R** is network input vector, comprised of *j x* values; *q* – the neuron number of the single hidden layer; *s* **w R** – all weights and network thresholds vector; *w – ij* weight entering the model nonlinearly between *j*-m input and *i*-m neuron of the hidden layer; *i v* – output layer neuron weight corresponding to the *i*-neuron of the hidden layer; 0 , *i b b* – thresholds of neurons of the hidden layer and output neuron; *f*σ – activation function (in our case the logistic sigmoid is used). ANN of this structure already has the universal approximation capability, in other words it gives the opportunity to approximate the arbitrary analog function with any given accuracy. The main stage of using ANN for resolving of practical issues is the neural network model training, which is the process of the network weight iterative adjustment on the basis of the learning set (sample)

*i j y vf b wx b*

*ii i* **x xR** *y ik* in order to minimize the network error – quality functional

1 ( ) ( ,) *k*

>

statistically distributed approximation algorithms may be used based on the back error

for non-numeric value *ys* equation (3) presents an arbitrary **xs** within the learning sample. We have *L*(*y*)*=R* with substitution of approximate solution (1) into equation (3), with *R* as equation residual. *R* is continuous function *R=f*(**w,x**), being a function of SLP inner

*i J Qf i*

propagation or the numerical methods of the differentiable function optimization.

**2. Neuronet's method of weighted residuals for computer simulation of**

*i i ij j*

**w wx** , , *i y y i i* – *i*-example error. For training purposes the

*q n*

0

**w w** , (2)

**w w** – ANN quality criterion as per the *i*-

*L y*() 0 (3)

**w x** , (1)

required accuracy.

, , , 1,..., *n*

training example; *f*

**hydrodynamics problems**

space by using the formula (Bishop, 1995):

( ,)

where **w** – ANN weight vector; 2 *Q f* ( ,) , *i f i*

Let us consider that a certain equation with exact solution *y*( ) *x*

parameters. Thus, ANN training under outlet functional is composed of inner parameters definition using trial solution (1) for meeting the equation (3) goal and its solution is realized through the corresponding modification of functional quality equation (2) training.

Usually total squared error at net outlets is presented as an objective function at neural net training and an argument is the difference between the resulted '*s*' net outlet and the real value that is known a priori. This approach to neural net utilization is generally applied to the problems of statistical set transformation along with definition of those function values unknown a priori (net outlet) from argument (net inlet). As for simulation issues, they refer to mathematical representation of the laws of physics, along with its modification to be applied practically. It is usually related to necessity for developing a digital description of the process to be modeled. Under such conditions we will have to exclude the a priori known computation result from the objective function and its functional task. Objective function during the known law simulation, therefore, shall only be defined by inlet data and law simulated:

$$E = \frac{1}{2} \sum\_{\mathbf{s}} \left( y^s - f\left(\mathbf{x}^s\right) \right)^2. \tag{4}$$

Use of neuronet's method of weighted residuals (NMWR) requires having preliminary systematic study for each specific case, aimed at: 1) defining the number of calculation nodes (i.e. the calculation grid size); 2) defining number of neurons within the network, required for obtaining proper approximation power; 3) choosing initial approximations for training neural network test solution ; 4) selecting additional criteria in the goal function for training procedure regularization in order to avoid possible solution non-uniformity; 5) analyzing the possibilities for applying multi-criteria optimization algorithms to search neural network solution parameters (provided that several optimization criteria are available).

Artificial neural network used for hydrodynamic processes studying is presented by two fundamentally different approaches. The first is the NMWR used for direct differential hydrodynamics equations solution. The NMWR description and its example realization for Navier-Stokes equations solution is presented in papers (Kretinin, 2006; Kretinin et al., 2008). These equations describe the 2D laminar isothermal flow of viscous incompressible liquid. In the paper (Stogney & Kretinin, 2005), the NMWR is used for simulating flows within a channel with permeable wall. Neural network solution results of hydrodynamic equations for the computational zone consisting of two sub-domains are presented below. One is rotating, while another is immobile. In this case, for NMWR algorithm realization specifying the conjugate conditions at the two sub-domains border is not required.

In the second approach, neural network structures are applied to computational experiment results approximation obtained by using traditional methods of computational hydrodynamics and for obtaining of hydrodynamic processes multifactor approximation models. This approach is illustrated by hydrodynamics processes neural network modeling in pipeline in the event of medium leakage through the wall hole.

#### **2.1 NMWR application: Preliminary studying**

There are specific ANN training programs such as STATISTICA NEURAL NETWORKS or NEURAL TOOLBOX in the medium of MATLAB, adjusting the parameters of the network

Neural Network Modeling of Hydrodynamics Processes 205

With IOSO NS 1.0, ANN internal parameter values were obtained with sum residual there were received the values of the internal parameters of the ANN, giving the sum

common vs

 do i=1,3 w(i)=vs(i) b(i)=vs(i+3) v(i)=vs(i+6) end do bv=vs(10) vyh=0. do i=1,3 t1(i)=w(i)\*t-b(i) q(i)=1./(1.+exp(-t1(i))) vyh=vyh+v(i)\*q(i)

 end do ynet=vyh-bv

*E* 0.000023 (fig. 2).

end

c b-thresholds of neurons

c bv-threshold of output neuron

c w-weights between neuron and input

c v-weights between neuron and output neuron

Fig. 2. Results of using IOSO NS 1.0 for the ANN training

to the known values of the objective function within the given points of its definitional domain. Using these packages in our case, therefore, does not seem possible. At the same time, many of optimization standard methods work well for ANN training, e.g. the conjugate gradients methods, or Newton, etc. To solve the issue of ANN training, we shall use the Russian program IOSO NS 1.0 (designed by prof. I.N. Egorov (Egorov et al., 1998), see www.IOSOTech.com) realizing the algorithm of indirect optimization method based on self-organizing. This program allows minimizing the mathematical model given algorithmically and presented as "black box", i.e. as external file module which scans its values from running variable file generated by optimization program, then calculates objective function value and records it in the output file, addressed in turn by optimization program. It is therefore sufficient for computer program forming, realizing calculations using the required neural network, where the input data will be network internal parameters (i.e. weights, thresholds); on the output, however, there'll be value of required equation sum residual based on accounting area free points. Let us suppose that the objective function 2 *y x* is determined within the interval 0;1 . It is necessary to define parameters of ANN perceptron type with one hidden layer, consisting of 3 neurons to draw the near-objective function with given accuracy, computed in 100 accounting points *i x* evenly portioned in determination field. Computer program for computing network sum residual depending on its parameters can be as follows (Fortran):

common vs

204 Hydrodynamics – Optimizing Methods and Tools

to the known values of the objective function within the given points of its definitional domain. Using these packages in our case, therefore, does not seem possible. At the same time, many of optimization standard methods work well for ANN training, e.g. the conjugate gradients methods, or Newton, etc. To solve the issue of ANN training, we shall use the Russian program IOSO NS 1.0 (designed by prof. I.N. Egorov (Egorov et al., 1998), see www.IOSOTech.com) realizing the algorithm of indirect optimization method based on self-organizing. This program allows minimizing the mathematical model given algorithmically and presented as "black box", i.e. as external file module which scans its values from running variable file generated by optimization program, then calculates objective function value and records it in the output file, addressed in turn by optimization program. It is therefore sufficient for computer program forming, realizing calculations using the required neural network, where the input data will be network internal parameters (i.e. weights, thresholds); on the output, however, there'll be value of required equation sum residual based on accounting area free points. Let us suppose that the objective function 2 *y x* is determined within the interval 0;1 . It is necessary to define parameters of ANN perceptron type with one hidden layer, consisting of 3 neurons to draw the near-objective function with given accuracy, computed in 100 accounting points *i x* evenly portioned in determination field. Computer program for computing network sum

residual depending on its parameters can be as follows (Fortran):

 dimension x(100),y(100) dimension vs(10) common vs

c 'inp'- file of input data,

del=del+(y(i)-x(i)\*\*2)\*\*2

 open(2,file='out') write(2,\*)del close(2) end

function ynet(t)

c sent to optimization program

 open(1,file='inp') read(1,\*)vs close(1)

 do i=1,100 x(i)=(i-1)/99. end do

 del=0. do i=1,100 y(i)=ynet(x(i))

end do

c vs- values of ANN internal parameters

c generated by optimization program

c calculation by subprogram ANN ynet c and finding of sum residual del

c 'out'-file of value of the minimization function ,

dimension vs(10),w(3),b(3),v(3),t1(3),q(3)

```
 do i=1,3
 w(i)=vs(i)
 b(i)=vs(i+3)
 v(i)=vs(i+6)
 end do
 bv=vs(10)
 vyh=0.
 do i=1,3
 t1(i)=w(i)*t-b(i)
 q(i)=1./(1.+exp(-t1(i)))
 vyh=vyh+v(i)*q(i)
 end do
 ynet=vyh-bv
 end
```
With IOSO NS 1.0, ANN internal parameter values were obtained with sum residual there were received the values of the internal parameters of the ANN, giving the sum *E* 0.000023 (fig. 2).

Fig. 2. Results of using IOSO NS 1.0 for the ANN training

Neural Network Modeling of Hydrodynamics Processes 207

*E y f w b bx*

,1 , *s ss s s*

Thereby, we have got all the components of the gradient of the objective function of minimization, comparatively which iterations will be consecutively realized in accordance

> **w w**

Parameter optimization of neural network trial solutions is achieved by applying several optimization strategies and by subsequently choosing the maximum effective one (see Cloete & Zurada, 2000). First strategy is to apply totality of effective gradient methods "starting" from various initial points. The other strategy is to apply structural-parametrical optimization to ANN training; this method is based on indirect statistic optimization method on self-organizing basis or parameter space research (see: Egorov et al., 1998;