Релятивистская разностная схема для расчёта динамики частиц

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

Рассматриваемая ниже схема возникла в ходе работ по транспортировке сильноточного электронного пучка в линейном индукционном ускорителе ЛИУ-5 [38–40]; физическое обоснование формул изложено там.

Основные величины

В начальный момент времени \(t\) заданы декартовы координаты \(i\)-й частицы \(\overrightarrow{r_{i}} = (x_i, y_i, z_i)\) и её начальный импульс \(\overrightarrow{p_{i}} = (p_{x_i}, p_{y_i}, p_{z_i})\). Задан минимальный шаг по времени \(\delta t\).

Для удобства обезразмерим физические величины, обозначенные дальше тильдой. \begin{equation} \widetilde{\overrightarrow{r_i}} = \frac{\overrightarrow{r_i}}{c \delta t}, \end{equation} \begin{equation} \widetilde{\overrightarrow{E_i}} = \frac{\overrightarrow{E_i}e}{mc} \delta t, \end{equation} \begin{equation} \widetilde{\overrightarrow{H_i}} = \frac{\overrightarrow{H_i}e}{mc} \delta t, \end{equation} где \(c\) — скорость света, \(e\) — заряд частицы, \(m\) — масса частицы.

Алгоритм

• Вычисление электрических полей \(\overrightarrow{E_i}\) в точках, занимаемых частицами.

• Приращение импульса. \begin{equation} \overrightarrow{p_i} = \overrightarrow{p_i} + 2\cdot\overrightarrow{E_i}, \end{equation} где коэффициент означает, что приращение импульса выполняется за весь интервал времени \(2\delta t\).

• По новым импульсам вычисляются новые скорости. \begin{equation} \overrightarrow{v_i} = \frac{\overrightarrow{p_i}}{\gamma_i}, \end{equation} где \(\gamma_i = \sqrt{1 + \overrightarrow{p_i}^2}.\)

• По новым скоростям вычисляются новые координаты. \begin{equation} \overrightarrow{r_i} = \overrightarrow{r_i} + 1\cdot\overrightarrow{v_i}, \end{equation} где коэффициент означает, что приращение координаты выполняется за интервал времени \(\delta t.\)

• Вычисление магнитных полей \(\overrightarrow{H_i}\) в новых точках, занимаемых частицами.

• Вычисляются новые значения скоростей, полученные после поворота в магнитном поле.

  \[ b_1 = 1 - \frac{H_i^2}{\gamma_i^2}, b_2 = 1 + \frac{H_i^2}{\gamma_i^2}, b_3 = 2\cdot \frac{\overrightarrow{v_i}\cdot\overrightarrow{H_i}}{\gamma_i}, \]

  \[ \overrightarrow{f_i} = 2\cdot \frac{\overrightarrow{v_i}\times \overrightarrow{H_i}}{\gamma_i}, \]

  \begin{equation} \overrightarrow{v_i} = \frac{\overrightarrow{v_i}b_1 + \overrightarrow{f_i} + \frac{\overrightarrow{H_i}}{\gamma_i}b_3}{b_2}. \end{equation}

Все члены здесь построены на векторе поворота \(\overrightarrow{t} = \overrightarrow{H_i}/\gamma_i\), поэтому в \(b_1\) и \(b_2\) входит его квадрат. Это проверка, заложенная в записи схемы: поворот, взятый в таком виде, сохраняет модуль скорости точно, тогда как любая другая степень \(\gamma\) делает выражение неоднородным, и частица начинает набирать или терять энергию в чисто магнитном поле, не совершающем над ней работы.

• По новым скоростям снова вычисляются новые координаты. \begin{equation} \overrightarrow{r_i} = \overrightarrow{r_i} + 1\cdot\overrightarrow{v_i}, \end{equation} где коэффициент означает, что приращение координаты выполняется за интервал времени \(\delta t.\)
• По новым скоростям вычисляются новые импульсы. \begin{equation} \overrightarrow{p_i} = \overrightarrow{v_i}\gamma_i. \end{equation}

Один цикл схемы, собранный из перечисленных шагов, выполняется за интервал времени \(2\delta t.\)

Python

REDPIC оформлен в виде библиотеки: пользователь описывает тракт и пучок и запускает трекинг, а разностная схема выполняется внутри.

import redpic as rp
import numpy as np
import holoviews as hv
hv.extension('matplotlib')

import warnings
warnings.filterwarnings('ignore')

Проверим установленную версию: расчёт, который потребуется воспроизвести через год, помечается версией кода, на которой он получен.

rp.__version__

Настройки графиков

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

%output size=100 backend='matplotlib' fig='png' dpi=300
%opts Curve Scatter [aspect=3 show_grid=True]
%opts Curve (linewidth=1 alpha=0.7 color='blue')
%opts Scatter (alpha=0.7 s=0.5)

Задаём параметры ускорительного тракта

Тракт, заданный границами по \(z\) и шагом сетки, описывается так же, как в KENV. Расчёт ведётся от 0.7 до 5 метров с шагом один сантиметр.

acc = rp.accelerator.Accelerator(0.7, 5, 0.01)

В тракте расположены два ускоряющих модуля, настроенных на −1.1 МВ/м.

#              Unique name,  z-position [m],  Ez [MV/m],  Ez(z) profile
acc.add_accel('Acc. 1',      4.096,          -1.1,         'Ez.dat')
acc.add_accel('Acc. 2',      5.944,          -1.1,         'Ez.dat')

Кроме того, в тракте расположены семь фокусирующих соленоидов. У первого, намотанного навстречу остальным, поле отрицательно.

#                 Unique name,  z-position [m],  Bz [T],  Bz(z) profile
acc.add_solenoid('Sol. 1',      0.450,          -0.0580,   'Bz.dat')
acc.add_solenoid('Sol. 2',      0.957,           0.0390,   'Bz.dat')
acc.add_solenoid('Sol. 3',      2.107,           0.0250,   'Bz.dat')
acc.add_solenoid('Sol. 4',      2.907,           0.0440,   'Bz.dat')
acc.add_solenoid('Sol. 5',      3.670,           0.0400,   'Bz.dat')
acc.add_solenoid('Sol. 6',      4.570,           0.0595,   'Bz.dat')
acc.add_solenoid('Sol. 7',      5.470,           0.0590,   'Bz.dat')

Три элемента из девяти расположены вне интервала счёта. Это Sol. 1 на 0.450 м, находящийся до его начала, а также Sol. 7 и Acc. 2 на 5.470 и 5.944 м, вышедшие за его конец. Профиль поля, заданный для каждого элемента, имеет конечную ширину, и своими краями он заходит внутрь окна. Однако полностью в расчёт попадает только один ускоряющий модуль из двух, и приведённые ниже графики, построенные на этом интервале, следует читать с поправкой: за \( z = 5 \) м расчёт не производится.

Метод compile() сшивает профили всех девяти элементов в две функции, \(E_z(z)\) и \(B_z(z)\), которые далее опрашиваются в каждой точке траектории каждой частицы.

acc.compile()

Опишем оси и построим оба поля вдоль тракта. Магнитное поле, измеряемое в теслах, переводится в гауссы умножением на \(10^4\).

dim_z  = hv.Dimension('z',  unit='m')
dim_Ez = hv.Dimension('Ez', unit='MV/m', label='$E_z$')
dim_Bz = hv.Dimension('Bz', unit='Gs', label='$B_z$')

z  = acc.z
z_Ez = hv.Curve((z, acc.Ez(z)), kdims=dim_z, vdims=dim_Ez)
z_Bz = hv.Curve((z, acc.Bz(z)*1e4), kdims=dim_z, vdims=dim_Bz)

(z_Ez + z_Bz).cols(1)

Собранную структуру полезно распечатать целиком: так проще заметить пропущенный или неверно расположенный элемент.

print(acc)

Accelerator structure.
	Solenoids:
	[ 0.45 m, -0.058 T, Bz.dat, Sol. 1, 0.0 m, 0.0 rad, 0.0 m, 0.0 rad] 
	[ 0.957 m, 0.039 T, Bz.dat, Sol. 2, 0.0 m, 0.0 rad, 0.0 m, 0.0 rad] 
	[ 2.107 m, 0.025 T, Bz.dat, Sol. 3, 0.0 m, 0.0 rad, 0.0 m, 0.0 rad] 
	[ 2.907 m, 0.044 T, Bz.dat, Sol. 4, 0.0 m, 0.0 rad, 0.0 m, 0.0 rad] 
	[ 3.67 m, 0.04 T, Bz.dat, Sol. 5, 0.0 m, 0.0 rad, 0.0 m, 0.0 rad] 
	[ 4.57 m, 0.0595 T, Bz.dat, Sol. 6, 0.0 m, 0.0 rad, 0.0 m, 0.0 rad] 
	[ 5.47 m, 0.059 T, Bz.dat, Sol. 7, 0.0 m, 0.0 rad, 0.0 m, 0.0 rad] 
	Accelerating modules:
	[ 4.096 m, -1.1 T, Ez.dat, Acc. 1, 0.0 m, 0.0 rad, 0.0 m, 0.0 rad] 
	[ 5.944 m, -1.1 T, Ez.dat, Acc. 2, 0.0 m, 0.0 rad, 0.0 m, 0.0 rad] 
	Quadrupoles:
	Correctors x:
	Correctors y:

У ускоряющих модулей поле напечатано в теслах, хотя задано оно было в МВ/м. Это недоработка метода __str__ в библиотеке, который выводит одну и ту же подпись для соленоидов и для ускоряющих секций. Числа верны, а единицы нет. Печать состояния объекта, написанная наспех, затем читается годами, и разбираться в чужих единицах, нигде не оговорённых, приходится по исходному коду.

Задаём параметры пучка

В отличие от KENV, где пучок описывался одной огибающей, здесь необходимо задать распределение, по которому разыгрываются макрочастицы; выбрано равномерное по радиусу распределение. Параметры пучка: электроны, энергия 1.32 МэВ, ток 500 А, радиус 48 мм, угловой разброс 70 мрад, нормализованный эмиттанс 200 мм·мрад. Продольный «радиус» 3.5 метра означает, что сгусток значительно длиннее самого тракта, и пучок здесь практически непрерывен, то есть неотличим от постоянного тока.

beam = rp.beam.RadialUniformBeam(
    type=rp.constants.electron, 
    energy = 1.32,          # MeV
    current = 0.5e3,  # A
    radius_x = 48e-3, # initial r (m)
    radius_y = 48e-3, # initial r (m)
    radius_z = 3.5,
    radius_xp = 2*35.0e-3,     # initial r' (rad)
    radius_yp = 2*35.0e-3,     # initial r' (rad)
    x  = 0.0e-3,   # horizontal centroid position (m)
    xp = 0.0e-3,     # horizontal centroid angle (rad)
    y = 0,          # vertical centroid position (m)
    normalized_emittance = 200e-6 # m*rad
)

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

beam.generate(10_000)

2023-08-31 21:35:33,649 - redpic.beam.base - INFO - Generate a radial-uniform beam with 10000 particles

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

beam.df

Опишем оси.

dim_x = hv.Dimension('x', unit='m', range=(-0.1, 0.1))
dim_y = hv.Dimension('y', unit='m', range=(-0.1, 0.1))
dim_z = hv.Dimension('z', unit='m', range=(acc.z_start, acc.z_stop))
dim_px = hv.Dimension('px', unit='MeV/c', label='$p_x$')
dim_py = hv.Dimension('py', unit='MeV/c', label='$p_y$')

и соберём четыре проекции, построенные по одному датафрейму: поперечное сечение \(x\)–\(y\), продольный вид \(z\)–\(x\) и два фазовых портрета, \(x\)–\(p_x\) и \(y\)–\(p_y\).

beam_x_y = hv.Scatter(beam.df, kdims=[dim_x, dim_y])
beam_z_x = hv.Scatter(beam.df, kdims=[dim_z, dim_x])
beam_x_px = hv.Scatter(beam.df, kdims=[dim_x, dim_px])
beam_y_py = hv.Scatter(beam.df, kdims=[dim_y, dim_py])

WARNING:param.Scatter: Chart elements should only be supplied a single kdim
2023-08-31 21:35:33,772 - param.Scatter - WARNING - Chart elements should only be supplied a single kdim
WARNING:param.Scatter: Chart elements should only be supplied a single kdim
2023-08-31 21:35:33,789 - param.Scatter - WARNING - Chart elements should only be supplied a single kdim
WARNING:param.Scatter: Chart elements should only be supplied a single kdim
2023-08-31 21:35:33,791 - param.Scatter - WARNING - Chart elements should only be supplied a single kdim
WARNING:param.Scatter: Chart elements should only be supplied a single kdim
2023-08-31 21:35:33,801 - param.Scatter - WARNING - Chart elements should only be supplied a single kdim

Предупреждения HoloViews можно не принимать во внимание: Scatter ожидает одну ключевую размерность, а передаются две; на изображение это не влияет.

(beam_x_y + beam_z_x + beam_x_px + beam_y_py).cols(2)

В поперечном сечении видно круглое пятно радиусом 48 мм. Фазовые портреты представляют собой наклонные полосы, а не облака; это видно и по числам в таблице, напечатанной сразу после генерации, где у первых же частиц \(p_x/x \approx 2.9\) для всех подряд. Пучок стартует расходящимся, поперечный импульс частицы почти пропорционален её координате, а толщина полосы соответствует эмиттансу, неустранимому никакой фокусировкой.

Параметры сгенерированного пучка распечатаем и сверим с заданными.

print(beam)

Beam parameters:
            Type	electron
            Distribution	radial-uniform
            Particles	10000
            Current	500.0 A
            Energy	1.32 MeV
            Total momentum	1.7582483378931433 MeV/c
            Rel. factor	3.583175582013202
            Radius x	48.0 mm
            Radius y	48.0 mm
            Radius z	3.5 m
            Radius x prime	70.0 mrad
            Radius y prime	70.0 mrad
            Horizontal centroid position	0.0 mm
            Vertical centroid position	0.0 mm
            Horizontal centroid angle	0.0 mrad
            Vertical centroid angle	0.0 mrad
            Normalized emittance x	200.0 mm*mrad
            Normalized emittance y	200.0 mm*mrad

Запуск моделирования

Полный импульс 1.7582 МэВ/c и \(\gamma = 3.58\) получены пересчётом заданной энергии и точно с ней согласуются.

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

rp_sim = rp.solver.Simulation(beam, acc)
rp_sim.track()

z = 4.98 m (99.5 %) 

Визуализация результатов моделирования

Расчёт дошёл до 4.98 метра, то есть до конца тракта. Результат rp_sim.result организован в виде словаря, где ключом служит координата вдоль тракта, а значением — датафрейм со всеми частицами, собранными в этом сечении. Это покадровая запись, хранимая целиком в памяти.

Напишем функцию, которая строит один кадр.

def plot(i):
    df = rp_sim.result[i]
    rp_z_x = hv.Scatter(df, kdims=[dim_z, dim_x], label='redpic')
    return rp_z_x

и объединим кадры в HoloMap — интерактивное изображение, управляемое ползунком по \(z\). Перемещение ползунка позволяет увидеть не только изменение размера пучка, но и перемешивание частиц внутри него и поведение его края.

items = [(i, plot(i)) for i in list(rp_sim.result.keys())]

holomap = hv.HoloMap(items, kdims = ['z'])

hv.output(holomap, widget_location='bottom')

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

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