Журнал вычислительной математики и математической физики, 2021, T. 61, № 11, стр. 1747-1758

Сравнение диссипативно-дисперсионных свойств компактных разностных схем для численного решения уравнения адвекции

Е. Н. Аристова^1, *, Г. О. Астафуров^1, **

¹ Институт прикладной математики им. М.В. Келдыша РАН
125047 Москва, Миусская пл., 4, Россия

^* E-mail: aristovaen@mail.ru
^** E-mail: astafurov.gleb@yandex.ru

Поступила в редакцию 10.11.2020
После доработки 14.02.2021
Принята к публикации 07.07.2021

DOI: 10.31857/S0044466921110028

Аннотация

В работе исследованы диссипативно-дисперсионные свойства эрмитовой характеристической схемы для решения одномерного уравнения адвекции. Метод основан на эрмитовой интерполяции, использующей не только значения функции в узлах, но также и значения пространственной производной функции в узлах. Используется вычисление производных на новом шаге по времени, обеспечивающее правильное перераспределение входящих потоков по выходным граням. Отметим, что схема строится в рамках одной ячейки, что позволяет ее отнести к классу бикомпактных схем. Восстановление производных на новом слое по времени производится с использованием интегрального среднего и формулы Эйлера–Маклорена. Проведен сравнительный анализ данной схемы с современными консервативными схемами, такими как бикомпактная схема Б.В. Рогова и схема В.М. Головизнина и Б.Н. Четверушкина. Показано, что эрмитова характеристическая схема обладает малой диссипацией и экстрамалой дисперсией для схем своего класса. Дисперсия эрмитовой характеристической схемы меньше дисперсии полудискретной бикомпактной схемы Рогова. В свою очередь, последняя схема при реализации метода трапеций аппроксимации по времени обладает нулевой диссипацией. Близкие идеи использования характеристических схем с дополнительным алгоритмом, обеспечивающим консервативность, использованы в наиболее простой в реализации схеме Головизнина–Четверушкина. Для сравнения выбраны схемы с компактным шаблоном и близкими чертами, используемыми для замыкания разностной схемы. Библ. 24. Фиг. 8.

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

1. ВВЕДЕНИЕ

Существует огромное число схем для численного решения уравнения переноса. Приложение этих схем к реальным задачам делится на два больших класса: решение нелинейных систем уравнений типа уравнений газовой динамики или решение уравнения переноса излучения и/или нейтральных частиц. Специфика каждого класса задач предъявляет свои требования к качеству используемых разностных схем. Для решения задач высокотемпературной радиационной газовой динамики требуются схемы, позволяющие строить аппроксимацию уравнения переноса излучения в пределах одной ячейки. Это значительно упрощает ситуацию около внешних границ и контактных разрывов решения. При этом желательно стоить схемы более высокого порядка, чем первый–второй. Схемы, которые имеют порядок аппроксимации более высокий, чем это позволяет тривиальный анализ размеров шаблона, называются компактными; схемы, в которых высокий порядок аппроксимации достигается на двухточечном шаблоне по каждой переменной, называются бикомпактными. В последнее десятилетие был предложен класс бикомпактных схем на основе использования метода прямых [1]–[8]. Пространственные производные дискретизуются, временные производные оставляются в дифференциальной форме. Полученная полудискретная бикомпактная схема может интегрироваться по времени любым методом, чаще всего используются диагонально-неявные методы Рунге–Кутты, например, третьего порядка аппроксимации. Для увеличения порядка аппроксимации по пространству используется расширение списка переменных и включение в него либо дополнительных переменных, например, интегрального среднего по ячейке, либо введение дополнительных внутренних точек шаблона. Схема консервативна, позволяет считать задачи методом бегущего счета или итерируемой факторизации в многомерных геометриях [9]. Метод не использует характеристические свойства уравнения переноса. Этот метод активно развивается, и на его основе получены прекрасные результаты в приложениях к задачам газовой динамики. С нашей точки зрения, этот метод обладает двумя недостатками. Первый – он работает на регулярных прямоугольных пространственных сетках, которые могут быть не слишком удобны при моделировании геометрии реальных технических устройств. Второй связан с особенностями решения уравнения переноса излучений и/или частиц. В этом случае уравнение переноса содержит недифференциальный член поглощения, так что точное решение уравнения вдоль каждой характеристики содержит экспоненциальный член. Развиваемый в настоящее время метод лебеговского осреднения уравнения переноса по энергии [10]–[15] обладает той особенностью, что сохраняет полный диапазон изменения коэффициента поглощения, который в центрах и в крыльях линий может отличаться на несколько порядков. Это означает, что при любом разумном выборе пространственной сетки в каких-то частях спектра ячейки будут оптически толстыми. Поэтому для нас важно, какая именно пространственная аппроксимация достигается в бикомпактных схемах. В [6] этот вопрос был исследован и показано, что пространственная аппроксимация бикомпактных схем четвертого порядка может быть получена для модельного уравнения Далквиста методом коллокации с функцией устойчивости

$R(z) = \frac{{1 + z{\text{/}}2 + {{z}^{2}}{\text{/}}12}}{{1 - z{\text{/}}2 + {{z}^{2}}{\text{/}}12}},$

где z – безразмерная оптическая толщина ячейки. Это означает, что данный способ пространственной аппроксимации не обладает “L-устойчивостью”, т.е. плохо передает экспоненциально затухающие решения в ситуации большой оптической толщины ячеек. Эти причины побудили вернуться к интерполяционно-характеристическим (или сеточно-характеристическим) схемам, в которых экспоненциальные зависимости решения могут быть частично учтены явно. Порядки аппроксимации по пространству и времени в характеристических схемах жестко связаны между собой и определяются точностью построения интерполянта для вычисления значения искомой величины в точке пересечения характеристики, выпущенной назад, до пересечения с границами ячейки, на которых решение известно. В [16] была предложена модификация известного CIP-метода. CIP (Cubic-Interpolated Pseudo-particle) метод и его модификация основаны на использовании эрмитовой интерполяции, требующей знания не только узловых значений переменных, но и узловых значений пространственной производной. Для данной схемы нельзя исследовать отдельно пространственную аппроксимацию и найти ее “функцию устойчивости” для уравнения Далквиста. Модификация CIP-метода основана на способе вычисления пространственных производных на новом временном шаге (т.е. на способе замыкания метода). В классическом CIP-методе вычисление производных на новом временном слое производится путем решения еще одного уравнения переноса, записанного для пространственных производных искомой функции (так называемых продолженных уравнений).

В данной работе исследуются диссипативно-дисперсионные свойства, основанные на эрмитовой интерполяции, и сравниваются с аналогичными свойствами бикомпактных схем, а также консервативной схемы второго порядка аппроксимации [17]. В последней схеме также используются переменные двух различных типов: узловые и консервативные. Это является общим для всех трех схем сравнения.

2. ПОСТРОЕНИЕ МОДИФИКАЦИИ CIP-СХЕМЫ

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

Пусть задано линейное однородное уравнение переноса (уравнение адвекции)

(1)

$Lu \equiv \frac{{\partial u(x,t)}}{{\partial t}} + a\frac{{\partial u(x,t)}}{{\partial x}} = 0.$

Дополним уравнение (1) начальными условиями

(2)

$u(x,0) = {{u}_{0}}(x),\quad 0 \leqslant x \leqslant X.$

В предположении положительности скорости переноса a краевые условия поставим на левой границе

(3)

$u(0,t) = \varphi (t),\quad 0 \leqslant t \leqslant T.$

Интерполяционно-характеристическая схема будет построена в рамках одной расчетной ячейки. Будем использовать аппроксимацию уравнения в рамках квадратной ячейки с вершинами (x_m, tⁿ), (${{x}_{m}}_{{ + 1}}$, tⁿ), (x_m, ${{t}^{n}}^{{ + 1}}$), (x_m+1, ${{t}^{n}}^{{ + 1}}$).

Для положительной скорости переноса a неизвестным сеточным значением искомой функции в каждой ячейке будет $y_{{m + 1}}^{{n + 1}}$. Выпустим из точки (${{x}_{m}}_{{ + 1}}$, ${{t}^{n}}^{{ + 1}}$) характеристику $x - at = \operatorname{const} $ назад до пересечения либо с нижней, либо с боковой гранью ячейки (см. фиг. 1). В первом случае число Куранта $\sigma = a\tau {\text{/}}h \leqslant 1$, во втором $\sigma \geqslant 1$. Обозначим координату пересечения обратной характеристики с границами ячейки либо x* ($\sigma \leqslant 1$), либо t* ($\sigma \geqslant 1$).

Фиг. 1.

Расчетная ячейка.

Для точного решения (1) значение из этой точки переносится без изменений в точку (${{x}_{m}}_{{ + 1}}$, ${{t}^{n}}^{{ + 1}}$). Таким образом, точность метода определяется точностью восстановления неузлового значения y*. В сеточно-характеристических методах для восстановления неузлового значения используется та или иная интерполяция.

Будем использовать кубическую интерполяцию Эрмита. Для случая $\sigma \leqslant 1$ пересечения с нижней гранью ячейки она в барицентрических пространственных координатах p и q на отрезке $[{{x}_{m}},{{x}_{{m + 1}}}]$ имеет вид

(4)

${{P}_{3}}(p,q) = {{H}^{R}}(p,q)y_{{m + 1}}^{n} + {{H}^{L}}(p,q)y_{m}^{n} + {{G}^{R}}(p,q)d_{{m + 1}}^{n}h + {{G}^{L}}(p,q)d_{m}^{n}h.$

Здесь $H$, $G$ – базисные функции Эрмита:

(5)

$\begin{gathered} p = (x{\text{*}} - {{x}_{m}}){\text{/}}h,\quad q = ({{x}_{{m + 1}}} - x*){\text{/}}h,\quad p + q = 1,\quad h = {{x}_{{m + 1}}} - {{x}_{m}}, \\ {{H}^{R}} = p(p + 2qp),\quad {{G}^{R}} = - p \cdot qp, \\ {{H}^{L}} = q(q + 2qp),\quad {{G}^{L}} = q \cdot qp, \\ \end{gathered} $

$y_{m}^{n}$ – сеточные значения искомой функции, а $d_{m}^{n}$ – сеточные значения пространственных производных искомой функции. Для случая $\sigma \geqslant 1$ точка начала характеристики, приходящей в точку (${{x}_{m}}_{{ + 1}}$, ${{t}^{n}}^{{ + 1}}$), попадает на левую границу, интерполяция на этой границе расчетной ячейки рассматривается аналогично с учетом производных по времени для искомой функции. Если в массиве данных хранятся только пространственные производные, производные по времени могут быть рассчитаны из уравнения переноса (1).

Как уже было указано, интерполяция вида (4), (5) была использована во многих работах (см., например, [18]–[24] и литературу в них). Отличие предлагаемого подхода от упомянутых работ заключается в способе определения величин $d_{m}^{{n + 1}}$, необходимых для возможности расчета величин на следующем временном слое. Как уже было сказано, в методе CIP стандартной процедурой является использование дифференциального продолжения уравнения (1), т.е. для пространственной производной также выписывается дифференциальное уравнение вида (1), которое должно решаться совместно с исходным уравнением переноса:

(1)'

$\frac{\partial }{{\partial x}}Lu \equiv \frac{\partial }{{\partial t}}\frac{{\partial u(x,t)}}{{\partial x}} + a\frac{\partial }{{\partial x}}\frac{{\partial u(x,t)}}{{\partial x}} = 0.$

В качестве начальных данных для вычисления $d_{{m + 1}}^{{n + 1}}$ в точке x* берется значение пространственной производной интерполянта (4). С нашей точки зрения, в приложении к задачам переноса излучений и частиц, у этого подхода есть два недостатка: во-первых, уравнение (1) допускает разрывные решения, и в этом случае производная уже должна рассматриваться как обобщенная функция, и, во-вторых, при рассмотрении неоднородного уравнения переноса с переменным коэффициентом поглощения продолженное уравнение будет содержать не только производную, но и саму функцию, что затрудняет построение схемы.

В предлагаемой модификации CIP-метода не используется дифференциальное продолжение уравнения (1'). Предположим, что в начальный момент времени мы знаем не только узловые значения функции ${{u}_{0}}({{x}_{m}})$, но и значения ее производной:

(6)

$d_{m}^{0} = u_{0}^{'}({{x}_{m}}).$

Для замыкания алгоритма нам необходима процедура получения сеточных значений $d_{{m + 1}}^{{n + 1}}$. Для граничного узла значение производной в случае краевых условий (3) может быть рассчитано из (1), (3) как $d_{0}^{{n + 1}} = - \dot {\varphi }({{t}^{{n + 1}}}){\text{/}}a$. Для однородного уравнения переноса производная $d_{{m + 1}}^{{n + 1}}$ может быть вычислена как производная интерполянта в точке x*: $d_{{m + 1}}^{{n + 1}} = P_{3}^{'}(p,q)$. Однако такой вариант вычисления производной не обобщается на неоднородное уравнение. Мы поступим иначе. Вычислим интегральное среднее на нижнем отрезке $[x*,{{x}_{{m + 1}}}]$ либо точно от интерполянта P₃, либо по формуле Симпсона. В силу характеристических свойств однородного уравнения (1) это интегральное среднее будет совпадать с интегральным средним по отрезку $[{{t}^{n}},{{t}^{{n + 1}}}]$ в точке ${{x}_{m}}_{{ + 1}}$ справа. Значение пространственной производной $d_{{m + 1}}^{n}$ может быть из уравнения (1) пересчитано в производную по времени: $g_{{m + 1}}^{n} = - ad_{{m + 1}}^{n}$. Тогда на правой границе нам известны два узловых значения $y_{{m + 1}}^{n}$, $y_{{m + 1}}^{{n + 1}}$, интегральное среднее ${{\bar {y}}_{{m + 1}}}$ и значение временной производной искомой функции $g_{{m + 1}}^{n}$. Эти данные позволяют получить значение $g_{{m + 1}}^{{n + 1}}$ из формулы Эйлера–Маклорена:

(7)

${{\bar {y}}_{{m + 1}}} \equiv \frac{1}{\tau }\int\limits_{{{t}^{n}}}^{{{t}^{{n + 1}}}} {u({{x}_{{m + 1}}},t)dt} = \frac{1}{2}(y_{{m + 1}}^{n} + y_{{m + 1}}^{{n + 1}}) - \frac{\tau }{{12}}(g_{{m + 1}}^{{n + 1}} - g_{{m + 1}}^{n}) + \frac{1}{{720}}u_{{tttt}}^{{(4)}}{{\tau }^{4}}.$

В свою очередь, величина $g_{{m + 1}}^{{n + 1}}$ может быть пересчитана в пространственную производную $d_{{m + 1}}^{{n + 1}} = - g_{{m + 1}}^{{n + 1}}{\text{/}}a$ из уравнения (1). Восстановление производных, исходя из сохранения интегральных средних, делает модификацию метода CIP консервативным для достаточно гладких решений (для разрывного начального профиля схема не является консервативной), а также позволяет обобщить схему на случай неоднородного уравнения переноса. Использование формулы Эйлера–Маклорена для замыкания системы уравнений идеологически очень близко к бикомпактным схемам Рогова.

В дальнейшем при ссылках на эту схему мы будем использовать обозначение CIP(3,3) как указание на порядки аппроксимации по времени и пространству.

3. КРАТКИЙ ВЫВОД ПОЛУДИСКРЕТНОЙ БИКОМПАКТНОЙ СХЕМЫ РОГОВА

Бикомпактная схема Рогова строится методом прямых. Полудискретная форма уравнений получается смешанным FD–FV-методом. Для повышения порядка аппроксимации рассматриваются не только узловые значения переменных, но и интегральные средние величины по пространственной ячейке. Интегрированием уравнения (1) по пространственной ячейке получается первое уравнение полудискретной формы схемы:

(8)

$\frac{\partial }{{\partial t}}{{\bar {y}}_{m}} + \frac{a}{h}({{y}_{{m + 1}}} - {{y}_{m}}) = 0.$

Интегрированием уравнения (1)' с использованием замыкающего соотношения формулы Эйлера–Маклорена (7) можно получить второе уравнение полудискретной системы:

(9)

$\frac{\partial }{{\partial t}}({{y}_{{m + 1}}} - {{y}_{m}}) + \frac{{6a}}{h}({{y}_{{m + 1}}} - 2{{\bar {y}}_{m}} + {{y}_{m}}) = 0.$

Система уравнений (8), (9) может быть проинтегрирована по времени любым из известных методов, обычно используются диагонально-неявные методы Рунге–Кутты третьего порядка аппроксимации. Методы третьего порядка аппроксимации сложны для диссипативно-дисперсионного анализа. Б.В. Рогов в [8] провел анализ полностью дискретных схем для неявной схемы Эйлера по времени и метода трапеций, который он называет методом типа Кранка–Николсон. Мы будем сравнивать диссипативно-дисперсионные свойства с реализацией метода трапеций по времени, а также с идеальными результатами полудискретной схемы (8), (9) .

Для реализации метода трапеций будем использовать ссылку на схему как BiC(4, 2) как указание на четвертый порядок аппроксимации по пространству и второй по времени, а на полудискретную схему как BiC(4, ∞).

4. КРАТКИЙ ВЫВОД КОНСЕРВАТИВНОЙ СХЕМЫ ГОЛОВИЗНИНА–ЧЕТВЕРУШКИНА

Консервативная схема Головизнина–Четверушкина также использует переменные двух видов: характеристические и консервативные. Узловые характеристические переменные вычисляются интерполяционно-характеристическим методом второго порядка (по схеме Бима–Уорминга) с использованием значения в полуцелом пространственном узле. Для обеспечения консервативности схемы введены консервативные переменные, отвечающие серединам сторон ячеек. На нижней стороне ячейки по известным значениям $y_{m}^{n}$, $y_{{m + 1/2}}^{n}$, $y_{{m + 1}}^{n}$ в трех точках строится интерполяция второго порядка, значение которой в точке пересечения с обратной характеристикой определяет $y_{{m + 1}}^{{n + 1}}$. Таким образом, с учетом $\sigma = a\tau {\text{/}}h$ получаем для характеристических переменных

(10)

$y_{{m + 1}}^{{n + 1}} = y_{{m + 1}}^{n} - 2\sigma (1.5y_{{m + 1}}^{n} - 2y_{{m + 1/2}}^{n} + 0.5y_{m}^{n}) + 2{{\sigma }^{2}}(y_{{m + 1}}^{n} - 2y_{{m + 1/2}}^{n} + y_{m}^{n}).$

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

$y_{{m + 1}}^{{n + 1/2}} = {{(y_{{m + 1}}^{n} + y_{{m + 1}}^{{n + 1}})} \mathord{\left/ {\vphantom {{(y_{{m + 1}}^{n} + y_{{m + 1}}^{{n + 1}})} 2}} \right. \kern-0em} 2},\quad y_{m}^{{n + 1/2}} = {{(y_{m}^{n} + y_{m}^{{n + 1}})} \mathord{\left/ {\vphantom {{(y_{m}^{n} + y_{m}^{{n + 1}})} 2}} \right. \kern-0em} 2}.$

Для расчета значений консервативных переменных в полуцелых узлах по пространству используется схема «крест» для уравнения переноса

(11)

$\frac{{y_{{m + 1/2}}^{{n + 1}} - y_{{m + 1/2}}^{n}}}{\tau } + a\frac{{y_{{m + 1}}^{{n + 1/2}} - y_{m}^{{n + 1/2}}}}{h} = 0.$

Исключая промежуточные величины, для консервативных переменных получаем

(12)

$y_{{m + 1/2}}^{{n + 1}} = y_{{m + 1/2}}^{n} - 0.5\sigma (y_{{m + 1}}^{n} + y_{{m + 1}}^{{n + 1}} - y_{m}^{n} - y_{m}^{{n + 1}}).$

Ссылки на схему (10), (12) будем оформлять в виде GC(2,2) как указание на порядки аппроксимации по пространству и по времени.

5. ДИССИПАТИВНО-ДИСПЕРСИОННЫЙ АНАЛИЗ СХЕМЫ СIP(3,3) ПРИ ЧИСЛАХ КУРАНТА, МЕНЬШИХ ЕДИНИЦЫ

В силу того, что и формула Симпсона, и формула Эйлера–Маклорена точны на многочленах третьей степени, то нетрудно показать, что для уравнения адвекции предлагаемый в разд. 2 способ вычисления пространственной производной на новом слое по времени эквивалентен переносу ее значения, получаемого дифференцированием интерполянта в точке x*, вдоль характеристики. Этот же способ вычисления производной для уравнения адвекции применяется для анализа и в классической схеме CIP (см. [21]). Проведем упрощенный (по сравнению с [21]) фурье-анализ полученной разностной схемы на диссипативно-дисперсионные свойства. Будем искать точное решение дифференциального уравнения адвекции в виде

(13)

$u(x,t) = {{e}^{{\lambda t}}}{{e}^{{ikx}}}.$

Подстановка фурье-гармоники в уравнение (1) для точного решения задачи дает связь

(14)

$\lambda = - ika.$

Каждая фурье-гармоника движется с одной и той же скоростью a, не изменяясь по амплитуде (λ – чисто комплексное). Для численного решения используем аналог (13) для узловых значений функции и узловых значений производной

(15)

$y_{m}^{n} = {{e}^{{\lambda n\tau }}}{{e}^{{ikmh}}},\quad d_{m}^{n} = \beta {{e}^{{\lambda n\tau }}}{{e}^{{ikmh}}}.$

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

(16)

$\begin{gathered} \frac{{\partial {{H}^{R}}}}{{\partial x}} = \frac{1}{h}\left( {\frac{{\partial {{H}^{R}}}}{{\partial p}} - \frac{{\partial {{H}^{R}}}}{{\partial q}}} \right) = \frac{{6pq}}{h},\quad \frac{{\partial {{H}^{L}}}}{{\partial x}} = \frac{1}{h}\left( {\frac{{\partial {{H}^{L}}}}{{\partial p}} - \frac{{\partial {{H}^{L}}}}{{\partial q}}} \right) = - \frac{{6pq}}{h}, \\ \frac{{\partial {{G}^{R}}}}{{\partial x}} = \frac{{p(1 - 3q)}}{h},\quad \frac{{\partial {{G}^{L}}}}{{\partial x}} = \frac{{q(1 - 3p)}}{h}. \\ \end{gathered} $

Точка пересечения характеристики с нижним слоем x* соответствует значению барицентрических координат $q = \sigma $, $p = 1 - \sigma $. Перенос значений функции и производной из этой точки в узел $({{t}^{{n + 1}}},{{x}_{{m + 1}}})$ эквивалентно системе

(17)

$\begin{gathered} y_{{m + 1}}^{{n + 1}} = p(p + 2pq)y_{{m + 1}}^{n} + q(q + 2pq)y_{m}^{n} - {{p}^{2}}qh \cdot d_{{m + 1}}^{n} + p{{q}^{2}}h \cdot d_{m}^{n}, \\ hd_{{m + 1}}^{{n + 1}} = 6pq \cdot y_{{m + 1}}^{n} - 6pq \cdot y_{m}^{n} + p(1 - 3q)h \cdot d_{{m + 1}}^{n} + q(1 - 3p)h \cdot d_{m}^{n}. \\ \end{gathered} $

Подстановка в эту систему гармоники (15) приводит к системе для коэффициента $\beta $:

(18)

$\begin{array}{*{20}{c}} {{{e}^{{\lambda \tau }}} = p(p + 2pq) + q(q + 2pq){{e}^{{ - ikh}}} - {{p}^{2}}qh\beta + p{{q}^{2}}h\beta {{e}^{{ - ikh}}},} \\ {h\beta {{e}^{{\lambda \tau }}} = 6pq - 6pq{{e}^{{ - ikh}}} + p(1 - 3q)h\beta + q(1 - 3p)h\beta {{e}^{{ - ikh}}},} \end{array}$

или

(19)

$\begin{array}{*{20}{c}} {h\beta [pq( - p + q{{e}^{{ - ikh}}})] + [p(p + 2pq) + q(q + 2pq){{e}^{{ - ikh}}} - {{e}^{{\lambda \tau }}}] = 0,} \\ {h\beta [p(1 - 3q) + q(1 - 3p){{e}^{{ - ikh}}} - {{e}^{{\lambda \tau }}}] + [6pq - 6pq{{e}^{{ - ikh}}}] = 0.} \end{array}$

Условие совместности системы приводит к квадратному уравнению относительно ${{e}^{{\lambda \tau }}}$:

(20)

$\begin{gathered} {{({{e}^{{\lambda \tau }}})}^{2}} - {{e}^{{\lambda \tau }}}[p(p + 2pq) + q(q + 2pq){{e}^{{ - ikh}}} + p(1 - 3q) + q(1 - 3p){{e}^{{ - ikh}}}] + \\ + \;[p(p + 2pq) + q(q + 2pq){{e}^{{ - ikh}}}][p(1 - 3q) + q(1 - 3p){{e}^{{ - ikh}}}] - 6{{p}^{2}}{{q}^{2}}[1 - {{e}^{{ - ikh}}}][p - q{{e}^{{ - ikh}}}] = 0. \\ \end{gathered} $

Это уравнение имеет два корня, корень, отвечающий дисперсионному соотношению для данной схемы, соответствует знаку “плюс” перед корнем из детерминанта. При $q = \sigma $, $p = 1 - \sigma $:

(21)

$\begin{gathered} \lambda = \frac{1}{\tau }\ln \left( {1 - 2\sigma + {{\sigma }^{3}} + \sigma {{e}^{{ - ikh}}}( - 1 + 3\sigma - {{\sigma }^{2}}) + {{e}^{{ - ikh}}}\sigma (1 - \sigma ) \times } \right. \\ \times \;\left. {\sqrt {{{\sigma }^{2}} - 4\sigma + 1 - 2{{e}^{{ikh}}}({{\sigma }^{2}} - \sigma - 5) + {{e}^{{2ikh}}}({{\sigma }^{2}} + 2\sigma - 2)} } \right). \\ \end{gathered} $

Разложение (21) в ряд Тейлора для длинноволновых гармоник $kh \ll 1$ дает

(22)

$\begin{gathered} \lambda = - iak + \frac{1}{{72}}ak(\sigma - 1)({{\sigma }^{2}} - \sigma + 1){{(kh)}^{3}} + \frac{i}{{540}}ak({{\sigma }^{2}} - 1)(2\sigma - 1)(\sigma - 2){{(kh)}^{4}} - \\ - \;\frac{1}{{648}}ak(\sigma - 1){{({{\sigma }^{2}} - \sigma + 1)}^{2}}{{(kh)}^{5}} + O)({{(kh)}^{6}}). \\ \end{gathered} $

Выражение (22) показывает, что эрмитова характеристическая схема действительно обладает третьим порядком аппроксимации, главный член погрешности по сравнению с (10) – диссипативный. На фиг. 2 представлена величина коэффициента при ${{(kh)}^{3}}$. Отрицательность коэффициента означает, что схема диссипативна (и, следовательно, устойчива).

Фиг. 2.

Коэффициент диссипации схемы в зависимости от числа Куранта при малых значениях безразмерного волнового числа kh.

Коэффициент при дисперсионном члене разложения представлен на фиг. 3. Для длинноволновых гармоник фазовая скорость переноса:

(23)

$a{\kern 1pt} * \approx a(1 - ({{\sigma }^{2}} - 1)(2\sigma - 1)(\sigma - 2){{(kh)}^{4}}{\text{/}}540).$

Там, где коэффициент отрицателен (для $0 < \sigma < 0.5$), дисперсия опережающая, так как эффективная скорость выше, чем a. И наоборот, там, где коэффициент положителен, дисперсия запаздывающая (для $0.5 < \sigma < 1$).

Фиг. 3.

Коэффициент при дисперсионном члене разложения (22) при различных числах Куранта при малых значениях безразмерного волнового числа kh.

6. СРАВНЕНИЕ ДИССИПАТИВНЫХ СВОЙСТВ РАЗНОСТНЫХ СХЕМ

На фиг. 4 представлены значения амплитудного множителя CIP(3,3)-схемы при значениях безразмерной величины $\phi = kh$ от 0 до $\pi $ и числе Куранта от 0 до 1. Паразитный корень приведен только для того, чтобы убедиться в устойчивости схемы.

Фиг. 4.

Временной множитель амплитуды $\left| {{{e}^{{\lambda \tau }}}} \right|$ в зависимости от числа Куранта $0 \leqslant \sigma \leqslant 1$ и безразмерного волнового числа $0 \leqslant kh \leqslant \pi $ для схемы CIP(3, 3). (a) Физический корень (шаг линий уровня 0. 01), (б) паразитный корень (шаг линий уровня 0.1).

Cхемы Рогова BIC(4,2) и BIC(4,∞) бездиссипативны в силу полной симметричности схемы. Поэтому рисунки $\left| {{{e}^{{\lambda \tau }}}} \right| = 1$ не приводятся. Для схемы Головизнина–Четверушкина GC(2,2) фурье-анализ дает дисперсионное соотношение

$\begin{gathered} {{e}^{{\lambda \tau }}} = \frac{1}{2}\left( {2 - 3\sigma + 2{{\sigma }^{3}} - \sigma \left( {{{e}^{{ - i\phi }}}(1 - 4\sigma + 2{{\sigma }^{2}})} \right.} \right. - [{{e}^{{ - i2\phi }}}{{(1 - 4\sigma + 2{{\sigma }^{2}})}^{2}} + \\ \left. {\left. { + \;{{e}^{{ - i\phi }}}(22 - 40\sigma + 8{{\sigma }^{2}} + 16{{\sigma }^{3}} - 8{{\sigma }^{4}}) - 7 + 16\sigma - 12{{\sigma }^{2}} + 4{{\sigma }^{4}}{{]}^{{0.5}}}} \right)} \right). \\ \end{gathered} $

На фиг. 5 представлены значения амплитудного множителя схемы GC(2,2) при значениях безразмерной величины $\phi = kh$ от 0 до $\pi $ и числе Куранта от 0 до 1.

Фиг. 5.

Временной множитель амплитуды $\left| {{{e}^{{\lambda \tau }}}} \right|$ в зависимости от числа Куранта $0 \leqslant \sigma \leqslant 1$ и безразмерного волнового числа $0 \leqslant kh \leqslant \pi $ для схемы GC(2, 2). (a) Физический корень (шаг линий уровня 0.05), (б) паразитный корень (шаг линий уровня 0.05).

Сравнение фиг. 4 и 5 показывает, что схема CIP(3,3) обладает достаточно малой диссипацией по сравнению со схемой GC(2,2). Эти схемы обладают разным порядком аппроксимации. Однако заметим, что схемы четного порядка аппроксимации обладают преимущественной дисперсионной ошибкой. Диссипативная ошибка появляется в следующем члене разложения $\lambda $ по степеням $\phi $. Схемы нечетного порядка аппроксимации обладают преимущественной диссипативной ошибкой. Поэтому сравнение схемы CIP(3,3) со схемой GC(2,2) по диссипации происходит в одном и том же порядке по $\phi $ (в третьем). Нетрудно видеть, что эрмитова схема обладает значительно меньшей диссипацией.

7. СРАВНЕНИЕ ДИСПЕРСИОННЫХ СВОЙСТВ КОНСЕРВАТИВНЫХ РАЗНОСТНЫХ СХЕМ

Бикомпактные схемы Рогова обладают численной дисперсией, поэтому приведем дисперсионные соотношения для полудискретной схемы BiC(4, ∞):

(24)

$\lambda = \frac{{a(3 + 3{{e}^{{i\phi }}} - {{{[ - 3 + 42{{e}^{{i\phi }}} - 3{{e}^{{i2\phi }}}]}}^{{0.5}}})}}{{h(1 - {{e}^{{i\phi }}})}},$

и схемы метода трапеций BiC(4, 2):

(25)

${{e}^{{\lambda \tau }}} = \frac{{\sin (\phi {\text{/}}2)(2 - 6{{\sigma }^{2}}) + i\sigma {{{[42 - 6\cos \phi ]}}^{{0.5}}}}}{{\sin (\phi {\text{/}}2)(2 + 6{{\sigma }^{2}}) - i6\sigma \cos (\phi {\text{/}}2)}}.$

Заметим, что для полудискретной бикомпактной схемы (24) у временного показателя пропадает зависимость от числа Куранта. На фиг. 6 приведены зависимости отношения численной скорости переноса гармоник к их точному значению $a{\text{*/}}a$ в зависимости от числа Куранта $\sigma $ и безразмерного волнового числа $\phi $.

Фиг. 6.

Зависимость отношения $a{\text{*/}}a$ скорости переноса гармоники с безразмерным волновым числом $0 \leqslant \phi \leqslant \pi $ к точной скорости переноса, равной a, при различных числах Куранта для эрмитовой схемы CIP(3, 3), бикомпактных схем Рогова BiC(4, ∞) и BiC(4, 2), и схемы Головизнина–Четверушкина GC(2, 2). (а) CIP(3, 3) (шаг линий уровня 0.005), (б) BiC(4, ∞), (в) BiC(4, 2) (шаг линий уровня 0.03), (г) GC(2, 2) (шаг линий уровня 0.05). Изменение направления оси по числу Куранта на Фиг. 6. сделано для лучшей видимости поведения функций.

Главные дисперсионные члены эрмитовой схемы и схемы Рогова имеют четвертый порядок аппроксимации по $\phi $, а схема Головизнина–Четверушкина – второй. Поэтому наиболее интересно сравнение по свойствам со схемой своего класса: с полудискретной схемой Рогова BiC(4, ∞). Эрмитова схема CIP(3,3) показывает меньшую дисперсию фурье-гармоник по сравнению с полудискретной бикомпактной схемой Рогова. Если отношение $a{\kern 1pt} *{\kern 1pt} {\text{/}}a$ для эрмитовой схемы для коротковолновых возмущений изменяется в пределах от 0.98 до 1.06, то для полудискретной бикомпактной схемы диапазон изменения от 1 до 1.1. Для полностью дискретных схем Рогова дисперсия может быть больше.

8. ДИССИПАТИВНО-ДИСПЕРСИОННЫЙ АНАЛИЗ СХЕМЫ ПРИ ЧИСЛАХ КУРАНТА, БОЛЬШИХ ЕДИНИЦЫ

С математической точки зрения в уравнении адвекции при постоянной скорости переноса переменные $x$ и $t$ равноправны. На этом можно было бы поставить точку, так как формально мы можем просто поменять местами переменные. Однако проведем полный анализ диссипативно-дисперсионных свойств и в этом случае. Интерполяция будет вестись по левому ребру ячейки

(26)

${{P}_{3}}(p,q) = {{H}^{U}}(p,q)y_{m}^{{n + 1}} + {{H}^{D}}(p,q)y_{m}^{n} + {{G}^{U}}(p,q)g_{m}^{{n + 1}}\tau + {{G}^{D}}(p,q)g_{m}^{n}\tau .$

Здесь $H$, $G$ – базисные функции Эрмита:

$\begin{gathered} p = (t{\text{*}} - {{t}^{n}}){\text{/}}\tau ,\quad q = ({{t}^{{n + 1}}} - t*){\text{/}}\tau ,\quad p + q = 1,\quad \tau = {{t}^{{n + 1}}} - {{t}^{n}}, \\ {{H}^{U}} = p(p + 2qp),\quad {{G}^{U}} = - p \cdot qp, \\ {{H}^{D}} = q(q + 2qp),\quad {{G}^{D}} = q \cdot qp. \\ \end{gathered} $(27)

В этом случае схема для уравнения адвекции также будет эквивалентна переносу из точки t* не только значения функции, но и значения ее временной производной. Аналогичный фурье-анализ типа (18) для функции и ее производной по времени

$\begin{gathered} {{e}^{{ikh}}} = p(p + 2qp) + {{e}^{{ - \lambda \tau }}}q(q + 2qp) + \tau \beta ( - {{p}^{2}}q) + \tau \beta {{e}^{{ - \lambda \tau }}}p{{q}^{2}}, \\ \tau \beta {{e}^{{ikh}}} = 6pq + {{e}^{{ - \lambda \tau }}}( - 6pq) + \tau \beta \cdot p(1 - 3q) + \tau \beta {{e}^{{ - \lambda \tau }}}q(1 - 3p) \\ \end{gathered} $

с учетом связи $p + q = 1$ приводит к квадратному уравнению

$\begin{gathered} ({{e}^{{2ikh}}} + {{(1 - q)}^{4}} - 2{{e}^{{ikh}}}(1 - 2q + {{q}^{3}})){{({{e}^{{\lambda \tau }}})}^{2}} - \\ - \;2q[1 - 2{{q}^{2}} + {{q}^{3}} - {{e}^{{ikh}}}(1 - 3q + {{q}^{2}})]({{e}^{{\lambda \tau }}}) + {{q}^{4}} = 0. \\ \end{gathered} $

Точка пересечения t* с левым ребром ячейки будет иметь координаты:

$q = \frac{{h{\text{/}}a}}{\tau } = {{\sigma }^{{ - 1}}},$

тогда корни этого квадратного уравнения

(28)

$\begin{gathered} {{e}^{{\lambda \tau }}} = q\left[ {1 - 2{{q}^{2}} + {{q}^{3}} - {{e}^{{ikh}}}(1 - 3q + {{q}^{2}}) \pm (1 - q) \times } \right. \\ \left. { \times \;\sqrt {1 + 2q - 2{{q}^{2}} + 2{{e}^{{ikh}}}( - 1 + q + 5{{q}^{2}}) + {{e}^{{2ikh}}}(1 - 4q + {{q}^{2}})} } \right] \times \\ \times \;{{\left( {{{{(1 - q)}}^{4}} - 2{{e}^{{ikh}}}(1 - 2{{q}^{2}} + {{q}^{3}}) + {{e}^{{2ikh}}}} \right)}^{{ - 1}}}. \\ \end{gathered} $

Выраженный через число Куранта физически осмысленный корень соответствует дисперсионному соотношению со знаком “+” перед корнем. Разложение в ряд Тейлора для длинноволновых гармоник дает дисперсионное соотношение:

$\begin{gathered} \lambda = - iak - \frac{1}{{72}}ak(\sigma - 1)({{\sigma }^{2}} - \sigma + 1){{(kh)}^{3}} + \frac{i}{{540}}ak({{\sigma }^{2}} - 1)(2\sigma - 1)(\sigma - 2){{(kh)}^{4}} - \\ - \;\frac{1}{{648}}ak(\sigma - 1){{({{\sigma }^{2}} - \sigma + 1)}^{2}}{{(kh)}^{5}} + O({{(kh)}^{6}}). \\ \end{gathered} $(29)

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

Фиг. 7.

Временной множитель амплитуды $\left| {{{e}^{{\lambda \tau }}}} \right|$ в зависимости от обратного числа Куранта $0 \leqslant q = {{\sigma }^{{ - 1}}} \leqslant 1$ и безразмерного волнового числа $0 \leqslant kh \leqslant \pi $ для CIP(3, 3). (a) Физический корень (шаг линий уровня 0.0625), (б) паразитный корень (шаг линий уровня 0.0625).

На фиг. 8 изображена зависимость отношения $a{\text{*/}}a$ в зависимости от угла и обратного числа Куранта для схемы CIP(3, 3) и BiC(4, 2). Чтобы избежать перехода с ветви на ветвь для аргумента от выражения (28), выражение (28) умножается на –1, после чего из аргумента вычитается π.

Фиг. 8.

Зависимость отношения $a{\text{*/}}a$ скорости переноса гармоники с безразмерным волновым числом $0 \leqslant \phi \leqslant \pi $ к точной скорости переноса, равной a, в зависимости от обратного числа Куранта для эрмитовой характеристической схемы и дискретной бикомпактной схемы Рогова. (а) CIP(3, 3) (шаг линий уровня 0.0625), (б) BiC(4, 2) (шаг линий уровня 0.0625).

Для полудискретной бикомпактной схемы график для чисел Куранта остается таким же, как это представлено на фиг. 6б. Однако расхождение по мере увеличения числа Куранта между идеальным дисперсионным соотношением и реальным значительно увеличивается, и не в пользу реального. Фигура 8 показывает, что дисперсия метода CIP достаточно мала и при числах Куранта, больших единицы.

9. ВЫВОДЫ

Все рассмотренные схемы построены на минимальном двухточечном шаблоне. Все они показывают неплохие диссипативно-дисперсионные свойства. Схема Головизнина–Четверушкина проста в реализации, что делает ее очень привлекательной в глазах вычислителей. Бикомпактные схемы при симметричной временной аппроксимации обладают нулевой диссипацией, что бесценно при рассмотрении задач аэроакустики. Предлагаемая в работе [16] модификация CIP -метода допускает обобщение на неоднородное уравнение переноса и неструктурированные сетки. Использование характеристик сближает модификацию CIP-метода со схемой Головизнина–Четверушкина, а использование формулы Эйлера–Маклорена – с бикомпактной схемой Рогова. Эрмитова характеристическая схема обладает достаточно малой диссипацией и экстрамалой дисперсией. Последнее свойство важно в задачах переноса излучений и частиц, поскольку значительно уменьшает амплитуду паразитных колебаний на границах раздела сред.

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

Рогов Б.В., Михайловская М.Н. Бикомпактные схемы четвертого порядка аппроксимации для гиперболических уравнений // Докл АН. 2010. Т. 430. № 4. С. 470–474.
Рогов Б.В., Михайловская М.Н. Монотонные бикомпактные схемы для линейного уравнения переноса // Матем. моделирование. 2011. Т. 23. № 6. С. 98–110.
Рогов Б.В., Михайловская М.Н. Монотонные компактные схемы бегущего счета для систем уравнений гиперболического типа // Ж. вычисл. матем. и матем. физ. 2012. Т. 52. № 4. С. 672–695.
Аристова Е.Н., Рогов Б.В. О реализации граничных условий в бикомпактных схемах для линейного уравнения переноса // Матем. моделирование. 2012. Т. 24. № 10. С. 3–14.
Aristova E.H., Rogov B.V. Bicompact scheme for the multidimensional stationary linear transport equation // Appl. Numer. Math. 2015. V. 93. P. 3–14.
Аристова Е.Н., Рогов Б.В., Чикиткин А.В. Оптимальная монотонизация высокоточной бикомпактной схемы для нестационарного многомерного уравнения переноса // Ж. вычисл. матем. и матем. физ. 2016. Т. 56. № 6. С. 973–988.
Chikitkin A.V., Rogov B.V. Family of central bicompact schemes with spectral resolution property for hyperbolic equations // Applied Numerical Mathematics. 2019. V. 142. P. 151–170.
Rogov B.V. Dispersive and dissipative properties of the fully discrete bicompact schemes of the fourth order of spatial approximation for hyperbolic equations // Applied Numerical Mathematics. 2019. V. 139. P. 136–155.
Рогов Б.В., Чикиткин А.В. О сходимости и точности метода итерируемой приближенной факторизации операторов многомерных высокоточных бикомпактных схем // Матем. моделирование. 2019. Т. 31. № 12 .С. 119–144.
Цветкова И.Л., Шильков А.В. Осреднение уравнения переноса в резонансно поглощающей среде // Матем. моделирование. 1989. Т. 1. № 1. С. 91–100.
Шильков А.В. Методы осреднения сечений и энергетического спектра в задачах переноса нейтронов // Матем. моделирование. 1991. Т. 3. № 2. С. 63–81.
Shilkov A.V. Generalized Multigroup Approximation and Lebesgue Averaging Method in Particle Transport Problems // Transp. Theory and Stat. Physics. 1994. V. 23. № 6. P. 781–814.
Шильков А.В., Герцев М.Н. Верификация метода лебеговского осреднения // Матем. моделирование. 2015. Т. 27. № 8. С. 13–31.
Аристова Е.Н., Герцев М.Н., Шильков А.В. Метод лебеговского осреднения в серийных расчетах атмосферной радиации // Ж. вычисл. матем. и матем. физ. 2017. Т. 57. № 6. С. 1033–1047.
Шильков А.В. Метод лебеговых моментов для решения уравнения переноса нейтронов // Матем. моделирование. 2020. Т. 32. № 5. С. 59–94.
Аристова Е.Н., Овчаров Г.И. Эрмитова характеристическая схема для неоднородного линейного уравнения переноса // Матем. моделирование. 2020. Т. 32. № 3. С. 3–18.
Головизнин В.М., Четверушкин Б.Н. Алгоритмы нового поколения в вычислительной гидродинамике // Ж. вычисл. матем. и матем. физ. 2018. Т. 58. № 8. С. 20–29.
Yabe T., Aoki T., Sakaguchi G. et al. The compact CIP (cubic-interpolated pseudo-particle) method as a general hyperbolic solver // Computers & Fluids. 1991. V. 19. № 3/4. P. 421–431.
Tsai T.-L., Chiang S.-W., Yang J.-G. Characteristics method with cubic–spline interpolation for open channel flow computation // Int. J. Numerical Methods in Fluids. 2004. V. 46. P. 663–683.
Yabe T., Xiao F., Utsumi T. The constrained interpolation profile method for multiphase analysis // J. Comput. 2001. V. 169. № 2. P. 556–593.
Aoki T. Stability and accuracy of the cubic interpolated propagation scheme // Comput. Phys. Commun. 1997. V. 101. № 1–2. P. 9–20.
Голубев В.И., Петров И.Б., Хохлов Н.И. Компактные сеточно-характеристические схемы повышенного порядка точности для трехмерного линейного уравнения переноса // Матем. моделирование. 2016. Т. 28. № 2. С. 123–132.
Фаворская А.В., Петров И.Б. Численное моделирование волновых процессов в скальных массивах сеточно-характеристическим методом // Матем. моделирование. 2018. Т. 30. № 3. С. 37–51.
Петров И.Б., Фаворская А.В., Хохлов Н.И. Сеточно-характеристический метод на системах вложенных иерархических сеток и его применение для исследования сейсмических волн // Ж. вычисл. матем. и матем. физ. 2017. Т. 57. № 11. С. 1804–1811.

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

Инструменты

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

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

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