ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2009, том 49, < 3, с. 549-566
УДК 519.634
ВОЗМОЖНОСТИ КВАЗИГАЗОДИНАМИЧЕСКОГО АЛГОРИТМА ДЛЯ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ТЕЧЕНИЙ
НЕВЯЗКОГО ГАЗА
© 2009 г. Т. Г. Елизарова, Е. В. Шильников
(125047 Москва, Миусская пл., 4а, ИММ РАН) e-mail: telizar@mail.ru; shiva@imamod.ru
Поступила в редакцию 12.08.2008 г. Переработанный вариант 16.09.2008 г.
Представлено решение десяти известных одномерных тестовых задач, отражающих характерные особенности нестационарных течений невязкого газа. Показано, что все тесты успешно решаются с помощью единого численного алгоритма, основанного на квазигазодинамической системе уравнений. Во всех случаях численное решение сходится к автомодельному решению при сгущении пространственной сетки. Библ. 15. Фиг. 27. Табл. 1.
Ключевые слова: квазигазодинамический алгоритм, уравнения Эйлера, одномерные течения.
ВВЕДЕНИЕ
Расчет нестационарных течений сжимаемого газа в одномерном по пространству приближении является наглядным способом оценки возможностей численных алгоритмов. В настоящее время имеется большое число примеров решения многомерных задач, которые также используются для сопоставления численных методов между собой. Однако в таких расчетах, как правило, точность решения определяется не только самим численным алгоритмом, но и структурой пространственной сетки, способом постановки граничных условий и прочими обстоятельствами, что затрудняет сопоставление самих численных алгоритмов. Кроме того, одномерные течения невязкого газа детально исследованы теоретически, а для течений типа распада разрыва имеются автомодельные решения уравнений Эйлера (см. [1], [2]), которые служат надежным эталоном проверки точности и сходимости численных решений.
Представительная система тестов собрана, например, в [3] и [4], где десять одномерных нестационарных течений невязкого газа рассчитываются на основе десяти численных алгоритмов. Эти алгоритмы достаточно полно представляют возможности конечно-разностного подхода к решению уравнений газовой динамики в эйлеровых переменных.
В данной работе на этих же тестах проверяются возможности численного алгоритма, основанного на квазигазодинамической (КГД) системе уравнений (см., например, [5]-[7]). КГД-алго-ритм и родственные ему кинетически-согласованные разностные схемы (см. [8]) успешно использовались для численного моделирования широкого круга течений вязкого сжимаемого газа. С их помощью рассчитывались двумерные и трехмерные нестационарные течения, но при этом тестовым расчетам одномерных задач уделялось недостаточное внимание. Настоящая работа восполняет этот пробел и демонстрирует применимость КГД-алгоритма для сквозного расчета разнообразных по своей природе нестационарных одномерных задач.
1. СИСТЕМА КВАЗИГАЗОДИНАМИЧЕСКИХ УРАВНЕНИЙ И ЧИСЛЕННЫЙ АЛГОРИТМ
Квазигазодинамические (КГД) уравнения для одномерного плоского течения тых обозначениях имеют вид
dP + j = 0,
д t д x
общеприня-(1)
Эр и + д ]ти + Эр _ ЭП д г д х дх дх'
дЕ + д)тН + дч _ ЭПи д г дх дх дх
(3)
Здесь Е и Н - полная энергия единицы объема и полная удельная энтальпия, которые вычисляются по формулам Е = ри2/2 + р/(у - 1) и Н = (Е + р)/р. Вектор плотности потока массы вычисляется как
)т _ Р(и - ),
где
Т д , 2 .
" _ рэх(ри + р)•
Компонента тензора вязких напряжений, входящая в систему уравнений (1)-(3), определяется таким образом:
„4 ди Г Эи дрЛ Г др диЛ
П _ з цэх + иТ1р и эх + др) + 4и дх + эхЛ •
Вектор теплового потока ч определяется по формуле
дТ
9 _-кэх "хри
и Э Г р
д Г1
_у - 1 д х V р у + рид х V Р
где р = рЯТ, ц - коэффициент динамической вязкости, к = цуЯ/((у - 1)Рг) - коэффициент теплопроводности, у - показатель адиабаты, Рг - число Прандтля, т = ц/(р Бе) - релаксационный параметр, имеющий размерность времени, Бе - число Шмидта.
Для удобства численного решения система уравнений (1)-(3) приводится к безразмерному виду с использованием базовых значений плотности р0, скорости звука с0 = ^КТ0 и длины Ь. Обез-размеривание не изменяет вид уравнений.
Введем равномерную сетку по координате х с шагом Н, и сетку по времени с шагом Дг. Значения всех газодинамических величин - скорости, плотности, давления - будем определять в узлах сетки. Значения потоков определяются в полуцелых узлах. Для решения задачи (1)-(3) используем явную по времени разностную схему следующего вида:
л Дг(.
Р; _ Р; - Н (]т, ; + 1/2 — ]п
-1/2
),
Дг
Р;и; _ Р;и; + н [(П" +1/2 - П;-1/2) - (р; +1/2 - р;-1/2) - От,; + 1/2и; + 1/2-]т,;-1/2и;-1/2)]>
(4)
(5)
Е; _ Е +дг
(П; +1/2и; +1/2 - П; -1/2и; -1/2
) - (+1/2 - -1/2) -
]т,; + 1/2 , ^ ч ] т,; - 1/2 / т-т ч
-(Е; + 1/2 + р; + 1/2) - —-(Е;-1/2 + р;-1/2)
; + 1/2
; - 1/2
2
р* _ (Г -1 >(Е
Дискретный аналог потока массы ]т имеет вид
]т,; +1/2 _ Р; +1/2(и; + 1/2 - +1/2) , где добавка к скорости вычисляется таким образом:
Т
; + 1/2
; + 1/2 1 2 2
7 ( Р; +1 и; + 1 + р; + 1- Р;и; - р;) •
Р; +1/2 Н
(6)
(7)
Таблица
Тест Pl ul Pl Pr Ur Pr fin
1 1 0.75 1 0.125 0 0.1 0.2
2 1 -2 0.4 1 2 0.4 0.15
3 1 1 10-6 1 -1 10-6 1
3a 1 -19.59745 1000 1 -19.59745 0.01 0.012
4 5.99924 19.5975 460.894 5.99924 -6.19633 46.095 0.035
5 1.4 0 1 1 0 1 2
6 1.4 0.1 1 1 0.1 1 2
7 0.1261192 8.9047029 782.92899 6.591493 2.2654207 3.1544874 0.0039
Дискретные выражения для П и q выписываются аналогично. Порядок точности разностной схемы (4)-(8) составляет O(h2 + At).
При численном решении уравнений Эйлера на основе системы (1)-(3) все диссипативные слагаемые, т.е. слагаемые с коэффициентами ц, к и т, рассматриваются как искусственные регуля-ризаторы. При этом релаксационный параметр и коэффициенты вязкости и теплопроводности связаны между собой и в безразмерном виде вычисляются по формулам
h о т р Sc
т = а -, ц = трSc, к = рг { _ 1), (9)
где а - числовой коэффициент, который, как правило, выбирался в пределах 0.2-0.7. В большинстве представленных далее расчетов полагаем Pr = 1 и Sc = 1.
Схема (4)-(9) формально имеет порядок O(ah + At). Приведенные далее расчеты подтверждают, что уменьшение коэффициента а в определенных пределах эквивалентно сгущению пространственной сетки в а раз.
Выписанная разностная схема (4)-(9) обладает условием устойчивости Куранта. Шаг по времени выбирается из соотношения
A t = в min f, ,h J, (10)
V| Щ + c J
где в - числовой коэффициент, который в большинстве расчетов составляет 0.1-0.4.
2. ЗАДАЧИ РИМАНА О РАСПАДЕ РАЗРЫВОВ
В этом разделе рассмотрены задачи о распаде разрывов, собранные в [3] и [4]. Эти задачи всесторонне отражают характерные и сложные для численного моделирования особенности нестационарных газодинамических течений. Начальные данные к задачам о распаде разрывов приведены в таблице в соответствии с обозначениями, принятыми в [3] и [4]. А именно, значения газодинамических величин слева от разрыва обозначены индексом L, справа - индексом R. Момент времени, для которого построены графики, указан в таблице и обозначен как tfin.
Граничные условия совпадают с соответствующими начальными условиями на концах расчетной области. Во всех вариантах расчетов у = 1.4, за исключением задачи Ноха (3), для которой Y = 5/3. Длина области расчета равна 1, разрыв расположен в точке 0.
Задачи 1, 2, 4, 5 и 6 впервые были предложены и численно решены в [9]. В [10] задача 2 решалась на основе КГД-алгоритма, представленного в виде схемы расщепления потока. Задача Ноха была впервые предложена и численно решена в [11].
Тест 1. Данная задача представляет собой один из вариантов задачи Сода. В образующемся течении имеются все характерные особенности, присущие сверхзвуковым течениям: звуковые точки на границах волны разрежения, контактный разрыв и ударная волна.
В рамках КГД-алгоритма устойчивое решение этой задачи обеспечивается выбором коэффициентов 0.2 < а < 0.5 и в = 0.4 в формулах (9) и (10), соответственно. Уменьшение коэффициента а
Фиг. 1.
Фиг. 2.
при расчете параметра регуляризации (9) приводит к уточнению численного решения. На фиг. 1 показаны распределения плотности, вычисленные для пространственной сетки с шагом к = = 0.0025. Сплошная линия - автомодельное решение. Здесь, также как и в последующих расчетах, уменьшение коэффициента а в два раза практически эквивалентно уменьшению шага пространственной сетки в два раза. Однако выбор а < 0.2 приводит к возникновению вычислительной неустойчивости, которая может быть устранена путем уменьшения шага по времени -уменьшения коэффициента |3. Сходимость численного решения к автомодельному решению задачи (сплошная линия) при сгущении пространственной сетки для а = 0.2 демонстрирует фиг. 2.
Тест 2. Течение здесь представляет собой две волны разрежения, разбегающиеся от центра области. Сложность численного решения этой задачи обусловлена тем, что в центре между разбегающимися потоками плотность газа и его давление очень малы, но в то же время внутренняя энергия 8 = р/[р(у - 1)] к нулю не стремится. Представляется, что никакая разностная схема в пе-
Фиг. 3.
Фиг. 4.
ременных Эйлера не описывает поведение внутренней энергии в этой задаче с высокой точностью (см., например, [3] и [4]).
Во всех расчетах этой задачи коэффициент в = 0.1. Увеличение шага по времени приводило к численной неустойчивости. Зависимость плотности и внутренней энергии газа от шага сетки и параметра регуляризации а приведены на фиг. 3 и 4. Линия 6 соответствует автомодельному решению. Где не оговорено специально, а = 0.5. Уменьшение параметра регуляризации, также как и уменьшение шага сетки, приводит к уточнению численного решения и его приближению к автомодельному решению. Это особенно видно на графиках внутренней энергии (фиг. 4). Результаты расчетов демонстрируют быструю сходимость распределения скорости по сетке.
Тест 3. Задача Ноха. Рассматриваемое в этой задаче течение представляет собой столкновение двух гиперзвуковых потоков холодного плотного газа, которое приводит к образованию двух расходящихся "бесконечно сильных" ударных волн, между которыми остается неподвижный газ с постоянными плотностью и давлением. Действительно, согласно начальным условиям
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.