Revision as of 12:27, 18 December 2020

Abstract

Among the main problems in the field of Computational Fluid Dynamics (CFD), there is the two-dimensional laminar flow in a transient regime of an incompressible fluid modeled by Navier-Stokes equations. Among the decoupled solutions for this equation system, that is, solutions for velocity regardless to pressure, there are the Projection methods, which separate the solution in three parts, solved in each time steps applying an iterative process regardless to the iterative process for the other variables. It may result, according to the Projection method applied, in two Reaction-Diffusion equations and one Poisson equation to be solved in each time step. This paper sought to develop an algorithm to solve the Navier-Stokes equation, applying the Finite Volume Method (FVM) with second order approximation scheme (CDS), beside a Projection method with incremental pressure-correction scheme, so that each Reaction-Diffusion and the Poisson equation are solved efficiently. Therefore, several solvers were tested for each equation, resulting in an algorithm with the combination that achieved the best result for each equation, with the preconditioned Conjugate Gradient method (PCG) with the Multigrid method (MG) and ILU solver (Incomplete LU factorization) being the methodology used in the whole problem solving process. The geometric Multigrid with V cycle, the correction scheme (CS), the full weighting restriction, the prolongation through bilinear interpolation and the maximum number of levels for the studied cases were utilized. The results achieved were satisfactory, since the proposed methodology accelerated the iterative process considerably in relation to the classical methods available in the literature.

Keywords: Navier-Stokes, multigrid, conjugated gradient, projection methods, ILU factorization

1. Introduction

In computational mathematics, determining the solution for linear or nonlinear equation systems is an important problem, since they often occur in the discretization process of mathematical models in Computational Fluid Dynamics (CFD). In order to achieve such solutions, iterative methods are widely applied. They can be categorized as linear, nonlinear, stationary and nonstationary. These approaches may result in difficulties related to the slow convergence of the iterative process applied [1]. In order to improve the convergence of the iterative methods, the Multigrid method has been proving very efficient [2-6]. Its philosophy is based on the application of several grids with several degrees of refinement, which are gone through during the iterative process. Nevertheless, the conditioning of the matrix for the linear system's coefficients may slow the iterative method's convergence. It may be accelerated by modifying this coefficient matrix in order to improve its conditioning. This process is called matrix preconditioning [7]. Among the preconditioning technics, one can mention the incomplete LU factorization (ILU), which consists in changing a matrix A into a product of a lower triangular matrix L and an upper triangular matrix U, attached to a matrix called residual R.

In the study of incompressible fluids in transient regime, the Navier-Stokes equations are well established by the conservation laws [8-13]. In order to know the state of a fluid, one has to determine the value of the variables that identify it through the time for each point in space occupied by the fluid. The variables that identify the state of an incompressible fluid in transient and isothermal regime are the velocity and pressure in each point in space.

In order to overcome the dependence between pressure and velocity in the Navier-Stokes equations, several numerical schemes have been developed, remarkably the Projection method [14], which split the solution in three parts, solved in each time step applying a local iterative process. These methods may be categorized in three classes: pressure-correction schemes, velocity-correction schemes and consistent splitting schemes.

In this paper, the problem of a two dimensional laminar flow in a transient regime of an incompressible fluid modeled by the Navier-Stokes equations with Dirichlet boundary conditions for the velocities was solved numerically. A Projection method with incremental pressure-correction scheme was applied, in which the strategies to speed the solutions of the linear systems from these equations were mixed, such as the preconditioned Krylov Subspace method (Conjugated Gradient method) with Multigrid and ILU solver. The behaviors of some computational parameters that measure the convergence ratio of the methods and norms regarding the errors in numerical solutions were analyzed, in order to produce the algorithm with the best methodology among the ones that were studied.

2. Mathematical model and solution method

The mathematical model in two dimensional Cartesian coordinates for the laminar flow of an incompressible fluid, isothermal and in a transient regime, given by the Navier-Stokes equations [16], in this paper, is represented by the continuity equation that depicts the physical principle of the conservation of mass, and by the equations of dimensionless movement in the directions x and y:

where x and y are the coordinates; t the temporal coordinate, u and v the velocity components of the vector ${\displaystyle {\bf {u}}}$ in the directions x and y, respectively; p the static fluid pressure and Re the Reynolds number associated to the viscous and the inertia forces. According to Hughes and Brighton, when Re << 1, the viscous forces prevail and when Re >> 1 the inertia forces prevail [17].

For this set of equations, the Taylor-Green Vortices problem [18] was solved, whose domain is given by ${\displaystyle \left\{\left(x,y\right)\in R^{2}:-\pi \leq x\leq \pi {\mbox{ and }}-\right.}$${\displaystyle \left.\pi \leq y\leq \pi \right\}}$ , and the analytical solution by:

where ${\displaystyle T_{i}}$ is the initial time and ${\displaystyle t>0}$.

The initial and boundary conditions are obtained directly from the analytical solution.

The methods for solution of the Navier-Stokes equations can be generally categorized in coupled methods and segregated methods [13,14].

Coupled methods seek to solve the complete equation system for each computational cycle, coupling the equations of continuity and movement conservation. This is the most immediate way to solve the Navier-Stokes equations, but it renders grater difficulties in its implementation and a higher computational cost, due to the strong influence of the non-linearity convective terms [13].

Segregated methods seek the decoupling of the equations, separating the nonlinear system in simpler problems, which can be solved sequentially. Among these methods, the Projection methods [14] stand out.

Basically, the main idea of the Projection method for the Navier-Stokes Eqs. (1) to (3) is to apply the movement equation to determine a temporary velocity field.

According to Ernst and Gander [19], Eqs. (5) and (6) are Helmholtz equations

where k and η are constants that depend on the studied problem and f the source term.

Equation (6) is common in CFD and describes stationary Reaction-Diffusion phenomena, being the model to be solved in order to obtain the provisional velocity field. Despite the similarities, Eq. (6) is a much greater challenge in the classic iterative methods application [19].

After having obtained this temporary field, for the two dimensional case, it can be shown that a system of three equations and three unknowns (u, v and p) can be obtained, where the pressure appears in only one of the equations. The role of pressure in incompressible flows is to make the velocity field satisfy the continuity equation. Therefore, it is necessary to calculate pressure in the following step. For this purpose, an elliptical Poisson equation is solved for the continuity equation to be satisfied and for the pressure to be determined [14,15].

The Poisson equation is a elliptic partial differential equation. with wide usage in incompressible fluid flows, described by [20]:

The importance of this equation for pressure is that it links the movement conservation and the continuity equations. Thus, the final velocity field and pressure are calculated.

3. Numerical models

3.1 Utilized grid

In this paper the Finite Volume method (FVM) with staggered grids was applied, as described in [21], also applied by [12], with the velocities placed on the faces and the pressure on the center of the volumes. Applying this strategy avoids the numerical instabilities on the pressure solution [6,22,23]. Such grid and its ordering are represented in Figure 1 (4 to 4 volumes).

Figure 1. Grid and lexicographic order for the pressure and velocities variables

In this grid, the volumes with dashed lines are the volumes with ghost volumes (which do not belong to the physical domain of the problem), and it may be noticed that in applying this kind of grid, the boundary conditions for the variable u are automatically prescribed in the east and west boundaries, while in the other boundaries some kind of procedure is necessary in order for the ghost volumes to consider the contour information. Here, linear extrapolation was applied, and a similar case for the v variable in the south and north contours.

3.2 Spatial discretization by finite volumes

Equations (1) a (3) can be integrated in each control volume described in Figure 2.

Figure 2. Control volumes

Thus, the mass conservation equations, given by Eq. (1), of linear momentum in the direction x, Eq. (2), and linear momentum in the direction y, Eq. (3), can be written as:

and

respectively.

The terms of Eqs. (8), (9) and (10) were discretized applying the FVM with staggered grids, where the integrals are calculated over the control volumes defined in Figure 2, where for a matter of simplicity, the same spatial refinement was applied in all directions, that is, ${\displaystyle h_{x}=h_{y}=h}$ (uniform grids).

In order to discretized the mass conservation, a volume centered in p was applied, Figure 2 (b), where there is:

The pressure terms were discretized applying the volumes centered in the velocity components defined in Figures 2(a) and (c), generating:

and

The same control volumes defined in Figures 2(a) and (c) were applied in the discretization of the diffusive terms, where one obtains:

and

Furthermore, applying control volumes centered in the velocity components, the advective terms were discretized, resulting in:

and

The velocities that appear in Eqs. (16) and (17) can be approximated by:

and

3.3 Projection methods

The Projection methods were introduced by [24,25], and they are characterized by the Navier-Stokes solution in two steps. Initially, an auxiliary velocity ${\displaystyle {\bf {u}}_{t}}$ is calculated ignoring the incompressibility condition given by Eq. (1) and the pressure term (considering ${\displaystyle p=0}$ in Eqs. (2) and (3)). In this step, a Reaction-Diffusion equation is solved by temporal discretization. In the second step, pressure is calculated by solving a Poisson equation. The most attractive property in the Projection methods is the fact that in each step a sequence of decoupled elliptical equations are solved for pressure and velocities, enabling large scale numerical simulations to be executed efficiently; however, it is not trivial to develop and analyze high order Projection methods [14].

There are several formulations for the Projection methods, differing in the way to calculate the auxiliary velocity ${\displaystyle {\bf {u}}_{t}}$ and the projection step. Each of these comes from the method of [24,25], which is a Projection method with a simpler correction scheme for the pressure, which is not yet a first order method (regarding time) for pressure on norm ${\displaystyle l_{2}}$ [14] due to the absence of the pressure gradient on the first step of the method. Goda observed that by adding the pressure gradient [26], Chorin's algorithm presents an improvement in accuracy, an idea applied by Van Kan to develop a incremental pressure-correction scheme of second order [27], which was improved by Timmemans et al. through a modification on the boundary conditions for pressure from the method mentioned previously [28]. This algorithm is known as “incremental correction for pressure on rotational form” and it will be adopted in this paper for the solution of the Navier-Stokes equations [29] . Such algorithm will be given below and it will be designated as Algorithm 1.

Furthermore, Guermond and Shen show that the incremental version on rotational form is a second order scheme for velocity on ${\displaystyle l_{2}}$-, ${\displaystyle l_{1}}$- and ${\displaystyle l_{\infty }}$-norms [29]. For pressure the scheme is a second order scheme for norms ${\displaystyle l_{2}}$- and ${\displaystyle l_{1}}$-norms and 3/2 order for ${\displaystyle l_{\infty }}$-norms. The method consists on the following steps:

First step:

where 'ut is the auxiliary velocities field, un is the velocities field on the time step n, ht is the temporal refinement, Pn is pressure on the time step n and bn the boundary conditions on time step n, ${\displaystyle \beta g=u\cdot \nabla u}$ with constant ${\displaystyle \beta }$ and ${\displaystyle \nu =1{\mbox{/Re}}}$.

Rearranging the terms on Eq. (21) , there

which has the form of a Reaction-Diffusion equation (Eq. (6)) for the variable ${\displaystyle {\bf {u}}_{t}}$, with:

and source term

Second step:

where ${\displaystyle {\phi }^{n+1}}$is the pressure correction.

If the time step ${\displaystyle h_{t}}$ is too large in relation to spatial refinement (${\displaystyle h_{x}}$ and ${\displaystyle h_{y}}$), the simulations can become unstable [30]. The criterion applied in this paper to avoid instabilities is given by [12,30]:

Since this paper considers ${\displaystyle h_{x}=h_{y}=h}$, the criterion can be rewritten as

In [30] other criteria are presented, which support the numerical stability. A common characteristic of the Projection methods is that they impose Neumann boundary conditions for the pressure on the contours [14]. Due to the boundary conditions, the following integrability condition has to be satisfied in order to assure the existence and uniqueness of the solution [5]:

For details on how Eq. (28) is satisfied for each time step, see the procedure proposed by [31].

4. Method for linear systems solution

The discretization for the partial differential equations to be solved results on the following algebraic equations system:

where b and ${\displaystyle f\in R^{n}}$ where b can be u, v or ${\displaystyle \phi }$ , according to steps 1 and 2 on Algorithm 1.

4.1 Gauss-Seidel method (GS)

Matrix A being divided in the form:

with non-singular P, and P is called preconditioning matrix or preconditioner.

To solve Eq. (29), there is the following method, which is called basic iterative method:

where the superscript of b is the representation of the iterations.

Considering the method on the following mode:

where ${\displaystyle B=P^{-1}N}$, ${\displaystyle N=P-A}$, and B is called the iteration matrix of the method given by Eq. (32).

If the diagonal entries of A are non-zeros, the unknown corresponding value can be isolated in each equation, resulting in the Gauss-Seidel method, given an initial estimate ${\displaystyle b^{(0)}}$:

4.2 Preconditioned Conjugated Gradient (PCG)

Descent methods are based on obtaining a solution b for a system of the type given by Eq. (29) by the minimization property of an appropriate norm for the error. From an initial estimate ${\displaystyle b^{(0)}}$, these methods determine a solution for the problem by diminishing the error.

The vectorial norm being given by ${\displaystyle {\Vert b\Vert }_{2}=\langle b,Sb\rangle }$, where ${\displaystyle S\in R^{nxn}}$ is a positive defined symmetrical matrix. Considering the residual equation ${\displaystyle Ae^{(k)}=r^{(k)}}$, it can be deduced that:

where ${\displaystyle R=A^{-T}SA^{-1}}$, which means that minimizing the error equals minimizing the residue in the same norm [7].

On the Descent method one starts from an initial estimate ${\displaystyle b^{(0)}}$ for the solution and a succession of vectors ${\displaystyle b^{1},\cdots ,b^{(k)},\cdots }$ is generated, moving successively towards diminishing the error, making this succession converge to b.

One of the ways to choose the best direction in each point ${\displaystyle b^{(k)}}$ would be to choose the direction of ${\displaystyle p^{(k)}}$ which is lined up with the direction of the function gradient, but in the opposite orientation (Gradients method).

Thus, the iterative expression of the Gradients methods is given by [7]:

where ${\displaystyle G_{k}=I-{\alpha }_{k}S}$ and ${\displaystyle c_{k}={\alpha }_{k}SA^{-1}f}$, ${\displaystyle {\alpha }_{k}}$ being a constant.

Conjugate Gradient method (CG) is an iterative descent method applied in cases where the iteration matrix A is a symmetric positive definite matrix. This method demands a more careful choice of the descent directions ${\displaystyle p^{(k)}}$.

In order to exemplify it, consider a two-dimensional problem. The lines ${\displaystyle {\Vert e^{\left(k\right)}\Vert }_{2}={\Vert r^{\left(k\right)}\Vert }_{2}=}$${\displaystyle {\mbox{constante}}}$, are concentric ellipses with the center in b, as illustrated in Figure 3.

Figure 3. Graphical representation of the Conjugated Gradient method

Supposing there is an initial estimate ${\displaystyle b^{(0)}}$, according to the Gradient method, one can start from this point and research the minimum of ${\displaystyle {\Vert e^{\left(k\right)}\Vert }_{2}}$, across the direction ${\displaystyle r^{(0)}}$, obtaining ${\displaystyle b^{(1)}}$. Observing Figure 3, it can be verified that the best progression of ${\displaystyle b^{(1)}}$ would not be the more inclined direction, but the direction ${\displaystyle p^{(1)}}$ which points to the center of the ellipse.

In fact, the new minimizing point ${\displaystyle b^{(2)}}$ would coincide with the exact solution of b, and thus the iterative process would end with just two iterations. Therefore a direction ${\displaystyle p^{(1)}}$ can be chosen as:

As it can be observed, complying with the orthogonality of the successive residues verified previously, there is

in such a way that

The sought directions are

where ${\displaystyle \omega }$ is given by

The iterations of the CG method satisfy the following error estimate [7]:

where ${\displaystyle cond_{2}A={\Vert A\Vert }_{2}\cdot {\Vert A^{-1}\Vert }_{2}}$.

Thus, it can be noticed that the CG method, theoretically, calculates the exact solution in n iterations, which is not verified in practice due to the application of this method in finite precision arithmetic’s and the introduction of rounding errors which destroy the exact conjugacy of the direction ${\displaystyle p^{(k)}}$ and the exact orthogonality of the residues ${\displaystyle r^{(k)}}$. Thus, this method takes a typical iterative form.

One way to make this iterative method faster consists in transforming the system ${\displaystyle Ab=f}$ into an equivalent system, where the matrix has a more favorable conditioning number. According to Eq. (41), if it is possible to reduce ${\displaystyle cond_{2}A}$, consequently a fast convergence is achieved.

The transformation of the pre- and post-multiplication of the matrix A by the matrixes P and Q being invertible by the following manner:

The system now corresponds, therefore, to ${\displaystyle {\overline {A}}{\overline {b}}={\overline {f}}}$, with ${\displaystyle {\overline {A}}=PAQ}$, ${\displaystyle {\overline {b}}=Q^{-1}b}$ and ${\displaystyle {\overline {f}}=Pf}$, the matrixes P and Q having to be chosen in the way that provides the best convergence properties to the matrix in relation to the original matrix A. This procedure is known as preconditioning of the original equations system, P and Q being known by the name of preconditioning matrixes.

Considering that in the CG method the transformed matrix has to keep being a symmetric positive definite matrix, the preconditioning will also have to preserve this property. In fact, it can be obtained making ${\displaystyle P=S}$ and ${\displaystyle Q=P^{T}=S^{T}}$, and in this manner, it can be concluded that ${\displaystyle {\overline {A}}{\overline {b}}={\overline {f}}}$, ${\displaystyle {\overline {A}}=SAS^{T}}$, ${\displaystyle {\overline {b}}=S^{-T}b}$ and ${\displaystyle {\overline {f}}=Sf}$.

Thus, it can be verified that ${\displaystyle {\overline {p}}^{\left(0\right)}={\overline {r}}^{\left(0\right)}=}$${\displaystyle {\overline {f}}-{\overline {A}}{\overline {b}}^{\left(0\right)}=}$${\displaystyle Sr^{\left(0\right)}.}$ On the other hand, choosing ${\displaystyle p^{\left(0\right)}=S^{T}{\overline {p}}^{\left(0\right)}=}$${\displaystyle S^{T}{\overline {r}}^{\left(0\right)}=S^{T}Sr^{\left(0\right)},}$ there is ${\displaystyle \langle {\overline {p}}^{\left(0\right)},{\overline {A}}{\overline {p}}^{\left(0\right)}\rangle =}$${\displaystyle \langle p^{\left(0\right)},Ap^{\left(0\right)}\rangle ,}$ the relation ${\displaystyle \langle {\overline {r}}^{\left(0\right)},{\overline {r}}^{\left(0\right)}\rangle =}$${\displaystyle \langle {\overline {r}}^{\left(0\right)},r^{\left(0\right)}\rangle }$ also being valid, with

and

From that, there is:

and

keeping these expressions for the following iterations.

Algorithm 2 presents the preconditioned CG method (PCG), where it can be noticed that at each iteration it is necessary to solve the auxiliary system:

and it is desirable that its resolution is not very laborious.

4.3 Incomplete LU decomposition (ILU)

The incomplete LU decomposition (ILU) consists in decomposing a matrix A in an incomplete manner, as the name suggests. Such decomposition is in the A = LU – R form, where L and U are, respectively, inferior and superior triangular matrixes A and R represents the residue or decomposition error. On the matrix R, conditions are imposed such as having null terms in determined positions in order to maintain certain characteristics of the original matrix A.

A general algorithm for ILU can be obtained applying the Gaussian Elimination [7] and having some elements in predetermined positions. For this, a set of points is determined ${\displaystyle P\subset \left\{\left(i,j\right)\vert i\not =j{\mbox{, com }}1\leq i,j\leq n\right\}}$ where the elements of the matrix will be null. In this case, there can be several variations of the ILU decomposition, for example: ILU(0), where the structure of the LU decomposition is identical to the structure of matrix A, that is, it does not allow extra fillings; ILU(1), where the filling of one additional diagonal is allowed in each LU factor, etc. In this paper, ILU(1) was chosen, which will be called 7-points ILU (since the stencil of the original problem has 5 points).

Thus, given the incomplete LU decomposition, it is solved in an iterative process composed of the following steps:

(1) ${\displaystyle Lz^{\left(k\right)}=r^{\left(k\right)}}$ to obtain ${\displaystyle z^{\left(k\right)}}$ , ${\displaystyle r^{\left(k\right)}}$ being the residue vector

(2) ${\displaystyle U{\delta }^{\left(k\right)}=z^{\left(k\right)}}$ to obtain ${\displaystyle {\delta }^{\left(k\right)}}$ .

Thus, the solution in the iteration k +1 is such that:

4.4 Multigrid method (MG)

The MG method [5,6,32] makes use of the error smoothing characteristics from the classic iterative methods: it occurs fast (in the initial iterations) for oscillatory components, while for smooth components, for a high number of iterations, such methods lose their efficiency. Thus, the MG method works with a system of auxiliary coarser grids (with a lower number of points) in which the error components are quickly smoothed, in order to return to the original grid. The information is transferred between grids through operators, called restriction operators (information from a fine to a following coarser grid), generically represented by ${\displaystyle {\left[I\right]}_{{\mbox{ }}h}^{{\mbox{ }}H}}$ and defined by

or prolongation ones (information from the coarse grid to a finer one), represented generically by ${\displaystyle {\left[I\right]}_{{\mbox{ }}H}^{{\mbox{ }}h}}$ and defined by

The restriction operator applied in this paper was the complete weighting operator and the prolongation operator was the bilinear interpolation [5,6,32,33].

In this paper, moreover, the V-Cycle was applied, for its computational efficiency. The equations of Poisson and Reaction-Diffusion being linear, the correction scheme (CS) was applied, which transfers only the residues to the coarser grids [5,6]. The coarsening ratio, considering uniform grids, is defined as r = H/h, where h represents the size of the spacing of the fine grid, ${\displaystyle {\Omega }^{h}}$ , (also called the dimension of the fine grid elements), and H the size at the element of the next coarser grid, ${\displaystyle {\Omega }^{H}}$ .

Figure 4 illustrates the cycle-V with CS for five grids. The superscripts h, 2h , 4h, ... indicate the grid where the vectors and matrixes are defined, the smoothing numbers (S) being achieved in the restriction process (R) (pre-smoothing) and prolongation (P) (post-smoothing) being represented, respectively, by ${\displaystyle {\nu }_{1}}$ and ${\displaystyle {\nu }_{2}}$.

Figure 4. Cycle-V for five grids

In this paper, the MG with GS and ILU solvers was applied. Such iterative methods were also applied as preconditioners for the PCG method, that is, they were used to solve the linear system given by Eq. (50) at each iteration of the CG method.

5. Results

5.1 Methodology and computational details

Computational tests were done to calculate the values of p, u and v in the discretized model problem, Eqs. (22) and (25), applying the following methodology:

1) Multigrid with Gauss-Seidel solver (MG-GS).

2) Multigrid with ILU solver (MG-ILU).

3) Preconditioned Conjugated Gradient with Multigrid and Gauss-Seidel solver (PCG-MG-GS).

4) Preconditioned Conjugated Gradient with Multigrid and ILU solver (PCG-MG-ILU).

In all simulations and for all the variables of interest, the stop criterion applied was the nondimensionalized residue ${\displaystyle l_{2}}$-norm based on the initial estimate, given by:

where ${\displaystyle r^{\left(k\right)}}$ is the residue on the current iteration, ${\displaystyle r^{\left(0\right)}}$ is the residue of the initial estimate, and ε is the tolerance.

In this paper the null vector was utilized as the initial estimate for all the variables of interest and with the use of double precision. The standard coarsening ratio was utilized, that is , ${\displaystyle r=2}$ [5], where the Multigrid method started from the finest grid and got to the coarsest grid possible, that is, the grid with ${\displaystyle N=2\times 2}$ volumes. The number of unknowns in the finest grids was: ${\displaystyle N=16\times }$16, 32${\displaystyle \times }$32, 64${\displaystyle \times }$64, 128${\displaystyle \times }$128, 256${\displaystyle \times }$256, 512${\displaystyle \times }$512 and 1024${\displaystyle \times }$1024. The number of pre- and post-smoothing was equals to ${\displaystyle \nu ={\nu }_{1}={\nu }_{2}=3}$.

The computer utilized for the simulations has an Intel Core i7-5500 processor, a 2.40 Ghz CPU, 16.0 Gb memory and a 64 bits Windows 10 operational system.

Four parameters were analyzed in order to compare the four methodologies:

• ${\displaystyle l_{\infty }}$-norm for numerical error (${\displaystyle ||E||_{\infty }}$);

• ${\displaystyle l_{2}}$-norm for nondimensionalized residue based on the initial estimate (Eq. (54));

• Empirical convergence factor, defined as [6]:

and

• CPU time.

5.2 Analysis of the pressure

This phase it is intended to determine the best method to solve the Poisson equation, which models pressure, determining a method to solve the velocities equations. In this case, MG-GS, MG-ILU, PCG-MG-GS and PCG-MG-ILU were utilized for pressure and MG-GS for the velocities.

The behavior of the ${\displaystyle l_{\infty }}$-norm was initially evaluated for the four studied methods, where it could be noticed that all of them presented a reduction in the numerical error norm with the increase in the size of the problem; a desirable characteristic in a convergent numerical method (Figure 5).

Figure 5. ${\displaystyle l_{\infty }}$-norm for the numerical error (${\displaystyle ||E||_{\infty }}$) versus size of problem (N) for p

Figure 6 presents the number of iterations necessary in order to reach the stop criterion for the problem with ${\displaystyle N'=1024\times 1024}$. Based on this figure, it can be noticed that PCG-MG-ILU reached the stop criterion in a lower number of iterations.

Figure 6. ${\displaystyle l_{2}}$-norm for the nondimensionalized residue by the norm of initial estimate versus number of iterations for the problem with ${\displaystyle N'=1024\times 1024}$ for p

Figure 7 illustrates the behavior of the empirical convergence factor with the increase of the size of the problem. Based on this figure, it can be claimed again that the PCG-MG-ILU methodology presented the best results, since it is the one that generated the best values for the convergence factors, that is, values near to zero [5,34].

Figure 7. Empirical convergence factor (${\displaystyle q^{(k)}}$) versus size of the problem (N) for p

Figure 8 illustrates the behavior of the CPU time (in seconds) of the simulations with the increase of the size of the problem. It can be noticed that the computational times became very close, with a slight advantage for the PCG-MG-ILU methodology, which can be confirmed in sequence.

Figure 8. CPU time (in seconds) versus size of the problem (N) for p

The speed-up is a measure utilized to determine the increase in velocity obtained during the execution of a program utilizing an algorithm “A” in relation to its execution utilizing an algorithm “B” [6]. It is given by:

Table 1 shows the Speed-up of the three methods MG-GS, MG-ILU and PCG-MG-GS in relation to the PCG-MG-ILU method.

According to the data presented in this table, it can be confirmed that the CPU times are close among themselves, but there is a slight advantage with the use of PCG-MG-ILU (Speed-up above 1), as described in [35].

It all leads us to conclude that the CGP-MG-ILU method has the potential to accelerate the convergence of the presented problem.

5.3 Analysis of the velocities

In this phase, it is intended to determine the best method to solve the Reaction-Diffusion equations, which model the velocities u and v, determining the best method for pressure presented in the previous section, which was PCG-MG-ILU. Since the results for the velocities u and v for all the parameters are similar, it was chose to display only the results for u.

Figure 9 shows the behavior of the numerical error versus the size of the problem, for the several solution methods. It can be noticed that the four studied methods presented an expected characteristic for the convergent iterative methods, that is, the decrease in the numerical error with the increase in the size of the problem. It can also be noticed that the presented error have values which are very close among themselves.

Figure 9. ${\displaystyle l_{\infty }}$-norm for the numerical error (${\displaystyle ||E||_{\infty }}$) versus size of the problem (N) for p

Figure 10 displays the behavior of ${\displaystyle _{2}}$-norm of the nondimensionalized residue by ${\displaystyle l_{2}}$-norm of the initial estimate in function of the number of iterations for the problem with ${\displaystyle N=1024\times 1024}$ volumes. Based on this figure, it can be noticed that the PCG-MG-ILU and the MG-ILU methods were the ones which achieved the stop criterion in the lowest number of iterations.

Figure 10. ${\displaystyle l_{2}}$-norm for the nondimensionalized residue by ${\displaystyle l_{2}}$-norm of the initial estimate versus number of iterations for the problem with ${\displaystyle N=1024\times 1024}$ for u

Figure 11 illustrates the behavior of the empirical convergence factor in function of the size of the problem. Based on this figure, it can be concluded that the PCG-MG-ILU method presented the best results, since it generated the best values for the convergence factors (closest values to zero), besides decreasing its value with the increase of the size of the problem.

Figure 11. Empirical convergence factor (${\displaystyle q^{(k)}}$) versus size of the problem (N) for u

Figure 12 shows the CPU times (in seconds) of the simulations with the increase of the size of the problem. Based on this figure, it can be noticed that the computational times were very close, with a slight advantage for the PCG-MG-ILU methodology, which can be confirmed in the results presented in the sequence.

Figure 12. CPU time (in seconds) versus size of the problem (N)

Table 2 presents the Speed-up (Eq.(56)) of the MG-GS, MG-ILU and PCG-MG-GS methods in relation to the PCG-MG-ILU method. In this table, it can be confirmed that the CPU times are close among themselves, but there is an advantage with the use of PCG-MG-ILU, that is, the method that most accelerated the convergence of the iterative process for the presented problem.

5.4 Algorithm

After the tests of the two previous phases (sections 5.2 and 5.3), whose objective was to identify the best methods for the solution of the equations for the variables p (Poisson equation) and u and v (Reaction-Diffusion equations), we get to the point where the best methods were applied simultaneously to solve the Navier-Stokes equations.

The results presented below are compared the classic Gauss-Seidel (GS) method in its singlegrid version (only grid method) so that the gains can be measured when using the proposed algorithm (adaptation of the Projection method presented in Algorithm 1), that is, using the PCG-MG-ILU solver in the solution of the linear systems given in steps 1 and 2 of Algorithm 1.

Figure 13 presents the behavior of ${\displaystyle l_{2}}$-norm of the nondimensionalized residue by the initial estimate in function of the number of iterations for the problem with ${\displaystyle N=128\times 128}$ volumes for the variables u, v and p. Based on this figure, it can be noticed that the PCG-MG-ILU method reached the stop criterion ε = 10^–10 in an extremely lower number of iterations, both for p and for u and v.

Figure 13. ${\displaystyle l_{2}}$-norm of the nondimensionalized residue by norm ${\displaystyle l_{2}}$ for the initial estimate versus number of iterations for the problem with ${\displaystyle N=128\times 128}$ for u, v and p

Figure 14 illustrates the behavior of the empirical convergence factor with the increase of the size of the problem. Based on this figure, it can be concluded that the PCG-MG-ILU method presented much better results for the convergence factors for the three variables of interest, comparing with the classic GS method.

Figure 14. Empirical convergence factor (${\displaystyle q^{(k)}}$) versus size of the problem (N) for u, v and p

Figure 15 contains the values of the CPU times (in seconds) of the simulations with the increase of the size of the problem. Based on this figure, it can be noticed that the PCG-MG-ILU methodology was the one that converged faster, being quite evident the difference in CPU time between the proposed methodology and the classic GS method.

Figure 15. CPU time (in seconds) versus size of the problem (N)

Table 3 presents the Speed-up (Eq.(56)) on the GS method in relation to the PCG-MG-ILU method. According to the data on this table, the methodology that uses the PCG-MG-ILU method proved much more efficient than the methodology that utilizes the GS method, reaching a 105 times higher velocity for the biggest problem compared in this table. It can also be noticed that the Speed-up increases as the size of the problem increases, which is a highly desirable property.

Another complementary study can be done through the comparative of complexities of algorithms utilizing the different methods. For complexity analysis of such algorithms, the exponent p was calculated; obtaining the by least square fitting, for the function given by