Журнал вычислительной математики и математической физики, 2023, T. 63, № 8, стр. 1354-1366
Моделирование эмиссионных процессов в сильных электромагнитных полях
Т. А. Кудряшова 1, *, С. В. Поляков 1, **, Н. И. Тарасов 1, ***
1 ИПМ им. М.В. Келдыша РАН
125047 Москва, Миусская пл., 4, Россия
* E-mail: kudryashova@imamod.ru
** E-mail: polyakov@imamod.ru
*** E-mail: nikita_tarasov@imamod.ru
Поступила в редакцию 25.02.2023
После доработки 14.03.2023
Принята к публикации 28.04.2023
- EDN: WSWCXJ
- DOI: 10.31857/S0044466923080100
Аннотация
Рассмотрена проблема расчета процессов электронной эмиссии с поверхности металлов при сильных электромагнитных полях с учетом релятивистских эффектов. Одним из методов моделирования в данной области является метод частиц, сочетающийся с сеточным расчетом полей на основе уравнений Максвелла. Подобные методики развиваются с 1960-х годов по настоящее время. При этом существующие подходы все еще имеют определенные ограничения. В настоящей работе для аксиально-симметричной геометрии, генерирующей системы, представлена новая численная методика моделирования процессов эмиссии электронов с поверхности металлических катодов. Методика использует представление крупных сглаженных гауссовых частиц и реализует расчеты электромагнитных полей на декартовых пространственных сетках. Программная реализация ориентирована на параллельные вычисления. Целью численных экспериментов было определение параметров электронной эмиссии. В качестве тестовых задач были выбраны диодные и триодные цилиндрические системы. В численных расчетах получены пространственно-временные характеристики релятивистских электронных пучков, порождаемых эмиссионными процессами, в том числе воспроизведен ток Чайлда–Ленгмюра. Разработанная численная методика подтвердила свою корректность и эффективность. Библ. 42. Фиг. 7.
ВВЕДЕНИЕ
Традиционно большой интерес представляют задачи о строении веществ и композитных материалов на их основе (см. [1, 2]). Один из экспериментальных способов исследования новых перспективных материалов основывается на воздействии пучка заряженных частиц на изучаемое вещество. В рамках этой технологии используется интенсивное импульсное воздействие на материал. В связи с этим разрабатываются генераторы сильноточных пучков и приборы на их основе (см., например, [3–6]). Наиболее востребованными являются коротко-импульсные генераторы релятивистских электронных пучков (РЭП) с длительностью импульса в диапазоне микро- и наносекунд (см. [7, 8]), которые не разрушают структуру материала.
Генерация РЭП часто основывается на явлении взрывной эмиссии электронов (ВЭЭ) (см. [9, 10]), которая существенно превосходит по энергии выхода как термоэлектронную (см. [11]), так и полевую автоэлектронную эмиссии (см. [12, 13]). Экспериментальному и теоретическому анализам этого явления посвящена обширная литература (см. [6]). Однако наряду с уже известными закономерностями остаются еще вопросы реализации конкретных приборов на основе ВЭЭ. Они решаются в том числе с помощью методов математического и компьютерного моделирования.
Для математического моделирования процессов, составляющих основу сильноточной электроники, с 70-х годов прошлого века создавались необходимые вычислительные основы (см., например, [14–18]) и соответствующие пакеты программ COMSOL (см. [19]), Ansys (см. [20]), КАРАТ (см. [21]), MEEP (см. [22]) и др. Общий тренд численного анализа данного класса задач базируется на применении сеточных методов контрольных объемов или конечных элементов в пространстве-времени FDTD (Finite Difference Time Domain) к решению уравнений Максвелла (см. [16, 21–24]). В дополнение применяются методы частиц в различных модификациях (см., например, [25–32]).
Современный технический потенциал, увеличение возможностей натурных и компьютерных экспериментов, а также широкие перспективы практического использования результатов фундаментальных исследований привели к бурному развитию сильноточного приборостроения. С внедрением этого оборудования в научных лабораториях и на производстве возник интерес к более детальному изучению эмиссионных процессов. В частности, это связано с использованием мощных электронных и ионных пучков в задачах напыления тонких пленок и установках нанолитографии.
Как известно, общей задачей практического применения эффекта эмиссии является создание источников электромагнитного излучения для технологических и хозяйственных целей. Среди типов источников различают электронные, ионные, оптические, смешанные и специальные. Одним из массово используемых механизмов эмиссии является электронная эмиссия, которая чаще всего принимает следующие формы:
• термоэлектронная эмиссия, при которой дополнительная энергия сообщается электронам в результате нагрева катода;
• фотоэлектронная эмиссия, при которой на поверхность катода воздействует электромагнитное излучение;
• полевая эмиссия, при которой сильное электрическое поле у поверхности катода создает силы, способствующие выходу электронов за его пределы;
• взрывная эмиссия – это электронно-ионный вариант в сильных полях, явление вырывания ионов с поверхности твердого тела под действие потоков электронов.
В ряде случаев различные механизмы эмиссии реализуются одновременно. Например, при действии тлеющего разряда на катод одновременно возникают четыре вида эмиссии: ионно-электронная, фотоэлектронная и эмиссия под действием быстрых невозбужденных атомов. Комбинация вторичной электронной эмиссии и эмиссии, возникающей в результате действия сильного электрического поля в тонком слое диэлектрика или полупроводника, называется эффектом Мальтера.
В настоящей работе исследуется взрывная эмиссия, сопровождаемая образованием плазмы и пучков релятивистских электронов (РЭП). Общей целью работы является разработка вычислительной технологии и комплекса параллельных программ для моделирования эмиссионных процессов взрывного типа, а в перспективе – численный анализ взаимодействия пучков частиц с различными средами и материалами. Технология включает разработку комплексной математической модели, множества численных методов, параллельных алгоритмов, комплекса расчетных программ для вычислительных кластеров и суперкомпьютеров, цифровой платформы для проведения детальных численных экспериментов. Для верификации разработанной численной методики использовались данные из [33–35].
1. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ
В математической модели необходимо иметь минимально два уровня описания – сплошная среда для полей (электрических, магнитных, тепловых) и описание в рамках моделей частиц. Мотивацией такого подхода является то, что модели сплошной среды позволяют рассматривать системы произвольного размера, которые актуальны для приложений, а модели частиц позволяют исследовать тонкие эффекты, связанные со свойствами сред и материалов.
Итоговым выбором для математического описания проблемы на первом этапе исследования явились: полная система уравнений Максвелла для электромагнитных полей, уравнение теплопроводности в твердых материалах, уравнения гидродинамики как альтернатива моделям частиц и описание крупных частиц, дополняемое при необходимости моделями молекулярной динамики.
В настоящей работе использованы уравнения Максвелла и специальная модель крупных сглаженных частиц. Последняя имеет тесные связи с методами частиц в ячейках (PiC – Particles in Cells) (см. [25–27]), крупных частиц (LPM – Large-Particle Method) (см. [28, 29]), облаков в ячейках (CIC – Cloud In Cell) (см. [30]), сглаженных частиц в гидродинамике (SPH, SPE – Smoothed-Particles in Hydrodynamics) и электродинамике (Smoothed-Particles in Electrodynamics) (см. [31, 32]). Поясним эти связи.
Во-первых, рассматриваемая нами проблема эмиссии связана с фактом рождения отдельных заряженных частиц (электронов). Появление этих частиц стимулируется средним полем, отнесенным к ячейке сетки, примыкающей к эмиссионной поверхности. Размер ячейки определяет количество рождаемых частиц. Первоначальная форма частиц также зависит от размера ячейки. В этом смысле применяется метод PiC.
Во-вторых, объединение отдельных заряженных частиц в макрочастицу (облако) и дальнейший анализ ее эволюции описываются уже методом крупных частиц, похожим на LPM. Однако специфика этого метода предполагает чередование эйлеровых и лагранжевых расчетных алгоритмов. Нами же используется несколько иной подход, в рамках которого макрочастицы движутся на фоне неподвижной эйлеровой сетки в условиях сохранения массы и заряда макрочастиц на этой сетке. Также нами предполагается возможность туннелирования одних частиц через другие, что не соответствует LPM, но находится в согласии с волновой природой зарядов.
В-третьих, форма частицы в свободной от границ зоне расчетной области имеет гауссов профиль, выходящий за рамки одной ячейки. Такой подход соответствует методу CIC. Однако при приближении к границам форма макрочастицы может искажаться (в методах CIC она обычно фиксирована). Например, при соприкосновении с металлической поверхностью включается механизм поглощения, и часть макрочастицы исчезает. При соприкосновении с диэлектрической поверхностью макрочастица может растекаться вдоль нее или частично отразиться. Таким образом, здесь присутствует аналогия с методами сглаженных частиц SPH и SPE.
Далее в работе обсуждаются основные компоненты модели, численные подходы для их реализации и результаты модельных расчетов.
1.1. Уравнения Максвелла
Основу модели составляют уравнения Максвелла, которые в СИ вместе с материальными уравнениями имеют вид
(1)
$\begin{gathered} \operatorname{div} {\mathbf{D}} = \rho ,\quad \operatorname{div} {\mathbf{B}} = 0, \hfill \\ \operatorname{rot} {\mathbf{E}} = - \frac{{\partial {\mathbf{B}}}}{{\partial t}},\quad \operatorname{rot} {\mathbf{H}} = {\mathbf{j}} + \frac{{\partial {\mathbf{D}}}}{{\partial t}}, \hfill \\ {\mathbf{D}} = {{\varepsilon }_{a}}{\mathbf{E}},\quad {\mathbf{B}} = {{\mu }_{a}}{\mathbf{H}}. \hfill \\ \end{gathered} $Уравнения (1) рассматриваются в области цилиндра $\Omega $, не занятой катодом ${{\Omega }_{C}}$ и анодом ${{\Omega }_{A}}$, т.е. в ${{\Omega }_{D}} = \Omega {\text{/}}\left( {{{\Omega }_{C}} \cup {{\Omega }_{A}}} \right)$. С учетом материальных уравнений и кусочно-постоянных $\varepsilon $ и $\mu $ из (1) можно получить
(1)'
$\begin{gathered} \operatorname{div} \left( {{{\varepsilon }_{a}}{\mathbf{E}}} \right) = \rho ,\quad \operatorname{div} {\mathbf{B}} = 0, \hfill \\ \operatorname{rot} {\mathbf{E}} = - \frac{{\partial {\mathbf{B}}}}{{\partial t}},\quad \operatorname{rot} \left( {\frac{1}{{{{\mu }_{a}}}}{\mathbf{B}}} \right) = {\mathbf{j}} + \frac{\partial }{{\partial t}}\left( {{{\varepsilon }_{a}}{\mathbf{E}}} \right). \hfill \\ \end{gathered} $1.2. Релятивистская динамика частиц
Как указано выше, нами используется специальный метод крупных частиц применительно к задачам эмиссионной электродинамики. В его рамках вместо классических уравнений неразрывности и импульса сплошной среды используются уравнения релятивистской динамики, записанные для отдельных заряженных частиц:
(2)
$\begin{gathered} \frac{{d{{{\mathbf{r}}}_{{\alpha ,k}}}}}{{dt}} = {{{\mathbf{v}}}_{{\alpha ,k}}},\quad \frac{{d{{{\mathbf{p}}}_{{\alpha ,k}}}}}{{dt}} = {{{\mathbf{F}}}_{{\alpha ,k}}} = {{q}_{{\alpha ,k}}}\left( {{\mathbf{E}} + \left[ {{{{\mathbf{v}}}_{{\alpha ,k}}} \times {\mathbf{B}}} \right]} \right), \hfill \\ {{{\mathbf{p}}}_{{\alpha ,k}}} = {{m}_{\alpha }}{{{\mathbf{v}}}_{{\alpha ,k}}}{{\gamma }_{{\alpha ,k}}},\quad {{\gamma }_{{\alpha ,k}}} = 1{\text{/}}\sqrt {1 - {{{\left( {{{v}_{{\alpha ,k}}}{\text{/}}c} \right)}}^{2}}} ,\quad k = 1,...,{{N}_{\alpha }}, \hfill \\ {{n}_{\alpha }} = \sum\limits_{k = 1}^{{{N}_{\alpha }}} {\delta \left( {{\mathbf{r}} - {{{\mathbf{r}}}_{{\alpha ,k}}}} \right)} ,\quad {{\rho }_{\alpha }} = \sum\limits_{k = 1}^{{{N}_{\alpha }}} {{{q}_{{\alpha ,k}}}\delta \left( {{\mathbf{r}} - {{{\mathbf{r}}}_{{\alpha ,k}}}} \right)} , \hfill \\ {{{\mathbf{v}}}_{\alpha }} = \frac{1}{{{{N}_{\alpha }}}}\sum\limits_{k = 1}^{{{N}_{\alpha }}} {{{{\mathbf{v}}}_{{\alpha ,k}}}} ,\quad {{{\mathbf{j}}}_{\alpha }} = \sum\limits_{k = 1}^{{{N}_{e}}} {{{q}_{{\alpha ,k}}}\delta \left( {{\mathbf{r}} - {{{\mathbf{r}}}_{{\alpha ,k}}}} \right){{{\mathbf{v}}}_{{\alpha ,k}}}} ,\quad {{{\mathbf{p}}}_{\alpha }} = \frac{1}{{{{N}_{\alpha }}}}\sum\limits_{k = 1}^{{{N}_{\alpha }}} {{{{\mathbf{p}}}_{{\alpha ,k}}}} . \hfill \\ \end{gathered} $

Уравнения (2) можно модифицировать в целях описания электродинамики крупных частиц. Для этого объединим несколько частиц сорта $\alpha $ в “крупную” частицу и обозначим число таких частиц как $N_{\alpha }^{{(p)}}$. Тогда вместо (2) получим
(2)'
$\begin{gathered} \frac{{d{\mathbf{r}}_{{\alpha ,l}}^{{(p)}}}}{{dt}} = {\mathbf{v}}_{{\alpha ,l}}^{{(p)}},\quad \frac{{d{\mathbf{p}}_{{\alpha ,l}}^{{(p)}}}}{{dt}} = {\mathbf{F}}_{{\alpha ,l}}^{{(p)}} = q_{{\alpha ,l}}^{{(p)}}\left( {{\mathbf{E}} + \left[ {{\mathbf{v}}_{{\alpha ,l}}^{{(p)}} \times {\mathbf{B}}} \right]} \right), \hfill \\ {\mathbf{p}}_{{\alpha ,l}}^{{(p)}} = m_{{\alpha ,l}}^{{(p)}}{\mathbf{v}}_{{\alpha ,l}}^{{(p)}}\gamma _{{\alpha ,l}}^{{(p)}},\quad \gamma _{{\alpha ,l}}^{{(p)}} = 1{\text{/}}\sqrt {1 - {{{\left( {v_{{\alpha ,l}}^{{(p)}}{\text{/}}c} \right)}}^{2}}} ,\quad m_{{\alpha ,l}}^{{(p)}} = k_{{\alpha ,l}}^{{(p)}}{{m}_{\alpha }},\quad q_{{\alpha ,l}}^{{(p)}} = k_{{\alpha ,l}}^{{(p)}}{{q}_{\alpha }},\quad l = 1,...,N_{\alpha }^{{(p)}}, \hfill \\ {{N}_{\alpha }} = \sum\limits_{l = 1}^{N_{\alpha }^{{(p)}}} {k_{{\alpha ,l}}^{{(p)}}} ,\quad {{n}_{\alpha }} \approx n_{\alpha }^{{(p)}} = \sum\limits_{l = 1}^{N_{\alpha }^{{(p)}}} {k_{{\alpha ,l}}^{{(p)}}\bar {\delta }\left( {{\mathbf{r}} - {\mathbf{r}}_{{\alpha ,l}}^{{(p)}}} \right)} ,\quad {{\rho }_{\alpha }} \approx \rho _{\alpha }^{{(p)}} = \sum\limits_{l = 1}^{N_{\alpha }^{{(p)}}} {q_{{\alpha ,l}}^{{(p)}}\bar {\delta }\left( {{\mathbf{r}} - {\mathbf{r}}_{{\alpha ,l}}^{{(p)}}} \right)} , \hfill \\ q_{{\alpha ,l}}^{{(p)}} = \pm ek_{{\alpha ,l}}^{{(p)}},\quad {{{\mathbf{v}}}_{\alpha }} \approx {\mathbf{v}}_{\alpha }^{{(p)}} = \frac{1}{{N_{\alpha }^{{(p)}}}}\sum\limits_{l = 1}^{N_{\alpha }^{{(p)}}} {{\mathbf{v}}_{{\alpha ,l}}^{{(p)}}} ,\quad {{{\mathbf{j}}}_{\alpha }} \approx {\mathbf{j}}_{\alpha }^{{(p)}} = \sum\limits_{l = 1}^{N_{\alpha }^{{(p)}}} {q_{{\alpha ,l}}^{{(p)}}\bar {\delta }\left( {{\mathbf{r}} - {\mathbf{r}}_{{\alpha ,l}}^{{(p)}}} \right){\mathbf{v}}_{{\alpha ,l}}^{{(p)}}} , \hfill \\ {{{\mathbf{p}}}_{\alpha }} \approx {\mathbf{p}}_{\alpha }^{{(p)}} = \frac{1}{{N_{\alpha }^{{(p)}}}}\sum\limits_{l = 1}^{N_{\alpha }^{{(p)}}} {{\mathbf{p}}_{{\alpha ,l}}^{{(p)}}} ,\quad {{\rho }^{{(p)}}} = \rho _{ - }^{{(p)}} + \rho _{ + }^{{(p)}},\quad {{{\mathbf{j}}}^{{(p)}}} = {\mathbf{j}}_{ - }^{{(p)}} + {\mathbf{j}}_{ + }^{{(p)}}. \hfill \\ \end{gathered} $Здесь и ниже будем считать, что частицы сорта представляют собой облака электронов, частицы сорта
представляют собой скопления положительно заряженных ионов.
1.3. Модификация уравнений Максвелла
Уравнения Максвелла в исходном виде не очень удобны для последующего численного решения сеточными методами. В связи с этим многие исследователи используют принцип суперпозиции полей для выделения отдельных их компонент. Наш подход не является исключением. В его рамках векторы электрического и магнитного полей представим в виде сумм квазистатических и нестационарных компонент:
(3)
${\mathbf{E}} = {\mathbf{\bar {E}}} + {\mathbf{\tilde {E}}},\quad {\mathbf{B}} = {\mathbf{\bar {B}}} + {\mathbf{\tilde {B}}},$(4)
$\operatorname{div} \left( {{{\varepsilon }_{a}}\nabla u} \right) = 0,\quad {\mathbf{\bar {E}}} = - \nabla u;\quad \operatorname{div} \left( {\frac{1}{{{{\mu }_{a}}}}\nabla \phi } \right) = 0,\quad {\mathbf{\bar {B}}} = - \nabla \phi .$(5)
$\begin{gathered} \operatorname{div} \left( {{{\varepsilon }_{a}}{\mathbf{\tilde {E}}}} \right) = {{\rho }^{{(p)}}},\quad \operatorname{div} {\mathbf{\tilde {B}}} = 0, \\ \frac{{\partial {\mathbf{\tilde {B}}}}}{{\partial t}} = - \operatorname{rot} {\mathbf{\tilde {E}}},\quad \operatorname{rot} \left( {\frac{1}{{{{\mu }_{a}}}}{\mathbf{\tilde {B}}}} \right) = {{{\mathbf{j}}}^{{(p)}}} + \frac{\partial }{{\partial t}}\left( {{{\varepsilon }_{a}}{\mathbf{\tilde {E}}}} \right). \\ \end{gathered} $Поле ${\mathbf{\tilde {E}}}$ представим в виде поля зарядов и волн:
Тогда уравнения (5) для слагаемых ${\mathbf{\tilde {E}}}$ примут вид(7)
$\begin{gathered} \operatorname{div} \left( {{{\varepsilon }_{a}}{{{{\mathbf{\tilde {E}}}}}^{{(p)}}}} \right) = {{\rho }^{{(p)}}},\quad \operatorname{div} \left( {{{\varepsilon }_{a}}{{{{\mathbf{\tilde {E}}}}}^{{(w)}}}} \right) = 0, \\ \frac{\partial }{{\partial t}}\left( {{{\varepsilon }_{a}}{{{{\mathbf{\tilde {E}}}}}^{{(w)}}}} \right) = {\text{rot}}\left( {\frac{1}{{{{\mu }_{a}}}}{\mathbf{\tilde {B}}}} \right) - {{{\mathbf{j}}}^{{(p)}}} - \frac{\partial }{{\partial t}}\left( {{{\varepsilon }_{a}}{{{{\mathbf{\tilde {E}}}}}^{{(p)}}}} \right). \\ \end{gathered} $1.4. Граничные и начальные условия
Системы уравнений (3)–(8) и (2)' дополняются граничными и начальными условиями, которые зависят от конкретной конструкции исследуемой технической системы. Будем считать, что в качестве таких систем выступают множества проводников в вакууме при наличии или отсутствии компактных слоев низкотемпературной квазинейтральной плазмы. При необходимости плазма задается системой зарядов, также описываемой в рамках вышеуказанного метода крупных частиц.
1.4.1. Граничные условия для поля. Граничные условия для потенциала электростатического поля задаются следующим образом:
(91)
$u = {{u}_{A}} = 0,\quad {\mathbf{r}} \in \partial {{\Omega }_{A}};\quad u = {{u}_{C}},\quad {\mathbf{r}} \in \partial {{\Omega }_{C}};\quad \left( {{\mathbf{\bar {E}}},{\mathbf{n}}} \right) = 0,\quad {\mathbf{r}} \in \partial {{\Omega }_{F}}.$Магнитостатическое поле не вычисляется, а задается некоторым распределением компонент во всей расчетной области:
(92)
${\mathbf{\bar {B}}} = {{\left( {{{B}_{r}},{{B}_{\varphi }},{{B}_{z}}} \right)}^{{\text{т}}}},\quad {\mathbf{r}} \in \Omega .$Граничные условия для зарядовой части динамического электрического поля имеют следующий вид:
(93)
${{u}^{{(p)}}} = 0,\quad {\mathbf{r}} \in \partial {{\Omega }_{A}} \cup \partial {{\Omega }_{C}};\quad \left( {{{{{\mathbf{\tilde {E}}}}}^{{(p)}}},{\mathbf{n}}} \right) = 0,\quad {\mathbf{r}} \in \partial {{\Omega }_{F}}.$Граничные условия для волновой части нестационарных полей являются следствием из общих граничных условий. Так, в приближении идеальной проводимости на металлических поверхностях последние задаются в виде
(94)
$\left[ {{\mathbf{E}} \times {\mathbf{n}}} \right] = 0,\quad \left( {{\mathbf{B}} \cdot {\mathbf{n}}} \right) = 0,\quad {\mathbf{r}} \in \partial {{\Omega }_{A}} \cup \partial {{\Omega }_{C}}.$Граничные условия на свободных границах имеют вид
(95)
$\left( {{\mathbf{E}} \cdot {\mathbf{n}}} \right) = 0,\quad \left( {{\mathbf{B}} \cdot {\mathbf{n}}} \right) = 0,\quad {\mathbf{r}} \in \partial {{\Omega }_{F}}.$Граничные условия на оси симметрии ($r = 0$) задаются следующим образом:
(96)
$\left[ {{\mathbf{E}} \times {{{\mathbf{l}}}_{z}}} \right] = 0,\quad \left( {{\mathbf{B}} \cdot {{{\mathbf{l}}}_{z}}} \right) = 0,\quad r = \sqrt {{{x}^{2}} + {{y}^{2}}} = 0.$1.4.2. Граничные условия для частиц.Граничные условия на эмиссионной поверхности для вновь появляющихся в результате эмиссии отрицательно заряженных частиц при использовании уравнений (2)' имеют вид
(10)
$k_{{ - ,l}}^{{(p)}} = \left[ {\frac{{\left| {\left( {{{\varepsilon }_{a}}{\mathbf{E}} \cdot {\mathbf{n}}} \right)} \right|}}{{\left| {{{q}_{ - }}} \right|}}\Delta {{S}_{C}}} \right],\quad {\mathbf{v}}_{{\alpha ,l}}^{{(p)}} = {\mathbf{\bar {v}}}_{{\alpha ,l}}^{{(p)}},\quad {\mathbf{p}}_{{\alpha ,l}}^{{(p)}} = {\mathbf{\bar {p}}}_{{\alpha ,l}}^{{(p)}},\quad l = N_{\alpha }^{{(p)}} + 1,...,N_{\alpha }^{{(p)}} + \bar {N}_{\alpha }^{{(p)}},\quad {\mathbf{r}} \in \partial {{\Omega }_{C}}.$Относительно выбранного граничного условия сделаем несколько замечаний.
Во-первых, в настоящей работе рассматривается автоэлектронная эмиссия при сильном самосогласованном электрическом поле в предположении идеальности эмиссионной поверхности и без учета тепловых и магнитных эффектов. Такое допущение справедливо в диапазоне электрических полей меньших, чем поле пробоя некоторых материалов, и при умеренных магнитных полях. В качестве примера реальных технических решений можно привести стальные катоды, покрытые слоем графита (см. [12, 13]).
Во-вторых, в этой части модели реализуется именно метод PiC, поскольку площадь $\Delta {{S}_{C}}$ определяется размерами ячейки используемой расчетной сетки. Эта площадь определяет число эмитированных электронов на каждом шаге по времени. Параметром модели является количество электронов, объединяемых в крупную частицу. В нашем подходе задается количество таких частиц, слетающих с поверхности одного сеточного элемента. Для уменьшения вычислительных затрат оно варьируется в диапазоне 1–10.
Граничные условия свободного выхода частиц из расчетной области реализуются с помощью геометрического анализа их положения и соответствующей процедуры удаления. В результате свободного выхода количество крупных частиц $N_{\alpha }^{{(p)}}$ в рассматриваемом объеме уменьшается.
1.4.3. Начальные условия. Будем считать, что в начальный момент времени выполнены следующие условия:
(11)
$\begin{gathered} {{{{\mathbf{\tilde {E}}}}}^{{(p)}}} = 0,\quad {{{{\mathbf{\tilde {E}}}}}^{{(w)}}} = 0,\quad {\mathbf{\tilde {B}}} = 0, \\ {\mathbf{v}}_{{\alpha ,l}}^{{(p)}} = {\mathbf{\bar {v}}}_{{\alpha ,l}}^{{(p)}},\quad l = 1,...,N_{\alpha }^{{(p)}},\quad \rho _{\alpha }^{{(p)}} = \bar {\rho }_{\alpha }^{{(p)}} \equiv \sum\limits_{l = 1}^{N_{\alpha }^{{(p)}}} {q_{{\alpha ,l}}^{{(p)}}{\mathbf{\bar {v}}}_{{\alpha ,l}}^{{(p)}}} . \\ \end{gathered} $Распределения частиц по пространству в начальный момент времени могут задаваться различными способами: равномерно или в соответствии с каким-либо физическим процессом. В рассматриваемом нами случае либо частицы отсутствуют (случай неразвитой эмиссии; $N_{\alpha }^{{(p)}} = 0$, $\rho _{\alpha }^{{(p)}} = 0$), либо в некоторой подобласти сформирована низкотемпературная квазинейтральная плазма, занимающая часть или весь расчетный объем. В этом случае $N_{\alpha }^{{(p)}} > 0$, а пространственные распределения отрицательно и положительно заряженных частиц обеспечивают квазинейтральность плазмы $\left( {\rho _{ - }^{{(p)}} + \rho _{ + }^{{(p)}} \approx 0} \right)$ в каждой точке расчетного объема.
Распределения скоростей частиц во втором случае выбираются в соответствии с температурной моделью плазмы. В частности, в рамках однотемпературной модели все частицы имеют модули скорости, соответствующие распределению Максвелла для заданной температуры, а направления скоростей обеспечивают равенство нулю суммарного импульса частиц:
В рамках двухтемпературной модели вышесказанное справедливо для каждой из подсистем частиц, и, следовательно,2. ЧИСЛЕННАЯ И ПРОГРАММНАЯ РЕАЛИЗАЦИИ МОДЕЛИ
2.1. Численная реализация
Численная реализация вышеуказанной модели проводилась для случая аксиальной симметрии в координатах (r, z). При этом предполагалось, что выполняется условие независимости полей от азимутальной переменной $\varphi $.
Сеточная модель заряженных частиц сорта $\alpha $ представляет собой аппроксимацию дельта-функции с модифицированным гауссовым профилем:
(12)
$\begin{gathered} \rho _{\alpha }^{{(p)}} = {{A}_{{p,\alpha }}}\exp \left[ { - {{R}^{2}}\left( {r,z;{{r}_{{p,\alpha }}},{{z}_{{p,\alpha }}}} \right){\text{/}}R_{p}^{2}} \right], \hfill \\ R\left( {r,z;{{r}_{{p,\alpha }}},{{z}_{{p,\alpha }}}} \right) = \sqrt {{{{\left( {r - {{r}_{{p,\alpha }}}} \right)}}^{2}} + {{{\left( {z - {{z}_{{p,\alpha }}}} \right)}}^{2}}} , \hfill \\ \end{gathered} $Численная реализация итоговых уравнений и дополнительных условий представляет собой совокупность четырех блоков.
Первый блок связан с построением расчетной сетки во всей области и численным решением уравнения Пуассона для электростатического потенциала $u$ (см. (4)). В качестве сетки используется либо ортогональная декартова сетка, либо криволинейная сетка, топологически эквивалентная декартовой. Аппроксимация уравнения Пуассона основана на известной схеме “крест” (см. [36]), модифицированной для случая (r,z)-геометрии и выбранного типа сеток. Решение дискретного аналога (4) проводится с помощью итераций. Ввиду возможной сингулярности расчетной области в качестве итерационного процесса применяется алгоритм Якоби из [36] с диагональной матрицей перехода.
Второй блок включает решение динамических уравнений Максвелла для расчета магнитного и электрического полей ${\mathbf{\tilde {B}}}$ и ${{{\mathbf{\tilde {E}}}}^{{(w)}}}$. Численная схема в этом случае реализована в рамках метода FDTD (см. [23, 24]). Фактически она объединяет методы расщепления по времени и контрольных объемов на случай согласованных векторных полей. Особенностью представляемой реализации также является использование криволинейных сеток, топологически эквивалентных декартовым. Также подчеркнем, что используемая аппроксимация уравнений Максвелла для магнитной индукции ${\mathbf{\tilde {B}}}$ обеспечивает точное выполнение условия соленоидальности поля $\operatorname{div} {\mathbf{\tilde {B}}} = 0$ на криволинейной сетке.
Третий блок включает процедуру интегрирования уравнений динамики частиц. Для этого используется известная симметричная схема по времени из [36]. Фактически она представляет собой неявную схему Адамса второго порядка точности, записанную в векторной форме.
Четвертый блок связан с решением квазистатического уравнения (5) для компоненты электрического поля ${{{\mathbf{\tilde {E}}}}^{{(p)}}}$ (также с помощью итераций на основе метода Якоби) и позволяет найти суммарные поля ${\mathbf{E}}$ и ${\mathbf{B}}$.
2.2. Программная реализация
На основе разработанной численной модели была создана параллельная программа. При ее реализации использовались языки программирования ANSI C и C++, а также стандарты параллельных вычислений MPI (см. [37]) и OpenMP (см. [38]). Методика распараллеливания базируется на методе разделения областей (DDM – Domain Decomposition Method) (см. [39, 40]) и алгоритмах динамической балансировки загрузки вычислителей (Dynamic Load Balancing) (см. [41, 42]).
Методика распараллеливания на основе DDT связана с разбиением декартовой цилиндрической сетки на компактные домены одинаковой мощности. Во время вычислений каждый домен прикрепляется к конкретному MPI-процессу, который последовательно выполняет все вычисления по четырем блокам динамического алгоритма.
В каждый фиксированный момент времени не все домены могут быть заняты частицами. Вследствие этого распределение сеточных доменов проводится только между MPI-процессами. Когда в конкретный домен попадает некоторое количество частиц, они обрабатываются группой параллельных потоков CPU, организованных с помощью стандарта OpenMP. При этом, если частиц немного, используется один поток. В противном случае добавляется такое число потоков, которое обеспечивает примерное равенство загрузки между MPI-процессами. Поскольку частицы динамически переходят из одного домена в другой (т.е. переходят в другой MPI-процесс), то возникает необходимость корректировки вычислительной нагрузки, которую реализует алгоритм динамической балансировки загрузки.
3. РЕЗУЛЬТАТЫ РАСЧЕТОВ
3.1. Эмиссия с боковой поверхности цилиндрического катода
Первая серия расчетов была посвящена численному анализу эмиссии с поверхности цилиндрического диода. Данная задача в условиях бесконечной протяженности поверхности диода имеет в стационарном варианте аналитическое решение, характеризующееся тем, что радиальная компонента плотности тока убывает пропорционально радиусу (${{j}_{r}} \cdot r = {\text{Const}}$) (см. [33, 34]). В проведенных численных экспериментах удалось воспроизвести это решение путем подбора настроечных параметров, а именно, шагов сеток по пространству и времени и константы эмиссии.
При локализации поверхности эмиссии до конечных размеров аналитическое решение не выражается в квадратурах. Однако в центре эмиссионной поверхности реализуется решение, близкое к аналитическому. В качестве примера приведем распределение частиц в зоне развития РЭП (фиг. 1), распределение потенциала электрического поля (фиг. 2) в момент времени 0.1325 нс, а также зависимость суммарного тока эмиссии от времени (фиг. 3). В расчетах шаги пространственной сетки были равны 0.125 см, шаг по времени – 0.125 пс. Напряжение на катоде было равно 511 кВ. Константа эмиссии ${{c}_{{ks}}}{\text{/}}c_{{ks}}^{0}$ равна $0.001$, где $c_{{ks}}^{0} = \frac{{{{\varepsilon }_{0}}{{E}_{n}}r_{n}^{2}}}{e} \approx 2.824 \times {{10}^{{11}}}$ – константа, получаемая при обезразмеривании граничного условия (10) (${{\varepsilon }_{0}} \approx 8.85 \times {{10}^{{ - 14}}}$ Ф/см – диэлектрическая проницаемость вакуума, ${{E}_{n}} = 511$ кВ/см, ${{r}_{n}} = 1$ см – нормировочные параметры, $e \approx 1.602 \times {{10}^{{ - 19}}}$ Кл – заряд электрона).
Фиг. 1.
Распределение частиц при эмиссии с поверхности цилиндрического диода в момент времени 0.1325 нс.

На фиг. 1 и 2 видна характерная структура РЭП и потенциала, соответствующая в центральном сечении z = 4 см аналитическому решению. Зависимость эмиссионного тока от времени на фиг. 3 иллюстрирует эффект начального скачка (поскольку потенциал был включен мгновенно) и последующий выход тока на более низкий стационарный уровень.
3.2. Эмиссия с торцевой поверхности цилиндрического катода
Вторая серия расчетов была посвящена подбору константы эмиссии в случае торцевого цилиндрического катода (см. фиг. 4). В качестве примера была рассмотрена ситуация, когда эмиссия идет с малого участка торца в целях создания трубчатого РЭП в условиях сильного однородного продольного магнитного поля. Физические параметры задачи: величина ${{B}_{z}}$ равна 1 Тл, напряжение на катоде 511 кВ; радиус катода равен 1 см, длина – 10 см; радиус анода – 2.72 см, длина – 30 см. Эмиссионный процесс запускается с помощью ТЕМ-волны, набегающей слева, с напряженностью поля 511 кВ/см. Расчеты проводились на сетке с шагами 0.05 см, шаг по времени составлял 0.5 пс.
В численных экспериментах варьировалась величина безразмерной константы эмиссии ${{c}_{{ks}}}$ относительно значения $c_{{ks}}^{0}$, составлявшей соответственно 0.004, 0.01, 0.02, 0.04, 0.06, 0.08. Результаты расчетов динамики полного тока эмиссии показаны на фиг. 5, кривые 1–6 соответствуют вышеуказанным значениям параметра ${{c}_{{ks}}}$. Они показывают, что полный ток примерно за 1 нс выходит на насыщение, величина тока пропорциональна ${{c}_{{ks}}}$.
Также величина тока зависит от амплитуды ТЕМ-волны. На фиг. 6 показаны кривые эмиссионного тока для амплитуды 511 кВ/см (кривые 1, 2) и 255.5 кВ/см (кривые 3, 4) для значений ${{c}_{{ks}}} = 0.02,\;0.04,\;0.08,\;0.16$ (кривые 1–4). При этом коэффициент эмиссии выбирается уже из другого диапазона.
Фиг. 6.
Зависимости эмиссионного тока от времени для различных значений амплитуды волны и константы эмиссии.

Заметим, что в натурном эксперименте амплитуда волны является контролируемым параметром. С другой стороны, величина ${{c}_{{ks}}}$ определяется материалом катода. Вносит свой вклад и геометрия эмиттеров. В связи с этим подбор константы эмиссии в численных расчетах подчиняется обычно каким-либо известным экспериментальным данным, например, данным о величине тока Чайлда–Ленгмюра (см. [33–35]).
3.3. Тестирование кода в условиях неоднородного магнитного поля
Третья серия расчетов была посвящена моделированию эмиссии в условиях неоднородного магнитного поля. Данная задача актуальна для управления направлением РЭП. При этом часто требуется обеспечить генерацию постоянного тока РЭП. В качестве примера была выбрана геометрия генерирующей системы, показанная на фиг. 7а. Там же показана итоговая квазистационарная конфигурация РЭП. На фиг. 7б,в приведены распределения радиальной и продольной компонент вектора магнитной индукции. Рассмотрен случай, когда генерирующая система частично помещена внутрь соленоида. Параметры соленоида: радиус 3.5 см, полудлина 12 см, амплитуда магнитного поля в точке генерации РЭП (z = 0, r = 2.3 см) составляла 1 Тл. Пучок генерировался ТЕМ-волной с амплитудой 511 кВ/см, проникающей слева через сечение над катодом. Целью численного эксперимента был подбор константы эмиссии, соответствующей току, наблюдаемому в натурном эксперименте.
Фиг. 7.
Расчетная геометрия и конфигурация РЭП (а), распределения радиальной (б) и продольной (в) компонент стационарного магнитного поля.

Численные расчеты проводились для различных значений константы эмиссии. В результате для реализации постоянного тока РЭП величиной 2 кА была определена константа ${{c}_{{ks}}} = 0.01$. При этом значении была подтверждена численно генерация постоянного тока РЭП в течение всего времени действия постоянного входного электромагнитного импульса.
ЗАКЛЮЧЕНИЕ
Рассмотрена проблема расчета процессов электронной эмиссии с поверхности металлов при сильных электромагнитных полях с учетом релятивистских эффектов. Для аксиально-симметричной геометрии генерирующей системы представлена новая численная методика решения задач эмиссии, сочетающая метод частиц и метод сеток. Методика использует представление крупных сглаженных гауссовых частиц и реализует расчеты электромагнитных полей на декартовых пространственных сетках. Программная реализация выполнена в параллельном варианте с использованием стандартов MPI и OpenMP. В численных экспериментах изучены характеристики электронной эмиссии. В качестве модельных геометрий были выбраны коаксиальный диод и цилиндрические системы типа анод–катод–коллектор. В численных расчетах были получены пространственно-временные характеристики релятивистских электронных пучков, порождаемых эмиссионными процессами. Предложенная численная методика подтвердила свои корректность и эффективность. Дальнейшее направление работы будет связано с рассмотрением технических систем реальной геометрии и уточнением механизма эмиссии на микроскопическом уровне.
Список литературы
Васильев В.В. Механика конструкций из композиционных материалов. М.: Машиностроение, 1988. С. 272.
Кербер М.Л. Полимерные композиционные материалы. Структура. Свойства. Технологии. СПб.: Профессия, 2008. С. 560.
Рухадзе A.A., Богданкевич Л.C., Росинский C.E. Физика сильноточных релятивистских электронных пучков. М.: Атомиздат, 1980. С. 167.
Миллер Р. Введение в физику сильноточных пучков заряженных частиц. М.: Мир, 1984. С. 432.
Бойко В.И., Евстигнеев В.В. Введение в физику взаимодействия сильноточных пучков заряженных частиц с веществом. М.: Энергоатомиздат, 1988. С. 136.
Месяц Г.А. Эктоны. Екатеринбург: УИФ “Наука”, 1993. С. 183.
Диденко А.Н., Юшков Ю.Г. Мощные СВЧ-импульсы наносекундной длительности. М.: Энергоатомиздат, 1984. С. 112.
Воронков С.Н., Лоза О.Т., Стрелков П.С. Ограничение длительности импульса излучения СВЧ генераторов на микросекундных РЭП // Физика плазмы. 1991. Т. 17. Вып. 6. С. 751–760.
Бугаев С.П., Литвинов Е.А., Месяц Г.А., Проскуровский Д.И. Взрывная эмиссия электронов // УФН 1975. Т. 115. С. 101–120.
Месяц Г.А. Взрывная электронная эмиссия. М.: Физматлит, 2011. С. 280.
Херинг К., Никольс М. Термоэлектронная эмиссия. М.: Изд-во иностр. лит., 1950. С. 196.
Шешин Е.П. Структура поверхности и автоэмиссионные свойства углеродных материалов. М.: Изд-во МФТИ, 2001. С. 288.
Иванов О.А., Лобаев М.А., Чернов В.В. и др. Экспериментальное исследование сильноточных катодов на основе алмазных пленок в составе мощного компрессора сверхвысокочастотных импульсов // Изв. вузов. Радиофизика. 2014. Т. LVII. № 10. С. 797–806.
Митра Р. (ред.) Вычислительные методы в электродинамике. М.: Мир, 1977. С. 485.
Birdsall C.K., Langdon A.B. Plasma Physics via Computer Simulation. New-York, McGraw-Hill book, 1985. P. 479.
Taflove Allen, Hagness Susan C. Computational Electrodynamics. The Finite-Difference Time-Domain Method. Third Ed. Artech House. 2005. P. 1038.
Inan U.S., Marshall R.A. Numerical Electromagnetics. The FDTD Method. Edinburgh, Cambridge (UK), Cambridge Univ. Press, 2011. P. 406.
Григорьев А.Д. Методы вычислительной электродинамики. М.: Физматлит, 2012. С. 432.
Программное обеспечение COMSOL Multiphysics URL: https://www.comsol.ru/comsol-multiphysics
Официальный сайт компании ANSYS Inc. URL: http://www.ansys.com/
Tarakanov V.P. User’s Manual for Code KARAT. Springfield, VA: Berkeley Res. VA, 1992. C. 262.
Программное обеспечение MEEP. URL: https://meep.readthedocs.io/en/latest/Materials/
Kane Yee. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media // IEEE Transact. Anten. Propagat. 1966. V. 14. №. 3. P. 302–307.
Benford J., Swegle J., Schamiloglu E. High Power Microwaves. Taylor & Francis, New York, 2nd ed. 2007. P. 1–12.
Харлоу Ф.Х. Численный метод частиц в ячейках для задач гидродинамики. Вычислительные методы в гидродинамике. М.: Мир, 1967. С. 316–342.
Дьяченко В.Ф. О расчетах задач бесстолкновительной плазмы // Ж. вычисл. матем. и матем. физ. 1985. Т. 25. № 4. С. 622–627.
Григорьев Ю.Н., Вшивков В.А., Федорук М.П. Численное моделирование методами частиц в ячейках. Новосибирск: Изд-во СО РАН, 2004. С. 358.
Белоцерковский О.М., Давыдов Ю.М. Метод крупных частиц в газовой динамике. М.: Наука, 1982. С. 370.
Садин Д.В. Эффективная реализация гибридного метода крупных частиц // Матем. моделирование. 2022. Т. 34. № 4. С. 113–127.
Birsdall C.K., Fuss D. Clouds-in-clouds, clouds-in-cells physics for many-body plasma simulation // J. Comp. Phys. 1969. V. 3. Iss. 4. P. 494–511.
Jianguo Wang, Dianhui Zhang, Chunliang Liu, Yongdong Li, Yue Wang, Hongguang Wang, Hailiang Qiao, Xiaoze Li. UNIPIC code for simulations of high power microwave devices // Phys. Plasm. 2009. № 16. P. 1–11.
Monaghan J.J. An introduction to SPH. // Comp. Phys. Comm. 1988. V. 48. P. 88–96.
Добрецов Л.Н. Электронная и ионная эмиссия. М.-Л., 1952. С. 312.
Райзер Ю.П. Физика газового разряда. М.: Наука, 1992. С. 536.
Lisovskiy V., Yegorenkov V. Validating the collision-dominated Child-Langmuir law for a dc discharge cathode sheath in an undergraduate laboratory // Europ. J. Phys. 2009. V. 30. № 6. P. 1345–1351.
Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1989. С. 432.
Official documentation and manuals on MPI. [Online]. Available from: http://mpi-forum.org/
Official documentation and manuals on OpenMP. [Online]. Available from: http://www.openmp.org, http://www.llnl.gov/computing/tutorials/openMP
Smith B.F. Domain Decomposition Methods for Partial Differential Equations / In: Keyes, D.E., Sameh, A., Venkatakrishnan, V. (eds) Parallel Numerical Algorithms. ICASE/LaRC Interdisciplinary Series in Science and Engineering, V. 4. Springer, Dordrecht, 1997. P. 225–243.
Dolean V., Jolivet P., Nataf F. An Introduction to Domain Decomposition Methods: algorithms, theory and parallel implementation. Master. France. 2015. P. 289. https://hal.science/cel-01100932v6
Alakeel A.A. Guide to Dynamic Load Balancing in Distributed Computer Systems // Inter. J. Comp. Sci. Network Security (IJCSNS). 2009. V. 10. № 6. P. 153–160.
Sanders P., Mehlhorn K., Dietzfelbinger M., Dementiev R. Sequential and parallel algorithms and data structures: the basic toolbox. Springer Nature, Cham (Switzerland), 2019. P. 516.
Дополнительные материалы отсутствуют.
Инструменты
Журнал вычислительной математики и математической физики