CC BY
40
УДК 519.63
НЕЯВНАЯ СХЕМА ПОВЫШЕННОЙ ТОЧНОСТИ С РАСЩЕПЛЕНИЕМ И ФАКТОРИЗАЦИЕЙ ДЛЯ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ НЕВЯЗКИХ ТЕЧЕНИЙ ГАЗА НА НЕСТРУКТУРИРОВАННЫХ СЕТКАХ
А.В. Русанов, Д.Ю. Косьянов
Институт проблем машиностроения им. А.Н. Подгорного НАН Украины, ул. Дм. Пожарского, 2/10, Харьков, 61046, Украина, e-mail: rusanov@ipmach.kharkov.ua, kosyanovdima@gmail.com
Аннотация. Рассмотрена новая неявная безытерационная конечно-объёмная схема с расщеплением и факторизацией для численного интегрирования системы уравнений Эйлера на неструктурированных сетках в задачах моделирования невязких течений сжимаемого газа. Представлены результаты численного моделирования обтекания аэродинамического профиля NACA0012 потоком сжимаемого газа.
Ключевые слова: неявная схема, расщепление, факторизация, метод контрольного объёма, невязкие течения газа, неструктурированная сетка.
1. Введение. Для численного решения начально-краевых задач газогидродинамики с помощью сеточных методов в областях сложной формы удобно использовать неструктурированные сетки. При моделировании стационарных физических полей с помощью методов установления по маршевой координате для ускорения сходимости численного решения целесообразно применять неявные схемы [2]. В работе [3] впервые для произвольной неструктурированной сетки разработана неявная безытерационная схема с расщеплением и факторизацией по координатным и характеристическим направлениям, использующая пространственные переключатели. В работах [4, 5] на её основе предложены методы повышенного порядка точности для численного моделирования невязких течений жидкости и газа. В настоящей работе представлена новая форма записи расщепления неявного оператора без применения переключателей, а также результаты численного моделирования невязкого трансзвукового течения газа возле профиля NACA0012.
2. Основные уравнения. Неявная схема с расщеплением и факторизацией. В
двумерном случае в декартовых координатах система уравнений Эйлера, описывающая течения невязкого сжимаемого совершенного газа при отсутствии массовых сил, может быть представлена в дивергентной форме:
Здесь - многосвязная расчётная область, декартовы координаты и время; Q,E,F
- вектор консервативных переменных и векторы потоков; р,п,ь,р - плотность, компоненты
dQ dE(Q) dF (Q)
0, (x,y) € ü, t> 0, Q = (p,pu,pv,ph)T , (1)
вектора скорости и статическое давление газа; 7, Н - показатель адиабаты и полная энергия единицы массы газа.
В момент I = 0 задано начальное условие <5 (ж, у, 0) = (р\ (ж, у), (ж, ¡/) е О = О и <90, а при £ > 0 - граничные условия, разновидность которых приведена в [2]. Дискретизация расчётной
ма
области выполняется с помощью согласованной неструктурированной сетки С= У Pi, где Pi -
i=1
ячейка сетки, N0 - количество ячеек сетки. Система уравнений (1) имеет гиперболический тип, гг г л дЕ г> дР
поэтому матрицы Якоои А=—— и В=—— векторов потоков можно привести к диагональному виду:
Ах = Р~1АРХ = <11щ(и—а,и,и,и+а), Ау = Ь^ВЬу = а,, у,у,у+а), а = \fypjp ■
Начально-краевая задача для (1) решается с применением метода контрольного объёма, в качестве которого выбирается ячейка сетки. После интегрирования по
Pí € 0,1 = и применения теоремы Грина система уравнений (1) принимает вид
= Кг{Я), = !я{х, у, Щ.хЛу , Щ<2) = -1 Епх + Рпу)с1Ь , (2)
где Li - граница i-й ячейки, состоящая из рёбер Гj, j=1,Nj.
На временной оси задаётся сетка узлов вида tn = At ■ n, где n = 0,1,... - номер временного слоя. Значение шага интегрирования вычисляется с помощью выражения Аt = max ц цд Ц);
i=TJJG у
в котором d = min min - минимальный линейный размер ячеек разностной сетки,
i=l,NGj=hNi
- длина ребра ячейки, а v - число Куранта, значение которого выбирается из условия удовлетворения устойчивости схемы.
В системе уравнений (2), записанной в момент времени t = tn+i, выполним аппроксимацию производной по времени трёхслойным разностным выражением вида
3 5Qi+l - 5Qi dQ
n+1
+ О (Ai2) , 5Q?+1 = Qf+1 -Qf, of ():il,)
2А <И
а также проведём линеаризацию векторов потоков:
Еп+1 » Еп + А|4=„ 5Яп+1, Рп+1 » Рп + В 1=п 5Яп+1,
5Яп+1 = Я{х, у, Ьп+1) - Я{х, у, Ьп) = Яп+1 - Яп.
После учёта описанных предположений, систему уравнений (2) с точностью до величины О (АЬ2) можно записать в виде
/1 О О /Ч
(Апх + Впу)\ь=ьп5дп+Ч,ь = + -АШ?, (Ц = - — . (3)
При моделировании стационарных решений методом установления значения 50>'1+1 и должны стремиться к 0 при увеличении номера временного слоя. В таком
случае порядок точности по пространству неявной схемы определяется порядком аппроксимации явного оператора Щ. Отсюда, интеграл в левой части (3) целесообразно вычислить приближенно с помощью формулы прямоугольников. Разложим матрицы Якоби на «положительную» и «отрицательную» части
Сф = С+ + С-, С± = ЬфЛ±£-1, Л± = 0.5Лф ± 0.5|ЛфОф=х=Л, Сф=у=В,
где |Л| - диагональная матрица, состоящая из модулей элементов матрицы Л. В результате, получим,
£(Лих + Впу) 1=гп5Яп+1йЬ £ |г;.|Нф,з( (с+ + с-) ¿дп+1), .
Г. 3 = 1 ф = х,у
Для определения значений С±£^п+1 на рёбрах ячейки воспользуемся информацией про области влияния дифференциального уравнения гиперболического типа. Слагаемые Сф+^п+1 вычислим приближенно по значениям вектора неизвестных в смежной ячейке, расположенной левее текущего ребра вдоль координатного направления ф, а слагаемые С—- по значениям в смежной ячейке, расположенной правее.
Для явной записи в вид схемы номера ячейки произвольной неструктурированной сетки предложено дополнительно раскладывать каждую из компонент вектора нормали ^-го ребра на «положительную» и «отрицательную» части
Пф,з = + п-3, п±3 = 0.5Пф,э ± 0.51Пф^ |, ф = х, у, что позволило записать следующие приближения (с первым порядком точности):
В представленных выражениях, уже в явном виде, содержатся номер г текущей и номер г3-смежной по ^-му ребру ячеек. С учётом всех промежуточных действий систему уравнений (3) можно записать в виде
/ + 1+ + 1" +В+ + = /", (4)
где I - единичный оператор; /п - правая часть;
(^п+1). = н Е ^КМ^1+* Е ^К^Ь/^Г1-
3=1 3=1
Неявный оператор (4) представлен в виде расщепления на локально-одномерные операторы, учитывающие координатные и характеристические направления. Т.к. =0(А1), то с ошибкой О (Д£2) его можно представить в факторизованном виде:
(I + Л+) (I + Л-) (I + В+) (I + В-) = /п
Разностная схема в дробных шагах имеет вид:
1 _ о
Шаг 0. /П = -^П + -АШП;
3 3
Шаг 1а. (/ + 1+) ^п+1/4 = /п;
Шаг 1Ь.
1 + А-)5Яп+2,А = 5Яп+1,А-
Шаг 2а. Шаг 2Ь.
Шаг 3.
/ + Б+Шп+3/4 = ^п+2/4;
1 + в-)бдп+1 = бдп+3/4;
Яп+1 = Яп + 5Яп+\
Для рёбер ячейки, лежащих на границе расчётной области, значение 6^ часто принимают равным 0. Способ эффективного решения систем уравнений для шагов 1а, 1Ь, 2а, 2Ь представлен в [3]. Для вычисления правой части на шаге 0 с высоким порядком точности по пространству используется кусочно-полиномиальная реконструкция функций с многочленом в форме Тейлора [4,5].
3. Численные результаты. Рассмотрена задача обтекания профиля МЛСЛ0012 (рис. 1 а) набегающим потоком невязкого сжимаемого газа с числом Маха Мж = 0.8 и углом аж = 1.25 [5]. Внешняя граница расчётной области удалена на расстояние 6-ти длин хорды профиля. Дискретизация выполнена с помощью «грубой» и «подробной» неструктурированных сеток, состоящих из 5106 (80 интервалов разбиения на каждой стороне профиля) и 20652 (160 интервалов разбиения) треугольных ячеек соответственно. Все результаты получены с применением кусочно-квадратичной реконструкции на дробном шаге 0. Число Куранта V равно 10.
а)
б)
)
в
Рис. 1. Разностная сетка и результаты расчёта: (а) - дискретизация расчётной области; (б) - изолинии коэффициента давления (подробная сетка); (в) -коэффициент давления на поверхности профиля.
На рис. 1 в прозрачными и закрашенными маркерами отмечены результаты на грубой и подробной сетках соответственно, а сплошной линией - результаты работы [5] (структурированная сетка с четырёхугольными ячейками; количество интервалов разбиения - 280). Полученные результаты хорошо соответствуют имеющимся представлениям о структуре трансзвукового течения со скачками уплотнения. Положение и интенсивность скачков одинаково хорошо моделируются как на нижней, так и на верхней поверхности профиля даже для грубой сетки. Нефизических осцилляций численного решения вблизи разрывов не наблюдается.
4. Выводы. В работе представлена новая неявная безытерационная конечно-объёмная схема повышенного порядка точности с расщеплением и факторизацией для численного интегрирования системы уравнений Эйлера на неструктурированных сетках при моделировании невязких течений сжимаемого газа. Результаты численного моделирования невязкого трансзвукового течения газа возле профиля КЛСЛ0012 показали высокую точность и устойчивость предложенного метода. Дальнейшие направления развития схемы связаны с задачами моделирования пространственных вязких нестационарных течений газа в областях сложной формы.
Литература
1. Русанов А.В., Ершов С.В. Математическое моделирование нестационарных газодинамических процессов в проточных частях турбомашин / Харьков: Ин-т пробл. машиностро-ен. НАН Украины, 2008. - 275 с.
2. Русанов А.В., Косьянов Д.Ю. Неявная схема для численного интегрирования уравнений гиперболического типа на неструктурированных сетках // Проблемы машиностроения. - 2010. - 13; 3. - С.30-37.
3. Русанов А.В., Косьянов Д.Ю., Косьянова А.И. Использование неявной схемы расщепления для численного моделирования трансзвуковых течений идеального газа // Вестник двигателестроения. - 2011. - 25; 2. - С.20-24.
4. Русанов А.В., Косьянов Д.Ю. Использование неявной схемы расщепления для моделирования течений невязкой несжимаемой жидкости // Восточно-Европейский журнал передовых технологий. - 2012. - 56; 2/4. - С.67-71.
5. Pulliam T.H., Barton J.T. Euler computations of AGARD working group 07. Airfoil test cases // AIAA Paper. - 1985. - 18. - P.1-34.
IMPLICIT HIGH-ORDER SCHEME WITH SPLITTING AND FACTORIZATION FOR NUMERICAL SIMULATION OF INVISCID IDEAL GAS FLOWS ON UNSTRUCTURED GRIDS
A.V. Rusanov, D.Yu. Kosyanov
A.N. Podgorny Institute for Mechanical Engineering problems NAS of Ukraine, Dm. Posharskoho St., 2/10, Kharkov, 61046, Ukraine, e-mail: rusanov@ipmach.kharkov.ua, kosyanovdima@gmail.com
Abstract. New implicit non-iterated finite-volume scheme with splitting and factorization for the numerical integration of the Euler equations on unstructured grids for modeling inviscid compressible gas flows is under consideration. Numerical results of compressible gas transonic flow over the airfoil NACA0012 simulation are presented.
Key words: implicit scheme, splitting, factorization, finite volume method, inviscid flows, unstructured grid.
