О применении WENO-схем к моделированию реагирующих газовых потоков

Аннотация

В работе исследованы WENO-схемы 5 и 7 порядка с оптимальными весовыми коэффициентами без учета индикаторов гладкости решения и с индикаторами гладкости. Показано, что при их применении к решению задач многокомпонентной газовой динамики с химическими реакциями, диффузией, вязкостью и теплопроводностью WENO- схема 5 порядка с оптимальными весовыми коэффициентами дает более точный результат при меньшей трудоемкости вычислений.

Полный текст

Введение. В настоящее время широкое распространение получило исследование газодинамических течений в задачах химической промышленности, поскольку технологам для получения целевых продуктов необходимо знать большое количество параметров проведения реакции. Наиболее значимым методом исследования стало математическое моделирование, которое позволяет рассмотреть поведение различных реакций в разных условиях без проведения лабораторных экспериментов. Аппарат математического моделирования расширяется и улучшается в прямой зависимости от появления новых задач: создаются новые математические модели, эффективные вычислительные алгоритмы, повышается точность расчетов. Настоящая работа направлена на исследование применения схем высокого порядка точности (WENO-схема 5 и 7 порядка) к решению задач теплопроводности и химических реакций.

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

$\frac{\partial U}{\partial t}+\frac{\partial \left(F\right(U)-H(U\left)\right)}{\partial x}=W\left(U\right).$

Векторы 𝑈, 𝐹(𝑈), 𝐻(𝑈), 𝑊(𝑈) имеют следующий вид:

$\begin{array}{l}U=\left(\begin{array}{l}p{Y}_m\\ pu\\ ph\end{array}\right),F\left(U\right)=\left(\begin{array}{l}pu{Y}_m\\ p{u}^2\\ phu\end{array}\right),H\left(U\right)=\left(\begin{array}{l}{J}_{mx}\\ {\tau }_{xx}\\ q_x\end{array}\right),W\left(U\right)=\left(\begin{array}{l}{R}_m\\ 0\\ 0\end{array}\right).\\ U=\left(\begin{array}{l}p{Y}_m\\ pu\\ ph\end{array}\right),F\left(U\right)=\left(\begin{array}{l}pu{Y}_m\\ p{u}^2\\ phu\end{array}\right),H\left(U\right)=\left(\begin{array}{l}{J}_{mx}\\ {\tau }_{xx}\\ q_x\end{array}\right),W\left(U\right)=\left(\begin{array}{l}{R}_m\\ 0\\ 0\end{array}\right).\end{array}$

В этой системе уравнений 𝑚 = 1, … , M, M - число компонент в газовой смеси, 𝜌 - плотность смеси, 𝑌_m - массовая доля -ой компоненты смеси, u  - скорость, h - энтальпия смеси, J_mx  - диффузионный поток, R_m  - скорость образования или расхода 𝑚-ой компоненты смеси, 𝜋  = 𝑝−𝑝₀ - динамическая составляющая давления, 𝑝₀ - термодинамическая составляющая давления, постоянная в области, 𝜏_xx – вязкий поток, 𝑞_𝑥 - потока тепла.

Система дополняется условием на дивергенцию вектора скорости:

$\nabla \cdot \overrightarrow{u}=\frac{1}{p{C}_{p}T}\left(\nabla \cdot \lambda \nabla T+\sum _{m}p{D}_{m,mix}\nabla Y_m\nabla h_m\right)+\frac{1}{p}\sum _m\frac{M_w}{M_{wm}}(\nabla \cdot p{D}_{m,mix}\nabla Y_m)+\frac{1}{p}\sum _m\left(\frac{M_w}{M_{wm}}-\frac{h_m}{C_{p}T}\right)R_m.$

Здесь C_𝑝 – теплоемкость смеси при постоянном давлении, T – температура смеси, 𝜆 – теплопроводность смеси, D_m,mix – коэффициент диффузии, ℎ_m – энтальпия образования компоненты смеси, M_w – молекулярная масса смеси, M_wm – молекулярная масса компоненты смеси. Более подробно математическая модель представлена в работе [1].

Для построения вычислительного алгоритма используем равномерную сетку отрезков:

${\Omega }_{\Delta x}=\left\{{\Delta }_i,i=\mathrm{1,}\mathrm{...},N,{\Delta }_i=\left[x_{i-1},x_i\right],\left|{\Delta }_i\right|=x_i-x_{i-1}=h_x,h_x{N}_x=L_x\right\}$

Численный алгоритм строится по схеме расщепления по физическим процессам [2]. На первом этапе в этом расщеплении решается система уравнений химической кинетики, на втором этапе интегрируются уравнения законов сохранения без учета давления, далее рассчитывается поле поправок к давлению из решения уравнения Пуассона, на последнем этапе корректируются поле давления и поле скорости [2]. Поскольку мы исследуем применение WENO-схем, приведем подробно алгоритм второго этапа. Для интегрирования законов сохранения используется разностная схема вида:

$\frac{U_i^{n+1}-U_i^n}{\Delta t}+\frac{F_{i+1/2}-F_{i-1/2}}{h_x}-\frac{H_{i+1/2}-H_{i-1/2}}{h_x}=0.$

Здесь H_𝑖+½, H_𝑖−½  – диффузионные, вязкие и тепловые потоки на границе 𝑖 и 𝑖 + 1, 𝑖 − 1 и 𝑖 ячейками соответственно, которые рассчитываются по схеме с центральными разностями, 𝐹𝑖+1/2, 𝐹𝑖−1/2  – конвективные потоки на границе 𝑖 и 𝑖 + 1, 𝑖 − 1 и 𝑖 ячейками соответственно, которые рассчитываются с использованием потоков Русанова [3]:

$\begin{array}{l}{F}_{i+1/2}=0.5\left(F\left(U_{i+1/2}^r\right)+F\left(U_{i+1/2}^l\right)-a\left(U_{i+1/2}^r-U_{i+1/2}^l\right)\right),\\ a=max\left(\left|u_{i+1/2}^r\right|,\left|u_{i+1/2}^l\right|\right),\end{array}$

где  и – значения вектора переменных U справа и слева от границы между i и i+1 ячейками, которые рассчитываем с использованием WENO схем 5 и 7 порядка точности [4; 5]. Их суть заключается в следующем. Для нахождения значения вектора переменных U справа и слева от границы между i и i+1 используется выражение:

$U_{i+1/2}=\sum _{v=1}^K{\Omega }^{\left(v\right)}{U}_{i+1/2}^{\left(v\right)}.$

Здесь K = 3 для случая схемы 5-го порядка, K = 4 для схемы 7-го порядка точности, Ω(v) – оптимальные весовые коэффициенты, $U_{i+1/2}^{\left(v\right)}$ – значение на границе, полученное на шаблоне $S=\left\{x_{i-v},\mathrm{...}\text{ , }{x}_i,\mathrm{...}\text{ , }{x}_{i-v+K-1}\right\}:$

$U_{i+1/2}^{\left(v\right)}=\sum _{p=0}^{K-1}{C}_{vp}{U}_{i-v+p}.$

В случае использования весовых коэффициентов с индикаторами гладкости значения вектора переменных U справа и слева от границы между i и i+1 находятся из выражения:

$U_{i+1/2}=\sum _{v=1}^K{\omega }^{\left(v\right)}{U}_{i+1/2}^{\left(v\right)},{\omega }^{\left(v\right)}=\frac{{\sigma }^{\left(v\right)}}{{\displaystyle {\sum }_{v=1}^K{\sigma }^{\left(v\right)}}},{\sigma }^{\left(v\right)}=\frac{{\Omega }^{\left(v\right)}}{{\left[\epsilon +I{S}^{\left(v\right)}\right]}^p}.$

Здесь IS^(v) – это индикаторы гладкости, 𝜀– некоторое малое число, вводимое, чтоб предотвратить деление на ноль. Величины 𝑐vp, Ω(v), IS^(v) для K = 3, K = 4 приведены в работах [4; 5].

Вычислительные эксперименты. Рассмотрим следующую одномерную постановку задачи. Размер области 0.2 м, шаг по пространству 2 ∙ 10⁻³ м, шаг по времени 1 ∙ 10⁻⁵ с, расчет ведем до 0.2 секунд. Принимаются следующие начальные данные: в области от 0 м до 0.1 м температура газа 800⁰С; в области от 0.1 м до 0.2 м температура газа 1200⁰С; скорость потока 0.1 м/c; давление 101325Па; состав газовой смеси – метан 100% (CH4). Граничные условия: на границе слева со скоростью 0.1 м/c втекает метан с температурой 800⁰С, на границе справа задаются условия вытекания. Такая постановка задачи принята для того, чтобы выяснить какая из схем будет лучше считать в области перепада температур. В этой же области начнут происходить интенсивные реакции с расходом метана и образования продуктов реакции.

Рассматриваем 3 варианта расчета значений вектора U на границе ячеек: 1) схемой первого порядка точности, т.е. $U_{i+1/2}^r=U_{i+1},U_{i+1/2}^l=U_i$; 2) WENO схемой 5-го порядка аппроксимации; 3) WENO-схемой 7-го порядка аппроксимации. На рисунках 1 и 2 представлены распределения температуры и метана. Из графиков можно сделать вывод, что схема первого порядка аппроксимации сглаживает решение в областях резкого изменения газодинамических параметров, разница температур составляет величину 10⁰С. Такой результат является достаточно большим расхождением в случае проведения лабораторных экспериментов. Расчет по схемам WENO 5-го и 7-го порядка практически совпадает, что говорит о преимуществе схемы WENO5, поскольку она использует более компактный шаблон. На рисунках 9 и 10 представлены расчеты по WENO-схемам с оптимальными весами и с расчетом индикаторов гладкости. Из графиков можно сделать вывод, что для рассматриваемых задач использование WENO-схем с оптимальными весами является преимущественным, поскольку графики полностью совпадают, а они являются менее трудоемкими. Полученный результат полного совпадения графиков можно объяснить преобладанием диффузионных процессов над конвективным переносом.

Рис. 1. Распределение температуры. Расчет по схемам 1-го порядка (синяя линия), WENO5 (красная линия), WENO7 (зеленая линия).

Рис. 2. Распределение метана. Расчет по схемам 1-го порядка (синяя линия), WENO5 (красная линия), WENO7 (зеленая линия).

Рис. 3. Распределение температуры. Расчет по схеме WENO5 (красная линия) и схеме WENO5 с оптимальными весами (бирюзовый пунктир).

Рис. 4. Распределение температуры. Расчет по схеме WENO7 (зеленая линия) и схеме WENO7 с оптимальными весами (синий пунктир).

Выводы. В работе проведено исследование применения WENO-схем 5-го и 7-го порядка точности с оптимальными весовыми коэффициентами и с весовыми коэффициентами, рассчитанными с использованием индикаторов гладкости к задачам многокомпонентных реагирующих течений с вязкостью, диффузией и теплопроводностью. Показано, что по сравнению со схемами первого порядка точности WENO-схемы меньше сглаживают решение в областях резкого изменения газодинамических параметров и концентраций компонент смеси. Сравнение результатов вычислительных экспериментов, в которых весовые коэффициенты в WENO-схеме принимались оптимальными и рассчитывались с индикаторами гладкости показало, что результаты расчетов совпадают. Таким образом, можно сделать вывод о преимуществе использования оптимальных весовых коэффициентов в WENO-схеме, поскольку в данном случае алгоритм является менее трудоемким.