Numerical analysis of slow steady and unsteady viscous flow by means of R-functions method Текст научной статьи по специальности «Математика»

МАТЕМАТИЧНЕ ТА КОМП'ЮТЕРНЕ МОДЕЛЮВАННЯ

МАТЕМАТИЧЕСКОЕ И КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ

MATHEMATICAL AND COMPUTER MODELING

UDC 004.93

NUMERICAL ANALYSIS OF SLOW STEADY AND UNSTEADY VISCOUS FLOW BY MEANS OF R-FUNCTIONS METHOD

Artiukh A. V. - PhD, Associate Professor of the Department of Applied Mathematics of Kharkiv National University of Radio Electronics, Kharkiv, Ukraine.

Lamtyugova S. N. - PhD, Associate Professor of the Department of Advanced Mathematics of O. M. Beketov National University of Urban Economy in Kharkiv, Kharkiv, Ukraine.

Sidorov M. V. - PhD, Associate Professor, Associate Professor of the Department of Applied Mathematics of Kharkiv National University of Radio Electronics, Kharkiv, Ukraine.

ABSTRACT

Context. This article is devoted to the linear problem of the steady and unsteady flows of a viscous incompressible fluid.

Objective. The purpose of this paper is to compare the previously developed methods of numerical analysis of the steady and unsteady flows of a viscous incompressible fluid.

Method. The flow of a viscous incompressible fluid can be described by the system of nonlinear Navier-Stokes equations. The variables of this system are velocity, pressure, density, volume forces, and fluid viscosity. Using the stream function, the Navier-Stokes equations can be transformed to the initial-boundary problem with the differential equation of the fourth order. To solve the problem the structural variational method of R-functions and Ritz method (steady problem) or Galerkin method (unsteady problem) are used. The R-functions method allows satisfying the boundary conditions accurately and transforming them to the homogeneous, which are the prerequisite for application of Ritz or Galerkin method. The problem transforms to the solving the system of linear algebraic equations or the system of ordinary differential equations for steady and unsteady flows respectively. The matrices elements are the scalar products in the norms of the corresponding differential operators. Numerical integration was made by means of Gaussian quadratures with 16 points. Solutions of the system of linear algebraic equations and the system of ordinary differential equations were found with the help of the Gauss method and the Runge-Kutta method with an automatic step-size control respectively. The existence of a unique solution of the problems is proved.

Results. The computational experiments for the problem of flows of a viscous incompressible fluid for the different rectangular domains carried out.

Conclusions. The conducted experiments have confirmed that the stream function, the flow velocity, and other flow characteristics are converging to the steady state when the time is increasing. This allows us to say that the obtained methods work as expected. The further research may be devoted to the comparison of the solution methods for the non-linear problems.

KEYWORDS: Navier-Stokes equations, steady flow, unsteady flow, viscous fluid, stream function, R-functions method, successive approximations method, Ritz method, Galerkin method.

ABBREVIATIONS

FDM is the Method of Finite Differences; FEM is the Method of Finite Elements; SIMPLE is the Semi-Implicit Method for Pressure Linked Equation;

SIMPLER is the SIMPLE Revised; RFM is the R-Functions Method.

NOMENCLATURE

x, y, z are variables of the Cartesian coordinate system;

t is a time;

H is a space region, a cylinder; dH is a boundary of H;

Q is a flow domain, a cross-section of H by a plane perpendicular to the axis Oz ;

dQ is a boundary of a streamlined body; Q is a closure of Q ; n is an outer normal to dQ; p is a density; p is a pressure;

v is a field of fluid velocities, v = (vx, vy, vz ) ;

X' "yi "z>

v0 is a field of fluid velocities at the initial instant of time, v0 = (v0, v0, v0);

a is a fluid velocity vector on the boundary 5S, a = (ax, ay, az );

F is a field of volume forces, F = (FX, Fy, Fz ) ;

Re is a Reynolds number;

v = Re-1 is a coefficient of viscosity;

y is a stream function, y = y(x, y, t) ;

y0 is a stream function at the initial instant of time,

yo =yo( x y);

EC is a gluing operator of boundary values; D1 is an operator, D1 = (V, V) ; ® is an undefined component of the solution structure;

ra = 0 is a normalized equation of dQ ; ra is a sufficiently smooth function describing the geometry of the domain Q ; V is a gradient operator;

A is a biharmonic operator; A is a Laplace operator;

2

A is an operator, A = A ;

Da is a domain of definition of an operator A ;

Ha is an energy space of an operator A ;

B is an operator, B = A ;

DB is a domain of definition of an operator B ;

HB is an energy space of an operator B ;

Cn (Q) is a space of continuous functions which are n-times continuously differentiable on Q ; W22 (Q) is a Sovolev space;

L2 (Q) is a Hilbert space of square-integrable on Q functions;

{xk } is a complete in L2 (Q) sequence of functions; {q>k } is a sequence of coordinate functions; N is a number of coordinate functions; ck (t ) is a sequence of unknown functions in Galerkin method;

M0(x0, y0) is a fixed point inside Q ;

M (x, y) is an arbitrary point inside Q ;

B6 ( x) is a Schoenberg spline of the sixth order.

INTRODUCTION

In this paper, we will address the problem of developing methods for the numerical simulation of the steady and unsteady flows of incompressible fluids. Indeed, the mathematical investigation of flows is one of the most complex and demanded problems nowadays. It takes place in chemistry, geophysics, biology, energetic thermal power engineering etc.

Usually, the fluid is considered as viscous Newtonian incompressible medium. This assumption allows us to express these flows by the Navier-Stokes differential equations, the solution complexity of which is very high due to its non-linear part. Therefore, we have an exact solution only for some particular cases, which are not suitable for the simulations of real processes. That is why it is quite important to have the numerical methods to get an approximation function for this problem.

The purpose of this paper is to develop numerical methods of getting an approximate solution of the Navier-Stokes equations by means of the R-functions, Ritz (steady flows) and Galerkin (unsteady flows) methods. The main advantage of these methods is that the solution of a differential equation can be obtained in an analytical form. This benefit allows one to investigate the solution more precisely.

1 PROBLEM STATEMENT

The unsteady creeping flow of a viscous incompressible fluid in the space region H is described by the well-known linearized system of equations (Stokes linearization) in natural «velocity-pressure» variables [1]

+ Fx

dvx 1 dp ----— + V ^ 2 d vx d2 vx +-x ~2 \ d vx +--

dt p ôx v cx2 dy2 dz2 J

dvy 1 dp ---— + V P & 5 vy a2vy 2vy ]

dt dx 2 V ' dy2 ' dz2 /

dvz 1 dp =---— + V 5 vz d2 vz + —Í d2 A d vz + —z

dt p dz V dx2 dy dz y

dvx+dvy+dvz= 0.

dx dy dz

+ Fy

+ Fz

(1)

The system of equations (1) should be supplemented with initial and boundary conditions:

vl t=o = v

Vl55 = a .

(2)

The conditions under which the problem (1), (2) has a unique solution are stated in [2].

Further let us consider plane-parallel flows, i.e. such ones for which the region H is cylindrical along the axis Oz, and the boundary and initial conditions do not depend on z . For such flows, one can introduce the stream function y = y(x, y, t) by the formulas

¥

vy =-

dx

(3)

Due to the introduction of the stream function we pass from the system (1) (suppose that the field F is potentially) to one fourth-order differential equation:

+vA2v = 0 in n.

dt

(4)

vx =

Based on the definition of the velocity vector v on

d S and at t = 0, one can obtain the values , —

dn

dQ

and t=0 [3].

Then for unsteady plane-parallel creeping flow of a viscous incompressible fluid, we can set the following initial-boundary value problem:

8(-Ay) + vA2y = 0, ( x, y) eQ , t > 0, (5)

dt

H=0 =^0(x, y), (x, y) eQ = QU dQ, (6) ^ = g0(s,t), 5 edQ , t > 0. (7)

Han = f0(s,t )

dn

dQ

The corresponding for (5)-(7) problem which describes stationary plane-parallel flows will have the form

h|5Q= f0( 5)

A> = 0, ( x, y) eQ

dn

= §0(5), 5 edQ .

(8) (9)

dQ

2 REVIEW OF THE LITERATURE

At present, the numerical mesh-based methods are widely used in different fields of science. Among them are methods of FDM [4, 5] and FEM [6-8], SIMPLE and SIMPLER [9] etc. The main advantage of these methods is the simplicity of implementation. However, when moving to a new domain, it is necessary to generate again and adjust the calculated grid under the task.

The R-functions method [10] of V. L. Rvachev, the Academician of Ukrainian National Academy of Sciences, allows satisfying the region's geometry precisely.

Lee Chern Earn and others [11] compared SIMPLE and SIMPLER algorithm in terms of their convergence rate, iteration number and computational time. They proved that SIMPLER algorithm is more efficient.

Gupta M. and Kalita J. [12] introduced a new methodology is based on a stream function-velocity formulation of the two-dimensional steady-state Navier-Stokes equations based on FDM. The main advantage of their formulation is that it avoids the difficulties associated with the computation of vorticity values.

Claeyssen J. [13] presented a velocity-pressure algorithm for incompressible viscous flow with a Neumann pressure boundary condition based on FDM. Simulations with the cavity problem were made for several Reynolds numbers.

Bognar G. and Csati Z. focused on the numerical investigation of the Navier-Stokes equation applying a spectral method [14]. The correctness of the MATLAB program was tested on the lid-driven cavity flow problem.

Bettaibi S. and Sedik E. [15] presented a hybrid scheme with multiple relaxation time Lattice Boltzmann Method to obtain the velocity field while the temperature field is deduced from energy balance equation by using

© Artiukh A. V., Lamtyugova S. N., Sidorov M. V., 2019 DOI 10.15588/1607-3274-2019-1-3

FDM. They considered a mixed convection heat transfer in two-dimensional lid driven rectangular cavity. The code has been validated for the classical lid-driven cavity. In addition Bettaibi S. and Sedik E. simulated the effect of Reynolds numbers on mixed convection fluid flow inside a lid-driven square cavity.

3 MATERIALS AND METHODS

Basic concepts of the R-functions method (RFM).

The R-functions were introduced by V.L. Rvachev, the academician of the National Academy of Sciences of Ukraine, in 1963 in his work [10]. Rvachev V.L. was able to notice among the functions of a continuous argument the existence of such functions (R-functions), which, in their properties, resemble the logic algebra functions. These functions form a set that has a non-empty intersection with a set of elementary functions, so that such operations as differentiation, integration, and other ones can be performed above the R-functions. Thus, each R-function was associated with a corresponding Boolean function, which further made it possible to use the developed apparatus for constructing the solution of the inverse problem of analytic geometry.

At the present time, a large number of R-functions systems were constructed [16]. In particular, the system R0 has the form

— I~2 2"

u =-u , u V 0 v = u + v + V u + v

V2 2

u + v . (10)

Such Boolean functions as negation, disjunction and conjunction correspond to the R-functions (10) (are the companion functions for them).

It is known [16], the inverse problem of analytic geometry consists in constructing for a given domain Q such a function ra(x, y) that, having the form of a single analytic expression, is zero on dQ, is strictly greater than zero in Q and negative outside Q. Consider, for example, the construction of such a function for a rectangular region.

The rectangular region is bounded by lines x = 0, x = a, y = 0, y = b . Let us introduce the supporting

regions: Z1 =j(x,y):—x(a-x) >0! is vertical band

between straight lines x = 0 1

and

x = a ;

Z2 = j(x,y): ^y(b - y) > 0 j is horizontal band between

straight lines y = 0 and y = b . Then it can be seen that the rectangular region has the form

Q = Z1 n Z2.

Considering Z1 and Z2 as two-valued predicates, we can write the predicate equation of the domain Q

az 2 = 1.

Substituting instead of Z and Z2 the left-hand sides of the inequalities defining them and Ro -conjunction a0 instead of the conjunction a , we obtain the function

®(x, y) =

1

a

- x(a - x)

A0

y(b - y)

=1 x(a - x) +1 y(b - y) -a b

-2 x2(a - x)2 +\y 2(b - y)2 a2 b2

(11)

ra( x, y) = 0 on dQ, ra(x,y) > 0 in Q , ra( x, y) < 0 outside Q .

(12)

90(s) =

90 1(s), s e 9Q1, 90, n(s), s edQn,

dQn are

where the sections of the boundary pair-wise distinct, do not have common interior points and 9Q = 5Q1 U-.3Qn . Further, let <i(x,y), i = 1,..., n ,

be such that <i |5Q = <o, i, and rai (x, y), i = 1,..., n , be

such that (x, y) = 0 on dQ and (x, y) > 0 in Q \ dQ,-. Then the function

9 =

91+.. + 9n_ " ®n Z9i n® j i=1 j=1 j *'

1 —+.. + 1 ®n nn in® j ,=1j=1 j *'

(13)

The constructed function satisfies the following properties

In addition, due to the introduction of constant multipliers into the left-hand sides of the inequalities defining the support domains, the function (12) has the normalization property: at regular points of 9Q the

5ra , , ,, equality — = -1 holds. dn

Applicable to the solution of boundary and initial-boundary value problems of mathematical physics, the R-functions method proposes to construct a so-called solution structure - a bundle of functions that depends on one or more indeterminate components, which for any choice of undefined components exactly satisfy all the boundary conditions of the problem. After constructing the solution structure, the indefinite component can be approximated accurately by some variational or projection method [17]. Thus, RFM is a mesh-free method, that allows to present an approximate solution of the problem in the form of a single analytical expression and accurately takes into account the geometry of the domain in which the task is considered.

The construction of the solution structure usually includes the stage of continuation of the boundary conditions inside the domain. For this, basically the following two approaches are used [17].

Suppose that at points of 9Q the function <o is given in the form

has the property 9|5Q =90. The formula (13) is called the «gluing» formula and is denoted by 9 = EC90 .

Another approach is connected with the continuation of differential operators defined on dQ inside the domain. Let ra = 0 be the normalized equation of the boundary dQ of the domain Q . Then the operator D1 at regular points of dQ satisfies the equality

-D1u\eQ = dn . In addition, the expression D1m has a

meaning everywhere in Q .

It is known [3; 17] that the boundary conditions (9) are satisfied by the bundle of functions

9 = f-®( g + DJ) + ra2®,

(14)

where f = EC/q , g = ECgo .

Method of numerical analysis of steady creeping flow. Using the solution structure (14) in the problem (8), (9) let us make the substitution

9 = 9 + u.

(15)

where 9 = f -ra(g + D1 f), u is a new unknown function. Then for the function u we will obtain a boundary value problem with homogeneous boundary conditions

2

A u = f , (x, y) e Q , =0,

I n du u = 0, — dn

(16) (17)

SQ

where f = -A 9.

Let us associate with the problem (16), (17) the operator A acting in by the rule Au = A2u on the

domain of definition

Da =\u eC4(Q)PIC1(Q):u

I = 0 —

Q ' dn

= 01

SQ

The operator A is positive definite on DA . Replenishing Da in the energy norm |UIIa = >/[u, u]a generated by the energy product [u, v]a = jj Au Avdxdy ,

Q

we obtain the energy space Ha of the operator A . Hence, the problem of finding a generalized solution of the problem (16), (17) will be equivalent to the problem of finding the minimum in HA of the energy functional [3; 18]

F[u] = |u||A - 2(u, f) = jj[(Au)2 + 2AuA|]dxdy . (18) Q

In addition, it follows from (14) and (15) that the structure of the solution of the boundary value problem (16), (17) has the form

u = ra2®.

(19)

N

®N = Z ck Tk . k=1

(20)

Then it follows from (19) that we seek the approximate solution of problem (16), (17) in the following form

N

N

u

N =ra Z ck Tk = Z ck 9k k=1 k=1

(21)

where 9k = Tk .

It is obviously that the sequence {9k } is a coordinate one [19]:

1) for anyone k 9k e Ha ;

2) for any N the elements 91, ..., 9N are linearly independent;

3) the system {9k } is complete in HA .

Then to find the coefficients c1, ..., cN in (21) we can use the Ritz method [18], according to which we obtain that c1 , . , cN is a solution of a system of the linear algebraic equations

N

Z ci [9i, 9j]A = -(A9 a9j), j = 1

, N,

(22)

i =1

where (A9, A9 j ) = JJA9 A9 jdxdy .

From the general theorems on the convergence of the Ritz method [18] the ensuing theorem follows.

Theorem 1. If 9 = f - ra(g + D1 f) e W22 (Q), then the

sequence yN = 9 + uN converges in the norm of W22 (Q) to a unique generalized solution of the problem (16), (17).

Method of numerical analysis of unsteady creeping flow.

In the problem (5)-(7) let us also make the substitution of the form (15) and for the function u we

obtain the initial problem with homogeneous boundary conditions

5(-Au). + vA2u = f , (x, y) eQ, t > 0, (23)

dt

u\t=0 = u0, ( x, y) eQ,

du

=0, t>0,

I n du

uL~> = 0, —

Q dn

(24)

(25)

dQ

Let us choose a complete in L2(Q) sequence of functions {Tk} and discover an approximation of the indeterminate component ® in the form

where f = -vA29 + dA^^, u0 = ^0 -9|t=0.

With the problem (23)-(25) for t e [0, T] (T > 0) let us associate the operators A and B acting in L2 (0, T; L2 (Q)) by the rules Au = A u , Bu = -Au on the domains of definition

Da =\u eC4(Q)nC1(Q):u|dQ = 0,^^

DB =\u eC2(Q)nC1(Q):u|„ = 0, —

dn

= 0!

dQ

= 0!

dQ

The operators A and B are positive definite on DA and DB respectively. Replenishing Da and DB in the energy norms ||u||a = <J[u, u]a , |u||B = V[u, u]B ,

generated by the energy products [u, v]a = jjAu Avdxdy ,

Q

[u, v]B = jj Vu Vvdxdy , we obtain the energy spaces Ha

Q

and HB .

Then the problem (23)-(25) can be written in the operator form

—Bu +vAu = f , t > 0, dt

u\t=0 = u0

(26) (27)

The structure of the solution of the initial-boundary value problem (23)-(25) also has the form (19), but since the problem is non-stationary, the approximation of the indeterminate component ® in contrast to (20) is sought in the form

N

®N = Z ck (t)Tk . k=1

(28)

Hence, we seek the approximate solution of the problem (23)-(25) in the form

N

uN = Z ck(t)9k k=1

(29)

where 9k =ra Tk are elements of coordinate sequence.

Q

Then to find the functions c1(t), ., cN (t) in (29) one can use the Galerkin method for non-stationary problems [20], according to which these functions satisfy the system of differential equations obtained from the condition of the orthogonality of the residual

Rn = —BuN +vAuN - f for equation (21) to the first

dt

N coordinate functions 91 , ., 9N :

N N

Z ci(t) [9i, 9j ]B + vZ ci(t) [9i, 9 j]A = i=1 i=1

= (f, 9j ), j = 1,..., N, (30)

where ( f, 9 j ) = jj(-vA9A9 j - V9 j j dxdy .

The initial conditions for the system (30) are obtained from the requirement of the orthogonality of the discrepancy rN = uN |f=0 - u0 of the initial condition (29)

to the first N coordinate functions 91 , ., 9N :

Cauchy problem (30), (31) was solved by the Runge-Kutta-Merson method with an automatic choice of the

integration step with accuracy 10-6.

In the stationary problem (8), (9), the boundary conditions were chosen in the form

IdQ

= 0.

dy dn

dQ

-1, 0,

y = b,

dQ \{ y = b},

which corresponds to the solid fixed walls of the cavity Q and to the movement of the upper cover to the left with unit velocity.

In the non-stationary problem (5)-(7) ( v = 1 ), the initial condition is chosen in the form

\t=0

= 0.

which corresponds to the beginning of the motion from the rest state, and the boundary conditions have the form

4rc= 0, -

dn

dQ

kf -1,

I 0,

y = b,

dQ \{ y = b},

N

ZCi (0)(9i, 9j) = (u0, 9j), j = 1,..., N . (31)

i=1

The unique solvability of the Cauchy problem (30), (31) for any N [21] follows from the positive definiteness of the operators A , B and the properties of the coordinate sequence. And from the theorems of convergence of the Galerkin method [20; 22] the next theorem follows.

Theorem 2. If u0 = y0 -9I =0 e L2(Q) and

f = -vA29 + dA9e L2(0, T; L2(Q)), then the sequence

dt

yN = 9 + uN converges in the norm of L2(0, T; W22 (Q)) to the unique generalized solution of the problem (5)-(7). 4 EXPERIMENTS A computational experiment was carried out for three rectangles regions with various aspect ratios: a = b = 1; a = 1, b = 2 and a = 1, b = 0.5 .

The undefined component in the solution structure (14) was approximated by an expression of the form

Nx +3 Ny +3

(

®N = Z Z CijB6

i=-3 j=-3

V hx

j ^ /

V hy

B6

ab

where hx = —, hy =- are steps of a splines grid by

x Nx y Ny V V 6 "

variables x and y accordingly.

Scalar products in systems (22) and (30), (31) were calculated with accuracy 10-6 according to the adaptive procedure of numerical integration based on the Gauss formula. The system (22) was solved by the Gauss method with the choice of the principal element, and the

which also corresponds to the solid impermeable walls of the cavity Q and the motion of the upper wall to the left

from the rest state with the speed e- -1.

For the stationary problem, in accordance with the «gluing» formula (13), we have:

f = 0, g = -

xy(a - x)

b - y + xy(a - x)

and for unsteady flow

f = 0, g =

(e 1 -1) xy(a - x) b - y + xy(a - x)

Then the problem solution structure for the domain has the form

y = ra

xy(a - x)

b - y + xy(a - x)

„-t\

+ ra ® (steady problem),

(1 - e )xy(a - x) 2.,. - . . .

y = ra---+ ra ® (unsteady problem),

b - y + xy(a - x)

where the function ra(x,y) is defined by (10).

For the pressure calculation, the following formula was used:

P(x,y) = j ^dx + dpdy + C =

./.. dx dy m 0m j

= VP J

M 0M

d y + d y dx2dy dy3

dx -

( £¡3 r3 \ d y + d y

dx3 dxdy2

dy + C.

5 RESULTS

Figure 1, 2 and 3 show the pressure contours of the obtained approximate solution for the rectangles a = b = 1; a = 1, b = 2 and a = 1, b = 0.5 respectively.

6 DISCUSSION

For small Reynolds numbers (v = 1 corresponds to Re = 1), the flow is symmetrical and slow. It can be seen from the figures 1-3 that with the increasing of time from 1 to 10, the pressure in the cavity has become steady. The pressure contours for unsteady flow are the same as for unsteady one when t = 10.

The maximum pressure of the flow is located at the top left corner of the cavity, and the minimum pressure occurs at the top right corner of the cavity. This behavior corresponds to the fluid movement from right to left on the lid. The pressure is decreases at the lower part of the cavity.

Thus, the stream function is an intermediary between the Navier-Stokes equation and obtaining the various characteristics of flow: speed, pressure etc. Other numerical results are shown in [23].

CONCLUSIONS

The solution methods for steady and unsteady flows in the driven cavity were further developed. These algorithms are based on Ritz or Galerkin method and the R-functions method, which allows satisfying the boundary conditions precisely. We compared steady and unsteady flow and it has been shown that unsteady flow

1.01

0.8

O. I

0.2

1.01

0.8

0.6

0.4

0.2h

0.01

\ N \ / 7 (

10 ) ntrv

-J 5 f -5

j f— \

\ 1 / \ 1 1

(1 i » 1 1 0

\ J

1.0 0.2 0.4 0.6 0.S 1.0 a

A 10 s \ N V 1 / / 'P- 20 [ o L

1 J f \ 1-

1 1 -S U

j v

1 Y > 1 1

J i /

\ /

1 1 ' 1 1 f]

ff №

0.0 0.2 0.4 0.6 0.S 1.0 c

becomes steady when the time increases. In addition, the various numerical experiments were carried out.

The scientific novelty of the research is the obtained methods do not change when the region is changed. The only thing one should do is to build a normalized boundary equation of a new flow domain. This benefit allows researches to get the solution faster and easier.

The practical significance is that the developed methods can be programmed easy and can be used for several practical applications. The only restriction is it must be possible to represent a boundary of flow domain by means of R-functions.

The prospects for further research are to investigate more complicated flow domains, for example, domains with inserts, to use other coordinate function with the aim of obtaining more precise solution and to reduce computation time with the help of parallel computing on CPU and GPU.

ACKNOWLEDGEMENTS The work was carried out at the department of Applied Mathematics at Kharkiv National University of Radio Electronics and department of Advanced Mathematics at O. M. Beketov National University of Urban Economy in Kharkiv within the framework of collaborative scientific research conducted by the departments.

1.0 F

0.8'

0.6

0.4

0.2

0.01

r#l H i(> \ N \ / / ( 7- fh -10 H

1 ) \ 1

MV 1

/ / v

\ I i v 1 h À

1 1 i —-

\

,1 0 1 1 31 1 0

[ i

1.0 0.2 0.4 0.6 0.S b 1.0

0.0 0.2 0.4 0.6 0.S 1.0 d

Figure 1 - Pressure contours for the rectangle a = b = 1 : a - t = 1, b - t = 3, c - t = 5, d - t = 10

0.0 0:2 0.4 O.G 0.S 1,0

a

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.1 0 Ci 0.S 1.0

C d

Figure 2 - Pressure contours for the rectangle a = 1, b = 2 :

a - t = 1, b - t = 3, c - t = 5, d - t = 10

(1.(11 L

u

u

(Mi

0.2

0.4 0.6

a

0.5F

0.4

0.3

0.2

0.1

o.oL

C±

<1.5

0 ;i

o.i

(I

oL

0,0

0.2

0.4

0.0

0.8

1.0 0.0

0.2

0.4 0.0

c d

Figure 3 - Pressure contours for the rectangle a = 1, b = 0.5 : a - t = 1, b - t = 3, c - t = 5, d - t = 10

o.s

1.0

3.

4.

REFERENCES

Landau L. D., Lifshits E. M. Teoreticheskaya fizika. V 10 t. T. VI. Gidrodinamika. Moscow, Fizmatlit, 2003, 736 p. Ladyizhenskaya O.A. Matematicheskie voprosyi dinamiki vyazkoy neszhimaemoy zhidkosti. Moscow, Nauka, 1970, 288 p.

Sidorov M.V. O postroenii struktur resheniy zadachi Stoksa, Radioelektronika i informatika, 2002, No. 3, pp. 52-54. Jiangfei L., Long J., Lian Y., Zhizhong F., Bo L., Wenxue C. Comparison of Finite Difference and Finite Volume Method for Numerical Simulation of Driven Cavity Flow Based on MAC, 2013 International Conference on Computational and Information Sciences, 2013, pp. 891894. DOI: 10.1109/ICCIS.2013.239. Bettaibi S., Kuznik F., Sediki E. Hybrid lattice Boltzmann finite difference simulation of mixed convection flows in a lid-driven square cavity, Physics Letters A, 2014, V. 378, 32-33, pp. 2429-2435. DOI:

10.1016/j.physleta.2014.06.032.

Wu Y., Wang H. Moving Mesh Finite Element Method for Unsteady Navier-Stokes Flow, Advances in Applied Mathematics and Mechanics, 2017, No. 9(3), pp. 742-756. DOI: 10.4208/aamm.20-16.m1457.

Li J., He Y., Chen Z. Performance of several stabilized finite element methods for the Stokes equations based on the lowest equal-order, Computing, 2009, No. 86, pp. 37-51. DOI: 10.1007/s00607-009-0064-5.

Berrone S., Marro M. Space-time adaptive simulations for unsteady Navier-Stokes problems, Computational Fluids, 2009, No. 38, pp. 1132-1144. DOI: 10.1016/j.compfluid.2008.11.004.

Patankar S. Numerical heat transfer and fluid flow. Washington, Hemisphere Pub. Corp, 1980. Rvachev V.L. Ob analiticheskom opisanii nekotoryih geometricheskih ob'ektov, Dokl. AN SSSR, 1963, Vol. 153, No. 4, pp. 765-768.

9.

10.

11. Earn L. C., Yen T. W., Ken T. L. The investigation on SIMPLE and SIMPLER algorithm through lid driven cavity Akademia Baru, 2017, No. 1, pp. 10-22.

12. Gupta M. M., Kalita J. C. A new paradigm for solving Navier-Stokes equations: Streamfunction-velocity formulation, Journal of Computational Physics, 2005, No. 207, pp. 52-68. DOI: 10.1016/j.jcp.2005.01.002.

13. Claeyssen J. R., Platte R. B., Bravo E. Simulation in primitive variables of incompressible flow with pressure Neumann condition, International Journal for Numerical Methods in Fluids, 1999, No. 30, pp. 1009-1026. DOI: 10.1002/(SICI)1097-0363(19990830)30:8<1009::AID-FLD876>3.0.C0;2-T.

14. Bognar G., Csâti Z. Spectral method for time dependent Navier-Stokes equations, Miskolc Mathematical Notes, 2016, No. 17, pp. 43-56. DOI: 10.18514/mmn.2016.1815.

15. Bettaibi S., Sediki E. Numerical simulation of mixed convection flows in lid-driven square cavity, Fluid Mechanics and Thermodynamics: proceedings of the 10th International Conference on Heat Transfer, 2014, pp. 967973.

16. Kravchenko V. F., Rvachev V. L. Algebra logiki, atomarnyie funktsii i veyvletyi v fiziche-skih prilozheniyah. Moscow, Fizmatlit, 2006, 416 p.

17. Rvachev V. L. Teoriya R-funktsiy i nekotoryie ee prilozheniya. Kiev, Nauk. dumka, 1982, 552 p.

18. Mihlin S.G. Variatsionnyie metodyi v matematicheskoy fizike. Moscow, Nauka, 1970, 512 p.

19. Kantorovich L. V., Kryilov V. I. Priblizhennyie metodyi vyisshego analiza. Leningrad, Fizmatgiz, 1962, 708 p.

20. Mihlin S. G. Chislennaya realizatsiya variatsionnyih metodov. Moscow, Nauka, 1966, 432 p.

21. Artiukh A. V., Sidorov M. V. Primenenie metodov R-funktsiy i Galerkina k raschetu ploskih nestatsionarnyih vyazkih techeniy, Visnik Zaporizkogo na-tsionalnogo

b

universitetu. Seriya: fiziko-matematichni nauki, 2011, No. 2. pp. 5-12.

22. Savula Ya. Chisloviy anallz zadach matematichnoi fIziki variatsiynimi metodami. Lviv, Vidavnichiy tsentr LNU Im. I. Franka, 2004, 224 p.

23. Artiukh A. V., Sidorov M. V. Issledovanie nestatsionarnyih ploskoparallelnyih techeniy vyazkoy neszhimaemoy zhidkosti (priblizhenie Stoksa) metodami R-funktsiy i Galerkina, Radioelektronika i informatika, 2011, No. 3 (54), pp. 16-21.

Received 16.10.2018.

Accepted 12.12.2018.

УДК 004.93

ЧИСЕЛЬНИЙ АНАЛ1З ПОВ1ЛЬНИХ СТАЦЮНАРНИХ ТА НЕСТАЦЮНАРНИХ В ЯЗКИХ ТЕЧ1Й МЕТОДОМ R-

ФУНКЦ1Й

Артюх А. В. - канд. фiз.-мат. наук, доцент кафедри прикладно! математики Харквського нацiонального ушверситету радiоелектронiки, Харкiв, Укра!на.

Ламтюгова С. М. - канд. фiз.-мат. наук, доцент кафедри вищо! математики Харкiвського нацiонального унiверситету мгського господарства iменi О. М. Бекетова, Харкв, Укра!на.

Сидоров М. В. - канд. фiз.-мат. наук, доцент, доцент кафедри прикладно! математики Харювського национального унiверситету радiоелектронiки, Харкв, Украша.

AНОТАЦIЯ

Актуальнiсть. У статтi розглядаеться лiнiйна стацiонарна i нестацiонарна задача течи в'язко! нестисливо! рiдини.

Мета роботи. Пор1вняння розроблених метод1в чисельного аналiзу стацiонарних та нестацiонарних задач течи в'язко! нестисливо! рщини.

Метод. Течiя в'язко! нестисливо! рiдини описуеться нелiнiйною системою Нав'е-Стокса вiдносно вектора швидкосп, тиску, густини, об'емних сил та в'язкосп рiдини. За допомогою введення функцi! течi! система зводиться до початково-крайово! задачi з нелшшним диференцiальним р1внянням четвертого порядку. Для розв'язання задачi використовуються структурно-варiацiйний метод R-функцiй i методи Р^ща та Гальоркiна для стацiонарно! та нестацюнарно! течi! вщповщно. Метод Я-функцш дозволяе точно врахувати крайовi умови та звести !х до однорщних, що е необидною умовою можливосл застосування методiв Рiтца та Гальоркша. У стацiонарному випадку задача зводиться до розв'язання системи лiнiйних алгебра!чних р1внянь, у нестацiонарному - системи звичайних диференщальних р1внянь. Елементами матриць е скалярнi добутки у нормах вщповщних диференщальних оператор1в. Чисельне штегрування виконувалось за допомогою квадратурних формул Гаусса з 16 вузлами, розв'язання системи лшшних алгебра!чних рiвнянь - методом Гаусса, системи звичайних диференщальних ршнянь - методом Рунге-Кутти з автоматичним вибором кроку штегрування. Доведено юнування единого розв'язку поставлених задач.

Результати. Обчислювальний експеримент проведено для задачi течi! в'язко! нестисливо! рiдини у рiзних прямокутних областях. Наведет графжи лшш р1вня тиску та виконано аналiз отриманих результатв, яю добре погоджуються з результатами, отриманими iншими авторами.

Висновки. Проведенi обчислювальш експерименти показали, що при збшьшенш часу функцiя течi!, вектор швидкосп та iншi характеристики течi! виходять у стацiонар, що пiдтверджуе працездатшсть метод1в. В подальшому плануеться поршняння нелшшних метод1в розв'язання нелiнiйних задач.

КЛЮЧОВ1 СЛОВА: р1вняння Нав'е-Стокса, стащонарна течш, нестацiонарна течiя, в'язка рщина, функщя течi!, метод R-функцiй, метод послщовних наближень, метод Рiтца, метод Гальоркша.

УДК 004.93

ЧИСЛЕННЫЙ АНАЛИЗ МЕДЛЕННЫХ СТАЦИОНАРНЫХ И НЕСТАЦИОНАРНЫХ ВЯЗКИХ ТЕЧЕНИЙ

МЕТОДОМ R-ФУНКЦИЙ

Артюх А. В. - канд. физ.-мат. наук, доцент кафедры прикладной математики Харьковского национального университета радиоэлектроники, Харьков, Украина.

Ламтюгова С. Н. - канд. физ.-мат. наук, доцент кафедры высшей математики Харьковского национального университета городского хозяйства имени А. Н. Бекетова, Харьков, Украина.

Сидоров М В. - канд. физ.-мат. наук, доцент, доцент кафедры прикладной математики Харьковского национального университета радиоэлектроники, Харьков, Украина.

AННОТАЦИЯ

Актуальность. В статье рассматривается линейная стационарная и нестационарная задача течения вязкой несжимаемой жидкости.

Цель работы. Сравнение разработанных методов численного анализа стационарных и нестационарных задач течения вязкой несжимаемой жидкости.

Метод. Течение вязкой несжимаемой жидкости описывается нелинейной системой Навье-Стокса относительно вектора скорости, давления, плотности, объемных сил и вязкости жидкости. С помощью введения функции тока система сводится к начально-краевой задаче с нелинейным дифференциальным уравнением четвертого порядка. Для решения задачи используются структурно-вариационный метод Я-функций и методы Ритца и Галеркина для стационарного и нестационарного течения соответственно. Метод Я-функций позволяет точно удовлетворить краевым условиям и свести их © АйшкИ А. V., Lamtyugova 8. К., Sidorov М. V., 2019 БО! 10.15588/1607-3274-2019-1-3

к однородным, что является необходимым условием применимости методов Ритца и Галеркина. В стационарном случае задача сводится к решению системы линейных алгебраических уравнений, в нестационарном - системы обыкновенных дифференциальных уравнений. Элементами матриц являются скалярные произведения в нормах соответствующих дифференциальных операторов. Численное интегрирование выполнялось с помощью квадратурных формул Гаусса с 16 узлами, решение системы линейных алгебраических уравнений - методом Гаусса, системы обыкновенных дифференциальных уравнений - методом Рунге-Кутты с автоматическим выбором шага интегрирования. Доказано существование единственного решения поставленных задач.

Результаты. Вычислительный эксперимент проведен для задачи течения вязкой несжимаемой жидкости в различных прямоугольных областях. Приведены графики линий уровня давления и выполнен анализ полученных результатов, которые хорошо согласовываются с результатами, полученными другими авторами.

Выводы. Проведенные вычислительные эксперименты показали, что при увеличении времени функция тока, вектор скорости и другие характеристики течения выходят в стационар, что подтверждает работоспособность методов. В дальнейшем планируется сравнение методов решения нелинейных задач.

КЛЮЧЕВЫЕ СЛОВА: уравнение Навье-Стокса, стационарное течение, нестационарное течение, вязкая жидкость, функция тока, метод R-функций, метод последовательных приближений, метод Ритца, метод Галеркина.

.MTEPATyPA / .HTEPATyPA

1. HaHgay H. fl. TeopeTHHecKaa $H3HKa. B 10 t. T. VI. rugpogMHaMHKa / H. fl. HaHgay, E. M. Hu^mu^ - M. : OroMainuT. - 2003. - 736 c.

2. HagtraeHCKaa O. A. MaTeMaTHHecKHe Bonpocti guHaMHKH b«3koh HecMMaeMoft rngKocTH / O. A. HagtraeHCKaa. -M. : HayKa, 1970. - 288 c.

3. CugopoB M. B. O nocipoeHHH cTpyKTyp pemeHHH 3agaHH CTOKca / M. B. CugopoB // Pagнoэneктpoннкa h HH^opMaTMKa. - 2002. - № 3. - C. 52-54.

4. Jiangfei L. Comparison of Finite Difference and Finite Volume Method for Numerical Simulation of Driven Cavity Flow Based on MAC / [L. Jiangfei, J. Long, Y. Lian et al.] // 2013 International Conference on Computational and Information Sciences, 2013. - P. 891-894. DOI: 10.1109/ICCIS.2013.239.

5. Bettaibi S. Hybrid lattice Boltzmann finite difference simulation of mixed convection flows in a lid-driven square cavity / S. Bettaibi, F. Kuznik, E. Sediki // Physics Letters A. - 2014. - V. 378, 32-33. - P. 2429-2435. DOI: 10.1016/j.physleta.2014.06.032.

6. Wu Y. Moving Mesh Finite Element Method for Unsteady Navier-Stokes Flow / Y. Wu, H. Wang // Advances in Applied Mathematics and Mechanics. - 2017. - № 9(3). -P. 742-756. DOI: 10.4208/aamm.20-16.m1457.

7. Li J. Performance of several stabilized finite element methods for the Stokes equations based on the lowest equal-order pairs / J. Li, Y. He, Z. Chen // Computing. - 2009. -№ 86. - P. 37-51. DOI: 10.1007/s00607-009-0064-5.

8. Berrone S. Space-time adaptive simulations for unsteady Navier-Stokes problems / S. Berrone, M. Marro // Computational Fluids. - 2009. - № 38. - P. 1132-1144. DOI: 10.1016/j.compfluid.2008.11.004.

9. Patankar S. Numerical heat transfer and fluid flow / S. Patankar. - Washington : Hemisphere Pub. Corp. - 1980.

10. PBaHeB B. H. O6 aHanHTHHecKoM onucaHHH HeKOToptix reoMeTpHHecKHX o6teKTOB / B. H. PBaHeB // flora. AH CCCP. - 1963. - T. 153, № 4. - C. 765-768.

11. Earn L. C. The investigation on SIMPLE and SIMPLER algorithm through lid driven cavity / L. C. Earn, T. W. Yen, T. L. Ken // Akademia Baru. - 2017. - № 1. - P. 10-22.

12. Gupta M. M. A new paradigm for solving Navier-Stokes equations: Streamfunction-velocity formulation / M. M. Gupta, J. C. Kalita // Journal of Computational

Physics. - 2005. - № 207. - P. 52-68. DOI: 10.1016/j.jcp.2005.01.002.

13. Claeyssen J. R. Simulation in primitive variables of incompressible flow with pressure Neumann condition / J. R. Claeyssen, R. B. Platte, E. Bravo. // International Journal for Numerical Methods in Fluids. - 1999. - № 30. -P. 1009-1026. DOI: 10.1002/(SICI)1097-0363(19990830)30:8<1009::AID-FLD876>3.0.CO;2-T.

14. Bognar G. Spectral method for time dependent Navier-Stokes equations / G. Bognar, Z. Csati // Miskolc Mathematical Notes. - 2016. - № 17. - P. 43-56. DOI: 10.18514/mmn.2016.1815.

15. Bettaibi S. Numerical simulation of mixed convection flows in lid-driven square cavity / S. Bettaibi, E. Sediki // Fluid Mechanics and Thermodynamics: proceedings of the 10th International Conference on Heat Transfer, 2014. - P. 967973.

16. KpaBHeHKO B. O. Anre6pa norHKH, aTOMapHtie ^yHKUHH h BeHBneTti b $H3HHecKHx npHnoxeHHax / B. O. KpaBHeHKo,

B.H. PBaHeB. - M. : OH3MaTnHT, 2006. - 416 c.

17. PBaHeB B. H. Teopua R-^yHKunft h HeKOToptie ee npHroxeHHfl / B. H. PBaHeB. - K. : HayK. gyMKa, 1982. -552 c.

18. MuxnuH C. r. BapнaцнoннвIe MeTogti b MaTeMamHecKOH $H3HKe / C. r. MuxnuH. - M. : HayKa, 1970. - 512 c.

19. KaHTopoBHH H. B. npH6nH:eHHBie MeTogti Bticmero aHanroa / H. B. KaHTopoBHH, B. H. KptinoB. - H.: OroMairro, 1962. - 708 c.

20. MuxnuH C. r. HucneHHaa peanнзaцнa вapнaцнoннвIX MeiogoB / C. r. MuxnuH. - M. : HayKa, 1966. - 432 c.

21. ApTrox A.B. npHMeHeHHe MeiogoB R-^yHKUHH h ranepKHHa K pacHeTy nnocKHx нecтaцнoнapнвIX ba3khx TeHeHHH / A. B. ApTrox, M. B. CugopoB // BicHHK 3anopi3tKoro Ha^oHantHoro ymBepcmeTy. Cepia: $i3HKo-MaTeMaTHHHi HayKH. - 2011. - № 2. - C. 5-12.

22. CaByna HucnoBHH aHani3 3agaH MaTeMaTHHHoi $i3HKH Bapia^HHHMH MeTogaMH / CaByna. - Hbbib : BugaBHHHHH цeнтp HHy iM. I. OpaHKa, 2004. - 224 c.

23. ApTrox A. B. HccnegoBaHue нecтaцнoнapнвIX nnocKonapannentHtix TeHeHHH b«3koh HecrnMaeMoft :ugKocTH (npu6nH:eHHe CToKca) MeTogaMH R-^yHKUHH h ranepKHHa / A. B. ApTrox, M. B. CugopoB // Pagнoэneктpoннкa h HH^opMaTHKa. - 2011. - № 3 (54). -

C. 16-21.