CC BY
21
Численное моделирование течений газа на основе квазигидродинамических уравнений
Т. Г. Елизарова1,fl, О. В. Булатов2
1 Институт математического моделирования РАН. Россия, 125047, Москва, Миусская пл., д. 4а.
2 Московский государственный университет имени М. В. Ломоносова, физический факультет,
кафедра математики. Россия, 119991, Москва, Ленинские горы, д. 1, стр. 2.
E-mail: а telizar@mail.ru
Статья поступила 10.06.2009, подписана в печать 08.07.2009
На примере семи известных тестовых задач, отражающих характерные особенности нестационарных течений невязкого газа, исследуются возможности квазигидродинамического алгоритма для численного решения системы уравнений Эйлера.
Ключевые слова: квазигидродинамический алгоритм, уравнения Эйлера, одномерные течения. УДК: 517.958:533.7. PACS: 47.ll.Df; 47.40.x.
Введение
Численный расчет нестационарных течений сжимаемого газа в одномерном по пространству приближении является известным способом оценки возможностей численных алгоритмов.
Одномерные по пространству течения идеального газа к настоящему времени хорошо изучены, а для задач типа распада разрыва имеются автомодельные решения уравнений Эйлера [1, 2], которые служат надежным эталоном для проверки точности и сходимости численных решений. Представительная система тестов собрана в работах [3, 4], где расчеты одномерных нестационарных течений идеального газа проведены на основе десяти различных численных алгоритмов. Эти алгоритмы достаточно полно представляют возможности конечно-разностного подхода к решению уравнений газовой динамики в эйлеровых переменных.
В настоящей работе на этих тестах проверяются возможности численного алгоритма, основанного на квазигидродинамической системе уравнений, построенной в работах Ю. В. Шеретова (см., напр., [5, 6]).
Квазигидродинамические уравнения, или уравнения Шеретова, могут рассматриваться как система, описывающая течения вязкого сжимаемого газа и обобщающая уравнения Навье-Стокса. По сравнению с системой Навье-Стокса система Шеретова включает в себя регуляризирующие добавки к уравнениям в виде вторых пространственных производных с малым коэффициентом в качестве множителя. Упрощение этих уравнений, полученное для случая течений вязкой несжимаемой жидкости, позволило построить эффективные численные алгоритмы решения нестационарных задач вынужденной и свободной конвекции. Однако использование уравнений Шеретова для расчета сжимаемых течений до настоящего времени практически не изучалось. Исключение составляет численное решение задачи о поршне [5], в которой демонстрируется работоспособность указанной модели.
Система квазигазодинамических уравнений, послужившая основой для системы Шеретова, позволила построить семейство численных алгоритмов, чрезвычайно эффективных для расчетов течений вязкого сжимаемо-
го газа [5-7]. Однако указанная система существенным образом опирается на уравнение состояния идеального газа, что является неадекватным для многих практических приложений. Поэтому представляется актуальным изучение возможностей численных алгоритмов расчета сжимаемых течений, основанных на системе Шеретова, для которой уравнение состояния имеет более общий вид.
1. Система уравнений и численный алгоритм
Квазигидродинамическая система уравнений Шеретова для одномерного плоского течения вязкого газа в общепринятых обозначениях имеет вид
9Р | d/m = Q dt дх ' дри | djmu | др _ ОТ
~дГ ' '
(1) (2) (3)
дх дх дх' дЕ д]тН <9^ ОТ и д1 дх дх дх ' Здесь Е и Я — полная энергия единицы объема и полная удельная энтальпия, которые вычисляются по формулам: Е = ри2/2 + ре и Н = (Е + р)/р. Вектор плотности потока массы вычисляется как
¡т = р(и - т), где добавка к скорости имеет вид
Ш р \идх дх) '
Компонента тензора вязких напряжений, входящая в систему уравнений (1)—(3), определяется в виде
„4 ди
Вектор теплового потока д равен
дТ
где р — коэффициент динамической вязкости, к = р"/Я/((1) Рг) — коэффициент теплопроводности, 7 — показатель адиабаты, Рг — число Прандтля,
т = р/(р Sc) — релаксационный параметр, имеющий размерность времени, Sc — число Шмидта.
Система уравнений (1)-(3) допускает замыкание вида
р = р(р,Т), е = е(р,Т). (4)
Однако в дальнейших расчетах в целях сопоставления с уже имеющимися численными результатами в качестве (4) мы будем использовать уравнение состояния идеального газа
Р
р = pRT, е ■■
Для удобства численного решения система уравнений (1)-(3) приводится к безразмерному виду с использованием базовых значений плотности р0, скорости звука со = и длины Ь. Обезразмеривание не
изменяет вид уравнений.
Введем равномерную сетку по координате х с шагом /г, и сетку по времени с шагом А1. Значения всех газодинамических величин — скорости, плотности, давления — будем определять в узлах сетки. Значения потоковых величин определяются в полуцелых узлах. Величины с «крышкой» относятся к верхнему временному слою. Для решения задачи (1)-(3) используем явную по времени разностную схему, впервые предложенную в [5] для численного решения задачи о поршне:
Pi— Pi--Jj-(jmi+1/2 — /mi—1/2).
P1U1 = P1U1 + — [(Пг+1/2 - П^1/2) - (pi+1/2 - pi-1/2)
— ijmi+l/2Щ+1/2 -/ш-1/2"|-1/2)]. At
(5)
(6)
Ei=Ei +
h
(<?i+ 1/2 /mi—1/2
(П(+1/2«г+1/2 - П(_i/2«i-l/2) -
- 9/-1/2) _ ( (Ei+l/2 + Pi+1/2)' \ Pi+1/2
Pi—1/2
№-1/2 + Pi-1/2)
(7)
Р1Щ 2
Дискретный аналог потока массы /т имеет вид
/т¿+1/2 = РН1/2(иН1/2 ~ 1/2). где добавка к скорости вычисляется как
* р1+\ р1
_ Ti+l/2 ( Щ+1
wi+1/2 — - [ Pi+1/2 Щ+1/2
Pi+1/2
h
h
(8)
(9)
Дискретные выражения для П и ^ выписываются аналогично. Порядок точности разностной схемы (5)-(9) составляет О (/г2 + ДО.
При численном решении уравнений Эйлера на основе системы (1)-(3) все диссипативные слагаемые, т.е. слагаемые с коэффициентами р, к и т, рассматриваются как искусственные регуляризаторы. При этом релаксационный параметр и коэффициенты вязкости и теплопроводности связаны между собой и в безразмерном виде вычисляются как
Н тр $>с
т = а—, р = тр Эс, к = ■=—-тт, (10)
с Рг(7 - 1)
где а > 0 — численный коэффициент. В представленных далее расчетах Рг = 1 и 5с = 1.
Схема (5)—(10) формально имеет порядок 0(оА + ДО. Приведенные далее расчеты подтверждают, что уменьшение коэффициента а в определенных пределах эквивалентно сгущению пространственной сетки в а раз.
Выписанная разностная схема (5)—(10) обладает условием устойчивости Куранта. Шаг по времени выбирается из соотношения
At = 0 min ( -с ) ,
(П)
где /3 > 0 — численный коэффициент, или число Куранта.
2. Задачи Римана о распаде разрывов
В этом разделе рассмотрены задачи о распаде разрывов, собранные в [3, 4]. Эти задачи всесторонне отражают характерные и сложные для численного моделирования особенности нестационарных газодинамических течений. Начальные данные к задачам о распаде разрывов приведены в таблице в соответствии с обозначениями, принятыми в [3, 4], а именно значения газодинамических величин слева от разрыва обозначены индексом Ь, справа — индексом Я. Момент времени, для которого построены графики, указан в таблице и обозначен как .
Граничные условия совпадают с соответствующими начальными условиями на концах расчетной области. Во всех вариантах расчетов 7= 1.4, за исключением задачи Ноха (3), для которой 7 = 5/3. Длина области расчета равна 1, разрыв расположен в точке 0. Тестирование квазигазодинамических уравнений на этих примерах представлено в [8].
Тест 1. В этой задаче имеются характерные особенности, присущие сверхзвуковым течениям, — волна разрежения, контактный разрыв и ударная волна.
Начальные условия для задач Римана
Тест 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^в 1 -1 10^B 1
За 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
Р
0.7 0.6 0.5 0.4 0.3 0.2 0.1
0
■ Тест 1 Тест 2 Тесты 5 и 6
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 а
б
1.0 ч0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1
■А = 0.0010 h = 0.0025 h = 0.0002 ■ автомодельное решение
-0.2 -0.1
0
0.1
0.2
0.3
0.4 0.5 х
Рис. 1. а — Тест 1. Условие устойчивости алгоритма для тестовых задач 1, 2 и 5, 6. б — Распределение плотности р. Зависимость решения от шага сетки, а = 0.4, /3 = 0.1. Пунктирная линия: /г = 0.001; штриховая: /г = 0.0025; сплошная: /г = 0.0002; серая — автомодельное решение
Рассмотренный алгоритм обеспечивает устойчивое решение этой задачи при согласованном выборе коэффициентов а и /3 в формулах (10) и (11) соответственно. На рис. 1 ,а показана зависимость максимального значения параметра /3, определяющего шаг по времени для устойчивого счета, от величины параметра регуляризации а для тестовой задачи 1 и последующих примеров.
Сходимость численного решения к автомодельному решению задачи (серая линия на этом и последующих рисунках) при сгущении пространственной сетки для а = 0.4, /3 = 0.1 демонстрирует рис. 1,6. Здесь величина параметра регуляризации а вдвое превышает соответствующее значение для квазигазодинамического алгоритма [8].
Тест 2. В этой задаче течение представляет собой две волны разрежения, разбегающиеся от центра области. Сложность численного решения этой задачи обусловлена тем, что в центре между разбегающимися потоками плотность газа и его давление очень малы, но в то же время внутренняя энергия е = р/{р{7—1)) к нулю не стремится. Следует отметить, что все известные разностные схемы в переменных Эйлера неудовлетворительно описывают поведение внутренней энергии в этой задаче (см., напр., [3, 4]).
На рис. 2 представлены распределения плотности и внутренней энергии в этой задаче для а = 0.5, /3 = 0.1. Видна уверенная сходимость решения к авто-
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 х
-0.5 -0.4 -0.3 -0.2 -0.1
Рис. 2. а — Тест 2. Распределение плотности р. Сходимость решения по сетке, а = 0.5, /3 = 0.1. б — Распределение внутренней энергии е. Сходимость по сетке. Сплошная линия: h = 0.005; штрихпунктир-ная: h = 0.002; штриховая: h = 0.001; пунктирная: h = 0.0001; серая — автомодельное решение
модельной зависимости при сгущении сетки для графика плотности и внутренней энергии. Однако, в отличие от результатов [8], даже в расчетах на самой подробной сетке в распределении внутренней энергии вблизи центра области сохраняется максимум, наличие которого противоречит автомодельному решению и представляет собой так называемый энтропийный след. Область устойчивого счета этой задачи представлена на рис. 1 ,а.
Тесты 3 (Задача Ноха) и За. Первый тест представляет собой столкновение двух гиперзвуковых потоков холодного плотного газа, которое приводит к образованию двух расходящихся «бесконечно сильных» ударных волн, между которыми остается неподвижный газ с постоянными плотностью и давлением. Действительно, согласно начальным условиям (таблица), скорость звука на невозмущенном фоне составляет с = л/^Pr/pr = 0.0013. Скорость распространения волн равна 1, т.е. число Маха Ма = ¿¿¿/с = 775. Известно, что в земных условиях максимально достижимое число Маха составляет порядка 30. Второй тест описывает газодинамическое течение типа сжатия газа в термоядерной мишени. Перепад давления Pl/Pr составляет 105, что соответствует перепаду температур того же порядка.
Разностная схема (5)—(10) не позволяет решить данные тестовые задачи. Отсюда следует, что имеющийся
в системе уравнении регуляризатор оказывается недостаточным для сглаживания возникающих в начальные моменты времени сверхсильных разрывов внутренней энергии.
Тест 4. В настоящей задаче рассматривается течение газа в виде двух расходящихся по газу ударных волн, между которыми располагается движущийся контактный разрыв.
Устойчивое решение этой задачи получается при параметрах расчета 0.3 ^ а ^ 0.8 и ¡3 ^ 0.01. Сходимость численного решения по сетке к автомодельным распределениям плотности и внутренней энергии показана на рис. 3. Полученное решение близко по точности к распределениям, приведенным в [8], но при его вычислении использовались большее значение коэффициента а = 0.7 и на порядок меньшее, чем в квазигазодинамическом алгоритме, число Куранта ¡3 = 0.01. Таким образом, в этом тестовом расчете разностная схема (5)—(10) оказывается менее точной
Р
35
30 25 20 15 10
-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 х
е
300
250
200
150
100
50
-0.3 -0.2 -0.1
0 0.1 0.2 0.3 0.4 0.5 х
Рис. 3. а — Тест 4. Распределение плотности р. Сходимость по сетке, б — Распределение внутренней энергии е. Сходимость по сетке. Штриховая линия: /г = 0.003; пунктирная: /г = 0.0015; сплошная: /г =0.00075; серая — автомодельное решение
и устоичивои, чем разностная схема, основанная на квазигазодинамических уравнениях.
Тест 5. Течение в этой задаче представляет собой неподвижный контактный разрыв. При отключении вязкости и теплопроводности (Se = 0) ширина контактного разрыва составляет один шаг сетки. Условие устойчивости этого алгоритма приведено на рис. 1 ,а. Рис. 4 демонстрирует распределение плотности в этой задаче для чисел Шмидта Se = 1, 0.1 и 0 для параметров расчета а = 0.5 и ¡3 = 0.5. Для Se = 0 там же приведена сходимость численного решения по сетке к автомодельному решению.
Р
1.5
h = 0.001, Se = 1 h = 0.001, Se = 0.1
1 ---h = 0.001, Se = 0
\ 1Л - А = 0.01, Se = 0
\ • 1.3
1.2
"I
, , , 1.0
-0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 0.10 *
Рис. 4. Тест 5. Распределение плотности р в неподвижном контактном разрыве. Зависимость решения от числа Шмидта Бс и шага пространственной сетки
Тест 6. Здесь течение представляет собой контактный разрыв, медленно движущийся по газу. Оптимальные параметры для устойчивого расчета такие же, как и в предыдущем тесте. В отличие от предыдущего случая значение Эс = 0 приводит к появлению осцилляций вблизи фронта волны. При 5с = 0.1, как и в задаче о неподвижном контактном разрыве, ширина разрыва составляет один шаг разностной сетки.
Результаты двух последних тестов близки по точности и устойчивости к данным, полученным на основе квазигазодинамического алгоритма [8]. Область устойчивости алгоритма для обоих тестов представлена на рис. 1 ,а.
Заключение
Проведенные расчеты показали, что система квазигидродинамических уравнений Шеретова, описывающая течения газа с уравнением состояния типа (4), позволяет сквозным образом, т. е. без выделения разрывов, рассчитывать нестационарные течения с контактными разрывами и ударными волнами небольшой интенсивности. При этом данный численный алгоритм обладает меньшей устойчивостью и точностью по сравнению с алгоритмами, основанными на квазигазодинамической системе уравнений, ограниченной использованием уравнения состояния идеального газа. Тем не менее рассмотренный алгоритм может иметь широкое применение при изучении трансзвуковых течений, встречающихся в технических приложениях, например течений
ТЕОРЕТИЧЕСКАЯ И МАТЕМАТИЧЕСКАЯ ФИЗИКА
33
в химических реакторах, численное моделирование которых, как правило, опирается на уравнение состояния неидеального газа.
Список литературы
1. Рождественский Б.Л., Яненко H.H. Системы квазилинейных уравнений. М., 1978.
2. Куликовский А.Г., Погорелое Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. М., 2001.
3. Liska R., Wendroff В. Comparison of several difference schemes on ID and 2D test problems for Euler equations. Technical report LA-UR-01-6225. LANL, Los Alamos, 2001.
4. Liska R., Wendroff В. // SIAM J. Sei. Comput. 2003. 25, N 3. P. 995.
5. Шеретов Ю.В. Математическое моделирование течений жидкости и газа на основе квазигидродинамических и квазигазодинамических уравнений. Тверь, 2000.
6. Елизарова Т.Г. Квазигазодинамические уравнения и методы расчета вязких течений. М., 2007.
7. Елизарова Т.Г., Соколова М.Е., Шеретов Ю.В. 11 Журн. вычисл. матем. и матем. физ. 2005. 45, № 3. С. 544.
8. Елизарова Т.Г., Шильников E.B. 11 Журн. вычисл. матем. и матем. физ. 2009. 49, № 3. С. 549.
Gas dynamic flow simulations based on quasi-hydrodynamic equation system T.G. Elizarova10, O.V. Bulatov2
1 Institute of Mathematical Modelling, Russian Academy of Sciences. Miusskaya sq. 4a, Moscow, 125047, Russia. 2Department of Mathematics, Faculty of Physics, M. V. Lomonosov Moscow State University, Moscow 119991, Russia.
E-mail: a telizar@mail.ru
Using sept well-known test problems reflecting the main features of unsteady inviscid gas flows, the possibilities of the quasihydrodynamic equation system for numerical simulation of the Euler equations are investigated.
Keywords: quasihydrodynamic algorithm, Euler equations, one-dimensional flows. PACS: 47.11.Df; 47.40.x. Received 10 June 2009.
English version: Moscow University Physics Bulletin 6(2009).
Сведения об авторах
1. Елизарова Татьяна Геннадьевна — докт. физ.-мат. наук, профессор, гл. науч. сотр.; e-mail: telizar@mail.ru.
2. Булатов Олег Витальевич — студент.
