Известия РАН. Механика жидкости и газа, 2020, № 3, стр. 125-140

НЕИЗОТЕРМИЧЕСКОЕ ТЕЧЕНИЕ РАЗРЕЖЕННОГО ГАЗА В ДЛИННОМ ЦИЛИНДРИЧЕСКОМ КАНАЛЕ ПРИ ПРОИЗВОЛЬНЫХ ПЕРЕПАДАХ ДАВЛЕНИЯ И ТЕМПЕРАТУРЫ

О. В. Гермидер a*, В. Н. Попов a**

a Северный (Арктический) федеральный университет им. М.В. Ломоносова
Архангельск, Россия

* E-mail: o.germider@narfu.ru
** E-mail: v.popov@narfu.ru

Поступила в редакцию 16.09.2019
После доработки 26.11.2019
Принята к публикации 17.12.2019

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

Аннотация

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

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

Исследования течений разреженного газа в микро- и наноканалах имеют существенное значение для большого числа практических приложений [13]. Корректное описание таких течений возможно на основе кинетического уравнения Больцмана или модельных кинетических уравнений [1]. Так, для цилиндрического канала с использованием S-модели кинетического уравнения Больцмана [4] выполнен ряд исследований [510]. В [57] численное решение задачи о стационарном неизотермическом течении разреженного газа в длинной трубе кругового сечения без учета концевых эффектов получено на основе линеаризованной S-модели. Распределение температуры и давления газа вдоль канала в [57] полагалось линейным. В полной постановке, без упрощающих предположений, задача для достаточно длинных труб решена в [810] прямым методом конечных разностей во всей области течения. В частности, в [9] проведен анализ влияния длины трубы на распределение макроскопических величин (плотности, температуры, массовой скорости и потока массы газа) вдоль ее оси при истечении газа из камеры высокого давления в вакуум, выполнено сравнение с результатами, полученными по методике, представленной в [1]. При использовании данной методики (см., например, [2] и [12]) концевыми эффектами, включая процессы, протекающие в резервуарах, пренебрегают. Поток массы газа определяют локально в каждом сечении трубы по линейной теории трубы бесконечной длины, а распределение давления вдоль трубы – интегрированием одномерного уравнения неразрывности с граничными условиями равенства значений давления соответствующим значениям давления в резервуарах, соединенных рассматриваемой трубой. При этом распределение температуры считается заданным граничными условиями. В [9] показано, что упрощенная постановка задачи без учета одного или обоих граничных условий на концах трубы и приближенный метод расчета потока массы газа применимы только для достаточно длинных труб и достаточно разреженных газов. Критерий применимости – отношение параметра разреженности к безразмерной длине канала.

Цель представленной работы заключается в решении в упрощенной постановке [1] полуаналитическим методом задачи о стационарном неизотермическом течении газа в длинной трубе кругового сечения при произвольных заданных значениях давления и температуры на концах трубы. В качестве основного уравнения в работе используется S-модель, а в качестве граничного условия на стенках канала – модель диффузного отражения. Методом коллокации для полиномов и рациональных функций Чебышева [13, 14] находится решение линейной задачи с диффузными граничными условиями в постановке [6], рассчитываются потоки массы газа и тепла, обусловленные градиентами давления и температуры газа в канале, по всему диапазону значений числа Кнудсена, выполняется сравнение с результатами [2, 57] и [11]. С использованием метода Рунге–Кутта численно находится решение нелинейной задачи распределения давления при заданном линейном распределении температуры, которое имеет место в случае, когда теплопроводность и толщина стенок не изменяются вдоль канала. Особое внимание уделяется восстановлению распределения давления газа в канале. Как следует из выполненных численных расчетов, предложенный в работе полуаналитический метод за счет предварительной проработки построения решения модельного уравнения позволяет существенно уменьшить объем вычислений.

1. ПОСТАНОВКА ЗАДАЧИ. КИНЕТИЧЕСКОЕ УРАВНЕНИЕ

Рассмотрим течение разреженного газа между двумя резервуарами, соединяющимися длинным цилиндрическим каналом радиусом $R{\text{'}}$. Давление и температура в первом и во втором резервуарах остаются постоянными и равными соответственно ${{p}_{1}}$, ${{T}_{1}}$ и ${{p}_{2}}$, ${{T}_{2}}$. Полагаем, что ${{p}_{2}} > {{p}_{1}}$ и ${{T}_{2}} > {{T}_{1}}$. Ось $z{\text{'}}$ направим вдоль оси канала. Будем рассматривать течение газа в средней части канала. Считаем, что длина канала $L{\text{'}} \gg R{\text{'}}$. Распределение температуры в канале задается температурой стенок [1]. Полагаем, что безразмерные градиенты давления и температуры газа являются малыми по абсолютной величине

(1.1)
${{G}_{p}} = \frac{{R{\text{'}}}}{{{{p}_{0}}}}\frac{{dp}}{{dz{\text{'}}}},\quad {\text{|}}{{G}_{p}}{\text{|}} \ll 1,\quad {{G}_{T}} = \frac{{R{\text{'}}}}{{{{T}_{0}}}}\frac{{dT}}{{dz{\text{'}}}},\quad {\text{|}}{{G}_{T}}{\text{|}} \ll 1$
где ${{p}_{0}}$, ${{T}_{0}}$ – давление и температура газа в некоторой точке, принятой за начало координат. В качестве размерного масштаба длины выбран радиус $R{\text{'}}$ цилиндрического канала. Далее безразмерные величины будем обозначать без штриха. В линейном приближении для функций давления и температуры газа в канале получаем следующие выражения

$p(z) = {{p}_{0}}(1 + {{G}_{p}}z),\quad T(z) = {{T}_{0}}(1 + {{G}_{T}}z)$

Состояние разреженного газа в точке, радиус-вектор которой ${\mathbf{r}}$ имеет координаты $\rho $, $\varphi $ и $z$ в цилидрической системе координат в конфигурационном пространстве, определяем функцией распределения молекул газа $f({\mathbf{r}},{\mathbf{C}})$,где ${\mathbf{C}} = {{\beta }^{{1/2}}}{\mathbf{v}}$ – безразмерная молекуляная скорость, β = = $m{\text{/}}(2{{k}_{B}}{{T}_{0}})$, ${{k}_{B}}$ – постоянная Больцмана, m – масса молекул газа. В пространстве скоростей также будем использовать цилиндрические координаты $({{C}_{ \bot }},\psi ,{{C}_{z}})$: ${{C}_{\rho }} = {{C}_{ \bot }}cos\psi $, ${{C}_{\varphi }} = {{C}_{ \bot }}sin\psi $, Cz. Так как безразмерные градиенты давления и температуры малы, то функция распределения $f({\mathbf{r}},{\mathbf{C}})$ в линейном приближении имеет вид [15]

(1.2)
$f({\mathbf{r}},{\mathbf{C}}) = {{f}_{0}}(C)\left( {1 + {{G}_{T}}\left( {{{C}^{2}} - \frac{5}{2}} \right)z + {{G}_{p}}z + h(\rho ,{\mathbf{C}})} \right)$
где ${{f}_{0}}(C) = {{n}_{0}}{{(\beta {\text{/}}\pi )}^{{3/2}}}exp( - {{C}^{2}})$ – абсолютный максвеллиан, n0 – концентрация газа в начале координат. Функцию возмущения $h(\rho ,{\mathbf{C}})$ представляем как сумму, слагаемые в которой связаны с ${{G}_{p}}$ и ${{G}_{T}}$

(1.3)
$h(\rho ,{\mathbf{C}}) = {{G}_{p}}{{h}_{1}}(\rho ,{\mathbf{C}}) + {{G}_{T}}{{h}_{2}}(\rho ,{\mathbf{C}})$

Введем, следуя [1] и [16], безразмерные компоненты массовой скорости газа и вектора потока тепла

(1.4)
${{U}_{z}} = {{\beta }^{{1/2}}}{{u}_{z}},\quad {{q}_{z}} = \frac{{{{\beta }^{{1/2}}}}}{{{{p}_{0}}}}q_{z}^{'}$

Здесь через ${{u}_{z}}$ и $q_{z}^{'}$ обозначены соответствующие размерные величины [16]

(1.5)
${{u}_{z}}({\mathbf{r}}) = \frac{1}{{n(z)}}\int {{{v}_{z}}f({\mathbf{r}},{\mathbf{C}}){{d}^{3}}{\mathbf{v}}} $
(1.6)
$q_{z}^{'}({\mathbf{r}}) = \frac{m}{2}\int {{\text{(}}{{v}_{z}} - {{u}_{z}}{\text{)(}}v - u{{{\text{)}}}^{2}}f({\mathbf{r}},{\mathbf{C}}){{d}^{3}}{\mathbf{v}}} $
где $n(z)$ – концентрация молекул газа.

Подставляя (1.2) в (1.5) и (1.6), для безразмерных компонент (1.4) с учетом (1.3) получаем

${{U}_{z}}(\rho ) = {{G}_{p}}{{U}_{{z,1}}}(\rho ) + {{G}_{T}}{{U}_{{z,2}}}(\rho ),\quad {{q}_{z}}(\rho ) = {{G}_{p}}{{q}_{{z,1}}}(\rho ) + {{G}_{T}}{{q}_{{z,2}}}(\rho )$
(1.7)
${{U}_{{z,i}}}(\rho ) = \frac{1}{{{{\pi }^{{3/2}}}}}\int {exp( - {{C}^{2}}){{C}_{z}}{{h}_{i}}(\rho ,{\mathbf{C}}){{d}^{3}}{\mathbf{C}}} ,\quad i = 1,2$
(1.8)
${{q}_{{z,i}}}(\rho ) = \frac{1}{{{{\pi }^{{3/2}}}}}\int {exp( - {{C}^{2}}){{C}_{z}}\left( {{{C}^{2}} - \frac{5}{2}} \right){{h}_{i}}(\rho ,{\mathbf{C}}){{d}^{3}}{\mathbf{C}}} ,\quad i = 1,2$

Здесь ${{U}_{{z,1}}}$ и ${{q}_{{z,1}}}$ определяют массовую скорость и компоненту вектора потока тепла вследствие наличия градиента давления при изотермическом течении газа в канале, а ${{U}_{{z,2}}}$ и ${{q}_{{z,2}}}$ представляют собой их величину при изобарическом течении, обусловленном градиентом температуры.

Наша цель на первом этапе исследования заключается в нахождении приведенных потоков массы газа и тепла, которые зависят от Gp и ${{G}_{T}}$

(1.9)
$\begin{gathered} {{J}_{M}} = \frac{{J_{M}^{'}}}{{\pi R{\kern 1pt} {{'}^{2}}{\kern 1pt} p(z)}}\sqrt {\frac{{2{{k}_{B}}T(z)}}{m}} = {{G}_{p}}{{J}_{{M,1}}} + {{G}_{T}}{{J}_{{M,2}}} \\ {{J}_{Q}} = \frac{{2J_{Q}^{'}}}{{\pi R{\kern 1pt} {{'}^{2}}{\kern 1pt} p(z)}}\sqrt {\frac{m}{{2{{k}_{B}}T(z)}}} = {{G}_{p}}{{J}_{{Q,1}}} + {{G}_{T}}{{J}_{{Q,2}}} \\ \end{gathered} $

Здесь штрихами обозначены размерные величины потоков. Коэффициенты ${{J}_{{M,i}}}$ и ${{J}_{{Q,i}}}$ (i = 1, 2) являются безразмерными коэффициентами пропорциональности между потоками через поперечное сечение канала и локальными градиентами давления Gp и температуры GT

(1.10)
${{J}_{{M,i}}} = 4\int\limits_0^1 \,{{U}_{{z,i}}}(\rho )\rho d\rho ,\quad {{J}_{{Q,i}}} = 4\int\limits_0^1 \,{{q}_{{z,i}}}(\rho )\rho d\rho ,\quad i = 1,2$

В качестве основного уравнения, описывающего кинетику процессов, будем использовать S-модель кинетического уравнения Больцмана в цилиндрической системе координат [6, 17]

(1.11)
$\begin{gathered} \left( {\frac{{\partial h}}{{\partial \rho }}cos\psi - \frac{{\partial h}}{{\partial \psi }}\frac{{sin\psi }}{\rho }} \right){{C}_{ \bot }} + {{G}_{T}}\left( {{{C}^{2}} - \frac{5}{2}} \right){{C}_{z}} + {{G}_{p}}{{C}_{z}} + \delta h(\rho ,{\mathbf{C}}) = \\ = \frac{\delta }{{{{\pi }^{{3/2}}}}}\int {K({\mathbf{C}},{\mathbf{C}}{\text{'}})exp( - {{C}^{{'2}}})h(\rho ,{\mathbf{C}}{\text{'}}){{d}^{3}}{\mathbf{C}}{\text{'}}} \\ \end{gathered} $
где $\delta = {\text{K}}{{{\text{n}}}^{{ - 1}}}$ – параметр разрежения газа, ${\text{Kn}} = {{l}_{g}}{\text{/}}R{\text{'}}$ – число Кнудсена, ${{l}_{g}} = {{\eta }_{g}}{{\beta }^{{ - 1/2}}}{\text{/}}p$ – средняя длина свободного пробега молекул газа, ${{\eta }_{g}}$ – динамическая вязкость газа, $K({\mathbf{C}},{\mathbf{C}}{\text{'}})$ – ядро этого уравнения [17]

$K({\mathbf{C}},{\mathbf{C}}{\text{'}}) = \frac{4}{{15}}{\mathbf{CC}}{\text{'}}\left( {{{C}^{{'2}}} - \frac{5}{2}} \right)\left( {{{C}^{2}} - \frac{5}{2}} \right) + 2{\mathbf{C}}{\text{'}}{\mathbf{C}} + \frac{2}{3}\left( {{{C}^{{'2}}} - \frac{3}{2}} \right)\left( {{{C}^{2}} - \frac{3}{2}} \right) + 1$

Умножим левую и правую части уравнения (1.11) на $C_{z}^{k}exp( - C_{z}^{2}){\text{/}}\sqrt \pi $ (k = 1, 3) и проинтегрируем по Cz от $ - \infty $ до $ + \infty $. Полученные таким образом уравнения с учетом (1.3), (1.7) и (1.8) сводятся к двум независимым системам уравнений. Для восстановления массовой скорости ${{U}_{{z,1}}}(\rho )$ и компоненты вектора потока тепла ${{q}_{{z,1}}}(\rho )$ система уравнений имеет вид

(1.12)
$\left( {\frac{{\partial {{Z}_{{1,1}}}}}{{\partial \rho }}\zeta + \frac{{\partial {{Z}_{{1,1}}}}}{{\partial \zeta }}\frac{{(1 - {{\zeta }^{2}})}}{\rho }} \right){{C}_{ \bot }} + \delta {{Z}_{{1,1}}}(\rho ,{{C}_{ \bot }},\zeta ) + \frac{1}{2} = \delta \left( {{{U}_{{z,1}}}(\rho ) + \frac{2}{{15}}{{q}_{{z,1}}}(\rho )(C_{ \bot }^{2} - 1)} \right)$
(1.13)
$\left( {\frac{{\partial {{Z}_{{2,1}}}}}{{\partial \rho }}\zeta + \frac{{\partial {{Z}_{{2,1}}}}}{{\partial \zeta }}\frac{{(1 - {{\zeta }^{2}})}}{\rho }} \right){{C}_{ \bot }} + \delta {{Z}_{{2,1}}}(\rho ,{{C}_{ \bot }},\zeta ) + \frac{3}{4} = \delta \left( {\frac{3}{2}{{U}_{{z,1}}}(\rho ) + \frac{1}{5}{{q}_{{z,1}}}(\rho )C_{ \bot }^{2}} \right)$
где $\zeta = cos\psi $ и

(1.14)
${{Z}_{{1,i}}}(\rho ,{{C}_{ \bot }},\zeta ) = \frac{1}{{\sqrt \pi }}\int\limits_{ - \infty }^{ + \infty } \,exp( - C_{z}^{2}){{C}_{z}}{{h}_{i}}(\rho ,{\mathbf{C}})d{{C}_{z}}$
(1.15)
${{Z}_{{2,i}}}(\rho ,{{C}_{ \bot }},\zeta ) = \frac{1}{{\sqrt \pi }}\int\limits_{ - \infty }^{ + \infty } \,exp( - C_{z}^{2})C_{z}^{3}{{h}_{i}}(\rho ,{\mathbf{C}})d{{C}_{z}},\quad i = 1,2$

Компоненты ${{U}_{{z,i}}}(\rho )$ и ${{q}_{{z,i}}}(\rho )$ (i = 1, 2) выражаются через введенные функции ${{Z}_{{k,i}}}(\rho ,{{C}_{ \bot }},\zeta )$ ($k,i = 1,2$) как

(1.16)
${{U}_{{z,i}}}(\rho ) = \frac{2}{\pi }\int\limits_0^{ + \infty } \,exp( - C_{ \bot }^{2}){{C}_{ \bot }}\int\limits_{ - 1}^1 \,\frac{1}{{\sqrt {1 - {{\zeta }^{2}}} }}{{Z}_{{1,i}}}(\rho ,{{C}_{ \bot }},\zeta )d\zeta d{{C}_{ \bot }}$
(1.17)
$\begin{gathered} {{q}_{{z,i}}}(\rho ) = \frac{2}{\pi }\int\limits_0^{ + \infty } \,exp( - C_{ \bot }^{2}){{C}_{ \bot }}\int\limits_{ - 1}^1 \,\frac{1}{{\sqrt {1 - {{\zeta }^{2}}} }}{{Z}_{{2,i}}}(\rho ,{{C}_{ \bot }},\zeta )d\zeta d{{C}_{ \bot }} - \frac{5}{2}{{U}_{{z,i}}}(\rho ) + \\ \, + \frac{2}{\pi }\int\limits_0^{ + \infty } \,exp( - C_{ \bot }^{2})C_{ \bot }^{3}\int\limits_{ - 1}^1 \,\frac{1}{{\sqrt {1 - {{\zeta }^{2}}} }}{{Z}_{{1,i}}}(\rho ,{{C}_{ \bot }},\zeta )d\zeta d{{C}_{ \bot }} \\ \end{gathered} $

Для восстановления ${{U}_{{z,2}}}(\rho )$ и ${{q}_{{z,2}}}(\rho )$ система уравнений имеет вид

$\left( {\frac{{\partial {{Z}_{{1,2}}}}}{{\partial \rho }}\zeta + \frac{{\partial {{Z}_{{1,2}}}}}{{\partial \zeta }}\frac{{(1 - {{\zeta }^{2}})}}{\rho }} \right){{C}_{ \bot }} + \delta {{Z}_{{1,2}}}(\rho ,{{C}_{ \bot }},\zeta ) + \frac{{C_{ \bot }^{2} - 1}}{2} = \delta \left( {{{U}_{{z,2}}}(\rho ) + \frac{2}{{15}}{{q}_{{z,2}}}(\rho )(C_{ \bot }^{2} - 1)} \right)$
$\left( {\frac{{\partial {{Z}_{{2,2}}}}}{{\partial \rho }}\zeta + \frac{{\partial {{Z}_{{2,2}}}}}{{\partial \zeta }}\frac{{(1 - {{\zeta }^{2}})}}{\rho }} \right){{C}_{ \bot }} + \delta {{Z}_{{2,2}}}(\rho ,{{C}_{ \bot }},\zeta ) + \frac{{3C_{ \bot }^{2}}}{4} = \delta \left( {\frac{3}{2}{{U}_{{z,2}}}(\rho ) + \frac{1}{5}{{q}_{{z,2}}}(\rho )C_{ \bot }^{2}} \right)$

В качестве граничного условия на обтекаемой газом поверхности Γ канала будем использовать модель диффузного отражения [16]

${{f}^{ + }}({{{\mathbf{r}}}_{\Gamma }},{\mathbf{C}}) = {{f}_{\Gamma }}({{{\mathbf{r}}}_{\Gamma }},{\mathbf{C}}),\quad {\mathbf{C}}{{{\mathbf{e}}}_{n}} > 0$

Здесь ${{f}^{ + }}({{{\mathbf{r}}}_{\Gamma }},{\mathbf{C}})$ – функция распределения молекул газа, отраженных от поверхности канала, en – вектор единичной нормали к этой поверхности, направленный в сторону газа; ${{f}_{\Gamma }}({{{\mathbf{r}}}_{\Gamma }},{\mathbf{C}})$ – локально равновесная функция распределения [18]

(1.18)
${{f}_{\Gamma }}({{{\mathbf{r}}}_{\Gamma }},{\mathbf{C}}) = {{n}_{\Gamma }}(z)\mathop {\left( {\frac{m}{{2\pi {{k}_{B}}{{T}_{\Gamma }}(z)}}} \right)}\nolimits^{3/2} exp\left( { - \frac{{{{T}_{0}}}}{{{{T}_{\Gamma }}(z)}}{{{\mathbf{C}}}^{2}}} \right)$
где ${{T}_{\Gamma }}(z)$ – температура газа на поверхности канала, ${{n}_{\Gamma }}(z)$ – концентрация молекул газа вблизи стенок канала. Ввиду малости величин ${{G}_{p}}$ и ${{G}_{T}}$ функцию (1.18) линеаризуем (1.9) относительно абсолютного максвеллиана ${{f}_{0}}(C)$ [1]. В результате, учитывая (1.2), (1.3), для функций ${{Z}_{{k,i}}}(\rho ,{{C}_{ \bot }},\zeta )$ (k, i = 1, 2) получаем следующие граничные условия

(1.19)
${{Z}_{{k,i}}}(1,{{C}_{ \bot }},\zeta ) = 0,\quad \zeta < 0,\quad i,k = 1,2$

2. РЕШЕНИЕ КРАЕВОЙ ЗАДАЧИ. ПОТОКИ ГАЗА ПРИ МАЛЫХ ПЕРЕПАДАХ ДАВЛЕНИЯ И ТЕМПЕРАТУРЫ

Найдем решение системы уравнений (1.12) и (1.13) с граничными условиями (1.19). Неизвестные функции ${{Z}_{{1,i}}}(\rho ,{{C}_{ \bot }},\zeta )$ (i = 1, 2) раскладываем в ряд по полиномам и рациональным функциям Чебышева

(2.1)
${{Z}_{{1,i}}}(\rho ,{{C}_{ \bot }},\zeta ) = \sum\limits_{{\mathbf{j}} = 0}^\infty \,a_{{\mathbf{j}}}^{{(i)}}{{T}_{{{{j}_{1}}}}}({{\rho }^{*}}){{T}_{{{{j}_{2}}}}}(\zeta ){{T}_{{{{j}_{3}}}}}(C_{ \bot }^{*})$
(2.2)
$\rho {\text{*}} = 2\rho - 1,\quad C_{ \bot }^{*} = \frac{{{{C}_{ \bot }} - 1}}{{{{C}_{ \bot }} + 1}}$
где ${\mathbf{j}} = ({{j}_{1}},{{j}_{2}},{{j}_{3}})$. Для определения ${{T}_{j}}(x)$ в разложении (2.1) применяем реккурентные соотношения [13]

${{T}_{0}}(x) = 1,\quad {{T}_{1}}(x) = x,\quad {{T}_{j}}(x) = 2x{{T}_{{j - 1}}}(x) - {{T}_{{j - 2}}}(x),\quad j \geqslant 2$

Ограничиваясь в разложении (2.1) членами с номерами ${{j}_{k}} \leqslant {{n}_{k}}$ ($k = \overline {1,3} $), получаем

(2.3)
${{Z}_{{1,i}}}(\rho ,{{C}_{ \bot }},\zeta ) = \sum\limits_{{\mathbf{j}} = 0}^{\mathbf{n}} \,{{a}_{{\mathbf{j}}}}{{T}_{{{{j}_{1}}}}}(\rho {\text{*}}){{T}_{{{{j}_{2}}}}}(\zeta ){{T}_{{{{j}_{3}}}}}(C_{ \bot }^{*}) = {{{\mathbf{T}}}_{1}}(\rho {\text{*}}){\mathbf{Q}}(\zeta ){\mathbf{W}}(C_{ \bot }^{*}){{{\mathbf{A}}}_{{\mathbf{i}}}},\quad i = 1,2$

Здесь ${\mathbf{n}} = ({{n}_{1}},{{n}_{2}},{{n}_{3}})$, $n_{k}^{'} = {{n}_{k}} + 1$ ($k = \overline {1,3} $)

(2.4)
${{{\mathbf{T}}}_{{\mathbf{k}}}} = {{{\mathbf{T}}}_{{1 \times n_{k}^{'}}}},\quad {{{\mathbf{T}}}_{{1 \times n}}}(x) = \left( {{{T}_{0}}(x){{T}_{1}}(x) \ldots {{T}_{n}}(x)} \right)$
${\mathbf{Q}}(\zeta )$ и ${\mathbf{W}}({{C}_{ \bot }})$ – блочно-диагональные матрицы
(2.5)
${\mathbf{Q}}(\zeta ) = {{\Phi }_{{n_{1}^{'} \times n_{1}^{'}n_{2}^{'}}}}(\zeta ),\quad {\mathbf{W}}(C_{ \bot }^{*}) = {{\Phi }_{{n_{1}^{'}n_{2}^{'} \times n_{1}^{'}n_{2}^{'}n_{3}^{'}}}}(C_{ \bot }^{*})$
(2.6)
${{\Phi }_{{m \times mn}}}(x) = {\text{diag}}\left( {{{{\mathbf{T}}}_{{1 \times n}}}(x),{{{\mathbf{T}}}_{{1 \times n}}}(x), \ldots {{{\mathbf{T}}}_{{1 \times n}}}(x)} \right)$
${{{\mathbf{A}}}_{{\mathbf{i}}}} = {{{\mathbf{A}}}_{{{\mathbf{i}}n_{1}^{'}n_{2}^{'}n_{3}^{'} \times 1}}}$ – неизвестная матрица-столбец ($i = 1,2$)

(2.7)
${{{\mathbf{A}}}_{{\mathbf{i}}}} = \mathop {(a_{{000}}^{{(i)}}a_{{001}}^{{(i)}} \ldots a_{{00{{n}_{3}}}}^{{(i)}} \ldots a_{{{{n}_{1}}{{n}_{2}}0}}^{{(i)}}a_{{{{n}_{1}}{{n}_{2}}1}}^{{(i)}} \ldots a_{{{{n}_{1}}{{n}_{2}}{{n}_{3}}}}^{{(i)}})}\nolimits^T $

Производные от ${{{\mathbf{T}}}_{1}}(\rho {\text{*}})$ по переменной ρ и ${{{\mathbf{T}}}_{2}}(\zeta )$ по ζ могут быть определены как

(2.8)
$\frac{{d{{{\mathbf{T}}}_{1}}(\rho {\text{*}})}}{{d\rho }} = 2{{{\mathbf{T}}}_{1}}(\rho {\text{*}}){{{\mathbf{J}}}_{1}},\quad \frac{{d{{{\mathbf{T}}}_{2}}(\zeta )}}{{d\zeta }} = {{{\mathbf{T}}}_{2}}(\zeta ){{{\mathbf{J}}}_{2}}$
где ${{{\mathbf{J}}}_{{\mathbf{k}}}} = {{{\mathbf{J}}}_{{n_{k}^{'} \times n_{k}^{'}}}}$ – матрица (k = 1, 2), ненулевые элементы которой вычисляются по формулам [19]

(2.9)
${{({{J}_{k}})}_{{i,j + 1}}} = \left\{ \begin{gathered} j,\quad i = 1\;\;(j\;{\text{нечет}}.) \hfill \\ 2j,\quad j > i - 1\;\;(i\;{\text{чет}}.,j\;{\text{чет}}.) \hfill \\ 2j,\quad j > i - 1,\quad i > 1\;\;(i\;{\text{нечет}}.,j\;{\text{нечет}}.) \hfill \\ \end{gathered} \right.$

Например, для ${{n}_{1}} = 5$ $(n_{1}^{'} = 6)$ согласно (2.9) получаем

${{{\mathbf{J}}}_{1}} = \left( {\begin{array}{*{20}{c}} 0&1&0&3&0&5 \\ 0&0&4&0&8&0 \\ 0&0&0&6&0&{10} \\ 0&0&0&0&8&0 \\ 0&0&0&0&0&{10} \\ 0&0&0&0&0&0 \end{array}} \right)$

Для компонент ${{U}_{{z,1}}}(\rho {\text{*}})$ и ${{q}_{{z,1}}}(\rho {\text{*}})$ подставкой (2.3) в (1.16) и (1.17) приходим к следующим выражениям

(2.10)
${{U}_{{z,1}}}(\rho {\text{*}}) = 4{{{\mathbf{P}}}_{1}}(\rho {\text{*}}){{{\mathbf{A}}}_{1}}$
(2.11)
${{q}_{{z,1}}}(\rho {\text{*}}) = 4\left( {{{{\mathbf{P}}}_{3}}(\rho {\text{*}}) - \frac{5}{2}{{{\mathbf{P}}}_{1}}(\rho {\text{*}})} \right){{{\mathbf{A}}}_{1}} + 4{{{\mathbf{P}}}_{1}}(\rho {\text{*}}){{{\mathbf{A}}}_{2}}$

Здесь Pi (i = 1, 3) – матрицы размером $1 \times n_{1}^{'}n_{2}^{'}n_{3}^{'}$, ненулевыые элементы которых определяются как

(2.12)
$\begin{gathered} {{({{P}_{i}})}_{{1,n_{2}^{'}n_{3}^{'}{{j}_{1}} + {{j}_{3}} + 1}}}(\rho {\text{*}}) = {{I}_{{i,{{j}_{3}}}}}{{({{T}_{1}})}_{{1,{{j}_{1}} + 1}}}(\rho {\text{*}}),\quad {{j}_{{1,3}}} = \overline {0,{{n}_{{1,3}}}} \\ 2{{I}_{{i,{{j}_{3}}}}} = \int\limits_0^{ + \infty } \,exp( - C_{ \bot }^{2})C_{ \bot }^{i}{{T}_{{{{j}_{3}}}}}(C_{ \bot }^{*})d{{C}_{ \bot }} \\ \end{gathered} $

Выражение (2.12) получено с учетом ортогональности полиномов $\{ {{T}_{{{{j}_{2}}}}}(\zeta )\} $ (${{j}_{2}} = \overline {0,{{n}_{2}}} $). Интегралы ${{I}_{{i,{{j}_{3}}}}}$ (i = 1, 2, ${{j}_{3}} = \overline {0,{{n}_{3}}} $) находим численно по формуле

(2.13)
${{I}_{{i,{{j}_{3}};n}}} = \frac{4}{n}{{\sum\limits_{k = 0}^n }^{{{\text{''}}}}}{{\tau }_{{ex}}}({{\hat {C}}_{{ \bot ,k}}},n){{w}_{i}}({{\hat {C}}_{{ \bot ,k}}}){{T}_{{{{j}_{3}}}}}({{\hat {C}}_{{ \bot ,k}}}),\quad {{\tau }_{{ex}}}(C_{ \bot }^{*},n) = {{\sum\limits_{j = 0}^{n/2} }^{{{\text{''}}}}}\frac{{{{T}_{j}}({{{\hat {C}}}_{ \bot }})}}{{1 - 4{{j}^{2}}}}$
${{w}_{i}}(C_{ \bot }^{*}) = \frac{{exp\left( { - \tfrac{{{{{(1 + C_{ \bot }^{*})}}^{2}}}}{{{{{(1 - C_{ \bot }^{*})}}^{2}}}}} \right){{{(1 + C_{ \bot }^{*})}}^{i}}}}{{{{{(1 - C_{ \bot }^{*})}}^{{i + 2}}}}}$
${{\hat {C}}_{{ \bot ,k}}} = cos\left( {\frac{{\pi (n - k)}}{n}} \right),\quad k = \overline {0,n} $
где под ${{\sum\limits_{k = 0}^n }^{{{\text{''}}}}}$ понимаем сумму, в которой первое и последнее слагаемые умножаются на 1/2.

В (2.13) n – четное. При n = 40 абсолютная разность значений для ${{I}_{{i,{{j}_{3}}}}}$ (i = 1, 3, ${{j}_{3}} = \overline {0,{{n}_{3}}} $), полученных с использованием (2.13) и методом Гаусса, не превышает 10–9.

Выражая из равенств (2.2) переменные $\rho $ и ${{C}_{ \bot }}$ через $\rho {\text{*}}$ и $C_{ \bot }^{*}$ и учитывая (2.3)–(2.11), преобразуем систему уравнений (1.12) и (1.13) в виде

(2.14)
${{\Lambda }_{1}}(\rho {\text{*}},\zeta ,C_{ \bot }^{*}){{{\mathbf{A}}}_{1}} + {{\Lambda }_{2}}(\rho {\text{*}},C_{ \bot }^{*}){{{\mathbf{A}}}_{2}} = - \frac{1}{2}$
(2.15)
${{\Delta }_{1}}(\rho {\text{*}},C_{ \bot }^{*}){{{\mathbf{A}}}_{1}} + {{\Delta }_{2}}(\rho {\text{*}},\zeta ,C_{ \bot }^{*}){{{\mathbf{A}}}_{2}} = - \frac{3}{4}$

Здесь матрицы ${{\Lambda }_{{\mathbf{i}}}}$ и ${{\Delta }_{{\mathbf{i}}}}$ (i = 1, 2) имеют размер $1 \times n_{1}^{'}n_{2}^{'}n_{3}^{'}$ и определяются следующим образом

$\begin{gathered} {{\Lambda }_{1}}(\rho {\text{*}},\zeta ,C_{ \bot }^{*}) = \left( {{{\Psi }_{1}}(\rho {\text{*}},\zeta ) + \frac{{C_{ \bot }^{*} + 1}}{{1 - C_{ \bot }^{*}}}{{\Psi }_{2}}(\rho {\text{*}},\zeta )} \right){\mathbf{W}}(C_{ \bot }^{*}) - \\ \, - 4\delta {{{\mathbf{P}}}_{1}}(\rho {\text{*}}) + \frac{4}{3}\delta \left( {\frac{{{{{(C_{ \bot }^{*} + 1)}}^{2}}}}{{{{{(1 - C_{ \bot }^{*})}}^{2}}}} - 1} \right)\left( {{{{\mathbf{P}}}_{1}}(\rho {\text{*}}) - \frac{2}{5}{{{\mathbf{P}}}_{3}}(\rho {\text{*}})} \right) \\ \end{gathered} $
${{\Psi }_{1}}(\rho {\text{*}},\zeta )\delta {{{\mathbf{T}}}_{1}}(\rho {\text{*}}){\mathbf{Q}}(\zeta ),\quad {{\Psi }_{2}}(\rho {\text{*}},\zeta ) = 2\zeta {{{\mathbf{T}}}_{1}}(\rho {\text{*}}){{{\mathbf{J}}}_{1}}{\mathbf{Q}}(\zeta ) + \frac{{2(1 - {{\zeta }^{2}})}}{{\rho {\text{*}} + 1}}{{{\mathbf{T}}}_{1}}(\rho {\text{*}}){\mathbf{Q}}{\text{'}}(\zeta )$
${\mathbf{Q'}}(\zeta ) = {\text{diag}}\left( {{{{\mathbf{T}}}_{2}}(\zeta ){{{\mathbf{J}}}_{2}},{{{\mathbf{T}}}_{2}}(\zeta ){{{\mathbf{J}}}_{2}}, \ldots {{{\mathbf{T}}}_{2}}(\zeta ){{{\mathbf{J}}}_{2}}} \right)$
${{\Lambda }_{2}}(\rho {\text{*}},C_{ \bot }^{*}) = \frac{8}{{15}}\delta \left( {1 - \frac{{{{{(C_{ \bot }^{*} + 1)}}^{2}}}}{{{{{(1 - C_{ \bot }^{*})}}^{2}}}}} \right){{{\mathbf{P}}}_{1}}(\rho {\text{*}})$
${{\Delta }_{1}}(\rho {\text{*}},C_{ \bot }^{*})\frac{{2\delta {{{(C_{ \bot }^{*} + 1)}}^{2}}}}{{{{{(1 - C_{ \bot }^{*})}}^{2}}}}\left( {{{{\mathbf{P}}}_{1}}(\rho {\text{*}}) - \frac{2}{5}{{{\mathbf{P}}}_{3}}(\rho {\text{*}})} \right) - 6\delta {{{\mathbf{P}}}_{1}}(\rho {\text{*}})$
${{\Delta }_{2}}(\rho {\text{*}},\zeta ,C_{ \bot }^{*}) = \left( {{{\Psi }_{1}}(\rho {\text{*}},\zeta ) + \frac{{C_{ \bot }^{*} + 1}}{{1 - C_{ \bot }^{*}}}{{\Psi }_{2}}(\rho {\text{*}},\zeta )} \right){\mathbf{W}}(C_{ \bot }^{*}) - \frac{{4\delta {{{(C_{ \bot }^{*} + 1)}}^{2}}}}{{5{{{(1 - C_{ \bot }^{*})}}^{2}}}}{{{\mathbf{P}}}_{1}}(\rho {\text{*}})$

Граничные условия (1.19) записываем в форме

(2.16)
$\Omega (1,\zeta ,C_{ \bot }^{*}){{{\mathbf{A}}}_{{\mathbf{j}}}} = 0,\quad {{\Omega }_{{1 \times n_{1}^{'}n_{2}^{'}n_{3}^{'}}}}(1,\zeta ,C_{ \bot }^{*}) = {{{\mathbf{T}}}_{1}}(1){\mathbf{Q}}(\zeta ){\mathbf{W}}(C_{ \bot }^{*}),\quad \zeta < 0,\quad j = 1,2$

Умножим уравнения (2.14), (2.15) и (2.16) на ${{w}_{i}}(C_{ \bot }^{*})$ ($i = \overline {0,{{n}_{3}}} $) и проинтегрируем по $C_{ \bot }^{*}$ от –1 до 1. В результате получаем

(2.17)
$\Lambda _{{1,{\mathbf{i}}}}^{'}(\rho {\text{*}},\zeta ){{{\mathbf{A}}}_{1}} + \Lambda _{{2,{\mathbf{i}}}}^{'}(\rho {\text{*}}){{{\mathbf{A}}}_{2}} = - \frac{1}{2}{{I}_{{i,0}}},\quad i = \overline {0,{{n}_{3}}} $
(2.18)
$\Delta _{{1,{\mathbf{i}}}}^{'}(\rho {\text{*}}){{{\mathbf{A}}}_{1}} + \Delta _{{2,{\mathbf{i}}}}^{'}(\rho {\text{*}},\zeta ){{{\mathbf{A}}}_{2}} = - \frac{3}{4}{{I}_{{i,0}}},\quad i = \overline {0,{{n}_{3}}} $
(2.19)
$\Omega _{{\mathbf{i}}}^{'}(1,\zeta ){{{\mathbf{A}}}_{{\mathbf{j}}}} = 0,\quad \Omega _{{\mathbf{i}}}^{'}(1,\zeta ) = {{{\mathbf{T}}}_{1}}(1){\mathbf{Q}}(\zeta ){\mathbf{W}}_{{\mathbf{i}}}^{'},\quad \zeta < 0,\quad j = 1,2;\quad i = \overline {0,{{n}_{3}}} $

Здесь

$\Lambda _{{1,{\mathbf{i}}}}^{'}(\rho {\text{*}},\zeta ) = {{\Psi }_{1}}(\rho {\text{*}},\zeta ){\mathbf{W}}_{{\mathbf{i}}}^{'} + {{\Psi }_{2}}(\rho {\text{*}},\zeta ){\mathbf{W}}_{{{\mathbf{i}} + 1}}^{'} + \frac{4}{3}\delta {{{\mathbf{P}}}_{1}}(\rho {\text{*}})({{I}_{{i + 2,0}}} - 4{{I}_{{i,0}}}) + \frac{8}{{15}}\delta {{{\mathbf{P}}}_{3}}(\rho {\text{*}})({{I}_{{i,0}}} - {{I}_{{i + 2,0}}})$
$\Lambda _{{2,{\mathbf{i}}}}^{'}(\rho {\text{*}}) = \frac{8}{{15}}\delta ({{I}_{{i,0}}} - {{I}_{{i + 2,0}}}){{{\mathbf{P}}}_{1}}(\rho {\text{*}})$
$\Delta _{{1,{\mathbf{i}}}}^{'}(\rho {\text{*}})2\delta {{{\mathbf{P}}}_{1}}(\rho {\text{*}})({{I}_{{i + 2,0}}} - 3{{I}_{{i,0}}}) - \frac{4}{5}\delta {{{\mathbf{P}}}_{3}}(\rho {\text{*}}){{I}_{{i + 2,0}}}$
$\Delta _{{1,{\mathbf{i}}}}^{'}(\rho {\text{*}},\zeta ) = {{\Psi }_{1}}(\rho {\text{*}},\zeta ){\mathbf{W}}_{{\mathbf{i}}}^{'} + {{\Psi }_{2}}(\rho {\text{*}},\zeta ){\mathbf{W}}_{{{\mathbf{i}} + 1}}^{'} - \frac{4}{5}\delta {{{\mathbf{P}}}_{1}}(\rho {\text{*}}){{I}_{{i + 2,0}}}$
${\mathbf{W}}_{{{\mathbf{i}}n_{1}^{'}n_{2}^{'} \times n_{1}^{'}n_{2}^{'}n_{3}^{'}}}^{'} = {\text{diag}}\left( {{{{\mathbf{I}}}_{{\mathbf{i}}}},{{{\mathbf{I}}}_{{\mathbf{i}}}} \ldots {{{\mathbf{I}}}_{{\mathbf{i}}}}} \right),\quad {{{\mathbf{I}}}_{{{\mathbf{i}}1 \times n_{{13}}^{'}}}} = \left( {{{I}_{{i,0}}},{{I}_{{i,2}}}, \ldots {{I}_{{i,{{n}_{3}}}}}} \right)$

В качестве точек коллокации в уравнениях (2.17)(2.19) будем использовать нули полиномов ${{T}_{{n_{i}^{'}}}}(x)$ ($i = 1,2$) на отрезке [–1; 1]

(2.20)
${{x}_{{i,k}}} = cos\left( {\frac{{\pi (2{{n}_{i}} - 2{{k}_{i}} + 1)}}{{2({{n}_{i}} + 1)}}} \right),\quad {{k}_{i}} = \overline {0,{{n}_{i}}} ,\quad i = 1,2$
где ${{x}_{{1,k}}} = \rho _{{{{k}_{1}}}}^{*}$, ${{x}_{{2,k}}} = {{\zeta }_{{{{k}_{2}}}}}$.

Подставляя (2.20) в (2.17), приходим к системе линейных $n_{1}^{'}n_{2}^{'}n_{3}^{'}$-уравнений, в которой заменяем уравнения с $\rho {\text{*}} = \rho _{{{{n}_{1}}}}^{*}$ и $\zeta = {{\zeta }_{{{{k}_{2}}}}}$ на уравнения, вытекающие из граничного условия (2.19) для Z1, 1 со значениями переменных ρ* = 1 и $\zeta = {{\zeta }_{{{{k}_{2}}}}}$ (${{k}_{2}} = \overline {0,k_{2}^{'}} $). Здесь $k_{2}^{'}$ – индекс такой, что ${{\zeta }_{{k_{2}^{'}}}} < 0$ и ${{\zeta }_{{k_{{2 + 1}}^{'}}}} \geqslant 0$, т.е. $k_{2}^{'} = {{n}_{2}}{\text{/}}2 - 1$, если n2 – четное, и $k_{2}^{'} = ({{n}_{2}} - 1){\text{/}}2$ иначе.

В результате получаем

${{{\mathbf{B}}}_{1}}{{{\mathbf{A}}}_{1}} + {{{\mathbf{B}}}_{2}}{{{\mathbf{A}}}_{2}} = {{{\mathbf{F}}}_{1}}$
где ${{B}_{j}}$ (j = 1, 2) – матрицы размером $n_{1}^{'}n_{2}^{'}n_{3}^{'} \times n_{1}^{'}n_{2}^{'}n_{3}^{'}$, а F1 – матрица-столбец ($n_{1}^{'}n_{2}^{'}n_{3}^{'} \times 1$)

${{{\mathbf{B}}}_{1}} = \left( {\begin{array}{*{20}{c}} {\Lambda _{{1,0}}^{'}(\rho _{0}^{*},{{\zeta }_{0}})} \\ {\Lambda _{{1,1}}^{'}(\rho _{0}^{*},{{\zeta }_{0}})} \\ \ldots \\ {\Lambda _{{1,{{{\mathbf{n}}}_{3}} - 1}}^{'}(\rho _{{{{n}_{1}} - 1}}^{*},{{\zeta }_{{{{n}_{2}}}}})} \\ {\Lambda _{{1,{{{\mathbf{n}}}_{3}}}}^{'}(\rho _{{{{n}_{1}} - 1}}^{*},{{\zeta }_{{{{n}_{2}}}}})} \\ \ldots \\ {\Omega _{0}^{'}(1,{{\zeta }_{0}})} \\ {\Omega _{1}^{'}(1,{{\zeta }_{0}})} \\ \ldots \\ {\Omega _{{{{{\mathbf{n}}}_{3}}}}^{'}(1,{{\zeta }_{{k_{2}^{'}}}})} \\ {\Lambda _{{1,0}}^{'}(\rho _{{{{n}_{1}}}}^{*},{{\zeta }_{{k_{2}^{'} + 1}}})} \\ \ldots \\ {\Lambda _{{1,{{{\mathbf{n}}}_{3}} - 1}}^{'}(\rho _{{{{n}_{1}}}}^{*},{{\zeta }_{{{{n}_{2}}}}})} \\ {\Lambda _{{1,{{{\mathbf{n}}}_{3}}}}^{'}(\rho _{{{{n}_{1}}}}^{*},{{\zeta }_{{{{n}_{2}}}}})} \end{array}} \right),\quad {{{\mathbf{B}}}_{2}} = \left( {\begin{array}{*{20}{c}} {\Lambda _{{2,0}}^{'}(\rho _{0}^{*},{{\zeta }_{0}})} \\ {\Lambda _{{2,1}}^{'}(\rho _{0}^{*},{{\zeta }_{0}})} \\ \ldots \\ {\Lambda _{{2,{{{\mathbf{n}}}_{3}} - 1}}^{'}(\rho _{{{{n}_{1}} - 1}}^{*},{{\zeta }_{{{{n}_{2}}}}})} \\ {\Lambda _{{2,{{{\mathbf{n}}}_{3}}}}^{'}(\rho _{{{{n}_{1}} - 1}}^{*},{{\zeta }_{{{{n}_{2}}}}})} \\ \ldots \\ 0 \\ 0 \\ \ldots \\ 0 \\ {\Lambda _{{2,0}}^{'}(\rho _{{{{n}_{1}}}}^{*},{{\zeta }_{{k_{2}^{'} + 1}}})} \\ \ldots \\ {\Lambda _{{2,{{{\mathbf{n}}}_{3}} - 1}}^{'}(\rho _{{{{n}_{1}}}}^{*},{{\zeta }_{{{{n}_{2}}}}})} \\ {\Lambda _{{2,{{{\mathbf{n}}}_{3}}}}^{'}(\rho _{{{{n}_{1}}}}^{*},{{\zeta }_{{{{n}_{2}}}}})} \end{array}} \right),\quad {\mathbf{F}} = - \frac{1}{2}\left( {\begin{array}{*{20}{c}} {{{I}_{{0,0}}}} \\ {{{I}_{{1,0}}}} \\ \ldots \\ {{{I}_{{{{n}_{3}} - 1,0}}}} \\ {{{I}_{{{{n}_{3}},0}}}} \\ \ldots \\ 0 \\ 0 \\ \ldots \\ 0 \\ {{{I}_{{0,0}}}} \\ \ldots \\ {{{I}_{{{{n}_{3}} - 1,0}}}} \\ {{{I}_{{{{n}_{3}},0}}}} \end{array}} \right)$

Аналогично, подставляя (2.20) в (2.18), приходим к системе линейных $n_{1}^{'}n_{2}^{'}n_{3}^{'}$-уравнений, в которой заменяем уравнения с $\rho {\text{*}} = \rho _{{{{n}_{1}}}}^{*}$ и $\zeta = {{\zeta }_{{{{k}_{2}}}}}$ на уравнения, вытекающие из граничного условия (2.19) для Z2, 1 со значениями переменных ρ* = 1 и $\zeta = {{\zeta }_{{{{k}_{2}}}}}$ (${{k}_{2}} = \overline {0,k_{2}^{'}} $). В результате получаем

(2.21)
${{{\mathbf{D}}}_{1}}{{{\mathbf{A}}}_{1}} + {{{\mathbf{D}}}_{2}}{{{\mathbf{A}}}_{2}} = {{{\mathbf{F}}}_{2}}$
где ${{{\mathbf{F}}}_{2}} = 3{\text{/}}2{{{\mathbf{F}}}_{1}}$, а матрицы D1 и D2 определяются аналогично B2 и B1, в которых $\Lambda _{{2,{\mathbf{i}}}}^{'}$ и $\Lambda _{{1,{\mathbf{i}}}}^{'}$ заменяем соответственно на $\Delta _{{1,{\mathbf{i}}}}^{'}$ и $\Delta _{{2,{\mathbf{i}}}}^{'}$ ($i = \overline {0,{{n}_{3}}} $).

Для нахождения A1 получаем уравнение

(2.22)
$({{{\mathbf{B}}}_{1}} - {{{\mathbf{B}}}_{2}}{\mathbf{D}}_{2}^{{ - 1}}{{{\mathbf{D}}}_{1}}){{{\mathbf{A}}}_{1}} = {{{\mathbf{F}}}_{1}} - {{{\mathbf{B}}}_{2}}{\mathbf{D}}_{2}^{{ - 1}}{{{\mathbf{F}}}_{2}}$
где ${\mathbf{D}}_{2}^{{ - 1}}$ – обратная матрица для D2.

Решение (2.22) находим LU-методом в системе компьютерной алгебры Maple. На основе полученных элементов матрицы A1 согласно (2.10) получаем выражение для ${{U}_{{z,1}}}(\rho {\text{*}})$. Подставляя (2.10) в (1.10), имеем

(2.23)
${{J}_{{M,1}}} = 4\int\limits_{ - 1}^1 \,{{{\mathbf{P}}}_{1}}(\rho {\text{*}})({{T}_{0}}(\rho {\text{*}}) + {{T}_{1}}(\rho {\text{*}}))d\rho {\text{*}}{{{\mathbf{A}}}_{1}} = 8\sum\limits_{{{i}_{3}} = 0}^{{{n}_{3}}} \,{{I}_{{1,{{i}_{3}}}}}\left( {\sum\limits_{{{i}_{1}} = 0}^{[{{n}_{1}}/2]} \,\frac{{a_{{2n_{2}^{'}n_{3}^{'}{{i}_{1}} + {{i}_{3}} + 1,1}}^{{(1)}}}}{{1 - 4i_{1}^{2}}} + \sum\limits_{{{i}_{1}} = 0}^{[n_{1}^{'}/2] - 1} \,\frac{{a_{{n_{2}^{'}n_{3}^{'}(2{{i}_{1}} + 1) + {{i}_{3}} + 1,1}}^{{(1)}}}}{{4 - {{{(2{{i}_{1}} + 1)}}^{2}}}}} \right)$

В табл. 1 приведены значения $ - {{J}_{{M,1}}}$, полученные по формуле (2.23) при $n{\text{*}} = {{n}_{1}} = 2{{n}_{{2,3}}}$ для $\delta < 0.5$ и при $n{\text{*}} = {{n}_{1}} = 2{{n}_{{2,3}}}$ для $\delta \geqslant 0.5$. Там же для сравнения приведены значения компонент ${{J}_{{M,1}}}$ из [2, 57] и [11]. Результаты в [2, 5, 6] получены на основе решения линеаризованного S-модели методами дискретных скоростей и ординат. В [7] применен спектральный метод с использованием кубических сплайнов (δ = 0.1, 0.5 и 1), в [11] – консервативный метод, представляющий собой неявную схему второго порядка на неструктурированных сетках. Видно, что предложенный метод уже при n* = 10 дает в целом приемлемые результаты для $\delta > 0.02$. В свободномолекулярном пределе (δ = 0) значение JM, 1 при n* = 10 равно $ - 1.5026$, а при $n{\text{*}} = 12$${{J}_{{M,1}}} = - 1.5045$, которое соответствует вычисленному аналитически $ - 8{\text{/}}(3\sqrt \pi )$ [1].

Таблица 1
δ JM, 1
(2.23) [2] [5] [6], [7] [11]
n* = 10 n* = 20
0.01 1.4861 1.4798 1.4781 1.4800 1.4770 1.4765
0.02 1.4710 1.4631 1.4617 1.4636 1.4616 1.4611
0.05 1.4386 1.4334 1.4330 1.4339 1.4334 1.4329
0.1 1.4097 1.4090 1.4086 1.4101 1.4090 1.4085
0.2 1.3892 1.3899 1.3891 1.3911 1.3893
0.5 1.3986 1.4008 1.3996 1.4011 1.4005 1.3998
1.0 1.4740 1.4768 1.4748 1.4758 1.4764 1.4731
2.0 1.6749 1.6784 1.6772 1.6799 1.6779 1.6757
5.0 2.3615 2.3664 2.3646 2.3666 2.3655 2.3619
10.0 3.5726 3.5778 3.5752 3.5749 3.5762 3.5674

Матрицу A2 находим из уравнения (2.21). На основе полученных элементов матриц ${{{\mathbf{A}}}_{{\mathbf{i}}}}$ ($i = 1,2$) согласно (2.11) восстанавливаем компоненту ${{q}_{{z,1}}}(\rho {\text{*}})$. Выражение (1.10) для ${{J}_{{Q,1}}}$ приводим к виду

(2.24)
$\begin{gathered} {{J}_{{Q,1}}} = 8\sum\limits_{{{i}_{3}} = 0}^{{{n}_{3}}} \,{{I}_{{3,{{i}_{3}}}}}\left( {\sum\limits_{{{i}_{1}} = 0}^{[{{n}_{1}}/2]} \,\frac{{a_{{2n_{2}^{'}n_{3}^{'}{{i}_{1}} + {{i}_{3}} + 1,1}}^{{(1)}}}}{{1 - 4i_{1}^{2}}} + \sum\limits_{{{i}_{1}} = 0}^{[n_{1}^{'}/2] - 1} \,\frac{{a_{{n_{2}^{'}n_{3}^{'}(2{{i}_{1}} + 1) + {{i}_{3}} + 1,1}}^{{(1)}}}}{{4 - {{{(2{{i}_{1}} + 1)}}^{2}}}}} \right) + \\ \, + 8\sum\limits_{{{i}_{3}} = 0}^{{{n}_{3}}} \,{{I}_{{1,{{i}_{3}}}}}\left( {\sum\limits_{{{i}_{1}} = 0}^{[{{n}_{1}}/2]} \,\frac{{a_{{2n_{2}^{'}n_{3}^{'}{{i}_{1}} + {{i}_{3}} + 1,1}}^{{(2)}}}}{{1 - 4i_{1}^{2}}} + \sum\limits_{{{i}_{1}} = 0}^{[n_{1}^{'}/2] - 1} \,\frac{{a_{{n_{2}^{'}n_{3}^{'}(2{{i}_{1}} + 1) + {{i}_{3}} + 1,1}}^{{(2)}}}}{{4 - {{{(2{{i}_{1}} + 1)}}^{2}}}}} \right) - \frac{5}{2}{{J}_{{M,1}}}. \\ \end{gathered} $

В табл. 2 приведены значения ${{J}_{{Q,1}}}$, полученные по формуле (2.24) при тех же значениях n, которые использованы ранее для ${{J}_{{M,1}}}$. Для нахождения компонент ${{U}_{{z,2}}}(\rho {\text{*}})$ и ${{q}_{{z,2}}}(\rho {\text{*}})$ и вычисления значений соответствующих приведенных потоков используем предложенный подход. В этом случае систему уравнений (1.14) и (1.15) с граничными условиями (1.19) приводим к виду

${{{\mathbf{B}}}_{1}}{\mathbf{A}}_{1}^{*} + {{{\mathbf{B}}}_{2}}{\mathbf{A}}_{2}^{*} = {\mathbf{F}}_{1}^{*}$
${{{\mathbf{D}}}_{1}}{\mathbf{A}}_{1}^{ * } + {{{\mathbf{D}}}_{2}}{\mathbf{A}}_{2}^{ * } = {\mathbf{F}}_{2}^{ * }$
где неизвестные матрицы ${\mathbf{A}}_{{\mathbf{k}}}^{*}$ (k = 1, 2)
${{Z}_{{2,j}}}(\rho ,{{C}_{ \bot }},\zeta ) = {{{\mathbf{T}}}_{1}}(\rho {\text{*}}){\mathbf{Q}}(\zeta ){\mathbf{W}}(C_{ \bot }^{*}){\mathbf{A}}_{{\mathbf{j}}}^{*},\quad j = 1,2$
а матрицы-столбцы ${\mathbf{F}}_{1}^{*}$ и ${\mathbf{F}}_{2}^{*}$ определяются аналогично ${{{\mathbf{F}}}_{1}}$ и ${{{\mathbf{F}}}_{2}}$, в которых ${{I}_{{i,0}}}$ заменяем соответственно на ${{I}_{{i + 2,0}}} - {{I}_{{i,0}}}$ и ${{I}_{{i + 2,0}}}$ ($i = \overline {0,{{n}_{3}}} $).

Таблица 2
δ JQ, 1 = JM, 2 JQ, 2
(2.23) [2] [5] [6], [7] [11], JM, 2 [11], JQ, 1 (2.23) [6], [7]
n* = 10 n* = 10 n* = 10 n* = 20
0.01 0.7301 0.7238 0.7223 0.7243 0.7210 0.7188 0.7207 3.3074 3.2912 3.2848
0.02 0.7112 0.7033 0.7022 0.7042 0.7020 0.6998 0.7016 3.2390 3.2197 3.2166
0.05 0.6678 0.6629 0.6627 0.6637 0.6630 0.6609 0.6626 3.0750 3.0636 3.0638
0.1 0.6214 0.6208 0.6469 0.6206 0.6209 0.6192 0.6189 0.6205 2.8799 2.8800
0.2 0.5663 0.5672 0.5896 0.5667 0.5653 0.5668 2.6159 2.6190
0.5 0.4777 0.4785 0.4953 0.4780 0.4784 0.4767 0.4780 2.1334 2.1363 2.1360
1.0 0.3961 0.3969 0.4087 0.3962 0.3968 0.3947 0.3946 1.6723 1.6749 1.6745
2.0 0.3022 0.3028 0.3095 0.3022 0.3027 0.3020 0.3020 1.1778 1.1797 1.1794
5.0 0.1757 0.1763 0.1781 0.1759 0.1762 0.1759 0.1759 0.6175 0.6186 0.6186
10.0 0.1012 0.1020 0.1024 0.1018 0.1020 0.1019 0.1023 0.3403 0.3411 0.3410

В табл. 2 приведены значения ${{J}_{{M,2}}}$ и ${{J}_{{Q,2}}}$ в зависимости от δ. Из табл. 1 и 2 видно, что наблюдается быстрая сходимость в среднем с увеличением ${{n}_{{1,2,3}}}$. Необходимо заметить, что совпадение значений ${{J}_{{M,2}}}$ и ${{J}_{{Q,1}}}$ подтверждается соотношением Онзагера [1] и служит дополнительным критерием точности полученных результатов.

3. ПОТОКИ ГАЗА ПРИ ПРОИЗВОЛЬНЫХ ПЕРЕПАДАХ ДАВЛЕНИЯ И ТЕМПЕРАТУРЫ

Здесь будет описан второй этап проблемы, а именно – будут получены массовые потоки как функции давлений ${{p}_{1}}$ и ${{p}_{2}}$ и температур ${{T}_{1}}$ и ${{T}_{2}}$ на концах канала. В случае малых перепадов температуры и давления на концах канала, распределение температуры и давления вдоль канала можно считать линейным [1]. При этом градиенты температуры и давления могут быть определены по формулам [1]

${{G}_{T}} = \frac{{{{T}_{2}} - {{T}_{1}}}}{{L{{T}_{{av}}}}},\quad {{G}_{p}} = \frac{{{{p}_{2}} - {{p}_{1}}}}{{L{{p}_{{av}}}}}$
где $L = L{\text{'/}}R{\text{'}}$, ${{T}_{{av}}} = ({{T}_{2}} + {{T}_{1}}){\text{/}}2$, ${{p}_{{av}}} = ({{p}_{2}} + {{p}_{1}}){\text{/}}2$ и перепады температуры и давления являются малыми: $({{T}_{2}} - {{T}_{1}}) \ll {{T}_{1}}$ и $({{p}_{2}} - {{p}_{1}}) \ll {{p}_{1}}$. В этом случае величина JM остается постоянной. Если отношения ${{T}_{2}}{\text{/}}{{T}_{1}}$ и ${{p}_{2}}{\text{/}}{{p}_{1}}$ являются большими, то распределение давления перестает быть линейным и происходит изменение величины JM вдоль канала.

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

(3.1)
$J_{M}^{*} = \frac{L}{{\pi {{R}^{{'2}}}{{p}_{1}}}}J_{M}^{'}\sqrt {\frac{{2{{k}_{B}}{{T}_{1}}}}{m}} $

Выражая $J_{M}^{'}$ из (1.9) и подставляя в (3.1), с учетом (1.1) получаем

(3.2)
$J_{M}^{*} = \frac{{L\sqrt {{{T}_{1}}} }}{{{{p}_{1}}}}\left( {\frac{1}{{\sqrt T }}{{J}_{{M,1}}}\frac{{dp}}{{dz}} + \frac{p}{{{{T}^{{3/2}}}}}{{J}_{{M,2}}}\frac{{dT}}{{dz}}} \right)$

Пусть $z{\text{*}} = z{\text{'/}}L{\text{'}} = z{\text{/}}L$, $T{\text{*}} = T{\text{/}}{{T}_{1}}$ и $p{\text{*}} = p{\text{/}}{{p}_{1}}$, тогда выражение (3.2) перепишем в виде

(3.3)
$J_{M}^{*} = \frac{1}{{\sqrt {T{\text{*}}} }}{{J}_{{M,1}}}\frac{{dp{\text{*}}}}{{dz{\text{*}}}} + \frac{{p{\text{*}}}}{{{{T}^{{*3/2}}}}}{{J}_{{M,2}}}\frac{{dT{\text{*}}}}{{dz{\text{*}}}}$

Необходимо заметить, что в отличие от $J_{M}^{*}$, постоянного вдоль канала, величины ${{J}_{{M,1}}}$ и ${{J}_{{M,2}}}$ зависят от z. Зависимость эта проявляется неявно, через параметр δ. Для модели жестких сфер имеем следующие соотношения [1]

(3.4)
$\delta = \frac{{p{{T}_{i}}}}{{{{p}_{i}}T}}{{\delta }_{i}},\quad {{\delta }_{i}} = \frac{{R{\text{'}}{{p}_{i}}}}{{{{\eta }_{g}}({{T}_{i}})}}\sqrt {\frac{m}{{2{{k}_{B}}{{T}_{i}}}}} ,\quad i = 1,2$

Рассмотрим изотермическое течение. В этом случае $T{\text{*}} = 1$ (${{T}_{1}} = {{T}_{2}}$) и уравнение (3.3) принимает вид

(3.5)
$J_{M}^{*}dz{\text{*}} = {{J}_{{M,1}}}dp{\text{*}}$

Выражая p* из (3.4) и подставляя в (3.3), интегрируем левые правые части (3.5). Принимая во внимание, что $J_{M}^{*}$ не зависит от $dz{\text{*}}$, получаем

(3.6)
$J_{M}^{*} = \int\limits_1^{{{p}_{2}}/{{p}_{1}}} \,{{J}_{{M,1}}}dp{\text{*}} = \frac{1}{{{{\delta }_{1}}}}\int\limits_{{{\delta }_{1}}}^{{{\delta }_{2}}} \,{{J}_{{M,1}}}d\delta $

Значения интеграла (3.6) находим численно по формуле

(3.7)
$J_{M}^{*} = \frac{{2({{\delta }_{2}} - {{\delta }_{1}})}}{{m{{\delta }_{1}}}}\sum\limits_{k = 0}^m {\text{''}}{{\tau }_{{ex}}}(\delta _{k}^{*},m){{J}_{{M,1}}}({{\delta }^{{(k)}}})$
(3.8)
${{\delta }^{{(k)}}} = \left( {\frac{{{{\delta }_{2}} - {{\delta }_{1}}}}{2}} \right)\delta _{k}^{*} + \frac{{{{\delta }_{1}} + {{\delta }_{2}}}}{2},\quad \delta _{k}^{*} = cos\left( {\frac{{\pi (m - k)}}{m}} \right),\quad k = \overline {0,m} $
где функция ${{\tau }_{{ex}}}(\delta _{k}^{*},m)$ определена выше (2.13). Значения ${{J}_{{M,1}}}({{\delta }^{{(k)}}})$ в (3.7) находим, используя разложение этой величины на отрезке $\delta \in [0,10]$ в ряд по полиномам Чебышева [13], где в качестве интерполяционных узлов при вычислении значений функции (2.23) выбираем точки (3.8). Как показывают расчеты, при числе членов ряда 33 абсолютная погрешность составляет порядка 10–3. В третьем, четвертом и пятом столбцах табл. 3 представлены значения $ - J_{M}^{*}$, полученные по формуле (3.7) при m = 2, 4 и 8 на отрезках [0.1, 1], [0.1, 10] и [1, 10]. В шестом столбце этой таблицы приведены значения, вычисленные по теореме Лагранжа о среднем значении интеграла $ - J_{M}^{*} = {{\kappa }_{p}}{{J}_{{M,1}}}(\delta {\text{'}})$, где ${{\kappa }_{p}} = ({{\delta }_{1}} - {{\delta }_{2}}){\text{/}}{{\delta }_{1}}$, $\delta {\text{'}} = ({{\delta }_{1}} + {{\delta }_{2}}){\text{/}}2$. В последнем столбце – значения $J_{M}^{*}$, полученные в [20] на основе БГК модели кинетического уравнения. Из табл. 3 видно, что наблюдается быстрая сходимость в среднем при небольших значениях m в (3.7), а промежуточное значение ${{J}_{{M,1}}}(\delta {\text{'}})$ при этом близко к среднему арифметическому, что соответствует предположению из [1].

Таблица 3
$J_{M}^{*}$ при T* = 1 $J_{M}^{*}$ при p* = 1 и $T_{2}^{*}$ = 3.8
δ1 δ2 (3.7) κpJM, 1(δ') [20] δ1 (3.11) TJM, 2(δ')
m = 2 m = 4 m = 8 m = 2 m = 4 m = 8
0.1 1.0 12.7720 12.7499 12.7494 12.6628 12.6270 0.1 0.6426 0.6375 0.6376 0.6007
1.0 10 22.4956 22.4719 22.4711 22.3691 22.3200 1.0 0.4579 0.4533 0.4535 0.4191
0.1 10 239.2577 237.5450 237.4764 235.4562 235.8180 10 0.1674 0.1630 0.1630 0.1376

Используя найденные значения параметра $J_{M}^{*}$ по формуле (3.6), получаем из (3.5) с учетом (3.4) дифференциальное уравнение для $p{\text{*}}(z{\text{*}})$

$\frac{{dp{\text{*}}}}{{dz{\text{*}}}} = \frac{{J_{M}^{*}}}{{{{J}_{{M,2}}}({{\delta }_{1}}p{\text{*}}(z{\text{*}}))}}$
с граничным условием $p{\text{*}}( - 1{\text{/}}2) = 1$.

На рис. 1 и 2 приведены полученные методом Рунге–Кутта распределения $p{\text{*}}(z{\text{*}})$ при значениях δ1 и δ2, указанных в табл. 3. Это соответственно кривые 1 и 2 для ${{\delta }_{1}} = 0.1$, 1 и ${{\delta }_{2}} = 1$, 10 на рис. 1 и кривая для ${{\delta }_{1}} = 0.1$ и ${{\delta }_{2}} = 10$ на рис. 2. Из рис. 1 и 2 видно, что распределение $p{\text{*}}(z{\text{*}})$ приближается к линейному при ${{\delta }_{1}} = 0.1$ и ${{\delta }_{2}} = 1$ (кривая 1).

Рис. 1.

Распределение давления $p{\text{*}}(z{\text{*}})$ при изотермическом режиме течения газа в канале: 1, 2$({{\delta }_{1}},{{\delta }_{2}}) = (0.1,1)$, $(1,10)$.

Рис. 2.

Распределение давления $p{\text{*}}(z{\text{*}})$ при изотермическом режиме течения газа в канале: $({{\delta }_{1}},{{\delta }_{2}}) = (0.1,10)$.

Рассмотрим изобарическое течение. В этом случае p* = 1 (${{p}_{1}} = {{p}_{2}}$) и уравнение (3.3) приводится к виду

(3.9)
$J_{M}^{*}dz{\text{*}} = T{{{\text{*}}}^{{ - 3/2}}}{{J}_{{M,2}}}dT{\text{*}}$

Интегрируем (3.9). Учитывая (3.4), получаем

(3.10)
$J_{M}^{*} = \int\limits_1^{{{T}_{2}}/{{T}_{1}}} \,T{{{\text{*}}}^{{ - 3/2}}}{{J}_{{M,2}}}dT{\text{*}} = - \delta _{1}^{{ - 1/2}}\int\limits_{{{\delta }_{1}}}^{{{\delta }_{2}}} \,{{\delta }^{{ - 1/2}}}{{J}_{{M,2}}}d\delta $

Значения интеграла (3.10) находим численно по формуле

(3.11)
где ${{\tau }_{{ex}}}(\delta _{k}^{*},m)$, $\delta _{k}^{*}$ и ${{\delta }^{{(k)}}}$ определяем согласно (2.13) и (3.8). В табл. 3 представлены значения $J_{M}^{*}$, вычисленные по формуле (3.11) при m = 2, 4 и 8 для $T_{2}^{*} = 3.8$ и ${{\delta }_{1}} = 0.1$, 1 и 10. Отношение температур $T_{2}^{*} = 3.8$ соответствует комнатной температуре ${{T}_{2}} = 293$ К, поддерживаемой во втором резервуаре, и температуре жидкого азота ${{T}_{1}} = 77.2$ К в первом резервуаре [2]. Интервалы для $\delta $ при указанных в табл. 3 значениях ${{\delta }_{1}}$ содержатся в [0, 10], поэтому как и в случае изотермического течения газа для интерполируемой функции ${{J}_{{M,2}}}(\delta )$ используем полиномы Чебышева. В последнем столбце табл. 3 приведены значения, вычисленные по теореме Лагранжа о среднем значении интеграла $J_{M}^{*} = ({{\kappa }_{T}}{{J}_{{M,2}}})(\delta {\text{'}})$, где ${{\kappa }_{T}}(\delta ) = ({{\delta }_{2}} - {{\delta }_{1}}){{(\delta {{\delta }_{1}})}^{{ - 1/2}}}$, $\delta {\text{'}} = ({{\delta }_{1}} + {{\delta }_{2}}){\text{/}}2$.

Используя найденные значения параметра $J_{M}^{*}$ в этом случае по формуле (3.10), получаем из (3.9) с учетом (3.4) дифференциальное уравнение для $T{\text{*}}(z{\text{*}})$

$\frac{{dT{\text{*}}}}{{dz{\text{*}}}} = \frac{{J_{M}^{*}T{{{\text{*}}}^{{3/2}}}(z{\text{*}})}}{{{{J}_{{M,1}}}\left( {\tfrac{{{{\delta }_{1}}}}{{T{\text{*}}(z{\text{*}})}}} \right)}}$
с граничным условием $T{\text{*}}( - 1{\text{/}}2) = 1$.

На рис. 3 приведены полученные методом Рунге–Кутта распределения $T{\text{*}}(z{\text{*}})$ при $T_{2}^{*} = 3.8$ и значениях ${{\delta }_{1}}$ из табл. 3. Это соответственно кривые 1–3. Из рис. 3 видно, что распределение $T{\text{*}}(z{\text{*}})$ приближается к линейному с ростом ${{\delta }_{1}}$ (кривая 3).

Рис. 3.

Распределение температуры $T{\text{*}}(z{\text{*}})$ для $T_{2}^{*} = 3.8$ при изобарическом режиме течения газа в канале: 1–3${{\delta }_{1}} = 0.1$, 1 и 10.

Рассмотрим течение при T* и p*, отличных от единицы. В этом случае полагаем, что теплопроводность и толщина стенок не изменяются вдоль канала. Тогда распределение температур можно считать линейным [2]: $T{\text{*}}(z{\text{*}}) = (T_{2}^{*} - 1)z{\text{*}} + (1 + T_{2}^{*}){\text{/}}2$. Распределение давления $p{\text{*}}(z{\text{*}})$ заранее неизвестно и должно быть найдено в результате решения уравнения (3.3), в котором $J_{M}^{*}$ является параметром. Значения $p{\text{*}}( - 1{\text{/}}2) = 1$ и $p{\text{*}}(1{\text{/}}2) = {{p}_{2}}{\text{/}}{{p}_{1}}$ определяют граничные условия. Пусть $\tau = \sqrt {T{\text{*}}} $, тогда уравнение (3.3) сводится к

(3.12)
$\frac{{dp{\text{*}}}}{{d\tau }} = 2\left( {\frac{{J_{M}^{*}{{\tau }^{2}}}}{{(T_{2}^{*} - 1){{J}_{{M,1}}}}} - \frac{{p{\text{*}}{{J}_{{M,2}}}}}{{\tau {{J}_{{M,1}}}}}} \right)$
с граничными условиями $p{\text{*}}(1) = 1$ и $p{\text{*}}(\sqrt {{{T}_{2}}{\text{/}}{{T}_{1}}} ) = {{p}_{2}}{\text{/}}p1$.

В свободномолекулярном режиме ($\delta = 0$) имеем ${{J}_{{M,2}}} = - {{J}_{{M,1}}}{\text{/}}2$, а общее решение дифференциального уравнения (3.12) имеет вид

(3.13)
$p{\text{*}}(\tau ) = \left( {\frac{{J_{M}^{*}{{\tau }^{2}}}}{{{{J}_{{M,1}}}(T_{2}^{*} - 1)}} + C} \right)\tau $
где C – константа интегрирования. Подставляя в (3.13) граничные условия $p{\text{*}}(1) = 1$ и $p{\text{*}}({{\tau }_{2}}) = p_{2}^{*}$, приходим к системе уравнений для определения C и $J_{M}^{*}$. В результате получаем

(3.14)
$J_{M}^{*} = \left( {\frac{{p_{2}^{*}}}{{{{\tau }_{2}}}} - 1} \right){{J}_{{M,1}}},\quad C = \frac{{p_{2}^{*} - \tau _{2}^{3}}}{{(1 - \tau _{2}^{2}){{\tau }_{2}}}}$

Для цилидрического канала в свободномолекулярном режиме ${{J}_{{M,1}}} = - 8{\text{/}}(3\sqrt \pi )$ = –1.5045. При $T_{2}^{*}$ = 3.8 для $p_{2}^{*} = 100$ и $p_{2}^{*} = 1$ значения $J_{M}^{*}$, найденные согласно (3.14), равны соответственно –75.6747 и 0.7327.

Подставляя (3.14) в (3.13) и учитывая, что $\tau = \sqrt {T{\text{*}}(z{\text{*}})} $ и $T{\text{*}}(z{\text{*}}) = (T_{2}^{*} - 1)z{\text{*}} + (1 + T_{2}^{*}){\text{/}}2$, получаем распределение давления в канале $p{\text{*}}(z{\text{*}})$

(3.15)
$p{\text{*}}(z{\text{*}}) = \left( {\frac{1}{2}\left( {\frac{{p_{2}^{*}}}{{\sqrt {T_{2}^{*}} }} - 1} \right)z{\text{*}} + \frac{1}{4}\left( {\frac{{p_{2}^{*}}}{{\sqrt {T_{2}^{*}} }} + 1} \right)} \right)\sqrt {4(T_{2}^{*} - 1)z{\text{*}} + 2(T_{2}^{*} + 1)} $

На рис. 4 и 5 штриховой линией (кривая 4) показаны распределения $p{\text{*}}(z{\text{*}})$, полученные по (3.15) при $T_{2}^{*} = 3.8$ для $p_{2}^{*} = 100$ и $p_{2}^{*} = 1$.

Рис. 4.

Распределение давления $p{\text{*}}(z{\text{*}})$ для $p_{2}^{*} = 100$ при $T_{2}^{*} = 3.8$: 1–4${{\delta }_{1}} = 0.1$, 1, 10 и 0.

Рис. 5.

Распределение давления $p{\text{*}}(z{\text{*}})$ для $p_{2}^{*} = 1$ при $T_{2}^{*} = 3.8$: 1–4${{\delta }_{1}} = 0.1$, 1, 10 и 0.

В промежуточном режиме ${{J}_{{M,1}}}$ и ${{J}_{{M,2}}}$ зависят от параметра δ, значения которого изменяются от δ1 до δ2. Учитывая, что $p{\text{*}} = p{\text{/}}{{p}_{1}}$, выражаем p* из (3.4) и дифференцируем полученное равенство. В результате

(3.16)
$p{\text{*}} = \frac{{\delta {{\tau }^{2}}}}{{{{\delta }_{1}}}},\quad \frac{{dp{\text{*}}}}{{d\tau }} = {{\delta }_{1}}\left( {\frac{{d\delta }}{{d\tau }}{{\tau }^{2}} + 2\delta \tau } \right)$

Подставляя (3.16) в (3.12), получаем дифференциальное уравнение для $\delta (\tau )$

(3.17)
$\frac{{d\delta }}{{d\tau }} = \frac{{2J_{M}^{*}{{\delta }_{1}}}}{{(T_{2}^{*} - 1){{J}_{{M,1}}}(\delta )}} - \frac{{2\delta }}{\tau }\left( {\frac{{{{J}_{{M,2}}}(\delta )}}{{{{J}_{{M,1}}}(\delta )}} + 1} \right)$
с граничным условием $\delta (1) = {{\delta }_{1}}$. Решение уравнения находим методом Рунге–Кутта, подбирая параметр $J_{M}^{*}$ так, чтобы $\delta (\sqrt {{{T}_{2}}/{{T}_{1}}} ) = {{\delta }_{2}}$ c абсолютной погрешностью не превышающей 10–4. В качестве начального значения $J_{M}^{*}$ используем решение уравнения
(3.18)
${{\delta }_{2}} - {{\delta }_{1}} = \frac{{2J_{M}^{*}{{\delta }_{1}}}}{{(T_{2}^{*} - 1){{J}_{{M,1}}}(\delta {\text{'}})}}({{\tau }_{2}} - {{\tau }_{1}}) - 2\delta {\text{'}}ln\left( {\frac{{{{\tau }_{2}}}}{{{{\tau }_{1}}}}} \right)\left( {\frac{{{{J}_{{M,2}}}(\delta {\text{'}})}}{{{{J}_{{M,1}}}(\delta {\text{'}})}} + 1} \right)$
где $\delta {\text{'}} = ({{\delta }_{1}} + {{\delta }_{2}}){\text{/}}2$.

В табл. 4 представлены значения $J_{M}^{*}$, полученные из (3.18) и (3.17) для $p_{2}^{*} = 100$ и $p_{2}^{*} = 1$ при $T_{2}^{*} = 3.8$ и ${{\delta }_{1}} = 0.1$, 1 и 10. Для интерполируемых функций ${{J}_{{M,1}}}(\delta )$ и ${{J}_{{M,2}}}(\delta )$ в (3.17) использованы полиномы Чебышева (m = 8). В последнем столбце табл. 4 приведены значения из [2].

Таблица 4
δ1 $J_{M}^{*}$ × 10–2 при $p_{2}^{*}$ = 100 и $T_{2}^{*}$ = 3.8 $J_{M}^{*}$ при $p_{2}^{*}$ = 1 и $T_{2}^{*}$ = 3.8
(3.18) (3.17) [2] (3.18) (3.17) [2]
0.1 0.8974 0.9712 0.971 0.5871 0.6340 0.6324
1.0 2.8557 3.8326 3.818 0.3389 0.4312 0.4315
10 22.6028 33.2470 32.82 –0.2360 0.1474 0.1496

Учитывая, что $\tau = \sqrt {T{\text{*}}(z{\text{*}})} $ и $T{\text{*}}(z{\text{*}}) = (T_{2}^{*} - 1)z{\text{*}} + (1 + T_{2}^{*}){\text{/}}2$, согласно (3.16) получаем распределение давления $p{\text{*}}(z{\text{*}})$ в канале. На рис. 4 и 5 приведены распределения $p{\text{*}}(z{\text{*}})$ для $p_{2}^{*} = 100$ и $p_{2}^{*} = 1$ при $T_{2}^{*} = 3.8$ и значениях δ1 из табл. 4.

ЗАКЛЮЧЕНИЕ

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

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

  1. Шарипов Ф.М., Селезнев В.Д. Движение разреженных газов в каналах и микроканалах. Екатеринбург. УрО РАН. 2008. 231 с.

  2. Graur I., Sharipov F. Non-isothermal flow of rarefied gas through a long pipe with elliptic cross section // Microfluid Nanofluid. 2009. V. 6. P. 267–275.

  3. Sharipov F. Benchmark problems in rarefied gas dynamics // Vacuum. Special Issue Vacuum Gas Dynamics: Theory, experiments and practical applications. 2012. V. 86. № 11. P. 1697–1700.

  4. Шахов Е.М. Об обобщении релаксационного кинетического уравнения Крука // Изв. АН СССР. МЖГ. 1968. № 5. С. 142–145.

  5. Sharipov F., Seleznev V. Data on internal rarefied gas flows // J. Phys. Chem. Ref. Data. 1998. V. 27. № 3. P. 657–706.

  6. Siewert C.E., Valougeorgis D. An analytical discrete-ordinates solution of the S-model kinetic equations for flow in a cylindrical tube // Journal of Quantitative Spectroscopy & Radiative Transfer. 2002. V. 72. P. 531–550.

  7. Kamphorst C.H., Rodrigues P., Barichello L.B. A Closed-Form Solution of a Kinetic Integral Equation for Rarefied Gas Flow in a Cylindrical Duct // Applied Mathematics. 2014. V. 5. P. 1516–1527.

  8. Шахов Е.М. Течение разреженного газа в трубе конечной длины // ЖВМ и МФ. 2000. Т. 40. № 4. С. 647–655.

  9. Титарев В.А., Шахов Е.М. Концевые эффекты при истечении разреженного газа через длинную трубу в вакуум // Изв. РАН. МЖГ. 2013. Т. 3. С. 146–155.

  10. Titarev V.A., Shakhov E.M. Computational study of rarefied gas flow through a long circular pipe into vacuum // Vacuum. 2012. № 6. P. 1709–1716.

  11. Титарев B.А., Шахов Е.М. Неизотермическое течение газа в длинном канале на основе кинетической S-модели // Ж. вычисл. матем. и матем. физ. 2010. Т. 50. № 12. С. 2246–2260.

  12. Ritos K., Lihnaropoulos Y., Naris S., Valougeorgis D. Study of the thermomolecular pressure difference phenomenon in thermal creep flows through microchannels of triangular and trapezoidal cross sections // 2nd Micro and Nano Flows Conference. Brunel University. West London. UK. 2009.

  13. Mason J., Handscomb D. Chebyshev polynomials. Florida. CRC Press. 2003.

  14. Boyd J.P. Chebyshev and Fourier spectral methods. Second Edition. Mineola. DOVER Publications. 2000.

  15. Germider O.V., Popov V.N. Nonisothermal Free-Molecular Flow of Gas in an Elliptic Channel with a Circular Cylindrical Element Inside // Technical Physics. 2019. V. 64. № 1. P. 19–23.

  16. Kogan M.N. Rarefied Gas Dynamics. New York: Plenum, 1969.

  17. Латышев А.В. Аналитические решения граничных задач для кинетических уравнений. М.: МГОУ, 2004. 286 с.

  18. Гермидер О.В., Попов В.Н. Потоки тепла и массы при неполной аккомодации молекул разреженного газа стенками эллиптического канала // Изв. РАН. МЖГ. 2017. № 5. С. 103–109.

  19. Tohidi E. Application of Chebyshev collocation method for solving two classes of non-classical parabolic PDEs // Ain Shams Engineering Journal. 2015. V. 6. № 230. P. 373–379.

  20. Graur I., Sharipov F. Gas flow through an elliptical tube over the whole range of the gas rarefaction // European Journal of Mechanics B Fluids. 2008. V. 27. P. 335–345.

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