\documentclass[11pt,reqno]{amsart}
\usepackage[cp1251]{inputenc}
\usepackage[english,russian]{babel}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage[colorlinks]{hyperref}
\usepackage{caption}
\usepackage{subcaption}
\usepackage{titlesec}
\titlespacing{\section}{0pt}{2.5ex}{1.5ex}
\titlespacing{\subsection}{0pt}{1.5ex}{1ex}
\titlespacing{\subsubsection}{0pt}{1.5ex}{1ex}
\titleformat{\section}{\large\bfseries\centering}{\thesection}{1em}{}
\titleformat{\subsection}[runin]{\bfseries}{\thesubsection.}{0.5em}{}[.\mbox{\ }]
\titleformat{\subsubsection}[runin]{\bfseries}{\thesubsubsection.}{0.4em}{}[.\mbox{\ }]


%%%%%---------- ORCiD
\newcommand\orcidicon[1]{\href{https://orcid.org/#1}{\includegraphics[scale=0.02]{orcid.pdf}}}


\def\udcs{???.??} %Here the author places classificators of the paper according to Russian classification system
\def\mscs{??X??} %Here the author places  classificators according to the AMS classification list.
\setcounter{page}{144}

\newtheorem{lemma}{Lemma}
\newtheorem{theorem}{Theorem}
\newtheorem{definition}{Definition}
\newtheorem{corollary}{Corollary}
\newtheorem{proposition}{Proposition}



\hypersetup{linkcolor=blue,filecolor=black,urlcolor=blue, citecolor=blue}

\def\logo{\textcolor{blue}{{\bf\Huge S\raisebox{0.2ex}{\hspace{0.55ex}\raisebox{0.05ex}e\hspace{-1.65ex}$\bigcirc$}MR}}}

\def\semrtop
     {
  \vbox{
     \noindent\logo \hspace{-3mm}  \parbox{12cm}{     \begin{center}
     {\LARGE \textsf{СИБИРСКИЕ ЭЛЕКТРОННЫЕ}} \\[2mm]
     {\LARGE \textsf{МАТЕМАТИЧЕСКИЕ ИЗВЕСТИЯ}} \\[2mm]
     {\Large \textsf{Siberian Electronic Mathematical Reports}} \\[1mm]
     {\Large\tt{\textcolor{blue}{http://semr.math.nsc.ru}}}\\[1mm]
     \raisebox{1ex}{ISSN }
     \vspace*{-3mm}
%     {\small Институт математики им. С.Л.Соболева СО РАН}\\[-1mm]
%     {\small Московский автомобильно-дорожный государственный технический университет (МАДИ)}
     \end{center}
}
     \vspace{-1mm}
     \noindent
      \begin{tabular}{c}
     \hphantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa!} \\
     \hline\hline
     \end{tabular}

     \vspace{3mm}
     {\flushleft\it Том 20, № 2, стр. 144--145 (2023) \hspace{\fill}{\rm\small УДК 519.63 \hspace{-4mm}}}
     \newline  %Volume, pages in Russian, filled by Editorial board
     	{\rm\small \textcolor{blue}{https://doi.org/10.33048/semi.2024.20.???}}\hphantom{aaaaaaaaaaaaaaaaaaaaaaaa\;}
     	{\rm\small MSC 93A30\ \ }
  }%\vbox
}


\begin{document}
\renewcommand{\refname}{References}
% \renewcommand{\proofname}{Proof.}
% \renewcommand{\figurename}{Рис.}

\thispagestyle{empty}

\title[Анализ модели переноса плазмы в винтовом магнитном поле]{\large Анализ модели переноса плазмы в винтовом магнитном поле}

\author[Г.Г. Лазарева, И.П. Оксогоева и А.В. Судников]{
{\bf \large \textsc{Г.Г. Лазарева\protect\orcidicon{0000-0001-8810-918X}, И.П. Оксогоева\protect\orcidicon{0000-0001-6000-188X}} и А.В. Судников}
\\ \\
\textit{Представлено} \textsc{}
%\textit{(заполняется редактором)}
}%

\address{Галина Геннадьевна Лазарева  
\newline\hphantom{iii} Московский государственный университет им. М.В. Ломоносова (МГУ),
\newline\hphantom{iii} 119991, Россия, Москва,
\newline\hphantom{iii} ул. Колмогорова, 1}%
\email{\textcolor{blue}{lazareva-gg@pfur.ru}}%

\address{Ирина Павловна Оксогоева
\newline\hphantom{iii} Российский университет дружбы народов имени Патриса Лумумбы (РУДН),
\newline\hphantom{iii} 117198, Россия, Москва,
\newline\hphantom{iii} ул. Миклухо-Маклая, 6}%
\email{\textcolor{blue}{oksogi@mail.ru}}%

\address{Судников Антон Вячеславович
\newline\hphantom{iii} Федеральное государственное бюджетное учреждение науки Институт ядерной физики им. Г.И. Будкера Сибирского отделения Российской академии наук (ИЯФ им. Г.И. Будкера СО РАН),
\newline\hphantom{iii} 630090, Россия, Новосибирск,
\newline\hphantom{iii} проспект Академика Лаврентьева, 11}%
\email{\textcolor{blue}{a.v.sudnikov@inp.nsk.su}}%


\thanks{\sc Lazareva G.G., Oksokoeva I.P., Sudnikov A.V.
Analysis of the plasma transfer model in a helical magnetic field}
\thanks{\copyright \ 2024 Лазарева Г.Г., Оксогоева И.П., Судников А.В}
\thanks{\rm Работа выполнена при поддержке Министерства науки и высшего образования Российской Федерации (мегагрант соглашение №075-15-2022-1115)}
\thanks{\it Поступила 1 ноября 2024 г., опубликована 31 декабря 2024 г.}%

\semrtop \vspace{1cm}
\maketitle {\small
\vspace{-12pt}
%\centerline{{\it Communicated by} {\sc P.P. Petrov}
%({\it filled by Editorial board})
%}
\begin{quote}
\noindent{\bf Abstract:} The paper presents the results of an analysis of a mathematical model of plasma transfer in a spiral open magnetic SMOLA trap created at the Budker Institute of Nuclear Physics of the Siberian Branch of the Russian Academy of Sciences. Plasma confinement in the trap is carried out by transmitting a pulse from a magnetic field with helical symmetry to a rotating plasma. The mathematical model is based on the stationary plasma transfer equation. The model introduced the consideration of the functional type of the coefficients of the model, set in the range corresponding to the physics of the confinement process. The dependences of the model parameters on the coordinates are given, at which a qualitative correspondence of the calculation to experimental data is observed. A mathematical model has been developed to predict the parameters of plasma confinement in designed trapes with a spiral magnetic field.

\medskip
\noindent{\bf Ключевые слова:} mathematical modeling, transfer equation, helical magnetic field
 \end{quote}
}

\bigskip

\section{Введение}
В ходе исследований термоядерного синтеза открытые магнитные системы для удержания плазмы рассматриваются в качестве возможных конфигураций для термоядерного реактора. В настоящее время достигнут большой прогресс в понимании физики открытых магнитных систем. В изучаемой ловушке формируются самоподдерживающиеся конфигурации плазмы и магнитного поля, способные существовать какое-то время, достаточно длительного по масштабам плазменных процессов. Установка Спиральная Магнитная Открытая Ловушка (СМОЛА) разработана и построена в 2017 году в Институте ядерной физики СО РАН им. Г.И. Будкера \cite{1}. Новая ловушка СМОЛА основана на использовании винтовой конфигурации магнитного поля для удержания плазмы в устройстве. Удержание плазмы осуществляется за счет тормозящей силы, возникающей при вращении плазмы в скрещенных электрическом и винтовом магнитном полях. Плазма удерживается в центральной части, которая работает как классическая газодинамическая ловушка. Роль предложенных секций со спиральным магнитным полем заключается в уменьшении потерь частиц вдоль магнитного поля и, следовательно, в уменьшении потерь энергии. Математическая модель переноса вещества в винтовом магнитном поле построена на основе системы уравнений магнитной гидродинамики (МГД) в цилиндрической системе координат с учетом симметрии \cite{2}. Рассмотрен стационарный случай, что обусловлено ходом экспериментов, когда после короткого этапа установления происходит интересующий экспериментаторов этап. Система уравнений МГД для данной задачи, включающая в себя поперечные переносы за счёт диффузии, изложена в безразмерном виде в \cite{3}.

Винтовые течения и течение жидкости в винтовом магнитном поле находятся на фронтире современной механики. В работе \cite{3} сделано обобщение решения В.А. Стеклова задачи о нестационарных винтовых течениях ньютоновской жидкости на случай жидкости второго порядка. В 2024 году на конференции Математические проблемы механики сплошных сред член-корреспондент РАН В.В. Пухначев показал известные примеры винтовых течений, представил плоские и вращательно симметричные аналоги винтовых течений, а также обобщение этих результатов на случай движения жидкости второго порядка \cite{4}. В работах \cite{5,6} приведено исследование устойчивости к малым возмущениям течения вязкой электропроводящей жидкости в трубе кольцевого сечения при наличии спирального магнитного поля в приближении малых магнитных чисел Рейнольдса. Использование симметрии приводит к зависимости параметров течения от радиуса и комбинации угла наклона и шага винта. В представленной модели параметр $\zeta$ является отношением скорости звука и продольной скорости движения магнитных возмущений при вращении плазмы в собственном амбиполярном электрическом поле, тем самым содержит информацию о геометрии магнитного поля. Кинетические модели процессов в установках, основанных на принципе многопробочного удержания, развиваются как в сотрудничестве ведущих физиков-теоретиков плазменных лабораторий ИЯФ СО РАН с вычислителями из ИВМиМГ СО РАН, так и самостоятельно. В работах \cite{7, 8, 9} на основе кинетической модели проведено исследование взаимодействия пучка ионов, инжектированного в осесимметричную магнитную ловушку под углом к центральной оси, с фоновой плазмой ловушки в двумерной постановке. Проведено моделирование истечения плазмы с максвелловским распределением по скоростям ионов из диамагнитной ловушки.  Наблюдалось раскручивание плазмы вокруг оси ловушки в направлении, совпадающем с направление циклотронного вращения в вакуумном магнитном поле. Вычислители из ИВМиМГ СО РАН работают с численными моделями, как кинетическими, так и гибридными (кинетическое описание ионной компоненты и гидродинамическое – для электронов), позволяющих учитывать большое разнообразие физических процессов с различными пространственно-временными масштабами, характерными для лабораторных экспериментов в линейных плазменных системах.  В работах \cite{10, 11} представлена кинетическая модель истечения плазмы из открытой ловушки, описывающая переход к газодинамическому режиму с учетом влияния амбиполярного потенциала. В результате расчетов обнаружен переходный режим истечения плазмы, в котором продольные потери происходят с образованием струи холодных ионов, и получено самосогласованное аналитическое решение для струи в ловушке с «плоским дном». Математическое моделирование в магнитогидродинамическом приближении проводилось в совместных работах сотрудников ИЯФ СО РАН и ИВТ СО РАН \cite{12}, где выполнена оценка влияния внешнего магнитного поля на динамику разлета плазменного облака. Традиционно лидером в разработке и развитии МГД моделей, программных комплексов и проведении вычислительных экспериментов для анализа процессов в существующих и перспективных магнитных ловушках для удержания плазмы, является Институт прикладной математики им. М.В. Келдыша РАН. Математические модели равновесных плазменных конфигураций в ловушках можно рассматривать в их распрямленных в цилиндр аналогах \cite{13}. Математический аппарат таких моделей — это двумерные краевые задачи с дифференциальным уравнением Грэда-Шафранова для функции магнитного потока. В ИПМ РАН в МГД приближении решаются задачи, близкие к предлагаемой в проекте. В работе \cite{14} представлена модель для изучения инжекции потока плотной плазмы в многопробочную ловушку, совмещенную с плазменным ускорителем.

Целью математического моделирования является демонстрация подавления потерь частиц и энергии из ловушки при включении секций с многопробочным магнитным полем. Точного аналитического решения для теории винтового удержания не построено. В настоящее время показано соответствие наблюдаемых результатов приближённым теоретическим оценкам. Винтовая конфигурация магнитного поля не дает возможности снизить размерность задачи без значительных упрощений, сделанных в работе \cite{2}. Результаты осуществляемого численного решения уравнения переноса плазмы востребованы как сложно вычисляемая оценка, учитывающая новую сложную геометрию магнитного поля. Рассматриваемая модель не требует больших вычислительных ресурсов, в отличие от моделей, построенных на решении системы уравнений МГД в полной трехмерной постановке \cite{15}. В перспективе возможен переход к пространственным моделям для исследования новых магнитных конфигураций многопробочных секций на основе уравнений МГД для сжимаемых сред. Например, комбинации однородного, спирального и гофрированного магнитных полей. Для построения полной математической пространственной модели сжимаемой среды в самосогласованных электромагнитных полях необходимо уточнить диапазон допустимых значений, в первую очередь коэффициента диффузии. Двумерные постановки для моделирования движения ионизированного газа в самосогласованных электромагнитных полях \cite{16, 17}, построенные с учетом геометрии задачи, имеют ряд преимуществ и широко используются. 

\section*{Постановка задачи}
Исследуется движение плазмы в поперечном сечении $[0, r_{max}]\times[0,z_{max}]$   центральной части установки (Рис. \ref{ris:im1a}), имеющей форму цилиндра радиуса $8$ см и длины $216$ см. Вещество входит в область удержания из источника плазмы через левый торец цилиндра и выходит в расширитель через правую границу. Расчетная область это поперечное сечение цилиндра в плоскости $(r, z)$.
\begin{figure}[h]
\centering
\begin{subfigure}[b]{0.49\textwidth}
\includegraphics[width=1\textwidth]{Fig1a.png}
\caption{}
\label{ris:im1a}
\end{subfigure}%
\begin{subfigure}[b]{0.49\textwidth}
\includegraphics[width=1\textwidth]{Fig1b.PNG}
\caption{}
\label{ris:im1b}
\end{subfigure}
\caption{Расчетная область в поперечном сечении центральной части ловушки (A). Экспериметальные значения концентрации плазмы торцах установки и их интерполяция полиномами (B): на входе (квадратики), на выходе (кружочки).}
\label{ris:im1}
\end{figure}
В статье \cite{2} получены выражения для компонентов радиального и продольного переноса частиц в винтовом магнитном поле. Система уравнений описывает динамику плазмы в МГД-приближении в аксиально-симметричной постановке. Различия в движении запертых и пролётных ионов учитываются в виде эффективной силы трения, зависящей от взаимной скорости компонент и доли запертых частиц. Продольная сила, действующая на плазму, возникает в результате взаимодействия радиального электрического тока захваченных ионов с азимутальной компонентой винтового магнитного поля. Учитывается диффузия плазмы поперёк магнитного поля. Исключение зависимых переменных сводит систему уравнений к уравнению неразрывности потока для концентрации плазмы: 
\begin{equation}
\begin{split}
&\frac{\partial}{\partial z}\frac{1}{\Lambda\kappa}\frac{\partial (Tu) }{\partial z}  + \frac{\partial}{\partial z}((1+\kappa)\zeta E u)
+\frac{\partial }{\partial r} Z^{-1}\zeta \frac{\partial (Tu) }{\partial z}  + \frac{\partial }{\partial r} D \frac{\partial u }{\partial r} = 0,
\\
&\frac{\partial u}{\partial r}(0,z)=0, \, u(r,0) = u_L(r), \, u(1,z) = 0, \, u(r,1)=u_R(r).
\end{split} \label{1}
\end{equation}
где $u$ -- концентрация вещества, $T=T_i + T_e$, $T_i$ и $T_e$ -- ионная и электронная температуры, $\Lambda$ - отношение длины системы к длине свободного пробега иона, $\kappa$ -- доля запертых частиц,  $l$ -- длина системы вдоль силовой линии, $Z_1$ -- среднее зарядовое число одного иона,  $D$ -- коэффициент диффузии в поперечном поле, $Е$ -- напряженность электрического поля. Параметр  $\zeta =c_s/V_z$ является отношением скорости звука $c_s=(T_e/M)^{1/2}$ и продольной скорости $V_z$ движения магнитных возмущений при вращении плазмы в собственном амбиполярном электрическом поле. В уравнении \eqref{1} физические величины обезразмерены на: $r_0=a$, \, $z_0=l$, \, $\phi_0=T_e/e$, \, $u_0=u_{max}$, \, $T_0=T_e$, где $a$ -- граница камеры.

В моделируемом эксперименте происходит установление параметров плазмы и электромагнитного поля в течении $40$ мс и затем наступает фаза, когда процесс стационарен в течении $120$ мс, после чего разряд отключается. Основная задача натурных и вычислительных экспериментов состоит в исследовании удержания плазмы, при которых все параметры постоянны. В безразмерных величинах область является единичным квадратом. В соответствии с экспериментом считаем, что вещество не достигает стенок ловушки, на оси $z$ ставится условие симметрии. На входе и выходе плазмы (Рис. \ref{ris:im1b}) задано граничное условие Дирихле:  $u(r,0)=u^{L}(r)$ и $u(r, z_{max})=u^{R}(r)$. 

Перенос излучения \cite{18} в модели \eqref{1} не учитывается. 
При характерной температуре электронов $T_e = 30$ эВ в водородной плазме с зарядовым числом (даже с учётом примесей) $Z \sim 1$  и плотностью $n$ = $10^{12}$ см$^{-3}$ как тормозным, так и рекомбинационным излучением плазмы можно пренебречь. Эмитированный с нагретого катода поток состоит из электронов, имеющих энергию, близкую к максимуму сечения ионизации. Для водорода линии излучения есть только у нейтральных атомов. В рассматриваемом случае степень ионизации велика $(>90\%)$. Выражение для тормозного излучения \cite{19} даёт значение около $(10^{-7}-10^{-6})$ Вт/см$^3$: $P_1=1.54 \cdot 10^{-32} \cdot Z^2 n_e$ $n_i \sqrt{T_e}$.

Мощность рекомбинационного излучения может быть оценена \cite{20} как $(10^{-7}
-10^{-5})$ Вт/см$^3$: $P_2=4.1 \cdot 10^{-31} {Z^4} n_e n_i \sqrt{T_e}$.

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


В работе \cite{23} было показано, что введение в коэффициенты \eqref{1} зависимости от осевой координаты позволяет получить распределение концентрации с соответствующим наблюдаемому в эксперименте характеру убывания по длине установки. Следующим шагом является детальная проработка функциональных зависимостей коэффициентов \eqref{1} от координат. Выделим постоянные параметры установки: $l =216$ см, $a =8$ см, $T_i=4$ Эв, $Z^{-1}\approx 1$, $\Lambda \approx 216$ см, $B=700$ Гс. Глубина гофрировки $R(r, R_m)$, определяющая долю запертых частиц $\kappa$, зависит только от радиуса установки. Для заданной средней глубины гофрировки $R_m=1.52$ выполняется $R(r, R_m)=2(R_m -1)(r/a)^{2}+1$. 

Представленные результаты получены в ходе определения характера убывания вдоль оси ловушки электронной температуры $T_e$, напряженности электрического поля $E$, доли запертых частиц $\kappa$ и параметра $\zeta$. Наибольший интерес вызывает коэффициент диффузии $D\in[0,1]$, который на данном этапе исследований считается постоянным. Для использования математической модели для более сложно устроенных перспективных установок, проектируемых в настоящее время, необходимо получить более точную оценку для величины коэффициента диффузии. Получить эту величину в ходе экспериментов невозможно. Знание вклада диффузии позволит построить математическую модель удержания плазмы в приближении сжимаемой среды с учетом самосогласованных электрического и магнитного полей. Такая модель должна учитывать зависимость всех параметров процесса от пространственных координат и не будет содержать упрощений, основанных на геометрии спирального магнитного поля.


\section{Метод решения}

Можно заметить, что задача \eqref{1} содержит параметр $\zeta(r) = 1/Ar$, полученный аппроксимацией экспериментальных данных, здесь $A=20$. На оси при $r=0$ параметр $\zeta$ ограничен значением при $r=\rho_B$, где  $\rho_B = V_Tmc/eB$ — ларморовский радиус. Для параметров установки СМОЛА $\rho_B \approx 0.3-0.4$ см. Это обусловлено тем, что ион движется по ларморовской орбите (вращается в равномерном магнитном поле, имеющем большую величину, чем спиральное), поэтому его радиальная координата осциллирует. Сильное однородное магнитное поле позволяет считать, что ионы движующиеся по орбитам перепендикулярны оси. В рассматриваемой модели все воздействия усредняются, иону приписывается координата центра окружности, по которой он движется в магнитном поле. То есть, координата равна нулю для ионов, которые облетают ось вокруг и находятся на расстоянии ларморовского радиуса от неё. Поэтому для исключения особенностей решения при расчетах в окрестности оси симметрии безразмерный параметр $\zeta$ задается следующим образом: 
\begin{equation}
\zeta(r) =
  \begin{cases}
    1/Ar, & \quad r>\rho_B, \\
    1/A\rho_B, & \quad  \rho_B \geq r \geq 0.
  \end{cases}
\label{2}
\end{equation}
Задача \eqref{1} была реализована методом установления и более экономичным методом Зейделя \cite{24} и аппроксимации смешанной производной на шаблоне \cite{25}. Наряду с другими преимуществами, метод Зейделя интересен удобством использования в цилиндрических координатах. Принцип выражения искомого элемента через соседние по схеме точек вида «крест» универсально и не зависит от выбора системы координат. Выбор граничного условия Дирихле, описывающего выход плазмы в расширитель через правую границу, позволяет контролировать верное изменение интегральной характеристики плотности вещества при изменении глубины гофрировки магнитного поля, диффузии и потенциала плазмы. Такой контроль постоянно осуществлялся при изменении коэффициентов уравнения.

\section{Результаты моделирования}

В расчётах \cite{26, 27, 28} наблюдается эффект пинчевания (уменьшения среднего радиуса) плазменной струи, проявляющийся и в эксперименте (Рис. \ref{ris:im2a}). 

  \begin{figure}[ht!]
  \centering
  \begin{subfigure}[b]{0.49\textwidth}
  \includegraphics[width=1\textwidth]{fig.PNG}
  \caption{}
  \label{ris:im2a}
  \end{subfigure}%
  \begin{subfigure}[b]{0.49\textwidth}
  \includegraphics[width=1\textwidth]{E.PNG}
  \caption{}
  \label{ris:im2b}
  \end{subfigure}
  \caption{ Расчетное распределение концентрации плазмы (A). Зависимость от радиуса установки напряженности $Е$ электрического поля (B): экспериментальные данные (квадратики) и интерполирующий полином (линия).}
\end{figure}
Функции параметров модели были модифицированы с использованием распределения по оси известных из эксперимента точек на оси симметрии \cite{23}. Характер убывания плотности плазмы вдоль длины центральной части ловушки при всех параметрах экспериментов имеет экспоненциальный вид. Фактически именно эта особенность распределения плотности плазмы выражает эффект удержания и очень важна в случае использования граничного условия Неймана на правой границе. В случае верного расчета убывания плотности по длине установки можно сделать априорную оценку размера расчетной области для задания граничного условия равенства нулю концентрации плазмы на бесконечности. Рассматриваемые экспериментальные данные на оси симметрии интерполированы функцией вида:
\begin{equation*}
u(0,z) = Ae^{-z/z_0}, \quad A=0.8, \quad z_0=0.38.
\end{equation*}
Необходимо отметить, что для всех экспериментов характер убывания остается практически неизменным при определении параметров разными способами: через точки на оси $z$, через точки максимальных значений плотности сечений, через значения плотности по сечениям. Следовательно, вид функции  $u(0,z)$ можно рассматривать как дополнительную информацию о процессе.

Измеряемый с помощью зондов потенциал плазмы в экспериментах на установке СМОЛА зависит от параметров эксперимента. Известно максимальное значение безразмерного потенциала для установок следующего поколения, в которых возможно использование принципа винтового удержания. Максимум $\phi$ варьируется в пределах от $2T_e/e$ до $3T_e/e$. Максимальное значение безразмерного потенциала $\phi$ для установок следующего поколения, в которых возможно использование принципа винтового удержания, также лежит в этих же пределах \cite{21, 22}. Экспериментально наблюдаемое распределение потенциала в центральной области плазмы (при значениях радиуса плазмы меньших $0.6$ в безразмерных величинах) близко к квадратичному, далее в периферийной области плазмы производная потенциала по радиусу снижается. Погрешность измерения потенциала в эксперименте составляет около $5\%$. 

Экспериментальные данные содержат таблицу значений напряженности электрического поля $E(r,0)$ на левой границе расчетной области. В модели учтено снижение величины напряженности электрического поля $Е$ с увеличением $z$ за счёт наличия поперечной проводимости плазмы как $E(r,z)=(1-z/2)E(r,0)$.  Распределение по радиусу $E(r,0)$ задается, как полином, который интерполирует экспериментальные данные (Рис. \ref{ris:im2b}):
\begin{equation*}
E(r,0)=0.155-0.71r+7.46r^2-3.82r^3-10.53r^4+7.55r^5.
\end{equation*}
Степень и коэффициенты полинома подобраны таким образом, чтобы его отклонение от экспериментально измеренных в референсном эксперименте значений было сопоставимым с экспериментальной погрешностью.

Рассмотрим вклад влияния производных в решение уравнения \eqref{1}. Вторая производная искомой функции по оси определяет требуемую вогнутость сечения вдоль длины установки (Рис. \ref{ris:im3a}). Первое слагаемое, содержащее вторую производную концентрации по $z$, включает коэффициенты: $\kappa$, $T$ и постоянную $\Lambda$. Если эти коэффициенты относительно малы, то влияние второй производной концентрации по $z$ на вогнутость решения задачи вдоль оси может оказаться не существенным. Первая производная искомой функции по оси, содержащаяся во втором слагаемом, принимает отрицательные значения, что соответствует убыванию решения по $z$. Третье слагаемое содержит смешанную производную, которая принимает сравнительно малые значения и отвечает за наблюдаемое сужение по радиусу поперечных сечений. Вторая производная по радиусу от искомой функции меняет знак в соответствии с поведением радиального распределения концентрации. Только первое слагаемое содержит производную, оказывающую влияние на вогнутость решения вдоль длины установки. Следовательно, для получения верным образом убывающей по $z$ плотности требуется, чтобы коэффициенты первого слагаемого были достаточно велики.

В модель введено убывание вдоль оси ловушки еще для двух параметров процесса. В допустимом диапазоне убывания можно записать функцию зависимости электронной температуры от координат в виде $T_e(r,z)=(1-5r^2/6)(1-z/4)$. Параметр $\zeta$ определяется через электронную температуру, убывающую по длине: $\zeta(r,z)=c_c/V_z \approx (48.7466 \sqrt{T_e})/Br$

Наибольшее влияние на характер решения оказывает доля запертых частиц $\kappa = 1- 1/R(R_m,r)$. Оценка показывает, что доля запертых частиц уменьшается в два раза по длине установки: $\kappa(r,z) = k(r,0)(1-z/2)$. 

\begin{figure}[ht!]
    \centering
    \begin{subfigure}[b]{0.49\textwidth}
      \includegraphics[width=1\textwidth]{Fig3a.PNG}
      \caption{}
      \label{ris:im3a}
      \end{subfigure}%
      \begin{subfigure}[b]{0.49\textwidth}
      \includegraphics[width=1\textwidth]{Fig3b.PNG}
      \caption{}
      \label{ris:im3b}
      \end{subfigure}\\
      \begin{subfigure}[b]{0.49\textwidth}
        \includegraphics[width=1\textwidth]{Fig3c.PNG}
        \caption{}
        \label{ris:im3c}
        \end{subfigure}%
        \begin{subfigure}[b]{0.49\textwidth}
        \includegraphics[width=1\textwidth]{Fig3d.PNG}
        \caption{}
        \label{ris:im3d}
        \end{subfigure}
    \caption{ Графики концентрации плазмы в сечении вдоль оси симметрии $u(0,z)$ (A), в радиальном сечении при $u(r,0.16)$ (B), $u(r,0.3)$ (C), $u(r,0.4)$ (D). Экспериментальные данные (квадратики), результаты расчетов при введении зависимости от $z$ функций  $E$ (1), $T$ (2) и $\kappa$ (3).}
    \label{ris:im3}
\end{figure}

\begin{figure}[ht!]
    \centering
    \begin{subfigure}[b]{0.49\textwidth}
      \includegraphics[width=1\textwidth]{figN.png}
      \caption{}
      \label{ris:im4a}
      \end{subfigure}%
      \begin{subfigure}[b]{0.49\textwidth}
      \includegraphics[width=1\textwidth]{FigN_b.PNG}
      \caption{}
      \label{ris:im4b}
      \end{subfigure}
    \caption{Распределение концентрации плазмы для постановки (3) при $z=0$ (1), $z=0.4$ (2), $z=1$ (3). }
    \label{ris:im4}
\end{figure}

Сравним результаты расчетов с учетом убывания параметров вдоль оси с результатами без такого учета и данными экспериментов в четырех сечениях (Рис. \ref{ris:im3}). На графиках показаны результаты трех новых расчетов, каждый из которых проведен при изменении только одного параметра. При изменении трех параметров происходит суммирование эффектов. График концентрации в первом сечении вдоль оси установки (Рис. \ref{ris:im3a}) показал, что только параметр $\kappa$ влияет на вогнутость решения вдоль оси. Сечения при $z=0.16$, $z=0.3$, $z=0.4$ имеют общий характер поведения решения (Рис. \ref{ris:im3b}, \ref{ris:im3c} и \ref{ris:im3d}). Выбор сечений обусловлен наличием данных новых экспериментов. Графики зависимости решения от радиуса показывают, что изменение $\kappa$ и $T$ дают меньшее и большее значение концентрации на оси при сохранении поведения решения.  Изменение напряженности электрического поля $E$ не влияет на решение. 

При этом уравнение переноса \eqref{1} реализует эффект удержания при единичных коэффициентах. Рассмотрим задачу в упрощенном виде:

\begin{equation}
\begin{split}
&\frac{\partial}{\partial z}\frac{\partial u }{\partial z}  + \frac{\partial u}{\partial z} 
+\frac{\partial }{\partial r} \frac{\partial u}{\partial z}  + \frac{\partial }{\partial r} \frac{\partial u }{\partial r} = 0,
\\
&\frac{\partial u}{\partial r}(0,z)=0, u(r,0) = u_L(r), u(1,z) = 0, \frac{\partial u}{\partial z}(r,1)=0.
\end{split}\label{3}
\end{equation}

Рачет показал, что \eqref{3} отражает не только эффект удержания, но и вогнутость решения вдоль оси $z$ (Рис \ref{ris:im4a}). Условие эллиптичности для задачи \eqref{3} выполняется безусловно.

Проверим выполнение условия эллиптичности задачи \eqref{1} в случае зависимости параметров только от радиуса.
Определитель матрицы главных коэффициентов будет выглядеть:
\begin{equation*}
\begin{split}
&\mathrm{det}\begin{pmatrix}
  \centering
  \frac{T(r)}{\Lambda \kappa (r)}&  \frac{T(r)\zeta (r)}{2} \\
  \frac{T(r)\zeta (r)}{2}& D
\end{pmatrix} = \\ 
&=\frac{T(r)D}{\Lambda\kappa (r)} - \frac{T^2(r)\zeta^2(r)}{4}=\frac{T(r)DR(R_m,r)}{\Lambda (R(R_m,r)-1)}-\frac{T(r)^2 48.746^2 T_e(r)}{4\cdot 600^2r^2}=
\\
&=\frac{D(1.13-5r^2/6)(1.04r^2 +a^2)-0.0197(1.13-5r^2/6)^2(1-5r^2/6)}{1.04\Lambda r^2}
\end{split}
\end{equation*}

Здесь $r\ne 0$ в силу \eqref{2}. Для выполнения условия эллиптичности:
\begin{equation*}
  \begin{split}
D(1.13-5r^2/6)(1.04r^2 +a^2)-0.0197(1.13-5r^2/6)^2(1-5r^2/6) \geq 
\\
0.296(0.562D-0.0197\cdot 0.296 \cdot 0.16=0.296(D-0.0016)>0
\end{split}
\end{equation*}
В таком случае для сохранения эллиптичности уравнения необходимо выполнение условия $D > 0.0016$.

Рассмотрим случай зависимости параметров от $(r,z)$.
\begin{equation*}
  \begin{split}
&\mathrm{det}\begin{pmatrix}
  \centering
  \frac{T(r,z)}{\Lambda \kappa (r,z)}&  \frac{T(r,z)\zeta (r)}{2} \\
  \frac{T(r,z)\zeta (r)}{2}& D
\end{pmatrix}  = \frac{T(r,z)D}{\Lambda\kappa (r,z)} - \frac{T^2(r,z)\zeta^2(r)}{4}=\\
&\frac{T(r,z)DR(R_m,r)}{\Lambda (1-0.5z)(R(R_m,r)-1)} -\frac{T^2(r,z) 48.746^2 T_e(r,z)}{4\cdot 600^2r^2}= (1-0.25z)(1.13-5r^2/6)\\ &\frac{D(1.04r^2 +a^2)-0.0197(1.13-5r^2/6)(1-5r^2/6)(1-0.25z)^2(1-0.5z)}{1.04\Lambda r^2(1-0.5z)}
\end{split}
\end{equation*}
Условие эллиптичности:

\begin{equation*}
\begin{split}
(1-0.25z)(1.13&-5r^2/6)(D(1.04r^2 +a^2) - \\
&-0.0197(1.13-5r^2/6)(1-5r^2/6)(1-0.25z)^2(1-0.5z)) \geq \\
&\geq 0.75\cdot 0.29(0.56D - 0.0197(0.296\cdot 0.16 \cdot 0.5625 \cdot 0.5) = \\&= 0.2175(D-0.00046)>0
\end{split}
\end{equation*}

В случае уточненных параметров для сохранения эллиптичности уравнения необходимо выполнение условия $D > 0.00046$.



\section{Заключение}

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


\bigskip
\begin{thebibliography}{1}


\bibitem{1} A.D. Beklemishev
\href{https://doi.org/10.13182/FST13-A16953}
{\it Helicoidal system for axial plasma pumping in linear traps }, Fusion Sci. Technol, {\bf 263}:1 (2013), 355.

\bibitem{2} A.D. Beklemishev
{\it Radial and axial transport in trap sections with helical corrugation}, AIP Conf. Proc., {\textbf 1771} (2016), 040006.

\bibitem{3} 
A.G. Petrova, V.V. Pukhnachev, O.A. Frolovskaya 
\href{https://doi.org/10.1134/S0081543823040156}{\it Exact Solutions of Second-Grade Fluid Equations}, Proc. Steklov Inst. Math., {\textbf 322} (2023), 173--187.

\bibitem{4}
В.В. Пухначев 
{\it Винтовые течения и их двумерные аналоги}, Всероссийская конференция Математические проблемы механики сплошных сред: тезисы докладов, Новосибирск: СО РАН (2024), 197.

\bibitem{5}A. V. Proskurin, A. M. Sagalakov,
\href{https://doi.org/10.1134/S0869864309030093}{\it Stability of the flow of weakly electro-conductive fluid in the presence of a spiral magnetic field}, Thermophysics and Aeromechanics, {\textbf 16}:3 (2009), 297-406.

\bibitem{6}  
А.В. Проскурин, 
{\it  Устойчивость магнитогидродинамических течений в каналах},Диссертация на соискание ученой степени доктора физико-математических наук. 1.1.9 – механика жидкости, газа и плазмы. (2023).

\bibitem{7} I.S. Chernoshtanov, I.G. Chernykh, G.I. Dudnikova, et al.,
\href{https://doi.org/10.1017/S0022377824000333}{\it . Effects observed in numerical simulation of high-beta plasma with hot ions in an axisymmetric mirror machine}, Journal of Plasma Physics, {\textbf 90}:2 (2024), 905900211.

\bibitem{8} T.V. Liseykina, V.A. Vshivkov, U.A. Kholiyarov
 \href{https://doi.org/10.1134/S1995080224010359}{\it An Efficient Algorithm for Calculating the Magnetic Field in a Cylindrical Plasma Trap}, Lobachevskii Journal of Mathematics, {\textbf 45}:1 (2024), 75-84.

\bibitem{9} T.V. Liseykina, G.I. Dudnikova, V.A. Vshivkov, et al. ,  
 \href{https://doi.org/0.14529/jsfi2303022023}{\it  MHD-PIC Supercomputer Simulation of Plasma Injection into Open Magnetic Trap }, Supercomputing Frontiers and Innovations, {\textbf 10}:3 (2023), 11-17.

\bibitem{10} D.I. Skovorodin, A.D. Beklemishev, 
\href{https://doi.org/10.13182/FST11-A11608}
{\it Numerical Kinetic Model of Axial Confinement in a Mirror Trap}, Fusion Science and Technology {\textbf 59}:1 (2011), 199-201.

\bibitem{11} D.I. Skovorodin, A.D. Beklemishev,
\href{https://doi.org/10.1134/s1063780x12030075}
{\it  Plasma Outflow from a Corrugated Trap in the Kinetic Regime },J. Appl. Mech. Tech. Phys., {\textbf 38}:3 (2012), 202-206.

\bibitem{12} V.T. Astrelin, V.M. Kovenya, T. V. Kozlinskaya,
 \href{https://doi.org/10.1007/s10808-006-0004-5}{\it  Numerical simulation of plasma motion in a magnetic field. Two-dimensional case. }, J. Appl. Mech. Tech. Phys., {\bf 48}:3 (2007), 401-411.

\bibitem{13} K.V. Brushlinskii, V.V. Kryuchenkov, E.V. Stepin,
 \href{https://doi.org/10.1134/S0081543823040053}{\it Mathematical Model of Equilibrium Plasma Configurations in Magnetic Traps and Their Stability Analysis }, Proc. Steklov Inst. Math., {\bf 322} (2023), 52-64. 

\bibitem{14} A.N. Kozlov, V.S. Konovalov,
\href{https://doi.org/10.1134/s2070048223040075}
{\it Empirical Stationary Condition of Two-Dimensional Flows of Ionizing Hydrogen in the Plasma Accelerator Channel },  Mathematical Models and Computer Simulations, {\bf 15} (2023), 630-642.  

\bibitem{15} А.Г. Куликовский, Г.А. Любимов, 
{\it Магнитная гидродинамика} М.: Физматгиз, (1962).

\bibitem{16} К.В. Брушлинский,
{\it Математические и вычислительные задачи магнитной газодинамики}, М.: Бином. Лаборатория знаний, (2009). 

\bibitem{17} К.В. Брушлинский, А.И. Морозов,   
{\it Расчет двумерных течений плазмы в каналах. Вопросы теории плазмы}. Под ред. М.А.Леонтовича. М.: Атомиздат (1974). 

\bibitem{18} Б.Н. Четверушкин, 
{\it Математическое моделирование задач динамики излучающего газа.}, М.: Наука, (1985), 304.

\bibitem{19} И.А. Котельников,
{\it Лекции по физике плазмы. Основы по физике плазмы}  Изательство Лань, 5 изд-е, {\bf 1} (2023), 400.

\bibitem{20} И.А. Котельников
{\it Лекции по физике плазмы. Магнитная гидродинамика} , 4 изд-е, {\bf 2} (2022), 448.
 
\bibitem{21} A. V. Sudnikov, A. D. Beklemishev, V. V. Postupaev, A. V. Burdakov, I. A. Ivanov, N. G. Vasilyeva, K. N. Kuklin, E. N. Sidorov,
\href{https://doi.org/10.1016/j.fusengdes.2017.09.005}
{\it SMOLA device for helical mirror concept exploration}  Fusion Eng. Des. , {\bf 122} (2017), 86-93.

\bibitem{22} D.I. Skovorodin, I.S. Chernoshtanov, V.H. Amirov, et al.
\href{https://doi.org/10.1134/S1063780X23600986}
{\it Gas-Dynamic Multiple-Mirror Trap GDMT} Plasma Physics Reports, {\bf 49}:9 (2023), 1039-1086.

\bibitem{23} Г. Г. Лазарева, И. П. Оксогоева, А. В. Судников  
{\it Влияние температуры на результат математического моделирования переноса плазмы в винтовом магнитном поле.}  Прикладная механика и техническая физика, {\bf 5} (2024).

\bibitem{24} А. А. Самарский, Е. С. Николаев,
{\it Методы решения сеточных уравнений.} Прикладная механика и техническая физика, М.: Наука, 1978. 

\bibitem{25} А. А. Самарский, В. И Мажукин., П. П. Матус, Г. И. Шишкин,
{\it Монотонные разностные схемы для уравнений со смешанными производными} Матем. Моделирование. {\bf 13} (2001). 

\bibitem{26} G.G. Lazareva, I.P. Oksogoeva, A.V. Sudnikov,
\href{https://doi.org/10.31857/S2686954323600945}
{\it Mathematical modeling of plasma transport in a helical magnetic field } Doklady Mathematics {\bf 108}:3 (2023), 533-538 

\bibitem{27} G.G. Lazareva, I.P. Oksogoeva, A.V. Sudnikov,
\href{https://doi.org/10.33048/SIBJIM.2023.26.405}
{\it Influence of mathematical model parameters on plasma transfer in a helical magnetic field} Journal of Applied and Industrial Mathematics {\bf 17}:4 (2023), 750-759. 

\bibitem{28} G.G. Lazareva, I.P. Oksogoeva, A.V. Sudnikov,
\href{https://doi.org/10.1134/S1995080222130248}
{\it Mathematical Modeling of Plasma Transport in a Helical Magnetic Field } Lobachevskii Journal of Mathematics {\bf 43}:10 (202), 2685-2691.

\end{thebibliography}

\end{document}
