Разработка неявной схемы для моделирования течений сжимаемого газа

Задачи газовой динамики постоянно возникают в самых разных областях науки и техники, что обуславливает потребность исследований в области механики газа. Наряду с теоретическим и экспериментальным подходами, в настоящее время, с учетом постоянного роста производительности компьютеров, вычислительная газодинамика становится все более эффективным инструментом исследований. Уравнения газодинамики нелинейны и многомерны – эти обстоятельства также говорят в пользу вычислительного подхода, поскольку теоретические исследования очень сложны. Вычислительная газовая динамика дает возможность проведения экспериментов, когда невозможно провести прямые измерения, и позволяет экономить средства, когда такие измерения слишком дороги. Интенсивное развитие и применение к расчету невязких сжимаемых течений идеального газа получили неявные разностные схемы, которые в отличие от явных схем снимаю жесткие ограничения на шаг интегрирования по времени. В код программы для моделирования задач гидродинамики cfd-2d [6] добавлена возможность расчета течения газа неявным конечно объемным методом. Приведен анализ решения типовых задач.

Основные уравнения. В основе математической модели теории движения газа и жидкости лежат законы сохранения массы, импульса и энергии. В декартовой системе координат (x, y) (в этой работе рассматривается только двумерные течения, при необходимости все рассуждения можно обобщить на трехмерный случай) течение сжимаемого идеального невязкого газа описывает уравнение Эйлера [3]:

д^'д^   ддЕ у

(1.1)

(1.2)

dt    дх    ду     ’ уравнение состояния представлено в форме:

V = £р(У - 1), где

U =

ри

/ Р\         p ри \         p ри \

Ри]          ри2+р              Рии

I ри 1, х( ) I рии Г у( ) I ри2 +р Г \ре/        \(Е+РрУи/        ХЕ+Рр^и/

(1.3)

Здесь t - время; р - плотность; p - давление; (u, v) - составляющие скорости в координатных направлениях (x, у) ; £ - внутренняя энергия единицы массы; Е = £ + ~~~ - полная энергия единицы массы; у = — - показатель адиабаты; С_р - теплоемкость при постоянном давлении;

c _v ...

C_v - теплоемкость при постоянном объеме.

Численный метод. В основу численного алгоритма положен метод конечных объемов. Рассматриваемый метод относится к семейству методов Годуновского типа и является консервативным. Методы Годуновского типа основаны на решении задачи Римана о распаде разрыва. Достаточно часто в области сложной геометрии не удается приемлемо построить не только структурированную сетку, но и многоблочную. Преимущество неструктурированных сеток перед регулярными заключается в большей гибкости при дискретизации физической области сложной формы [1]. В данной работе выбрана неструктурная треугольная сетка, удовлетворяющая принципу триангуляции Делоне. В таких сетках треугольники построены так, что в круг, описанный около любого треугольника, не попадает ни одного узла, отличного от вершин указанного треугольника [2, 7]. Треугольные ячейки сетки, построенной таким образом, выступают в роли непересекающихся контрольных объемов.

Система (1.1) записывается в виде интегральных законов сохранения по контрольному объему A j c границей д Aj , ориентация которой задается внешней единичной нормалью п:

^dSP dwEdS^. (2.1)

Ai dt JJAi v 7

Затем, применив формулу Остроградского-Гаусса, связывающую интеграл по объему с поверхностным интегралом, получим:

?UdS + ^.Fdl = 0,                               (2.2)

Д; dt         "dbi            , отсюда, е dUi _

S:~

чД^

Запишем (2.3) в полудискретном виде:

U^-U? ^T i

— 1ФЯЛ F • ndl.

S i

(2.3)

(2.4)

Обозначим через F/-+1 = (f^1, F—'1) поток через грань с номером j, контрольного объема Д[ на временном слое (n+1); Т[ - шаг по времени в Д^ контрольном объеме; St -площадь i-ого контрольного объема.

Контурный интеграл из соотношения (2.4) может быть аппроксимирован:

яг¹ = Z j—i hl [С¹ • Пх,+ ф» • nyj

(2.5)

где lti- длина j-ой грани треугольника Д^; п^ = (пх.., Пу.^ - единичная нормаль к j-ой грани i-ой ячейки.

Потоки через грань контрольного объема вычисляются на основе приближенного решения задачи Римана [3] методом Лакса-Фридрихса [5].

Таким образом, из (2.4) и (2.5) получаем:

п+1_пп

'   / S = —     1ц [f^1 • n-. + F - • Пу).         (2.6)

Рассмотрим разложение потоковой функции в ряд Тейлора:

С1 = F-^-'1 = ^ + (%)"   (U"' — u."),      (2.7)

^U — ^U i]

где   U- - осредненное по Роу [4] значение между U" и и"; и— = U$ - такое переобозначение, когда k-ая ячейка является соседней с i-ой ячейкой, и для них обоих j-ая граница i-ой ячейки является общей.

Матрица Якоби (^~) обладает набором действительных собственных значений и векторов, и поэтому может быть представлена в виде:

А = RAL = R(A+ + A~)L = А+ + А-, А+ = RA+L,

А ^- = RA - L,                            (2.8)

где A = diagfa], A+= 1(A + |A|), A- = 1(A — |A|), |A| = diag{\ki\}.

Опираясь на этот факт, для того чтобы сконструировать противопоточную разностную схему в дельта-форме, представим разложение (2.7) в виде:

С1 = Ъ + (Х) w+ - W ■ (Х)R+1 -Щ(2.9)

Таким образом, получим неявную разностную схему:

■■ ■ $ =1 ^{i [}Х⁺ (X ,) -"' ⁺(X» ^Л^ ⁺

- . i: Хл+(X)Л;  '+(X,) 4ЙН| - °.

Для стационарного случая вычисления могут проводиться с собственным шагом по времени для каждой ячейки, в качестве параметров используя условие Куранта-Фридрихса-

Леви. Тогда,

Т; = CFL • v3 г X.

(2.11)

Для нестационарного случая, Т ;

^‘ j=i []- iJ 'P(.^A iJ ^)]

где р(А) - спектральный радиус матрицы Якоби A- ( ^“ ).

для всех ячеек остаются постоянными.

Классическая реализация метода С. К. Годунова предполагает, что в каждой из ячеек значение газодинамических параметров U постоянно. Для двумерного случая с треугольной сеткой это можно представить как область, состоящую из треугольных ступенек, где высота ступеньки определяется значением U (понятно, что U это вектор, так что реально получается 4 области). Для повышения порядка точности по пространству применяются методы более интеллектуальной, чем кусочно-постоянная, реконструкции значений в ячейках.

В данной работе применяется, как кусочно-постоянная реконструкция, так и кусочнолинейная.

Кусочно-линейное распределение сеточной функции U ищется в виде:

U(x,y) = U; + (х-х^ + (у-у№,            (2.12)

где а. ; и Р ; - некоторые коэффициенты; (х ; , у ; )- координаты центра масс ячейки с номером i.

Формула (2.12) определяет, как могут быть найдены значения газодинамических параметров на границе i-ого треугольника. На настоящий момент не получены достаточные условия, обеспечивающие ограниченность всех вариаций численного решения на произвольной сетке, т. е. применимость каждого конкретного метода реконструкции должна быть проверена [3].

Роль коэффициентов а ; и Р ; играет знaчениеgrad и(х,у)в точке (х ; , у ; ).

Градиент вычисляется в общем случае из теоремы о градиенте:

VUdS^, Udl/VU*1^ Udi.             (2.13)

Л ;           ^ТЛ^    ,       s ^ГЛУ

Для решения линейной системы уравнений, возникшей из-за (2.10), в данной работе применялся метод Гаусса-Зейделя.

Численный эксперимент. Возможности реализованной неявной схемы в рамках пакета программ cfd-2d демонстрируются на примере решения задач обтекания профиля NACA0012 невязким сжимаемым газом и задачи о сверхзвуковом обтекании клина.

Обтекание профиля NACA0012. Симметричный профиль NACA0012 обтекается потоком невязкого газа с числом Маха M=0.7 под углом атаки 1.489 градусов, при давлении 46066.16 Па и температуре 248 К. Используя эти данные, вычисляются остальные параметры исходя из термодинамических соотношений.

В качестве сравнительной характеристики между результатом численного метода и ранее проведенным экспериментом рассматривается коэффициент давления по поверхности профиля.

Ср = -^2—,                                 (3.1.1)

(P w ^_M)/2

где р ^, и ^ - параметры набегающего потока.

Рис. 1. Распределение коэффициентов давления

Распределение коэффициентов давления показано на (рис. 1) и демонстрирует хорошее совпадение с экспериментом (синий маркер – эксперимент; красный – результаты расчета).

Сверхзвуковое обтекание клина. Рассчитывается задача о сверхзвуковом обтекании клина с углом в 10 градусов в плоском канале (рис. 2), при числе Маха M=2, при давлении P = 101325Па и температуре 300К. Производится анализ системы скачков уплотнения возникающих при обтекании клина и многократного отражения начального скачка от стенок канала.

Рис. 2. Сверхзвуковое обтекание клина с углом в 10 градусов в плоском канале

На (рис. 3) и (рис. 4) представлено распределение полей числа Маха и давления.

Рис. 3 Распределение полей числа Маха

Рис.4 Распределение полей давления