2. Математическая модель и численный метод

2.1. Математическая модель

2.1.1 Основные уравнения и преобразование координат

С целью применения метода искусственной сжимаемости и использования итерационного процесса на каждом шаге по физическому времени осредненные по Рейнольдсу уравнения Навье-Стокса для трехмерных течений сжимаемого газа могут быть записаны в тензорной форме следующим образом [12], [13], [21]:

, (2.1.1)

, , (2.1.2)

где r - плотность среды,

xk - декартовы координаты,

uk - компоненты скорости в декартовой системе координат,

p - статическое давление,

t - физическое время,

t - псевдо-время,

- декартовы компоненты тензора эффективных напряжений,

- вязкие напряжения,

- напряжения Рейнольдса;

b - определяемый пользователем параметр искусственной сжимаемости.

Целью итерационного процесса является получение нулевого значения производной по псевдо-времени на каждом шаге по физическому времени. При расчете стационарного течения производная по физическому времени приравнивается к нулю. Система уравнений (2.1.1), (2.1.2) дополняется уравнением энергии в форме баланса энтальпии.

Для определения тензора вязких напряжений применяется модель ньютоновской жидкости

. (2.1.3)

где m mol - молекулярная вязкость.

Уравнения для характеристик турбулентности могут быть выражены в обобщенной форме закона сохранения скалярной величины

, (2.1.4)

где F - некоторый параметр турбулентности,

, - компоненты диффузионного потока F,

- компоненты турбулентного потока F,

- объемный источник F.

Полный диффузионный поток скаляра F описывается моделью градиентного типа:

, (2.1.5)

где - эффективный коэффициент диффузии скалярной величины.

Для областей сложной геометрии используемые уравнения преобразуются с помощью обобщенных криволинейных координат [12]. Независимые переменные, преобразующие декартовы координаты в обобщенные криволинейные, вводятся следующим образом:

.

Якобиан преобразования определяется как

,

Метрические коэффициенты могут быть вычислены, например, следующим образом

,

где якобиан обратного преобразования представляется в виде

Полезно отметить, что вектор имеет следующие декартовы компоненты

является нормалью к поверхности . Кроме того, если вычислительная область разделена на однородные ячейки с размерами , длина этого вектора равна площади поверхности соответствующей грани ячейки в физическом пространстве.

Применение преобразования координат к уравнению неразрывности приводит к следующему:

(2.1.6)

где контрвариантные компоненты потока массы, Gm , определяются как

Результатом преобразования уравнения сохранения количества движения является

. (2.1.7)

Уравнение переноса скаляра после преобразования примет вид

. (2.1.8)

Систему уравнений неразрывности, количества движения и переноса скаляра удобно представить в общей форме

где - вектор переменных,

- конвективный поток,

- диффузионный поток, и

- источниковый член.

Консервативная форма используемых уравнений позволяет использовать интегральную формулировку уравнений сохранения в пределах каждого конечного объема внутри расчетной области.

2.1.2 Моделирование турбулентности

Используемые уравнения. Принята концепция изотропной турбулентной вязкости, основанная на приближении Буссинеска. Напряжения Рейнольдса считаются пропорциональными осредненным по времени скоростям деформаций. Введением коэффициента турбулентной вязкости , данная концепция выражается как

. (2.1.9)

Коэффициент эффективной вязкости вычисляется как сумма молекулярного и турбулентного коэффициентов

В данной работе для определения коэффициента турбулентной вязкости используется двухпараметрическая высокорейнольдсовая k-e модель турбулентности со следующими уравнениями и замыкающими соотношениями [17]

(2.1.10)

(2.1.11)

где - квадрат осредненного тензора напряжений. Константами модели являются

. (2.1.12)

Как показали Като и Лаундер [15], k-e модель в стандартной формулировке переоценивает величину генерации кинетической энергии турбулентности в областях с малыми скоростями сдвиговых деформаций. Для устранения этого эффекта Като и Лаундер переопределили генерационный член через произведение вращательного и деформационного сомножителей

, (2.1.13)

где - модуль осредненной завихренности.

Пристенные функции. Используемые уравнения высокорейнольдсовой k-e модели применяются в комбинации с методикой пристенных функций [17]. Согласно методу пристенных функций рекомендуется располагать ближайшую к стенке точку сетки на расстоянии, удовлетворяющем следующему условию: , где - координата стенки. В случае сложного пространственного течения удовлетворить этому условию во всей расчетной области трудно. Для устранения этой проблемы используются модифицированные пристенные функции, предложенные в работе [16].

Согласно [16] используемое в пристенных функциях распределение скорости вблизи стенки описывается следующей трехслойной аппроксимацией

(2.1.14)

(2.1.15)

(2.1.16)

где определяется равенством значений, заданных (2.1.15) и (2.1.16).

Стандартные пристенные функции для кинетической энергии турбулентности и скорости ее диссипации определяются через

(2.1.17)

где - динамическая скорость. Модификация стандартных пристенных функций состоит в добавлении демпфирующего коэффициента, D, в (2.1.17):

(2.1.18)

где n=0.35 - эмпирическая константа. Подстановка (2.1.18) в соотношение для турбулентной вязкости

(2.1.19)

показывает, что коэффициент D можно рассматривать как поправочный коэффициент. Для определения D используются следующие выражения

(2.1.20)

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

2.2. Основы численной схемы и характеристика программного комплекса SINF

2.2.1 Пространственная дискретизация

Уравнения движения, неразрывности и переноса турбулентных характеристик дискретизируются в пространстве по методу контрольных объемов с определением искомых величин в центрах контрольных объемов (ячеек). Разностные соотношения, описывающие законы сохранения, получаются интегрированием дифференциальных уравнений по контрольным объемам с использованием теоремы Гаусса. Для кубической вычислительной ячейки с размерами интегрирование уравнений, записанных в обобщенных координатах, приводит к получению таких же выражений, как в методе конечных разностей.

Вычисление вклада конвективных слагаемых. Рассмотрим дискретизацию конвективной части невязки. Дискретизация второго порядка имеет вид:

(2.2.1)

Индексы i, j, k отвечают независимым переменным x 1, x 2, x 3. Дробный индекс обозначает, что величина вычисляется на грани ячейки. Будем рассматривать только дискретизацию потока , остальные потоки, и , вычисляются аналогично.

К потоку через грань ячейки добавляется слагаемое, связанное с численной диссипацией

. (2.2.2)

Первое слагаемое в (2.2.2) получено путем линейной интерполяции потоков в соседних ячейках

, , (2.2.3)

где - расстояние между центрами граней и . Для определения численной диссипации используется QUICK схема второго порядка Леонарда [18] (см. также [12], [13]); для уравнений переноса скаляра в модели турбулентности используется верхнепоточная схема первого порядка.

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

, (2.2.4)

Применяя изотропную модель эффективной вязкости, выражение для вычисления напряжения можно записать в форме

. (2.2.5)

Производные компонент скорости в направлении нормали к поверхности грани и в двух касательных направлениях вычисляется с использованием десятиточечного шаблона [13]. Диффузионные потоки скаляров вычисляются таким же способом. Поток через грань определяется как

, (2.2.6)

где производные снова вычисляются по десятиточечному шаблону .

2.2.2. Численная реализация

Вышеописанная численная схема реализована в программном комплексе SINF. Разработка кода начата в 1992 году [23]. Вычисления базируются на пространственной дискретизации второго порядка точности с заданием величин в центрах ячеек структурированных сеток, вписанных в произвольные границы расчетной области. Первая версия кода ограничивалась использованием одноблочных сеток. В 1997-1998 разработана и протестирована модернизированная многоблочная стационарно/нестационарная версия программного комплекса SINF. В 1998-1999 возможности программы расширены за счет моделирования трансзвуковых течений.

Трехслойная схема второго порядка применяется для продвижения по физическому времени в совокупности с эффективными неявными методами, реализующими итерационный процесс эволюции по псевдо-времени. Для вычисления конвективных потоков возможно использование верхнепоточных схем высокого порядка. Для моделирования турбулентности применимы как высоко-, так и низкорейнольдсовые модели турбулентности. В высокорейнольдсовых моделях используются модифицированные пристенные функции, которые обеспечивают слабую чувствительность результатов к положению первого узла сетки вблизи стенки. Набор вариантов граничных условий, задаваемых на сегментированной поверхности, позволяет решать широкий диапазон задач фундаментальной и промышленной аэродинамики. Некоторые подробности и результаты применения описанного программного обеспечения перечислены в списке ссылок [23], [20], [2], [7], [16].