Scientific journal
Modern problems of science and education
ISSN 2070-7428
"Перечень" ВАК
ИФ РИНЦ = 1,006

THE SOFTWARE LOGOS. CALCULATION METHOD FOR VISCOUS COMPRESSIBLE GAS FLOWS ON A BLOCK-STRUCTURED MESHES

Veselova E.A. 1 Zhalnin R.V. 1 Deryugin Yu.N. 1 Zelenskiy D.K. 1 Kozelkov A.S. 1 Struchkov A.V. 1
1 Russian Federal Nuclear Center – VNIIEF
1508 KB
The paper presents the methodology for calculating three-dimensional viscous gas in a block- structured Euler grids. The method is applied for internal and external aerodynamics. The mathematical model described the laminar and turbulent flow of compressible flows bases on the Navier-Stokes equations and one- and two-parameter turbulence models. The original equations are written in Cartesian coordinates. Both explicit and implicit methods are constructed for the integration of time-dependent equations. The explicit time integration techniques base on the method of Runge-Kutta method. To construct the implicit finite difference techniques equations are written in a delta form. To approximate the derivatives describing convective transport in the implicit operator we use one-sided differences depending on the eigenvalues of the Jacobi matrix. The terms describing the viscous effects are approximated by central differences. To improve the approximation accuracy of terms describing convective transport reconstruction is used based on the linear and quadratic distribution of the parameters in cells of the difference grid satisfying monotony. For calculation of the flow on the faces the following algorithms are used: the solution of the Riemann’s problem, the Roe’s algorithm, the Harten’s method and the method of Chakravarti-Osher. The techniques implemented as part of the software LOGOS [3]. The structure of the block-structured grid data and the interaction of the database elements are described. Parallelization of numerical algorithms is given. Opportunities of methods are illustrated by calculating a number of test problems.
high order of accuracy method
block-structured meshes

Одним из важных приложений пакета программ ЛОГОС является решение задач аэродинамики. Для их решения были созданы вычислительные алгоритмы, основанные на неструктурированных сетках, методе конечного объема и явных неявных разностных схемах [3; 4] имеющих в общем случае первый порядок аппроксимации. Создание алгоритмов на неструктурированных сетках определялось тем, что сеточные генераторы для областей сложной конфигурации создавали сеточные модели, ячейками которых являются произвольные многогранники. В то же время существует широкий класс задач аэродинамики, где можно достаточно оперативно построить блочно-структурированные сеточные модели. На таких сетках развиты алгоритмы повышенного порядка аппроксимации. В данной работе приводятся вычислительные методики, созданные в рамках пакета программ ЛОГОС, применительно к расчету трехмерных вязких течений на блочно структурированных эйлеровых сетках.

Математическая модель. В консервативной форме уравнения Навье-Стокса, описывающие нестационарное трехмерное течение вязкого сжимаемого газа в декартовых координатах, записываются в виде:

. (1)

Здесь вектор консервативных переменных, вектор конвективных потоков, вектор вязких потоков, источниковый член, который учитывает неинерциальность системы отсчета, связанную с влиянием кориолисовой и центробежной сил.

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

Для описания турбулентных течений уравнения (1) выражаются через осредненные по Фавру параметры и дополняются уравнениями модели турбулентности, а вместо молекулярных коэффициентов переноса в уравнениях Навье-Стокса используются их эффективные значения.

В качестве начальных условий задаются распределения всех параметров потока в момент времени . На твердых стенках ставятся граничное условие непротекания для нормальной компоненты скорости и условие прилипание для касательной компоненты. Для температуры на твердой стенке могут быть выставлены различные граничные условия, в частности отсутствие теплового потока в стенку. На входной и выходной границах расчетной области могут быть поставлены различные граничные условия. Некоторые границы области могут быть объявлены свободными границами и границами периодичности потока.

Постановка задачи осуществляется в размерных переменных, а расчет проводится в безразмерных переменных. При обезразмеривании определяются характерные (критериальные) числа потока, параметры потока «нормализуются», изменяясь в обозначенных пределах, и это позволяет объединять невязки, возникающие в результате дискретизации исходных уравнений.

Сетка. Построение вычислительных алгоритмов основано на использовании блочно- структурированных сеток. Ячейки (контрольные объемы) такой сетки являются криволинейными шестигранниками, покрывающими расчетную область без зазоров и наложений. Гранями ячеек являются криволинейные поверхности, натянутые на четыре точки (узлы сетки). Отдельный блок блочно-структурированной сетки топологически эквивалентен кубу. Организация структуры данных для хранения информации о блочно- структурированной сетке сделана по аналогии организации структуры данных для представления неструктурированной сетки [3]. Это позволило практически полностью использовать все расчетные алгоритмы, созданные на неструктурированных сетках, изменив только некоторые алгоритмы, связанные с расчетом конвективных потоков на блочно-структурированных сетках.

Структура данных состоит из списков, описывающих ячейки, грани и узлы сетки. В списке ячеек сетки для каждой ячейки указывается номера граней , образующих ячейку, координаты геометрического центра ячейки, объем ячейки , а также список определяющих расчетных параметров, которые относятся на каждый момент времени к центру ячейки:

.

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

Все грани сортируются по типам. В списке сначала указываются внутренние грани. Затем указываются внешние грани, которые сортируются по типам граничных условий.

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

Расчетные методики. Для построения дискретной модели система уравнений (1) интегрируется по объему ячейки ограниченной поверхностью и, применяя теорему Гаусса-Остроградского, записывается в виде:

. (2)

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

, (3)

где вектор невязки, определяемый соотношением:

. (4)

Интегрирование разностных уравнений по времени осуществляется явными и неявными методами.

Явная дискретизация уравнения (3) основана на методах Рунге-Кутты 1-4-го порядка. В частности, используется следующая форма метода Рунге-Кутты 4-го порядка точности [11]:

(5)

При решении нестационарных задач явным методом интегрирование проводится с временным шагом, определяемым по оценке шага в каждом узле расчетной сетки:

, (6)

где число Куранта.

Шаг по времени в ячейке определяется из оценки спектрального радиуса якобианов дискретных операторов невязких и вязких потоков:

, (7)

где спектральный радиус матрицы , а спектральный радиус матрицы . Значения спектральных радиусов определяются по формулам:

, (8)

. (9)

При расчете стационарных течений интегрирование каждой ячейки проводится со своим шагом , определяемым по формуле (7).

Построение неявной разностной схемы основано на записи разностных уравнений в дельта-форме. Для этого представим неявную аппроксимацию в следующем виде:

, (10)

где:

, (11)

. (12)

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

, (13)

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

. (14)

С учетом (10)-(14) неявная разностная схема записывается в так называемой дельта-форме:

(15)

где параметр для схемы первого порядка точности имеет значение , для схемы второго порядка точности .

При решении нестационарных задач неявным методом используется процедура с введением шага по псевдовремени. Для этого в уравнение (2) вводится производная по псевдовремени

, (16)

где физическое время, псевдовремя. При производная по псевдовремени исчезает, а решение при этом соответствует физическому времени. Решение на новом временном слое определяется в процессе внутренних итераций по псевдовремени. Дискретизация уравнения (16) проводится аналогично аппроксимации уравнения (2) с тем лишь отличием, что производная по псевдовремени аппроксимируется с первым порядком, а для аппроксимации производной по физическому времени используется схема первого либо второго порядка точности:

(17)

где - приращение на итерациях, счетчик внутренних итераций по псевдовремени. Во внутренних итерациях и полагаются постоянными. Шаг по физическому времени выбирается исходя из требований точности, а шаг по псевдо-времени определяется условием Куранта-Фридрихса-Леви.

Расчет стационарных течений основан на методе установления, в котором уравнения (15) интегрируются в каждой ячейке со своим шагом по времени , который определяется по формуле (7) с коэффициентом запаса .

При использовании уравнений в дельта-форме граничные условия разрешаются явным методом, а для приращений искомых функций на границах области ставятся условия:

. (18)

Решение стационарных и нестационарных задач неявными методами находится с использованием расщепления по физическим процессам. Полная система уравнений расщепляется на систему уравнений Навье-Стокса, уравнения неразрывности компонент смеси и уравнения моделей турбулентности. Решение систем линейных алгебраических уравнений основано на использовании решателей библиотеки PMLP [1].

Точность получаемых решений во многом зависит от способа вычисления конвективных потоков. Использование блочно-структурированных сеток позволяет использовать реконструкции решения на гранях ячеек повышенного порядка точности и методы решения задачи Римана, учитывающие эти реконструкции [10; 13]. Опишем применяемые реконструкции для определения параметров на гранях ячеек. Пусть, как показано на рисунке 1, рассматривается сеточная линия , проходящая через центры ячеек и грань .

ξ, j 

Рисунок 1 – Шаблон, используемый для определения переменных слева (qL) и справа (qR) от грани j-1/2

С использованием введенных обозначений формулы реконструкции для определения значения параметров непосредственно слева и справа от рассматриваемой грани имеют вид, представленный в таблицах 1 и 2, где введено обозначение:

Таблица 1. Реконструкция переменных до пятого порядка без ограничителей

Параметры слева от грани (qL)

Параметры справа от грани (qR)

 

Таблица 2. Реконструкция переменных с ограничителями

Тип ограничителя

Параметры слева от грани (qL)

Параметры справа от грани (qR)

VAN ALBADA

VAN ALBADA LIGHT

MINMOD

ROE-DIFFUS

(KOLGAN)

В таблицах введены следующие обозначения:

, , . (19)

Параметр s может принимать значения s= -1 для реконструкции 2-го порядка и s=1\3 для реконструкции 3-го порядка. Слева от грани j-1/2 = 1+1/2 приграничной ячейки j=2 и справа от грани j+1/2 = N-1/2 приграничной ячейки j=N-1 вместо формул, приведенных в таблицах 1 и 2, используется следующая интерполяция 2-го порядка точности:

(20)

Реконструкция без ограничителей применяется для дозвуковых низкоскоростных течений. Интерполяции с ограничителями используется в задачах, где возможны сильные скачки уплотнения. Реконструкция может проводиться для вектора консервативных переменных , вектора примитивных переменных , вектора неконсервативных переменных , в том случае, когда используется предобуславливание уравнений при моделировании низкоскоростных течений и относительно инвариантов :

где матрица из левых собственных векторов матрицы Якоби .

Для определения конвективных потоков на грани в явном операторе используются следующие алгоритмы: решение задачи Римана (метод Годунова ) [2], схема Рое [12], метод Хартена [9], метод Чакраварти-Ошера [7] и гибридный метод на основе схемы Рое. В этом методе конвективные потоки определяются по следующей формуле:

. (21)

Здесь и весовые коэффициенты, а .

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

(22)

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

В неявной схеме аппроксимация конвективных потоков основана на вычислении односторонних потоков в зависимости от собственных чисел матрицы Якоби:

(23)

где матрицы и определяются как

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

Неявная аппроксимация диффузионного потока записывается в виде:

, (24)

где матрица имеет вид

. (25)

Программная реализация. Пакет программ ЛОГОС написан на языке программирования С++ с использованием библиотек межпроцессорных обменов MPI.

Пакет ЛОГОС ориентирован на использование неструктурированных пространственных сеток, состоящих из произвольных многогранников. Для представления расчётной сетки используется так называемая face-based модель памяти, в которой ведущую роль в топологии связей объектов сетки играют грани расчётных ячеек. Объект данных «Грань» содержит информацию об индексах вершин расчётной сетки, которыми она образована, и индексы двух ячеек, для которых грань является общей (рисунок 2).

Рисунок 2 – Структура сеточных данных

Блочно-структурированные расчётные сетки являются частным случаем неструктурированных сеток. Для реализации методов повышенного порядка точности на блочно-структурированных сетках используется модель памяти для неструктурированных сеток с введением в структуру данных двух дополнительных массивов для граней ячеек – информации о противоположных гранях.

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

В случае если сетка является блочно-структурированной, для каждой грани определяются индексы противоположных граней для ячейки С1 и ячейки С2. Используя эту информацию, можно получить индексы N ячеек в направлении ячеек С1 и С2.

Для проведения расчётов на многопроцессорных ЭВМ выполняется декомпозиция сеточной модели по процессорам. Для каждого фрагмента расчётной сетки формируется несколько слоёв обменных ячеек, как показано на рисунок 3.

Рисунок 3 – Организация межпроцессорных обменов

Результаты расчетов. Созданные методики были проверены на решении ряда задач, имеющих известное решение. Первой является задача расчета стационарного сверхзвукового обтекания клина с числом Маха На рисунке 4 представлена геометрия задачи. На рисунке 5 дано распределение полного давления, полученное по базовой схеме: метод Роу с линейной реконструкцией.

Рисунок 4 – Геометрия задачи

Рисунок 5 – Распределение полного давления, полученное по методу Роу с линейной реконструкцией

На рисунке 6 представлены зависимости полного давления от продольной координаты вдоль сечения y=0.85, z=0.5, полученные по методу Роу с различными реконструкциями.

Рисунок 6 – Зависимости полного давления от продольной координаты вдоль сечения y=0.85, z=0.5

На рисунке 7 представлено распределение давления, полученное по различным методам: метод Роу с линейной реконструкцией, метод Хартена, метод Ошера-Соломона, метод Чакраварти-Ошера. На рисунке 8 даны зависимости давления от продольной координаты вдоль сечения y=0.85, z=0.5, полученные по различным методам.

а) метод Роу с линейной реконструкцией

б) метод Хартена

в) метод Ошера-Соломона

г) метод Чакраварти-Ошера

Рисунок 7 – Распределение давления, полученное по различным методам.

Рисунок 8 – Зависимости давления от продольной координаты вдоль сечения y=0.85, z=0.5

Во второй задаче рассматривается дозвуковое обтекание профиля NACA0012 невязким потоком M=0.15 под углами атаки α=0° и α=5° [5]. На рисунке 9 приведено поле давления для двух углов атаки.

а) α=0°

б) α=5°

Рисунок 9 – Распределение давления для двух углов атаки

Расчеты проводятся по методу Роу с линейной реконструкцией и методу Роу с реконструкцией Ван Лира. Значения силы лобового сопротивления Fx и подъемной силы Fy для угла атаки α=0° приведены в таблице 3. Значения силы лобового сопротивления Fx для угла атаки α=5° приведены в таблице 4.

Таблица 3. Значения силы Fx и Fy для угла атаки α=0°

Теория

Fx = 0

Fy = 0

Линейная реконструкция

Fx = 0.0094

Fy = 0.0015

Реконструкция Ван Лира

Fx = 0.0015

Fy = 0.0015

Таблица 4. Значения силы Fx для угла атаки α=5°

Теория

Fx = 0

Линейная реконструкция

Fx = 0.036

Реконструкция Ван Лира

Fx = 0.0017

В третьей задаче рассматриваются обтекания профиля NACA0012 вязким потоком M=0.7 под углами атаки α=1.489°, α=3.0462°, α=3.9897°, и α=4.841°. Расчеты проводятся по методу Роу с линейной реконструкцией, с реконструкцией Колгана и с квадратичной реконструкцией. На рисунке 10 приведены поля давления для четырех углов атаки.

а) α=1.489°

б) α=3.0462°

в) α=3.9897°

г) α=4.841°

Рисунок 10 – Распределение давления для четырех углов атаки

На рисунке 11 представлена зависимость коэффициента подъемной силы от угла атаки, а зависимость коэффициента силы сопротивления от коэффициента подъемной силы показана на рисунке 12. Для сравнения на рисунках 10-11 приведены экспериментальные данные [12].

1 – эксперимент;

2 – метод Роу с линейной реконструкцией;

3 – метод Роу с реконструкцией Колгана;

4 – метод Роу с квадратичной реконструкцией.

Рисунок 11 – Зависимость коэффициента подъемной силы от угла атаки

1 – эксперимент;

2 – метод Роу с линейной реконструкцией;

3 – метод Роу с реконструкция Колгана;

4 – метод Роу с квадратичной реконструкцией.

Рисунок 12 – Зависимость коэффициента силы сопротивления от коэффициента подъемной силы

В четвертой задаче рассчитывается обтекание четырехзвенного профиля NASA с вязким сжимаемым газом с числом Маха M=0.201 под нулевым углом атаки [6]. Расчеты проводятся по методу Роу с различными реконструкциями. На рисунке 13 приведено распределение числа Маха. На рисунке 14 представлены значения коэффициента давления вдоль поверхности профиля для расчетов в сравнении с экспериментальными данными.

Рисунок 13 – Распределение числа Маха

Рисунок 14 – Значения коэффициента давления вдоль поверхности профиля

Пятая задача - задача о сверхзвуковом ламинарном течении в каверне. Постановка задачи взята из работы [14]. Рассчитываются колебания в трансзвуковой открытой ламинарной полости под действием распространяющихся ударных волн. Была выбрана конфигурация с отношением длины к высоте L / D = 1.5, числом Маха M∞ = 1,5 и ламинарным входящим граничным слоем толщиной δ = D / 5. В расчетную область входит 1,5D в направлении вдоль течения и 3D в перпендикулярном направлении. Начало оси находится на переднем краю полости. На входной границе скорость свободного течения U∞, число Маха M∞, статическая температура T∞ и статическое давление P∞ составляют соответственно 514,67 м/сек, 1,5, 293 К и 1,013 гПа. Профиль входящего граничного слоя задается на основании отношений Крокко-Буземана (Crocco-Busemann). Молекулярная вязкость удовлетворяет закону Сазерленда (Sutherland). В воздухе число Рейнольдса для свободного потока составляет Re = 1,73 x 104. Мы использовали равномерную сетку 600 х 240 (80 точек на длину D). Геометрия задачи показана на рисунке 15. Расчеты проводятся по схеме Годунова с линейной реконструкцией и реконструкцией пятого порядка аппроксимации.

Рисунок 15 – Геометрия задачи

Анализ результатов показывает, что в полости образуются регулярные вихри. Под действием вихрей, развивающихся в сдвиговом слое, образуется ударная волна на правой стенке каверны. Вдали от края эта ударная волна является стационарной, но рядом со сдвиговым слоем образуется маховская конфигурация, обеспечивающая вихревое движение. Под действием вихря маховская конфигурации выталкивается из каверны и распадается на две волны. Одна волна распространяется по течению, а вторая против течения в каверне. После того как волна отразится от левой стенки, начинается формирование новых вихрей. Приход вихря на нижний угол вызывает образование вихревых структур, перемещающихся внутрь полости. В результате чего возникают высокочастотные колебания. Для иллюстрации картины течения на рисунках 16 и 17 показаны поля давлений на два последовательных момента времени.

Рисунок 16 – Распределение давления на момент времени t=0.001

Рисунок 17 – Распределение давления на момент времени t=0.0028

В точке на дне каверны записывались пульсации давления. На рисунке 18 показаны пульсации давления, полученные в [14] по схемам седьмого и третьего порядка (слева) и по схемам второго и пятого порядка (справа), рассчитанные по пакету программ ЛОГОС. В работе [14] расчет начинался с установившегося решения, полученного по схеме первого порядка. В разработке по пакету программ ЛОГОС расчет проводился с начального постоянного распределения скорости. Поэтому в начальный момент времени имеются значительные различия в пульсациях давления.

Рисунок 18 – Пульсации давления в точке на дне каверны для схем седьмого и третьего порядка (слева) и схем второго и пятого порядка (справа)

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

Заключение. В рамках пакета программ ЛОГОС созданы модули расчета задач аэродинамики на блочно-структурированных сетках. В этих модулях реализованы одномерные реконструкции с повышенной разрешающей способностью, на основе которых определяются входные параметры для решения задачи Римана. Для расчета конвективных потоков реализованы ряд методов: метод Годунова, метод Рое, метод Хартена, метод Чакраварти-Ошера и гибридный метод на основе схемы Рое. Методики проверены на ряде задач, решенных по другим методикам. Показано, что схемы высокого порядка аппроксимации могут быть использованы для решения задач аэроакустики.

Рецензенты:

Щенников В.Н., д.ф.-м.н., профессор, заведующий кафедрой дифференциальных уравнений факультета математики и информационных технологий ФГБОУ ВПО «Мордовский государственный университет им. Н.П. Огарёва», г. Саранск.

Малыханов Ю.Б., д.ф.-м.н., профессор, профессор кафедры физики и методики обучения физике физико-математического факультета ФГБОУ ВПО «Мордовский государственный педагогический институт имени М.Е. Евсевьева», г. Саранск.