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

Уравнения, описывающие волны в трубах с упругими стенками, и численные методы с низкой схемной диссипацией

И. Б. Бахолдин^*

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

Поступила в редакцию 25.06.2019
После доработки 22.01.2020
Принята к публикации 10.03.2020

DOI: 10.31857/S0044466920070030

Аннотация

Рассматриваются уравнения трубы с упругими стенками: труба с контролируемым давлением, труба, наполненная жидкостью, труба с газом. Для описания стенок трубы используются полная модель мембраны и нелинейная теория гиперупругих материалов. Решается задача о распаде произвольного разрыва, решения подтверждают теорию обратимых структур разрывов. Дисперсия коротких волн для данных уравнений исчезает, поэтому допускается включение и диссипативных структур разрывов. В связи со сложным характером уравнений развиваются общие численные методы. Рассматривается применение центрированной трехслойной схемы типа крест и схем, основанных на аппроксимации временных производных по методу Рунге–Кутты различного порядка. Разрабатывается методика коррекции схем на основе метода Рунге–Кутты посредством добавления диссипативных членов. Методы третьего и четвертого порядков в скалярном случае коррекции не требуют. Анализируется возможность использования членов с производными высокого порядка для расчета решений, в которых одновременно присутствуют диссипативные и недиссипативные разрывы. Библ. 18. Фиг. 5.

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

1. ВВЕДЕНИЕ

Исследование распространения волн в трубах важно для технических приложений, а также для изучения процессов в биообъектах. Наиболее простые варианты моделей основаны на уравнениях типа уравнения Кортевега–де Вриза и Буссинеска [1]. Здесь исследуются уравнения с более точным учетом упругих свойств стенок трубы. Основные уравнения были предложены в [2], решения типа уединенных волн для них аналитически исследовались в [3]–[5], а в [6] были проведены расчеты эволюции уединенных волн и распада произвольного разрыва. Изучался и более простой вариант этих уравнений – модель с контролируемым давлением [7], [8]. В [8] было предложено усовершенствовать уравнения путем учета жесткости стенок на изгиб с использованием модели тонкой пластины. В [6] были также учтены сжимаемость материала, диссипация, предложена модель трубы с газом. Ранее была разработана теория бездиссипативных и слабодиссипативных структур разрывов [9]–[11], одним из элементов которой является методика расчетов уравнений с применением численных схем без диссипации или со слабой диссипацией. Проведение расчетов для волн в трубах показало, что данная методика требует дальнейшего развития, в связи с чем было проведено исследование некоторых численных схем. Целью настоящей работы являются анализ свойств уравнений, описывающих волны в трубах, и разработка численных методов с низкой схемной диссипацией, применимых для широкого круга недиссипативных и диссипативных уравнений. В разд. 2 рассматриваются уравнения, описывающие волны в трубах с упругими стенками. Анализируются свойства дисперсионных соотношений $\omega = \omega (\kappa )$, $\omega $ – циклическая частота, $\kappa $ – волновое число, для этих уравнений. Анализ дисперсионных соотношений необходим в связи с тем, что каждому пересечению дисперсионной ветви в начале координат при решении задачи о распаде произвольного разрыва соответствует структура разрыва или центрированная простая волна. В разд. 3 описывается общая методика расчета с помощью численных схем с низкой схемной диссипацией. Проводится анализ схем с применением спектрального метода исследования устойчивости для ряда модельных уравнений с высшими производными. В разд. 4 даются примеры решения задачи о распаде произвольного разрыва для уравнений, описывающих волны в трубах.

2. ОСНОВНЫЕ УРАВНЕНИЯ И ИХ МОДИФИКАЦИИ

2.1. Уравнения для случая контролируемого давления

Уравнения движения волн в трубе с упругими несжимаемыми осесимметричными стенками и фиксированным внутренним и внешним давлением имеют вид [7]

(2.1)

$\left( {R{{\sigma }_{1}}\frac{{z{\text{'}}}}{{\lambda _{1}^{2}}}} \right){\text{'}} - P{\kern 1pt} {\text{*}}rr{\text{'}} = \rho R\ddot {z},$

(2.2)

$\begin{gathered} \left( {R{{\sigma }_{1}}\frac{{r{\text{'}}}}{{\lambda _{1}^{2}}}} \right){\text{'}} - \frac{{{{\sigma }_{2}}}}{{{{\lambda }_{2}}}} + P{\kern 1pt} {\text{*}}z{\text{'}} = \rho R\ddot {r}; \\ {{\lambda }_{1}} = \sqrt {{{r}^{{'2}}} + {{z}^{{'2}}}} ,\quad {{\lambda }_{2}} = \frac{r}{R},\quad {{\lambda }_{3}} = \frac{h}{H},\quad P{\kern 1pt} {\text{*}} = \frac{P}{H},\quad {{\sigma }_{i}} = {{\lambda }_{i}}\frac{{\partial W}}{{\partial {{\lambda }_{i}}}} - p,\quad i = 1,2,3. \\ \end{gathered} $

Штрихом обозначено дифференцирование по переменной $Z$, являющейся начальной лагранжевой пространственной координатой вдоль трубы, точкой обозначено дифференцирование по времени $t$, $H$ и $h$ – толщина стенки трубы в ненапряженном и напряженном состоянии, $W = W({{\lambda }_{1}},{{\lambda }_{2}},{{\lambda }_{3}})$ – упругий потенциал, $p$ – давление, связанное с несжимаемостью материала трубы, аналогичное давлению в несжимаемой жидкости, параметр $P$ – разность между внутренним и внешним давлением, параметр $\rho $ – плотность материала трубы на единицу площади. Для описания упругих свойств стенки трубы используется полная мембранная модель, неизвестные $z$ и $r$ задают поверхность трубы в цилиндрической системе координат, ось $z$ этой системы совпадает с осью трубы, ${{\lambda }_{i}}$ – главные удлинения, а ${{\sigma }_{i}}$ – главные напряжения, компоненты тензора напряжений. Индексы 1, 2, 3 соответствуют широтному (круговому), меридианальному (касательному в плоскости, проходящей через ось вращения) и ортогональному направлению деформируемой поверхности. При отсутствии нагрузки $z = Z$, $r = R$, $h = H$. Модель с контролируемым давлением пригодна в случае заполнения трубы газом пренебрежимо малой плотности. Давление поддерживается за счет перетекания в трубу среды из резервуара большого объема или работы компрессора.

Материал трубы считается несжимаемым, поэтому

${{\lambda }_{1}}{{\lambda }_{2}}{{\lambda }_{3}} = 1,\quad {{\sigma }_{i}} = {{\lambda }_{i}}{{\hat {W}}_{i}},\quad \hat {W}({{\lambda }_{1}},{{\lambda }_{2}}) = W({{\lambda }_{1}},{{\lambda }_{2}},1{\text{/}}({{\lambda }_{1}},{{\lambda }_{2}})),\quad i = 1,2.$

Здесь и далее используется обозначение ${{\hat {W}}_{i}} = \partial{ \hat {W}}{\text{/}}\partial {{\lambda }_{i}}$. В случае сжимаемости также можно найти зависимости ${{\sigma }_{i}} = {{\sigma }_{i}}({{\lambda }_{1}},{{\lambda }_{2}})$, $i = 1,2$ [6].

Потенциал в численных расчетах, проведенных далее, соответствует модели материала Гента, считающейся пригодной для резиноподобных материалов [7]:

$W = - \tilde {\mu }{{J}_{m}}ln\left( {1 - \frac{{\lambda _{1}^{2} + \lambda _{2}^{2} + \lambda _{3}^{2} - 3}}{{{{J}_{m}}}}} \right).$

Возможное использование иного потенциала или учет сжимаемости материала не меняет методику расчета.

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

2.2. Учет жесткости на изгиб

Учет жесткости повышает точность модели, кроме того, он необходим в тех случаях, когда стенка трубы не растянута в продольном направлении, а сжата, поскольку уравнения становятся некорректными. Учет делается на основе уравнения Жермен–Лагранжа для тонкой пластины. Вводится дополнительное “изгибное давление” $P = P + {{P}_{b}}$ [8]

$P \to P + {{P}_{b}},\quad {{P}_{b}} = - \cos \alpha \left[ {\frac{{{{\partial }^{2}}}}{{\partial {{x}^{2}}}}\left( {D\frac{{{{\partial }^{2}}}}{{\partial {{x}^{2}}}}r} \right) - \operatorname{tg} \alpha \frac{{{{\partial }^{2}}}}{{\partial {{x}^{2}}}}\left( {D\frac{{{{\partial }^{2}}}}{{\partial {{x}^{2}}}}z} \right)} \right],\quad \frac{\partial }{{\partial x}} = \frac{{\cos \alpha }}{{z{\text{'}}}}\frac{\partial }{{\partial Z}},\quad \operatorname{tg} \alpha = \frac{{r{\text{'}}}}{{z{\text{'}}}}.$

Здесь $\alpha (Z)$ – угол наклона по отношению к оси $Z$ касательной к кривой, образованной пересечением поверхности трубы с плоскостью, проходящей через ее ось, $x$ – координата в локальной системе координат, ось которой направлена по касательной к рассматриваемой кривой в некоторой точке $({{r}_{0}}(Z),{{z}_{0}}(Z))$. Если геометрическая нелинейность не учитывается, то ${{P}_{b}} = - Dr{\text{''''}}$. Поэтому в случае упрощенного учета жесткости на изгиб в левую часть уравнения (2.2) только добавляется член $ - br{\text{''''}}$, где $b = DR{\text{/}}H = \mu {{H}^{2}}R{\text{/}}3$, $\mu $ – параметр Ламе.

2.3. Уравнения для жидкости

В случае, когда труба заполнена инерционной несжимаемой жидкостью, давление в трубе будет не постоянным, следует добавить уравнения движения жидкости [3]

(2.3)

$\dot {r}z{\text{'}} - r{\text{'}}\dot {z} + vr{\text{'}} + \frac{1}{2}rv{\text{'}} = 0,$

(2.4)

${{\rho }_{f}}(\dot {v}z{\text{'}} - v{\text{'}}\dot {z} + vv{\text{'}}) + P{\text{'}} = 0.$

Здесь $v$ – скорость жидкости, ${{\rho }_{f}}$ – ее плотность. Дисперсионное соотношение содержит четыре дисперсионные ветви, пересекающие начало координат, которые можно сопоставить с продольными упругими и поперечными волнами.

Для применения численных методов удобно исключить давление, выразив его из уравнения (2.2) и сделав последующие исключения временных производных $\dot {r}$ и $\ddot {z}$ с помощью уравнений (2.3) и (2.1) [6]

$P{\kern 1pt} {\text{*}} + P_{b}^{*} = \frac{{P + {{P}_{b}}}}{H} = \frac{{ - 1}}{\xi }\left[ {\frac{{r{\text{'}}\rho R}}{{r\mathop {z'}\nolimits^2 }}\dot {v} + \frac{{\rho R}}{{2\mathop {z'}\nolimits^2 }}\dot {v}{\text{'}} + \mathcal{P}} \right],\quad \xi = 1 + \frac{{\mathop r\nolimits^{'2} }}{{\mathop z\nolimits^{'2} }},$

$\begin{gathered} \mathcal{P} = \frac{1}{{rz{\text{'}}}}\left[ {\left( {R{{\sigma }_{1}}\frac{{r{\text{'}}}}{{\lambda _{1}^{2}}}} \right) - \frac{{r{\text{'}}}}{{z{\text{'}}}}\left( {R{{\sigma }_{1}}\frac{{z{\text{'}}}}{{\lambda _{1}^{2}}}} \right) - \frac{{{{\sigma }_{2}}}}{{{{\lambda }_{2}}}}} \right] - \frac{{\rho R}}{{r{{z}^{{'2}}}}}(q - v)\left( {\frac{{z{\text{'}}q - vr{\text{'}} - \tfrac{1}{2}r{v}{\kern 1pt} '}}{{z'}}} \right) - \\ \, - \frac{{\rho R}}{{r{{z}^{{'3}}}}}\left( {r{\text{'}}q - vr{\text{'}} - \frac{1}{2}rv{\text{'}}} \right)\left( {\frac{1}{2}v{\text{'}} - q{\text{'}}} \right). \\ \end{gathered} $

Уравнения для расчета течений в трубе, заполненной жидкостью, теперь таковы, что в каждом уравнении есть временные и пространственно-временные производные только от одной переменной [6]:

(2.5)

$\begin{gathered} \left[ {\rho _{f}^{*}z{\text{'}} - \left( {\frac{{r{\text{'}}\rho {{R}^{2}}}}{{r{{z}^{{'2}}}\xi }}} \right)} \right]\dot {v} - \left[ {\frac{{r{\text{'}}\rho {{R}^{2}}}}{{r{{z}^{{'2}}}\xi }} + \left( {\frac{{\rho {{R}^{2}}}}{{2{{z}^{{'2}}}\xi }}} \right)} \right]\dot {v}{\text{'}} - \frac{{\rho {{R}^{2}}}}{{2{{z}^{{'2}}}\xi }}\dot {v}{\text{''}} = \\ \, = \rho _{f}^{*}(v{\text{'}}q - vv{\text{'}}) + \left[ { - P_{b}^{*} + \frac{\mathcal{P}}{\xi }} \right] - \{ {{c}_{{v{\text{4}}}}}v{\text{''''}}\} ,\quad \rho _{f}^{*} = \frac{{{{\rho }_{f}}}}{H}, \\ \end{gathered} $

(2.6)

$\dot {r}z{\text{'}} - r{\text{'}}q + vr{\text{'}} + \frac{1}{2}rv{\text{'}} = 0 + \{ {{c}_{{r6}}}r{\text{''''''}}z{\text{'}}\} ,$

(2.7)

$\dot {z} = q - \{ {{c}_{{z4}}}z{\text{''''}}\} ,\quad \rho R\dot {q} = \left( {R{{\sigma }_{1}}\frac{{z{\text{'}}}}{{\lambda _{1}^{2}}}} \right) - (P{\kern 1pt} {\text{*}} + P_{b}^{*})rr{\text{'}}.$

В фигурных скобках здесь добавлены некоторые члены, иногда полезные для улучшения свойств численной схемы при применении метода Рунге–Кутты второго порядка (см разд. 3, 4).

2.4. Уравнения для газа

Анализ методики вывода уравнений для трубы с жидкостью [2] показывает, что в случае заполнения трубы инерционным сжимаемым газом нужно модифицировать уравнение сохранения массы (2.3) [6]:

(2.8)

$({{\rho }_{f}}z{\text{'}} - {{\rho }_{f}}\dot {z}){{r}^{2}} + 2{{\rho }_{f}}r(\dot {r}z{\text{'}} - r{\text{'}}\dot {z}) + ({{\rho }_{f}}v{{r}^{2}}){\text{'}} = 0.$

В изотермическом и адиабатическом случае имеется баротропия, добавляется уравнение состояния вида $P = P({{\rho }_{f}})$, $P{\kern 1pt} {\text{*}} = (P - {{P}_{e}}){\text{/}}H$, ${{P}_{e}}$ – внешнее давление. В этой модели (2.1), (2.2), (2.4), (2.8) шесть дисперсионных ветвей: в основном продольных упругих волн, в основном поперечных упругих волн и газодинамические ветви. Ветви, связанные с распространением в основном поперечных упругих волн, не пересекают начало координат. Такие ветви есть и в случае модели с контролируемым давлением, поэтому такую модель следует рассматривать как упрощенную модель трубы, заполненной газом. Заметим при этом, что дисперсионное соотношение для ветвей продольных упругих волн в модели с контролируемым давлением все же можно получить предельным переходом ${{\rho }_{f}} \to 0$ из дисперсионного соотношения для модели трубы с жидкостью, найденного в [4]. В общем случае потребуются еще включение уравнения энергии, учет термодинамики газа и уравнение состояния с зависимостью давления от температуры.

2.5. Учет диссипации

В случае применения для описания материала трубы вязкоупругой модели типа Кельвина–Фойтха [12] логично добавить вязкие напряжения, например, так [6]:

${{\sigma }_{i}} \to {{\sigma }_{i}} + {{\sigma }_{{vi}}},\quad {{\sigma }_{{vi}}} = {{\nu }_{{vi}}}{{\dot {\lambda }}_{i}},\quad i = 1,2,3.$

Коэффициенты вязкости могут зависеть от деформации. Можно показать, что при некотором упрощенном подходе [6] вязкость материала трубы в случае ее заполнения жидкостью можно учесть, включив в левую часть уравнения (2.5) член ${{c}_{{v2}}}v{\text{''}}$. Таким же способом учитывается и вязкость наполнителя. Можно учитывать трение о стенки, включив в правую часть уравнения (2.1) член $rF{\text{/}}H$, где $F = F(v)$ – сила трения на единицу длины, обычно линейная функция, а в левую часть уравнения (2.4) – член $ - 2F{\text{/}}r$, величина $r{\text{'}}$ предполагается малой [6].

2.6. Упрощенные уравнения

В теории решений с разрывами [12] используется понятие упрощенные уравнения для обозначения гиперболических уравнений, в решениях которых имеются разрывы. Уравнения с дисперсией или вязкостью, в решениях которых этим разрывам соответствуют непрерывные переходы, структуры разрывов, называются полными уравнениями. Упрощенные уравнения для стенок трубы получаются путем удаления старших производных из уравнения для $r$ и исключения давления из уравнения для $z$

(2.9)

$\dot {u} = q{\text{'}};\quad \rho R\dot {q} = \left( {R{{\sigma }_{1}}\frac{u}{{\lambda _{1}^{2}}}} \right){\text{'}} - P{\kern 1pt} {\text{*}}r{\text{'}},\quad {{\lambda }_{1}} = u = z{\text{'}},\quad P{\kern 1pt} {\text{*}} = \frac{{{{\sigma }_{2}}}}{{{{\lambda }_{2}}ru}}.$

Уравнения для жидкости и газа упрощения не требуют и добавляются к (2.9).

Таким образом, распространение волн в трубах описывается уравнениями с производными различных порядков. Уравнения обратимы по пространству и времени. Во всех этих уравнениях скорость распространения волн конечна, но только последние из представленных уравнений в корректном случае являются гиперболическими. Остальные – уравнения с дисперсией, исчезающей при $\kappa \to \infty $ (кроме уравнений с учетом жесткости на изгиб или с учетом вязкости).

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

Ранее была разработана теория обратимых и слабодиссипативных разрывов [9]– [11]. Одна из целей данной работы – приложение этой теории к уравнениям, описывающим волны в трубах с упругими стенками. Согласно теории, решения задачи о распаде произвольного разрыва в недиссипативном случае содержат расширяющиеся со временем области однородных участков, волновых зон и центрированных простых волн. Эти области могут разделяться обратимыми структурами разрывов (кинками), представляющими локализованные переходы между однородными, периодическими или стохастическими состояниями. Упорядоченные волновые зоны описываются усредненными уравнениями, а центрированные простые волны – упрощенными уравнениями. В слабодиссипативном случае волновые зоны также описываются усредненными уравнениями, но со временем перестают расширяться. Решения в недиссипативном и слабодиссипативном случае качественно отличаются. В связи с этим требуется применение недиссипативных численных схем либо схем с пренебрежимо малой схемной диссипацией. Схемная обратимая (недиссипативная) дисперсия при достаточно малых пространственных и временных шагах не оказывает влияния на качественный вид решений, не меняются усредненные параметры упорядоченных волновых зон и типы структур разрывов, поскольку она суммируется с дисперсией исследуемых уравнений. Качественные схемные эффекты, связанные со схемной дисперсией, могут проявляться на грубых сетках, но они устраняются путем уменьшения пространственного шага. По возможности использовались явные схемы с центральными пространственными разностями. В случае наличия в уравнениях членов со смешанными пространственно-временными производными первого и второго порядка по времени только эти члены аппроксимировались неявным образом, что позволяло применять метод прогонки. Для моделирования решений задачи Коши посредством расчетов в конечной области на достаточно удаленных границах ставились условия постоянства рассчитываемых неизвестных в стольких краевых узлах, сколько необходимо для реализации численной схемы. При необходимости область расчета расширялась. Иногда вводились диссипативные зоны вблизи границ для поглощения приходящих волн. В случае наличия в системе уравнений производных второго порядка по времени система преобразовывалась к уравнениям с производными первого порядка.

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

(3.1)

$\frac{{\partial u}}{{\partial t}} = \mathcal{L}u,\quad \mathcal{L} = \sum\limits_{j = 1}^M \,{{a}_{j}}{{\mathcal{L}}_{j}},\quad {{\mathcal{L}}_{j}} = \frac{{{{\partial }^{j}}u}}{{\partial {{x}^{j}}}},\quad \omega = \sum\limits_{j = 1}^M \,{{a}_{j}}{{i}^{{j - 1}}}{{\kappa }^{j}}.$

Производным нечетного порядка соответствуют обратимые дисперсионные члены, а четным – диссипативные. Диссипативный член корректен, если ${{a}_{{4m - 2}}} > 0$, а ${{a}_{{4m}}} < 0$, $m = 1,2,...$. В соответствии с методикой определения корректности уравнений по дисперсионному соотношению [13], если диссипативный член со старшей производной является некорректным диссипативным членом, то уравнение становится некорректным, в иных случаях включение некорректных диссипативных членов может вести только к неустойчивости для некоторых диапазонов волн.

Для пространственных производных применяются центральные разности, аппроксимация второго порядка точности:

${{u}_{x}} \to \frac{{u_{{k + 1}}^{n} - u_{{k - 1}}^{n}}}{{2\tilde {h}}},\quad {{u}_{{xx}}} \to \frac{{u_{{k + 1}}^{n} + u_{{k - 1}}^{n} - 2u_{k}^{n}}}{{{{{\tilde {h}}}^{2}}}},\quad {{u}_{{xxx}}} \to \frac{{u_{{k + 2}}^{n} - 2u_{{k + 1}}^{n} + 2u_{{k - 1}}^{n} - u_{{k - 2}}^{n}}}{{2{{{\tilde {h}}}^{3}}}},$

${{u}_{{xxxx}}} \to \frac{{u_{{k + 2}}^{n} - 4u_{{k + 1}}^{n} + 6u_{k}^{n} - 4u_{{k - 1}}^{n} + u_{{k - 2}}^{n}}}{{{{{\tilde {h}}}^{4}}}},$

${{u}_{{xxxxx}}} \to \frac{{u_{{k + 3}}^{n} - 4u_{{k + 2}}^{n} + 5u_{{k + 1}}^{n} - 5u_{{k - 1}}^{n} + 4u_{{k - 2}}^{n} - u_{{k - 3}}^{n}}}{{2{{{\tilde {h}}}^{5}}}},...$

При использовании спектрального метода исследования устойчивости делается подстановка $u_{k}^{n} = {{q}^{n}}expik\phi $, при $\left| q \right| > 1$ имеет место рост возмущений.

3.1. Центрированная трехслойная схема

Используется следующая аппроксимация производной по времени: $(u_{k}^{{n + 1}} - u_{k}^{{n - 1}}){\text{/}}(2\tau )$. Схема в устойчивых случаях для недиссипативных уравнений (3.1) не обладает схемной диссипацией, ${\text{|}}q{\text{|}} = 1$. Спектральный метод показывает, что для уравнения (3.1) условие устойчивости $\tau < c{{\tilde {h}}^{M}}$, $c$ – некоторая постоянная, в случае, если все производные нечетные, что соответствует скорости роста $\omega (\kappa )$ при $\kappa \to \infty $. В случае наличия четных производных для обеспечения устойчивости приходится применять иные схемы. Можно делать аппроксимацию по нижнему слою. Если член со старшей производной диссипативный, то условие устойчивости при этом также $\tau < c{{\tilde {h}}^{M}}$, что соответствует скорости убывания амплитуды волн при $\kappa \to \infty $. Можно также применять схему типа схемы Дюфорта–Франкеля для уравнения теплопроводности (шаблон ромб), где в разностной аппроксимации центральный член $u_{k}^{n}$ заменяется полусуммой $(u_{k}^{{n + 1}} + u_{k}^{{n - 1}}){\text{/}}2$, но при этом вместо ${{\partial }^{{2j}}}u{\text{/}}\partial {{x}^{{2j}}}$ аппроксимируется ${{\partial }^{{2j}}}u{\text{/}}\partial {{x}^{{2j}}} + {{\tau }^{2}}{\text{/}}{{\tilde {h}}^{{2j}}}{{\partial }^{2}}u{\text{/}}\partial {{t}^{2}}$. Такая аппроксимация не создает ограничений по устойчивости, но для достижения точности величину $\tau $ необходимо брать достаточно малой. В случае производной второго порядка удобно применять также неявную схему и метод прогонки. Численный эксперимент показывает, что и в случае не скалярных уравнений при выборе условия устойчивости можно ориентироваться на скорость роста $\omega (\kappa )$ или на скорость убывания амплитуды волн. Центрированная трехслойная схема успешно применялась для многих бездиссипативных систем уравнений [9] и, в частности, для случая уравнений трубы с контролируемым давлением [8]. В случае применения для уравнений трубы, заполненной жидкостью или газом, обнаружилась сильная краевая неустойчивость. Она может быть устранена посредством введения вблизи границы расчетного объема в уравнении для $v$ диссипативных зон c добавлением члена ${{c}_{{v2}}}v{\text{''}}$ с аппроксимацией по шаблону ромб. Коэффициент при диссипативном члене должен быть тем больше, чем меньше пространственный шаг, т.е. вблизи границы фактически решается уравнение с производной второго порядка по времени. Ввиду неудобства такого подхода и невозможности решения краевых задач в этих случаях применялись двухслойные схемы с аппроксимацией временных производных по методу Рунге–Кутты. Число ветвей величины $q$ у трехслойной схемы в два раза больше числа ветвей дисперсионного соотношения рассчитываемых уравнений, фактически при расчетах задается в два раза больше начальных данных. Не всегда все ветви $q$ аппроксимируют дисперсионные ветви рассчитываемого уравнения. Это ведет к неустойчивости, например, для уравнения теплопроводности. Возможно, граничная неустойчивость при трехслойной схеме связана с недостаточностью числа граничных условий. Она наблюдалась и при расчете уравнений электронной магнитной гидродинамики плазмы, устранялась преобразованием уравнений от вида с пространственными производными первого порядка к виду с производными третьего порядка [9], [14].

3.2. Схемы на основе аппроксимации временных производных методом Рунге–Кутты

Имеются различные многошаговые реализации методов Рунге–Кутты высокого порядка, см., например, [15]. Будем считать, что не используются неявные аппроксимации, что в результирующей формуле метода порядок полинома от $\tau $ и число шагов совпадают с порядком метода и что методы получены классическим способом приравнивания коэффициентов соответствующего порядка аппроксимации к нулю [16]. Тогда в случае линейного уравнения одношаговое операторное представление линейной схемы с аппроксимацией временной производной по методу Рунге–Кутты имеет вид разложения $expL$ с точностью до соответствующего порядка:

(3.2)

${{u}^{{n + 1}}} = \left( {1 + L + \frac{{{{L}^{2}}}}{{2!}} + \frac{{{{L}^{3}}}}{{3!}} + \frac{{{{L}^{4}}}}{{4!}} + ...} \right){{u}^{n}},$

где $L{\text{/}}\tau $ – разностный оператор, аппроксимирующий дифференциальный оператор $\mathcal{L}$. Пусть $F(L)$ – результат подстановки ${{e}^{{ik\phi }}}$ в пространственно обратимый оператор, например,

$\begin{gathered} F({{L}_{1}}) = (\tau {\text{/}}\tilde {h})isin\phi ,\quad F({{L}_{2}}) = - 4(\tau {\text{/}}{{{\tilde {h}}}^{2}})si{{n}^{2}}(\phi {\text{/}}2),\quad F({{L}_{3}}) = (\tau {\text{/}}{{{\tilde {h}}}^{3}})i(sin2\phi - 2sin\phi ),\; \ldots \\ F({{L}_{{2m - 1}}}) = (\tau {\text{/}}{{{\tilde {h}}}^{{2m - 1}}})i{\kern 1pt} \sum\limits_{j = 1}^m {{{{\tilde {c}}}_{{mj}}}sinj\phi } ,\quad {{{\tilde {c}}}_{{mj}}} \in R,\quad F({{L}_{{2m}}}) = {{( - 1)}^{m}}{{2}^{{2m}}}(\tau {\text{/}}{{{\tilde {h}}}^{{2m}}})si{{n}^{{2m}}}(\phi {\text{/}}2). \\ \end{gathered} $

Справедливы алгебраические законы: $F(A + B) = F(A) + F(B)$, $F(AB) = F(A)F(B)$. Для недиссипативных операторов

$\begin{gathered} {\text{|}}q{\kern 1pt} {{{\text{|}}}^{2}} = \mathop {\left[ {1 + F\left( {\frac{{{{L}^{2}}}}{{2!}} + \frac{{{{L}^{4}}}}{{4!}} + ...} \right)} \right]}\nolimits^2 + \mathop {\left[ {F\left( {L + \frac{{{{L}^{3}}}}{{3!}} + \frac{{{{L}^{5}}}}{{5!}} + ...} \right){\text{/}}i} \right]}\nolimits^2 , \\ {\text{|}}q{\kern 1pt} {{{\text{|}}}^{2}} = F\mathop {\left( {1 + \frac{{{{L}^{2}}}}{{2!}} + \frac{{{{L}^{4}}}}{{4!}} + ...} \right)}\nolimits^2 - F\left( {\mathop {\left( {L + \frac{{{{L}^{3}}}}{{3!}} + \frac{{{{L}^{5}}}}{{5!}} + ...} \right)}\nolimits^2 } \right). \\ \end{gathered} $

Перейдем к неотрицательной функции $G$, $G({{L}^{2}}) \equiv - F({{L}^{2}})$, $G({{L}^{{2p}}}) = {{( - 1)}^{{p + 1}}}G{{({{L}^{2}})}^{p}}$

$\begin{gathered} {\text{|}}q{\kern 1pt} {{{\text{|}}}^{2}} = 1 - \frac{{2G\left( {{{L}^{2}}} \right)}}{{2!}} + \frac{{G{{{({{L}^{2}})}}^{2}}}}{{2{{!}^{2}}}} + \frac{{2G{{{({{L}^{2}})}}^{2}}}}{{4!}} - \frac{{G{{{({{L}^{2}})}}^{3}}}}{{4!}} + \frac{{G{{{({{L}^{2}})}}^{4}}}}{{4{{!}^{2}}}} + \ldots \\ + \;G({{L}^{2}}) - \frac{{2G{{{({{L}^{2}})}}^{2}}}}{{3!}} + \frac{{G{{{({{L}^{2}})}}^{3}}}}{{3{{!}^{2}}}} + \frac{{2G{{{({{L}^{2}})}}^{3}}}}{{5!}} - \frac{{2G{{{({{L}^{2}})}}^{4}}}}{{3!5!}} + \frac{{G{{{({{L}^{2}})}}^{5}}}}{{5{{!}^{2}}}} + \ldots \\ \end{gathered} $

Приведем результаты вычисления величины $\left| q \right|$ для методов Рунге–Кутты 1, 2, 3, 4, 5 порядка:

$\begin{gathered} {\text{|}}q{\kern 1pt} {{{\text{|}}}^{2}} = 1 + G({{L}^{2}});\quad {\text{|}}q{\kern 1pt} {{{\text{|}}}^{2}} = 1 + \frac{{G{{{({{L}^{2}})}}^{2}}}}{4};\quad {\text{|}}q{\kern 1pt} {{{\text{|}}}^{2}} = 1 - G{{({{L}^{2}})}^{2}}\left( {\frac{1}{{12}} - \frac{{G({{L}^{2}})}}{{3{\kern 1pt} {{!}^{2}}}}} \right); \\ {\text{|}}q{\kern 1pt} {{{\text{|}}}^{2}} = 1 - G{{({{L}^{2}})}^{3}}\left( {\frac{1}{{72}} - \frac{{G({{L}^{2}})}}{{4{{!}^{2}}}}} \right);\quad {\text{|}}q{\kern 1pt} {{{\text{|}}}^{2}} = 1 + G{{({{L}^{2}})}^{3}}\left[ {\frac{1}{{360}} - G({{L}^{2}})\left( {\frac{1}{{960}} - \frac{{G({{L}^{2}})}}{{5{\kern 1pt} {{!}^{2}}}}} \right)} \right]. \\ \end{gathered} $

В общем случае ${\text{|}}q{\kern 1pt} {{{\text{|}}}^{2}} = 1 + G{{({{L}^{2}})}^{s}}({{C}_{*}} + O(G({{L}^{2}})))$. Параметр $s$, определяющий скорость роста или убывания возмущений, равен 1 для метода первого порядка, равен 2 для методов второго и третьего порядка и равен 3 для методов четвертого и пятого порядка, для методов шестого и седьмого порядков он равен 4, и т.д. Для методов 1–5 порядка ${{C}_{*}} = 1;\;1{\text{/}}4;\; - {\kern 1pt} 1{\text{/}}12;\; - {\kern 1pt} 1{\text{/}}72;\;1{\text{/}}360$, ${{C}_{*}} > 0$ для методов 6–7 порядка. В случае нелинейных уравнений явные методы Рунге–Кутты выше четвертого порядка содержат большее число шагов, чем порядок метода [16], поэтому при анализе линеаризованных вариантов этих методов в (3.2) добавляются более высокие степени оператора $L$ с иными коэффициентами, значения ${{C}_{ * }}$ другие. Согласно спектральному методу численная схема устойчива, если ${\text{|}}q{\kern 1pt} {\text{|}} < 1 + c\tau $, $c > 0$ [17]. Поэтому рассматриваемые схемы устойчивы для $\mathcal{L} = {{a}_{{2m - 1}}}{{\mathcal{L}}_{{2m - 1}}}$ при условии $\tau < c{{\tilde {h}}^{{2s(2m - 1)/(2s - 1)}}}$, $c$ – произвольная константа. Но в случае ${{C}_{*}} > 0$ надежный расчет с незначительным ростом возмущений при фиксированном $\tau $ возможен только на конечном интервале $0 < t < T$, $T \sim {{\tau }^{{ - (2s - 1)}}}$, в дальнейшем произойдет аварийная остановка из-за роста возмущений. Продолжительность этого интервала резко растет с увеличением порядка метода. Схема первого порядка мало эффективна из-за обременительного условия устойчивости и быстрого роста возмущений, тем не менее она применялась для контроля при расчетах электронной магнитной гидродинамики [9]. Результаты применения схемы второго порядка на практике мало отличаются от результатов применения трехслойной схемы в устойчивых случаях. Для методов третьего и четвертого порядка роста возмущений нет при естественных условиях устойчивости $\tau < c{{\tilde {h}}^{{2m - 1}}}$; $c = {\text{|}}{{C}_{*}}{{a}_{{2m - 1}}}{{{\text{|}}}^{{ - 1/2}}}$ при $m = 1$. Можно найти диссипативный член с низшей производной дифференциального приближения схемы методом восстановления по дисперсионному соотношению ${\text{|}}q{\kern 1pt} {\text{|}} \approx 1 + {{C}_{*}}{{({{C}_{m}}\tau {{a}_{{2m - 1}}}{{\phi }^{{2m - 1}}}{\text{/}}{{\tilde {h}}^{{2m - 1}}})}^{{2s}}}{\text{/}}2$ при $\phi \to 0$; ${{C}_{1}} = 1$, ${{C}_{2}} = - ({{2}^{3}} - 2){\text{/}}3!$, ${{C}_{3}} = ({{3}^{5}} - 4 \cdot {{2}^{5}} + 5){\text{/}}5!$, … . Такой показатель роста или убывания можно получить, если добавить диссипативный член

${{c}_{d}}{{\partial }^{{2s(2m - 1)}}}u{\text{/}}{{\partial }^{{2s(2m - 1)}}}x,\quad {{c}_{d}} = {{C}_{*}}{{( - 1)}^{{s(2m - 1) + 1}}}{{({{C}_{m}}\tau {{a}_{{2m - 1}}})}^{{2s}}}{\text{/}}(2\tau ),\quad T \sim {\text{|}}{{c}_{d}}{\kern 1pt} {{{\text{|}}}^{{ - 1}}}.$

Такой же результат можно получить непосредственно из дифференциального приближения разностной схемы, но это более трудоемко. Здесь подразумевается, что дифференциальное приближение приводится к виду (3.1) эквивалентными преобразованиями с сохранением порядка аппроксимации при условии $\tau \sim \tilde {h}$.

В случае роста возмущений численную схему можно улучшить, добавив к $\mathcal{L}$ корректный диссипативный оператор ${{\mathcal{L}}_{d}} = {{a}_{{2j}}}{{\mathcal{L}}_{{2j}}}$. При этом рост или убывание возмущений определяются в основном знаком выражения, ${{C}_{*}}G{{({{L}^{2}})}^{s}} + 2F({{L}_{d}})$, если ${{a}_{{2m - 1}}}\tau {\text{/}}{{\tilde {h}}^{{2m - 1}}} \ll 1$ и ${{a}_{{2j}}}\tau {\text{/}}{{\tilde {h}}^{{2j}}} \ll 1$. Метод Рунге–Кутты обычно реализуется как многошаговый. Для ускорения счета можно ограничиться добавлением диссипативного оператора только на заключительном шаге. Коррекция для всех волн, а не только коротких, возможна при $j \leqslant s(2m - 1)$, оптимальный вариант с полной коррекцией и наименьшим значением коэффициента достигается при $j = s(2m - 1)$. В этом случае при выборе коэффициента ${{a}_{{2j}}}$ можно ориентироваться на величину ${{c}_{d}}$. В общем случае коэффициент ${{a}_{{2j}}}$ должен иметь порядок ${{\tau }^{{2s - 1}}}{{\tilde {h}}^{{2j - 2s(2m - 1)}}}$. Если оператор $\mathcal{L}$ содержит несколько членов, то при выборе порядка производных в диссипативном операторе придется ориентироваться на член с низшей производной или добавлять несколько корректирующих членов. Здесь есть существенное отличие от коррекции дифференциального уравнения корректность определяется диссипативным членом с производной наиболее высокого порядка и откорректировать уравнение добавлением члена с производной более низкого порядка невозможно. При применении методов с аппроксимацией по времени высокого порядка даже в случае наличия роста возмущений этот рост медленный и коррекция численной схемы может не требоваться.

В случае, если уравнение корректно и содержит только диссипативные члены, условие устойчивости при применении этой схемы $\tau < c{{\tilde {h}}^{M}}$, $M$ – порядок производной старшего диссипативного члена, $c$ – некоторая константа.

Для методов первого и второго порядков было показано, что сделанные оценки порядка роста возмущений, условия устойчивости и порядка корректирующего диссипативного члена справедливы и для обобщенного волнового уравнения

$\frac{{{{\partial }^{2}}u}}{{\partial {{t}^{2}}}} = {{( - 1)}^{{m + 1}}}{{c}_{{2m}}}\frac{{{{\partial }^{2}}u}}{{\partial {{x}^{{2m}}}}},\quad {{c}_{{2m}}} > 0;\quad {{\omega }^{2}} = {{c}_{{2m}}}{{\kappa }^{{2m}}},\quad m = 1,2,3,...{\kern 1pt} {\kern 1pt} ,$

при замене $2m - 1$ на $m$. Уравнения (2.2) и (2.1) аналогичны этому уравнению. Для расчетов используется форма записи с производными по времени первого порядка:

(3.3)

$\begin{gathered} \frac{{\partial u}}{{\partial t}} = {{b}_{l}}\frac{{{{\partial }^{l}}v}}{{\partial {{x}^{l}}}} - \left\{ {{{{( - 1)}}^{j}}{{a}_{{2j}}}\frac{{{{\partial }^{{2j}}}u}}{{\partial {{x}^{{2j}}}}}} \right\},\quad \frac{{\partial v}}{{\partial t}} = b_{{2m - l}}^{*}\frac{{{{\partial }^{{2m - l}}}u}}{{\partial {{x}^{2}}}} - \left\{ {{{{( - 1)}}^{j}}a_{{2j}}^{*}\frac{{{{\partial }^{{2j}}}v}}{{\partial {{x}^{{2j}}}}}} \right\}, \\ {{( - 1)}^{{m + 1}}}{{c}_{{2m}}} = {{b}_{l}}b_{{2 - l}}^{*},\quad {{a}_{{2j}}} > 0,\quad a_{{2j}}^{*} > 0,\quad l = 0,1,...,2m,\quad j = 1,2,...\;. \\ \end{gathered} $

При $m = 1$, $l = 1$ имеем аналог газодинамических уравнений. Добавочные корректирующие диссипативные члены даны в фигурных скобках.

В связи с тем, что число ветвей $q$ совпадает с числом ветвей дисперсионного соотношения можно предполагать, что в случае наличия аппроксимации эти ветви аппроксимируют ветви дисперсионного соотношения и рассматривать схемы с применением методов Рунге–Кутты как универсальные, в том числе и для уравнений с диссипацией.

Схема типа Лакса–Вендроффа в версии Бернстайна [18] отличается от схемы с аппроксимацией по методу Рунге–Кутты второго порядка только тем, что предиктор в них вычисляется с применением схемы Лакса, в которой используется пространственное усреднение. Для уравнения переноса одношаговое представление этой схемы представляет собой фактически схему с аппроксимацией временной производной первого порядка с добавленным вязким членом c пространственной производной второго порядка и коэффициентом порядка $\tau $. То есть это и есть скорректированная схема. Анализ дифференциальной формы для уравнения переноса также показывает, что в схемах типа Лакса–Вендроффа остаточные диссипативные члены корректны, описываются производными четвертого порядка, а их коэффициенты имеют порядок $O(\tau \mathop {\tilde {h}}\nolimits^2 ,{{\tau }^{3}})$. Уменьшать пространственный шаг сложнее, чем уменьшать временной шаг, поэтому схемная вязкость для таких схем оказывается существенно большей, чем при использовании метода Рунге–Кутты. При наличии в исходных уравнениях производных второго порядка, диссипативных членов в случае схем типа Лакса–Вендроффа, как и в случае трехслойной схемы, требуются специальные аппроксимации. При применении схемы Лакса–Ведроффа для уравнений волн в трубах были замечены проблемы с устойчивостью.

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

3.3. Методы контроля

Дисперсия схемы может вызвать “лишние” волны, но их легко отличить от настоящих путем изменения шага сетки. В случае существования классического решения “лишние” волны легко устраняются уменьшением пространственного шага. У нелинейных гиперболических уравнений дисперсии нет, классического решения часто не существует и устранить схемные волны таким способом нельзя. Длина волны “лишних” волн определяется шагом сетки, это 3–12 шагов в зависимости от вида уравнений и используемой численной схемы, т.е. с уменьшением шага пропорциональным образом уменьшается и длина волны. Амплитуда таких волн при наличии классического решения убывает экспоненциальным образом при уменьшении шага, поэтому от них легко избавиться. Длина “настоящих” волн для сходящихся схем при уменьшении шага стремится к некоторой конечной величине, к конечным величинам стремятся также амплитуда волны и средний уровень. Это позволяет предположить сходимость для огибающей усредненных уравнений упорядоченных волновых зон [9]. Удовлетворительное приближение, без возникновения выраженных схемных волн, обычно достигается примерно начиная от 20 шагов. При непродолжительных по времени расчетах не сложно проверить и сходимость поточечно. Существуют еще и методы проверки на основе аналитических исследований дисперсионных систем [9]–[11].

3.4. Применение высших производных для получения решений с разрывами

Расчеты магнитной гидродинамики показали, что рассматриваемая схема на основе метода Рунге–Кутты позволяет рассчитывать и нелинейные гиперболические уравнения, в частности уравнения (2.9), без дисперсии, только придется добавить вязкость. Если ее не будет, то схема вместо диссипативной структуры разрыва может автоматически вставить недиссипативную структуру, которая качественно по виду огибающей не отличается от физических структур дисперсионной системы электронной магнитной гидродинамики, но длина волн в ней зависит от шага сетки согласно п. 3.3. При этом было обнаружено, что по крайней мере для разрывов малой амплитуды граничные условия на разрывах не меняются при использовании схемной дисперсии, реальной дисперсии или реальной диссипации. Уравнения, приведенные в разд. 1, относятся к типу уравнений с фазовой скоростью, стремящейся к постоянному значению при $\kappa \to \infty $, т.е. с исчезающей дисперсией при стремлении длины волны к нулю, исключение – уравнения с учетом жесткости на изгиб. Именно уравнения такого типа наиболее адекватно описывают реальные физические явления. Можно предположить, что дисперсия в случае сильной нелинейности не способна устранить опрокидывание волны, т.е. для разрывов малой амплитуды возникает недиссипативная структура, а для большой – опрокидывание волны (данная ситуация не встречалась в случае уравнений Кортевега–де Вриза, нелинейного уравнения Шрёдингера и им подобных уравнений, у которых дисперсия растет при уменьшении длины волны). Для уравнений трубы опрокидывание упругих волн является вполне физическим допустимым [6], в случае контролируемого давления оно описывается теми же уравнениями (2.2), (2.1), при наличии жесткости на изгиб уравнения при этом будут корректными, но для расчета понадобятся специальные численные методы, поскольку в нескольких точках производные обращаются в бесконечность. Опрокидывание, как и в случае получения обобщенных решений нелинейных гиперболических уравнений, можно устранить, добавив дополнительные члены с диссипацией, коэффициенты при которых стремятся к нулю при уменьшении $\tilde {h}$. В случае выявления при расчете уравнений трубы возникновения коротких волн, длина которых существенно уменьшается при уменьшении шага сетки, и возможного последующего вычислительного блоуапа, предлагается добавлять диссипацию, блоуапа можно также избежать, добавив дисперсионный член, например, включив жесткость на изгиб упрощенным способом, см. разд. 1. Такой же подход можно применять, если в уравнениях дисперсия или вязкость слишком малы и не представляется возможным сделать достаточно малый пространственный шаг. При этом в случае решения задачи о распаде разрыва искаженным будет только участок структуры разрыва, ассоциированной с той дисперсионной ветвью, для которой возникает опрокидывание, что позволит исследовать оставшиеся структуры разрыва. Если же требуется исследование именно этой структуры, то коэффициенты при добавочных членах должны быть физически обоснованными и постоянными, при выборе пространственного шага для устранения схемных эффектов нужно руководствоваться п. 2.3.

Вязкость, описываемая производными высокого порядка, обладает сильной избирательностью по длине волны. Поэтому она позволяет выбирать коэффициент при добавленном члене малым и проводить расчеты, в которых в течение длительного времени присутствуют одновременно диссипативные и бездиссипативные структуры разрывов. Так, для уравнений типа нелинейного уравнения переноса ${{u}_{t}} + f{{(u)}_{x}} = 0$, исходя из сохранения сеточной длины (числа узлов) структуры разрыва при изменении пространственного шага и временного шага при условии $\tau = c\tilde {h}$ необходимый порядок коэффициента вязкого или невязкого добавленного члена нужно брать равным ${{\tilde {h}}^{{l - 1}}}$, где $l$ – порядок производной. При использовании нелинейной вязкости Неймана и ее обобщений с высшими производными величина $m$ может быть нечетной, учитывается суммарный порядок производных, так, члену ${{({\text{|}}{{u}_{x}}{\text{|}}{{u}_{x}})}_{x}}$ соответствует $l = 3$. В случае добавки дополнительного дисперсионного члена в волновое уравнение ${{u}_{{tt}}} = f{{(u)}_{{xx}}}$ порядок коэффициента целесообразно брать равным $\mathop {\tilde {h}}\nolimits^{l - 2} $, сохранение сеточной длины волны при $\tau = c\tilde {h}$ достигается в случае применения для аппроксимации ${{u}_{{tt}}}$ трехслойной схемы или же в стационарном случае при применении двухслойной схемы к уравнениям типа (3.3). Поскольку коэффициенты определяются шагом сетки, то такие добавки можно охарактеризовать как искусственную схемную вязкость и дисперсию. Такие добавки не требуют изменения естественного условия устойчивости $\tau < c\tilde {h}$. В случае бездиссипативных добавок с наиболее низким порядком производных $l = 3$ для уравнения переноса и $l = 4$ для волнового уравнения, коэффициенты имеют порядок ${{\tilde {h}}^{2}}$, эти добавки аналогичны низшим членам дифференциального приближения, т.е. порядок аппроксимации схемы не изменился. В случае использования производных более высокого порядка можно и не понижать порядок коэффициентов (такой подход к коррекции схемы предлагался в [6]), но тогда условие устойчивости станет более обременительным. Применение добавок с высшими производными позволяет также рассчитывать решения со структурами слабых разрывов (разрывов производных).

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

Часто применяемые для гиперболических уравнений типа газовой динамики численные схемы с использованием коррекции на основе принципа максимума вместо вязкости не применимы в случае недиссипативных дисперсионных систем, появление новых волн – естественное свойство недиссипативных систем с дисперсией общего вида $\partial \omega {\text{/}}\partial \kappa \ne 0$, $\omega \in R$, $\kappa \in R$, принцип максимума для них не действует.

4. ЗАДАЧА О РАСПАДЕ ПРОИЗВОЛЬНОГО РАЗРЫВА

Были проведены серии расчетов задачи о распаде произвольного разрыва. Для случая контролируемого давления и трубы с жидкостью они подробно описаны в [8] и [6]. В качестве начальных данных для $r$ и $z{\text{'}}$ бралась сглаженная ступенька. В случае заполнения жидкостью или газом начальные скорость и плотность считались не зависящими от $Z$. Полученные результаты соответствуют теории обратимых структур разрывов. Общие данные: $R = 1$, $H = 1$ (расчет делается с использованием безразмерных неизвестных, величина данного параметра становится существенна только при вычислении коэффициента жесткости на изгиб [6], если жесткость учитывается), $\rho = 1$, $\tilde {\mu } = 2$, ${{J}_{m}} = 30$, $z{\text{'}} = 1.1$ при $Z \to - \infty $.

На фиг. 1 показан пример решения задачи о распаде разрыва для трубы с контролируемым давлением. Решалась система уравнений (2.1), (2.2), $t = 3000$. Использована центрированная трехслойная схема. Однородные участки разделены участками с в основном продольными упругими волнами, бегущими влево, и с в основном продольными упругими волны, бегущими вправо. Переходных элементов, связанных с поперечными волнами, нет, поскольку ветвь поперечных волн не пересекает начало координат и возможны поперечные колебания трубы как единого целого. Показаны графики решений с различными радиусами деформированной трубы с правой стороны. Имеются расширяющиеся солитонные (при $t \to \infty $ волна на границе волновой зоны стремится к уединенной волне) структуры разрывов разного типа с излучением волн в разных направлениях, иначе говоря с положительной и отрицательной дисперсией. При непрерывном изменении радиуса справа выявляется нестационарная структура промежуточного типа с двумя волновыми зонами, представленная в [8].

Фиг. 1.

На фиг. 2 показаны графики решений задачи о распаде произвольного разрыва для трубы, заполненной жидкостью при разных значениях плотности, ${{\rho }_{f}} = 10$ – кривая 1, ${{\rho }_{f}} = 0.25$ – кривая 2, с применением метода Рунге–Кутты второго порядка с коррекцией, $t = 1200$, $v = {{v}_{0}} = 0.4$ при $t = 0$, ${{c}_{{v4}}} = {{c}_{{z4}}} = {{10}^{{ - 5}}}$, $\tilde {h} = 0.05$. Решалась система уравнений (2.5)–(2.7). Неявное разностное уравнение разрешалось методом прогонки. Заметим, что сходная система уравнений с пространственно-временными производными была и при расчетах электронной магнитной гидродинамики плазмы [14]. Наблюдаются (слева направо) в основном продольные упругие волны, бегущие влево, классический кинк и в основном поперечные упругие волны, бегущие вправо. Ветвь поперечных волн не пересекает начало координат, колебания трубы как единого целого невозможны. Обнаружено, что при уменьшении плотности жидкости выход на решения с однородными участками замедляется или не происходит, т.е. упрощенные уравнения плохо описывают эволюцию решения при малых значениях плотности. Заметим, что это не наблюдается для трубы с газом. Упрощенные уравнения для трубы с жидкостью можно вывести непосредственно из уравнений (2.5)–(2.7), отбросив все высшие производные. Но при ${{\rho }_{f}} \to 0$ в уравнении (2.5) члены с высшими производными очевидно начинают играть существенную роль, и модель не переходит в модель с контролируемым давлением. Заметим также, что при ${{\rho }_{f}} = 10$ структура разрыва для продольных волн, бегущих вправо, фактически рассчитана как диссипативная (сравним добавочные коэффициенты с рекомендацией п. 2.4: ${{c}_{4}} = \tilde {C}{{\tilde {h}}^{3}}$, $\tilde {C} = 0.08$), поскольку она не расширяется со временем, но структура поперечных волн рассчитана как бездиссипативная.

Фиг. 2.

На фиг. 3 показан пример расчета для трубы, наполненной газом, с помощью схемы на основе метода Рунге–Кутты четвертого порядка (метод (c) из [15], расчет на основе метода третьего порядка (b) также проводился, существенных отличий нет), $t = 5500$, ${{v}_{0}} = - 0.25$. Приведены графики: $r$ – кривая 1, $z{\text{'}}$ – кривая 2, $v$ – кривая 3, график скорости участка газодинамических волн показан в рамке. Решалась система (2.1), (2.2), (2.4), (2.8), изотермическая модель $P = C{{\rho }_{f}}$, ${{P}_{e}} = 0$, $C = 10$, $H = 0.05$, $\tilde {h} = 0.1$. Существенного роста амплитуды возмущений при таком длительном времени расчета обнаружено не было. Наблюдается (слева направо) в основном продольные упругие волны, бегущие влево (простая волна), газодинамические волны, бегущие влево (предположительно простая волна), классический кинк, газодинамические волны, бегущие вправо (солитонная структура), продольные упругие волны, бегущие вправо (простая волна). Переходных элементов, соответствующих в основном поперечным волнам, как и в случае контролируемого давления, нет, поскольку соответствующая дисперсионная ветвь не пересекает начало координат и возможны поперечные колебания трубы как единого целого.

Фиг. 3.

Были проведены исследования на основе предложенной в п. 2.4 методики диссипативных и дисперсионных добавок с коэффициентами, определяемыми шагом сетки.

На фиг. 4 приведены результаты расчета для ${{v}_{0}} = - 1$, $t = 3500$, использован метод Рунге–Кутты второго порядка с добавленной вязкостью во все уравнения ${{c}_{{r6}}} = {{c}_{{z6}}} = {{c}_{{v6}}} = {{c}_{{\rho f6}}} = 0.00001$ (добавка по образцу уравнений (2.6) и (2.7)), $\tilde {h} = 0.1$. Добавка согласуется с рекомендацией п. 3.4 для расчета разрывов в случае аналогов уравнения переноса (${{c}_{6}} = \tilde {C}{{\tilde {h}}^{5}}$, $\tilde {C} = 1$). Расчет без вязкости здесь осуществить не удается, в том числе и с применением метода Рунге–Кутты четвертого порядка. При отсутствии диссипативного члена спустя некоторое время от начала расчета на некотором участке наблюдаются рост и несимметричное укручение волн (до этого момента имеется поточечная сходимость) с последующим их дроблением на более короткие волны. При уменьшении пространственного шага длина этих волн уменьшается. Это можно расценивать как результат действия схемной дисперсии в случае, когда происходит опрокидывание волны. На фиг. 4 наблюдаются (слева направо) в основном продольная структура, бегущая влево, газодинамические волны, бегущие влево (солитонная структура), классический кинк, газодинамические волны, бегущие вправо (солитонная структура), продольные упругие волны, бегущие вправо (простая волна). Между структурами продольных и газодинамическими волн имеется участок с фактически последовательностью из семи уединеннных волн, структуру разрыва не образующий. Участок графика скорости, представляющий наибольший интерес, на фиг. 4 показан в рамке. Структура разрыва продольных упругих волн слева здесь фактически рассчитана как диссипативная, занимающая примерно 6 ячеек, но она излучает волны малой амплитуды. Уменьшение $\tilde {h}$ в два раза при сохранении $\tilde {C}$ показало, что структура продолжает занимать 6 ячеек. Амплитуда и длина излучаемой волны при этом не меняются, эти величины не меняются и при изменении $\tilde {C}$ в диапазоне $4 - 0.25$. Можно предположить наличие в области продольных волн предельного обобщенного решения, получаемого при $\tilde {h} \to 0$.

Фиг. 4.

В случае $C = 0.1$, ${{v}_{0}} = - 0.5$ проводился аналогичный расчет с методом Рунге–Кутты второго порядка с добавленной вязкостью, фиг. 5а, $t = 2000$, а также методом Рунге–Кутты четвертого порядка с добавкой дисперсионного члена для устранения блоуапа продольных волн с $b = 0.01$, согласно рекомендациям п. 3.4 для волнового уравнения ($b = {{c}_{4}} = \tilde {C}{{\tilde {h}}^{2}}$, $\tilde {C} = 1$). На фиг. 5б показан фрагмент зоны продольных волн для сравнения результатов расчетов с вязкой (жирной линией) и дисперсионной добавкой (тонкой линией). В зоне продольных волн имеются участок центрированной простой волны и две структуры разрыва, одна из которых догоняет другую (выявляется при более длительном счете). Пример расчета трубы с жидкостью и с добавленной жесткостью при $b = 0.1$ есть в [6]. Добавление жесткости на изгиб порождает излучение дополнительной короткой волны, приводящей к решению хаотического типа, в случае [6] длина волны удовлетворяет условиям контроля п. 2.3, т.е. рассматривается физическая структура. В данном случае, как и следовало ожидать, ее длина сопоставима с шагом сетки, порядка 12 шагов, полученную структуру можно рассматривать как схемную. Короткая дополнительная волна накладывается на длинную, аналогичную наблюдаемой в расчете с диссипацией, на начальном участке структура упорядоченная. Расчет показал, что отличие от варианта с диссипативной добавкой проявилось только в области структуры продольных волн. Это было использовано для исследования других структур, газодинамической и кинка, для проверки отсутствия воздействия на них добавленных диссипативных членов. Данный расчет с низкой скоростью звука специально проводился для проверки возможности сращивания кинка и газодинамической структуры, сращивания обнаружено не было (см. график для области кинка на фиг. 5в рамке).

Фиг. 5.

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

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

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

Волобуев А.Н. Течение жидкости в трубках с эластичными стенками // Успехи физ. наук. 1995. Т. 165. № 2. С. 177–186.
Epstein M., Johnston J.W. On the exact speed and amplitude of solitary waves in fluid-filled elastic tubes // Proc. Roy. Soc. London. A. 2001. V. 457. I. 2009. P. 1195–1213.
Fu Y.B., Il’ichev A. Solitary waves in fluid-filled elastic tubes: existence, persistence, and the role of axial displacement // IMA J. Appl. Math. 2010. V. 75. № 2. P. 257–268.
Il’ichev A.T., Fu I.B. Stability of aneurism solutions in fluid-filled elastic membrane tube // Acta Mech. Sinica. 2012. V. 28. № 4. P. 1209–1218.
Il’ichev A.T., Fu I.B. Localized standing waves in a hyperelastic membrane tube and their stabilization by a mean flow // Math. Mech. Solids. 2015. V. 20. I. 10. P. 1198–2014.
Бахолдин И.Б. Применение теории обратимых разрывов для исследования уравнений, описывающих волны в трубах с упругими стенками // Прикл. матем. и механ. 2017. Т. 81. Вып. 4.
Fu Y.B., Pearce S.P. Characterization and stability of localized bulging/necking in inflated membrane tubes // IMA J. Appl. Math. 2010. V. 75. № 4. P. 581–602.
Бахолдин И.Б. Численное исследование уединенных волн и структур разрывов в трубах с контролируемым давлением // Ж. вычисл. матем. и матем. физ. 2015. Т. 55. № 11. С. 120–136.
Бахолдин И.Б. Бездиссипативные разрывы в механике сплошной среды. М.: Физматлит, 2004. 318 с.
Бахолдин И.Б. Стационарные и нестационарные структуры разрывов для моделей, описываемых обобщенным уравнением Кортевега–Бюргерса // Прикл. матем. и механ. 2011. Т. 75. Вып. 2. С. 271–302.
Бахолдин И.Б. Теория и классификация обратимых структур разрывов в моделях гидродинамического типа // Прикл. матем. и механ. 2014. Т. 78. Вып. 6. С. 833–852.
Куликовский А.Г., Свешникова Е.И. Нелинейные волны в упругих средах. М.: Моск. лицей, 1998. 412 с.
Куликовский А.Г. Об устойчивости однородных состояний // Прикл. матем. и механ. 1966. Т. 30. Вып. 1. С. 148–153.
Бахолдин И.Б., Егорова Е.Р., Исследование магнитозвуковых уединенных волн для уравнений электронной магнитной гидродинамики // Ж. вычисл. матем. и матем. физ. 2011. Т. 51. № 3. С. 515–528.
Korn G.A., Korn T.M. Mathematical handbook for scientists and engineers. N.Y.: McGraw-Hill Book Companey. 1968; Корн Г., Корн Т. Справочник по математике для научных работников и инженеров. М.: Наука, 1978. 832 с.
Альшина Е.А., Закс Е.М., Калиткин Н.Н. Оптимальные схемы Рунге–Кутты с первого по шестой порядок точности // Ж. вычисл. матем. и матем. физ. 2008. Т. 48. № 3. С. 418–429.
Годунов С.К., Рябенький В.С. Разностные схемы. М.: Наука, 1977. 439 с.
Roach P. J. Computational Fluid Dynamics. Albuquerque, NM: Hermosa, 1976; Роуч П. Дж. Вычислительная гидродинамика. М.: Мир, 1980. 616 с.

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

Инструменты

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

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

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