Полуаналитический метод расчета течения вязкой несжимаемой жидкости в канале постоянной ширины Текст научной статьи по специальности «Математика»

УДК 532.5: 519.6

ПОЛУАНАЛИТИЧЕСКИЙ МЕТОД РАСЧЕТА ТЕЧЕНИЯ ВЯЗКОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ В КАНАЛЕ ПОСТОЯННОЙ ШИРИНЫ

© 2014 г. М.А. Сумбатян, В.В. Абрамов

Сумбатян Межлум Альбертович - доктор физико-математических наук, профессор, заведующий кафедрой теоретической и компьютерной гидроаэродинамики, факультет математики, механики и компьютерных наук, Южный федеральный университет, ул. Мильчакова, 8а, г. Ростов н/Д, 344090, e-mail: sumbat@math.sfedu.ru.

Абрамов Владимир Владимирович - аспирант, кафедра теоретической и компьютерной гидроаэродинамики, факультет математики, механики и компьютерных наук, Южный федеральный университет, ул. Мильчакова, 8а, г. Ростов н/Д, 344090, e-mail: abram5189@yandex.ru.

Sumbatyan Mezhlum Albertovich - Doctor of Physical and Mathematical Science, Professor, Head of the Department of Theoretical and Computer Hydroaerodynamics, Faculty of Mathematics, Mechanics and Computer Sciences, Southern Federal University, Milchakov St., 8a, Rostov-on-Don, Russia, 344090, e-mail: sumbat@math.sfedu.ru.

Abramov Vladimir Vladimirovich - Post-Graduate Student, Department of Theoretical and Computer Hydroaerodynamics, Faculty of Mathematics, Mechanics and Computer Sciences, Southern Federal University, Milchakov St., 8a, Rostov-on-Don, Russia, 344090, e-mail: abram5189@yandex. ru.

Предлагается новый численно-аналитический метод решения двумерной задачи о течении вязкой несжимаемой жидкости в канале постоянной ширины. Метод основан на неявной разностной схеме по времени для уравнений Навье-Стокса, который приводит к итерационному процессу, где на каждом временном шаге решаются эллиптические краевые задачи для функции тока и завихренности. Эти решения строятся через функции Грина, полученные в виде Фурье-разложений по поперечной координате. Граничные условия приводят к интегральным уравнениям с разностным ядром для функции завихренности на стенках канала, что достигается с помощью быстрого преобразования Фурье. Приведен пример расчета в ламинарном режиме.

Ключевые слова: функция тока, завихренность, уравнения Навье-Стокса, численно-аналитический метод, неявная конечно-разностная схема, функция Грина, быстрое преобразование Фурье.

In the paper we propose a new numerical and analytical method to solve a two-dimensional problem on a flow of viscous incompressible fluid in a channel of constant width. The method is based on an implicit finite-difference scheme in time for the Navier-Stokes equations, which reduces the problem to an iteration process, where at each temporal step there are solved elliptic boundary value problems for stream function and vorticity. These solutions are constructed in terms of Green's function obtained as the Fourier expansions over transverse coordinate. The boundary conditions reduce to some integral equations with convolution kernels for vorticity over the faces of the channel, which are further solved by fast Fourier transform. An example of the calculations is demonstrated in the laminar regime.

Keywords: stream function, vorticity, Navier-Stokes equations, numerical and analytical method, implicit finite-difference scheme, Green's function, fast Fourier transform.

Математическая постановка задачи

Рассмотрим классическую задачу об однородном потоке несжимаемой вязкой жидкости в канале постоянной ширины Ь в двумерном случае. Уравнения движения Навье-Стокса в терминах «завихренность» -«функция тока» имеют вид

дС дс дс

— = -и — -V— + vAC , С = Ду,

dt дх ду

(1)

где компоненты вектора скорости частиц жидкости и(х,уД), у(х,уД) выражаются в виде

^ _ ду ду Q _ди dv ду ' дх ' ду дх

(2)

Q = \ u(x,y)dy = уу=ь -уу=0.

0

Уравнения (1), (2) и условия прилипания жидкости на стенках канала дают

дС

д/ С = Ау

-vAC= /, f =

ду дАу ду дАу

дх ду ду дх

у у=0 = 0 , у у=b = Q

ду ¥

у=0

ду ¥

= 0

(3)

(4)

(5)

у=b

Координата х направлена вдоль нижней границы канала х е (-), а вертикальная координата у -поперек канала у е (0;Ь). Считаем, что известен расход жидкости через нормальное сечение канала, который не зависит от времени и координаты х :

В такой постановке к рассматриваемой задаче применялись различные численные методы [1, 2], которые допустимы как в ламинарном, так и в турбулентном режиме. Обзор современных численных методов можно найти в [3, 4].

Рассматриваемая задача (5) по переменной I является классической задачей Коши, если считать известными начальные условия для скоростей. Обычно

в развитом потоке такие условия редко бывают заданы заранее, однако численные эксперименты показывают, что выбор тех или иных начальных условий не меняет качественного характера течения. Известно [1, 2], что устойчивое численное решение можно получить различными итерационными процессами по времени, если использовать конечно-разностную аппроксимацию производной по времени в варианте «назад». При этом простейшей схемой является схема Эйлера

,5t,

(Сn n-l)

е

где е - шаг по времени. При-

Ay(n)-вД2у(п) = g„-1,

(6)

gn-i=c(n-1)+en-i.

Эта схема обладает линейной скоростью сходимости с уменьшением временного шага в. Заметим, что в (6) £ =УТ> 0 .

Таким образом, задача сводится к следующему виду:

C(n)-BAC(n) = gn

(n) _,

C(n) =Av(n),

V(n)

y=0

5v(n)

5у

= 0, V(n)

_öy(n)

y=b

= Q,

(7)

(8)

y=0

= 0.

y=b

Построение функций Грина

Для реализации метода необходимо построить функции Грина для функции тока и завихренности. Выделим в канале прямоугольную область Б шириной Ь и длиной Ь . Характерный размер Ь выбирается из условия однородности течения жидкости в канале, т.е. однородность течения можно понимать как его периодичность вдоль потока в достаточном удалении от начала канала с некоторым периодом Ь . Заметим, что тогда Ь/Ь >> 1. Для построения требуемых функций Грина рассмотрим два следующих дифференциальных уравнения с соответствующими граничными условиями:

AGv (#, л,х, У) = 5(# - х)5(л - У):

О^(#,0,х,y) = О^(#Ах,y) = 0 ,

Gw(0,q, х, y) = Gv (М, х, y),

SGV SGV --(0,Л,х,y) = —г(Мх,y) ,

(9)

S#

S#

Ос (#, л, х, y) -еДОс (#, л, х, y) = 5(# - х)5(л - y), (10) Ос (#,0, х, у) = Ос (#, ь, х, у) = 0, Gc (0, л, х, y) = Gc (L, л, х, у),

SG

SG

С

(0,Л,х,у) =—тСМхУ),

дд дд

где д е (0;Ь) и л е (0;Ь), а ху следует понимать как числа, которые изменяются в пределах х е (0;Ь) и У е(0;Ь).

Функции Грина (9), (10) строятся методом разделения переменных и имеют вид

» ът(Ьту)ъЩЬтл) еЬт+ еЬт(Ь-д-х)

°у=2--:——-, (11)

менение схем более точного порядка не меняет сути метода, при этом малый параметр £ > 0 (см. ниже) принимает несколько иной вид в зависимости от используемой схемы. Подстановка этого выражения в (3), (4) с выбором функции /, содержащей нелинейные члены на предыдущем временном слое, а члены со старшей четвертой производной от функции у -на новом временном слое, приводит к эллиптической краевой задаче 4-го порядка относительно функции тока:

m=1

, 7Ш

bm = — > b

Gc =-Z

1 - ebmL

Sin(bmУ)sin(bmЛ)

Am\#-х| + ßAm (.^\#-х

b£Ä„

1 - eAmL

Лт = ^ Ьт + 1 .

Заметим, что изложенные построения родственны методу, примененному в [5].

Построение решений для функции тока и завихренности

С целью сокращения выкладок продемонстрируем алгоритм построения решения только для функции тока, решение для завихренности строится аналогично. Умножим уравнение (9) на функцию тока и вычтем из него уравнение (7), умноженное на функцию Грина для функции тока. Проинтегрировав данное уравнение по области Б , применив теорему Грина и свойство свертки для дельта-функции, получим

J

SD

V(n)(M)G(#;л;х;y)- Gw(#;л;х;у)

Зу

(n)

Sn

Sn

-(#;л)

Lb

х а=у(п) (х; у)-цс(п) (д, (д, л, х, уадал. (12)

00

Применяя к левой части (12) условия (8) и однородные граничные условия для функции Грина, получим

ь осу (1 /—у (д, ь, х, у)ад =

0 дЛ

= у(п) (х; у) - л С(п) (д, (д, Л, х, у)адал. 00

Подставим (11) в полученное выражение. Рассмотрим его левую часть

Lb

(n)

L SGV

Qi—V(#,b,х,y)d# = --Q Z

0 5л П m=1

-Q f ny

-Q » (-1)m Sin(bmy)

m

-b

Таким образом, окончательно получаем V(n) (х; У) = Z V(mn) (х) Sin(bmy) + ^ ,

m=1 b

(13)

n

от

m=1

X

Y^Cx) = "

2b т (1 - e~bmL )

0

bm (|X-L) +g"bm

с m)(№

2 L

С(тП)(4) = "Г JC(n)(#,7)sin(bm7)d7. b 0

Рассуждая аналогично, получим формулы для завихренности

ад

C(n)(x,7) = ZCin)(x)sin(bmy) , (14)

сти)( x)=■

bÄm (1 - ekm )

L [(n) (#,0)- (-1)m [(#, b)\eÄmx + eÄm (L-li-xl)

1_L [

i kmL-, Je

bkm (1 - eXmL )0

eÄm\Z-A + ek (L-|i-x| )

?{mn-1)(4)d4:

g(m"-1)(4) = J g ("-1)(4,7)sin(bm^.

0

Итак, мы построили точные аналитические решения для функции тока и завихренности в виде (13), (14). Для нахождения оставшихся неизвестных значений завихренности на стенках канала подставим (14) в (13). После рутинных преобразований получаем

Y(n)(x, У) = I

1

(15)

2bebmlm (1-e~bmL )(1 - ekL )' L b

j gmn-1) wm (x, --—ТУ x

0 m 2bbm^m(1 -e m )(1 -eÄmL)

x L [(n) (r,0) - (-1) m [(г, bJ (x, \)d\ j sin(bmy) + ^

x, О =

L

= J 0

ebm (|x|-L) + e~bm x

Xm M + ^m (L-|x-[|)'

V {H-1][e

km + bm

XmL _ 1 bm(\-xl-L) + e-bmx\ 1 +

+ 1 - e-bmL ] [eÄm(-l^-x|+L) + ekm\r-x\ | у +

V {И -1] [

km - bm

Xm (-\\-x|+L) + eXm\-x\ | +

+ [eÄmL - 1] [ebm(\-x|-L) + e-bm\T-x\ I j _

;2 /,2 r "m

Xm - bm

-k

I kL 1 (\-x\-L) , \-x\

\e m -11 le ^ + e m

- ^ [1 - e bmL ] [ek(-|\-xl+L) + ^^m^x\ ] j.

Удовлетворим для функции (14) условию прилипания на стенках канала (5) и, производя элементарные преобразования, придем к следующей системе двух интегральных уравнений:

b2m

L ад

0 m=1 bX 2m (1 - e~bm )(1 - ek 2mL)

; [[(n)(\,0) -С(Х,Ъ)]/2m

x[- (\,0)-С(Т,Ъ)]720m (x,x)dx = 1

L ад

= JI

)m:

0 m=1 bek 2m (1 - e b2mL )(1 - eX2mL)

x g2m 1}(T)J2m (x, \)d\ +

2£ Ъ '

y2m+1

L ад

J I --■

0 m=0 bk2m+1 (1 - e ~b2m+1L )(1 - ehm+1L )

[(й) (\,0) + [(Г,b)] J2m+1 (x, Г)^Г =

= L

£ ад 1

= | ^ _1_х

0ш=0 ЬеХ2т+1 (1 - е"Ь2™+1£ )(1 - 2™ +1Ь)

х Я 2т~+1(т) 3 2ш+1( х>

Данные интегральные уравнения являются уравнениями Фредгольма первого рода с непрерывным разностным ядром. Решать их будем с помощью классического метода коллокаций, который переводит каждое из этих уравнений в систему линейных алгебраических уравнений (СЛАУ) размером N х N, где N - количество узлов сетки, взятых вдоль одной из стенок канала. При этом матрицы соответствующих СЛАУ имеют теплицеву структуру [6] и допускают решение с помощью быстрого преобразования Фурье [7].

Вычисление правой части для эллиптической краевой задачи

Обсудим алгоритм расчета функции яп-1 (6), известной из предыдущего шага:

d# = (16)

gn-1 =Ду(n-1) +

Г (п-1) 5Ау(п-1) (п-1) 5Ау(п-1)" дх ду ду дх ^

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

5 KY

dxk

-(x, У) =I

1

ril 2bekmbm (1 - e-bmL)(1 -ekmL)

:L gm--2)(\)5 kj°m( x,r)dr -

dxK

2bXmbm (1 - e bmL)(1 -eXmL)'

Lr Л5кJ0

J[[(n)(r,0) - (-1) m[(r, b)]J (x,T)dT +

dxK

2beXmbm (1 - e -bmL )(1 - ekL )

{gmn-2) (x) - bme [[(n) (x,0) - (-1)m [(x,b)]jj

1

L

x

X

e

m=1

X

x

x

0

e

ад

+

b

m

2

X

0

1

X

X

дк J 0 дк J 0 lim -m (х,т) - lim -m (х,т)

т^х-0 дх

к =1,2,3,

к

т^х+О дх'

>Sin(bmy) ,

/ J

О к I0

где производные-т(х, х) вычисляются явно из (16).

дх

Можно показать, что 1т(х, х) е С [0; Ь], поэтому разность пределов в (17) обращается в нуль. Таким образом, (17) окончательно может быть записано в

о у(и-1)

виде

к = 1,2,3, где

дх

(n-1)

ш дkw

Чх, У) = Z Vm (х) Sin(bmy) ,

m

=1 дх'

як (n-1)

д Vm

дхк

(х) =•

-bAmbm (1 - e ~bmL )(1 - eAmL )

L д к T 0

xJ gmn--)(T)-f- (х,т)Ст-

дх

0

-bAmbm (1 - e~bmL )(1 - eA mL )

X„ L4

l[c(

)

В итоге

x l[C(n)(T,0) - (-1)m 0

к(т, b)J

д к T 0 д J m

дхк

(х, т)с1т\.

g^-1)(#) =Z

д-V(n-1)

dVm (#)-bmvmn-1)(#)

дх

2

:Jsin(baл)sin(bmл)dл +eZ Zb

sv

(n-1)

к=1m=1

дх

-(#) X

Я- (n-1)

5 -Vm- (#) - bm v mn-1)(#)

-

b

(18)

X J cos(bm л) sin(bk л) sin(ba л)Сл -

0

-ez Z bV-1^)

к=1m=1

Xn-1)

- bm

дх3

дх

J cos(bkЛ) sin(bm л) sin(b ал)Сл -

eQ,

X Z

m=1

где

3,.,(n-1)

5 V

дх

-(#) - b-

Sy1

(n-1)

дх

b [0,m Фа

Jsin(bаЛ)sin(bmЛ)CЛ = j ^ m = а , 0 l-'

b

J лС =

0

J Sin(bаЛ)sin(bшЛ)CЛ,

(19)

0, т + к-аФ 0, т - к - а Ф 0, а + т - к Ф 0

Ь

—, а = т + к, а = к - т 4

Ь ,

--, а = т - к.

4

Таким образом, с учетом (19) правая часть выражения (18) представляется обычным однократным рядом.

Тестовый пример

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

v = -

15Q

11561b

1-550144y9 + 1859891-by^ -

- 34307072Ь2у7 + 101839360 Ь3у6 - 97143368 Ь4у5 + 3 5

12880357,5 4 3468401, 6 3 180975 , 7 2

+-Ь у--Ь у +-Ь у

2 3 2

Сходимость на первых трех шагах итерации показана на рисунке.

Изменение профиля скорости в зависимости от итерации по времени

Данные результаты получены для профиля поперечной скорости внутри канала со следующей геометрией: Ь = 3, Ь = 1, с расходом жидкости через нормальное сечение трубы Q = 1, количеством узлов сетки N = 120 вдоль нижней стенки канала и с помощью 35 членов ряда для функции тока. При этом вязкость жидкости V = 1, шаг по времени в = 0,05.

X

1

X

X

да

X

m=1

0

X

X

0

b

b

ЗО

m

0

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

Литература

1. Роуч П.Дж. Вычислительная гидродинамика. М., 1980.

616 ^

2. Флетчер К. Вычислительные методы в динамике жид-

костей. Т. 1.1. М., 1991. 552 ^

Поступила в редакцию_

3. Lesieur M., Metais O., Comte P. Large-eddy simulations of

turbulence. Cambridge, 2005. 219 p.

4. Meseguer A., Trefethen L.N. Linearized pipe flow to Reyn-

olds number 107 // J. Comput. Physics. 2002. Vol. 186. P. 178-197.

5. Сумбатян М.А., Троян Э.А. Об одном эффективном ме-

тоде вычисления звуковых полей, генерируемых турбулентными потоками жидкости // Изв. вузов. Сев.-Кавк. регион. Естеств. науки. 2009. № 4. С. 27-31.

6. Воеводин В.В., Тыртышников Е.Е. Вычислительные

процессы с теплицевыми матрицами. М., 1987. 320 с.

7. Van Loan C. Computational frameworks for the Fast Fourier

Transform. Philadelphia, 1992. 341 p.

5 сентября 2013 г.