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

ВВЕДЕНИЕ

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

Два основных подхода к моделированию неравновесных сред развиваются в настоящее время. Первый подход основан на построении конечно-разностных или конечно-элементных схем для уравнений сплошной среды [2], [3]. Второй подход заключается в переходе к рассмотрению динамической системы путем построения обобщенных решений уравнений сплошной среды. Следует отметить, что методы моделирования динамических систем сходятся медленно, но обеспечивают эффективное распараллеливание вычислений.

Наиболее известным динамическим методом является метод частиц (см. [4], [5]). Он разработан для моделирования бесстолкновительной плазмы как метод решения квазилинейной системы квазилинейных электродинамических и однородных кинетических уравнений Власова. Рассматривается движение заряженных частиц под действием самосогласованного электромагнитного поля. Переход к моделированию динамической системы заключается в построении обобщенного решения кинетического уравнения посредством так называемой “дельта-подстановки”. Возможность такой подстановки обусловлена тем, что уравнение Власова допускает “микроскопические решения” в виде набора дельта-функций от фазовых переменных системы многих частиц.

Решение некоторых актуальных практических задач потребовало использования метода частиц для моделирования переноса электронов в рассеивающей среде. Электроны высокой энергии ионизируют [6] и возбуждают [7] атомы рассеивающей среды, генерируют тормозное излучение, испытывают упругое рассеяние [8]. При ударной ионизации образуются вторичные электроны [9]. Распространение электронов формирует пространственный заряд и самосогласованное электромагнитное поле [10]. Влияние поля на движение электронов приводит к неустойчивости потока [11]. Процессы переноса и рассеяния интенсивного электронного потока описывает квазилинейное неоднородное кинетическое уравнение с интегралом столкновений. Электронный пучок с интенсивностью, достаточной для формирования самосогласованного поля, создает, например, ускоритель КАЛЬМАР, эксплуатируемый в Национальном исследовательском центре “Курчатовский институт” [12].

Применимость метода частиц для численного решения квазилинейного неоднородного кинетического уравнения с интегралом столкновений необходимо обосновать построением аналога “дельта-подстановки”. Обобщенное решение неоднородного уравнения Власова в однородной рассеивающей среде построено в [13]. Но это не исчерпывает все практически важные задачи. Например, в работе [12] рассматривается эксперимент, в котором пучок электронов образует анодную плазму при рассеянии в преграде из эпоксидной смолы. Как сечения рассеяния электронов, так и коэффициенты уравнений Максвелла для самосогласованного электромагнитного поля испытывают разрывы на границе преграды. В результате испытывают разрывы компоненты силы Лоренца и ядро оператора рассеяния в кинетическом уравнении. Неоднородность имеет место во всех физических экспериментах и практических приложениях.

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

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

Рассмотрим функцию распределения электронов $f = f\left( {t,{\mathbf{r}},{\mathbf{p}}} \right)$ в фазовом пространстве $\left( {{\mathbf{r}},{\mathbf{p}}} \right) = \mathbb{R}_{{\mathbf{r}}}^{3} \times \mathbb{R}_{{\mathbf{p}}}^{3}$ координат ${\mathbf{r}} = \left( {x,y,z} \right)$ и импульсов ${\mathbf{p}} = \left( {{{p}_{x}},{{p}_{y}},{{p}_{z}}} \right)$. Данная функция удовлетворяет уравнению

где $t$ – лабораторное время, ${\mathbf{v}}$ – скорость электрона, ${\mathbf{E}} = {\mathbf{E}}\left( {t,{\mathbf{r}}} \right)$ и ${\mathbf{H}} = {\mathbf{H}}\left( {t,{\mathbf{r}}} \right)$ – напряженности электрического и магнитного полей, ${\text{di}}{{{\text{v}}}_{{\mathbf{r}}}}$ и ${\text{di}}{{{\text{v}}}_{{\mathbf{p}}}}$ – дивергенции в пространствах координат и импульсов, e – заряд электрона, ${\mathbf{\beta }} = {{\mathbf{v}} \mathord{\left/ {\vphantom {{\mathbf{v}} c}} \right. \kern-0em} c}$, $с$ – скорость света в вакууме, ${{\sigma }^{t}}$ – полное макроскопическое сечение рассеяния электронов, $\sigma \left( {{\mathbf{p}},{\mathbf{p}}{\kern 1pt} '} \right)$ – дифференциальное макроскопическое сечение рассеяния электронов, ${\mathbf{p}}{\kern 1pt} '$ и ${\mathbf{p}}$ – импульсы электрона до и после рассеяния, соответственно, ${\mathbf{j}}$ – плотность тока электронов. Внешний источник электронов задается функцией $Q = Q\left( {t,{\mathbf{r}},{\mathbf{p}}} \right)$ в правой части уравнения (1.1).

Задача Коши для уравнений (1.1), (1.2), как показано в [13], имеет решение в пространстве финитных обобщенных функций [14], определенных на основном пространстве финитных бесконечно дифференцируемых функций $\varphi = \varphi \left( {{\mathbf{r}},{\mathbf{p}}} \right)$ в фазовом пространстве $\left( {{\mathbf{r}},{\mathbf{p}}} \right)$. Обобщенные функции зависят от параметра $t \geqslant 0$. Начальные условия предполагаются однородными. Плотность тока ${\mathbf{j}} = {\mathbf{j}}\left( {t,{\mathbf{r}}} \right)$ в уравнениях Максвелла (1.2) определена как действие финитной обобщенной функции ${\mathbf{v}}f\left( {t,{\mathbf{r}},{\mathbf{p}}} \right)$ на бесконечно дифференцируемую функцию $W\left( {{\mathbf{r}},\alpha ,\Delta } \right)$, $\alpha \in \mathbb{R}_{r}^{3}$, $\Delta > 0$, удовлетворяющую условиям:

Функция $W\left( {{\mathbf{r}},{\mathbf{\alpha }},\Delta } \right)$ не удовлетворяет свойству финитности в пространстве $\left( {{\mathbf{r}},{\mathbf{p}}} \right)$. Тем не менее плотности тока и заряда электронов $n = n\left( {t,{\mathbf{r}}} \right)$ определяются корректно именно по причине финитности $f\left( {t,{\mathbf{r}},{\mathbf{p}}} \right)$:

В рамках определения (1.3) плотность тока и, как следствие, напряженности электрического и магнитного полей, вычисляемые решением уравнений Максвелла (1.2), являются бесконечно дифференцируемыми функциями координат. Тогда слагаемое ${\text{di}}{{{\text{v}}}_{{\mathbf{r}}}}\left( {{\mathbf{v}}f} \right) + e{{\operatorname{div} }_{{\mathbf{p}}}}\left[ {\left( {{\mathbf{E + }}\left[ {{\mathbf{\beta ,H}}} \right]} \right)f} \right]$ в левой части уравнения (1.1) есть финитная обобщенная функция [13].

В [13] также показано, что интеграл столкновений в уравнении (1.1) может быть представлен как свертка финитной обобщенной функции $\sigma \left( {{\mathbf{p}},{\mathbf{p}}{\kern 1pt} '} \right)\delta \left( {{\mathbf{r}} - {\mathbf{r}}{\kern 1pt} '} \right)$ и $f = f\left( {t,{\mathbf{r}},{\mathbf{p}}} \right)$. Внешний источник электронов в правой части уравнения (1.1) задается финитной обобщенной функцией $Q = Q\left( {t,{\mathbf{r}},{\mathbf{p}}} \right)$ переменных ${\mathbf{r}},\;{\mathbf{p}}$.

Неоднородность среды приводит к следующим проблемам. Величины $\sigma $ и ${{\sigma }^{t}}$ становятся функциями $\sigma = \sigma \left( {t,{\mathbf{r}},{\mathbf{p}},{\mathbf{p}}{\kern 1pt} '} \right)$, ${{\sigma }^{t}}\left( {t,{\mathbf{r}},p} \right)$. По переменным ${\mathbf{r}}$ эти функции испытывают разрывы I рода. Уравнения (1.2) для электромагнитного поля в вакуумоподобной среде превращаются в уравнения Максвелла для сплошной среды:

где $ \in = \in {\kern 1pt} \left( {\mathbf{r}} \right)$ – диэлектрическая проницаемость, $\xi = \xi \left( {\mathbf{r}} \right)$ – собственная проводимость, $\mu = \mu \left( {\mathbf{r}} \right)$ – магнитная проницаемость.

Возникают сложности при обосновании численного алгоритма решения уравнения (1.1) посредством доказательства его эквивалентности интегральному уравнению:

где функции ${{{\mathbf{r}}}^{s}} = {{{\mathbf{r}}}^{s}}\left( {t,\tilde {t},{\mathbf{\tilde {r}}},{\mathbf{\tilde {p}}}} \right)$, ${{{\mathbf{p}}}^{s}} = {{{\mathbf{p}}}^{s}}\left( {t,\tilde {t},{\mathbf{\tilde {r}}},{\mathbf{\tilde {p}}}} \right)$ являются решениями уравнений движения: с начальными условиями ${{\left. {{{{\mathbf{r}}}^{s}}} \right|}_{{t = \tilde {t}}}} = {\mathbf{\tilde {r}}}$, ${{\left. {{{{\mathbf{p}}}^{s}}} \right|}_{{t = \tilde {t}}}} = {\mathbf{\tilde {p}}}$. Подынтегральное выражение в показателе экспоненты имеет вид: где ${{{\mathbf{r}}}^{{s{\kern 1pt} '}}} = {{{\mathbf{r}}}^{s}}\left( {t{\kern 1pt} ',\tilde {t},{\mathbf{\tilde {r}}},{\mathbf{\tilde {p}}}} \right)$, ${{{\mathbf{p}}}^{{s{\kern 1pt} '}}} = {{{\mathbf{p}}}^{s}}\left( {t{\kern 1pt} ',\tilde {t},{\mathbf{\tilde {r}}},{\mathbf{\tilde {p}}}} \right)$.

Эквивалентность уравнений (1.1) и (1.5) рассмотрена в [13] для бесконечно дифференцируемых в фазовом пространстве $\left( {{\mathbf{r}},{\mathbf{p}}} \right)$ коэффициентов $\sigma = \sigma \left( {t,{\mathbf{r}},{\mathbf{p}},{\mathbf{p}}{\kern 1pt} '} \right),$ ${{\sigma }^{t}}\left( {t,{\mathbf{r}},p} \right)$, ${\mathbf{E}} = {\mathbf{E}}\left( {t,{\mathbf{r}}} \right)$ и ${\mathbf{H}} = {\mathbf{H}}\left( {t,{\mathbf{r}}} \right)$. Рассмотрим действие результата подстановки (1.5) в (1.1) на функцию из основного пространства. Получаем невязку следующей структуры [15]:

Если функции $\sigma \left( {t,{\mathbf{r}},{\mathbf{p}},{\mathbf{p}}{\kern 1pt} '} \right)$, ${{\sigma }^{t}}\left( {t,{\mathbf{r}},p} \right)$, ${\mathbf{E}}\left( {t,{\mathbf{r}}} \right)$ и ${\mathbf{H}}\left( {t,{\mathbf{r}}} \right)$ бесконечно дифференцируемы, то интегрирование по ${\mathbf{r}}$ и ${\mathbf{p}}$ обращает невязку в ноль. В противном случае представленное рассмотрение теряет смысл.

В актуальных задачах зависимость функций $\sigma \left( {t,{\mathbf{r}},{\mathbf{p}},{\mathbf{p}}{\kern 1pt} '} \right)$, ${{\sigma }^{t}}\left( {t,{\mathbf{r}},p} \right)$ от пространственных переменных определяется, главным образом, распределением плотности. Невозмущенное распределение задается кусочно-постоянной функцией координат, т.е. испытывающей разрывы I рода. Если преграда разрушается ионизирующим излучением, то плотность определяется решением уравнений сплошной среды. В этом случае распределение плотности может иметь разрывы I рода и быть дифференцируемым вне точек разрыва. Аналогично, нормальная компонента электрического поля испытывает разрыв I рода на границах сред с разными диэлектрическими проницаемостями, сохраняя дифференцируемость вне областей разрывов.

2. ФУНКЦИЯ РАСПРЕДЕЛЕНИЯ КАК ПЛОТНОСТЬ МЕРЫ

Рассмотрим функцию Хэвисайда $\Theta \left( {{{{\mathbf{r}}}^{s}} - {\mathbf{r}}} \right)\Theta \left( {{{{\mathbf{p}}}^{s}} - {\mathbf{p}}} \right)$, определяющую меру Дирака в фазовом пространстве [16]. Эта мера сконцентрирована в точке $\left( {{{{\mathbf{r}}}^{s}},{{{\mathbf{p}}}^{s}}} \right)$. Она принимает значение 1 для каждого открытого множества евклидова пространства $\left( {{\mathbf{r}},{\mathbf{p}}} \right)$, содержащего точку $\left( {{{{\mathbf{r}}}^{s}},{{{\mathbf{p}}}^{s}}} \right)$. Для открытого множества, которое не содержит точку$\left( {{{{\mathbf{r}}}^{s}},{{{\mathbf{p}}}^{s}}} \right)$, мера Дирака равна нулю. В этой интерпретации соотношение (1.8) является интегралом функции

по мере Дирака, определенной распределением $\Theta \left( {{{{\mathbf{r}}}^{s}} - {\mathbf{r}}} \right)\Theta \left( {{{{\mathbf{p}}}^{s}} - {\mathbf{p}}} \right)$.

Области разрыва коэффициентов являются поверхностями в пространстве $\mathbb{R}_{r}^{3}$, поэтому их мера равна нулю. Ограниченность и финитность $\varphi \left( {{\mathbf{r,p}}} \right)$, а также компонент ее градиента, наряду с ограниченностью остальных функций в подынтегральном выражении, обеспечивают конечность интеграла в бесконечных пределах и его равенство нулю. В результате уравнения (1.1) и (1.5) эквивалентны. В этой интерпретации функцию распределения электронов, являющуюся решением уравнений (1.1) или (1.5), можно рассматривать как плотность меры Дирака, определенной на открытых множествах в пространстве $\left( {{\mathbf{r}},{\mathbf{p}}} \right)$. Например, при отсутствии рассеяния решение уравнения (1.1) имеет вид

Плотность (2.1) соответствует мере Дирака, определяемой соотношением

Из (2.1) и (2.2) следует $dF = fd{\mathbf{r}}d{\mathbf{p}}$.

Требование бесконечной дифференцируемости $W\left( {{\mathbf{r}},{\mathbf{\alpha }},\Delta } \right)$ в такой интерпретации становится излишним. Соотношения (1.3) преобразуются к виду

3. СГЛАЖИВАНИЕ КОЭФФИЦИЕНТОВ

Рассмотрим геометрическую модель рассеивающей среды в трехмерном евклидовом пространстве ${{\mathbb{R}}^{3}}$. Пусть она состоит из М элементов с заданными физическими свойствами. Простые линейно связные области воздуха или вакуума внутри и в окрестности преграды также рассматриваются как отдельные элементы. Поставим в соответствие каждому элементу с номером $m$, $m = 1,...,M$, максимальное открытое множество ${{D}_{m}}$ точек $\mathbb{R}_{r}^{3}$ внутри него. Символом $\partial {{D}_{m}}$ обозначим границу множества ${{D}_{m}}$, $m = 1,...,M$. Пусть ${{\bar {D}}_{m}}$ есть замыкание множества ${{D}_{m}}$. Обозначим ${{D}_{0}} = \mathbb{R}_{r}^{3}{{\backslash }}\bigcup\nolimits_{m = 1}^M {{{{\bar {D}}}_{m}}} $.

Уравнения (1.1), (1.2) содержат неоднородные нестационарные коэффициенты. Будем рассматривать эти коэффициенты как ограниченные функции в пространстве $\mathbb{R}_{r}^{3}$, дифференцируемые на каждом из множеств ${{D}_{m}}$ и испытывающие разрыв I рода на границах $\partial {{D}_{m}}$.

Существует бесконечно дифференцируемая функция ${{h}_{m}}\left( {\mathbf{r}} \right)$ из основного пространства, равная единице для ${\mathbf{r}} \in {{D}_{m}}$, нулю для ${\mathbf{r}} \in \mathbb{R}_{r}^{3}{{\backslash }}{{\bar {D}}_{m}}$ и принимающая значения между нулем и единицей для ${\mathbf{r}} \in {{\bar {D}}_{m}} \cap \left( {\mathbb{R}_{r}^{3}{{\backslash }}{{D}_{m}}} \right)$. Данное утверждение доказывается аналогично лемме, приведенной в [14, с. 62]. Действительно, следуя [14], обозначим через $V_{m}^{\delta }$ объединение открытых шаров радиуса $\delta $ с центрами в точках множества $\mathbb{R}_{r}^{3}{{\backslash }}{{D}_{m}}$. Данное объединение называют дельта-окрестностью множества $\mathbb{R}_{r}^{3}{{\backslash }}{{D}_{m}}$. Число $\delta $ выберем настолько малым, чтобы выполнялось условие $V_{m}^{{2\delta }} \subset \mathbb{R}_{r}^{3}{{\backslash }}{{\bar {D}}_{m}}$. Введем функцию

где $\left| {\mathbf{r}} \right| = r$, ${{C}_{\delta }} = \int_{r \leqslant \delta } {\exp \left\{ { - {{{{\delta }^{2}}} \mathord{\left/ {\vphantom {{{{\delta }^{2}}} {\left( {{{\delta }^{2}} - {{r}^{2}}} \right)}}} \right. \kern-0em} {\left( {{{\delta }^{2}} - {{r}^{2}}} \right)}}} \right\}d{\mathbf{r}}} $, причем $\int_{{{\mathbb{R}}^{3}}} {\psi \left( {{\mathbf{r}},\delta } \right)d{\mathbf{r}}} = 1$.

Рассмотрим функцию

Подынтегральная функция в (3.1) принимает ненулевое значение ${{C}_{\delta }}\exp \left\{ { - {{{{\delta }^{2}}} \mathord{\left/ {\vphantom {{{{\delta }^{2}}} {\left( {{{\delta }^{2}} - {{{\left( {r - r{\kern 1pt} '} \right)}}^{2}}} \right)}}} \right. \kern-0em} {\left( {{{\delta }^{2}} - {{{\left( {r - r{\kern 1pt} '} \right)}}^{2}}} \right)}}} \right\}$ в шаре радиуса $\delta $ с центром в точке ${\mathbf{r}}$. Данный шар является носителем функции $\psi \left( {{\mathbf{r}} - {\mathbf{r}}{\kern 1pt} ',\delta } \right)$. Если ${\mathbf{r}} \in \mathbb{R}_{r}^{3}{{\backslash }}{{D}_{m}}$, то носитель целиком оказывается в $V_{m}^{\delta }$. Интегрирование по $V_{m}^{\delta }$ сводится к интегрированию по носителю и $\int_{{{V}_{\delta }}} {\psi \left( {{\mathbf{r}} - {\mathbf{r}}{\kern 1pt} ',\delta } \right)d{\mathbf{r}}{\kern 1pt} '} = 1$. Функция ${{h}_{m}}\left( {\mathbf{r}} \right)$ в этом случае равна нулю. Если ${\mathbf{r}} \in {{D}_{m}}$, то в силу условия $V_{m}^{{2\delta }} \subset \mathbb{R}_{r}^{3}{{\backslash }}{{\bar {D}}_{m}}$ носитель подынтегральной функции не попадает в область интегрирования $V_{m}^{\delta }$ и интеграл обращается в ноль. Функция $h_{m}^{\delta }\left( {\mathbf{r}} \right)$ в этом случае равна единице. Если ${\mathbf{r}} \in {{\bar {D}}_{m}} \cap \left( {\mathbb{R}_{r}^{3}{{\backslash }}{{D}_{m}}} \right)$, то в область интегрирования попадает часть носителя подынтегральной функции. Тогда функция $h_{m}^{\delta }\left( {\mathbf{r}} \right)$ принимает значение от нуля до единицы. Функция $h_{m}^{\delta }\left( {\mathbf{r}} \right)$ бесконечно дифференцируема также, как и $\psi \left( {{\mathbf{r}},\delta } \right)$.

Рассмотрим функции

Функции (3.2) во всем пространстве $\mathbb{R}_{r}^{3}$ обладают теми же свойствами гладкости, что и функции ${\mathbf{E}}\left( {t,{\mathbf{r}}} \right)$, ${\mathbf{H}}\left( {t,{\mathbf{r}}} \right)$, $\sigma _{{}}^{t}\left( {t,{\mathbf{r}}} \right)$ и $\sigma \left( {t,{\mathbf{r}},{\mathbf{p}},{\mathbf{p}}{\kern 1pt} '} \right)$ на множествах ${{D}_{m}}$. На границах $\partial {{D}_{m}}$ функции (3.2) обращаются в ноль.

Заменим разрывные коэффициенты в уравнениях (1.1), (1.5) соответствующими функциями из набора (3.2):

где функции ${\mathbf{r}}_{\delta }^{s} = {\mathbf{r}}_{\delta }^{s}\left( {t,\tilde {t},{\mathbf{\tilde {r}}},{\mathbf{\tilde {p}}}} \right)$, ${\mathbf{p}}_{\delta }^{s} = {\mathbf{p}}_{\delta }^{s}\left( {t,\tilde {t},{\mathbf{\tilde {r}}},{\mathbf{\tilde {p}}}} \right)$ являются решениями уравнений движения: с начальными условиями ${{\left. {{{{\mathbf{r}}}^{s}}} \right|}_{{t = \tilde {t}}}} = {\mathbf{\tilde {r}}}$, ${{\left. {{{{\mathbf{p}}}^{s}}} \right|}_{{t = \tilde {t}}}} = {\mathbf{\tilde {p}}}$.

Подставляя (3.4) в (3.3), получаем невязку

Обозначим символом $V_{\Gamma }^{\delta }$ объединение открытых шаров радиуса $\delta $ с центрами в точках множества $\Gamma = \bigcup\nolimits_{m = 0}^M {\partial {{D}_{m}}} $. Данное объединение есть дельта-окрестность множества $\Gamma = \bigcup\nolimits_{m = 0}^M {\partial {{D}_{m}}} $. Коэффициенты уравнений (1.1) и (3.3) различаются только для ${\mathbf{r}} \in V_{\Gamma }^{\delta }$. В этой области функции (3.2) при каждом ${\mathbf{r}}$ принимают значения, меньшие по модулю, чем ${\mathbf{E}}\left( {t,{\mathbf{r}}} \right)$, ${\mathbf{H}}\left( {t,{\mathbf{r}}} \right)$, $\sigma _{{}}^{t}\left( {t,{\mathbf{r}}} \right)$и $\sigma \left( {t,{\mathbf{r}},{\mathbf{p}},{\mathbf{p}}{\kern 1pt} '} \right)$, а при ${\mathbf{r}} \in \Gamma $ обращаются в ноль. Из этого очевидно, что различие решений уравнений (1.1) и (3.3) будет максимальным, если при ${\mathbf{r}} \in V_{\Gamma }^{\delta }$ функции (3.2) положить равными нулю. Это будет означать, что в дельта-окрестности границы $V_{\Gamma }^{\delta }$ электроны переносятся свободно, т.е. не поглощаясь, не рассеиваясь и не изменяя импульс под действием электромагнитного поля.

Рассмотрим погрешности, вносимые заменой коэффициентов уравнения (1.1) функциями (3.2). Анализ выполним отдельно для сечений в интеграле столкновений, и напряженностей электромагнитного поля.

Пусть функции $\sigma _{{}}^{t}\left( {t,{\mathbf{r}}} \right)$и $\sigma \left( {t,{\mathbf{r}},{\mathbf{p}},{\mathbf{p}}{\kern 1pt} '} \right)$ в уравнениях (1.1) и (1.5) заменены на функции $\sigma _{\delta }^{t}\left( {t,{\mathbf{r}}} \right)$ и ${{\sigma }^{\delta }}\left( {t,{\mathbf{r}},{\mathbf{p}},{\mathbf{p}}{\kern 1pt} '} \right)$. Решения уравнений (1.5) и (3.4) могут быть построены методом последовательных поколений [17]. Например, для уравнения типа (3.4) данное решение имеет вид

где

Рассмотрим действие разности этих решений

на функцию из основного пространства:

В силу положительности функций $\sigma _{\delta }^{t}$ и $\sigma _{{}}^{\delta }$ для неотрицательных $\varphi \left( {{\mathbf{r}},{\mathbf{p}}} \right)$ справедливы следующие неравенства:

Интегрирование по переменной ${\mathbf{\tilde {r}}}$ в (3.8) ведется по области $V_{\Gamma }^{\delta }$. Сравним (3.8) с действием (3.7) на функцию из основного пространства. Действия решения (3.7) и его погрешности (3.8) будут сравнимы, если объемы носителя основной функции и его пересечения с множества $V_{\Gamma }^{\delta }$ окажутся сравнимыми величинами. Таким образом, погрешность, вносимая заменой коэффициентов уравнения (1.1) функциями (3.5), будет пропорциональна объему области пересечения $V_{\Gamma }^{\delta }$ с носителем функции $W\left( {{\mathbf{r}},{\mathbf{\alpha }},\Delta } \right)$.

Погрешность, вносимую сглаживанием функций ${\mathbf{E}}\left( {t,{\mathbf{r}}} \right)$ и ${\mathbf{H}}\left( {t,{\mathbf{r}}} \right)$, достаточно рассмотреть для уравнений (1.1) и (1.5) без учета столкновений. В соответствии с формулой (2.1) для этого необходимо сравнить величины $\int_0^t {d\tilde {t}} \int {d{\mathbf{\tilde {r}}}} \int {d{\mathbf{\tilde {p}}}} \,\varphi \left( {{{{\mathbf{r}}}^{s}},{{{\mathbf{p}}}^{s}}} \right)\,Q\left( {\tilde {t},{\mathbf{\tilde {r}}},{\mathbf{\tilde {p}}}} \right)$ и $\int_0^t {d\tilde {t}} \int {d{\mathbf{\tilde {r}}}} \int {d{\mathbf{\tilde {p}}}} \,\varphi \left( {{\mathbf{r}}_{\delta }^{s},{\mathbf{p}}_{\delta }^{s}} \right)\,Q\left( {\tilde {t},{\mathbf{\tilde {r}}},{\mathbf{\tilde {p}}}} \right)$.

Для функций $\Delta {\mathbf{r}}_{{}}^{s} \equiv {\mathbf{r}}_{\delta }^{s} - {\mathbf{r}}_{{}}^{s}$, $\Delta {\mathbf{v}}_{{}}^{s} \equiv {\mathbf{v}}_{\delta }^{s} - {\mathbf{v}}_{{}}^{s}$, $\Delta {\mathbf{p}}_{{}}^{s} \equiv {\mathbf{p}}_{\delta }^{s} - {\mathbf{p}}_{{}}^{s}$ справедливы уравнения

Рассмотрим погрешность

Пользуясь гладкостью функции $\varphi \left( {{\mathbf{r}}_{\delta }^{s},{\mathbf{p}}_{\delta }^{s}} \right)$, разложим ее в ряд Тейлора с центром в точке $\left( {{\mathbf{r}}_{{}}^{s},{\mathbf{p}}_{{}}^{s}} \right)$ по всем аргументам, сохраняя нулевой и первый члены разложения:

Рассмотрим следующую оценку:

Тогда

Формула (3.9) показывает, что погрешность определяется отношением длин траекторий в фазовом пространстве и их возмущений, вносимых сглаживанием компонент электромагнитного поля. Положим при решении уравнений (3.9), что ${\mathbf{E}}_{{}}^{\delta }\left( {t,{\mathbf{r}}} \right) = 0$, ${\mathbf{H}}_{{}}^{\delta }\left( {t,{\mathbf{r}}} \right) = 0$ в области $V_{\Gamma }^{\delta }$. Получим, что погрешность пропорциональна сумме двух величин. Первая величина, ${{\left| {\Delta {\mathbf{r}}_{\delta }^{s}} \right|} \mathord{\left/ {\vphantom {{\left| {\Delta {\mathbf{r}}_{\delta }^{s}} \right|} {\left| {{\mathbf{r}}_{{}}^{s}} \right|}}} \right. \kern-0em} {\left| {{\mathbf{r}}_{{}}^{s}} \right|}}$, есть отношение длины пересечения траектории частицы в пространстве $\mathbb{R}_{r}^{3}$ с областью $V_{\Gamma }^{\delta }$ к длине самой траектории. Вторая величина, ${{\left| {\Delta {\mathbf{p}}_{\delta }^{s}} \right|} \mathord{\left/ {\vphantom {{\left| {\Delta {\mathbf{p}}_{\delta }^{s}} \right|} {\left| {{\mathbf{p}}_{{}}^{s}} \right|}}} \right. \kern-0em} {\left| {{\mathbf{p}}_{{}}^{s}} \right|}}$, есть отношение изменения импульса, накапливаемого частицей при прохождении области $V_{\Gamma }^{\delta }$, к характерной величине импульса.

4. ХАРАКТЕРНЫЕ РАЗМЕРЫ ЗАДАЧИ

Полученные оценки имеют простой физический смысл. Для его выявления рассмотрим связь $\Delta {\mathbf{r}}_{\delta }^{s}$ и $\Delta {\mathbf{p}}_{\delta }^{s}$ с характерными частотами задачи (1.1), (1.2).

Обобщенное решение уравнения (1.1) сводит задачу к моделированию движения частиц с координатами ${\mathbf{r}} = {{{\mathbf{r}}}^{s}}\left( {t,\tilde {t},{\mathbf{\tilde {r}}},{\mathbf{\tilde {p}}}} \right)$, ${\mathbf{p}} = {{{\mathbf{p}}}^{s}}\left( {t,\tilde {t},{\mathbf{\tilde {r}}},{\mathbf{\tilde {p}}}} \right)$ и вероятностью $\exp \left\{ { - \int_{\tilde {t}}^t {dt{\kern 1pt} '\sigma _{{}}^{t}{{v}^{{s{\kern 1pt} '}}}} } \right\}$ испытать столкновение к моменту времени t.

Пусть в момент времени t частица оказалась на границе области $V_{m}^{\delta }$, а в момент времени $t + {{t}_{\delta }}$ ее покинула. Импульс частицы изменялся при этом за счет взаимодействия с электромагнитным полем. Времена характерного изменения самосогласованного электромагнитного поля и импульсов генерирующих его заряженных частиц совпадают. Рассмотрим характерные частоты самосогласованного электромагнитного поля, описываемого уравнениями (1.1), (1.2). Для этого построим момент уравнения (1.1), соответствующий плотности тока – действие на функцию ${\mathbf{v}}\,W\left( {{\mathbf{r}},\alpha ,\Delta } \right)$. Для интеграла столкновений используем приближение времени релаксации к равновесному распределению. Пренебрегая градиентом давления в левой части уравнения, получаем

где ${{\tau }_{g}}$ – время релаксации, $\gamma = \sqrt {1 + {{{\left( {{{{{p}_{c}}} \mathord{\left/ {\vphantom {{{{p}_{c}}} {mc}}} \right. \kern-0em} {mc}}} \right)}}^{2}}} $, ${{p}_{c}}$ – характерный импульс частиц, например, средний по распределению $Q\left( {t,{\mathbf{r,p}}} \right),$ m – масса электрона.

Первые два слагаемых (4.1) в совокупности с локально плоским уравнением Максвелла ${{\partial {\mathbf{E}}} \mathord{\left/ {\vphantom {{\partial {\mathbf{E}}} {\partial t}}} \right. \kern-0em} {\partial t}} + 4\pi {\mathbf{j}} = 0$ определяют частоту плазменных колебаний $\omega _{p}^{2} = {{4\pi {{e}^{2}}n} \mathord{\left/ {\vphantom {{4\pi {{e}^{2}}n} {\left( {m{{\gamma }^{3}}} \right)}}} \right. \kern-0em} {\left( {m{{\gamma }^{3}}} \right)}}$. Первое и третье слагаемые (4.1) определяют ларморовскую частоту $\omega _{L}^{{}} = {{eH} \mathord{\left/ {\vphantom {{eH} {\left( {mc\gamma } \right)}}} \right. \kern-0em} {\left( {mc\gamma } \right)}}$.

Первое и четвертое слагаемые (4.1) определяют характерное время релаксации импульса за счет столкновений ${{\tau }_{g}}$. Одновременно второе и четвертое слагаемые определяют характерное время затухания электрического поля за счет компенсации стороннего тока быстрых частиц дрейфом термализованных электронов ${{\tau }_{\sigma }}{{ = {{\gamma }^{3}}m} \mathord{\left/ {\vphantom {{ = {{\gamma }^{3}}m} {{{\tau }_{g}}{{e}^{2}}n}}} \right. \kern-0em} {{{\tau }_{g}}{{e}^{2}}n}}$. Величина ${1 \mathord{\left/ {\vphantom {1 {{{\tau }_{\sigma }}}}} \right. \kern-0em} {{{\tau }_{\sigma }}}}$ есть электронная электропроводность ионизованной среды. Здесь следует отметить, что величина ${{\tau }_{g}}$ может существенно превосходить время жизни до первого столкновения ${{\left( {{{\sigma }^{t}}v} \right)}^{{ - 1}}}$. Все дифференциальные сечения рассеяния электронов устроены так, что наиболее вероятными являются малые изменения импульса при столкновениях [7]–[9]. Поэтому время релаксации импульса можно оценить снизу величиной транспортного времени релаксации

где $\sigma _{{el}}^{{}}$ – дифференциальное сечение упругого рассеяния электронов на угол $\vartheta $.

Плазменная и ларморовская частоты, электропроводность и транспортное время релаксации импульса являются классическими фундаментальными параметрами ионизованной среды [18, 19]. Представленный анализ лишь показывает их связь с задачей (1.1), (1.2). Из анализа непосредственно следуют условия малости импульса, приобретаемого электроном в области $V_{m}^{\delta }$:

Первые три из неравенств (4.2) влекут очевидные условия малости возмущения, вносимого сглаживанием коэффициентов уравнения (1.1):

где ${{r}_{D}} = {{v\left( {{{p}_{c}}} \right)} \mathord{\left/ {\vphantom {{v\left( {{{p}_{c}}} \right)} {{{\omega }_{p}}}}} \right. \kern-0em} {{{\omega }_{p}}}}$ – радиус Дебая, ${{r}_{L}} = {{v\left( {{{p}_{c}}} \right)} \mathord{\left/ {\vphantom {{v\left( {{{p}_{c}}} \right)} {{{\omega }_{L}}}}} \right. \kern-0em} {{{\omega }_{L}}}}$ – радиус Лармора, λ_tr = ${{\lambda }_{{tr}}} = {{\left( {v\left( {{{p}_{c}}} \right)\int_{ - 1}^1 {\sigma \left( {{\mathbf{p}}{\kern 1pt} ',{\mathbf{p}}} \right)\left( {1 - \cos \vartheta } \right)d{\mathbf{p}}{\kern 1pt} '} } \right)}^{{ - 1}}}$ – транспортная длина пробега.

Последнее из соотношений (4.2) имеет нетривиальное следствие. В соответствии с уравнениями Максвелла (1.2) время затухания электрического поля за счет тока проводимости определяется величиной ${{4\pi } \mathord{\left/ {\vphantom {{4\pi } {{{\tau }_{\sigma }}}}} \right. \kern-0em} {{{\tau }_{\sigma }}}} = {{4\pi {{e}^{2}}n{{\tau }_{g}}} \mathord{\left/ {\vphantom {{4\pi {{e}^{2}}n{{\tau }_{g}}} {m{{\gamma }^{3}}}}} \right. \kern-0em} {m{{\gamma }^{3}}}} = \omega _{p}^{2}{{\tau }_{g}}$. Тогда рассматриваемое соотношение принимает вид $\omega _{p}^{{}}{{\tau }_{g}}\omega _{p}^{{}}{{\tau }_{\delta }} \ll 1$.

С учетом первого из соотношений (4.2) данное условие не допускает, чтобы величина $\omega _{p}^{{}}{{\tau }_{\delta }}$ существенно превысила единицу. Формально это означает, что столкновениями можно пренебречь по сравнению с плазменными колебаниями, т.е. справедливо приближение бесстолкновительной плазмы. Но сравнение тока быстрых электронов с дрейфовым током термолизованных электронов предполагает приближение проводимости. Это сугубо столкновительный эффект. В этом приближении термализованные электроны под действием электрического поля приходят в движение мгновенно. Следует отметить, что в ряде задач вторичные электроны низкой энергии, образующиеся при ударной ионизации среды, удобно рассматривать как электроны проводимости [20]. Для данного приближения рассмотрим характерный пространственный параметр – глубину скин-слоя [21] ${{d}_{\sigma }} = с\sqrt {{2 \mathord{\left/ {\vphantom {2 {\omega _{p}^{3}{{\tau }_{g}}}}} \right. \kern-0em} {\omega _{p}^{3}{{\tau }_{g}}}}} $, где предполагается, что электрическое поле имеет характерную частоту, совпадающую с плазменной. Если $\omega _{p}^{{}}{{\tau }_{g}}$ стремится к бесконечности, то глубина скин-слоя стремится к нулю. При этом очевидно, что электромагнитное поле в приближении проводимости не должно искажаться существенно областью $V_{m}^{\delta }$. Это означает, что должно выполняться условие ${{d}_{\sigma }} \gg \delta $, что действительно невозможно, если $\omega _{p}^{{}}{{\tau }_{\delta }} \gg 1$.

5. ОСОБЕННОСТИ ЧИСЛЕННОГО АЛГОРИТМА

Численное решение уравнений (1.1)–(1.4) осуществляется с помощью алгоритмов, подробно представленных в работах [13], [15], [20], [22–26]. Например, в работе [22] изложен алгоритм численного решения уравнений движения частиц и интерполяции создаваемого ими тока в узлы разностной сетки. На данной сетке аппроксимируются уравнения Максвелла, разностная схема решения которых представлена в работе [23]. Рассмотрим некоторые элементы алгоритмов, позволившие избежать проблем с разрывами свойств рассеивающей среды.

Первым элементом алгоритма является разностная сетка [23]. Пусть Ω – прямоугольная подобласть пространства $\mathbb{R}_{r}^{3}$, в которой рассматривается численное решение начально-краевой задачи для уравнений (1.1)–(1.4). Введем в Ω прямоугольную декартову разностную сетку для разностного решения уравнений Максвелла. По переменной х сетка имеет вид

По переменным $(y,z)$ сетка вводится аналогично. Узлы сетки разместим так, чтобы разрывы сечений в интеграле столкновений кинетического уравнения (1.1) и электрофизических параметров в уравнениях Максвелла (1.4) оказались на поверхностях $x = {{x}_{i}}$, $y = {{y}_{j}}$ и $z = {{z}_{k}}$. Таким образом, расчетная область делится на декартовы прямоугольные ячейки. Значения коэффициентов уравнений Максвелла задаются в точках сетки с дробными индексами. Эти точки совпадают с центрами ячеек. Компонента напряженности электрического поля, нормальная к поверхности разрыва диэлектрической проницаемости, может претерпевать разрыв при переходе через эту поверхность. Касательная к этой поверхности компонента напряженности электрического поля при этом непрерывна. Поэтому сеточные компоненты электрического поля определяются в середине соответствующих ребер указанных выше прямоугольных ячеек.

Будем рассматривать ячейки разностной сетки как множества ${{D}_{m}}$, $m = 1,...,M$, определенные в предыдущем разделе. Введем величину

В однородном материале вероятность испытать столкновение за время $\tau $ определяется по формуле $p = \exp \left\{ { - \tau {{\sigma }^{t}}v} \right\}$. Вероятность $\xi = 1 - p$ не испытать столкновение за время $\tau $ является случайной величиной, равномерно распределенными на отрезке $\left[ {0,1} \right]$. Тогда время жизни $\tau $ до первого столкновения определяется по формуле $\tau = - {{\log \left( \xi \right)} \mathord{\left/ {\vphantom {{\log \left( \xi \right)} {{{\sigma }^{t}}v}}} \right. \kern-0em} {{{\sigma }^{t}}v}}$.

Рассмотрим координату ${\mathbf{r}} = {{{\mathbf{r}}}^{s}}\left( {{{t}_{n}},\tilde {t},{\mathbf{\tilde {r}}},{\mathbf{\tilde {p}}}} \right)$ и импульс ${\mathbf{p}} = {{{\mathbf{p}}}^{s}}\left( {{{t}_{n}},\tilde {t},{\mathbf{\tilde {r}}},{\mathbf{\tilde {p}}}} \right)$ частицы в момент времени $t = {{t}_{n}}$. Будем считать, что частица в этот момент времени находится в ячейке с номером $\left( {i,j,k} \right)$, соответствующей множеству ${{D}_{m}}$. Рассмотрим луч с началом в точке ${\mathbf{r}} = {{{\mathbf{r}}}^{s}}\left( {{{t}_{n}},\tilde {t},{\mathbf{\tilde {r}}},{\mathbf{\tilde {p}}}} \right)$ в направлении вектора ${\mathbf{p}} = {{{\mathbf{p}}}^{s}}\left( {{{t}_{n}},\tilde {t},{\mathbf{\tilde {r}}},{\mathbf{\tilde {p}}}} \right)$. Рассмотрим отрезки, соединяющие точку ${\mathbf{r}} = {{{\mathbf{r}}}^{s}}\left( {{{t}_{n}},\tilde {t},{\mathbf{\tilde {r}}},{\mathbf{\tilde {p}}}} \right)$ с ближайшей (№ 1) и наиболее (№ 2) удаленной из точек области $V_{m}^{\delta }$ на луче. Обозначим через ${{\Delta }_{\delta }}$ – время, потребное частице, чтобы преодолеть отрезок между точками ${\mathbf{r}} = {{{\mathbf{r}}}^{s}}\left( {{{t}_{n}},\tilde {t},{\mathbf{\tilde {r}}},{\mathbf{\tilde {p}}}} \right)$ и № 1. Рассмотрим $T = \min \left\{ {\Delta ,{{\Delta }_{\delta }},\tau } \right\}$. Вычислим ${\mathbf{r}} = {{{\mathbf{r}}}^{s}}\left( {{{t}_{n}} + T,\tilde {t},{\mathbf{\tilde {r}}},{\mathbf{\tilde {p}}}} \right)$ и ${\mathbf{p}} = {{{\mathbf{p}}}^{s}}\left( {{{t}_{n}} + T,\tilde {t},{\mathbf{\tilde {r}}},{\mathbf{\tilde {p}}}} \right)$. Для этого используется центрированная по времени разностная схема, рассмотренная в [22]. В конце вычисления пересчитывается время жизни частицы на данном временном шаге: $\Delta = \Delta - T$.

Если оказалось, что $T = \tau $, то рассматривается столкновение. Если частица поглощается, то из рассмотрения она выбывает. Если частица рассеивается, то ее новый импульс определяется статистически. Если в столкновении рождается новая частица, то она включается в рассмотрение с параметрами, также определенными статистически [13], [22].

Если оказалось, что $T = {{\Delta }_{\delta }}$, то частица поглощается в точке № 1 и генерируется в точке № 2 с тем же импульсом.

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

ЗАКЛЮЧЕНИЕ

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

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

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