CC BY
0_ВЕСТНИК ПЕРМСКОГО УНИВЕРСИТЕТА_
2023 • Математика. Механика. Информатика • Вып. 4(63)
Научная статья УДК 519.6: 532.5
DOI: 10.17072/1993-0550-2023-4-12-21
N-кратное расщепление явной разностной схемы для уравнения вихря в вязкой несжимаемой жидкости
Н. К. Волосова1, К. А. Волосов2, А. К. Волосова2, М. И. Карлов3, Д. Ф. Пастухов4, Ю. Ф. Пастухов4
Московский государственный технический университет (МГТУ) им. Н.Э. Баумана, Москва, Россия 2 Российский университет транспорта (МИИТ), Москва, Россия ^Московский физико-технический университет (МФТИ), Москва, Россия 4Полоцкий государственный университет, Новополоцк, Республика Беларусь
Автор, ответственный за переписку: Дмитрий Феликсович Пастухов, dmitrij.pastuhov@mail.ru
Аннотация. В работе впервые рассматривается возможность п-кратного(п=100,200) расщепления явной разностной схемы для уравнения вихря в системе уравнений гидродинамической задачи в прямоугольной каверне с вязкой несжимаемой жидкостью и с числом Рейнольдса Re=1000. Предложенный в работе алгоритм позволяет значительно увеличить максимальный временной шаг за одну итерацию общей задачи и уменьшить в десятки раз общее время расчета. Алгоритм расщепления для явной разностной схемы уравнения вихря эффективен в случае, если время, затраченное программой на цикл расщепления во много раз меньше времени решения общей задач на одну итерацию. Численно показано, что качественно решение без расщепления совпадает с решением расщепленной схемы (совпадение в пяти значащих цифрах). При этом решение задачи без расщепления не является полностью установившимся (постоянны во времени первые пять значащих цифры после 400000 итераций). Численно показано, что двухслойная и трехслойная явные разностные схемы имеют установившиеся решения с совпадением полей в 11-12 значащих знаках в каждом узле расчетной сетки (скорости, вихря, функции тока) после 21000 итераций.
Ключевые слова: численные методы; гидродинамика; метод расщепления; устойчивость; разностные схемы; вихрь; функция тока
Для цитирования: ВолосоваН.К., ВолосовК.А., ВолосоваА.К., КарловМ.И., ПастуховД.Ф., ПастуховЮ.Ф. N-кратное расщепление явной разностной схемы для уравнения вихря в вязкой несжимаемой жидкости // Вестник Пермского университета. Математика. Механика. Информатика. 2023. Вып. 4(63). С. 12-21. DOI: 10.17072/1993-0550-2023-4-12-21.
Статья поступила в редакцию 30.09.2023; одобрена после рецензирования 01.11.2023; принята к публикации 28.11.2023.
Research article
Explicit Difference Scheme A-fold Splitting For the Vortex Equation in a Viscous Incompressible Fluid
N.K. Volosova1, K.A. Volosov2, A.K. Volosova2, M.I. Karlov3, D.F. Pastuhov4, Yu.F. Pastuhov4
1Bauman Moscow State Technical University (BMSTU), Moscow, Russia 2Russian University of Transport (RUT MIIT), Moscow, Russia 3Moscow University of Physics and Technology (MIPT), Moscow, Russia 4Polotsk State University, Novopolotsk, Republic of Belarus
Corresponding author: Dmitriy F. Pastukhov, dmitrij.pastuhov@mail.ru
Эта работа О 2023 Волосова Н.К., Волосов К.А., Волосова А.К., Карлов М.И., Пастухов Д.Ф., Пастухов Ю.Ф. под лицензией СС BY 4.0. Чтобы просмотреть копию этой лицензии, посетите http://creativecommons.Org/licenses/by/4.0/
Abstract. This work is the first to consider the possibility of N-fold (n=100,200) splitting of an explicit difference scheme for the vortex equation in the system of equations of a hydrodynamic problem in a rectangular cavity with a viscous incompressible fluid and with the Reynolds number Re=1000. The algorithm proposed in the work allows us to significantly increase the maximum time step per iteration of the general problem and reduce the total calculation time by tens to hundreds of times. The splitting algorithm for the vortex equation explicit difference scheme is effective if the time spent by the program on the splitting cycle is many times less than the general problem on one iterationsolving time. It is shown numerically that the solution without splitting qualitatively coincides with the solution of the split circuit (match to five significant figures). In this case, the solution to the problem without splitting is not completely steady (the first five significant digits are constant in time after 400000 iterations). It is shown numerically that two-layer and three-layer explicit difference schemes have steady-state solutions with fields matching in 11-12 significant signs at each node of the computational grid (velocity, vortex, stream function) after 21000 iterations.
Keywords: numerical methods; hydrodynamics; splitting method, stability; difference schemes; vortex; stream function
For citation: Volosova N.K., Volosov K.A., Volosova A.K., Karlov M.I., Pastuhov D.F., Pastuhov Yu.F. Explicit Difference Scheme N-fold Splitting for the Vortex Equation in a Viscous Incompressible Fluid. Bulletin of Perm University. Mathematics. Mechanics. Computer Science. 2023;4(63):12-21. (In Russ.). DOI: 10.17072/19930550-2023-4-12-21.
The article was submitted 30.09.2023; approved after reviewing 01.11.2023; acceptedfor publication 28.11.2023.
Введение
Рассматривается гидродинамическая задача для вязкой несжимаемой жидкости в прямоугольной каверне с подвижной крышкой. Данная тестовая задача является полигоном для апробации новых вычислительных методов и алгоритмов, так как имеет в прямоугольной области две особые точки поля скорости [1], [2], [3]. При движении точки наблюдения по левой боковой стороне прямоугольника вверх в угол значение скорости равно нулю. Справа от угла на верхнем отрезке скорость скачком меняется с нуля до единицы и направлена вместе с крышкой вправо. Аналогичная ситуация происходит в верхнем правом углу. В работе предложен алгоритм «-кратного расщепления явной разностной схемы для уравнения вихря.
В работе Р.П. Федоренко описан [4, с. 137] метод расщепления уравнения в частных производных по физическим направлениям (разделение по слагаемым в общем уравнении). В методе расщепления на каждом дробном временном интервале т1 =т0/ п использовалось одно и то же уравнение вихря в общем виде, а число дробных шагов (кратность расщепления) было не «=2-3 как в работе [4], а «=100, 200.
Постановка задачи
Рассмотрим классическую гидродинамическую задачу в прямоугольной области (каверне) с системой уравнений в частных производных, начальными и краевыми условиями для физических полей [1].
Обозначим (и(х, у), у(х, у)) вектор скорости жидкой частицы (на твердой границе -боковых отрезках и на нижнем отрезке прямоугольной каверны скорость равна нулю - условие прилипания частиц жидкости). Также нормальная компонента скорости равна нулю на всей прямоугольной границе. Начало системы координат расположим в нижнем левом углу прямоугольника, направим ось у-вверх, ось х-вправо. Ширину прямоугольной каверны обозначим Ь, высоту буквой Н.
В гидродинамической задаче в закрытой каверне подвижная верхняя крышка перемещается вправо с постоянной скоростью итах характерные масштабы длины Ь, времени Ь
скорости
функции тока Lum
вихря
L
, числа Рейнольдса Re . Введем без-
размерные переменные: х - горизонтальная координата, у - вертикальная координата, —, w - безразмерные функции тока и вихря соответственно, (и, у) - вектор безразмерной скорости, I - безразмерное время по формулам
0 < X = X < 1, 0 < У = У < к = Н — = — тах = Ьитах Ь Ь Ь —т
u =-
u - V — w -, V =-w =-
u
Wmax L
t = - ,T = —,Re = ^ T w „ v
u
max
u
m ax
u
m ax
uu
max max
С безразмерными переменными и функциями запишем систему уравнений гидродинамики [1], [2]:
Wxx + W-y-y = -w(xy), 0< x = - < 1, 0<y <К
w = vx - uy,
u = wy ;v = -wx
w't + u ■ Wx + v • w
— (wxx + Wyy 10 < t =
Re V p T
(1)
uixn , y m ) = -U0 ixn ) ~Jr I Slnl
ym LJ
2k
xn = nK y m = mh2 , h1 = — , h2 = —
(2)
n = 0, n, m = 0, n2, n = П = 100
где профиль (3) горизонтальной компоненты скорости на верхнем отрезке каверны
(у = к = 1) имел вид симметричной трапеции [2] без особых точек поля скорости:
l(x, k) = щ (x) =
x „ n 1
— ,0< x<t,t = -^ = — 20
n,
1,t< x < 1 -T,
(3)
1 - x
1 -t< x< 1,
^ = 0, V = 0, щ = 0, щ = 1
1г 1г 1г, 1г\г,
Здесь Г1 - объединение боковых сторон и нижнего отрезка, Г\Г1 - верхний отрезок прямоугольника Г. Первым в системе (1) следует уравнение Пуассона для неизвестной функции тока и известной функции вихря. Двумерное уравнение Пуассона на прямоугольнике решается в матричном виде за конечное число арифметических действий с четвертым или шестым порядком погрешности [5, 6]. Далее мы опустим черту сверху над безразмерными функциями, временем и координатами.
Вторая строка системы (1) - функция вихря вычисляется через координатные производные поля скорости. Третья строка - компоненты скорости - вычисляются как частные производные от функции тока. Четвертая строка - уравнение динамики вихря, которое в системе уравнений (1) единственно явно зависит от времени. Слева стоит полная (конвективная) производная по времени. На границе прямоугольника отсутствует вертикальная компонента скорости, а горизонтальная равна нулю на нижнем отрезке и боковых сторонах, а на верхнем отрезке она равна единице.
Кроме двух упомянутых особых точек поля скорости для тестирования алгоритма использовалось сильно нестационарное и завихренное начальное поле скоростей. Оно задавалось нами [2] на равномерной прямоугольной сетке по формуле (2):
Аналитически метод «-кратного расщепления уравнения вихря внутри одного временного интервала можно записать в виде
w
,k+((i+1)/n) _ w—+(i/n)
к к+(i/n) , k k+(i/n) _
T0/ n
■ + un ■ wyK""' + v" ■ W
y
5 (
wkx+(i/n) + wky+(i/n) ), i = 0, n -
), i = 0, n -1 . (4)
Система рекуррентных уравнений (4) для вихря с замороженным полем скорости (щк (х, у), V1' (х, у)) состоит из п дробных шагов
г = 0, п—1, верхний индекс / - указывает дробный слой времени в (4), индекс к- номер целого слоя времени в уравнении вихря в системе (1). Поля скорости, функции тока меняются последовательно, согласно системе (1), в которой поле вихря имеет уже не дробные индексы
, „к+ г / п п «
^ , а целые ^ . В цикле с системой уравнений (4) изменяется только поле вихря м>к+г / п, г = 0, п — 1. Поле скорости скачком изменяется, когда переменная внешнего цикла увеличивается на единицу от к до к+1в системе уравнений (1).
Идея решения системы уравнений (4) заключается в уменьшении накопления ошибки округления и времени вычислений. Дифференциальные операторы по координате в (4) аппроксимированы, как и все уравнения системы (1), а также граничные условия с точностью о(г4), а по времени с точностью О(г).
Предположим, что система уравнений (1) спектрально устойчива с максимальным временным шагом равным т0, а система уравнений (4)
с максимальным временным шагом т0 / п .
Таким образом, за время т0 / п, решая п раз уравнение (4), мы получаем уже скачок по времени т0 (в п раз больший, чем последовательное решение системы уравнений (1)) и уменьшаем ошибку округления, так как внутри цикла системы (4) не решаем другие уравнения системы (1).
T
n
n
2
Второе важное преимущество схемы п-кратного расщепления уравнения вихря (4) заключается в достижении установившихся конечных полей скорости, вихря и функции тока (15 значащих цифр в каждой точке поля не меняются во времени после 21000 итераций).
Уравнение (4) линейно относительно координатных производных w1X, w1y, wXX, w1yy . Аппроксимируя каждую производную с порядком о(к4), получим невязку уравнения (4) по координате с той же точностью. Для первой производной на симметричном шаблоне методом неопределенных коэффициентов [5], [6] имеем
Wx(0) =1 (21" w_l- )] + о(к4) . (5)
1 ( 5 4
к \ 3
Получим аналогичную формулу для второй производной на симметричном шаблоне [5], [6]:
Wxx(0) = ,2 1 „ "0
1 (-1"0 + 4+ "-1)-^("2 + "-2)| + 0(к') . (6)
™ 1 (5 5 1 7 1 1 I /
х(0) = —| 7"-1 --"0--"1 +-"2--"3 + — "4 |+ 0(к
к2 I 2 0 3
Я + - + ^)-- ^ + "-2 )| = ^
12
к4"хб)(#2), (10)
#1,#2 е(- 2к,2к).
Таким образом, для уравнения вихря в системе уравнений (1) во внутренних узлах с учетом формул (9), (10) можно найти невязку (11) на временном слое к с шагом т0:
"к+2 ик.("х (0) - ± к4 |+
+ ^'("у(0)к4| = ^("х*х + "Уу -
к-("хб)(#з) + "уб)(#4))|,#1,#2,#3,#4 е(-2к,2к)<
К = 2 к4(и^х5)(^1) + ^(Ф
7 4
+ 9(^ё ("хб)(^з) + о(20 + к4). (11)
Порядок аппроксимации невязки
Формулы (5), (6), можно использовать на внутренних узлах, удаленных от сторон прямоугольника не менее чем на два координатных шага сетки.
На первом прямоугольном контуре с узлами, удаленными на один шаг от сторон прямоугольника, получим первую производную (методом м. н. к.) на шаблоне со смещенным центром [5], [6] (индекс -1 имеет узлы сетки на границе прямоугольника):
/т 1 ( 1 5 3 1 1 I п1, Л " (0) = —I-—- —^ +—Щ-—^ +— Щ 1 + 0(к ) . (7) х к{ 4 -1 6 0 2 1 2 2 12 3) у '
Аналогично, для второй производной получим формулу (8) на шаблоне со смещенным центром (узел с нулевым индексом расположен на первом прямоугольном контуре, а узел с индексом -1 расположен на границе прямоугольника). Назовем прямоугольник и узлы на нем нулевым прямоугольным контуром [5], [6], [1], [2].
Я = о(т0+ к4)
сохраняется и на первом прямо-
Разложив формулы (5), (6) в ряд Тейлора на центральном шаблоне, получим невязку (центр шаблона с координатой х=0):
угольном контуре согласно формулам (7) и (8).
Рассмотрим спектральную устойчивость разностного уравнения вихря в системе (1) в пределе при быстром движении жидкости (число Рейнольдса Ке = 1000 ) и при медленном вязком течении, например, движении крови в капиллярах Яе ^ 0 )[7]. Приведем определение спектральной устойчивости из работы [6, с. 125].
Определение [6], [4]. Если при заданном законе стремлении 2, к ^ 0 существует постоянная 0 < с < да , такая, что для всех р, — е [0,2;?] справедливо неравенство \Ар—)< ехр(ст), то
спектральный признак выполнен.
В первом случае в (1) при Яе ^-да пренебрежем диффузией в уравнении вихря. Тогда с центральным шаблоном на внутренних узлах сетки запишем разностное уравнение переноса в системе уравнений (1):
к+1 к
" -" к к к к Г, 20
-+ ик ' + ук ' = 0, ^ = —
х у 1 к
2
+ ^ (ик 1"к + ук I"к ) = 0, (12)
т,п т,п 1 \ т,п х т,п т,п у т,п/ > х*-*"/
1 (2 (
к 13
("1 -"-1)-1 ("2-"-2)]="х(0)-1к^а (9) I=2-"т,п-1)-1 (
" , о - " о
1 т,п+2 т,п-2
I
1 k 2 i k k \ 1 i k k
l wk = — lw, — w, I--lw„ — wk „
y m,n 3 \ m+1,n m—1,n / ^^ \ m+2,n m—2,n
).
Подставляя в уравнение (12) функцию возмущения [4] м>ктп =^к (ф,^)ехр ('(пф+т^)\
г = 4-1, [0,2ж] получим спектральное уравнение (13):
./2max = max^ | 4 — Vl^
xxe[0,1] 3
/m ax
max( /1
ma^ 2max ) = /1 m
Тогда из (13) получим абсолютное значение:
H = J1+U.JС)+vinfУ)) ^ 1+¿M/mi и1+2/maZ?:
^ =1+Ц<п I —1sin(2^)|+vm,n I jsin(^)—1sin(2^) II.
Рассмотрим задачу на максимум:
/ (ф) =
/ (ф) =
4 1
-sm(iP) — - sm(2p) 3 6
|4 —
cos (Ф ) = ■!
4 ±д/1 — sin2(p)
5
<- .
3
к задаче
мума:
/(x) = xI 4 + x2 max ^ / (x) = 0
/ (x)=3
4 + -n/T— x2 + x
^ — 2x ^
241—:
= 0 »
//
WI—x2 + 1 — 2x2
= 0 ^ 4x4 +12 x2 —15 = 0 ^
= _3W2i = — 3 + ^ <x= JV6 — 3 < 1
/-=as! (41=3 (4 1
x=J V6—3
V6—3
4 + J5 —46
i 1,37222 ;
b) во втором случае
/(x) = x (4 — 4\—x£ max ^ / (x) = 0
/(x) =
1 — x 2 —:
^ — 2x ^
241—x
= 0 »
W1 — x2 —
x2 —1 + 2x2
41—x
= 0 ^ 4x +12x2 —15 = 0 .
Но x2 =-3 + V6.W1 - x2 = 1 - 2 x2 = 4 - 246 < 0
2
и локальных точек экстремума нет x e (0,1):
^ max, фе [0,2^], (14)
Обозначив x = |sin(p)|,x е[0,1], перейдем
/(x) = x(4±41 — x2 j^max, xе[0,1] : a) запишем необходимое условие экстре-
= 1 + 2/^А = 1 + 2/Цу = 1 + Т ^2/тт-±j .
Сравнивая последнее выражение с определением спектрального признака устойчивости получим [5], [6] постоянную c
Щ < ex?(сто) и 1 + сто ° с = 2/1x^1 = const° Т = const.
При медленном вязком течении [7] Re ^ 0 можно пренебречь конвективным переносом, но сохранить слагаемое с диффузией вихря.
Разностная схема для уравнения вихря в (1) на внутренних узлах сетки примет вид
wk+1 — wk 1
= ^ ^ + wyy) z2 =
т Rev xx yyp 2 Reh2'
wk+1 = wk + 7 f— 5wk +—(wk ,+ wk ,)—L(wk + wk )+
" " т a2 I m,n 3 Wn+1 m,n—1jy^ V m,n+2 T ''m,n—2 Г
' i k k V 1 / k k
+ 3K+1,n + Wm—1,n )^(Wm+2,n + Wm —2,n
^ )-1 W
)|. (15)
Подставим в разностную схему (15) функцию возмущения вида [4]:
wkmn = I(р,у)exp(i(np + mу/)),i = V-I,р,у e [0,2^]
получим
8
спектральное 1
уравнение.
Х = 1 + z21— 5+-(cos(i>)+cos(//))—(cos(2ф)+cos(2^)) |. (16)
3 6
Введем вспомогательные переменные: x = cos(p), y = cos(y), x, y e [-1,1],cos(2p) = -1 + 2cos2(p) . Перепишем уравнение (16) в виде
8
к(х, у) = 1 + 72 5 + -(х + у) — - (х2 + у2 — 1) . (17)
Исследуем на экстремум выражение, стоящее в круглых скобках формулы (17):
я(х,у)=—5+3( х+у)—3( х 2+у 2—1) ^ ^ 8 2
Ях (х, у) = 0 «3 — 3 х = 0 х = 4, 8 2
Яу (х, у) = 0 —- у = 0 у = 4 .
2
x
2
3
3
т
0
2
2
3
2
2
Но точка локального экстремума (4,4) й [-1,1] х [-1,1] функции g(х, у). Из симметрии функционала также следует равенство координат в точке экстремума. Поэтому достаточно вычислить функцию g(х, у) в точках (-1, -1) и (1, 1).
8 1 32
g (-1,-1) = -5 + -(-1 -1) - ^(1 +1 -1) = - 32 3 3 3
8 1
g(1,1) = -5 + -(1 +1)-^(1 +1 -1) = 0 .
точно
32
Для спектральной устойчивости доста-чтобы \Л(х, y)| < 1Vx, y е [-1,1], то
1 - — Z2 = 1 + Z2g(-1,-1) = -1 < Л(X, y) < 1 = 1 + Z2g(1,1),
32 6 3
— z, = 2 « z, = — = —, 3 2 2 32 16
0Д 16
(18)
0,— 16
Или z, =
2 Re h2
3 3
< — orn< —Re h2 (*). 16 0 16
[1], [2] в системе (19). На границе прямоугольника в (19) значения функции вихря связаны с граничными значениями функции тока и компонентами скорости [1], [2].
Таким образом, если уравнение переноса
спектрально устойчиво для любого закона: f
= const, то разностное уравнение (15) с диффузионным членом спектрально устойчиво только на отрезке z2 е
1„к+1 е../+1 +1 +1
k+1 'У m,0 - Щ m,1 + У m,2 V m,0 W m,0 =-:--3-,
1 k+1 _o k+1 , k+1 k+1
I у m,n Оу m,n -1 + у m,n-2 ^ ^ V m,nl
k+1 'у m,n
W m,n.
2h
m = 1, n, -1
k+1 k+1 k+1 k+1 k+1 >W 0,n - 1,n +У 2,n ,u 0,n
W 0,n =-;--+ 3-
2h
k+1 k+1 k+1 k+1
k+1 n2,n - оу n2-1,n + у n2-2,n ,u n2,n
2h2
--3-
n = 1, n -1,
k+1 k+1 k+1 k+1 k+1 W 0,1 + W 1,0 k+1 W n2 ,n1 -1 + W n2 -1,n W 0,0 =-, W n, ,n, =-
'(19)
k+1 k+1 k+1 k+1 k+1 W 0,n1 -1 + W 1,n1 k+1 W n2,1 + W n2 -1,0
W 0,n1 = -
2
-, W n2,0 =-
2
В формуле (*) шаг 20 максимальный (соответствует верхней границе спектральной устойчивости) в задаче (1). На практике реальный шаг времени разностной схемы
21 = 20 / п много меньше 20 из-за решения других уравнений системы в (1) и 2 особых точек поля скорости.
Обозначим реальный шаг
2 3
2, = —,Н2, <— Яек2 = 2„ 1 п 1 16 0
Видно, что при достаточно большой кратности расщепления п явной разностной схемы (4) с шагом Т1 всегда можно выполнить
3
условие п^ < — Яе к2 .
Отметим, что поиск новых значений вихря в цикле (4) проводится только для внутренних узлов, то есть для узлов первого прямоугольного контура и узлов, расположенных дальше от границы прямоугольника.
Граничные значения вихря со вторым порядком погрешности приведены в работах
Оценим время решения задачи (1) в трех случаях.
В первом случае (алгоритм А) система уравнений (1) решается последовательно, алгоритм расщепления (« = 1) в системе уравнений (4) не используется, временной шаг (т1=ть/«=ть).
Во втором случае (алгоритм В) уравнение вихря (4) решается на двухслойной схеме с кратностью расщепления « = 100. В алгоритме С « = 200 и трехслойная схема по времени. Программы алгоритмов А, В, С запускались одновременно и заметного временного различия для достижения выбранного числа итераций (внешней переменной цикла к) не обнаружено (из-за простоты явной разностной схемы (4)). Но если для наступления стационарного режима в алгоритме В было достаточно 1,5 часа работы программы, то для алгоритма А стационарный режим не был достигнут после двух суток работы программы. Тогда (та/тв) ~ (3/16)/(1/300) = 56 (рис. 2).
На рис. 1 показано поле линий тока для решения задачи с кратностью « =1 (отсутствие расщепления). Видно, что спустя 6000 шагов времени поле линий тока находится далеко от стационарного состояния. Центральный вихрь первого порядка обменивается моментом импульса с вихрями второго порядка, расположенных в нижних правом и левом углах прямоугольной каверны.
h
h
z, е
2
2
2
г.
В работе [3] авторы А.А. Фомин и Л.Н. Фомина также получили численно стационарное решение задачи (1) и вихри трех порядков.
В работе Фоминых большее разрешение сетки 2500x2500 (Яе>10000), число ячеек в нашей сетке 100x100 (и Яе=1000 - в 10 раз меньше). Хотя нестационарный вихрь третьего порядка обнаружен нами на рис. 1.
Рис. 1. Поле линий тока для решения задачи (1) без расщепления («=1) Яе=1000.
Шаг времени
т = -
h2Re 300
n = n = 100,
A
- . i- —
■
-
-
/0)
I ь
h = 1/100 в момент t = 6000 т
B
На рис. 2 показаны поля линий тока в предельных стационарных состояниях.
Рис. 2(А) соответствует решению задачи (1) без расщепления.
Рис. 2(В) соответствует решению задачи (4) с расщеплением «=100 с двухслойной аппроксимацией производной по времени
и первым порядком погреш-
wk+(i+1)/n _ wk+'/n т0/ n
ности О(т1), стационарное состояние достигается спустя N=21000 итераций, после чего max и min значения функции тока сохраняются с 15 значащими цифрами.
Третий рисунок (рис. 2(С)) соответствует стационарному решению с трехслойной
wk+(i+1)/n _ wk+(i-1)/n
временной производной
2т0 / n
ОТ2)
и с
со вторым порядком погрешности кратностью расщепления «=200.
Поле тока становится стационарным спустя N=43000 итераций.
ч ..... — "л
/ \
a / /
/ / V ; i
/ i / /
i \ J /J
N 4 \ "").....,/ f\
_ \ ■■■■ / \
- ......
ч ч Щ
Ж
i i
C
Рис. 2. Установившееся поле линий тока для задачи (1), число Рейнольдса Яе=1000, « - кратность расщепления, 2 = г0 / п .
Рис. 2(А). 2=2= ,* = 400000 2,п = 1;
300
3
Рис. 2(B). т0 = — h1 Re, t = 21000 т, n = 100;
16 3
Рис. 2(C). т0 = — h2 Re, t = 43000 т0, n = 200 16
и
Общими показателями решения (1) являются максимальное и минимальное значение поля функции тока, поскольку все остальные поля можно получить через дифференциальные операторы от поля функции тока (как видно из задачи (1)).
Рассмотрим max и min значения функции тока в предельном (стационарном) состоянии решения задачи (1) для трех случаев:
Таблица 1. Сравнение значащих цифр в стационарном поле функции тока в задаче (1) для различных алгоритмов
Алгоритм ^max ^mm
A) n=1 1,71588 10"3 -0,1183079
B)n=100 1,715838733007 10-3 -0,1183038468
C)n=200 1,715838733009 10-3 -0,1183038468
Примечание:
А)п=1 (2-х слойная схема), N=400000; B)п=100 (2-х слойная схема), N=21000; ^п=200 (3-слой-ная схема), N=43000
Таблица 2. Сравнение значащих цифр в нескольких точках стационарного поля горизонтальной компоненты скорости в задаче (1) для трех алгоритмов
Ал-
го- A:n=1 B:n=100 C:n=200
рит
м
Г1 3,25838 3,258397414826 3,258397414825
9 10-2 9 10-2 7 10-2
Г2 -9,772 -9,77226561 -9,77226561
15 10-2 55239 10-2 55180 10-2
Г3 -1,732 -1,73241760 -1,73241760
42 10"2 2921899 10"2 2922225 10"2
Анализ ушп, Ушах показывает, что два разных алгоритма: В - двухслойная схема по времени с расщеплением п=100; С - трехслойная временная схема с п=200 дают одинаковые предельные решения, в которых стационарные Ушп, Ушах совпадают в 11-12 значащих цифрах (табл. 1). То же можно сказать о горизонтальной компоненте и стационарного поля скорости в нескольких, но одних и тех же узлах сетки - совпадение в 11-12 значащих цифрах (табл. 2). Поэтому рис. 2(В) и рис. 2(С) неразличимы графически.
Решение без расщепления (алгоритм А) и решения с расщеплением (алгоритмы В, С) совпадают в 5 значащих цифрах.
Основные полученные результаты:
1) Впервые предложен алгоритм (4) п-кратного расщепления (п=100, 200), явной разностной схемы для уравнения вихря в десятки раз ускоряющий время решения гидродинамической задачи.
2) В предельных случаях с большими и малыми числами Рейнольдса показано, что для спектральной устойчивости разностной схемы уравнения вихря следует выбирать закон [5], [6] стремления к нулю т0, h.
т0 = о(и2\г0, И ^ 0,птх < — И2 Яе = г0 .
3) Решения с расщеплением в системе
3
уравнений (1) и шагом времени т0 =—h Re
16
соответствуют верхней границе спектральной устойчивости. Спектральная устойчивость с т0 подтверждена численно.
4) Предложена двухслойная временная схема для уравнения вихря с аппроксимацией û(T + h4 ).
Список источников
1. Salih A. Streamfunction - Vorticity Formulation // Department of Aerospace Engineering Indian Institute of Space Science and Technology, Thiruvananthapuram-Mach 2013. p.10.
2. Сборник статей по гидродинамике / Н.К. Волосова, К.А. Волосов, А.К. Волосова [и др.]. 2-е изд. М.: МИИТ, 2023. 231 с. EDN: UDVEDI.
3. Фомин А.А., Фомина Л.Н. Численное моделирование течения жидкости в плоской каверне при больших числах Рейнольдса // Вычислительная механика сплошных сред. 2014. Т. 7, № 4. С. 363-377.
4. Федоренко Р.П. Введение в вычислительную физику: учеб. пособие для вузов / Р.П. Федоренко; Федоренко Р.П.; под ред. и с доп. А.И. Лобанова. 2-е изд., испр. и доп. Долгопрудный (Моск. обл.): Изд. дом Интеллект, 2008. 503 с. (Физтеховский учебник). ISBN 978-5-91559-011-2. EDN: QJUAEP.
5. Бахвалов, Н. С. Численные методы: учеб. пособие для студ. физ.-мат. специальностей вузов / Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков; Мос. гос. ун-т им. М.В. Ломоносова. 7-е изд. М.: Бином. Лаб. знаний, 2011. 636 с. (Классический университетский учебник). ISBN 978-5-9963-0449-3. EDN: QJXMXL
6. Бахвалов Н.С., Лапин А.В., Чижонков Е.В. Численные методы в задачах и упражнениях. М.: БИНОМ, 2010. 240 с.
7. Волосов К.А., Вдовина Е.К., Пугина Л.В. Моделирование "пульсирующих" режимов динамики свертывания крови. Математическое моделирование. 2014. Т 26, № 12. С.14-32.
References
1. Salih A. Streamfunction - Vorticity Formulation // Department of Aerospace Engineering Indian Institute of Space Science and Technology, Thiruvananthapuram-Mach 2013. P. 10.
2. Sbornik statej po gidrodinamike / N.K. Voloso-va, K.A. Volosov, A.K. Volosova [i dr.]. 2-e izd. M.: MIIT; 2023. 231 p. EDN: UDVEDI. (In Russ.).
3. Fomin A.A., Fomina L.N. CHislennoe modeli-rovanie techeniya zhidkosti v ploskoj kaverne pri bol'shih chislah Rejnol'dsa. Vychis-litel'naya mekhanika sploshnyh sred. 2014;(7(4):363-377. (In Russ.).
4. Fedorenko R.P. Vvedenie v vychislitel'nuyu fiziku: ucheb. posobie dlya vuzov / R.P. Fedorenko; Fedorenko R.P.; pod red. i s dop. A.I. Lobanova. 2-e izd., ispr. i dop. Dolgo-prudnyj (Mosk. obl.): Izd. dom Intellekt; 2008. 503 p. (Fiztekhovskij uchebnik). ISBN 978-591559-011-2. EDN: QJUAEP. (In Russ.).
5. Bahvalov N.S. CHislennye metody: ucheb. posobie dlya stud. fiz.-mat. special'nostej vuzov. Mos. gos. un-t im. M.V. Lomono-sova. 7-e izd. Moskva: Binom. Lab. Znanij; 2011. 636 p. (Klassicheskij universitetskij uchebnik). ISBN 978-5-9963-0449-3. EDN: QJXMXL. (In Russ.).
6. Bahvalov N.S., Lapin A.V., CHizhonkov E.V. CHislennye metody v zadachah i uprazh-neniyah. M.: BINOM; 2010. 240 p. (In Russ.).
7. Volosov K.A., Vdovina E.K., Pugina L.V. Modelirovanie "pul'sirubshchih" rezhimov dinamiki svyortyvaniya krovi. Matemati-cheskoe modelirovanie. 2014;(26(12)):14-32. (In Russ.).
Информация об авторах:
Наталья Константиновна Волосова - аспирант МГТУ им. Н. Э. Баумана (105005, Россия, г. Москва, 2-я Бауманская ул., д. 5, стр. 1), navalosova@yandex.ru;
Константин Александрович Волосов - доктор физико-математических наук, профессор кафедры прикладной математики Российского университета транспорта (127994, ГСП-4, Россия, г. Москва, ул. Образцова, д. 9, стр. 9), konstantinvolosov@yandex.ru, АиШогГО 128228; Александра Константиновна Волосова - кандидат физико-математических наук, начальник аналитического отдела ООО "Трамплин" Российского университета транспорта (127994, ГСП-4, Россия, г. Москва, ул. Образцова, д. 9, стр. 9), alya01@yandex.ru, АиШогГО 607500;
Михаил Иванович Карлов - кандидат физико-математических наук, доцент кафедры высшей математики Московского физико-технического университета (МФТИ) (141701, Россия, Московская область, г. Долгопрудный, Институтский пер., 9.), karlov@shade.msu.ru, А^^гГО 14680;
Дмитрий Феликсович Пастухов - кандидат физико-математических наук, доцент кафедры технологий программирования Полоцкого государственного университета (211440, Республика Беларусь, Витебская обл., г. Новополоцк, ул. Блохина, 29), dmitrij.pastuhov@mail.ru, Аи-^гГО 405101;
Юрий Феликсович Пастухов - кандидат физико-математических наук, доцент кафедры технологий программирования Полоцкого государственного университета (211440, Республика Беларусь, Витебская обл., г. Новополоцк, ул. Блохина, 29), pulsar1900@mail.ru, А^^гГО 405109.
Information about the authors:
Natalya K. Volosova - Post-graduate Student of Bauman Moscow State Technical University (5-1, 2nd Baumanskaya St., Moscow, Russia, 105005), navalosova@yandex.ru;
Konstantin A. Volosov - Doctor of Physical and Mathematical Sciences, Professor of the Department of Applied Mathematics of the Russian University of Transport (9-9, Obraztsova St., Moscow, GSP-4, Russia, 127994), konstantinvolosov@yandex.ru, AuthorlD 128228;
Aleksandra K. Volosova - Candidate of Physical and Mathematical Sciences, Chief Analytical Department "Tramplin" LLC, Russian University of Transport (9-9, Obraztsova St., Moscow, GSP-4, Russia, 127994), alya01@yandex.ru, AuthorlD 607500;
Mikhail I. Karlov - Candidate of Physical and Mathematical Sciences, Associate Professor of the Department of Higher Mathematics, Moscow University of Physics and Technology (9, Institutskiy per., Dolgoprudny, Moscow region, Russia, 141701), karlov@shade.msu.ru, AuthorlD 14680;
Dmitriy F. Pastukhov - Candidate of Physical and Mathematical Sciences, Associate Professor of Polotsk State University (29, Blokhin St., Novopolotsk, Vitebsk Region, Republic of Belarus, 211440), dmitrij.pastuhov@mail.ru, AuthorlD 405101;
Yuriy F. Pastukhov - Candidate of Physical and Mathematical Sciences, Associate Professor of Polotsk State University (29, Blokhin St., Novopolotsk, Vitebsk Region, Republic of Belarus, 211440), pulsar1900@mail.ru, AuthorlD 405109.
