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

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

М. В. Абгарян 1, А. М. Бишаев 2*, В. А. Рыков 3

1 МАИ
125993 Москва, Волоколамское шоссе, 4, Россия

2 МФТИ
141701 М.о., Долгопрудный, Институтский пер. 9, Россия

3 ФИЦ ВЦ РАН
119333 Москва, ул. Вавилова, 44, кор. 2, Россия

* E-mail: bishaev@mail.ru

Поступила в редакцию 15.02.2020
После доработки 12.03.2020
Принята к публикации 09.04.2020

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

Аннотация

Исследуется построенная ранее кинетическая модель для описания поведения неидеального газа. Проводится более полная оценка безразмерных параметров, показывающих, когда необходимо учитывать неидеальность газа. Обнаружено, что в интеграле столкновений для частиц в связанных состояниях можно провести интегрирование по пространству скоростей, что позволяет существенно упростить исходную систему уравнений и доказать Н-теорему. Полученная система приводится к безразмерному виду. Для ее решения предлагается консервативная численная схема. Библ. 11. Фиг. 2.

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

ВВЕДЕНИЕ

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

1. ОСНОВНЫЕ УРАВНЕНИЯ

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

$U(r) = \left[ \begin{gathered} + \infty ,\quad 0 \leqslant r \leqslant d, \hfill \\ - {{U}_{0}}{{\left( {\frac{d}{r}} \right)}^{4}},\quad r > d. \hfill \\ \end{gathered} \right.$

Это один из типов потенциала Сазерленда. Ясно, что молекулы, двигающиеся со скоростями $\vec {\xi }$ и ${{\vec {\xi }}_{1}}$ такие, что $\frac{{m{{g}^{2}}}}{4} + U(r) < 0$, $g = \left| {{{{\vec {\xi }}}_{1}} - \vec {\xi }} \right|$, образуют связанные состояния и, в отличие от случая, когда потенциал взаимодействия сводится только к отталкиванию, никогда не разойдутся. Система кинетических уравнений, о которой будет идти речь, была получена методом Боголюбова из BBGKY цепочки с учетом образования связанных состояний. В этом случае для функции распределения $f(t,\vec {x},\vec {\xi })$ получается следующее уравнение:

$\frac{{\partial f}}{{\partial t}} + {{\xi }^{i}}\frac{{\partial f}}{{\partial {{x}^{i}}}} = {{J}_{B}} - {{J}_{{{\text{del}}}}},\quad {{J}_{B}} = {{d}^{2}}\int {d{{\xi }_{1}}} \left( {\int\limits_0^{ + \infty } {\int\limits_{}^{2\pi } {g(f{\text{'}}f_{1}^{'} - f{{f}_{1}})} } bdbd\theta } \right),$
где
${{J}_{{{\text{del}}}}} = 2{{\chi }_{1}}{{d}^{2}}\int\limits_{g \leqslant {{\chi }_{1}}} {d\vec {g}} \int\limits_0^\pi {\int\limits_0^{2\pi } {\left| {\cos \omega } \right|} {\kern 1pt} [{{F}_{2}}]} \sin \omega d\omega d\theta ,\quad {{\chi }_{1}} = \sqrt {\frac{{{{U}_{0}}}}{m}} ,$
$[{{F}_{2}}] = {{F}_{{2f}}} - {{F}_{{2d}}}$ – разрыв двухчастичной функции распределения ${{F}_{2}}(t,\vec {x},\vec {w},\vec {r},\vec {g})$ на поверхности $\gamma :\left\{ {\frac{m}{4}{{g}^{2}} - {{U}_{0}}{{{\left( {\frac{d}{r}} \right)}}^{4}} = 0} \right\}$. Приведенное уравнение не является замкнутым, поэтому для построения замкнутой системы кинетических уравнений в [1] вводится ряд предположений. Предполагается, что внутри поверхности $\gamma $ молекулы находятся в связанном состоянии и их двухчастичная функция распределения представляется в виде На внешнюю поверхность $\gamma $ приходят молекулы в свободном состоянии. Так как двухчастичная функция сохраняется на поверхности $\gamma $ и на бесконечности имеет место условие молекулярного хаоса, то получаем
${{F}_{{2f}}} = {{\left. {f(t,\vec {x},\vec {\xi })f(t,\vec {x},\vec {\xi } + \vec {g})} \right|}_{{g = 0}}} = {{f}^{2}}(t,\vec {x},\vec {\xi }) = {{f}^{2}}(t,\vec {x},\vec {w}).$
Последняя формула получилась из-за того, что скорость центра масс молекул есть $\vec {w} = \vec {\xi } - \frac{{\vec {g}}}{2}$.

Чтобы получить замкнутую систему кинетических уравнений, вводятся функции $s(t,\vec {x},\vec {w})$ и $h(t,\vec {x},\vec {w})$:

$s = s(t,\vec {x},\vec {w}) = \iint\limits_{{{\Omega }_{{{\text{del}}}}}} {{{{\bar {F}}}_{{2d}}}d\vec {r}d\vec {g},\quad h = h(t,\vec {x},\vec {w}) = \iint\limits_{{{\Omega }_{{{\text{del}}}}}} {\left( {\frac{{m{{g}^{2}}}}{4} - {{U}_{0}}{{{\left( {\frac{d}{r}} \right)}}^{4}}} \right){{{\bar {F}}}_{{2d}}}}}d\vec {r}d\vec {g} \leqslant 0,$
где ${{\Omega }_{{{\text{del}}}}} = \left\{ {0 \leqslant g \leqslant 2\frac{{{{\chi }_{1}}}}{{{{r}^{2}}}},\;d \leqslant r < + \infty } \right\}$ – область фазового пространства, в которой молекулы находятся в связанном состоянии (см. [1]). Нетрудно видеть, что $s$ $\Delta \vec {x}\Delta \vec {w}$ есть число молекул, находящихся в связанном состоянии, чьи координаты и скорости находятся в соответствующем элементе фазового пространства, а $h\Delta \vec {x}\Delta \vec {w}$ – их энергия в системе координат, связанной с центром масс этих молекул, который движется со скоростью $\vec {w}$. Непосредственно из определения функций $s(t,\vec {x},\vec {w})$ и $h(t,\vec {x},\vec {w})$ следует, что
$s(t,\vec {x},\vec {w}) = 16{{\pi }^{2}}{{f}_{0}}(t,\vec {x},\vec {w}){{d}^{3}}{{Z}_{s}}(\bar {\chi })\chi _{1}^{{3/2}},\quad h(t,\vec {x},\vec {w}) = 16{{\pi }^{2}}{{f}_{0}}(t,\vec {x},\vec {w}){{d}^{3}}{{U}_{0}}{{Z}_{h}}(\bar {\chi })\chi _{1}^{{3/2}},$
где
$\begin{gathered} {{Z}_{s}}(\bar {\chi }) = \frac{{16}}{{105}}\left( {\frac{{35}}{6} - \bar {\chi } + \frac{2}{{11}}{{{\bar {\chi }}}^{2}} - \frac{4}{{135}}{{{\bar {\chi }}}^{3}} + \; \ldots } \right), \\ {{Z}_{h}}(\bar {\chi }) = \frac{{16}}{{105}}\left( {1 - \frac{4}{{11}}\bar {\chi } + \frac{4}{{15}}{{{\bar {\chi }}}^{2}} - \frac{{32}}{{99 \cdot 19}}{{{\bar {\chi }}}^{3}} + \; \ldots } \right),\quad \bar {\chi } = \frac{{{{U}_{0}}}}{{k\varphi }}. \\ \end{gathered} $
Тогда в [1] получена система кинетических уравнений, описывающих поведение $f(t,\vec {x},\vec {\xi })$, функции распределения молекул в свободном состоянии, функций $s(t,\vec {x},\vec {w})$ и $h(t,\vec {x},\vec {w})$. Она имеет следующий вид:

$\frac{{\partial f}}{{\partial t}} + {{\xi }^{i}}\frac{{\partial f}}{{\partial {{x}^{i}}}} = {{J}_{B}} - {{J}_{{{\text{del}}}}},\quad \frac{{Ds}}{{Dt}} = {{J}_{{{\text{del}}}}},\quad \frac{{Dh}}{{Dt}} = {{A}_{d}},\quad {{J}_{B}} = {{d}^{2}}\int {d{{{\vec {\xi }}}_{1}}} \left( {\int\limits_0^\infty {\int\limits_0^{2\pi } g } (f{\kern 1pt} {\text{'}}f_{1}^{'} - f{{f}_{1}})bdbd\theta } \right),$
(1.1)
$\frac{D}{{Dt}} = \frac{\partial }{{\partial t}} + {{w}^{i}}\frac{\partial }{{\partial {{x}^{i}}}} + {{a}^{j}}\frac{\partial }{{\partial {{w}^{j}}}},\quad {{J}_{{{\text{del}}}}} = 2{{\chi }_{1}}{{d}^{2}}\int\limits_{g \leqslant 2{{\chi }_{1}}} {d\vec {g}\left( {\int\limits_0^{2\pi } {\int\limits_0^\pi {\left| {\cos \omega } \right|({{f}^{2}} - {{f}_{0}})\sin \omega d\omega d\theta } } } \right)} ,$
${{A}_{d}} = - 2{{\chi }_{1}}{{d}^{2}}\int\limits_{g \leqslant 2{{\chi }_{1}}} {\frac{m}{4}{{g}^{2}}d\vec {g}\int\limits_0^{2\pi } {\int\limits_0^\pi {\left| {\cos \omega } \right|({{f}^{2}} - {{f}_{0}})\sin \omega d\omega d\theta } } } ,\quad \vec {a} = \frac{{4\pi {{U}_{0}}{{d}^{3}}\operatorname{grad} (n)}}{{3m}}.$

В [1] вводятся следующие макропараметры от фигурирующих в (1.1) функций распределения:

${{n}_{f}} = \int {fd\vec {\xi },\quad {{n}_{d}} = \int {sd\vec {w}} } ,$
являющиеся соответственно числовыми плотностями молекул в свободном и связанном состояниях. Тогда $n = {{n}_{f}} + {{n}_{d}}$ есть числовая плотность газа. Определяется макроскопическая скорость неидеального газа:
$\vec {u} = \frac{1}{n}\left( {\int {\vec {\xi }fd\vec {\xi }} + \int {\vec {w}sd\vec {w}} } \right).$
Вводя аналоги тепловых скоростей в виде ${{\vec {c}}_{f}} = \vec {\xi } - \vec {u}$, ${{\vec {c}}_{d}} = \vec {w} - \vec {u}$, определяется тензор напряжений в неидеальном газе:
${{P}^{{ij}}} = m\int {c_{f}^{i}c_{f}^{j}fd\vec {\xi }} + m\int {c_{d}^{i}c_{d}^{j}sd\vec {w}} $
и давление в нем
$3P = {{P}^{{11}}} + {{P}^{{22}}} + {{P}^{{33}}} = m\left( {\int {\vec {c}_{v}^{2}fd\vec {\xi } + \int {\vec {c}_{d}^{2}sd\vec {w}} } } \right).$
В [1] отмечено, что такое определение давления решает поставленную в [2] проблему определения с точки зрения кинетической теории давления в неидеальном газе. Также вводятся следующие величины:
$\frac{3}{2}nkT = \int {m\frac{{\vec {c}_{f}^{2}}}{2}} fd\vec {\xi } + \int {\left( {m\frac{{\vec {c}_{d}^{2}}}{2}s - \frac{1}{2}h} \right)} d\vec {w},$
$\vec {Q} = \int {m{{{\vec {c}}}_{f}}\frac{{\vec {c}_{f}^{2}}}{2}} fd\vec {\xi } + \int {} {{\vec {c}}_{d}}\left( {m\frac{{\vec {c}_{d}^{2}}}{2}s - \frac{1}{2}h} \right)d\vec {w},$
являющиеся аналогом температуры неидеального газа и вектором теплового потока в нем.

В [1] из (1.1) были получены уравнения (сохранения), связывающие описанные выше макропараметры. Эти уравнения имеют следующий вид:

$\frac{{\partial n}}{{\partial t}} + \frac{{\partial (n{{u}_{i}})}}{{\partial {{x}_{i}}}} = 0,\quad i = 1,\;2,\;3,\quad \frac{{\partial (\rho {{u}^{j}})}}{{\partial t}} + \frac{{\partial (\rho {{u}^{i}}{{u}^{i}})}}{{\partial {{x}_{i}}}} + \frac{{\partial {{P}^{{ij}}}}}{{\partial {{x}_{i}}}} - {{n}_{d}}{{a}_{j}} = 0,\quad j = 1,\;2,\;3,$
$\frac{\partial }{{\partial t}}\left( {\frac{3}{2}nkT{\text{*}} + \frac{\partial }{{\partial t}}\left( {\frac{3}{2}nkT + \rho {{{\vec {u}}}^{2}}{\text{/}}2} \right)} \right) + \frac{\partial }{{\partial {{x}_{i}}}}\left( {{{u}^{i}}\left( {\frac{3}{2}nT + \rho {{{\vec {u}}}^{2}}{\text{/}}2} \right)} \right) + \frac{{\partial ({{u}_{k}}{{P}^{{ki}}})}}{{\partial {{x}_{i}}}} + \frac{{\partial {{Q}_{i}}}}{{\partial {{x}_{i}}}} - {{n}_{d}}u_{d}^{i}{{a}_{i}} = 0.$
Используя первое уравнение, три последующих уравнения можно привести к уравнениям вида

$mn\frac{{\partial {{u}^{j}}}}{{\partial t}} + mn{{u}^{i}}\frac{{\partial {{u}^{j}}}}{{\partial {{x}_{i}}}} + \frac{{\partial {{P}^{{ij}}}}}{{\partial {{x}_{i}}}} - {{n}_{d}}{{a}_{j}} = 0,\quad j = 1,\;2,\;3.$

Эти инвариантные относительно преобразования Галилея уравнения иногда называют уравнениями движения. Используя их и уравнение неразрывности, можно преобразовать последнее уравнение к инвариантному относительно преобразования Галилея виду

$\frac{3}{2}kn\left( {\frac{{\partial T}}{{\partial t}} + {{u}^{i}}\frac{{\partial T}}{{\partial {{x}_{i}}}}} \right) + {{P}^{{ki}}}\frac{{\partial {{u}_{k}}}}{{\partial {{x}_{i}}}} + \frac{{\partial {{Q}_{i}}}}{{\partial {{x}_{i}}}} - n(u_{d}^{i} - {{u}^{i}}){{a}_{i}} = 0.$

По своей сути это уравнение есть полученное из кинетической теории первое начало термодинамики для неидеального газа. А.А. Власов в одной из своих работ [3] заметил, что, с точки зрения кинетической теории, уравнение сохранения энергии получено только для идеального газа, так что последнее уравнение есть решение проблемы, поставленной А.А. Власовым.

Из приведенного выше определения макропараметров, характеризующих неидеальный газ, можно получить, что

${{\left( {P - \frac{1}{3}\int h d\vec {w}} \right)} \mathord{\left/ {\vphantom {{\left( {P - \frac{1}{3}\int h d\vec {w}} \right)} n}} \right. \kern-0em} n} = kT,$
т.е. уравнение состояния неидеального газа. Если отождествить $1{\text{/}}n$ с объемом, занятым газом, а $\frac{1}{3}\int {hd\vec {w}} $ с потенциальной энергией взаимодействия молекул, то последнее уравнение есть уравнение состояния, по форме совпадающее с уравнением состояния Ван-Дер-Ваальса (см. [2]).

Используя это обстоятельство, в [4] удалось получить выражение критической температуры для газа с рассматриваемым потенциалом взаимодействия. Получилось, что ${{T}_{k}} = \frac{8}{{135}}{{U}_{0}}$.

2. ПРИВЕДЕНИЕ СИСТЕМЫ К БЕЗРАЗМЕРНОМУ ВИДУ

В выражениях для ${{J}_{{{\text{del}}}}}$ и ${{A}_{d}}$ можно выполнить интегрирование. В случае использования потенциала межмолекулярного взаимодействия Сазерленда, который используется в данной работе, это делается особенно просто. Тогда

(2.1)
${{J}_{{{\text{del}}}}} = \frac{{128}}{3}{{\pi }^{2}}{{d}^{2}}\chi _{1}^{4}({{f}^{2}} - {{f}_{0}}),\quad {{A}_{d}} = - \frac{{128}}{5}{{\pi }^{2}}{{d}^{2}}{{U}_{0}}\chi _{1}^{4}({{f}^{2}} - {{f}_{0}}).$

Комбинируя второе и третье уравнение (1.1), получаем

$\frac{{D\left( {\frac{3}{5}{{U}_{0}}s + h} \right)}}{{Dt}} = 0.$
Отсюда $\frac{3}{5}{{U}_{0}}s + h = C$. Ясно, что если нет молекул в связанном состоянии, то нет и их энергии, поэтому имеем

(2.2)
$\frac{3}{5}{{U}_{0}}s(t,\vec {x},\vec {w}) = - h(t,\vec {x},\vec {w}).$

Из (2.2) следует, что

(2.3)
$\frac{{{{Z}_{h}}(\bar {\chi })}}{{{{Z}_{s}}(\bar {\chi })}} = - \frac{3}{5}.$
Уравнение (2.3) есть уравнение для определения $\varphi $. Из (2.3) следует, что введенная ранее величина $\varphi $ есть константа. Если $\varphi $ придать физический смысл, то она есть энергия, приходящаяся на одну молекулу в связанном состоянии. Тогда получившийся результат объясняется тем, что движение молекул в связанном состоянии происходит в $d$-масштабе, и их равновесное состояние достигается за временной промежуток, существенно меньший временного масштаба ${{t}_{\lambda }}$. Поэтому их энергия, приходящаяся на одну частицу, определяется только потенциалом их взаимодействия между собой, что является аргументом в пользу соотношения (2.2) (не исключено, что соотношение, определяющее $\varphi $, должно находиться при помощи квантовой механики).

Выше было показано, что фигурирующая в определении величина $\varphi (t,\vec {x},\vec {w})$ есть константа. Это позволяет установить, что для неидеального газа, который описывается системой (1.1), имеет место Н-теорема.

Определим Н-функцию молекул в связанном состоянии в виде

${{H}_{d}} = \int {S(t,\vec {w})d\vec {w}} ,$
где

$S(t,\vec {w}) = \int\limits_{{{\Omega }_{{{\text{del}}}}}} {{{{\bar {F}}}_{d}}ln{{{\bar {F}}}_{d}}d\vec {r}d\vec {g} = } s\ln {{f}_{0}} + \frac{h}{{k\varphi }} = s\left( {\ln {{f}_{0}} - \frac{{3{{U}_{0}}}}{{5k\varphi }}} \right).$

Это выражение получено, если учесть (2.2). Принимая во внимание (2.1), можно получить, что

(2.4)
${{\dot {H}}_{d}} = \frac{k}{2}\int {\frac{{\partial S}}{{\partial t}}d\vec {w} = \frac{k}{2}} \int {\left( {\ln {{f}_{0}} + 1 - \frac{{3{{U}_{0}}}}{{5k\varphi }}} \right)\frac{{\partial s}}{{\partial t}}d\vec {w}} = r\frac{k}{2}\int {\left( {\ln {{f}_{0}} + 1 - \frac{{3{{U}_{0}}}}{{5k\varphi }}} \right)} ({{f}^{2}}(t,\vec {x},\vec {\xi }) - {{f}_{0}}(t,\vec {x},\vec {w}))d\vec {w},$
где $r = \frac{{128}}{3}{{d}^{2}}\chi _{1}^{4}$. Функцию молекул в свободном состоянии Нf введем в виде
${{H}_{f}} = k\int {f\ln \left( {f{{e}^{{ - \left( {\frac{1}{2} + \frac{{3{{U}_{0}}}}{{10k\varphi }}} \right)}}}} \right)d\vec {\xi }} .$
Тогда получим
${{\dot {H}}_{f}} = \frac{k}{2}\int {\frac{{\partial f}}{{\partial t}}\left( {1 - \frac{{3{{U}_{0}}}}{{5k\varphi }} + \ln {{f}^{2}}} \right)d\vec {\xi }} = {{\bar {J}}_{B}} - r\frac{k}{2}\int {\left( {\ln {{f}^{2}} + 1 - \frac{{3{{U}_{0}}}}{{5k\varphi }}} \right)} ({{f}^{2}}(t,\vec {x},\vec {\xi }) - {{f}_{0}}(t,\vec {x},\vec {w}))d\vec {w},$
где
${{\bar {J}}_{B}} = k\int\limits_D {d\vec {\xi }\iiint {\ln \frac{{f{{f}_{1}}}}{{f{\kern 1pt} {\text{'}}{{{f'}}_{1}}}}(f{\kern 1pt} {\text{'}}f_{1}^{'} - f{{f}_{1}})bdbd\varepsilon d\vec {\xi } \leqslant 0}} .$
Откуда имеем
${{\dot {H}}_{f}} + {{\dot {H}}_{d}} = {{\bar {J}}_{B}} + r\frac{k}{2}\int {\ln \frac{{{{f}_{0}}}}{{{{f}^{2}}}}} ({{f}^{2}}(t,\vec {x},\vec {\xi }) - {{f}_{0}}(t,\vec {x},\vec {w}))d\vec {w} \leqslant 0,$
что доказывает Н-теорему.

Приведем уравнения (1.1) к безразмерному виду. Пусть ${{T}_{0}}$ есть характерное значение температуры, а $L$ – масштабный размер длины рассматриваемой задачи. Как обычно, штрихи будут соответствовать безразмерным переменным. Тогда получим

${{x}_{i}} = Lx_{i}^{'},\quad {{\xi }_{i}} = {{\xi }_{0}}\xi _{i}^{'},\quad {{w}_{i}} = {{\xi }_{0}}w_{i}^{'},\quad t = \frac{L}{{{{\xi }_{0}}}}t{\text{'}},\quad {{\xi }_{0}} = \sqrt {\frac{{2k{{T}_{0}}}}{m}} .$
Полагая ${{n}_{0}}$ характерным значением плотности, получим

$f = \frac{{{{n}_{0}}}}{{\xi _{0}^{3}}}f{\kern 1pt} {\text{'}},\quad {{f}_{0}} = \frac{{n_{0}^{2}}}{{\xi _{0}^{6}}}f_{0}^{'},\quad s = s{\text{'}} \cdot {{s}_{0}} = 16{{\pi }^{2}}{{d}^{3}}{{Z}_{s}}(\bar {\chi })\chi _{1}^{{3/2}}\frac{{n_{0}^{2}}}{{\xi _{0}^{6}}}f_{0}^{'}.$

Отсюда следует, что

${{s}_{0}} = 16{{\pi }^{2}}{{d}^{3}}{{Z}_{s}}(\bar {\chi })\chi _{1}^{{3/2}}\frac{{n_{0}^{2}}}{{\xi _{0}^{6}}},\quad s{\text{'}} = f_{0}^{'}.$

Во введенных выше безразмерных переменных, штрихи у которых опущены, система (1.1) будет иметь следующий вид:

(2.5)
$\begin{gathered} \frac{{\partial f}}{{\partial t}} + {{\xi }^{i}}\frac{{\partial f}}{{\partial {{x}^{i}}}} = \frac{1}{{Kn}}\left( {{{J}_{B}} - \frac{{16{{\pi }^{2}}}}{3}{{\chi }^{2}}({{f}^{2}} - {{f}_{o}})} \right),\quad Kn = \frac{1}{{{{n}_{0}}{{d}^{2}}L}},\quad \chi = \frac{{{{U}_{o}}}}{{k{{T}_{0}}}}, \\ \frac{{\partial {{f}_{0}}}}{{\partial t}} + {{w}^{i}}\frac{{\partial {{f}_{0}}}}{{\partial {{x}^{i}}}} - \frac{{2\pi \varepsilon }}{{3Kn}}\chi \frac{{\partial n}}{{\partial {{x}_{j}}}}\frac{{\partial {{f}_{0}}}}{{\partial {{w}^{j}}}} = \frac{{4\sqrt 2 }}{{3\varepsilon Kn}}\sqrt \chi ({{f}^{2}} - {{f}_{0}}),\quad \varepsilon = {{n}_{0}}{{d}^{3}},\quad h = - {{f}_{0}}. \\ \end{gathered} $
Последнее соотношение (2.5) получается из (2.3), если характерное значение $h$ положить в виде ${{h}_{0}} = \frac{3}{5}{{U}_{0}}{{s}_{0}}$.

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

$n = \int {fd\vec {\xi } + 4\sqrt 2 {{\pi }^{2}}\varepsilon {{\chi }^{{3/2}}}\int {{{f}_{0}}d\vec {w}} } ,\quad n\vec {u} = \int {\vec {\xi }fd\vec {\xi } + 4\sqrt 2 {{\pi }^{2}}\varepsilon {{\chi }^{{3/2}}}\int {\vec {w}{{f}_{0}}d\vec {w}} } .$
Если принять за характерные значения тензора напряжений величину $\bar {P} = m{{n}_{0}}\xi _{0}^{2}$, температуры ${{T}_{0}}$, а теплового потока $\bar {Q} = \frac{1}{2}m{{n}_{0}}\xi _{0}^{3}$, то их безразмерные значения будут вычисляться по следующим формулам:

${{P}^{{ij}}} = \int {c_{f}^{i}c_{f}^{j}fd\vec {\xi }} + 4\sqrt 2 {{\pi }^{2}}\varepsilon {{\chi }^{{3/2}}}\int {c_{d}^{i}c_{d}^{j}{{f}_{0}}d\vec {w}} $,      $\frac{3}{2}nT = \int {\vec {c}_{f}^{2}} fd\vec {\xi } + \int {\left( {\vec {c}_{d}^{2}{{f}_{0}} - \frac{3}{{10}}\chi {{f}_{0}}} \right)} d\vec {w}$, $\vec {Q} = \int {{{{\vec {c}}}_{f}}\vec {c}_{f}^{2}} fd\vec {\xi } + $ $ + \;\int {{{{\vec {c}}}_{d}}} \left( {\vec {c}_{d}^{2}{{f}_{0}} - \frac{3}{{10}}\chi {{f}_{0}}} \right)d\vec {w}$.

Появившийся (2.5) параметр $\varepsilon $ есть внепорядковый член асимптотического разложения, с помощью которого была получена система (1.1), см. [1]. Однако следует заметить, что он присутствует только в комбинации с параметром $\chi $. Из второго уравнения (2.5) можно получить, что

$\frac{{\sqrt \chi }}{{Kn}}\sqrt \chi ({{f}^{2}} - {{f}_{0}}) = \frac{3}{{4\sqrt 2 }}\varepsilon \frac{{D{{f}_{0}}}}{{Dt}}.$

Тогда правая часть первого уравнения (2.5) будет иметь вид

$\frac{{{{J}_{B}}}}{{Kn}} - \frac{{4{{\pi }^{2}}}}{{\sqrt 2 }}\varepsilon {{\chi }^{{3/2}}}\frac{{D{{f}_{0}}}}{{Dt}}.$
Отсюда неидеальность газа наступает при $\varepsilon {{\chi }^{{3/2}}} \approx 1$. В [1] показано, что, если вместо длины свободного пробега взять величину $l = n_{0}^{{ - 1/3}}$, то, как это показано в [1], интеграл столкновений будет иметь порядок ${{\varepsilon }^{{2/3}}}$. Тогда неидеальность будет наступать при ${{\varepsilon }^{{1/3}}}{{\chi }^{{3/2}}} \approx 1$, т.е. при существенно меньших значениях параметра $\chi $.

3. ПОСТРОЕНИЕ ЧИСЛЕННОГО МЕТОДА

При построении численной схемы решения системы (2.1) предполагается использовать метод расщепления по физическим параметрам. Этот метод был предложен в [5]. В [6] были предложены разные модификации этого метода, а также показано, что он обладает вторым порядком точности по временному шагу $\Delta t$. Применительно к системе (2.1) можно предложить следующую схему метода расщепления.

Пусть при $t = {{t}_{i}}$ известны $f({{t}_{i}},\vec {x},\vec {\xi }) = {{f}_{i}}$, ${{f}_{0}}({{t}_{i}},\vec {x},\vec {w}) = f_{0}^{i}$ и $n({{t}_{i}},\vec {x}) = {{n}_{i}}$. На промежутке $\left[ {{{t}_{i}},{{t}_{i}} + \frac{{\Delta t}}{2}} \right]$ предполагается решать следующую задачу:

(3.1)
$\begin{gathered} \frac{{\partial f}}{{\partial t}} = \frac{1}{{Kn}}\left( {{{J}_{B}} - \frac{{16{{\pi }^{2}}}}{3}{{\chi }^{2}}({{f}^{2}} - {{f}_{0}})} \right),\quad f({{t}_{i}}\vec {x},\vec {\xi }) = {{f}_{i}},\quad {{\alpha }_{1}} = \frac{{2\pi \varepsilon }}{{3Kn}}\chi , \\ \frac{{\partial {{f}_{0}}}}{{\partial t}} - {{\alpha }_{1}}\frac{{\partial {{n}_{i}}}}{{\partial {{x}_{j}}}}\frac{{\partial {{f}_{0}}}}{{\partial {{w}^{j}}}} = \frac{{{{\alpha }_{2}}}}{\varepsilon }({{f}^{2}} - {{f}_{0}}),\quad {{f}_{0}}({{t}_{i}},\vec {x},\vec {w}) = f_{0}^{i},\quad {{\alpha }_{2}} = \frac{{4\sqrt 2 }}{{3Kn}}. \\ \end{gathered} $

Решение второго уравнения (3.1) можно записать в виде

$\begin{gathered} {{f}_{0}}(t,\vec {x},\vec {w}) = f_{0}^{i}(\vec {x},\vec {w} + {{\alpha }_{1}}\operatorname{grad} {{n}_{i}}(t - {{t}_{i}})){\kern 1pt} {\text{exp}}\left\{ { - \frac{\alpha }{\varepsilon }(t - {{t}_{i}})} \right\} + \frac{\alpha }{\varepsilon }\int\limits_{{{t}_{i}}}^t {f_{i}^{2}(s,\vec {x},\vec {w} + } \\ + \;{{\alpha }_{1}}\operatorname{grad} {{n}_{i}}(t - s)){\text{exp}}\left\{ { - \frac{\alpha }{\varepsilon }(t - s)} \right\}ds,\quad t \in \left[ {{{t}_{i}},{{t}_{i}} + \frac{{\Delta t}}{2}} \right]. \\ \end{gathered} $

Используя схему первого порядка точности, на первом этапе метода расщепления получаем

(3.2)
$\begin{gathered} f_{0}^{{i + 1/2}} - f_{i}^{2} = \left( {f_{o}^{i}(\vec {x},\vec {w} + {{\alpha }_{1}}\operatorname{grad} {{n}_{i}}\frac{{\Delta t}}{2}} \right) - f_{i}^{2}(\vec {x},\vec {w})){\kern 1pt} {\text{exp}}\left\{ { - \frac{{\alpha \Delta t}}{{2\varepsilon }}} \right\}, \\ {{f}_{{i + \frac{1}{2}}}} = {{J}_{B}}({{f}_{i}})\frac{{\Delta t}}{{2Kn}} + \frac{{16{{\pi }^{2}}}}{3}{{\chi }^{2}}\left( {f_{0}^{{i + 1/2}} - f_{i}^{2}} \right). \\ \end{gathered} $
Следующий этап метода расщепления – это разлет, т.е. на $\left[ {{{t}_{i}},{{t}_{i}} + \Delta t} \right]$ решаем систему:
$\frac{{\partial f}}{{\partial t}} + {{\xi }^{i}}\frac{{\partial f}}{{\partial {{x}^{i}}}} = 0,\quad \frac{{\partial {{f}_{0}}}}{{\partial t}} + {{w}^{i}}\frac{{\partial {{f}_{0}}}}{{\partial {{x}^{i}}}} = 0$
с начальными условиями $f({{t}_{i}},\vec {x},\vec {\xi }) = {{f}_{{i + 1/2}}}$, ${{f}_{0}}({{t}_{i}},\vec {x},\vec {w}) = f_{0}^{{i + 1/2}}$. Тогда получим
${{\tilde {f}}_{{i + 1/2}}} = {{f}_{{i + 1/2}}}(\vec {x} - \vec {\xi }\Delta t),\quad \tilde {f}_{0}^{{i + 1/2}} = f_{0}^{{i + 1/2}}(\vec {x} - \vec {w}\Delta t).$
Вычислив по функциям ${{\tilde {f}}_{{i + 1/2}}}$ и $\tilde {f}_{0}^{{i + 1/2}}$ значение плотности ${{n}_{{i + 1/2}}}$, на $\left[ {{{t}_{i}} + \frac{{\Delta t}}{2},\;{{t}_{i}} + \Delta t} \right]$ осуществим снова этап релаксации с начальными функциями распределения ${{\tilde {f}}_{{i + 1/2}}}$, $\tilde {f}_{0}^{{i + 1/2}}$.

Выше был описан вариант схемы метода расщепления, который может быть применен для численного решения системы (3.1). В отличие от случаев, рассмотренных в [5] и [6], система (3.1) содержит силовой член. Его действие в предложенной схеме, как и в [7], учитывается на этапе релаксации. Если следовать [6], то при построении метода расщепления в зависимости от специфики задачи можно менять местами этапы релаксации и свободно молекулярного разлета. Из анализа значений безразмерных величин в (3.1) следует, что $\frac{\chi }{\varepsilon } \gg 1$. Тогда ${{f}_{0}}$, как это следует из (3.2), будет быстро релаксировать к ${{f}^{2}}$, поэтому правую часть в первом уравнении (3.1) можно считать зависящей от предыдущего момента времени. В зависимости от величин параметров ${\text{Kn}}$, $\chi $, $\varepsilon $ в рамках метода расщепления по физическим параметрам возможны различные вариации.

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

(3.3)
${{J}_{B}} = \int {\int\limits_0^\infty {\int\limits_0^{2\pi } {g(f{\text{'}}f_{1}^{'} - f{{f}_{1}})bdbd\theta d{{{\vec {\xi }}}_{1}}} } } .$
Этот интеграл есть вклад в изменение функции распределения молекул в свободном состоянии, который возникает в результате взаимодействия (столкновения) молекул, двигающихся со скоростями $\vec {\xi }$ и ${{\vec {\xi }}_{1}}$ соответственно. Так как сама система уравнений (1.1) выведена в предположении, что $\frac{d}{\lambda } \ll 1$, считается, что молекулы взаимодействуют мгновенно, приобретая при этом скорости $\vec {\xi }{\text{'}}$ и $\vec {\xi }_{1}^{'}$. В (3.3) имеем
$f = f(t,\vec {\xi }),\quad {{f}_{1}} = f(t,{{\vec {\xi }}_{1}}),\quad f{\kern 1pt} {\text{'}} = f(t,\vec {\xi }{\text{'}}),\quad f_{1}^{'} = f(t,\vec {\xi }_{1}^{'}).$
Конкретное выражение скоростей молекул после столкновения получается при решении задачи о взаимодействии молекул, имеющих до столкновения скорости $\vec {\xi }$ и ${{\vec {\xi }}_{1}}$. Подробное решение этой задачи для молекул, которые отталкиваются как точечные центры, приведено в [8]. Там получено, что
(3.4)
$\vec {\xi }{\text{'}} = \vec {\xi } + (\vec {g} \times \vec {k})\vec {k},\quad \vec {\xi }_{1}^{'} = {{\vec {\xi }}_{1}} - (\vec {g} \times \vec {k})\vec {k},$
где единичный вектор $\vec {k} = \vec {k}(b,g,\varepsilon )$ определяется в ходе решения задачи о столкновении этих молекул. Рассмотрим решение этой задачи для введенного ранее потенциала притяжения молекул.

Будем называть молекулы, двигающиеся со скоростью $\vec {\xi }$, $\xi $-молекулы. Удобно (см. [8]) решение поставленной выше задачи проводить в системе координат, связанной с $\xi $-молекулой и двигающейся со скоростью

$\vec {w} = \frac{{\vec {\xi } + {{{\vec {\xi }}}_{1}}}}{2} = \frac{{\vec {\xi }{\text{'}} + \vec {\xi }_{1}^{'}}}{2},$
которая есть скорость центра масс $\xi $ и $\xi $1 молекул. Так как взаимодействие молекул происходит в $d$-масштабе, то задача сводится к определению движения материальной точки, массы $\mu = \frac{m}{2}$ и имеющей на бесконечно удаленном расстоянии от центра скорость $\vec {g} = {{\vec {\xi }}_{1}} - \vec {\xi }$ в центральном поле с потенциалом
$U(r) = \left[ \begin{gathered} + \infty ,\quad 0 \leqslant r \leqslant d, \hfill \\ - {{U}_{0}}{{\left( {\frac{d}{r}} \right)}^{4}},\quad r > d. \hfill \\ \end{gathered} \right.$
Из курса механики известно, что при движении материальной точки в центральном поле сохраняется ее момент количества движения $\vec {L} = \mu [\vec {r} \times \vec {g}]$, где $\vec {r}$и $\vec {g}$ суть радиус-вектор точки и ее скорость соответственно. Из последнего соотношения следует, что движение рассматриваемой точки будет все время происходить в плоскости, перпендикулярной вектору $\vec {L}$ (плоскость Р на фиг. 1). Поэтому интеграл столкновений в уравнении Больцмана вычисляется в цилиндрической системе координат, в которой ось $z$ выбирается коллинеарной вектору $\vec {g}$, а значение полярного угла $\varepsilon $ определяет плоскость, где происходит движение.

Фиг. 1.

В плоскости Р вводится полярная система координат, где за полярную ось принимается ось $z$. Пусть $\vec {\tilde {g}} = \{ {{\tilde {g}}_{r}},{{\tilde {g}}_{\theta }}\} $ – скорость материальной точки, двигающейся относительно силового центра О, а $\;r$ и $\theta $ – ее полярные радиус и угол.

Из условий сохранения величины $\left| {\vec {L}} \right|$ и полной энергии имеем

(3.5)
$\tilde {g}r\sin \theta = {{\tilde {g}}_{\theta }}r = gb = g{\text{'}}b{\text{'}},\quad \tilde {g}_{r}^{2} + \tilde {g}_{\theta }^{2} - 2\frac{\chi }{{{{r}^{4}}}} = {{g}^{2}} = {{(g{\text{'}})}^{2}}.$

В (3.5) значения $b$, $b{\text{'}}$ – прицельные параметры – минимальное расстояние между центром О и движущейся материальной точкой, какое было бы, если бы точка двигалась с постоянной скоростью $g$ или $g{\text{'}}$. Соответственно $\vec {g}$ и $\vec {g}{\text{'}}$ – скорости материальной точки до и после “столкновения” ($r = \infty $). Изображенный на фиг. 2 единичный вектор $\vec {k}$ есть вектор, коллинеарный радиус-вектору, соединяющему центр О с двигающейся относительно него материальной точкой, когда она находится на минимальном от него расстоянии ${{r}_{m}}$. Из (3.5) следует, что $g = g{\text{'}}$, $b = b{\text{'}}$. Таким образом, движение точки относительно центра О будет симметрично относительно оси, которая задается вектором $\vec {k}$, поэтому будем иметь $\vec {g}{\text{'}} = \vec {g} + 2(\vec {g} \times \vec {k})\vec {k}$. Из последней формулы и получается (3.4) – выражение для скоростей молекул после столкновения. Сам же вектор $\vec {k}$ имеет следующие компоненты:

$\vec {k} = \{ \sin \theta \cos \varepsilon ,\;\sin \theta \sin \varepsilon ,\;\cos \theta \} $.
Фиг. 2.

В случае межмолекулярного потенциала взаимодействия, который используется в данной работе, ${{r}_{m}}$ – минимальное расстояние между силовым центром и движущейся точкой в случае, когда ${{b}^{4}} > \frac{{8\chi }}{{{{g}^{2}}}}$, и определяется формулой

$r_{m}^{2} = \frac{{{{b}^{2}} + \sqrt {{{b}^{4}} - 8\frac{\chi }{{{{g}^{2}}}}} }}{2}.$

Этот случай соответствует фиг. 2а. Если ${{b}^{4}} < \frac{{8\chi }}{{{{g}^{2}}}}$, то движущаяся точка достигает границы сферы с радиусом 1 и упруго отражается от нее (см. фиг. 2б). Движение материальной точки остается симметричным относительно вектора $\vec {k}$ и в этом случае. Поэтому имеем

(3.6)
$\theta = \int\limits_0^{1/{{r}_{m}}} {\frac{{bdu}}{{\sqrt {\frac{{2\chi }}{{{{g}^{2}}}}{{u}^{4}} - {{b}^{2}}{{u}^{2}} + 1} }}} ,\quad {{r}_{m}} = \left[ \begin{gathered} \sqrt {\frac{{{{b}^{2}} + \sqrt {{{b}^{4}} - 8\frac{\chi }{{{{g}^{2}}}}} }}{2}} ,\quad {{b}^{4}} - 8\frac{\chi }{{{{g}^{2}}}} > 0, \hfill \\ 1,\quad {{b}^{4}} - 8\frac{\chi }{{{{g}^{2}}}} < 0. \hfill \\ \end{gathered} \right.$

Заметим, что при ${{b}^{4}} - 8\frac{\chi }{{{{g}^{2}}}} = 0$ равенство (3.6) будет следующим:

$\theta = \int\limits_0^{\sqrt 2 /b} {\frac{{bdu}}{{{{{\left( {\frac{{{{b}^{2}}{{u}^{2}}}}{2}} \right)}}^{2}} - 1}}} .$
Этот интеграл расходится в $\sqrt 2 {\text{/}}b$. Это соответствует тому, что движущаяся точка будет вращаться вокруг силового центра по круговой орбите с радиусом, равным 1. В [9] указывалось, что замкнутые траектории в случае межмолекулярного потенциала притяжения могут быть и при положительном значении полной энергии взаимодействующих частиц. Там же говорилось, что эти траектории неустойчивые. Это становится понятно из приведенного выше анализа. Действительно, если ${{b}^{4}} - 8\frac{\chi }{{{{g}^{2}}}} \approx 0$, то в случае положительного значения этой величины картинку (см. фиг. 2а), причем величина $\theta $ может быть больше, чем $2\pi $. В случае же отрицательного значения ${{b}^{4}} - 8\frac{\chi }{{{{g}^{2}}}}$ движение будет происходить в соответствии с фиг. 2б, т.е. движущаяся точка отразится от сферы и уйдет на бесконечность.

На этапе релаксации всегда будет решаться первое уравнение (3.1), куда входит пятикратный интеграл столкновений ${{J}_{B}}$. Как уже было отмечено выше, его вычисление является наиболее сложной задачей при построении численной схемы решения рассматриваемых кинетических уравнений. В данной статье предлагается следующая схема вычисления интеграла столкновений.

Сам интеграл ${{J}_{B}}$ можно переписать в виде

(3.7)
${{J}_{B}} = \int {I(\vec {\xi },{{{\vec {\xi }}}_{1}})d{{{\vec {\xi }}}_{1}}} ,$
где

$I(\vec {\xi },{{\vec {\xi }}_{1}}) = \int\limits_0^\infty {\int\limits_{}^{2\pi } {(f{\kern 1pt} {\text{'}}f{\kern 1pt} {\text{'}} - f{{f}_{1}})g} } bdbd\varepsilon .$

Интеграл (3.7) и все макропараметры от функций распределения $f$ и $s$ будут вычисляться с помощью квадратурных формул типа Гаусса (см. [10]). В частности,

$\int\limits_{ = \infty }^{ + \infty } {{{e}^{{ - {{x}^{2}}}}}G(x)dx \cong \sum\limits_{k = 1}^n {{{C}_{k}}G({{x}_{k}})} } ,$
где ${{x}_{k}}$ суть нули полинома Эрмита ${{H}_{n}}(x)$.

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

Тогда если $\psi (\vec {\xi }) = \psi ({{\xi }_{1}},{{\xi }_{2}},{{\xi }_{3}})$ есть функция от скоростей молекулы, то имеем

$\int {\psi (\vec {\xi })f(t,\vec {\xi })d\vec {\xi } \cong \sum\limits_{{{k}_{1}} = 1}^n {\sum\limits_{{{k}_{2}} = 1}^n {\sum\limits_{{{k}_{3}} = 1}^n {{{C}_{{{{k}_{1}}}}}{{C}_{{{{k}_{2}}}}}{{C}_{{{{k}_{3}}}}}} } } } \psi \left( {{{\xi }_{{{{k}_{1}}}}},{{\xi }_{{{{k}_{2}}}}},{{\xi }_{{{{k}_{3}}}}}} \right)\varphi \left( {t,{{\xi }_{{{{k}_{1}}}}},{{\xi }_{{{{k}_{2}}}}},{{\xi }_{{{{k}_{3}}}}}} \right) = \sum\limits_{\vec {k} = 1}^n {{{{\vec {C}}}_{k}}\psi ({{{\vec {\xi }}}_{k}})\varphi (t,{{{\vec {\xi }}}_{k}})} ,$
где $\varphi (t,\vec {\xi })$ определяется соотношением

$f(t,\vec {\xi }) = {{e}^{{ - {{{\vec {v}}}^{2}}}}}\varphi (t,\vec {\xi }),\quad {{\vec {v}}^{2}} = \xi _{1}^{2} + \xi _{2}^{2} + \xi _{3}^{2}.$

Подставив его в (3.1), получим уравнение, которое будет определять $\varphi (t,\vec {\xi })$.

Численный метод решения системы (3.1) состоит в следующем. В скоростном пространстве фиксируются точки ${{\vec {\xi }}_{k}} = \{ {{\xi }_{{{{k}_{1}}}}},{{\xi }_{{{{k}_{2}}}}},{{\xi }_{{{{k}_{3}}}}}\} $, ${{k}_{i}} = 1,\; \ldots ,\;n$, $i = 1,\;2,\;3$, являющиеся нулями полинома Эрмита ${{H}_{n}}(x)$. Из (3.4) следует, что ${{(\operatorname{v} _{1}^{'})}^{2}} + {{(\operatorname{v} ')}^{2}} = {\text{v}}_{1}^{2} + {{{\text{v}}}^{2}}$. Тогда получим, что

$\frac{{\partial \varphi (t,{{{\vec {\xi }}}_{k}})}}{{\partial t}} = \int {{{I}_{k}}} \left| {{{{\vec {\xi }}}_{1}} - \vec {\xi }} \right|d{{\vec {\xi }}_{{_{1}}}} + {{Q}_{k}},$
где
${{I}_{k}} = I({{\vec {\xi }}_{k}},{{\vec {\xi }}_{1}}) = \int\limits_0^\infty {\int\limits_0^{2\pi } {{{e}^{{ - \operatorname{v} _{1}^{2}}}}(\varphi _{1}^{'}\varphi {\text{'}} - {{\varphi }_{1}}{{\varphi }_{k}})\left| {{{{\vec {\xi }}}_{1}} - {{{\vec {\xi }}}_{k}}} \right|bdbd\varepsilon } } ,\quad {{\varphi }_{k}} = \varphi (t,{{\vec {\xi }}_{k}}),$
а
${{Q}_{k}} = \frac{{16{{\pi }^{2}}}}{{3Kn}}{{\chi }^{2}}({{e}^{{ - {{{v}}^{2}}}}}\varphi _{k}^{2} - {{f}_{0}}),$
где ${{f}_{0}}$ определяется при решении второго уравнения (3.1). Положим, что
(3.8)
$\int {{{I}_{k}}({{{\vec {\xi }}}_{k}},{{{\vec {\xi }}}_{1}})\left| {{{{\vec {\xi }}}_{1}} - \vec {\xi }} \right|d{{{\vec {\xi }}}_{{_{1}}}}} \cong \sum\limits_{\vec {l}}^n {{{{\vec {C}}}_{{\vec {l}}}}\int\limits_0^\infty {\int\limits_0^{2\pi } {\varphi _{{1lk}}^{'}\varphi _{{lk}}^{'} - {{\varphi }_{l}}{{\varphi }_{k}}){{g}_{{lk}}}bdbd\varepsilon } } } ,$
где $\varphi _{{1lk}}^{'} = \varphi (t,\vec {\xi }_{{1lk}}^{'})$, $\varphi _{{lk}}^{'} = \varphi (t,{{\vec {\xi '}}_{{lk}}})$, ${{\varphi }_{l}} = \varphi (t,{{\vec {\xi }}_{l}})$, ${{\vec {g}}_{{lk}}} = {{\vec {\xi }}_{l}} - {{\vec {\xi }}_{k}}$, ${{g}_{{lk}}} = \left| {{{{\vec {g}}}_{{lk}}}} \right|$, $\vec {\xi }_{{1/k}}^{'} = {{\vec {\xi }}_{l}} - ({{\vec {g}}_{{lk}}} \cdot \vec {k})\vec {k}$, $\vec {\xi }_{{lk}}^{'} = {{\vec {\xi }}_{k}} + ({{\vec {g}}_{{lk}}} \cdot \vec {k})\vec {k}$, $l = 1,\;2,\; \ldots ,\;n$.

Заменив правую часть уравнения для $\varphi (t,{{\vec {\xi }}_{k}})$ представлением (3.8), получим дифференциальное уравнение для определения функции $\varphi (t,{{\vec {\xi }}_{k}})$. Полученный при его решении набор функций $\varphi (t,{{\vec {\xi }}_{k}})$, $k = 1,\;2, \ldots ,\;n$, есть решение поставленной задачи.

Известно (см. [8]), что интеграл столкновений Больцмана ${{J}_{B}}(\vec {\xi })$ обладает тем свойством, что

(3.9)
$\int {\psi (\vec {\xi }){{J}_{B}}(\vec {\xi })d\vec {\xi } = \int {\int {I(\vec {\xi },{{{\vec {\xi }}}_{1}})\psi (\vec {\xi })d{{{\vec {\xi }}}_{1}}d\vec {\xi }} } } = - \frac{1}{4}\int {\int {I(\vec {\xi })g(\psi (\vec {\xi }_{1}^{'} + \psi (\vec {\xi }{\text{'}}) - \psi ({{{\vec {\xi }}}_{1}}) - \psi (\vec {\xi }))} } d{{\vec {\xi }}_{1}}d\vec {\xi }.$

В работе [5] была предложена численная схема вычисления интеграла столкновений Больцмана, которая обладала свойством (3.9) для $\psi (\vec {\xi }) = 1$, ${{\xi }_{i}}$, ${{\vec {\xi }}^{2}}$, $i = 1,\;2,\;3$. Численные схемы решения уравнения Больцмана, которые обладали таким свойством, автор [5] назвал консервативными. Ясно, что если первые тринадцать моментов от функции распределения вычислять при помощи таких схем, то разностные численные схемы, аппроксимирующие моментные уравнения для этих моментов, будут удовлетворять уравнениям сохранения.

Надо отметить, что в настоящее время консервативные схемы широко используются при решении кинетических уравнений. Так, в работе [11] консервативная численная схема решения была получена для $S$-модельного кинетического уравнения.

Покажем, что предложенный выше метод вычисления ${{J}_{B}}$ приводит к консервативной численной схеме. Рассмотрим равенство

$\int {\psi (\vec {\xi }){{J}_{B}}(\vec {\xi })d\vec {\xi } = } \int {\int {\psi (\vec {\xi })I(\vec {\xi },{{{\vec {\xi }}}_{1}})d{{{\vec {\xi }}}_{1}}d\vec {\xi }} } ,$
где

$I(\vec {\xi },{{\vec {\xi }}_{1}}) = \int\limits_0^\infty {\int\limits_0^{2\pi } {{{e}^{{ - {{\operatorname{v} }^{2}} - \operatorname{v} _{1}^{2}}}}(\varphi _{1}^{'}\varphi {\text{'}} - {{\varphi }_{1}}\varphi )gbdbd\varepsilon } } .$

Сделаем в интеграле замену переменных интегрирования $\vec {\xi } \to {{\vec {\xi }}_{1}}$, ${{\vec {\xi }}_{1}} \to \vec {\xi }$. Так как $I(\vec {\xi },{{\vec {\xi }}_{1}}) = I({{\vec {\xi }}_{1}},\vec {\xi })$, то получим, что

(3.10)
$Q = \int {\psi (\vec {\xi }){{J}_{B}}(\vec {\xi })d\vec {\xi } = } \frac{1}{2}\int {\int {I(\vec {\xi },{{{\vec {\xi }}}_{1}})(\psi (\vec {\xi }) + \psi ({{{\vec {\xi }}}_{1}}))d{{{\vec {\xi }}}_{1}}d\vec {\xi }} } .$

Схема вычисления $Q$ будет следующей:

$Q \cong \sum\limits_{\vec {k}}^n {\sum\limits_{\vec {l}}^n {{{{\vec {C}}}_{{\vec {k}}}}{{{\vec {C}}}_{{\vec {1}}}}\psi ({{{\vec {\xi }}}_{k}})} } \bar {I}({{\vec {\xi }}_{k}},{{\vec {\xi }}_{l}}),\quad {\text{где}}\quad \bar {I}({{\vec {\xi }}_{k}},{{\vec {\xi }}_{l}}) = \int\limits_0^\infty {\int\limits_0^{2\pi } {(\varphi _{{1lk}}^{'}\varphi {{{_{{lk}}^{'}}}_{{}}} - {{\varphi }_{l}}{{\varphi }_{k}}){{g}_{{lk}}}bdbd\varepsilon } } .$
Меняя в последней формуле индекс $\vec {k}$ на $\vec {l}$, а $\vec {l}$на $\vec {k}$, получаем разностный аналог написанной выше формулы:
$Q \cong \frac{1}{2}\sum\limits_{\vec {l}}^n {\sum\limits_{\vec {k}}^n {{{{\vec {C}}}_{{\vec {k}}}}{{{\vec {C}}}_{{\vec {1}}}}(\psi ({{{\vec {\xi }}}_{k}}) + \psi ({{{\vec {\xi }}}_{l}}))} } \bar {I}({{\vec {\xi }}_{k}},{{\vec {\xi }}_{l}}).$
В интеграле (3.10) сделаем замену переменных интегрирования $\vec {\xi } = \vec {\xi }{\text{'}}$, ${{\vec {\xi }}_{1}} = \vec {\xi }_{1}^{'}$. Так как $\vec {g} = - \vec {g}{\text{'}} = \vec {\xi }_{1}^{'} - \vec {\xi }{\text{'}}$, а $\vec {\xi } = \vec {\xi }{\text{'}} - (\vec {g} \cdot \vec {k})\vec {k} = \vec {\xi }{\text{'}} + (\vec {g}{\text{'}} \cdot \vec {k})\vec {k}$, ${{\vec {\xi }}_{1}} = \vec {\xi }_{1}^{'} + (\vec {g} \cdot \vec {k})\vec {k} = \vec {\xi }_{1}^{'} - (\vec {g}{\text{'}} \cdot \vec {k})\vec {k}$, и $d\vec {\xi }{\text{'}}d\vec {\xi }_{1}^{'} = d\vec {\xi }d{{\vec {\xi }}_{1}}$ (см. [7]), то
$Q = \frac{1}{2}\int {\int {(\psi (\vec {\xi }{\text{'}}) + \psi (\vec {\xi }_{1}^{'}))} } I(\vec {\xi },{{\vec {\xi }}_{1}})d{{\vec {\xi }}_{1}}d\vec {\xi } \cong \frac{1}{2}\sum\limits_{\vec {l}}^n {\sum\limits_{\vec {k}}^n {{{{\vec {C}}}_{{\vec {k}}}}{{{\vec {C}}}_{{\vec {1}}}}(\psi (\vec {\xi }_{{kl}}^{'}) + \psi (\vec {\xi }_{{1kl}}^{'}))} } \bar {I}({{\vec {\xi }}_{k}},{{\vec {\xi }}_{l}}).$
Отсюда имеем

$\begin{gathered} \int {\int {I(\vec {\xi },{{{\vec {\xi }}}_{1}})\psi (\vec {\xi })d{{{\vec {\xi }}}_{1}}d\vec {\xi }} } \cong \sum\limits_{\vec {k}}^n {\sum\limits_{\vec {l}}^n {{{{\vec {C}}}_{{\vec {k}}}}{{{\vec {C}}}_{{\vec {1}}}}\psi ({{{\vec {\xi }}}_{k}})} } \bar {I}({{{\vec {\xi }}}_{k}},{{{\vec {\xi }}}_{l}}) = \\ = \; - \frac{1}{4}\sum\limits_{\vec {l}}^n {\sum\limits_{\vec {k}}^n {{{{\vec {C}}}_{{\vec {k}}}}{{{\vec {C}}}_{{\vec {1}}}}(\psi (\vec {\xi }_{{kl}}^{'}) + \psi (\vec {\xi }_{{1kl}}^{'}) - \psi ({{{\vec {\xi }}}_{k}}) + \psi ({{{\vec {\xi }}}_{l}}))} } \bar {I}({{{\vec {\xi }}}_{k}},{{{\vec {\xi }}}_{l}}). \\ \end{gathered} $

Получаем, что предложенная выше схема вычисления интеграла столкновений обладает такими же свойствами, что и классический интеграл столкновений Больцмана. Это дает возможность построить консервативную численную схему решения системы (3.1).

ЗАКЛЮЧЕНИЕ

В статье приведен окончательный вид кинетических уравнений, описывающих поведение неидеального газа и ее безразмерный вид. Получившиеся кинетические уравнения оказались не сложнее уравнения Больцмана, что позволило предложить численный метод для их решения. Этот метод без труда может быть реализован на имеющейся вычислительной базе.

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

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

  1. Бишаев А.М., Рыков В.А. Построение системы кинетических уравнений для неидеального газа // Теплофиз. высоких температур. 2017. Т. 5б. № 21. С. 31–43.

  2. Сивухин Д.В. Общий курс физики. Т. II. Термодинамика и молекулярная физика. М.: Наука, 1975. 552 с.

  3. Власов А.А. Нелокальная статистическая механика. М.: Наука, 1976. 264 с.

  4. Bishaev A.M., Rykov V.A., Abgaryan M.V. The H-theorem and equation of state for kinetic model of imperfect gas // J. Phys. 2018. Conf. Ser. P. 1–5.

  5. Черемисин Ф.Г. Решение уравнения Больцмана при движении газа с большими скоростями // Ж. вычисл. матем. и матем. физ. 2006. Т. 46. № 2. С. 329–343.

  6. Ларина И.Н., Рыков В.А. Численное решение уравнения Больцмана методом симметричного расщепления // Ж. вычисл. матем. и матем. физ. 2003. Т. 43. № 4. С. 601–6013.

  7. Абгарян М.В., Бишаев А.М. Модернизация метода расщепления для решения системы кинетических уравнений, описывающих поведение струи разреженной плазмы // Ж. вычисл. матем. и матем. физ. 2018. Т. 58. № 7. С. 1132–1146.

  8. Коган М.Н. Динамика разреженных газов. М.: Наука, 1967. 440 с.

  9. Рудяк В.Я. Статистическая теория диссипативных процессов в газах и жидкостях. М.: Наука, 1987. 270 с.

  10. Никифоров А.Ф., Уваров В.Б. Специальные функции математической физики. М.: Наука, 1978. 320 с.

  11. Титарев В.А., Шахов Е.М. Численный анализ винтового движения Куэтта разреженного газа между коаксиальными цилиндрами // Ж. вычисл. матем. и матем. физ. 2006. Т. 46. № 3. С. 527–535.

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