CC BY
37
САЛЬНИКОВ В. Д.
РАЗРАБОТКА НЕЯВНОЙ СХЕМЫ ДЛЯ МОДЕЛИРОВАНИЯ ТЕЧЕНИЙ СЖИМАЕМОГО ГАЗА
Аннотация. В статье описана неявная противопоточная конечно объемная схема в дельта-форме для моделирования течений идеального сжимаемого газа. В код программы для моделирования задач гидродинамики cfd-2d, разрабатываемой на кафедре прикладной математики, дифференциальных уравнений и теоретической механики, добавлена возможность расчета течения газа неявным конечно объемным методом. Приведен анализ решения типовых задач.
Ключевые слова: газовая динамика, численное моделирование, уравнение Эйлера.
SALNIKOV V. D.
CREATING IMPLICIT SCHEME FOR COMPRESSIBLE GAS FLOWS SIMULATION
Abstact. The article describes the implicit finite volume upwind scheme in delta-form for calculation of an ideal compressible gas flows. The ability to calculate an implicit finite volume upwind scheme in delta-form for calculation of an ideal compressible gas flows was added to the program for simulation of hydrodynamics cfd-2d, developed by the Chair of Applied Mathematics, Differential Equations and Theoretical Mechanics. The article includes an analysis of typical problems.
Keywords: gas dynamics, numerical simulation, Euler equation.
Задачи газовой динамики постоянно возникают в самых разных областях науки и техники, что обуславливает потребность исследований в области механики газа. Наряду с теоретическим и экспериментальным подходами, в настоящее время, с учетом постоянного роста производительности компьютеров, вычислительная газодинамика становится все более эффективным инструментом исследований. Уравнения газодинамики нелинейны и многомерны - эти обстоятельства также говорят в пользу вычислительного подхода, поскольку теоретические исследования очень сложны. Вычислительная газовая динамика дает возможность проведения экспериментов, когда невозможно провести прямые измерения, и позволяет экономить средства, когда такие измерения слишком дороги. Интенсивное развитие и применение к расчету невязких сжимаемых течений идеального газа получили неявные разностные схемы, которые в отличие от явных схем снимаю жесткие ограничения на шаг интегрирования по времени. В код программы для моделирования задач
гидродинамики [6] добавлена возможность расчета течения газа неявным конечно
объемным методом. Приведен анализ решения типовых задач.
Основные уравнения. В основе математической модели теории движения газа и жидкости лежат законы сохранения массы, импульса и энергии. В декартовой системе координат (х, у) (в этой работе рассматривается только двумерные течения, при необходимости все рассуждения можно обобщить на трехмерный случай) течение сжимаемого идеального невязкого газа описывает уравнение Эйлера [3]:
(1.1)
"у
уравнение состояния представлено в форме:
^ + ^ + ^ = о
д1 дх ду '
Р = £Р(У - 1)
(12)
где
и =
/Р
ри рр
\рЕ,
/
рх(и) =
ри ри2 + р
/
риу \(рЕ + р)и/
Ру(и) =
ру
рир рр2 + р
\(рЕ + р)у)
(13)
Здесь время; р - плотность; р - давление; (и, \) - составляющие скорости в координатных
и2+У2
--полная энергия
направлениях (х, у); г - внутренняя энергия единицы массы; Е = е + ■
единицы массы; у = — - показатель адиабаты; Ср - теплоемкость при постоянном давлении;
Су
Су - теплоемкость при постоянном объеме.
Численный метод. В основу численного алгоритма положен метод конечных объемов. Рассматриваемый метод относится к семейству методов Годуновского типа и является консервативным. Методы Годуновского типа основаны на решении задачи Римана о распаде разрыва. Достаточно часто в области сложной геометрии не удается приемлемо построить не только структурированную сетку, но и многоблочную. Преимущество неструктурированных сеток перед регулярными заключается в большей гибкости при дискретизации физической области сложной формы [1]. В данной работе выбрана неструктурная треугольная сетка, удовлетворяющая принципу триангуляции Делоне. В таких сетках треугольники построены так, что в круг, описанный около любого треугольника, не попадает ни одного узла, отличного от вершин указанного треугольника [2, 7]. Треугольные ячейки сетки, построенной таким образом, выступают в роли непересекающихся контрольных объемов.
Система (1.1) записывается в виде интегральных законов сохранения по контрольному объему с границей дА^, ориентация которой задается внешней единичной нормалью п:
ГГ ГГ (НУР йБ = 0.
(2.1)
Затем, применив формулу Остроградского-Гаусса, связывающую интеграл по объему
с поверхностным интегралом, получим:
I ТТ & , _
¡¡А.%(К + &АМ = 0, (2.2)
отсюда,
^^На/^ (23)
Запишем (2.3) в полудискретном виде:
^-= -±ф пй1. (2.4)
Обозначим через Р]ц+1 = (рху^1, Ру^1) поток через грань с номером у', контрольного
объема А I на временном слое (п+1);т1- шаг по времени в А I контрольном объеме; -площадь г-ого контрольного объема.
Контурный интеграл из соотношения (2.4) может быть аппроксимирован:
= 13=11ц [г?*1 • пх.. + Р^1 • пу.], (2.5)
где - длинау'-ой грани треугольника А[; щу = (щх^.,Пу.^ - единичная нормаль к j-ой грани 1-ой ячейки.
Потоки через грань контрольного объема вычисляются на основе приближенного решения задачи Римана [3] методом Лакса-Фридрихса [5]. Таким образом, из (2.4) и (2.5) получаем:
ип+1_,,п г -|
= -13=11ц ^ • Пхц + Р^1 • Пу.]. (2.6)
Рассмотрим разложение потоковой функции в ряд Тейлора:
кГ = Рх(и%+1) = Рх1 + (и^1 - (2.7)
где и™ - осредненное по Роу [4] значение между и™ и и ¡у; иу = иЦ - такое переобозначение, когда к-ая ячейка является соседней с /-ой ячейкой, и для них обоих у'-ая граница ьой ячейки является общей.
Матрица Якоби обладает набором действительных собственных значений и
векторов, и поэтому может быть представлена в виде:
А = ЯЛЬ = И(Л+ + Л-)Ь = А+ + А-, А+ = ЕЛ+Ь,
А- = ЕЛ-Ь, (2.8)
1 1 где Л = <Иад{\}, Л+ = ~(Л + |Л|), Л-= -(Л- |Л|), \Л\ = <ИадШ}.
Опираясь на этот факт, для того чтобы сконструировать противопоточную разностную схему в дельта-форме, представим разложение (2.7) в виде:
Р^1 = + №) (Оп+1 — Оп) + ) (0Г)+1 — Оп). хч хч \ди;(П,1Лу 1 ^ \ди/(ПЛ)^ I] ч;
(пЦ)
(2.9)
Таким образом, получим неявную разностную схему:
Аи,
1 $ + Е/=1 к]
п
№'■ + №)+ АО^ + ^У
Ч\ ХЧ Vди '(п,1]) 1 \dUJrnii
(п, I ])
АО?
п+1
+13=11
п
Уц
+ №) АО¡1+1 + №) АО?
УЧ \ди)(пП) 1 \ди)(пцл 1]
п+1
+
= 0.
(2.10)
*ди'(пЦ) ^ ди / (п,1 ])
Для стационарного случая вычисления могут проводиться с собственным шагом по времени для каждой ячейки, в качестве параметров используя условие Куранта-Фридрихса-Леви. Тогда,
т- = ГР1 ■ -
1 = Г Я=1[кГР(Аф]'
(2.11)
(дР\
где p(A) - спектральный радиус матрицы Якоби A= Для нестационарного случая, т£
для всех ячеек остаются постоянными.
Классическая реализация метода С. К. Годунова предполагает, что в каждой из ячеек значение газодинамических параметров и постоянно. Для двумерного случая с треугольной сеткой это можно представить как область, состоящую из треугольных ступенек, где высота ступеньки определяется значением и (понятно, что и это вектор, так что реально получается 4 области). Для повышения порядка точности по пространству применяются методы более интеллектуальной, чем кусочно-постоянная, реконструкции значений в ячейках.
В данной работе применяется, как кусочно-постоянная реконструкция, так и кусочно-линейная.
Кусочно-линейное распределение сеточной функции и ищется в виде:
О(х,у) = О[ + (х — х^щ + (у-
(2.12)
где а.1 и - некоторые коэффициенты; (х^, у{)- координаты центра масс ячейки с номером ь Формула (2.12) определяет, как могут быть найдены значения газодинамических параметров на границе ього треугольника. На настоящий момент не получены достаточные условия, обеспечивающие ограниченность всех вариаций численного решения на произвольной сетке, т. е. применимость каждого конкретного метода реконструкции должна быть проверена [3].
Роль коэффициентов а^ и ^ играетзначениедгай О(х,у)в точке (х^ у{). Градиент вычисляется в общем случае из теоремы о градиенте:
Ма^ О dS = фдAUdl,VU~1JдAUdL
дА
(2.13)
Для решения линейной системы уравнений, возникшей из-за (2.10), в данной работе применялся метод Гаусса-Зейделя.
Численный эксперимент. Возможности реализованной неявной схемы в рамках пакета программ еГё-2ё демонстрируются на примере решения задач обтекания профиля NACA0012 невязким сжимаемым газом и задачи о сверхзвуковом обтекании клина.
Обтекание профиля ЫАСА0012. Симметричный профиль КАСА0012 обтекается потоком невязкого газа с числом Маха М=0.7 под углом атаки 1.489 градусов, при давлении 46066.16 Па и температуре 248 К. Используя эти данные, вычисляются остальные параметры исходя из термодинамических соотношений.
В качестве сравнительной характеристики между результатом численного метода и ранее проведенным экспериментом рассматривается коэффициент давления по поверхности профиля.
Ср=-—(3.1.1) где рт, и22- параметры набегающего потока.
Рис. 1 . Распределение коэффициентов давления Распределение коэффициентов давления показано на (рис. 1) и демонстрирует хорошее совпадение с экспериментом (синий маркер - эксперимент; красный - результаты расчета).
Сверхзвуковое обтекание клина. Рассчитывается задача о сверхзвуковом обтекании клина с углом в 10 градусов в плоском канале (рис. 2), при числе Маха М=2, при давлении Р = 101325Па и температуре 300К. Производится анализ системы скачков уплотнения возникающих при обтекании клина и многократного отражения начального скачка от стенок канала.
Рис. 2. Сверхзвуковое обтекание клина с углом в 10 градусов в плоском канале На (рис. 3) и (рис. 4) представлено распределение полей числа Маха и давления.
Рис. 3 Распределение полей числа Маха
Рис.4 Распределение полей давления 6
ЛИТЕРАТУРА
1. Волков К. Н., Емельянов В. Н. Течения и теплообмен в каналах и вращающихся плоскостях. - М.: Физматлит, 2010. - 486 с.
2. Елизарова Т. Г. Лекции. Математическая модель и численные методы в динамике жидкости и газа. Подходы, основанные на системах квазигазодинамиеских и квазигидродинамических уравнений. - М.: Физ. факультет МГУ, 2005. - 224 с.
3. А. Г. Куликовский, Н. В. Погорелов, А. Ю. Семёнов. Математические вопросы численного решения гиперболических систем уравнений. - М.: Физматлит, 2001. - 608 с.
4. Roe P. L. Approximate Rieman Solvers, Parameters Vectors and Difference Schemes // Journal of Computations Phys., 1981. - Vol. 43. - pp. 357-378.
5. Toro E. F. Riemann solvers and numerical methods for fluid dynamics. - Verlag: Springer, 1999. - 624 p.
6. Пакетпрограмм для моделирования турбулентных течений сжимаемого газа. -[Электронныйресурс]. - Режим доступа: http://code.google.com/pZcfd-2d/
7. Программа генерации сетки. - [Электронный ресурс]. - Режим доступа: http://www.cs.cmu.edu/~quake/triangle.html
