Журнал вычислительной математики и математической физики, 2019, T. 59, № 3, стр. 481-493

ВВЕДЕНИЕ

Квазигазо/квазигидродинамические (КГД/КГДД) системы уравнений можно рассматривать как регуляризованные системы уравнений газовой динамики. С ними связан класс явных двухслойных по времени и симметричных по пространству разностных схем, достаточно активно используемых при численном решении многих задач [1]–[3]. Баротропные КГД/КГДД системы были введены и изучены в [4]–[7]. Соответствующие разностные схемы применялись при расчете задач мелкой воды, двухкомпонентных течений с поверхностными эффектами, астрофизических задач и др. [8]–[11].

Вопросы теории указанных разностных схем, включая практически важный анализ условий их устойчивости в зависимости от параметров схем, пока разработаны недостаточно. Для 1D разностных схем, линеаризованных на постоянном решении с нулевой скоростью ${{u}_{ * }} = 0$, достаточные условия ${{L}^{2}}$-диссипативности масштабированных решений были получены в [3], [12] энергетическим методом. Более точные критерии (необходимые и достаточные условия) в баротропном случае при ${{u}_{ * }} = 0$ выведены недавно в [13] спектральным методом.

В настоящей работе в баротропной постановке впервые анализируется случай ${{u}_{ * }} \ne 0$ и выводятся как необходимые, так и достаточные условия ${{L}^{2}}$-диссипативности решений задачи Коши в зависимости от числа Маха M. Эти условия различаются между собой не более чем в 2 раза и содержат существенную для практики информацию. В том числе выяснено, что при использовании КГД системы можно обеспечить независимость условия на число Куранта от ${\text{M}}$. Вывод условий базируется на новых результатах для абстрактных разностных схем с регуляризацией и спектральном подходе (см. [14]). Эти результаты существенно развивают известные для линеаризованной схемы Лакса–Вендроффа в [15]. Они применимы также для анализа схем изучаемого типа для полной 1D системы уравнений газовой динамики [16].

Приводятся также результаты численных экспериментов с использованием КГД системы в связи с анализом применимости найденных условий в нелинейной постановке при решении задачи Римана в изэнтропическом случае при ${\text{M}}$ = 2,4,6,8. Расчеты выполнены для трех схем – схемы стандартного типа и двух схем с энергетически диссипативными дискретизациями по пространству из [17]. Подчеркнем, что линеаризация всех трех схем на постоянном решении одинаковая. Однако в нелинейной постановке свойства устойчивости последних двух схем оказываются существенно лучше, чем у схемы стандартного типа, тем более с ростом ${\text{M}}$, и эти свойства в определенной степени соответствуют найденным условиям.

1. УСЛОВИЯ ${{L}^{2}}$-ДИССИПАТИВНОСТИ АБСТРАКТНОЙ ЯВНОЙ ДВУХСЛОЙНОЙ СХЕМЫ С РЕГУЛЯРИЗАЦИЕЙ

Пусть ${{\omega }_{h}}$ – равномерная сетка на $\mathbb{R}$ с узлами ${{x}_{k}} = kh$, $k \in \mathbb{Z}$ и шагом $h = X{\text{/}}N$, а $\omega _{h}^{ * }$ – вспомогательная сетка с узлами ${{x}_{{k - 1/2}}} = (k - 0.5)h$, $k \in \mathbb{Z}$. Пусть $\mathop {\bar {\omega }}\nolimits^{\Delta t} $ – равномерная сетка по $t$ на $[0,\infty )$ с узлами ${{t}_{m}} = m\Delta t$, $m \geqslant 0$, и шагом $\Delta t > 0$. Введем сеточные операторы

Пусть $H$ – гильбертово пространство вектор-функций со значениями из ${{\mathbb{C}}^{n}}$, заданных и суммируемых в квадрате на ${{\omega }_{h}}$, со скалярным произведением ${{({\mathbf{v}},{\mathbf{y}})}_{H}} = h\sum\nolimits_{k = - \infty }^{ + \infty } {{{{({{{\mathbf{v}}}_{k}},{{{\mathbf{y}}}_{k}})}}_{{{{\mathbb{C}}^{n}}}}}} $.

Рассмотрим сначала абстрактную явную двухслойную по времени и симметричную трехточечную по пространству линейную разностную схему

где ${{{\mathbf{y}}}^{m}} \in H$ при $m \geqslant 0$, $A = A{\text{*}}$ и $B = B{\text{*}}$ – матрицы (вязких и конвективных слагаемых) порядка $n$, ${{c}_{0}} > 0$ – характерная скорость, ${{\tau }_{ * }} > 0$ – параметр регуляризации. При $A = {{B}^{2}}$, ${{\tau }_{ * }} = \Delta t{\text{/}}2$ устойчивость подобной линеаризованной схемы Лакса–Вендроффа изучена в [15].

Поставим вопрос об условиях справедливости оценки

Нетрудно убедиться, что она эквивалентна свойству ${{L}^{2}}$–диссипативности Перепишем схему (1) в рекуррентном виде: Пусть $\Delta t$ и ${{\tau }_{ * }}$ задаются формулами с параметрами $\tilde {\beta } > 0$ и $\alpha > 0$. Цель дальнейшего анализа – вывести условия на $\tilde {\beta }$ в зависимости от $\alpha $, обеспечивающие выполнение оценки (2).

Для этого в соответствии с [14] рассмотрим частные решения вида ${\mathbf{y}}_{k}^{m} = {{e}^{{{\mathbf{i}}k\xi }}}{{{\mathbf{v}}}^{m}}(\xi )$, $k \in \mathbb{Z}$, $m \geqslant 0$, где ${\mathbf{i}}$ – мнимая единица и $0 \leqslant \xi \leqslant 2\pi $ – параметр. Их подстановка в (3) с учетом (4) дает

Здесь матрица $G$ – символ оператора ${\mathbf{A}}$, матрица $I$ – единичная, $\sigma = si{{n}^{2}}\tfrac{\xi }{2} \in [0,1]$ и берутся оба знака $ + $ и $ - $, отвечающие соответственно $0 \leqslant \xi \leqslant \pi $ и $\pi < \xi \leqslant 2\pi $.

Следующая теорема вытекает из [16]. Для полноты дадим ее с доказательством.

Теорема 1. 1. Выполнение матричного неравенства

необходимо и достаточно для справедливости оценки (2). Здесь $[A,B] = AB - BA$ – коммутатор матриц $A$ и $B$, и матрица ${\mathbf{i}}[A,B]$ эрмитова.

2. Выполнение матричных неравенств

необходимо, а выполнение неравенствадостаточно для справедливости оценки (2).

Доказательство. 1. По аналогии с леммой 1 из [13] доказывается, что выполнение спектрального неравенства

необходимо и достаточно для справедливости оценки (2); здесь и ниже ${{\lambda }_{{max}}}(C)$ – максимальное из собственных значений матрицы $C$ (когда они вещественные).

Это неравенство эквивалентно матричному неравенству $G{\text{*}}G \leqslant I$, т.е. неравенству

Для $F = {{F}_{R}} + {\mathbf{i}}{{F}_{I}}$, где матрицы ${{F}_{R}}$ и ${{F}_{I}}$ – эрмитовы матрицы, оно принимает вид В силу последней формулы (5) после сокращения на $8\sigma \alpha $ при $\sigma \ne 0$ оно переходит в (6).

2. По непрерывности неравенство (6) выполняется и при $\sigma = 0$. При $\sigma = 1,0$ оно дает

Эти неравенства можно переписать в виде (7).

Далее, для любого ${\mathbf{z}} \in H$ справедливы соотношения

Поэтому при $0 \leqslant \sigma \leqslant 1$ верно матричное неравенство (при $\sigma = 0,1$ оно очевидно). При любом $\varepsilon > 0$ оно влечет достаточное условие При $\varepsilon = (1 - \sigma ){\text{/}}\sigma $ (для $0 < \sigma < 1$) получаем (8).

Замечание 1. 1. Для $A \geqslant 0$ первое неравенство (7) эквивалентно неравенству $\tilde {\beta } \leqslant 1{\text{/}}[2\alpha {{\lambda }_{{max}}}(A)].$ Второе неравенство (7) влечет свойство $A \geqslant 0$.

2. Если $A = {{B}^{2}} + D$, где $D$ – диагональная матрица с элементами ${{d}_{{ii}}} \geqslant 0,\;1 \leqslant i \leqslant N$, причем ${{d}_{{{{i}_{0}}{{i}_{0}}}}} = 0$, а ${{a}_{{{{i}_{0}}{{i}_{0}}}}} \ne 0$ для некоторого ${{i}_{0}}$, то второе неравенство (7) сводится просто к числовому неравенству $\tilde {\beta } \leqslant 2\alpha $.

3. Если ${{B}^{2}} \leqslant A$, то неравенство (8) выполнено при $\tilde {\beta } \leqslant 1{\text{/}}\left[ {2\alpha {{\lambda }_{{max}}}(A) + {{{(2\alpha )}}^{{ - 1}}}} \right].$

Замечание 2. Неравенства (7) с удвоенными левыми частями

влекут неравенство (8) и поэтому служат более простым достаточным условием.

2. РАЗНОСТНЫЕ СХЕМЫ С РЕГУЛЯРИЗАЦИЕЙ ДЛЯ 1D УРАВНЕНИЙ БАРОТРОПНОЙ ГАЗОВОЙ ДИНАМИКИ, ИХ ЛИНЕАРИЗАЦИЯ И L²–ДИССИПАТИВНОСТЬ

1D баротропная КГД система уравнений состоит из регуляризованных уравнений баланса массы и импульса

Здесь $\rho > 0$, $u$ – плотность и скорость, $p(\rho )$ – давление с $p \in {{C}^{1}}(0, + \infty )$ и $p{\kern 1pt} ' > 0$, $j$, $\Pi $ – регуляризованные поток массы и вязкое напряжение, $\rho w$, $\rho \hat {w}$ – регуляризующие импульсы, $\tau = \tau (\rho ,u) > 0$ – параметр регуляризации, $\mu {{\partial }_{x}}u$ – вязкое напряжение Навье–Стокса, где $\mu = \mu (\rho ,u) > 0$ пропорционален коэффициенту вязкости.

КГДД система проще и отличается тем, что в ней $\rho w = \rho \hat {w} = \tau [\rho u{{\partial }_{x}}u + p{\kern 1pt} '(\rho )]$, $\Pi = \mu {{\partial }_{x}}u + u\rho \hat {w}$.

Обе системы переходят в баротропные системы уравнений газовой динамики при $\tau = \mu = 0$ и Навье–Стокса вязкого сжимаемого газа при $\tau = 0$, $\mu > 0$.

Сначала рассмотрим явную двухслойную по времени и симметричную по пространству стандартного типа аппроксимацию КГД уравнений (11)–(13):

Основные искомые функции $\rho > 0$, $u$, а также $\tau $ определены на ${{\omega }_{h}}$, а $j,w,\hat {w},\Pi ,\mu $ – на $\omega _{h}^{*}$.

В [17] были построены две нестандартные пространственные аппроксимации КГД уравнений (11)–(13), диссипативные по энергии (их обобщения на многомерный случай даны в [18]). В сочетании с явным методом Эйлера по времени первая из них (схема А) имеет вид

Здесь ${\text{h}}(\rho ) = \int_{{{r}_{0}}}^\rho \tfrac{{p{\kern 1pt} '(r)}}{r}dr$ (с некоторым ${{r}_{0}} > 0$) – энтальпия, причем ${\text{h'}}(\rho ) = p{\kern 1pt} '(\rho ){\text{/}}\rho $. В изэнтропическом случае $p(\rho ) = {{p}_{1}}{{\rho }^{\gamma }}$ с ${{p}_{1}} > 0$, $\gamma > 1$ и можно взять ${{r}_{0}} = 0$; тогда ${\text{h}}(\rho ) = \gamma p(\rho ){\text{/}}[(\gamma - 1)\rho ]$ и ${\text{h'}}(\rho ) = \gamma p(\rho ){\text{/}}{{\rho }^{2}}$.

Обратим внимание на нестандартные зависящие от ${\text{h}}(\rho )$ дискретизации ${{\partial }_{x}}p(\rho )$ в (17), (19) и $\tau {{\partial }_{x}}(\rho u)$ в (18), (19) согласно (20).

Схема Б со второй из аппроксимаций имеет вид

Здесь, например, $\tau = \tau (s\rho ,su)$ на $\omega _{h}^{*}$.

В схеме Б наряду с обычными усреднениями $s\rho $, $su$ применяются нестандартные усреднения $\rho $ и $p{\kern 1pt} '(\rho )$, зависящие от ${\text{h}}$:

где ${{\rho }_{{ \pm ,k - 1/2}}} = {{\rho }_{{k - 1/2 \pm 1/2}}}$. Кроме того, $g(\alpha ;\beta )$ – разделенная разность функции $g \in {{C}^{1}}(0, + \infty )$:

В изэнтропическом случае введенные усреднения принимают вид

В частном случае $\gamma = 2$ (случай уравнений мелкой воды) они становятся стандартными: ${{s}_{p}}\rho = s\rho $, $\widetilde {p{\kern 1pt} '(\rho )} = 2{{p}_{1}}s\rho $. На практике во избежание потери точности при малых ${\text{|}}({{\rho }_{ + }}{\text{/}}{{\rho }_{ - }}) - 1{\text{|}} < {{\varepsilon }_{0}}$ вместо указанных формул используются квадратурные формулы для их интегральных представлений, например, формула парабол. В случае усреднений (21) эта формула принимает вид который нестандартен, но несущественно сложнее простейших усреднений.

Схему (14)–(16) линеаризуем на постоянном решении ${{\rho }_{ * }} \equiv {\text{const}} > 0$, ${{u}_{ * }} \equiv {\text{const}}$. Пусть ${{c}_{ * }} = \sqrt {p{\kern 1pt} '({{\rho }_{ * }})} $ – фоновая скорость звука, а ${{\tau }_{ * }} = \tau ({{\rho }_{ * }},{{u}_{ * }})$, ${{\mu }_{ * }} = \mu ({{\rho }_{ * }},{{u}_{ * }})$. Запишем решение схемы в виде $\rho = {{\rho }_{ * }}(1 + \tilde {\rho })$, $u = {{u}_{ * }} + {{c}_{ * }}\tilde {u}$, отбросим члены второго порядка малости по отношению к масштабированным (безразмерным) возмущениям $\tilde {\rho }$, $\tilde {u}$ и получим линеаризованную схему

Пусть коэффициент вязкости задается обычной КГД-формулой $\mu = {{\alpha }_{S}}\tau (\rho ,u)\rho p{\kern 1pt} '(\rho )$, где ${{\alpha }_{S}} \geqslant 0$ – параметр. Тогда ${{\mu }_{ * }} = {{\alpha }_{S}}{{\tau }_{ * }}{{\rho }_{ * }}c_{ * }^{2}$. Положив ${{c}_{0}} = {{c}_{ * }}$, можно схему (22), (23) переписать в векторном виде (1) с где ${\text{M}} = {{u}_{ * }}{\text{/}}{{c}_{ * }}$ (при этом ${\text{|M|}}$ – число Маха). Полученная схема имеет прямое отношение к линеаризованной баротропной КГД системе уравнений (см. [4], [6]).

Линеаризации на постоянном решении схем А и Б имеют такой же вид.

Отметим, что первая формула (4) на практике обычно берется в виде

где $\beta $ – число Куранта. Вторую формулу (4) перепишем в модифицированном виде: Подобная формула использовалась в [11].

Теорема 2. Рассмотрим схему (22), (23) , т.е. (1), (4), (24).

1. Для нее необходимое условие ${{L}^{2}}$-диссипативности (7) принимает вид

где

В случае ${{\alpha }_{S}} = 0$ правая часть этого условия максимальна, а само условие становится не только необходимым, но и достаточным, и упрощается до

2. Достаточное условие ${{L}^{2}}$-диссипативности (8) выполнено при

Доказательство. Матрицы $A$, $B$ в (24) вещественны и симметричны. При этом

Поэтому с помощью замечаний 1.1, 1.2 заключаем, что необходимое условие (7) принимает вид (25), а с помощью замечания 1.3 получаем, что достаточное условие (8) выполнено при условии (27).

Кроме того, при ${{\alpha }_{S}} = 0$ имеем ${{B}^{2}} = A$, $[A,B] = 0$ и, значит, неравенство (6) сводится к

Оно эквивалентно двум неравенствам: $\tilde {\beta } \leqslant 1{\text{/}}[2\alpha {{\lambda }_{{{\text{max}}}}}(A)]$ с ${{\lambda }_{{{\text{max}}}}}(A) = {{(\left| {\text{M}} \right| + 1)}^{2}}$ и $\tilde {\beta } \leqslant 2\alpha $, т.е. неравенству (26).

Обратим внимание на то, что при ${{\alpha }_{S}} > 0$ величина ${{\widehat \alpha }_{{{\text{opt}}}}}$ возрастает по ${\text{|M|}}$ (как нетрудно проверить) и стремится к $1{\text{/}}2$ при $\left| {\text{M}} \right| \to \infty $. Поэтому член $\widehat \alpha _{{{\text{opt}}}}^{2}{\text{/}}\widehat \alpha $ в условии (25) и правая часть условия (27) возрастают по $\left| {\text{M}} \right|$ и стремятся к $1{\text{/}}(4\widehat \alpha )$ и $1{\text{/}}[2\widehat \alpha + {{(2\widehat \alpha )}^{{ - 1}}}]$ соответственно при $\left| {\text{M}} \right| \to \infty $, и в пределе при $\left| {\text{M}} \right| \to \infty $ условие (25) снова переходит в (26).

Очевидно, что правые части условий (25) и (27) максимальны при $\widehat \alpha = {{\widehat \alpha }_{{{\text{opt}}}}}$ и равны соответственно $2{{\widehat \alpha }_{{{\text{opt}}}}}$ и ${{\widehat \alpha }_{{{\text{opt}}}}}$; при этом

Существенно, что ${{\widehat \alpha }_{{{\text{opt}}}}}$ лежит в сегменте $[1{\text{/}}(2\sqrt {{{\alpha }_{S}} + 1} ),1{\text{/}}2]$, не зависящем от ${\text{|M|}}$.

Если вместо указанной выше формулы ${{\tau }_{ * }} = \widehat \alpha h{\text{/}}({{c}_{ * }} + \left| {{{u}_{ * }}} \right|)$ использовать стандартную ${{\tau }_{ * }} = \alpha h{\text{/}}{{c}_{ * }}$, то в условиях (25) и (27) следует заменить $\widehat \alpha $ на $\alpha ({\text{|M|}} + 1)$. При этом при ${{\alpha }_{S}} > 0$ значение ${{\alpha }_{{{\text{opt}}}}} = {{\widehat \alpha }_{{{\text{opt}}}}}{\text{/}}(\left| {\text{M}} \right| + 1)$ станет убывать по ${\text{|M|}}$ и стремиться к $0$ при $\left| {\text{M}} \right| \to \infty $, что неудобно на практике.

В типичном случае ${{\alpha }_{S}} = 1$ правые части условий (25) и (27) (верхняя и нижняя поверхности) представлены на фиг. 1а. Результаты замены в них $\widehat \alpha $ на $\alpha (\left| {\text{M}} \right| + 1)$ представлены на фиг. 1б; наблюдаются негативные изменения в графиках.

При ${\text{M}} = 0$, т.е. ${{u}_{ * }} = 0$, условие (25) переходит в

полученное в [13] с помощью прямого анализа условия (9) и являющееся критерием (и необходимым, и достаточным условием).

Обратимся к КГДД системе. С учетом указанных выше упрощений ей отвечает упрощенная линеаризованная схема (22), (23) , а именно, схема (1) с теми же ${\mathbf{y}}$, $B$ и иной $A$:

Теорема 3. Рассмотрим схему (1), (4), (28) .

1. Для нее необходимое условие ${{L}^{2}}$-диссипативности (7) принимает вид

где $q({{\alpha }_{S}},{\text{M}}) = [{{({{{\text{M}}}^{2}} - 1)}^{2}} + {{\alpha }_{S}}({{{\text{M}}}^{2}} + 1)]{\text{/}}2$ и

2. Достаточное условие ${{L}^{2}}$-диссипативности (10) получается из (29) в результате замены $\beta $ на $2\beta $.

Доказательство. Матрицы $A$, $B$ в (28) вещественны и симметричны. Легко проверяются свойство $A \geqslant 0$ и формула (30). С помощью замечания 1.1 первое неравенство (7) принимает вид $\tilde {\beta } \leqslant 1{\text{/}}[2\alpha {{\lambda }_{{{\text{max}}}}}(A)]$.

Анализ второго неравенства (7) сложнее. Положив $\zeta = \tilde {\beta }{\text{/}}(2\alpha )$, перепишем его в виде

Это неравенство эквивалентно неотрицательности диагональных элементов и определителя матрицы $A - \zeta {{B}^{2}}$, т.е. неравенствам При ${\text{|M|}} = 1$ и ${{\alpha }_{S}} > 0$ последнее неравенство означает, что $\zeta \leqslant 1{\text{/}}2$, а это следует из (31).

Пусть ${\text{|M|}} \ne 1$. Заметим, что

Поэтому с учетом (31) неравенство (32) означает, что где ${{\zeta }_{ - }}$ и ${{\zeta }_{{\text{в }}}}$ – меньший из корней и вершина квадратного трехчлена ${{p}_{2}}(\zeta )$ в (32). Отсюда и вытекает неравенство (29).

Если вместо указанной выше формулы ${{\tau }_{ * }} = \widehat \alpha h{\text{/}}({{c}_{ * }} + {\text{|}}{{u}_{ * }}{\text{|}})$ использовать стандартную ${{\tau }_{ * }} = \alpha h{\text{/}}{{c}_{ * }}$, то опять в условии (29) следует заменить $\widehat \alpha $ на $\alpha ({\text{|M|}} + 1)$.

В типичном случае ${{\alpha }_{S}} = 1$ правая часть условия (29) (верхняя поверхность) и его половина – в соответствии с теоремой 3, п. 2 (нижняя поверхность) представлены на фиг. 2а. Результаты замены в них $\widehat \alpha $ на $\alpha (\left| {\text{M}} \right| + 1)$ представлены на фиг. 2б; возникающие различия опять существенны. Графики достаточно быстро стремятся к 0 с ростом ${\text{|M|}}$; это вполне соответствует тому, что КГДД система изначально не предназначена для расчета течений с большими числами Маха. Сами же полученные выше результаты представляют интерес и при умеренных числах Маха.

При ${\text{M}} = 0$, т.е. ${{u}_{ * }} = 0$, условие (29) переходит в условие

полученное в [13] с помощью прямого анализа условия (9) и являющееся критерием.

Отметим, что при ${{\alpha }_{S}} = 0$ схема (1), (4), (28) неустойчива (правая часть (29) равна 0) в отличие от схемы (1), (24) .

3. ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ И АНАЛИЗ УСТОЙЧИВОСТИ В НЕЛИНЕЙНОЙ ПОСТАНОВКЕ

Численные эксперименты выполнены для модельных уравнения состояния $p(\rho ) = {{\rho }^{\gamma }}$ с $\gamma = 5{\text{/}}3$ и задачи Римана с начальными данными

где U = $\sqrt \gamma {{{\text{M}}}_{0}}$, M₀ – число Маха при x < 0. Задача описывает набегание движущейся массы газа на покоящуюся вдвое менее плотную. В течении между ними возникает слой газа со значением плотности больше 1 и скоростью из (0, U); числа Маха лежат от 0 до M₀.

Задача решалась с применением КГД-системы для ${{M}_{0}} = 2,4,6,8$ на отрезке $[ - X{\text{/}}2,X{\text{/}}2]$ с $X = 0.6$ и шагом $h = 0.002$ при $0 \leqslant t \leqslant {{t}_{{{\text{ф и н }}}}}$ с ${{t}_{{{\text{ф и н }}}}} = 0.05$. Брались модифицированная формула $\tau = \widehat \alpha h{\text{/}}({\text{|}}u{\text{|}} + {{c}_{s}}(\rho ))$, где ${{c}_{s}}(\rho ) = \sqrt {p{\kern 1pt} '(\rho )} $ – скорость звука, и ${{\alpha }_{S}} = 1$.

Результаты расчетов отражены на фиг. 3–5 для схемы стандартного типа и схем А и Б соответственно. На них пунктирной и сплошной линией изображены графики правых частей необходимых условий (25) при ${\text{M}} = {{{\text{M}}}_{0}}$ и (26) соответственно. Различие между ними невелико, уменьшается с ростом ${{{\text{M}}}_{0}}$ и визуально исчезает при ${{{\text{M}}}_{0}} = 8$.

Тестировались значения $\widehat \alpha = 0.3,0.4,\; \ldots ,0.9$ и $\beta = 0.1kmin\{ 2\widehat \alpha ,1{\text{/}}(2\widehat \alpha )\} $, $k = 1,2,\; \ldots ,\;10$ (и $k = 11$ для схем А и Б) в долях от его максимального значения в условии (26).

Для анализа качества численных решений было решено использовать величину ${{\varepsilon }_{V}}$ – отклонение их относительной вариации от 1 (на момент времени $t = {{t}_{{{\text{ф и н }}}}}$):

Здесь $V(y)$ – вариация функции $y$, определенной на сетке ${{x}_{i}} = - X{\text{/}}2 + ih$, $0 \leqslant i \leqslant N$, а ${{V}_{{ * \rho }}} = {\text{|}}{{\rho }_{{{\text{max}}}}} - 1{\text{|}} + \;{\text{|}}{{\rho }_{{{\text{max}}}}} - 0.5{\text{|}}$ с ${{\rho }_{{{\text{max}}}}} = ma{{x}_{{0 \leqslant i \leqslant N}}}{{\rho }_{i}}$ и ${{V}_{{ * u}}} = U$ – “точные” вариации $\rho $ и $u$.

На фигурах точки с ${{\varepsilon }_{V}} \in [0,0.1]$, $(0.1,0.2],(0.2,\infty ]$ помечены закрашенными, полузакрашенными, незакрашенными кружками (окружностями) соответственно; здесь ${{\varepsilon }_{V}} = \infty $ означает, что расчет не был завершен из-за переполнения или появления отрицательных значений $\rho $. Условно можно считать, что они отвечают устойчивым, промежуточным, неустойчивым решениям соответственно, поскольку графики первых содержат только незначительные одиночные впадинки/пички, к тому же (как было проверено) не растущие во времени и тем самым вполне приемлемые на практике (в сложных расчетах), а графики последних – наоборот, уже заметные осцилляции. Забегая вперед, отметим, что довольно типичная ситуация отражена на фиг. 6. Подчеркнем, что при подсчете вариации учитывается амплитуда каждой впадинки/пичка, причем дважды, так что значение ${{\varepsilon }_{V}} = 0.1$ не так уж велико.

Для схемы стандартного типа (фиг. 3) устойчивые результаты намного (кратно) хуже отвечающих необходимому условию (26), причем соответствие ухудшается с ростом ${{{\text{M}}}_{0}}$ и для ${{{\text{M}}}_{0}} = 8$ становится минимальным (или вовсе исчезает) в зависимости от $\widehat \alpha $. При этом при $\widehat \alpha \geqslant 0.5$ соответствие лучше, чем при $\widehat \alpha = 0.4,0.3$.

Для схем А и особенно Б результаты существенно лучше. Для схемы А (фиг. 4) устойчивые результаты хорошо соответствуют необходимому условию (26) при $\widehat \alpha \geqslant 0.6$ для ${{{\text{M}}}_{0}} = 2,4,6$ (также при $\widehat \alpha = 0.5$ для ${{{\text{M}}}_{0}} = 6$). Для ${{{\text{M}}}_{0}} = 8$ соответствие хуже. Соответствие также намного хуже при $\widehat \alpha = 0.5,0.4$ и особенно 0.3.

Для схемы Б (фиг. 5) устойчивые результаты очень хорошо соответствуют необходимому условию (26) при $\widehat \alpha \geqslant 0.6$ для всех ${{{\text{M}}}_{0}} = 2,4,6,8$ (также при $\widehat \alpha = 0.5$ для ${{{\text{M}}}_{0}} = 2,6$). Но это соответствие резко падает с уменьшением $\widehat \alpha $ от 0.5 и становится плохим при $\widehat \alpha = 0.3$.

Довольно типичные численные решения с ${{\varepsilon }_{V}} \approx 0.1$ и ${{\varepsilon }_{V}} \approx 0.2$ для ${{{\text{M}}}_{0}} = 4$ при $\widehat \alpha = 0.4$ для $t = {{t}_{{{\text{ф и н }}}}}$ по схеме Б показаны на фиг. 6. В случае ${{\varepsilon }_{V}} \approx 0.1$ заметны незначительные одиночные впадинки/пички, а в случае ${{\varepsilon }_{V}} \approx 0.2$ наблюдаются уже заметные осцилляции.

Дополнительно для примера в табл. 1 и 2 приведены значения ${{\varepsilon }_{V}}$ для ${{{\text{M}}}_{0}} = 8$ при различных $\widehat \alpha $ и $\beta = 0.1kmin\{ 2\widehat \alpha ,1{\text{/}}(2\widehat \alpha )\} $ для схем А и Б соответственно. Прочеркам соответствует разрушение численных решений к моменту $t = {{t}_{{{\text{ф и н }}}}}$. Это более детальная информация, чем на фиг. 4г и 5г соответственно. Обратим внимание на то, что ${{\varepsilon }_{V}}$ не убывает по $\beta $ при каждом $\widehat \alpha $ – это естественно. При этом при $\widehat \alpha = 0.6,0.7,0.8,0.9$ происходит резкий финальный скачок от малых значений ${{\varepsilon }_{V}}$ к большим (для схемы Б) или к ситуации с разрушением решения (для схемы А и частично Б). Кроме того, ${{\varepsilon }_{V}}$ быстро убывает по $\widehat \alpha $ при каждом $k$ (с частичным исключением $k = 0.9$ в табл. 2), что довольно неожиданно. В том числе при каждом $k$ значения ${{\varepsilon }_{V}}$ при $\widehat \alpha = 0.7,0.8,0.9$ на несколько порядков меньше, чем при меньших $\widehat \alpha $. Поэтому такие значения $\widehat \alpha $ с успехом могут использоваться на практике.

Полученные результаты служат в пользу применения схем А и особенно Б с энергетически диссипативными дискретизациями по пространству.