CC BY
36
Вычислительные технологии
Том 7, № 5, 2002
СРАВНЕНИЕ НЕЯВНЫХ СХЕМ РУНГЕ-КУТТЫ
ТРЕТЬЕГО ПОРЯДКА *
В. И. Пинчуков Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: pinchvi@ict.nsc.ru
Implicit third order Runge — Kutta schemes are constructed for multidimensional transfer equation with diffusion and for compressible gas equations. Linear stability is proved for any time step. Two and three steps schemes are considered. Various sets of coefficients are founded to provide absolute stability of every type of scheme. The results of numerical calculations are presented.
В последние годы в различных областях вычислительной математики, в частности в вычислительной аэродинамике, быстро развиваются схемы высоких порядков. Весьма плодотворным оказался ENO-подход (см., например, обзор [1]), основанный на конструировании схем высоких порядков посредством выбора формул с наименьшими осцилляционны-ми свойствами из класса формул одного порядка для определения каких-либо элементов схем, как правило, потоков через границы элементарной разностной ячейки. В качестве примера упомянем адаптивные схемы высоких порядков типа Рунге — Кутты и Адамса [1, гл. 16, 2]. Следует отметить, что ввиду переменчивости шаблона этих схем линейный анализ устойчивости их сопряжен с техническими сложностями и обычно не проводится. Однако для некоторых задач линейная устойчивость схемы весьма важна, в ее отсутствие возможно, по нашему мнению, появление численных решений, не сходящихся к стационарным при интегрировании на больших временных интервалах. Примеры построения явных линейно устойчивых схем третьего порядка содержатся в классических работах [3, 4].
Во многих задачах предпочтительными являются неявные схемы. В настоящей работе предлагаются неявные схемы, устойчивые в приближении замороженных коэффициентов для любых временных шагов. Проводится сравнение эффективности как разных типов схем, так и схем одного типа с разными коэффициентами.
1. Схемы для диффузионного уравнения переноса
Анализ схем произведем на примере многомерного диффузионного уравнения переноса
* Работа поддержана Американским фондом гражданских исследований и развития для независимых государств, образованных на территории бывшего Советского Союза (СКОР), грант Б,М1-2324-МО-02. © В. И. Пинчуков, 2002.
m
m
(1)
Введем обозначения и = ±и(...,хг ± Ахг, хг,...), А±г = /Ахг и =
т т
агАХг (1-/6) — оператор конвективного переноса 4-го порядка; Ш = ^ А+ /2 — г=1 г=1
оператор диффузии 2-го порядка с множителем 1/2, введенным для удобства выкладок;
т
Л2к — семейство вспомогательных операторов Я2к = (-1)к ^ (А- А+ )2к.
г=1
Из анализа устойчивости следует возможный вид двухшаговой схемы (3/2 + т^/2 + £т4Я4 - Ат2Я2 - тЖк)(/* - /п)/т + (^1 - 2Ш)/п = -2т3рЯ4/п, (2) (1 + ^т4Я4 - ПТ2Я2 - тШ)(/п+1 - Л/т + (тпЯ2 + Ш)(/* - /п)3/2 +
+ (^1 - 2Ш)(/п/4 + /*3/4) = -2т3^Д4/га. (3)
Второе слагаемое в уравнении (3) введено для компенсации членов низкого порядка в стабилизирующем операторе. Если отказаться от требования второго порядка по времени для диффузии, то это слагаемое можно модифицировать. Итак, рассмотрим также схему с другим представлением вязких слагаемых
(3/2 + т^/2 + £т4Я4 - Ат2Я2 - тШк)(/* - /п)/т + (^ - 2Ш)/п = -2т^4/п, (4) (1 + ^т4Я4 - пт2Я2 - тШ)(/п+1 - /п)/т + (тп^2 - Ш)(/* - /п)3/2 +
+ (^1 - 2Ш)(/п/4 + /*3/4) = -2т3^Д4/п. (5)
Изменение знака перед Ш во втором слагаемом уравнения (3) приводит к понижению порядка по времени до первого для вязких членов. Если эти вязкие члены имеют смысл искусственной диффузии, то это понижение не является существенным для точности схемы, в то же время, как показали тестовые расчеты течений газовой динамики (см. ниже), позволяет использовать больший шаг по времени.
Параметры схемы определялись в процессе анализа устойчивости ее на основе метода Фурье. Требование ограниченности множителя перехода схемы эквивалентно неравенству на степенную функцию от четырех, как мы увидим, независимых переменных. В процессе рассмотрения различных асимптотических случаев определялась часть коэффициентов, при этом нужное неравенство упрощалось и становилось возможным обеспечить его выполнение за счет выбора оставшихся коэффициентов. Таким образом были найдены наборы коэффициентов (6) и (9), а также параметрические семейства (10), (11) и (12), (13). Другим способом определения устойчивой схемы было задание коэффициентов или соотношений между ними, за исключением одного, например р, и отыскание этого последнего методом перебора, так чтобы обеспечить выполнение условия | р |< 1 + е. Поскольку для нулевого волнового вектора р =1 и возможны ошибки округления, следует выбирать е > 0. Обычно использовалось значение е = 0.0001. Так найдены наборы (7), (8).
В результате можно сформулировать теорему:
Теорема 1. Схема (2), (3) устойчива при любых значениях т, если ее коэффициенты
принимают значения
П = 6/48, А = п3/2, к = 3/2, £ = р3/2, ^ = р = 63/72, (6)
А = 0, п = 6/18, к = 3/2, £ = р3/2, ^ = р = 63/86, (7)
А = -6/18, п = 0, к = 3/2, £ = р3/2, ^ = р = 63/513, (8)
А = 0, п = 0, к = 1, £ = р3/2, ^ = р = 63/135. (9)
Схемы с указанными значениями параметров не обладают свойством полной аппроксимации, т.е. стационарное решение, полученное с их помощью, зависит от временного шага, что неудобно при их использовании как итерационных методов. Представляет интерес отыскание наборов параметров, в которых ^ = 0 и которые обеспечивают абсолютную устойчивость схемы. Такие наборы были найдены, причем найдены параметрические семейства схем. Особенно обширное семейство коэффициентов существует в отсутствие вязких слагаемых. Далее будет доказана
Теорема 2. Схема (2), (3) при V = 0 абсолютно устойчива, если ее параметры принимают значения
П = -2А - 6/12 - г, к = 3/2, £ = р3/2, ^ = 0, (10)
р = 6[(6/12 + г)(2А + 6/12 + г) + А2/0.75]2/[4ф + 6/12)], г> 0. (11)
Итак, есть два свободных параметра — А и г. Следует отметить, что первая из формул (10) имеет следствием соотношение п + 2А < 0, таким образом, хотя бы один из коэффициентов п или А отрицателен. Это нестандартная черта данной схемы, так как на неявный слой принято выносить положительные операторы. Тем не менее экспериментальные расчеты продемонстрировали эффективность схемы.
При обобщении ее на случай уравнения переноса с диффузией с целью упрощения анализа устойчивости семейство коэффициентов было выбрано более узким. В частности, параметр А принят равным А = -(г+6/12)3/4. Это значение соответствует экстремуму функции (11), оно минимизирует коэффициент р. Коэффициент А, как можно видеть, отрицателен, соответствующий коэффициент п, вычисленный по формуле (10), п = (г + 6/12)/2, положителен. Следует отметить, что дополнительные, обусловленные вязкими членами требования устойчивости на коэффициенты схемы выполнены с большим запасом, как будет видно из нижеследующего анализа. Поэтому, возможно, существует и более широкое семейство устойчивых схем для уравнения переноса с диффузией.
Итак, верна
Теорема 3. Схема (2), (3) устойчива при любых значениях т, если ее параметры принимают значения
п =(6/12 + г)/2, А = -(6/12 + г)3/4, к = 3/2, £ = р3/2, ^ = 0, (12)
р = 6(6/12 + г)2/16 тах[(6/12 + г)/(4г), 1] г> 0. (13)
Первая ветвь функции (13) соответствует функции (11).
Мы ограничимся здесь лишь доказательством теорем 2, 3. Проверка их вычислением на компьютере множителя перехода (15) требует много арифметических операций, так как к перебору четырех независимых переменных Р^ Р2, Р4, Ш в данном случае добавляется необходимость перебора свободных параметров г, А. Теорема 1 с формулами (6) доказана
в [1, 5].
Исключая промежуточную функцию f *, схему (2), (3) можно свести к уравнению
(3/2 + т#1/2 + £т4Я4 + Ат2Я2 - тШк)(1 + рт4Я4 - пт^2 - тШ)ап+1 - Л/т+ + (3/2 + т2АЯ2 + т2п^23/2 - т#1 /4 + £т4Я4 - тШк)(^1 - 2Ш + 2т3^4)Г = 0. (14)
Для исследования устойчивости применим метод Фурье. Представим решение в виде г = рп ехр (Е ОНц^ , 0 < а < 2п, ц =1, ..., /г, ] =
Введя вспомогательные величины Pk, k = 1, 2, 4, связанные с операторами (rfcRk) и даваемые формулами
m m m
Pi = ^2 Oi sin ai (1+qi/6), P2 = ^ ofqi, P4 = ^ afqj2, 07 = a¡ т/Axi, q¡ = 2- 2cos ai, i=i i=i i=i
m
и V = dlrql/(2Ax¿2) > 0 — спектральный образ оператора (-tW), из уравнения (14) i=i
легко получить множитель перехода для схемы (2), (3):
р =1 - [(3/2 + ЛP2 + nP¿3/2 + + Vk)(2V + 2^4) + P2/4+
+jPi(3/2 + ЛP2 + nP¿3/2 + £P4 + VK - V/2 - ^4/2)]/
/(3/2 + jPi/2 + £P4 + ЛP2 + tVk)/ü, h = 1 + nP2 + <P4 + V. (15)
Отметим соотношения [5] между величинами P-, P2, P4, которые потребуются при исследовании устойчивости:
Pi2 < bP2, b = m(1 + 7/243), P22 < bP4. (16)
Подставим в (15) вместо параметра £ его значение £ = <3/2. Легко получить, что при V = 0 и ^ = 0 условие устойчивости | р |< 1 подстановкой р из (15) сводится к неравенству
P?[h/2 - (3/2 + ЛP2 + P4<3/2 + nP23/2)]2+
+ [h(3/2 + ЛP2 + P4<3/2) - P2/4]2 < h2[(3/2 + ЛP2 + P4<3/2)2 + P2/4].
Если перенести все слагаемые в левую часть, то получим квадратный трехчлен относительно аргумента P-j2. Поскольку коэффициент при старшей степени этого аргумента положителен, квадратный трехчлен максимален при минимальном и максимальном значениях этого аргумента, т. е. при P-j2 = 0 и ввиду (16) Pj2 = bP2. Подставив первое значение, имеем в левой и правой частях последнего неравенства одинаковое выражение, т. е. оно верно. Подставив Pj2 = bP2, имеем
[h(3/2 + ЛP2 + P4<3/2) - bP2/4]2 - h2(3/2 + ЛP2 + P4<3/2)2+ +bP2[h/2 - (3/2 + ЛP2 + P4<3/2 + nP23/2)]2 - h2bP2/4 < 0.
Раскрывая разности квадратов по формуле x2 - y2 = (x - y)(x + y) и подставляя выражение для h, h = 1 + nP2 + <P4, умножая на -1, имеем
[(1 + nP2 + <P4)(3 + 2ЛP2 + 3P4<) - bP2/4]bP2/4-
-bP2/4 (1 + 2ЛP2 + P4< + nP2)(3 + 2ЛP2 + 3P4< + 3nP2) > 0. (17)
Отбрасывая неотрицательный множитель bP2/4 и сокращая одинаковые слагаемые, получим
(1 + 2ЛP2 + P4< + nP2)3nP2 + 2ЛP2(3 + 2ЛP2 + 3P4<) + bP2/4 > 0.
Сокращая Р2 > 0, подставляя выражение (10) для п, имеем
Р4<3(Ь/12 + г) - (Ь/12 + г)Р23(2Л + Ь/12 + г) - 4Л2Р2 + 3г > 0.
Поскольку коэффициент при Р4 положителен, можно заменить Р4 меньшим в силу второго из неравенств (16) значением Р|/Ь. В результате в левой части последнего неравенства получаем квадратный трехчлен относительно переменной Р2. Легко проверить, что выражение (11) для < делает равным нулю дискриминант этого трехчлена, т.е. это неравенство выполнено. Таким образом, теорема 2 доказана.
Доказательство теоремы 3 аналогично доказательству теоремы 2. Легко получить, что условие устойчивости | р |< 1 подстановкой р из (15) и заменами к = 3/2, £ = <3/2, ^ = 0 сводится к неравенству
[к(3/2 + ЛР2 + Р4<3/2 + V3/2) - V(3 + 2ЛР2 + 3Р4< + 3У + 3^) - р2/4]2+ +^>/2 - (3/2 + ЛР2 + Р4<3/2 + п^23/2 - V)]2 <
< к2 [(3/2 + ЛР2 + Р4<3/2 + V3/2)2 + Р2/4]. (18)
Если перенести все слагаемые в левую часть, то получим квадратный трехчлен относительно аргумента Р-р. Поскольку коэффициент при старшей степени этого аргумента положителен, квадратный трехчлен максимален при минимальном и максимальном значениях этого аргумента, т.е. при Р1 = 0 и ввиду (16) Р\ = ЬР2. Подставив в (18) первое значение, имеем
[к(3/2 + ЛР2 + Р4<3/2 + V3/2) - V(3 + 2ЛР2 + 3Р4< + 3V + 3пР2)]2--к2(3/2 + ЛР2 + Р4<3/2 + V3/2)2 < 0.
Раскрывая разность квадратов и подставляя в полученное неравенство выражение для к, к = 1 + пР2 + <Р4 + V, после преобразований запишем
- V (3 + 2ЛР2 + 3Р4< + 3V + 3пР2) х
х [(1 + пР2 + <Р4)(3 + 2ЛР2 + 3Р4<) + (1 + <Р4^] < 0,
куда подставим выражения (12) для параметров Л, п, после чего получим
(1 + Р4< + V)[(1 + Р2^/2 + <Р4)(1 - Р2^/2 + Р4<) + (1 + <Р4) V] < 0,
где d = Ь/12 + г > 0. Из фигурирующих здесь выражений в круглых скобках лишь третье может быть отрицательным. Покажем, что имеет место соотношение
1 - P2d/2 + Р4< > 0. (19)
Заменим в нем Р4 на меньшее в силу (16) выражение Р|/Ь. Чтобы получившийся квадратный трехчлен относительно Р^ был неотрицательным, достаточно неотрицательности его дискриминанта, что дается условием < > М2/16. Это выполнено для величины даваемой формулой (13). Следовательно, (19) верно. Поскольку все выражения в квадратных скобках предпоследнего неравенства неотрицательны, оно выполнено.
Таким образом, исходное неравенство (18) верно при р2 — 0. Доказав его при р2 — 6Р2, мы обоснуем абсолютную устойчивость схемы. Итак, рассмотрим получающееся при указанной замене неравенство
[Л(3/2 + ЛР2 + Р^3/2 + V3/2) - V(3 + 2АР2 + + ЗУ + 3пР) - 6Р2/4]2+ +6Р2[й/2 - (3/2 + ЛР2 + Р4^3/2 + ^3/2 + V)]2--Л2[(3/2 + ЛР2 + Р4^3/2 + V3/2)2 + 6Р2/4] < 0. (20)
Представляя в (20) разности квадратов по формуле х2 — у2 — (х — у)(х + у), подставляя выражение для Л, умножая на —1 и группируя слагаемые по степеням параметра V, можно получить
9V3(1 + ^4) + ^2[(1 + пР + ^Р4)(1 + 2/3 ЛР2 + Р^) - ЬР2/12+
+ (1 + 2/3 ЛР2 + Р^ + Пр2)(1 + ^4)] + +V [(3 + 2ЛР2 + 3Р^ + 3пРг)(1 + Пр + ¥>Р0(3 + 2ЛР2 + 3Р^)--ЬР2/4(2 + 6ЛР2 + 2Р^ + 5^)] + + (ЬР2/4)[(1 + пР + ^Р0(3 + 2ЛР2 + 3Р^) - &Р2/4--(3 + 2ЛР2 + 3Р^ + 3пРг)(1 + 2ЛР2 + Р^ + ПР2)] > 0.
Последнее слагаемое равно левой части уже доказанного неравенства (17), т.е. оно неотрицательно. Очевидно, что первое слагаемое так же неотрицательно. Таким образом, достаточно доказать неотрицательность выражений в квадратных скобках при второй и первой степенях V. Используя формулы (12) для параметров Л, п, можно свести исследование неотрицательности выражения при второй степени к доказательству неравенства
2 + 4Р4^ + 2Р42^2 - ^2Р|/4 - 6Р2/12 > 0.
Заменим Р4 на минимально возможное согласно (16) значение Р|/6, уменьшив левую часть этого неравенства. Получим
(р/6 - ^2/16)4Р22 + (2 + 2Р24^2/62 - Р^) > 0.
Неотрицательность выражения в первой скобке сразу следует из определения (13) для
Выражение во второй скобке имеет единственный экстремум при Р2 — 6/(^296)(1/3), который соответствует минимальному значению этого выражения. Нетрудно проверить, это минимальное значение положительно, если выполнено соотношение ^ > 63/(32)2/3(1/2). Можно показать, что если в первой ветви формулы (13) выбрать значение свободного параметра г, минимизирующее эту функцию, т.е. положить г — 6/24, то нужное соотношение будет выполнено даже в этом неблагоприятном случае. Таким образом, коэффициент при второй степени V неотрицателен.
Аналогично неотрицательность выражения при первой степени эквивалентна неравенству
(1 + Р^)(1 + ^2/2 + ^£0(1 - ^2/2 + Р4¥>) - (1 - ^2 + Р4^)6Р2/18 > 0.
Увеличим вычитаемое, заменив ЬР2/18 на ЬР2/18 + г/1.5 = ^/1.5 и (1 — ^Р2 + Р4р) на (1 — ^Р2/2 + Р4р). Неотрицательный в силу соотношения (19) множитель (1 — ^Р2/2 + Р4р) можно сократить, в итоге получаем неравенство
(1 + Р4^)(1 + ^/2 + ^Р4) — ^2/3 > 0.
Заменим множитель (1 + ^Р2/2 + Р4р) на Р2^, лишь уменьшив в силу соотношения (19) левую часть неравенства. Неотрицательный множитель Р2^ можно сократить. В итоге получаем очевидное соотношение
1 + Р4Р — 2/3 > 0.
Таким образом, теорема доказана.
Первый шаг схемы (2), (3) реализует вычисление промежуточного значения а* с третьим порядком точности. Для достижения этой цели можно использовать два шага, образующие схему типа предиктор — корректор. В итоге можно получить следующий возможный вид трехшаговой схемы Рунге — Кутты:
а * — Г)3/(2т) + = 2ЖГ, (21)
(3/2 + Ат — ктЖ )(а ** — Л/т + Я1(Г + а *)/2 = (22)
(1 + ^т4Я4 + Ст 3Я2 W — тЖ )(Г+1 — Л/т + ^(^/4 + а **3/4) = 2ЖЛ (23)
Здесь Л1, Л2, Л4, Ж — разностные операторы, определенные выше.
Первая неявная трехшаговая схема третьего порядка построена в [5,6]. Ее коэффициенты даются формулами (24). В [7] обоснована абсолютная устойчивость схемы с меньшими коэффициентами стабилизирующего оператора, приведенными в формулах (25). Итак, верна
Теорема 4. Схема (21)-(23) устойчива при любых значениях т, если ее параметры принимают значения
А = С = т250/243/4, р = т3250/243/36, к = 3/2, (24)
А = С = т250/243/6, р = т3250/243/36, к = 1, (25)
где т — число пространственных измерений.
Наряду с аналитическим исследованием устойчивости производилось прямое вычисление коэффициента перехода, возникающего при использовании метода Фурье. Для приведенных выше параметров схемы максимальное значение коэффициента перехода равно единице с точностью по меньшей мере семи значащих цифр. Следует отметить, значение параметра р однозначно определяется в анализе устойчивости и уменьшено быть не может. В то же время значения (24), (25) параметров А и С, по-видимому, зависят от найденного конкретного доказательства устойчивости и могут быть завышенными. В пробных вычислениях коэффициента перехода получено, что при одновременном уменьшении параметров А и С от значения, данного в формулах (24), до значения А = С = 0.004 Ь (Ь = т 250/243 = 2 х 250/243) этот коэффициент практически не меняется, потом в диапазоне от 0.004 Ь до 0 он возрастает до значения 1.00017. Таким образом, при 0 < А, С < 0.004 Ь схема (21)-(23) имеет слабо растущие гармоники, при А, С > 0.004 Ь схема абсолютно устойчива. Однако в численных расчетах (см. далее) оптимальная сходимость наблюдалась при А = С = 0.
2. Двухшаговая схема для уравнений идеального газа
Запишем уравнения сжимаемого газа в консервативной форме:
и + ^ + Су = 0. (26)
Обобщение схемы (2), (3) произведено с упрощенным представлением искусственной диффузии в стабилизирующем операторе (весовой коэффициент принят равным единице):
3 т (sx)-- т2ar2x + £т4r4x - twx
2 2 DU4 7
X
X
3 T DG(Sy)- 1 ЛуSy - T2AR2y + £т4R4y - tWy
2
(Uk - u* )3T + Rn
2 2 DU
= (W - 2^т3Д4* + W - )U*, (27)
[1 - ПТ2R2x + <T4R4x - TWX] [1 - ПТ2R2y + <T4R4y - tWy] (Uk+1 - un)/T+ +(TnR2x + Wx + TnR2y + Wy)(U* - Un)3/2 + (Rn + R*3)/4 =
= (Wx - 2^т3R4x + Wy - 2^T3R4y)(U™fc + 3Uik)/4. (28)
Здесь A, n, < — параметры, заданные в формулах (6) - (13); R2x = Д+(аж)2Д-, ax = c+|u|, c — скорость звука, величины с индексом y вычисляются аналогично;
Ri = + ЛуGifc, Лх = (1 - №/6)ДХ, Лу = (1 - ; (29)
Sy — матрицы, диагонализирующие матрицы DF/DU, DG/DU: Sx(DF/DU)(Sx)-1 = Lx, Lx — диагональная матрица; Wx, Wy — нелинейные адаптивные операторы искусственной диффузии, переключающиеся в зависимости от гладкости решения (диффузия третьего порядка на "гладких" решениях и первого возле ударных волн):
WxUik = Д-И+^к - i+e^WU*, e(4x) = ax7/8, e(2x) = ax(1 - T)/2, (30)
0 < 7ifc < 1 — коэффициенты локальной гладкости (приводятся ниже), аналогично по адаптивным формулам строятся операторы R4x, R4y:
R4x = Д- [Д+(«Х*)47ifcД- - (aX+i/2*)4(1 - Y+i/2*)4/Дх2]Д+. (31)
Формула (29) аппроксимирует нелинейное выражение + Gy с четвертым порядком. Тем не менее с учетом искусственной диффузии схема будет иметь лишь третий порядок по пространственным переменным. Вес этих диффузионных слагаемых в стабилизирующем операторе повышен до единицы, что не влияет на порядок аппроксимации схемы, в то же время, согласно тестовым расчетам, существенно расширяет область чисел Куранта, при которых схема работоспособна.
Возвращаясь к структуре искусственной диффузии, можно отметить, что она строится на основе классического подхода комбинирования диффузии высокого порядка, предназначенной для лучшей сходимости к стационарным решениям, и диффузии низкого порядка, предназначенной для подавления осцилляций на разрывах. Хотя этот подход разрабатывается давно, построение конкретных процедур, обеспечивающих хорошее качество как гладких, так и разрывных решений, представляет проблему. Здесь при их построении
используются следующие принципы: а) консервативность; б) сохранение монотонности решений для уравнений, включающих лишь производную по времени и диффузионное слагаемое; в) обеспечение диссипации энергии согласно анализу [8]. Для реализации этого подхода величины 0 < 7^ < 1 вычисляются по давлению Р^ на основе формул
Тгк = max[0, шт(ад+, ад-, 1) — с]/(1 — с), (32)
= (*+*- + вР)/[(г+)2 + 10-9], д- = (г+г- + е2)/[(г-)2 + 10-9],
--рП 1Л~ > 1 -- V И '
= 6+Pik, г- = 5ГР^
и являются индикаторами его "гладкости". Параметры а и с позволяют регулировать свойства диффузии, параметр ер предназначен для предотвращения понижения порядка диффузии на гладких экстремумах или участках почти постоянного решения. По результатам пробных расчетов были выбраны следующие значения эмпирических величин: с = 0.4, а = 1.3, ер = 2(Ртах — Рт1П)//, I — количество узлов сетки по переменной х. Отметим, что для рассмотренных формул на гладких решениях имеем 7 = 1 и реализуется режим диссипации высокого порядка. В то же время вблизи разрывов коэффициенты 7 уменьшаются и включается более грубая диссипация второго порядка. Численные расчеты показали, что это обеспечивает как приемлемый уровень монотонности схем для уравнений сжимаемого газа, так и хорошие свойства сходимости к стационарному решению.
3. Трехшаговая схема для уравнений идеального газа
Обобщение схемы (21) - (23) осуществляется с упрощенным представлением искусственной диффузии в стабилизирующем операторе:
3
(ик — ик)- + яп = (Ж* + Жу )ик, (33)
3 3
2 + т3 лдж — тЖ*
3 3
2 + т3л^у — тЖу
2
те — ик) - + (ЯП + Я*)/2 =(Ж* + Жу)ип, (34)
[1 + <Т4Я4* + £т3дж — тЖ^ [1 + <т4Я% + £т3фу — тЖ^ (ип+ 1 — ик)/т+
+(Я? + Я**3)/4=(ЖЖ + Жу )Цк. (35)
Здесь < = т3250/243/36 = (250/243)2/9; параметры Л, £ приведены в формулах (24), (25); Я4ж определяется формулой (32); Я 1 приведен в формуле (29); Жу — нелинейные адаптивные операторы искусственной диффузии, определенные формулой (31).
Стабилизирующие операторы ф*, Qy соответствуют слагаемому Я2Ж в схеме (21)-(23). Стандартным рецептом замены такого произведения при факторизации стабилизирующего оператора является замена стабилизирующего оператора на мажорирующий его более простой, например, по формуле Я2Ж = (Я2ж + Я2у)(Жг + Жу) ^ (Я| + Ж2)/2 ^ ЬЯ4/2 + Ж2 + Ж2. По нашему мнению, эта замена может оказаться слишком грубой в случае очень больших или очень малых значений параметра, отвечающего соотношению конвективных и вязких слагаемых решаемого уравнения. Ориентируясь на использование данной схемы для расчета течений вязкого газа, где упомянутый выше параметр — число Рейнольдса — принимает, как правило, большие значения, рассмотрим здесь другой способ. В двумерном случае имеем
Я2 V = (Я2* + Я2у )(Ж* + Шу) = (^5+5- + Г2у 5+5-)^ж5+5- + Wy 5+5-),
где г2х = а^/Дж2, г2у = а2/Ду2, их = и^/Дж2, = и2/Ду2 — коэффициенты операторов. Если раскрыть это произведение и заменить операторы смешанных производных по формуле ^ + (^+^-)2]/2, в итоге получим
^ (^+^Х")2(г2х^х + Г2х^у/2 + Г2уих/2) + (¿+^-)2(г2у+ Г2у^х/2 + Г2х^у/2).
В рамках этого подхода для вычисления коэффициентов стабилизирующих операторов используются соотношения
^х = Д-[(ах)2 сх + (ах)У/2 + (ау )У/2]т/2к Д-Д+Д-,
х
е(2х) +4е(4х), сУ = 0.5ау.
С
Здесь параметры сх и су соответствуют редукции адаптивных пятиточечных операторов искусственной диффузии Ш и V к трехточечному виду при их представлении в параметр сУ вычисляется приближенно, что позволяет избежать трудоемкого вычисления величин е(2у) и е(4у) в программном модуле для обращения стабилизирующего оператора с индексом "ж", оператор с индексом "у" вычисляется аналогично.
4. Примеры расчетов
Продемонстрируем работоспособность схемы вначале на задачах идеального газа. Наибольшие проблемы с обеспечением сходимости по времени возникают при расчете разрывных решений уравнений газовой динамики. Поэтому для иллюстрации качества построенного численного метода рассмотрим результаты расчета двумерной задачи распада разрыва. Пусть в начальный момент заданы два параллельных полубесконечных потока с параметрами Р =1, р =1, и = 2.4(1.4)1/2 (при у > 0) и Р = 0.25, р = 0.5, и = 4(1.4)1/2 (при у < 0). Отметим, для этих параметров вниз распространяется ударная волна, далее контактный разрыв, и выше — веер волн разрежения. Рассматривается область 0 < ж < 1, -0.3 - 0.3ж < у < 0.3 + 0.3ж. Используется сетка из 45 х 60 узлов.
Эта задача решалась с помощью схемы (27), (28) с разными наборами параметров, заданными в формулах (6)-(9), и с набором, полученным с помощью формул (12), (13) при подстановке свободного параметра г = 6/24,
п = 6/16, А = -36/32, £ = р3/2, ^ = 0, р = 63/1024, (36)
где 6 = 2 х 250/243. Использовалась также модификация схемы (27), (28), аналогичная схеме (4), (5). Как оказалось, она работоспособна при несколько больших числах Куранта. При небольших числах Куранта (0.7, 1) результаты работы схемы с разными наборами коэффициентов оказались практически идентичными. С ростом чисел Куранта различие увеличивается, и при К = 4 для достижения невязки порядка 10-12 схемой с разными наборами параметров требуются разные количества итераций, различающиеся до 20%. Оптимальным оказался набор (8), близкие результаты получены также при параметрах (9) и (36). Поскольку формулы (36) обусловливают выполнение свойства полной аппроксимации, т. е. независимость решения, полученного методом установления, от шага по времени, и обеспечивают скорость сходимости по времени, очень близкую к оптимальной, они приняты как предпочтительные и дальнейшие результаты получены с их использованием.
На рис. 1 слева показана динамика сходимости схемы (27), (28) при разных числах Куранта, справа — для ее модификации, аналогичной (4), (5). Изображается эволюция
невязки плотности Я = ^шах^ | рП+1 /Рп — 1 I /т• Использование шага по времени, соответствующего следующему в этой последовательности числу Куранта (4 для данных слева и 5.6 справа), вызывает отсутствие сходимости. Оптимальная скорость сходимости наблюдается при числах Куранта 2 •••2.8, она близка к оптимальной скорости сходимости простейшей одношаговой максимально неявной схемы с такой же пространственной дискретизацией, т. е. она весьма высока.
-2
-4
-10
-12
; 1 «К=0.7 ^
МОВВйФл, +0 0 4 +°А4 по т4 а и ф ■ +0 П 4 □ К=1 4 К= 1.4 -°К=2
- П 0 Ч++ 0 V 0 Д ++, 0 0 д "„ \ ос сч
о 4 0 ' : : : **ч ч Х>:
4 : 0 4 П Ф
: 0 % 0 А оп
г 0 По 4 0
0 ; °ооос □п : 4 П
200
Рис. 1. Динамика сходимости: слева — данные для схемы (27), (28), справа — для ее модификации по типу (4), (5).
Рис. 2. Динамика сходимости трехшаговой схемы (33)-(35): слева — число Куранта 1.4, схемные параметры (37)-(40), справа — результаты для параметров (40) при разных числах Куранта.
На рис. 2 изображается динамика сходимости схемы (33)-(35). Слева приведены данные, полученные при К = 1.4 (К — число Куранта) для ^ = 23 х 250/243/36 и следующих значений параметров £, Л:
Л = £ = 6/4, (37)
Л = £ = 6/6, (38)
Л = £ = 6 0.004, (39)
Л = £ = 0, (40)
6 =2 х 250/243. Как видим, уменьшение параметров £, Л вызывает ускорение сходимости. Оптимально нулевое значение этих параметров. Как указывалось ранее, в приближении замороженных коэффициентов при нулевом значении параметров £, Л существуют слабо растущие гармоники. Поэтому естественно полагать, что оптимален выбор формулы (39), которая, как и формулы (37) и (38), обеспечивает абсолютную устойчивость схемы, но порождает наименьший коэффициент в стабилизирующем операторе. Однако результаты оказались другими.
Для оптимального значения £ = Л = 0 на рис. 2 справа приводится динамика сходимости при разных числах Куранта. Оптимальным является число Куранта 1.4, которое меньше оптимального числа Куранта схемы (27), (28) и скорость сходимости при котором в два раза меньше оптимальной скорости сходимости этой схемы.
Преимущество схем высоких порядков наглядно проявляется при расчетах гладких течений на больших интервалах по времени. С помощью схемы Рунге — Кутты третьего порядка была просчитана задача о пассивной конвекции вихря в однородном потоке, имеющая точное решение. Пусть в начальный момент времени течение описывается формулами
s = P/pY = 1, T = P/p =1 - e2 1 e1-r2, б = 5, r2 = x2 + y2,
87П2
u =1 - e(1-r2)/V/(2n), v =1 + е(1-г2)/2хб/(2п), -5 < x < 5, -5 < y < 5.
Точное решение заключается в смещении исходных распределений в соответствии со средними значениями компонент скорости u = 1, v = 1. Численное интегрирование проводилось на сетке 81 х 81 до момента времени t = 100, когда исходный вихрь переместился на 10 размеров расчетной области по каждой переменной. В процессе численного интегрирования расчетная область периодически сдвигалась вслед за перемещением вихря.
На рис. 3 приведены результаты моделирования конвекции вихря схемой (27), (28), на рис. 4 — схемой (33)-(35). Погрешность приблизительно одинакова во всех случаях, что говорит о малом вкладе в суммарную погрешность ее временной составляющей. В [1, гл. 17] приведены аналогичные данные для схем второго — пятого порядков. Погрешность схем второго порядка значительно больше, и решение не является приемлемым, расчеты по схемам 4-5-го порядков позволяют получить очень хорошее решение.
Итак, предложенные схемы являются удобным и достаточно эффективным инструментом решения задач газовой динамики. Затраты на реализацию стабилизирующих операторов относительно невелики и не превышают 25% обших затрат для схемы (27), (28) и 7% для схемы (33)-(35). В то же время схемы позволяют использовать несколько больший, чем явные схемы, шаг по времени, они сохраняют консервативный характер, что является проблемой для многих неявных схем, имеют неплохую скорость сходимости по времени за счет эффективной искусственной вязкости и повышенный порядок точности для гладких решений. Двухшаговая схема наиболее эффективна по всем критериям.
Рис. 3. Задача о пассивной конвекции вихря, плотность: треугольники — число Куранта 0.7, прямоугольники — число Куранта 1.1, крестики — точное решение.
Рис. 4. Задача о пассивной конвекции вихря, плотность: треугольники — число Куранта 1.1, прямоугольники — число Куранта 2.2, крестики — точное решение.
Список литературы
[1] Пинчуков В. И., Шу Ч.-В. Численные методы высоких порядков для задач аэрогидродинамики. Новосибирск: Изд-во СО РАН, 2000. 232 с.
[2] Shu C.-W. Total-variation-diminishing time discretizations // SIAM J. Sci. Stat. Comput. 1988. Vol. 9, No. 6. P. 1073-1084.
[3] Русанов В. В. Разностные схемы третьего порядка точности для сквозного расчета разрывных решений // Докл. АН СССР. 1968. Т. 180, №6. С. 1303-1305.
[4] Burstein S. Z., Mirin A. A. Third order difference methods for hyperbolic equations // J. Comput. Phys. 1970. Vol. 5, No. 3. P. 547-571.
[5] Пинчуков В. И. О неявных схемах типа Рунге — Кутты третьего порядка по времени // Вычисл. технологии. 1999. Т. 4, №2. С. 59-73.
[6] Пинчуков В. И. Абсолютно устойчивые схемы Рунге — Кутты третьего порядка // Журн. вычисл. математики и мат. физики. 1999. Т. 39, №11. С. 1855-1868.
[7] Пинчуков В. И. О численном решении уравнений сжимаемого газа неявной схемой Рунге — Кутты третьего порядка // Там же. 2002. Т. 42, №6. (В печати).
[8] Пинчуков В. И. Энтропийные и квазиэнтропийные условия для схем высоких порядков аппроксимации // Вычисл. технологии. 1997. Т. 2, №6. С. 71-87.
Поступила в редакцию 2 апреля 2002 г., в переработанном виде — 20 июня 2002 г.
