Доклады Российской академии наук. Математика, информатика, процессы управления, 2020, T. 492, № 1, стр. 79-84

1. В работах [1, 2] было показано, что NFC (Nonlinear Flux Correction) схемы сквозного счета, к которым относятся MUSCL-схемы, TVD-схемы, WENO-схемы, DG-схемы, САВARET-схемы (ссылки на эти схемы см., например, в [2]), имеют не более чем первый порядок локальной сходимости в областях влияния ударных волн. Как указано в [1–4], причина падения порядка сходимости до первого связана с тем, что в NFC-схемах коррекция численных потоков приводит к снижению их гладкости и, как следствие, к снижению порядка аппроксимации ε-условий Гюгонио на фронтах ударных волн [1].

В работе [3] был предложен метод построения комбинированных схем сквозного счета, которые монотонно локализуют фронты ударных волн и одновременно имеют повышенную точность в областях гладкости рассчитываемых обобщенных решений. Решение в комбинированной схеме строится путем коррекции решения базисной немонотонной схемы, которое находится во всей расчетной области и имеет повышенный порядок сходимости в областях влияния ударных волн, где решение является гладким. В областях больших градиентов, где решение базисной схемы имеет нефизические осцилляции, оно корректируется путем численного решения внутренних начально-краевых задач по одной из NFC-схем. Комбинированная схема [3] состояла из разнородных схем: в качестве базисной использовалась неявная симметричная компактная схема третьего порядка слабой аппроксимации, а в качестве внутренней NFC‑схемы – монотонизированный вариант явной схемы CABARET второго порядка точности на гладких решениях. Недостаток, связанный с указанной разнородностью схем, был устранен в работе [4], где базисная и внутренняя схемы были немонотонным и монотонизированным вариантами разрывного метода Галёркина. В работах [3, 4] свойства комбинированных схем демонстрировались на расчете модельной одномерной задачи Коши, в которой решение содержит одну изолированную ударную волну. Применение комбинированной схемы [3] к случаю взаимодействующих одномерных ударных волн показывает [5], что реализация комбинированной схемы этого типа значительно усложняется по сравнению со случаем одной ударной волны из-за слияния и разъединения расчетных областей для внутренних начально-краевых задач. Сложность реализации комбинированной схемы [3] увеличивается в случае многомерных задач. Поэтому в настоящей работе предлагается принципиально новый метод построения комбинированной схемы сквозного счета с вышеуказанными свойствами. В этом методе построения не требуется решать внутренние начально-краевые задачи, поэтому он является существенно более простым по сравнению с методом построения комбинированных схем в работах [3–5].

2. Для квазилинейной гиперболической системы законов сохранения

рассмотрим задачу Коши на временном интервале $(0,T]$ с периодическими начальными условиями где символ ${{\partial }_{x}} \equiv \frac{\partial }{{\partial x}}$; ${\mathbf{u}}(x,t)$ – искомая, а ${\mathbf{f}}({\mathbf{u}})$ – заданная вектор-функции из m компонент; X > 0 – период.

Полудискретная бикомпактная схема четвертого порядка аппроксимации для системы уравнений (1) имеет вид (см. [6])

где N – число ячеек на отрезке оси Ox с длиной X. Далее будем рассматривать отрезок $[0,X]$. Разностные операторы $A_{0}^{x}$, $\Lambda _{1}^{x}$, $\Lambda _{2}^{x}$ определяются для произвольной сеточной функции U так: где ${{h}_{{x,j + 1/2}}} = {{x}_{{j + 1}}} - {{x}_{j}}$ – в общем случае переменный шаг сетки по x, обозначаемый далее как h. Система двух обыкновенных дифференциальных уравнений (ОДУ) (3) решается относительно двух сеточных функций ${{{\mathbf{u}}}_{j}}(t)$ и ${{{\mathbf{u}}}_{{j + 1/2}}}(t)$, определенных в целых x_j и полуцелых ${{x}_{{j + 1/2}}} = \frac{{{{x}_{j}} + {{x}_{{j + 1}}}}}{2}$  узлах.

Начальные условия (2) для функций ${{{\mathbf{u}}}_{j}}(t)$ и ${{{\mathbf{u}}}_{{j + 1/2}}}(t)$ будут следующими:

Граничное условие периодичности для сеточной функции ${{{\mathbf{u}}}_{j}}(t)$ имеет вид

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

3. Полностью дискретные бикомпактные схемы получаются путем применения какого-либо устойчивого численного метода интегрирования по времени к системе ОДУ (3). В данной работе рассмотрены две такие схемы.

Первая из схем, далее обозначаемая как DIRK1B4, получается путем интегрирования уравнений (3) по неявному методу Эйлера. Она имеет погрешность аппроксимации $O(\tau ,{{h}^{4}})$, является абсолютно устойчивой и консервативной [8]. В работе [8] исследованы дисперсионные и диссипативные свойства, а также свойства монотонности схемы DIRK1B4. Аналитически показано, что эта схема для линейного скалярного уравнения переноса является монотонной по Годунову при числе Куранта $\kappa \geqslant 0.25$. Отметим, что коэффициент схемной вязкости у данной схемы пропорционален $\kappa $ и положителен при $\kappa > 0$. Численный эксперимент [8] показал, что схема DIRK1B4 для квазилинейного уравнения Хопфа монотонна, если локальное число Куранта $\kappa \geqslant 0.15$.

Численные решения второй схемы, обозначаемой как Rich2B4, получаются уточнением разностных решений схемы DIRK1B4 со вторым порядком точности по времени методом экстраполяции Ричардсона. Для этого проводятся два вычисления во всей расчетной области по схеме DIRK1B4 на сетках с временными шагами $\tau $ и $\frac{\tau }{2}$. При сгущении пространственно-временной сетки поддерживается постоянство величины отношения шагов r = $\frac{\tau }{h}$. Решение экстраполяционной бикомпактной схемы Rich2B4 в момент времени ${{t}_{n}} = n{\kern 1pt} \tau $ рассчитывается по формуле

где первый параметр сеточной функции ${\mathbf{u}}_{k}^{n}$ указывает на схему, а второй – на временной шаг сетки, на которой определялось решение.

4. В качестве конкретной гиперболической системы (1) возьмем систему уравнений первого приближения теории мелкой воды в случае прямоугольного русла без учета донного трения, для которой

где H, Q и ${v}$ = $\frac{Q}{H}$ – глубина, расход и горизонтальная скорость жидкости, g – ускорение свободного падения (в расчетах полагалось, что g = 10). Периодические начальные условия (2) для системы (1), (5) задаются так: где a = 2, b = 10, X = 10.

Результаты расчета задачи Коши (1), (5), (6) для момента времени T = 1 представлены на рис. 1, 3, 4 (на отрезке $[0,X]$) и на рис. 2 (в окрестности фронта ударной волны).

На рис. 1 приведены профили глубины H и величины ${{\left( {{{H}_{x}}} \right)}_{k}} = 2\frac{{{{H}_{{k + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}} - {{H}_{k}}}}{h}$ как в целых k = j, так и полуцелых $k = j + \frac{1}{2}$ узлах на отрезке $[0,X]$. Эти профили получены по схеме DIRK1B4 на мелкой сетке с N = 2000. Значения величины H, рассчитанной по той же схеме на сетках с N = 2000 и N = 100, а также по схеме Rich2B4 на сетке с N = = 100, показаны в окрестности ударной волны на рис. 2. Укажем на то, схема DIRK1B4 выдает профили H и H_x, свободные от нефизических осцилляций, в отличие от схемы Rich2B4, порождающей осцилляции вблизи ударной волны.

Рассмотрим вопрос о поведении погрешности решения схемы Rich2B4 при сгущении сетки. В силу неизвестности аналитического решения вместо него для вычисления погрешностей берется решение, рассчитанное по схеме Русанова [9] третьего порядка аппроксимации на очень мелкой сетке с h = 10⁻⁴ и r = 0.05. Решение этой схемы при T = 1 содержит нефизические осцилляции, локализованные в малой окрестности ударной волны, что демонстрирует рис. 2. Однако этот дефект решения схемы Русанова приемлем для оценки погрешности решений схемы Rich2B4, рассчитанных на сетках с шагом h > 10^–3. На рис. 3 показаны распределения абсолютных погрешностей err(H) глубины H, полученные по схеме Rich2B4 на последовательности сгущающихся сеток с N = 500, 1000, 2000 и 4000 при отношении шагов r = 0.05.

Кроме того, по err(H) были оценены порядки p локальной сходимости глубины H:

На рис. 1 точками в каждом двадцатом целом узле показано распределение p для N = 2000. Из данных на рис. 1 и 3 следует, что при сгущении сетки локальные абсолютные погрешности H убывают со вторым порядком всюду за исключением тех областей, где ошибка резко возрастает или убывает. Очевидно, резкое возрастание имеет место в окрестности ударной волны. Из рис. 1 видно, что резкое убывание происходит в окрестностях точек экстремума H_x, т.е. там, где H ведет себя как линейная функция x. Такое поведение ошибки характерно для схем второго порядка точности и выше.

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

Решение комбинированной схемы CombinB4 зададим по формулам:

где u_s – s-я компонента искомой вектор-функции u, C > 0 – задаваемый параметр комбинированной схемы, Ω_N – пространственная сетка на отрезке длины периода X, состоящая из N ячеек. В формулах (7), (8) u_s, $u_{s}^{A}$, $u_{s}^{B}$ – это сеточные функции, у которых опущен нижний индекс k = j, $j + \frac{1}{2}$ и верхний индекс n, поскольку эти формулы задают локальные в пространстве и во времени связи между указанными сеточными функциями. Согласно формулам (7), (8) базисной схемой в предложенной комбинированной схеме является бикомпактная схема Rich2B4, а решение внутренней схемы, область определения которой расположена в окрестности фронта ударной волны и задается неравенством $\left| {{{D}_{s}}} \right| > C{{N}_{s}}$, является нелинейной комбинацией решений бикомпактных схем DIRK1B4 и Rich2B4.

На рис. 2 вблизи фронта ударной волны показаны профили глубины H, посчитанные по комбинированной бикомпактной схеме CombinB4, по схеме WENO5 и по комбинированной схеме CombinC3 [3] на сетке с N = 100. У схемы CombinB4 параметр C = 10^–2. Схема WENO5 имеет пятый порядок аппроксимации по пространству и третий по времени. У комбинированной схемы CombinC3 [3] базисной схемой является компактная схема третьего порядка классической и слабой аппроксимации, а внутренней схемой – монотонизированная версия CABARET-схемы второго порядка классической аппроксимации. Из данных на рис. 2, во-первых, видно, что решение схемы CombinB4 не имеет осцилляций. Во-вторых, схема CombinB4 размазывает фронт ударной волны больше, чем схема CombinC3, но меньше, чем схема WENO5. Другими словами, уровень диссипации схемы CombinB4 находится между уровнями диссипации схем WENO5 и CombinC3.

На рис. 4 показаны распределения относительных погрешностей инвариантов Римана ${{w}_{{1,2}}} = {v} \mp 2c$ для схемы CombinB4 на последовательности сгущающихся сеток. Величина c = $\sqrt {gH} $ есть скорость распространения малых возмущений. Из данных на этом рисунке следует, что локальные относительные погрешности инвариантов w₁, w₂ убывают при сгущении сетки со вторым порядком подобно локальным абсолютным погрешностям H, см. рис. 3. На рис. 4 дается сравнение погрешностей ${{w}_{{1,2}}}$, полученных по схемам CombinB4 (штриховые линии) и WENO5 (точки) на сетке с N = 2000. Видно, что в области влияния нестационарной ударной волны (в следе волны [3]), т.е. на отрезке $x \in [4,9]$, погрешность схемы CombinB4 в 10–15 раз меньше, чем у схемы WENO5. Отметим, что благодаря простоте построения комбинированной бикомпактной схемы ее можно использовать для расчетов многомерных задач с взаимодействующими ударными волнами без принципиальных затруднений.