CC BY
49
УДК 53-73
А.В. Сафронов1, Ю.В. Фомин2
1 Центральный научно-исследовательский институт машиностроения 2 Московский физико-технический институт (Национальный исследовательский университет)
Метод численного решения уравнений газодинамики с помощью соотношений
на разрывах
Представлена экономичная схема Годунова типа для расчёта сложных течений газа на основе аппроксимации потоков на границе ячеек разностной сетки из приближённого безытерационного решения модельной задачи распада газодинамического разрыва с помощью соотношений на разрывах в массовых переменных с максимальной оценкой скоростей волн. Приведены результаты тестовых расчётов одномерных и двумерных задач в широком диапазоне изменения параметров в сравнении с расчётами методом Годунова.
Ключевые слова: схема Годунова, сложные разрывные течения газа, аппроксимации потоков на границе ячеек, задача распада газодинамического разрыва, максимальная оценка скоростей волн.
I. Введение
Схема С.К. Годунова [1] численного решения уравнений газодинамики без сомнения является наилучшей с точки зрения точности и надежности при сквозном расчёте сложных разрывных течений газа, включающих скачки уплотнения различной интенсивности, зоны разрежения и контактные разрывы. Метод Годунова основан на аппроксимации потоков на границах ячеек разностной сетки с помощью точного решении задачи Римана распада газодинамического разрыва. Произвольный разрыв газа распадается на три волны: на левую волну, контактный разрыв и правую волну. Левые и правые волны могут быть в зависимости от перепада давления как волнами разрежения, так и волнами сжатия (скачками). Точное решение задачи Римана требует трудоемкого решения нелинейной системы уравнений методом итераций даже в случае совершенного газа. Кроме этого, в ряде случаев, например в расчёте магнитогидро-динамических течений, решении задач с химическими реакциями в многокомпонентных потоках с переменными теплофизическими свойствами, в многофазных средах и др., точное решение задачи Римана получить затруднительно. В связи с этим широкое распространение имеют более экономичные схемы численного решения уравнений газодинамики на основе приближённого решения задачи распада разрыва с возможностями сквозного расчёта разрывных течений, приближающимися к схеме Годунова.
Рассматриваются базовые схемы первого порядка аппроксимации. Подробный обзор имеющихся схем, основанных на различных способах решения задачи Римана (риман-солверы), и их сравнение с другими подходами можно найти в работах [2-4].
II. Краткий анализ риман-солверов
1. Схема С.К. Годунова. В методе Годунова [1] аппроксимация параметров на границах ячеек разностной сетки осуществляется с помощью точного решении задачи Римана распада произвольного разрыва с параметрами, равными параметрам газа в соседних ячейках сетки. Точное решение этой задачи требует решения нелинейной системы уравнений методом итераций. В работе [1] предложен итерационный алгоритм решения этой задачи методом Ньютона, в качестве начального приближения используется акустическое (линеаризованное) приближение. Обзор способов получения точных решений задачи Римана можно найти в работе [4].
2. Схема PVRS-TSRS-TSRS. (Primitive Variable Riemann Solvers — Two-Rarefaction Riemann Solvers — Two-Shock Riemann Solvers) [2]. Схема основана на сложении трёх приближённых решений задачи Римана — линеаризованном решении (акустическое приближение), полученное давление из которого затем уточняется в приближении с двумя волнами разрежения или с двумя скачками. Практически эта схема является первой итерацией схемы Годунова, поэтому имеет осцилляции на разрывах, что затрудняет получение решения в сложных случаях.
3. Сеточно-характеристические схемы (СХМ) типа схем KUP (Курант и др.), Холодова, Рое [10] применяют для аппроксимации параметров на границах ячеек разностной сетки соотношения на характеристиках. Аппроксимация в рамках данных схем интерпретируется как линеаризованное решение задачи о распаде разрыва в конфигурации с двумя скачками. Схема Роу уникальна тем, что при определённом осреднении параметров поток на границе ячеек соответствует приближённому решению задачи о распаде разрыва с выпол-
нением соотношений на разрывах. СХМ хорошо зарекомендовали себя при расчёте течений с ударными волнами и контактными разрывами, однако имеют известные энтропийные проблемы (появление нефизических скачков) в расчёте зон разрежения при смене знака характеристик. Для устранения этих эффектов СХМ применяются процедуры «энтропийной» коррекции, обзор которых можно найти в [2, 3].
4. Двухволновая схема HLL (Harten, Lax, Leer) и её модификации с добавлением третьей волны HLLC [2]. Схема HLL основана на двухволновом приближении, без рассмотрения контактного разрыва. В схеме HLL предложен простой, но эффективный способ выбора скоростей этих волн по максимальным наклонам характеристик в соседних ячейках разностной сетки. Этот способ исключил энтропийные проблемы при смене знака характеристик. Однако не учёт условий на контактном разрыве, как показано, например, в работе [2], приводит к «размазыванию» контактного разрыва в расчёте методом установления. Для устранения этого недостатка в схеме HLLC [2] предложено несколько способов введения контактного разрыва.
5. Схема Ошера. В этой схеме задача Рима-на решается в изоэнтропическом приближении при замене ударной волны пучком сходящихся характеристик. Схема не имеет осцилляций на разрывах (сильный разрыв «размазывается» на 2-3 ячейки сетки) и проблем прохождения звуковой точки. Однако в связи с допущением о постоянстве энтропии эта схема имеет ограничение на область применимости в случае интенсивных скачков, кроме того, объём вычислений приближается к точному решению задачи Римана. К тому же схема требует аналитических соотношений на характеристиках, что затрудняет её обобщение на случай переменных термодинамических свойств газа.
Схемы п. 2-5 основаны на различных способах приближённого решения задачи Римана, в работе [5] предложен новый способ приближённого решения этой задачи на основе выполнения соотношений на разрывах.
III. Разностный метод для уравнений газодинамики из соотношений на разрывах в массовых переменных
Рассмотрим описание схемы на примере двумерных уравнений гидродинамики в расщеплении по оси x:
Ut + Fx = 0, (1)
U = {ppupvpE }T,
F = {pu,pu2 + P,puv,puE + Pu}T,
где p — плотность газа, u — компонента скорости по оси x (нормальная к грани ячейки), v — ком-
понента скорости, касательная к грани ячеики), Р — давление, Е = Р/— 1) + V?/2 — полная внутренняя энергия на единицу объёма, 7 — показатель адиабаты.
Разностную схему запишем в консервативном виде:
= (2)
где п — номер шага по времени с интервалом ДЬ, г — номер ячеИки сетки по оси х с разбиением Дх.
В методе Годунова [1] поток на границе ячеек -^+1/2 вычисляется из решения задачи распада произвольного разрыва с состояниями, соответствующими параметрам газа в соседних ячеИках сетки. При этом схема устоИчива при условии Куранта: СЕЬ = и>тахД£/Дх < 1, где штах — максимальная скорость распространения возмущениИ.
Предложенный в работе [5] метод состоит в расчёте потока на границе ячеек путём приближённого решения задачи Римана с использованием соотношениИ на поверхности разрывов:
AF = wAU,
(3)
где и> — скорость распространения разрыва. Символ Д обозначает разность значениИ соответствующей величины на разрыве.
1+1/2 в)
Рис. 1. (а) Левая ударная волна. (Ь) Левая волна разрежения. (с) Распад разрыва
Основная идея метода [5] состоит в применении выражениИ законов сохранения массы, им-
пульса и энергии на поверхности газодинамического разрыва (3) в виде, соответствующем уравнениям (1) в лагранжевых (массовых) переменных:
Дш = 0, (4)
шДи + ДР = 0, шДЕ + Д(Ри) = 0, шД^ = 0,
(5)
(6) (7)
где ш = р(п — и>) — массовая скорость (плотность потока газа) через поверхность разрыва. На контактном разрыве ш = 0.
В дальнейшем условимся считать ш положительным.
Рассмотрим подробнее «левую» волну (или волну, обращенную влево), поток газа через неё вдоль оси х направлен слева направо (рис. 1а, 1Ь).
Рассмотрим некоторые свойства (см., например, [6]) перехода из состояния газа до разрыва «1» в состояние после, которое обозначим параметрами без индекса.
Введём обозначение с — скорость звука.
В случае ударной волны ^В) (схема на рис. 1а) справедливы следующие неравенства:
рС1 > ш > Р1С1, П1 — С1 > IV > и — с.
Подставляя максимальное значение скорости распространения YB в выражение для массовой скорости, получим
рс> Р1(и1 — и + с) > ш > Р1С1.
В случае, если левая волна является волной разрежения, максимальное значение скорости распространения возмущения определяется характеристикой и1 — с1 , а максимальное значение массовой скорости определяется по состоянию «1» (рис. 1Ь). Таким образом, максимальное значение массовой скорости левой волны определяется по формулам
штаХ = тах(р1С1,рс),
или
штах = Р1(п1 — тт(п1 — с1,и — с)). (8)
Для правой волны оценки аналогичны. В дальнейшем будем использовать максимальные оценки изменения массовой скорости на разрыве.
С помощью уравнений (4)-(7) решаем задачу Римана с распадом на левую волну, контактный разрыв и правую волну. Ниже обозначим индексами: 1 — параметры в левой ячейке сетки, 2 — в правой, 3 — между левой волной и контактным разрывом, 4 — между контактным разрывом и правой волной, как показано на схеме течения рис. 1с.
Решение имеет вид
из = П4
№2 = и2 + ш2/р2, (10)
П2Ш2 + П1Ш1 — Р2 + Р1
ш1 + ш2
(11)
Рз = Р4 = Р =
Р2 /77-1 + Р\ 2 ~ ~ 1)
ш1 + ш2
■№1 = и1 — ш1/р1,
(9)
(12)
1/рз = 1/р1 + (и — и\)/ш1, (13)
Ез = Е1 — (Ри — —Рщ1)/ш1, (14) 1/р4 = 1/р2 — (и — и1)/ш2, (15)
Е4 = Е2 — (Ри — —Р2П2)/ш2, (16)
Уз = г>1,г>4 = У2.
Здесь индексы соответствуют номеру зоны на схеме течения рис. 1с.
Когда скорость распространения контактного разрыва и ^ 0: при скорости левой волны №1 ^ 0 параметры на границе ячеек равны параметрам в зоне 1, а при №1 < 0 — параметрам в зоне 3. В случае и < 0: при скорости правой волны №2 ^ 0 параметры на границе ячеек равны параметрам в зоне 2, а при №2 > 0 — параметрам в зоне 4.
Характерная особенность предложенной схемы заключается в том, что, поскольку на контактном разрыве ш = 0, решение задачи зависит от двух параметров — массовых скоростей на левой и правой волне ш1 и ш2. Причём в сеточно-харак-теристических схемах, так же, как в «звуковом» приближении, соотношения (5)-(7) так или иначе рассматриваются в линеаризованном виде.
На основе оценок (8) в работе [7] рекомендован следующий алгоритм выбора массовых скоростей в задаче распада разрыва (рис. 1с), отвечающий условию неубывания энтропии:
1) вычисляются массовые скорости левой и правой волны по максимальным наклонам характеристик:
ш1 = р1 тах(с1,с2 + и1 — и2),
ш2 = р2 тах(с2,с1 + и1 — и2), (17)
этот выбор обеспечивает выполнение энтропийного условия ш1 ^ р1с1, ш2 ^ р2с2 в зонах разрежения, а также условия №2 > №1, при этом скорость волны определяется характеристикой перед волной;
2) в случае конфигурации с ударными волнами, когда оценки по массовой скорости (8) превышают оценку по наклонам характеристик:
ш1 = р2с2,при^с2 > р1 (с2 + и1 — и2) > р1 с1;
ш2 = р1с1,при1 с1 > р2 (с1 + и1 — и2) > р2 с2. (18)
С выбором скоростей (12)-(17) схема положительно определена по плотности (13), (15) и энергии
(14), (16)[7].
В случае интенсивных ударных волн Р ^ ж, (и1 — и2) ^ ж соотношения (17) при вычислении давления по формуле (12) приводят к квадратичной схемной вязкости
Р
Р1Р2 Р1+Р2
(и1 — и2)2. Это позволяет вести расчёт
и
—>
интенсивных ударных волн с шагом по времени, близким к условию Куранта [4].
В случае сильного разрежения для параметра P, вычисленного по формуле (12), энтропийно обосновано [7] применение малого ограничения: P = max(P,fix), где fix = 1.E — 6.
Таким образом, предлагаемое приближённое решение задачи Римана из соотношений на разрывах вектора потока (3) заключается в формулах (9-22), по которым вычисляются «простые» переменные и по ним в свою очередь вектор потока на границе ячеек ^¿+1/2. Полученный поток используется затем в разностном уравнении (2) аналогично схеме Годунова.
Представленная схема хорошо зарекомендовала себя при решении различных одномерных и пространственных задач и обобщена на сверхзвуковой стационарный случай [8].
IV. Результаты тестовых расчётов
В качестве иллюстрации ниже представлены расчёты нестационарных тестовых задач, собранных в работе [2], и расчёты одномерных течений в сопле методом установления. Результаты расчётов по новой схеме сопоставляются с результатами расчётов по схеме Годунова.
Исходные данные задач [2] приведены в табл. 1, где xo — начальное положение разрыва, £ — время. Начальные условия: в момент времени £ = 0 при х < хо параметры газа задавались равными параметрам с индексом 1, при х > хо — параметрам с индексом 2. Число разбиений N = 50, Дх = 1/50, СЕЬ = 0,9.
Результаты расчётов представлены на рис. 2, 3, 4, 5 и 6, на которых приведены точное решение задач и результаты расчётов с 1-м порядком по новой схеме, условное название схема «С», и схеме Годунова.
Таблица 1
№ rl ul PI r2 u2 Р2 хО t
Тест 1 1 0,75 1 0,125 0 0,1 0,3 0,2
Тест 2 1 -2 0,4 1 2 0,4 0,5 0,15
Тест 3 1 0 1000 1 0 0,01 0,5 0,012
Тест 4 5,99924 19,5975 460,894 5,99242 -6,19633 46,095 0,4 0,035
Тест 5 1 -19,59745 1000 1 -19,59745 0,01 0,8 0,012
Из рисунков видно, что представленная схема не имеет осцилляций на разрывах, нет проблем прохождения звуковой точки и не требуется снижение числа Куранта на сильных разрывах. В зонах разрежения и при ударных волнах слабой и «средней» интенсивности результаты очень близки к схеме Годунова, сильный разрыв «размазывается» на 2-- 3 ячейки сетки, в то время как в схеме Годунова на 1--2.
Заметим, что сеточно-характеристическая схема Рое [9] не проходит тест 2, см. [2].
С целью сравнения сходимости новой схемы и схемы Годунова проведены расчёты одномерного течения в сопле методом установления.
Геометрия сопла задавалась следующим образом: у = ут(х/5 — 1)2 + 1, где
( 1,964, 0 <х < 5, Ут = \ 0,6875, 5 < х < 10.
Расчёты проводились по схемам 1-го и 2-го порядков аппроксимации. Повышение порядка схем проводилось методом [10] с ограничителем тштс^
Расчёты проводились при числе СЕЬ = 0,99.
На рис. 7, 8 приведены результаты расчёта течения в сопле без скачка. Число ячеек сетки
N = 50, Дх = 10/50. Исходные параметры в начальном и выходном сечениях сопла приведены в табл. 2. Рассматривались два режима течения в сопле без скачка и со скачком.
На рисунках показаны распределения в ячейках сетки вдоль сопла числа Маха, плотности и давления, а также зависимость логарифма параметра сходимости log (rms) от числа итераций с шагом 100.
Параметр сходимости rms определялся как отношение максимального по сетке приращения плотности по времени на n-м шаге счета по времени ((5рП) = рП+1 — Рп ) к максимальному по сетке приращению плотности по времени на 1-м шаге счета:
rms = max(5prn)/max(^p').
Общее число итераций принималось 6000.
На рис. 7 приведён расчёт по схемам 1-го порядка, на рис. 8 — по схемам 2-го порядка.
Из данных рис. 7, 8 видно, что на этом тесте данные расчётов и сходимость схем практически совпадают.
Характер зависимостей 7, 8 практически не меняется при изменении числа ячеек сетки.
Таблица 2
Плотность р Давление Р Число Маха М
входное сечение 1,55 1,32 0,22
выход без скачка 0,363 0,173 2
выход со скачком 1,363 1,13 0,404
Рис. 8. Течение в сопле 1М = 50. 2-й порядок
0 2000 4000 6000 iter
Рис. 10. Течение в сопле со скачком 1М = 50. 2-й порядок
0 2000 4000 6000 iter
Рис. 12. Течение в сопле со скачком 1М = 52. 2-й порядок
На рис. 9, 10, 11, 12 приведены данные расчётов течения в сопле со скачком.
На рис. 9, 10 показаны схемы 1-го и 2-го порядка при числе узлов N = 50. На рис. 11, 12 N = 52. Разное число разбиений N обеспечивает различное положение скачка относительно ячеек сетки. При дальнейшем изменении N зависимости 9-12 практически не меняются.
Из рис. 9-12 видно, что сходимость предложенной схемы при течении со скачком близка к схеме Годунова при несколько большем на 1 ячейку сетки) «размазывании» разрыва.
Аналогичные результаты имеют решения пространственных задач.
Так, на рис. 13 показаны сравнительные результаты расчётов в виде изолиний чисел Маха при взаимодействии затопленной осесимметрич-ной струи с преградой. Параметры на срезе сопла: число Маха М а = 2,5, отношение давления к внешнему Ра /Ре = 15, 7 = 1,4, температура торможения равна внешней температуре. На рис. 14 приведены распределения давления и числа Маха на оси струи и вдоль преграды. Расчёты проводились по схемам 2-го порядка, на квадратной сетке 60 х 105 ячеек. Граничные условия: на оси симметрии и на стенке — условия отражения, на срезе сопла — параметры истечения, на свободных границах параметры определелись из решения задачи распада разрыва с начальными данными, соответствующими параметрам внешней среды и параметрам в ближайшей ячейке сетки. В обоих слу-
чаях через 5000 итераций параметр тшз имел порядок 10~6. Область течения включает зоны разряжения, скачки, контактный разрыв и переход от дозвукового течения к сверхзвуковому (вдоль преграды), которые рассчитываются сквозным образом. Как видно, результаты расчётов предложенным методом и методом Годунова согласуются между собой. В этих расчётах время счета по предложенной схеме со 2-м порядком было в 2 раза меньше, чем по схеме Годунова, а с 1-м — в 3,5 раза. Следует отметить, что при расчёте со 2-м порядком, кроме решения задачи Римана, требуются вычисления на промежуточном шаге с процедурами экстраполяции параметров, объём которых одинаков для рассматриваемых подходов.
Выводы. Приведены результаты тестирования новой экономичной схемы расчёта сложных разрывных течений газа на основе приближённого решения задачи Римана с помощью соотношений на разрывах в массовых переменных с максимальной локальной оценкой скоростей волн.
Предложенный метод не имеет проблем, присущих другим схемам на основе приближённого решения задачи Римана. Представленная схема не имеет осцилляций на разрывах, нет проблем прохождения «звуковой точки» и не требуется снижение числа Куранта на сильных разрывах.
Результаты численных расчётов зон разрежения и скачков близки к схеме Годунова при сокращении времени счета типовых двумерных задач в 2-- 3 раза.
Рис. 13. Изолинии числа Маха. M = 2,52, n = 15. Сетка 60 х 105. 2-й порядок
итерации
Рис. 14. Ма = 2.52, ва = 7°, п = 15
Литература
1. Годунов С.К. Разностный метод численного расчёта разрывных решений уравнений гидродинамики // Матем. сб. — 1959. — Т. 47, вып. 3.
2. Toro E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics. — Springer-Verlag. Second Edition, June 1999.
3. Куликовский А.Г., Погорелое Н. В, Семенов А.Ю. Математические вопросы численного ре-
шения гиперболических систем уравнений. — М.: Физматлит, 2001.
4. Чарахчьян А.А. Об алгоритмах расчёта распада разрыва для схемы С.К. Годунова // ЖВМ и МФ. — 2000. — Т. 4, № 5.
5. Сафронов А.В. Разностный метод решения нестационарных уравнений газодинамики на основе соотношений на разрывах // Космонавтика и ракетостроение. — 2006. — № 2 (43). — С. 152-158.
6. Овсянников Л.В. Лекции по основам газовой динамики. — М.: Наука. 1981.
7. Сафронов А.В. Кинетические схемы для уравнений газодинамики / / Вычислительные ме-
тоды и программирование. — 2009. — Т. 10. — С. 62-74.
8. Сафронов А.В. Разностный метод для уравнений газодинамики из соотношений на разрывах // Математическое моделирование. — 2008. — Т. 20, № 2. — С. 76-84.
9. Roe P.L. Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes //J. Comput. Phis. — 1981. — N 43.
10. Родионов А.В. Повышение порядка аппроксимации схемы С.К. Годунова // ЖВМ и МФ. — 1987. — Т. 27.
Поступила в редакцию 14.09.2009.
