Журнал вычислительной математики и математической физики, 2020, T. 60, № 5, стр. 884-899

1. ВВЕДЕНИЕ

В классической работе [1] было показано, что среди линейных разностных схем для линейного уравнения переноса нет монотонных схем с минимальным порядком аппроксимации выше первого по всем независимым переменным (здесь и далее под “аппроксимацией” мы будем понимать классическую аппроксимацию на гладких решениях, если явно не указано иное). Вследствие этого “запрета Годунова” разностные схемы повышенного порядка аппроксимации для гиперболических систем законов сохранения при сквозном счете генерируют осцилляции в областях высоких градиентов искомых функций. В литературе предложены различные способы подавления этих осцилляций или, другими словами, монотонизации численного решения. Одним из распространенных способов монотонизации является нелинейная коррекция численных потоков. Такая коррекция используется в MUSCL-схемах [2], TVD-схемах [3], WENO-схемах [4], NED-схемах [5], DG-схемах [6], САВARET-схемах [7]. Далее эти схемы будем сокращенно называть NFC-схемами (Nonlinear Flux Correction). Наряду с ограничителями (limiters) численных потоков для уменьшения осцилляций около разрывов в схемах высокого порядка аппроксимации применяются искусственная диссипация [8]–[11], численные фильтры [12]–[15], а также различные процедуры гибридизации [16]–[20], в которых обычно комбинируются монотонизированная и немонотонная схемы.

Однако коррекция потоков около сильных разрывов решения в схемах высокого порядка приводит к нежелательным эффектам. В работах [21]–[23] было показано, что NFC-схемы имеют не более чем первый порядок локальной сходимости в областях влияния ударных волн. Также снижение порядка сходимости NFC-схем до первого имеет место и для интегралов от искомых функций, рассчитанных по квадратурной формуле повышенного порядка точности на интервалах, одна из границ которых находится в области влияния ударной волны [24]–[29]. Такую сходимость интегралов от сеточных функций авторы работ [24]–[29] называют интегральной сходимостью. Как указано в [24]–[29], причина падения порядка сходимости до первого связана с тем, что в NFC-схемах коррекция численных потоков приводит к снижению их гладкости и, как следствие, снижению порядка аппроксимации ε-условий Гюгонио на фронтах ударных волн [30]. О связи гладкости численных потоков и порядка сходимости разностных схем говорят также и результаты работ [24], [27], в которых установлено, что немонотонная схема Русанова [31] и немонотонная компактная схема из [23] с достаточно гладкими численными потоками имеют не менее чем второй порядок интегральной сходимости в областях, содержащих сильные разрывы. В отличие от монотонизированных NFC-схем указанные немонотонные схемы обладают повышенным порядком сходимости в областях влияния ударных волн. Отметим, что эти немонотонные схемы генерируют ложные осцилляции, локализованные около ударных волн.

В работах [32]–[34] были предложены центральные компактные схемы на неразнесенных декартовых сетках для решения квазилинейных гиперболических систем законов сохранения. Схемы получены методом прямых и имеют четвертый порядок пространственной аппроксимации. Шаблон вышеупомянутых компактных схем по каждой пространственной переменной состоит из двух целых узлов, по этой причине эти схемы стали называться бикомпактными. Для значений искомых функций в дополнительных полуцелых узлах шаблона выводятся специальные уравнения – следствия исходных гиперболических уравнений. В итоге дискретные уравнения схем имеют первый разностный порядок по пространственным переменным для искомых сеточных функций в целых узлах [35]. Поэтому решение этих уравнений требует обращения только блочно-двухдиагональных матриц, хотя компактные схемы [32]–[34] – неявные, абсолютно устойчивые и могут быть использованы в широком диапазоне локальных чисел Куранта. Особенностью этих схем является то, что они сохраняют порядок пространственной аппроксимации на произвольно неравномерных сетках из целых узлов. Кроме того, бикомпактные схемы имеют улучшенные дисперсионные свойства по сравнению с классическими трехточечными компактными схемами [35] того же порядка аппроксимации. Отметим, что полудискретные центральные бикомпактные схемы являются бездиссипативными. Диссипативные свойства полностью дискретных бикомпактных схем определяются диссипацией методов интегрирования по времени полудискретных схем. В качестве методов интегрирования по времени в [32]–[34] предлагается использовать экономичные неявные методы типа диагонально-неявных методов Рунге–Кутты (DIRK) [36]. В работах [37], [38] предложены обобщения этих схем на случай многомерных гиперболических уравнений, а в работе [39] – обобщения на случай произвольного четного (четвертого, шестого восьмого и так далее) порядка аппроксимации пространственных производных.

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

2. БИКОМПАКТНЫЕ СХЕМЫ БЕЗ МОНОТОНИЗАЦИИ

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

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

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

где 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. Система (3) двух обыкновенных дифференциальных уравнений (ОДУ) является системой уравнений для двух сеточных функций ${{{\mathbf{u}}}_{j}}(t)$ и ${{{\mathbf{u}}}_{{j + 1/2}}}(t)$, определенных в целых ${{x}_{j}}$ и полуцелых ${{x}_{{j + 1/2}}} = {{({{x}_{j}} + {{x}_{{j + 1}}})} \mathord{\left/ {\vphantom {{({{x}_{j}} + {{x}_{{j + 1}}})} 2}} \right. \kern-0em} 2}$ узлах, причем разностный порядок системы (3) для сеточной функции ${{{\mathbf{u}}}_{j}}(t)$ равен единице, а для сеточной функции ${{{\mathbf{u}}}_{{j + 1/2}}}(t)$ – нулю [35]. Отметим, что понятие разностного порядка для одного разностного уравнения относительно одной сеточной функции можно найти, например, в монографии [40]. Обобщение этого понятия в случае системы разностных уравнений для нескольких сеточных функций дано в [35]. Разностный порядок дискретного уравнения для сеточной функции равен числу граничных условий, необходимых для корректной постановки разностной задачи.

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

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

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

Полностью дискретные бикомпактные схемы получаются путем применения какого-либо устойчивого численного метода интегрирования по времени системы ОДУ полудискретной схемы (3) . В данном разделе будут рассмотрены свойства разностных решений двух полностью дискретных бикомпактных схем на равномерных сетках с постоянным пространственным шагом ${{h}_{{x,j + 1/2}}} = h$ и постоянным временным шагом $\tau $.

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

Разностные решения второй схемы, обозначаемой как Rich2B4, получаются уточнением разностных решений базовой бикомпактной схемы со вторым порядком точности по времени методом экстраполяции Ричардсона [42], [43]. Проводятся два расчета по базовой бикомпактной схеме на сетках с временными шагами $\tau $ и $\tau {\text{/}}2$. При этом пространственный шаг сетки h берется таким, чтобы вклад ошибки пространственной дискретизации в погрешность численного решения оказался малым по сравнению со вкладом ошибки дискретизации по времени. При сгущении пространственно-временной сетки поддерживается постоянство величины отношения шагов $r = {\tau \mathord{\left/ {\vphantom {\tau h}} \right. \kern-0em} h}$. Решение экстраполяционной бикомпактной схемы Rich2B4 в момент времени ${{t}^{n}} = n{\kern 1pt} \tau $ рассчитывается по следующей формуле:

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

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

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

где H, Q и $u = {Q \mathord{\left/ {\vphantom {Q H}} \right. \kern-0em} H}$ – глубина, расход и горизонтальная скорость жидкости, g – ускорение свободного падения (в расчeтах полагалось, что g = 10). Периодические начальные условия (2) для системы (1), (5) задаются в виде где $a = 2$, $b = 10$, $X = 10$.

Результаты расчета задачи Коши (1), (5), (6) для моментов времени $T = 0.5$ и $T = 1.0$ представлены на фиг. 1–4 и на фиг. 5–10 соответственно. Из данных на фиг. 1 и 5 видно, что при $t \approx 0.5$ и $x \approx 6$ формируется прерывная волна, которая движется в положительном направлении оси Ox и достигает положения $x \approx 9$ в момент времени $t = 1.0$.

В согласии с анализом монотонности [33], [34] базовая бикомпактная схема DIRK1B4 при определенных шагах сетки не генерировала нефизические осцилляции при расчете профилей H, Q. На фиг. 1, 5 и 6 показаны профили глубины H, рассчитанные по схеме DIRK1B4. Также схема DIRK1B4 не генерирует нефизических осцилляций при расчете первых разностных производных величин H, Q. На фиг. 7a и 7б представлено распределение величины ${{\left( {{{H}_{x}}} \right)}_{k}} = 2{{\left( {{{H}_{{k + 1/2}}} - {{H}_{k}}} \right)} \mathord{\left/ {\vphantom {{\left( {{{H}_{{k + 1/2}}} - {{H}_{k}}} \right)} h}} \right. \kern-0em} h}$ вдоль оси Oх как в целых $k = j$, так и полуцелых $k = j + 1{\text{/}}2$ узлах на момент времени $T = 1.0$. При этом на фиг. 7a показана верхняя часть данного распределения при значениях ${{\left( {{{H}_{x}}} \right)}_{k}} > 0$, а на фиг. 7б – распределение в окрестности ударной волны.

На фиг. 2–4, 6 и 8–10 представлены результаты расчетов по экстраполяционной бикомпактной схеме Rich2B4 на последовательности сгущающихся пространственных сеток с числом ячеек $N = 1250$, 2500, 5000 и 10 000 при отношении шагов $r = 0.05$. Расчеты показали, что схема Rich2B4 не генерирует нефизические осцилляции при $T = 0.5$, а при $T = 1.0$ эта схема генерирует нефизические осцилляции, локализованные около ударной волны, как это видно из фиг. 6. Для исследования вопроса о поведении погрешности решения схемы Rich2B4 при сгущении сетки желательно иметь аналитическое решение нелинейной задачи Коши (1), (5), (6). Поскольку такое решение не известно, в качестве него будем рассматривать решение, рассчитанное по схеме Русанова [31] на очень мелкой сетке с пространственным шагом h = 10^–4 и отношением шагов $r = 0.05$. Схема Русанова позволяет получить решение задачи Коши (1), (5), (6) без ложных осцилляций для момента времени $T = 0.5$. Решение этой схемы для $T = 1.0$ содержит нефизические осцилляции, локализованные в малой окрестности ударной волны, как это видно из фиг. 6. Однако этот дефект решения схемы Русанова приемлем для оценки погрешности решений бикомпактных схем, рассчитанных на сетках с шагом $h \geqslant {{10}^{{ - 3}}}$. На фиг. 2–4 и 8–10 показаны распределения ошибок решений схемы Rich2B4 для глубины H, а также для инвариантов Римана ${{w}_{1}},\;{{w}_{2}}$ (см. [29])

гиперболической системы уравнений (1), (5). В формулах (7) величина $с = \sqrt {gH} $ есть скорость распространения малых возмущений. Из данных на фиг. 2–4 и 8–10 следует, что локальная ошибка численных решений для H, ${{w}_{1}},\;{{w}_{2}}$ при сгущении сетки убывает со вторым порядком везде за исключением тех областей, где ошибки резко возрастают или резко падают. Заметим, что ошибка численных решений резко возрастает в окрестности ударной волны. На фиг. 1 и 5 точками в каждом пятидесятом целом узле показаны порядки p локальной сходимости глубины H, определенные согласно формуле по абсолютным погрешностям ${\text{err}}(H)$ расчетов H на двух сетках при $N = 5000$. Из анализа данных на фиг. 7a и 8, соответствующих расчетам на сетке с $N = 10\,000$, видно, что резкое падение ошибки в расчете величины H происходит в окрестностях точек экстремума первой разностной производной ${{H}_{x}}$, т.е. в окрестности точек, где решение ведет себя как линейная функция x. Такое поведение ошибки характерно для схем второго и более высокого порядка точности.

3. МОНОТОНИЗИРОВАННАЯ БИКОМПАКТНАЯ СХЕМА

Рассмотрим бикомпактную схему SDIRK3B4, получающуюся путем применения L-устойчивого однократно диагонально-неявного метода Рунге–Кутты третьего порядка точности с тремя неявными стадиями SDIRK3 [45] к интегрированию системы ОДУ полудискретной бикомпактной схемы (3) . Вследствие “запрета Годунова” схема SDIRK3B4 является немонотонной и порождает нефизические осцилляции в областях высоких градиентов решения. Мы будем применять процедуру консервативной монотонизации, предложенную в работе [46], в процессе счета по схеме SDIRK3B4 на каждом временном слое для уменьшения нефизических осцилляций.

На фиг. 11–14 и 15–18 представлены результаты расчетов задачи Коши (1), (5), (6) по монотонизированной бикомпактной схеме SDIRK3B4 для моментов времени $T = 0.5$ и $T = 1.0$ соответственно. Расчеты выполнены на последовательности сгущающихся сеток с числом ячеек $N = 1250$, 2500, 5000 и 10 000 при отношении шагов $r = 0.07$. Профили глубины H, вычисленные по монотонизированной схеме SDIRK3B4, не имеют заметных осцилляций и представлены на фиг. 6, 11 и 15.

Важно отметить, что задача Коши для гиперболической системы уравнений (1), (5) с периодическими условиями (6) описывает докритическое течение жидкости c горизонтальной скоростью

для которой характеристические скорости ${{\lambda }_{1}} = u - c$, ${{\lambda }_{2}} = u + c$ [29] инвариантов Римана ${{w}_{1}},{{w}_{2}}$ удовлетворяют неравенствам

Неравенства (9) указывают на то, что информация о величинах ${{w}_{1}}$ и ${{w}_{2}}$ распространяется в отрицательном и положительном направлениях оси Ox соответственно. Ниже будет показано, как неравенства (9) сказываются на профилях величин ${{w}_{1}}$ и ${{w}_{2}}$, а также на профиле величины

Профили абсолютных ошибок расчета глубины H и инвариантов Римана ${{w}_{1}},{{w}_{2}}$ показаны на фиг. 12–14 на момент времени $T = 0.5$, когда формируется нестационарная ударная волна, движущаяся вправо. На ударной волне генерируются осцилляции величины ${{w}_{1}}$, которые переносятся за ее фронт в согласии с неравенством ${{\lambda }_{1}} < 0$. Из фиг. 13 видны области распространения осцилляций малой амплитуды для разных шагов сетки. Поскольку величина H зависит от величины ${{w}_{1}}$ согласно (10), то на фиг. 12 также можно наблюдать осцилляции за ударной волной. В то же время на фиг. 14 осцилляции за ударной волной отсутствуют, что согласуется с неравенством ${{\lambda }_{2}} > 0$. На фиг. 11 в каждом пятидесятом целом узле точками показаны порядки p локальной сходимости величины H, рассчитанные по формуле (8) при $N = 5000$. Видно, что порядок этой сходимости близок к трем везде, за исключением зон, где погрешность расчетов резко уменьшается. Подобный порядок сходимости имеет место и для инвариантов ${{w}_{1}},\;{{w}_{2}}$, что видно из фиг. 13–14.

На профилях абсолютных ошибок H и инвариантов Римана ${{w}_{1}},\;{{w}_{2}}$, рассчитанных по монотонизированной схеме SDIRK3B4 и соответствующих моменту времени $T = 1.0$, осцилляции малой амплитуды за ударной волной наблюдаются для H, ${{w}_{1}}$ и отсутствуют для ${{w}_{2}}$. Причем по сравнению со случаем $T = 0.5$ эти осцилляции имеют бóльшую амплитуду. Чтобы яснее видеть динамику убывания ошибок со сгущением сетки, разностные решения для H и ${{w}_{1}}$ были усреднены по пространству в области $x \in [6,8.9]$ за ударной волной. Для осреднения решения $U({{x}_{j}},T)$ использовалась следующая формула [29]:

где $U = H,{{w}_{1}}$, $l = 5$, $\varepsilon = 5h$, а интеграл рассчитывается по квадратурной формуле трапеций.

Профили ошибок H и инвариантов Римана ${{w}_{1}},\;{{w}_{2}}$, представленные на фиг. 16–18, показывают следующую динамику при сгущении сетки. В области влияния $x \in [4,9]$ нестационарной ударной волны (то есть в следе волны [25]) убывание ошибок для монотонизированной схемы SDIRK3B4 при сгущении сетки соответствует не более чем первому порядку локальной сходимости. Вне области влияния порядок локальной сходимости близок к третьему порядку. На фиг. 15 в каждом пятидесятом целом узле точками показаны порядки p локальной сходимости величины H, рассчитанные по формуле (8) при $N = 5000$. Таким образом, процедура монотонизации [46] численного решения приводит к снижению порядка локальной сходимости схемы в следе нестационарной ударной волны.

4. ЗАКЛЮЧЕНИЕ

Исследована точность двух бикомпактных схем повышенного порядка классической аппроксимации по времени и четвертого по пространству при расчете разрывных решений квазилинейной гиперболической системы законов сохранения. Первая из изученных бикомпактных схем получена экстраполяцией по Ричардсону бикомпактной схемы первого порядка аппроксимации по времени и имеет второй порядок точности по времени на гладких решениях. Согласно “запрету Годунова”, она является немонотонной. Отметим, что экстраполяция в этой схеме реализуется пассивно, то есть после завершения расчетов по монотонной схеме первого порядка по времени на двух вложенных сетках и только для тех моментов времени, для которых требуется найти решение и исследовать его точность. В результате такой экстраполяции получается численное решение повышенной точности с локализованной немонотонностью небольшой амплитуды. Вторая схема получена путем применения процедуры консервативной монотонизации [46] к бикомпактной схеме третьего порядка аппроксимации по времени. В качестве разрывного решения рассмотрено решение задачи Коши с нестационарной ударной волной для гиперболической системы уравнений теории мелкой воды. В результате численных экспериментов показано, что первая бикомпактная схема имеет второй порядок сходимости в области влияния нестационарной ударной волны. Вторая бикомпактная схема, несмотря на высокую точность на гладких решениях, снижает свой порядок сходимости до первого порядка в области влияния ударной волны.

Таким образом, проведенное в настоящей работе исследование подтверждает вывод работ [23], [24], [29] о невозможности одновременно с высокой точностью монотонно локализовать двигающиеся с переменной скоростью сильные разрывы и сохранить повышенный порядок сходимости в областях их влияния.

Авторы выражают благодарность Н.А. Хандеевой за любезно предоставленные результаты расчетов по схеме Русанова.