1. На главную
  2. Электронные версии
  3. Программирование
  4. № 1, 2023

Программирование, 2023, № 1, стр. 32-37

ИССЛЕДОВАНИЯ РАЗНОСТНЫХ СХЕМ ДЛЯ ДВУМЕРНЫХ УРАВНЕНИЙ НАВЬЕ–СТОКСА АЛГОРИТМАМИ КОМПЬЮТЕРНОЙ АЛГЕБРЫ

Ю. А. Блинков a, b, c, *, А. Ю. Ребрина d, **

a Саратовский национальный исследовательский государственный университет имени Н.Г. Чернышевского
410012 Саратов, ул. Астраханская, д. 83, Россия

b Российский университет дружбы народов (РУДН)
117198 Москва, ул. Миклухо-Маклая, д. 6, Россия

c Объединенный институт ядерных исследований
141980 Дубна, ул. Жолио-Кюри, д. 6, Россия

d Саратовский государственный технический университет имени Гагарина Ю.А.
410054 Саратов, ул. Политехническая, д. 77, Россия

* E-mail: blinkovua@info.sgu.ru
** E-mail: anrebrina@yandex.ru

Поступила в редакцию 01.07.2022
После доработки 23.08.2022
Принята к публикации 07.09.2022

Полный текст (PDF)

Аннотация

На основе алгоритма построения базисов Грёбнера рассмотрен класс совместных разностных схем для уравнений Навье–Стокса несжимаемой жидкости в физических переменных и их дифференциальные приближения. Представлены результаты исследования первых дифференциальных приближений этих схем, выполненные авторскими программами, реализованными в системе компьютерной алгебре SymPy. Для рассмотренных разностных схем показана квадратичная зависимость погрешности рассмотренных разностных схем для больших чисел Рейнольдса и обратная пропорциональная для ползущих течений.

1. ВВЕДЕНИЕ

Найти точные решения систем уравнений в частных производных удается лишь в весьма редких случаях, поэтому обычно такие системы в механике и физике приходится решать численно, то есть путем их замены дискретными аналогами. Наиболее известными методами дискретизации и численного решения являются методы конечных элементов, конечных объемов и конечных разностей. Последний метод исторически является первым [1] и основан на замене дифференциальных уравнений разностными, определенными на выбранной сетке решений. Чтобы построить численное решение, разработанное конечно-разностное приближение уравнений в частных производных дополняется соответствующей дискретизацией начальных или/и граничных условий. При этом близость и сходство численного и точного решений системы уравнений в частных производных в решающей степени зависит от качества разностной схемы.

В работах [2, 3] разработан и реализован средствами компьютерной алгебры Maple алгоритм построения консервативных разностных схем. После задания искомой системы в виде условий связи разностных функций и их разностных производных в виде интегральных соотношений и выбора соответствующего упорядочения путем вычисления разностного базиса Грёбнера, который адаптирует к разностному случаю известный алгоритм для исследования полиномиальных систем [4]. В результате применения алгоритма построения базиса Грёбнера разностного идеала для дифференциальных уравнений с производными выше первого порядка строится разностная схема в виде соотношения только на сами функции. При это проходит проверка на совместность полученной разностной схемы и определяется ее размерность на пространство решений.

Одним из алгоритмических подходов для исследования разностных схем является построение первого дифференциального приближения [5]. На его основе и с использованием базисов Грёбнера построены дифференциальные приближения для уравнений эволюционного типа [6] и для систем уравнений [7]. В работе [8] рассмотрены вопросы вычисления первого дифференциального приближения в системах компьютерной алгебры. Если разностная схема не совместна или имеет меньшее пространство решений по сравнению с исходной системой дифференциальных уравнений, то при построении первого дифференциального приближения будут получены лишние уравнения и вычисления можно остановить. Эта проверка значительно проще, чем проверка совместности в разностном случае как по требуемому времени, так и по памяти и легко может быть реализована обычными средствами систем компьютерной алгебры.

В данной работе предлагается в рамках универсального алгоритма для коммутативной алгебры средствами компьютерной алгебры проверить совместность исходной системы дифференциальных уравнений и аппроксимирующей ее разностной схемой, прямо и в рамках первого дифференциального приближения. Хотя алгоритм построения базисов Грёбнера и встроен в большинство систем компьютерной алгебры, он как правило имеет ограничения на работу с дифференциальными уравнениями, разностными схемами и формально бесконечными рядами Тейлора при построении первого дифференциального приближения [9–11].

В качестве модельной системы для реализованных в системе компьютерной алгебры SymPy будет использована система двухмерных уравнений Навье–Стокса несжимаемой жидкости. В работах [12–17] были построены и численно исследованы разностные схемы на некоторых точных решениях.

2. НЕПРЕРЫВНЫЙ СЛУЧАЙ

Запишем систему двухмерных уравнений Навье–Стокса несжимаемой жидкости в декартовой системе координат

(2.1)
$\begin{gathered} {{F}^{1}}:\;\;{{u}_{x}} + {{v}_{y}} = 0, \\ {{F}^{2}}:\;\;{{u}_{t}} + u{{u}_{x}} + v{{u}_{y}} + {{p}_{x}} - \frac{1}{{{\text{Re}}}}\Delta u = 0, \\ {{F}^{3}}:\;\;{{v}_{t}} + u{{v}_{x}} + v{{v}_{y}} + {{p}_{y}} - \frac{1}{{{\text{Re}}}}\Delta v = 0. \\ \end{gathered} $

Пусть $p \succ u \succ v$ и $t \succ x \succ y$. Введем допустимое упорядочение по POT (позиция старше терма) таким способом, что любая производная по t старше любой другой производной.

Тогда лидирующими членами для системы (2.1) будут соответственно ${{u}_{x}},{{u}_{t}},{{{v}}_{t}}$ и можно будет построить следующий $S$-полином:

(2.2)
$\begin{gathered} (F_{x}^{2} - F_{t}^{1}) + F_{y}^{3} + \frac{1}{{{\text{Re}}}}\Delta {{F}^{1}} = \\ \, = {{\left( {u{{u}_{x}} + v{{u}_{y}}} \right)}_{x}} + {{\left( {u{{v}_{x}} + v{{v}_{y}}} \right)}_{y}} + \Delta p. \\ \end{gathered} $

В результате получим инволютивную систему:

${{F}^{1}}{\text{ : }}{{u}_{x}} + {{v}_{y}} = 0,$
(2.3)
$\begin{gathered} {{F}^{2}}:{{u}_{t}} + u{{u}_{x}} + v{{u}_{y}} + {{p}_{x}} - \frac{1}{{{\text{Re}}}}\Delta u = 0, \\ {{F}^{3}}:{{v}_{t}} + u{{v}_{x}} + v{{v}_{y}} + {{p}_{y}} - \frac{1}{{{\text{Re}}}}\Delta v = 0, \\ \end{gathered} $
${{F}^{4}}:{{\left( {u{{u}_{x}} + v{{u}_{y}}} \right)}_{x}} + {{\left( {u{{v}_{x}} + v{{v}_{y}}} \right)}_{y}} + \Delta p = 0,$
где лидирующими членами будут соответственно ${{u}_{x}},{{u}_{t}},{{{v}}_{t}},{{p}_{{xx}}}$. Эту систему можно переписать в виде законов сохранения:

${{F}^{1}}:{\text{div}}(u,v) = 0,$
(2.4)
$\begin{gathered} {{F}^{2}}:{{u}_{t}} + {\text{div}}\left( {{{u}^{2}} + p - \frac{1}{{{\text{Re}}}}{{u}_{x}},vu - \frac{1}{{{\text{Re}}}}{{u}_{y}}} \right) = 0, \\ {{F}^{3}}:{{v}_{t}} + {\text{div}}\left( {uv - \frac{1}{{{\text{Re}}}}{{v}_{x}},{{v}^{2}} + p - \frac{1}{{{\text{Re}}}}{{v}_{y}}} \right) = 0, \\ \end{gathered} $
${{F}^{4}}:{\text{div}}\left( {u{{u}_{x}} + v{{u}_{y}} + {{p}_{x}},{\kern 1pt} u{{v}_{x}} + v{{v}_{y}} + {{p}_{y}}} \right) = 0.$

Понятие инволютивности было введено около ста лет тому назад Картаном (Cartan) [18] при исследовании уравнений типа Пфаффа в полных дифференциалах. Инволютивная система содержит в себе все условия интегрируемости и дифференциальные продолжения системы не дают новых условий совместности. Пополнение системы ее условиями интегрируемости называется замыканием. Подход Картана был обобщен Келером (Kähler) [19] к произвольным системам внешних дифференциальных уравнений. Рикье (Riquier) [20], для исследования решений дифференциальных уравнений в виде формальных рядов, предложил полное упорядочение для частных производных. Используя упорядочение он выделил часть производных, называемых главными, относительно которых можно разрешить систему ДУЧП. Оставшиеся производные, называемые параметрическими, задают произвол в решении и влияют на постановку начальных условий. В результате Рикье была построена теория содержащая теорему Коши–Ковалевской как частный случай.

3. КОНЕЧНО-РАЗНОСТНАЯ АППРОКСИМАЦИЯ

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

1. Для дискретной производной D не выполняется правило Лейбница: ${\text{D}}(ab) \ne {\text{D}}(a)b + a{\text{D}}(b)$.

2. Если в непрерывном случае допустимое упорядочение было связанно с производными, то в дискретном случае оно привязано к дискретным координатам n, j, k (им соответствуют t, x, y) самой функции. Другими словами, в непрерывном случае упорядочение допустимо относительно производных, а в дискретном допустимо относительно сдвигов.

3. Сдвиги могут быть вперед и могут быть назад, а базисы Грёбнера могут работать только с положительными, поэтому желательно сдвинуть первоначальные разностные соотношения вперед для упрощения итогового базиса Грёбнера.

В статье [12], чтобы обойти пункт 2 для получения конечно разностной аппроксимации системы (2.1) использовали дискретные операторы, которые обладали коммутативностью и позволяли применить теорию базисов Грёбнера, а в статье [13] для вычисления разностных базисов Грёбнера использовали их специализированную реализацию средствами компьютерной алгебры.

3.1. MAC схема

Метод маркеров и частиц (MAC) [21] представляет очень известный вариант явной разностной схемы на разнесенной сетке с использованием осреднений для компонент скорости:

$\begin{gathered} \frac{{u_{{j + 1/2,k}}^{n} - u_{{j - 1/2,k}}^{n}}}{h} + \frac{{v_{{j,k + 1/2}}^{n} - v_{{j,k - 1/2}}^{n}}}{h} = 0, \\ u_{{j + 1/2,k}}^{n} = F_{{j + 1/2,k}}^{n} - \frac{\tau }{h}(p_{{j + 1,k}}^{{n + 1}} - p_{{j,k}}^{{n + 1}}), \\ \end{gathered} $
(3.1)
$\begin{gathered} v_{{j,k + 1/2}}^{n} = G_{{j,k + 1/2}}^{n} - \frac{\tau }{h}(p_{{j,k + 1}}^{{n + 1}} - p_{{j,k}}^{{n + 1}}), \\ \frac{{p_{{j + 1,k}}^{{n + 1}} - 2p_{{j,k}}^{{n + 1}} + p_{{j - 1,k}}^{{n + 1}}}}{{{{h}^{2}}}} + \frac{{p_{{j,k + 1}}^{{n + 1}} - 2p_{{j,k}}^{{n + 1}} + p_{{j,k - 1}}^{{n + 1}}}}{{{{h}^{2}}}} = \\ \, = \frac{1}{\tau }\left( {\frac{{F_{{j + 1/2,k}}^{n} - F_{{j - 1/2,k}}^{n}}}{h} + \frac{{G_{{j,k + 1/2}}^{n} - G_{{j,k - 1/2}}^{n}}}{h}} \right), \\ \end{gathered} $
$\begin{gathered} F_{{j + 1/2,k}}^{n} = \ldots , \\ G_{{j,k + 1/2}}^{n} = \ldots . \\ \end{gathered} $

Уравнение Пуассона решается на каждом шаге по времени и после того, как определено значение $p_{{j,k}}^{{n + 1}}$, подстановка этого значения в остальные уравнения системы (3.1) позволяет определить компоненты скорости на временном слое n + 1.

Также стоит отметить, что фактически в статье [21] осуществлена проверка совместности системы (3.1) путем построения s-полинома. Как дословно пишут авторы, что дискретная производная по времени от уравнения неразрывности равна нулю. Удивительно, что одновременно, в том же году, была построена теория базисов Грёбнера [22] Бруно Бухбергером.

3.2. Обычная разностная схема

Поскольку при построении базиса Грёбнера нельзя работать с отрицательными сдвигами, рассмотрим разностные соотношения в точке $j + l$, $k + l$, где $l > 0$ достаточно большое число, чтобы избежать отрицательных сдвигов. Запишем соотношения для уравнения неразрывности, например, в точке $j + l + m,$ $k + l$ и $m > 0$, чтобы получить более компактный базис Грёбнера. В результате при любом выборе разностных производных ${{D}_{1}},\;{{D}_{2}},\;{{D}_{3}},\;\Delta $ и осреднения по времени ${{I}_{t}}$ (для построения явной, полунеявной или чисто неявной схемы) следующая система будет представлять собой разностный базис Грёбнера

$\begin{gathered} {{F}^{1}}:{{D}_{{1x}}}u_{{j + l + m,{\kern 1pt} k + l}}^{n} + {{D}_{{1y}}}v_{{j + l + m,{\kern 1pt} k + l}}^{n} = 0, \\ {{F}^{2}}:\frac{{u_{{j + l,{\kern 1pt} k + l}}^{{n + 1}} - u_{{j + l,{\kern 1pt} k + l}}^{n}}}{\tau } + {{I}_{t}}\left( {u_{{j + l,{\kern 1pt} k + l}}^{n}{{D}_{{2x}}}u_{{j + l,{\kern 1pt} k + l}}^{n}{{ + }_{{_{{_{{_{{_{{_{{}}}}}}}}}}}}}} \right. \\ \, + v_{{j + l,{\kern 1pt} k + l}}^{n}{{D}_{{2y}}}u_{{j + l,{\kern 1pt} k + l}}^{n} + {{D}_{{3x}}}p_{{j + l,{\kern 1pt} k + l}}^{n} - \left. {\frac{1}{{{\text{Re}}}}\Delta u_{{j + l,{\kern 1pt} k + l}}^{n}} \right) = 0, \\ \end{gathered} $
(3.2)
$\begin{gathered} {{F}^{3}}:\frac{{v_{{j + l,{\kern 1pt} k + l}}^{{n + 1}} - v_{{j + l,{\kern 1pt} k + l}}^{n}}}{\tau } + {{I}_{t}}\left( {u_{{j + l,{\kern 1pt} k + l}}^{n}{{D}_{{2x}}}v_{{j + l,{\kern 1pt} k + l}}^{n}{{ + }_{{_{{_{{_{{_{{_{{}}}}}}}}}}}}}} \right. \\ \, + v_{{j + l,{\kern 1pt} k + l}}^{n}{{D}_{{2y}}}v_{{j + l,{\kern 1pt} k + l}}^{n} + \\ \, + {{D}_{{3y}}}p_{{j + l,{\kern 1pt} k + l}}^{n} - \left. {\frac{1}{{{\text{Re}}}}\Delta v_{{j + l,{\kern 1pt} k + l}}^{n}} \right) = 0, \\ \end{gathered} $
$\begin{gathered} {{F}^{4}}:{{D}_{{1x}}}(u_{{j + l,{\kern 1pt} k + l}}^{n}{{D}_{{2x}}}u_{{j + l,{\kern 1pt} k + l}}^{n} + v_{{j + l,{\kern 1pt} k + l}}^{n}{{D}_{{2y}}}u_{{j + l,{\kern 1pt} k + l}}^{n}) + \\ \, + {{D}_{{1y}}}(u_{{j + l,{\kern 1pt} k + l}}^{n}{{D}_{{2x}}}v_{{j + l,{\kern 1pt} k + l}}^{n} + v_{{j + l,{\kern 1pt} k + l}}^{n}{{D}_{{2y}}}v_{{j + l,{\kern 1pt} k + l}}^{n}) + \\ \, + ({{D}_{{1x}}}{{D}_{{3x}}}p_{{j + l,{\kern 1pt} k + l}}^{n} + {{D}_{{1y}}}{{D}_{{3y}}}p_{{j + l,{\kern 1pt} k + l}}^{n}) = 0. \\ \end{gathered} $

Для проверки совместности системы (3.2) построим следующий $S$-полином, аналогичный (2.2) для непрерывного случая:

(3.3)
$\begin{gathered} {{D}_{{1x}}}{{F}^{2}} + {{D}_{{1y}}}{{F}^{3}} - \frac{{({{D}_{{1x}}}u_{{j + l,{\kern 1pt} k + l}}^{{n + 1}} + {{D}_{{1y}}}v_{{j + l,{\kern 1pt} k + l}}^{{n + 1}})}}{\tau } - \\ \, - \frac{{({{D}_{{1x}}}u_{{j + l,{\kern 1pt} k + l}}^{n}\, + \,{{D}_{{1y}}}v_{{j + l,{\kern 1pt} k + l}}^{n})}}{\tau }\, + \,{{I}_{t}}\left( {\frac{1}{{{\text{Re}}}}\Delta {{F}^{1}}} \right)\, - \,{{I}_{t}}({{F}^{4}}). \\ \end{gathered} $

Система (3.2) совместна, если этот многочлен тождественно равен нулю. В отличии от схемы (3.1) , схема (3.2) используя уравнения F2, F3 строит значение для скоростей на следующем временном слое. Уравнение неразрывности F1 выполнено автоматически. Уравнения F4 служит только для определения давления для следующего временного слоя.

3.3. Разностная схема в дивергентном виде

Аналогично, при любом выборе разностных производных ${{D}_{1}},\;{{D}_{2}},\;{{D}_{3}},\;\Delta $ и осреднения по времени ${{I}_{t}}$ следующая система будет представлять собой разностный базис Грёбнера

${{F}^{1}}:{{D}_{{1x}}}u_{{j + l + m,{\kern 1pt} k + l}}^{n} + {{D}_{{1y}}}v_{{j + l + m,{\kern 1pt} k + l}}^{n} = 0,$
$\begin{gathered} {{F}^{2}}:\frac{{u_{{j + l,{\kern 1pt} k + l}}^{{n + 1}} - u_{{j + l,{\kern 1pt} k + l}}^{n}}}{\tau } + {{I}_{t}}\left( {{{D}_{{2x}}}(u{{{_{{j + l,{\kern 1pt} k + l}}^{n}}}^{2}}){{ + }_{{_{{_{{_{{_{{_{{_{{}}}}}}}}}}}}}}}} \right. \\ \, + {{D}_{{2y}}}(v_{{j + l,{\kern 1pt} k + l}}^{n}u_{{j + l,{\kern 1pt} k + l}}^{n}) + {{D}_{{3x}}}p_{{j + l,{\kern 1pt} k + l}}^{n} - \\ \left. {\, - \frac{1}{{{\text{Re}}}}\Delta u_{{j + l,{\kern 1pt} k + l}}^{n}} \right) = 0, \\ \end{gathered} $
(3.4)
$\begin{gathered} {{F}^{3}}:\frac{{v_{{j + l,{\kern 1pt} k + l}}^{{n + 1}} - v_{{j + l,{\kern 1pt} k + l}}^{n}}}{\tau } + \\ \, + {{I}_{t}}\left( {{{D}_{{2x}}}(u_{{j + l,{\kern 1pt} k + l}}^{n}v_{{j + l,{\kern 1pt} k + l}}^{n}) + {{D}_{{2y}}}(v{{{_{{j + l,{\kern 1pt} k + l}}^{n}}}^{2}}){{ + }_{{_{{_{{_{{_{{_{{}}}}}}}}}}}}}} \right. \\ \left. {\, + {{D}_{{3y}}}p_{{j + l,{\kern 1pt} k + l}}^{n} - \frac{1}{{{\text{Re}}}}\Delta v_{{j + l,{\kern 1pt} k + l}}^{n}} \right) = 0, \\ \end{gathered} $
$\begin{gathered} {{F}^{4}}:{{D}_{{1x}}}({{D}_{{2x}}}(u{{_{{j + l,{\kern 1pt} k + l}}^{n}}^{2}}) + {{D}_{{2y}}}(v_{{j + l,{\kern 1pt} k + l}}^{n}u_{{j + l,{\kern 1pt} k + l}}^{n})) + \\ \, + {{D}_{{1y}}}({{D}_{{2x}}}(u_{{j + l,{\kern 1pt} k + l}}^{n}v_{{j + l,{\kern 1pt} k + l}}^{n}) + {{D}_{{2y}}}(v{{_{{j + l,{\kern 1pt} k + l}}^{n}}^{2}})) + \\ \, + ({{D}_{{1x}}}{{D}_{{3x}}}p_{{j + l,{\kern 1pt} k + l}}^{n} + {{D}_{{1y}}}{{D}_{{3y}}}p_{{j + l,{\kern 1pt} k + l}}^{n}) = 0. \\ \end{gathered} $

Для проверки совместности системы (3.4) можно построить $S$-полином, в точности совпадающий с (3.3).

Аналогично схеме (3.2) для схемы (3.4) используя уравнения F2, F3 строится значение для скоростей на следующем временном слое. Уравнение неразрывности ${{F}^{1}}$ выполнено автоматически, а уравнения F4 служит только для определения давления для слоя n + 1.

4. ДИФФЕРЕНЦИАЛЬНЫЕ ПРИБЛИЖЕНИЯ

Построение первого дифференциального приближения было сделано с оригинальными программами авторов статьи в системе алгебры SymPy и может быть скачана https://github. com/blinkovua/sharing-blinkov/blob/ master/NavierStokes2D.ipynb.

Рассмотрим сначала аппроксимацию уравнений (3.1) до первых ненулевых членов:

${{F}^{1}}:{{u}_{x}} + {{v}_{y}} + \tau ({{u}_{{tx}}}{\text{/}}2 + \ldots ) + {{h}^{2}}({{u}_{{xxx}}}{\text{/}}24 + \ldots ),$
$\begin{gathered} {{F}^{2}}:{{p}_{x}} + 2u{{u}_{x}} + u{{v}_{y}} + {{u}_{t}} + {{u}_{y}}v + ( - {{u}_{{xx}}} - {{u}_{{yy}}}){\text{/Re}} + \\ \, + \tau ({{p}_{{tx}}} + \ldots )\, + \,{{h}^{2}}(({{p}_{{xxx}}}{\text{/}}6 + \ldots )\, - \,(5{{u}_{{xxxx}}}{\text{/}}24 + \ldots ){\text{/Re}}), \\ \end{gathered} $
$\begin{gathered} {{F}^{3}}:{{p}_{y}} + u{{v}_{x}} + {{u}_{x}}v + 2v{{v}_{y}} + {{v}_{t}} + ( - {{v}_{{xx}}} - {{v}_{{yy}}}){\text{/Re}} + \\ + \tau ({{p}_{{ty}}} + \ldots ) + {{h}^{2}}(({{p}_{{yyy}}}{\text{/}}6 + \ldots ) - ({{v}_{{xxxx}}}{\text{/}}12 + \ldots ){\text{/Re}})), \\ \end{gathered} $
$\begin{gathered} {{F}^{4}}:{{p}_{{xx}}} + {{p}_{{yy}}} + 2u{{u}_{{xx}}} + 2u{{v}_{{xy}}} + 2{{u}_{{xy}}}v + 2u_{x}^{2} + \\ \, + 2{{u}_{x}}{{v}_{y}} + 2{{u}_{y}}{{v}_{x}} + 2v{{v}_{{yy}}} + 2v_{y}^{2} + \\ \, + ( - {{u}_{{xxx}}} - {{u}_{{xyy}}} - {{v}_{{xxy}}} - {{v}_{{yyy}}}){\text{/Re}} + \\ \, + \tau ({{p}_{{txx}}} + \ldots ) + {{h}^{2}}({{p}_{{xxxx}}}{\text{/}}12 + \ldots ). \\ \end{gathered} $

Некоторая несимметрия в разложениях для F2, F3 связана с порядком переменных. Преобразовав эти уравнения с учетом исходных уравнений (2.1), получим первое дифференциальное приближение для (2.3):

${{F}^{1}}:\;\;{{u}_{x}} + {{v}_{y}} + {{h}^{2}}({\text{Re}}\left( { - {{p}_{{yy}}}{\text{/}}24 + \ldots } \right)) + {{v}_{{yyy}}}{\text{/}}12),$
$\begin{gathered} {{F}^{2}}:\;\;{{p}_{x}} + 2u{{u}_{x}} + u{{v}_{y}} + {{u}_{t}} + {{u}_{y}}v + ( - {{u}_{{xx}}} - {{u}_{{yy}}}){\text{/Re}} + \\ \, + \tau ({{p}_{{tx}}} + \ldots ) + {{h}^{2}}({\text{R}}{{{\text{e}}}^{2}}{{v}^{2}}({{p}_{x}}{\text{/}}12 + \ldots ) + \\ \, + {\text{Re}}( - {{p}_{{tx}}}{\text{/}}12 + \ldots ) + ({{p}_{{xyy}}}{\text{/}}12 + \ldots ) - {{v}_{{xyyy}}}{\text{/}}(6{\text{Re}})), \\ \end{gathered} $
$\begin{gathered} {{F}^{3}}:\;\;{{p}_{y}} + u{{v}_{x}} + {{u}_{x}}v + 2v{{v}_{y}} + {{v}_{t}} + ( - {{v}_{{xx}}} - {{v}_{{yy}}}){\text{/Re}} + \\ \, + \tau (({{p}_{{ty}}} + \ldots )) + {{h}^{2}}({\text{R}}{{{\text{e}}}^{2}}{{u}^{2}}({{p}_{y}}{\text{/}}12 + \ldots ) + \\ \, + {\text{Re}}({{p}_{{ty}}}{\text{/}}12 + \ldots ) + (5{{p}_{{yyy}}}{\text{/}}24 + \ldots ) - {{v}_{{yyyy}}}{\text{/}}(6{\text{Re}})) \\ \end{gathered} $
$\begin{gathered} {{F}^{4}}:\;\;{{p}_{{xx}}} + {{p}_{{yy}}} + 2u{{u}_{{xx}}} + 2u{{v}_{{xy}}} + 2{{u}_{{xy}}}v + 2u_{x}^{2} + \\ \, + 2{{u}_{x}}{{v}_{y}} + 2{{u}_{y}}{{v}_{x}} + 2v{{v}_{{yy}}} + 2v_{y}^{2} + \\ \, + ( - {{u}_{{xxx}}} - {{u}_{{xyy}}} - {{v}_{{xxy}}} - {{v}_{{yyy}}}){\text{/Re}} + \tau ( - 2{{u}_{{ty}}}{{v}_{x}} + \ldots ) + \\ \, + {{h}^{2}}({\text{R}}{{{\text{e}}}^{2}}(19{{p}_{x}}v{{v}_{x}}{\text{/}}24 + \ldots ) - \\ \, - {\text{Re}}( - 7{{p}_{{xyy}}}u{\text{/}}24 + \ldots )({{p}_{{yyyy}}}{\text{/}}6 \ldots )) \\ \end{gathered} $

Можно сделать следующие выводы.

• Для ${{F}^{1}}$ нет зависимости от $\tau $, а вот зависимость от h2 имеет зависимость ${{h}^{2}}({\text{Re}}( \ldots ) + ( \ldots ))$. Это означает, что для достижения той же точности для ${{F}^{1}}$ при увеличении числа ${\text{Re}} \gg 1$ необходимо в квадратичный корень раз уменьшить шаг h.

• Многочлены ${{F}^{2}},{{F}^{3}},{{F}^{4}}$ зависят от $\tau $, но не зависит от ${\text{Re}}$. Для ${{F}^{2}},{{F}^{3}}$ зависимость от h2 имеет вид ${{h}^{2}}({\text{R}}{{{\text{e}}}^{2}}( \ldots ) + {\text{Re}}( \ldots ) + ( \ldots )) + ( \ldots )){\text{/Re}}$. Для ${{F}^{4}}$ зависимость от h2 имеет зависимость ${{h}^{2}}({\text{R}}{{{\text{e}}}^{2}}( \ldots )$ + + Re(...) + (...)). Поэтому при ${\text{Re}} \gg 1$ для достижения той же точности для ${{F}^{2}},{{F}^{3}},{{F}^{4}}$ при увеличении числа Re необходимо линейно уменьшить шаг $h$. В случае же ползущих течений, когда ${\text{Re}} \ll 1$, необходимо для достижения той же точности для ${{F}^{4}}$ в квадратичный корень раз уменьшить шаг $h$.

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

(4.1)
$\begin{gathered} {{D}_{x}}g = \frac{{g_{{j + 1,k}}^{n} - g_{{j - 1,k}}^{n}}}{{2h}},\quad {{D}_{y}}g = \frac{{g_{{j,k + 1}}^{n} - g_{{j,k - 1}}^{n}}}{{2h}}, \\ {{\Delta }_{1}} = \frac{{g_{{j + 1,k + 1}}^{n} + g_{{j + 1,k}}^{n} + g_{{j,k - 1}}^{n} + g_{{j - 1,k - 1}}^{n} - 4g_{{j,k}}^{n}}}{{{{h}^{2}}}} \\ {{\Delta }_{2}} = \frac{{g_{{j + 2,k + 2}}^{n} + g_{{j + 2,k}}^{n} + g_{{j,k - 2}}^{n} + g_{{j - 2,k - 2}}^{n} - 4g_{{j,k}}^{n}}}{{{{h}^{2}}}} \\ \end{gathered} $

Последние соотношение в (4.1) на ${{\Delta }_{2}}$ использовано, поскольку шаблон для ${{F}^{4}}$ получается минимум размером 3 × 3.

Для получения явных, полунеявных и неявных схем можно использовать следующие операторы

(4.2)
$I{{t}_{0}}g = g_{{j,k}}^{n},\quad I{{t}_{1}}g = \frac{{g_{{j,k}}^{{n + 1}} + g_{{j,k}}^{n}}}{2},\quad I{{t}_{2}}g = g_{{j,k}}^{{n + 1}}.$

Этим путем из (3.2), (3.4) получены 12 различных разностных схем, для них выполнена проверка равенства нулю S-полинома. Выводы, сделанные из первого дифференциального приближения для (3.1), справедливы и по схемам (3.2), (3.4) , за одним исключением: для всех 12 схем для ${{F}^{4}}$ не зависит от τ. Вероятно, это связано с тем, что для этих схем в отличии от (3.1) для следующего временного слоя сначала вычисляются компоненты скорости $u,\;v$, а уже потом определяется давление p.

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

5. ОЦЕНКА КАЧЕСТВА РАЗНОСТНЫХ СХЕМ НА ТОЧНЫХ РЕШЕНИЯХ

Первое дифференциальное приближение при подстановке в него точного решения позволяет оценить саму схему без ее программирования и проведение вычисленных экспериментов для ее проверки. В двумерном случае известно несколько точных решений системы Навье–Стокса (2.1), описанные в [23] и [24]. Приведем здесь первое из них:

(5.1)
$\begin{gathered} u = - {{e}^{{ - 2t/{\text{Re}}}}}\cos (x)\sin (y), \\ v = {{e}^{{ - 2t/{\text{Re}}}}}\sin (x)\cos (y), \\ p = - {{e}^{{ - 4t/{\text{Re}}}}}(\cos (2x) + \cos (2y)){\text{/}}4. \\ \end{gathered} $

Подстановка этого решения в первое дифференциальное приближение схемы MAC (3.1), приводит к следующему. Первое дифференциальное приближение аппроксимации уравнения неразрывности ${{F}^{1}}$ обращается тождественно в нуль. Первые дифференциальные приближения ${{F}^{2}},{{F}^{3}}$ имеет вид

$\tau ( \ldots ){\text{/R}}{{{\text{e}}}^{2}} + {{h}^{2}}(( \ldots ) + ( \ldots ){\text{/Re}}).$

Подстановка решения (5.1) в первое дифференциальное приближение полученных выше 12 разностных схем дает различные выражения для ${{F}^{2}},{{F}^{3}}$, вид которых зависит от выбора осреднения по времени (4.2). Для схемы (3.2) :

(5.2)
$\begin{gathered} I{{t}_{0}}:\tau ( \ldots ){\text{/R}}{{{\text{e}}}^{2}} + {{h}^{2}}\left( {( \ldots ) + ( \ldots ){\text{/Re}}} \right), \\ I{{t}_{1}}:{{h}^{2}}\left( {( \ldots ) + ( \ldots ){\text{/Re}}} \right), \\ I{{t}_{2}}:\tau ( \ldots ){\text{/R}}{{{\text{e}}}^{2}} + {{h}^{2}}(( \ldots ) + ( \ldots ){\text{/Re}} + \\ \, + ( \ldots ){\text{/R}}{{{\text{e}}}^{2}} + ( \ldots ){\text{/R}}{{{\text{e}}}^{3}}), \\ \end{gathered} $
а для схемы (3.4) :

(5.3)
$\begin{gathered} I{{t}_{0}}:\tau ( \ldots ){\text{/R}}{{{\text{e}}}^{2}} + {{h}^{2}}\left( {( \ldots ){\text{/Re}}} \right), \\ I{{t}_{1}}:{{h}^{2}}\left( \ldots \right){\text{/Re}}, \\ I{{t}_{2}}:\tau ( \ldots ){\text{/R}}{{{\text{e}}}^{2}} + {{h}^{2}}(( \ldots ){\text{/Re}} + ( \ldots ){\text{/R}}{{{\text{e}}}^{3}}). \\ \end{gathered} $

Таким образом, лучшее поведение демонстрируют осреднение для полунеявной схемы.

Особо отметим поведение первого дифференциального приближения для ${{F}^{4}}$ на решении (5.1). Для схемы MAC (3.1):

$\begin{gathered} - \tau (4\left( {\cos \left( {2x} \right) + \cos \left( {2y} \right)} \right){{e}^{{ - \frac{{4t}}{{{\text{Re}}}}}}}){\text{/Re}} + \\ \, + {{h}^{2}}((12\mathop {\sin }\nolimits^2 \left( x \right){{\sin }^{2}}\left( y \right) - 8\mathop {\sin }\nolimits^2 \left( x \right) - \\ \, - 8{{\sin }^{2}}\left( y \right) + 5){\text{/}}4 - \\ \, - 16(\cos \left( {2x} \right) + \cos \left( {2y} \right)){{e}^{{ - \frac{{4t}}{{{\text{Re}}}}}}}{\text{/R}}{{{\text{e}}}^{2}}). \\ \end{gathered} $

Для схемы (3.2) при ${{I}_{t}}$, равным $I{{t}_{0}}$ или $I{{t}_{1}}$, первое дифференциальное приближение для F4 равно:

${{h}^{2}}(({{\sin }^{2}}\left( x \right) + {{\sin }^{2}}\left( y \right) - 1){{e}^{{ - \frac{{4t}}{{{\text{Re}}}}}}}).$

Для схемы (3.2) с вариантами с It2 и для всех вариантов (3.4) первое дифференциальное приближение для F4 равно нулю.

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

Применение алгоритма базисов Грёбнера позволяет проверить на совместность исходную систему, равно как и аппроксимирующую ее разностную схему в первом дифференциальном приближении, что особенно важно, если задача является вычислительно трудной, как например, схема MAC (3.1). Первое дифференциальное приближение позволяет получить важную информацию о разностной схеме. В частности, для рассмотренных в статье разностных схем, аппроксимирующих плоские уравнения Навье–Стокса, это позволило определить их погрешность в зависимости от числа Рейнольдса Re.

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

  1. Самарский А.А. Теория разностных схем. 3-е изд., испр. – М. Наука. 1989. 616 с.

  2. Блинков Ю.А., Мозжилкин В.В. Генерация разностных схем для уравнения Бюргерса построением базисов Грёбнера // Программирование. 2006. № 2. С. 71–74.

  3. Gerdt V.P., Blinkov Yu.A., Mozzhilkin V.V. Gröbner Bases and Generation of Difference Schemes for Partial Differential Equations // SIGMA. 2006. № 2(051). 26 p.

  4. Buchberger B. Gröbner bases: an Buchberger algorithmic method in polynomial ideal theory // Recent Trends in Multidimensional System Theory / Ed. by N.K. Bose. V. 6. Reidel, Dordrecht, 1985. P. 184–232.

  5. Шокин Ю.И., Яненко Н.Н. Метод дифференциального приближения. Применение к газовой динамике. Новосибирск: Наука: Сиб. отд-ние, 1985. 364 с.

  6. Блинков Ю.А., Гердт В.П., Маринов К.Б. Дискретизация квазилинейных эволюционных уравнений методами компьютерной алгебры // Программирование. 2017. № 2. С. 28–34.

  7. Zhang X., Gerdt V.P., Blinkov Y.A. Algebraic Construction of a Strongly Consistent, Permutationally Symmetric and Conservative Difference Scheme for 3D Steady Stokes Flow // Symmetry. 2019. № 11(269). 15 p.

  8. Блинков А.Ю., Малых М.Д., Севастьянов Л.А. О дифференциальных приближениях разностных схем // Изв. Сарат. ун-та. Нов. cер. Сер. Математика. Механика. Информатика. 2021. № С. 472–488.

  9. Gerdt V.P., Robertz D. Consistency of Finite Difference Approximations for Linear PDE Systems and its Algorithmic Verification / In: S. Watt (ed.). Proceedings of ISSAC 2010 P. 53–59.

  10. Gerdt V.P. Consistency Analysis of Finite Difference Approximations to PDE Systems / Mathematical Modelling in Computational Physics 2011. LNCS. V. 7125. 2012. P. 28–42.

  11. Scala R.L. Gröbner bases and gradings for partial difference ideals // Math. Comput. 2015. № 84. P. 959–985.

  12. Gerdt V.P., Blinkov Yu.A. Involution and Difference Schemes for the Navier–Stokes Equations / In: Gerdt V.P., Mayr E.W., Vorozhtsov E.V. (eds.) CASC 2009. LNCS. 2009. V. 5743. P. 94–105.

  13. Amodio P., Blinkov Yu.A., Gerdt V.P., La Scala R. On Consistency of Finite Difference Approximations to the Navier-Stokes Equations / In: Gerdt V.P., Koepf W., Mayr E.W., Vorozhtsov E.V. (eds.) CASC 2013. LNCS. 2013. V. 8136. P. 46–60.

  14. Amodio P., Blinkov Yu.A., Gerdt V.P., Scala R.La. Algebraic construction and numerical behavior of a new s-consistent difference scheme for the 2D Navier–Stokes equations // Appl. Math. and Comput. 2017. № 314. P. 408–421.

  15. Blinkov Yu.A., Gerdt V.P., Lyakhov D.A., Michels D.L. A Strongly Consistent Finite Difference Scheme for Steady Stokes Flow and its Modified Equations / In: Gerdt V.P., Koepf W., Mayr E.W., Vorozhtsov E.V. (eds.) CASC 2018. LNCS. 2018. V. 11077. P. 67–81.

  16. Michels D.L., Gerdt V.P., Blinkov Y.A., Lyakhov D.A. On the Consistency Analysis of Finite Difference Approximations // Journal of Mathematical Sciences (United States). 2019. № 5. P. 665–677.

  17. Gerdt V.P., Robertz D., Blinkov Yu.A. Strong Consistency and Thomas Decomposition of Finite Difference Approximations to Systems of Partial Differential Equations / eprint arXiv:2009.01731. 2019. 48 p.

  18. Cartan E. Sur certaines expressions différentielles et le problème de Pfaff // Annales scientifiques de l’École Normale Supérieure Sér. 3. 1899. № 16. P. 239–332.

  19. Kähler E. Einführung in die Theorie der Systeme von Differentialgleichungen. Hamburger mathematische Einzelschriften 16. Leipzig: B.G. Teubner. 1934. Vol IV. 80 p.

  20. Riquier C. Les Systèmes d’Equations aux Dérivées Partielles. Mémorial Sci. Math. XXXII, Gauthier-Villars, Paris. 1910.

  21. Harlow F.H., Welch J.E. Numerical Calculation of Time-Dependent Viscous Incompressible Flow of Fluid with Free // Phys. Fluids. 1965. № 8. P. 2182–2189.

  22. Buchberger B. Ein Algorithmus zum Auffinden der Basiselemente des Restklassenrings nach einem nulldimensionalen Polynomideal: Ph.D. thesis / Universiät Innsbruck. 1965.

  23. Taylor G.I. On the decay of vortices in a viscous fluid // Philosophical Magazine. 1923. V. 46. P. 671–674.

  24. Kovasznay L.I.G. Laminar flow behind a two-dimensional grid // Mathematical Proceedings of the Cambridge Philosophical Society. 1948. V. 44. № 1. P. 58–62.

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

Инструменты

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

Программирование

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