Физика плазмы, 2021, T. 47, № 1, стр. 29-39

Плазмофизический код DOUBLE-MC: моделирование потоков атомов, выходящих из плазмы

М. И. Миронов a*, Ф. В. Чернышев a**, В. И. Афанасьев a***, А. Д. Мельник a****, А. С. Наволоцкий a*****, В. Г. Несеневич a******, М. П. Петров a, С. Я. Петров a

a Физико-технический институт им. А.Ф. Иоффе РАН
С.-Петербург, Россия

* E-mail: maxim@npd.ioffe.ru
** E-mail: fvc@npd.ioffe.ru
*** E-mail: val@npd.ioffe.ru
**** E-mail: amelnik@npd.ioffe.ru
***** E-mail: anavolotsky@mail.ioffe.ru
****** E-mail: vnesenevich@npd.ioffe.ru

Поступила в редакцию 03.05.2020
После доработки 15.07.2020
Принята к публикации 17.07.2020

Полный текст (PDF)

Аннотация

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

Ключевые слова: моделирование, потоки атомов, корпускулярная диагностика, токамак, нейтрализационная мишень, ионная температура

1. ВВЕДЕНИЕ

Моделирование потоков атомов, испускаемых высокотемпературной плазмой, можно использовать в качестве инструмента, помогающего интерпретировать данные корпускулярной диагностики (КД) плазмы.

Метод КД [1, 2] основан на регистрации и анализе потоков атомов, выходящих из плазмы, что осуществляется при помощи специальных энерго-масс-спектрометров – так называемых анализаторов атомов [37]. Потоки атомов могут возникать в результате нейтрализации плазменных ионов при их взаимодействии с частицами основной плазмы (пассивная диагностика), а также с атомами пучков, инжектируемых в плазму (активная диагностика) [8, 9]. Нейтрализация ионов происходит благодаря процессам перезарядки и рекомбинации, которые протекают практически без изменения энергии частиц. Это обстоятельство позволяет, производя анализ покинувших плазму атомов, получать информацию о функции распределения ионов по энергии, а также определять величину таких важных параметров, как ионная температура и компонентный состав плазмы. Между тем на практике не всегда удается установить прямую связь между энергетическими распределениями атомов, выходящих из плазмы, и соответствующими распределениями плазменных ионов. В случае пассивной диагностики это связано прежде всего с тем, что измерения производятся вдоль линии наблюдения, пересекающей различные области плазмы с сильно отличающимися параметрами. В результате, в итоговое энергетическое распределение атомов, регистрируемых анализаторами, вносит вклад множество локальных распределений ионов, не соответствующее какому-либо конкретному месту в плазме. Тем не менее существует подход, позволяющий использовать эти измерения для получения информации о ионах, находящихся в центральной области плазмы. Для этого при анализе энергетических распределений атомов следует исключить из рассмотрения низкие энергии, относящиеся преимущественно к периферийной плазме [2]. Несмотря на то, что этот подход часто используется на практике для оценки центральной ионной температуры, следует помнить, что он применим только для плазмы невысокой концентрации и установок небольшого размера. Это ограничение связано с эффектом экранировки центра плазмы, который заключается в сильном ослаблении потоков атомов в результате их прохождения через периферийные области. Считается, что применение активной диагностики может упростить задачу нахождения локальных функций распределения. Между тем это справедливо только в случае достаточно интенсивной инжекции, когда преобладающим источником атомов является нейтрализация ионов при взаимодействии с пучком, и пассивной составляющей потока атомов по сравнению с активной можно пренебречь. Это позволяет считать, что энергетическое распределение атомов соответствует распределению ионов в месте пересечения линии наблюдения анализатора и траектории пучка. В то же время, если пассивная составляющая потока атомов сопоставима с активной, то такого допущения сделать нельзя. В этом случае создаются дополнительные трудности для интерпретации данных, и связать энергетическое распределение атомов с определенным местом в плазме может оказаться проблематичным.

Для установления связи между энергетическими распределениями ионов и атомов в случаях, когда не работают отмеченные выше упрощающие допущения, может быть использовано компьютерное моделирование. Оно может быть выполнено как с помощью уже имеющихся, так и с помощью специально созданных для этой цели компьютерных кодов. Как правило, существующие коды, которые можно использовать для подобных расчетов, представляют собой сложные математические программы, моделирующие различные процессы, протекающие в высокотемпературной плазме. Чаще всего расчет потоков атомов, выходящих из плазмы, осуществляется в так называемых транспортных кодах [10]. При этом, как правило, подобный расчет является лишь вспомогательной частью этих кодов, нацеленных на решение более общих задач баланса и переноса частиц и энергии в плазме. К сожалению, в большинстве случаев эти математические коды способны производить только общую оценку количества частиц и энергии, уносимых из плазмы потоками атомов. В их задачу не входит моделирование энергетических распределений атомов, испускаемых плазмой в определенном направлении, что можно было бы использовать для интерпретации данных КД. Однако в транспортных кодах заложена возможность расчета такого важного параметра, необходимого для моделирования потоков атомов, как пространственное распределение в плазме концентрации остаточных атомов (нейтралов). Важность информации о плазменных нейтралах для КД обусловлена тем, что в случае пассивной диагностики потоки атомов, выходящие из плазмы, возникают главным образом в результате реакции перезарядки ионов плазмы на этих нейтралах. Экспериментальное определение профиля нейтралов очень затруднительно, а его расчет является достаточно сложной задачей, требующей моделирования проникновения пристеночных атомов в центр плазмы. Это обстоятельство объясняет применение транспортных кодов для подобных расчетов, результаты которых используются для последующего моделирования потоков, регистрируемых анализаторами атомов [11]. Расчет пространственного распределения нейтралов может быть выполнен различными методами. Наиболее известным является метод последовательных приближений, основанный на решении кинетического уравнения для нейтралов и учитывающий поправки, связанные с многократной эстафетной перезарядкой пристеночных атомов [12, 13]. Также в последнее время в связи с развитием компьютерной техники широкое распространение получили алгоритмы, основанные на методе Монте-Карло. Эти алгоритмы позволяют производить расчет проникновения нейтралов без привлечения сложного математического аппарата и учета поправок, связанных с многократной перезарядкой, а также позволяют легко включать в рассмотрение большое количество процессов взаимодействия между частицами.

Необходимо отметить, что существующие транспортные коды, которые можно применять для расчета концентрации нейтралов при интерпретации данных КД, являются очень громоздкими и не предоставляют всей необходимой информации. В то же время относительно небольшой самодостаточный код DOUBLE-MC, которому посвящена настоящая статья, является чрезвычайно востребованным и успешно используется в экспериментах на многих плазменных установках [1420]. Данный код, моделирующий потоки атомов, испускаемые плазмой тороидальных установок, был создан в ФТИ им. А.Ф. Иоффе специально для интерпретации результатов измерений КД. Он является модифицированной версией кода DOUBLE [21], разработанного ранее в ФТИ. Новизна кода DOUBLE-MC по сравнению с его предыдущей версией заключается в том, что в нем используется трехмерная постановка задачи, что позволяет более точно моделировать геометрию эксперимента, а также в использовании метода Монте-Карло, что позволило включить в рассмотрение большее число процессов атомных взаимодействий. Код DOUBLE-MC также с успехом применяется на стадии разработки комплексов КД как для имеющихся плазменных установок, не обладающих подобным диагностическим оборудованием, так и для проектируемых установок, в том числе для токамака-реактора ИТЭР [22].

Изложение материала статьи организовано следующим образом. Разд. 2 посвящен физическим основам, заложенным в код DOUBLE-MC. В этом разделе проанализированы основные элементарные процессы, приводящие к возникновению потоков атомов и к их ослаблению при выходе из плазмы. Даны библиографические ссылки на данные по сечениям процессов атомных столкновений, учитываемых в коде. В разд. 3 приведена структура кода и изложен алгоритм моделирования. Продемонстрировано практическое применение кода для ряда плазменных установок. Произведено сравнение результатов расчета с экспериментальными данными, полученными с помощью КД. В разд. 4 представлены выводы и приведен перечень плазменных установок, на которых код прошел апробацию.

2. ЭЛЕМЕНТАРНЫЕ ПРОЦЕССЫ, УЧИТЫВАЮЩИЕСЯ В КОДЕ DOUBLE-MC

Для расчета потока атомов, покидающих плазму, в коде DOUBLE-MC используется следующее выражение, представляющее собой интеграл вдоль линии наблюдения:

(1)
$\begin{gathered} \Gamma _{0}^{{}}(E) = \frac{1}{{4\pi }}\int\limits_L^{} {\,{{n}_{i}}(x)\sum\limits_j {\left[ {n_{j}^{0}(x){{{\left\langle {\sigma _{j}^{0}v} \right\rangle }}_{{{{v}_{j}}}}}} \right]} } \times \\ \, \times {{f}_{i}}(E,x)\mu (E,x)dx. \\ \end{gathered} $

Здесь E – энергия частиц (выходящих из плазмы атомов и соответствующих им ионов плазмы); ni(x) – концентрация ионов основной плазмы; $n_{j}^{0}(x)$ – концентрация доноров j, дающих вклад в процесс нейтрализации; ${{\left\langle {\sigma _{j}^{0}v} \right\rangle }_{{{{v}_{j}}}}}$ – скоростной коэффициент соответствующей реакции, усредненный по распределению относительной скорости ионов и доноров; fi(E, x) – функция распределения ионов по энергии; μ(E, x) – фактор ослабления потока атомов в результате процессов, приводящих к повторному превращению атомов в ионы; x – координата вдоль линии наблюдения. Донорами могут быть как частицы основной плазмы, так и частицы пучка в случае нейтральной инжекции. Под частицами пучка подразумеваются как атомы самого пучка, так и гало пучка, то есть вторичные атомы и водородоподобные ионы, образующиеся в результате взаимодействия атомов пучка с основной плазмой. Суммирование производится по основным процессам, приводящим к нейтрализации ионов плазмы. Интеграл берется по хорде L, являющейся пересечением линии наблюдения и объема, занимаемого плазмой. Функция распределения ионов по энергии fi(E) имеет изотропный характер и в самом простом случае представляет из себя максвелловское распределение. Также имеется возможность задавать более сложную форму энергетического распределения ионов, обладающих анизотропией.

Выражение для фактора ослабления потока атомов имеет вид

(2)
$\mu (E,x) = \exp \left\{ { - \int\limits_L^{} {\sum\limits_k {\left[ {n_{k}^{i}(l)\frac{{{{{\left\langle {\sigma _{k}^{i}{v}} \right\rangle }}_{{{{{v}}_{k}}}}}}}{{v}}} \right]} } dl} \right\},$
где $n_{k}^{i}(l)$ – концентрация соответствующего компонента плазмы, дающего вклад в процесс k, приводящий к потере атома; ${{\left\langle {\sigma _{k}^{i}{v}} \right\rangle }_{{{{{v}}_{k}}}}}$ – скоростной коэффициент, усредненный по распределению относительной скорости испускаемых атомов и компонента плазмы k; $v$ – скорость частиц (выходящих из плазмы атомов и соответствующих им ионов плазмы); l – координата вдоль линии наблюдения. Суммирование производится по основным процессам, которые вносят вклад в ослабление потока выходящих атомов. Интеграл берется по участку хорды L от координаты x, соответствующей месту образования атомов в плазме, до координаты на ее границе, соответствующей месту выхода атомов из плазмы.

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

Перечень атомных взаимодействий, учитываемых в коде DOUBLE-MC на различных стадиях расчета, приведен в табл. 1. В первом столбце таблицы находится полный список процессов, распределенных по группам согласно типу взаимодействия. Во втором столбце приведены ссылки на соответствующие библиографические источники с указанием раздела или номера страницы. В последних пяти столбцах отмечены стадии расчета, на которых принимались во внимание представленные процессы.

Таблица 1.

Перечень атомных взаимодействий, учитываемых в коде DOUBLE-MC при расчете: А – пространственного распределения водородных и гелиевых нейтралов; Б – образования гало пучка; В – ионизационного баланса водородоподобных ионов примесей; Г – нейтрализации ионов плазмы; Д – ослабления выходящего потока атомов и пучков атомов, инжектируемых в плазму

Тип взаимодействия/Формула, описывающая процесс [Источник], Раздел/страница Стадия расчета, на которой используется
А Б В Г Д
Перезарядка
  H+ + H0 → H0 + H+    [23], A-22  
  He2+ + H0 → He+ + H+    [23], A-88      
  H+ + He+ → H0 + He2+    [23], A-54        
  He2+ + He+ → He+ + He+2    [23], A-84        
  He2+ + He0 → He0 + He2+    [23], A-104  
  H+ + He0 → H0 + He+    [23], A-32    
  He+ + H0 → He0 + H+    [23], A-62        
  H+ + Be3+ → H0 + Be4+    [24], 2801        
  H+ + C5+ → H0 + C6+    [25], 2905        
  C6+ + H0 → C5+ + H+    [26], 1-42      
Рекомбинация
  H+ + e → H0 + γ    [27], 470      
  He2+ + e → He+ + γ    [27], 470        
  Be4+ + e → Be3+ + γ    [27], 470        
  C6+ + e → C5+ + γ    [27], 470        
Ионизация
  e + H0 → H+ + 2e    [28], (1.2.1.)    
  H+ + H0 → H+ + H+ + e    [23], D-6    
  He2+ + H0 → He2+ + H+ + e    [23], D-96      
  C6+ + H0 → C6+ + H++ e    [26], (2-28)    
  e + He0 → He+ + 2e    [29], 5    
  e + He+ → He2+ + 2e    [30], 2377      
  e + Be3+ → Be4+ + 2e    [30], 2377        
  e + C5+ → C6+ + 2e    [30], 2377      
  H+ + He+ → H+ + He2+ + e    [23], D-50        
  H+ + He0 → H+ + He+ + e    [23], D-24      
  He2+ + He0 → He2+ + He+ +e    [23], D-114      

В качестве примечания к таблице отметим, что в коде не рассматриваются атомные процессы с участием возбужденных состояний частиц, так как их учет при расчете величины потока атомов для типичных плазменных параметров приводит к несущественным поправкам, Эти поправки, как правило, находятся на уровне не выше 20–30%, что сопоставимо с неопределенностями сечений основных атомных процессов, участвующих в формировании потока нейтральных частиц. Также на данный момент не рассчитывается образование гало пучков атомов, состоящее из водородоподобных ионов бериллия Be3+.

3. ВЫЧИСЛИТЕЛЬНЫЙ АЛГОРИТМ КОДА

Блок-схема кода DOUBLE-MC представлена на рис. 1. Код содержит несколько последовательно исполняемых этапов. На первом этапе происходит чтение начальных параметров с их последующей верификацией, которая заключается в проверке на квазинейтральность моделируемой плазмы. Далее следует подготовительный этап, во время которого создаются вспомогательные массивы данных, необходимые для расчета магнитных полей. На следующем этапе рассчитывается нейтрализационная мишень, т.е. концентрация доноров, дающих вклад в процесс нейтрализации при формировании выходящего из плазмы потока атомов. Вначале производится моделирование проникновения нейтральных частиц с границы плазмы и образование пассивной мишени. Затем осуществляется вычисление активной мишени, которое включают в себя моделирование проникновения нейтрального пучка и образования локальных нейтрального и ионного гало. Далее рассчитывается мишень, образующаяся в результате рекомбинационных процессов. На следующем этапе работы кода происходит расчет энергетического распределения выходящих атомов. На этом этапе также выполняется вычисление вспомогательных параметров. На заключительном этапе производится сохранение результатов расчета в выходных файлах. Остановимся более подробно на основных моментах работы кода на перечисленных выше этапах.

Рис. 1.

Блок-схема кода DOUBLE-MC.

3.1. Чтение начальных параметров

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

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

(3)
$Y(t) = \left( {{{Y}_{0}} - {{Y}_{1}}} \right){{\left( {1 - {{t}^{{Q1}}}} \right)}^{{Q2}}} + {{Y}_{1}},\quad t = {r \mathord{\left/ {\vphantom {r a}} \right. \kern-0em} a},$
где t – значение приведенного малого радиуса r; Y0, Y1 – значение соответствующего параметра в центре и на границе плазмы; Q1, Q2 – степенные показатели; a – малый радиус плазмы. Допускается задавать профили с помощью комбинации двух таких функций, что может быть использовано для моделирования плазмы с транспортными барьерами. Также имеется возможность задавать тороидальное вращение плазмы.

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

(4)
${{f}_{i}}\left( E \right) = {{f}_{T}}\left( E \right) = \frac{{2\pi }}{{\sqrt {{{{\left( {\pi {{T}_{i}}} \right)}}^{3}}} }}\sqrt E \exp \left( { - \frac{E}{{{{T}_{i}}}}} \right),$
с температурой Ti. Кроме того, существует возможность задавать надтепловую часть функции распределения ионов по энергии. Такая потребность возникает, например, при применении дополнительных методов нагрева плазмы или при образовании в плазме продуктов термоядерного синтеза. В этих случаях распределение ионов имеет комплексную форму, состоящую из тепловой максвелловской части и надтеплового “хвоста” частиц высоких энергий:
(5)
${{f}_{i}}(E) = {{f}_{T}}(E) + {{({{n}_{S}}} \mathord{\left/ {\vphantom {{({{n}_{S}}} {{{n}_{i}}}}} \right. \kern-0em} {{{n}_{i}}}}){{f}_{S}}(E),$
где nS – концентрация надтеплового компонента в плазме,  fS(E) – надтепловой хвост. Надтепловой хвост, в отличие от тепловой части распределения, может быть неизотропным и иметь сложную форму. Предусмотрена возможность задания надтеплового хвоста в виде экспоненциально спадающей с ростом энергии функции, либо в виде зависимости другого произвольного вида.

3.2. Подготовка массивов данных

На этом этапе происходит создание массивов, содержащих градиенты магнитного потока, что может осуществляться двумя способами. Первым способом является формирование массивов на основе карты распределения магнитного потока ψ(R, Z), которая создается внешними программами расчета равновесия плазмы, например, такими как EFIT [31]. Другой вариант предусматривает использование упрощенной конфигурации полоидального поля, задаваемой в виде вложенных магнитных поверхностей, имеющих форму эллипса. В этом случае расчет траекторий заряженных частиц, о котором будет говориться ниже, не производится, а конфигурация полоидального поля используется исключительно для задания распределения плазменных параметров.

3.3. Расчет нейтрализационной мишени

Наиболее трудоемкой частью кода, занимающей основное время вычислений, является этап, связанный с расчетом многокомпонентной пространственно неоднородной нейтрализационной мишени. На этом этапе применяются два различных метода. Расчет той составляющей мишени, которая образуется в результате проникновения в плазму пристеночных атомов, а также атомов, поступающих в плазму в результате нейтральной инжекции, выполняется с помощью метода Монте-Карло [32]. Для расчета мишени, создаваемой в результате рекомбинационных процессов, используются соотношения, полученные из модели ионизационного баланса в корональном приближении [33]. Остановимся более подробно на особенностях расчета с использованием этих двух методов.

Метод Монте-Карло. Расчет проникновения атомов в плазму начинается с запуска пробных частиц с последующим отслеживанием их движения в плазме. Схема расчета одинакова как для атомов, проникающих с границы плазмы, так и для атомов, инжектируемых в плазму в виде пучков. С помощью метода Монте-Карло разыгрываются скорости пробной частицы и частиц плазмы, затем разыгрывается механизм их взаимодействия. После каждого взаимодействия происходит расчет траектории отслеживаемой частицы и разыгрываются параметры следующего взаимодействия. При этом принимается следующая модель плазмы и перемещения в ней атомов и ионов. Плазма считается тороидально-симметричной. Гофрировка магнитного поля, присущая установкам типа токамак, не учитывается, так как она не оказывает существенного влияния на полную концентрацию частиц мишени. При расчете используется цилиндрическая система координат (R, θ, Z), в которых координата R соответствует значению большого радиуса тороидальной плазменной установки, координата θ соответствует тороидальному углу, а координата Z направлена вдоль вертикальной оси установки.

Для нейтральных частиц уравнения движения сводятся к прямолинейному распространению с постоянной скоростью:

(6)
${{\dot {v}}_{R}} = Rv_{\theta }^{2},\quad {{\dot {v}}_{\theta }} = - \frac{1}{R}2vR{{v}_{\theta }},\quad {{\dot {v}}_{Z}} = 0,$
где ${{{v}}_{R}},\;{{{v}}_{\theta }},\;{{{v}}_{Z}}$ – составляющие скорости частицы в системе координат (R, θ, Z), $v$ – полная скорость частицы.

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

(7)
${{{\mathbf{v}}}_{g}} = {{v}_{{||}}}{\mathbf{b}} + \frac{{mv_{ \bot }^{2}{\text{/}}2}}{{q{{B}^{2}}}}({\mathbf{b}} \times \nabla B) + \frac{{mv_{{||}}^{2}}}{{q{{B}^{2}}}}({\mathbf{b}} \times ({\mathbf{b}} \cdot \nabla ){\mathbf{B}}),$
где ${{{v}}_{{||}}},\;{{{v}}_{ \bot }}$ – составляющие скорости частицы вдоль и поперек направления магнитного поля; m, q – масса и заряд частицы; b, B, B – единичный и полный вектор магнитного поля, а также его абсолютное значение. Полоидальные компоненты магнитного поля определяются из карты распределения полоидального магнитного потока ψ(RZ), о способах задания которого говорилось выше,
(8)
${{B}_{R}} = - \frac{1}{R}\frac{{\partial \psi }}{{\partial Z}},\begin{array}{*{20}{c}} {} \end{array}{{B}_{z}} = \frac{1}{R}\frac{{\partial \psi }}{{\partial R}},$
где BR, Bz – радиальная, и вертикальная компоненты магнитного поля. Тороидальная составляющая поля считается спадающей обратно пропорционально величине большого радиуса:
(9)
${{B}_{\theta }} = \frac{{{{B}_{0}}{{R}_{0}}}}{R},$
где B0, R0 – значения тороидального магнитного поля и большого радиуса в центре камеры установки.

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

При проведении расчета используется следующий подход. Для розыгрыша начальных координат, направления и скорости пробных частиц используется штатная процедура компилятора Intel Fortran для генерирования ряда псевдослучайных вещественных чисел в интервале от 0 до 1 с повторяемостью порядка 1018. Подробности устройства алгоритма этого числового генератора изложены в [36, 37]. Этот генератор можно использовать непосредственно для задания равномерно распределенных величин. Однако данный генератор также можно использовать для розыгрыша произвольно распределенных величин с помощью так называемых преобразователей закона случайных чисел. В коде DOUBLE-MC для этих целей применяется метод обратных функций, аналогичный описанному в [38].

При расчете область пространства, занимаемая плазмой, разбивается на ячейки, образующие трехмерный массив. Ячейки имеют фиксированный размер по каждой из координат (R, θ, Z). Во время своего движения в плазме пробные частицы переходят из одной ячейки массива в другую. Для расчета траекторий движения частиц используется схема Рунге–Кутты второго порядка. Применение схем более высокого порядка неоправданно, так как для конкретной задачи не дает выигрыша в точности расчета, но значительно увеличивает время вычислений. На каждом шаге вдоль траектории движения частицы рассчитываются ее весовой коэффициент и вероятность ее преобразования. Весовой коэффициент представляет собой суммарную вероятность того, что частица достигнет рассматриваемой точки в пространстве без полной ионизации, и вычисляется по формуле, аналогичной соотношению (2) для ослабления выходящего потока атомов. Использование в расчетах весовых коэффициентов позволяет в большинстве случаев избежать вычислений для ионизовавшихся частиц, которые не вносят вклад в формирование мишени, а сосредоточиться на отслеживании нейтральных частиц. Под преобразованием подразумеваются в основном те атомные процессы, при которых первоначальная частица может стать полностью ионизованной и выбыть из рассмотрения, а вновь образованный атом или водородоподобный ион займет ее место. В результате произойдет изменение (преобразование) кинетических свойств пробной частицы. Примером такого процесса может служить перезарядка атома водорода на ионе водорода. Если во взаимодействие вступают частицы, отличающиеся атомным номером и массой, то при преобразовании также может произойти изменение типа пробной частицы. Примером такого процесса может служить перезарядка атома гелия на ионе водорода. Вероятность преобразования для каждого возможного процесса определяется как $p = n\left\langle {\sigma {v}} \right\rangle \Delta t$, где n – концентрация частиц мишени, на которых происходит преобразование, $\left\langle {\sigma {v}} \right\rangle $ – скорость реакции, усредненная по локальному энергетическому распределению частиц мишени, а Δt – время, которое тратится на прохождение пробной частицей одного шага вдоль траектории. После того как все вероятности преобразований определены, происходит розыгрыш преобразования пробной частицы. В результате она либо не претерпит изменения, либо будет преобразована в другую частицу. В любом случае, весовой коэффициент частицы будет постоянно уменьшаться и при достижении определенного минимального значения, которое является входным параметром кода, она будет выключена из рассмотрения, и в расчет будет запущена новая пробная частица.

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

На каждом шаге вычислений вклад Δn в результирующую концентрацию мишени будет пропорционален времени Δt, проведенному частицей внутри ячейки, и может быть выражен как

(10)
$\Delta n = \frac{{{{K}_{{norm}}}}}{{NV}}w{{w}_{b}}\Delta t.$

Здесь Knorm – нормировочный множитель, учитывающий начальную концентрацию пробных частиц, N – общее число пробных частиц, V – объем ячейки, w – весовой коэффициент частицы, wb – коэффициент ветвления. Движение пробных частиц отслеживается вплоть до их потери на границе плазмы или достижения минимально допустимого значения весового коэффициента частицы. В результате запуска заданного количества частиц, число которых для плазмы больших установок может достигать 107, формируется трехмерный массив пространственного распределения нейтрализационной мишени.

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

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

(11)
$n_{1}^{0} = \frac{{{{n}_{e}}n_{1}^{i}{{{\left\langle {\sigma v} \right\rangle }}_{{rec}}}\left( {\left( {n_{1}^{i} + n_{2}^{i}} \right)\left\langle {\sigma v} \right\rangle _{{cx}}^{{2 \to 1}} + \Sigma _{{loss}}^{2}} \right)}}{{n_{1}^{i}\left\langle {\sigma v} \right\rangle _{{cx}}^{{2 \to 1}}\Sigma _{{loss}}^{1} + n_{2}^{i}\left\langle {\sigma v} \right\rangle _{{cx}}^{{1 \to 2}}\Sigma _{{loss}}^{2} + \Sigma _{{loss}}^{1}\Sigma _{{loss}}^{2}}},$
(12)
$n_{2}^{0} = \frac{{{{n}_{e}}n_{2}^{i}{{{\left\langle {\sigma v} \right\rangle }}_{{rec}}}\left( {\left( {n_{1}^{i} + n_{2}^{i}} \right)\left\langle {\sigma v} \right\rangle _{{cx}}^{{1 \to 2}} + \Sigma _{{loss}}^{1}} \right)}}{{n_{1}^{i}\left\langle {\sigma v} \right\rangle _{{cx}}^{{2 \to 1}}\Sigma _{{loss}}^{1} + n_{2}^{i}\left\langle {\sigma v} \right\rangle _{{cx}}^{{1 \to 2}}\Sigma _{{loss}}^{2} + \Sigma _{{loss}}^{1}\Sigma _{{loss}}^{2}}},$
(13)
$n_{{H{{e}^{ + }}}}^{{}} = {{n}_{{H{{e}^{{2 + }}}}}} \cdot \frac{{{{{\left\langle {\sigma v} \right\rangle }}_{{rec}}}{{n}_{e}} + {{{\left\langle {\sigma v} \right\rangle }}_{{CX}}}\left( {n_{1}^{0} + n_{2}^{0}} \right)}}{{\Sigma _{{loss}}^{{H{{e}^{ + }}}}}},$
где $n_{1}^{0}$, $n_{2}^{0}$, nHe+ – концентрация атомов водородных компонентов и ионов Не+; ne, $n_{1}^{i}$, $n_{2}^{i}$, nHe2+ – концентрация электронов, ионов водородных компонентов и ионов Не2+; ${{\left\langle {\sigma v} \right\rangle }_{{rec}}}$, ${{\left\langle {\sigma v} \right\rangle }_{{cx}}}$ – скоростные коэффициенты реакций рекомбинации и перезарядки, усредненные по относительной скорости взаимодействующих частиц (т.е. число соответствующих актов в единице объема за единицу времени); Σloss – скорость полных потерь в результате ионизации соответствующих частиц (т.е. число потерь в единицу времени). Индексами 1 → 2 и 2 → 1 обозначены процессы перезарядки, приводящие к превращению частиц первого водородного компонента во второй и обратно.

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

На рис. 2 приведен пример результатов расчета нейтрализационной мишени, выполненного для проектируемого термоядерного реактора ИТЭР [33]. Показаны концентрации отдельных компонентов, составляющих общую мишень для нейтрализации ионов плазмы. При этом формирование водородных и гелиевого компонентов происходит за счет проникновения атомов с края плазмы и рекомбинации с электронами, а концентрация водородоподобных ионов примесей рассчитана в приближении корональной модели. Что касается объемной рекомбинации, вклад которой в нейтрализацию тепловых частиц может оказаться существенным для центральных областей плазмы, то обычно для плазмы ИТЭР он составляет около 20% [22].

Рис. 2.

Радиальное распределение концентрации компонентов мишени, на которых происходит нейтрализация ионов изотопов водорода в установке ИТЭР: водородоподобных ионов углерода C5+, бериллия Be3+ и гелия He+, а также нейтралов дейтерия и трития D0+T0.

3.4. Расчет энергетического распределения выходящих атомов и вычисление вспомогательных параметров

На следующем этапе моделирования происходит расчет энергетического распределения выходящих из плазмы атомов в соответствии с формулой (1).

Как говорилось выше, энергетические распределения, рассчитанные с помощью кода, можно использовать в целях интерпретации данных КД. Для этого необходимо произвести сравнение энергетических распределений, полученных в результате расчета, с распределениями, полученными в ходе эксперимента. Сравнение этих распределений позволяет сделать выводы о правильности задания некоторых параметров, которые вводятся в код на начальном этапе. Таким образом, возникает возможность оценить отклонение тех параметров, величина которых неизвестна из эксперимента, от их реальных значений. К таким параметрам можно отнести концентрацию нейтралов на границе плазмы, ионную температуру, скорость тороидального вращения плазмы, концентрацию надтеплового ионного компонента и др. После коррекции можно повторить расчет и, выполнив несколько раз подобную итерационную процедуру, определить значения этих параметров. Еще раз отметим, что такой подход не претендует на однозначность, но позволяет произвести оценку важных плазменных параметров. На рис. 3 приведено энергетическое распределение атомов, рассчитанное с помощью кода DOUBLE-MC для токамака COMPASS [40]. На рисунке также показан экспериментальный энергетический спектр, полученный с помощью КД. Метод определения ионной температуры по наклону экспериментального спектра, который часто применяется на практике, для случая, представленного на рисунке, дал оценку ионной температуры на уровне 340 эВ. Применение кода позволило путем итераций оценить значение центральной ионной температуры, которое оказалось равным 400 эВ. На рис. 4 показан результат коррекции значения ионной температуры, измеренной на токамаке Глобус-М [41]. Расчет, выполненный с помощью кода DOUBLE-MC, показал, что центральная ионная температура могла достигать величины около 750 эВ, в то время как значение температуры, полученное по наклону экспериментального энергетического спектра, соответствовало уровню 650 эВ.

Рис. 3.

Энергетическое распределение атомов дейтерия, выходящих из плазмы токамака COMPASS: кружки – экспериментальные данные, полученные с помощью КД; сплошная кривая – результат расчета, выполненного с помощью кода DOUBLE-MC в предположении о значении центральной ионной температуры 400 эВ.

Рис. 4.

Временная зависимость ионной температуры на токамаке Глобус-М во время нейтральной инжекции: кружки и сплошная кривая – температура, определенная по наклону экспериментального энергетического спектра; штриховая кривая – результат коррекции, выполненной с помощью кода DOUBLE-MC; ВЧ – период действия нейтральной инжекции.

Помимо энергетического распределения вылетающих из плазмы атомов определяются некоторые вспомогательные характеристики, помогающие разобраться в особенностях образования их потока. Одной из вспомогательных характеристик является так называемая функция светимости, которая пропорциональна вероятности рождения и выхода атома с определенной энергией из данной точки плазмы. Анализ формы функции светимости позволяет получить представления о том, из каких слоев плазмы поступают атомы с той или иной энергией. Пример расчета функций светимости для случая плазмы установки ИТЭР приведен на рис. 5 [33]. Как видно из рисунка, атомы, имеющие более низкую энергию, рождаются намного ближе к внешней границе плазмы, чем атомы с высокой энергией. Это объясняется эффектом экранировки плазмы, который сильнее сказывается на ослаблении потоков атомов низких энергий.

Рис. 5.

Функции светимости атомов дейтерия (штриховые кривые) и трития (сплошные кривые) для разных энергий в плазме установки ИТЭР в зависимости от малого радиуса плазмы. Линия наблюдения лежит в центральной плоскости установки и направлена перпендикулярно плазменному шнуру. Горизонтальной стрелкой указано направление выхода атомов.

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

4. ЗАКЛЮЧЕНИЕ

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

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

Приведем перечень плазменных установок, на которых успешно применялся код DOUBLE-MC:

– TCV (CRPP, Лозанна, Швейцария) [14],

– JET (Culham Science Centre, Абингдон, Великобритания) [15, 21]*,

– Глобус-М (ФТИ им. А.Ф. Иоффе, СПб, Россия) [16, 41],

– ФТ-2 (ФТИ им. А.Ф. Иоффе, СПб, Россия) [17],

– NSTX (PPPL, Принстон, США) [18],

– COMPASS (IPP, Прага, Чехия) [19, 40],

– TJ-II (CIEMAT, Мадрид, Испания) [20]* (в случаях, отмеченных звездочкой *, использовалась предыдущая версия кода).

Также код DOUBLE-MC использовался для оценки величины сигнала диагностики по потокам атомов на установке – ИТЭР (CEA Cadarache, Кадараш, Франция) [22, 33].

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

  1. Афросимов В.В., Гладковский И.П., Гордеев Ю.С., Калинкевич И.Ф., Федоренко Н.В. // ЖТФ. 1960. Т. 30. С. 1456.

  2. Петров М.П. // Физика плазмы. 1976. Т. 2. С. 371.

  3. Афросимов В.В., Березовский Е.Л., Гладковский И.П., Кисляков А.И., Петров М.П., Садовников В.А. // ЖТФ. 1975. Т. 45. С. 56.

  4. Gladkovsky I.P., Izvozchikov A.B., Petrov M.P. // Nucl. Instr. and Meth. 1980. V. 175. P. 441.

  5. Извозчиков А.Б., Петров М.П., Петров С.Я., Чернышев Ф.В., Шустов И.В. // ЖТФ. 1992. Т. 62. Вып. 2. С. 157.

  6. Hayashi K., Hashimoto K., Yamato H., Takeuchi H., Miura Y., Nishitani T., Shiho M., Maeda H. // Rev. Sci. Instrum. 1985. V. 56. P. 359. https://doi.org/10.1063/1.1138303

  7. Davis S.L., Mueller D., Keane C.J. // Rev. Sci. Instrum. 1983. V. 54. P. 315. https://doi.org/10.1063/1.1137366

  8. Афросимов В.В., Иванов Б.А., Кисляков А.И., Пет-ров М.П. // ЖТФ. 1966. Т. 31. С. 89.

  9. Кисляков А.И., Крупник Л.И. // Физика плазмы. 1981. Т. 7. С. 866.

  10. Ongena J.P.H.E., Voitsekhovitch I., Edvard M., McCu-ne D. // Fusion Sci. Tech. 2012. V. 61. P. 180. https://doi.org/10.13182/FST12-A13505

  11. Akers R.J., Appel L.C., Carolan P.G., Conway N.J., Counsell G.F., Cox M., Gee S.J., Gryaznevich M.P., Martin R., Morris A.W., Nightingale M.P.S., Sykes A., Mironov M., Walsh M.J. // Nucl. Fusion. 2002. V. 42. P. 122. https://doi.org/10.1088/0029-5515/42/2/302

  12. Dnestrovskij Yu.N., Lysenko S.E., Kislyakov A.I. // -Nucl. Fusion. 1979. V. 19. P. 293. https://doi.org/10.1088/0029-5515/19/3/002

  13. Днестровский Ю.Н., Костомаров Д.П. Математическое моделирование плазмы. М.: Наука, 1982.

  14. Karpushov A.N., Duval B.P., Schlatter Ch., Afana-syev  V.I., Chernyshev F.V. // Proc. 32nd EPS Conf. on Plasma Phys., Tarragona, 2005. ECA. V. 29C. P–4.097.

  15. Mironov M.I., Afanasyev V.I., Murari A., Santala M., Beaumont P. and JET-EFDA contributors // Plasma Phys. Control. Fusion. 2010. V. 52. 105008. https://doi.org/10.1088/0741-3335/52/10/105008

  16. Gusev V.K., Aleksandrov S.E., Alimov V.Kh., Arkhipov I.I., Ayushin B.B., Barsukov A.G., Ber B.Ya., Chernyshev F.V., Chugunov I.N., Dech A.V., Golant V.E., Gorodetsky A.E., Dyachenko V.V., Kochergin M.M., Kurskiev G.S., Khit-rov S.A., Khromov N.A., Lebedev V.M., Leonov V.M., Litunovstky N.V., Mazul I.V., Minaev V.B., Mineev A.B., Mironov M.I., Miroshnikov I.V., Mukhin E.E., Nikola-ev Yu.A., Novokhatsky A.N., Panasenkov A.A., Patrov M.I., Petrov M.P., Petrov Yu.V., Podushnikova K.A., Rozhansky V.A., Rozhdestvensky V.V., Sakharov N.V., Shcherbinin O.N., Senichenkov I.Yu., Shevelev A.E., Suhov E.V., Trapesnikova I.N., Terukov E.I., Tilinin G.N., Tolstya-kov S.Yu., Varfolomeev V.I., Voronin A.V., Zakharov A.P., Zalavutdinov R.Kh., Yagnov V.A., Kuznetsov E.A., Zhi-lin E.G. // Nucl. Fusion. 2009. V.49. ArtNo: 104021. https://doi.org/10.1088/0029-5515/49/10/104021

  17. Lashkul S.I., Altukhov A.B., Gurchenko A.D., Gusa-kov E.Z., Dyachenko V.V., Esipov L.A., Kouprienko D.V., Chernyshev F.V., Mironov M.I. // Proc. 46th EPS Conf. on Plasma Phys., Milan, 2019. ECA. V. 43B. P4.1080.

  18. https://nstx.pppl.gov/DragNDrop/NSTX_Meetings/Monday_Physics_Meetings/2005/08-01-05/DOUBLE%20Analysis_Medley.pdf

  19. Mitosinkova K., Melnik A., Tomes M., Stockel J., Janky F., Komm M., Imrisek M., Hacek P., Varju J., Weinzettl V. // Proc. 1st EPS Conf. on Plasma Diagnostics, Frascati, 2015. P. 074. https://doi.org/10.22323/1.240.0074

  20. Balbin R., Tabares F., Tribaldos V., Petrov S. and TJ-II Team. // Proc. 30th EPS Conf. on Contr. Fusion and Plasma Phys., St. Petersburg, 2003. ECA. V. 27A. P-1.23.

  21. Afanasyev V.I., Gondhalekar A., Kislyakov A.I. // JET preprint, JET-R-(00)04. Luxemburg: ECSC/EEC/ EURATOM, 1999.

  22. Afanasyev V.I., Mironov M.I., Nesenevich V.G., Petrov M.P., Petrov S.Ya. // Plasma Phys. Control. Fusion. 2013. V. 55. 045008. https://doi.org/10.1088/0741-3335/55/4/045008

  23. Atomic Data for Fusion, V. 1: Collisions of H, H2, He, and Li Atoms and Ions with Atoms and Molecules / Ed. by C.F. Barnett. 1990. ORNL-6086.

  24. Ermolaev A.M., Korotkov A.A. // J. Phys. B: At. Mol. Opt. Phys. 1996. V. 29. P. 2797. https://doi.org/10.1088/0953-4075/29/13/015

  25. Winter T.G. // Phys. Rev. A. 1997. V. 56. P. 2903. https://doi.org/10.1103/PhysRevA.56.2903

  26. Atomic Data for Fusion. V. 5: Collisions of Carbon and Oxygen Ions with Electrons, H, H2 and He / Ed. by R.A. Phaneuf, R.K. Janev, M.S. Pindzola. 1987. O-RNL-6090.

  27. Verner D.A., Ferland G.J. // Astr. J. Suppl. Series. 1996. V. 103. P. 467.

  28. Janev R.K., Smith J.J. Atomic and Plasma-Material Data for Fusion (Supplement to the Journal Nucl. Fusion) / Ed. by R.K. Janev, C. Bobeldijk. Vienna: IAEA, 1993. V. 4.

  29. Ralchenko Yu.,V., Janev R.K., Kato T., Fursa D.V., Bray I., de Heer F.J. NIFS report, NIFS-DATA-059. Toki, 2000.

  30. Aichele K., Hartenfeller U., Hathiramaniet D., Hof-mann G., Schafer V., Steidl M., Stenke M., Salzborn E., Pattard T., Rost J.M. // J. Phys. B: At. Mol. Opt. Phys. 1998. V. 31. P. 2369. https://doi.org/10.1088/0953-4075/31/10/023

  31. Lao L.L., Ferron J.R., Groebner R.J., Howl W., John H.St., Strait E.J., Taylor T.S. // Nucl. Fusion. 1990. V. 30. P. 1035. https://doi.org/10.1088/0029-5515/30/6/006

  32. Hammersley J.M., Handscomb D.C. Monte Carlo Methods. Norwich: Fletcher & Son Ltd, 1964.

  33. Afanasyev V.I., Chernyshev F.V., Kislyakov A.I., Kozlovski S.S., Lyublin B.V., Mironov M.I., Melnik A.D., Nesenevich V.G., Petrov M.P., Petrov S.Ya. // Nucl. Instrum. Methods Phys. Res. A. 2010. V. 621. P. 456. https://doi.org/10.1016/j.nima.2010.06.201

  34. Голант В.Е., Жилинский А.П., Сахаров С.А. Основы физики плазмы. М.: Атомиздат, 1977.

  35. Сорокина Е.А., Ильгисонис В.И. // Физика плазмы. 2012. Т. 38. С. 307.

  36. Lècuyer P. // Communications of the ACM. 1988. V. 31. P. 742. https://doi.org/10.1145/62959.62969

  37. Bratley P., Fox B.L., Schrage L.E. A Guide to Simulation. 2nd Edition. NY: Springer-Verlag, 1983.

  38. Соболь И.М. Численные методы Моте-Карло. М.: Наука, 1973.

  39. Leonov V.M., Zhogolev V.E. // Plasma Phys. Control. Fusion. 2005. V. 47. P. 903. https://doi.org/10.1088/0741-3335/47/6/011

  40. Мельник А.Д., Митошинкова К., Томеш М., Чернышев Ф.В., Штекел Я. // В сб. Тезисы докладов XVI Всероссийской конф. по диагностике высокотемпературной плазмы, Звенигород, 2015. С. 104.

  41. Чернышев Ф.В., Афанасьев В.И., Гусев В.К., Ива-нов А.Е., Курскиев Г.С., Мельник А.Д., Минаев В.Б., Миронов М.И., Несеневич В.Г., Патров М.И., Петров М.П., Петров С.Я., Петров Ю.В., Сахаров Н.В., Толстяков С.Ю. // Физика плазмы. 2011. Т. 37. С. 595.

Дополнительные материалы отсутствуют.