Numerical modeling of plasma dynamics and neutron generation in Z-pinch at the Angara-5-1 facility

Толық мәтін

1. Введение

Исследования генерации нейтронов при сжатии Z -пинчей продолжаются уже давно, начиная с самых первых экспериментов на дейтериевых Z -пинчах, проведенных в 1950-е годы [1–3]. Эти исследования переросли в исследования генерации нейтронов на плазменном фокусе [4, 5], который представляет собой нецилиндрический Z -пинч. Однако уже в самых первых экспериментах с Z -пинчами было выяснено, что получаемые в них нейтроны не являются термоядерными, а генерируются за счет столкновений ионов, ускоренных до энергий, намного превышающих тепловые, с ионами плазмы пинча. Этот механизм генерации нейтронов получил название ускорительного или пучково-мишенного (beam-target). Новым шагом в исследованиях генерации нейтронов на Z -пинчах стали эксперименты на быстрых ускорителях, начиная с экспериментов на установке “Ангара-5-1” [6], в которых был получен значительный нейтронный выход, достигавший 8∙10 ¹¹ ДД- нейтронов за импульс при уровне токов 2 –2.5 МА. В начале 2000-х проводились эксперименты на установке Z [7–9], где при уровне токов 18 МА был получен нейтронный выход 4∙10 ¹³ ДД- нейтронов и где полученные данные указывали на преимущественно термоядерный механизм генерации нейтронов. Помимо заметного укорочения импульса тока существенным отличием экспериментов на установках “Ангара-5-1” и Z от классического Z -пинча является импульсный напуск газа в разрядный промежуток, позволяющий создать неоднородное распределение плотности дейтерия в узле нагрузки и сосредоточить его основную массу вблизи оси, где происходит фокусировка плазменной оболочки. Следует отметить, что импульсный напуск газа применялся также и в плазменном фокусе [10–12], где применение импульсного напуска в качестве основного или дополнительного источника дейтерия в комбинации со стационарным напуском позволило увеличить выход нейтронов из плазменного фокуса по сравнению с работой в режиме стационарного напуска газа.

Несмотря на долгую историю Z -пинчевых и плазмо-фокусных исследований, полной ясности и окончательной модели генерации нейтронов, которая бы объясняла все экспериментальные данные и позволяла бы уверенно предсказывать результаты экспериментов, до сих пор не существует. Для моделирования динамики токово -плазменной оболочки и нейтронного выхода в плазменном фокусе успешно используется модель Ли [13, 14], однако она является феноменологической и хорошо адаптирована только для плазменного фокуса мейзеровского типа [5]. Включение кинетического описания движения заряженных частиц (например, в приближении частиц в ячейке) в области фокусировки токово -плазменной оболочки и сопряжение этого описания с МГД-подходом для всей остальной области плазменного фокуса [15–17] кажется перспективным подходом, но является очень трудоемким и затратным с вычислительной точки зрения. Поэтому представляет интерес провести расчеты работы быстрого Z -пинча с учетом импульсного напуска газа, с помощью двумерной МГД- модели [18], включающей генерацию термоядерных и ускорительных нейтронов. Наша работа посвящена расчетному исследованию экспериментов [ 6]. Такие исследования могут помочь при планировании соответствующих экспериментов на существующих [ 6, 19] или перспективных быстрых ускорителях.

2. Описание физической модели

В представленных расчетах использовались уравнения двумерной магнитной гидродинамики в цилиндрических координатах r, z с азимутальным магнитным полем, с учетом диффузии магнитного поля, потерь энергии на излучение и теплопроводности [18].

Электронное и ионное давления определялись уравнением состояния идеального газа. Для определения связи между температурой электронов, внутренней энергией и степенью ионизации использовалось уравнение Саха.

Из классических кинетических коэффициентов в наших расчетах мы учитывали только электрическое сопротивление плазмы, электронную и ионную теплопроводности, при этом в этих кинетических коэффициентах мы учитывали замагниченность электронов и ионов. В наших расчетах не учитываются некоторые двухжидкостные эффекты замагниченной плазмы, которые могут быть важны для расчетов магнитотепловых процессов при магнитоинерциальном удержании [20] и для описания динамики малоплотных Z -пинчей. В их числе явления переноса, представленные эффектом Холла, эффектами Эттинсгаузена, Нернста и Ледюка–Риги в обобщенном законе Ома и в выражении для теплового потока и описывающиеся с помощью классических кинетических коэффициентов [21, 22]. Для расчетов Z -пинчей их роль может быть более важной в области перетяжек, где погонная плотность плазмы уменьшается. Поскольку в области перетяжек плазменное β ~ 1, а плазма замагничена, в качестве критерия степени “гидродинамичности” описания плазмы можно использовать ионное погонное число [23] ${\Pi }_i\equiv N{e}^2/m_i{c}^2$, где $N=\underset{0}{\overset{\infty }{\int }}n\cdot 2\pi rdr$. При больших ${\Pi }_i$ роль этих эффектов должна быть мала. В нашем случае на стадии образования перетяжки (см. рис. 6, 10, время t = 38.742 мкс) ${\Pi }_i$ ~ 30, достаточно велико, чтобы этими эффектами можно было пренебречь. На стадии развала перетяжки (рис. 6, 10, время t = 38.744 мкс) ионное погонное число становится ${\Pi }_i$ ~ 1, однако в это время здесь действует аномальное сопротивление, эффективные турбулентные частоты столкновений сильно превышают классические кулоновские и классические кинетические коэффициенты здесь неприменимы. Поэтому описание плазмы перетяжки в это время является в значительной степени феноменологическим.

 

Рис. 6. Распределение плотности плазмы в различные моменты времени.

 

Для расчета проводимости плазмы использовались классические спитцеровские формулы с учетом рассеяния электронов на нейтралах [21]. Также учитывалось возможное появление в плазме аномального сопротивления вследствие развития нижнегибридной дрейфовой неустойчивости [24], которое включалось при условии замагниченности ионов и электронов: (ωτ) _e ≥ 1 и (ωτ) _i ≥ 1. Соответствующие аномальные поправки учитывались также при расчете электронной теплопроводности плазмы. Проведенные ранее расчеты показали, что учет аномального сопротивления является ключевым фактором в генерации нейтронов по ускорительному механизму, поскольку именно оно приводит к образованию электрических полей на оси, достаточных для ускорения ионов до необходимых энергий. Модели проводимости плазмы, основанные на развитии разных неустойчивостей, могут давать различную величину аномального сопротивления [25–27]; в данной работе использовалась модель аномального сопротивления [24], которая дает хорошее согласие с экспериментальными данными по величине нейтронного выхода, форме и положению источника нейтронов [18]. К тому же результаты полностью кинетического моделирования плазменного фокуса [28] также указывают на развитие нижнегибридной неустойчивости.

Нейтронный выход вычислялся как сумма вкладов термоядерного и ускорительного механизмов генерации нейтронов. При этом для оценок нейтронного выхода по ускорительному механизму предполагалось, что вблизи оси, на расстояниях от нее меньше трех ларморовских радиусов ускоренных ионов, часть ионов имеет возможность ускоряться в электрическом поле вдоль оси, одновременно испытывая силу торможения. Величина этой силы в зависимости от энергии иона, электронной и ионной температур плазмы вычислялась по формулам, представленным в [29, 30]. Считалось, что ускоренные ионы замещали в приосевой области α _i ~ 15% максимальной (до текущего значения z) плотности тока j _z [30]. Для начала ускорения ионов величина ускоряющей силы со стороны электрического поля вдоль оси должна была превышать силу торможения. Такие условия могут возникнуть, если электрические поля обуславливаются аномальным сопротивлением, возникающим в плазме в момент ее пинчевания на оси камеры.

При расчете нейтронного выхода по термоядерному и ускорительному механизмам дополнительно насчитывалось распределение термоядерных нейтронов по температуре ионов и распределение ускорительных нейтронов по энергии ускоренных ионов в пучке. При вычислении нейтронного выхода по ускорительному механизму использовались данные по сечению ДД-реакции из библиотеки ENDF/B-VII [31].

3. Постановка задачи и схема расчета

На рис. 1 представлена прямоугольная область (в координатах z, r), для которой проводились расчеты, размеры указаны в сантиметрах.

Область, выделенная белым цветом на рис. 1, имитирует внутреннюю, вакуумную область узла нагрузки установки “Ангара-5-1”, геометрия которого взята из статьи [6]. Отрезок А (z = 0) – анодный электрод, К — катодный электрод. Область С — сопло с выходным диаметром 2 см, предназначенное для напуска дейтерия в разрядный промежуток, область КЛ – быстродействующий дейтериевый клапан. В начальный момент времени в расчетах область заполнялась разреженным дейтерием с молярной массой А = 4 г/моль и фоновой “вакуумной” плотностью ρ _vac = (2–8)∙10 ^–11 г/см ^3, начальная температура газа T ₀ = 300 К (0.0259 эВ).

 

Рис. 1. Расчетная геометрия узла нагрузки установки “Ангара-5-1” (А — анод; К — катод; С — сопло; КЛ – быстродействующий клапан; B ₀ – значение магнитного поля на границе; V _in – скорость дейтерия на входе).

 

Область, выделенная на рис. 1 серым цветом, имитирует вещество катодного электрода. В этой области задавалась плотность твердого вещества 8.9 г/см ^3, коэффициент теплопроводности 4 Вт/(см∙К) и проводимость 1.1∙10 ²⁵ См/см. Начальная температура плотного вещества T ₀ = 0.0259 эВ.

Границы прямоугольной расчетной области Z = 0, Z = 10.2, R = 6 считались холодными, идеально проводящими стенками с нулевым значением нормальной компоненты скорости везде, кроме отрезка Z = 10.2, 0 ≤ R ≤ 0.76 (выделен синим на рис. 1) и R = 6, 0 ≤ Z ≤ 1 (выделен красным на рис. 1). На отрезке Z = 10.2, 0 ≤ R ≤ 0.76 задавалось граничное условие по газодинамике “скорость на входе (V _in)” для создания массового потока газа для имитации работы быстродействующего клапана. На отрезке R = 6, 0 ≤ Z ≤ 1 задавалось значение граничного магнитного поля HB ₀ (t), которое определяется по значению тока из осциллограммы, приведенной в статье [6]. Граница R = 0 – ось симметрии. На границах между катодным электродом и внутренней областью узла нагрузки нормальная компонента скорости считалась нулевой для исключения движения вещества катодного электрода во внутреннюю область узла нагрузки.

Расчеты проводились в два этапа. На первом, газодинамическом этапе моделировалось заполнение внутренней области концентратора молекулярным дейтерием (D ₂₎ с начальной температурой T ₀ = 300 К через быстродействующий клапан и сверхзвуковое сопло. Скорость дейтерия на входе составляла V _in = –10 ⁵ см/с в соответствии с оценками скорости струи в [6], плотность на входе варьировалась в диапазоне ρ _in = (1.7–6.8)∙10 ^{– 4} г/см ^3. Расчеты проводились с использованием газодинамической части программы с учетом теплопроводности. При расчете теплопроводности дейтерия использовались коэффициенты теплопроводности для молекулярного водорода из [32, 33]. Расчеты на первом этапе проводились в однотемпературном приближении, распределения магнитного и электрического полей, а также проводимость, джоулев нагрев и потери энергии на излучение не рассчитывались.

На втором этапе с помощью МГД-расчета моделировалась динамика плазмы под действием тока, создаваемого генератором импульсов тока установки “Ангара-5-1”. При переходе от газодинамического к МГД-расчету выполнялась повторная инициализация начальных данных. При этом полученное на первом этапе расчета распределение плотности и скорости газа использовалось в качестве начального для МГД-расчета, начальная температура газа принималась равной T _i0 = T _e0 = 5800 К (0.5 эВ), для определения давления дейтерия использовалось уравнение состояния идеального одноатомного газа с атомной массой А = 2 г/моль, остальные параметры рассчитывались в соответствии с формулами, используемыми в МГД-расчете [18].

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

$t_{\text{mhd}}=t_{\text{nozzle}}+t_{\text{delay}},$ (1)

где t _nozzle – время достижения в узкой части сопла заданной плотности, t _delay – время задержки.

Плотность дейтерия в разрядном промежутке в момент перехода к МГД-расчету и сам момент перехода t _mhd подбирались так, чтобы наилучшим образом обеспечить согласование с экспериментальными данными, приведенными в [6]. Плотность дейтерия в разрядном промежутке варьировалась путем изменения плотности на входе ρ _in, которая подбирались таким образом, чтобы при скорости дейтерия на входе V _in = –10 ⁵ см/с обеспечить достижение погонного числа частиц на срезе сопла N ₀ ≈ 2∙10 ¹⁸ см ^–1 за ≈ 30 мкс и перепад плотности от катода к аноду примерно на порядок величины [6]. В ремя перехода к МГД-расчету t _mhd варьировалось таким образом, чтобы обеспечить согласование расчетной и экспериментальной зависимости напряжения между анодом и катодом от времени на растущей части кривой. В процессе расчетов данные параметры варьировались также с целью исследования их влияния на величину нейтронного выхода. В ремя t _mhd изменялось в диапазоне от 27.5 до 38.9 мкс.

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

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

Газодинамический этап расчета проводится по явной схеме. Для аппроксимации дифференциальных уравнений используются разностная схема “против потока” [34]. Расчет электронной и ионной теплопроводности проводился по неявной схеме методом последовательных прогонок по столбцам и строкам с итерациями [35].

При расчете распределения магнитного и электрического полей также использовалась неявная схема с прогонками. При этом расчет компонент магнитного поля, связанных с альфвеновской скоростью, был совмещен с расчетом электрических полей (см. работы [30, 18, 36, 37]). Такой подход позволяет значительно уменьшить ограничение на величину счетного шага по времени, возникающего при расчете по явным схемам из-за наличия областей с низкой плотностью.

Расчеты проводились на квадратных сетках с размером ячейки h = 0.08 см, 0.04 см, 0.02 см и 0.01 см. Для ускорения счета на многопроцессорных системах с общей памятью в данной программе применялось распараллеливание кода средствами программного интерфейса OpenMP [38].

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

В данном разделе представлены результаты газодинамического этапа расчетов, на котором моделировалось заполнение дейтерием внутренней области узла нагрузки через быстродействующий клапан и сверхзвуковое сопло. На рис. 2 приведено распределение плотности инжектируемого в узел нагрузки дейтерия в различные моменты времени, полученное в расчете на сетке h = 0.01 см. Представленные ниже результаты получены при следующих параметрах: время перехода к МГД-расчету t _mhd = 38.64 мкс, плотность на входе ρ _in = 3.0∙10 ^–4 г/см ³ и соответствуют МГД-расчету, в котором достигнут максимальный нейтронный выход.

 

Рис. 2. Распределение плотности дейтерия в различные моменты времени.

 

Рис. 2 показывает, что к 5 мкс дейтерий достигает узкой части сопла и в ней формируется струя, которая к 10 мкс достигает среза сопла и начинает заполнять разрядный промежуток между анодным и катодным электродами. К 15 мкс струя достигает анодного электрода и к 20 мкс начинает распространяться радиально от оси вдоль анода; при этом формируется отраженная от анода ударная волна, которая отчетливо видна, начиная с момента времени 25 мкс. Момент времени 38.64 мкс соответствует переходу к МГД-расчету.

В целом из рис. 2 видно, что на приведенные моменты времени заполнение внутренней области узла нагрузки дейтерием является существенно неоднородным. Степень этой неоднородности демонстрирует также рис. 3, где приведены профили плотности дейтерия вдоль оси на участке 0 ≤ z ≤4 см, соответствующие моментам времени, для которых представлены распределения на рис. 2, начиная с 20 мкс.

 

Рис. 3. Профили плотности дейтерия вдоль оси в различные моменты времени.

 

Согласно данным статьи [6] погонное число частиц на срезе сопла достигает значения N ₀ ≈ 2∙10 ¹⁸ см ^–1 (соответствующая плотность 4.2∙10 ^–6 г/см ³⁾ за время ≈ 30 мкс после срабатывания клапана, а количество газа у катода может превышать количество газа у анода примерно на порядок величины. Согласно рис. 3 на момент времени 30 мкс на срезе сопла (z = 3 см) плотность равна 5.3∙10 ^–6 г/см ^3, что соответствует погонному числу частиц 2.5∙10 ¹⁸ см ^{– 1,} а отношение значения плотности на срезе сопла к ее минимальному значению вблизи анода составляет 16.8¹. Это отношение уменьшается со временем, для момента 35 мкс оно равно 11, а для момента 38.64 мкс, когда происходит переход к МГД-расчету, оно равно 8.7.

Для расчетов на сетках с h = 0.08 –0.01 см со значением плотности на входе ρ _in в диапазоне от 1.7∙10 ^–4 до 6.8∙10 ^–4 г/см ³ время достижения заданной плотности в узкой части сопла изменялось в диапазоне от 8 до 12.8 мкс. Для расчета с ρ _in = 3.4∙10 ^–4 г/см ³ это время для заданной плотности 4.2∙10 ^–6 г/см ³ составило t _nozzle = 12.7 мкс. Таким образом, в момент перехода к МГД-расчету t _mhd = 38.9 мкс используемое в экспериментах время задержки составило, согласно (1), t _delay = 26.2 мкс.

5. Результаты МГД-расчетов и сравнение с экспериментом

В данном разделе представлены результаты МГД-расчета, на котором моделировалась динамика плазмы под действием тока, создаваемого генератором импульсов тока установки “Ангара-5-1”.

На рис. 4 представлена используемая в МГД-расчете зависимость тока генератора установки “Ангара-5-1” от времени I (t), взятая из статьи [6], а также ее аппроксимация полиномом 10-й степени I _fit (t), и производная тока dI _fit / dt, полученная путем дифференцирования по времени аппроксимирующего полинома I _fit (t).

 

Рис. 4. Зависимости тока и производной тока генератора установки “Ангара-5-1” от времени.

 

На рис. 5 представлены экспериментальная (U _exp) и сглаженная расчетная (U _calc) зависимости напряжения в генераторе от времени. Экспериментальная зависимость взята из [6]. Зависимость U _calc вычислялась как сумма напряжения между анодом и катодом U _MHD и добавки, связанной с индуктивностью линии U _L. В МГД-расчете вычислялось напряжение U _MHD путем интегрирования осевой компоненты напряженности электрического поля вдоль ячеек, на которых задавалось граничное значение магнитного поля $B_0$ (см. рис. 1). Добавочное напряжение U _L = L ∙ dI _fit / dt возникает из-за ненулевой индуктивности передающей линии на участке от узла нагрузки до места установки измерителя напряжения. Индуктивность L подбиралась таким образом, чтобы обеспечить наилучшее согласование экспериментального и расчетного импульса напряжения на растущем участке кривой. Найденное значение индуктивности L = 17 нГн. Далее, поскольку расчетная зависимость U _calc (t) на временах пинчевания плазмы содержала большие высокочастотные осцилляции плазмы (возможно, связанные с численной методикой), а в реальности измерения должны сглаживать эти осцилляции, представленная на рис. 5 кривая U _calc (t) была сглажена с характерным “приборным” временем τ = 5 нс.

 

Рис. 5. Экспериментальная и расчетная зависимости напряжения от времени.

 

Сравнивая кривые, представленные на рис. 5, можно отметить, что не удалось достичь полного согласования экспериментальной и расчетной зависимостей напряжения. Подъем расчетной кривой напряжения в начале в момент времени 38.64–38.65 мкс связан с поведением производной тока, которое вызвано перегибом экспериментального тока в моменте 38.65 мкс (см. рис. 4). Но наиболее важное отличие расчетной кривой от экспериментальной состоит в том, что у расчетной кривой имеется пик (двойной) напряжения, связанный с пинчеванием (уходом тока на малый радиус) и соответствующим ростом индуктивности, что обычно имеет место в Z-пичах и плазменном фокусе. У экспериментальной же кривой такого пика нет или же он распределен на большие времена (и, значит, пинчевание происходит более медленно). Подчеркнем, что в наших расчетах мы варьировали плотность газа в широких пределах и варьировали распределение газа в момент включения тока, и даже таким образом не смогли добиться согласования расчетной и экспериментальной кривых напряжения. Можно указать некоторые причины, возможно, влияющие на отличие этих кривых. Во-первых, различие между экспериментальной и расчетной кривыми напряжения можно объяснить неучетом в расчете трехмерных эффектов, которые могут иметь место в эксперименте. В расчетах плазменного фокуса и МАГО получено хорошее согласие описания динамики плазмы в 2D-моделировании с экспериментом [18, 20, 30] и трехмерность не оказывала заметного влияния. Однако в этих экспериментах проводилась тщательная подготовка внутренней поверхности вакуумной камеры, которая обрабатывалась “тренировочными” импульсами, что уменьшало развитие азимутально несимметричных возмущений (несимметричных пробойных явлений, филаментации тока) в плазменной оболочке. Чем тщательнее велась эта подготовка, тем лучше было согласие 2D-расчетов с экспериментом. При работе систем с импульсным напуском газа проводить такие “тренировочные” импульсы затруднительно. Возможно, кроме того, что причина размытия пика напряжения на экспериментальной кривой связана с тем, что измерение коротких импульсов напряжения в эксперименте является очень сложной задачей и при этом в экспериментах иногда наблюдается растягивание «хвостов» сигнала.

На рис. 6 приведено распределение плотности плазмы на различные моменты времени до и после образования пинча. Из распределений, представленных на рис. 6, видно, что с момента времени 38.69 мкс начинается сгребание магнитным полем дейтериевой плазмы в направлении оси системы. В момент времени 38.74 мкс на оси формируется сжимающаяся плазменная перетяжка – пинч, которая к моменту времени 38.742 мкс достигает диаметра 0.2 см и значения плотности 1.9∙10 ^–4 г/см ^3, что в примерно в 18 раз превышает плотность плазмы на срезе сопла и в 160 раз превышает плотность плазмы у анода в момент перехода к МГД-расчету. В момент времени 38.744 мкс происходит разрушение плазменной перетяжки, а затем к моменту 38.748 мкс перетяжка замыкается вновь.

На рис. 7 приведено распределение температуры ионов плазмы на различные моменты времени до и после образования пинча. Распределения температуры, представленные на рис. 7, показывают, что в расчетах присутствует значительный разогрев вакуумных областей позади плазменной оболочки, со значением температуры в диапазоне от 0.05 до 1 кэВ (значение 1 кэВ соответствует искусственно введенному в программу ограничению на температуру вакуума). Температура ионов в плазменной оболочке находится в диапазоне от 3 до 5 кэВ, а в момент образования перетяжки 38.74‒38.742 мкс возрастает до ~20 кэВ. В момент развала перетяжки температура ионов в образовавшемся зазоре превышает значение 1 МэВ.

 

Рис. 7. Распределение температуры ионов плазмы в различные моменты времени.

 

На рис. 8 приведено распределение осевой компоненты напряженности электрического поля E _z на различные моменты времени до и после образования пинча. На рисунке в момент времени 38.744 мкс на врезке с увеличением показана область 0.5 ≤ z ≤ 1.5, 0 ≤ r ≤ 0.6, в которой происходит развал перетяжки и формирование зазора. В этой области на оси осевая компонента напряженности электрического поля E _z достигает величины примерно 10 МВ/см, что при длине зазора порядка 2 мм может давать напряжение на зазоре до 2 МВ и, соответственно, энергию ускоренных в этом зазоре ионов до 2 МэВ. Эти ионы, согласно используемой нами модели ускорительного механизма [18, 30], вызывают генерацию нейтронов, взаимодействуя с ионами плазмы, которые в нашей модели мы считали покоящимися. На моменты 38.74‒38.744 мкс, в которые происходит сжатие и разрушение плазменной перетяжки, достигнутые в узком канале наибольшие значения E _z порядка 10 МВ/см близки к результатам кинетического моделирования плазменного фокуса [17].

 

Рис. 8. Распределение осевой компоненты напряженности электрического поля E _z в различные моменты времени.

 

На рис. 9 приведен полученный в МГД-расчете график зависимости напряжения на оси от времени. Из рисунка видно, что напряжение на оси достигает максимума 3.3 МВ в момент времени 38.745 мкс. Полуширина импульса напряжения составляет 1.1‒1.8 нс. Временной интервал 38.744‒38.748 мкс, на который приходится основной пик напряжения на оси, соответствует разрушению перетяжки (см. рис. 6). Замыкание перетяжки приводит далее к заметному снижению напряжения на оси.

 

Рис. 9. Зависимость напряжения на оси от времени.

 

На рис. 10 приведено распределение величины Br на различные моменты времени до и после образования пинча. На увеличенных фрагментах, представленных на рисунке, видно, что ток Br, несущий магнитное поле, достаточно близко подходит к оси; так на момент времени 38.742 мкс, соответствующий образованию плазменной перетяжки, он достигает радиуса 0.1 см (при этом магнитное поле достигает величины B ≈ 4 МГс), а в момент 38.744 мкс, соответствующий сбросу массы из перетяжки, этот ток проникает еще глубже – до радиуса 0.02–0.03 см (магнитное поле достигает величины B ≈ 5.5–6 МГс). В последующие моменты времени наблюдается уменьшение магнитного поля в области пинчевания.

 

Рис. 10. Распределение величины B r в различные моменты времени.

 

Интегральный выход нейтронов в расчете составил 1.5∙10 ¹¹ ДД-нейтронов за импульс, при этом на термоядерную составляющую нейтронного выхода приходилось 9∙10 ¹⁰ ДД-нейтронов, а на ускорительную, соответственно, 6∙10 ¹⁰ ДД-нейтронов. Такое соотношение термоядерных и ускорительных нейтронов существенно отличается от результатов расчетов нейтронного выхода в плазменном фокусе, где обычно ускорительных нейтронов в 100–1000 раз больше, чем термоядерных [18, 30].

На рис. 11 приведено полученное в расчете пространственное распределение выхода термоядерных и ускорительных нейтронов из единицы объема по расчетной области. Области, выделенные на рис. 11а, б рамкой, ниже даны с увеличением. Из рис. 11в, г видно, что длина области генерации ускорительных нейтронов вдоль оси z составляет примерно 0.7 см и существенно больше осевого размера области генерации термоядерных нейтронов, которая составляет примерно 0.2 см. При этом радиусы области генерации близки и составляют примерно 0.03 см для ускорительных нейтронов и 0.02 см для термоядерных нейтронов.

 

Рис. 11. Пространственное распределение нейтронного выхода из единицы объема : термоядерные нейтроны (а); ускорительные нейтроны (б); термоядерные нейтроны (увеличена область 0 ≤ r ≤ 0.1, 0 ≤ z ≤ 1.5) (в); ускорительные нейтроны (увеличена область 0 ≤ r ≤ 0.1, 0.5 ≤ z ≤ 2) (г).

 

На рис. 12 приведены интенсивности нейтронного выхода для термоядерных (dN _терм / dt) и ускорительных нейтронов (dN _уск / dt) в зависимости от времени. Графики, представленные на рисунке, показывают, что положения максимумов интенсивности нейтронного выхода практически совпадают и соответствуют моменту развала плазменной перетяжки пинча (см. рис. 6) и максимуму импульса напряжения на оси (см. рис. 9). При этом максимальное значение dN _терм / dt примерно в 2 раза превышает максимальное значение dN _уск / dt. Полуширина интенсивности термоядерного нейтронного выхода dN _терм / dt составляет 1 нс, а для интенсивности ускорительного нейтронного выхода dN _уск / dt – 3 нс.

 

Рис. 12. Интенсивности нейтронного выхода для термоядерных (dN _терм / dt) и ускорительных нейтронов (dN _уск / dt) в зависимости от времени.

 

На рис. 13 представлено распределение нейтронного выхода термоядерных нейтронов по температуре ионов. Это распределение имеет два локальных максимума в области 5.8 кэВ и 10.8 кэВ, при этом 90% нейтронного выхода (от 5% до 95% от полного выхода) лежит в диапазоне температур 3–42 кэВ.

 

Рис. 13. Распределение нейтронного выхода термоядерных нейтронов по температуре ионов.

 

На рис. 14 представлено распределение нейтронного выхода ускорительных нейтронов по энергии налетающих ионов. Это распределение имеет основной пик в диапазоне энергий 40–300 кэВ и высокоэнергетический “хвост” 300–1200 кэВ. При этом 90% нейтронного выхода (от 5% до 95% от полного выхода) лежит в диапазоне 55–900 кэВ. Особенностью распределения выхода ускорительных нейтронов по энергии ионов пучка является его существенная немонотонность с большим количеством локальных минимумов и максимумов. Полученный в МГД-расчете диапазон энергий ускоренных ионов, приводящих к генерации ускорительных нейтронов, хотя и весьма широк, в целом согласуется с оценками, приведенными в [6], где данный диапазон энергий составляет 200–500 кэВ.

 

Рис. 14. Распределение нейтронного выхода ускорительных нейтронов по энергии налетающих ионов.

 

В процессе проведения МГД-расчетов было исследовано влияние размера ячейки расчетной сетки на газодинамические параметры расчета и величину нейтронного выхода. В качестве газодинамических параметров выбрано отношение значения плотности на срезе сопла к ее минимальному значению вблизи анода в момент перехода к МГД-расчету, а также время перехода от газодинамического этапа к МГД-расчету t _mhd. Результаты расчетов представлены в табл. 1. В расчетах, представленных в табл. 1, для каждого h подбирали время перехода к МГД-расчету t _mhd для наилучшего согласования с данными [ 6], плотность на входе ρ _in = 3.38∙10 ^–4 г/см ³ и остальные параметры не менялись.

 

Таблица 1. Влияние размера ячейки сетки на результаты расчета.

h, см	ρ max / ρ m in	t mhd, м кс	N терм, нейтр. ×10 10	N уск, н ейтр. × 10 10	N, нейт р. × 10 10
0.08	11.3	27.5	3.87	3.56	7.43
0.04	13.5	32.6	3.65	3.84	7.49
0.02	9.8	36.7	5.79	4.81	10.6
0.01	8.5	38.9	5.84	3.28	9.12

Примечание : ρ _max / ρ _min – отношение значения плотности на срезе сопла к ее минимальному значению вблизи анода в момент перехода к МГД-расчету; N _терм – выход термоядерных нейтронов; N _уск – выход ускорительных нейтронов; N – полный выход нейтронов.

 

По результатам, представленным в табл. 1 можно отметить, что сходимость по времени t _mhd проявляется достаточно медленно. Такое поведение может быть объяснено использованием газодинамической схемы первого порядка точности по пространству и недостаточным количеством ячеек в области сужения сопла, так для сетки с h = 0.08 см на сопло приходится всего 4 ячейки. Нейтронный же выход по термоядерному и ускорительному каналам ведет себя с измельчением сетки немонотонно, что свидетельствует о неустойчивости процессов пинчевания. Надо сказать, что немонотонное поведение нейтронного выхода, свидетельствующее о неустойчивости, обнаруживается также при изменении других параметров.

В процессе проведения МГД-расчетов было исследовано влияние времени перехода от газодинамического этапа к МГД-расчету t _mhd и начальной плотности газа в струе в момент перехода к МГД-расчету на величину нейтронного выхода. Вариация времени перехода от газодинамического этапа к МГД-расчету t _mhd достигалась путем изменения в программе задержки t _delay. Вариация начальной плотности струи в момент перехода к МГД-расчету достигалась путем изменения плотности на входе ρ _in при фиксированной задержке t _delay. При таком изменении параметров t _mhd изменяется слабо и при переходе к МГД-расчету заполнение плазмой разрядного промежутка меняется незначительно, а плотность возрастает в примерно одинаковое число раз во всех точках расчетной области². В расчетах на сетке h = 0.02 и 0.01 см параметры варьировались в следующих диапазонах: плотность на входе ρ _in от 1.7∙10 ^–4 до 6.8∙10 ^–4 г/см ^3, время перехода к МГД-расчету t _mhd от 34.7 до 38.9 мкс. Результаты расчетов нейтронного выхода на сетках h = 0.02 и 0.01 см в зависимости от времени t _mhd и плотности на входе ρ _in приведены в табл. 2. Остальные параметры, используемые в МГД-расчете, не менялись.

 

Таблица 2. Зависимость выхода нейтронов от параметров t _mhd и ρ _in.

h, см	ρ in, г/ см 3 × 10 –4	t mhd, м кс	N терм, нейтр. × 10 10	N уск, н ейтр. × 10 10	N, нейт р. × 10 10
0.02	3.38	34.7	3.97	3.67	7.64
0.02	3.38	36.4	6.75	5.15	11.9
0.02	3.38	37.0	6.33	4.75	11.08
0.02	3.38	37.4	7.52	5.16	12.68
0.02	3.38	38.7	5.66	2.86	8.52
0.01	3.38	38.5	3.70	2.04	5.74
0.01	3.38	38.9	5.84	3.28	9.12
0.01	1.69	38.5	2.75	2.58	5.33
0.01	2.37	38.6	5.79	3.72	9.51
0.01	2.70	38.6	5.36	3.34	8.7
0.01	3.04	38.6	8.95	5.95	14.9
0.01	3.38	38.7	5.76	4.52	10.28
0.01	4.06	38.7	3.61	1.79	5.4
0.01	6.76	38.8	1.97	1.75	3.72

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

 

Величина интегрального нейтронного выхода при изменении перечисленных выше параметров в указанных диапазонах варьировалась от 3.7∙10 ¹⁰ до 1.5∙10 ¹¹ ДД-нейтронов за импульс (см. табл. 2). Соотношение термоядерного и ускорительного выхода нейтронов в проведенных расчетах варьировалось в диапазоне от 1.1 до 2, при этом в большинстве расчетов термоядерный выход превышал ускорительный. Все остальные результаты расчетов при варьировании параметров газового напуска в указанных пределах менялись не сильно.

Нейтронный выход в экспериментах [6] варьировался в диапазоне от 5∙10 ¹⁰ до 8∙10 ¹¹ ДД-нейтронов за импульс в зависимости от задержки между сигналом датчика давления газа в сопле и моментом пуска генератора. Следует отметить, что в расчетах так и не удалось достичь экспериментального нейтронного выхода на уровне 8∙10 ¹¹ ДД-нейтронов за импульс [6], также как и совпадения расчетной и экспериментальной кривых напряжения, несмотря на то, что плотность газа варьировалась в достаточно широких пределах, что может свидетельствовать о необходимости корректировки модели ускорительного механизма [18, 30] применительно к расчетам систем с импульсным напуском газа. Кроме того, как показали проведенные расчеты, выход нейтронов достаточно чувствителен к изменению варьируемых параметров плотности в времени перехода к МГД-расчету. Так при изменении параметров t _mhd и ρ _in в указанных диапазонах отношение максимального выхода к минимальному ~4 (см. табл. 2), хотя общая картина течения плазмы в узле нагрузки при этом сохраняется. Таким образом, даже небольшое отличие в начальных данных от экспериментов [6] или неточность в подобранных параметрах плотности и времени перехода к МГД-расчету, описывающих начальное распределение плазмы в узле нагрузки, могут приводить к отличию выхода нейтронов в МГД-расчетах от результатов, полученных в экспериментах [6].

6. Заключение

Проведены двумерные МГД-расчеты динамики быстрого Z -пинча на установке “Ангара-5-1” с учетом импульсного напуска газа и генерации нейтронов с учетом термоядерного и ускорительного механизмов. Расчеты проводились в два этапа: на первом, газодинамическом этапе моделировалось заполнение разрядного промежутка дейтерием, на втором этапе происходил переход к МГД-расчету и моделировалась динамика плазмы под действием тока, создаваемого генератором импульсов тока установки “Ангара-5-1”. В расчете подобраны параметры – плотность газа на входе в клапан и время перехода к МГД-расчету, которые позволили удовлетворительно описать имеющиеся экспериментальные данные по распределению плотности газа в разрядном промежутке и напряжению на установке “Ангара-5-1” [ 6]. В результате МГД-расчетов получены следующие параметры пинча: диаметр ~0.2 см, сжатие по плотности более 160, температура в пинче ~20 кэВ, магнитное поле на радиусе пинча ~4–6 МГс. Размер области генерации нейтронов z ×r составляет ~ 0.7×0.03 см для ускорительных нейтронов и ~ 0.2×0.02 см для термоядерных нейтронов.

В проведенных расчетах исследовано влияние времени перехода от газодинамического этапа к МГД-расчету и начальной плотности газа в струе в момент перехода к МГД-расчету на величину нейтронного выхода. Расчеты показали, что при вариации времени перехода к МГД-расчету в диапазоне от 34.7 до 38.9 мкс и плотности на входе от 1.7∙10 ^–4 до 6.8∙10 ^–4 г/см ^3, расчетный нейтронный выход находится в диапазоне (0.4–1.5)∙10 ¹¹ ДД-нейтронов за импульс. При этом в экспериментах выход варьировался в диапазоне (0.5–8)∙10 ¹¹ ДД-нейтронов за импульс в зависимости от задержки между сигналом датчика давления газа в сопле и моментом пуска генератора. Следует отметить, что в расчетах так и не удалось достичь экспериментального нейтронного выхода на уровне 8∙10 ¹¹ ДД-нейтронов за импульс, что может свидетельствовать о необходимости корректировки модели ускорительного механизма применительно к расчетам систем с импульсным напуском газа.

Максимальный выход нейтронов в расчете составил 1.5∙10 ¹¹ ДД-нейтронов за импульс, при этом на термоядерную составляющую нейтронного выхода приходилось 9∙10 ¹⁰ ДД-нейтронов, а на ускорительную, соответственно, 6∙10 ¹⁰ ДД-нейтронов. Такое соотношение термоядерных и ускорительных нейтронов существенно отличается от результатов расчетов нейтронного выхода в плазменном фокусе, где обычно ускорительных нейтронов в 100–1000 раз больше, чем термоядерных. Анализ распределения ускорительного выхода нейтронов по энергии налетающих ионов позволил выявить основной пик в диапазоне энергий 40–300 кэВ и высокоэнергетический “хвост” 300–1200 кэВ. При этом 90% нейтронного выхода лежит в диапазоне 55–900 кэВ. Полученный в расчете диапазон энергий ускоренных ионов, приводящих к генерации ускорительных нейтронов, хотя и весьма широк, в целом согласуется с оценками, приведенными в [6], где данный диапазон энергий составляет 200–500 кэВ.

Авторы выражают благодарность участникам семинара ТРИНИТИ 22.03.2023, и персонально В.П. Смирнову, Е.В. Грабовскому и В.Е. Черковцу, за обсуждение возможности проведения исследований, представленных в настоящей работе.

 

¹ В процессе проведения расчетов плотность дейтерия на входе ρ _in подбиралась таким образом, чтобы наилучшим образом обеспечить согласование с экспериментальными параметрами дейтериевой струи, приведенными в [6], т.е. обеспечить достижение значения погонного числа частиц N ₀ ≈ 2∙10 ¹⁸ см ^–1 за время ≈ 30 мкс и перепад плотности от анода к катоду ~10, однако из-за ограничений по времени это делалось лишь для расчетов на грубых сетках с размерами ячейки h = 0.08–0.04 см. При переходе на подробные сетки с размерами ячейки h = 0.02–0.01 см, значение плотности на срезе сопла и перепад плотности от катода к аноду несколько изменились.

² Небольшое изменение t _mhd и не строго одинаковое распределение плотности по объему обусловлено учетом теплопроводности в газодинамическом расчете.

Авторлар туралы

S. Garanin

Russian Federal Nuclear Center—All-Russian Scientific Research Institute of Experimental Physics

Email: vyudolinskij@vniief.ru
Ресей, Sarov

V. Dolinskii

Russian Federal Nuclear Center—All-Russian Scientific Research Institute of Experimental Physics

Хат алмасуға жауапты Автор.
Email: vyudolinskij@vniief.ru
Ресей, Sarov

Әдебиет тізімі

Қосымша файлдар