CC BY
41
УДК 532.5:519.63
С. И. Мартыненко
СОВЕРШЕНСТВОВАНИЕ ВЫЧИСЛИТЕЛЬНЫХ АЛГОРИТМОВ ДЛЯ РЕШЕНИЯ УРАВНЕНИЙ НАВЬЕ-СТОКСА НА СТРУКТУРИРОВАННЫХ СЕТКАХ
Предложен новый подход к вычислению давления при решении полных уравнений Навье-Стокса в переменных скорость-давление на структурированных сетках. В его основу положено использование интегральных форм уравнения неразрывности (уравнений постоянства массового расхода) и декомпозиции давления, исходя из которых формулируется вспомогательная задача. Ее решение близко к решению полных уравнений Навье-Стокса при меньшем объеме вычислений. Приведены методика построения вспомогательной задачи, описание вычислительного алгоритма и результаты решения тестовых задач.
Двумерные стационарные уравнения Навье-Стокса в переменных скорость-давление в координатах X, У для несжимаемой среды могут быть записаны в следующей безразмерной форме:
— уравнение неразрывности
ди дь ...
дх ду
— уравнение движения по координате X
д (и2) д (ьи) др 1 (д2и д2и\
дх dy dx Re \ дх2 ду — уравнение движения по координате Y
д (uv) д (v2) dp 1 (d2v д2у\
( ) + ^ = -дР + тН ^ + ^ . (3)
дх дуу дуу Re \ дх2 ду2 /
Численное решение уравнений Навье-Стокса при больших значениях числа Рейнольдса связано с многочисленными трудностями, наиболее значимые из которых нелинейность уравнений движения, которая вызывает необходимость принятия специальных мер по обеспечению сходимости итераций по нелинейности; сингулярная возмущен-ность уравнений движения, связанная с наличием малого параметра (1/Re) при старших производных. Как правило, решения уравнений Навье-Стокса обладают особенностями типа пограничного слоя, для разрешения которых необходимо сгущение сетки. Построение вычислительных сеток в областях со сложной геометрией, адаптированных к особенностям численного решения, — трудоемкая задача, которая
еще не получила однозначного решения; уравнения Навье-Стокса, в конечном итоге, сводятся к системе линейных алгебраических уравнений высокого порядка, однако выбор наиболее эффективного метода ее решения не является очевидным. Кроме того, нелинейность уравнений Навье-Стокса крайне затрудняет их теоретическое исследование и анализ вычислительных алгоритмов. Поэтому методы численного решения уравнений Навье-Стокса носят эвристический характер и приближенные представления об их свойствах могут быть получены только при решении тестовых задач.
Структура уравнений Навье-Стокса осложняет построение эффективных вычислительных алгоритмов при любых значениях числа Рей-нольдса. Построим сетку и осуществим на ней аппроксимацию уравнений (1)-(3). Результирующую линеаризованную систему разностных уравнений можно записать в матричной форме
ai3\ fu\ /ЬЛ
7 J(pJ = Ы ■ (4)
Очевидно, что полученная система линейных алгебраических уравнений (4) не может быть решена итерационными методами типа Яко-би или Зейделя из-за наличия нулевого блока на главной диагонали (a33 = 0), обусловленного отсутствием слагаемого, связанного с давлением, в уравнении неразрывности (1).
Первоначально для решения системы (4) применяли сегрегированные (раздельные) алгоритмы, в которых сначала вычислялись компоненты скорости из решения системы
'яц ai2\ /u\ = /bi - ai3p &21 0,22/ Vv) \Ь2 - Я23Р; а затем корректировалось давление из условия, что компоненты скорости будут удовлетворять разностному аналогу уравнения неразрывности a31u + a32v = 0. В настоящее время наибольшее распространение получил метод SIMPLE, в котором для поправок к давлению (5p) предложено уравнение
д25,0 д25,0 ди dv
—- +--- =--1--■
дх2 ду2 дх ду
Сходимость метода SIMPLE достигается посредством применения последовательной нижней релаксации, т.е. давление пересчитывается по соотношению p := p + ap5p (равенство в смысле присвоения), где ap есть параметр нижней релаксации [1].
К неустранимым недостаткам метода SIMPLE и его модификаций следует отнести дополнительные граничные условия для давления,
p := p - ßp\— + —
которые отсутствуют в постановке задачи. Поэтому в приложениях часто используется метод искусственной сжимаемости [2], в котором для вычисления давления предложено уравнение
/ди дь^ дх ду
где вР — параметр нижней релаксации. Метод искусственной сжимаемости исключительно просто реализуется, поскольку для отыскания давления не требуется решения дополнительных краевых задач.
Однако сегрегированное (раздельное) отыскание компонент скорости и давления затрудняет сходимость итераций по нелинейности при решении уравнений движения. Поэтому сегрегированные алгоритмы применяются, как правило, к решению нестационарных уравнений Навье-Стокса, а стационарные решения отыскиваются посредством счета на установление.
В первых совместных алгоритмах для решения системы (4) применяли прямые методы [3]. Совместные алгоритмы выглядят предпочтительнее сегрегированных, поскольку полнее сохраняется взаимосвязь между компонентами скорости и давлением в итерационном процессе. Однако объем вычислений при решении системы (4) прямыми методами оказывается слишком большим, особенно если для аппроксимации уравнений Навье-Стокса используются мелкие сетки. Позднее были предложены другие методы, которые еще не получили широкого распространения [4].
Сегрегированные алгоритмы проще совместных, что особенно заметно при решении уравнений Навье-Стокса для сжимаемых сред. В последнем случае не только конвективные, но и диффузионные члены уравнений движения — нелинейные, и разностный аналог уравнений Навье-Стокса трудно записать в виде системы (4). Однако скорость сходимости совместных алгоритмов выше, чем сегрегированных.
В настоящей статье предлагается вычислительный алгоритм, который занимает промежуточное место между сегрегированными и совместными алгоритмами и объединяет их преимущества. Сначала рассмотрим течение между параллельными пластинами, которое описывается уравнениями Навье-Стокса в приближении пограничного слоя:
— уравнение неразрывности
ди ^ дь ^ дх ду
— уравнение движения по координате X
д(и2) ^ д(ьи) йр ^ 1 д2и_ дх ду йх Яе ду2'
— уравнение постоянства массового расхода 1 1
J и(х,у) йу = У и(0,у) йу . (7)
0 0
Уравнения (5)-(7) получены из (1)-(3) в предположении, что изменение давления поперек потока равно нулю (р'у = 0) и |иуу | ^ |иХх Принятые допущения позволяют решать уравнения (5)-(7) маршевым методом по направлению X (по потоку). Положим, что в (г — 1)-м сечении, перпендикулярном направлению течения, известно распределение скорости. Тогда линеаризованный разностный аналог уравнения движения (6) в г-м сечении может быть записан в виде
А, Щ]-1 + В, щ, + С, 1 = Б, + Рг 1 , (8)
^х
где кх есть шаг сетки по направлению х. Численное решение уравнений (5)-(7) проводится в следующем порядке: 1) решение разностного уравнения (8) методом прогонки по мере сходимости итераций Ньютона по нелинейности в конвективном члене; 2) получение компоненты скорости V из разностного аналога уравнения неразрывности (5); 3) проверка условия постоянства массового расхода (7), коррекция давления рг и возврат к п. 1, если необходимо; 4) переход к (г + 1)-му сечению.
Коррекция давления основана на численном решении уравнения
1 1
Р(Рг) = ! и(хг/у,'Рг) йу — J и(0,у) йу , 00 которое осуществляется методом секущих [5, 6].
К очевидным достоинствам приведенного выше алгоритма следует отнести: высокую вычислительную эффективность определения давления (давление вблизи г-го сечения изменяется почти линейно, а метод секущих является точным для линейных задач, поэтому требуется всего несколько итераций метода секущих для подбора значения давления рг таким образом, что погрешность разностного аналога уравнения постоянства массового расхода (7) станет сопоставимой с погрешностями округления); согласованное вычисление компонент скорости и давления.
Предлагаемый алгоритм решения полных уравнений Навье-Стокса основан на неформальном разделении расчета компонент скорости и давления с привлечением методов, используемых ранее для решения уравнений Навье-Стокса в приближении пограничного слоя (5)-(7).
Вспомогательная задача. Формулировку вспомогательной задачи проиллюстрируем на примере стационарного течения несжимаемой среды в каверне Н х Н с крышкой, движущейся со скоростью и,ю. Положим, что плотность р и вязкость ^ среды постоянные. В данном случае число Рейнольдса в уравнениях Навье-Стокса (1)-(3) определяется так:
^ ришн
Ке =-.
На стенках и крышке каверны граничные условия для компонент скорости заданы в виде
и(х, 0) = п(0, у) = и(1, у) = 0; и(х, 1) = 1; у(х, 0) = г>(0, у) = у(1, у) = у(х, 1) = 0.
Теперь сформулируем вспомогательную задачу для расчета течения в каверне, которая по вычислительным затратам была бы сравнима с решением уравнений Навье-Стокса в приближении пограничного слоя, но решение которой близко к решению исходной задачи (1)-(3). Для этого проинтегрируем уравнение неразрывности (1) по контрольным объемам У1 и У2, показанным на рис. 1, и получим два уравнения постоянства массового расхода
1
J и(х,У) ¿у = 0; (9)
о
1
/ф,у) ¿х = о. овд
о
Далее представим давление в виде
р(х,у)= рх(х) + ру (у) + рху(х,у), (11)
где верхние индексы х, у и ху показывают зависимость того или иного слагаемого от соответствующих координат. Во вспомогательной зада-
че вместо одного уравнения неразрывности (1) используются уравнения постоянства массового расхода (9) и (10). Поскольку из (11) непосредственно следует
др (1рх дрху др дх (х дх ду
dpy dpxy dy dy
то вспомогательная задача принимает следующий вид:
— уравнение движения по координате X и уравнение постоянства массового расхода (9)
д (и2) д (vu)
дх
дУ
dpx dx
др
xy
дх
1
+ Re
д2и д2и дх2 ду2
(12)
и(х, у) dy = 0,
v о
— уравнение движения по координате У и уравнение постоянства массового расхода (10)
д (uv) д(v2)
дх
дУ
dpy dy
др
xy
1 i д2,о д2,о Re \ дх2 дy2
(13)
v(x, y) dx = 0,
где квадратные скобки означают, что производные (рху)'х и (рху)'у являются фиксированными (т.е. вычисляются по значениям рху, взятым с предыдущей итерации), а фигурные скобки означают, что уравнения движения решаются согласованно с соответствующими уравнениями постоянства массового расхода.
Отметим, что во вспомогательной задаче (12) и (13) уравнения постоянства массового расхода (9) и (10) рассматриваются как априорная информация об искомом решении, использование которой позволяет уменьшить объем вычислительной работы. Представление давления в форме (11) неоднозначно, поэтому по отдельности каждое из слагаемых рх(х), ру (у) и рху (х, у) лишено физического смысла. В частности, если известно значение давления в одной точке (например, в начале координат р0 = р(0, 0)), то это значение должно быть некоторым образом перераспределено между слагаемыми рх(х), ру(у) и рху(х,у) (например, рх(0) = 0, ру (0) = 0 и рху (0, 0) = р0); очевидно, что погрешности вычисления слагаемых рх(х) и ру(у) взаимно независимы.
Численное решение вспомогательной задачи (12)-(13) осуществляется при помощи метода Зейделя с блочным упорядочиванием неизвестных. Запишем линеаризованный разностный аналог уравнения
i
i
движения по X в (12) в виде
а ^^ л e л s л n л p т i _1
Aij ui_lj + Aij ui+lj + Aij uij_ 1 + Aij uij+1 + Aij = bij + 7 ■
hx
„X „X
Верхними индексами ы, е, в, п и р обозначены соответственно западная, восточная, южная, северная и центральная точки шаблона или грани контрольного объема [1]. Тогда итерация метода Зейделя (д) с блочным упорядочиванием неизвестных по столбцам примет вид
^^ ^^ ^^ ^^ ^^ ррХ _ррХ
Ащ + Аец иг+1_1 + АЦ иц-1 + АЦ и(Ц+1 + АЦ иц — Ь'Ц + Т
Тх
и сводится к решению системы линейных алгебраических уравнений с трехдиагональной матрицей коэффициентов
А5 ц(,1) + Ар и(я) + Ап и(<1) -
- Ьц + _ АЦ _ АЩ. и^ (14)
Тх
методом прогонки.
Нетрудно видеть, что вычислительная процедура для решения (12) аналогична решению уравнений Навье-Стокса в приближении пограничного слоя (5)-(7): сначала система (14) решается методом прогонки по мере сходимости итераций по нелинейности, далее рХ подбирается методом секущих таким образом, чтобы компонента скорости и удовлетворяла разностному аналогу уравнения постоянства массового расхода (9). Численное решение уравнений (13) проводится при помощи метода Зейделя с блочным упорядочиванием по строкам.
Поскольку вспомогательная задача (12) и (13) предназначена для быстрого отыскания "части давления" (т.е. рХ(х) + ру (у)), то необходимо выполнить только несколько итераций метода Зейделя, которые назовем Р-итерациями. Потом значения рХ(х) и ру(у) фиксируются и осуществляется переход к исходной задаче (1)-(3), в которой уравнения движения с учетом (11) принимают вид
д (и2) д (vu) дх ду
dpx dx
dpxy 1 (д2и д2и .
+ ^ тт-2 + ^ ; (15)
дх R^ \ дх2 ду
д (uv) д(v2)
дх
ду
dpy dy
др
xy i /д 2v д 2v ду R^ \ дх2 ду2
(16)
В традиционных методах решения уравнений Навье-Стокса значительную трудность представляет отыскание давления р(х,у); теперь, благодаря вспомогательной задаче (12) и (13), основная трудность состоит в отыскании "части давления", а именно слагаемого рху(х,у).
Очевидно, что эффективность предлагаемого подхода зависит от типа течения. В случае
m
min ^min
dpx
dx
, min
j
dpy
dV
^ max max
\ ij
dp
xy
dx
, max
ij
dp
xy
dV
который характерен для гидродинамических задач с выделенным направлением течения среды, уменьшение объема вычислительной работы будет максимальным.
Описание вычислительного алгоритма. Исходная задача, состоящая из уравнения неразрывности (1) и уравнений движения (15), (16), может быть решена как посредством сегрегированного, так и совместного алгоритмов. В данной работе для наглядности и простоты применялся сегрегированный алгоритм.
Аппроксимация уравнений Навье-Стокса осуществлялась на равномерной разнесенной сетке (Их + 1) х (Иу + 1) с шагами Нх = 1/Их и Ну = 1/Иу. Соответствующие разностные уравнения принимают вид: — разностный аналог уравнения неразрывности
ui+1j uij vij+1 vij
hx
+
hy
= 0;
(17)
— разностный аналог уравнения движения по координате X
(U2)e - (U2)w (vu)n - (vu)
hx
+
hy
-т1 + ¿А''и<>; <18)
— разностный аналог уравнения движения по координате Y
(uv)e - (uv)w , (v2)n - (V2)
hx
+
hy
jj + , <19)
где Ah — разностный оператор Лапласа, имеющий вид
д , Фг~1] - 2Фгз + Фг+1] . Фг]-1 - Цг] + ^j+1 Ahtyij = -гт;--+
hl
h2 hy
Нижние индексы е, ы, в и п в уравнениях (18) и (19) обозначают значения переменных соответственно на восточной, западной, южной и северной гранях контрольных объемов. Значения компонент скорости на гранях контрольного объема аппроксимировались при помощи взвешенных разностей против потока; линеаризация разностных уравнений проводилась методом Ньютона. Введем обозначение
М(ф) =
д (иф) d(vy) 1 f д2ф д 2ф
+
дх ' ду Ие \ дх2 ду2
Классический сегрегированный алгоритм (КСА) для решения стационарных уравнений Навье-Стокса представим в виде следующей последовательности действий.
s
s
1. Задание начальных значений и(0), -у(0) и р(0).
2. Задание пр = 0.
3. Получение начального приближения к решению уравнений Навье-Стокса из решения системы
(0) „(0)
Mft(u) = Mh(v) =
Pij - Pi-
i—lj
hx (0) _ (0) pij - pij- l hy '
4. Задание np = np + 1.
(np)
5. Коррекция давления p,
6. Решение уравнений движения (18) и (19)
ij
Mh (u) =
Mh (v) =
P(np) _ P(np) Vij Vi— lj
hx :
p(nP) _ p(nP)
Vij V ij—l
hy '
(20)
7. Проверка сходимости, если необходимо — возврат к п. 4.
Далее совокупность действий п. 4-7 будет называться итерацией по давлению, а np — соответствующим счетчиком итераций.
Модифицированный сегрегированный алгоритм (МСА), с включением вспомогательной задачи, представим в виде следующей последовательности действий.
1. Задание начальных значений u(0), v(0), (pxy )(0), (px)(0) и (py )(0).
2. Задание np = 0.
3. Получение начального приближения (P-итерации)
px px (Pxy )(0) _ (Pxy )(0)
Mh(u) = - Px- Px—1 (P )ij (P )i—lj
(0)
hx
u(x,y) dy = ...;
hx
Mh(v) i
pj - рУ—i (pxyj - (pxy)g—1
hy
hy
v(x, y) dx = ...
0
4. Задание np = np + 1;
5. Коррекция давления (pxy)(jlp).
(21)
l
6. P-итерации
Mh(u) = i
/ u(x,y)dy =
nx ^x (pxy )(np) __ (pxy )(np)
Pi Pi-1
hx
hx
v 0
<
Mh(v) 1
У У
Py - Py-
1 (pxy j - (Pxy )j-1
hy
hy
v(x, y) dx = ...
7. Решение уравнений движения (15) и (16)
Mh(u) =
Mh(v) =
pxx - PX-1 (Pxy j - (Pxy №
1j
hx
hx
j - Py-1 (Pxy j - (Pxy)j-1
hy
hy
8. Проверка сходимости, если необходимо — возврат к п. 4.
Нетрудно видеть, что вспомогательная задача в МСА выполняет роль предиктора, а основная — корректора.
Проверка сходимости состоит в вычислении погрешности разностного аналога уравнения неразрывности (17); критерий останова итераций по давлению задан в виде
Ruv = max
ij
ui+1j uij + vij+1 vij
hx
hy
< £ .
(22)
где е — малый параметр, принимаемый равным 10 5.
В данной работе коррекция давления осуществлялась при помощи метода искусственной сжимаемости, т.е.
(Пр)
Pj
P ( ПР-1) D Pij - D
ui+1j - u hx
•ij + vij+1 vij hy
в КС А;
(Pxy j := (P
xy)(np-1) _
)ij
Ш
ui+1j - u hx
4] + vij+1 vij
hy
в МСА,
где ш — параметр нижней релаксации.
Способ аппроксимации уравнений постоянства массового расхода (9) и (10) однозначно определяется аппроксимацией уравнения нераз-
рывности, поэтому из (17) непосредственно следует
N
0 = / u(xi, y) dy ~ Uij и 0 = / v(x, yj) dx ~ V v
0 j=1 0 т.е. интегралы приближенно вычисляются при помощи формулы средних значений [5].
Отдельные компоненты данных алгоритмов рассмотрены более подробно при анализе результатов вычислительных экспериментов. Заметим, что P-итерации в МСА можно выполнять не на каждой итерации по давлению. Более того, в любой момент можно перейти от МСА к КСА путем отказа от P-итераций.
Рассмотренные выше алгоритмы могут применяться и для решения нестационарных уравнений Навье-Стокса. В этом случае приведенная последовательность действий выполняется на каждом временном слое.
Вычислительные эксперименты включают расчеты течений в каверне с движущейся крышкой и обтекание ступеньки.
Течение в каверне с движущейся крышкой характеризуется отсутствием выделенного направления движения среды, поэтому предлагаемый подход в этом случае будет наименее эффективен. Рассмотрим простейший случай, когда Re = 100. При малых числах Рейнольдса, как правило, не возникает проблем, связанных со сходимостью итераций по нелинейности и решением уравнений движения (18) и (19).
Выберем нулевое начальное значение компонент скорости и давления u(0) = 0, v(0) = 0 и p(0) = 0. В этом случае компонента скорости v точно удовлетворяет уравнению движения по координате Y и граничным условиям, поэтому в КСА получение начального приближения к решению уравнений Навье-Стокса (20) сводится к решению уравнения
д (u2) 1 (д2и д2и\
( ) = тН ^ + . (23)
i=l
dx Re \ dx2 dy2 Аналогичная задача (21) в МСА имеет вид
д (и2) dpx 1 / д2 и д2 и
дх dx Re \ дх2 ду2
1
J и(х,у) dy = 0.
о
Обе задачи в конечном итоге сводятся к системам вида Ax = b, которые решались при помощи метода Зейделя с блочным упорядочиванием неизвестных. В методологических целях критерий
(24)
i
Рис. 2. Распределение компоненты скорости u(0,5,y) в вертикальном сечении, проходящем через середину дна каверны, при Re = 100:
+ данные работы [7];--расчет; о —
КСА; • — МСА
останова итераций задан в виде
\\Ах — Ь\\ _7 -—- < 10 . Расчет течения
\\о\\
при Ие = 100 проводился на равномерной разнесенной сетке 101x101
(К = ку = 1/100).
На рис. 2 проиллюстрированы преимущества использования уравнений постоянства массового расхода (9) и (10) в качестве априорной информации. На рисунке сплошной линией показано полученное решение уравнений Навье-Стокса в переменных скорость-давление, а маркерами "+" — в переменных вихрь-функция тока [7]. Начальное приближение в КСА, полученное из (23), близко к решению только в непосредственной близости от движущейся крышки (у > 0,7). Начальное приближение в МСА, полученное из (24), близко к решению во всем вертикальном сечении.
Естественно ожидать, что скорость сходимости МСА будет выше, чем КСА. На рис. 3 показано уменьшение погрешности уравнения неразрывности Ка'" (22) в процессе выполнения итераций по давлению (пр). Использование вспомогательной задачи (12) и (13) в МСА привело к уменьшению числа итераций по давлению приблизительно в полтора раза.
Следующий этап тестирования МСА — расчет течения в каверне при Ие = 1000. В этом случае итерации по нелинейности будут расходиться при нулевом начальном приближении, поэтому необходимо принять дополнительные меры для обеспечения сходимости. Принят наиболее распространенный способ получения стационарного решения подобных задач путем счета на установление. Этим же показана применимость МСА к решению нестационарных уравнений Навье-Стокса. Расчет течения при Ие = 1000 проводили на равномерной сетке 301 х 301.
Рис. 3. Зависимость погрешности разностного аналога уравнения неразрывности В""" от числа итераций по давлению пр:
О — КСА; • — МСА
Рис. 4. Распределение компоненты скорости u в вертикальном сечении, проходящем через середину дна каверны (а), и распределение компоненты скорости v в горизонтальном сечении, проходящем через середину боковых стенок каверны (б), при Re = 1000 (+ — данные работы [7])
На рис. 4 сопоставлены компоненты скорости и и v в характерных сечениях каверны с результатами решений уравнений Навье-Стокса в переменных вихрь-функция тока [7]. Изолинии функции тока и изобары показаны на рис. 5.
Число выполненных итераций по давлению (np), необходимое для достижения критерия останова (22) на первых двадцати временных слоях n, показано на рис. 6. В МСА одна P-итерация выполнялась после каждой сотой итерации по давлению, поэтому вычислительная стоимость каждой итерации по давлению в МСА и КСА практически одинакова. Как и при решении стационарных уравнений Навье-Стокса, вспомогательная задача позволила сократить объем вычислений приблизительно в полтора раза.
Рис. 5. Изолинии (а) функции тока (Яе = 1000) и изобары (б) (Яе = 1000, р(0,0) = 0) при течении в каверне
Поскольку в основе вспомогательной задачи лежат упрощения, характерные для уравнений Навье-Стокса в приближении пограничного слоя или длинного канала, то уменьшение объема вычислений в МСА по сравнению с КСА на примере расчета течения в каверне является минимальным.
Задача об обтекании ступеньки при Re = 800 часто применяется для тестирования вычислительных алгоритмов решения уравнений Навье-Стокса. Расчетная схема задачи показана на рис. 7. Рассматривается стационарное течение среды с постоянными свойствами в канале длиной L и высотой H. На входе в канал расположена ступенька размера H/2 х Н/2. Предполагается, что на входе в канал компонента скорости u распределена по параболическому закону, а v = 0. При больших числах Рейнольдса образуются два вихря — непосредственно за ступенькой и на верхней стенке канала.
Запишем интегральные формы уравнения неразрывности для дан-
ной задачи:
H H
Ju(x,y) dy = Ju(-H/2, y) dy при - H ^ x ^ 0;
H/2 H/2
H H
J u(x,y) dy = Ju(-H/2,y) dy при 0 <x ^ L;
0 H/2
L У
J v(x,y) dx = — J u(L,() d( при 0 ^ y ^ H/2; 00
0 4 8 12 16 П
Рис. 6. Зависимость числа выполненных итераций пр по давлению, необходимых для достижения критерия останова, от числа п временных слоев
а ~TL ¡Т X TR 3
'/////
_— —
V " -----------
/ / / / / / / X
п L
Рис. 7. Расчетная схема задачи об обтекании ступеньки
L У У
J v(x,y) dx=-J u(L,£) d£ + j u(-H/2,£) d£ при H/2 <y < H.
-H/2 0 H/2
Для удобства доопределим компоненты скорости, полагая их равными нулю в ступеньке, т.е.
u(x,y)= v(x, y) = 0 при - H/2 < x < 0 и 0 < y < H/2 .
Тогда уравнения постоянства массового расхода упрощаются и принимают вид
H H
J u(0,y) dy = J u(x,y) dy; 00
L У
I = - / (,L,) - ,0,)) dy.
00 В качестве примера выполнен расчет обтекания ступеньки при Re = 800, L = 14H. Расчет проводился на равномерной разнесенной сетке 101 х 1401 (hx = hy = 1/100). Получено, что при более высоких значениях числа Рейнольдса уже не существует стационарного течения.
Изолинии функции тока показаны на рис. 8, а в таблице приводится сопоставление полученных результатов1 с ранее опубликованными данными.
Таблица
Сравнение результатов моделирования течения при обтекании ступеньки
(Re = 800)
Литературный источник 1Б It wT XTL XTR Число узлов
[8] 6,0150 5,6600 - 4,8200 10,4800
[9] 6,1000 5,6300 - 4,8500 10,4800 129681
[10] 6,0820 5,6260 - 4,8388 10,4648 245760
[10] 6,1000 5,6300 - 4,8600 10,4900 ^ 8000
[11] 6,0964 5,6251 - 4,8534 10,4785 3737
Данная работа 6,1000 5,6300 0,28 4,8400 10,4700 141501
Сравнение времени решения задачи об обтекании ступеньки показывает, что объем вычислений в МСА при L/H = 14 приблизительно в 300 раз меньше, чем в КСА. Причина этого понятна из рис. 9, где
1В настоящей работе границей вихрей считались узлы, в которых сеточная функция иу меняет знак, т.е. координаты точек отрыва и присоединения хТь, хтя и ¡в определялись с точностью до шага сетки Нх = 1/100.
Рис. 8. Изолинии функции тока при обтекании ступеньки (Re = 800)
1.0 -|
0.8 -
О 2 4 6 8 10 12 14
Рис. 9. Изобары функции тока при обтекании ступеньки (Re = 800)
показаны изобары. Очевидно, что, за исключением небольшой подобласти, изобары — почти вертикальные линии, т.е. давление в данной задаче изменяется преимущественно в направлении потока. Поэтому "градиент давления" (px)X, отыскиваемый при помощи вспомогательной задачи, будет доминирующим. Преобладание (px)X приводит к существенному снижению трудоемкости вычислений для отыскания слагаемого pxy в (11). Заметим, что при увеличении отношения L/H будет наблюдаться дальнейший рост эффективности МСА по сравнению с КСА.
Выводы. Предложенный алгоритм для решения уравнений Навье-Стокса в переменных скорость-давление на структурированных сетках занимает промежуточное место между сегрегированными и совместными алгоритмами, поскольку компоненты скорости и "часть давления" px + py всегда отыскиваются совместно при решении вспомогательной задачи.
Организация вычислений по схеме предиктор-корректор, где в качестве предиктора используется вспомогательная задача, позволяет уменьшить объем вычислений, необходимых для решения уравнений Навье-Стокса в переменных скорость-давление на структурированных сетках. Уменьшение трудоемкости вычислений особенно заметно при моделировании течений с выделенным направлением течения среды.
Решение вспомогательной задачи не требует глобальной линеаризации матрицы коэффициентов результирующей системы алгебраических уравнений. Поэтому данный метод может быть легко обобщен для моделирования течений сжимаемых сред, теплофизические свойства которых далеки от свойств совершенного газа.
Дальнейшее повышение эффективности алгоритмов численного решения уравнений Навье-Стокса возможно только методами прикладной математики. В частности, используемый в данной работе метод искусственной сжимаемости имеет малую скорость сходимости. Однако вычислительные эксперименты показывают, что на первых итерациях по давлению скорость сходимости выше, т.е. данный метод обладает сглаживающими свойствами. По-видимому, многосеточные методы наиболее перспективны для создания высокоэффективных алгоритмов решения уравнений Навье-Стокса. В этом случае вспомогательная задача служит источником дополнительного сглаживания на самой мелкой сетке.
СПИСОК ЛИТЕРАТУРЫ
1. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости: Пер. с англ. - М.: Энергоатомиздат, 1984. - 152 с.
2. C h o г i n A. J. A numerical method for solving incompressible viscous flow problems // J. Comput. Phys. - 1967. - V. 2. - P. 12-26.
3. Vanka S. P.,Leaf G. K. An efficient finite-difference calculation procedure for multi-dimensional fluid flows // AIAA Paper. - 1984. - № 1244. - P. 12.
4. Ольшанский М. А. Лекции и упражнения по многосеточным методам. -М.: Физматлит, 2005. - 168 с.
5. К а л и т к и н Н. Н. Численные методы. - M.: Наука, 1978. - 512 c.
6. B r i l e y W. R. Numerical method for predicting three-dimensional steady viscous flow in ducts // J. Comp. Phys. - 1974. - V. 14. - P. 8-28.
7. G h i a U.,Ghia K. N.,Shin C. T. High-Re solutions for incompressible flow using the navier-stokes equations and a multigrid method // J. Comp. Physics. - 1982. -V. 48. - P. 387-411.
8. B arton I. E. The entrance effect of laminar flow over a backward-facing step geometry // International Journal for Numerical Methods in Fluids. - 1997. - V. 25. - P. 633-644.
9.Gartling D. A test problem for outflow boundary conditions flow over a backward-facing step // International Journal for Numerical Methods in Fluids. -1990. -V. 11. - P. 953-967.
10. G r e s h o P. M., G a r 11 i n g D. K., T o r c z y n s k i J. R., C l i f f e K. A., Winters K. H., G a r r a 11 T. G., S p e n c e A., G o o d r i c h J. W. Is a steady viscous incompressible two-dimensional flow over a backward-facing step at Re = 800 stable? // International Journal for Numerical Methods in Fluids. - 1993. -V. 17.-P. 501-541.
11. Keskar J.,Lin D. A. Computation of laminar backward-facing step flow at Re = 800 with a spectral domain decomposition method // International Journal for Numerical Methods in Fluids. - 1999. - V. 29. - P. 411-427.
Статья поступила в редакцию Сергей Иванович Мартыненко — окончил МВТУ им. Н.Э. Баумана в 1988 г., канд. физ.-мат. наук, научный сотрудник отдела химмотологии и спецдвигателей ФГУП ЦИАМ им. П.И. Баранова.
S.I. Martynenko graduated from the Bauman Moscow Higher Technical School in 1988. Ph. D. (Phys.-Math.), researcher of department for chemistry of motor combustive-lubricating materials and special motors of the Federal State Unitary Enterprise "Central Institute of Aviation Motors n.a. P.I. Baranov".
