CC BY
9
DOI: 10.17516/1997-1397-2021-14-4-519-527 УДК 517.9
The Splitting Algorithm in Finite Volume Method for Numericai Solving of Navier-Stokes Equations of Viscous Incompressible Fluids
Viktor M. Kovenya* Daniel Tarraf
Novosibirsk State University Novosibirsk, Russian Federation
Received 10.03.2021, received in revised form 05.04.2021, accepted 20.05.2021 Abstract. For the numerical solution of the Navier-Stokes equations, written in an integral form, an implicit of finite-volume algorithm is proposed, which is a generalization of previously proposed differences schemes. Using the integral form of equations allowed to ensure its conservatism, and the technology of splitting — the economy of the algorithm. The numerical test of the algorithm on the exact solution, in problems about the viscosity flow in the cavern with a moving lid and the current of the heated walls of the channel, confirmed the sufficient accuracy of the algorithm and its effectiveness. The work is presented in the issue of the memory of Prof. Yu. Ya. Belov.
Keywords: Navier-Stokes equations, viscous flows, finite-volume method, splitting algorithms. Citation: V.M.Kovenya, D.Tarraf, The Splitting Algorithm in Finite Volume Method for Numericai Solving of Navier-Stokes Equations of Viscous Incompressible Fluids, J. Sib. Fed. Univ. Math. Phys., 2021, 14(4), 519-527. DOI: 10.17516/1997-1397-2021-14-4-519-527.
Introduction
The Navier-Stokes equations of a viscous incompressible fluid are the basic model for solving various classes of problems in hydromechanics [1,2]. They are nonlinear and their solutions contain areas of high gradients, boundary layers, separation zones, etc., which imposes additional difficulties in their study. Therefore, the problem of constructing economical numerical algorithms for solving the Navier-Stokes equations is still relevant today. Some approaches for constructing finite-difference and finite-volumes schemes are given, for example, in [3-11]. When solving multidimensional problems, including those in curvilinear and multiply connected domains, the method of finite volumes [2, 8, 9], based on the approximation of equations in integral form, may turn out to be more convenient. It has the property of conservatism, the approximation of the equations in it, is constructed for each cell, the shape of which is easier to adapt to the boundaries of the region.
The use of explicit schemes in solving the Navier-Stokes equations leads to large expenditures of computer resources, especially in the multidimensional case, due to strict restrictions on the ratio of the temporal and spatial steps of the grid. Implicit unfactorized algorithms are also uneconomical due to the need to invert large matrices (see, for example, [5,11]). An alternative to this approach is the splitting and factorization methods [3], which make it possible to reduce the solution of a multidimensional problem to the solution of its one-dimensional analogs or simpler problems. In [11], a difference scheme was proposed for solving the Navier-Stokes equations for a viscous incompressible fluid in physical variables "velocity, pressure", based on the method
* [email protected] © Siberian Federal University. All rights reserved
of splitting into physical processes and spatial directions. This made it possible to simplify the implementation of the algorithm.
Below we generalize it to the finite volume method. The properties of the algorithm in terms of the accuracy of calculations and the rate of convergence are investigated when finding a stationary solution by the established method. The algorithm was tested on the solution of the Poisson flow problem, which has an exact solution, on the problems of fluid flow in a square cavity with a moving cover and flow in a rectangle with heated walls. The results obtained illustrate the capabilities of the proposed method and allow us to conclude about its effectiveness.
1. Initial equations. Algorithm for solving the Navier-Stokes equations
When studying the flows of a viscous incompressible fluid, taking into account the effects of heat conduction, one usually uses models described by the Navier-Stokes equations of an incompressible fluid, supplemented by the heat conduction equation. Let us consider them in gas-dynamic variables "velocity-pressure" in the form of a system of integral conservation laws in gas-dynamic variables [1]:
M d
d_ dt
UdV + / (Wn)ds = i FdV, 1/ j s jv
i TdV + <£ (WTn)ds = 0, jv j S
(1)
+ f (WTn)ds = 0, (2)
>v J S
where V the volume of the computational domain, S its boundary, n the external normal to the boundary area, U, T the vector of the sought functions and temperature, W, WT are matrices composed of columns of flows at the boundaries of the volume, n, k the coefficients of viscosity and thermal conductivity, g the acceleration coefficient, and F the force of gravity. The algorithm will be presented using the example of two-dimensional equations in Cartesian coordinates written in dimensionless form in the absence of external forces. Then:
p vi v2 0 0 0 0
U= vi , W = 2 i vf + p — ai v v2 — a2 , F = 0 , M = 0 1 0
V°2_ _ v i vf — a f vf + p — a|_ d 0 0 1
(3)
dvi
dT
= , aj = —, WT = [viT — af v2T — af] , d = agT, j = const.
dxj'
Let's set the grid step t = T/N, where N is the number of time steps. U, T the functions will be set at the nodes i,j of the cell, and the flows W at the boundaries of the cells at the nodes i ± 1/2, j and i,j ± 1/2 (Fig. 1).
We introduce the averaging of the sought functions over the elementary volume
Vi,j = u,
U i. = 1
u
U du, Ti. = 1 u
Tdu,
F. . = 1
i j
u
F du,
and we approximate the integral operators in cells by grid operators by the formulas
d f Un+i — Un
— U du « u-
dt. v t
, <f (WSn)dS « Q = ^ Am(WmS),
JS m=i
S(Wtn)ds « ^ Am (WTmS)
a
j
Fig. 1.
where
SW
v\S\ + V2S2 (vl + p - <7i)Si + (V1V2 - al)s2 (v 1V2 - a 2)Si + (v2 + P - a|)S2
Am(WS)m = [Sm+1 /2(Wm+1 + Wm) - Sm-1 /2(Wm + Wm-1 )]/2
is the flows through opposite faces of a cell, m = 1 corresponding to index i, m = 2 index j. We construct an algorithm for solving the system of equations (1), (2), at first for equations (1), assuming that the temperature value is known. Finite-dimensional scheme with weights
M-
U
n+1
+ Un
1
+ -
Am(aW 1 + I3W nm)Sn
m= 1
(4)
approximates the original equations (1), (3) with order O(r2 + h?) for a > 0.5 and at a = 0, it is nonlinear. Here h = w1/2. Operators alm on the boundaries contain directional derivatives with respect to xm, that cannot coincide with the direction of the cell faces, so we introduce parameterization Xi = Xi(qj), qj = qj(xi) where 0 ^ qj ^ 1. Then:
d
d
— = Y^ zl — zl = — jj
dx
1=1
dqi''
dxn
n>zl
dvi
^ j dqi
l=1 l
(5)
and in new variables aj contain derivatives with respect to normals q1 and tangents q2 to the cell boundaries w. We approximate them at the nodes i ± 1/2 or j ± 1/2 by symmetric operators
aj = M E
zm Amvi, Am±1/2v = ±(vm±1j - vm) and linearize the vector Wn+1 with respect to
m=1
U and known values m, z,
l: m'
W
n+1
Wn + T* r m. 1 '
dWm dUn ~3U dt
+ O(T 2) = W m + TBr
U
n+1
Un
+...,
where Bm
0 S1 S2
51 V - tmk 0
52 0 V - tmk
tmk = zmSkAk, V = v1S1 + v2S2 the projection
k=1
of the velocity vector V times the normal to the face area. To construct economical algorithms, we introduce an operator Bm in which only derivatives with respect to the normals are stored in
a
j
m
the coefficients tmk, and represent it in the form of a splitting into physical processes:
0 0 0 0 Si Sf'
Bm = Bim + B2m, Bim = 0 V — tm 0 , Bfm = Si 0 0
0 0 V — tm Sf 0 0
t = nZmS A
the finite-volume scheme:
— Un+i — Un C-
(6)
- - nn + Fn, C = M + — ]T Am(Bnm + Bnm), № №
m=1
2
Qn =J2Am (WnSn)
m=1
linear, but approximates the original equations with order 0(t + h?). Note that the matrix M is degenerated, which does not allow standard factorization methods. However, there is a special splitting of the operator C in the physical process and spatial directions (see [11]), in which it can be factorized in the form:
C = C + O(T), C = TT (I + dAmBnm) • (M + d V AmBnm), d = —
u
Then the scheme of approximate factorization:
U n+i — U 1
C-=--nn + Fn
t u
(7)
or it equivalent scheme in fractional steps:
e = — - nn + Fn, (I + dAiB^i)Cn+i/3 = c, (I + d,A2B^2)Cn+2/3 = C+i/3, (8)
(M + d J2 Am^nm)^n+1 = £n+2/3, un+1 = Un + tC+1
m= 1
Where £ = (£p, £1,£2)T the residuals of the solution at fractional steps approximate Eqs. (1), (2) with the same order as (6). The values £n are computed explicitly. At the first (m = 1) and second (m = 2) fractional steps of the equation scheme
[1 + dAm(V - tm)]$+m/3 = $+(m-1)/3(l = 1, 2)
are solved by scalar sweeps independently for each component of the velocity residual £n+m/3 and £p+m/3 = £pn .At the third fractional step of scheme (8), the system of equations is solvec
d[A1S1tn+1 + A2S2&+1] = $+2/3, en+1 = tn+2/3 - dA1S1^n+1,
(9)
£p , S1 = "<"1"1sp
C2+1 = n2/3 - dA2S2£np+1.
Eliminating the velocity components from the continuity equation, we obtain the equation for it
AC+i = f
where A = An + A22, Amm = AmSmAmSm, f
m=1
/d - C/d2. The solution to
the Poisson equation can be obtained by various iterative algorithms [12], for example, by the iterative approximate factorization scheme
(I - toAii)(I - T0A22 W+1 - ?)/to = A(^+1/3y - f
or an equivalent scheme in fractional steps
(I + T0Au)ni =A(C+1/3y - f, (I - T0A22)n = ni, C+1 = C +
Ton.
Realized with fractional steps, also by scalar sweeps. The solution is carried out until the iterations converge, i.e., until the condition A£v - f = O(T0h?) is met in all cells. Then, from (9), the new values of the velocity residuals are clearly found. The new values of the functions U are explicitly calculated from (8) and, if necessary, the calculation process is repeated.
2. Algorithm for solving the heat equation
We approximate the integral operators in (2) by the grid operators
d r Tn+1 _ Tn
— Tdw « w-,
dt J v t '
¿(Wtn)ds « ^ Am[(v1 - kaf)S1 + (V2T - ka3)S2]T
m=1
and, like scheme (7), we consider the finite-volume scheme of approximate factorization
2 Tn+1 _ Tn 1 2
T[[I + dAm(V - tm)]- = - -J^AmWn (11)
J-J- T W ^
m=1 k=1
or the equivalent of a scheme in fractional steps
en = --£m=1 Amwn, [i+dA1(v - t1)]en+1/2 = en, , s
w (12)
[I + dA2(V - t2)]^n+1 = ^n+1/2, Tn+1 = Tn + T^n+1.
2
Here tm = k(zYiS1+zJ^lS2)Am, a'im = k J2 zkmAkT, = [VT-(a:^S1+a^S2)]. It approximates
k=1
the thermal conductivity equation (2), (5) with order 0(t2+h2) when a = 0.5+0(t) . The values ^n are calculated explicitly, then the equations are solved at fractional steps by scalar sweeps in each direction. The new temperature values are calculated explicitly from the last equation in scheme (12). This completes one step of the calculations and, if necessary, the process continues to find a solution at subsequent points in time.
3. Examples of numerical calculations
The proposed algorithm was tested on a number of problems. Testing was carried out on three problems, the numerical solution of which was obtained by different authors and different numerical algorithms (see, for example, [5,7-9]). This made it possible to compare and evaluate the properties of the algorithm. When carrying out numerical calculations, equations (2), (3)
were reduced to dimensionless form [1]. This led to the appearance in the equations of the di-mensionless parameters of the Reynolds numbers Re, Rayleigh Ra and Prandtl Pr, respectively: Re = 1Ra = ag, Pr = ¡i/k. Then the characteristic size of length L « 1, speed |v| = 1, and pressure p = 1 (it is set up to a constant).
In the first problem, the stationary flow of a viscous incompressible fluid in a channel was investigated in the framework of the Navier Stokes model (1) where F = 0. Its solution is reduced to the Poiseuille flow, the exact solution of which is: v2 = 0, v1 = 1 — x\, p = 1 — 2^x1.
In the computational domain, a square grid with a number of cells J = I x I was used. The velocity and pressure v2 = 0, v1 = 1 — x\, p = 1 were set at the channel inlet, and v1 = v2 = 0 for the adhesion conditions on the channel walls. At the initial moment of time, constant values of velocity and pressure were set inside the region. The stationary solution of the problem was found by the establishment method. The establishment criterion was set in the form:
max lpn+1 — pnl < K(rw), K « 0.1 — 1.0
for all interior points. Since the sought functions are specified at the centers of the cells, and the flows are determined at the boundaries of the cells, the implementation of the algorithm requires the introduction of near-boundary dummy cells and the specification of functions in them. On the upper and lower walls of the channel at dummy points, the velocity components are set equal in magnitude, but opposite in sign, and at the input, their values are set equal to those at the input. To evaluate the accuracy of the algorithm and to estimate the rate of convergence of the solution to the stationary one, we performed calculations on grids with different numbers of cells. Tab. 1 shows estimates of the errors of solutions
Table 1
T hi = h2 Ap Avi Av2
0.1 0.1 0.008817 0.009947 0.006017
0.05 0.05 0.002204 0.002487 0.001505
0.025 0.025 0.000551 0.000622 0.000376
0.0125 0.0125 0.000138 0.000155 0.000094
As follows from the calculation results, an increase in the number of cells (decrease in grid steps) by 2 times in each direction leads to a decrease in the error by 4 times, which confirms the second order of accuracy of the algorithm. The number of iterations before stopping depends on the initial guess and the number of nodes J. Their typical number is given in Tab. 2 on a J = 80 x 80 grid at various values of the viscosity coefficients
Table 2
0.01 0.025 0.001
iterations 1111 2731 3652
In the second problem, the fluid flow in a square cavity with a moving cover in the absence of gravity was studied. No-slip conditions vi = v2 =0 were set on the stationary walls of the channel. At the initial moment of time v = 0. At t > 0, the lid begins to move at a constant speed (Fig. 2).
The stationary solution of the problem was found by the establishment method. The calculations were carried out on a sequence of grids at various values of the Reynolds number Re.
Fig. 3 shows the distribution of the longitudinal and transverse velocity components on various grids at Re = 3 • 103.
A V|=1
v, =0 v2=0
p = 0 -►
Fig. 2
Fig. 3
The convergence of solutions is observed with an increase in the number of grid steps. Its further refinement practically did not lead to differences in values on the 81 x 81 grid. At numbers Re > 102, a vortex appears in the cavity, the center of which is shifted to the right, and two small vortices at the corners of the cavity. With an increase in Re, the angular vortices increase, their intensity increases, which follows from theoretical estimates and calculations using other algorithms [11,13-15]. A typical flow pattern is shown in Fig. 4.
Comparison of the results obtained with the calculations in [5, 7-9] shows the visual coincidence of the flow fields.
In the third problem, some results of calculations of fluid flows in a closed flat cavity with heating of one of the sides are presented. The system of Navier-Stokes equations (1) was supplemented by the equation for temperature (2), and a term of the form d = RaT was added to the equation of motion. At t = 0, the liquid was assumed to be stationary, and the adhesion conditions were set at the boundaries of the region. On the left and right walls of the region, T =1 and T = 0 respectively, and on the upper and lower walls, according to a linear law T =1 — xi. Due to the temperature difference in the region, a rotational motion of the liquid occurs, its intensity is determined by the numbers Re and Ra. The numerical solution of the problem was found according to schemes (7), (11) on various grids. The calculations for various values of the parameters of the problem and comparison with the calculations [13-15] showed the identity of the solutions obtained. For example, in Fig. 5 streamlines for Ra = 1 Re = 3 • 102(a) and 103(6) are shown.
A large central region of the circulation flow ("central vortex") and secondary "corner vortices" in the lateral part of the cavity are distinguished. Note the increase in the vortex velocity with the increase in the Re numbers. A change in the Pr numbers for fixed Re and Ra also leads to a significant rearrangement of the flow pattern, and a change in the Ra numbers in a wide range of parameters has little effect on convection, as follows from calculations by other authors (see [11, 13-15]).
Fig. 5
Conclusion
The paper proposes a generalization of the finite difference splitting scheme for the numerical solution of the Navier-Stokes equations of a viscous incompressible fluid to the finite volume method. It has been tested for solving a number of problems (Poiseuille flows, in a cavity with a moving cover, and in a square region with heated sides). The performed comparisons in terms of the accuracy of the algorithm and the rate of convergence when finding a stationary solution by the establishmed method showed the efficiency of the algorithm and its sufficient accuracy.
References
[1] L.G.Loytsyansky, Mechanics of liquid and gas, Moscow, Nauka, 1978 (in Russian).
[2] P.J.Roache, Computational Fluid Dynamics, Hermosa Publishers, 1976.
[3] N.N.Yanenko, The method of fractional steps for solving multidimensional problems of mathematical physics, Novosibirsk, Science. Sib. Branch, 1967 (in Russian).
[4] S.K.Godunov, A.V.Zabrodin, M.Ya.Ivanov et al., Numerical solution of multidimensional problems of gas dynamics, Moscow, Nauka, 1976 (in Russian).
[5] A.A.Samarskiy, Difference Scheme Theory, Moscow, Nauka, 1989 (in Russian).
[6] A.I.Tolstykh, Compact difference schemes and their applications to the problems of aerohy-drodynamics, Moscow, Nauka, 1990 (in Russian).
[7] R.J. Le Veque, Finite volume methods for hyperbolic problems, Cambridge University Press, Cambridge, 2002.
[8] P.J.Roache, Computational Fluid Dynamics, Hermosa Publishers, 1976.
[9] R.J. Le Veque, Finite Volume methods for Hyperbolic Problems, Cambridge University Press, Cambridge, 2002
10] J.B.Vos, A.Rizzi, D.Darrac, E.H.Hirschel, Navier-Stokes solvers in European aircraft design, Progress in Aerospace Sciences, 38(2002), no. 8, 601-697.
DOI: 10.1016/S0376-0421(02)00050-7
11] V.M.Kovenya, Algorithms of optimal splitting in computational fluids dynamics, Matem. Mod., 16(2004), no. 6, 23-27 (in Russian).
12] A.A.Samarskiy, E.S.Nikolaev, Methods for solving grid equations, Moscow, Nauka, 1978 (in Russian).
13] U.Ghia, K.N.Ghia, C.T.Shin, High-Re Solutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method, J. Comput. Phys., 48(1982), 387-411.
DOI: 10.1016/0021-9991(82)90058-4
14] S.E.Rogers, D.Kwak, An upwind differencing scheme for the incompressible Navier-Stokes equations, Applied Numerical Mathematics, 8(1991), no. 1, 43-64.
15] M.Hafez, Numerical simulation of incompressible, World Scientific, 2002.
Алгоритм расщепления в методе конечных объемов для численного решения уравнений Навье-Стокса вязкой несжимаемой жидкости
Виктор М. Ковеня Даниэль Тарраф
Новосибирский государственный университет Новосибирск, Российская Федерация
Аннотация. Для численного решения уравнений Навье-Стокса, записанных в интегральной форме, предложен неявный конечно-объемный алгоритм, являющийся обобщением предложенных ранее разносных схем. Использование интегральной формы уравнений позволило обеспечить его консервативность, а технологии расщепления — экономичность алгоритма. Проведена численная апробация алгоритма на точном решении, в задачах о течении жидкости в каверне с движущейся крышкой и течении с подогревом стенок канала, подтвердившая достаточную точность алгоритма и его эффективность. Работа представлена в выпуск памяти профессора Ю. Я. Белова.
Ключевые слова: уравнения Навье-Стокса, вязкие течения, конечно-объемный метод, алгоритмы расщепления.
