Физика плазмы, 2023, T. 49, № 12, стр. 1357-1374

3.1. Уравнение баланса сил для равновесия со свободной границей

Равновесие плазмы в токамаке описывается с помощью уравнения баланса сил, которое можно свести к уравнению в частных производных второго порядка в цилиндрических координатах ($R,Z,\phi $) – к уравнению Грэда–Шафранова. А учитывая симметрию плазмы по тороидальному углу $\phi $, можно написать [31]

где $\psi \left( {R,Z} \right)$ – поток полоидального магнитного поля (нормирован на $2\pi $), шафрановский оператор – $\Delta {\kern 1pt} * = {{R}^{2}}\nabla \frac{\nabla }{{{{R}^{2}}}}$, ${{J}_{\phi }}\left( {R,Z} \right)$ – плотность тороидального тока, которая включает в себя токи, текущие в обмотках магнитной системы и в плазме, $~{{J}_{{\phi ,cond}}}\left( {R,Z} \right)$ – плотность тороидального тока в проводниках. Для условий магнитной системы Глобус-М2 предполагается, что форма полоидальных катушек прямоугольная, где плотность тока распределена равномерно по сечению обмотки: где ${{N}_{{cond}}}$ – количество обмоток полоидального поля, а ${{J}_{{cond,k}}}$, ${{I}_{{cond,k}}}$, ${{S}_{k}}$ и ${{{{\Omega }}}_{{cond,k}}}$ – это плотность тороидального тока, тороидальный ток, площадь полоидального сечения проводника и область пространства k-го проводника соответственно, ${{J}_{{\phi ,pl}}}\left( {R,Z} \right)$ – плотность тороидального тока, имеющая сильную зависимость от потока полоидального магнитного поля, что привносит нелинейность в уравнение Грэда–Шафранова где $g{\text{'}} = \partial g{\text{/}}\partial \psi $ – производная по магнитному потоку, $p\left( \psi \right)$ – давление плазмы, предполагающееся изотропным, $F\left( \psi \right) = R{{B}_{\phi }}$ – функция, определяющая форму тороидального поля.

Задание зависимостей $p{\text{'}}\left( \psi \right)$ и $F\left( \psi \right)F(\psi ){\text{'}}$ осуществляется обычно в виде параметрических функций или полиномов. В коде pyGSS они представлены в двух вариантах, первый из которых задан аналогично используемой зависимости на сферическом токамаке NSTX [32] (вариант А)

где функции заданы от нормализованного потока полоидального магнитного поля ${{\psi }_{a}}$ – поток на магнитной оси, ${{\psi }_{b}}$ – поток на границе. $~{{A}_{0}}$, ${{A}_{1}}$, ${{B}_{1}}$ – свободные параметры, которые подбираются в ходе вычислений под экспериментальные измерения. При этом в данной параметрической функции предполагается, что давление плоское на магнитной оси – $p{\kern 1pt} '\left( 0 \right) = 0$ и что $\left( {FF{\kern 1pt} '} \right){\kern 1pt} '$ = = 0 на границе. Для удобства и повышения стабильности вычислений параметры ${{A}_{0}}$ и ${{A}_{1}}$ выражены через ток плазмы а также через полоидальное бета, определенное как

Из выражений (4), (6), (7), а также выражения для плотности тороидального тока в (1), можно составить систему уравнений, а из решения набор свободных параметров ${{A}_{0}}$, ${{A}_{1}}$, ${{B}_{1}}$ превращается в набор ${{I}_{p}}$, ${{\beta }_{p}}$, ${{B}_{1}}$. Отсюда можно исключить ${{I}_{p}}$, так как он измеряется во всех экспериментах, и, таким образом, свободными параметрами равновесия остаются ${{\beta }_{p}}$ и ${{B}_{1}}$.

Вторым является классический случай (вариант Б)

где свободными параметрами являются ${{C}_{0}}$, ${{C}_{1}}$, ${{\alpha }_{m}}$, ${{\alpha }_{n}}$. Используя плазменный ток и полоидальное бета, как было сделано выше, можно уменьшить количество неизвестных параметров до трех: ${{\beta }_{p}}$, ${{\alpha }_{m}}$, ${{\alpha }_{n}}$.

Во всех расчетах данной статьи использовался вариант А, сравнение двух типов параметрических функций приведено в разделе 0.

3.2. Решение уравнения Грэда–Шафранова

Уравнение Грэда–Шафранова можно решить, как двумерное уравнение Пуассона в тороидальной геометрии или методом функций Грина. В коде pyGSS реализован последний способ. Используя равномерную и прямоугольную вычислительную сетку в координатах $R,Z$ уравнение Грэда–Шафранова (2) можно переписать следующим образом:

где $G\left( {R,Z,R{\kern 1pt} ',Z{\kern 1pt} '} \right)$ – четырехмерная функция Грина, с помощью которой можно получить поток полоидального магнитного поля в точке с координатами $\left( {R,Z} \right)$ от единичного токового элемента, находящегося в точке с координатами $\left( {R{\kern 1pt} ',Z{\kern 1pt} '} \right)$. Данная функция определена как где $K\left( k \right)$ и $E\left( k \right)$ – эллиптические функции первого и второго порядка соответственно.

Уравнение Грэда–Шафранова обладает сильной нелинейностью, поскольку правая часть уравнения (2), т. е. плотность тороидального тока ${{J}_{\phi }}\left( {R,Z} \right)$ зависит от потока полоидального магнитного поля, являющегося искомой функцией. По этой причине, решение уравнения Грэда–Шафранова осуществляется итерационным методом, известным как метод Пикара [33]

Здесь n – номер итерации Пикарда. Таким образом, задавая некоторое начальное распределение потока поля ${{\psi }^{{\left( 0 \right)}}}\left( {R,Z} \right)$ (в коде pyGSS используется Гауссова функция), можно вычислить плотность тороидального магнитного тока J_ϕ(ψ⁽⁰⁾(R, Z)), которая интегрируется с функцией Грина $G\left( {R,Z,R',Z'} \right)$, и дает новое значение потока ${{\psi }^{{\left( 1 \right)}}}\left( {R,Z} \right)$, который передается в следующую итерацию и т.д. Такая процедура продолжается до выполнения условия сходимости – ${{a}_{{tol}}} < {{m}_{1}}$ и ${{r}_{{tol}}} < {{m}_{2}}$, где абсолютная ошибка определяется как

а относительная ошибка равна

С целью оптимизации, перед началом решения уравнения Грэда–Шафранова осуществляется отдельный расчет вакуумного полоидального магнитного потока, создаваемого катушками с помощью предварительно рассчитанных и сохраненных функций Грина. Для этого площадь катушки разбивается на 10⁴ точечных элементов, для каждого высчитывается функция Грина $G\left( {R,Z,{{R}_{i}},{{Z}_{i}}} \right)$ и усредняется. Полученная матрица умножается на количество обмоток ${{n}_{{coil}}}$ и на ток ${{I}_{{coil}}}$ получая, таким образом, поток вакуумного полоидального поля ${{\psi }_{{vac}}}\left( {R,Z} \right)$

3.3. Поиск границы плазмы

Определение границы плазмы требуется на каждом этапе Пикардовских итераций, так как $p{\kern 1pt} '\left( \psi \right)$ и $F\left( \psi \right)F\left( \psi \right){\kern 1pt} '$ должны быть заданы как функции полоидального магнитного потока $\psi $, и они должны быть обнулены вне плазмы.

Поиск границы плазмы начинается с нахождения точек нуль-поля с ${{\left| {\nabla \psi } \right|}^{2}} = 0$, из них можно отличить магнитную ось ($S\left( {R,Z} \right) > 0$) и Х-точку ($S\left( {R,Z} \right) < 0$). Поиск таких точек осуществляется с помощью 2-мерной интерполяции кубическим сплайном [22]

Для получения границы плазмы, происходит поиск контура с потоком одной из Х-точек, внутри которой линии замкнуты.

3.4. Стабилизация решения уравнения Грэда–Шафранова

Численная вертикальная неустойчивость – общая проблема процесса решения уравнения Грэда–Шафранова, когда все токи и параметры, описывающие плазму зафиксированы, а при наличии обмоток полоидального поля, в которых токи направлены в разные стороны, плазма смещается наверх и дальнейшее решение уравнения невозможно. Пример эволюции вертикального положения магнитной оси во время Пикардовских итераций без стабилизации приведен на рис. 2а, где ось неконтролируемо двигается вверх, а после 40-й итерации плазма оказывается вне границ вычислительной сетки и решение уравнения теряет смысл. В связи с этим, в коде pyGSS представлены два метода стабилизации решения.

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

где $\Delta {{\psi }_{{loop}}}({{R}_{j}},{{Z}_{j}})\, = \,{{\psi }_{{loop,sim}}}({{R}_{j}},{{Z}_{j}})\, - \,{{\psi }_{{loop,exp}}}({{R}_{j}},{{Z}_{j}})$ – разница между экспериментальным измерением ${{\psi }_{{loop,exp}}}~$ и рассчитанным значением потока ${{\psi }_{{loop,sim}}}$, $\Delta {{I}_{{wall}}}$ – изменение тока по камере, $\gamma $ – тихоновский параметр для регуляризации. $\sum\nolimits_{i = 1}^{{{N}_{{wall}}}} {G({{R}_{j}},{{Z}_{j}},{{R}_{i}}} $, Z_i)ΔI_wall.i – изменение вклада каждого i-го элемента камеры в поток j-й магнитной петли. Таким образом, на каждом этапе пикардовских итераций осуществляется решение уравнения (16) с помощью тихоновской регуляризации [34] и определяется изменение тока по камере, тогда на следующей итерации ток по камере равен $I_{{wall,i}}^{{\left( n \right)}} = I_{{wall,i}}^{{\left( {n - 1} \right)}} + \Delta {{I}_{{wall,i}}}$. В случае сильного смещения плазменного шнура по вертикали, измерения магнитных петель начинают резко отличаться от рассчитанных кодом, в результате чего, вычисленные по формуле (16) токи по камере стабилизируют плазменный шнур. Пример использования приведен на рис. 2б, где представлено изменение вертикального положения магнитной оси и суммарного тока по камере в ходе пикардовских итераций. Чтобы сделать рассчитанный поток близким к измеренному с помощью магнитных петель, в уравнение добавлена сумма разностей рассчитанного тока по камере ${{I}_{{wall,i}}}$ и экспериментального ${{I}_{{wall - exp,i}}}$ (процедура вычисления описана в разделе 1).

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

где $C$ – множитель плазменного тока ${{I}_{p}}$, $Z_{{mag}}^{{\left( n \right)}}$ – вертикальная координата магнитной оси на $n$-й итерации, ${{Z}_{{target}}}$ – задаваемая пользователем желаемая вертикальная координата магнитной оси. Таким образом, если вертикальная неустойчивость возникает, $Z_{{mag}}^{{\left( n \right)}}$ значительно отклоняется от ${{Z}_{{target}}}$, возникает стабилизирующий ток, который возвращает плазму вниз. Если подобрать ${{Z}_{{target}}}$ близким к ${{Z}_{{mag}}}$, то в конце пикардовских итераций, стабилизирующий ток ${{I}_{{feedback}}} \approx 0$, и не оказывает значительного влияния на глобальное равновесие. К недостаткам такого метода относится необходимость подбора ${{Z}_{{target}}}$ и множителя $C$, чтобы минимизировать стабилизирующий ток. Пример использования приведен на рис. 2в, где представлено изменение вертикального положения магнитной оси и тока в виртуальной обмотке в ходе пикардовских итераций.

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

3.5. Поиск свободных параметров равновесия

Поиск свободных параметров ${{\beta }_{p}}$ и ${{B}_{1}}$ или ${{\beta }_{p}}$, ${{\alpha }_{m}}$, ${{\alpha }_{n}}$ происходит за счет сравнения конечного результата решения уравнения Грэда–Шафранова с экспериментальными данными – измерениями магнитных петель. Плазменный ток и токи в обмотках измеряются в каждом разряде и заданы как неизменяющиеся входные параметры. Сигнал в петле складывается из трех слагаемых

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

Подбор осуществляется следующим образом: подпрограмма, решающая уравнение Грэда–Шафранова, имеет на входе свободные параметры равновесия (а также данные о катушках, токах, функциях Грина и координатах, которые хранятся в файле), а на выходе коэффициент MAPE (средняя относительная погрешность в процентах)

Для минимизации MAPE и подбора свободных параметров равновесия используется метод Байесовской оптимизации [36]. Такой метод выбран был в связи с существованием большого количества несходящихся решений уравнения Грэда–Шафранова и наличием локальных минимумов. В таком случае, обычные методы, например, градиентного спуска, могут попасть в область с неплавным изменением коэффициента MAPE и в результате, выдать неправильный результат. На рис. 3а, приведена погрешность измерений магнитных петель в виде параметра MAPE в процентах в зависимости от искомых в процессе реконструкции равновесия величин – ${{\beta }_{p}}$ и ${{B}_{1}}$. Звездой обозначено наиболее оптимальное значение ${{\beta }_{p}}$ и ${{B}_{1}}$ с MAPE ≈ 0.5%. Как видно, при низких ${{\beta }_{p}}$ и ${{B}_{1}}$ параметр MAPE плавно меняется и падает с приближением к минимуму, который оказывается близок к границе области с несходящимися решениями, где резко потом возрастает, что дополнительно может помешать нестохастическим методам оптимизации. На рис. 3б дополнительно приведена зависимость MAPE от ${{\beta }_{p}}$ при фиксированном ${{B}_{1}} = 2.2$, где видно монотонное падение погрешности при изменении бета ${{\beta }_{p}}$ от 0 до 0.7, а на рис. 3в представлена зависимость MAPE от ${{B}_{1}}$ при ${{\beta }_{p}} = 0.57$.

Для самого наилучшего варианта свободные параметры ${{\beta }_{p}}$ и ${{B}_{1}}$ сохраняются вместе с прочими необходимыми интегральным и производными величинами, а равновесие выводится в geqdsk [37] или divgeo [13] формате.

4.1. Граница плазмы и измерения магнитных петель

Для демонстрации работы кода был выбран разряд #42368 (тороидальное магнитное поле B_T = 0.8 T, ток плазмы I_p = 350 кА), осциллограммы которого изображены на рис. 4. Нагрев плазмы осуществлялся двумя инжекторами НИ-1 и НИ-2, создающими пучки дейтерия с энергией E_b = 28 кэВ и мощностью P_b = 500 кВт (НИ-1), и с энергией E_b = 46 кэВ, мощностью P_b = 925 кВт (НИ-2). Инжекция НИ-2 начиналась на фазе роста тока (150 мс) и длилась 90 мс, инжекция НИ‑1 началась позже, на стадии плато тока при 180 мс и продолжалась 40 мс.

На рис. 1 представлены результаты работы кода pyGSS для стационарной фазы t = 200 мс: распределение полоидального магнитного потока и граница плазмы с нижним положением Х-точки. Для сравнения, здесь также приведена граница плазмы, полученная кодом PET и методом токовых колец. Как видно из рисунка, размеры и конфигурация поперечного сечения плазмы находятся в хорошем соответствии друг с другом при использовании трех разных кодов.

На рис. 5 представлено распределение измеренного магнитными петлями потока полоидального магнитного поля в зависимости от номера датчика (расположение петель изображено на рис. 1), как видно, измерения согласуются с расчетом. При этом вклад от токов по камере оказывается малым, меньше 1 мВб/рад. На рис. 6 приведено изменение абсолютной ${{a}_{{tol}}}$ и относительной ${{r}_{{tol}}}$ ошибки в ходе Пикардовских итераций, как видно, за почти 100 итераций решение уравнения Грэда–Шафранова плавно сходится до предела в 10^–7. При этом, абсолютная ошибка сходится быстрее, чем относительная, опережая последнюю примерно на 40 итераций. Был проведен тест на сеточную сходимость и результаты расчетов идентичны при разрешениях сетки 65 × 65 (основной вариант) в сравнении с 128 × 128, разница вычисленной границы плазмы, энергозапаса, внутренней индуктивности и т.д. не превышает 10^–4.

На данный момент для реконструкции одной магнитной конфигурации требуется от 20 до 100 шагов Байесовской оптимизации, для каждой примерно 100 итераций методом Пикарда. Время расчета одного разряда, разбитого на примерно 100 временных моментов с шагом 1 мс составляет около 40 мин на 24-ядерном процессоре Intel X-eon Gold.

4.2. Сравнение с методом токовых колец

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

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

Таким образом, для корректного сравнения данных вычисления равновесия методом токовых колец и кода pyGSS подходит только граница плазмы, а также следующие геометрические параметры: малый радиус ${{a}_{0}} = \left( {{{R}_{{max}}} - {{R}_{{min}}}} \right){\text{/}}2$, большой радиус ${{R}_{0}} = \left( {{{R}_{{max}}} + {{R}_{{min}}}} \right){\text{/}}2$, вытянутость κ = = $({{Z}_{{max}}} - {{Z}_{{min}}}){\text{/}}2a$ и средняя треугольность δ = = $({{\delta }_{{upper}}} + {{\delta }_{{lower}}}){\text{/}}2$, где верхняя треугольность вычисляется как ${{\delta }_{{upper}}} = ({{R}_{0}} - {{R}_{{upper}}}){\text{/}}2a$, а нижняя δ_tower = $({{R}_{0}} - {{R}_{{lower}}}){\text{/}}2a$, при этом ${{R}_{{upper}}}$ равен большому радиусу самой верхней точки на границе ${{Z}_{{max}}}$, а ${{R}_{{lower}}}$ соответственно большому радиусу точки ${{Z}_{{min}}}$.

На рис. 7 представлены вышеописанные геометрические параметры в зависимости от времени для разряда #42368 (левый столбец), а в правом столбце представлено сравнение на диаграмме для серии экспериментов (разряды: 41 585, 41 615, 41 629, 41 644, 41 645, 41 649, 41 665, 41 666, 42 119, 42 123, 42 343, 42 368), полученных при магнитном поле B_T = 0.8–0.9 Tл, токе плазмы I_p = 300–400 кА, с инжекцией дейтерия или водорода НИ-1 и НИ-2 (E_b = 27–28 кэВ и 40–47 кэВ, суммарная мощность P_b = 1200–1400 кВт). Как видно, малый, большой радиус и вытянутость хорошо соответствуют друг другу с MAPE < 2%, но треугольность плазмы имеет отличия с MAPE = 9%. При этом, наблюдается резкий скачок треугольности с 0.3 до 0.4 во время работы НИ-1+2. Также наблюдаются отличия на фазе роста тока, в это время сильны токи по камере, которые затрудняют реконструкцию равновесия. Визуальное сравнение границы плазмы для одного разряда на рис. 1 также показывает согласованность результатов расчета границы плазмы.

4.3. Сравнение с измерениями положения ноги сепаратрисы с помощью тепловизора

В сферическом токамаке Глобус-М2 установлена быстрая тепловизионная камера, линия наблюдения которой направлена вертикально на нижние диверторные пластины, что позволяет увидеть распределение тепла на них (рис. 8а). К-детектор изготовлен из теллурида кадмия-ртути, что обеспечивает работу в диапазоне длин волн 3.5–4.7 мкм, размер кадра тепловизора 320 × × 256  пикселей в полнокадровом режиме, что соответствует пространственному разрешению ~1.6 мм/пиксель с частотой ~220 Гц.

Для определения положения внешней ноги сепаратрисы используется координата максимума температуры диверторной пластины, которая может отличаться от реального положения из-за процессов диффузии тепла у графитовой плитки, однако погрешность измерения относительно мала из-за узости области с большим потоком тепла на дивертор (менее 3 см, рис. 8б).

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

На рис. 8в представлена эволюция во времени положения ноги сепаратрисы на стороне со слабым полем в разряде #41649, близком по параметрам к #42368, по данным вычислений кодом равновесия pyGSS и измерениям тепловизора, а на рис. 8г представлена диаграмма сравнения для серии разрядов. Как видно, результаты реконструкции равновесия согласуются с экспериментальными данными, полученными оптическим методом, точность реконструкции составляет меньше MAPE = 1%.

4.4. Сравнение с измерениями диамагнитной петли

Глобус-М2 оснащен диамагнитной петлей, расположенной по полоидальному обходу внутри камеры для измерения диамагнитного сигнала, являющегося разностью тороидального потока магнитного поля с плазмой и в ее отсутствии

Данные измерения широко используются в экспериментах для определения энергосодержания и полоидального бета ${{\beta }_{p}}$. Но использование данных измерений сопряжено с трудностями в связи с малостью сигнала на фоне большого потока тороидального поля (доля сигнала 10^–4), в связи с многочисленными наводками, значительным влиянием вибраций, влиянием токов в катушках полоидального и тороидального поля и т.д.

Диамагнитный сигнал не используется непосредственно для реконструкции равновесия кодом pyGSS, но может быть сравнен с сигналом, полученным с помощью кода. На рис. 9а изображена зависимость от времени диамагнитного сигнала, полученного экспериментально и с помощью кодов равновесия pyGSS и PET для разряда #42368. На фазе роста тока плазмы с 140 до 170 мс наблюдаются отличия примерно в 0.5 мВб, вызванные низкой точностью реконструкции равновесия вне плато тока (MAPE = 5–10%). Но в случае квазистационарной фазы разряда (170–240 мс) измерения находятся в хорошем соответствии, как и на диаграмме с сравнением большой статистики экспериментальных измерений и данных реконструкции (рис. 9б), MAPE = 8.6%.

4.5. Сравнение с результатами работы кода равновесия PET

Для реконструкции равновесия на токамаке Глобус-М2 также используется код со свободной границей PET [24], написанный на языке FORTRAN. В коде PET решается уравнение Грэда–Шафранова, но в отличии от кода pyGSS используется иная формула для распределения производной давления $p{\kern 1pt} '\left( \psi \right)$ и функции $F\left( \psi \right)F\left( \psi \right){\kern 1pt} '$

где параметры ${{\alpha }_{0}}$ и ${{b}_{0}}$ определяются через полоидальный ${{\beta }_{p}}$ и ток плазмы ${{I}_{p}}$, а коэффициенты ${{\alpha }_{1}}$ и ${{\alpha }_{2}}$ подбираются пользователем. Ток плазмы и ток в полоидальных катушках заданы на основе экспериментальных измерений.

Реконструкция равновесия в коде осуществляется через подбор неизвестных величин, т.е. полоидального ${{\beta }_{p}}$ и коэффициентов ${{\alpha }_{1}}$ и ${{\alpha }_{2}}$ через сравнение результатов расчета c границей плазмы, определенной методом токовых колец и измеренным сигналом диамагнитной петли. Как показано на рис. 9а, реконструкция равновесия вышеописанным кодом вполне может описать измерения, получаемые диамагнитной петлей.

Результаты реконструкции равновесия кодом PET были использованы для валидации расчетов кодом pyGSS на примере двух величин: полоидального бета ${{\beta }_{p}}$ (7) и для нормализованной внутренней индуктивности

Сравнение результатов, полученных двумя равновесными кодами, представлено на рис. 10 в зависимости от времени для разряда #42368, а также для большой статистики экспериментов. Как видно, данные находятся в хорошем соответствии друг с другом, точность для всех величин в виде коэффициента MAPE немного больше 11%. При этом внутренняя индуктивность, рассчитанная кодом PET, немного больше, чем по коду pyGSS, что говорит о более пикированном профиле плотности тока в PET, вызванном разницей в параметрических функциях, используемых в расчете. На рис. 1 можно увидеть сравнение границы плазмы из PET, близкое к границе из метода токовых колец, и также находящееся в согласии с результатами pyGSS.

4.6. Вычисление энергозапаса

Одним из методов вычисления энергосодержания плазмы является прямое интегрирование равновесного давления плазмы по объему, которое можно получить из кодов равновесия pyGSS и PET

Стоит отметить, что энергозапас, полученный таким методом, включает в себя сумму тепловой энергии (от электронов ${{W}_{e}}$ и от ионов ${{W}_{i}}$), а также быстрых частиц ${{W}_{f}}$ (вкладом от вращения плазмы в данном случае пренебрежем)

На рис. 11а представлено сравнение энергозапаса, полученного кодам pyGSS и PET, точность совпадения результатов вычислений обоих кодов – MAPE = 11.4%. Оба метода используют различные измерения данные магнитных измерений (PET – диамагнитные измерения, а pyGSS – данные магнитных петель), результаты обработки которых, в целом находятся в согласии друг с другом.

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

Вычисление электронного энергозапаса проводилось на основе равновесия pyGSS, а также с помощью профилей электронной плотности и температуры, полученных диагностикой ТР. Для каждой точки ТР производился поиск соответствующего ей нормализованного полоидального потока (5) с последующей интерполяцией профиля по всему объему плазмы. Таким образом, можно получить распределение температуры и концентрации электронов по всему объему плазмы – ${{n}_{e}}(R,Z)$ и ${{T}_{e}}(R,Z)$. Электронный энергозапас можно получить интегрированием по объему распределения температуры и концентрации

Ионный энергозапас был получен приблизительно с помощью отношения максимальной ионной и электронной температур

Результаты расчета кинетического энергозапаса в зависимости от времени представлены на рис. 11б для разряда #42 368, дополнительно, добавлены результаты вычислений кодами PET и pyGSS, как видно все три метода находятся в согласии друг с другом, с учетом погрешности. Измерения CXRS возможны только во время работы пучка нейтральных частиц НИ-1, выступающего здесь также в роли диагностического пучка, поэтому полный кинетический энергозапас доступен только с 182 по 217 мс.

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

4.7. Сравнение использования разных параметрических функций для реконструкции равновесия

При реконструкции равновесия возникает вопрос об использовании разных параметрических функций, которые описывают вид зависимостей $p{\kern 1pt} {\text{'}}\left( \psi \right)$ и $F\left( \psi \right)F\left( \psi \right){\kern 1pt} '$. Как говорилось в разделе 0, в коде pyGSS используются 2 варианта таких функций.

На рис. 12а представлено сравнение границы плазмы и структуры магнитных поверхностей для двух вариантов функций: вариант А – фиолетовый, формула (4), Б – черный, формула (8). Как видно, граница плазмы хорошо совпадает для двух типов функций, при этом статистический анализ на серии экспериментов показал, что малый и большой радиус также совпадают, как и вытянутость, с MAPE < 1%, различия могут наблюдаться только для треугольности, где отличия могут достигать MAPE = 10%. Магнитные поверхности, находящиеся близко к сепаратрисе, находятся в согласии друг с другом, но, начиная с малого радиуса в 0.15 м наблюдаются небольшие отличия, магнитная ось также немного отличается в двух вариантах расчета на 1.2 см, как и полоидальное бета ${{\beta }_{p}}$ и энергозапас также имеют небольшие различия. На рис. 12б и 12в наблюдается общее очень небольшое смещение в сторону реконструкции вариантом Б, показанное красной линией. Следовательно, вариант Б дает чуть большие значения, чем А, но разница может наблюдаться только для случаев с большим энергозапасом, выше 12 кДж. Отличия ${{\beta }_{p}}$ и энергозапаса меньше MAPE < 7%.

Нормализованная внутренняя индуктивность ${{l}_{i}}$ оказывается выше для варианта А, чем в Б, что видно на рис. 12д, что говорит о более пикированном профиле плотности тороидального тока в А. Данный факт подтверждается на рис. 12е, где приведена плотность тороидального тока (формула (3)) в экваториальной плоскости. Плотность тока А оказывается смещена в сторону сильного магнитного поля, а вариант Б в сторону слабого, что и сказывается на форме магнитных поверхностей на рис. 12а. При этом такие существенные различия в форме тока оказывают малое влияние на интегральные параметры и форму границы плазмы.

Параметрическая функция типа А, дает более плоский профиль давления, чем Б, что видно на рис. 12г. Для сравнения, приводятся результаты измерения электронного давления диагностикой ТР (красная область) и ионного давления с помощью диагностики CXRS (синяя область) в экваториальной плоскости

Ионная концентрация, необходимая для расчета ионного давления оценивалась из уравнения квазинейтральности и эффективного заряда, принятого здесь для простоты одинаковым по всему объему плазмы и равным ${{Z}_{{eff}}} = 2$. В данном случае, наиболее близким к суммарному кинетическому давлению (зеленая область), оказалась оценка вариантом функции Б. Но, вместе с этим, разница в энергозапасе, составляет лишь 500 Дж на фоне полного значения примерно 15 кДж, поскольку давление в центре плазмы, где наблюдаются наибольшие различия, дают очень малый вклад в энергозапас, из-за малого объема охватываемого магнитными поверхностями.

Как подтверждается вышеприведенными результатами, реконструкция равновесия с использованием различных параметрических функций не дает принципиально иных результатов. Различия наблюдаются только в виде разной формы профиля давления и тока, что дает немного различающиеся значения полоидального бета ${{\beta }_{p}}$, нормализованной внутренней индуктивности ${{l}_{i}}$ и энергозапаса, отличия меньше MAPE < 7%, а также немного разной формы магнитных поверхностей в центре плазмы.

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