ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2015, том 55, № 6, с. 1058-1075
УДК 519.634
ПРЕДОБУСЛОВЛИВАНИЕ УРАВНЕНИЙ ГАЗОВОЙ ДИНАМИКИ В РАСЧЕТАХ ТЕЧЕНИЙ СЖИМАЕМОГО ГАЗА ПРИ МАЛЫХ ЧИСЛАХ МАХА
© 2015 г. К. Н. Волков, А. Г. Карпенко
(190005Санкт-Петербург, ул. 1-я Красноармейская, 1, Балтийский гос. техн. ун-т "Военмех")
e-mail: dsci@mail.ru Поступила в редакцию 24.07.2013 г.
Рассматриваются особенности моделирования низкоскоростных течений невязкого и вязкого сжимаемого газа и обсуждается конечно-объемная дискретизация уравнений газовой динамики при малых числах Маха на неструктурированных сетках. Для ускорения сходимости метода установления и повышения точности стационарного решения применяется метод предобусловливания, основанный на использовании физических переменных. Обсуждаются структура матрицы предобусловливания и диагонализация якобиана предобусловленной системы уравнений. Возможности подхода демонстрируются на примере решения ряда модельных задач газовой динамики в широком диапазоне чисел Маха. Библ. 34. Фиг. 9.
Ключевые слова: предобусловливание, газовая динамика, неструктурированная сетка, метод конечных объемов, уравнения Эйлера и Навье—Стокса.
DOI: 10.7868/S0044466915060150
1. ВВЕДЕНИЕ
Течения жидкости и газа описываются системой уравнений, включающей в себя уравнение неразрывности, уравнение изменения количества движения, уравнение изменения энергии и уравнение состояния (см. [1]). Уравнение изменения количества движения имеет различный вид в модели невязкой и вязкой среды (уравнения Эйлера для невязких течений и уравнения Навье— Стокса для вязких потоков). При необходимости уравнения Эйлера и Навье—Стокса дополняются уравнениями модели турбулентности, уравнениями химической кинетики и другими соотношениями.
Для интегрирования уравнений, описывающих течения несжимаемой жидкости, широкое применение находят метод искусственной сжимаемости (Artificial compressibility method, Pseudo-compressibility method, см. [2]), связанный с введением производной по времени от давления в уравнение неразрывности, метод проекции (Projection method, см. [3]), основанный на принципе расщепления по физическим процессам, а также методы, использующие процедуру коррекции давления (Pressure projection method, см. [4]). Общей чертой методов, основанных на процедуре коррекции давления, является формулировка разностной схемы относительно приращений искомых функций и решение уравнения Пуассона для поправки давления на каждом шаге во времени. При использовании неявных разностных схем широкое применение находят метод переменных направлений (Alternating direction implicit, ADI, см. [5]), методы релаксационного типа (см. [6]), симметричный метод Гаусса—Зейделя с LU-факторизацией (Lower-upper symmetric Gauss—Seidel, LU—SGS, см. [7]) и другие. Обзор и сравнение различных подходов дается в [8].
Численные методы решения уравнений сжимаемого газа, хорошо работающие при умеренных дозвуковых и сверхзвуковых скоростях потока, оказываются неэффективными или непригодными для расчета течений с малыми числами Маха (при М < 0.2, см. [1]), что проявляется в ухудшении сходимости процесса установления и падении точности получаемых стационарных решений (см. [9]—[12]). Замедление сходимости метода установления объясняется возрастающей при М —»- 0 жесткостью уравнений Эйлера и Навье—Стокса (особенность проявляется на дифференциальном уровне), описывающих течения сжимаемого газа, которая характеризуется отношением максимального и минимального собственных значений якобиана (отношение мак-
симальной и минимальной скорости распространения возмущений). Шаг интегрирования по времени определяется скоростью наиболее быстрой волны (акустические волны, X = \u + c\), а время достижения стационарного состояния зависит от скорости наиболее медленной волны (конвективные волны, X = \и\). При решении вязких задач и в расчетах турбулентных течений на растянутых сетках в пограничном слое шаг интегрирования по времени ограничивается акустическими модами решения и шагом сетки в направлении, нормальном к стенке (см. [13], [14]).
Численное моделирование течений при малых числах Маха проводится в рамках уравнений Эйлера или Навье—Стокса, записанных в приближении несжимаемой жидкости, с использованием методов, разработанных для этого класса уравнений. При М < 0.3 модель несжимаемой среды оказывается достаточно точным приближением, давая погрешность около 5%. Необходимость использования полных уравнений Эйлера или Навье—Стокса возникает при моделировании высокоскоростных течений с обширными подобластями низких скоростей (см. [15]) (например течений с зонами торможения потока и рециркуляционными зонами, внутренних течений в диффузорах с дозвуковой скоростью на входе) и низкоскоростных течений с изменениями плотности и температуры в результате подвода тепла (например свободно-конвективных течений).
Переход к предельной при М —► 0 форме уравнений Навье—Стокса для гипозвуковых неизотермических течений вязкого газа позволяет частично устранить трудности, возникающие при расчете этих течений в рамках полной системы уравнений Навье—Стокса (см. [16]).
Распространенным способом устранения вычислительных трудностей при М —»- 0 является применение различных методов предобусловливания (preconditioning) исходных уравнений, направленных на выравнивание порядков собственных значений якобиана при всех М < 1 (см. [17]—[23]). На дифференциальном уровне предобусловливание модифицирует члены с производной по времени в уравнениях изменения количества движения. При установлении решение модифицированной (предобусловленной) системы совпадает с решением исходной системы уравнений. Для нахождения нестационарного решения задачи применяется двойная процедура интегрирования по времени (dual time-stepping, см. [24]).
Предобусловливание находит также широкое применение для ускорения сходимости итерационных методов решения систем разностных уравнений, порожденных конечно-разностной или конечно-объемной дискретизацией уравнений Эйлера и Навье—Стокса (жесткость проявляется на матричном уровне).
Предобусловливание позволяет устранить жесткость исходной системы уравнений и ускорить сходимость метода установления (см. [9]—[12]), а также восстановить точность расчета дозвуковых течений путем использования модифицированной дискретизации конвективных потоков в предобусловленных уравнениях (см. [25], [26]). В общем случае предобусловливание приводит к изменению формы основных уравнений и свойств разностных схем из-за введения искусственной вязкости, а также ставит вопрос о применимости граничных условий. Точность предобусловленных разностных схем ухудшается при увеличении числа Маха. Теоретические вопросы, связанные с предобусловливанием уравнений Эйлера и Навье—Стокса при малых числах Маха, обсуждаются в [8], [20], [27], [28], а сравнение различных подходов приводится в [29].
На практике широкое распространение получили методы, развитые в [9]—[11] и [17]—[19]. Область применения метода из [9]—[11] ограничивается центрально-разностными схемами, которые хорошо работают при М < 1, но становятся диссипативными при моделировании сверхзвуковых течений. Метод из [17]—[19] может быть без особых трудностей применен к противо-поточным разностным схемам и получил широкое распространение при решении задач внешней газовой динамики (см. [13], [14]). Применяется также метод локального предобусловливания уравнений Эйлера и Навье—Стокса (см. [30], [31]), в котором переход на модифицированные уравнения осуществляется в зависимости от локального числа Маха (внешние течения) или локального поля давления (внутренние течения). Во многих случаях методы предобусловливания используются вместе с другими методами ускорения сходимости (см. [15]), такими как метод сглаживания невязки (Residual smoothing) и многосеточными методами (Multigrid method).
В то время как численные расчеты обычно производятся на основе уравнений, записанных в консервативных переменных, для построения матрицы предобусловливания применяются физические переменные, что позволяет упростить построение матрицы предобусловливания (см. [20]). B подходе из [20] в качестве зависимой переменной используется энтропия (симмет-ризованные переменные), а в подходах из [9]—[11], [16] — температура (физические переменные). В [9], [10] используется матрица предобусловливания, которая видоизменяет только урав-
нение энергии. Метод, реализованный в [11], предназначен для моделирования вязких течений. В [32] предобусловливание реализуется таким образом, чтобы оптимизировать скорости распространения волн во всем диапазоне чисел Маха (оптимальное число обусловленности).
В данной работе развивается метод предобусловливания уравнений Эйлера и Навье—Стокса, использование которого позволяет построить универсальную вычислительную процедуру расчета течений невязкого и вязкого сжимаемого газа в широком диапазоне чисел Маха (от существенно дозвуковых до транс- и сверхзвуковых скоростей потока). Для построения матрицы предобусловливания применяется подход, предложенный в [16] и реализованный в одномерном случае в [1]. Подход основывается на использовании физических переменных (в качестве одной из физических переменных принимается температура). Особенности реализации подхода состоят в форе записи потоков, расчете диссипативного члена при нахождении потоков через грани контрольного объема и представлении матриц при диагонализации якобиана невязких потоков предобусловленной системы уравнений. Дается компактная форма записи диссипативного члена в разностной схеме расчета потоков. Возможности подхода демонстрируются на примере решения ряда модельных задач внутренней газовой динамики.
2. РАСЧЕТЫ ПРИ МАЛЫХ ЧИСЛАХ МАХА Рассмотрим линеаризованные уравнения Эйлера
ди+А ди = 0,
д? дх
где А = д¥/ди — якобиан. Собственные значения якобиана таковы: = и, Х2 = и — с, Х3 = и + с. Число обусловленности матрицы А имеет вид
к(А) = е(А)е(А-1),
где р — спектральный радиус. При к > 1 матрица А является плохо обусловленной, что приводит к трудностям решения системы разностных уравн
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.