Conservative semi-Lagrangian method for continuity equation with various time steps in different parts of computational domain

Full Text

Введение

В настоящее время задача полноценного моделирования течения газа [1] из сопла ракеты с описанием всех возникающих физических процессов полностью не решена. Основным подходом для поиска решения дифференциальных уравнений, используемых при физико-математическом моделировании течения газа, является использование численных методов. Одним из препятствий в этом подходе является ограниченная вычислительная мощность вычислительных систем. Решение этой сложности состоит в увеличении мощности вычислительных систем и снижении вычислительной сложности используемых численных методов. Наиболее точно течение газа описывается системой дифференциальных уравнений Навье – Стоксах [2–6]. Эта система включает в себя уравнение неразрывности. Каждое уравнение из системы уравнений Навье – Стокса может быть рассмотрено в терминах конвективных и диффузионных слагаемых [7], а также правой части. Разработанный метод может быть использован для аппроксимации конвективных слагаемых системы уравнений Навье – Стокса. Остальные слагаемые из системы уравнений Навье – Стокса могут аппроксимироваться, например, в соответствии с технологией метода конечных элементов [2] или метода конечных разностей [8]. Идея использовать полулагранжевы методы (ранее эти методы назывались методами характеристик или методами траекторий) для описания течения газа появилась еще во второй половине XX в. [9–11]. В дальнейшем полулагранжевы методы начали использоваться при решение многих других задач [12–14]. Современные версии полулагранжевых методов являются консервативными [12; 15; 16]. При решении уравнения неразрывности под консервативностью, как правило, понимается выполнение закона сохранения для численного решения задачи в норме пространства L₁. [5; 15]. С точки зрения физики выполнение этого закона сохранения для численного решения уравнения неразрывности означат выполнение закона сохранения массы [1]. Для снижения вычислительной сложности расчетов строятся методы, использующие неравномерную (по пространству) вычислительную сетку [17; 18]. Это позволяет сократить время расчетов или использовать более мелкую сетку для поиска более точного решения. Несмотря на то, что при построении консервативных полулагранжевых методов дифференциальное уравнение неразрывности заменяется интегральным тождеством [15], использование неравномерной сетки, как правило, влечет за собой невозможность достижения локального закона сохранения для окрестностей узлов вычислительной сетки. Поэтому, чтобы добиться выполнения закона сохранения, используются поправочные коэффициенты, которые привносят в алгоритм дополнительную искусственную вязкость. Это существенно снижает точность вычислений в частях вычислительной области с высоким уровнем изменения численногорешения.

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

Постановка задачи

На множестве $D=\left[\mathrm{0,}T\right]\times \left[\mathrm{0,1}\right] , T>0 ,$рассмотрим одномерное уравнение неразрывности с ненулевой правой частью

$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u\right)}{\partial x}=f(t,x).$                                                                                                                    (1) 

Здесь $\rho \left(t,x\right)$ – искомая функция, а функции $u\left(t,x\right)$ и  $f\left(t,x\right)$ известны. Полагаем, что функции $\rho \left(t,x\right)$ и $u\left(t,x\right)$ достаточно гладкие. Для определенности, будем считать, что справедливы соотношения

$u(t,x)\left|{}_{x=0}\right.\ge 0 , u(t,x)\left|{}_{x=1}\right.\ge 0 .$                                                                                                   (2)

Для искомой функции $\rho \left(t,x\right)$ зададим начальное условие

$\rho (t,x)\left|{}_{t=0}\right.={\rho }_{\text{init}}\left(x\right)$,                                                                                                                         (3)

и граничное условие

$\rho (t,x)\left|{}_{x=0}={\rho }_{\text{left}}\left(t\right)\right..$                                                                                                                            (4)

Здесь ${\rho }_{\text{init}}\left(x\right), {\rho }_{\text{left}}\left(t\right)$,  – известные функции.

Множество значений $\mathrm{\Omega }=\left[0,1\right]$ переменной х разобьём на два подмножества ${\mathrm{\Omega }}^{\tau }=\left[0;\mathrm{0,4}\right]\cup \left[\mathrm{0,6};1\right]$ и ${\mathrm{\Omega }}^{2\tau }=\left(\mathrm{0,4};\mathrm{0,6}\right)$. На множестве $\left[\mathrm{0,}T\right]\times {\mathrm{\Omega }}^{\tau }$введем равномерную сетку с шагом $h=\frac{{\displaystyle 1}}{{\displaystyle N}}$, $N=10\cdot 2^n$, n – целое положительное число, по оси $Ox$ и шагом $\tau =\frac{{\displaystyle 1}}{{\displaystyle M}}, M -$четное положительное число, по оси Ot в следующем виде:

${\overline{D}}^{\tau }=\left\{(t_m,x_i):t_m=m\tau ,x_i=ih;i=\mathrm{0,...,}{N}_1,N_2\mathrm{,...,}N;m=\mathrm{0,...,}M\right\}$,

где  $N_1:N_{1}h=\mathrm{0,4}$  и  $N_2:N_{2}h=\mathrm{0,6}$.

На множестве $\left[\mathrm{0,}T\right]\times {\mathrm{\Omega }}^{2\tau }$ введем сетку с шагом $h=\frac{{\displaystyle 1}}{{\displaystyle N}}$ по оси $Ox$ и шагом $2\tau$ по оси $Ot$ в виде

${\overline{D}}^{2\tau }=\left\{(t_{2m},x_j):t_{2m}=2m\tau ,x_j=jh;j=N_1+\mathrm{1,...,}{N}_2-1;m=\mathrm{0,...,}M/2\right\}.$

Численное решение ${\rho }^h\left(t,x\right)$ задачи (1)–(4) будем искать во всех узлах сетки ${\overline{D}}^{h,\text{full}}={\overline{D}}^{\tau }\cup {\overline{D}}^{2\tau }$ (рис. 1).

 

Рис. 1. Объединенная вычислительная сетка

Fig. 1. Combined computational grid

 

Для этого на каждом слое по времени $t = t_m$ в каждой пространственной окрестности ${\mathrm{\Omega }}_k=\left[x_k-h/2,x_k+h/2\right]\cap \left[\mathrm{0,1}\right]$ узла $(t_m,x_k)$ сетки ${\overline{D}}^{h,\text{full}}$ искомую функцию ${\rho }^h\left(t,x\right)$ определим в виде кусочно-постоянной функции

${\rho }^h\left(t_m,x\right)={\rho }^h\left(t_m,x_k\right)\forall x\in {\mathrm{\Omega }}_k.$                                                                                                     (5)

В дальнейшем для удобства используем обозначения $g(t_m,x_k)=g_k^m$ и ${\rho }^h\left(t_m,x_k\right)={\rho }_k^{h,m}$. Для упрощения теоретических выкладок полагаем, что $\tau =с h$, где $c$ – некоторая константа такая, что

$с\le \frac{1}{4\max \left|u\left(t,x\right)\right|}$.                                                                                                                            (6)

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

Интегральное тождество во внутренних узлах сеток

Рассмотрим узел $\left(t_m,x_k\right), m>0; k>0$ при $k\ne N_1,k\ne N_1+1$ и $k\ne N_2-1, k\ne N_2$. В этом случае узел $\left(t_m,x_k\right)$ является внутренним узлом сетки ${\overline{D}}^{\tau }$ или сетки ${\overline{D}}^{2\tau }$. Обозначим $\left(t_m,x_l\right)$ и $\left(t_m,x_r\right)$ – две точки, которые являются границами множества ${\mathrm{\Omega }}_k$ на слое по времени $t_m$. Для каждой из этих точек в плоскости $\left(t,x\right)$ построим обратно по времени со слоя $t=t_m$ на слой $t=t_{m-s}$ характеристическую траекторию $C_q\left(t,\tilde{x}\left(t\right)\right), q=l,r$,(1). Здесь $s=1$, если $\left(t_m,x_k\right)\in {\overline{D}}^{\tau }$, и $s=2$, если $\left(t_m,x_k\right)\in {\overline{D}}^{2\tau }$ (рис. 2). Каждая траектория является решением задачи Коши для обыкновенного дифференциального уравнения

$\frac{d{\tilde{x}}_q}{dt}=u\left(t,{\tilde{x}}_q\right),\left[t\in t_{m-s},t_m\right],$                                                                                                          (7)

с начальными данными

${\tilde{x}}_q\left(t_m\right)=x_q, q=l,r.$                                                                                                                          (8)

В силу ограничения (6) несложно показать, что ${\tilde{x}}_l\left(t\right)\in \left(x_{k-1},x_k\right)$ и ${\tilde{x}}_r\left(t\right)\in \left(x_k,x_{k+1}\right)$ (рис. 2). Координаты траектории $C_q\left(t,\tilde{x}\left(t\right)\right)$ при $t=t_{m-s}$обозначим в виде $A_{k,q}^{m-s}\left(A_{k,q}^{t,m-s},A_{k,q}^{x,m-s}\right)=C_q\left(t_{m-s},{\tilde{x}}_q\left(t_{m-s}\right)\right), q=l,r$. Полагаем, что траектории $C_l$ и $C_r$ не пересекаются, поэтому $A_{k,l}^{x,m-s}<A_{k,r}^{x,m-s}$. Обозначим $G$ – множество, ограниченное сверху отрезком $t_m\times \left[x_l,x_r\right]$, траекториями $C_l$ и $C_r$ слева и справа, соответственно, и отрезком  снизу, как показано на рис. 2. В соответствии с теоремой Гаусса – Остроградского справедлива теорема 1.

 

Рис. 2. Характеристические траектории, выпущенные из границ окрестности

Fig. 2. Trajectories issued out from boundary of grid node neighborhood

 

Теорема 1. Решение задачи (1)–(4) удовлетворяет интегральному уравнению

$\underset{{\mathrm{\Omega }}_k}{\int }\rho \left(t_m,x\right)dx=\underset{A_{k,l}^{x,m-s}}{\overset{A_{k,r}^{x,m-s}}{\int }}\rho \left(t_{m-s},x\right)dx+\underset{G}{\int }f\left(t,x\right)dtdx$.                                                                           (9)

Построение численного решения во внутренних узлах сеток

Чтобы построить численный метод, аппроксимируем каждое слагаемое из тождества (9). Раскладывая функцию $\rho \left(t_m,x\right)$ в ряд Тейлора в окрестности точки $x_k$ по оси $Ox$, видно, что

$\underset{{\mathrm{\Omega }}_k}{\int }\rho \left(t_m,x\right)dx=\rho \left(t_m,x_k\right)h+\frac{{\partial }^2\rho }{\partial x^2}\left(t_m,x_k\right)\frac{h^3}{24}+O\left(h^5\right)$.                                                                      (10)

Для аппроксимации остальных слагаемых из (9) нам потребуется вычислить координаты точек $A_{k,q}^{m-s}\left(A_{k,q}^{t,m-s},A_{k,q}^{x,m-s}\right), q=l,r$, . Для этого мы решим задачу Коши (7)–(8) численным методом Ньютона

$A_{k,q}^{h,x,m-s}=x_q-s\tau u\left(t_m,x_q\right)\text{, }q=l,r$,                                                                                                    (11)

где $x_l=x_k-h/2 , x_r=x_k+h/2 .$ Таким образом, каждую траекторию мы аппроксимируем прямой линией, как показано на рис. 3. Определим величину ${\delta }_k^{m-s}$ такой, что

$\underset{A_{k,l}^{x,m-s}}{\overset{A_{k,r}^{x,m-s}}{\int }}\rho \left(t_{m-s},x\right)dx=\underset{A_{k,l}^{h,x,m-s}}{\overset{A_{k,r}^{h,x,m-s}}{\int }}\rho \left(t_{m-s},x\right)dx+{\delta }_k^{m-s}$.                                                                              (12)

 

Рис. 3. Аппроксимация траекторий движения прямыми линиями

Fig. 3.Trajectories approximation by straight lines

 

Справедлива следующая теорема.

 

Теорема 2. Для приближенного решения (11) задачи (7)–(8) верно

 

${\delta }_k^{m-s}=O\left(h{\tau }^2+{\tau }^3\right), q=l,r$ .                                                                                                           (13)

Доказательство.

Не вдаваясь во все точные детали доказательства, представим его в коротком виде. Разложим функцию ${\tilde{x}}_q\left(t\right)$ в окрестности точки $t=t_m$ в ряд Тейлора

${\tilde{x}}_q\left(t_m-s\tau \right)={\tilde{x}}_q\left(t_m\right)-s\tau \frac{d{\tilde{x}}_q}{dt}\left(t_m\right)+\frac{d^2{\tilde{x}}_q}{d{t}^2}\left({\theta }_q\right)\frac{{\left(s\tau \right)}^2}{2},$  где  ${\theta }_q\in \left[t_m-s\tau ,t_m\right]$.

С учетом тождества $\frac{d{\tilde{x}}_q}{dt}\left(t_m\right)=u\left(t_m,x_q\right)$, а также с учетом малости величин $\tau$ и $h$ понятно, что знак величины $A_{k,q}^{x,m-s}-A_{k,q}^{h,x,m-s}$ определяется знаком величины $\frac{d^2{\tilde{x}}_q}{d{t}^2}\left(t_m\right)$.

Если величины $\frac{d^2{\tilde{x}}_l}{d{t}^2}\left(t_m\right)$ и $\frac{d^2{\tilde{x}}_r}{d{t}^2}\left(t_m\right)$ имеют одинаковый знак, например, положительный (см. рис. 3), то

${\delta }_k^{m-s}=\underset{A_{k,r}^{h,x,m-s}}{\overset{A_{k,r}^{x,m-s}}{\int }}\rho \left(t_{m-s},x\right)dx-\underset{A_{k,l}^{h,x,m-s}}{\overset{A_{k,l}^{x,m-s}}{\int }}\rho \left(t_{m-s},x\right)dx$.

В этом случае несложно показать, что ${\delta }_k^{m-s}=O\left(h{\tau }^2+{\tau }^3\right)$. Если величины $\frac{d^2{\tilde{x}}_l}{d{t}^2}\left(t_m\right)$ и $\frac{d^2{\tilde{x}}_r}{d{t}^2}\left(t_m\right)$ имеют разные знаки, то существует некоторая точка ${\tilde{x}}_k^{*}\in \left[x_l,x_r\right]$ такая, что $\frac{d^2{\tilde{x}}_k^{*}}{d{t}^2}\left(t_m\right)=0$. Тогда, используя разложение в ряд Тейлора, можно показать, что

$\underset{A_{k,l}^{h,x,m-s}}{\overset{A_{k,l}^{x,m-s}}{\int }}\rho \left(t_{m-s},x\right)dx=O\left(h{\tau }^2+{\tau }^3\right)$   и   $\underset{A_{k,r}^{h,x,m-s}}{\overset{A_{k,r}^{x,m-s}}{\int }}\rho \left(t_{m-s},x\right)dx=O\left(h{\tau }^2+{\tau }^3\right)$.

Тогда из (12) вытекает, что ${\delta }_k^{m-s}=O\left(h{\tau }^2+{\tau }^3\right)$. Теорема 2 доказана.

Аналогично теореме 2, из разложения в ряд Тейлора функции $f\left(t,x\right)$ можно показать, что при использовании соотношения $\tau =сh$ справедливо

$\underset{G}{\int }f\left(t,x\right)dxdt=\underset{t_m-s\tau }{\overset{t_m}{\int }}dt\underset{{\tilde{x}}_l\left(t\right)}{\overset{{\tilde{x}}_r\left(t\right)}{\int }}f\left(t,x\right)dx=s\tau \text{meas}\left({\mathrm{\Omega }}_k\right)f\left(t_m,x_i\right)+O\left(h^3\right)$,                                                  (14)

где $\text{meas}\left({\mathrm{\Omega }}_k\right)$ – длинна окрестности ${\mathrm{\Omega }}_k$, которая определяется в виде

$\text{mes}\left({\mathrm{\Omega }}_k\right)=\left\{\begin{array}{c}h,\text{ если }k=\mathrm{1,...,}N-1;\\ h/\mathrm{2,}\text{ если }k=\mathrm{0,}k=N.\end{array}\right.$

Подставив (10), (12), (14) в (9), мы получим формулу для вычисления численного решения  ${\rho }_k^{h,m}$

${\rho }_k^{h,m}=\frac{1}{\text{meas}\left({\mathrm{\Omega }}_k\right)}\underset{A_{k,l}^{h,x,m-s}}{\overset{A_{k,r}^{h,x\mathrm{,,}m-s}}{\int }}{\rho }^h\left(t_{m-s},x\right)dx+s\tau f\left(t_m,x_i\right)$.                                                                      (15)

Напомним, что функция ${\rho }^h\left(t_{m-s},x\right)$ является кусочно-постоянной, поэтому ее интеграл легко вычисляется аналитически. Исходя из построения алгоритма, формула (15) применима для всех узлов сетки ${\overline{D}}^{h,\text{full}}$ при нечетном числе $m$. Если же $m$ четное число, то формулу (15) можно применять .

$\forall k>0:k\ne N_1,k\ne N_1+\mathrm{1,}k\ne N_2-\mathrm{1,}k\ne N_2$

Построение численного решения на границах между сетками

Рассмотрим случай границы между сетками ${\overline{D}}^{\tau }$ и ${\overline{D}}^{2\tau }$ при $x=\mathrm{0,6}$. На границе $x=\mathrm{0,4}$ алгоритм строится аналогично. Определим способ вычисления численного решения для ${\rho }_{N_2}^{h,m}$ и ${\rho }_{N_2-1}^{h,m}$ при четном числе $m$. В зависимости от знака величины $u\left(t_m,x_{N_2}-h/2\right)$ возможны три случая. В первом случае, когда $u\left(t_m,x_{N_2}-h/2\right)=0$, для вычисления ${\rho }_{N_2-1}^{h,m}$,  ${\rho }_{N_2}^{h,m}$ следует использовать формулу (15).

Пусть $u\left(t_m,x_{N_2}-h/2\right)>0$, тогда $A_{N_2,l}^{h,x,m-1}<x_{N_2}-h/2$, как это показано на рис. 4. В этом случае величина ${\rho }_{N_2-1}^{h,m}$ определяется по формуле (15) при $s = 2$. Чтобы вычислить значение ${\rho }_{N_2}^{h,m}$, также используем формулу (15) при $s = 1$ и получим

${\rho }_{N_2}^{h,m}=\frac{1}{h}\underset{A_{N_2,l}^{h,x,m-1}}{\overset{A_{N_2,r}^{h,x,m-1}}{\int }}{\rho }^h\left(t_{m-1},x\right)dx+\tau f\left(t_m,x_{N_2}\right)$.                                                                                     (16)

 

Рис. 4. Случай положительного значения функции скорости на границе двух сеток

Fig. 4. The case of positive value of velocity function on boundary

 

Функция ${\rho }^h\left(t_{m-1},x\right)$ не определена (ранее не вычислена) на всем множестве $\left[A_{N_2,l}^{h,x,m-1},A_{N_2,r}^{h,x,m-1}\right]$, поэтому вычисление этого интеграла невозможно. Множество $\left[A_{N_2,l}^{h,x,m-1},A_{N_2,r}^{h,x,m-1}\right]$ разобьем на два отрезка $\left[x_{N_2}-h/\mathrm{2,}{A}_{N_2,r}^{h,x,m-1}\right]$, на котором функция ${\rho }^h\left(t_{m-1},x\right)$ определена и оставшийся отрезок $\left[A_{N_2,l}^{h,x,m-1},x_{N_2}-h/2\right]$ (рис. 4). Однако, по аналогии с теоремой 1, можно показать, что

$\underset{A_{N_2,l}^{h,x,m-1}}{\overset{x_{N_2}-h/2}{\int }}\rho \left(t_{m-1},x\right)dx=\underset{A_l^{m-1}}{\overset{A_r^{m-2}}{\int }}\rho \left(t_{m-2},x\right)dx+\underset{G_{N_2}^{m-1}}{\int }f\left(t,x\right)dtdx$.                                                          (17)

Здесь величина $A_r^{m-1}\approx A_{N_2,l}^{h,x,m-2}$, $A_l^{m-2}\approx A_{N_2-\mathrm{1,}r}^{h,x,m-2}$. В итоге подставим (17) в (16) и получим

${\rho }_{N_2}^{h,m}=\frac{1}{h}\underset{x_{N_2}-h/2}{\overset{A_{N_2,r}^{h,x,m-1}}{\int }}{\rho }^h\left(t_{m-1},x\right)dx+\frac{1}{h}\underset{A_{N_2-\mathrm{1,}r}^{h,x,m-2}}{\overset{A_{N_2,l}^{h,x,m-2}}{\int }}{\rho }^h\left(t_{m-2},x\right)dx+\tau f\left(t_m,x_{N_2}\right)$.                                          (18)

Рассмотрим второй случай, когда $u\left(t_m,x_{N_2}-h/2\right)<0$ и $A_{N_2,l}^{h,x,m-1}>x_{N_2}-h/2$ (рис. 5). В этом случае величина ${\rho }_{N_2}^{h,m}$ определяется по формуле (15) при $s=1$. Чтобы вычислить значение ${\rho }_{N_2-1}^{h,m}$, используем формулу (15) при $s=1$ и получим

${\rho }_{N_2-1}^{h,m}=\frac{1}{h}\underset{A_{N_2-\mathrm{1,}l}^{h,x,m-1}}{\overset{A_{N_2-\mathrm{1,}r}^{h,x,m-1}}{\int }}{\rho }^h\left(t_{m-1},x\right)dx+\tau f\left(t_m,x_{N_2}\right)$.                                                                                 (19)

 

Рис. 5. Случай положительного значения функции скорости на границе двух сеток

Fig. 5. The case of negative value of velocity function on boundary

 

Отрезок $\left[A_{N_2-\mathrm{1,}l}^{h,x,m-1},A_{N_2-\mathrm{1,}r}^{h,x,m-1}\right]$ разложим на две части: $\left[A_{N_2-\mathrm{1,}l}^{h,x,m-1},x_{N_2-1}+h/2\right]$ и $\left[x_{N_2-1}+h/\mathrm{2,}{A}_{N_2-\mathrm{1,}r}^{h,x,m-1}\right]$. Значение функции  на отрезке $\left[x_{N_2-1}+h/\mathrm{2,}{A}_{N_2-\mathrm{1,}r}^{h,x,m-1}\right]$ определено. Вновь, по аналогии с теоремой 1, можно показать, что

$\underset{A_{N_2-\mathrm{1,}l}^{h,x,m-1}}{\overset{x_{N_2-1}+h/2}{\int }}\rho \left(t_{m-1},x\right)dx=\underset{{\tilde{A}}_l^{m-2}}{\overset{{\tilde{A}}_r^{m-2}}{\int }}\rho \left(t_{m-2},x\right)dx+\underset{{\tilde{G}}_{N_2}^{m-1}}{\int }f\left(t,x\right)dtdx$.                                                       (20)

Причем, ${\tilde{A}}_l^{m-2}\approx A_{N_2-\mathrm{1,}l}^{h,x,m-2}$, ${\tilde{A}}_r^{m-1}\approx A_{N_2,l}^{h,x,m-2}$. Используя (19) и (20) определим

${\rho }_{N_2-1}^{h,m}=\frac{1}{h}\underset{x_{N_2-1}+h/2}{\overset{A_{N_2-\mathrm{1,}r}^{h,x,m-1}}{\int }}{\rho }^h\left(t_{m-1},x\right)dx+\frac{1}{h}\underset{A_{N_2-\mathrm{1,}l}^{h,x,m-2}}{\overset{A_{N_2,l}^{h,x,m-2}}{\int }}{\rho }^h\left(t_{m-2},x\right)dx+2\tau f\left(t_m,x_{N_2-1}\right)$.                             (21)

Закон сохранения и сходимость

Проверку выполнения закона сохранения выполним только на четных слоях по времени, поскольку только на этих временных слоях численное решение ${\rho }^h$ определено $\forall x\in \left[\mathrm{0,1}\right]$. Классический закон сохранения массы при переходе с текущего слоя по времени на следующий слой с учетом границ втекания, вытекания и функции правой части (внутренних источников) имеет следующий вид [1]:

$\underset{0}{\overset{1}{\int }}\rho \left(t_m,x\right)dx=\underset{0}{\overset{1}{\int }}{\rho }^h\left(t_{m-2},x\right)dx+\underset{t_{m-2}}{\overset{t_m}{\int }}\rho \left(t\mathrm{,0}\right)u\left(t\mathrm{,0}\right)dt-\underset{t_{m-2}}{\overset{t_m}{\int }}\rho \left(t\mathrm{,1}\right)u\left(t\mathrm{,1}\right)dt+\underset{t_{m-2}}{\overset{t_m}{\int }}dt\underset{0}{\overset{1}{\int }}f\left(t,x\right)dx.$    (22)

Покажем, что тождество (22) справедливо для построенного численного решения ${\rho }^h\left(t,x\right)$.

В соответствии с (15) выпишем формулу для вычисления численного решение на слое $t=t_{m-1}$ в узлах сетки  ${\overline{D}}^{\tau }. \forall k=\mathrm{1,...,}{N}_1,N_2\mathrm{,...,}N$:

${\rho }_k^{h,m-1}=\frac{1}{\text{meas}\left({\Omega }_k\right)}\underset{A_{k,l}^{h,x,m-2}}{\overset{A_{k,r}^{h,x\mathrm{,,}m-2}}{\int }}{\rho }^h\left(t_{m-2},x\right)dx+\tau f\left(t_{m-1},x_k\right)$.                                                             (23)

Просуммируем (23) по всем $k=\mathrm{1,...,}{N}_1,N_2\mathrm{,...,}N$. С учетом того, что ${\rho }^h\left(t,x\right)$ является кусочно-постоянной функцией в каждой окрестности узла $x_k$, мы получим

$\underset{h/2}{\overset{0.4+h/2}{\int }}{\rho }^h\left(t_{m-1},x\right)dx+\underset{0.6-h/2}{\overset{1}{\int }}{\rho }^h\left(t_{m-1},x\right)dx=$

$=\underset{A_{\mathrm{1,}l}^{h,x,m-2}}{\overset{A_{N_1,r}^{h,x\mathrm{,,}m-2}}{\int }}{\rho }^h\left(t_{m-2},x\right)dx+\underset{A_{N_2,l}^{h,x,m-2}}{\overset{A_{N,r}^{h,x\mathrm{,,}m-2}}{\int }}{\rho }^h\left(t_{m-2},x\right)dx+\tau \underset{h/2}{\overset{0.4+h/2}{\int }}{f}^h\left(t_{m-1},x_k\right)dx+\tau \underset{0.6-h/2}{\overset{1}{\int }}{f}^h\left(t_{m-1},x\right)dx.$  (24)

Аналогично теореме 1, можно показать, что на промежутке времени $\left[t_{m-2},t_{m-1}\right]$ справедливы тождества, учитывающие втекание и вытекание субстанции

$\underset{0}{\overset{h/2}{\int }}\rho \left(t_{m-1},x\right)dx=\underset{t_{m-2}}{\overset{t_{m-1}}{\int }}\rho \left(t\mathrm{,0}\right)u\left(t\mathrm{,0}\right)dt+\underset{0}{\overset{A_{\mathrm{0,}r}^{x,m-2}}{\int }}\rho \left(t_{m-2},x\right)dx+\underset{t_{m-2}}{\overset{t_{m-1}}{\int }}dt\underset{0}{\overset{{\tilde{x}}_r\left(t\right)}{\int }}f\left(t,x\right)dx,$

$\underset{A_{N,r}^{x,m-2}}{\overset{1}{\int }}\rho \left(t_{m-2},x\right)dx=\underset{t_{m-2}}{\overset{t_{m-1}}{\int }}\rho \left(t\mathrm{,1}\right)u\left(t\mathrm{,1}\right)dt+\underset{t_{m-2}}{\overset{t_{m-1}}{\int }}dt\underset{{\tilde{x}}_l\left(t\right)}{\overset{1}{\int }}f\left(t,x\right)dx.$                                                         (25)

С учетом (25) определим ${\rho }^h\left(t_{m-1},x_0\right)$ и ${\rho }^h\left(t\mathrm{,1}\right)u^h\left(t\mathrm{,1}\right)$ такими, что

$\frac{h}{2}{\rho }^h\left(t_{m-1},x_0\right)=\underset{t_{m-2}}{\overset{t_{m-1}}{\int }}{\rho }^h\left(t\mathrm{,0}\right)u^h\left(t\mathrm{,0}\right)dt+\underset{0}{\overset{A_{\mathrm{0,}r}^{h,x,m-2}}{\int }}{\rho }^h\left(t_{m-2},x\right)dx+\tau \frac{h}{2}f\left(t_{m-1},x_0\right),$  

$\underset{t_{m-2}}{\overset{t_{m-1}}{\int }}{\rho }^h\left(t\mathrm{,1}\right)u^h\left(t\mathrm{,1}\right)dt=\underset{A_{N,r}^{h,x,m-2}}{\overset{1}{\int }}{\rho }^h\left(t_{m-2},x\right)dx.$                                                                                 (26)

Добавим тождества (26) в (24) и получим

$\underset{0}{\overset{0.4+h/2}{\int }}{\rho }^h\left(t_{m-1},x\right)dx+\underset{0.6-h/2}{\overset{1}{\int }}{\rho }^h\left(t_{m-1},x\right)dx=$

$=\underset{t_{m-2}}{\overset{t_{m-1}}{\int }}{\rho }^h\left(t\mathrm{,0}\right)u^h\left(t\mathrm{,0}\right)dt+\underset{0}{\overset{A_{N_1,r}^{h,x\mathrm{,,}m-2}}{\int }}{\rho }^h\left(t_{m-2},x\right)dx+\underset{A_{N_2,l}^{h,x,m-2}}{\overset{1}{\int }}{\rho }^h\left(t_{m-2},x\right)dx-$

$-\underset{t_{m-2}}{\overset{t_{m-1}}{\int }}{\rho }^h\left(t\mathrm{,1}\right)u^h\left(t\mathrm{,1}\right)dt+\tau \underset{0}{\overset{0.4+h/2}{\int }}{f}^h\left(t_{m-1},x_k\right)dx+\tau \underset{0.6-h/2}{\overset{1}{\int }}{f}^h\left(t_{m-1},x\right)dx.$

С одной стороны, значения функции $\rho (t,x)$ на границе $x=0$ нам известны, а с другой – для выполнения закона сохранения мы использовали ${\rho }^h\left(t_{m-1}\mathrm{,0}\right)$ из первого соотношения в (26). Чтобы одновременно добиться выполнения сохранения и использовать точное решение на границе втекания, определим величину

$E\left(t_{k-1}\mathrm{,0}\right)=\rho \left(t_{k-1}\mathrm{,0}\right)-{\rho }^h\left(t_{k-1}\mathrm{,0}\right).$                                                                                                (27)

Тогда использование точного решение на границе втекания приводит закон сохранения к виду

$\underset{0}{\overset{0.4+h/2}{\int }}{\rho }^h\left(t_{m-1},x\right)dx+\underset{0.6-h/2}{\overset{1}{\int }}{\rho }^h\left(t_{m-1},x\right)dx=$

$=\underset{t_{m-2}}{\overset{t_{m-1}}{\int }}{\rho }^h\left(t\mathrm{,0}\right)u^h\left(t\mathrm{,0}\right)dt+\underset{0}{\overset{A_{N_1,r}^{h,x\mathrm{,,}m-2}}{\int }}{\rho }^h\left(t_{m-2},x\right)dx+\underset{A_{N_2,l}^{h,x,m-2}}{\overset{1}{\int }}{\rho }^h\left(t_{m-2},x\right)dx+\frac{h}{2}E\left(t_{m-1}\mathrm{,0}\right)-$

$-\underset{t_{m-2}}{\overset{t_{m-1}}{\int }}{\rho }^h\left(t\mathrm{,1}\right)u^h\left(t\mathrm{,1}\right)dt+\tau \underset{0}{\overset{0.4+h/2}{\int }}{f}^h\left(t_{m-1},x_k\right)dx+\tau \underset{0.6-h/2}{\overset{1}{\int }}{f}^h\left(t_{m-1},x\right)dx.$                                        (28)

Проведя аналогичные выкладки для численного решения на слое $t=t_m$ в узлах сетки  ${\overline{D}}^{\tau }$, $\forall k=\mathrm{1,...,}{N}_1,N_2\mathrm{,...,}N$ мы получим формулу аналогичную (28). Для внутренних узлов ${\overline{D}}^{2\tau }$ при $k=N_1+\mathrm{2,...,}{N}_2-2$ после суммирования (15), получим

$\underset{N_1+2-h/2}{\overset{N_2-2+h/2}{\int }}{\rho }^h\left(t_m,x\right)dx=\underset{A_{N_1,l}^{h,x,m-2}}{\overset{A_{N_2,r}^{h,x,m-2}}{\int }}{\rho }^h\left(t_{m-2},x\right)dx+2\tau \underset{N_1+2-h/2}{\overset{N_2-2+h/2}{\int }}{f}^h\left(t_m,x\right)dx$.                                         (29)

С учетом специфики построения формул для вычисления численного решения на четном слое $t=t_m$ для $k=N_1,N_1+\mathrm{1,}{N}_2-\mathrm{1,}{N}_2$ и используя (28)–(29), получим

$\underset{0}{\overset{1}{\int }}{\rho }^h\left(t_m,x\right)dx=\underset{t_{m-2}}{\overset{t_{m-1}}{\int }}{\rho }^h\left(t\mathrm{,0}\right)u^h\left(t\mathrm{,0}\right)dt+\underset{0}{\overset{1}{\int }}{\rho }^h\left(t_{m-2},x\right)dx+$

$+\frac{h}{2}E\left(t_m\mathrm{,0}\right)+\frac{h}{2}E\left(t_{m-1}\mathrm{,0}\right)-\underset{t_{m-2}}{\overset{t_m}{\int }}{\rho }^h\left(t\mathrm{,1}\right)u^h\left(t\mathrm{,1}\right)dt+\underset{t_{m-2}}{\overset{t_m}{\int }}dt\underset{0}{\overset{1}{\int }}{f}^h\left(t,x\right)dx.$                                         (30)

Тождество (30) представляет собой закон сохранения для искомой функции ${\rho }^h$ при переходе с временного слоя $t=t_{m-2}$ на слой $t=t_m$ с учетом втекания через границу $x=0$, вытекания через границу $x=1$, внутренних источников, описываемых функций $f\left(t,x\right)$ и корректировки закона сохранения в виду использования точного решения на границе втекания.

Для проверки сходимости численного решения к точному решению задачи (1)–(4) используем норму пространства L₁:

$er{r}_n\left(t_m\right)={||{\rho }^h\left(t_m,x\right)-\rho \left(t_m,x\right)||}_{L_1}=\sum _{k=0}^N\text{meas}\left({\Omega }_k\right)\left|{\rho }_k^{h,m}-{\rho }_k^m\right|$.                         (31)

С учетом выполнения закона сохранения и наличия аппроксимации интегрального тождества (9) с точностью $O\left(h^3\right)$ при выполнении равенства $\tau =сh$ несложно показать [13], что численное решение задачи (1)–(4) сходится к точному решению с первым порядком точности.

Результаты вычислительного эксперимента

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

$\rho \left(t,x\right)=\mathrm{1,1}+\sin \left(t x\right), u\left(t,x\right)=\mathrm{0,6}+\mathrm{0,5}\cos \left(\text{2}\pi x\right).$                                                                         (32)

Используя функции (32), были заданы функции ${\rho }_{\text{init}}\left(x\right)$, ${\rho }_{\text{left}}\left(t\right)$ и $f\left(t,x\right)$. После этого была сформулирована задача (1)–(4). Решение задачи искалось на серии сеток при $N=10\cdot 2^n$, $n=\mathrm{1,...,6}$. Для удовлетворения ограничения (6) использовалось тождество $\tau =\mathrm{0,2}h$. Для оценки погрешности $er{r}_n\left(t_m\right)$ численного решения используем норму (31). Порядок сходимости численного решения к точному решению оценивается параметром ${\text{conv}}_n={\log }_2\left(er{r}_{n-1}\left(T\right)/er{r}_n\left(T\right)\right)$. Результаты расчетов приведены в табл. 1. Из таблицы видно, что алгоритм имеет первый порядок сходимости.

 

Таблица 1. О порядке сходимости разработанного метода

n	N	M	errn(T)	convn
1	20	100	0,0234467
2	40	200	0,0103564	1,17646
3	80	400	0,0048117	1,10544
4	160	800	0,0023115	1,05774
5	320	1600	0,0011317	1,03027
6	640	3200	0,0005598	1,01542

 

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

— масса газа на слое по времени  $t=t_m$

${\text{mass}}_m=\underset{0}{\overset{1}{\int }}{\rho }^h\left(t_m,x\right)dx=\sum _{k=0}^N\text{meas}\left({\mathrm{\Omega }}_k\right){\rho }_k^{h,m}$;

— масса жидкости, которая втекла через границу втекания $x=0$ за время $\left[t_{m-2},t_m\right]$,

${\text{input}}_m=\frac{\tau }{2}{\rho }_{\text{left}}\left(t_{m-2}\right)u\left(t_{m-2}\mathrm{,0}\right)+\tau {\rho }_{\text{left}}\left(t_{m-1}\right)u\left(t_{m-1}\mathrm{,0}\right)+\frac{\tau }{2}{\rho }_{\text{left}}\left(t_m\right)u\left(t_m\mathrm{,0}\right)$;

— масса жидкости, которая вытекла через границу вытекания $x=1$ за время $\left[t_{m-2},t_m\right]$,

${\text{output}}_m=\underset{A_{N,r}^{h,x,m-2}}{\overset{1}{\int }}{\rho }^h\left(t_{m-2},x\right)dx+\underset{A_{N,r}^{h,x,m-1}}{\overset{1}{\int }}{\rho }^h\left(t_{m-1},x\right)dx$;

— масса жидкости, которая втекла из внутренних источников за время $\left[t_{m-2},t_m\right]$,

${\text{source}}_m=\tau \sum _{q=т-1}^m\sum _{k=1}^{N_1-1}\text{meas}\left({\Omega }_k\right)f_k^q+\tau \sum _{q=т-1}^m\sum _{k=N_2}^N\text{meas}\left({\Omega }_k\right)f_k^q+2\tau h\sum _{k=N_1+1}^{N_2-1}{f}_k^m$;

— коэффициент корректировки в связи с использованием точного решения на границе втекания

${\text{corr}}_{\text{m}}=\frac{h}{2}\rho \left(t_{k-1}\mathrm{,0}\right)+\frac{h}{2}\rho \left(t_k\mathrm{,0}\right)-{\text{input}}_m$;

— коэффициент, показывающий погрешность выполнения закона сохранения,

${\text{conlow}}_m={\text{mass}}_m-\left({\text{mass}}_{m-2}{\text{+input}}_m-{\text{output}}_m+{\text{source}}_m+{\text{corr}}_{\text{m}}\right)$.

Результаты расчетов для $N = 20$ представлены в табл. 2. Как видно из расчетов, закон сохранения массы выполняется при переходе с каждого четного слоя по времени на следующий четный слой с точностью до представления нуля в машинной арифметике.

 

Таблица 2. О выполнении закона сохранения 

m	massm	massm-2	inputm	outputm	sourcem	corrm	conlowm
2	1,124607	1,10	0,0605	0,060827	0,026348	–0,001415	3,13E-15
4	1,149764	1,124607	0,0605	0,062812	0,029021	–0,001552	1,88E-15
6	1,174822	1,149764	0,0605	0,065404	0,031652	–0,001690	2,06E-15
8	1,199616	1,174822	0,0605	0,068115	0,034235	–0,001827	–4,20E-15
10	1,224083	1,199616	0,0605	0,070831	0,036763	–0,001965	8,02E-17
…
94	1,481967	1,463581	0,0605	0,099126	0,060626	–0,003615	6,94E-15
96	1,499675	1,481967	0,0605	0,100956	0,061916	–0,003752	1,39E-15
98	1,516679	1,499675	0,0605	0,102682	0,063077	–0,003891	–2,49E-15
100	1,532958	1,516679	0,0605	0,1043	0,064106	–0,004027	3,51E-15

 

Заключение

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

Благодарности. Работа поддержана Красноярским математическим центром, финансируемым Минобрнауки РФ в рамках мероприятий по созданию и развитию региональных НОМЦ (соглашение 075-02-2023-912).

Acknowledgment. This work is supported by the Krasnoyarsk Mathematical Center and financed by the Ministry of Science and Higher Education of the Russian Federation in the framework of the establishment and development of regional Centers for Mathematics Research and Education (Agreement No. 075-02-2023-912)

About the authors

Aleksandr V. Vyatkin

Institute of Computational Modeling SB RAS

Author for correspondence.
Email: vyatkin@icm.krasn.ru

Cand. Sc., scientific researcher

Russian Federation, 50/44, Akademgorodok St., Krasnoyarsk, 660036

Aleksand V. Maltsev

Siberian Federal University

Email: maltsev@mail.ru

graduate student

Russian Federation, 79, Svobodny Av., Krasnoyarsk, 660041

References

Supplementary files