Физика плазмы, 2020, T. 46, № 10, стр. 890-903

1. ВВЕДЕНИЕ

В середине 1950-х гг. Т.И. Филиппова и Н.В. Филиппов и несколько позже Д. Мейзер независимо предложили, исследовали и опубликовали модификации двух цилиндрических коаксиальных газоразрядных камер с металлическими стенками [1, 2]. В таких камерах в начальной стадии высоковольтного разряда в разреженном дейтерии формируется вокруг изолятора токовая плазменная оболочка (ТПО), которая затем ускоряется магнитным полем разрядного тока. При схождении оболочки к оси камеры происходит частичное сгребание и сжатие нейтрального газа. В итоге при схлопывании ТПО в приосевой области формируется плотный плазменный фокус (ПФ) с высокой плотностью энергии. В начале 1960-x гг. во ВНИИЭФ по предложению Н.Г. Макеева начались исследования новой оригинальной конструкции сферических газоразрядных камер с ПФ. Камеры получили название СФК [3, 4]. Анод и катод камеры имеют сферическую форму и образовывают компактный протяженный газоразрядный промежуток, по которому движется ТПО, плавно изменяя направление движения на 180° в течение разряда.

Разряды типа плазменного фокуса являются частным случаем Z-пинча [5–8], отличаясь от него формой плазменной камеры, которая дает возможность разделить по характерным временам стадию движения ТПО по объему камеры и стадию пинчевания, и, соответственно, использовать для получения пинча более медленные источники тока. Самые первые эксперименты 1950‑х годов проводились с Z-пинчами и исследователи надеялись в них получить нагрев плазмы и зажечь термоядерные реакции, предполагая, что в каждом элементе объема устанавливается максвелловское распределение ионов, а столкновения ионов между собой приводят к реакциям термоядерного синтеза и генерации нейтронов. Этот механизм генерации нейтронов получил название термоядерного. Однако уже в самых первых экспериментах с Z-пинчами было выяснено, что получаемые в них нейтроны не являются термоядерными, а генерируются за счет столкновений ионов, ускоренных до энергий, намного превышающих тепловые, с ионами плазмы пинча. Этот механизм генерации нейтронов получил название ускорительного или пучково-мишенного (beam-target).

Следует сказать, что несмотря на долгую историю Z-пинчевых и плазмо-фокусных исследований, полной ясности и окончательной модели генерации нейтронов, которая бы объясняла все экспериментальные данные и позволяла бы уверенно предсказывать результаты экспериментов, до сих пор не существует. Важным этапом этих исследований, было исследование развития сосисочной неустойчивости. Ссылки на посвященные этому вопросу работы можно найти в обзорах [5, 6], а далее этот вопрос изучался в [9–12]. Активная дискуссия, посвященная вопросу механизмов генерации нейтронов, была в конце 1980‑х–начале 1990-х гг. представлена в обзорных статьях в журнале “Физика плазмы” [13–15]. Важную роль в вопросе исследования механизма нейтронной генерации, сыграло получение высокого нейтронного выхода в экспериментах с 100 нс генераторами тока: установка Ангара-5-1 [16] и установка S-300 [17]. Влияние размеров и времен плазмо-фокусного пинча [18], так же, как и уровень тока [19] влияют на генерацию нейтронов и выяснение этой взаимосвязи также может помочь исследованиям механизма генерации нейтронов. Нельзя не упомянуть об экспериментальных исследованиях на установке PF-1000 (см., например, [20]), в которых также исследовался механизм нейтронной генерации, и, в частности, рассматривалось влияние аномального сопротивления на разрыв тока в пинче. Важная информация о механизме нейтронной генерации может содержаться в экспериментальных данных по анизотропии нейтронного выхода [21]. В последние годы продолжается активное изучение механизмов нейтронной генерации в Z-пинчах и плазменном фокусе. Предлагаются гипотезы о спектре нейтронов из плазмы Z-пинча [22]. В работах [23–25], где в Z-пинче зарегистрированы ионы, ускоренные до десятков МэВ, обсуждаются возможные механизмы генерации столь энергичных ионов и их влияние на нейтронный выход.

В настоящее время во ВНИИЭФ проводятся исследования сильноточных разрядов типа ПФ в камерах СФК. Они нацелены на увеличение выхода нейтронного излучения из области ПФ. Эти исследования имеют большое научное и прикладное значение. Они расширяют возможности применения источников нейтронного излучения на основе ПФ для ядерно-физических, нейтронографических и других прикладных задач.

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

С целью расчетного исследования динамики ТПО в камерах СФК и изучения влияния геометрии камеры и параметров электрического контура на величину нейтронного выхода был разработан и реализован двумерный (в эйлеровых цилиндрических координатах) магнитогидродинамический код, в который был добавлен расчет нейтронного выхода с учетом термоядерного и ускорительного механизма генерации нейтронов в соответствии с общепринятой терминологией. Конкретно для наших расчетов под этими терминами мы будем иметь в виду следующее. Поскольку в наших МГД-расчетах мы предполагаем распределение ионов в каждом элементе объема равновесным, что, как можно показать, выполняется для основной массы плазмы, мы будем использовать термин “термоядерный механизм” при расчетах генерации нейтронов с использованием максвелловского распределения ионов. Кроме того, мы на фоне МГД-расчетов используем модель ускорения ионов в электрическом поле, которая объясняет многие экспериментальные факты, связанные с генерацией нейтронов, называя этот механизм генерации нейтронов ускорительным.

Двумерное МГД-моделирование широко применяется для описания различных процессов, происходящих в камерах ПФ [5, 6, 26], однако, в рамках МГД-приближения и в предположении локального термоядерного механизма генерации нейтронов не удается описать наблюдаемый нейтронный выход, поскольку нейтронный выход в камерах ПФ в основном определяется ускорительным механизмом. Включение кинетического описания движения заряженных частиц (например, в приближении частиц в ячейке) в области фокусировки токово-плазменной оболочки и сопряжение этого описания с МГД-подходом для всей остальной области ПФ (например, в [27, 28]) кажется перспективным подходом, но является очень трудоемким и затратным с вычислительной точки зрения, и пока еще не доведен до того уровня, когда им можно пользоваться для описания генерации нейтронов. В нашей работе мы продолжаем подход, предложенный в [29], в котором на фоне МГД-расчета учитывается движение ионных пучков, ускоряемых напряжением, вызванным в области фокусировки аномальным сопротивлением. Столкновение ускоренных ионов с плотной плазмой в области фокусировки приводит к генерации нейтронов и в [29] показано, что такой подход хорошо описывает экспериментальные результаты по генерации нейтронов для ряда геометрий плазменного фокуса.

2. ФИЗИЧЕСКАЯ МОДЕЛЬ. ОСНОВНЫЕ УРАВНЕНИЯ И ФОРМУЛЫ

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

где ρ – массовая плотность; $\vec {v}({{\nu }_{r}},{{{v}}_{z}})$ – ионная скорость и ее компоненты; $M = \rho {{v}_{r}}$ и $N = \rho {{v}_{z}}$ – r и z-компоненты плотности импульса; P = P_e + P_i – суммарное давление, P_e – давление электронов, P_i – давление ионов; H – магнитное поле, которое в наших задачах имеет только φ-ю компоненту; ε_e, ε_i – плотности внутренней электронной и ионной энергии; ${{\varepsilon }_{{kin}}} = \rho (\nu _{r}^{2} + \nu _{z}^{2}){\text{/}}2$ – плотность кинетической энергии; ε_sum = ε_i + ε_e + ε_kin – сумма плотностей ионной, электронной и кинетической энергии; E_r, E_z – r и z-компоненты напряженности электрического поля; ${{\vec {Q}}_{e}}$, ${{\vec {Q}}_{i}}$, – электронная и ионная плотности потоков тепла; ${{Q}_{{joule}}} = \sigma (E_{r}^{2} + E_{z}^{2})$ – джоулев нагрев плазмы; Q_ei – поток энергии от ионов к электронам; S_rad – потери на излучение; c – скорость света в вакууме; σ – проводимость плазмы; T_e, T_i – температуры электронов и ионов; ${{j}_{z}} = \frac{c}{{4\pi }}\frac{1}{r}\frac{\partial }{{\partial r}}(rH)$, ${{j}_{r}} = - \frac{c}{{4\pi }}\frac{{\partial H}}{{\partial z}}$ – компоненты плотности тока.

При вычислении P_i(ρ, T_i) и P_e(ρ, T_e) использовалось уравнение состояния идеального одноатомного газа, для электронной компоненты плазмы использовалось уравнение Саха для определения связи между температурой T_e, внутренней энергией и степенью ионизации α. P_i(ρ, T_i) и ε_i включали вклад нейтральной компоненты.

Для расчета проводимости σ использовались классические формулы с учетом рассеяния электронов на нейтралах

где L_e – электронный кулоновский логарифм, $k = \frac{{1 + 3.92 \cdot \left( {\omega \tau } \right)_{e}^{2} + 0.265 \cdot \left( {\omega \tau } \right)_{e}^{4}}}{{1 + 4.33 \cdot \left( {\omega \tau } \right)_{e}^{2} + 0.524 \cdot \left( {\omega \tau } \right)_{e}^{4}}}$ – коэффициент замагничивания проводимости, (ωτ)_e – электронная замагниченность [30]. Степень ионизации α рассчитывалась в результате решения уравнения Саха [31] и уравнения связи между удельной внутренней энергией и температурой электронов с учетом вклада энергии диссоциации где I = 13.6 эВ, R = 9.65, $\alpha {\kern 1pt} * = \left\{ \begin{gathered} \frac{{{{\varphi }_{0}}}}{{{{\alpha }_{0}}}}\alpha ,\quad \alpha < {{\alpha }_{0}} \hfill \\ {{\varphi }_{0}},\quad \alpha \geqslant {{\alpha }_{0}} \hfill \\ \end{gathered} \right.$, φ₀ = = 2.28 эВ, α₀ = 0.02. В формулах (2) σ имеет размерность с^–1·10⁷, температура электронов T_e выражается в кэВ. Кроме того, в случае если (ωτ)_e > 1 и (ωτ)_i > 1 ((ωτ)_i – замагниченность ионов), проводимость вычисляется с учетом поправок на появление аномального сопротивления плазмы вследствие развития нижнегибридной дрейфовой неустойчивости [32]: где β – отношение токовой скорости к тепловой скорости ионов, магнитное поле H измеряется в Э × 10⁷, A – атомный вес иона (г/моль), скорость света c в см/с × 10⁷.

Нейтронный выход вычислялся как сумма вкладов термоядерного и ускорительного механизмов генерации нейтронов:

где $\left\langle {\sigma \nu } \right\rangle $ – интенсивность ДД или ДТ реакции; n – плотность ядер плазмы, на которых идет реакция; α_i – доля тока, переносимого ионами; j_z – осевая составляющая плотности тока; σ(ε_i) – сечение ДД или ДТ реакции, зависящее от энергии ускоренного иона ε_i.

При этом для оценок нейтронного выхода по ускорительному механизму предполагалось, что вблизи оси, на расстояниях от нее меньше трех ларморовских радиусов ускоренных ионов, часть ионов имеет возможность ускоряться в электрическом поле вдоль оси, одновременно испытывая силу торможения F. Величина этой силы F в зависимости от энергии иона, электронной и ионной температур плазмы вычислялась по формулам, представленным в [29, 33]. Считалось, что ускоренные ионы замещали в приосевой области α_i ~ 15% максимальной (до текущего значения z) плотности тока j_z [29]. Расчет энергии пучка ионов проводился по строкам расчетных ячеек в направлении от анода к катоду и начинался в той ячейке, в которой сила, действующая на ион со стороны электрического поля, впервые превышала минимальное значение силы торможения ионов F со стороны плазмы. Для безразмерной функции торможения значение минимума достигается при энергии дейтрона ${{\varepsilon }_{0}} = {{\left( {\frac{{9\pi }}{{\mu {{\lambda }^{2}}}}} \right)}^{{\frac{1}{3}}}}{{T}_{e}}$, где μ – отношение массы электрона к массе иона, λ – отношение электронного кулоновского логарифма к ионному кулоновскому логарифму, T_e – температура электронов [10]. Мы считали, что в каждой строке, в той ячейке, в которой начинается расчет энергии пучка ионов, начальное значение энергии равно ε₀. В дальнейшем ионы могут увеличивать свою энергии, ускоряясь в электрических полях E_z, возникающих при вытекании магнитного потока на ось, т.е.

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

Следует заметить, что мы пренебрегали вкладом ускоренных ионов в баланс импульса (закон Ома) и энергии (джоулево тепловыделение), хотя их вклад в плотность тока в приосевой области мог достигать α ~ 15%, не такой уж малой величины. В оправдание можно сказать, что полный баланс энергии и импульса из-за появления ионных пучков в системе не меняется, поскольку уменьшение джоулевого тепловыделения из-за затрат на разгон ионов в области больших электрических полей компенсируется тепловыделением при торможении этих ионов в плотной и относительно холодной плазме по пути их движения от анода к катоду. По существу, часть джоулева тепла выделяется и распределяется по более широкой области плотной плазмы. Если говорить о полной энергии пучка E_b, которая перераспределяется таким образом, то по оценкам она составляет ~0.1% тепловой энергии плазмы. Поскольку скорости ускоренных ионов ${{{v}}_{i}}$ велики по сравнению с тепловыми скоростями плазмы, то характерный импульс ${{P}_{b}}\sim {{E}_{b}}{\text{/}}{{{v}}_{i}}$, приобретаемый ими, относительно совсем невелик по сравнению с импульсом, характерным для динамики плазмы в это время, и составляет примерно 3 × 10^–5 последнего. Представляется, что перераспределение по пространству очень малых величин порядка энергии и импульса пучка не должно сколь-нибудь заметно сказываться на экспериментальных результатах, получаемых в плазменном фокусе.

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

Металлические границы в камере считаются жесткими, холодными, идеально проводящими стенками. На изоляторе задается значение граничного магнитного поля, которое определяется по значению тока в цепи разряда. Ток в цепи определяется из решения электротехнического уравнения:

где L – индуктивность контура, $U(t) = {{U}_{0}} - $ $ - \;\frac{1}{C}\int_0^t {I(t)dt} $, U₀ – зарядное напряжение батареи, C – емкость батареи, R – сопротивление контура, U_PF – напряжение на входе в камеру, которое определяется по результатам МГД-расчета.

3. СХЕМА РАСЧЕТА. ОСОБЕННОСТИ ЧИСЛЕННОЙ МОДЕЛИ И ВЛИЯНИЕ “ВАКУУМНОЙ ПЛОТНОСТИ”

МГД-расчет состоит из следующих этапов:

• газодинамический этап счета, на котором определяются стандартные газодинамические величины: ионная скорость, плотность, электронная и ионная энергии;

• расчет проводимости, распределение магнитного и электрического полей;

• расчет обмена энергией между электронной и ионной компонентой, а также расчет электронной и ионной теплопроводностей, потерь энергии на излучение и джоулева нагрева плазмы;

• расчет степени ионизации, электронной и ионной температуры и давления;

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

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

При расчете по явным схемам есть ограничение Куранта на величину счетного шага по времени:

где Δl – линейный размер ячейки, q = 0.25–0.75, ${{c}_{{ms}}} = \sqrt {c_{s}^{2} + c_{A}^{2}} $, c_s – скорость звука в веществе, ${{c}_{A}} = \frac{H}{{\sqrt {4\pi \rho } }}$ – альфвеновская скорость. При достаточно низких плотностях газа в разреженной области за оболочкой (ρ < ρ₀/10³) величина c_A становится велика, и временной шаг уменьшается значительно. От данного ограничения можно избавиться, если совместить расчет компонент ионной скорости с расчетом электрических полей и использовать неявную схему для определения электрических полей E_r и E_z как это сделано в работе [29]: где $\kappa = \frac{{{{c}^{2}}}}{{4\pi \sigma }}$ – коэффициент диффузии магнитного поля.

Используя разностные выражения для E_r, E_z и ∂H/∂t (3), можно получить трехточечную диффузионную разностную систему для E_r и E_z, но для коэффициента диффузии $\kappa + c_{A}^{2}\Delta t$, которая решается методом последовательных прогонок по столбцам и строкам с итерациями. Сходимость итераций определяется по значению H. Данный подход позволяет значительно снизить влияние плотности газа в разреженной области на величину счетного шага и проводить расчет со значением плотности до 6-и порядков меньшим, чем ее начальное значение.

В проведенных расчетах было обнаружено, что за ТПО в “вакуумной” области происходит сильный разогрев плазмы (по-видимому, из-за нефизических счетных эффектов), настолько сильный, что это приводит к тепловому давлению, сравнимому с магнитным, и, в результате, к частичной экранировке области пинча от тока, протекающего по этой малоплотной плазме. В результате пинч не добирает тока и в нем не создаются условия для генерации нейтронов по ускорительному механизму. Повышение напряжения на оси камеры и соответственно величины нейтронного выхода в расчете было достигнуто путем ограничения температуры ионов и электронов до 1 кэВ в низкоплотной области: там, где плотность в 10 и более раз меньше начальной плотности невозмущенного газа. Меньшая температура в вакуумной области способствовала большему сжатию плазмы в момент фокусировки, т.к. область с меньшей температурой (и давлением) слабее экранировала магнитное поле. После введения ограничения на температуру плазмы в низкоплотной области удалось значительно улучшить согласование экспериментального и расчетного выхода нейтронов.

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

Расчет с учетом аномального сопротивления позволяет получить значительный прирост напряжения и, соответственно, электрического поля на оси камеры, однако из-за резкого включения аномального сопротивления в счете коэффициент диффузии для магнитного поля резко возрастал, что приводило к численной неустойчивости. Для устранения численной неустойчивости были введены ограничение на уменьшение проводимости плазмы за один счетный шаг на уровне 5%, а также разработана модель “плавного включения” аномального сопротивления.

Один из вопросов моделирования течений плазмы в камере для плазменного фокуса в эйлеровых координатах состоит в том, каким образом проводить расчет в областях, в которых уже не осталось вещества (вакуумных областях). Такие области возникают за плазменной оболочкой в процессе ее движения к оси камеры. Для того чтобы преодолеть эту трудность, приходится полагать, что плотность не может быть меньше какого-то минимального значения ρ_min = ρ₀/10^m. При малых m = 2–3 плотность газа за плазменной оболочкой достаточно велика, в результате происходит шунтирование тока в области позади оболочки, и процесс фокусировки плазмы на оси получается менее выраженным. При m > 6 начинает уменьшаться счетный шаг по времени из-за требования устойчивости расчета.

Для изучения влияния величины плотности остаточного газа на процесс фокусировки плазмы была проведена серия расчетов с разными значениями m в диапазоне от 2 до 6.

Расчеты проводились для сферической плазменной камеры с диаметром анода 60 мм и диаметром катода 120 мм. Геометрия камеры представлена на рис. 1.

Расчетная область представляет собой прямоугольник. Расчетная сетка имеет размер ячеек 0.1 мм. Внутренняя область камеры заполнена ДТ-смесью. Плотность газа ρ₀ = 4.8 × 10^–6 г/см³ (18 мм рт. ст.). Начальная температура ионов и электронов T_e = T_i = 11600 К (1 эВ). В ячейках на границе с изолятором задается граничное значение магнитного поля H₀. Для расчета H₀ находилась величина тока в цепи из уравнения (2). Параметры электрического контура для уравнения (2) принимались следующими: C = 69.5 мкФ, L = = 27 нГн, R = 3 мОм, U₀ = 21 кВ. При расчете граничных условий нормальная компонента скорости на границе с электродами считалась равной нулю. Теплопроводность и радиационные потери энергии, а также диффузия магнитного поля в плазме в расчете не учитывались (проводимость плазмы считалась бесконечной). Полученные в МГД-расчетах зависимости напряжения на входе в камеру (U_вход на рис. 1) от времени для различных значений минимальной плотности газа за плазменной оболочкой представлены на рис. 2.

Исходя из закономерностей, полученных в расчетах, можно сделать вывод, что при больших значениях остаточной плотности действительно происходит шунтирование тока в вакуумной области и экранировка магнитного поля, в результате фокусировка плазмы слабо выражена и пик напряжения практически отсутствует, что не согласуется с экспериментальными данными. При малых значениях остаточной плотности (m ≥ 5) отличие результатов расчета для различных значений m становится незначительным и существенно улучшается согласие с экспериментом, поэтому во всех последующих расчетах плотность остаточного газа принималось равной ρ_min = ρ₀/10⁵.

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

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

4. РЕЗУЛЬТАТЫ РАСЧЕТОВ И СРАВНЕНИЕ С ЭКСПЕРИМЕНТОМ

На рис. 3 приведены напряжения на входе в камеру (U_вход) и на оси камеры (U_ось), полученные в МГД-расчете с учетом аномального сопротивления плазмы. Расчет выполнен для сферической камеры с диаметром катода 340 мм и диаметром анода 180 мм, используемой в составе установки КПУ-200 [36]. Напряжение зарядки конденсаторной батареи в расчете составляло 30 кВ. Пик напряжения на входе в камеру формируется в момент схлопывания плазменной оболочки и образования ПФ. Пик напряжения на оси соответствует моменту развала пинча и выхода магнитного потока на ось камеры. Длительность пика напряжения на полувысоте на входе в камеру составляет примерно 130 нс, напряжение на оси камеры имеет 2 пика с длительностями на полувысоте 20 и 15 нс. Амплитуда напряжения на оси достигает 600 кВ, что примерно в 4 раза превышает амплитуду напряжения на входе в камеру (150 кВ). При этом напряжение на входе в 5 раз превышает напряжение на конденсаторной батарее. Такое значительное повышение напряжения может привести к генерации энергичных ионов. В момент обострения напряжения в ПФ происходит ускорение ионов вдоль оси камеры и генерация ускорительных нейтронов. Максимум напряжения на оси соответствует максимальной интенсивности генерации ускорительных нейтронов. Надо сказать, что в недавних работах [24, 25, 37–39] на разных установках зарегистрированы ионы, ускоренные до нескольких единиц и даже до десятков МэВ. Эти энергии превышают в десятки (и даже сотню) раз те энергии, которые могут обеспечивать начальные напряжения на генераторах тока, используемых в экспериментах. Возможно, полученные в наших расчетах напряженности электрического поля, которые с учетом их зависимости от координат и времени, могут разгонять ионы до больших энергий, помогут понять механизм генерации столь энергичных ионов.

Для проверки корректности разработанной модели было проведено сравнение результатов расчета с экспериментальными данными по току, производной тока, напряжению на камере и величине нейтронного выхода. На рис. 4 приведено сравнение экспериментальных и расчетных кривых тока, напряжения и производной тока для установок ПФ с камерой диаметром 120 мм (рис. 4а, 4б) и камерой диаметром 340 мм (рис. 4в). Здесь следует отметить, что для согласования с экспериментом времени возникновения особенности и, соответственно, времени фокусировки плазменной оболочки приходилось полагать, что в камере диаметром 340 мм начальная плотность в 1.75 раза меньше реальной. Причина отличия расчетного и экспериментального времени фокусировки, по-видимому, связана с тем, что в эксперименте вследствие несимметричности ТПО по азимутальному направлению возникает волокнистая структура оболочки, и часть газа просачивается между волокнами и не вовлекается в процесс сгребания плазмы к оси. Полностью учесть такие эффекты в двумерном расчете не представляется возможным, поэтому расчет проводился с уменьшенным значением начального давления газа [40, 41]. Можно упомянуть также некоторые эффекты, которые потенциально могли бы приводить к более быстрому приходу ТПО к оси. Это, например, отрыв плазмы от анода, обусловленный приносом к аноду магнитного потока благодаря эффекту Холла [42–44 см. также 12]. Однако наши попытки учесть этот эффект в расчетах с помощью значительного уменьшения начальной плотности газа в прианодной области на масштабах ~c/ω_pi (ω_pi – плазменная ионная частота) [12, 43], также, впрочем, как попытки учета этого эффекта в камере МАГО [12] (МАГО/MTF, МАГнитное Обжатие/Magnetized Target Fusion, – альтернативный подход к решению проблемы управляемого термоядерного синтеза, активно исследовался в 1980-х–2000-х гг.), включающие даже и прямой учет эффекта Холла, не привели к сколь-нибудь ощутимому уменьшению времени прихода ТПО к оси. Причинами тому, по-видимому, является малость масштаба ~c/ω_pi по сравнению с характерными масштабами задачи, двумерное гидродинамическое перемешивание этой зоны с бо́льшими масштабами и относительно низкие температуры, и, следовательно, высокая столкновительность плазмы. Следует упомянуть также смытие приэлектродного материала с анода, которое играет важную роль в камере МАГО [12], но не исследовано для плазменного фокуса, и этот эффект также может в экспериментах подавлять отрыв плазмы от анода. Возможным эффектом, приводящим к сбросу части массы плазмы при ускорении ТПО и его движении к оси, могла бы быть рэлей-тейлоровская (РТ) неустойчивость, развитие которой на коротких длинах волн может приводить к сбросу части массы плазмы, и, соответственно к дополнительному ускорению остальной массы, также как это происходит при ускорении лайнеров [45]. Однако для осесиммметричных возмущений этот эффект вполне моделируется [45] двумерными МГД-расчетами. Что же касается неосесимметричных возмущений, развивающихся при РТ-неустойчивости, вполне возможно, что их развитие приводит к тому, что часть массы плазмы остается невовлеченной в движение и вынуждает уменьшать массу плазмы в двумерных расчетах. Правда остается вопрос, связанный с тем, что развитие РТ-неустойчивости должно быть подобно для малых и больших камер, а для описания экспериментов нам приходится уменьшать массу плазмы только для большой камеры.

Для параметров работы установки ПФ, приведенных на рис. 4а, экспериментальный нейтронный выход составил (3.0 ± 0.6) × 10⁹ ДД-нейтронов. В расчете выход составил 2.8 × 10⁹ нейтронов. Для камеры диаметром 340 мм при зарядке конденсаторной батареи до 30 кВ нейтронный выход составил (6.2 ± 0.9) × 10¹¹ ДД-нейтронов. В расчете нейтронный выход составил 2.0 × 10¹¹ нейтронов, что в 3.1 раза меньше экспериментального.

Систематическое занижение результатов расчета нейтронного выхода для камеры с внешним диаметром 340 мм имело место для всего диапазона энергии, для которого проводилось моделирование и имелись надежные экспериментальные данные по нейтронному выходу. Отличие расчетного выхода от экспериментального частично можно объяснить наличием анизотропии нейтронного излучения для ДД-реакции. При измерении нейтронного выхода детекторы устанавливались близко к оси разрядной камеры, при этом при определении интегрального выхода поправка на анизотропию не вводилась. Но согласно данным, представленным в [46], отличие дифференциального сечения выхода нейтронов вперед от среднего по телесному углу может находиться в диапазоне от 1.9 до 3.5 для энергий налетающих дейтронов 150–1000 КэВ. Если предположить, что в наших экспериментах энергия налетающих дейтронов была ~150 КэВ (см. ниже обсуждение величины характерных энергий налетающих дейтронов), то можно ожидать, что отношение экспериментальных и расчетных нейтронных выходов при учете анизотропии нейтронного выхода уменьшится с 3.1 до 1.6.

На рис. 5 представлены распределения плотности (рис. 5а) и температуры ионов (рис. 5б) и электронов плазмы (рис. 5в), полученные в МГД-расчетах для камеры диаметром 340 мм. На рис. 5а отчетливо прослеживается, что при движении ТПО часть массы газа прижимается к катоду. Можно обратить внимание на то, что граница ТПО на рис. 5а довольно гладкая, хотя априори могла бы быть искажена из-за развития РТ-неустойчивости. Эта загадка является общей для многих течений в условиях нелинейного развития РТ-неустойчивости: известно, что “пузыри” (bubbles) в РТ-неустойчивости имеют гладкую форму, это демонстрируют и эксперименты (см., например, [45, 47]) и расчеты (см., например, [45, 48]). Аналогичная гладкость границы плазмы имеется и в расчетах камеры МАГО [12]. Удовлетворительного количественного объяснения этого феномена, насколько мы знаем, не существует. Можно предположить, что хотя, несомненно, на малых масштабах развитие РТ-неустойчивости должно приводить к неоднородной турбулентной форме границы, на больших масштабах, сравнимых с полным масштабом задачи, двумерные эффекты (растягивание масштабов в пузыре и шир (shear) скорости) подавляют развитие неустойчивости, что и наблюдается в экспериментах и расчетах.

Можно отметить, что форма и положение плазменной оболочки, полученные в МГД-расчетах, хорошо соотносятся с экспериментальными оптическими снимками ТПО, полученными в ряде работ. Представленный на рис. 5 момент времени t = 5.42 мкс соответствует пику напряжения на входе в камеру и минимальному радиусу пинча. В этот момент плотность плазмы в области пинча достигает своего максимального значения. Момент времени t = 5.43 мкс соответствует пику напряжения на оси камеры (рис. 4) и максимальной интенсивности генерации нейтронов. В этот момент тонкая перетяжка пинча разрушается и начинается стадия его распада.

Помимо величины интегрального выхода нейтронов, представляет интерес также распределение генерируемых в ПФ нейтронов по энергии, т.е. энергетический спектр нейтронов. В ряде работ [3, 49, 50] были представлены результаты экспериментальных исследований спектрального состава ДТ-нейтронов плазменного фокуса в направлениях 0°, 90° и 180° относительно оси газоразрядной камеры. В данной работе с использованием заложенной в МГД-расчет модели генерации нейтронов были проведены расчетные оценки спектра ДТ-нейтронов для направлений 0°, 90° и 180° относительно оси камеры. В ходе расчета генерации нейтронов по ускорительному механизму было получено распределение нейтронного выхода по энергии налетающих ионов дейтерия или трития (рис. 6). Спектр, представленный на рис. 6, соответствует МГД-расчету для плазменной камеры диаметром 340 мм при наполнении смесью дейтерия и трития. Напряжение зарядки конденсаторной батареи составляло 30 кВ, начальное давление газа в камере – 12 мм рт.ст. Интегральный выход в расчете составил 9.1 × 10¹² ДТ-нейтронов. Характерная величина энергии ионов, при которой реализуется максимум генерации нейтронов ~70 кэВ, примерно соответствует средней энергии ускоренных ионов 50–70 кэВ, полученной с помощью измерений пространственно-энергетической анизотропии нейтронного излучения [19].

Далее был осуществлен переход от распределения нейтронного выхода по энергии налетающих ионов f(E_i) к распределению нейтронного выхода по энергии нейтрона f(E_n) (т.е. искомого спектра) для фиксированного угла вылета нейтрона 0°, 90° и 180° к оси разрядной камеры. Такой переход был осуществлен путем умножения полученного из МГД-расчета распределения f(E_i) на производную ${{\left( {\frac{{d{{E}_{i}}}}{{d{{E}_{n}}}}} \right)}_{{_{{{{\psi }_{0}}}}}}}$, вычисленную для трех фиксированных значений углов: ψ₀ = 0°, ψ₀ = 90° и ψ₀ = 180°:

При вычислении производной ${{\left( {\frac{{d{{E}_{i}}}}{{d{{E}_{n}}}}} \right)}_{{_{{{{\psi }_{0}}}}}}}$ мы использовали формулу связи между энергией нейтрона и энергией налетающего иона в виде:

где E_n – энергия нейтрона; E_i1 – энергия налетающего иона (в наших расчетных оценках налетает ион дейтерия); m_n – масса нейтрона; m_i1 – масса налетающего иона; m_i2 – масса покоящегося иона (в наших расчетных оценках покоится ион трития); Θ – угол разлета продуктов реакции в системе центра масс; $\beta = \frac{{({{m}_{{i1}}} + {{m}_{{i2}}}){{m}_{{i2}}}{{m}_{\alpha }}}}{{({{m}_{{i3}}} + {{m}_{\alpha }}){{m}_{{i1}}}{{m}_{n}}}}$ + + $\frac{{{{{({{m}_{{i1}}} + {{m}_{{i2}}})}}^{2}}{{m}_{\alpha }}}}{{({{m}_{n}} + {{m}_{\alpha }}){{m}_{{i1}}}{{m}_{\alpha }}}}\frac{Q}{{{{E}_{{i1}}}}}$; m_α – масса ядра ${}_{2}^{4}{\text{Не}}$; Q = = 17.6 МэВ – энергия ДТ-реакции. Для перехода от угла Θ к углу вылета нейтрона в лабораторной системе отсчета ψ использовалось соотношение

При получении расчетных оценок спектра ДТ‑нейтронов предполагалось, что сечение ДТ‑реакции является изотропным в системе центра масс продуктов реакции, независимо от энергии. Согласно экспериментальным данным [51] это предположение справедливо для относительно невысоких энергий налетающих ионов (<200 кэВ). Кроме того, используемые соотношения (8) и (9) получены из рассмотрения кинематики ДТ-реакции в нерелятивистском случае. Для нейтрона с энергией 15 МэВ отличие в скорости при расчете по релятивистской и по нерелятивистской формуле составляет ~1%. Поэтому в приведенных оценках для упрощения использовались нерелятивистские формулы для энергии и импульса.

На рис. 7 приведено сравнение расчетного спектра нейтронов с экспериментальными данными [50]. Для того чтобы сравниться с экспериментом, расчетные спектры были нормированы так, чтобы максимальное значение соответствовало единице.

По результатам расчетов можно отметить хорошее согласие с экспериментом по положению максимума спектральных распределений для всех трех направлений вылета нейтрона. Но ширина расчетного спектра для направления вылета нейтрона 90° получилась значительно уже экспериментального спектра для данного направления. Уширение спектра в направлении 90° в эксперименте можно объяснить наличием радиальной составляющей скорости движения налетающего иона, которое может быть вызвано отклонением движущегося вблизи оси иона в магнитном поле, что не учитывалось в нашей модели. Это отклонение может также привести и к дополнительному уширению спектра нейтронов под углами 0° и 180°. Кроме того к уширению спектра может также привести взаимодействие нейтронов с элементами конструкции разрядной камеры, что также не учитывалось в расчете.

5. СКЕЙЛИНГИ И ВОЗМОЖНОСТИ ПОВЫШЕНИЯ НЕЙТРОННОГО ВЫХОДА

При практическом использовании устройств на основе ПФ актуальной задачей является повышение нейтронного выхода. С использованием разработанной программы была проведена серия расчетов для двух сферических камер ПФ диаметром 120 мм и 340 мм.

Задачей этих расчетов было исследование возможности повышения выхода нейтронов при увеличении энергии, запасаемой в конденсаторной батарее установки. Для установки, работающей с камерой 120 мм, были проведены расчеты в диапазоне энергии от 8 до 29 кДж. Для установки, работающей с камерой 340 мм, расчеты были проведены в диапазоне энергии от 45 до 245 кДж. Для согласования электрического контура установки и плазменной нейтронной камеры для каждой точки по энергии были проведены расчеты с различным начальным давлением дейтерия. Для камеры диаметром 120 мм расчеты проводились в диапазоне от 4 до 30 мм рт. ст. (рис. 8). Для камеры диаметром 340 мм расчеты проводились в диапазоне от 18 до 30 мм рт. ст.

Оптимальной по давлению считалась точка, при которой нейтронный выход для данной энергии максимален. Расчеты были выполнены для наполнения камеры дейтерием, при расчете нейтронного выхода по ускорительному механизму использовались сечения ДД-реакции [52].

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

Исходя из анализа полученной зависимости, можно отметить, что полученная в расчете степень для скейлинга 5.3 несколько больше экспериментально наблюдаемой величины, которая находится обычно в диапазоне от 3 до 5. Наибольший нейтронный выход 3.4 × 10¹⁰ нейтронов достигается при токе ~1.0 МА. При этом энергия, запасенная в конденсаторной батарее установки равна 29.4 кДж. При переходе на ДТ-наполнение можно ожидать выход ~1.7 × 10¹² нейтронов.

На рис. 10 приведены расчетные и экспериментальные токовые скейлинги для установки ПФ, работающей с камерой диаметром 340 мм.

Можно отметить, что расчетный скейлинг лежит ниже экспериментального и имеет меньший показатель степени, т.е. выход растет медленнее по сравнению с экспериментом при увеличении амплитуды разрядного тока. В диапазоне тока от 0.8 до 2.0 МА это отличие изменяется от 2 до 5 раз. Как отмечалось выше, с учетом анизотропии нейтронного излучения для ДД-реакции, которая может находиться в диапазоне от 1.9 до 3.5 для выхода вперед по отношению к среднему выходу в 4π, экспериментальный выход может быть скорректирован и приближен к расчетному. Если продолжить экспериментальный токовый скейлинг (который, с учетом сказанного, должен все же быть скорректирован) в область более высоких токов, то при токе ~2.5 МА достигается выход ~3 × 10¹² ДД-нейтронов или ~1.5 × 10¹⁴ нейтронов при переходе на ДТ-наполнение. Экстраполяция расчетного скейлинга в область токов ~2.5 МА дает меньшую величину нейтронного выхода: ~5 × × 10¹¹ ДД-нейтронов или ~2.5 × 10¹³ ДТ-нейтронов.

При расчетах также был исследован вклад термоядерного механизма генерации нейтронов в полный нейтронный выход. На рис. 11 представлена доля термоядерных нейтронов в зависимости от энергии, рассчитанная для установок, работающих с камерой 120 мм и 340 мм в диапазоне энергий от 8 до 250 кДж.

По результатам расчетов видно, что доля термоядерных нейтронов возрастает с увеличением энергии, запасенной в конденсаторной батарее установке от ~4 × 10^–4 до ~2 × 10^–2. В ходе МГД-расчетов работы установки ПФ с камерой 340 мм при увеличении емкости конденсаторной батареи до 600 мкФ и напряжении зарядки 40 кВ доля термоядерных нейтронов составила ~4 × 10^–2 при энергии 480 кДж.

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

Разработана двумерная двухтемпературная МГД-модель расчета течений плазмы в камерах ПФ с учетом генерации нейтронного излучения. С использованием МГД-модели была разработана методика численного моделирования процессов ускоренного движения плазменной оболочки, а также расчета нейтронного выхода. В расчетах исследовано влияние диффузии магнитного поля, теплопроводности плазмы и аномального сопротивления на динамику плазменной оболочки. Изучено влияние минимального значения плотности остаточного газа в вакуумной области. Показано, что удовлетворительное согласие с экспериментом достигается при плотностях остаточного газа ≤10^–5 от начальной плотности.

Сравнение МГД-расчетов с экспериментом по электротехническим параметрам, таким как производная тока, ток и напряжение на входе в камеру показало адекватность разработанной модели и возможность использования программы для расчетного исследования динамики плазменной оболочки в камерах ПФ. В то же время для больших камер для согласования расчетного времени сжатия с экспериментальным приходится предполагать, что в наших двумерных расчетах значительная часть плазмы не вовлекается в движение. Возможной причиной этого может быть то, что в эксперименте вследствие несимметричности ТПО по азимутальному направлению возникает волокнистая структура оболочки, и часть газа просачивается между волокнами и не вовлекается в процесс сгребания плазмы к оси и, может быть, остается в виде слабоионизованных азимутально-несимметричных “облаков”. Во всяком случае следует признать, что в настоящее время наша двумерная МГД-модель не является полностью согласованной и для больших камер для лучшего согласования времен сжатия требует введения дополнительного эмпирического параметра – доли вовлеченной в движение плазмы.

При расчете нейтронного выхода с моделью ускорительного механизма генерации нейтронов было получено удовлетворительное согласие по величине интегрального выхода нейтронов.

Оценен вклад термоядерных нейтронов в полный нейтронный выход. По результатам проведенных расчетов доля термоядерных нейтронов составляла от ~4 × 10^–4 до ~4 × 10^–2 от полного выхода. Таким образом, при данных энергиях термоядерные нейтроны вносят незначительный вклад в полный выход. Тем не менее, если экстраполировать полученную зависимость в область больших энергий, то при энергиях >1.5 МДж, вклад термоядерных нейтронов в полный нейтронный выход составит >10%.

В ходе проведенных расчетов было получен энергетический спектр нейтронов в направлении 0°, 90° и 180° относительно оси разрядной камеры. Полученные результаты были сопоставлены с экспериментальными результатами по измерению спектра ДТ-нейтронов из плазменного фокуса в указанных направлениях. Получено хорошее согласие с экспериментом по положению максимума спектра для всех трех направлений. Ширина экспериментального спектра для направления 90° значительно превосходит ширину распределения, полученного в расчете, что может быть связано с наличием радиальной составляющей скорости движения налетающего иона и взаимодействием нейтронов с элементами конструкции установки. Эти факторы не учитывались в проведенных расчетах.

Удовлетворительное согласие расчетов с экспериментом по величине интегрального выхода нейтронов позволило провести расчетные оценки по возможности увеличения нейтронного выхода на двух установках ПФ со сферическими камерами с внешним диаметром 120 и 340 мм. Расчеты показали, что при увеличении энергии, запасаемой в конденсаторной батарее установок, нейтронный выход растет в соответствии с известным скейлингом Y ~ I ^3–5. Однако, в случае с камерой 340 мм показатель степени для скейлинга, полученного из МГД-расчетов, отличается от экспериментального (3.19 и 4.16 соответственно). В то же время для установки с малой камерой диаметром 120 мм этот показатель несколько выше (5.29).