CC BY
88
Вычислительные технологии
Том 6, № 2, 2001
ПРИМЕНЕНИЕ АЛГОРИТМОВ РАСЩЕПЛЕНИЯ В МЕТОДЕ КОНЕЧНЫХ ОБЪЕМОВ *
В. М. КОВЕНЯ Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: kovenya@ict.nsc.ru
The splitting ideology for the efficient solution of Euler equations in the integral form, which are approximated by finite volume method, is applied. This ideology was developed earlier for solving aerohydrodynamics problems. The constructed algorithms are economical and characterized by the minimum dissipation.
Введение
Аппроксимация дифференциальных уравнений Эйлера или Навье — Стокса разностными схемами широко используется при решении различных классов задач аэро- и гидродинамики (см., например, [1-10]). Разработанные алгоритмы: явные и неявные разностные схемы, их комбинации, схемы повышенного порядка с различными формами монотонизации решения и созданные на их основе комплексы и пакеты прикладных программ — позволяют решать двумерные и некоторые трехмерные задачи с достаточной точностью на современных ЭВМ. Однако для ряда задач — задач в многосвязных областях, в областях со сложной геометрией — более эффективными могут оказаться подходы, использующие аппроксимацию расчетных областей различными структурами регулярных и нерегулярных сеток, что позволяет более естественно и просто проводить сгущение расчетных ячеек в отдельных подобластях, например, в областях больших градиентов. При этом исходные уравнения удобнее выбрать в интегральной форме и аппроксимировать их на расчетной сетке, состоящей из ячеек различной конфигурации (треугольных, четырех- или шестиугольных ячеек Дирихле и т. д.), наиболее удобных при решении конкретного класса задач. Аппроксимация уравнений в интегральной форме может оказаться предпочтительнее и на регулярных сетках, например, при решении задач, содержащих разрывы различного рода. Такой подход, носящий название интегро-интерполяционного метода или метода контрольного объема, использовался начиная с 60-х годов при решении двумерных задач [1-5, 9, 10]. Сегодня его называют методом конечных объемов (МКО). Аппроксимация исходных уравнений в интегральной форме на основе МКО приводит к системе линейных или нелинейных алгебраических уравнений большой размерности, как и при аппроксимации дифференциальных уравнений. Их решение представляет самостоятельную проблему
* Работа выполнена в рамках Федеральной целевой программы "Интеграция" (проект №274), программы "Университеты России" и при финансовой поддержке Российского фонда фундаментальных исследований, гранты №99-01-00619, №99-07-90418.
© В.М. Ковеня, 2001.
и находится, как правило, или итерационными методами, или на основе алгоритмов расщепления и приближенной факторизации. В настоящей работе идеология расщепления, развитая в работах [6, 7], применена для построения экономичных алгоритмов МКО.
1. Алгоритмы МКО для уравнений газовой динамики
Представим уравнения газовой динамики в виде интегральных законов сохранения массы, импульса и энергии для произвольного фиксированного объема П с границей ¿П (для простоты изложения рассматривался двумерный случай):
д Г
- шп +
W ¿8 = 0,
(1.1)
до
где
и =
Р
рУ 1 РУ 2
V Е )
w =
Е = р(е + у), у2 = ь( +
( рv \
рь^ + рв1
рV2V + рв2
\ v(E + р) /
v¿ — составляющие вектора скорости v в декартовых координатах хт (т = 1, 2), в1, в2 — базис декартовой системы координат; ¿8 = ий-8 — элемент поверхности, умноженный на единичную внешнюю нормаль п к ней. Для замыкания исходных уравнений (1.1) задано уравнение состояния р = р(р, в). Введем в рассмотрение четырехугольную ячейку (см. рисунок).
Значения сеточных функций будем определять в центре ячейки г,], а потоки на границах — в узлах г ± 1/2,] и г,] ± 1/2. Объем ячейки П обозначим через V, а площадь грани — через Sm±1/2. Введем среднее значение сеточных функций и по объему:
и = и- = V
ШП,
(1.2)
до
а значение контурного интеграла в (1.1) аппроксимируем оператором И, по формуле
Wd8 « и = (^)г+1/2 - (^)г-1/2 + (^),-+1/2 - (£1/2.
(1.3)
до
Разностные потоки И получены на основе симметричной аппроксимации контурного интеграла с вторым порядком точности. Заметим, что £¿+1/2 — £¿-1/2 + £7+1/2 + £7-1/2 = 0 как
следствие интегрального равенства у) = 0. С учетом (1.2), (1.3) рассмотрим разностную
йП
схему с весами
ттга+1 _ ттп Vи-— + аИ,п+1 + (1 _ а)И,п = 0.
т
(1.4)
Здесь и в дальнейшем, где это возможно, индексы г,] опущены. Схема (1.4) аппроксимирует исходные уравнения с порядком 0(тт + Л2) и при а = 0 нелинейна, а т = 2 при а = 0.5 и т =1 при а = 0.5, Л = тах(Лт), Лт = У/$т.
Линеаризуем векторы ип+1 и Wn+1 относительно вектора искомых функций 7
ип+1 = ип + АпДГ + 0(т2), Wn+1 = ^)п + ВпДР + 0(т2),
(1.5)
. 5и - дW
где А = ; В = "дГ;
дгп
7п+1 _ 7п. Вид матриц А, В- зависит от вида искомой функ-
ции 7 В частности, если в качестве искомых функций выбраны плотность, компоненты скорости и давление, то
А
( 1 0 0
У1 Р 0
У2 0 Р
V2
\ 2 РУ1 РУ2
0 0 0 1
\
( г
В
РП1
РП2 р(г + ^2^2)
2
-г
'42
7 _ 1 / 0
П1
П2 7 ,
43
/
7 _ 1
( Р \
У1
У2
\Р/
(1.6)
где г = Уп = г>1П1 + ^2^2; 642 =
7
7 . ,
р+Р1 и1
П1 + ; 643 = +
1 р + р ^ — + ) п2; — проекции вектора нормали к поверхности на координатные оси хг (г = 1, 2). Вид матриц А и В- приведен для уравнения состояния р = (7 _ 1)ре. С учетом (1.5) разностная схема (1.4) может быть преобразована к виду
(УАп + таМп) ДР = —тЯп
(1.7)
Схема (1.7) линейна и аппроксимирует исходные уравнения (1.1) с тем же порядком, что и исходная схема (1.4). Здесь
Мп
(5Вп)г+1/2 _ (5Вп)г_1/2 + (^+1/2 _ (ЯВ*),-^.
(1.8)
В силу симметричной аппроксимации контурного интеграла в (1.1) разностная схема (1.7) немонотонна. Для получения монотонной схемы заменим симметричную аппроксимацию интеграла (1.8) на его несимметричную противопотоковую аппроксимацию первого порядка.
Пусть ^ = Д7п. Тогда справедливо приближенное равенство
2
V
Мф = (£ВТ^)г+1/2 — (£Вга^)г-1/2 + (£Вга— (£Вгаф),--1/2 -- В? [(£ф)т/2 — (£р)г-1/2] + Щ [(£^+1/2 — (£^),-1/2] - (вВПАх + В?А2)ф, (1.9)
где
В А = В Г £7т-1/2А-т при (уга)7т > 0,
В тАт В т 1 о Л ~
I £т + 1/2А+т При (уга)7т < 0,
т =1, 2, а = г, ]2 = ]. Матрицы 51 и 52 взяты в целых узлах ячейки как их средние значения по граням, т.е. вт = (в^т+1/2 + 57т-1/2)/2, а разностные операторы А± определены по формуле А^т фт = ±(фт — Значение проекции нормальной компоненты скорости (ьп)7т на координатное направление хт определяется так же, как среднее ее значение в целом узле:
¿т = (уга)7т = (¿т+1/2 — ¿т-1/2) /2,
где ¿т±1/2 = (ь1п1 + ь2п2)7т±1/2; пт — внешние нормали к поверхности элемента, а знак — взят с учетом направления нормалей на противоположных гранях. Согласно введенным упрощениям и обозначениям (1.9), стабилизирующий оператор схемы (1.7) может быть представлен в виде
+ таМ - + га(Б?А1 + Б^) = ( I + та ^ В£Ат ) , (1.10)
т=1
где
1
Вт = V (А )"
/ ¿т 0 \
0 ¿т 0 гт
0 0 / г2
тт
V 0 ^т ^т ¿т /
д7 = рп7-; г7 = (1/р)п; 7 = рс2п (] = 1, 2); с2 = 7р/р — квадрат скорости звука. Тогда
разностная схема
2
I + та£ -тА^ АР = — V (А-1)пИп (1.11)
т=1 /
аппроксимирует исходные уравнения (1.1) с порядком 0(тт + Л2) и линейна. Ее решение может быть получено методом матричной прогонки.
Приближенно факторизуя оператор I + та ^т=1 ВтАт, получим разностную схему
2 »п+1 £ п 1
П(1 + та -тАт)-— = — V (А-1)пИп (1.12)
т=1
или эквивалентную ей схему в дробных шагах
Г = — -1(А-1)пКп,
(I + та В?А1) |п+1/2 = Г, (I + таВ?А2) Г+1 = |п+1/2,
£ П+1 = £ п + т ^П+1
аппроксимирующую исходные уравнения (1.1) с порядком 0(тт + тЛ + Л2). Как следует из вида матричных операторов ВПД^, разностная схема реализуется на дробных шагах векторными прогонками, требующими обращения матриц размером 4 х 4 в каждом узле расчетной ячейки. Для f = и схема (1.11) при аппроксимации уравнений МКО по структуре близка к разностной схеме Макдональда — Брили [11], а для f = и — к схеме Бима — Уорминга [12], предложенным для решения уравнений в дифференциальной форме.
Замечание 1. Известно, что введение операторов расщепления или приближенной факторизации для построения экономичных разностных алгоритмов приводит к схемам, свойства которых отличаются от нефакторизованных. Это связано с появлением дополнительных членов в стабилизирующем операторе при его факторизации. Так, например, в схеме (1.12) стабилизирующий оператор
П (I + тавтДт) = I +5] Вт + Т2а2ВПД1ВПД2
т=1 т=1
содержит дополнительные члены Q = т2а2ВПД1ВПД2 второго порядка малости, играющие роль диссипативных членов, отсутствующих в исходной нефакторизованной схеме. Это может приводить к повышению погрешности при решении нестационарных задач или замедлению скорости сходимости при нахождении стационарного решения методом установления. Подобная ситуация возникает и при увеличении числа дробных шагов в схемах расщепления и схемах приближенной факторизации. Это "плата" за эффективность алгоритма. При построении таких схем необходимо оценивать и минимизировать влияние расщепления или приближенной факторизации.
2
2
1.1. Алгоритм расщепления по физическим процессам и пространственным направлениям
Для получения экономичных разностных схем в работах [6, 7] введено расщепление исходных операторов по физическим процессам таким образом, чтобы каждый оператор описывал элементарный физический процесс. Это позволило построить разностные схемы, реализуемые на дробных шагах скалярными прогонками, что делало эти алгоритмы экономичными. Как отмечалось в замечании 1, введение операторов расщепления приводит к появлению в разностной схеме дополнительных членов. В работах [13, 14] предложены оптимальные формы расщепления по физическим процессам, при которых влияние расщепления было минимальным, т. е. дополнительные члены появились лишь в уравнении энергии. Подобные алгоритмы могут быть построены и при аппроксимации уравнений в интегральной форме. Следуя [14], введем расщепление операторов Вт в виде
Вт Дт
5>т д™
10 0 0 0 \ ( ¿тД^т ?тДТт 0
0 0 0 гт Д±т + 0 ¿т Д ^т 0 0
0 0 0 гт Д±т 0 0 ¿т Д ^т 0
1° /1 Д /2 Д ¿т Д ^т / 0 0 0 0
Заметим, что подобно [6, 7] аппроксимация конвективных членов в операторах ВтДт выбрана с учетом знака скорости, а для членов с давлением в уравнениях движения — по сопряженным к ним формулам согласно (1.9). Приближенно факторизуя оператор
2
(I + таВтАт) — (I + та—тАт)(/ + та—тАт) в (1.12), рассмотрим разностную схему
П
т=1
£ п+1 £ п 1
(I + та—т Ат) (I + та—т Ат)-— = — -(А-1)пИп
т V
(1.13)
или эквивалентную ей схему в дробных шагах
Г = — ^(А-1)?:^,
(/ + та—1 А1) |п+1/4 = Г, (/ + таВ12А1) |п+1/2 = Г+1/4, (1 + таВ21А2) |п+3/4 = Г+1/2, (1 + таВ22А2) Г+1 = Г+3/4,
£п+1 = + т|п+1. (1.14)
Как и в схеме (1.12), матричные операторы —т в (1.13) и (1.14) аппроксимированы на п-м временном слое. Разностная схема (1.13) аппроксимирует исходные уравнения с тем же порядком, что и (1.11), но реализуется на дробных шагах скалярными прогонками. Действительно, значения компонент вектора вычисляются по явным формулам из первого уравнения схемы (1.14). На шаге п + 1/4 разностные уравнения
¿п+1/4 _ ¿п
¿р ¿р )
п+1/4 2
= С
¿П+1/4 + таг1 А±1Срп+1/4 = СП, СП+1/4 + таг2А±1Срп+1/4 = Срп+1/4 + та*1 А^рп+1/4 + та (V А^+1/4 + /2А^
^п+1/4 ^п+1/4
решаются в такой последовательности: исключая и С2 из последнего уравнения,
п+1/4 1 2
приходим к уравнению относительно
[1 + ат^А^ — т2а2 (/1А^А^ + /2А^1г2А±1)] ¿рп+1/4 = = С? — та (ГА^ + /2АТ1СП) .
Его решение может быть получено скалярной прогонкой, так как оператор А^1гА±1 есть сеточная аппроксимация второй производной на трехточечном шаблоне. Затем явно вы-
>-п+1/4 ^п+1/4
числяются значения и С2 из второго и третьего уравнений схемы.
На втором дробном шаге
С+1/2 + та (VАТ1СП+1/2 + д2А^2П+1/2) = С
С п+1/4
)
Сп+1/2 + та*1 А^
п+1/2 _ ¿.п+1/4
п+ ¿1
Сп+1/2 + та^2АТ1Сп+1/2 = Сп+1/4,
С п+1/2 = ¿п+1/4
2
разностные уравнения могут быть решены независимо друг от друга по неявной схеме бегущего счета или скалярными прогонками. Реализация третьего и четвертого дробных шагов подобна рассмотренным выше. Наконец, из последнего уравнения (1.14) явно определяется новое значение £п+1. В линейном приближении схема (1.13) безусловно устойчива при а > 0.5. Таким образом, реализация предложенной разностной схемы в МКО сводится к скалярным прогонкам, что делает этот подход экономичным.
Подобным образом разностные схемы, основанные на аппроксимации уравнений в интегральной форме МКО, могут быть построены и для трехмерного случая. Отметим, что все схемы приближенной факторизации в случае трех пространственных переменных теряют свойство безусловной устойчивости (см., например, [15]). Изложенные выше алгоритмы описаны для случая четырехгранных ячеек. Очевидно, что аналогично могут быть построены схемы МКО и для других форм ячеек, наиболее удобных при аппроксимации областей со сложными криволинейными границами. Отметим, что в МКО использование ячеек различной формы не приводит к понижению порядка аппроксимации или усложнению аппроксимации при сохранении второго порядка, как в разностных схемах на неравномерной сетке. Не требуется и преобразования координат для сгущения ячеек в подобластях. Разбиение расчетной области на ячейки может быть проведено до начала или в процессе вычислений. В качестве недостатка при аппроксимации МКО отметим усложнение расчетных формул, полученных при усреднении функций по ячейке.
1.2. Алгоритм расщепления по физическим процессам
С учетом введенных упрощений (1.10) нефакторизованная базовая схема (1.11) может быть представлена в виде
[I + таЬ]
£п+1 _ £п
т
— ^(А-1 )пИп,
(1.15)
где
ь = —пА1 +
( ¿1А1 + ¿2А2 0
0 0
( г 52 0
0 г 0 г1
0 0 ь г2
0 к 12 1)
А1 + д21А2 ¿1А1 + ¿2 А2 0
11А1 + /1А2
д2А1 + ^2 0
¿1А1 + ¿2 А2 12 А1 + д2А2
0
г1А1 + г1А2 г2А1 + г2А2
\
¿1А1 + ¿2А2 у
Она аппроксимирует исходные уравнения с порядком О(тт+тЛ+Л2), а при установлении — стационарные уравнения в консервативной форме с порядком О (Л,2).
Введем расщепление оператора ь по физическим процессам, позволяющее свести решение системы уравнений к независимому решению отдельных уравнений, т. е. представим оператор ь в виде суммы
Ь = Ь1 + ь
2=
0 0 0 0 ( г 51 52 0
0 0 0 г1 + 0 £ 0 0
0 0 0 г2 0 0 £ 0
0 11 12 * / 0 0 0 0
Приближенно факторизуя оператор (I + таЬ) , рассмотрим разностную схему
_ гге 1
(I + та^) (I + таЬ2)-= _-(А-1)пЯга (1.16)
т V
или эквивалентную ей схему в дробных шагах
( А-1 )"
/•га V ) тэга
I =--—К ,
(I + таЬ1) |п+1/2 = Г, (I + таЬ2) |п+1 = |п+1/2,
Р+1 = р + т|п+1. (1.17)
Опишем ее реализацию. Нулевой шаг схемы реализуется явно по формуле
Iп = _ ^(А-1)^,
где К и А определены в (1.3) и (1.6). На первом дробном шаге п +1/2 решение системы разностных уравнений
¿«.+1/2 _ ¿те
¿р ¿р >
¿"+1/2 + таг1^2 = ¿"+1/2, ¿2га+1/2 + таг2£"+1/2 = ¿2га+1/2,
¿Г1/2 + та(/^+1/2 + / 2 ¿П+1/2 + ^+1/2) = ¿п (1.18)
^«,+1/2 ^«,+1/2 /., 0\ после исключения ¿1 и ¿2 из последнего уравнения системы (1.18) сводится к решению уравнения
[I + та£ _ т2а2 (/1Г1 + /2Г2)] £рп+1/2 = _ та (/^ + (1.19)
¿п+1/2 ¿га+1/2
и явному вычислению оставшихся компонент ¿1 и ¿2 . Как следует из вида разностных операторов /г и гг, уравнение (1.19) содержит аппроксимацию вторых производных (смешанных и повторных) по каждому направлению х1 и х2 и его решение может быть получено по какой-либо итерационной схеме. На втором дробном шаге п + 1 все уравнения
¿П+1 + та^1 + ^¿Г1 + ^¿Г1) = ¿п+1/2, ¿П+1+та^ег1 = ¿п+1/2,
¿П+1 + та^еп+1 = еп+1/2,
¿П+1 = еп+1/2 (1.2°)
решаются независимо друг от друга также по какой-либо итерационной схеме. Наконец, новые значения функций Р+1 находятся из последнего уравнения схемы в дробных шагах (1.17) по формулам Р+1 = Р + т|"'+1. Таким образом, разностная схема (1.16), основанная на расщеплении операторов по физическим процессам, позволяет свести решение системы уравнений к решению отдельных уравнений.
1.3. Диссипативные свойства алгоритмов
Проведем анализ диссипативных свойств рассмотренных выше схем приближенной факторизации: с расщеплением операторов по пространственным направлениям (схема (1.12)), расщеплением по пространственным направлениям и физическим процессам (1.13) и с расщеплением по физическим процессам (1.16). Введем следующие обозначения:
Г
т
таг! Д
±т:
Гг
та1тД^т, Т^
1 + та£тД™, То = Т? + Т2 - 1, где ко-
эффициенты ¿т, ^, г^, определены в (1.9), (1.10), а знаки в операторах выбраны в зависимости от знака скорости ¿т. С учетом введенных обозначений стабилизирующий оператор Со = I + таЬ в нефакторизованной схеме (1.15) может быть представлен в виде
Со
( То з? + 32 3? + 32 о о
V 0
д
и содержит члены типа 1 + та—— или та
То о
о То
¿? + ¿2 ¿2 + ¿2 д
о
В? + Л + Л
То
\
дх
Стабилизирующий оператор С? ризации (1.12) представим в виде
/Т?Т2 Т?32 + Т231 Т?32 + Т23?
П (I + таВПДт) для схемы приближенной факто-
т=?
С? =
о Т?Т2 + в? ¿2
V
о о
т?ь? + т^т?
Т?Т2 + ^^ Т?Ь2 + Т2Ь2
3?д? + з?^2 \ я2т? + д?Т2
т?^2 + Т2Д2 Т?Т2 + + Л^? /
Как легко видеть, диссипативная матрица Д? = С? — С0 содержит дополнительные члены порядка 0(т2), т.е. в каждом разностном уравнении при временных производных содер-
дд
22
жатся диссипативные члены т а
^ •
джт дх
Для разностной схемы (1.13) стабилизирующий оператор имеет вид
о о о о \
о (Т2 — 1)л?^2 Л? ¿2(Т2 — 1) о
С2 = С?+ о (Т2 — 1№? Л2 ¿2(Т2 — 1) о
\о т?ь?(Т2—1)+Т2Ь?(Т? — 1) т?¿2(Т2—1)+Т2Ь2(Т?—1) (т? — 1)(ь?л?+¿2я2)/
и содержит все дополнительные члены, как и в операторе С?, и, кроме того, дополнительные члены третьего порядка малости. Наконец, для схемы приближенной факторизации (1.16) с расщеплением по физическим процессам стабилизирующий оператор С3 имеет вид
о о о о
о о о о
о о о о
V о Ь? + ¿2 Ь + ¿2 0 у
Сз = (I + та£?)(1 + та^2) = Со + (¿о — 1)
и содержит дополнительные члены порядка 0(т2) лишь в последней строке, т.е. диссипативные члены присутствуют только в разностном уравнении энергии.
2
На основании вышеизложенного можно сделать следующие выводы: схема расщепления по физическим процессам (1.16) по своим свойствам наиболее близка к нефактори-зованной схеме. Ее реализация на каждом дробном шаге требует итераций для каждого уравнения, а не для системы уравнений, как в нефакторизованной схеме (1.15). Схемы приближенной факторизации (1.12) и (1.13) близки по своим свойствам, но значительно отличаются от нефакторизованной схемы (1.15), так как содержат в каждом уравнении дополнительные члены порядка 0(т2) и 0(т3). Заметим, что схема (1.13) реализуется на дробных шагах скалярными прогонками, а схема (1.12) — векторными.
Замечание 2. Еще большее отличие в диссипативных операторах получается при рассмотрении схем приближенной факторизации для трехмерных уравнений. Наряду с дополнительными членами порядка 0(т2), в построенных схемах за счет расщепления и приближенной факторизации возникают члены порядка 0(т3). Напомним, что в трехмерном случае схемы приближенной факторизации теряют свойство безусловной устойчивости. Можно ожидать, что среди рассмотренных схем наиболее экономичной по числу итераций будет схема (1.16) с расщеплением операторов по физическим процессам.
2. Алгоритмы МКО для уравнений Эйлера несжимаемой жидкости
Представим уравнения Эйлера в виде интегральных законов сохранения для произвольного фиксированного объема П с границей ^П:
м-д у Шу + у
п да
0,
(2.1)
где
и
VI
ТО"
( V \
VIV + рб1 ^2 V + рв2
\ ^ + рез )
М =
0000 0 10 0 0010 0001
Все обозначения аналогичны приведенным в п. 1. Сеточные функции будем определять в центре ячейки ^'1,^2, ^'з, а значения потоков через границу — в дробных ячейках ^±1/2 (т = 1, 2, 3). Объем ячейки обозначим через V, а сеточные функции определим как среднее значение по объему:
V £ ШП. (2.2)
и
Ил:шз
Ниже, где это возможно, сеточные индексы опущены.
Контурные интегралы по ячейкам, как и в п. 1, аппроксимируем со вторым порядком точности по формулам
¿П
Я =¿[(5^т + 1/2 - (^т-1/2]
т=1
( Яо \
Я1 Я2
Я3
(2.3)
Очевидно, что выполняется условие ^ ^(5?т+1/2 — ^?-т-1/2) = 0 как следствие интегрального
/т=1
= 0. С учетом введенных аппроксимаций (2.2), (2.3) рассмотрим разност-
5
ную схему с весами
Ип+1 — ип УМ-+ аК^1 + (1 — а)Яга = 0,
т
аппроксимирующую уравнения (2.1) с порядком 0(тг + Л2), где I = 2 при при а = 0.5, Л, = шах(У/$).
Линеаризуем вектор Wra+1 относительно вектора искомых функций и:
а
(2.4) 0.5 и I = 1
W
+1
Wra + В? ди + 0(т2), ди
и +1 — и
где
В1
/0 Щ Щ Щ \
_ щ £ + ^1п2 ^1п3
5 И Щ ^2^1 £ + ^2^2 ^2^3
\ Щ ^3^1 £ + ШП2 £ + ^зПз )
щ — проекции вектора нормали на координатные оси ж^; £ = гп = г1щ1 + г2щ2 + г3щ3 — проекция вектора скорости на внешнюю нормаль щ. Тогда с учетом линеаризации разностная
схема (2.4) может быть представлена в виде
МДИ + ^ N1 = — У Ип.
(2.5)
Здесь N1 = Ет=1 [(ЗДДИ)^+1/2 — (ЗДДИ)^-1/2].
Схема (2.5) аппроксимирует уравнения (2.1) с тем же порядком, что и базовая схема (2.4). Представим матричный оператор В1 в виде
в = в + в,
0 =
0 П1 П2 П3 \ 0 0 0 0
П1 £ 0 0 + 0 г1П1 г1П2 г1П3
П2 0 £ 0 0 г2П1 г2П2 г2П3
\ Щ 0 0 £ / 0 г3П1 г3П2 г3П3 )
(2.6)
Тогда N1 = N+N0, N0 = ^[(£Во)ДИ);т+1/2 — (5(Во)ДИ)^т_1/2]. Так как ДИп = Ип+1 — Ип
т=1 3
то N0 = — N = 5][(5ВоИп+1)^т+1/2 — (5ВоИп+1)^т_1/2] — ^[(5ВоИп)^+1/2
т=1
т=1
N0 = £
(5В0Ип)^т-1/2]. Но вектор N0, согласно (2.6), можно записать в виде
0
+1/2 — (г1^Ьт-1/2 (г2^Ьт+1/2 — (г2^Ьт_1/2
V (г3^т+1/2 — (г3^т_1/2 У как следствие равенства нулю сеточного аналога
3
До = ^][(^т+1/2 — (^т_1/2] = 0
т=1
0 0
г1До 0
г2До 0
\ г3До / 0
(2.7)
т=1
3
3
интегрального уравнения неразрывности у) = 0. Заметим, что приближенное равенство (2.7) выполнено с погрешностью 0(Л2). Тогда с учетом упрощений (2.6) разностная схема (2.5) примет вид
МДИ + ^ N = - V я, (2.8)
где
3
N = ^][(££ДИ)^+1/2 - (£ВДИ)з„—1/2].
т=1
Как и схема (2.5), она аппроксимирует уравнение (2.1) с погрешностью 0(т1 + Л2), линейна относительно Ига+1 и немонотонна в силу симметричной аппроксимации контурного интеграла в (2.1).
Как и в п. 1, проведем замену симметричных аппроксимаций в И. на несимметричные с первым порядком по пространству. Пусть ф = ДИП. Тогда
Еф = [(£^„+1/2 - (£Вф)з„—1/2]
' Зт
т=1
~ ^ Вт [(£Ф)З„+1/2 - (£Ф)Зт-1/^ = ^ ВтДтФ, (2^
т= 1 т= 1
где
В Д = В Г £Зт —1/2 Д—т при ¿т > 0,
втДт вт 10 Л „ ^ п
I £Зт+1/2Д+т при ¿т < 0;
Вт = (Взт+1/2 + Взт—1/2)/2; ¿т = (¿т+1/2 - ¿т—1/2)/2, а разностные операторы Д± определены по формулам Д^Зт = ±(1 - (I — единичный оператор, Т^т — оператор сдвига в направлении хт). При усреднении величин следует учитывать знак вектора внешней нормали к поверхности граней.
Согласно (2.9), разностная схема (2.8) после упрощений может быть представлена в виде
м+^ Е вт д т) Дип = - у (2.10)
т=1 /
Она линейна и аппроксимирует уравнения (2.1) с порядком 0(т1 + тЛ + Л2). В силу вырожденности матрицы М (уравнение неразрывности стационарно) непосредственно разрешить разностные уравнения (2.10) не удается, поэтому обычно применяется либо метод расщепления [6, 15, 18], либо вводится слабая сжимаемость уравнений (см., например, [20]), после чего модифицированные уравнения становятся системой Коши — Ковалевской [20, 21].
2.1. Алгоритмы расщепления в модели искусственной сжимаемости
Следуя [20], введем в уравнение неразрывности дополнительные члены, моделирующие слабую сжимаемость жидкости:
[р^У + /(УпЫЗ = 0 ] /
П ¿П
3
3
(- — малый параметр). Тогда исходные уравнения Эйлера можно представить в виде
д
Мо
'д*
ШУ + ф Wds = 0.
(2.11)
Здесь
/ р \
И = г>1 ; Мо =
(- 0 0 0
р
0 1 0 0
0 0 1 0
0 0 0 1
Повторяя все выкладки, проведенные при построении разностных схем для уравнений Эйлера (2.1), рассмотрим разностную схему типа (2.10):
м0 + ^¿втА
■п
тАт
т=1
аи" = - у и"
(2.12)
аппроксимирующую уравнение (2.11) с порядком 0(тг + тЛ + Л2) и линейную относительно вектора И. Она может быть реализована матричной прогонкой, что делает данный подход неэкономичным. В работах [22-24] для решения стационарных и нестационарных уравнений Навье — Стокса предложены численные алгоритмы повышенного порядка, использующие неявный метод конечных объемов и модель искусственной сжимаемости с различными формами их реализации итерационными схемами.
Построим разностные схемы, использующие расщепление оператора в (2.12) по пространственным направлениям. Приближенно факторизуя стабилизирующий оператор
3
3
Мп + V ^ БтЛт ~ МоП I + -^М0"1ВтДт
то—т ~ М0
т=1 ,7=1
с,
рассмотрим разностную схему
и"+1 И" 1
С И-— = -1 И"
т V
(2.13)
или эквивалентную ей схему в дробных шагах
¿" _ 1 тэ"
* = - Vи ,
Мо + ^угВГА^ *"+1/3 =
(Мо + *"+2/3 = Мо*"+1/3,
(Мо + А3) = Мо*"+1/3, И"+1 = И" + т*"+1.
Схема (2.13) аппроксимирует уравнения (2.11) с тем же порядком, что и базовая схе-
та
ма (2.12) и, как следует из структуры матричных операторов Мо + -^ВтАт (см. (2.6)),
на каждом дробном шаге реализуется векторными прогонками вдоль каждого пространственного направления хт.
Рассмотрим дальнейшее расщепление одномерных операторов Вт по физическим про-
Вт = вт+В
тт
Вт в виде
0 00 0
0 ¿т 0 0
0 0 ¿т 0
0 00 ¿т У
+
( 0 (щ)
Ы (П2)т \
Ыт 0 0 0
Ыт 0 0 0
\ (Пз)т 0 0 0 /
(2.14)
и приближенно факторизуем стабилизирующий оператор в схеме (2.12) по формулам
.... - - хх ч г+1ВтДт) (/+-угм—1ВтДт) = с.
З=1
Тогда разностная схема
33
:>>т+Вт )Дт« Мех;
т=1 З=
Ип+1 Ип 1
(7 И-— = -1 И
т V
или эквивалентная ей разностная схема в дробных шагах
/•п 1 тэ п
* = - VИ ,
(2.15)
(ме + ^^аВ! Д^ *п+1/6 = *п,
Ме + *п+2/6 = Ме*п+1/6,
(Ме + уВЗДз) *п+5/6 = Ме*п+4/6,
Ме + ^В32Д^ *п+1 = Ме*п+5/6
V та
V
Ип+1 = Ип + т *п+1 (2.16)
аппроксимирует исходные уравнения (2.11) с порядком 0(тг + тЛ + Л2) и реализуется скалярными прогонками. Действительно, на нечетных шагах (1-м, 3-м и 5-м) система разностных уравнений
£
£
п+(2т—1)/6 . ™ ,п д ¿п+(2т—1)/6 _ ¿п+(т—1)/3
п+(2т—1)/6 = £п+(т-1)/3
п
О" = 1, 2, 3)
решается независимо для каждой компоненты вектора £п+(2т—1)/6 скалярными прогонками или по неявной схеме бегущего счета.
На четных шагах решается система уравнений
С+т/3 + (щ )Дт6
З=1
п+т/3 = £п+(2т—1)/6
£п+т/3 + ^(ПЗ)тДт£Гт/3 = езп+(2т—1)/6 О' = 1, 2, 3).
2
Исключая £п+т/3 Из первого уравнения системы, приходим к разностному уравнению относительно £р:
I -
т 2а2рп V 2е
N
^ у(п7)тАт(п3
3=1
лп+т/3 _ лп+(2т-1)/6
тар'
^ ^ (nj )т 3 = 1
п+(2т-1)/6
тЧ3 •
Его решение может быть получено скалярными прогонками вдоль каждого пространственного направления хт. После вычисления £р явно определяются значения функции . Наконец, из последнего уравнения системы (2.16) явно находятся значения функций и на новом временном слое. Таким образом, разностная схема (2.16) реализуется скалярными прогонками, их число равно Ид (И — размерность задачи по пространству, д — число уравнений). Для рассмотренного трехмерного случая Ид _ 12.
Подобно схеме расщепления, рассмотренной в п. 1.2, строится схема с расщеплением по физическим процессам для уравнений Эйлера слабосжимаемой жидкости.
Замечание 3:
1. Для получения решений нестационарных уравнений Эйлера (2.11) с достаточной точностью необходимо провести серию расчетов при различных значениях е.
2. При решении стационарных задач методом установления необходимо добиваться достижения сходимости к стационарному решению с высокой точностью, чтобы избежать влияния на него искусственной сжимаемости. Это может привести к большому числу итераций.
3
т
2.2. Алгоритмы расщепления в уравнениях Эйлера
Как отмечалось выше, в силу вырожденности матрицы М не удается разрешить разностную схему (2.10). В работах [16-24] для решения уравнений Эйлера и Навье — Стокса несжимаемой жидкости рассматривались специальные алгоритмы расщепления с явной аппроксимацией конвективных (и вязких) членов, неявной аппроксимацией давления и уравнения неразрывности. Построенные схемы являлись условно устойчивыми. Полностью неявная схема типа стабилизирующей поправки предложена в работе [8]. Ниже рассмотрены полностью неявные разностные схемы с различной формой расщепления уравнений.
Представим нефакторизованную разностную схему (2.10) в виде
[М + та(В1 + В2)]
Ип+1 - ип
т
- У *
(2.17)
соответствующем расщеплению матричных операторов по физическим процессам, где
13
у ^ ^ВтАт
т=1
Вт дт
1 3 1 3 у вт дт+у
т=1
Вт Дт
В1 + В2,
( 0
о
о
о о \ о о
т=1
о о
о о
о
о
^т
Вт2 Дт
/ 0 (п1)тДТт (П-2)тДТт (^3)тДТт \
(П1)тД±т 0 0 0
(п2)тД±т 0 0 0
\ (п3)тД±т 0 0 0 )
Приближенно заменим стабилизирующий оператор схемы (2.17) на факторизованный:
(2.18)
М + та(В1 + В2) « (I + таВ1 )(М + таВ2) = М + та(В1 + В2) + т2а2(В1В2) = Сь так как В1М = В1. Разностная схема
(1
Ип+1 - ип
-V и
(2.19)
или эквивалентная ей схема в дробных шагах
*п = - У Ип, (/ + таВ1) *п+1/2 = *п,
(М + таВ2) *п+1 = *п+1/2, Ип+1 = Ип + т*п+1 (2.20)
аппроксимирует исходные уравнения Эйлера (2.1) с тем же порядком 0(тг + Л2), что и базовая нефакторизованная схема (2.10).
Остановимся на ее реализации. Нулевой шаг *п = - у И" определяется явно. Полагаем, что на п-м временном слое уравнение неразрывности выполнено, т. е. ф = 0 и его
сеточный аналог Д = ^ [(£гп)зт+1/2 - (£гп)зт —1/2] = 0. Тогда
т=1
£1 £2 V £3
1
V
( Де \
Д1 Д2
\Д3/
1
V
0
Д1 Д2 V Д3 у
где ДЗ- — аппроксимации потоков через границу (см. уравнение (2.3)). На первом дробном шаге схемы (2.20) система разностных уравнений
¿п+1/2 _ ¿т = п
£е = £е = 0,
Г1/2 + у Г сд*тГ1/2 = Г1/2 с? = 1, 2, з)
З=1
может быть решена независимо для каждой компоненты £п+1/2 одним из известных методов. На втором дробном шаге разностные уравнения
^тГ1 = £
33
п+1 _ ¿-п+1/2
З т т З е
Г1 + (пз)кДкГ1 = £зп+1/2 (? = 1,2,3)
т=1З=1
3
(2.21)
к=1
3
0
решаются в такой последовательности: исключая .+1 из разностного аналога уравнения неразрывности, приходим к уравнению относительно £П+1:
33 2 2 3 3 / 3
^ Е Е^")тДт^Г1/2 = ^ Е Е Е^^ 1 • (2-22)
Уравнение (2.22) — разностный аналог уравнения Пуассона, содержащего смешанные производные и известную правую часть. Его решение может быть получено известными методами (см., например [9, 19, 24]), после чего явно определяется значение .+1. Заметим, что разностный аналог стационарного уравнения неразрывности (первое уравнение системы (2.21)) тождественно удовлетворяется на слоях п и п +1 в силу соотношения
3 3 3 3
Е Е(П)тДт< = Е Е(П)тДт^™+1,
т=1 7=1 т=1 7=1
так как £П+1 = — . Наконец, из последнего уравнения схемы (2.20) явно вычисляются значения функций на п + 1 слое. Таким образом, реализация схемы (2.19) сводится к решению отдельных уравнений для компонент скорости и решению уравнения Пуассона для давления (точнее, для их невязок). Заметим, что предложенная схема близка к конечно-разностной неявной схеме типа стабилизирующей поправки, рассмотренной в [8] для решения уравнений Навье — Стокса несжимаемой жидкости в дифференциальной форме.
Схема приближенной факторизации с расщеплением операторов по физическим процессам и пространственным переменным строится аналогично схеме (2.19). Заменяя стабилизирующий оператор схемы (2.17) по формуле
М + та(В1 + В2) « П (7 + Т^т) (м + у В^) = с
.7 = 1
2,
получим разностную схему
ип+1 — т тп 1 С-= — - Яп (2.23)
или эквивалентную ей схему в дробных шагах
¿п_ 1 "О5
1 = — V к
I + Г^в1) |п+1/6 = м + £в?) Г+2/6 = Г+1/6,
+ ^ ^п+5/6 = ^п+4/6
М + ^^ав!) Г+1 = Г+5/6,
ип+1 = ип + т^п+1. (2.24)
3
Как следует из вида разностных операторов В^, разностная схема (2.24) на дробных шагах реализуется скалярными прогонками. Покажем это на примере первых двух шагов. На первом шаге система разностных уравнений
¿П+1/6 _ ¿и _ о _ п
е™+1/6 + Т%Л1Г1/6 — з с? — 1,2, з)
та
V
^га+1/6
решается независимо для каждой компоненты ^ скалярными прогонками или по неявной схеме бегущего счета при несимметричной противопотоковой аппроксимации конвективных членов. На втором дробном шаге система уравнений
3
уЕ к )1Л1Г1/3 — Г1/6,
3 = 1
С+1/3 + та к-)1Л1еРи+1/3 — с+1/6 с? — 1,2, з) (2.25)
^га+1/3
после исключения из первого уравнения приводится к виду
3
та £М1Л1 к-)1Л1^г1/3 — )1Л1е™+1/6.
3=1 3=1
Его решение может быть получено методом прогонки, после чего явно находятся £П+1 из последних уравнений системы (2.25). По второму и третьему пространственным направлениям реализация схемы подобна рассмотренным выше для 1-го и 2-го дробных шагов. Наконец, новые значения функций ии+1 вычисляются явно из последнего шага схемы (2.24).
2.3. Диссипативные свойства алгоритмов
Диссипативные свойства разностных конечно-объемных схем (2.13), (2.15) для численного решения уравнений (2.12) подобны свойствам схем для решения уравнений Эйлера сжимаемого газа (см. п. 1.3). Остановимся на диссипативных свойствах схемы приближенной факторизации с расщеплением по физическим процессам (2.19). Введем обозна-
3 3 3
та та - та V--«., ,
чения: Ц — — ^ ¿тЛ^т, ¿0 — 1 + V, Ь — — ^(п^)тЛтт, Ь — — ^(п^)тЛтт, где
т=1 _ т=1 т=1
верхний знак в операторах Ь и Ь^ выбирается при ¿т > 0 и нижний — при ¿т < 0. С учетом введенных обозначений стабилизирующий оператор нефакторизованной схемы (2.17) представляется в виде
Ь \
0. 0.
¿0
в схеме приближенной факторизации с расщеплением по физическим процессам (2.19) в уравнениях движения возникают дополнительные члены порядка 0(т2). Действительно,
С1 — (I + таВ1) (М + таВ2) — С0 + т2а2В^ — С + Б,
С0 — М + та(В1 + В2)
( 0 Ь1 Ь
Ь ¿0 о
Ь ¿0 0
V Ь 0 0
0 \
0
0 .
0
Если рассмотреть разностную схему (2.19), (2.20) с заменой второго дробного шага на первый, а первого на второй, то диссипативная матрица равна
¿2^ ¿3^ \ 00,
00, 00
т. е. дополнительные члены возникают лишь в уравнении неразрывности. Но они равны
3
. — "П) = 0,
. =1
и эта схема совпадает с нефакторизованной схемой (2.17). Как легко показать в конечно-объемной разностной схеме (2.23), диссипативные члены второго и более высоких порядков возникают во всех уравнениях. Замечание 4.
1. В настоящей работе построены разностные схемы на основе МКО для численного решения уравнений Эйлера сжимаемого газа и несжимаемой жидкости. Подобные схемы могут быть обобщены и для уравнений Навье — Стокса сжимаемого теплопроводного газа и несжимаемой жидкости. При этом вязкие члены в уравнениях Навье — Стокса могут аппроксимироваться на дробных шагах, например, совместно с конвективными членами.
2. Построенные алгоритмы (2.13), (2.15), (2.23) ориентированы на регулярные сетки. Алгоритмы с расщеплением по физическим процессам наряду с регулярными сетками могут применяться и при решении уравнений Эйлера и Навье — Стокса на нерегулярных неструктурированных сетках, так как их реализация сводится к последовательному решению отдельных уравнений для компонент скоростей и давления. Эти алгоритмы могут служить базовыми для решения задач на многопроцессорных ЭВМ.
Список литературы
[1] Флетчер К. Вычислительные методы в динамике жидкостей. Т. 1, 2. М.: Мир, 1991.
[2] Андерсон Д., Таннехил Дж., Флетчер Р. Вычислительная гидродинамика и теплообмен. Т. 1, 2. М: Мир, 1991. 123 с.
[3] Численные методы решения многомерных задач газовой динамики / С. К. Годунов, А. В. Забродин, М.Я. Иванов и др. М.: Наука, 1989. 614 с.
[4] Роуч П. Вычислительная гидродинамика. М.: Мир, 1980. 616 с.
[5] Пинчуков В. И., Шу Ч.-В. Численные методы высших порядков для задач аэрогидродинамики. Новосибирск: Изд-во СО РАН, 2000. 231 с.
[6] КовЕня В. М., Яненко Н. Н. Метод расщепления в задачах газовой динамики. Новосибирск: Наука, 1981. 304 с.
Б
/000 д! 1 0 0 № 0 0 \ д!3 00
Б = т 2а2В2В1
0 ад
0 0
0 0
0 0
[7] КовЕня В. М., Тлрнлвский Г. А., Черный С. Г. Применение метода расщепления в задачах аэродинамики. Новосибирск: Наука, 1990. 243 с.
[8] Толстых А. И. Компактные разностные схемы и их приложения к проблемам аэродинамики. М.: Наука, 1990. 230 с.
[9] Самарский А. А., Николаев Е. С. Методы решения сеточных уравнений. М.: Наука, 1978. 501 с.
[10] Самарский А. А. Теория разностных схем. М.: Наука, 1989. 614 с.
[11] Briley W. R., McDonald H. Solution of the 3D compressible Navier — Stokes equations by an implicit technique // Lect Notes in Phys. 1975. Vol. 35.
[12] Beam R. M., Warming R. F. An implicit finite-difference algorithm for hyperbolic systems in conservation law form //J. Comput. Phys., 1976. Vol. 22. P. 87-108.
[13] КовЕня В. М., Лебедев А. С. Модификация метода расщепления для конструирования экономичных разностных схем // Журн. вычисл. математики и мат. физики. 1994. Т. 34, №6. С. 886-897.
[14] КовЕня В. М. Схемы расщепления в методе конечных объемов // Журн. вычисл. математики и мат. физики. 2001. Т. 41, №1.
[15] КовЕня В. М. Методы вычислений (дополнительные главы): Уч. пособие. Новосибирск: НГУ, 1995. 92 с.
[16] БЕлоцЕрковский О.М., Гущин В. А., ЩЕнников В. В. Метод расщепления в применении к решению задач динамики вязкой несжимаемой жидкости // Журн. вычисл. математики и мат. физики. 1975. Т. 15, №1. С. 197-207.
[17] БЕлоцЕрковский О. М. Численное моделирование в механике сплошных сред. М.: Наука, 1984. 520 с.
[18] Марчук Г. И. Методы расщепления и переменных направлений. М.: Отдел вычисл. математики АН СССР. 1986. 334 с.
[19] Марчук Г. И. Методы вычислительной математики. М.: Наука, 1980. 536 с.
[20] ЯнЕнко Н. Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, 1967. 196 с.
[21] Chorin A. J. A numerical method for solving incompessible viscons fow problems //J. Comput. Phys. 1967. №2. P. 12-26.
[22] Пейре Р., Тейлор Т. Д. Вычислительные методы в задачах механики жидкости. Л.: Гидрометеоиздат, 1986. 320 с.
[23] Грязин Ю.А., Черный С. Г., Шаров С. П. Численное моделирование течений несжимаемой жидкости на основе метода искусственной сжимаемости // Вычисл. технологии. 1995. Т. 4, №13. С. 180-203.
[24] Грязин Ю.А., Черный С. Г., Шаров С. П. Об использовании методов типа попеременно-треугольных решений неявных разностных схем для трехмерных уравнений динамики несжимаемой жидкости // Там же. С. 306-320.
Поступила в редакцию 4 декабря 2000 г.
