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

ВВЕДЕНИЕ

В последнее время существенное развитие получил класс высокоточных бикомпактных схем для уравнений гиперболического типа [1]–[8]. Эти схемы выводятся при помощи метода прямых и интегро-интерполяционного метода. В каждой ячейке искомое решение приближается финитным полиномом, построенным только по данным из этой ячейки, причем совокупность этих полиномов образует как минимум непрерывную функцию во всей расчетной области. В свою очередь, чтобы полином имел высокую степень (т.е. чтобы соответствующая бикомпактная схема была высокого порядка аппроксимации), вводятся дополнительные сеточные функции, определяемые либо в уже имеющихся целых узлах, либо в новых полуцелых узлах. По своему смыслу эти вспомогательные функции могут выражать саму искомую функцию, ее производные или первообразную. Для их отыскания используются дискретизированные дифференциальные следствия исходного уравнения в частных производных. С одной стороны, в результате перечисленных действий получается компактная конечно-разностная аппроксимация пространственных производных на шаблоне, который по каждому направлению включает не более двух целых узлов сетки, чем и объясняется название бикомпактных схем. С другой стороны, бикомпактные схемы в некотором понимании близки [5] к конечно-элементным схемам типа Галеркина, однако, схемы этих классов имеют разные, не тождественные уравнения.

Главным достоинством бикомпактных схем является сочетание следующих свойств: (а) спектральное разрешение, лучшее по сравнению с классическими компактными схемами равного порядка аппроксимации [4]; (б) совпадение числа граничных условий в дискретной (разностной) и дифференциальной задачах; (в) неявность (т.е. хорошая устойчивость); (г) низкая вычислительная трудоемкость, близкая к явным схемам. На примере практически важных задач о распространении детонационных волн в идеальном (невязком) газе показано преимущество бикомпактных схем в сравнении с другими широко применяющимися в этой области схемами [7], [8].

Хорошее спектральное разрешение бикомпактных схем делает желательным их применение для расчетов турбулентных течений вязкой жидкости. Учесть вязкость можно, применяя, например, расщепление по физическим процессам и далее какие-либо известные схемы для аппроксимации параболической части нестационарной задачи. Тем не менее более естественным было бы включить эту часть однородным образом, в рамках того же численного подхода, что применяется в бикомпактных схемах. Здесь возникает новая проблема: бикомпактные схемы разработаны лишь для одномерных уравнений параболического типа [9]. Целью настоящей работы является решение этой проблемы применительно к линейному многомерному уравнению конвекции-диффузии.

Приложение компактных схем, в том числе мультиоператорных, к уравнениям параболического типа подробно рассматривается в монографии [10], а конечно-элементных схем Галеркина (преимущественно локальных разрывных схем Галеркина, LDG) – в работах [11], [12] (см. также ссылки в этих трех работах). Для сравнения, бикомпактные схемы отличаются от первых схем меньшим пространственным шаблоном, а от последних – меньшим количеством дополнительных зависимых переменных и возможностью использовать неявную аппроксимацию по времени без существенных потерь в эффективности (экономичности) счета. Отметим дополнительно, что для преодоления высокой вычислительной трудоемкости неявных DG-схем сейчас используются два подхода, представляющие самостоятельный интерес: явно-неявная аппроксимация по времени [13], [14] и гиперболизация либо на уровне схемы [15], либо на уровне исходных дифференциальных уравнений [16], [17].

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

1. ПОЛУДИСКРЕТНЫЕ БИКОМПАКТНЫЕ СХЕМЫ

Рассмотрим сперва простейшее одномерное линейное уравнение диффузии:

где $u = u(x,t)$ – искомая скалярная функция, $\nu = {\text{const}} > 0$ – заданный коэффициент, символ ${{\partial }_{x}} \equiv \partial {\text{/}}\partial x$. Не конкретизируя пока постановку начальных и граничных условий, обсудим вопрос построения бикомпактных схем для уравнения (1) и его более сложных версий.

Введем на отрезке $x \in [0,{{x}_{{max}}}]$ (в общем случае, неравномерную) сетку

Ради краткости записи индекс $j + 1{\text{/}}2$ у шага ${{h}_{x}}$ в дальнейшем опускается всюду, где это допустимо. Наряду с функцией $u(x,t)$ мы будем использовать ее производную по $x$: Определим на $\Omega $ две сеточные функции: ${{u}_{j}}(t)$ и ${{v}_{{x,j}}}(t)$, аппроксимирующие в узлах $\Omega $ функции $u(x,t)$ и ${{v}_{x}}(x,t)$ соответственно.

Полудискретная бикомпактная схема четвертого порядка аппроксимации для уравнения (1) известна (см. [9]), она имеет следующий вид:

где индекс $j = \overline {0,{{N}_{x}} - 1} $. Схема (3) выводится методом прямых и интегро-интерполяционным методом. Уравнение (1) и соотношение (2) при фиксированном $t = {\text{const}}$ осредняются по ячейке ${{K}_{{j + 1/2}}} = [{{x}_{j}},{{x}_{{j + 1}}}]$ (т.е. интегрируются по $x$ от ${{x}_{j}}$ до ${{x}_{{j + 1}}}$ и делятся на ${{h}_{x}}$), после чего интегралы либо вычисляются точно по формуле Ньютона–Лейбница, либо приближенно по квадратурной формуле Эйлера–Маклорена с отбрасыванием остаточного члена $O(h_{x}^{4})$ (${{h}_{x}} \to 0$). При применении этой квадратуры к левой части соотношения (2) возникает производная ${{\partial }_{x}}{{v}_{x}}$, которая заменяется на ${{\nu }^{{ - 1}}}{{\partial }_{t}}u$ при помощи уравнения (1).

Обратим внимание на последнее действие, а именно, на замену производной $\partial {{v}_{x}}{\text{/}}\partial x$ при дискретизации (2) на ${{\nu }^{{ - 1}}}{{\partial }_{t}}u$. Ясно, что сделать аналогичное при переменном коэффициенте $\nu $ или для двумерного уравнения диффузии (даже линейного) весьма проблематично: в первом случае придется отказаться от консервативной аппроксимации и привлечь производную ${{\partial }_{x}}\nu $ (которая вообще может не существовать), а во втором случае помешает производная по $y$.

Найдем другой, более универсальный подход к конструированию бикомпактных схем для уравнения диффузии. Рассмотрим уравнение (1) с добавлением линейного конвективного члена и источника:

где ${{c}_{x}} = {\text{const}}$ – заданный коэффициент, $f(x,t)$ – заданная функция. Переход от уравнения диффузии (1) к уравнению конвекции-диффузии (4) обусловлен возможностью испытать получаемые бикомпактные схемы на более интересных вычислительных примерах (см. разд. 3).

Построим в ячейке ${{K}_{{j + 1/2}}}$ по значениям ${{u}_{j}}$, ${{u}_{{j + 1}}}$, ${{v}_{{x,j}}}$, ${{v}_{{x,j + 1}}}$ кубический интерполяционный полином Эрмита:

где ${{x}_{{j + 1/2}}} = ({{x}_{j}} + {{x}_{{j + 1}}}){\text{/}}2$, а коэффициенты ${{a}_{m}} = {{a}_{m}}(t)$, $m = \overline {0,3} $, даются выражениями: Если функция $u(x,t)$ достаточно гладкая, то полином (5) аппроксимирует ее в ячейке ${{K}_{{j + 1/2}}}$ с погрешностью $O(h_{x}^{4})$. Далее, возьмем уравнение (4) и его дифференциальное следствие (т.е. производную по $x$) и, вновь следуя методу прямых и интегро-интерполяционному методу, при фиксированном $t = {\text{const}}$ осредним эти уравнения по ячейке ${{K}_{{j + 1/2}}}$. Искомую функцию $u(x,t)$ при этом заменим на полином (5) с погрешностью $O(h_{x}^{4})$. Приведем выражения для интегралов, появляющихся в процессе осреднения уравнений (4), (6) и замены $u$ на ${{u}_{h}}$: Заметим, что формула для интеграла от ${{u}_{h}}$ повторяет квадратуру Эйлера–Маклорена, чего и следовало ожидать. В итоге после выполнения всех действий мы приходим к полудискретной бикомпактной схеме четвертого порядка аппроксимации для уравнения (4): где Нетрудно проверить, что схема (8) в точности переходит в схему (3) при ${{c}_{x}} = 0$, $f(x,t) \equiv 0$. Укажем на то, что вышеизложенный метод построения бикомпактных схем для уравнения теплопроводности по своей идее схож с аналогичным методом [18] построения схем этого типа уже для гиперболических уравнений.

Отметим, что формулы (7) и схема (8) могут быть записаны более компактно при помощи сеточных операторов. Определим такие операторы:

где ${{U}_{j}}$ – произвольная сеточная функция. Оператор ${{P}^{x}}$ заменяет значение сеточной функции в узле на значение соответствующей ей сеточной функции-производной по $x$ (в смысле интерполяции Эрмита) в этом же узле. Например, действие ${{P}^{x}}$ на ${{u}_{j}}$ дает ${{v}_{{x,j}}}$. Сразу уточним, что в настоящей работе этот оператор будет применяться не только к ${{u}_{j}}$, но и функциям типа производной от ${{u}_{j}}$ по времени, приращения ${{u}_{j}}$ за шаг по времени или за итерацию; в этом случае результат действия ${{P}^{x}}$ будет означать соответственно производную от ${{{v}}_{{x,j}}}$ по времени, приращение ${{v}_{{x,j}}}$ за шаг по времени или за итерацию. Итак, с применением операторов (9) формулы (7) приобретают вид а схема (8) имеет вид

Покажем теперь, каким именно образом новый метод построения бикомпактных схем является более универсальным по сравнению с [9].

Начнем с переменного коэффициента $\nu $. В этом случае диффузионный член в уравнении (4) примет вид:

Его осреднение по ячейке выполняется точно по формуле Ньютона–Лейбница Если коэффициент $\nu (x)$ терпит разрыв на одной из границ ячейки ${{K}_{{j + 1/2}}}$, его можно определить там, например, как полусумму пределов слева и справа. При дискретизации дифференциального следствия необходимо будет вычислить интегральное среднее уже от производной $\mathcal{D}u$ по $x$: Здесь функция ${{{v}}_{x}}$ заменяется на полином ${{{v}}_{{x,h}}}$ – производную по $x$ от полинома (5), а функция $\nu (x)$ – либо на (а) интерполяционный полином второй степени, построенный по данным в точках ${{x}_{j}}$, ${{x}_{{j + 1/2}}}$, ${{x}_{{j + 1}}}$, либо на (б) производную от кубического интерполяционного полинома Эрмита для первообразной $\nu (x)$. В любом случае после замены интеграл преобразуется в алгебраическое выражение, содержащее сеточные значения $\nu (x)$ и, быть может, значения ее первообразной (но не производной), а также значения сеточных функций ${{u}_{j}}$ и ${{{v}}_{{x,j}}}$.

Перейдем к обобщению полудискретной схемы (8) и ее частного случая (3) на многомерные варианты уравнения (4), что и составляет основную цель и предмет настоящей работы.

Рассмотрим двумерную версию (4):

где $u = u(x,y,t)$ – искомая функция, ${{c}_{x}} = {\text{const}}$, ${{c}_{y}} = {\text{const}}$, $\nu = {\text{const}} > 0$ – заданные коэффициенты, $f(x,y,t)$ – заданная функция.

По аналогии с сеткой ${{\Omega }_{x}}$ введем на отрезке $[0,{{y}_{{max}}}]$ сетку ${{\Omega }_{y}}$:

Таким образом, мы получаем сетку $\Omega = {{\Omega }_{x}} \times {{\Omega }_{y}}$ в замкнутой расчетной области $\bar {D}$. В отличие от одномерного случая, помимо функции $u(x,y,t)$ мы будем работать не с одной, а с тремя дополнительными функциями: Функциям $u(x,y,t)$, ${{{v}}_{x}}(x,y,t)$, ${{{v}}_{y}}(x,y,t)$, ${{{v}}_{{xy}}}(x,y,t)$ поставим в соответствие аппроксимирующие их сеточные функции ${{u}_{{j,k}}}(t)$, ${{{v}}_{{x,j,k}}}(t)$, ${{{v}}_{{y,j,k}}}(t)$, ${{{v}}_{{xy,j,k}}}(t)$, заданные на $\Omega $. Аналогично оператором (9) определим сеточные операторы с верхним индексом “$y$”. Операторы с верхним индексом “$x$” действуют только по первому индексу сеточной функции, а операторы с верхним индексом “$y$” – только по второму, например: где $C = (j + 1{\text{/}}2,k + 1{\text{/}}2)$ – мультииндекс. Поясним принцип действия операторов ${{P}^{x}}$, ${{P}^{y}}$: Ясно, что любой оператор с индексом “$x$” коммутирует с любым оператором с индексом “$y$”.

Полудискретная бикомпактная схема четвертого порядка аппроксимации для уравнения (12) строится аналогично одномерному случаю и тому, как это делается для уравнений гиперболического типа. Уравнение (12) и три его дифференциальных следствия:

(соответственно, производные по $x$, по $y$, смешанная производная по $x,y$) при фиксированном $t = {\text{const}}$ осредняются по прямоугольной ячейке Дальше действовать можно двумя способами. Первый способ повторяет рассуждения в одномерном случае: в ячейке ${{K}_{C}}$ строится двумерная версия полинома (5): где ${{y}_{{k + 1/2}}} = ({{y}_{k}} + {{y}_{{k + 1}}}){\text{/}}2$. Коэффициенты ${{a}_{m}} = {{a}_{m}}(t)$, $m = \overline {0,15} $, полинома (13) определяются значениями сеточных функций ${{u}_{{j,k}}}$, ${{{v}}_{{x,j,k}}}$, ${{{v}}_{{y,j,k}}}$, ${{{v}}_{{xy,j,k}}}$ в вершинах (углах) ячейки ${{K}_{C}}$. Отметим, что выражения для коэффициентов полинома легче получить не из решения системы линейных уравнений, а по аналогии с одномерным случаем, используя введенные выше сеточные операторы. Затем полином (13) подставляется вместо функции $u(x,y,t)$ в появляющиеся интегралы. Второй способ состоит в том, чтобы свести кратные интегралы к повторным и затем аппроксимировать одномерные интегралы по формулам (10). Оба способа приводят к одной и той же искомой схеме: где индексы $j = \overline {0,{{N}_{x}} - 1} $, $k = \overline {0,{{N}_{y}} - 1} $, а члены $\Phi _{C}^{1}(t), \ldots ,\Phi _{C}^{4}(t)$ выражаются как

Добавим, что интерполяция вида (13) называется бикубической (в трехмерном случае – трикубической). Каждая из совокупностей полиномов (5), (13) образует во всей соответствующей ей расчетной области гладкую кусочно-полиномиальную функцию, или сплайн гладкости 1. В связи с этим необходимо упомянуть работу [19], где уравнения Навье–Стокса для несжимаемой жидкости решались численно при помощи конечно-элементной схемы Галеркина, основанной на глобально гладкой аппроксимации точного решения трикубическим сплайном. Подход настоящей работы принципиально отличается от подхода [19] тем, что носители полиномов (5), (13) всегда занимают лишь одну ячейку сетки, а не несколько, как в [19].

Обратимся к вопросу о постановке начальных и граничных условий.

Корректная смешанная (начально-краевая) задача для одномерного уравнения (4) обязана включать в себя начальное условие

и два граничных условия, которые мы рассмотрим в общем виде: где ${{\alpha }_{1}},\;{{\alpha }_{2}},\;{{\beta }_{1}},\;{{\beta }_{2}}$ – заданные постоянные параметры. Если функции ${{u}_{*}}(x)$, ${{\mu }_{{1,2}}}(t)$ гладкие, то условия (15), (16) дискретизируются в виде (штрих сверху обозначает производную). Благодаря дифференцированию по времени в граничных условиях (18) индекс системы уравнений (11), (18) равен нулю (см. [20]), так как в ней отсутствуют алгебраические связи. В свою очередь, нулевой индекс позволяет избежать потенциального понижения порядка точности по времени при интегрировании системы (11), (18) высокоточными методами Рунге–Кутты.

Разумно проверить, является ли вообще одномерная схема (11) , дополненная условиями (17) и (18), корректной в смысле баланса числа уравнений и неизвестных. Сетка $\Omega $ состоит из ${{N}_{x}} + 1$ узлов. На ней определены две искомые сеточные функции ${{u}_{j}}(t)$ и ${{v}_{{x,j}}}(t)$, их начальные значения известны во всех узлах $\Omega $ из условия (17). Значит, число неизвестных (функций времени) равно $2{{N}_{x}} + 2$. На каждую ячейку ${{K}_{{j + 1/2}}}$, $j = \overline {0,{{N}_{x}} - 1} $, приходится по два уравнения схемы (11) , следовательно, общее число уравнений схемы составляет $2{{N}_{x}}$. Кроме того, мы имеем два граничных условия (18), что дает еще два уравнения. Всего уравнений получается $2{{N}_{x}} + 2$, столько же, сколько и неизвестных, т.е. схема (11) , дополненная условиями (17) и (18), действительно корректна.

Перейдем к двумерному уравнению (12). Корректная смешанная задача для него обязана включать в себя начальное условие

и граничные условия на всей $\partial D$, которые мы, как и в предыдущем одномерном случае, рассмотрим в общем виде: Здесь ${{\alpha }_{1}}, \ldots ,{{\alpha }_{4}},$ ${{\beta }_{1}}, \ldots ,{{\beta }_{4}}$ – заданные постоянные параметры. Предполагая, что функция ${{u}_{*}}(x,y)$ по крайней мере дважды непрерывно дифференцируемая, мы можем ясным образом дискретизировать начальное условие (19): При дискретизации граничных условий (20) мы не только проецируем их на сетку и дифференцируем по времени, но и дополнительно привлекаем их дифференциальные следствия – производные вдоль границы: Конечно, задавать дискретные граничные условия в виде (22) можно только для дважды непрерывно дифференцируемых функций ${{\mu }_{{1,2}}}(y,t)$, ${{\mu }_{{3,4}}}(x,t)$.

Проверим схему (14) , дополненную условиями (21) и (22) на предмет совпадения числа уравнений и неизвестных. Сетка $\Omega $ состоит из $({{N}_{x}} + 1)({{N}_{y}} + 1)$ узлов, на ней определены четыре искомые сеточные функции ${{u}_{{j,k}}}(t)$, ${{v}_{{x,j,k}}}(t)$, ${{v}_{{y,j,k}}}(t)$, ${{{v}}_{{xy,j,k}}}(t)$, начальные значения для которых известны из условия (21) во всех узлах $\Omega $. Таким образом, у нас есть $4({{N}_{x}} + 1)({{N}_{y}} + 1)$ неизвестных (функций времени). Теперь посчитаем число уравнений. Число ячеек сетки равно ${{N}_{x}}{{N}_{y}}$, на каждую из них приходится четыре уравнения схемы (14) , т.е. на всю сетку получается $4{{N}_{x}}{{N}_{y}}$ уравнений схемы. Граничные условия (22) дают еще

соотношений. Суммируя число уравнений схемы и число дискретных соотношений на границе, получаем т.е. число уравнений, равное числу неизвестных, что говорит в этом смысле о корректности схемы (14) , дополненной условиями (21) и (22).

Что касается трехмерного случая, то он рассматривается совершенно аналогично двумерному, поэтому мы не будем останавливаться на нем подробно. Приведем лишь само аппроксимируемое уравнение:

а также полудискретную бикомпактную схему для него (без вывода): где мультииндекс $C = (j + 1{\text{/}}2,k + 1{\text{/}}2,l + 1{\text{/}}2)$, индексы $j = \overline {0,{{N}_{x}} - 1} $, $k = \overline {0,{{N}_{y}} - 1} $, $l = \overline {0,{{N}_{z}} - 1} $, а члены $\Phi _{C}^{1}(t), \ldots ,\Phi _{C}^{8}(t)$ выражаются как

2. ПОЛНОСТЬЮ ДИСКРЕТНЫЕ БИКОМПАКТНЫЕ СХЕМЫ И МЕТОДЫ ИХ РЕАЛИЗАЦИИ

Полностью дискретные бикомпактные схемы получаются в результате применения к полудискретным схемам (11), (14), (23) тех или иных численных методов интегрирования по времени, например, методов Рунге–Кутты. В этой работе мы используем два таких метода: неявный метод Эйлера и $L$-устойчивый жестко-точный диагонально-неявный метод Рунге–Кутты (DIRK) второго порядка со следующей таблицей Бутчера [21, теорема 5]:

Бикомпактные схемы с аппроксимацией по времени неявным методом Эйлера мы будем сокращенно называть DIRK1B4 (независимо от числа пространственных измерений), а бикомпактные схемы с аппроксимацией по времени методом (24) – как SDIRK2B4. Погрешность аппроксимации у DIRK1B4 равна $O(\tau ,{{h}^{4}})$ ($\tau ,h \to 0$), у SDIRK2B4 – $O({{\tau }^{2}},{{h}^{4}})$, где $\tau > 0$ – шаг по времени, $h$ – максимальный пространственный шаг.

Рассмотрим произвольную полностью дискретную бикомпактную схему, интегрирование по времени в которой осуществляется многостадийным DIRK-методом, например, SDIRK2B4. Ясно, что реализация такой схемы на одном временном шаге сводится к композиции схем DIRK1B4 с правильно подобранными шагами по времени и промежуточными начальными и граничными условиями, подобно тому, как реализация многостадийных DIRK-методов сводится к композиции неявных методов Эйлера. Это обстоятельство делает достаточным обсуждение только реализации схем(ы) DIRK1B4.

Прежде всего опишем экономичный метод реализации одномерной схемы DIRK1B4. Он понадобится далее для описания реализации уже многомерных схем DIRK1B4.

Одномерная схема DIRK1B4 имеет вид

где $u_{j}^{{n + 1}}$ – неизвестная сеточная функция на временном слое ${{t}^{{n + 1}}}$ ($n = 0,1,2, \ldots $). Уравнения (25) дополняются граничными условиями, являющимися результатом применения неявного метода Эйлера к условиям (18):

Перепишем уравнения схемы (25) и граничные условия (26) в стандартном виде двухточечного векторного уравнения и граничных условий для него, см. [22]:

где Параметры $\kappa _{x}^{d}$ и $\kappa _{x}^{c}$ суть диффузионное и конвективное числа Куранта соответственно:

Отметим, что каждая из матриц $P_{{j + 1}}^{{ - 1}}{{Q}_{j}}$, $Q_{j}^{{ - 1}}{{P}_{{j + 1}}}$ при любых $\kappa _{x}^{d} > 0$, $\kappa _{x}^{c}$ имеет по два вещественных собственных значения, одно из которых больше 1, а другое – меньше 1, причем произведение этих значений в точности равно 1. Следовательно, стандартный метод двухточечной матричной прогонки [22] непригоден для решения уравнения (27): независимо от направления прогонки многократное перемножение матриц $P_{{j + 1}}^{{ - 1}}{{Q}_{j}}$ или $Q_{j}^{{ - 1}}{{P}_{{j + 1}}}$ при больших ${{N}_{x}}$ будет порождать очень большие числа, которые потенциально приведут к переполнению при счете на компьютере. Необходимо понимать, однако, что этот вывод говорит лишь о неустойчивости стандартного метода двухточечной матричной прогонки [22] при решении уравнения (27), сама же схема (25) абсолютно устойчива.

Предлагается решать уравнение (27) иначе. Во-первых, сделаем замену переменной:

введя также параметры отрицательные при любых $\kappa _{x}^{d} > 0$, $\kappa _{x}^{c}$. Во-вторых, умножим обе части (27) на матрицу В результате этих действий уравнение (27) запишется в виде Благодаря специальному выбору матриц $L$ и ${{H}_{{j + 1}}}$ матрицы ${{H}_{{j + 1}}}{{P}_{{j + 1}}}L$ и ${{H}_{{j + 1}}}{{Q}_{j}}L$ являются диагональными, и векторное уравнение (31) в действительности представляет из себя совокупность двух независимых скалярных уравнений: где $W_{j}^{{(1)}}$ и $W_{j}^{{(2)}}$ – компоненты вектора ${{W}_{j}}$, $F_{{j + 1}}^{{(1)}}$ и $F_{{j + 1}}^{{(2)}}$ – компоненты вектора ${{F}_{{j + 1}}}$. Первое уравнение в (32) решается слева-направо от граничного значения $W_{0}^{{(1)}}$, а второе – справа-налево от граничного значения $W_{{{{N}_{x}}}}^{{(2)}}$. Вводя функцию и величины мы приходим к прогоночным формулам для $W_{{j + 1}}^{{(1)}}$ и $W_{j}^{{(2)}}$: Дополнительный нижний индекс $j + 1{\text{/}}2$ у параметров ${{z}_{1}}$ и ${{z}_{2}}$ подчеркивает их зависимость от пространственного шага ${{h}_{x}} = {{h}_{{x,j + 1/2}}}$, который, как указывалось в начале разд. 1, может быть переменным. Добавим, что параметры ${{\lambda }_{1}}$ и ${{\lambda }_{2}}$, а также матрица $L$ не зависят от ${{h}_{x}}$, поэтому снабжать их индексом не нужно. Граничные же значения $W_{0}^{{(1)}}$, $W_{{{{N}_{x}}}}^{{(2)}}$ вычисляются из системы линейных уравнений, вытекающей из формул (28) и (35):

Таким образом, мы построили метод решения уравнения (27) и, тем самым, метод реализации одномерной схемы DIRK1B4 (25): он состоит из формул (29), (30), (33)–(36). Этот метод является безусловно устойчивым в силу неравенства $R(z) < 1$, выполненного при всех $z < 0$, и отрицательности параметров ${{z}_{{1,2}}}$ при любых $\kappa _{x}^{d} > 0$, $\kappa _{x}^{c}$.

Из формулы (35) очевидно, что предлагаемый метод реализации одномерной бикомпактной схемы DIRK1B4 принадлежит семейству методов прогонки. Это говорит об экономичности этой схемы и других бикомпактных схем на ее основе для уравнения (4). Более того, несмотря на векторный характер искомого численного решения (напомним, что неизвестных сеточных функций две, это $u_{j}^{{n + 1}}$ и $v_{{x,j}}^{{n + 1}}$), оно вычисляется посредством независимых скалярных прогонок, а не матричной прогонки, что обеспечивает еще большую эффективность счета.

Полезно сравнить одномерную бикомпактную схему DIRK1B4 для уравнения (4) с ее ближайшим аналогом – разрывной схемой Галеркина (DG-схемой) того же четвертого порядка аппроксимации по $x$, такой же дискретизацией по $t$ неявным методом Эйлера, для того же уравнения (4). Чтобы достичь четвертого порядка, такой DG-схеме понадобится аппроксимировать точное решение полиномом третьей степени в каждой ячейке [11]. В силу разрывности этой аппроксимации общее число неизвестных на каждом слое составит $4{{N}_{x}}$ у DG-схемы против почти двукратно меньшего числа $2({{N}_{x}} + 1)$ у DIRK1B4. Один шаг по времени в этой DG-схеме сводится к решению системы линейных алгебраических уравнений с блочно-трехдиагональной матрицей. Наилучшим прямым методом решения этой системы уравнений является трехточечная матричная прогонка, которая в случае данной DG-схемы требует обращения и перемножения матриц размера $4 \times 4$. Один шаг по времени в схеме DIRK1B4 сводится к двум независимым двухточечным скалярным прогонкам, что немногим сложнее двух шагов по времени по какой-либо простейшей явной схеме. Отсюда следует вывод, что схема DIRK1B4 (а) оперирует с меньшим числом неизвестных (т.е. требует меньших объемов машинной памяти) и (б) существенно быстрее совершает каждый шаг по времени.

Обратимся к двумерной схеме DIRK1B4:

Для реализации схемы (37) предлагается использовать итерационный метод, основанный на приближенной факторизации [23], [24] ее уравнений. Ранее этот метод успешно применялся для реализации многомерных бикомпактных схем для линейных и квазилинейных уравнений гиперболического типа [25], [26].

Пусть $\mathop u\limits^s $ – приближение номер $s$ к искомому решению ${{u}^{{n + 1}}}$ ($s = 0,1,2, \ldots $), причем начальное приближение $\mathop u\limits^0 = {{u}^{n}}$. Определим сеточные операторы

Заметим, что операторы $B_{1}^{x}(\tau )$ и $B_{2}^{x}(\tau )$ совпадают с операторами в скобках в левых частях одномерной схемы DIRK1B4 (25). Тогда уравнения метода итерируемой приближенной факторизации для схемы (37) , записанные относительно приращений $\mathop {\delta u}\limits^{s + 1} = \mathop u\limits^{s + 1} \; - \mathop u\limits^s $, принимают вид где Нетрудно видеть, что нахождение $\mathop {\delta u}\limits^{s + 1} $ из уравнений (39) сводится к последовательному решению четырех серий одномерных сеточных уравнений типа (25): где $U_{{j + 1/2,k}}^{{(1)}}$, $U_{{j + 1/2,k}}^{{(2)}}$ – вспомогательные сеточные функции, определенные при $j = \overline {0,{{N}_{x}} - 1} $, $k = \overline {0,{{N}_{y}}} $. Уравнения (40) дополняются граничными условиями а уравнения (41) – граничными условиями Уравнения (40), (41) решаются методом, построенным в первой половине этого раздела для реализации одномерной схемы DIRK1B4. Итерационный процесс останавливается, как только выполнено неравенство где $\varepsilon > 0$ – наперед заданная относительная точность вычислений.

Основной объем вычислений на одной итерации метода (39) состоит из одномерных скалярных двухточечных прогонок, т.е. из простейших явных действий. Следовательно, вопрос об экономичности метода (39) сводится к тому, сколько итераций придется выполнить, прежде чем наступит сходимость с точностью $\varepsilon $. Очевидно, число итераций до сходимости во многом определяется постановкой дифференциальной задачи, поэтому разумно исследовать вычислительную эффективность метода (39) практически на конкретных примерах, что и делается в следующем разделе.

Для реализации трехмерной схемы DIRK1B4 также предлагается использовать метод итерируемой приближенной факторизации. Уравнения этого метода в трехмерном случае могут быть без труда выписаны, исходя из полудискретной схемы (23) по аналогии с формулами (37) и (39). Мы не будем этого делать, так как эти уравнения громоздки и не вносят чего-то принципиально нового по сравнению с вышесказанным, кроме того, примеры расчетов в разд. 3 ограничиваются одномерным и двумерным случаями.

3. СЕТОЧНАЯ СХОДИМОСТЬ

Испытаем бикомпактные схемы DIRK1B4 и SDIRK2B4 на смешанных задачах для одномерного уравнения (4) и двумерного уравнения (12). Для каждого из уравнений рассмотрим по две задачи: одной стационарной и одной нестационарной, что позволит продемонстрировать точность бикомпактных схем отдельно по пространству и по времени.

Задача 1. Исследуем сходимость схемы DIRK1B4 на стационарном (предельном) решении уравнения (4):

где ${\mathbf{R}} = {{c}_{x}}{\text{/}}\nu $. На этот параметр можно смотреть как на некий аналог числа Рейнольдса. Решение (42) получается при следующих входных данных: Начальное условие можно взять любым; пусть это будет простейший линейный профиль, соединяющий граничные условия: Заметим, что при ${\mathbf{R}} \gg 1$ решение (42) содержит пограничный слой у правой границы расчетной области, т.е. у точки $x = 1$.

Для полноты проверки схемы DIRK1B4 и метода ее реализации, построенного в разд. 2, имеет смысл провести вычисления в широком диапазоне значений параметра ${\mathbf{R}}$. Полагая ${{c}_{x}} = 1$, возьмем $\nu = {{10}^{{ - 4}}}$ (${\mathbf{R}} = {{10}^{4}}$ есть погранслой очень малой толщины), $\nu = 1$ (${\mathbf{R}} = 1$), $\nu = {{10}^{4}}$ (${\mathbf{R}} = {{10}^{{ - 4}}}$). Для расчетов с ${\mathbf{R}} = 1,\;{{10}^{{ - 4}}}$ выберем равномерные сетки с узлами ${{x}_{j}} = j{{h}_{x}}$, ${{h}_{x}} = 1{\text{/}}{{N}_{x}}$, а для расчетов с ${\mathbf{R}} = {{10}^{4}}$ – неравномерные сетки с узлами

Число ячеек сетки ${{N}_{x}}$ возьмем равным $25,50, \ldots ,800$. В случае ${\mathbf{R}} = {{10}^{4}}$ при ${{N}_{x}} = 25$ на пограничный слой приходятся 7 ячеек. Сгущение сеток проводится при постоянном конвективном числе Куранта $\kappa _{x}^{c} = 10$ (для неравномерной сетки (43) оно равно 10 при максимальном ${{h}_{x}}$), что означает $\tau = O(h)$, $\kappa _{x}^{d} \to \infty $ при $h \to 0$. Счет по времени прекращается, если выполнено неравенство

На фиг. 1 в двойном логарифмическом масштабе изображены графики абсолютных погрешностей $u$ в норме ${{L}_{\infty }}$ как функций ${{N}_{x}}$. Из приведенных данных ясно, что схема DIRK1B4 действительно имеет четвертый порядок точности по пространству, и она не испытывает каких-либо проблем из-за очень малых или очень больших значений ${\mathbf{R}}$ или значительной неравномерности сетки (43).

Задача 2. Изучим сходимость схемы SDIRK2B4 на нестационарном решении уравнения (4) с $f(x,t) \equiv 0$:

Параметры возьмем следующими: ${{x}_{{max}}} = 2$, ${{c}_{x}} = 10$, $\nu = 0.01$ (${\mathbf{R}} = {{10}^{3}}$). Начальное и граничные условия для уравнения (4) получаются путем взятия решения (44) при $t = 0$, $x = 0$, $x = 2$ (такая смешанная задача моделирует задачу Коши на всей оси $Ox$).

Расчеты проведем до момента времени $t = 0.1$ на последовательности равномерных сеток с ${{N}_{x}} = 25,50, \ldots ,800$. Сгущение сетки при этом осуществим двумя способами: при $\kappa _{x}^{c} = 0.5 = {\text{const}}$ ($\tau = O(h)$, $\kappa _{x}^{d} \to \infty $) и $\kappa _{x}^{c} = 12.5{\text{/}}{{N}_{x}} = O(h)$ ($\tau = O({{h}^{2}})$, $\kappa _{x}^{d} = {\text{const}}$).

На фиг. 2 показаны графики абсолютных погрешностей $u$ в норме ${{L}_{\infty }}$ как функций ${{N}_{x}}$. Погрешности посчитаны в конечный момент времени $t = 0.1$. Хорошо видно, что при постоянном конвективном числе Куранта схема SDIRK2B4 сходится со вторым порядком, а при постоянном диффузионном числе Куранта – с четвертым порядком, как того и следовало ожидать.

Задача 3. Вернемся к схеме DIRK1B4 и проверим ее сходимость на стационарном решении теперь уже двумерного уравнения (12) с $f(x,y,t) \equiv 0$:

Параметры расчетной области ${{x}_{{max}}} = {{y}_{{max}}} = 1$. Для постановки граничных условий берутся значения функции (45) на границе $\partial D$ расчетной области $D$. Начальное условие можно выбрать любым, пусть это будет тождественный нуль: $u{{{\text{|}}}_{{t = 0}}} \equiv 0$. Для двумерных задач определим параметр ${\mathbf{R}}$ в виде Решение (45) не содержит каких-либо особенностей (например, погранслоев) ни в пределе при ${\mathbf{R}} \to 0$, ни в пределе при ${\mathbf{R}} \to \infty $.

По аналогии с задачей 1 зафиксируем коэффициенты ${{c}_{x}} = 1$, ${{c}_{y}} = 2$ и проведем расчеты при разных значениях коэффициента $\nu $: $\nu = \sqrt 5 \times {{10}^{{ - 4}}}$ (${\mathbf{R}} = {{10}^{4}}$), $\nu = \sqrt 5 $ (${\mathbf{R}} = 1$), $\nu = \sqrt 5 \times {{10}^{4}}$ (${\mathbf{R}} = {{10}^{{ - 4}}}$). Несмотря на формальное отсутствие у решения особенностей при малых и больших ${\mathbf{R}}$, варьировать $\nu $ в таком широком диапазоне не бессмысленно: тем самым мы испытаем метод итерируемой приближенной факторизации при сильном различии коэффициентов в уравнениях схемы. Для всех ${\mathbf{R}}$ вычисления выполним на равномерных сетках; при ${\mathbf{R}} = {{10}^{4}}$ используем сетки с ${{N}_{x}} = {{N}_{y}} = 50$, 100, 200, 400, а при ${\mathbf{R}} = 1,\;{{10}^{{ - 4}}}$ – сетки с ${{N}_{x}} = {{N}_{y}} = 5$, 10, 20, 40 (и дополнительно 80 для ${\mathbf{R}} = 1$). Сгущение сеток делается при постоянном конвективном числе Куранта, независимо от значения ${\mathbf{R}}$; $\kappa _{y}^{c} = 2\kappa _{x}^{c} = 1,\;0.1,\;{{10}^{{ - 5}}}$ при ${\mathbf{R}} = {{10}^{4}},\;1,\;{{10}^{{ - 4}}}$ соответственно. Параметр $\varepsilon = {{10}^{{ - 10}}}$, счет по времени останавливается при выполнении неравенства

На фиг. 3 изображены графики абсолютных погрешностей $u$ в норме ${{L}_{\infty }}$ как функций ${{N}_{x}}$. Очевидно, во всех случаях схема DIRK1B4 сходится с четвертым порядком точности (для ${\mathbf{R}} = {{10}^{{ - 4}}}$ он даже несколько больше четырех). Из фиг. 1 становится понятным выбор столь грубых сеток для ${\mathbf{R}} = 1,\;{{10}^{{ - 4}}}$ – уже при ${{N}_{x}} = {{N}_{y}} = 40{\kern 1pt} - {\kern 1pt} 80$ точность схемы приближается к точности итераций по факторизации.

Важно также привести и проанализировать данные о сходимости итераций по факторизации. Независимо от значения ${\mathbf{R}}$, число итераций до сходимости ${{s}_{\varepsilon }}$ меняется при увеличении $n$ таким образом: при $n = 1$ (первый шаг по времени) число ${{s}_{\varepsilon }}$ очень велико, затем оно резко падает до нескольких десятков итераций и начинает монотонно убывать, и после некоторого $n$ становится меньше десяти и быстро выходит на константу в три итерации. Для ${\mathbf{R}} = {{10}^{4}}$ самое первое большое значение ${{s}_{\varepsilon }}$ составляет примерно 750, а меньше десяти число ${{s}_{\varepsilon }}$ становится примерно при $n = 50$. Для ${\mathbf{R}} = 1,\;{{10}^{{ - 4}}}$ первое ${{s}_{\varepsilon }}$ равно примерно 2000, меньше десяти ${{s}_{\varepsilon }}$ становится примерно при $n = 10$. Такой характер зависимости ${{s}_{\varepsilon }}(n)$ объясняется тем, что (а) начальное приближение на первом шаге по времени совпадает с достаточно грубо выбранным начальным условием, которое даже “не попадает” в граничные условия, и (б) с ростом $n$ численное решение устанавливается из-за стационарности постановки задачи.

Задача 4 является двумерной версией задачи 2. Рассматривается сходимость схемы SDIRK2B4 на нестационарном решении уравнения (12) с $f(x,y,t) \equiv 0$:

Параметры ${{x}_{{max}}} = {{y}_{{max}}} = 2$, ${{c}_{x}} = {{c}_{y}} = 10$, $\nu = 0.01$ (${\mathbf{R}} = \sqrt 2 \times {{10}^{3}}$). Начальное и граничные условия для уравнения (12), подобно задаче 2, ставятся путем взятия точного решения (46) при $t = 0$ и на $\partial D$ (получающаяся смешанная задача моделирует задачу Коши на всей плоскости $Oxy$).

Расчеты выполним до момента времени $t = 0.1$ на последовательности равномерных сеток с ${{N}_{x}} = {{N}_{y}} = 50$, 100, 200, 400. Сгущение сетки при этом осуществим двумя способами: при $\kappa _{{x,y}}^{c} = 0.5 = {\text{const}}$ ($\tau = O(h)$, $\kappa _{{x,y}}^{d} \to \infty $) и $\kappa _{{x,y}}^{c} = 12.5{\text{/}}{{N}_{x}} = O(h)$ ($\tau = O({{h}^{2}})$, $\kappa _{{x,y}}^{d} = {\text{const}}$). Как и в задаче 3, относительная точность итераций по факторизации $\varepsilon = {{10}^{{ - 10}}}$.

На фиг. 4 приведены графики абсолютных погрешностей $u$ в норме ${{L}_{\infty }}$ как функций ${{N}_{x}}$. Нетрудно видеть, что апостериорные порядки сходимости схемы SDIRK2B4 совпадают с теоретическими, вторым при $\kappa _{{x,y}}^{c} = {\text{const}}$ и четвертым при $\kappa _{{x,y}}^{c} = O(h)$. Поскольку решение гладкое и мало меняется за один шаг по времени, метод итерируемой факторизации сходится очень быстро: На всех слоях число ${{s}_{\varepsilon }} = 6$, 5, 4, 3 (среднее по стадиям метода Рунге–Кутты) при ${{N}_{x}} = 50$, 100, 200, 400 соответственно.

ЗАКЛЮЧЕНИЕ

В работе предложены новые бикомпактные схемы для линейного многомерного уравнения конвекции-диффузии. Эти схемы получены при помощи метода прямых, интегроинтерполяционного метода, а также би- и трикубической интерполяции Эрмита искомой функции в ячейке сетки. Для интегрирования по времени в этих схемах рекомендуется использовать диагонально-неявные методы Рунге–Кутты. Построенные бикомпактные схемы абсолютно устойчивы, консервативны, имеют четвертый порядок аппроксимации по пространству. Для реализации этих схем разработан итерационный метод, основанный на приближенной факторизации их многомерных уравнений. Каждая итерация этого метода не является вычислительно трудоемкой, так как сводится к совокупности одномерных скалярных двухточечных прогонок.

На одномерных и двумерных, стационарных и нестационарных смешанных задачах проверена сеточная сходимость бикомпактных схем с аппроксимацией по времени неявным методом Эйлера и $L$-устойчивым жестко-точным диагонально-неявным методом Рунге–Кутты второго порядка. Расчеты проведены в широком диапазоне отношений конвективных скоростей и коэффициента диффузии, от ${{10}^{4}}$ до ${{10}^{{ - 4}}}$. В одной из задач для разрешения очень тонкого пограничного слоя использована существенно неравномерная сетка. На стационарных задачах продемонстрирован четвертый порядок точности по пространству у выбранных бикомпактных схем. На нестационарных задачах бикомпактная схема второго порядка аппроксимации по времени показала второй порядок точности по времени при постоянном конвективном числе Куранта и четвертый при постоянном диффузионном числе Куранта. Число итераций до сходимости по факторизации в нестационарной двумерной задаче на не слишком грубых сетках оказалось равным 3–5, что вкупе с вычислительной простотой каждой итерации говорит о высокой экономичности предлагаемых бикомпактных схем.

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