Эффективность неявных схем Рунге-Кутты четвертого порядка в задачах газовой динамики Текст научной статьи по специальности «Физика»

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

Том 8, № 6, 2003

ЭФФЕКТИВНОСТЬ НЕЯВНЫХ СХЕМ РУНГЕ - КУТТЫ ЧЕТВЕРТОГО ПОРЯДКА В ЗАДАЧАХ ГАЗОВОЙ ДИНАМИКИ*

В. И. Пинчуков

Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: pinchvi@net.ict.nsc.ru

Two types of implicit fourth order Runge — Kutta schemes are constructed for both multidimensional transfer equation with diffusion and for compressible flow equations. The sets of coefficients providing absolute stability of schemes are found. Adaptive artificial diffusion is used to provide convergence in time and to damp oscillations near shocks. Two types of absolutely stable fourth order schemes are compared. The results of test calculations for the case of compressible flows illustrating the schemes effectiveness are presented.

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

Первые схемы повышенного (третьего) порядка типа Рунге — Кутты предложены в работах Русанова [1], Бернстейна — Мирина [2], Кутлера — Ломекса [3]. Современные явные адаптивные ENO-схемы третьего и более высоких порядков сконструированы, например, в [4, 5] (см. также обзор [6]). В некоторых задачах предпочтительны неявные методы. В [7] предложены компактные неявные схемы высоких порядков. Абсолютно устойчивые схемы Рунге — Кутты третьего порядка развиты в [6, 8, 9]. В [10, 11] предложены неявные схемы четвертого порядка двух типов, доказана их абсолютная устойчивость. Здесь сравнивается работоспособность этих схем в задачах динамики сжимаемого газа.

""Исследования, произведенные В. И. Пинчуковым и описанные в этой статье, стали возможными отчасти благодаря поддержке Американского фонда гражданских исследований и развития для независимых государств, образованных на территории бывшего Советского Союза (CRDF) (грант № RM1-2324-NO-02).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2003.

1. Схемы для уравнения переноса с диффузией

Рассмотрим т-мерное уравнение

О т п

ди ^^ ди

а1 дж, г=1 1

д*

т п2

, д2и ¿^«/т^, а/

1=1

' джг2:

еопэ^ а = еопэ^ I = 1,...,т.

(1)

Введем разностный оператор конвекции

Д1м = ^ а, ДМ 1

г=1

г "г

и

т о

£ а, дХ/ + О ,

г=1 г

Л, = тахт=1 (Дж,), разностный оператор диффузии

т т д2 и (

уи = £ а/ДГД++и = £ дХи + О (Ь2 г=1 г=1 1

(2)

(3)

и вспомогательные операторы, используемые при построении схемы,

Я2

£а2Д+Д,

^4

Еа4 (Д++ДГЛ2

г=1

г=1

(4)

Здесь и далее используются обозначения 8±и = ± [и (ж1,..., ж, ± Дж,,..., жт) — и (ж1,..., жт)], Д±и = #±и/Дж,.

Анализ устойчивости различных вариантов аппроксимации позволил построить нижеследующую четырехшаговую схему:

и(1) - ип

т/2

+ ^ип = Уип

(5)

(1 — Лт2^)

и(2) - и(1)

т/2

т

— V — т>УЯ2

и(2) — ип

+ Д1 (и(1) — ип) = 0;

(6)

и(3) — ип

+ ^и(2) = Уип

(7)

(1 + Ст4^4 — ПтУ + т 3^УР2)

ип+1 - ип

+ 2пУ (и(1) — ип) +

ип и(1) и(2) и(3)

Т + 1Г + - + -

Уи(2).

В результате анализа устойчивости [10] методом Фурье доказана

Теорема 1. Схема (5)-(8) при V = 0 абсолютно устойчива, если

Л = в/6, С = в3/72, в = т 250/243.

(8)

(9)

т

6

Для обобщения схемы на случай ненулевой диффузии далее с помощью ЭВМ был проведен поиск значений параметров ^ и р (параметр п принимал значения 1/2 и 1/3, найденные в процессе рассмотрения асимптотических случаев), гарантирующих ограниченность по модулю коэффициента перехода схемы единицей при любом т. В результате перебора нескольких десятков вариантов получены следующие наборы параметров, обеспечивающие абсолютную устойчивость схемы:

^ = в/12, р = в/22, п = 1/3

(10)

и

^ = в/15, р = в/18, П =1/2-

(11)

Выше приведены четырехшаговые схемы Рунге — Кутты. Позже были найдены [11] трехшаговые схемы, также обладающие абсолютной устойчивостью для многомерного уравнения переноса:

тт

1 + 4 д. - 2 у

и(1) - ип

т/2

+ Д.ип = Уип;

(12)

(1 + ^т4Д4 - 2V) ц(2) т иП + Д.и(1) = Уип

(13)

1 ^ 4п ттЛ ип+. - ип п /ип 2и(1) и(2)', тгп 1 + Ст4Д4 - 2 V) - + Д. I — + + — 1 = Уига

т

6

3

6

(14)

Теорема 2. [11]. Схема (12)-(14) устойчива при любом шаге по времени т, если

С = в3/384 = (т250/243)3 /384.

(15)

2. Схемы для уравнений сжимаемого газа

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

Ц + Fx + Су = 0-Обобщение схемы (5)-(8) имеет вид

(Ц. - Цк) / (т/2) + дп = {у: + V;) ип;

(16)

1 - Ат2Д2х - 2Ух*п + 2т3^у;*пдпх 1 - Ат2Д2у - -у;п + 2т3^у;*пд

2

п 2у

и(2) - ип

т/2

+2АтДп (и(1) - Цп) + д11) = (Ухп + V")

(1) 1

у ) игк1

(17)

U(3) - Un

+ r12) = (Vj2) + Vy(2) ) Uffc;

(18)

1+ Ст4RÎÏ - 2Vj*2) + т32^(**2)r2ï

x

x

1 + Ст 4r42) - Tv(*2) + т 32^уу(**2)я2У)

U^1 _ иn

*2) (U(1) - Un) + ( ^ + ^ + ^ + R

(1)

6

Ri.

3

¡>(2)

(3)'

3

6

(v(2) + vy(2)) u (2).

(19)

Здесь Л, С, — параметры, заданные в формулах (9)-(11); R2x = А+ (а2) А-, а2 c + |u|; c — скорость звука; величины с индексом y вычисляются аналогично;

Ri = Л2^ + Л;Gifc, Л2 = АХ (1 - №/6), Л; = А; (1 - ¿+¿"/6) ;

(20)

Т/Х, ТУ — нелинейные адаптивные операторы искусственной диффузии, переключающиеся в зависимости от гладкости решения (диффузия пятого порядка на "гладких" решениях и первого возле ударных волн):

VxU

2Uifc

А-

Д2х) _ ¿+л-^(6х) ¿+Л-

"î+1/2fc ¿ж ¿2 ci+1/2fc¿ж ¿2

¿+U

ifc )

e(62) = а276, б(22) = а2 (1 - y) 16b,

(21)

где, как правило, используется значение Ь = 1/32, 0 < 7^ < 1 — индикаторы локальной гладкости, вычисляются по давлению на основе формул

Tifc

1 - 1 -

max

о)]:

(22)

[min (agjfc, 1) - c] / (1 - c) , gifc = (z z+ + ep|) /{max (z+) , (z ) + 10 9}

= ¿-Ptfc, z+ = 5+Pi.

jfc •

Константы а и с позволяют регулировать "чувствительность" величин 7^ к значениям разностной функции Р^ в узлах сетки, параметр е предназначен для предотвращения понижения порядка диффузии на гладких экстремумах или участках почти постоянного решения. По результатам пробных расчетов были выбраны следующие значения эмпирических величин: с = 0.4, а = 1.3, е = 10-3. Близкие формулы для коэффициентов 7 использовались при решении двумерного уравнения переноса схемой четвертого порядка аппроксимации [10, 11] и уравнений газовой динамики схемами различных порядков [12].

Как видно, искусственная диффузия носит скалярный характер, что позволяет удобно представлять ее в стабилизирующем операторе. Поскольку представленные выше диффузионные операторы являются семиточечными, в стабилизирующих операторах схемы (16)-(19) диффузионные члены приходится урезать для приведения их к пятиточечному или трехточечному, если они умножаются дополнительно на разностные операторы, виду:

v;Uifc

А-

(22) + 4^(б2) 5+5-i+1/2fc + i+1 /2fc¿2

¿+Uifc,

z

V** и,

гк

А:

_(2х) ~г+1/2к:

е;, , + 16е.

(6х) г+1/2к

и,

гк •

Адаптивная диффузия (21) переключается в зависимости от локальных индикаторов "гладкости", она увеличивает диссипацию возле ударных волн. Структура формулы (21) такова, что обеспечивается наличие стоков некоторой "разностной энергии" [14]. Опыт расчетов задач аэродинамики говорит о том, что это способствует подавлению осцилля-ций возле скачков и улучшению сходимости по времени. Формула (22) обусловливает для слабых разрывов квадратичную зависимость коэффициента вязкости низкого порядка от градиента давления, в отличие от ранее использованных формул [10-12], где имела место линейная зависимость. В результате улучшилось разрешение слабых разрывов. Такая же квадратичная зависимость имеет место в составной искусственной вязкости [13].

Операторы Я4х, Я4у также строятся по адаптивным формулам:

Я

4х —

А:

А+ «к)4 7гкА- - (яХ+1/2к) 4 (1 - 7г+1/2к) 4/АХ

А+.

(23)

Формулы (20) аппроксимируют нелинейное выражение Ех + Оу с четвертым порядком. При расчете течений с сильными ударными волнами, например при расчете маховского отражения ударной волны от плоскости (см. далее), вместо формул (20) целесообразно применение формул с элементами адаптации:

ЯП — Ах^к + Лу ,

Лх — АХ (1 - 5+5-/6), Лу — А0 (1 - 5+5-/6) •

Их применение позволяет избежать появления отрицательных значений давления и плотности и других аварийных ситуаций.

Следует отметить принципиальное отличие трехшаговых и четырехшаговых схем, описанных в предыдущем разделе. Использование факторизации стабилизирующего оператора, существенно облегчающей реализацию схем в многомерных задачах, оставляет неизменным порядок точности четырехшаговой схемы и снижает его до третьего в случае трехшаговой схемы. Это объясняется тем, что стабилизирующий оператор первого шага этой схемы (12)-(14) содержит недиффузионные слагаемые порядка О (т), которые при факторизации породят аппроксимационную погрешность в уравнении (12) порядка О (т)2, что понизит итоговый порядок схемы до О (т)3. Поэтому для сохранения порядка аппроксимации первый шаг этой схемы заменяется на два, имеющих форму начальной стадии итерационного процесса для построения нефакторизованного разностного уравнения. Для наших целей проводить его до конца нет необходимости. В итоге имеем следующее обобщение схемы (12)-(14):

т БЕ

1 + (Бх) 4 Биу '

-1

ЛхЯх - - у;

1 + тБЕ (5 у)

4 В\]У '

-1 Лу Бу - - Уу*

и(1) - ип

т/2

+

+ЯП — (У;п + УП) ип;

(24)

т пр т

1+4 Ш (БХ)-1 Л;Бх - 2У;*

т пр т

1 + (Бу )-1 Лу Бу - -Уу

4 Пи

2

и(2) и(1)

т/2

+

+

и(1) — и(п)

т/2

+

ЯП + Я(11) - (Ужп + Ууп) ип - (Уж(1) + Уу(1)) и(1) /2 — 0;

(25)

*

(l + cr4r4x - 2 v;) (1 + £r4ry - 2v;)--— + r12) = (Vxn + v;) uu; (26)

(1 + ct4R4x - 2v;) (l + CT4R4y - 2v;) uu+_-uu

2-J — +

(du 2 R(2) R(3) \

Rt+2rT+Rh) = (VU+VU) i/u (27)

Здесь Sx, Sy — матрицы, диагонализирующие матрицы

DF/DU, DG/DU : Sx(DF/DU)(Sx)-1 = Lx, Lx — диагональная матрица, остальные величины определены выше.

3. Результаты расчетов

Сначала сформулируем основные итоги расчетов, которые далее частично проиллюстрируем примерами.

1. Схема (24)-(27) с векторными прогонками на первом и втором шагах существенно скорее сходится к стационарным решениям (до двух раз) и допускает расчеты при больших числах Куранта, чем схема со скалярными прогонками на всех шагах (16)-(19).

2. Схема (24)-(27) при нулевой искусственной вязкости b = 0 допускает расчет лишь непрерывных течений. При появлении разрывов приходится увеличивать этот коэффициент, причем значение b = 1/32 обеспечивает расчет почти любых разрывов.

Продемонстрируем работоспособность схемы на задачах идеального газа. Наибольшие проблемы с обеспечением сходимости по времени возникают при расчете разрывных решений уравнений газовой динамики. Поэтому для иллюстрации качества построенного численного метода рассмотрим результаты расчета двумерной задачи распада разрыва. Пусть в начальный момент заданы два параллельных полубесконечных потока с P =1, р = 1, u = 2.4 (1.4)1/2 (при y > 0) и P = 0.25, р = 0.5, u = 4(1.4)1/2 ( при y < 0). Отметим, для этих параметров ударная волна распространяется вниз, далее — контактный разрыв и выше — веер волн разрежения. Рассматривается область 0 < x < 1, -0.3 - 0.3x < y < 0.3 + 0.3x. Используется сетка из 45 х 60 узлов.

На рис. 1 показана динамика сходимости при разных числах Куранта для схемы (16)-(19) слева и для схемы (24)-(27) справа. Схема (16)—(19) использовалась с набором параметров (9), (11). При замене формул (11) на формулы (10) результаты меняются незначительно. Изображается эволюция невязки плотности R = log10 [max^ | pU+VpUc - 1 I /т]. Использование шага по времени, соответствующего следующему в этой последовательности числу Куранта, -5.6 слева и 8 справа — вызывает отсутствие сходимости. Оптимальная скорость сходимости наблюдается при числах Куранта 2... 2.8, причем скорость сходимости схемы (24)-(27) оказалась лишь несколько хуже (требуется в 1.5 раза больше итераций) оптимальной скорости сходимости двухшаговой схемы Рунге — Кутты третьего порядка [9].

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

Рис. 1. Динамика сходимости: слева — для схемы (16)—(19); справа — для схемы (24)—(27).

-4 -2 0 2 4 -4 -2 0 2 4

Рис. 2. Задача о пассивной конвекции вихря, плотность: слева — для схемы (16)-(19); справа — для схемы (24)-(27).

имеющая точное решение. Пусть в начальный момент времени течение описывается формулами

s = P/p7 = 1, T = P/p = 1 - e2e1-r2, e = 5, r2 = x2 + y2,

87П2

u = 1 - e(l-r2)/2ye/ (2n), v = 1+ e(l-r2)/2xe/ (2n), -5 < x < 5, -5 < y < 5.

Точное решение заключается в смещении исходных распределений в соответствии со средними значениями компонент скорости u = 1, v = 1. Численное интегрирование проводилось на сетке 81 х 81 до момента времени t = 100, когда исходный вихрь переместился на 10 размеров расчетной области по каждой переменной. В процессе численного интегрирования расчетная область периодически сдвигалась вслед за перемещением вихря.

На рис. 2 приведены результаты моделирования конвекции вихря схемой (16)—(19) слева и для схемы (24)-(27) справа. Изображены профили плотности в сечении y = 0. Крестики соответствуют точному решению, треугольники — расчетным данным для числа Куранта, равного 1.2, прямоугольники — для K = 2.3, кружки — для K = 3.5. Согласие численного и точного решений хорошее для K = 1.2, значительно ухудшается при увеличении числа Куранта доК = 2.3 иК = 3.5, причем решение для схемы (16)-(19) при К = 3.5 не приводится, поскольку оно совершенно не передает точное решение. Таким образом, схема (16)-(19) существенно уступает в точности. Для объяснения этого факта следует указать различие в коэффициентах этих схем. В стабилизирующих операторах обеих схем слагаемое наиболее высокой (четвертой) степени по локальному числу Куранта (£т4R4) содержит коэффициент £ = ß3/384 для схемы (24)-(27) и значительно больший коэффициент £ = ß3/72 для схемы (16)-(19). Кроме того, последняя схема содержит также быстрорастущие с увеличением числа Куранта слагаемые типа т3^VP2, отсутствующие в схеме (24)-(27). По-видимому, это и объясняет различие в точности при больших числах Куранта. Отметим, что погрешность неявных схем Рунге — Кутты третьего порядка [9] при K = 1 .2 значительно больше, чем погрешность схем четвертого, однако с увеличением шага по времени она растет не столь быстро. В [6, гл. 17] приведены аналогичные данные для явных ENO-схем второго-пятого порядков. Погрешность схем второго порядка значительна и решение не является приемлемым, расчеты по схемам 4-5-го порядков позволяют получить очень хорошее решение.

Моделирование нестационарных разрывных течений представляет трудную задачу для существующих неявных схем. Рассматриваемые схемы на завершающем шаге имеют скалярные консервативные стабилизирующие операторы, что позволяет использовать их для этой цели. Рассмотрим результаты решения задачи о взаимодействии вихря с ударной волной. Используется вычислительная область [0,1] х [0,1]. Стационарный скачок расположен посередине области нормально к оси x. Параметры слева от скачка принимают значения (р, u, v, P) = (1,1.1 -y/Y, 0,1), y — показатель адиабаты. Небольшой вихрь расположен в потоке левее скачка с центром в точке (xc,yc) = (0.25, 0.5). Вихрь задается возмущением скорости (u, v) и температуры (T = P/p) по отношению к среднему потоку:

u = 1 + ea(l-r2|rl) е (y - yc) /re, v = -ea(l-r2|rl) б (x - Xc) /re,

(y - 1) e2a(1-r2Irl) 2

T = P/р =1 - (Y 1) e-е2, P/pY = 1,

4aY

б

Рис. 3. Взаимодействие вихря и ударной волны, изолинии давления: сетка 251 х 100, 1.2 < K < 1.4, £ = 0.05, 0.20, 0.35 (а — в соответственно).

п

(л

Рис. 4. Маховское отражение ударной волны от плоскости, 0.4 < K < 0.8.

где г = у (х - хс)2 + (у - ус)2; е имеет смысл интенсивности вихря; а управляет скоростью убывания возмущения; ге — критический радиус, внутри которого вихрь имеет максимальную интенсивность. В расчетах выбрано е = 0.3, ге = 0.05 и а = 0.204. Используется сетка 251 х 100, которая равномерна по у, но сгущается по х возле скачка. Результаты (изолинии давления) показаны на рис. 3 для разных моментов времени.

Выше рассмотрена задача с нестационарными скачками небольшой интенсивности. Работоспособность схем для расчета течений с сильными разрывами проиллюстрируем на задаче о маховском отражении скачка от плоскости. Используется прямоугольная вычислительная область с равномерной сеткой. Стенка расположена внизу вычислительной области начиная с точки х = 1/6. Перемещающийся вправо скачок с числом Маха = 10 расположен первоначально в точке х = 1/6, у = 0 и имеет угол 60° с осью х. Для границы внизу области точные параметры за скачком приняты на участке от х = 0 до х = 1/6, отражающие граничные условия используются для остальной части границы. На границе сверху значения потока соответствуют точному движению скачка с числом Маха, равным 10. На рис. 4 приведены изолинии давления на момент времени I = 0.2 для разного количества узлов сетки. Ударные волны переданы правильно в расчетах, однако для разрешения сложных структур возле точки ветвления необходимо большое число узлов сетки, и результаты близки к правильным лишь на сетке 541 х 135 (см. [6]).

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

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

[1] РУСАНОВ В.В. Разностные схемы третьего порядка точности для сквозного расчета

разрывных решений // Докл. АН СССР. 1968. Т. 180, № 6. С. 1303-1305.

[2] Burstein S.Z., Mirin A.A. Third order difference methods for hyperbolic equations // J. Comput. Phys. 1970. Vol. 5, N 3. P. 547-571.

[3] Kutler P., Lomax H., Warming R. Computation of Space Shuttle Flow Fields Using Noncentered Finite-difference Schemes. AIAA Paper. 1972. N 72-193. 25 p.

[4] Shu C.-W. Total-variation-diminishing time discretizations // SIAM J. Sci. Stat. Comput. 1988. Vol. 9, N 6. P. 1073-1084.

[5] Osher S., Shu C.-W. Efficient implementation of essentially nonoscillatory shock-capturring schemes // J. Comput. Phys. 1988. Vol. 77, N 2. P. 439-471.

[6] Пинчуков В.И., Чи-ВАнг Шу. Численные методы высоких порядков для задач аэрогидродинамики. Новосибирск: Изд-во Сиб. отд-ния РАН, 2000. 232 с.

[7] Толстых А.И. Компактные разностные схемы и их приложения к проблемам аэрогидродинамики. М.: Наука, 1990. 230 с.

[8] Пинчуков В.И. Абсолютно устойчивые схемы Рунге — Кутты третьего порядка // Журн. вычисл. математики и мат. физики. 1999 Т. 38, № 11. С. 1855-1868.

[9] Пинчуков В.И. Сравнение неявных схем Рунге — Кутты третьего порядка // Вы-числ. технологии. 2002. Т. 7, № 5. С. 44-57.

[10] Пинчуков В.И. О неявных абсолютно устойчивых схемах Рунге — Кутты четвертого порядка // Вычисл. технологии. 2002. Т. 7, № 1. C. 96-105.

[11] Пинчуков В.И. Трех- и четырехшаговые неявные абсолютно устойчивые схемы Рунге — Кутты четвертого порядка // Журн. вычисл. математики и мат. физики. (В печати).

[12] Пинчуков В.И. О численном исследовании трансзвуковых турбулентных течений возле крыла неявными схемами высоких порядков // Вычисл. технологии. 2001. T. 6, № 2. С 110-121.

[13] Jameson A., Schmidt W., Turcel E. Numerical Solution of the Euler Equations by Finite-volume Method Using Runge — Kutta Time Stepping Schemes // AIAA Paper. 1981. P. 81-1259.

[14] Пинчуков В.И. Нелинейные сеточные фильтры и их использование в схемах высоких порядков для задач аэродинамики // Мат. моделирование. 1999. Т. 10, № 11. С. 111-115.

Поступила в редакцию 26 июня 2003 г.