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

Abstract

Рассматривается новая тестовая задача для одномерных уравнений газовой динамики. Начальные данные в этой задаче представляют собой периодическую гладкую волну. За конечное время в течении газа формируются ударные волны. Исследуется сеточная сходимость двух линейных схем третьего порядка аппроксимации: бикомпактной схемы и схемы Русанова. Демонстрируется, что в области влияния ударной волны обе схемы имеют лишь первый порядок интегральной сходимости. В то же время при расчете уравнений изоэнтропической газовой динамики выбранные схемы сходятся не менее чем со вторым порядком. Библ. 38. Фиг. 6.

Full Text

Введение

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

Одной из основных характеристик численной схемы является порядок сходимости, или порядок точности. Высокий порядок точности (второй и выше) позволяет добиться заданной погрешности искомого решения на менее подробной сетке, т.е. с меньшими затратами машинного времени и памяти (см. [1]). Это особенно актуально при решении многомерных нестационарных задач.

Другой важной характеристикой численной схемы является наличие или отсутствие свойства монотонности (см. [2]). В монотонной схеме исключены ложные осцилляции численного решения в окрестностях сильных разрывов точного решения. Осцилляции такого рода вредны не только тем, что, собственно, искажают физику моделируемого процесса, но и тем, что могут привести к потере устойчивости счета и его аварийному завершению.

В известной работе С. К. Годунова [3] была доказана следующая теорема: для линейного уравнения переноса не существует линейных монотонных схем с порядком аппроксимации выше первого. В монографии [4] была высказана любопытная трактовка этой теоремы: из трех свойств схемы для линейного уравнения переноса — высокого порядка аппроксимации, монотонности и линейности — одновременно иметь место могут только два.

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

Наиболее очевидным, если не единственным, способом получить монотонную схему высокого порядка точности является применение нелинейных аппроксимаций, зависящих от самого численного решения. За последние полвека было предложено множество методов их построения: ограничители наклонов/потоков (см. [6–8]), реконструкции ENO/WENO (см. [9–11]), численные фильтры (см. [12–14]), адаптивные искусственные вязкости (см. [15–17]) и т.д. В их основе лежит одна и та же идея: автоматически (в рамках подхода сквозного счета) выделить окрестности сильных разрывов решения и в этих окрестностях ввести регуляризующие поправки в аппроксимацию формально высокого порядка. Общепринятая гипотеза состоит в том, что порядок точности в областях гладкости точного решения должен оставаться высоким, так как коррекция аппроксимации проводится локально, в малых окрестностях сильных разрывов (где нет и не ожидается классической аппроксимации высокого порядка).

Эта гипотеза была опровергнута в работах [18–28]. На контрпримерах было показано, что все нелинейные монотонные схемы высокого порядка аппроксимации имеют в действительности не более чем первый порядок точности в областях влияния ударных волн. Это явление объясняется так. Невязка имеет высокий порядок по шагу сетки благодаря высокочастотным гармоникам численного решения. Но именно эти гармоники искажаются либо подавляются вблизи фронтов ударных волн (ради монотонности схемы), причем это происходит непрерывно, при каждом переходе со слоя на слой.

Чтобы решить данную проблему, в работах [29–32] был разработан принципиально новый класс схем, названных комбинированными. В общих чертах идея комбинированных схем из [29–31] сводится к следующему. Исходная начально-краевая задача решается целиком (во всей пространственно-временной расчетной области) с помощью линейной схемы высокого порядка аппроксимации, или базисной схемы. Решение базисной схемы используется для постановки внутренних начально-краевых задач в (подвижных) окрестностях ударных волн. Для решения этих внутренних задач применяется нелинейная монотонная схема, или внутренняя схема. В интересующий момент времени решение внутренней схемы используется, чтобы заменить (скорректировать) решение базисной схемы вблизи ударных волн. Поскольку эта процедура имеет характер постобработки, высокий порядок точности в областях влияния ударных волн должен сохраняться.

Теория комбинированных схем опирается на гипотезу, что линейные схемы высокого порядка аппроксимации (базисные схемы) могут иметь высокий порядок точности в областях влияния ударных волн. Результаты работ [22–27, 29–31] как будто бы свидетельствуют в пользу этого. Однако нельзя не отметить, что в перечисленных работах рассматривалась исключительно система уравнений мелкой воды.

Цель настоящей работы — показать, что эта гипотеза неверна в общем случае; что даже для линейных схем высокого порядка аппроксимации, вообще говоря, гарантирован не более чем первый порядок точности в областях влияния ударных волн. Для этого мы рассмотрим более сложную гиперболическую систему уравнений — систему уравнений газовой динамики. Данный негативный результат создает мотивацию для переосмысления понятия высокоточной схемы и теории комбинированных схем.

Работа структурирована следующим образом. В разд. 1 строится тестовая задача Коши для одномерных уравнений Эйлера. В разд. 2 приводятся численные схемы. В разд. 3 описываются критерии оценки их точности. В разд. 4, 5 обсуждаются результаты расчетов.

1. Тестовая задача для уравнений газовой динамики

Рассмотрим одномерные уравнения Эйлера для совершенного газа:

${\partial }_{t}u+{\partial }_{x}f\left(u\right)=\mathrm{0,}u=\left[\begin{array}{c}\rho \\ \rho v\\ E\end{array}\right],f\left(u\right)=\left[\begin{array}{c}\rho v\\ \rho v^2+p\\ v(E+p)\end{array}\right],E=\frac{p}{\gamma -1}+\frac{\rho v^2}{2}.$ (1)

Здесь x — пространственная координата $\left({\partial }_x\equiv \partial /\partial x\right),$ t — время, $u=u(x,t)$ — искомый вектор консервативных переменных, f(u) — вектор потоков, r — плотность, v — скорость, p — давление, E — удельная полная энергия, $\gamma =const$ — показатель адиабаты. Выпишем собственные значения матрицы Якоби $A\left(u\right)={\partial }_{u}f\left(u\right)$ (см. [5]):

${\lambda }_1=v-c,{\lambda }_2=v,{\lambda }_3=v+c,$

где $c=c(p,\rho )=\sqrt{\gamma p/\rho }$ — адиабатическая скорость звука. Система (1) всегда имеет как минимум один инвариант Римана — удельную энтропию газа:

$w_2=s=s(p,\rho )=c_v\ln \frac{p}{{\rho }^{\gamma }}+const,$

где c_v = const — удельная теплоемкость при постоянном объеме. Если течение является изоэнтропическим (s = const), то существуют еще два инварианта Римана:

$w_{\mathrm{1,3}}=v\mp \frac{2c}{\gamma -1}.$

(Верно и обратное: если эти функции — инварианты Римана, то течение изоэнтропическое.)

Поставим задачу Коши для уравнений (1). Начальные условия для нее мы зададим по аналогии с задачей Коши для уравнений мелкой воды, сформулированной в работе [23]:

$\begin{array}{c}{w}_1\left(x\mathrm{,0}\right)=-b_0,\\ w_2\left(x\mathrm{,0}\right)=s_0=s(p_0,{\rho }_0),\\ w_3\left(x\mathrm{,0}\right)=2a_0\sin \left(\frac{2\pi x}{X}+\frac{\pi }{4}\right)+b_0,\end{array}$(2)

где p₀, r₀, s₀, X, a₀, b₀ — постоянные параметры. Отличие от начальных условий из [23] состоит в наличии инварианта w₂ = s, который полагается тождественно постоянным (в противном случае полный набор инвариантов Римана не существует глобально и формула (2) теряет смысл). Перепишем начальные условия (2) для инвариантов Римана в виде начальных условий для примитивных переменных (v,r,p):

$\begin{array}{c}v\left(x\mathrm{,0}\right)=a_0\sin \left(\frac{2\pi x}{X}+\frac{\pi }{4}\right),\\ \frac{\rho \left(x\mathrm{,0}\right)}{{\rho }_0}={\left[\frac{\gamma -1}{2c_0}\left(v\left(x\mathrm{,0}\right)+b_0\right)\right]}^{\frac{2}{\gamma -1}},\\ \frac{p\left(x\mathrm{,0}\right)}{p_0}={\left(\frac{\rho \left(x\mathrm{,0}\right)}{{\rho }_0}\right)}^{\gamma },\end{array}$ (3)

где $c_0=\sqrt{\gamma p_0/{\rho }_0}.$ Далее значения параметров будут браться следующими:

$\gamma =\mathrm{1.4,}{p}_0=\mathrm{1,}{\rho }_0=\mathrm{1,}{s}_0=\mathrm{1,}X=\mathrm{10,}{a}_0=\mathrm{1,}{b}_0=10.$ (4)

Начальные данные (3) при параметрах (4) изображены на фиг. 1.

Фиг. 1. Скорость, плотность и давление газа в начальный момент времени (t = 0).

Ясно, что начальные данные (3) представляют собой периодическое гладкое изоэнтропическое течение. С течением времени решение задачи (1), (3) ведет себя как опрокидывающаяся прогрессивная волна. Она движется слева направо (по оси x); те ее участки, где градиенты параметров газа положительны, растягиваются, а те участки, где эти градиенты отрицательны, напротив, сокращаются. Приблизительно в момент времени t = 1.3 происходит градиентная катастрофа и формируются периодические ударные волны, которые затем распространяются с числом Маха, колеблющимся в диапазоне от 1.2 до 1.4. Каждая из ударных волн имеет область влияния, которая примыкает к фронту слева и расширяется.

2. Численные схемы

Рассмотрим две линейные численные схемы высокого порядка аппроксимации для уравнений (1): бикомпактную схему (см. [33; 34]) и схему Русанова (см. [35; 36]). Как уже пояснялось во введении, линейными мы называем такие численные схемы, в которых аппроксимация не зависит от самого численного решения. В разд. 5 мы обсудим, является ли выбор именно этих двух схем принципиальным.

Бикомпактная схема

Поскольку решение задачи Коши (1), (3) является периодическим, мы имеем право сузить вещественную ось до конечного интервала (0, X) и поставить на его границах периодические условия.

Зададим шаги h, t > 0. На отрезке [0, X] построим сетку

$\begin{array}{c}{\omega }_h={{\omega }_h}^{\text{'}}\cup {{\omega }_h}^{\text{'}\text{}\text{'}},\\ {{\omega }_h}^{\text{'}}={\left\{x_j=jh\right\}}_{j=0}^N,\\ {{\omega }_h}^{\text{'}\text{}\text{'}}={\left\{x_{j+1\text{}/2}=(j+1\text{}/2)h\right\}}_{j=0}^{N-1},\end{array}$

а ось времени разобьем на слои $t_n=n\tau$ $\left(n=\mathrm{0,1,}\dots \right).$ Определим на ${\omega }_h$ сеточную функцию $v_i^n$ $\left(i=\mathrm{0,1}\text{}/\mathrm{2,1,}\dots ,N\right),$ аппроксимирующую точное решение $u(x_i,t_n)$ уравнений (1). Пользуясь глобальным расщеплением Лакса–Фридрихса, введем векторы потоков

$f^{\pm }\left(u\right)=\frac{1}{2}f\left(u\right)\pm C_{2}u.$

Параметр $C_2>0$ подбирается так, чтобы матрицы Якоби $A^{\pm }\left(u\right)={\partial }_u{f}^{\pm }\left(u\right)$ были положительно/отрицательно определены на некотором релевантном множестве значений решения. Например, C₂ можно пересчитывать перед каждым переходом со слоя на слой по формуле

$C_2=(0.5+\delta )\underset{i}{max}{\left(\right|v|+c)}_i^n.$

Число $\delta >0$ — это постоянный коэффициент запаса.

Простейшая чисто неявная бикомпактная схема для уравнений (1) имеет вид

$\left\{\begin{array}{l}{A}_0^x\frac{v_{j+1\text{}/2}^{*}-v_{j+1\text{}/2}^n}{\tau }+{\Lambda }_1^x{f}^{+}\text{}\left(v_{j+1\text{}/2}^{*}\right)=\mathrm{0,}\\ {\Lambda }_1^x\frac{v_{j+1\text{}/2}^{*}-v_{j+1\text{}/2}^n}{\tau }+{\Lambda }_2^x{f}^{+}\text{}\left(v_{j+1\text{}/2}^{*}\right)=\mathrm{0,}\end{array}\right.\left\{\begin{array}{l}{A}_0^x\frac{v_{j+1\text{}/2}^{n+1}-v_{j+1\text{}/2}^{*}}{\tau }+{\Lambda }_1^x{f}^{-}\text{}\left(v_{j+1\text{}/2}^{n+1}\right)=\mathrm{0,}\\ {\Lambda }_1^x\frac{v_{j+1\text{}/2}^{n+1}-v_{j+1\text{}/2}^{*}}{\tau }+{\Lambda }_2^x{f}^{-}\text{}\left(v_{j+1\text{}/2}^{n+1}\right)=0.\end{array}\right.$(5)

Сеточные операторы $A_0^x$, ${\Lambda }_1^x$, ${\Lambda }_2^x$ определяются формулами

$\begin{array}{c}{A}_0^x{v}_{j+1\text{}/2}^n=\frac{v_j^n+4v_{j+1\text{}/2}^n+v_{j+1}^n}{6},\\ {\Lambda }_1^x{v}_{j+1\text{}/2}^n=\frac{v_{j+1}^n-v_j^n}{h},\\ {\Lambda }_2^x{v}_{j+1\text{}/2}^n=\frac{4\left(v_j^n-2v_{j+1\text{}/2}^n+v_{j+1}^n\right)}{h^2}.\end{array}$

Индекс j в уравнениях (5) пробегает только целые значения от 0 до N – 1. Нетрудно видеть, что в схеме (5) для дискретизации по времени применяется неявный метод Эйлера. Погрешность аппроксимации схемы (5) (на гладких решениях) составляет $O(h^4+\tau ).$ Схема (5) устойчива при любых $h,\tau >0;$ в случае линейного уравнения переноса в отсутствие расщепления потоков она монотонна при числе Куранта $\kappa ⩾0.25$ (см. [33]).

Удобно поставить в соответствие частям схемы (5) операторы послойного перехода $S_1^{\pm }\left(\tau \right)$ и записать ее в краткой операторной форме:

$v^{n+1}=S_1^{-}\left(\tau \right)v^{*}=S_1^{-}\left(\tau \right)S_1^{+}\left(\tau \right)v^n.$ (6)

Далее в разд. 4, 5 мы будем исследовать точность иной бикомпактной схемы, имеющей повышенный порядок аппроксимации по времени. Для дискретизации по времени в ней применяется трехстадийный L-устойчивый однократно диагонально-неявный метод Рунге–Кутты из [37]. Запишем эту бикомпактную схему в форме, аналогичной (6):

$v^{n+1}=S_3^{-}\left(\tau \right)v^{*}=S_3^{-}\left(\tau \right)S_3^{+}\left(\tau \right)v^n,$ (7)

где оператор послойного перехода $S_3^{-}\left(\tau \right)$ является трехкратной композицией оператора $S_1^{-}(\cdot )$ (формулы для оператора $S_3^{+}\left(\tau \right)$ выписываются аналогично):

$\begin{array}{c}{v}^{\left(1\right)}=S_1^{-}\left(a\tau \right)v^{*},\\ v^{\left(2\right)}=S_1^{-}\left(a\tau \right)\left(D_{21}{v}^{*}+D_{22}{v}^{\left(1\right)}\right),\\ v^{n+1}=S_1^{-}\left(a\tau \right)\left(D_{31}{v}^{*}+D_{32}{v}^{\left(1\right)}+D_{33}{v}^{\left(2\right)}\right).\end{array}$ (8)

Приведем числовые значения коэффициентов a и D_ij в формулах (8) (с двойной точностью):

$\begin{array}{c}a=0.4358665215084579;\\ \begin{array}{l}{D}_{21}=\mathrm{0.3528598198604800,}\\ D_{22}=1-D_{21};\end{array}\\ \begin{array}{c}{D}_{31}=-\mathrm{1.250979895056064,}\\ D_{32}=\mathrm{3.729329662444573,}\\ D_{33}=1-D_{31}-D_{32}.\end{array}\end{array}$

Погрешность аппроксимации схемы (7) составляет $O(h^4+{\tau }^3).$

Обратим внимание на любопытную деталь: справедливость теоремы Годунова для схемы (7) обеспечивается отрицательностью коэффициента D₃₁. Если бы все коэффициенты D_ij были неотрицательны, то мы получили бы линейную монотонную схему с порядком аппроксимации выше первого (оператор послойного перехода $S_1^{\pm }\left(\tau \right)$ сохраняет монотонность при достаточно больших t), а это невозможно. Это, кстати, объясняет, почему условие монотонности схемы (5) имеет вид $\tau ⩾qh,$ где q не зависит от h. Если бы такого ограничения снизу не было, то, положив $\tau =O\left(h^k\right),$ где k > 1, мы получили бы линейную монотонную схему с порядком аппроксимации выше первого, что, как уже говорилось, невозможно.

Схема Русанова

Уравнения этой схемы имеют вид

$\begin{array}{c}{v}_{j+1\text{}/2}^{\left(1\right)}=\frac{v_j^n+v_{j+1}^n}{2}-\frac{R}{3}\left(f_{j+1}^n-f_j^n\right),v_j^{\left(2\right)}=v_j^n-\frac{2R}{3}\left(f_{j+1\text{}/2}^{\left(1\right)}-f_{j-1\text{}/2}^{\left(1\right)}\right),\\ \frac{v_j^{n+1}-v_j^n}{\tau }+\frac{7}{12}\frac{f_{j+1}^n-f_{j-1}^n}{2h}-\frac{1}{3}\frac{f_{j+2}^n-f_{j-2}^n}{4h}+\frac{3}{4}\frac{f_{j+1}^{\left(2\right)}-f_{j-1}^{\left(2\right)}}{2h}=D_j^n,\end{array}$ (9)

где

$\begin{array}{c}R=\frac{\tau }{h},f_j^n=f\left(v_j^n\right),f_{j+1\text{}/2}^{\left(1\right)}=f\left(v_{j+1\text{}/2}^{\left(1\right)}\right),f_j^{\left(2\right)}=f\left(v_j^{\left(2\right)}\right),\\ D_j^n=-\frac{C{h}^3}{24R}\frac{v_{j+2}^n-4v_{j+1}^n+6v_j^n-4v_{j-1}^n+v_{j-2}^n}{h^4}.\end{array}$ 

Член $D_j^n$ — это искусственная вязкость четвертого порядка дивергентности (т.е. пропорциональная четвертой производной решения). Коэффициент вязкости равен $C{h}^3/\left(24R\right)=O\left(h^3\right)$, где C — настраиваемый постоянный параметр. Погрешность аппроксимации схемы (9) составляет $O(h^3+{\tau }^3).$ Схема (9) устойчива при выполнении двух условий:

$\kappa <1\mathrm{и}{\kappa }^2(4-{\kappa }^2)⩽C⩽3.$

Подчеркнем одно важное отличие между схемой Русанова и бикомпактными схемами. В схеме Русанова используются сеточные функции $v_{j+1\text{}/2}^{\left(1\right)},$ $f_{j+1\text{}/2}^{\left(1\right)},$ определенные на ${{\omega }_h}^{\text{'}\text{}\text{'}}.$ Однако значения этих функций носят пассивный (вспомогательный, промежуточный) характер и исключаются в процессе перехода со слоя на слой. В конечном итоге схема Русанова работает с сеточной функцией $v_j^n,$ определенной на ${{\omega }^{\text{'}}}_h.$ Бикомпактные схемы же работают с сеточной функцией $v_i^n,$ определенной на всей ${\omega }_h={{\omega }_h}^{\text{'}}\cup {{\omega }_h}^{\text{'}\text{}\text{'}}.$ В бикомпактных схемах значения решения в полуцелых узлах носят активный характер и выступают как дополнительные степени свободы.

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

3. Критерии оценки точности

Для оценки точности численных схем, описанных в разд. 2, мы будем применять критерии, которые были предложены в работах [22; 23].

Критерий 1: порядок интегральной сходимости

Зафиксируем момент времени t. Введем первообразную точного решения уравнений (1):

$U(x,t)=\underset{0}{\overset{x}{\int }}u(\xi ,t)d\xi .$ (10)

Рассмотрим некоторое численное решение на слое $t_n=t$ (не имеет значения, определено оно на ${\omega }_h$ или $\left.{{\omega }^{\text{'}}}_h\subset {\omega }_h\right).$ Пользуясь значениями этого численного решения, построим на каждом отрезке ${\Delta }_{j+1}=[x_j,x_{j+1}]$ интерполяционный полином как минимум первой степени. Совокупность этих полиномов на всем отрезке [0, X] образует функцию $v_h(x,t).$ Введем с ее помощью дискретный аналог первообразной (10):

$V_h(x_j,t)=\underset{0}{\overset{x_j}{\int }}{v}_h(\xi ,t)d\xi ,x_j\in {{\omega }^{\text{'}}}_h.$(11)

Под интегральной сходимостью решения понимается сходимость численной первообразной (11) к точной первообразной (10) при сгущении сетки по всем направлениям. Порядком интегральной сходимости называется сеточная функция

$k_{h_1}(x_j,t)={\log }_r\frac{\left|V_{h_1}\right(x_j,t)-U(x_j,t\left)\right|}{\left|V_{h_2}\right(x_j,t)-U(x_j,t\left)\right|},x_j\in {{\omega }^{\text{'}}}_{h_1}.$(12)

Здесь $\left|u\right|=\sqrt{u^{T}u}$ — это евклидова норма вектора, а шаги h₁ и h₂ связаны соотношением $h_2=h_1/r,$ где $r>1$ — целое число. Поскольку точное решение уравнений (1), как правило, неизвестно, порядок интегральной сходимости оценивается методом Рунге:

$k_{h_1}(x_j,t)={\log }_r\frac{\left|V_{h_1}\right(x_j,t)-V_{h_2}(x_j,t\left)\right|}{\left|V_{h_2}\right(x_j,t)-V_{h_3}(x_j,t\left)\right|},x_j\in {{\omega }^{\text{'}}}_{h_1},$ (13)

где $h_3=h_2/r.$

Из формул (12), (13) ясно, что выбор точки x = 0 в качестве постоянного нижнего предела интегрирования не вредит общности наших рассуждений. Выбор иной точки или, что то же самое, ненулевой постоянной интегрирования не влияет на значения $k_{h_1}(x_j,t),$ так как эта постоянная уничтожается при вычислении разностей между первообразными. Аналогично выбор в качестве пределов интегрирования $x,X$ вместо $\mathrm{0,}x$ (как в [22; 23]) не влияет на значения $k_{h_1}(x_j,t),$ так как разница между этими двумя интегралами состоит в знаке и аддитивной постоянной; последняя уничтожается при вычислении разностей, а первый устраняется евклидовой нормой.

Приведем формулы расчета численной первообразной (11) для схем из разд. 2.

Для дискретизации по пространству в бикомпактных схемах применяется интегро-интерполяционный метод. Чтобы достичь четвертого порядка аппроксимации, в каждой ячейке сетки, т.е. на отрезке ${\Delta }_{j+1},$ решение приближается интерполяционным полиномом второй степени (см. [34]). Этот полином строится по трем значениям: $v_j^n,$ $v_{j+1\text{}/2}^n,$ $v_{j+1}^n.$ Весьма разумно задать функцию $v_h(x,t)$ на ${\Delta }_{j+1}$ равной этому полиному. Отсюда мы немедленно получаем формулу для численной первообразной:

$\begin{array}{c}{V}_h(x_{j+1},t)=\sum _{m=0}^j\left(A_0^x{v}_{m+1\text{}/2}^n\right)h=\\ =\sum _{m=0}^j\frac{v_m^n+4v_{m+1\text{}/2}^n+v_{m+1}^n}{6}h.\end{array}$

Обратимся к схеме Русанова. Функцию $v_h(x,t)$ на ${\Delta }_{j+1}$ можно задать, например, в виде интерполяционного полинома первой степени, построенного по двум значениям: $v_j^n$ и $v_{j+1}^n.$ Такой подход применялся в работах [22; 23]. Однако у этого подхода есть недостаток: поскольку мы фактически пользуемся квадратурой трапеций, порядок интегральной сходимости ограничен двумя, и мы не можем наблюдать совпадение эмпирического (апостериорного) и теоретического порядков. Предлагается повысить точность расчета численной первообразной. Для этого функция $v_h(x,t)$ на ${\Delta }_{j+1}$ задается как интерполяционный полином второй степени, построенный по трем значениям: $v_{j-1}^n,$ $v_j^n,$ $v_{j+1}^n.$ Такой выбор интерполянта приводит к формуле

$V_h(x_{j+1},t)=\sum _{m=0}^j\frac{5v_{m+1}^n+8v_m^n-v_{m-1}^n}{12}h\left(v_{-1}^n=v_{N-1}^n\right).$

Критерий 2: относительные погрешности инвариантов Римана

Посмотрим на инварианты Римана как на функции от вектора консервативных переменных: $w_i=w_i\left(u\right).$ Относительной погрешностью инварианта Римана $w_i\left(u\right)$ (i = 1,2,3) называется сеточная функция

$rerr\left(w_i\right)(x_j,t)=\frac{\left|w_i\text{}\left(v(x_j,t)\right)-w_i\text{}\left(u(x_j,t)\right)\right|}{\left|w_i\text{}\left(u(x_j,t)\right)\right|},x_j\in {{\omega }_h}^{\text{'}},$ (14)

где $t_n=t$ — наперед заданный фиксированный момент времени. Поскольку точное решение $u(x,t)$ уравнений (1), вообще говоря, неизвестно, вместо него разумно использовать численное решение, полученное на очень подробной сетке с шагом $h_{*}\ll h.$

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

Отметим, что порядок (13) и относительную погрешность (14) можно было бы определить более типичным способом — не как сеточные функции, а как числа, заменив модули на сеточные нормы. Однако рассмотрение этих величин в виде пространственных распределений позволяет проводить более детальный анализ точности тех или иных численных схем.

4. Сеточная сходимость

Проведем расчеты задачи (1), (3) на сетках с шагами $h=\mathrm{0.01,0.005,0.0025}$ (число ячеек $N=\mathrm{1000,2000,4000}$ соответственно). Для бикомпактной схемы (7) возьмем $R=\tau /h=1/7.8$ (число Куранта $\left.\kappa \approx 0.5\right),$ $\delta =\mathrm{1.5.}$ Для схемы Русанова (9) возьмем $R=1/6.4$ $\left(\kappa \approx 0.5\right),$ $C=\mathrm{2.5.}$

На фиг. 2 представлены графики плотности газа и порядка интегральной сходимости при $t=1.3$ и $t=\mathrm{2.5,}$ полученных по схемам (7), (9) на сетке с шагом $h=\mathrm{0.01.}$ В момент времени $t=\mathrm{1.3,}$ когда происходит градиентная катастрофа и ударная волна только зарождается, обе схемы демонстрируют сходимость с третьим порядком (совпадает с теоретическим). Это ожидаемый результат, т.к. область влияния ударной волны только начинает расширяться.

Фиг. 2. Плотность газа (сплошные кривые) и порядок интегральной сходимости (кружки) при t = 1.3 и t = 2,5, полученные по бикомпактной схеме (а) и схеме Русанова (б) на сетке с шагом h = 0.01.

К моменту времени $t=2.5$ эта область уже сформировалась, она занимает интервал $\left(\mathrm{4.3,9.1}\right).$ На нем обе схемы теряют свою высокую точность и сходятся всего лишь с первым порядком.

Проанализируем относительные погрешности инвариантов Римана при $t=2.5$ на той же сетке. Интересно добавить к сравнению нелинейную монотонизированную бикомпактную схему на основе (7) (см. [34]). Для этой схемы параметр монотонизации $C_1=\mathrm{5,}$ монотонной схемой-партнером полагается схема “явный уголок”. Графики погрешностей изображены на фиг. 3. В качестве эталонного решения бралось решение схемы Русанова на сетке с шагом $h_{*}={10}^{-4}$ $\left(N_{*}={10}^5\right).$ Бикомпактная схема и схема Русанова имеют примерно одинаковую точность (первая немного точнее). Погрешности у монотонизированной бикомпактной схемы по сравнению с ними больше на полтора-два порядка (равномерно). Отметим, что средний уровень погрешностей инвариантов w₁, w₃ в области влияния существенно выше, чем вне ее. Это вполне отвечает снижению порядка интегральной сходимости в этой области. Для инварианта w₂ (удельной энтропии) картина несколько иная: средний уровень погрешности заметно выше в подобласти области влияния, на интервале (7.4, 9.1). На нем же, кстати, наблюдается дополнительное ухудшение точности вычисления инвариантов w₁, w₃. Как видно на фиг. 4, в этой подобласти локализовано возмущение энтропии относительно ее начального (фонового) значения. Этому возмущению соответствует участок с изломом на профиле плотности (см. фиг. 2).

Фиг. 3. Относительные погрешности инвариантов Римана при $t=2.5$ , полученные по бикомпактной схеме (сплошные кривые), монотонизированной бикомпактной схеме (штриховые кривые) и схеме Русанова (кружки) на сетке с шагом h = 0.01.

Фиг. 4. Удельная энтропия газа при t = 2.5, рассчитанная по бикомпактной схеме на сетке с шагом h = 0.01.

Итак, мы построили тестовую задачу для уравнений газовой динамики и показали, что на решении этой задачи выбранные нами линейные численные схемы высокого порядка аппроксимации сходятся всего лишь с первым порядком (который достигается в области влияния ударной волны). На первый взгляд, данный результат явно противоречит тому, что было установлено в работах [22–27; 29–31]. В них рассматривалась аналогичная тестовая задача для уравнений мелкой воды, проводились расчеты по тем же самым схемам либо другим линейным схемам, показывалось, что порядок интегральной сходимости не ниже второго. Этот результат распространялся с уравнений мелкой воды вообще на любые гиперболические системы уравнений. Казалось бы, такая индукция верна, поскольку дело шло об анализе аппроксимаций, которые записаны для произвольных векторов $u$ и $f\left(u\right)$ и сами по себе от них не зависят. Однако свойства гиперболической системы уравнений могут оказывать значительное влияние на характер сеточной сходимости линейных численных схем.

5. Сеточная сходимость в случае перманентно изоэнтропического течения

Проиллюстрируем тезис, что характер сеточной сходимости линейных численных схем зависит от свойств решаемой гиперболической системы уравнений. Мы не будем возвращаться к уравнениям мелкой воды (уже рассмотренным в [22–31]) и менять тем самым физическую модель. Вместо этого мы останемся в рамках модели невязкого совершенного газа и упростим уравнения газовой динамики (1) так, чтобы они стали похожи на уравнения мелкой воды.

Предположим, что течение газа является изоэнтропическим перманентно, т.е. во все моменты времени. Отсюда следует равенство

$p=p_0{\left(\frac{\rho }{{\rho }_0}\right)}^{\gamma }\equiv \mathcal{P}\left(\rho \right),$ (15)

справедливое для всех (x, t), т.е. глобальная адиабата Пуассона. Интеграл (15) позволяет исключить уравнение энергии, и от уравнений (1) остаются первые два уравнения — уравнение массы и уравнение импульса:

$\begin{array}{c}{\partial }_{t}u+{\partial }_{x}f\left(u\right)=\mathrm{0,}\\ u=\left[\begin{array}{c}\rho \\ \rho v\end{array}\right],f\left(u\right)=\left[\begin{array}{c}\rho v\\ \rho v^2+\mathcal{P}\left(\rho \right)\end{array}\right].\end{array}$ (16)

Можно показать, что собственные значения l_i и отвечающие им инварианты Римана w_i (i = 1.2) для системы (16) имеют вид

${\lambda }_{\mathrm{1,2}}=v\mp c,w_{\mathrm{1,2}}=v\mp \frac{2c}{\gamma -1},$

где $c=c\left(\rho \right)=\sqrt{\gamma \mathcal{P}\left(\rho \right)/\rho }$ — скорость звука. Уравнения (16) похожи на уравнения мелкой воды в том смысле, что вторые получаются из первых заменой плотности газа r на высоту жидкости H и члена $\mathcal{P}\left(\rho \right)$ на член $\frac{1}{2}g{H}^2$ (где $g=const$ — ускорение свободного падения).

Поясним постановку задачи Коши для уравнений (16). Начальные условия для скорости и плотности газа берутся из формул (3); начальное условие для давления газа является излишним, т.к. оно имелось в виду при написании формулы (15), т.е. было учтено заранее.

Проведем расчеты тестовой задачи для уравнений (16) по тем же численным схемам, на тех же сетках, при тех же параметрах.

Графики плотности газа и порядка интегральной сходимости при t = 1.3 и t = 2.5 приведены на фиг. 5. Сразу же обращает на себя внимание тот факт, что и бикомпактная схема, и схема Русанова в области влияния ударной волны сходятся со вторым, а не первым порядком. Например, схема Русанова сходилась именно так на решении аналогичной тестовой задачи для уравнений мелкой воды (см. [23]). Заметим, что профили плотности при t = 2.5 тоже выглядят немного иначе: отсутствует участок с изломом.

Фиг. 5. Плотность газа (сплошные кривые) и порядок интегральной сходимости (кружки) при t = 1.3 и t = 2.5, полученные по бикомпактной схеме (а) и схеме Русанова (б) на сетке с шагом h = 0.01.

На фиг. 6 представлены графики относительных погрешностей инвариантов Римана в момент времени t = 2.5. Перечислим основные отличия от результатов в случае полной системы уравнений газовой динамики (1). Точность вычисления инварианта w₁ в области влияния ударной волны для бикомпактной схемы примерно на порядок-полтора хуже, чем для схемы Русанова. Погрешности w₁ вне области влияния и погрешности w₂ во всей расчетной области для этих схем примерно одинаковы. Разница в точности между бикомпактной схемой и ее монотонизированной версией оказывается еще более глубокой: для инварианта w₂ она составляет целых три порядка.

Фиг. 6. Относительные погрешности инвариантов Римана при t = 2.5, полученные по бикомпактной схеме (сплошные кривые), монотонизированной бикомпактной схеме (штриховые кривые) и схеме Русанова (кружки) на сетке с шагом h = 0.01.

Рассмотренный нами контрпример (тестовая задача Коши) выявил, что одна и та же линейная численная схема может сходиться в области влияния ударной волны то с первым, то со вторым порядком, в зависимости от того, какая гиперболическая система решается. Это явление нельзя списать на специфику аппроксимации, например, бикомпактной схемы, т.к. схема Русанова, будучи схемой другого класса, обнаружила совершенно такое же поведение. В связи с этим выбор именно бикомпактной схемы и схемы Русанова не является принципиальным. Понижение порядка не может быть объяснено и тем, что сходимость якобы не установилась и действительный порядок еще не достигнут; мы специально брали достаточно подробные сетки (об этом свидетельствует низкий уровень погрешностей вне области влияния: от 10^–10 до 10^–8). Значит, причиной вариаций в порядке сходимости является сама гиперболическая система уравнений. Но в чем отличие между системами (1) и (16)?

Известно, что любая гиперболическая система, состоящая из двух уравнений, всегда может быть записана в инвариантах (см. [38]). Иначе обстоит дело для гиперболических систем из трех или более уравнений: инвариантов Римана для них в общем случае нет. В зависимости от постановки задачи инварианты Римана могут существовать в течение некоторого времени во всей расчетной области, а затем не существовать в некоторой подобласти. Например, это происходит в тестовой задаче (1), (3): одновременно с формированием ударной волны возникает (и далее расширяется) зона за ее фронтом, где $s=w_2\ne const$ и, как следствие, не существуют инварианты w₁, w₃.

Возможность записи в инвариантах объединяет уравнения мелкой воды и уравнения изоэнтропической газовой динамики и отличает обе этих системы от полных уравнений газовой динамики. По-видимому, имеется некая связь между существованием инвариантов Римана и характером сеточной сходимости линейных численных схем. Анализ этой связи составляет предмет для дальнейших исследований.

Заключение

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

Исследована сеточная сходимость бикомпактной схемы и схемы Русанова на решении данной задачи. Обе схемы линейные, немонотонные, высокого третьего порядка аппроксимации. Порядок сходимости определялся поточечно (в узлах, локально) путем сравнения дискретных первообразных численных решений, полученных на вложенных сетках. Установлено, что выбранные схемы сходятся с третьим порядком до момента образования ударной волны и с первым порядком — после этого момента. Понижение порядка сходимости происходит в области влияния ударной волны. Однако, несмотря на формальное совпадение порядков аппроксимации, немонотонная бикомпактная схема оказывается на полтора-два порядка точнее, чем ее монотонизированная версия, если судить по относительным погрешностям инвариантов Римана.

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