Журнал вычислительной математики и математической физики, 2023, T. 63, № 12, стр. 1984-1992

ВВЕДЕНИЕ

Наиболее полное решение задачи о струе, выходящей из стационарного плазменного двигателя (СПД), рассматривалось в работах [1] и [2]. Геометрические аспекты постановки этой задачи представлены на фиг. 1.

Из кольцевого отверстия двигателя, который моделируется параллелепипедом AB-CDA₁B₁C₁D₁, выходят в окружающее пространство ионы и нейтралы. Электроны испускаются катодом (он изображен на фиг. 1) и поступают через кольцевое отверстие в двигатель. В названных выше работах приводятся результаты ранее сделанных оценок, которые показывают, что движущаяся среда есть квазинейтральная плазма и для ее адекватного описания необходим кинетический подход.

В [1] и [2] для описания движения ионов и нейтралов использовался кинетический подход, т.е. вводились функции распределения ионов и нейтралов и формулировались кинетические уравнения, откуда они находились. Для описания электронной компоненты использовались ее уравнения движения, откуда с помощью выдвинутой А.И. Морозовым гипотезы “термализованного потенциала” (см. [3]) получалось выражение для электрического поля. Это позволяло замкнуть систему кинетических уравнений для ионов и нейтралов.

В [1] был построен численный метод для решения получившейся системы кинетических уравнений и, используя его, сделать расчеты, результаты которых приведены в [2]. На фиг. 1 параллелепипед KLMNK₁L₁M₁N₁ обозначает область (расчетную), где осуществлялись расчеты. Они производились как вверх, так и вниз по потоку.

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

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

Пусть ${{f}_{i}} = {{f}_{i}}(t,{\mathbf{x}},{\mathbf{\xi }})$ – функция распределения ионов, где ${\mathbf{\xi }}$ – их скорость, а$\;{{D}_{i}}$ ($\;{\mathbf{\xi }} \in {{D}_{i}}$) – скоростное пространство ионов. Соответственно, ${{f}_{e}} = {{f}_{e}}(t,{\mathbf{x}},{\mathbf{v}})$ – функция распределения электронов, а ${\mathbf{v}},\;{{D}_{e}}$ суть скорость и скоростное пространство электронов. Аналогично ${{f}_{n}} = {{f}_{n}}(t,{\mathbf{x}},{\mathbf{w}})$ является функцией распределения нейтралов, тогда как w – скорость нейтралов с ${\mathbf{w}} \in {{D}_{n}}$ их скоростным пространством. Предположим, что так введенные функции распределения подчиняются следующей системе кинетических уравнений:

В (1) электрическое поле ${\mathbf{E}} = \left\{ {{{E}^{i}}} \right\}$ предполагается потенциальным, т.е. ${{E}_{k}} = - \frac{{\partial \varphi }}{{\partial {{x}_{k}}}}$, $k = 1,2,3$, где $\varphi = \varphi (t,{\mathbf{x}})$ есть потенциал электрического поля, для определения которого используется уравнение Максвелла e, как в (1), так и в (2) есть модуль заряда электрона, а m_i и m_e – массы иона и электрона соответственно. Фигурирующие в (2) ${{n}_{e}}$ и ${{n}_{i}}$ суть плотности электронов и ионов соответственно, и определяются следующими формулами: Плотность тока определяется в виде суть плотности ионного и соответственно электронного токов. Величины ${{{\mathbf{u}}}_{i}} = \frac{{{{{\mathbf{j}}}_{i}}}}{{{{n}_{i}}}}$, ${{{\mathbf{u}}}_{e}} = \frac{{{{{\mathbf{j}}}_{e}}}}{{{{n}_{e}}}}$ являются макроскопическими скоростями ионов и электронов. Выражения определяют такие величины, как температуры ионов и электронов. Правые части первого и третьего уравнения (1) суть интегралы столкновений ионов с нейтралами и нейтралами с ионами. В цитируемых выше работах указывалось, что основной вклад в них дает резонансная перезарядка (см. [4]) – взаимодействие, заключающееся в том, что ион отбирает электрон у нейтрала, становясь при этом нейтралом; а нейтрал, отдавая электрон, становится ионом. В [1] использовались для интегралов столкновений следующие выражения ${{J}_{{in}}} = {{\nu }_{{in}}}{{f}_{n}} - {{\nu }_{{ni}}}{{f}_{i}}$, ${{J}_{{ni}}} = {{\nu }_{{ni}}}{{f}_{i}} - {{\nu }_{{in}}}{{f}_{n}}$, где ${{\nu }_{{in}}} = {{n}_{i}}\nu $, ${{\nu }_{{ni}}} = {{n}_{n}}\nu $, а $\nu $ есть функция от макропараметров ионов и нейтралов (ее выражение имеется в [1]). Эти выражения будут использоваться в данной работе. В [2], применяя процедуру феноменологического вывода уравнения Больцмана, были получены более точные выражения для частот столкновений. Как показали тестовые расчеты, результаты численных решений систем кинетических уравнений в диапазоне используемых входных данных, слабо зависели от вида частоты в интегралах столкновений.

Предполагается, что электроны в струе двигаются в бесстолкновительном режиме. Отметим сразу, что в формулах (1), (2) и тех, которые будут приведены ниже, по повторяющимся индексам производится суммирование. Стоящие внизу индексы i, e, n используются исключительно для обозначения принадлежности соответствующей величины к ионам, электронам или нейтралам. Отметим также, что аналогичный подход к описанию плазмы был предложен А.А. Власовым в [5].

Приведем уравнения (1) и (2) к безразмерному виду. За масштаб длины примем полуширину квадрата ABCD – a (см. фиг. 1). За масштабные значения скоростных пространств примем следующие величины:

где ${{U}_{0}}$ – разрядное напряжение (≈300 В), $T_{e}^{0},\;T_{n}^{0}$ – суть характерные значения температуры электронов и нейтралов соответственно. Так как ${{m}_{e}} \ll {{m}_{i}}$, то, как правило, ${{{v}}_{0}} > {{\xi }_{0}}$, ${{t}_{0}} = \frac{a}{{{{\xi }_{0}}}}$, хотя в нашем случае $\frac{{{{\xi }_{0}}}}{{{{{v}}_{0}}}} \approx {{10}^{{ - 2}}}$. Такой выбор временного масштаба связан с тем, что в задаче о струе, выходящей из СПД, основной интерес представляет поведение ионов. Электрическое поле в плазме (самосогласованное) определяется движением ее компонент, поэтому за его масштаб принимается величина ${{\varphi }_{0}} = ~~\frac{{kT_{e}^{{\text{0}}}}}{e}$, ${{E}_{0}} = \frac{{kT_{e}^{{\text{0}}}}}{{ae}}$.$~~$Как обычно, величины где $T_{i}^{0}$ – характерное значение температуры ионов на выходе из отверстия, принимаются за характерные значения соответствующих функций распределения. Здесь ${{n}_{0}}$ – масштабное значение плотности тонов и электронов. Она принимается одной и той же для обоих компонент. В безразмерном виде система (1) будет следующей: В (3) ${\text{S}}{{{\text{h}}}_{e}} = \frac{{{{\xi }_{0}}}}{{{{\operatorname{v} }_{0}}}} = \sqrt {\frac{{{{m}_{e}}{{U}_{0}}}}{{{{m}_{i}}T_{e}^{0}}}} < 1$ – есть аналог числа Струхаля для электронов, а ${\text{S}}{{{\text{h}}}_{n}} = \frac{{{{\xi }_{0}}}}{{{{w}_{0}}}}$ для нейтралов, $B = \frac{{eT_{e}^{0}}}{{2{{U}_{0}}}} < 1$. ${\text{K}}{{{\text{n}}}_{{in}}}$, ${\text{K}}{{{\text{n}}}_{{ni}}}$ – числа Кнудсена ион-нейтрального и нейтрал-ионного взаимодействий, а стоящие перед ними члены суть безразмерные значения интегралов столкновений. Они определяются в [1] и [2].

Безразмерный вид уравнения (2) задается формулой (5)

Причем $\varepsilon = \frac{{{{r}_{d}}}}{a},$ где ${{r}_{d}} = \sqrt {\frac{{kT_{e}^{0}}}{{4\pi {{e}^{2}}{{n}_{0}}}}} $ есть длина Дебайя или дебаевский радиус. В задаче о струе его величина составляет где-то 3 × 10^–3, т.е. является малой.

2. МЕТОД РЕШЕНИЯ

Малость $\varepsilon $ (дебаевского фактора) создает сложность для нахождения численного решения (3), (4). Она заключается в том, что из (4) следует, что, если электрическое поле в струе есть $O(1)$, то

Соотношение (5) называется условием квазинейтральности. Его обычно принимают за определение плазмы – как среды, состоящей из заряженных частиц, чья объемная плотность заряда равна нулю (как в магнитной гидродинамике (см. [6])) или “близка “ к нулю.

Экспериментальные данные свидетельствуют, что в плазме, если исключить дебаевские слои, нет сильных электрических полей $\left( {O\left( {\frac{1}{{{{\varepsilon }^{2}}}}} \right)} \right)$. Тогда из уравнения Максвелла (4) не определяется электрическое поле (оно дает только условие квазинейтральности), и поэтому возникает задача определения электрического поля в ходе совместного решения макроскопических уравнений сохранения, описывающих плазменную среду, и уравнений Максвелла.

Упомянутая выше, гипотеза “термализованного потенциала” является частным случаем способа определения электрического поля, когда предполагают, что $\varphi = f({{n}^{i}})$ (см. [7]).

В некоторых работах проблему определения электрического поля решают так: полагают, что ${{n}_{i}} = {{n}_{e}}$. Подставляя это соотношение в (4), получают, что $\Delta \varphi = 0$. Откуда и получают электрическое поле. В основе такого подхода лежит непонимание того, что (5) есть не изначально заданный факт, а есть результат движения плазменной среды, подчиняющейся уравнениям сохранения и уравнениям Максвелла, в ходе которого получается ${{n}^{i}} - {{n}^{e}} = {{\varepsilon }^{2}}g(t,{\mathbf{x}})$. Подставляя этот результат в (4), получим $\Delta \varphi = g(t,{\mathbf{x}})$. Это соотношение и будет определять электрическое поле в плазме, которое называется самосогласованным.

Целью данной работы является получение конкретного выражения $g(t,{\mathbf{x}})$ в написанном выше соотношении. Из системы кинетических уравнений (1) получим уравнение неразрывности и уравнения сохранения импульса для электронной и ионной компонент. Для этого проинтегрируем первое и второе уравнение (1) каждое по своему пространству скоростей, а затем умножим первое ${\mathbf{\xi }}$, второе на v, и снова проинтегрируем каждое по своему скоростному пространству (см. [8]). Получим

где

Обозначим через N = (${{n}^{i}} - {{n}^{e}}$) и вычтем из уравнения неразрывности для ионов уравнение неразрывности для электронов. Получим $~\frac{{\partial N}}{{\partial t}} + \frac{{\partial {{j}^{k}}}}{{\partial {{x}^{k}}}} = 0$. Вычтем из уравнений сохранения импульса ионов соответствующие уравнения для электронов. После соответствующей манипуляции, которая заключается в переносе соответствующих членов из одной части уравнения в другую и вычитанию величины ${{n}_{n}}\nu j_{e}^{s}$ из обеих частей уравнения, получим

Возьмем дивергенцию от левой и правой частей (7). Воспользовавшись тем, что $\operatorname{div} {\mathbf{E}} = 4\pi eN$, а $\operatorname{div} {\mathbf{\vec {j}}} = - \frac{{\partial N}}{{\partial t}}$, после переноса в правую часть некоторых членов получим где Приведем полученную выше формулу к безразмерному виду. Будем использовать тот факт, что фигурирующая в (1) величина ${{n}_{n}}\nu $ есть (см. [1]) ${{n}_{n}}\nu = n_{n}^{0}{{\xi }_{0}}\sigma {{\bar {n}}_{n}}\bar {\nu }$, где $\sigma $ есть сечение столкновений резонансной перезарядки (в расчетах $\sigma = {{10}^{{ - 14}}}$ см²), а величины с чертой наверху обозначают безразмерные переменные. Чтобы иметь возможность анализа наиболее быстро протекающих процессов, происходящих в плазме, за временной масштаб примем самое короткое время – ${{\bar {t}}_{0}} = \frac{a}{{{{\operatorname{v} }_{0}}}}$. Рассмотрим подробнее переход к безразмерным переменным в левой части уравнения. Имеем где ${\text{K}}{{{\text{n}}}_{d}}$ есть число Кнудсена, определенное по дебаевскому радиусу ${{r}_{d}}$: Если, как это принято, опустить штрихи в безразмерных переменных и индекс у величины ${{\bar {n}}^{d}}$, то приведенное выше уравнение в безразмерном виде будет иметь следующий вид: где В последней формуле ${\mathbf{j}} = {\text{S}}{{{\text{h}}}_{e}}{{{\mathbf{j}}}_{i}} - {{{\mathbf{j}}}_{e}},$ ${\text{K}}{{{\text{n}}}_{{in}}} = \frac{1}{{n_{n}^{0}\sigma a}}$.

Наличие малого параметра $\varepsilon $ в левой части (8) указывает на то, что макроскопическое движение заряженных компонент плазмы будет происходить в двух временных масштабах: а именно – в режиме “ быстрого времени”, когда $\frac{{\partial N}}{{\partial t}} = O\left( {\frac{1}{\varepsilon }} \right)$, и в масштабе, когда $\frac{{{{\partial }^{2}}N}}{{\partial {{t}^{2}}}} \approx \frac{{\partial N}}{{\partial t}} = O(1)$. Для исследования возникающего движения воспользуемся разработанными в [9] асимптотическими методами.

Введем “быстрое время” формулой $\bar {t} = \frac{{(t - {{t}_{0}})}}{\varepsilon }$. Нетрудно видеть, что по порядку величины это – время, которое требуется электронам, чтобы переместиться на расстояние порядка дебаевского радиуса. Будем искать решение (8) в виде разложения по параметру $\varepsilon $:

Подставим (9) в уравнение Учитывая, что получим

Откуда получим, что

и Пусть $N\left( {{{t}_{0}}} \right) \ne 0$, что означает отсутствие квазинейтральности в момент времени ${{t}_{0}}$. Тогда Не ограничиваясь общностью, положим ${{N}_{1}}(\bar {t},{\mathbf{x}}) = 0$, а ${{N}_{2}}(\bar {t},{\mathbf{x}}) = \frac{2}{{\bar {n}\left( {{{t}_{0}},{\mathbf{x}}} \right)}}Q\left( {{{t}_{0}},{\mathbf{x}}} \right)$. Тогда с точностью до $O({{\varepsilon }^{3}})$ решение (10) записывается в виде Из (12) следует, что если в какой-то момент времени нарушается квазинейтральность, то в режиме быстрого времени квазинейтральность устанавливается или в режиме обычной экспоненциальной релаксации.

При получении решения (9) в режиме времени t решение также раскладывается в ряд по параметру $\varepsilon $. Опуская несложные выкладки, можно получить, что в этом случае

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

Из (13) следует, что, если вдруг в плазме нарушится квазинейтральность, т.е. $N\left( {{{t}_{0}}} \right) \ne O({{\varepsilon }^{2}})$, то в течение времени $\Delta t = (t - {{t}_{0}}) \approx \varepsilon $ возникнут процессы, которые приведут к установлению квазинейтральности.

Полученный результат является обобщением результатов, имеющихся в [7]. Там, когда рассматривался вопрос, связанный с затуханием Ландау, после ряда допущений было получено, что если квазинейтральности нет, то возникают колебания с частотой ${{\omega }^{2}} = \frac{{{{n}^{e}}}}{{4\pi {{e}^{2}}{{m}^{e}}}}$. В (8) именно этот член в безразмерной форме определяет частоту возникающих колебаний. В [7] также отмечается, что возникающие колебания будут затухающими, и на примере диффузионного приближения показывается, что если $N({{t}_{0}}) \ne 0$, то $N(t)$ будет за время $\Delta t = (t - {{t}_{0}}) \gg \varepsilon $ будет приближаться к квазинейтральному режиму. В нашем случае этот результат также получается. Действительно, если в (11) $\frac{{{{\alpha }^{2}}}}{4} - \frac{{\bar {n}}}{2} > 0$, то решение определяется формулой (12), и ${{\alpha }_{2}}$ будет соответствовать полученной в [7] константе в показателе экспоненты. Идея проведенного выше анализа возникла именно тогда, когда авторы ознакомились с проведенным в [7] анализом диффузионного приближения. Собственно, проведенный выше анализ есть обобщение изложенного в [7] метода на систему уравнений сохранения, которые были получены из приведенных кинетических уравнений. Нетрудно видеть, что фигурирующая в (9) $\alpha $ пропорциональна числу столкновений, которые испытают носители заряда на длине Дебайя в единицу времени, поэтому ее величина будет определяться числом ${\text{K}}{{{\text{n}}}_{d}}$.

Следует отметить, что (9) было получено без условия потенциальности электрического поля. Кинетическое уравнение для электронов в случае учета ионизации будет другим. Это приведет к тому, что изменятся выражения для частот столкновений в левой части уравнения, а также $Q(t,{\mathbf{x}})$. Но тип уравнения (8) не изменится, поэтому выше приведенный результат имеет достаточно общий характер.

Можно предложить естественную схему численного решения (3), (4). Пусть в момент времени ${{t}_{i}}$ известны функции распределения и электрическое поле. Методом, описанным в [1], можно найти все функции распределения в момент времени $[{{t}_{i}},{{t}_{i}} + \Delta t]$, а значит, все макропараметры, в том числе $N({{t}_{i}} + \Delta t)$.

Потенциал электрического поля определится тогда из уравнения Пуассона

Так и делается во всех работах, где решалась задача о движении плазмы в канале ускорителя. Только в этих работах $N({{t}_{i}} + \Delta t,{\mathbf{x}})$ находилась методом статистического моделирования (метод Берда). Для получения решения уравнения Пуассона (14) необходимо корректно поставить граничные условия, которые ниоткуда не следуют (речь идет о граничных условиях для потенциала электрического поля в плазме). В упомянутых выше работах для решения этой проблемы приходилось идти на ряд существенно ограничивающих полученные результаты допущений. Например, предполагалось, что электрическое поле где z – длина канала, а φ – полярный угол (движение в канале предполагалось осесимметричным), а до этого делались абсурдные предположения, что моделирование происходит для плазмы с другой диэлектрической проницаемостью. Нужно отметить также, что величина шага интегрирования должна выбираться так, чтобы ${\text{min}}\{ \Delta x,\Delta y,\Delta z\} < {{\varepsilon }^{2}}$. Последнюю проблему удается как-то решить, так как за последнее время существенно возросли возможности ЭВМ. Следует также отметить, что в описанной выше схеме никак не учитывается процесс установления квазинейтральности в плазме, а это может привести к возникновению в численном решении незатухающих плазменных колебаний.

Чтобы ввести в численную схему учет процесса установления квазинейтральности, на промежутке $[{{t}_{i}},{{t}_{i}} + \Delta t]$ будем решать уравнение (8). Для получения его решения запишем явную схему, т.е.

Решение уравнения (15) есть $N({{t}_{i}} + \Delta t,{\mathbf{x}}) = g({{t}_{i}},{\mathbf{x}}) + {{\varepsilon }^{2}}Q({{t}_{i}},{\mathbf{x}})$, где Отсутствие числа Струхаля Sh_e объясняется тем, что шаг $\Delta t$ берется в режиме решения кинетических уравнений.

Подставляя найденное решение в правую часть (14), получим $\Delta \varphi = \frac{{g({{t}_{i}},{\mathbf{x}})}}{{{{\varepsilon }^{2}}}} + Q({{t}_{i}},{\mathbf{x}})$. Откуда находится

В (16) интегрирование ведется по всей счетной области V, а параметр ${{\varphi }_{{oot}}}$ есть значение потенциала электрического поля, в котором находится двигатель. Предполагается, что он известен, и удовлетворяет уравнению $\Delta {{\varphi }_{{oot}}} = 0$. Такой способ действия решает проблему граничных условий, которые должны быть установлены, чтобы получить решение (14). С этим электрическим полем на шаге промежутка $[{{t}_{i}},{{t}_{i}} + \Delta t]$ находится решение системы кинетических уравнений. По найденным функциям распределения находятся все макропараметры в момент времени ${{t}_{i}} + \Delta t$, и процесс, если нужно, повторяется.

ЗАКЛЮЧЕНИЕ

Применение полностью кинетического подхода к исследованию плазменной струи позволило получить систему уравнений сохранения, которым подчиняется это движение. Проведенный анализ уравнений сохранения позволил изучить процесс установления квазинейтральности в движущейся плазменной среде, что позволило построить численный метод совместного решения кинетических уравнений и уравнения Максвелла. Этот метод основывается на том, что некоторые макропараметры находятся из решения уравнений сохранения, замкнутых с помощью соответствующих функций распределения. Такой метод решения кинетических уравнений был разработан в [10] и показал свою эффективность. Подводя итог приведенному выше анализу, надо отметить, что он позволил вскрыть механизм установления квазинейтральности в плазме и установить временные масштабы этого процесса, что является основным результатом данной статьи.

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