SIMULATION OF GAS FLOWS AT LOW MACH NUMBERS ON THE BASIS OF QUASIGASDYNAMIC SYSTEM OF EQUATIONS
Journal: VESTNIK OF ASTRAKHAN STATE TECHNICAL UNIVERSITY. SERIES: MANAGEMENT, COMPUTER SCIENCE AND INFORMATICS
Volume 2014 № 3
, 2014
Rubrics: MATHEMATICAL MODELING
Authors:
1.
Keldysh Institute of Applied Mathematics of Russian Academy of Sciences
2.
Keldysh Institute of Applied Mathematics of Russian Academy of Sciences
Type:
Article
Pages:
from 113
to 124
Status:
Published
Received:
05.07.2014
Accepted:
25.07.2014
Published:
25.07.2014
Subject area:
VAC 05.12.2013 Системы, сети и устройства телекоммуникаций
VAC 05.13.01 Системный анализ, управление и обработка информации (по отраслям)
VAC 05.13.06 Автоматизация и управление технологическими процессами и производствами (по отраслям)
VAC 05.13.10 Управление в социальных и экономических системах
VAC 05.13.18 Математическое моделирование, численные методы и комплексы программ
VAC 05.13.19 Методы и системы защиты информации, информационная безопасность
UDK 519.63
GRNTI 20.01 Общие вопросы информатики
GRNTI 28.01 Общие вопросы кибернетики
GRNTI 49.01 Общие вопросы связи
GRNTI 50.01 Общие вопросы автоматики и вычислительной техники
GRNTI 82.01 Общие вопросы организации и управления
VAC 05.13.01 Системный анализ, управление и обработка информации (по отраслям)
VAC 05.13.06 Автоматизация и управление технологическими процессами и производствами (по отраслям)
VAC 05.13.10 Управление в социальных и экономических системах
VAC 05.13.18 Математическое моделирование, численные методы и комплексы программ
VAC 05.13.19 Методы и системы защиты информации, информационная безопасность
UDK 519.63
GRNTI 20.01 Общие вопросы информатики
GRNTI 28.01 Общие вопросы кибернетики
GRNTI 49.01 Общие вопросы связи
GRNTI 50.01 Общие вопросы автоматики и вычислительной техники
GRNTI 82.01 Общие вопросы организации и управления
Language:
Russian
Keywords:
slightly-compressible gas flows, low Mach number, quasigasdynamic system of equations, decomposed pressure, pressure increments
Abstract (English):
A new approach to the simulation of slow slightly-compressible gas flows based on the system of equations distinct from the Navier - Stokes equations is under consideration. A peculiarity of substantially subsonic gas flows consists in the fact that on the one hand they are similar to liquid flows as the speed of small perturbations propagation is much higher than the gas velocity, but on the other hand, they show the property of medium compressibility. While Mach number tends to zero pressure-density coupling is violated. In the considered problems, gas is slightly-compressible, i. e. strong density drops are not observed even in the case, when the medium compressibility can not be neglected. At the same time, the pressure changes more perceptibly. Unlike the traditional approach to computation of gasdynamic flows in such cases, the pressure instead of the density should be considered as the primary dependent variable. In the paper, the procedure of pressure decomposition into two components - a volume-averaged part and a dynamic part - is proposed. To derive the dimensionless form of the equations, different reference values for these components are used that allows eliminating singularities at Mach number tending to zero. The notions of increments of volume-average and dynamic pressure parts as the difference between the current value of the corresponding component and its value on the previous time level are introduced. A special computational algorithm including the solution of elliptic equation for the increment of the dynamic component is developed. The test predictions indicate high accuracy and robustness of the proposed methods.
A new approach to the simulation of slow slightly-compressible gas flows based on the system of equations distinct from the Navier - Stokes equations is under consideration. A peculiarity of substantially subsonic gas flows consists in the fact that on the one hand they are similar to liquid flows as the speed of small perturbations propagation is much higher than the gas velocity, but on the other hand, they show the property of medium compressibility. While Mach number tends to zero pressure-density coupling is violated. In the considered problems, gas is slightly-compressible, i. e. strong density drops are not observed even in the case, when the medium compressibility can not be neglected. At the same time, the pressure changes more perceptibly. Unlike the traditional approach to computation of gasdynamic flows in such cases, the pressure instead of the density should be considered as the primary dependent variable. In the paper, the procedure of pressure decomposition into two components - a volume-averaged part and a dynamic part - is proposed. To derive the dimensionless form of the equations, different reference values for these components are used that allows eliminating singularities at Mach number tending to zero. The notions of increments of volume-average and dynamic pressure parts as the difference between the current value of the corresponding component and its value on the previous time level are introduced. A special computational algorithm including the solution of elliptic equation for the increment of the dynamic component is developed. The test predictions indicate high accuracy and robustness of the proposed methods.
Keywords:
slightly-compressible gas flows, low Mach number, quasigasdynamic system of equations, decomposed pressure, pressure increments
slightly-compressible gas flows, low Mach number, quasigasdynamic system of equations, decomposed pressure, pressure increments
Введение Предмет исследования - проблема численного моделирования течений вязкого сжимаемого газа при относительно небольших по сравнению со скоростью звука газодинамических скоростях, то есть при малых числах Маха - . Медленные течения газа являются основой технологических процессов в некоторых промышленных установках. Описание многих процессов, связанных с моделированием климата и переноса загрязнений воздушными массами, также опирается на исследование течений газа при небольших значениях скорости. Именно поэтому расчеты подобного рода течений представляют большой практический интерес. Особенность течений газа при малых числах Маха состоит в том, что они, с одной стороны, близки к течениям жидкости - в них скорость распространения малых возмущений намного больше скорости газа, а с другой - в них проявляется свойство сжимаемости среды. Численное моделирование таких течений сталкивается с многочисленными трудностями, которые возникают, когда сжимаемость значительно проявляется даже при существенно дозвуковой скорости течения. В качестве примеров такого рода задач в первую очередь можно назвать задачи свободной и смешанной конвекции при большой разнице максимальной и минимальной температуры. Математическая модель и вычислительный алгоритм должны учитывать особенности поведения давления при малых числах Маха, иначе возможно существенное снижение точности. При стремлении числа Маха к нулю нарушается взаимосвязь давления и плотности. В рассматриваемых задачах газ является слабосжимаемым, т. е., вообще говоря, не наблюдается сильных перепадов плотности даже в тех случаях, когда нельзя пренебречь сжимаемостью среды. При этом давление изменяется более ощутимо. Именно поэтому, в отличие от традиционного подхода к расчету газодинамических течений, давление, а не плотность должно рассматриваться в качестве основной зависимой переменной при разработке алгоритма. Реализация подобного подхода была осуществлена в [1] для полных уравнений Навье - Стокса. Ниже будет рассмотрена реализация этого подхода для другой модели - квазигазодинамической системы уравнений [2]. Квазигазодинамическая система уравнений с расщепленным давлением Квазигазодинамическая (КГД) система уравнений является дифференциальным приближением для кинетически согласованных разностных схем (КСРС) [2] в том случае, когда время пролета молекулой газа расчетной ячейки равно времени между столкновением частиц t. Кинетически согласованные разностные схемы были получены путем осреднения по скоростям молекул разностной аппроксимации кинетического уравнения Больцмана для одночастичной функции распределения. Расчеты на основе этой модели дают практически те же результаты, что и при использовании уравнений Навье - Стокса в области справедливости последних. В [2] было показано, что КГД-система совпадает с системой уравнений Навье - Стокса с точностью до где - газодинамическое время; - время прохождения звуком исследуемого объема. Выпишем КГД-систему для пространственно двумерного случая: (1) (2) (3) (4) Здесь r - плотность; u и v - горизонтальная и вертикальная компоненты скорости соответственно; P - давление; - полная энергия; e - внутренняя энергия. Параметр t = 2m/P, где m - коэффициент динамической вязкости. Правые части уравнений КГД-системы вносят дополнительную вязкость и дают новые вычислительные возможности по сравнению с уравнениями Навье - Стокса. Практика показывает, что при проведении расчетов членами со смешанными производными в правых частях уравнений (1)-(4) можно пренебречь. В случае дозвуковых течений, характеризующихся малыми числами Маха, решение этой системы сталкивается с теми же трудностями, что и традиционные подходы. В первую очередь сохраняется требование к шагу по времени. Это требование делает крайне неэффективным применение явных схем, которые являются предпочтительными, если в дальнейшем предполагается распараллеливание алгоритма. Однако, ввиду характерного поведения давления при малых числах Маха, даже применение неявных схем может не дать требуемой точности, а следовательно, скорость их сходимости может существенно замедлиться, что иногда на порядок увеличивает расчетное время. Физические особенности течений газов при малых числах Маха иллюстрирует рис. 1. Рис. 1. Особенности поведения давления при малых числах Маха Давление в каждой точке области колеблется около некоторой средней величины, постоянной для всей области. При характерное значение P0 среднего давления становится много больше характерного значения отклонений давления от среднего значения. Поле скорости, зависящее от градиента давления, определяется этими отклонениями. При переходе к безразмерным величинам используется следующее масштабирование: (5) и т. д. Здесь знаком штриха помечены безразмерные величины, а индексом 0 - характерные значения. При обычном способе обезразмеривания, не учитывающем описанной особенности поведения давления, общее давление обезразмеривается на указанное среднее значение P0: P = P0p', где P - размерное, а p' - безразмерное давление. Если представленное таким образом давление вместе со всеми остальными величинами (5) подставить в уравнения движения (2) и (3), то перед слагаемым grad p' появится множитель 1/M2. При проведении расчетов ошибки округления, возникающие при вычислении grad p', будут делиться на M2. При это приведет к значительным погрешностям вычислений. Для того чтобы правильно обезразмерить уравнения КГД-системы, представим исходное размерное давление в виде суммы двух составляющих - средней (осредненной по объему - ) и динамической (П) частей: . Среднее давление имеет термодинамический смысл и зависит только от времени. Динамическое давление зависит как от пространственных переменных, так и от времени. Эти слагаемые имеют разные характерные значения: обезразмеривается на среднюю величину , а - на , где и - характерные значения плотности и скорости соответственно. Можно показать, что такой подход к получению безразмерной формы уравнений позволяет избежать множителей типа 1/M2, следовательно, в уравнениях движения не появляется особенностей. Безразмерное давление имеет вид . (6) Здесь учитывается, что число Маха , где g - отношение удельных теплоемкостей. Так как была введена дополнительная переменная, для замыкания системы уравнений вводится интегральное соотношение, которое определяет : , (7) где VW - объем области течения. Если проинтегрировать (6) по области течения, то, учитывая определение (7), можно получить эквивалентное соотношение, которое более удобно при построении численных алгоритмов: В результате система КГД-уравнений с расщепленным давлением в безразмерном виде выглядит следующим образом (для краткости приведен пространственно одномерный случай): . (8) . (9) (10) Уравнения состояния: , (11) . (12) Уравнение энергии (10) записано относительно внутренней энергии. Система с расщеплённым давлением (8)-(10) отличается от исходной КГД-системы (1)-(4). Некоторые слагаемые в уравнениях (8)-(10) вместо общего давления содержат только его динамическую часть. Нет сомножителей типа 1/M2, следовательно, при стремлении числа Маха к нулю в уравнениях не появляется особенностей. Более того, несжимаемые жидкости, для которых M = 0, могут быть приняты к рассмотрению в качестве особого (предельного) случая в рамках данной модели, поскольку числа Маха не только нет в знаменателе, но при равенстве числа Маха нулю в правых частях уравнений сохраняются ненулевые слагаемые. Важно заметить, что первое уравнение состояния из (11) следует решать относительно плотности, которая определяется через давление, а не наоборот. При традиционном подходе к расчету газодинамических течений плотность определяется из уравнения неразрывности (8), а давление - через плотность из уравнения состояния. Однако для слабосжимаемых газов необходим другой порядок вычислений, поскольку, как уже отмечалось выше, в данном случае нарушается взаимозависимость давления и плотности. Вследствие этого для получения точного решения нужна специальная процедура для вычисления давления. При построении численного алгоритма будем опираться на уже разработанные методы расчета течений несжимаемой жидкости - методы коррекции давления (или предиктор-корректор). Поправки давления С целью определения давления введем поправки средней и динамической частей давления: , где n - номер слоя по времени. Учитывая эти поправки, для полного давления на текущем временном слое (n + 1) можно написать уравнение . (13) Для того чтобы выполнялось интегральное соотношение (12) для Pn и Pn+1 поправка динамической части давления должна удовлетворять тому же соотношению: (14) Если известны некоторые промежуточные значения поправок и такие, что для них верно (15) но не выполняется интегральное соотношение (14): , то с их помощью можно получить окончательные значения поправок. Для этого осредним по области течения равенство , (16) учтем (14) и получим выражение для поправки среднего давления: (17) Подставим (17) в (16) и найдем выражение для поправки динамического давления: (18) Таким образом, на каждом слое по времени сначала находим промежуточные значения поправок давления, потом окончательные значения поправок и уже затем - полное давление по формуле (13). Разностные уравнения для промежуточных значений поправок можно вывести путем алгебраических преобразований дискретных аналогов уравнений движения (9), неразрывности (8) и состояния (11). В двумерном случае уравнения движения в разностном виде для некоторого промежуточного временного слоя, который условно обозначим (n + 1/2), можно записать следующим образом (для конечных разностей использованы стандартные обозначения [3]): (19) (20) Аналогичные уравнения можно записать для (n + 1)-го слоя, где значения P в левой части берутся тоже с (n + 1)-го слоя: (21) (22) Вычтем (19) из (21), а (20) из (22) и получим выражения для (ru)n+1 и (rv) n + 1. Подставим эти выражения в разностную аппроксимацию уравнения неразрывности (8). Учтем, что на основании уравнения состояния (11) и выражения (15) (23) В результате получим уравнение для промежуточного значения поправки (dP)*: (24) Дифференциальное уравнение, аппроксимируемое данным дискретным уравнением, меняет свой тип при M = 0. При М ≠ 0 это уравнение гиперболического типа, а при M = 0 - эллиптического. В задаче с граничными условиями первого рода для скорости при M = 0 для (dP)* получим уравнение Пуассона с граничными условиями Неймана. Уравнение для промежуточного значения поправки среднего давления можно получить, осреднив по области течения разностную аппроксимацию уравнения неразрывности (8) с учетом дискретного уравнения состояния (23): (25) Алгоритм вычислений Для расчета медленных течений газа в рамках предложенной выше новой модели разработан следующий вычислительный алгоритм. На каждом слое по времени (n + 1 - номер текущего слоя по времени) вычисления происходят в два этапа: 1. Итерационный процесс для нахождения промежуточных значений искомых величин. 2. Процедура коррекции. В ходе итерационного процесса (k + 1 - номер текущей итерации) находим: 1) внутреннюю энергию e k + 1 из дискретного аналога уравнения энергии (10); 2) импульс (ru) k + 1 и (rv) k + 1 из разностных уравнений движения (19), (20); 3) температуру T k + 1 из дискретного аналога уравнения состояния для энергии (12); 4) промежуточные значения поправок средней и динамической частей давления и (dP)* из уравнений (25) и (24), используя вычисленные в ходе текущей итерации (этапы 2 и 3) значения импульса и температуры; 5) плотность r k + 1 из разностного уравнения состояния (23) через давление с предыдущего слоя по времени и вычисленные в ходе текущей итерации (этапы 3 и 4) значения температуры и промежуточных поправок давления. Для достижения заданной точности всех искомых величин требуется всего несколько итераций. После завершения итерационного процесса проводится коррекция, которая состоит в вычислении: a) окончательных значений поправок и dP согласно выражениям (17) и (18); б) компонент скорости на текущем слое через найденные в ходе итерационного процесса величины: в) всех компонент давления: Затем происходит переход на новый слой по времени. Задачи расчета течений вязкого сжимаемого газа при малых числах Маха требуют подробных пространственно-временных разностных сеток. Эти задачи связаны с обработкой на каждом временном слое больших объемов данных при довольно мелком шаге по времени и потому являются весьма трудоемкими. Успешному решению таких сложных задач способствует использование современных вычислительных систем сверхвысокой производительности. Предложенный алгоритм в целом легко адаптируется к многопроцессорной архитектуре, поскольку в основном используются явные численные методы. Некоторые трудности вызывает лишь параллельная реализация методов решения разностного эллиптического уравнения (24) для определения промежуточного значения поправки динамической части давления. Следует отметить, что основные вычислительные затраты связаны именно с этим уравнением. Однако использование современных экономичных методов решения систем линейных уравнений [4] позволяет сократить общее время счета и повысить эффективность расчетов. Примеры тестовых расчетов Классической тестовой для методов, используемых при моделировании динамики слабосжимаемого газа и несжимаемой жидкости, является задача о течении в каверне с подвижной верхней крышкой. Для нее в [5] при исследовании динамики несжимаемой жидкости были получены значения параметров течения на основе высокоточных расчетов системы уравнений Навье - Стокса. Результаты, полученные на основе КГД-модели с расщепленным давлением, сравнивались с эталоном из [5]. В начальный момент времени плотность и температура постоянны по всей области, компоненты скорости равны нулю, за исключением горизонтальной составляющей скорости на верхней границе («крышке»), которая равна 1,0 и постоянна во времени. Для всех остальных величин заданы условия второго рода: равенство нулю производной по нормали. В качестве расчетной области использовался единичный квадрат. Расчеты проводились на квадратных сетках при числе Маха M = 0,1 и числах Рейнольдса Re = 100 или 400. На рис. 2, а, б изображены линии тока установившегося течения для Re = 100 и Re = 400 соответственно, полученные с использованием вышеприведенного алгоритма. Линии тока имеют характерный для рассматриваемой задачи вид. а б Рис. 2. Линии уровня функции тока: а - Re = 100, M = 0,1, сетка 162 × 162; б - Re = 400, M = 0,1, сетка 162 × 162 На рис. 3 и 4 показаны профили скорости в центральных вертикальном и горизонтальном сечениях каверны для разных чисел Рейнольдса, расчетная сетка 162 × 162 узла. Сплошная линия, помеченная цифрой 2, соответствует результатам решения системы уравнений (8)-(12) с помощью рассмотренного выше алгоритма. Точками, помеченными цифрой 1, представлены данные из [5]. Наблюдается хорошее совпадение результатов вычислительных экспериментов с эталоном, что позволяет сделать вывод о возможности применения предложенного подхода для достаточно точного моделирования течений газа при малых числах Маха. Рис. 3. Профили компонент скорости в вертикальном и горизонтальном сечениях каверны: Re = 100, M = 0,1, сетка 162 × 162 Рис. 4. Профили компонент скорости в вертикальном и горизонтальном сечениях каверны: Re = 400, M = 0,1, сетка 162 × 162 Следующей тестовой задачей является задача о естественной конвекции в квадратной неоднородно нагретой каверне: вертикальные стенки нагреты до разной температуры, а горизонтальные изолированы от тепла. Обе компоненты скорости на границах равны нулю. В расчётах учитывается сила тяжести, поэтому следует явно учитывать гидростатическую компоненту давления. Иными словами, размерное давление состоит из трёх составляющих: В этом случае в правых частях уравнений сохранения импульса (9) и энергии (10) появляются соответствующие дополнительные члены. Эта задача является традиционным тестом для численных алгоритмов решения задач конвекции. В частности, результаты расчётов, представленные выше, показали хорошее совпадение с результатами высокоточных расчётов, описанных в [6] и полученных при использовании алгоритма, основанного на приближении Буссинеска. На рис. 5 представлены полученные нами параметры установившегося течения. а б в г Рис. 5. Параметры установившегося течения: а - линии тока; б - поле температуры; в - горизонтальная скорость; г - вертикальная скорость Расчеты производились при числе Рэлея и числе Прандтля . Решение прикладной задачи о течении в химическом реакторе В качестве прикладной задачи рассматривается течение газа в горизонтальном химическом реакторе в поле силы тяжести. Реактор представляет собой двумерный канал, в котором течёт смесь газов (метан, кислород, аргон), вступающих в химическую реакцию с платиновыми вставками на поверхности реактора (рис. 6). В рамках нашего исследования рассматриваются только газодинамические аспекты процесса. Вычислительная область - прямоугольник. Через левую границу втекает входной поток фиксированной температуры - около 300 К. На нижней и верхней границах для скорости задаются условия прилипания. Справа поток свободно вытекает из реактора. Температура в начальный момент везде постоянна, кроме двух симметричных сегментов на верхней и нижней границах, имеющих постоянную температуру 1200 К. Остальные части верхней и нижней границ являются термоизоляторами. Рис. 6. Постановка задачи о течении в химическом реакторе На рис. 7 и 8 показаны поле температуры и векторы скорости установившегося течения, полученные в результате расчётов. Рис. 7. Поле температуры Рис. 8. Векторы скорости Температура минимальна на входе и максимальна на нагретых сегментах. Благодаря влиянию силы тяжести в левом верхнем углу возникает зона рециркуляции (рис. 8), вследствие чего поле температуры становится несимметричным. Заключение Таким образом, нами подробно рассмотрен один из возможных подходов к моделированию течений слабосжимаемого газа при малых числах Маха (вплоть до граничного случая: М = 0) на основе квазигазодинамической системы уравнений с расщепленным давлением. Этот подход проиллюстрирован результатами тестовых расчетов, однако он может быть с успехом применен и для решения сложных задач, представляющих практический интерес, например задач охлаждения в турбомашиностроении, задач выращивания кристаллов в газоэпитаксиальных реакторах в микроэлектронной промышленности, задач горения топлива в различного рода камерах сгорания и других задач о течениях реагирующих смесей газов с сильным выделением тепла, а также для расчетов газодинамических течений в областях сложной формы, например в порах грунтов.
1. Churbanov A. G. A pressure-based algorithm to solve the full Navier - Stokes equations at low Mach number / A. G. Churbanov, A. N. Pavlov // Proceedings of the 4th European Computational Fluid Dynamics’98. 1998. Vol. 1, part 2. Chichester, United Kingdom, John Wiley and Sons, 1998. P. 894-899.
2. Chetverushkin B. N. Kineticheskie shemy i kvazigazodinamicheskaya sistema uravneniy / B. N. Chetverushkin. M.: MAKS Press, 2004. 332 s.
3. Samarskiy A. A. Teoriya raznostnyh shem / A. A. Samarskiy. M.: Nauka, 1989. 616 s.
4. Saad Y. Iterative methods for sparse linear systems / Y. Saad. SIAM, Philadelphia, 2003. 528 p.
5. Ghia V. High-Re solutions for incompressible flow using the Navier - Stokes equations and a multi-grid method' / V. Ghia, K. N. Ghia, C. T. Shin // J. Comp. Phys. 1982. Vol. 48. P. 387-411.
6. De Vahl Davis G. Natural convection of air in a square cavity: a bench mark numerical solution / G. de Vahl Davis // Int. J. for numerical methods in fluids. 1983. Vol. 3. P. 249-264.
