Алгоритмы расщепления в методе конечных объемов Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 20, № 5, 2015

Алгоритмы расщепления в методе конечных объемов

В.М. Ковеня1'2, П. В. Бабинцев2'*

1 Институт вычислительных технологий СО РАН, Новосибирск, Россия 2Новосибирский национальный исследовательский государственный университет, Россия *Контактный e-mail: [email protected]

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

Ключевые слова: уравнения Эйлера и Навье — Стокса, метод расщепления, метод конечных объемов, регулярное и маховское отражение волн.

Введение

Уравнения Навье — Стокса сжимаемого теплопроводного газа являются наиболее общей моделью для описания различных классов задач в аэро- и гидродинамике. При различных предположениях о характере течений из них могут быть получены уравнения Эйлера, пограничного слоя, уравнения вязкого и ударного слоя, параболические уравнения и т. д. Поэтому остается актуальной задача построения эффективных численных алгоритмов полных уравнений Навье — Стокса. К настоящему времени литература по разработке и использованию различных классов алгоритмов методов конечных разностей (МКР) и конечных объемов (МКО) для численного решения уравнений Эйлера и Навье — Стокса составляет тысячи публикаций, поэтому дать достаточно полный обзор алгоритмов практически невозможно. Некоторые подходы описаны, например, в [1-17]. Использование конечно-разностных схем в криволинейных областях приводит к необходимости введения преобразований координат, являющихся отображением исходной расчетной области в прямоугольную. Как следствие, это значительно усложняет уравнения и соответственно численные алгоритмы, поэтому в последние два десятилетия особое внимание уделялось построению алгоритмов МКО [9-15].

При разработке алгоритмов МКО основные проблемы сводятся к аппроксимации исходных областей дискретными ячейками различных типов и интегральных операторов исходных уравнений, а также поиску способов решения полученных систем алгебраических уравнений. Отметим, что современные алгоритмы МКО основаны главным образом на развитии метода Годунова с целью повышения порядка аппроксимации и получения монотонных решений [2].

© ИВТ СО РАН, 2015

При решении многомерных задач использование явных алгоритмов не всегда целесообразно в силу жестких ограничений на соотношение временного и пространственных шагов сетки. Особенно это важно при нахождении стационарного решения на основе принципа установления. Неявные алгоритмы свободны от этих ограничений (или они значительно слабее), но при решении многомерных задач их прямое использование приводит к необходимости обращения матриц большой размерности, требующих значительных вычислительных ресурсов при реализации алгоритмов.

Для построения эффективных алгоритмов чаще всего используются подходы с применением ЬИ-факторизации [9-18], позволяющей упростить решение и свести его к схемам типа бегущего счета. Однако это требует обращения матриц размера т х т в каждом узле расчетной ячейки, где т — число уравнений исходной системы. Ниже используется другой подход, основанный на методе расщепления по физическим процессам [7, 8].

Идеология метода расщепления состоит в сведении решения многомерной задачи к последовательному (или параллельному) эффективному решению одномерных задач или задач более простой структуры. Так как расщепление исходных уравнений может быть достигнуто различными способами, как следствие, это может приводить к различным классам численных алгоритмов, построенных на основе расщепления. В работах [7, 8, 20] метод расщепления использовался при построении экономичных численных алгоритмов МКР решения уравнений Эйлера и Навье — Стокса сжимаемого теплопроводного газа. В настоящей статье этот подход был перенесен на метод конечных объемов при аппроксимации уравнений Эйлера и Навье — Стокса, записанных в интегральной форме.

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

1. Основы алгоритма расщепления в методе конечных объемов

Приведем уравнения Эйлера в интегральной форме в декартовых координатах [7, 18]

Здесь V — замкнутая область произвольного объема; дУ — замкнутая поверхность объема; = п • dS — элемент поверхности в, умноженный на единичную нормаль п к ней; и — вектор искомых функций, Ш — матрица, состоящая из векторов столбцов потоков через границы поверхности. Ниже при построении и реализации алгоритма для

(1)

простоты изложения ограничимся рассмотрением двумерных уравнений. Обобщение алгоритмов на трехмерные уравнения подобно двумерному случаю. Тогда в системе уравнений (1) вектор и и матрица Ш соответственно равны

и

рЬх

РУ2

V Е }

( РУ1 РЬ2 \

2

Ш

РЬ2 + Р РУ\Ь2 РУ\У2 РУ2 + р

\ Ны НУ2

Н = Е+ р,

/

где р, р, V = (ь 1, у2)т — плотность, давление и вектор скорости; Е = р(е + V2/2) — полная энергия (е — внутренняя энергия); V2 = + у"2. Для определенности зададим уравнение состояния идеального газа в виде р = (7 — 1) ре.

При построении численного алгоритма, основанного на методе конечных объемов, воспользуемся идеологией расщепления [1, 7, 8]. Представим матрицу потоков Ш в виде суммы двух матриц

Ш = Ш + Ш2,

(

ру 2

РУ2 \ рУх Ь2

2

РЬ1У2 РУ2

\ ру2уг/2 ру2У2/2

( О

Шо

о \

о

Н Н 2

о

Н = Р + ре.

Первая из них содержит конвективные члены во всех уравнениях, а вторая Ш2 — оставшиеся члены (члены с давлением и внутренней энергией). Расщепление (3) аналогично расщеплению по физическим процессам для уравнений, записанных в дифференциальной форме [7, 8]. Тогда расщепленную систему уравнений Эйлера можно представить в виде системы двух уравнений

2 -Ц И(У + У адя = 0, 1 И(У + У Щ^Б = 0, (4)

V дУ V дУ

слабо аппроксимирующих исходную систему уравнений (1), (2). В скалярном виде эти уравнения принимают вид

2у %(У=0,

У дУ

1 г д г

2 + (рУ = 0, У дУ

1 д

2 д1р1'2(^ + (pv2v)dS = 0,

У дУ

1 г д

1 д

1 д

- ^(У = 0, 2} д1 , У

"У -¡Р^У + !(рег)(Б = 0,

У дУ

2У -¡Р^У + У (ре2) (Б = 0,

У дУ

1 г д

- I —Е(У + J (ру2/2^)(Б = 0 ^ 2] -¡¡(У = 0, 2] -¡Е(У + ] (Ре+Р^)(Б = 0.

дУ V V дУ

2 д

V

Здесь ег = (1, 0), е2 = (0,1) — базис декартовой системы координат.

1 Г дЕ Г

Расщепленное уравнение энергии - + (РV /2)v • ¿8 = 0 эквивалентно урав-

V дУ

1 Г др

нению для давления — -^rdV = 0. Действительно, уравнение энергии в дифференци-

2} т

V

1 дЕ 2 я

альном виде -—- = ¿2 т:—(р^2/2)г^ = 0 после его дифференцирования и исключения из

2 оь ох^

него плотности и скорости эквивалентно уравнению для давления в дифференциальной др [ др

— = 0 или соответственно в интегральной -¿тйУ = 0 форме. Таким образом, рас-дЪ ] дЪ

V

щепление по физическим процессам позволило упростить уравнения (4), фактически уменьшив число решаемых уравнений для каждой подсистемы на одно уравнение.

Исходные уравнения (1)-(3) будем решать в расчетной области V х (0,Т). Во временном интервале (0,Т) введем равномерную сетку с шагом т = 1/И, где N — число временных шагов, а в замкнутой расчетной области V — регулярную сетку, состоящую из ячеек ^ в виде четырехугольников (рис. 1). В каждой ячейки может быть введена параметризация по формуле Хг = х^г), где ^ изменяют свои значения в единичном квадрате.

Введем усреднение искомых функций по элементарной ячейке ш = и аппроксимируем интегральные операторы ^ ^ б!8 (I = 1, 2) на границах ячеек по следующим

дУ

формулам [18]:

= 1/ = / Щ ^ = А*(БЩ)<^ + а*,

у эч

А] (8ШгЬ = (8^)^+1/2 - (8ШгЬ-1/2. (5)

Здесь (8^)г±1/2^, (8^)г^±1/2 — аппроксимация потоков в (5) через грани Б ячейки , на которых задаются потоки. С учетом введенных обозначений конечно-объемная схема предиктор-корректор:

Рис. 1. Расчетная сетка

И тп+1/4 и™

и и , 0п+1/4 „ Ш--+ "ль = 0,

га ш

~ ~ , г^п+2/4 Ш--+

та 2Н

И»г+2/4 _ цИ-1/4

а

п«+1 — ТТ'-

Ши и + "Т1/2 = 0, (6)

где " = "ш + а интегральные выражения приближены по формулам (5), аппроксимирует исходные уравнения с порядком 0(тт + к2). Здесь к = у/Ш, а т = 2, если весовой множитель а = 0.5 + 0(т), иначе т =1. По построению конечно-объемная схема (6) консервативна для каждой ячейки и при а = 0 является нелинейной. Остановимся на ее реализации. На первом дробном шаге схемы (6) систему уравнений для отдельной ячейки представим в виде

Ип+1/4 тт п гпп+1/4 „п

Т Т . Лп+1/4 „ Р -Р „ ^

ш--+ " Л, =0, Ш-= 0, (7)

та 1П та

где

и = (р,рУ1,рУ2)т, = V ■ ТТ, 1/4 = Ai{Sv ■ И у*1/4 + Д, (5у ■ И )п+1/4

г,]

Введение расщепления (4) позволило представить уравнения (7) в виде расщепленных уравнений движения и неразрывности, а уравнение энергии свести к уравнению для давления, из которого следует, что п 1/4 = п, т. е. давление переносится с предыдущего временного слоя. Для решения оставшихся уравнений введем линеаризацию вектора потоков по формуле (V ■ И)п+1/4 ~ vn ■ Ип+1/4. Тогда вместо (7) получаем конечно-объемную схему

Ип+1/4 _ "ип

-+ Дг(Svn ■ Ип+1/4) + Д]( Svn ■ Ип+1/4) = 0.

а

Известно, что симметричная аппроксимация интегральных операторов может приводить к немонотонным решениям, поэтому на этапе предиктор может быть рассмотрена их аппроксимация с первым порядком несимметричными операторами (с учетом знака проекции вектора скорости на соответствующую нормаль к границе S ячейки Уг,]) по формуле

1 п ( $>Я)т±1/2 = 22рт±1/2 [( Я т± 1 + Ят) Т (vn)(Яm±1 - Ят)]

где т = I,а Я = vn ■ Ип+1/4. Схема (8) может быть реализована различными способами: векторными прогонками для каждой компоненты вектора, итерационными методами (например, с использованием ЬИ-факторизации) или введением расщепления по пространственным направлениям. Остановимся на последнем подходе. Представим схему (8) в виде расщепления по пространственным направлениям:

Ип+1/8 - Ип ~ , 1/8 -+ Д( Svn ■ ип+1/8) = 0,

а

"Ип+1/4 _ "Ип+1/8

-+ д л Svn ■ Ип+1/4) = 0.

а

Решение уравнений (9) на каждом дробном шаге может быть получено независимо для всех компонент скорости и плотности скалярными прогонками по каждому пространственному направлению по схемам первого (8) или второго (7) порядков аппроксимации потоков.

На втором дробном шаге схемы (6) система уравнений для отдельной ячейки может быть представлена в виде

ип+2/4 _ ип+1/4

та

+ п£2/4 = 0, *

рП+2/4 - рП+1/4

та

:ю)

откуда следует, что рп+2/4 = рп+1/4. Здесь

(

и

РЬ1 РУ2

2

Р V2

+ Р^

\

\1 - 1 "2 /

2

П2Н = ^ Дг№г), Л 1=1

р 0 0 р

Ьру1 Ьрь2

1

(1 - 1)'

(11)

Относительно вектора и схема (10) нелинейна. В качестве искомых функций выберем вектор ¥ = (рг>1, ру2,р)т и относительно него линеаризуем векторы "О и ЛЛ по формулам

ии П+2 / 4 = "П+1/4 + ДП(£П+2/4 — £П+1/4) + 0(т 2) = дп ■ £П+2/4 + 2),

лЛП+2/4 = Л+ Вп(Ип+2/4 - ип+1/4) + 0(т2) = Сп ■ Г+2/4 + 0(т2).

'12)

Здесь И = А ■ ¥, С = В ■ А, А = дИ/д£, В = <9Л¥/дИ, ру1 \ /10 0

П1

¥ = | рУ2 Р

А

01

0

С

У1 У2 1/(7 - 1)

00 0 0 2

Ьп1 Ьп2 Ь\

Так как конечно-объемная схема (10) аппроксимирует систему уравнений с первым порядком по времени, без потери точности матричные операторы А, В на дробных временных шагах могут быть выбраны на п-м шаге. Тогда в силу соотношений (10) из (12) получим линейную схему

р+2/4 - ¥п+1/4 "-+ (А )п П 21г 1

0.

'13)

та

Пусть С = А-1 ■ С. Тогда

(А-')пд Г4;

гп+2/4 = с п¥п+2/4) + (3 с пГ

и с учетом соотношений (12) схема (13) примет вид

(А-1)пПп+2/4 + О(Н), (14)

^п+2/4 _ £п+1 /4

ш-+ [Д(БСп ■ ¥п+2/4) + д.(^Сп ■ ¥п+2/4)] = 0.

та

Введем расщепление по пространственным направлениям:

¥п+3/8 _ ¥п+2/8

ш-

та

+ Дг(БСп ■ Р+2/8) = 0,

0

Ь

£п+1/2 _ £п+3/8

3

та

+ А3 ( ЬСп • fп+1/2) = 0.

:15)

Для повышения устойчивости алгоритма на дробных шагах будем использовать несимметричную аппроксимацию потоков через границы ячейки (8). На дробных шагах схема (15) реализуется скалярными прогонками. Покажем это на примере третьего дробного шага этапа предиктор, например Ь = (ьгПп1 + ^П^)^ > 0. В скалярной форме на слое (п + 3)/8 уравнения (15) имеют вид

7 п+3/8 _ ч ,п+1/4

ирП+1/4^-+5+1 Аг+1^Г/8 = 0,

а

_ п+1/4

п+1/4и 2 и2 . с д п+3/8 п

эа+1/4--+5+2 Аг+1/2'Рг+ = 0,

иР | ~+2— г+1/21:-

а

пп+3/8 _ 'п+1/4

Р Р по Л Г п+3/8 . п+3/8\ д п+3/8

где А±1/2Д.,- = ±(- Д.,-); Ь± = Ь±1/2^п4 (I = 1, 2) — нормаль вектора скорости к ребру Ь±-\_/2,]. Исключая значение давления, получим для него уравнение

сти к ребру Ь±]_/2,]. Исключая значение компонент скоростей у™+3/8 из уравнения для

1 + Н Ьг-1/2^ Аг-1/2 - Ь-^-1/2 ^Ь+А+ф-

-Ь-2А-1/2 -п+/4Ь+2Аг+1/2

п+3/8 = УгА = 1Р

где г = та/и; ¡Р = рп+1/4 - гЬг-1/2,зА—1/2(п1у'1+1/2 + п2Уп+1/2)^. Его решение находится трехточечной скалярной прогонкой, после чего явно определяются новые значения компонент скоростей по формулам

1 п+3/8 = 1 П+1/4 г Ь А п+3/8 „п+3/8 = 1п+1/4 г Ь А п+3/8

го = г/рп+1/4.

Напомним, что плотность вычислялась явно, ее значение переносилось с предыдущего временного слоя.

На четвертом дробном шаге схема (15) реализуется аналогично: скалярной прогонкой в направлении х2 находится значение давления рп+1/2, после чего явно вычисляют-

„ п+1/2 п+1/2 Тт

ся новые значения компонент скоростей 1 и 2 . На этом решение уравнений на этапе предиктор заканчивается.

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

ТТп+1 _ ип

и -^ + АДЯЖ )п+1/2 + А, (Б1У )п+1/2 = 0. (16)

Симметричная аппроксимация уравнений на этапе корректор может приводить к ос-цилляциям решения. Для их устранения вводится монотонизирующий оператор сглаживания по каждому пространственному направлению второго порядка малости по формулам

Дк9к = 9к+1/2 - дк-1/2 - Sgn(F) (£к+1/2Дк+1/29к - £к-1/2Дк-1/29к) . Здесь V — скорость; £к±1/2 = (¿к + £к±\)/2, Дк±1/29к = ±(#к±1 - 9к), |Дк+1/2 9к - Дк-1/29к1

£к

|Дк+1/2^к | + |Дк-1/2#к |

или 0, если |Дк+1/2^к | + | Дк-1/2 9к | = 0.

При £к = 1 схема (16) является противопотоковой схемой первого порядка. Предложенный алгоритм экономичен по числу операций на узел сетки, обладает свойством консервативности для каждой расчетной ячейки и, как показали расчеты, устойчив в широком диапазоне параметров сетки.

2. Алгоритм расщепления для уравнений Навье — Стокса

Рассмотрим уравнения Навье — Стокса вязкого сжимаемого теплопроводного газа при отсутствии внешних сил в виде системы интегральных законов сохранения [7, 17]:

^ [ и ¿V + [ ШйЪ = 0,

:17)

у

ЗУ

где И — вектор искомых функций, а Ш — матрица, составленная из столбцов потоков через границу элементарного объема, т. е.

и

( р \

РЬ1 РУ2

V Е

(

w

рУ1

ру"2 + р - а\ 2

рьхУ2 - о2

Р^2 \

ру\Ю2 - а2 ру% + р - а^ 22

(Е + Р)ы - (Е + Р)Ь2 - Я2 -Е VIа12

\ 1=1 1=1

:18)

а элементы тензора вязких напряжений а* и тепловые потоки Qj равны

дь1

(Ул

дх1

2р—--+ А div V, &2 = 2р—--+ А div V,

дх

2

<г1 = а1 = р

(дь1 дьЛ

(ей + а^

+ тН , ^ = кдТ/дх "2 дх1 1 j ■

9У1 ду2 2

Здесь div V = —--+ ——; А = р--р (р, р — коэффициенты первой и второй вязкости);

ох1 ох2 3

к — коэффициент теплопроводности. Для замыкания системы (17) заданы уравнения состояния в виде

Р = Р(Р,е) = (1 - 1)pe, е = е(Т ) = Т (19)

и зависимости коэффициентов вязкости и теплопроводности от температуры [7].

При вычислении потоков на гранях ячейки для уравнений Навье — Стокса (17) и (19), наряду с газодинамическими переменными, необходимо аппроксимировать члены а* и Я 7, содержащие производные от скорости и температуры на границах ячейки, не совпадающих с осями декартовых координат. Поэтому введем невырожденное преобразование системы координат

41 = 41 (Х1,Х2), I = 1, 2, отображающее ячейку в единичный квадрат. Тогда

д дxj

д

д

ддг 7 дд2

дх7

Коэффициенты г1, и производные д/д^ будем аппроксимировать симметричными опе-

раторами:

/ % ддЛ

I дх\ дх2)

Л^

г+1/2

И

г+1/2

(ее у

^ )г+1/2, , дх2)

После преобразований матрица Ш примет вид

г+1/2

Л®

И

г+1/2

(БХ1, БХ2 )ъ+1/2.

(

W

рУ1

ру"2 + р - а\ 2

рЬхЬ" - и2

рУ2 рЬхЬ" - (т\ ру" + р - т"

(Е + р)ьх - Яг - Е Чт[ (Е + р)ь" - Я2 - Е ^т" 1=1 1=1

где

1 1 л 1ч ду1 2 2 л 2 \ ду2 2 1 9У2

о\ = + А^1) —, т2 - ' ^ -2----1

дд1

(2^ + ^2)—, у = рг2—.

дц

ддг

То

2

од2

Введем по аналогии с расщеплением уравнений Эйлера расщепление по физическим процессам:

Ш = Щ + W2,

( ру\ РУ2 \

ру\ 21

РЩ -

2

рЬ1Ь2 - У2

РУ2 РУ1У2 - Т2 рь% - У2

2

V2 т 2 у2 т 2

р—У1 - Я1 - Е V1т[ р—Ъ2 - Я2 - Е VIт2 2 1=1 2

2=

0

1=1

о

р о

0 р \ (ер + р)г>1 (ер + р)У2 )

Матрица включает конвективные и диссипативные члены (вязкие члены в уравнениях движения и члены с теплопроводностью в уравнении энергии), а Ш2 — члены с давлением, как и при расщеплении уравнений Эйлера (2).

1

2

2

7

7

В [7] в исходных уравнениях при расщеплении по физическим процессам вводилось расщепление векторов потоков на три вида: инерционные или конвективные, обусловленные градиентом давления и обусловленные диссипативными эффектами вязкости и теплопроводности. Там же отмечалось, что введение расщепления по физическим процессам и пространственным направлениям для уравнений Навье — Стокса, записанных

в дифференциальной форме, не позволяет свести их к решению одномерных задач, так

к

как элементы тензора напряжений а к содержат производные по каждому направлению. Поэтому вязкие и диссипативные члены сохраним в операторе Ш1. Тогда уравнения (18) и (19) могут быть представлены в виде расщепленной системы, слабо аппроксимирующей уравнения Навье — Стокса [1, 7]:

1 ж/ш"+/

У дУ

Ш^Б = 0,

21/ш"+/

V дУ

ВДБ = 0.

(20)

Для построения экономичных схем, реализуемых скалярными прогонками, в операторе Ш1 сохраним лишь часть вязких членов, перенеся остальные в вектор правых частей, т. е. перепишем уравнения (20) в виде

^ I И ¿V + I

2Ы} У

V дУ

Л¥ 1 -¿Б

Е • (!Б,

(21)

дУ

[и(IV + I 2Ы} У

V дУ

= 0,

(22)

где

Л¥ 1

рУ1

оу2 Р19Ь1 д у2

РУ1У2 - й —

рщ

д 1 дТ

2 д 1

рУ 1У2 - К* — 1 д 2

0,2 Р2^

д Т

\

ру1- у2/2 - Кх— рУ2 • У2/2 - К2-ТГ-дЯ! дд2

2

Е = Ш- Л¥

1;

£ = (2^ + Хг\), ^ = (2^2 + Хг1), С2

к

В качестве искомых функций выберем вектор ¥ = (р, ру1, ру2, рТ)т или ¥ = (р, ру1, ру2, р)т. Тогда с учетом (19) получим Е = рТ + ру2/2 и система уравнений (22) может быть переписана в виде

^ /и(IV +/Л¥ •¿Б = Еь (23)

(1 = 1, 2).

2 д

Л¥

РУ1

У

РУ 2

дУ

2 ^дъ 1 .2дУ1

дя1 дд2

1 д 2 2 2

ру IV2 - й ТГ" РЩ - & -7Г-

2 д 1 2 2 д 2

К1

дТ д д1

К2

дТ дд2

Е = /Е-Л8 -1 ж/

дУ У

0 0 0

V РЪ2/2 )

¿V.

Аппроксимируем интегральные операторы у сеточными операторами:

дУ

Пш = У ^ (1Б = )+ Л^^)г,7 = ПцН + ПцН,

дУ

где Лт(8И^)т = (8И^)т+1/2 - (8Жг)т-1/2, а (8Жг)т±1/2 — потоки через грани 5 ячейки по направлению т. Тогда схема предиктор-корректор

Сп+1/8 _

йо-

та

+ П

п+1/8 11Н

0, ш-

£п+2/8 _ £п+1/8

та

+ П

п+2/8 12Н

Ш

с п+3/8 _ с п+2/8 с п+1/2 _ с п+3/8

-+ = о, и-+ пп+У2 = 0,

та

та

ь22к

ип+1 — ип

ш-

+ пп+1/2

0

(24)

с расщеплением операторов по физическим процессам и пространственным направлениям аппроксимирует уравнения Навье — Стокса (17) с порядком 0(т2 + к2) при а = 0.5 + О(т). Так как на этапе предиктор уравнения аппроксимируются с первым порядком по времени, вектор Е1 без потери порядка аппроксимации исключается из первого дробного шага.

Остановимся на реализации схемы (24). На первом дробном шаге после линеаризации оператора Пп+У8 подобно (7) система уравнений может быть представлена в виде

с п+1/8 с п

ш-+ Лг(Д1:Т+1/8) = 0,

та

(25)

где с

( р \

РЬ1 РУ2

\ рт )

Я1

0

1

0 и - Лг+1/2 —

V

0 0

0 0

11

^1 - ^2$1Лг+1/2 —

0

0 0

-г+1/2 — !

и = У1Б1 + У2Б2; Лг±1/2!г7 = ±(¡1+1,7 - Л^ Лг= ¡г+1/2,7 - !г-1/2,7\ $I — проекция границы в ячейки на ось х\. Как следует из вида оператора К1, система уравнений (25) для каждой компоненты вектора с может быть представлена в виде

та

1+— Л, ш

(

и - г5*1 Лг+1/2 —

Г+1/4 = Г, и = (* 1, Ь, и, 0), г = (0, £,«)

0

0

т. е. все уравнения решаются независимо друг от друга скалярными трехточечными прогонками. На втором дробном шаге схемы (24) уравнения

сп+2/4 _ сп+1/4

и-+ Аг(П2Г+2/4) = 0

а

решаются скалярными прогонками по второму пространственному направлению подобно первому дробному шагу. На третьем и четвертом дробных шагах на этапе предиктор реализация схем аналогична схеме (15), так как матрица Ш2 совпадает для уравнений Эйлера и Навье — Стокса.

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

Предложенный алгоритм обобщается на пространственные уравнения по аналогии с двумерным случаем: вводится расщепление по физическим процессам и для расщепленных уравнений — по пространственным направлениям. Особенностью алгоритма является возможность реализации уравнений на дробных шагах на этапе предиктор скалярными прогонками, т. е. он экономичен по числу операций на расчетную ячейку, в отличие от схем с расщеплением по пространственным направлениям, которые реализуются векторными прогонками. На регулярной равномерной сетке данный алгоритм эквивалентен разностной схеме с расщеплением по физическим процессам и пространственным направлениям [8] и, как следствие, обладает теми же свойствами — он безусловно устойчив в двумерном случае и условно устойчив в трехмерном, как и все схемы, приведенные к каноническому виду [7]. Наиболее эффективным представляется использование данного типа алгоритмов решения методом установления, когда итерационный параметр может быть выбран лишь из условия наиболее быстрой сходимости.

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

ип+1/4 ип

Т Т , г^п+1/4 „

и--+ = 0,

та 1П

Ип+2/4 _ ип+1/4

' Г1,2Н

и----+ пп+2/4 = 0,

а

ип+1 - ип .,+1/2 и-+ ^п+ = 0.

т п

На каждом дробном шаге этапа предиктор алгоритм можно реализовать с использованием, например, ЬИ-факторизации по схемам бегущего счета для каждого расщепленного уравнения [9-17]. Такой подход не требует обращения матриц, так как матрица П^ диагональная, а система уравнений для матрицы П2^ после исключения компонент скоростей сводится к неявному решению уравнения для давления. Его решение может быть получено также методом ЬИ-факторизации, а значения скоростей затем могут быть восстановлены по явным формулам.

3. Примеры численных расчетов 3.1. Тестирование алгоритма

Алгоритм протестирован на двух задачах, имеющих точное решение: задаче о распаде произвольного разрыва и задаче о течении в канале переменного сечения [19, 21]. В первой задаче расчетная область имела длину Ь = 10. В начальный момент слева и справа от координаты х = 0 задавались различные значения плотности и давления

(р,у,р)тх<0 = (1.0,0,1.0), (Р,у,р)тх>0 = (0.125,0, 0.1).

Решение уравнений газовой динамики находилось при £ > 0. Расчеты проводились в приближении одномерных уравнений Эйлера по конечно-объемной схеме (8) с первым и вторым порядками аппроксимации до момента времени £ = 2.5 [8, 19, 21]. Шаги сетки к = Ь/(Н - 1) варьировались при различном числе узлов N = 101, 201, 401.

На рис. 2 приведено распределение продольной компоненты скорости и плотности, полученное в расчетах по схемам первого и второго порядков аппроксимации (крестики) и из точного решения (сплошная линия) на сетке N = 401. Более детально различие в расчетах по схемам первого (*) и второго (+) порядков аппроксимации видно на рис. 2, в для части расчетной области. Таким образом, предложенный алгоритм решения уравнений Эйлера верно передает скорости фронтов, размазывает ударную волну в схеме первого порядка на 8-9, а в схеме второго порядка на 5-6 ячеек, а контактный разрыв — на большее число ячеек.

а б

в

Р

1.5 2.5 3.5 4.5 х

Рис. 2. Распределение продольной компоненты скорости (а) и плотности (б), полученное в расчетах по схемам первого и второго порядков аппроксимации (крестики) и из точного решения (сплошная кривая)

Во второй задаче исследовалось течение в сужающемся канале (рис. 3). На левой границе расчетной области задавался равномерный сверхзвуковой поток с параметрами

р = 1, V! = 1, = 0, р = 1/[(7 - 1)М2], Т = 1/[7(7 - 1)М2]

на верхней границе Г2 — условия непротекания V • п|г = 0 при течении невязкого газа или условия прилипания и тепловой изоляции дТ/дп = 0 для вязкого газа, на нижней границе — условия симметрии, на выходной границе Г3 — мягкие условия д¥/дп = 0, где ¥ = (р,у1 ,у2 ,р). При £ = 0 внутри области задавались значения набегающего потока, как на входной границе Г1. Все приведенные ниже результаты численных расчетов получены по схеме предиктор-корректор второго порядка аппроксимации (12), (18), (19) на различных сетках, что позволило оценить точность алгоритма и его эффективность.

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

е = тах

£П+Т _

Тт ¥7

<Н2, т ^ Н,

где Т « 30 ^ 40, так как в расчетах по схеме второго порядка со сглаживанием (19) классический критерий установления (Т = 1) не всегда выполнялся из-за появления осцилляций невязок решения в схемах с симметричной аппроксимацией первых производных.

Расчеты в сужающемся канале выполнены при различных параметрах течения и углах наклона клина. В первой серии расчетов исследовался регулярный режим отражения. На рис. 4 приведено распределение плотности и температуры в установившемся

Рис. 3. Расчетная сетка в канале

0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

0.5

1.5

2.5

3.5

4.5

Рис. 4. Распределение плотности р (а) и температуры Т (б) в канале

Рис. 5. Распределение плотности: точное решение — сплошная кривая; численные решения — штриховая кривая — сетка (1) 300 х 150, штрихпунктирная — сетка (2) 200 х 100 ячеек

0.5

1.5

2.5

3.5 4 4.5

5.5

-I-

. I.

-и

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1.1

Рис. 6. Распределение продольной компоненты скорости и (а) и плотности р (б) в канале

решении при М =3 и угле наклона tg(а) = 0.28, полученное на сетке 300 х 150 ячеек при т = 0.01. Угол наклона падающей и отраженной волн для данных параметров равен 3 = 32° и совпадает с точным решением [22-24]. Точное и численное решения в продольном сечении х2 = 0.21 (для 0.7 < х < 1.7) приведены на рис. 5 при М = 3, tg(a) = 0.28 на двух сетках. Типичное число итераций до сходимости для этих режимов отражения составило 7000. Из рис. 5 можно видеть полное совпадение значений параметров и положений падающей и отраженной волн, численное решение размазано на 4-5 ячеек. Расчеты по схеме первого порядка приводили к более сильному размазыванию ударных волн.

Изменение угла наклона образуюшей стенки канала может приводить к изменению режима течения. На рис. 6 приведено распределение продольной компоненты скорости и плотности при угле наклона tg(a) = 0.36 на сетке 200 х 150. Этот режим течения соответствует маховской конфигурации отражения.

Сравнение параметров течения, полученных в настоящем расчете, с данными расчетов [22, 23] говорит о практическом совпадении картины течения и углов наклона падающих и отраженных скачков, а также волн разрежения.

3.2. О режиме двойного решения

Из теоретических и экспериментальных данных [22, 24] известно, что при больших значениях числа Маха (М > 2.2) в одной и той же области изменения углов клина и, как следствие, углов наклона падающие волны при отражении от клина могут иметь две конфигурации — регулярную и маховскую. Так, в диапазоне а^ < а < а^ при углах меньше а^ существует только регулярная конфигурация, при углах больше а в — только маховская, а внутри этого интервала появляется область существования двух волновых конфигураций. Возникает вопрос: от чего зависит появление той или иной конфигурации в области двойного решения? Обычно полагают, что в этом случае регулярная конфигурация неустойчива к возмущениям, которые могут возникать в течении. Этот эффект и был проверен в расчетах течений для М = 3 при различных формах возмущения начальных данных. Моделировалось течение вязкого газа при Ив = 103. Расчеты выполнялись на различных сетках. На теле задавались условия прилипания и тепловой изоляции, остальные условия подобны рассмотренным выше. Известно [22-24], что для данного числа Маха двойное решение находится в интервале 37.3 < а < 39.5°. Угол наклона задавался таким образом, чтобы ударная волна падала на нижнюю грань симметрии под углом в а = 38°. В первом расчете в начальный момент внутри расчетной области задавались те же газодинамические параметры, как и на входной границе, что соответствовало случаю мгновенного помещения тела в сверхзвуковой поток. В итоге реализовалась регулярная конфигурация ударной волны, отраженной от плоскости симметрии (рис. 7).

Следующий расчет проводился с теми же граничными условиями и углом наклона клина, что и в первом случае, но с другими начальными данными. В начальный момент времени внутренняя область заполнена покоящимся газом с давлением р^ = 0.03р^. Таким образом, имеется разрыв газодинамических параметров на границах области. На рис. 8 представлено распределение продольной компоненты скорости и плотности при установлении. Из-за присутствия эффектов вязкости фронты ударных волн немного смазаны, но тройная маховская конфигурация сохраняется. Расчеты проведены на сетках из 200 х 100 и 300 х 150 ячеек. Картины течений для обеих сеток совпадают. Так как в начальный момент времени на входе в канал существует разрыв газодинамиче-

-0.2 0 0.2 0 А 0.6 0.8 1

1.2

0.5

1.5

2.5

3.5

4.5

5.5

Рис. 7. Распределение продольной компоненты скорости и (а) и плотности р (б) падающей волны регулярной конфигурации. Сетка 200 х 100

Рис. 8. Распределение продольной компоненты скорости и (а) и плотности р (б) волны махов-ской конфигурации. Сетка 300 х 150

ских параметров, в течение времени он будет генерировать две ударные волны внутрь канала. Отразившись от клина, первая ударная волна приобретает маховскую конфигурацию, как показано на рис. 8. Последующие расчеты подтвердили, что образование маховского отражения с помощью разрыва газодинамических параметров в области двойного решения зависит от относительного положения клина и плоскости симметрии, а также начального значения давления во внутренней области. Результаты проведенных расчетов качественно совпадают с данными, полученными экспериментально [22, 23]. Наиболее быстрая сходимость результатов расчетов получена при г ~ Н3/2, что дает возможность для последующей оценки эффективности алгоритма.

Замечание 2. Возможности распараллеливания предложенного алгоритма исследованы лишь на этапе предиктор, так как на этапе корректор алгоритм явный и его распараллеливание производится стандартным образом. Введенное специальное расщепление на этапе предиктор позволило свести решение системы уравнений к их независимому решению для каждого расщепленного уравнения (оно может быть получено скалярными прогонками по каждому пространственному направлению). При распараллеливании на каждом этапе расщепления расчетная область разделялась на равные области по направлениям и каждое уравнение решалось независимо от других. В данной работе не был проведен анализ эффективности распараллеливания алгоритма в полном объеме, это планируется в дальнейшем. Выполнена только серия расчетов с разным числом вычислительных ядер, которая показала близкое к линейному ускорение на дробных шагах, зависящее от числа вычислительных узлов и элементарных объемов. Так, например, на сетке из 250 х 100 ячеек время расчета при использовании одного вычислительного ядра в 1.8 раза превышает время расчета при использовании двух ядер, но при сравнении времени вычисления на двух и четырех вычислительных ядрах относительное ускорение составило порядка 1.7. Таким образом, расчеты показали возможность распараллеливания предложенного алгоритма МКО.

Заключение

Для численного решения уравнений Эйлера и Навье — Стокса, записанных в интегральной форме, проведено обобщение алгоритма типа предиктор-корректор на метод конеч-

ных объемов с использованием расщепления операторов по физическим процессам и пространственным направлениям в случае регулярных структурированных сеток. Реализация алгоритма сводится к скалярным прогонкам по каждому направлению на этапе предиктор и явным вычислениям на этапе корректор.

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

Благодарности. Работа выполнена при финансовой поддержке РФФИ (грант № 1401-00191).

Список литературы / References

[1] Яненко Н.Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука. Сиб. отд-ние, 1967. 197 с.

Yanenko, N.N. The method of fractional steps the solution of problems of mathematical physics in several variables. Berlin: Springer-Verlag, 1971. 158 p.

[2] Годунов С.К., Забродин А.В., Иванов М.Я., Крайко А.Н., Прокопов Г.П. Численное решение многомерных задач газовой динамики. М.: Наука, 1976. 400 с. Godunov, S.K., Zabrodin, A.V., Ivanov, M.Y., Krayko, A.N., Prokopov, G.P.

Numerical solution of multi-dimensional problems in gas dynamics. Moscow: Nauka, 1976. 400 p. (in Russ.)

[3] Флетчер К. Вычислительные методы в динамике жидкостей: В 2 т. М.: Мир, 1991. 504 с. Fletcher, C.A.J. Computational techniques for fluid dynamics. Berlin: Springer-Verlag, 1988. 409 p.

[4] Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. М.: Физматлит, 2001. 608 с. Kulikovskii, A.G., Pogorelov, N.V., Semenov, A.Y. Mathematical problems in the numerical solution of hyperbolic systems. Moscow: Fizmatlit, 2001. 608 p. (in Russ.)

[5] Yamomoto, S., Daiguji, H. Higher-order accurate upwind schemes for solving the compressible Euler and Navier —Stokes equations // Computer and Fluids. 1993. Vol. 22. P. 259-270.

[6] Vos, J.B., Rizzi, A., Darrac, D., Hirschel, E.H. Navier — Stokes solvers in European aircraft design // Progress in Aerospace Sciences. 2002. Vol. 38. P. 601-697.

[7] Ковеня В.М., Яненко Н.Н. Метод расщепления в задачах газовой динамики. Новосибирск: Наука. Сиб. отд-ние, 1981. 304 с.

Kovenya, V.M., Yanenko, N.N. The splitting method in problems of gas dynamics. Novosibirsk: Nauka. Sibirskoe otdelenie, 1981. 304 p. (in Russ.)

[8] Ковеня В.М. Алгоритмы расщепления при решении многомерных задач аэрогидродинамики. Новосибирск: Изд-во СО РАН, 2014. 280 с.

Kovenya, V.M. Splitting algorithms for solving multidimensional problems of aerodynamics. Novosibirsk: Izdatel'stvo SO RAN, 2014. 280 p. (in Russ.)

[9] Le Veque, R.J. Finite volume methods for hyperbolic problems. Cambridge: Cambridge University Press, 2002. 580 p.

[10] Meister, А., Sonar, Th. Finite-volume schemes for compressible fluid flow // Surveys on Mathematics for Industry. 1998. Vol. 8. P. 1-36.

[11] Toro E.F. Riemann solvers and numerical methods for fluid dynamics. A practical introduction. 2nd Edition. Berlin: Springer-Verlag, 2009. 724 p.

[12] Morton, K.W., Sonar, Th. Finite volume methods for hyperbolic conservation laws // Acta Numerica. 2007. Vol. 16. P. 155-238.

[13] Titarev, V.A., Toro, E.F. ADER: Arbitrary high order Godunov approach // Journal Scientic Computing. 2002. Vol. 17, No. 4. P. 609-618.

[14] Remaki, L., Hassan, O., Morgan, K. Aerodynamic computations using a finite volume method with an HLLC numerical flux function // Mathematical Modelling of Natural Phenomena. 2010. Vol. 10, No. 5. P. 1-20.

[15] Bijl, H., Carpenter, M.H., Vatsa, V.N., Kennedy, C.A. Implicit time integration schemes for the unsteady compressible Navier —Stokes equations: Laminar flow // Journal of Computational Physics. 2002. Vol. 179. P. 313-329.

[16] Ketcheson, D.I., Macdonald, C.B., Gottlieb, S. Optimal implicit strong stability preserving Runge —Kutta methods // Applied Numerical Mathematics. 2009. Vol. 59. P. 373-392.

[17] Чирков Д.В., Черный С.Г. Неявный метод численного моделирования пространственных течений вязкого газа // Вычисл. технологии. 2003. Т. 8, № 1. С. 66-83.

Chirkov, D.V., Cherny, S.G. Implicit method of numerical modeling of spatial fluids of viscous gas // Computational Technologies. 2003. Vol. 8, No. 1. P. 66-83. (in Russ.)

[18] Черный С.Г., Чирков Д.В., Лапин В.Н., Скороспелов В.А., Шаров С.В. Численное моделирование течений в турбомашинах. Новосибирск: Наука, 2006. 202 с. Cherny, S.G., Chirkov, D.V., Lapin, V.N., Skorospelov, V.A., Sharov, S.V. Numerical simulation of flows in turbomachinery. Novosibirsk: Nauka, 2006. 202 p. (in Russ.)

[19] Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1978. 736 с. Loitsiansky, L.G. Fluid mechanics. Moscow: Nauka, 1978. 736 p. (in Russ.)

[20] Ковеня В.М., Слюняев А.Ю. Алгоритмы расщепления при решении уравнений На-вье —Стокса // Журнал вычисл. математики и мат. физики. 2009. Т. 49, № 4. C. 700-714. Kovenya, V.M., Slyunyaev, A.Y. Splitting algorithms for solving Navier — Stokes equations // Computational Mathematic and Mathematical Physics. 2009. Vol. 49, No. 4. P. 700-714. (in Russ.)

[21] Лебедев А.С., Черный С.Г. Практикум по численному решению уравнений в частных производных: Учеб. пособие. Новосибирск: НГУ, 2000. 136 p.

Lebedev, A.S., Cherny, S.G. Workshop on the numerical solution of partial differential equations: Tutorial. Novosibirsk: NGU, 2000. 136 p. (in Russ.)

[22] Ivanov, M.S., Bendor, G., Elperin, T., Kudryavtsev, A.N., Khotyanovsky, D.V.

Flow-Mach-number-variation induced hysteresis in steady flow shock wave reflections // American Institute of Aeronautics and Astronautics Journal. 2001. Vol. 39, No. 5. P. 972-974.

[23] Ivanov, M.S., Vandromme, D., Fomin, V.M., Kudryavtsev, A.N., Hadjadj, A., Khotyanovsky, D.V. Transition between regular and Mach reflection of shock waves: new numerical and experimental results // Shock Waves. 2001. Vol. 11, No. 3. P. 197-207.

[24] Von Neumann, J. Oblique reflection of shock waves. Oxford: Pergamon Press, 1963. Vol. 6. P. 238-299.

Поступила в 'редакцию 6 июля 2015 г., с доработки — 18 сентября 2015 г.

84

B. M. KoBeH3,n. B. Ba6HHu,eB

Splitting algorithms in finite volume method

Kovenya, Viktor M.1,2, Babintsev, Pavel V.2'*

institute of Computational Technologies SB RAS, Novosibirsk, 630090, Russia 2Novosibirsk State University, Novosibirsk, 630090, Russia * Corresponding author: Babintsev, Pavel V., e-mail: [email protected]

We propose a class of the finite volume predictor-corrector algorithms for the numerical solution of the Euler and Navier —Stokes equations presented in the integral form in the case of a compressible heat-conducting gas. In the predictor step we have introduced a special form of splitting of the governing equations. In the approximation of the Euler equations, the vector flows across the borders of a cell is split into two vectors: the first vector contains only the convective terms, and the second one contains the terms with pressure and internal energy, which is an analogue of splitting into physical processes in the differential form. This kind of splitting is equivalent to representation of the original equations in the form of two systems that approximate the original Euler equation. Upon the introduction of this splitting, we are able to reduce the solution of equations in fractional steps to effective scalar sweeps. Approximation of the Navier — Stokes was introduced by analogy with the splitting of the Euler equations: flux vector is the sum of two vectors, the first of which includes convective and dissipative terms of the equations, and the second term includes pressure in the same way as it was done for the splitting of the Euler equations. In the corrector step we deal with the complete system of equations in the conservative form using the explicit scheme. The algorithm is tested on the two problems with the exact solution: the Riemann problem and the flow in a channel with variable cross-section. For the flow of gas in the convergent channel with regular and irregular reflection effects we have numerically confirmed the existence of the set of dual solutions, depending on the initial data.

Keywords: Euler and Navier — Stokes equations, splitting method, finite-volume method, regular and Mach reflection of waves.

Acknowledgements. This work was financially supported by RFBR (grant No. 1401-00191).

Received 6 July 2015 Received in revised form 18 September 2015

© ICT SB RAS, 2015