**1. Introduction**

Mathematical modeling of hydrodynamics is the base for research of various natural phenomena, technological processes, and environmental problems. The main equations describing this problem are the Navier–Stokes equations. Development and research of effective numerical algorithms for solving these equations and their practical realization is an actual task. The use of the MGM for the numerical solution of the Navier–Stokes equations describing the motion of an incompressible viscous fluid is discussed. Currently, various discretization methods for the corresponding differential model are known. However, with any choice of the discretizing method, the problem of constructing effective methods for solving large systems of algebraic equations—to which the discrete model is reduced—arises. This problem is especially relevant in the nonstationary case, when multiple solutions of the systems of algebraic equations are required at each discrete time step.

To discretize the system of two-dimensional Navier–Stokes equations on regular grids, we use the finite difference method. The equations are considered in the natural variables "velocity-pressure":

$$\frac{\partial \mathbf{V}}{\partial t} + (\mathbf{V} \cdot \nabla)\mathbf{V} = -\nabla P + \nu \Delta \mathbf{V}, \quad \operatorname{div} \mathbf{V} = 0,\tag{1}$$

where *P*/*ρ* is replaced by *P* (i.e., *ρ* is normalized at 1), *P* is the static pressure, **V** is the velocity vector, and *ν* is the kinematic viscosity coefficient. At the initial moment of time and at the boundary of the domain, the initial and boundary conditions are set, respectively.

Flow simulating is accompanied by a number of mathematical difficulties. One of the problems in solving this system is the nonlinearity associated with convective terms in the equations, which can lead to the appearance of oscillations of the solution in regions with large gradients. The main efforts of the researchers were directed at overcoming the difficulties associated with the nonlinearity of the Navier–Stokes system of equations.

One of the most time-consuming stages of the computational procedure is finding the solution of the system of linear algebraic equations (SLAE). Modern application packages usually use the linearization of the original equations, and Krylov subspace methods are used to solve the resulting SLAEs. Despite the fact that these methods have proven themselves well, they have some problems in cases of significant nonsymmetry of the SLAEs—associated, for example, with variable coefficients in differential equations or using complex numerical boundary conditions. For time discretization of the unsteady problem, we use an implicit difference scheme. Here, we do not specifically consider the stages of discretization and linearization of the Navier–Stokes equations, but focus on solving SLAEs. Given that the SLAEs resulting from the use of the implicit time schemes have a large dimension and a sparse nonsymmetric matrix, we propose using the MGM to solve them.

Thus, we consider the iterative solution of the large sparse SLAE

$$Av = b, \quad v, b \in \mathbb{C}^n,\tag{2}$$

where *<sup>A</sup>* <sup>∈</sup> <sup>C</sup>*n*×*<sup>n</sup>* is a non-Hermitian and positive definite matrix.

Naturally, the matrix *A* can be split as

$$A = A\_0 + A\_{1\prime} \tag{3}$$

where

$$A\_0 = \frac{1}{2}(A + A^\*), \quad A\_1 = \frac{1}{2}(A - A^\*) \tag{4}$$

and *A*∗ denotes the conjugate transpose of the matrix *A*. Positive definiteness of the matrix *A* means that for all *<sup>x</sup>* <sup>∈</sup> <sup>C</sup>*<sup>n</sup>* \ {0}, *<sup>x</sup>*∗*A*0*<sup>x</sup>* <sup>&</sup>gt; 0. Here, *<sup>x</sup>*<sup>∗</sup> denotes the conjugate transpose of the complex vector *x*. Let in some matrix norm ||| · |||, |||*A*0||| << |||*A*1|||, then the matrix *A* is called a strongly non-Hermitian one. This situation occurs in many real applications, such as the discretization of the Navier–Stokes equations.

The Hermitian and skew-Hermitian splitting (HSS) iteration methods, based on HS splitting (3) and (4), for solving large sparse non-Hermitian positive definite SLAE were firstly proposed in [1]. The HSS iteration method has been widely developed in [2–5] and others.

Then, we can split the skew-Hermitian part *<sup>A</sup>*<sup>1</sup> of the matrix *<sup>A</sup>* <sup>∈</sup> <sup>C</sup>*n*×*<sup>n</sup>* into

$$A\_1 = \mathcal{K}\_L + \mathcal{K}\_{lI\prime} \tag{5}$$

where *KL* and *KU* are the strictly lower and the strictly upper triangular parts of *A*1, respectively. Obviously, that *KL* = −*K*<sup>∗</sup> *U*.

Based on the splitting (3)–(5) in [6–8] classes of skew-Hermitian triangular splitting (STS), iteration methods for solving SLAE (2) have been proposed. The triangular operator of the STS uses only the skew-Hermitian part of the coefficient matrix *A*. These methods have been further developed in [9–12].

The use of the multigrid method (MGM) with the STS-based smoothers for solving convection–diffusion problems has been studied in [13]. The convergence of the MGM with the STS-based smoothers has also been proved in this research. The local Fourier analysis of the MGM with the triangular skew-symmetric smoothers has been performed in [14]. The results of numerical experiments for convection–diffusion problems with large Peclet numbers by the geometric MGM have been presented in both researches.

In [15], it was shown that the MGMs with the HSS-based smoothers converge uniformly for second-order nonselfadjoint elliptic boundary value problems. This happens if the mesh size of the coarsest grid is sufficiently small, but independent of the number of the multigrid levels.
