Журнал вычислительной математики и математической физики, 2019, T. 59, № 6, стр. 913-919

ВВЕДЕНИЕ

Обширный набор методов разного порядка аппроксимации, применяющийся в настоящее время для численного решения задач CFD, представляет дополнительные возможности для оценки погрешности расчета. Здесь мы рассмотрим определение окрестности численного решения, в котором находится истинное решение (мы будем обозначать это как “захват истинного решения”), основанное на использовании ансамбля расчетов, полученных схемами разного порядка аппроксимации. Обычно методы оценки погрешности расчета основаны на некоторых нормах. В данной работе основное внимание уделено оценке расстояния между решениями с помощью метрики, которая может быть не связана с нормой. На языке норм соответствующий подход представлен в работах [1]–[3].

1. ЗАХВАТ ИСТИННОГО РЕШЕНИЯ НА АНСАМБЛЕ ЧИСЛЕННЫХ РЕШЕНИЙ

Рассмотрим ансамбль численных решений, полученных схемами разного порядка аппроксимации на одной и той же сетке. Обозначим численное решение как вектор ${{u}^{{(i)}}} \in {{R}^{N}}$ ($i$ – номер схемы, $N$ – число переменных), а значения неизвестного истинного решения в узлах сетки как вектор $\tilde {u} \in {{R}^{N}}$. Используется некоторая (произвольная) метрика $d(u,\text{v})$ в пространстве решений. Неизвестное отклонение истинного решения в узлах сетки от численного решения (величина ошибки расчета) обозначено как $d({{u}^{{(k)}}},\tilde {u}) = {{\delta }_{{0,k}}}$ (например, $d({{u}^{{(k)}}},\tilde {u}) = {{\left\| {{{u}^{{(k)}}} - \tilde {u}} \right\|}_{{{{L}_{2}}}}}$).

Следующая теорема может быть сформулирована для численных решений ${{u}^{{(1)}}}$ и ${{u}^{{(2)}}}$, для которых априори известно соотношение величин ошибок ${{\delta }_{{0,1}}} \geqslant 2{{\delta }_{{0,2}}}$.

Теорема 1. Пусть расстояние ${{\delta }_{{1,2}}}$ между двумя численными решениями ${{u}^{{(1)}}} \in {{R}^{N}}$ и ${{u}^{{(2)}}} \in {{R}^{N}}$ известно из расчета, а также априори известно соотношение величин ошибок (расстояний между точным и численным решением)

тогда точное решение находится внутри гиперсферы радиуса ${{\delta }_{{1,2}}}$ с центром в более точном решении ${{u}^{{(2)}}}$:

Доказательство. Неравенство треугольника [4] для расстояний ${{\delta }_{{0,1}}}$, ${{\delta }_{{1,2}}}$, ${{\delta }_{{0,2}}}$ между точками ${{u}^{{(1)}}}$, ${{u}^{{(2)}}}$ и $\tilde {u}$ может быть представлено как ${{\delta }_{{0,1}}} \leqslant {{\delta }_{{1,2}}} + {{\delta }_{{0,2}}}$ и переформулировано как ${{\delta }_{{0,1}}} - {{\delta }_{{0,2}}} \leqslant {{\delta }_{{12}}}$. Учитывая (1) в форме ${{\delta }_{{0,1}}} - {{\delta }_{{0,2}}} \geqslant {{\delta }_{{0,2}}}$, получим ${{\delta }_{{0,2}}} \leqslant {{\delta }_{{0,1}}} - {{\delta }_{{0,2}}} \leqslant {{\delta }_{{12}}}$ и желаемое выражение ${{\delta }_{{0,2}}} \leqslant {{\delta }_{{1,2}}}$.

2. АНАЛИЗ СООТНОШЕНИЯ ПОГРЕШНОСТЕЙ НА АНСАМБЛЕ ЧИСЛЕННЫХ РАСЧЕТОВ

Несмотря на естественное мнение о том, что схемы более высокого порядка аппроксимации более точны, очевидной слабостью Теоремы 1 является предположение о существовании решений с упорядоченной по величине погрешностью. Поэтому целесообразно рассмотреть возможности проверки данной упорядоченности в численных экспериментах. Простейший анализ показывает (а расчеты подтверждают), что в некоторых случаях набор расстояний между решениями ${{\delta }_{{i,j}}}$ позволяет выделить решения, расположенные поблизости друг от друга, и удаленные решения. Например, если одно из решений ${{u}^{{(1)}}}$ существенно менее точно, чем остальные (${{\delta }_{{0,1}}} \gg {{\delta }_{{0,i}}}$), то набор расстояний ${{\delta }_{{i,j}}}$ расщепляется на

кластер “неточных” решений с большой величиной ${{\delta }_{{1,j}}}$ и

кластер “точных” решений ${{\delta }_{{i,j}}}$ ($i \ne 1$).

Это вызвано асимптотиками ${{\delta }_{{1,j}}}{\text{/}}{{\delta }_{{0,1}}} \to 1$ и ${{\delta }_{{i,j}}}\;(i \ne 1){\text{/}}{{\delta }_{{0,1}}}\sim ({{\delta }_{{0,i}}} + {{\delta }_{{0,j}}}){\text{/}}{{\delta }_{{0,1}}} \to 0$ при ${{\delta }_{{0,i}}}{\text{/}}{{\delta }_{{0,1}}} \to 0$.

Расслоение набора расстояний между приближенными решениями на кластеры свидетельствует об упорядочении погрешностей. Сравним набор расстояний ${{\delta }_{{1,j}}}$ и ${{\delta }_{{i,j}}}$, где ${{u}^{{(1)}}}$ – решение с максимальной погрешностью, ${{u}^{{(i)}}}$ – одно из более точных решений (мы будем искать истинное решение в его окрестности), обозначим $\delta _{{0,i}}^{{\max }}$ максимальную ошибку в подмножестве более точных решений. Максимальное значение из ${{\delta }_{{i,j}}}$ ($i \ne 1$) (расстояние от нуля до максимального значения в кластере “точных” решений) обозначим как верхнюю границу первого кластера ${{d}_{1}}$ (набора расстояний между “точными” решениями), минимальное из ${{\delta }_{{1,j}}}$ обозначим как нижнюю границу второго кластера ${{d}_{2}}$ (расстояний между “точными” решениями и наиболее неточным).

Можно сформулировать следующее (полу)эвристическое правило, обеспечивающее применимость Теоремы 1.

Критерий 1. Если расстояние между кластерами больше размера кластера точных решений ${{d}_{2}} - {{d}_{1}} > {{d}_{1}}$, тогда истинное решение расположено внутри гиперсферы радиуса ${{\delta }_{{1,i}}}$ с центром в ${{u}^{{(i)}}}$: ${{\delta }_{{0,i}}} \leqslant {{\delta }_{{1,i}}}$, где ${{u}^{{(i)}}}$ принадлежит кластеру более точных решений, а ${{u}^{{(1)}}}$ – максимально неточное решение.

Это предположение основано на допущении о том, что размер кластера точных решений (расстояние от нуля до границы кластера) равен удвоенной максимальной ошибке в этом кластере ${{d}_{1}} = 2\delta _{{0,i}}^{{\max }}$, а кластер неточных решений принадлежит интервалу $({{\delta }_{{0,1}}} - \delta _{{0,i}}^{{\max }},{{\delta }_{{0,1}}} + \delta _{{0,i}}^{{\max }})$ и ${{d}_{2}} = {{\delta }_{{0,1}}} - \delta _{{0,i}}^{{\max }}$. Обе эти оценки завышены и соответствуют коллинеарным векторам ошибки. Если принять их верными, соотношение размера точного кластера и расстояния между кластерами (${{d}_{2}} - {{d}_{1}} > {{d}_{1}}$, ${{d}_{2}} > 2{{d}_{1}}$) будет ${{\delta }_{{0,1}}} - \delta _{{0,i}}^{{\max }} > 4\delta _{{0,i}}^{{\max }}$. Это ведет к соотношению ${{\delta }_{{0,1}}} > 5\delta _{{0,i}}^{{\max }}$, соответственно справедливо условие (1) Теоремы 1 ${{\delta }_{{0,1}}} > 2\delta _{{0,i}}^{{\max }}$.

Этот критерий не является строгим, тем не менее, численные эксперименты для двумерных невязких сверхзвуковых течений указывают на его практическую применимость. Возможно, это связано с наличием некоторого “запаса” (${{\delta }_{{0,1}}} > 5\delta _{{0,i}}^{{\max }}$).

Теорема 1 и критерий 1 сформулированы для расстояний, заданных некоторой метрикой. Представляет интерес, какая из метрик лучше подходит для оценки погрешности расчета. Мы сравним метрики, порожденные нормами ${{L}_{1}}$ и ${{L}_{2}}$, а также метрику, имитирующую относительную погрешность (далее ее обозначаем как REM) и метрику IMED [5].

Рассмотрим двумерные уравнения Эйлера, имеющие четыре компоненты ${{u}^{{(i)}}} = \{ {{\rho }^{{(i)}}},{{U}^{{(i)}}},{{V}^{{(i)}}},{{e}^{{(i)}}}\} $. Для метрик, порожденных нормой ${{L}_{2}}$, расстояние между решениями имеет вид

Для проблем вычислительной газодинамики вектор решения содержит элементы с различным физическим смыслом, такие как плотность, компоненты скорости. Их прямое сопоставление, строго говоря, может быть проблематичным. Поэтому параллельно выражению (3) расстояние между решениями рассчитывалось с использованием нормализованного выражения

которое имитирует относительную погрешность (REM). Выражение (4) соответствует метрике Махаланобиса [6]. Это расстояние определено метрическим тензором (глобальным) с диагональной матрицей ${{M}_{{j,k}}}$. С учетом представления $M = A{\text{*}}A$ (верного для симметричной положительно определенной матрицы), можно выразить $(\Delta u_{{}}^{{(i)}},M\Delta u_{{}}^{{(i)}}) = {{(\Delta u_{{}}^{{(i)}},A{\text{*}}A\Delta u_{{}}^{{(i)}})}^{{1/2}}}$ = ${{(A\Delta u_{{}}^{{(i)}},A\Delta u_{{}}^{{(i)}})}^{{1/2}}} = {{(\Delta z_{{}}^{{(i)}},\Delta z_{{}}^{{(i)}})}^{{1/2}}}$, что соответствует норме ${{L}_{2}}$ в трансформированном пространстве $Au_{{}}^{{(i)}}$.

Вышеописанные метрики чувствительны к малым вариациям поля течения таким, как, например, сдвиг положения ударной волны на одну ячейку. При этом два численных решения, порожденных таким сдвигом, воспринимаются как удаленные и описывающие разные течения. Однако с практической точки зрения такие течения идентичны. Поэтому эти расстояния не позволяют оценить структурную близость между решениями. Эвклидова метрика, модифицированная для анализа изображений (IMage Euclidean Distance (IMED)), более удобна с этой точки зрения [5]. Она может быть описана метрическим тензором

Величина $\left| {{{P}_{i}} - {{P}_{j}}} \right|$ – это расстояние между узлами ${{P}_{i}}$ и ${{P}_{j}}$ на сетке. Для нашей двумерной задачи, расстояние между решениями оценивалось с помощью выражения (представленного здесь только для плотности)

Это расстояние соответствует усредненной по пространству погрешности. IMED тоже является метрикой типа Махалонобиса, в некотором базисе сводящейся к ${{L}_{2}}$. Асимптотически IMED стремится к ${{L}_{2}}$ при $\sigma \to 0$.

Метрики, имеющие физический смысл, иллюстративную мощь и потенциал для сравнения течений, не исчерпываются указанными выше, поэтому поиск оптимальной метрики представляет существенный интерес.

3. ТЕСТОВЫЕ ЗАДАЧИ

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

Здесь предполагается суммирование по повторяющимся индексам, $i,k = 1,2$, ${{U}^{1}} = U$, ${{U}^{2}} = V$ компоненты скорости, ${{h}_{0}} = ({{U}^{2}} + {{V}^{2}}){\text{/}}2 + h$, $h = \frac{\gamma }{{\gamma - 1}}\frac{P}{\rho } = \gamma e$, $e = \frac{{RT}}{{\gamma - 1}}$, $E = \left( {e + \frac{1}{2}({{U}^{2}} + {{V}^{2}})} \right)$ энтальпии и энергии, $P = \rho RT$ – уравнение состояния, $\gamma = {{C}_{p}}{\text{/}}{{C}_{\text{v}}}$.

В качестве тестовых примеров рассмотрены: косой скачок уплотнения, взаимодействие ударных волн I типа по классификации Edney, взаимодействие ударных волн VI типа по классификации Edney [7].

Расчеты проводились при ${{C}_{p}}{\text{/}}{{C}_{\text{v}}} = 1.4$, в диапазоне чисел Маха $M = 3{\kern 1pt} - {\kern 1pt} 5$ и углов отклонения потока $\alpha = 10{\kern 1pt} - {\kern 1pt} 30^\circ $. Некоторые типичные результаты представлены ниже для иллюстрации. На фиг. 1 представлено распределение плотности для конфигурации скачков VI типа по Edney ($M = 4$, два последовательных отклонения потока на углы ${{\alpha }_{1}} = 10^\circ $ и ${{\alpha }_{2}} = 15^\circ $). Поле течения определяется сливающимися скачками уплотнения, линией контактного разрыва и веером волн разрежения. Фиг. 2 представляет изолинии плотности для конфигурации I типа по Edney ($M = 3$, углы отклонения потока ${{\alpha }_{1}} = 20^\circ $ и ${{\alpha }_{2}} = 15^\circ $). Основными элементами являются пересекающиеся скачки и линия контактного разрыва.

Для этих задач достаточно просто построить аналитические решения, проекции которых на расчетную сетку мы будем рассматривать далее как “истинные” решения. Это позволяет тестировать работоспособность предлагаемого метода оценки погрешности и захвата решения. Если истинное решение $\tilde {u} \in {{R}^{N}}$ находится в гиперсфере $d(\tilde {u},{{u}^{{(i)}}}) \leqslant d({{u}^{{(1)}}},{{u}^{{(i)}}})$, мы считаем, что метод продемонстрировал работоспособность. В расчетах мы сравним метрики, порожденные нормами ${{L}_{2}}$ и ${{L}_{1}}$, метрики типа Махалонобиса [6], одна из которых описывает нормированную ошибку, другая [5] обеспечивает некоторое размывание ошибки.

4. РЕЗУЛЬТАТЫ ЧИСЛЕННЫХ ЭКСПЕРИМЕНТОВ

В работе анализировался ансамбль расчетов поля течения, полученных разными методами. Использовалась схема 1 порядка точности типа CIR (Courant Isaacson Rees) [8] в варианте [9], обозначенная как $S1$, схема второго порядка точности, основанная на методе MUSCL [10], использующая на границах AUFS [11] и обозначенная как $S2$, схема 3 порядка точности – модификация метода Chakravarthy-Osher $S3$ [12, 13] и схема 4 порядка точности [14], обозначенная как $S4$. Расчеты проводились на равномерных сетках, содержащих $100 \times 100$, $200 \times 200$ и $400 \times 400$ пространственных узлов.

Методы $S1$, $S2$, $S3$, $S4$ (1, 2, 3 и 4 номинальных (заявленных разработчиками) порядков аппроксимации) демонстрировали порядок сходимости чуть ниже $n = 1{\text{/}}2$ в норме ${{L}_{2}}$. В норме ${{L}_{1}}$ те же самые расчеты продемонстрировали порядок сходимости чуть выше $n = 1{\text{/}}2$.

После получения набора численных решений мы проверяем выполнение Критерия 1 и, далее, проверяем захват истинного решения. Захват считаем успешным, если оценка погрешности больше истинной погрешности, полученной при сравнении с аналитическим решением $d(\tilde {u},{{u}^{{(i)}}}) \leqslant d({{u}^{{(1)}}},{{u}^{{(i)}}})$. Например, в ${{L}_{2}}$ ${{\left\| {{{u}^{{(i)}}} - {{u}^{{(k)}}}} \right\|}_{{{{L}_{2}}}}}$ больше истинной погрешности ${{\left\| {{{u}^{{(k)}}} - \tilde {u}} \right\|}_{{{{L}_{2}}}}}$.

Численные эксперименты позволяют заключить, что схема $S1$ (“неточная”) и схемы $S2$, $S3$, $S4$ (“точные”) позволяют найти окрестность численного решения, содержащую истинное решение. Сравнение схем $S2$, $S3$, $S4$ не позволяет выделить кластеры и захватить точное решение. Погрешности этих схем близки по величине и расщепления на кластеры не происходит.

Численные эксперименты, проведенные для течений, соответствующих уединенному косому скачку уплотнения, взаимодействию скачков I и VI типов по Edney демонстрируют возможность захвата истинного решения, если Критерий 1 удовлетворен.

Например, фиг. 3 показывает, что набор расстояний между численными решениями ${{\left\| {d{{u}_{{i,k}}}} \right\|}_{{{{L}_{2}}}}} = {{\left\| {{{u}^{{(i)}}} - {{u}^{{(k)}}}} \right\|}_{{{{L}_{2}}}}}$ распадается на два кластера для схем взаимодействия скачков VI типа по Edney (с “неточной” схемой $S1$, сетка $100 \times 100$). Это позволяет захватить точное решение, что продемонстрировано на фиг. 4. Нужно обратить внимание, что данные достаточно громоздки и для удобства их представления на фиг. 3, 4 по обеим осям расположена норма погрешности, хотя анализируемые данные одномерны.

Численные эксперименты показали, что метрика, порожденная ${{L}_{1}}$, обеспечивала успешный захват во всех случаях, когда выполняется Критерий 1. Захват точного решения в норме ${{L}_{2}}$ получается не всегда. Однако нарушение критерия захвата не очень велико. Максимальная относительная величина нарушения условия захвата, наблюдавшаяся в расчетах, около 0.1. Метрика REM (4) дает некоторое улучшение в сравнении с ${{L}_{2}}$. Метрика IMED [5] дает с точки зрения успешности захвата истинного решения результаты, сходные с ${{L}_{1}}$.

В качестве примера в табл. 1 представлены результаты оценок погрешности расчета $d({{u}^{{(i)}}},\tilde {u})$ и ${{\delta }_{{1,i}}}$ для различных метрик с точки зрения захвата истинного решения для теста Edney-I ($M = 3$, углы отклонения потока ${{\alpha }_{1}} = 20^\circ $ и ${{\alpha }_{2}} = 15^\circ $, сетка $400 \times 400$).

5. ОБСУЖДЕНИЕ

Повсеместно используемый контроль сходимости по сетке основывается на правиле Рунге [15]. По этому правилу, если разность двух приближенных решений на грубой сетке ${{T}_{h}}$ и на мелкой сетке ${{T}_{{h,ref}}}$ с шагом сетки ${{h}_{{ref}}}$ мала, тогда ${{u}_{{h,ref}}}$ и ${{u}_{h}}$, вероятно, близки к точному решению. С точки зрения практики хотелось бы иметь оценку типа $\left\| {{{u}_{h}} - \tilde {u}} \right\| \leqslant \delta $ с вычислимой $\delta $. Достаточно близко к этому идеалу подходит метод Ричардсона [16], который позволяет определить уточненное решение и оценку погрешности в том случае, если во всем поле существует единый порядок точности. Для этого используется набор решений на последовательно измельчаемых сетках. К сожалению, во многих задачах вычислительной аэрогазодинамики порядок ошибки на разных элементах течения различен, что существенно затрудняет или исключает возможность применения метода Ричардсона.

В данной работе мы рассматриваем некоторую альтернативу методу Ричардсона и правилу Рунге. Она основана на ансамбле расчетов, выполненных на одной и той же сетке с помощью численных методов разного порядка аппроксимации. Величина погрешности $\delta ({{u}_{h}},\tilde {u}) \leqslant \delta $ оценивается в некоторой метрике (или в норме $\left\| {{{u}_{h}} - \tilde {u}} \right\| \leqslant \delta $). Расчеты можно прекращать при достижении заданной точности. Рассматриваемый метод может быть использован в форме постпроцессора аналогично экстраполяции Ричардсона. Однако он не нуждается в измельчении сетки и может применяться вне асимптотического диапазона. Полуэмпирическая природа Критерия 1 является основным недостатком рассматриваемого метода. Поэтому он не может заменить стандартный подход, основанный на мельчении сетки, однако может дополнить его вычислительно дешевым односеточным алгоритмом.

ЗАКЛЮЧЕНИЕ

Набор решений, полученных различными солверами на одной и той же сетке, содержит информацию о погрешности расчета и расположении истинного решения.

Если набор расстояний между численными решениями распадается на кластеры, соответствующие “точным” и “неточным” схемам, а расстояние между кластерами больше размера кластера “точных” решений (Критерий 1), то истинное решение расположено в гиперсфере с центром в более точном решении и с радиусом, равным расстоянию между “точным” и “неточным” решениями.

Вычислительные эксперименты подтвердили применимость данного эвристического правила для двумерных сверхзвуковых течений, описываемых уравнениями Эйлера и содержащих скачки уплотнения. Успех в значительной степени зависит от выбора метрики.

Метрика, порожденная нормой ${{L}_{1}}$, обеспечивает оценку погрешности и захват истинного решения во всех случаях, когда выполнен Критерий 1.

Метрика, порожденная ${{L}_{2}}$, обеспечивает оценку погрешности не во всех случаях, хотя наблюдавшиеся нарушения не велики.

Метрика, имитирующая относительную погрешность, позволяет захват истинного решения в некоторых случаях, когда захват на основе ${{L}_{2}}$ не удается.

Метрика IMED [5] демонстрирует характеристики, сходные с ${{L}_{1}}$.