Физика плазмы, 2022, T. 48, № 11, стр. 993-1004

1. ВВЕДЕНИЕ

Для крупных открытых магнитных ловушек, используемых в исследованиях УТС, нагрев плазмы с субтермоядерными параметрами электромагнитным излучением в условиях электронного циклотронного резонанса (ЭЦР) впервые был реализован на установке TMX-U, в которой ЭЦ-нагрев позволил достичь наибольшей на тот момент температуры электронов 280 эВ [1]. Однако в конце 1980-х программа развития открытых ловушек для УТС в Ливерморе была закрыта, и до недавнего времени применение ЭЦР-нагрева плазмы в открытых магнитных ловушках было ограничено или задачами МГД-стабилизации плазмы достаточно низкой плотности [2–4], или компактными лабораторными установками, например, ЭЦР-источники плазмы или ионов. Всплеск интереса именно к нагреву тепловой фракции электронов в условиях ЭЦ-резонанса возник после серии успешных экспериментов по комбинированному нагреву плазмы нейтральными пучками и СВЧ-излучением на российской установке ГДЛ в ИЯФ СО РАН, в которых было продемонстрировано удержание плазмы с рекордными параметрами [5–8]. В результате этих исследований была достаточно убедительно продемонстрирована перспективность использования магнитных ловушек открытого типа с осесимметричной конфигурацией в качестве мощных нейтронных источников для термоядерных приложений [9–12].

Реализация эффективного ЭЦР-нагрева плотной плазмы в большой открытой ловушке потребовала пересмотра сложившихся представлений как о физике взаимодействия СВЧ-излучения с плазмой в режиме ЭЦР [13], так и о последующем транспорте энергии и стабилизации плазменного шнура [8], поскольку ни одна из используемых на практике схем нагрева тороидальной плазмы в этом случае не эффективна. Важную роль при этом играют методы моделирования распространения и поглощения электромагнитных волн в неоднородной плазме осесимметричной ловушки. До недавнего времени такое моделирование проводилось исключительно в рамках приближения геометрической оптики методом трассировки лучей, например, именно таким способом была обоснована основная схема ЭЦР-нагрева, используемая в ГДЛ [13–15]. В этой схеме возникают области отражения и поглощения волн, в которых среда перестает быть плавно-неоднородной в масштабах длины волны, поэтому детальное понимание физики происходящих процессов требовало выхода за рамки приближения геометрической оптики.

Основные эффекты, приводящие к нарушению приближения геометрической оптики, связаны с пространственной дисперсией в сильно неоднородной области резонансного поглощения, с дифракционным расплыванием волнового пучка вне этой области и с формированием каустик в окрестности точек полного внутреннего отражения. Прямой учет этих эффектов в рамках уравнений Максвелла сложен, в первую очередь, из-за малой длины волны излучения по сравнению с размерами установки¹¹ и, как правило, нецелесообразен, поскольку СВЧ-излучение распространяется в плазме в виде хорошо сформированных квазиоптических (но не обязательно гауссовых) пучков. Хорошей альтернативой этому подходу является последовательный квазиоптический подход, основанный на асимптотическом разложении уравнений Максвелла в параксиальном приближении в окрестности выделенной геометро-оптической моды [16, 17]. Этот метод был адаптирован для описания распространения и поглощения волновых пучков в открытых ловушках [18, 19]. Заметим, что специфика СВЧ-нагрева в современных прямых ловушках стимулировала дальнейшее развитие квазиоптической теории, изначально предложенной для тороидальных систем. Это связано с тем, что в прямой ловушке магнитное поле в основном меняется вдоль своего направления, а в тороидальной – поперек. Поэтому приближения и численные приемы, которые хорошо проявили себя при моделировании распространения волн в тороидальных установках [20], не сработали в случае прямой ловушки. Это потребовало создания уточненной теоретической модели [21] и существенной модификации численных алгоритмов, что привело в итоге к созданию абсолютно нового кода QOOT [18]. Данный код применялся для моделирования не только ЭЦР-нагрева плазмы в установке ГДЛ [18, 19], но и коллективного рассеяния микроволнового излучения в нерезонансной плазме на этой же установке [22–24]. Во всех приложениях кода пространственные распределения концентрации и температуры электронов в установке предполагались заданными. Таким образом, решалась несамосогласованная задача – по известным из эксперимента квазистационарным профилям плазмы анализировалось распределение СВЧ-поля и удельный энерговклад за счет диссипации волн в плазме.

В данной статье сделан естественный следующий шаг в развитии квазиоптического моделирования – реализован самосогласованный расчет распределения температуры электронов, неизвестного заранее и устанавливающегося в результате баланса ЭЦР-нагрева и потерь энергии. При этом мы воспользовались тем обстоятельством, что, в отличие от тороидальных ловушек, основной канал потерь энергии для установок типа ГДЛ хорошо изучен – это газодинамический вынос мишенной плазмы вдоль силовых линий магнитного поля [25]. В результате мы получили относительно простой инструмент, позволяющий описывать нестационарную динамику электронной фракции в процессе СВЧ-нагрева. Учитывая импульсный характер разрядов в большинстве крупных современных установок, этот инструмент был давно востребован, однако, на наш взгляд, не мог быть адекватно реализован в рамках лучевых моделей. На основании нашего опыта это связано с тем, что картина лучей даже в статическом случае может быть довольно запутанной и нерегулярной, учет изменения среды за счет ЭЦР-нагрева привносит дополнительные нефизические эффекты, связанные с заменой непрерывного распределения поля в пространстве конечным числом бесконечно тонких лучей. С другой стороны, квазиоптический метод, с самого начала ориентированный на расчет плавного непрерывного распределения комплексной амплитуды волнового поля, представляется вполне адекватным для описания сложной самосогласованной динамики поля и плазмы, сопровождаемой, например, развитием “перегревной” неустойчивости, обсуждаемой в конце этой статьи.

Поскольку такого рода вопросы возникали к нашим предыдущим публикациям²², сразу подчеркнем, что при интенсивном СВЧ-нагреве мы не учитываем возможное отклонение функции распределения электронов по скоростям от максвелловского распределения. Этот подход оправдан тем обстоятельством, что в наиболее важном практическом случае нагрева достаточно плотной плазмы энергия СВЧ-поля вкладывается в тепловую компоненту, для которой кулоновское взаимодействие эффективно поддерживает распределение электронов по скоростям, близкое к максвелловскому. В частности, такие распределения, формирующиеся сразу после ЭЦР-нагрева, были прямо зарегистрированы методом лазерного рассеяния на установке ГДЛ [5–7]. Кинетические эффекты, связанные с модификацией функции распределения резонансных электронов, проявляются при нагреве более разреженной плазмы, например, в установке GAMMA-10 [3] или на стадии ЭЦР-пробоя в ГДЛ [26]. Такие эффекты описываются в рамках кинетического уравнения типа Фоккера–Планка для функции распределения электронов, учитывающего совместное действие кулоновских соударений и СВЧ-поля в квазилинейном приближении [27]. Для описания распределения СВЧ-поля в разреженной плазме обычно достаточно геометро-оптического приближения, точнее уравнения переноса СВЧ-излучения вдоль луча с коэффициентом поглощения, отвечающим возмущенной функции распределения электронов.

Статья организована следующим образом. В разд. 2 предложена транспортная модель, описывающая баланс энергии электронов при ЭЦР-нагреве плазмы в газодинамической ловушке. В разд. 3 мы кратко напоминаем основные идеи квазиоптического подхода и обсуждаем детали, связанные с одновременным численным решением квазиоптического уравнения для СВЧ-поля и уравнения баланса энергии для температуры электронов. В разд. 4 обсуждаются первые результаты моделирования перспективных сценариев ЭЦР-нагрева плазмы на установке ГДЛ. В разд. 5 сформулированы основные результаты.

2. БАЛАНС ЭНЕРГИИ ПРИ ЭЦР-НАГРЕВЕ ПЛАЗМЫ В ГАЗОДИНАМИЧЕСКОЙ ЛОВУШКЕ

Открытые магнитные ловушки в зависимости от физических задач и диапазонов параметров плазмы могут работать в различных режимах удержания плазмы [25]. С точки зрения авторов статьи, для термоядерных приложений наиболее перспективным и интересным является так называемый режим газодинамического удержания, реализованный, в частности, в аксиально-симметричной ловушке ГДЛ в ИЯФ СО РАН. В плазме газодинамической магнитной ловушки формируются две фракции – энергичные ионы, возникающие вследствие процессов перезарядки и ионизации вводимых в плазму пучков нейтральных частиц с энергиями порядка десятков кэВ, и так называемая мишенная плазма с умеренными “дотермоядерными” параметрами: температура электронов и ионов порядка сотен эВ. Высокоэнергичные частицы, вводимые в магнитную ловушку под углом к магнитной оси, удерживаются в её центральной секции адиабатически за счет сохранения магнитного момента. Высокоэнергичные ионы испытывают кулоновские соударения с электронами плазмы и постоянно теряют на них свою энергию, соударения с ионами мишенной плазмы приводят в основном к изотропизации распределения энергичных ионов по скоростям. Мишенная плазма, основное предназначение которой состоит в подавлении высокочастотных неустойчивостей фракции энергичных частиц, удерживается в газодинамическом режиме.

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

где ${{T}_{e}}$ и ${{T}_{i}}$ – температуры электронов и ионов, соответственно, Z и M – средний заряд и масса ионов, для простоты при моделировании ЭЦР-разрядов ионы будем считать холодными. Связанная с этим потоком мощность продольных потерь энергии электронной компоненты плазмы ${{\Delta }}P$ из магнитной силовой трубки, отвечающей в центральном сечении ловушки элементу площади $\Delta S$, равна

Множитель два учитывает наличие двух пробок, ограничивающих силовую трубку, $\Delta S{\text{/}}R$ – сечение силовой трубки в магнитной пробке, R = $ = {{B}_{{\max }}}{\text{/}}{{B}_{{\min }}}$ – пробочное отношение, A ~ $\sim \ln \sqrt {M{\text{/}}m} $ – безразмерный коэффициент, определяющий энергию на одну электрон-ионную пару, вылетающую из ловушки в расширитель, m – масса электронов, ${{n}_{e}}$ – концентрация электронов в центре ловушки, больцмановский множитель $\exp \left( { - \Delta \phi {\text{/}}{{T}_{e}}} \right) \approx 1{\text{/}}\sqrt {\text{e}} $ учитывает отношение концентраций электронов в пробке и в центре ловушки и соответствует перепаду амбиполярного потенциала $\Delta \phi \approx \tfrac{1}{2}Mv_{s}^{2}$ (при условии $Z{{T}_{e}} \gg {{T}_{i}}$), необходимому на ускорение ионов до ионно-звуковой скорости (здесь е – основание натурального логарифма). Последние теоретические и экспериментальные исследования физики расширителя дают для установок масштаба ГДЛ значение $A \approx 7$ в дейтериевой плазме [28–30]. Заметим, что формирование продольного амбиполярного термобарьера полностью подавляет потери энергии, связанные с классической “спитцеровской” электронной теплопроводностью [5], при этом для мишенной плазмы внутри ловушки хорошо выполняется изотермическое приближение. Рассматривая динамику электронной температуры, разумно предполагать, что во всем объеме плазменного шнура вплоть до лимитера конвективные газодинамические потери (2) доминируют по сравнению поперечным переносом [8].

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

• Основная энергия вносится захваченными в плазме пучками атомов в виде кинетической энергии высокоэнергичных ионов.

• Высокоэнергичные ионы испытывают силу кулоновского трения (collisional drag force) в мишенной плазме, что приводит к нагреву электронов.

• При включении дополнительного ЭЦР-нагрева плазмы внешним СВЧ-излучением энергия вкладывается непосредственно в тепловые электроны.

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

• Принимая специальные меры, все остальные каналы, включая поперечный транспорт энергии за счет развития МГД-неустойчивостей, удается подавить настолько, что их вкладом можно пренебречь на фоне мощного продольного выноса [9–11].

Для установки ГДЛ эта схема была неоднократно подтверждена в экспериментах [5–8, 30].

Отвечающее данной схеме уравнение баланса энергии, запасенной в равновесных электронах в заданной силовой трубке $\Delta S$, можно записать так

Первое слагаемое в правой части описывает передачу энергии от энергичных ионов за счет кулоновских соударений, второе – определяет нагрев электронов в результате поглощения СВЧ-волн, третье – соответствует газодинамическим потерям (2). Здесь ${{\nu }_{{i{\text{e}}}}}$ – транспортная частота кулоновских соударений пробного энергичного иона на электронах мишенной плазмы, $\left\langle {{{n}_{{hi}}}} \right\rangle $ – усредненная по длине силовой трубки концентрация энергичных ионов, ${{E}_{{hi}}}$ – средняя энергия энергичных ионов, ${{P}_{{ECR}}} = {{Q}_{{abs}}}{\text{/}}\Delta S$ – поглощаемая СВЧ-мощность, приходящаяся на единицу площади $\Delta S$ центрального поперечного сечения, то есть “поверхностная плотность” вкладываемой СВЧ-мощности, $L = {{B}_{{{\text{min}}}}}\int {dl{\text{/}}B} $ приведенная длина силовой трубки.

Радиальная неоднородность амбиполярного потенциала приводит к E × B-вращению плазмы с азимутальной дрейфовой скоростью, пропорциональной градиенту электронной температуры $\nabla {{T}_{e}}$ [31]. Характерный период вращения плазмы в приосевой области ГДЛ около 10 мкс [8]; на больших временах можно говорить о формировании магнитных (точнее, магнитодрейфовых) поверхностей в аксиально-симметричной ловушке. При этом все величины усредняются по силовой трубке за счет быстрого движения частиц вдоль магнитного поля и по всем азимутам из-за вращения плазмы. Поэтому уравнение (3) следует сразу рассматривать для совокупности силовых линий, образующих магнитную поверхность с меткой ρ – все величины в (3) являются функциями ρ и времени t. В качестве метки магнитной поверхности $\rho $ возьмем ее радиус в центральном сечении, тогда элемент площади в (3) можно определить как $\Delta S = 2\pi \rho \,\Delta \rho $. Будем также считать, что вся СВЧ‑мощность, независимо от того, в каком месте она была вложена, усредняется по всей магнитной поверхности ρ. В плотной субтермоядерной плазме газодинамической ловушки эффекты пространственной дисперсии и резонансной диссипации заметно влияют на распространение и поглощение греющего СВЧ-поля, поэтому при расчете профиля ${{P}_{{ECR}}}(\rho ,t)$ мы не можем воспользоваться приближением заданного поля, и вынуждены решать самосогласованную задачу, сформулированную ниже.

Уравнение (3) можно привести к более простому виду, если предположить, что концентрация электронов, устанавливающаяся в основном за счет баланса продольных потерь и процессов ионизации, не меняется в процессе ЭЦР-нагрева. Это упрощающее предположение подтверждается сравнением полученных в ходе экспериментов на ГДЛ-профилей концентрации плазмы в отсутствии и при наличии дополнительного ЭЦР-нагрева, а также измеренной динамикой концентрации плазмы сразу после выключения дополнительного ЭЦР-нагрева [6, 8]. В результате мы получаем следующее уравнение:

Для вычисления $\left\langle {{{n}_{{hi}}}} \right\rangle $ необходимо знать функцию распределения энергичных ионов, которая может быть восстановлена, например, путем решения уравнения Фоккера–Планка [32]. Мы используем другой прием, основанный на экспериментальных данных. Допустим, что нам известен радиальный профиль распределения электронной температуры $T_{e}^{*}(\rho )$ в стационарном разряде, поддерживаемом только атомными пучками без дополнительного ЭЦР-нагрева. В этом случае из (4) следует равенство

где “*” означает стационарные значения величин до включения ЭЦР-нагрева. Считая, что параметры распределения энергичных ионов не изменились при включении ЭЦР-нагрева, пренебрегая температурой мишенных ионов по сравнению с электронами и учитывая зависимость частоты соударений от температуры электронов ${{\nu }_{{ie}}}{\text{/}}{{n}_{e}} \propto T_{e}^{{ - 3/2}}$, из уравнений (4) и (5) получаем

Это окончательное уравнение, описывающее изменение температуры электронов ${{T}_{e}}(\rho ,t)$ на магнитной поверхности ρ; ${{n}_{e}}(\rho )$ и $T_{e}^{*}(\rho ) = {{T}_{e}}(\rho ,0)$ считаются заданными из эксперимента или по результатам других моделей; профиль поглощаемой мощности ${{P}_{{ECR}}}(\rho ,t)$ рассчитывается с помощью квазиоптического кода QOOT для заданного мгновенного профиля ${{T}_{e}}(\rho ,t)$. Таким образом, уравнения (6), формально полученные независимо для каждой магнитной поверхности, оказываются связанными через общий профиль поглощения СВЧ-мощности, который в свою очередь, сложным нелокальным образом зависит от профиля температуры.

3. ОБЪЕДИНЕНИЕ КВАЗИОПТИЧЕСКОЙ И ТРАНСПОРТНОЙ МОДЕЛЕЙ

Самосогласованный квазиоптический расчет был реализован в новой версии кода QOOT–SC (Quasi-Optics for Open Traps – Self-Consistent). Здесь мы обсудим некоторые технические детали, связанные с совместным численным решением уравнения баланса энергии электронов и квазиоптического уравнения для СВЧ-поля.

Расчет профиля поглощенной мощности производится на основе квазиоптической модели, детально описанной в [18, 19]. В качестве эволюционной координаты мы используем координату z вдоль оси ловушки, в качестве поперечных координат – декартову систему ${{{\mathbf{r}}}_{ \bot }} = (x,y)$ в плоскости поперек оси ловушки. Поиск распределения волнового поля в объеме ловушки сводится к решению эволюционного уравнения вида

для скалярной комплексной амплитуды $u({{{\mathbf{r}}}_{ \bot }},z)$ с нелокальным эволюционным оператором, где ${{k}_{0}} = \omega {\text{/}}c$, а ${{\hat {H}}_{H}}$ и $i{{\hat {G}}^{2}}$ представляют собой, соответственно, эрмитову и антиэрмитову (положительно определенную для диссипативной среды) части, действие которых определено как

В параксиальном приближении интегральное ядро оператора может быть восстановлено из решения ${{k}_{z}}({{{\mathbf{r}}}_{ \bot }},{{{\mathbf{k}}}_{ \bot }},z)$ локально-однородного геометрооптического дисперсионного соотношения для выбранной моды в теплой магнитоактивной плазме: с точностью до аддитивной постоянной $H = {{k}_{z}}({{{\mathbf{r}}}_{ \bot }},{{{\mathbf{k}}}_{ \bot }},z){\text{/}}{{k}_{0}}$ (подробнее см. формулы (18) и (19) в [18]). В диссипативной среде мнимая часть этого решения определяет антиэрмитову часть эволюционного оператора, которая в свою очередь определяет объемную плотность поглощаемой в среде СВЧ-мощности $\Pi = 2{{k}_{0}}\operatorname{Re} (u{\kern 1pt} *\hat {G}[\hat {G}[u]])$.

Как уже отмечалось, при ЭЦР-нагреве плазмы в магнитной ловушке вкладываемая мощность быстро перераспределяется по магнитным поверхностям, поэтому в уравнение баланса энергии входит одномерное распределение мощности по магнитным поверхностям, определяемое как

где $l(\rho ) = \iint {\delta \left( {\rho - \rho (x,y,0)} \right)dxdy}$ определяет эффективный периметр магнитной поверхности в центральном сечении. В длинной аксиально-симметричной ловушке типа ГДЛ магнитные поверхности задаются приближенными уравнениями, учитывающими сохранение потока (вмороженность) магнитного поля

Здесь $z = 0$ отвечает центру ловушки (минимуму магнитного поля), ρ имеет смысл эффективного радиуса магнитной поверхности в центральном сечении. Для этого случая $l(\rho ) = 2\pi \rho $.

Исходный квазиоптический код QOOT позволяет рассчитать и визуализировать решение уравнения (7), а также построить профиль вкладываемой СВЧ-мощности (10) для открытых магнитных конфигураций с теплой (нерелятивистской) плазмой. На входе задаются:

• магнитная конфигурация ${\mathbf{B}}(x,y,z)$ и $\rho (x,y,z)$,

• распределение поля в квазиоптическом пучке на входе в плазму,

• радиальные профили температуры ${{T}_{e}}(\rho )$ и концентрации ${{n}_{e}}(\rho )$ электронов.

На выходе получаются

• радиальный профиль плотности поглощаемой мощности ${{P}_{{ECR}}}(\rho )$,

• распределение интегральной интенсивности поля, например, $J(y,z) = \int {{\text{|}}u{{{\text{|}}}^{2}}} dx$ для визуализации поля в плоскости $(y,z)$.

В расчетах используется гладкая непрерывная функция $\Pi (x,y,z)$, полученная при помощи билинейной интерполяции значений поглощенной мощности на внутренней вычислительной сетке. При этом интеграл (10) берется аналитически, что дает гладкое распределение ${{P}_{{ECR}}}(\rho )$ поглощенной мощности по магнитным поверхностям.

При построении решения уравнения баланса энергии электронов (6) квазиоптический код QOOT приходится вызывать на каждом шаге по времени в соответствии с меняющимся профилем температуры ${{T}_{e}}(\rho ,t)$. Учитывая сложность вычисления обратного оператора для правой части, мы применили явную разностную схему по времени

где $R[{{T}_{e}}(\rho ,{{t}_{i}})]$ обозначает правую часть уравне-ния (6), рассчитанную с помощью квазиоптического кода в момент времени ${{t}_{i}}$. Уравнение (12) позволяет найти ${{T}_{e}}(\rho ,{{t}_{{i + 1}}})$ по ${{T}_{e}}(\rho ,{{t}_{i}})$, сделав одну итерацию квазиоптического расчета ${{P}_{{ECR}}}(\rho )$. Шаг по времени используется адаптивный, при котором максимальное по $\rho $ изменение температуры в ходе итераций не превышало бы определенного порогового значения $\Delta {{T}_{{\max }}}$, а сам шаг ${{t}_{{i + 1}}} - {{t}_{i}} < \Delta {{t}_{{\max }}}$ ограничен некоторым естественным масштабом времени, известным из эксперимента. Этот подход позволяет воспроизвести процесс установления профиля температуры, но при этом требует большого времени для расчета − 20–80 отдельных вызовов кода QOOT для приведенных ниже примеров.

Для более оперативного анализа и проверки результатов динамических расчетов было отдельно реализовано вычисление стационарного распределения ${{T}_{e}}(\rho )$, устанавливающегося при $t \to \infty $. Расчет самосогласованного профиля температуры, обеспечивающего $R\,[{{T}_{e}}(\rho )] = 0$, проводится по следующей итеративной схеме:

Иными словами, новое значение температуры ${{T}_{e}}$ для каждой магнитной поверхности находится как решение статического уравнения баланса температуры с профилем вкладываемой мощности ${{P}_{{ECR}}}(\rho )$, подсчитанным на предыдущем шаге. Эта последовательность повторяется до тех пор, пока среднее квадратичное изменение температуры за одну итерацию не окажется меньше выбранного порогового значения $\Delta {{\bar {T}}_{{\max }}}$. Этот подход требует меньшего числа вызовов QOOT по сравнению с динамическим расчетом (10−30 вызовов для примеров, приведенных ниже). При большом числе итераций (шагов по времени) решение, полученное с помощью алгоритма (12), сходится к решению стационарной задачи (13).

В реализованной версии кода QOOT–SC профиль концентрации ${{n}_{e}}(\rho )$ и начальный профиль температуры электронов $T_{e}^{*}(\rho )$ задаются аналитическими выражениями. Эволюция профиля температуры в ходе нагрева описывается на эквидистантной сетке ${{\rho }_{i}}$ по радиальной координате. Для реализации гладкого пересчета температуры на узлы внутренней двумерной сетки по (x, y), используемой в квазиоптическом блоке QOOT, значение температуры между узлами внешней сетки ${{\rho }_{i}}$ интерполируется с помощью кубического сплайна.

4. МОДЕЛИРОВАНИЕ РЕЖИМОВ “ШИРОКОГО” И “УЗКОГО” ЭЦР-НАГРЕВА НА УСТАНОВКЕ ГДЛ

В этом разделе мы представим первый пример применения новых возможностей интегрированного кода QOOT–SC для физического моделирования сложных режимов СВЧ-нагрева плазмы на установке ГДЛ. В экспериментах на этой установке было обнаружено, что при различных, но весьма близких конфигурациях магнитного поля разряд в условиях комбинированного ЭЦР-нагрева плазмы и нагрева нейтральными пучками может протекать по двум совершенно различным сценариям [5–7]. Реализовывался либо “узкий профиль нагрева” – неустойчивый режим с рекордно высокой температурой электронов на оси ловушки, либо “широкий профиль нагрева” – относительно устойчивый режим с умеренно повышенной электронной температурой в протяженной в радиальном направлении области плазменного шнура (см. рис. 1). Предпринятое ранее моделирование профилей поглощения СВЧ-мощности не позволяло объяснить драматическую разницу в результатах экспериментов. Новая гибридная модель позволила впервые воспроизвести бифуркацию режимов с узким и широким нагревом в численном эксперименте.

Аналитические формулы для профилей концентрации и начальной температуры электронов, интерполирующие экспериментальные данные для установки ГДЛ имеют вид

Параметры этих распределений, соответствующие приведенным на рис. 1 профилям до включения ЭЦР-нагрева, следующие:

• для режима “широкого” нагрева: ${{n}_{0}} = 1.44 \times {{10}^{{13}}}$ см^–3, ${{\rho }_{0}} = 15.8\,$ см, ${{T}_{0}} = 120$ эВ, $a = 3.5$, $b = 2.9$, $с = 0.38$, $\Delta {{\rho }_{1}} = 9$ см, $\Delta {{\rho }_{2}} = 19$ см;

• для режима “узкого” нагрева: ${{n}_{0}} = 0.9 \times $ × 10¹³ см^–3, ${{\rho }_{0}} = 9.5$ см, ${{T}_{0}} = 170$ эВ, $a = 3.5$, $b = 0.42$, $c = 0.4$, $\Delta {{\rho }_{1}} = 3.3$ см, $\Delta {{\rho }_{2}} = 17.8$ см.

Магнитная конфигурация установки ГДЛ моделировалась системой соосных круглых катушек, параметры которых для двух моделируемых режимов приведены в приложении.

Переключение между “широким” и “узким” нагревом контролировалось в эксперименте током через две выделенные “конусные” катушки ($z = \pm 309$ см), расположенные ближе всего к области ЭЦР-поглощения. На рис. 2 приведено взаимное положение входного порта и поверхностей циклотронного резонанса для режимов “широкого” и “узкого” нагрева в боковой проекции ГДЛ. Для сравнения приведена также визуализация квазиоптических пучков, соответствующих экспериментальным профилям температуры и плотности электронов до включения ЭЦР-нагрева, приведенным на рис. 1. На рис. 3 изображены соответствующие этим профилям распределения плотности поглощаемой СВЧ-мощности по радиусам для двух режимов нагрева. Профили в целом подобны, главная особенность заключается в том, что в режиме с “широким” нагревом заметный энерговклад приходится на область 5–10 см, что не наблюдается в режиме с “узким” нагревом.

На рис. 4 приведены результаты динамического моделирования для режимов “широкого” и “узкого” ЭЦР-нагрева с вводимой мощностью 400 кВт. Построены радиальные профили температуры электронов в последовательные моменты времени (верхние графики) и зависимость температуры от времени в нескольких выделенных точках по радиусу (нижние графики). Эволюционный подход позволяет воспроизвести процесс установления самосогласованного профиля, но при этом требует большого времени для расчета. Расчеты выполнены для $\Delta {{T}_{{\max }}} = 50$ эВ и $\Delta {{t}_{{\max }}} = 0.05$ мс на эквидистантной пространственной сетке, состоящей из 200 точек по $\rho $ в диапазоне от 0 до 14 см. Для выхода счета на стационарное решение требовалось 20 вызовов кода QOOT для режима “широкого” нагрева и 80 вызовов кода QOOT для режима “узкого” нагрева, при этом физическое время расчетов на одном ядре рабочей станции (Intel Core i7-4770k/3.5GHz) занимало 15 и 60 суток, соответственно. Для контроля проводились независимые расчеты стационарного самосогласованного профиля температуры, которые требовали от 10 до 30 итераций для точности $\Delta {{\bar {T}}_{{\max }}} = 1$ эВ.

Для режима “широкого” нагрева плазмы начальный профиль электронной температуры эволюционирует достаточно предсказуемо: за время около 0.3 мс устанавливается новый более широкий профиль электронной температуры. Расчеты показали, что конечный профиль электронной температуры практически не зависит от формы начального профиля (в пределах разброса экспериментальных данных порядка 50–100 эВ). Это связано с тем, что, как показывает наше численное моделирование, асимптотически устойчивое стационарное решение уравнения баланса энергии (6) отвечает условию ${{T}_{e}} \gg T_{e}^{*}$ в центральной приосевой области ловушки, поэтому член с $T_{e}^{*}$ можно не учитывать. Поскольку $T_{e}^{*}$ в нашей модели совпадает с начальным распределением электронной температуры, приходим к выводу, что самосогласованный стационарный профиль температуры определяется только параметрами вводимого СВЧ-излучения и пространственным распределением плотности плазмы и магнитного поля.

Результаты моделирования для “узкого” нагрева оказались неожиданными. Несмотря на близкие к предыдущему случаю распределения поглощенной СВЧ-мощности на первом шаге итерации и близкие предсказываемые установившиеся значения температуры в центре ловушки, эволюция профиля электронной температуры существенно отличается. Процесс установления разбивается на два этапа. Вначале узкая приосевая часть плазменного шнура перегревается до температур, превышающих 1 кэВ за счет специфической неустойчивости: повышение электронной температуры, ведет к увеличению поглощения, а увеличение поглощения, в свою очередь, приводит к повышению температуры. При этом на оси возникает конечная производная $d{{T}_{e}}{\text{/}}d\rho $. Это является нефизической особенностью модели, возникающей из-за отсутствия радиального транспорта, к которой (к особенности) электродинамический расчет не чувствителен. В центральной области плазменного шнура происходит формирование плато с характерной шириной около 6 см. Затем, по мере прогрева более отдаленных от оси областей плазменного шнура, доля мощности, приходящаяся на центральную часть, падает. В результате приосевая плазма начинает остывать с формированием временно немонотонного радиального профиля температуры. По мере выхода на стационарное решение плато сдвигается к периферии плазмы, формируя достаточно сложный профиль температуры с резким перепадом.

На рис. 2 видно, что в магнитной конфигурации с “узким” нагревом ЭЦР-область расположена дальше от порта ввода СВЧ-мощности, чем в конфигурации с “широким” нагревом. По всей видимости, триггером перехода в режим с перегревной неустойчивостью, определяющий “узкий” нагрев, является удаление области циклотронного резонанса в сторону от СВЧ-пучка. Это приводит к тому, что зависимость коэффициента ЭЦР-поглощения от температуры электронов в начале эволюции профиля оказывается более резкой, поэтому более холодная периферия плазмы прогревается слабее. Когда разница между энерговыделением в центральной и периферийной областях плазмы превышает некоторое пороговое значение, развивается режим с “узким” нагревом. В процессе самосогласованной эволюции периферия плазмы медленно прогревается и начинает поглощать все больше СВЧ-мощности. Это приводит к недостатку СВЧ-мощности в центральной области плазмы и ее остыванию. Это особенно хорошо видно на зависимостях температуры от времени в центре и на периферии плазмы, приведенных на нижних графиках на рис. 4.

При сопоставлении результатов моделирования с данными эксперимента на установке ГДЛ, приведенными на рис. 1, следует учитывать, что выход на стационарный профиль электронной температуры был возможен только для режима с “широким” нагревом. Режим с “узким” нагревом в экспериментах, отраженных на рис. 1, всегда заканчивался срывом до достижения стационарного профиля электронной температуры [6] (частично задача стабилизации таких разрядов была решена в [8]). Причиной срыва является высокий перепад электронной температуры между периферией и центральной частью плазмы, который, согласно существующим представлениям, должен разрушать вихревое удержание, стабилизирующее МГД-неустойчивости в обсуждаемых экспериментах. Наблюдаемое при этом охлаждение центральной плазмы вовремя ЭЦР-нагрева (см. рис. 1 в [5], рис. 6 в [6]) ранее всегда интерпретировалось как результат ухудшения удержания плазмы перед срывом. Проведенное в настоящей работе моделирование показывает, что немонотонный ход электронной температуры в центре может происходить и по сугубо электродинамическим причинам, не зависящим от условий срыва плазмы, − из-за увеличения поглощения СВЧ-мощности на периферии плазмы. Понимание физики этих процессов очень важно с практической точки зрения, поскольку реализация рассмотренных нами динамических режимов позволяет достичь в экспериментах рекордных пиковых температур электронов на оси плазмы.

В целом, представленная здесь модель склонна завышать значение электронной температуры в центре плазменного шнура и занижать его на периферии. Особенно сильно это проявляется для режима “узкого” нагрева, в которых сразу после резкого перепада температуры формируется область пониженной по сравнению с начальным уровнем температуры (см. верхний правый график на рис. 4). Это противоречит имеющимся экспериментальным данным, но это настоящий эффект в рамках нашей электродинамической модели. Для диссипативной среды уравнения (7)–(10) гарантируют знакоопределенность интегрального поглощения $\int {{{P}_{{ECR}}}\,\rho d\rho } > 0$, но на отдельных магнитных поверхностях, примыкающих к области ЭЦР, может выполняться ${{P}_{{ECR}}}(\rho ) < 0$ за счет пространственной дисперсии поглощения. Поэтому в процессе ЭЦР-нагрева на отдельных магнитных поверхностях возможно охлаждение. Как и токамаках, в открытых ловушках этот тонкий эффект, скорее всего, подавляется радиальным транспортом энергии электронов.

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

В результате объединения транспортного и квазиоптического подходов впервые продемонстрирована возможность самосогласованного моделирования динамических процессов установления радиальных профилей электронной температуры при электронно-циклотронном резонансном нагреве высокотемпературной плазмы в открытой магнитной ловушке с необходимой для интерпретации современного эксперимента точностью. Простейший вариант транспортной модели, учитывающий только продольные потери энергии в газодинамической режиме удержания, изначально рассматривался нами как тест развиваемой численной квазиоптической модели на устойчивость и отсутствие нефизических артефактов, которые мы наблюдали при попытке реализации аналогичных расчетов с помощью многолучевой геометрической оптики. Однако на примере задачи моделирования бифуркации режимов “узкого” и “широкого” СВЧ-нагрева плазмы в установке ГДЛ мы убедились, что рассмотренная простейшая модель на качественном уровне воспроизводит довольно тонкую бифуркацию режимов “узкого” и “широкого” ЭЦР-нагрева на установке ГДЛ. Этот эффект, экспериментально открытый в 2014 г., не воспроизводился стандартными методами моделирования.

В настоящей работе мы не ставили своей задачей продемонстрировать количественную верификацию кода QOOT–SC на данных ГДЛ. Такая верификация требует проведения отдельной экспериментальной кампании с использованием новой многоканальной системы томсоновского рассеяния для измерения профилей электронной температуры и концентрации плазмы [33]. Результаты предварительного моделирования говорят о целесообразности включения радиального транспорта энергии электронов в развиваемую модель. Несмотря на незначительный вклад радиального транспорта по сравнению с продольными потерями энергии, радиальный перенос может устранить большинство недостатков приведенного выше моделирования, такие как острые профили температуры на оси, формирующиеся при быстром СВЧ-нагреве, или охлаждение плазмы на периферии. При необходимости, возможности разработанного гибридного кода могут быть расширены за счет включения эволюции плотности плазмы в транспортную модель. Электродинамическая часть при этом останется неизменной при условии, что основная часть СВЧ‑мощности вкладывается в тепловые электроны. Обобщение квазиоптического подхода на случай плотной плазмы с немаксвелловской функцией распределения является открытой, но не очень актуальной проблемой для современных открытых ловушек с плотной плазмой.

Обнаруженное в модельных расчетах явление формирования в открытой ловушке сложного профиля электронной температуры с резким перепадом, разделяющим центральную и периферийную область плазмы, с сугубо электродинамической природой (поперечный транспорт энергии и частиц в нашей модели отсутствует) заслуживает более подробного исследования. На наш взгляд, это наряду с МГД-турбулентностью еще один возможный способ самоорганизации плазмы. Теорию этого нового явления мы планируем опубликовать отдельно [34].

Авторы выражают благодарность коллегам из ИЯФ СО РАН за доступ к экспериментальным данным ГДЛ и полезное обсуждение работы на семинаре института. Работа выполнена в рамках государственного задания ИПФ РАН (тема № 0030-2021-0002) и ФТИ РАН (по договору между ФТИ РАН и ИПФ РАН № ОК44-2-21 от 06.08.2021).