Моделирование зоны возмущения разреженной многокомпонентной низкотемпературной плазмы неподвижным симметричным телом

Полный текст

Физические условия и предварительное обсуждение модели

Ионизированные многокомпонентные среды, включая разреженные, их свойства и взаимодействие с различными телами и системами традиционно привлекают пристальное внимание исследователей. Хотя многие приложения низкотемпературной плазмы относительно недавно возникли в таких областях, как медицина, экология, связь, энергетика, физика твердого тела и технологии [1–10], позиции аэрокосмической тематики в подобных исследованиях остаются весьма значительными [11–22].  

Вместе с тем в задачах теплообмена разреженная низкотемпературная плазма практически не рассматривается, что обусловлено спецификой данного состояния вещества и особенностями его описания. Прежде всего, исследование разреженных сред требует отказа от традиционных для специалистов по теплообмену макроскопических уравнений, которые принципиально неприменимы к газам в непосредственной близости от взаимодействующих с ними поверхностей. Макроскопическое описание предполагает изотропию распределения частиц по скоростям, однако это условие нарушается в кнудсеновском пристеночном слое, где динамику газа можно описывать только на основе кинетического подхода [23, 24]. Кроме того, в ионизованных газах с высокой концентрацией заряженных частиц существенное влияние на динамику среды оказывает электромагнитное поле, описание которого нельзя рассматривать изолированно. Наконец, плазма в общем случае представляет собой настолько сложную и многопараметрическую систему, что для неё практически не существует универсальной и одновременно эффективной математической модели [25, 26]. Таким образом, решение задач теплообмена в такой среде сохраняет актуальность и в настоящее время.  

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

Отметим, что подобное специфическое состояние ионизированного газа помимо возможности его лабораторной реализации наблюдается в космической плазме так называемой переходной зоны — обширной области, простирающейся от внешних границ ионосферы (порядка 2000 км) до плазмопаузы (3–5 радиусов Земли) и далее, где летательные аппараты можно считать квазипокоящимися телами по отношению к частицам [27, 28]. В этой области степень ионизации плазмы возрастает с удалением от Земли от значений порядка $10^{4}$, эффективный свободный пробег ионов превышает $10^{4}$ м, а дебаевский радиус варьируется от нескольких сантиметров до величин, сопоставимых с размерами летательных аппаратов или их ключевых конструкционных элементов. Последнего условия, как правило, достаточно для решения задач диагностики в зондовых экспериментах, проводимых для тел, размеры которых значительно меньше, чем у летательных аппаратов [29]. Магнитное поле Земли на таких расстояниях существенно слабее, чем на поверхности; ларморовский радиус электронов увеличивается с ростом расстояния от нескольких метров до нескольких десятков метров, а у ионов он всегда значительно больше. В таких условиях в области возмущения толщиной до нескольких десятков дебаевских радиусов взаимодействием ионов с магнитным полем Земли можно уверенно пренебречь. 

При наличии в газе свободных электронов абсорбирующие поверхности тел обычно приобретают отрицательный плавающий потенциал [27], обеспечивающий баланс потоков заряженных частиц разного знака на поверхности. Уже в этих условиях положительные ионы преобладают в пристеночной области, поскольку скорости электронов, как правило, значительно выше. Расчет плавающего потенциала представляет собой отдельную обратную задачу, однако его элементарные оценки можно найти в различных источниках (см., например, [30]). Эти оценки показывают, что в рассматриваемой области космической плазмы с температурой до $10^{4}$ K абсолютные значения плавающего потенциала невелики и не превышают нескольких вольт.

При задании на поверхности тела потенциалов, сопоставимых с плавающим, электрические токи в возмущенной зоне оказываются незначительными, а их плотность, согласно расчетам, не превышает десятка мА/м$^{2}$. Следовательно, собственным магнитным полем ионизированного газа в данном случае также можно пренебречь, ограничившись рассмотрением лишь квазистационарного самосогласованного электрического поля.

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

Возможная анизотропия распределения ионов разреженного газа по скоростям в окрестности тела требует их описания на основе бесстолкновительного кинетического уравнения Власова, предложенного в работе [31]. С тех пор накоплен значительный объем исследований, посвященных этому уравнению и связанным с ним задачам, что привело к формированию целого научного направления. Хотя само кинетическое уравнение продолжает носить имя А. А. Власова [25, 32], это название часто используется в контексте различных его модификаций и уточнений.

Основная система уравнений, рассматриваемая в данной работе применительно к плазме с самосогласованным электростатическим полем, называется системой Власова—Пуассона, хотя это наименование встречается и в гравитационной теории. Уравнение Власова успешно применяется для решения физических задач в различных областях, оно всесторонне исследовано математически, для него известны различные методы решения, найдены и изучены важные аналитические решения. Обзоры по данной тематике представлены в относительно недавних работах [33, 34].

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

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

В невозмущенной плазме концентрации и заряды различных частиц удовлетворяют условию квазинейтральности $\sum_{\alpha}e_{\alpha}n_{\alpha}=0$, где индекс $\alpha$ обозначает сорт (тип) частиц с зарядом $e_{\alpha}$ и концентрацией $n_{\alpha}$. Следовательно, на начальном этапе формирования возмущенной зоны макроскопическое электрическое поле создается исключительно внесенным телом и может рассчитываться, как в вакууме.

Для компонентов невозмущенного равновесного газа справедливо максвелловское распределение по скоростям, которое принимается в качестве начального приближения. Однако при моделировании это асимптотическое внешнее граничное условие приходится переносить из «бесконечности» в точки $\boldsymbol{r}_{\infty}$ внешней границы $\partial_{e}V$ конечной зоны возмущения $V$. Если такой перенос для макроскопической плотности заряда вполне корректен (поскольку она действительно становится пренебрежимо малой на расстоянии нескольких дебаевских радиусов от тела), то предположение о сохранении на внешней границе возмущенной зоны распределения частиц по скоростям, идентичного невозмущенной плазме, является некорректным. На границе $\partial_{e}V$ следует учитывать ненулевую среднюю скорость компонент, необходимую для обеспечения непрерывности полного тока проводимости. Таким образом, на внешней границе возмущенной зоны исчезает лишь макроскопическое электрическое поле, тогда как у функций распределения ионов устраняются проявления анизотропии, вызванные воздействием поверхности тела $\partial_{b}V$.

Наконец, систему уравнений задачи целесообразно представить в безразмерной форме. Для перехода к безразмерным величинам удобно использовать естественную систему масштабов, не приводящую к появлению дополнительных коэффициентов в уравнениях. Принимая ионы первого типа в качестве наиболее тяжелых и медленных частиц, для многокомпонентного газа в системе СИ такую систему масштабов можно записать в виде
\[\begin{gather}
M_{n}=n_{1\infty},\quad M_{T}=T_{1\infty}, 
\quad 
M_{L}=R_{D1\infty}=\Bigl( \frac{\varepsilon_{0} k_{B} T_{1\infty}}{e^{2}n_{1\infty}} \Bigr)^{1/2},
\\ 
\tag{1}
M_{\varphi}=\frac{k_{B}T_{1\infty}}{e},\quad 
M_{v}=\Bigl( \frac{2k_{B}T_{1\infty}}{m_{1}} \Bigr)^{1/2},
\\
M_{t}=\frac{M_{L}}{M_{v}},\quad M_{E}=\frac{M_{\varphi}}{M_{L}},\quad M_{f}=\frac{M_{n}}{M_{v}^{3}},
\end{gather}\]
где $M_n$ —  масштаб концентрации частиц; $M_T$ — масштаб температуры; $M_L$  — масштаб длины; $M_\varphi$ — масштаб потенциала  электрического поля; $M_v$ — масштаб скорости; $M_t$ — масштаб времени; $M_E$ — масштаб напряженности электрического поля; $M_f$ — масштаб функции распределения; $n_1$ — концентрация ионов первого типа; $T_1$ — температура ионов первого типа; $R_{D1}$ — дебаевский радиус  для ионов первого типа; $e$ — элементарный заряд;  $k_{B}$ и $\varepsilon_{0}$ — постоянная Больцмана и электрическая постоянная; $m_1$ — масса иона первого типа. Индекс $\infty$ относится к параметрам невозмущенной плазмы.

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

1. Неголономные координаты и математическая формулировка задачи

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

В монографии [35] исследована связь структуры конвективного фазового оператора кинетического уравнения с метрикой пространства, доказана его ковариантность и установлены правила преобразования в криволинейные координаты евклидова пространства. Однако в этой классической работе А. А. Власова рассматривался переход от декартовых фазовых координат лишь к произвольным голономным криволинейным координатам $u_{k}$, $k=1,\dots,6$, для которых почти всюду выполняется условие взаимно однозначного непрерывного соответствия с декартовыми координатами. В монографии показано, что когда функция распределения не зависит от части переменных $u_{k}$, $k=m+1,\dots,6$, кинетическое уравнение принимает вид
\[\begin{equation} 
\tag{2}
\frac{\partial (J^{\ast}f )}{\partial t} + 
\sum_{k=1}^{m} \frac{\partial}{\partial u_{k}}  (\dot{u}_{k}J^{\ast}f ) = 0, 
\quad 
J^{\ast} = \int_{\Omega_{m+1,\dots,6}} J(u_{1},\dots,u_{6})d\Omega_{m+1,\dots,6},
\end{equation}\]
где $J^{\ast}$ — якобиан преобразования, проинтегрированный по области $\Omega_{m+1,\dots,6}$ полного изменения переменных $u_{k}$, $k=m+1,\dots,6$. Уравнение (2) формально совпадает с уравнением неразрывности для декартовых координат, выражая условие баланса вещества с плотностью $J^{\ast}f$ в криволинейном фазовом пространстве.

Однако голономные координаты не всегда обеспечивают минимизацию области решения. В частности, при наличии очевидной угловой и трансверсальной симметрии, характерной для неподвижных сферических и цилиндрических тел, фазовое пространство становится более локализованным при использовании координат $r$, $v$, $\mu=\boldsymbol{r\cdot v}/(rv)$ [29]. Эти координаты являются неголономными, так как их связь с компонентами трансверсальной скорости не является взаимно однозначной, несмотря на непрерывность. Следовательно, переход к указанным неголономным переменным требует дополнительных пояснений.

Во-первых, отметим, что преобразование кинетического уравнения к переменным $r$, $v$, $\mu$ в случае центральной и осевой симметрии задачи при отсутствии поля приводит к исчезновению в фазовом конвективном операторе компонент, содержащих производные по $v$. Следовательно, соответствующие преобразования следует выполнять только при наличии поля.

Во-вторых, поскольку форма конвективного фазового оператора в голономных координатах установлена в работе [35], преобразование координат можно проводить исходя из формы уравнения Власова для сферических и цилиндрических координат, полученных в указанной работе. Обозначив через $Z$ зарядовое число частиц и используя масштабы (1), запишем эти исходные формы для дальнейших преобразований.

В случае цилиндрических координат имеем
\[\begin{multline}
\tag{3}
\frac{Df}{dt} = \frac{\partial f}{\partial t} + u^{\alpha} \frac{\partial f}{\partial x^{\alpha}} - 
\Gamma^{m}_{\alpha\beta}u^{\alpha}u^{\beta} \frac{\partial f}{\partial u^{m}} = {}
\\
{}=
\frac{\partial f}{\partial t} + v_{r}\frac{\partial f}{\partial r} + \frac{v_{\phi}}{r} \frac{\partial f}{\partial \phi} + v_{z}\frac{\partial f}{\partial z} + {}
\\
{}+ \Bigl(\frac{v_{\phi}^{2}}{r} + \frac{ZE_{r}}{2}\Bigr) 
\frac{\partial f}{\partial v_{r}} + 
\Bigl(-\frac{v_{r}v_{\phi}}{r} + \frac{ZE_{\phi}}{2}\Bigr) \frac{\partial f}{\partial v_{\phi}} + 
\dot{v}_{z}\frac{\partial f}{\partial v_{z}}.
\end{multline}\]

В случае сферических координат имеем
\[\begin{multline}
\tag{4}
\frac{Df}{dt} = \frac{\partial f}{\partial t} + u^{\alpha}\frac{\partial f}{\partial x^{\alpha}} - \Gamma^{m}_{\alpha\beta}u^{\alpha}u^{\beta}\frac{\partial f}{\partial u^{m}} = {}
\\
{}= 
\frac{\partial f}{\partial t} + v_{r}\frac{\partial f}{\partial r} + \frac{v_{\theta}}{r}\frac{\partial f}{\partial \theta} + \frac{v_{\phi}}{r\sin\theta}\frac{\partial f}{\partial \phi} +   \Bigl(\frac{v_{\theta}^{2}+v_{\phi}^{2}}{r} + \frac{ZE_{r}}{2}\Bigr)\frac{\partial f}{\partial v_{r}} + {}
\\
{} +
\Bigl(\frac{v_{\phi}^{2}}{r}ctg\theta - \frac{v_{r}v_{\theta}}{r} + \frac{ZE_{\theta}}{2}\Bigr)\frac{\partial f}{\partial v_{\theta}}  
 + \Bigl(-\frac{v_{r}v_{\phi}}{r} - \frac{v_{\theta}v_{\phi}}{r}ctg\theta + \frac{ZE_{\phi}}{2}\Bigr)\frac{\partial f}{\partial v_{\phi}}.
\end{multline}\]

Здесь $\Gamma^{m}_{\alpha\beta}$ — символы Кристоффеля; по повторяющимся верхним и нижним индексам производится суммирование. Формы уравнений (3) и (4) допускают наличие электрического поля произвольной ориентации, а смысл входящих в них переменных ясен из контекста.

Рассмотрим сначала цилиндрический оператор (3), предполагая осевую симметрию в физическом пространстве и дополнительную однородность вдоль оси $Oz$, совпадающей с осью цилиндра. В этом случае поле, функция распределения и конвективный фазовый оператор не зависят от координат $\phi$ и $z$, однако сохраняют зависимость от радиальной и трансверсальной скоростей.

Из (3) получаем
\[\begin{equation}
\tag{5}
\frac{Df}{dt} = \frac{\partial f}{\partial t} + v_{r}\frac{\partial f}{\partial r} + 
\Bigl(\frac{v_{\phi}^{2}}{r} + \frac{ZE_{r}}{2}\Bigr)\frac{\partial f}{\partial v_{r}}
- \frac{v_{r}v_{\phi}}{r}\frac{\partial f}{\partial v_{\phi}}.
\end{equation}\]

Связь между компонентами скорости в цилиндрических и неголономных координатах задается следующими соотношениями:
\[\begin{equation}
\tag{6}
v_{r} = v\mu, \quad v_{\phi} = \pm \sqrt{1-\mu^{2}}, \quad 
v = \sqrt{v_{r}^{2}+v_{\phi}^{2}}, \quad \mu = \frac{v_{r}}{\sqrt{v_{r}^{2}+v_{\phi}^{2}}},
\end{equation}\]
причем именно многозначность одного из выражений в (6) обусловливает неголономность соответствующих преобразований. 

Учитывая, что
\[\begin{equation*}
\frac{\partial}{\partial v_{r}} = \frac{\partial \mu}{\partial v_{r}}\frac{\partial}{\partial \mu} + \frac{\partial v}{\partial v_{r}}\frac{\partial}{\partial v}, 
\quad
\frac{\partial}{\partial v_{\phi}} = \frac{\partial \mu}{\partial v_{\phi}}\frac{\partial}{\partial \mu} + \frac{\partial v}{\partial v_{\phi}}\frac{\partial}{\partial v},
\end{equation*}\]
и вычисляя частные производные, из (6) получаем
\[\begin{equation}\tag{7}
\frac{\partial}{\partial v_{r}} = \frac{1-\mu^{2}}{v}\frac{\partial}{\partial \mu} + \mu\frac{\partial}{\partial v}, 
\quad
\frac{\partial}{\partial v_{\phi}} = \mp\frac{\mu\sqrt{1-\mu^{2}}}{v}\frac{\partial}{\partial \mu} \pm \sqrt{1-\mu^{2}}\frac{\partial}{\partial v}.
\end{equation}\]

Подстановка (7) в (5) дает окончательный вид кинетического фазового оператора в неголономных координатах:
\[\begin{equation}
\tag{8}
\frac{Df}{dt} = \frac{\partial f}{\partial t} + v\mu\frac{\partial f}{\partial r} + 
\frac{ZE_{r}\mu}{2}\frac{\partial f}{\partial v}
+ (1-\mu^{2})\Bigl(\frac{v}{r} + \frac{ZE_{r}}{2v}\Bigr)\frac{\partial f}{\partial \mu}.
\end{equation}\]

В случае неподвижного сферического тела функция распределения частиц зависит лишь от $t$, $r$, $v_r$, $v_\tau$ и обладает трансверсальной симметрией (трансверсальные компоненты векторов обозначены индексом $\tau$). Следовательно, конвективный фазовый оператор в неголономных координатах сохраняет аналогичный вид. Данный факт легко устанавливается посредством преобразований, аналогичных приведённым выше.

С учётом симметрии и уравнения (8) безразмерная система принимает следующий вид:
\[
\frac{\partial f_\alpha}{\partial t} + v\mu\frac{\partial f_\alpha}{\partial r} + 
\frac{ZE_r\mu}{2}\frac{\partial f_\alpha}{\partial v}
+ (1-\mu^2)\Bigl(\frac{v}{r}+\frac{ZE_r}{2v}\Bigr)\frac{\partial f_\alpha}{\partial\mu} = 0,
\]

\[
\beta_\alpha = m_\alpha/T_{\alpha\infty}, \quad
\alpha=1,\ldots, N-1, \quad t > 0, \quad v \in (0,\infty), \quad r \in [r_b, r_\infty],
\]

\[
\frac{\partial^2\varphi}{\partial r^2} + \frac{\gamma}{r}\frac{\partial\varphi}{\partial r} = n_e - \sum_{\alpha=1}^{N-1} Z_\alpha n_\alpha, \quad
E_r = -\frac{\partial\varphi}{\partial r}, \quad
\gamma = \begin{cases}
1, & \text{цилиндр},\\
2, & \text{шар},
\end{cases}
\]

\[
f_\alpha(0, r, v, \mu) = n_{\alpha\infty}(\beta_\alpha/\pi)^{3/2}\exp(-\beta_\alpha v^2), \quad
\varphi(0, r) = \varphi_1 r_1/r,
\]

\[\begin{equation} 
\tag{9}
f_\alpha(t, r_b, v, \mu) = 0, \quad \mu<0, \quad
\varphi(t,r_b) = \varphi_1, \quad
\varphi(t,r_\infty) = 0,
\end{equation}\]

\[
f_\alpha(t,r_e,v,\mu) = n_{\alpha\infty}(\beta_\alpha/\pi)^{3/2}
\exp [-\beta_\alpha (v^2-2vv_{\alpha\infty}+v_{\alpha\infty}^2 ) ],
\]

\[
n_\alpha(t,r) = 2\pi\int_0^\infty v^2\int_{-1}^1 f_\alpha(t,r,v,\mu) d\mu dv,
\]

\[
T_\alpha(t,r) = \frac{2\pi}{n_\alpha}\int_0^\infty v^2\int_{-1}^1 
v_{T,\alpha}^2(t,r,v,\mu)f_\alpha(t,r,v,\mu) d\mu dv,
\]

\[
q_{\alpha,r}(t,r) = \frac{2\pi}{n_\alpha}\int_0^\infty v^3\int_{-1}^1 
\mu v_{T,\alpha}^2(t,r,v,\mu)f_\alpha(t,r,v,\mu)d\mu dv,
\]

\[
\langle v_{\alpha,\tau}\rangle(t,r) = 0, \quad
\langle v_{\alpha,r}\rangle(t,r) = \frac{2\pi}{n_\alpha}\int_0^\infty v^3\int_{-1}^1 
\mu f_\alpha(t,r,v,\mu)d\mu dv,
\]

\[
v_{T,\alpha}^2(t,r,v,\mu) = (v\mu-\langle v_{\alpha,r}(t,r)\rangle)^2 + v^2(1-\mu^2),
\]
где $N$ — количество компонент газа; массы частиц и температуры компонент плазмы приведены в безразмерных единицах. Индекс $\alpha$, как и ранее, обозначает тип частиц (при этом $\alpha=N$ заменён на $e$ для электронов). Индекс $b$ соответствует поверхности тела $\partial_b V$, $\infty$ — внешней границе $\partial_e V$ возмущённой зоны, радиус которой задаётся a priori; $T$ обозначает тепловую скорость частиц. Угловые скобки $\langle{}\cdot{}\rangle$ означают усреднение по ансамблю частиц.

Для электронов можно использовать безразмерное распределение Больцмана
\[\begin{equation} 
\tag{10}
n_{e}(t,r) = n_{e,\infty} \exp\bigl(\varepsilon \varphi(t,r)\bigr)
\end{equation}\]
либо модифицированное распределение [27]:
\[\begin{equation}
\tag{11}
n_{e}(t,r) = \frac12 n_{e,\infty} \exp\bigl(\varepsilon \varphi(t,r)\bigr)  
\biggl\{ 
1 + p + (1 - q)\sqrt{1 - \frac{r_{1}^{2}}{r^{2}}}
\exp\Bigl[\frac{\varepsilon r_{b}^{2}(\varphi(t,r) - \varphi_{b})}{r^{2} - r_{b}^{2}}\Bigr]
\biggr\}, 
\end{equation}\]
где 
\[
p = \Phi\Bigl(\sqrt{\varepsilon(\varphi(t,r) - \varphi_{b})}\Bigr), 
\quad
q = \Phi\biggl(\sqrt{\frac{\varepsilon r^{2}(\varphi(t,r) - \varphi_{b})}{r^{2} - r_{b}^{2}}}\biggr),
\]
$\Phi({}\cdot{})$ — функция ошибок, а $\varepsilon = T_{1,\infty}/T_{e,\infty}$.

Важно отметить, что дивергентная форма (2) кинетического оператора (8) позволяет применять для его аппроксимации как консервативные вычислительные схемы [36], так и имитационные динамические методы, основанные на моделировании переноса вещества в фазовом пространстве [37, 38]. 

Учитывая целесообразность использования дивергентной формы кинетического уравнения в задаче (9), отметим зависимость проинтегрированного якобиана преобразования $J^*$ от геометрии тела:
\[\begin{equation}\tag{12}
J^{\ast} = \begin{cases}
4\pi rv/\sqrt{1 - \mu^{2}}, & \text{цилиндр},\\
8\pi^{2}r^{2}v^{2}, & \text{шар},
\end{cases}
\end{equation}\]
где интегрирование выполняется по всем углам физического пространства, включая азимутальный угол для сферической системы или полярный угол для цилиндрической системы в пространстве скоростей. В обоих случаях пространство скоростей двумерно, что ограничивает контроль лишь радиальной компонентой скорости и модулем её трансверсальной компоненты.

2. Численное решение и некоторые особенности реализации метода

Для решения кинетического уравнения задачи (9) использовался метод крупных частиц с постоянным форм-фактором [38], реализованный на равномерной неоднородной сетке криволинейного фазового пространства $\{r, v, \mu\}$ с шагами $h_{r}$, $h_{v}$, $h_{\mu}$. Метод предполагает, что последовательно рассматриваемые ячейки сетки фазового пространства (крупные частицы) объемом $h_{r}h_{v}h_{\mu}$ (якобиан (12) учитывается в фазовой плотности вещества согласно (2)) сначала перемещаются как жесткие объекты в соответствии со скоростями, рассчитанными для их центров, и шагом метода по времени, определяемым из условия устойчивости. Затем вещество смещенной крупной частицы в соответствии с выбранным форм-фактором перераспределяется между ячейками разностной сетки, с которыми у нее имеется непустое пересечение. При постоянном форм-факторе вещество частицы равномерно распределено по ее объему, а сам метод крупных частиц эквивалентен явной схеме дифференцирования против потока, имеющей первый порядок точности аппроксимации по шагам и условие устойчивости Куранта [36–38]:  
\[\begin{equation}
\tag{13}
\Delta t < \frac{1}{\sqrt{K}}\min_{1 \leqslant k \leqslant K} \Bigl( \frac{h_{k}}{\max |\dot{u}_{k}|} \Bigr),
\end{equation}\]
где $K$ — размерность пространства задачи, а экстремумы определяются по всем ячейкам области решения. С физической точки зрения условие (13) не позволяет крупным частицам в течение одного шага метода по времени сместиться дальше, чем на одну ячейку сетки, по любому из координатных направлений.  

Для расчета макропараметров ионов использовались квадратуры интегралов, входящих в систему (9). При этом проблему вычисления несобственных интегралов по $v \in (0,\infty)$ в (9) можно решить, например, используя адаптивную подвижную по $v$ сетку [29], отслеживающую смещение по этой переменной носителя функции распределения. Тогда бесконечный интервал интегрирования можно заменить подвижным интервалом фиксированной конечной ширины, задаваемой a priori, гарантированно перекрывающим основную часть носителя. Также несобственный интеграл можно представить в виде суммы двух интегралов по единичному промежутку с переменными интегрирования $v$ и $1/v$. Подобный прием, хотя и несколько более ресурсоемкий, успешно применялся автором для расчета несобственных интегралов в задачах переноса излучения [39].  

В данной работе использовался первый подход. Для уменьшения контролируемого диапазона переменной $v$ в условиях смещения носителя функции распределения в фазовом пространстве задачи её сетка была выполнена неоднородной и подвижной по $v$. На каждом временном шаге метода первоначальные положения крупных частиц ассоциировались с центрами фазовых ячеек $r_{i}=r_{b}+h_{r}(i-1/2)$, $\mu_{k}=-1+h_{\mu}(k-1/2)$, $v_{j,k,i}=v_{\min\,k,i}+h_{v}(j-1/2)$, где $i=1,\ldots, N_{r}$, $j=1,\ldots, N_{v}$, $k=1,\ldots, N_{\mu}$. Необходимая подвижность сетки обеспечивалась использованием переменной нижней границы $v=v_{\min\,k,i}$, рассчитываемой для каждого $(r_{i},\mu_{k})$ — «пучка» крупных частиц. Для её вычисления использовался следующий алгоритм.

Определим сеточную функцию $Q_{j,k,i}=Q(r_{i},v_{j,k,i},\mu_{k})$, где $Q(r,v,\mu) =J^{\ast} (r,v,\mu)f(r,v,\mu)h_{r}h_{v}h_{\mu}$, задающую количество частиц в фазовых ячейках. Здесь и далее для краткости опускается временной аргумент. В качестве примера рассмотрим случай шара. Пусть $V(r,h_{r})=4\pi r^{2}h_{r}$ — объём шарового слоя толщины $h_{r}$. Тогда значения $Q_{1,k,i}$ и $Q_{N_{v},k,i}$, полученные после каждого этапа перемещения крупных частиц и перераспределения их содержимого, использовались для определения параметра $\eta_{k,i}$ направления смещения носителя функции распределения, параметра величины сдвига $s_{i,k}$ сетки, минимальной скорости $v_{\min\,k,i}$ и пересчёта сеточной функции $Q$ при сдвиге сетки следующим образом (для краткости несущественные для понимания сути преобразований индексы $i,k$ в последующих соотношениях также опущены):
\[\begin{gather} 
s=\frac{|Q_{N_{v}}-Q_{1}|}{\eta(Q_{2}+Q_{N_{v}})+(1-\eta)(Q_{1}+Q_{N_{v}-1})},\quad  
\eta=\Theta(Q_{N_{v}}-Q_{1}),
\\
\begin{array}{l}
v_{\min}:=\eta(v_{\min}+h_{v}s)+(1-\eta)\max(0,v_{\min}-h_{v}s),
\\
Q_{1}:=(1-s)Q_{1}+s \eta Q_{2} , 
\\ 
Q_{N_{v}}:=(1-s)Q_{N_{v}}+s(1-\eta)Q_{N_{v}-1} ,
\\
Q_{j}:=(1-s)Q_{j}+ s[\eta Q_{j+1} +(1-\eta)Q_{j-1}],\quad j=2 ,\ldots, N_{v}-1. 
\end{array}
\end{gather}\]
Здесь $\Theta({}\cdot{})$ — функция Хэвисайда, $:=$ —  символ присваивания. Введение подвижных сеток давно и эффективно используется при решении подобных кинетических задач и позволяет, в частности, по модулю скорости рассматривать относительно небольшой диапазон $\Delta v$ значений $v$. В проведенных расчетах он нигде не превышал 5 (в единицах тепловой скорости наиболее тяжелых ионов невозмущенной плазмы).

Для оценки величины средней скорости $v_{\infty}$ перетекания ионов через внешнюю границу возмущённой зоны запишем уравнение баланса (неразрывности) полного числа частиц $N$ для предпоследнего радиального слоя толщиной $h_{r}$, середина которого имеет радиальную координату $r_{N_{r}-1}=r_{e}-3h_{r}/2$. В качестве примера вновь рассмотрим случай сферической геометрии. Поскольку объём такого слоя с центральным радиусом $r$ на используемой сетке определяется выражением $V(r, h_{r})=4\pi r^{2}h_{r}$, получаем
\[\begin{multline*} 
0 = \Bigl[ \frac{\partial N}{\partial t} + \nabla \cdot (N\langle\boldsymbol{v}\rangle) \Bigr]_{N_{r}-1}^{k} 
= 
\Bigl[ \frac{\partial N}{\partial t} + \frac{\partial}{\partial r}(N\langle v_{r}\rangle) \Bigr]_{N_{r}-1}^{k} = {}
\\
{}= 
4\pi r_{N_{r}-1}^{2} h_{r} \Bigl( \frac{\partial n}{\partial t} \Bigr)_{N_{r}-1}^{k} + 
4\pi h_{r} \Bigl[ \frac{\partial}{\partial r}(N\langle v_{r}\rangle) \Bigr]_{N_{r}-1}^{k} = {}
\\
{}= 4\pi r_{N_{r}-1}^{2}h_{r} \Bigl( \frac{n^{k} - n^{k-1}}{\Delta t} \Bigr)_{N_{r}-1}  
+ 2\pi \bigl( r_{N_{r}}^{2}n_{N_{r}}^{k}\langle v_{r}\rangle_{N_{r}}^{k} - r_{N_{r}-2}^{2}n_{N_{r}-2}^{k}\langle v_{r}\rangle_{N_{r}-2}^{k} \bigr),
\end{multline*}\]
откуда следует искомая оценка:
\[
v_{\infty}^{k} = \langle v_{r}\rangle_{N_{r}}^{k}  
= \frac{r_{N_{r}-2}^{2}}{r_{N_{r}}^{2}} \frac{n_{N_{r}-2}^{k}}{n_{N_{r}}^{k}} \langle v_{r}\rangle_{N_{r}-2}^{k} - 
2h_{r}\frac{r_{N_{r}-1}^{2}}{r_{N_{r}}^{2}} \frac{1}{n_{N_{r}}^{k}} \Bigl( \frac{n^{k} - n^{k-1}}{\Delta t} \Bigr)_{N_{r}-1}.
\]

В приведённых соотношениях верхние индексы соответствуют временному слою (случаи квадратов величин понятны из контекста), нижние индексы обозначают пространственную координату. Пространственная производная аппроксимируется симметричной схемой, временная производная — методом «дифференцирования назад».

Свойства квазилинейных эллиптических операторов и подходы к численному решению нелинейного уравнения Пуассона с правой частью, содержащей распределения типа (10) и (11), подробно исследованы в работе [40]. В указанной работе помимо классического метода установления применимого для случая экспоненциальной нелинейности (распределение Больцмана (10)), предложен новый метод разбиения области решения, основанный на принципе сжимающих отображений и теореме Банаха о неподвижной точке. Данный метод использует связь между нормой эллиптического оператора перехода в процессе простых итераций и диаметром области решения нелинейной задачи Дирихле. 

Метод разбиения области существенно расширяет класс допустимых функций плотности в уравнении Пуассона, включая распределения вида (11). Оба рассмотренных метода обладают локальным характером и могут быть эффективно использованы для вычисления электрического потенциала при численном решении задачи (9).

3. Результаты моделирования и их обсуждение

Описанный в предыдущем разделе метод решения задачи (9) был реализован в среде программирования Matlab. В данном разделе представлены и анализируются результаты, демонстрирующие как адекватность разработанного вычислительного инструмента, так и эффективность кинетического подхода в целом. Исследование проводилось для трёхкомпонентной плазмы со следующими безразмерными параметрами: $r_{b}=7$, $r_{\infty}=25$, $\varphi_{1}=-6$, $T_{1,2,\infty}=T_{e,\infty}=1$, $Z_{1}=1$, $Z_{2}=-1$, $m_{1,2}=1$, $n_{1,\infty}=1$, $n_{2,\infty}=n_{e,\infty}=0.5$. В расчётах использовалось распределение Больцмана (10) для электронов, причём их температура считалась постоянной во всей расчётной области.

Для численного решения в фазовом пространстве неголономных криволинейных координат $(r, v, \mu)$ применялась равномерная подвижная сетка с размерностями $N_{r}=45$, $N_{v}=30$, $N_{\mu}=20$. Ширина контролируемой области по скоростной переменной составляла $\Delta v=4.5$.

Из рис. 1, a видно, что электрическое поле быстро локализуется в окрестности заряженного тела: радиальная компонента напряжённости $E_{r}$ становится практически постоянной малой отрицательной величиной в большей части расчётной области. Такая конфигурация поля обеспечивает приток положительных ионов из внешней зоны к поверхности тела, компенсирующий их поглощение. При этом, как показано на рис. 1, b большая часть возмущённой зоны становится квазинейтральной, что подтверждает корректность выбора границы расчётной области $r_{\infty}$.

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

Рис. 1. Эволюция профиля радиальной напряженности электрического поля $E_{r}$ в ходе развития возмущенной зоны (a) и стационарные профили концентрации электрических зарядов (b). Здесь и далее нумерация кривых соответствует таблице
[Figure 1. Evolution of the radial electric field intensity $E_{r}$ profile during the development of the disturbed zone (a) and stationary profiles of electric charge concentration (b). Here and elsewhere, the curve numbering corresponds to the table]

Моделирование продемонстрировало, что в стационарном режиме электрический потенциал $\varphi$ изменяется монотонно. Однако данное поведение наблюдается не на всех стадиях процесса. На начальном этапе формирования возмущённой области в пристеночный слой поступает избыточное количество положительных ионов, которые не успевают немедленно достичь поверхности тела. В этот период в возмущённой зоне формируются слой положительного объёмного заряда и соответствующий потенциальный барьер, снижающий интенсивность притока положительных ионов. Наличие данного барьера можно идентифицировать при детальном анализе кривой 1 на рис. 1, a, где сохраняется область с положительными значениями радиальной компоненты напряжённости электрического поля $E_r$.

Поведение частиц с различным знаком заряда в возмущённой зоне имеет принципиальные различия. Особый интерес представляют распределения фазовой плотности $Q=J^*f$ для положительных и отрицательных ионов, приведенные на рис. 2  для стационарного состояния вблизи поверхности тела. Для визуализации использовалось фазовое пространство с голономными координатами — радиальной $v_r$ и трансверсальной $v_\tau$ скоростью, что обеспечивает наглядность интерпретации результатов. Переход к переменной $v_\tau$ осуществлён с учётом осевой симметрии трансверсального движения. 

Следует отметить, что все функции распределения $f_{\pm}(r,v_r,v_\tau)$ обладают трансверсальной симметрией. Однако данный эффект, как для $Q$, так и для $f$, необходимо специально учитывать в графических процедурах, поскольку его отображение в неголономных координатах принципиально невозможно.

Рис. 2. Вид стационарной фазовой плотности $Q=J^{\ast}f$ распределения положительных (a) и отрицательных (b) ионов в пристеночном слое тела
[Figure 2. The type of the stationary phase density $Q=J^{\ast}f$ distribution for positive (a) and negative (b) ions in the near-wall layer of the body]

В процессе развития возмущенной области наряду со смещением носителя плотности распределения, отражающего ускорение положительных ионов (причем его величина возрастает по мере приближения к поверхности тела), происходит формирование и усиление группировок частиц, обладающих относительно высокими трансверсальными скоростями \( v_{\tau} \), что отчетливо видно на рис. 2, a. Как будет показано далее, именно такие частицы обуславливают повышение температуры положительных ионов в тонком приповерхностном слое газа. Данный механизм аналогичен процессу разогрева вещества в аккреционной области, наблюдаемому в окрестности черных дыр. В этом смысле поглощающая все падающие частицы сфера демонстрирует сходное поведение.

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

Следует подчеркнуть, что на рисунках представлена именно фазовая плотность \( Q = J^{*}f \) распределения частиц, а не их функция распределения, на которой, в частности, отсутствует видимое разделение потоков падающих на поверхность и отраженных отрицательных ионов. Как видно, кинетическое распределение частиц в пристеночной области является существенно неравновесным и принципиально различается для частиц с разным знаком заряда.

Указанные процессы в фазовом пространстве подтверждаются эволюцией профилей температуры и радиальной плотности тепловых потоков ионов разного знака, представленных на рис. 3.  

Рис. 3. Поведение радиальных профилей тепловых характеристик возмущенной зоны в ходе ее развития: a и b — температура $T_{+}$ и радиальная плотность теплового потока положительных ионов $q_{r+}$; c и d — температура $T_{-}$ и радиальная плотность теплового потока отрицательных ионов $q_{r-}$
[Figure 3. Behavior of radial profiles of thermal characteristics of the disturbed zone during its development: a and b — temperature $T_{+}$ and radial heat flux density of positive ions $q_{r+}$; c and d — temperature $T_{-}$ and radial heat flux density of negative ions $q_{r-}$]

Особого внимания заслуживает поведение профиля температуры положительных ионов \( T_{+} \) в возмущенной области (рис.3, a). Появление в плазме поглощающей отрицательно заряженной поверхности первоначально приводит к «вымыванию» из пристеночной области наиболее энергичных частиц данного типа, вызывая их охлаждение. Однако последующая динамика \( T_{+} \) демонстрирует интенсивный разогрев ионов в приповерхностном слое. Поскольку, как показывают расчеты, этот процесс не сопровождается увеличением радиального ионного тока (являющегося мерой интенсивности радиального движения частиц), рост температуры может быть объяснен лишь накоплением в пристеночной области частиц с высокой угловой скоростью, удерживаемых вблизи сферической поверхности балансом электрических и центробежных сил. Вследствие трансверсальной симметрии задачи для шара средняя орбитальная скорость ионов равна нулю, поэтому накопление положительных частиц с высокой вращательной скоростью приводит к повышению температуры. 

Фазовые процессы также влияют на теплоперенос, определяя вид профилей плотности радиального теплового потока положительных ионов \( q_{r,+} \) в возмущенной зоне (рис. 3, b). В области слабого поля и относительно низкой температуры плотность теплового потока близка к нулю. Основные значения потока сосредоточены в непосредственной близости от поверхности тела и формируются при относительно высоких значениях \( T_{+} \) и \( E \). Как видно из рисунков, пространственная локализация величин \( q_{r,+} \) и \( E \) максимальна у поверхности тела по сравнению с другими характеристиками. При этом усиление теплового потока в приповерхностной области может быть весьма существенным и способно влиять на физические свойства поверхности.

Тепловое поведение отрицательных ионов (рис. 3, c, d), хотя и менее значимо для поверхности тела в рассматриваемом случае, также представляет интерес и связано с фазовыми процессами. Частицы с более высокой энергией первыми покидают возмущенную зону либо замедляются полем. Этот результат ожидаем, однако расчеты показывают, что процесс происходит исключительно быстро, несмотря на равенство масс положительных и отрицательных ионов. Удаление отрицательных ионов носит лавинообразный характер, о чем свидетельствуют наличие и динамика одиночного импульса на радиальных профилях их тока проводимости и теплового потока. Тепловой поток, представленный на рис. 3, d, фактически повторяет поведение тока проводимости и демонстрирует тепловую волну, связанную с выносом отрицательных частиц из возмущенной зоны. Хотя температура отрицательных ионов (рис. 3, c) эволюционирует сложным образом, ее профиль со временем становится монотонным, как и у положительных частиц. Однако, в отличие от положительных ионов, отрицательные частицы вблизи тела, как и ожидалось, охлаждаются.

Заключение

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

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

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

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

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

Конкурирующие интересы. Автор заявляет об отсутствии конфликта интересов в отношении данной публикации. 
Авторский вклад и ответственность. Автор несет полную ответственность за содержание и подготовку окончательной версии рукописи. Окончательная версия рукописи одобрена автором.
Финансирование. Исследование выполнено за счет бюджетного финансирования.
Благодарность. Автор выражает благодарность рецензентам за внимательное рассмотрение работы, ценные замечания и конструктивные предложения.

Об авторах

Валерий Вениаминович Черепанов

Автор, ответственный за переписку.
Email: vvcherepanov@yandex.ru
ORCID iD: 0000-0003-2733-752X
SPIN-код: 3740-3292
Scopus Author ID: 14419630000
ResearcherId: R-7589-2016
https://www.mathnet.ru/rus/person63230

доктор технических наук, доцент; профессор; кафедра физики

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

Дополнительные файлы

Доп. файлы

Действие

1. JATS XML

Скачать

2. Рис. 1. Эволюция профиля радиальной напряженности электрического поля $E_{r}$ в ходе развития возмущенной зоны (a) и стационарные профили концентрации электрических зарядов (b). Здесь и далее нумерация кривых соответствует таблице

Скачать (127KB)

Метаданные

3. Рис. 2. Вид стационарной фазовой плотности $Q=J^{\ast}f$ распределения положительных (a) и отрицательных (b) ионов в пристеночном слое тела

Скачать (270KB)

Метаданные

4. Рис. 3. Поведение радиальных профилей тепловых характеристик возмущенной зоны в ходе ее развития: a и b — температура $T_{+}$ и радиальная плотность теплового потока положительных ионов $q_{r+}$; c и d — температура $T_{-}$ и радиальная плотность теплового потока отрицательных ионов $q_{r-}$

Скачать (343KB)

Метаданные