Журнал вычислительной математики и математической физики, 2022, T. 62, № 8, стр. 1428-1444

Экспоненциальные многосеточные методы сквозного счета для стационарных сжимаемых потоков с ударными волнами

Шу-Джи Ли^*

Пекинский исследовательский центр вычислительных наук
Пекин, Китай

Поступила в редакцию 10.10.2021
После доработки 21.01.2022
Принята к публикации 11.04.2022

EDN: BRYNKK
DOI: 10.31857/S0044466922080087

Аннотация

Предлагается надежная и эффективная экспоненциальная многосеточная схема для численного моделирования стационарных сжимаемых течений. В алгоритме, основанном на экспоненциальной схеме интегрирования по времени c глобальной связью, обеспечивается подавление ошибки для ускорения сходимости к устойчивому состоянию, а пространственные моды ошибок высокой частоты и высокого порядка сглаживаются с помощью $s$-ступенчатого предобусловленного метода Рунге–Кутты. Показано, что полученная экспоненциальная многосеточная схема эффективна для гладких течений и может стабилизировать сквозной счет для ударных волн средней амплитуды без лимитеров и без добавления искусственной диссипации. Библ. 21. Фиг. 16.

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

1. ВВЕДЕНИЕ

Важным требованием вычислительной гидродинамики является возможность точного моделирования устойчивых течений, например, при обтекании тела, чтобы можно было оценить основные рабочие параметры, например, коэффициенты подъемной силы и сопротивления. По сравнению с традиционными методами второго порядка, методы высокого порядка, такие как разрывные методы Галеркина, могут предложить преимущества в точности, но недостаток производительности и проблемы с устойчивостью. Для снижения вычислительных затрат разрывные методы Галеркина могут быть успешно ускорены с помощью $p$-многосеточных методов, однако ранее такой подход применялся только для задач с гладкими течениями (см. [1]–[5]). Еще не проводилось исследование $p$-многосеточных методов в задачах сквозного счета ударной волны. В настоящей работе рассматриваются задачи, включающие как гладкие течения, так и ударные волны, и исследуется возможность применения $p$-многосеточной методологии к сквозному счету ударной волны. Учитывая вопросы надежности и эффективности методов интегрирования по времени многосеточных сглаживающих методов, выбрана схема экспоненциального временного интегратора с предиктором-корректором (PCEXP) (см. [6]–[10]), которая демонстрирует некоторые преимущества в плане точности и эффективности как для установившихся, так и для нестационарных потоков. Схема PCEXP используется в рамках $p$-многосеточного метода, который включает в себя, в частности, метод экспоненциального интегрирования по времени и $s$-ступенчатый предобусловленный метод Рунге–Кутты. Полученный алгоритм тестируется как эффективный способ повышения реализуемости разрывных методов Галеркина высокого порядка для моделирования стационарных сжимаемых течений.

Работа организована следующим образом. В разд. 2 представлен вариант многосеточного алгоритма, который объединяет два отдельных метода в V-цикле $p$-многосеточного метода. В разд. 3 представлена пространственная дискретизация с помощью модального разрывного метода Галеркина высокого порядка. В разд. 4 обсуждается вычисление шагов по времени в $p$-многосеточном методе. В разд. 5 представлены численные результаты для гладких и разрывных течений. Наконец, в разд. 6 сделаны выводы по данной работе.

2. ЭКСПОНЕНЦИАЛЬНЫЙ МНОГОСЕТОЧНЫЙ МЕТОД

В этом разделе представлен многосеточный метод произвольного порядка ($p$-го порядка). Метод объединяет два отдельных алгоритма: алгоритм экспоненциального интегрирования по времени и $s$-стадийный предобусловленный метод Рунге–Кутты. Эти два алгоритма представлены по отдельности, а затем интегрированы в целостный мнoгосеточный метод с использованием V-цикла.

2.1. Экспоненциальное интегрирование по времени

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

(1)

$\frac{{{\text{d}}{\mathbf{u}}}}{{{\text{d}}t}} = {\mathbf{R}}({\mathbf{u}}),$

где ${\mathbf{u}} = {\mathbf{u}}(t) \in {{\mathbb{R}}^{K}}$ – вектор переменных решений и ${\mathbf{R}}({\mathbf{u}}) \in {{\mathbb{R}}^{K}}$ – член правой части, вычисленный с помощью произвольной дискретизации в пространстве. Без потери общности мы рассматриваем ${\mathbf{u}}(t)$ в одном шаге по времени, т.е. $t \in [{{t}_{n}},{{t}_{{n + 1}}}]$.

Разделение правой части приводит к другому выражению:

(2)

$\frac{{{\text{d}}{\mathbf{u}}}}{{{\text{d}}t}} = {{{\mathbf{J}}}_{n}}{\mathbf{u}} + {\mathbf{N}}({\mathbf{u}}),$

где индекс $n$ указывает на значение, вычисленное при $t = {{t}_{n}}$; ${{{\mathbf{J}}}_{n}}$ – матрица Якоби, образованная глобальной дискретизированной системой, размерность матрицы Якоби для глобальной невязки при этом равна суммарному числу степеней свободы глобальной системы, а именно, ${{\left. {{{{\mathbf{J}}}_{n}} = \partial {\mathbf{R}}({\mathbf{u}}){\text{/}}\partial {\mathbf{u}}} \right|}_{{t = {{t}_{n}}}}} = \partial {\mathbf{R}}({{{\mathbf{u}}}_{n}}){\text{/}}\partial {\mathbf{u}};$ ${\mathbf{N}}({\mathbf{u}}) = {\mathbf{R}}({\mathbf{u}}) - {{{\mathbf{J}}}_{n}}{\mathbf{u}}$ обозначает невязку, которая в общем случае является нелинейной. Уравнение (2) допускает следующее формальное решение:

(3)

${{{\mathbf{u}}}_{{n + 1}}} = \exp (\Delta t{{{\mathbf{J}}}_{n}}){{{\mathbf{u}}}_{n}} + \int\limits_0^{\Delta t} {\exp \left( {(\Delta t - \tau ){{{\mathbf{J}}}_{n}}} \right)} {\mathbf{N}}({\mathbf{u}}({{t}_{n}} + \tau )){\text{d}}\tau ,$

где матричная экспонента определяется как

(4)

$\exp ( - t{{{\mathbf{J}}}_{n}}) = \sum\limits_{m = 0}^\infty {\frac{{{{{\left( { - t{{{\mathbf{J}}}_{n}}} \right)}}^{m}}}}{{m!}}} .$

В формуле (3) жесткий линейный компонент и нелинейный интегральный компонент могут быть оценены отдельно. В результате линейный компонент может быть выражен аналитически для некоторых специализированных уравнений, но нелинейный компонент обычно аппроксимируется численно. Недавно была разработана двухступенчатая экспоненциальная схема PCEXP (см. [10]) для эффективного вычисления различных областей течения как для установившихся, так и для нестационарных течений. Но для установившихся течений схема экспоненциального интегрирования по времени первого порядка (EXP1), полученная путем аппроксимации первого порядка интеграла по времени нелинейного члена (см. [6], [8], [11]), оказывается более эффективной, несмотря на ее неспособность разрешать нестационарные течения. Схема EXP1 имеет две эквивалентных формы:

(5)

${{{\mathbf{u}}}_{{n + 1}}} = \exp (\Delta t{{{\mathbf{J}}}_{n}}){{{\mathbf{u}}}_{n}} + \Delta t{{\Phi }_{1}}(\Delta t{{{\mathbf{J}}}_{n}}){{{\mathbf{N}}}_{n}},$

(6)

${{{\mathbf{u}}}_{{n + 1}}} = {{{\mathbf{u}}}_{n}} + \Delta t{{\Phi }_{1}}(\Delta t{{{\mathbf{J}}}_{n}}){{{\mathbf{R}}}_{n}},$

где

(7)

${{\Phi }_{1}}(\Delta t{\mathbf{J}}): = \frac{{{{{\mathbf{J}}}^{{ - 1}}}}}{{\Delta t}}\left[ {\exp \left( {\Delta t{\mathbf{J}}} \right) - {\mathbf{I}}} \right],$

а ${\mathbf{I}}$ есть $K \times K$-единичная матрица. Далее используется формула (6).

Физическая природа таких типов экспоненциальных схем опирается на свойство глобальной связи через глобальную матрицу Якоби ${\mathbf{J}}$, так что информация о переносе в течении может транслироваться на всю вычислительную область без ограничения на глобальное число Куранта–Фридрихса–Леви (CFL). Поэтому экспоненциальные схемы ведут себя как полностью неявные методы, но зависят только от текущего решения, т.е. подобны явному методу, как показано в уравнении (6). Поэтому схема EXP1 с такой вычислительной особенностью используется в рамках $p$-многосеточного метода для расчетов установившегося потока.

2.2. Реализация EXP1 с помощью метода Крылова

Реализация схем экспоненциального интегрирования по времени требует вычислений матрично-векторных произведений, в частности, произведения экспоненциальных функций матрицы Якоби на вектор, например, ${{\Phi }_{1}}(\Delta t{{{\mathbf{J}}}_{n}}){\mathbf{N}}$ в (6). Для их эффективного вычисления может быть использован метод Крылова (см. [12], [13]). Рассмотрим $m$-мерное подпространство Крылова

(8)

${{\mathbb{K}}_{m}}({\mathbf{J}},{\mathbf{N}}) = {\text{span}}\left\{ {{\mathbf{N}},{\mathbf{JN}},{{{\mathbf{J}}}^{2}}{\mathbf{N}},\; \ldots ,\;{{{\mathbf{J}}}^{{m - 1}}}{\mathbf{N}}} \right\}.$

Ортогональная базисная матрица ${{{\mathbf{V}}}_{m}}: = ({{{\mathbf{v}}}_{1}},\; \ldots ,\;{{{\mathbf{v}}}_{m}}) \in {{\mathbb{R}}^{{K \times m}}}$ удовлетворяет так называемому разложению Арнольди (см. [13]):

(9)

${\mathbf{J}}{{{\mathbf{V}}}_{m}} = {{{\mathbf{V}}}_{{m + 1}}}{{{\mathbf{\tilde {H}}}}_{m}},$

где ${{{\mathbf{V}}}_{{m + 1}}}: = ({{{\mathbf{v}}}_{1}},\; \ldots ,\;{{{\mathbf{v}}}_{m}},{{{\mathbf{v}}}_{{m + 1}}}) = ({{{\mathbf{V}}}_{m}},{{{\mathbf{v}}}_{{m + 1}}}) \in {{\mathbb{R}}^{{K \times (m + 1)}}}$. Матрица ${{{\mathbf{\tilde {H}}}}_{m}}$ является $(m + 1) \times m$ верхней матрицей Гессенберга. Тогда уравнение (9) становится

(10)

${\mathbf{J}}{{{\mathbf{V}}}_{m}} = {{{\mathbf{V}}}_{m}}{{{\mathbf{H}}}_{m}} + {{h}_{{m + 1,m}}}{{{\mathbf{v}}}_{{m + 1}}}{\mathbf{e}}_{m}^{{\mathbf{T}}}.$

Поскольку ${\mathbf{V}}_{m}^{{\mathbf{T}}}{{{\mathbf{V}}}_{m}} = {\mathbf{I}}$,

(11)

${{{\mathbf{H}}}_{m}} = {\mathbf{V}}_{m}^{{\mathbf{T}}}{\mathbf{J}}{{{\mathbf{V}}}_{m}}.$

Таким образом, ${{{\mathbf{H}}}_{m}}$ является проекцией линейного преобразования ${\mathbf{J}}$ на подпространство ${{\mathbb{K}}_{m}}$ с базисом ${{\mathbb{V}}_{m}}$. Поскольку ${{{\mathbf{V}}}_{m}}{\mathbf{V}}_{m}^{{\mathbf{T}}} \ne {\mathbf{I}}$, уравнение (11) приводит к приближению

(12)

${\mathbf{J}} \approx {{{\mathbf{V}}}_{m}}{\mathbf{V}}_{m}^{{\mathbf{T}}}{\mathbf{J}}{{{\mathbf{V}}}_{m}}{\mathbf{V}}_{m}^{{\mathbf{T}}} = {{{\mathbf{V}}}_{m}}{{{\mathbf{H}}}_{m}}{\mathbf{V}}_{m}^{{\mathbf{T}}},$

и $\exp ({\mathbf{J}})$ можно аппроксимировать с помощью $\exp ({{{\mathbf{V}}}_{m}}{{{\mathbf{H}}}_{m}}{\mathbf{V}}_{m}^{{\mathbf{T}}})$ следующим образом:

(13)

$\exp ({\mathbf{J}}){\mathbf{N}} \approx \exp \left( {{{{\mathbf{V}}}_{m}}{{{\mathbf{H}}}_{m}}{\mathbf{V}}_{m}^{{\mathbf{T}}}} \right){\mathbf{N}} = {{{\mathbf{V}}}_{m}}\exp ({{{\mathbf{H}}}_{m}}){\mathbf{V}}_{m}^{{\mathbf{T}}}{\mathbf{N}}.$

Первый вектор-столбец ${{{\mathbf{V}}}_{m}}$ равен ${{{\mathbf{v}}}_{1}} = {{\mathbf{N}} \mathord{\left/ {\vphantom {{\mathbf{N}} {{{{\left\| {\mathbf{N}} \right\|}}_{2}}}}} \right. \kern-0em} {{{{\left\| {\mathbf{N}} \right\|}}_{2}}}}$ и ${\mathbf{V}}_{m}^{{\mathbf{T}}}{\mathbf{N}} = {{\left\| {\mathbf{N}} \right\|}_{2}}{{{\mathbf{e}}}_{1}}$, таким образом, (13) становится

(14)

$\exp ({\mathbf{J}}){\mathbf{N}} \approx {{\left\| {\mathbf{N}} \right\|}_{2}}{{{\mathbf{V}}}_{m}}\exp ({{{\mathbf{H}}}_{m}}){{{\mathbf{e}}}_{1}}.$

Следовательно, ${{\Phi }_{1}}$ можно аппроксимировать как

(15)

${{\Phi }_{1}}(\Delta t{\mathbf{J}}){\mathbf{N}} = \frac{1}{{\Delta t}}\int\limits_0^{\Delta t} {\exp } ((\Delta t - \tau ){\mathbf{J}}){\mathbf{N}}{\text{d}}\tau \approx \frac{1}{{\Delta t}}\int\limits_0^{\Delta t} {{{{\left\| {\mathbf{N}} \right\|}}_{2}}} {{{\mathbf{V}}}_{m}}\exp \left( {(\Delta t - \tau ){{{\mathbf{H}}}_{m}}} \right){{{\mathbf{e}}}_{1}}{\text{d}}\tau .$

В общем случае размерность подпространства Крылова $m$ выбирается намного меньше размерности ${\mathbf{J}}$, поэтому ${{{\mathbf{H}}}_{m}} \in {{\mathbb{R}}^{{m \times m}}}$ легко оценить. Таким образом, ${{\Phi }_{1}}$ вычисляется как

(16)

${{\Phi }_{1}}(\Delta t{\mathbf{J}}){\mathbf{N}} \approx \frac{1}{{\Delta t}}{{\left\| {\mathbf{N}} \right\|}_{2}}{{{\mathbf{V}}}_{m}}\int\limits_0^{\Delta t} {\exp \left( {(\Delta t - \tau ){{{\mathbf{H}}}_{m}}} \right){{{\mathbf{e}}}_{1}}{\text{d}}\tau } = \frac{1}{{\Delta t}}{{\left\| {\mathbf{N}} \right\|}_{2}}{{{\mathbf{V}}}_{m}}{\mathbf{H}}_{m}^{{ - 1}}\left[ {\exp (\Delta t{{{\mathbf{H}}}_{m}}) - {\mathbf{I}}} \right]{{{\mathbf{e}}}_{1}},$

где экспонента матрицы $\exp (\Delta t{{{\mathbf{H}}}_{m}})$ может быть эффективно вычислена рациональным приближением Чебышёва (ср. [12], [13]).

2.3. Предобусловленный метод Рунге–Кутты

Рассмотрим $s$-стадийный предобусловленный метод Рунге–Кутты (PRK) в следующей форме:

${{{\mathbf{u}}}^{{(0)}}} = {{{\mathbf{u}}}^{n}},$

(17)

${{{\mathbf{u}}}^{{(k)}}} = {{{\mathbf{u}}}^{n}} + {{\beta }_{k}}{{{\mathbf{P}}}^{{ - 1}}}\left( {{{{\mathbf{u}}}^{n}}} \right){\mathbf{R}}\left( {{{{\mathbf{u}}}^{{(k - 1)}}}} \right),\quad k = 1,\;2,\; \ldots ,\;s,$

${{{\mathbf{u}}}^{{n + 1}}} = {{{\mathbf{u}}}^{{(s)}}},$

где ${{\beta }_{k}} = 1{\text{/}}(s - k + 1)$, $P$ берется как диагональная часть матрицы Якоби невязки ${\mathbf{J}} = \partial R{\text{/}}\partial u$, которая распределяет информацию о распространении волны по элементам сетки. В данной работе для всех тестовых случаев используется $s = 4$.

Физическая природа метода PRK может быть интерпретирована с двух различных точек зрения. Во-первых, мы рассматриваем пространственную дискретизацию первого порядка методом конечных объемов или, во-вторых, разрывным методом Галеркина для $i$-го элемента, окруженного смежными ячейками $j$ ($1 \leqslant j \leqslant N$) с площадью разделяющей поверхности, равной ${{{\text{S}}}_{{ij}}}$. Понятие матричного шага по времени можно раскрыть, рассмотрев пространственную невязку, полученную с помощью какой-то схемы на основе распада разрыва Римана. Без потери общности мы выбрали схему Рое (см. [15]) в качестве примера для лучшего представления связи между скалярным и матричным временными шагами. Отметим, что в практических расчетах конвективный поток вычисляется римановым решателем SLAU2 (см. [16]) из соображений надежности:

(18)

${{V}_{i}}\frac{{\Delta {{u}_{i}}}}{{\Delta t}} = {{R}_{i}} = \sum\limits_{j = 1}^N {\frac{1}{2}} \left[ {F({{u}_{i}}) + F({{u}_{j}})} \right]{{n}_{{ij}}}{{{\text{S}}}_{{ij}}} + \frac{1}{2}\left| {A_{{ij}}^{n}} \right|\left( {{{u}_{i}} - {{u}_{j}}} \right){{{\text{S}}}_{{ij}}}.$

Диагональная часть матрицы Якоби $P$ глобальной невязки определяется как

(19)

${{P}_{i}} = \frac{{\partial {{R}_{i}}}}{{\partial {{u}_{i}}}} \approx \sum\limits_{j = 1}^N \frac{1}{2}\left| {A_{{ij}}^{n}} \right|{{S}_{{ij}}}.$

Определяя матрицу $\Delta {\mathbf{t}}$ как

(20)

$\Delta {\mathbf{t}} = {{V}_{i}}{{P}^{{ - 1}}} = \frac{{{{V}_{i}}}}{{\sum\limits_{j = 1}^N {\frac{1}{2}\left| {A_{{ij}}^{n}} \right|{{{\text{S}}}_{{ij}}}} }},$

можно раскрыть связь между матрицей $\Delta {\mathbf{t}}$ и традиционным скалярным определением временного шага, рассмотрев постоянную в ячейках скалярную аппроксимацию спектрального радиуса $\lambda _{{ij}}^{{{\text{max}}}}$ к $\left| {A_{{ij}}^{n}} \right|$, т.е.

(21)

$\Delta t = \frac{{2{{V}_{i}}}}{{\sum\limits_{j = 1}^N {\lambda _{{ij}}^{{{\text{max}}}}} {{{\text{S}}}_{{ij}}}}}\;\mathop \to \limits^{1D} \;\frac{{\Delta {{x}_{i}}}}{{\lambda _{i}^{{{\text{max}}}}}}.$

Поэтому ${{P}^{{ - 1}}}$ эквивалентно матричному шагу по времени, что согласуется с обычным определением шага по времени в скалярном случае. Как было продемонстрировано в [17], эта матрица является предобусловливателем, который может обеспечить эффективную кластеризацию конвективных собственных значений и существенное улучшение сходимости временного шага метода Рунге–Кутты. В данной работе, в отличие от [17], используется точный способ оценки временных шагов матрицы с точной матрицей Якоби (ср. (8)), так что все эффекты жесткости пространственной дискретизации и граничных условий могут быть точно учтены. Поскольку традиционного скалярного физического шага по времени не существует, схема PRK не имеет временного порядка точности. С другой точки зрения, схему PRK можно рассматривать как упрощенную полностью неявную схему или как неявно-явные методы Рунге–Кутты, игнорируя недиагональные члены матрицы Якоби. Таким образом, это локализованный метод с малым объемом памяти. Для повышения надежности мы рекомендуем увеличить диагональное доминирование до $P$, а именно, $P = \partial R{\text{/}}\partial u + I{\text{/}}\delta \tau $, где $\delta \tau $ можно рассматривать как псевдо-временной шаг, вычисляемый в (33).

2.4. V-цикл $p$p-многосеточного метода

В этой многосеточной схеме $p$-многосеточный метод применяется к пространственной дискретизации с порядком аппроксимации $p$. Таким образом, количество уровней многосеточного метода равно пространственному порядку аппроксимации $p$. Временной шаг применяется к каждому уровню многосеточного V-цикла, где решение сглаживается в течение одного шага по времени. Весь V-цикл состоит из операторов ограничения и продолжения, которые также применяются к одному шагу по времени. Для сглаживания решения можно использовать либо явные методы Рунге–Кутты, либо PRK. Однако они оказываются неэффективными для устранения низкочастотных мод ошибок при низких порядках точности. В данной работе рассматривается схема EXP1 с учетом ее численных особенностей. В отличие от явного сглаживателя Рунге–Кутты, который дает слабый эффект затухания только локально (точечно), схема EXP1 относится к методам глобальной связи, что позволяет использовать большие шаги по времени с сильным эффектом затухания для всех частотных мод (см. [8]).

В экспоненциальном $p$-многосеточном методе (eMG) схема EXP1 используется на уровне точности $p = 0$, а метод PRK – на уровне точности $p > 0$. При многосеточном сглаживании используется V-цикл $p$-многосеточного процесса, где рекурсивно используется двухуровневый алгоритм. Для иллюстрации алгоритма рассмотрим нелинейную задачу $A\left( {{{u}^{p}}} \right) = {{f}^{p}}$, где ${{u}^{p}}$ – вектор решения, $A\left( {{{u}^{p}}} \right)$ – нелинейный оператор, ${{f}^{p}}$ – член правой части, а $p$ обозначает уровень точности разрывного метода Галеркина. Пусть ${{v}^{p}}$ – приближение к вектору решения. Тогда $A\left( {{{u}^{p}}} \right)$ и ${{f}^{p}}$ получены пространственной дискретизацией разрывного метода Галеркина в сочетании со схемами EXP1 или PRK. Мы используем схему EXP1 в качестве сглаживающей при $p = 0$ и PRK при $p > 0$. Невязка $r\left( {{{v}^{p}}} \right)$ определяется как

$r\left( {{{v}^{p}}} \right) = {{f}^{p}} - {{A}^{p}}\left( {{{v}^{p}}} \right).$

В рамках eMG решение на уровне $p - 1$ используется для корректировки решения на уровне $p$ следующим образом.

1. Проводится шаг по времени со схемой PRK на самом высоком уровне точности ${{p}_{{\max }}}$.

2. Решение и остаток $p$ ограничиваются уровнем $p - 1$ ($1 \leqslant p \leqslant {{p}_{{\max }}}$):

(22)

$v_{0}^{{p - 1}} = \mathbb{R}_{p}^{{p - 1}}{{v}^{p}},\quad {{r}^{{p - 1}}} = \mathbb{R}_{p}^{{p - 1}}{{r}^{p}}({{v}^{p}}),$

где $\mathbb{R}_{p}^{{p - 1}}$ – оператор ограничения с уровня $p$ на уровень $p - 1$.

3. Вычисляется форсирующий член для уровня $p - 1$:

(23)

${{s}^{{p - 1}}} = {{A}^{{p - 1}}}\left( {v_{0}^{{p - 1}}} \right) - {{r}^{{p - 1}}}.$

4. Решение сглаживается с помощью схемы PRK на уровне $p - 1$, но переключается на использование схемы EXP1 на самом низком уровне точности $p = 0$:

(24)

${{A}^{{p - 1}}}\left( {{{v}^{{p - 1}}}} \right) = \mathbb{R}_{p}^{{p - 1}}{{f}^{p}} + {{s}^{{p - 1}}}.$

5. Оценивается ошибка уровня $p - 1$:

(25)

${{e}^{{p - 1}}} = {{v}^{{p - 1}}} - v_{0}^{{p - 1}}.$

6. Продолжается ошибка $p - 1$ и корректируется аппроксимация уровня $p$:

(26)

${{v}^{p}} = {{v}^{p}} + \mathbb{P}_{{p - 1}}^{p}{{e}^{{p - 1}}},$

где $\mathbb{P}_{{p - 1}}^{p}$ – оператор продолжения.

V-Цикл $p$-многосеточного метода проводит попеременное глобальное и локальное сглаживание, которое оказывается эффективным как для гладких, так и разрывных полей потока.

2.5. Сквозной счет

Установлено, что eMG обладает естественным механизмом стабилизации вычислений, связанных с разрешением ударной волны, так что возможно применение методов, не использующих лимитеры на разрывах. В eMG два источника численного демпфирования обеспечиваются EXP1 и $p$-многосеточными схемами. Обе схемы могут проводить большую численную диссипацию вдали от устойчивого состоянии, повышая устойчивость сквозного счета. Диссипация постепенно исчезает при приближении к стационарному состоянию решения, так что сходящиеся решения могут сохранять очень острые профили ударной волны. Такое поведение численного метода исследуется в разд. 5 для трансзвуковых и сверхзвуковых областей течения, в которых могут быть получены очень резкие профили разрывов.

3. ПРОСТРАНСТВЕННАЯ ДИСКРЕТИЗАЦИЯ

В данной работе eMG применяется для решения трехмерных сжимаемых уравнений Эйлера, дискретизированных модальным разрывным методом Галеркина (см. [6]–[10]).

Рассмотрим трехмерные сжимаемые уравнения Эйлера

(27)

$\frac{{\partial {\mathbf{U}}}}{{\partial t}} + \nabla \cdot {\mathbf{F}} = 0.$

Здесь ${\mathbf{U}}$ – вектор консервативных переменных, а ${\mathbf{F}}$ – конвективный поток

(28)

${\mathbf{U}} = \left( {\begin{array}{*{20}{c}} \rho \\ {\rho v} \\ {\rho E} \end{array}} \right),\quad {\mathbf{F}}\left( {\begin{array}{*{20}{c}} {\rho v} \\ {\rho v \otimes v + p{\mathbf{I}}} \\ {\rho Hv} \end{array}} \right),$

где $v = (u,v,w{{)}^{{\mathbf{T}}}}$ – вектор абсолютной скорости; $\rho $, $p$ и $e$ – плотность потока, давление и удельную внутреннюю энергию; $E = e + \frac{1}{2}{{\left\| v \right\|}^{2}}$ и $H = E + p{\text{/}}\rho $ – полная энергия и полная энтальпия соответственно; ${\mathbf{I}}$ – тождественный тензор; а давление $p$ задается уравнением состояния идеального газа

(29)

$p = \rho \left( {\gamma - 1} \right)e,$

где $\gamma = 7{\text{/}}5$ – отношение удельных теплоемкостей для идеального газа.

3.1. Модальный разрывный метод Галеркина

Рассматривая вычислительную область $\Omega $, разделенную на множество непересекающихся элементов произвольной формы, модальный разрывный метод Галеркина (см. [8]) ищет аппроксимацию ${{{\mathbf{U}}}_{h}}$ в каждом элементе $E \in \Omega $ с помощью конечномерного пространства полиномов ${{P}^{p}}$ порядка $p$ в пространстве разрывных конечных элементов:

(30)

${{\mathbb{V}}_{h}}: = \left\{ {{{\psi }_{i}} \in {{L}^{2}}(\Omega )\,:\,{{{\left. {{{\psi }_{i}}} \right|}}_{E}} \in {{P}^{p}}(\Omega ),\;\forall E \in \Omega } \right\}.$

Численное решение ${{{\mathbf{U}}}_{h}}$ может быть аппроксимировано в пространстве конечных элементов ${{\mathbb{V}}_{h}}$:

(31)

${{{\mathbf{U}}}_{h}}(x,t) = \sum\limits_{j = 1}^n {{{{\mathbf{u}}}_{j}}} (t){{\psi }_{j}}(x).$

В слабой формулировке уравнения Эйлера (27) на элементе $E$ можно выписать следующие соотношения:

(32)

$\int\limits_E {{{\psi }_{i}}} {{\psi }_{j}}dx\frac{{d{{{\mathbf{u}}}_{j}}}}{{dt}} = - \int\limits_{\partial E} {{{\psi }_{i}}} \widetilde {\mathbf{F}} \cdot \hat {n}d\sigma + \int\limits_E {\mathbf{F}} \cdot \nabla {{\psi }_{i}}dx: = {{{\mathbf{R}}}_{i}},$

где $\hat {n}$ – внешний единичный вектор элемента поверхности $\sigma $ ячейки сетки $E$, ${\mathbf{\tilde {F}}}$ оценивается решателем Римана (см. [18]), и в данной формуле используется суммирование Эйнштейна. Для ортонормального базисного набора $\{ {{\psi }_{i}}\} $ член в левой части уравнения (32) становится диагональным, поэтому система не отличается от стандартной формы обыкновенных дифференциальных уравнений (1), что позволяет избежать решения линейной системы, требуемой неортогональным базисом. Еще более важно, что использование ортогонального базиса дает более точные решения, особенно для методов высокого порядка, например, $p = 6$.

4. СТРАТЕГИЯ ВЫБОРА ШАГА ВО ВРЕМЕНИ

В этом разделе обсуждается стратегия выбора шага по времени в рамках метода eMG. Необходимо вычислить два различных временных шага. Один для временного шага PRK $\delta \tau $, а другой – для сглаживания EXP1, последний эмпирически выбран как $({{p}_{{\max }}} + 1){\kern 1pt} \delta \tau $. Таким образом, необходимо определить только $\delta \tau $:

(33)

$\delta \tau = \frac{{\operatorname{CFL} {{h}_{{{\text{3}}D}}}}}{{\left( {2p + 1} \right)\left( {\left\| v \right\| + c} \right)}},\quad {{h}_{{{\text{3}}D}}}: = 2d\frac{{\left| E \right|}}{{\left| {\partial E} \right|}},$

где $p$ – уровень точности, $v$ – вектор скорости в центре ячейки, $c$ – скорость звука, $d$ – пространственная размерность, $\left| E \right|$ и $\left| {\partial E} \right|$ – объем ячейки и площадь поверхности границы $E$ соответственно; ${{h}_{{{\text{3D}}}}}$ – характерная длина ячейки в трехмерном пространстве, которая определяется отношением объема клетки и площади поверхности. В решателе HA3D квазидвумерные вычисления реализуются путем экструзии двумерной сетки в трехмерную с помощью одного слоя ячеек. Мы используем ${{h}_{{{\text{2D}}}}}$ вместо ${{h}_{{{\text{3d}}}}}$, чтобы устранить влияние измерения $z$ на получение истинно двумерного временного шага. Учитывая размер ячейки $\Delta z$ в направлении $z$, ${{h}_{{{\text{2D}}}}}$ определяется как

(34)

$\frac{2}{{{{h}_{{{\text{2D}}}}}}} = \frac{3}{{{{h}_{{{\text{3D}}}}}}} - \frac{1}{{\Delta z}}.$

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

(35a)

${\text{CF}}{{{\text{L}}}_{n}} = \min \left\{ {{\text{CF}}{{{\text{L}}}_{{\max }}},\;\max \left[ {\left\| {R({{\rho }_{n}})} \right\|_{2}^{{ - 1}},\;1 + \frac{{(n - 1)}}{{(2p + 1)}}} \right]} \right\},$

(35b)

${{\left\| {R({{\rho }_{n}})} \right\|}_{2}}: = \frac{1}{{\left| \Omega \right|}}{{\left[ {\int\limits_\Omega R {{{({{\rho }_{n}})}}^{2}}dx} \right]}^{{1/2}}},$

где $R({{\rho }_{n}})$ – невязка уравнения для плотности, ${\text{CF}}{{{\text{L}}}_{{{\text{max}}}}}$ – определяемое пользователем максимальное число CFL, $n$ – число итераций для неявного метода или многосеточных циклов (MG-циклов) для схемы eMG, $p$ – пространственный порядок точности. Такая переменная стратегия эволюции CFL обеспечивает надежный запуск кода и хорошую общую вычислительную эффективность на практике. Во всех рассмотренных тестовых случаях верхнее граничное число ${\text{CF}}{{{\text{L}}}_{{\max }}}$ принимается равным ${{10}^{3}}$ для неявного метода Эйлера и ${{10}^{2}}$ для метода eMG.

5. ЧИСЛЕННЫЕ РЕЗУЛЬТАТЫ

Численные результаты представлены для расчетов как гладких, так и разрывных течений. Производительность и сходимость системы eMG исследуются для двух случаев гладкого потока: дозвукового обтекания кругового цилиндра и сферы. Метод eMG сравнивается с методом GMRES c быстрым полностью неявным ILU-предобусловливателем относительно процессорного времени и использования памяти. Оба метода используют одну и ту же размерность подпространства Крылова 30 и критерий останова ${{10}^{{ - 5}}}$ (подробнее см. [13], [19]). Для неявного интегрирования по времени используется неявный метод Эйлера первого порядка. Для сквозного счета рассматриваются сверхзвуковое сопло Лаваля при числе Маха 3 и трансзвуковое крыло M6 при числе Маха 0.8395. Исследуются надежность и устойчивость сходимости метода eMG.

5.1. Гладкое обтекание кругового цилиндра

Представим результаты, полученные для обтекания кругового цилиндра при числе Маха $0.3$. Цилиндр имеет радиус $1$ и окружен круговой расчетной областью радиуса $15$, как показано на фиг. 1. Квази-2D сетка с $128 \times 32 = 4096$ квадратичными изогнутыми гексаэдральными элементами создается путем экструзии 2D сетки на один слой сеток. Окончательная 3D сетка используется решателем потока HA3D разрывного метода Галеркина произвольного высокого порядка (см. [6]–[10]). Общее число степеней свободы достигает 81 920 при $p = 3$.

Фиг. 1.

Обтекание кругового цилиндра в квази-2D: 128 $ \times $ 32 = 4.096 квадратичных криволинейных элементов.

На фиг. 2 ${{L}_{2}}$-норма невязки по плотности $R({{\rho }_{n}})$ построена в зависимости от количества циклов многосеточного метода (MG-циклов) с использованием схемы eMG. Показано, что скорость сходимости не зависит от пространственного порядка точности $p$, т.е. $p$-независима. Результаты вычислены с помощью быстрого неявного метода GMRES с предобусловливателем ILU и приведены на фиг. 3, где показаны истории сходимости неявного метода с различной пространственной точностью. Видна быстрая квадратичная ньютоновская сходимость, однако скорость сходимости зависит от пространственного порядка точности $p$. Схема eMG сравнивается с полностью неявным методом относительно процессорного времени на фиг. 4, где процессорное время нормировано на время схемы eMG. Видно, что неявный метод (IMP) быстрее для случаев $p = 1,2$, но медленнее, чем схема eMG для случая $p = 3$. Метод eMG, по крайней мере, сравним с неявным методом по общей производительности.

Фиг. 2.

$p$-Независимые сходимости в методе eMG.

Фиг. 3.

Истории сходимости неявного метода с различной пространственной точностью.

Фиг. 4.

Сравнение производительности между методом eMG и неявным методом.

5.2. Гладкое обтекание сферы

Вычислительная эффективность схемы eMG исследуется для трехмерного потока, обтекающего сферу с числом Маха 0.3, радиусом сферы $1$ и радиусом внешней границы $5$. На поверхности сферы задано граничное условие стенки скольжения, а на внешней границе используется характерное граничное условие дальнего поля с инвариантами Римана. Сетка соблюдает симметрии потока в горизонтальной и вертикальной плоскостях, на которые накладывается граничное условие симметрии. Созданная криволинейная сетка состоит из 9778 тетраэдров и 4248 призм, всего 14 026 ячеек. На фиг. 5 показан крупный план сетки сферы и контура скорости, рассчитанного по схеме eMG при $p = 3$.

Фиг. 5.

Контуры течения, рассчитанного для обтекания сферы при ${\text{Ma}} = 0.3$, $p = 3$ с помощью метода eMG.

На фиг. 6 приведены истории сходимости метода eMG для пространственного порядка точности $p$ до $3$. Вновь наблюдается $p$-независимая сходимость, что аналогично двумерному случаю. На фиг. 7 показаны истории сходимости IMP с отсчетами многосеточных циклов, которые зависят от $p$. На фиг. 8 сравниваются оба метода, измеренные в процессорном времени. Видно, что хотя IMP является быстрым по количеству итераций, вычислительные затраты на каждую итерацию относительно высоки, поэтому общее время процессора ухудшается. При использовании пространственных схем высокого порядка вместе с неявным методом глобальная матрица Якоби высокого порядка занимает большой объем памяти. Наиболее значительная часть использования памяти методами eMG и IMP оценивается следующим образом:

${{{\text{M}}}_{{{\text{emg}}}}} = {\text{NE}}\left[ {\frac{5}{3}(p + 1)(p + 2)(p + 3) + 150} \right],$

${{{\text{M}}}_{{{\text{imp}}}}} = 6{\text{NE}}{{\left[ {\frac{5}{6}(p + 1)(p + 2)(p + 3)} \right]}^{2}}.$

Для задач размером до ${\text{NE}} = {{10}^{5}}$ элементов с пространственной точностью четвертого порядка полностью неявный метод требует $45$ Гб памяти только для хранения матрицы Якоби. Для сравнения eMG требует только $0.03$ Гб памяти для хранения векторов решений плюс матрицу якобиана первого порядка при том же размере задачи. Таким образом, метод eMG является более экономичным, чем полностью неявный метод для моделирования высокого порядка с ограниченными вычислительными ресурсами.

Фиг. 6.

$p$-Независимая сходимость метода eMG.

Фиг. 7.

История сходимости неявного метода с различной пространственной точностью.

Фиг. 8.

Сравнение производительности метода eMG и неявного метода.

5.3. Сверхзвуковое сопло Лаваля

Для того чтобы проверить надежность метода eMG для сквозного счета ударной волны, было рассчитано сверхзвуковое сопло Лаваля при скорости 3 Маха. Геометрия сопла симметрична относительно оси $y$, как показано на фиг. 9. Верхний и нижний профили даны аналитически:

(36)

$y = \pm \left( { - \frac{{{{x}^{3}}}}{{405}} + \frac{{{{x}^{2}}}}{{18}} - \frac{x}{3} + 1} \right),\quad x \in [0,\;8].$

Фиг. 9.

Сверхзвуковое сопло Лаваля: конфигурация сопла.

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

Для этого случая нас интересуют поведение сходимости в стационарном режиме и устойчивость решателя. Мы обнаружили, что схема eMG может повысить устойчивость решателя при расчете ударных волн средней силы, т.е. дозвуковых и сверхзвуковых полей потока. Для контроля устойчивости вычислений в качестве переменной для отслеживания выбрано глобальное максимальное значение числа Маха, исходя из того факта, что решатель HA3D обычно не разваливается, если максимальное число Маха не претерпевает значительных резких изменений. В частности, история сходимости на фиг. 10а,в показывает, что eMG с разрывным методом Галеркина с $p = 1$ столкнулся с низким риском неустойчивости (небольшое повышение числа Маха) только на начальном этапе вычислений, где итерационное решение, стартовавшее из начального однородного потока, подверглось сильному нелинейному перестроению за счет взаимодействия с границами. После этого невязка устойчиво сходится к устойчивому состоянию, как показано на фиг. 10б, г. Контуры давления для решений второго и третьего порядка показаны на фиг. 11, где хорошо видны неограниченные ударные волны с резкими профилями. Однако на фиг. 10а, в показано, что решение $p = 2$ страдает от сильной неустойчивости, так что максимальное число Маха достигало совершенно нефизического значения $12$. К счастью, eMG все еще поддерживает стабильность решателя без развала счета.

Фиг. 10.

История сходимости для сверхзвукового сопла Лаваля при числе Маха 3: глобальный максимум числа Маха (а), (в), невязка по плотности (б), (г); $p = 1$ (а), (б), $p = 2$ (в), (г).

Фиг. 11.

Контуры давления для сверхзвукового сопла Лаваля при числе Маха 3: $p = 1$ (а), $p = 2$ (б).

Для решения $p = 2$ сходимость по невязке демонстрирует периодические осцилляции, следовательно, поле течения переходит в нестационарное периодическое состояние. Нестационарность, вероятно, вызвана переносом возмущений, которые происходят от колебаний разрывов. Затем колебания переносятся вниз по течению и проходят весь путь до выхода из сопла, как показано на фиг. 11б. При ближайшем рассмотрении фиг. 12 видно, что основная ударная структура была хорошо разрешена, и ее ширина находится в пределах одного элемента. Решения довольно резкие, без эффектов сглаживания разрывов. Более того, ударные волны решения $p = 2$ почти в два раза уже, чем у решения $p = 1$.

Фиг. 12.

Крупный план контуров давления и разрешение сетки: $p = 1$ (а), $p = 2$ (б). Наблюдается подсеточное разрешение для нескольких отражающихся скачков.

Эти результаты показывают, что схема eMG действительно может стабилизировать сквозной счет, в то время как использование традиционных методов без ограничения разрывов обычно приводит к расходящимся решениям (см., например, [20]). Из результатов, показанных на фиг. 13, следует, что профили ударной волны высокого порядка, выделенные вдоль оси $x$, находятся в пределах ширины одной ячейки, в то время как решение первого порядка ($p = 0$) не захватывает никаких разрывов при таком разрешении сетки. Отметим, что такая аппроксимация разрывов высокого порядка не противоречит теореме Годунова. Эта теорема лишь утверждает, что линейные численные схемы решения дифференциальных уравнений, обладающие свойством не порождать новых экстремумов, могут быть не более чем первого порядка точности. Поэтому нелинейные схемы могут создавать профили скачков высокого порядка, как, например, схема WENO. Хотя профили скачков высокого порядка являются резкими, решения демонстрируют колебания вблизи разрывов, испытывая так называемый феномен Гиббса. Использование ограничения или добавление искусственной диссипации для подавления осцилляций также заслуживает будущего исследования, несмотря на то что они все еще являются открытыми проблемами для специалистов по сквозному счету высокого порядка.

Фиг. 13.

Распределения числа Маха, полученные вдоль оси $x$ сопла Лаваля: y = 0 (а), y = 1/5 (б).

5.4. Трансзвуковое крыло ONREA M6

Крыло ONERA M6 рассматривается для оценки надежности сквозного счета в 3D. Эта конфигурация была использована в качестве эталонного случая (см. [21]). Решения по дозвуковому обтеканию были рассчитаны при числе Маха $0.8394$ и угле атаки $\alpha = 3.06^\circ $. Стационарные решения получены с использованием схемы по времени eMG и пространственной дискретизации разрывным методом Галеркина с $p = 1$.

На фиг. 14 показано установившееся течение, полученное на неструктурированной криволинейной сетке. В решении четко видна ударная волна $\lambda $, образованная двумя внутренними ударными волнами, которые сливаются вместе, образуя одну сильную ударную волну в области законцовки крыла. Подсеточное разрешение ударной волны можно наблюдать в центральной области верхней поверхности, где ячейки на поверхности очень крупные. Коэффициенты давления на поверхности крыла сравниваются с экспериментальными результатами в четырех различных пролетных точках вдоль крыла, а именно, $\eta = 0.44$, $0.65$, $0.90$, $0.99$ (фиг. 15). Видно, что при всех расположениях крыла по пролету согласие вычисленного решения с экспериментальными данными достаточно хорошее. История сходимости полученного решения, которая монотонно сходится к устойчивому состоянию, показана на фиг. 16б. Максимальное число Маха также не подвержено резким изменениям, как показано на фиг. 16а. Таким образом, 3D-расчет, полученный с помощью схемы eMG, является эффективным и надежным.

Фиг. 14.

Контур давления и разрешение сетки: можно наблюдать подсеточное разрешение $\lambda $ ударной волны.

Фиг. 15.

Крыло ONREA M6: сравнение с экспериментальными данными в четырех точках пролета.

Фиг. 16.

История сходимости для крыла ONREA M6: глобальное максимальное число Маха (а), невязка плотности (б).

6. ВЫВОДЫ

Схема экспоненциального интегрирования по времени первого порядка EXP1 была использована для повышения реализуемости решателя на основе разрывных конечных элементов Галеркина высокого порядка HA3D для расчета стационарных сжимаемых потоков. Метод eMG был разработан для сбалансированного рассмотрения эффективности, устойчивости и использования памяти. Обсуждены алгоритмы и физическая природа методов. Представлены результаты численных расчетов как гладких, так и разрывных течений. Для вычислений гладких течений метод eMG показывает $p$-независимую скорость сходимости для порядка аппроксимации $p = 1{\kern 1pt} - {\kern 1pt} 3$, и он в 20 раз быстрее, чем ILU-GMRES при расчете трехмерного течения около сферы. По сравнению с полностью неявным методом, метод eMG использует меньше памяти и достигает значительной вычислительной эффективности.

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

Список литературы

Helenbrook B., Mavriplis D., Atkins H. Analysis of p-multigrid for continuous and discontinuous finite element discretizations. 16th AIAA Comput. Fluid Dynam. Conf. AIAA-2003-3989, 2003.
Helenbrook B.T., Atkins H.L. Application of p-multigrid to discontinuous Galerkin formulations of the Poisson equation // AIAA J. 2006. V. 44. № 3. P. 566–575.
Fidkowski K.J., Oliver T.A., Lu J., et al. p-Multigrid solution of high-order discontinuous Galerkin discretizations of the compressible Navier-Stokes equations // J. Comput. Phys. 2005. V. 207. № 1. P. 92–113.
Atkins H., Helenbrook B. Numerical evaluation of p-multigrid method for the solution of discontinuous Galerkin discretizations of diffusive equations. 17th AIAA Comput. Fluid Dynam. Conf. AIAA-2005-5110, 2005.
Bassi F., Franchina N., Ghidoni A., et al. Spectral p-multigrid discontinuous Galerkin solution of the Navier–Stokes equations // Inter. J. Numer. Meth. Fluid. 2011. V. 67. № 11. P. 1540–1558.
Li S.-J., Wang Z.J., Ju L., Luo L.-S. Explicit large time stepping with a second-order exponential time integrator scheme for unsteady and steady flows. 55th AIAA Aerospace Sci. Meet. AIAA-2017-0753, 2017.
Li S.-J., Wang Z.J., Ju L., Luo L.-S. Fast time integration of Navier–Stokes equations with an exponential-integrator scheme. 2018 AIAA Aerospace Sci. Meet. AIAA-2018-0369, 2018.
Li S.-J., Luo L.-S., Wang Z.J., Ju L. An exponential time-integrator scheme for steady and unsteady inviscid flows // J. Comput. Phys. 2018. V. 365. P. 206–225.
Li S.-J. Mesh curving and refinement based on cubic Bézier surface for high-order discontinuous Galerkin methods // Comput. Math. and Math. Phys. 2019. V. 59. №. P. 2080–2092.
Li S.-J., Ju L., Hang S. Adaptive exponential time integration of the Navier–Stokes equations. AIAA-2020-2033, 2020.
Caliari M., Ostermann A. Implementation of exponential Rosenbrock-type integrators // Appl. Numer. Math. 2009. V. 59. № 3. P. 568–582.
Tokman M., Loffeld J. Efficient design of exponential-Krylov integrators for large scale computing // Procedia Comput. Sci. 2010. V. 1. № 1. P. 229–237.
Saad Y. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems // SIAM J. Sci. Stat. Comput. 1986. V. 7. № 3. P. 856–869.
Saad Y. Analysis of some Krylov subspace approximations to the matrix exponential operator // SIAM J. Numer. Anal. 1992. V. 29. № 1. P. 209–228.
Moler C. Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later // SIAM J. Numer. Anal. 2003. V. 45. № 1. P. 3–49.
Roe P. Approximate Riemann solvers, parameter vectors, and difference schemes // J. Comput. Phys. 1981. V. 43. № 2. P. 357–372.
Kitamura K., Shima E. Towards shock-stable and accurate hypersonic heating computations: A new pressure flux for AUSM-family schemes // J. Comput. Phys. 2013. V. 245. P. 62–83.
Pierce N., Giles M. Preconditioning compressible flow calculations on stretched meshes. 34th Aerospace Sci. Meet. and Exhibit. AIAA-1996-889, 1996.
Toro E. Riemann solvers and numerical methods for fluid dynamics. Springer, 1999.
Persson P.-O., Peraire J. Sub-cell shock capturing for discontinuous Galerkin methods. 44th AIAA Aerospace Sci. Meet. and Exhibit. AIAA-2006-112, 2006.
Schmitt V., Charpin F. Pressure distributions on the ONERA-M6-Wing at transonic mach numbers. Experimental data base for computer program assessment. Rep. Fluid Dynam. Panel Work. Group 04, AGARD-AR 138, May 1979.

Дополнительные материалы отсутствуют.

Инструменты

предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал