Новый эффективный алгоритм расчета траектории заряда и его тестирование на точных решениях Текст научной статьи по специальности «Физика»

УДК 533.95

П. В. Сецко, М. Н. Мельник, О. В. Мингалёв, И.В. Мингалёв

НОВЫЙ ЭФФЕКТИВНЫЙ АЛГОРИТМ РАСЧЕТА ТРАЕКТОРИИ ЗАРЯДА И ЕГО ТЕСТИРОВАНИЕ НА ТОЧНЫХ РЕШЕНИЯХ

Аннотация

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

Ключевые слова:

расчет траектории заряженной частицы, ускорение заряженных частиц.

P. V. Setsko, M. N. Melnik, O. V. Mingalev, I. V. Mingalev

NEW EFFICIENT ALGORITHM FOR CALCULATION OF THE CHARGED PARTICLES TRAJECTORIES AND ITS TESTING ON THE EXACT SOLUTION

Abstract

The paper describes an algorithm for calculating of the charged particles trajectories having a 2nd order accuracy. The algorithm makes the most efficient to calculate a large number of trajectories and the ability to simulate the acceleration of charged particles to relativistic energies.

Keywords:

the calculation of the trajectory of the charged particle, the acceleration of charged particles.

Введение

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

Оптимальное решение проблемы состоит в создании гибкого алгоритма расчета траектории заряда, который в области нерелятивистской энергии использует оптимизированный для классического приближения алгоритм с достаточно большим шагом по времени, а при приближении к релятивистской энергии переходит на расчет траектории по неявной схеме 2-го порядка точности с мелким шагом по времени. Первый из указанных алгоритмов был разработан нами более 15 лет назад, широко используется с 2000-го года, был всесторонне протестирован и на практике доказал свою правильность и эффективность (см., например, [1-6]). Этот алгоритм основан на точном решении системы уравнений движения заряда Ньютона — Лоренца для нерелятивистского случая в постоянных (по пространству и времени) электрическом и магнитном полях, имеет нулевую фазовую ошибку и снимает условие необходимости достаточно большого числа шагов по времени на период циклотронного вращения частиц. В случае плавно изменяющихся полей и достаточно сильного магнитного поля алгоритм допускает шаг по времени в десятки или более периодов циклотронного вращения, что позволяет получить значительную экономию вычислительных ресурсов.

В работе предложена новая неявная схема 2-го порядка точности для численного решения задачи Коши для системы уравнений движения заряда Ньютона — Лоренца в релятивистском случае. Для ее тестирования получено точное решение задачи Коши в случае постоянных и параллельных друг другу магнитного и электрического полей в явной векторной форме и для произвольных начальных условий. Отметим, что это решение обсуждается в монографии [7], но там рассматривался только частный случай, когда начальная скорость заряда ортогональна магнитному полю, а также была выбрана параметризация, которая дает для решения только неявное соотношение и не позволяет получить решение в явной форме. Тестирование показало хорошую точность предложенной схемы.

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

Алгоритм расчета траектории частиц

Для расчета траектории частицы с зарядом е и массой покоя т0

необходимо численно решать задачу Коши для системы уравнений Ньютона — Лоренца, которая в системе единиц СИ в релятивистском случае имеет вид:

¿х (г) т0У

-Г~ = v (0 , Р(v )= , 0 . 2 ,

= * (£ (х (г), г И V (г) хВ (х (г), г)])

где с — скорость света в вакууме. В классическом случае |у|/с □ 1 эта система может быть представлена в виде

(1)

^ = V С )' ^ = £ (* <х ('),'И V (') «В (' О, ')])■ (2)

Начальные условия имеют вид:

х ) = х0, V ) = V0. (3)

Кинетическая энергия в электрон-вольтах Ж( V ) связана с модулем

скорости V =| V | формулами

„ {

1

W(|v|) =

mc

-1

eW

eW

v m0c

+ 2

eW

v m0c

+1

(4)

и в классическом пределе переходит в кинетическую энергию для классического

/ \ 1 2 v

приближения: W( v ) ^— mV ПРИ--^ 0. Отметим, что из уравнений (1)

2e с

dw( |v (t) |) . . . ,

и (4) вытекает уравнение ——-—= (E (x (t), t); v (t)).

В ходе расчета траектории для экономии вычислительных ресурсов при значении параметра ( = |v |/c < (0 = 0.02 (которому соответствует

кинетическая энергия W ~ 188 кэВ) выполняется временной шаг в рамках

классической системы уравнений (2), а при ( > (0 = 0.02 выполняется

временной шаг в рамках релятивистской системы уравнений (1).

Для сокращения вычислительных затрат и удобства системы уравнений (1) и (2)

нужно представить в оптимальной безразмерной (нормализованной) форме. Для

этого надо выбрать масштабы e = |e|, m = m0, V = c, а также произвольное значение какого-либо одного из следующих масштабов: t, x, В, EE, а остальные из них определить из соотношений x = ct, t = m0 j(eB) , EE = cB . В

результате, если обозначить безразмерный заряд как e' = sign (e) = ±1, и по традиции обозначить безразмерный импульс через u :

u = p' = p/(m0c) = v/Vl-|vР/c2 = v'/yj 1 -|vf ,

то получим следующие безразмерные формы систем (1) и (2) (далее штрихи опущены):

dx (t)

= v (t), u (v ) = -

ф-|v|2 '

^ = e'( E ( x (t), t) + [ v (t) XB ( X (t), t)]),

= v (t), d^l = e' (E (x (t), t) + [ v (t) XB (x (t), t)]).

(5)

(6)

Для численного решения задачи Коши (6), (3) используется разработанный авторами и показавший очень высокую эффективность алгоритм

e

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

Решение задачи Коши (6), (3) для последующего изложения удобно обозначить как

*(0 = x(t, t0, x0, v0 ), v (t) = v(t, t0, x0, v0 ),

(7)

то есть эти векторные функции удовлетворяют следующим уравнениям и начальным условиям:

ax(t, t0, x0, v0)

dt

V (t, t0, x0, v0)

av(t,t0,x0,v0) , . . ,

(, ^^ = e'E(x(t,t0,x0,v0),t)■

at

+ e

v{

(t, t0, x0, v0) XB (X(t, t0, x0, v0), t)

X(t0, t0, x0, v0 ) = x0, V(t0, t0, x0, v0 ) = v0.

Разработанный авторами алгоритм имеет 2-й порядок точности, допускает переменный шаг по времени, а также явный и неявный варианты и основан на точном решении задачи (6), (3) в случае постоянных полей Е = const, В = const. Это решение с использованием стандартные обозначений:

B Г E х b 1 . . .

b = —, vE =i-*, ®c = B , E | = (E; b) b,

B ' E — ' c 1 1, 1 1 ( ; ) ,

^ =(v; b) b, vL = v - v - v£ , т = t -10 можно представить в виде

2

x(t, t0, x0, v0) = x0 + (v£ + vc0) т + T— e'EL +

> х vL ],

(8)

+

— sin (юст)- e

"(l- cos (®CT))

V(t, t°, x°, v°) = vE + + те'Е, + v°L cos (а>ст) - e' [b x v°L ] sin (а>ст).

(9)

Рассмотрим алгоритм в терминах послойного перехода г ^ г = г +т в полях Е(х, г) и В(х, г) , заданных в любой точке во все моменты времени

г > г0 . Обозначим /12 = /0 + г/2, для функций вида _/(г) и Е(х, г) введем

обозначения /а=/(га) и Еа=Е(х(га),га) , а = 0,1,! I. Стартовая итерация

1.1. Находим поля Е°= Е(х°,г0) , В°= В(х0,г0) , затем по ним и формулам (8), находим Ь° , V® , = (е°-,Ъ° , уп° = ,

1.2. По первой формуле в (9) находим х12(0) = х( г12,г0, х0, V0)

,[хVI]

V!=V0 - *0 - ^

в постоянных полях Е0, В0 :

12(0) 0/0 0\77 гт~"0 ■

1 ( ' = х + (V0 + V|0 )- +—е Е0 + —81 2 8 1 а

00 С 7

ч 2 ,

с

1 - С08

С ^

V 2 „

Находим поля е1/2(0)=Е (х12(0), г12), В1/2(0)=В (х12(0), г12), а затем

по ним и формулам (8), находим Ь

1/2(0)

1/2(0)

а

12(0)

е1

=( ; ЬУ2(0)) Ь^2(0) ,,12(0^^.^2(0^ Л12(0) ¿/2(0) = V0 - ^^ _

= ( Е

= (V0; Ь12(0))

Ь1

= V - V,

1.3. По первой формуле в (9) находим х1(0)=х(г1, г0, х0, V0 ) в постоянных полях Е12(0), В12(0):

х1(К) = х0 +(v12C") + ^1/2(К) )т + — е е1{2(к) +

,1/2(к)

+

-8т (®12(К)г)

[ ь

2

1/2(к) ^,1/2«) ^

X V

1/2(к)°"Ч "с V 1/2(к)

сс с

(1 - С08 (С)),

(10)

где К = 0 . Находим поля Е1(0)= Е (х1(0), ), В1(0)= В (х1(0), ^ ) . 1.4. Находим при К = 0 поля

е1/%+1) = i (е0 + е1(к)) , в12(к+1) = 1 (в0 + в1(К)) ,

далее по ним и формулам в (8), находим Ь

12(к+1) х> 12(К+1) ^112(К+1) /7У2(К+1)

Е

, с

Е1

1/2(к+1 = (V0. ь12(К+1) ) Ь12(К+1) ^

1/2(к+1^, 0_,1/2(К+1^, 1/2(К+1)

= V - V

- V,

. Затем по

формуле (10) при к = 1 находим х1(1)=х(/1, г0, х0, V0) в постоянных полях Е1/2(1), В^2(1). Далее находим поля Е1(1) = Е (х1(1), г1) , В1(1) = В (х1(1), г1), а также

параметр точности 81 =

х х ^

. Если достигнута заданная точность 8 :

Е

ь

с

V

< 5, то процесс закончен. В противном случае выполняется обычная итерация.

II. Обычная итерация К ^ К +1. По уже найденным х1((К), E

гК*)

B

■1(К)

E

выполняя действия, описанные в пункте I.4, находим х1К+!)

.1(К+1) ^1(К)

1(к+1)

X

— X

B , а также параметр относительном точности о ± = параметр сходимости ук+1 = 5г+1 — 5К . Итерации выполняются до достижения заданной относительной точности 5^+1 < 5. При этом если условие сходимости итерационного процесса YK+\ = 5^+1 — 5 < 0 на некоторой итерации нарушается или количество итераций превышает заданное максимально допустимое число к0, то текущий шаг по времени т уменьшается в два раза —

т ^ т/2 и описанный процесс начинается сначала.

III. Блок завершения шага. Выполняем присваивание x1 = x ,

E1 = E^+1), B1 = B^+1). Затем находим поля

E12 = - (E0 + E1 ) , B12 = - (B0 + B1 ) ,

далее по ним и формулам в (7) находим b1/2, v]l2, со12, E/2,

12 ( 0 »1/2 \ ,1/2 .,1/2 .,0 1/2 12 г, „ , ,оч

vr' =(v ; b ) b' , v[ = v — vj — v^ . Затем по второй формуле в (8) находим v1 = V(t, t0, x0, v0 ) в постоянных полях E12, B12 :

V1 = vf + vV2 + ^ ^ 2 + v12 c o s ( COlJ 2 Г ) — e '

-,1/2 v,,1/2'

sin (с/Г).

На этом шаг по времени закончен.

Для численного решения задачи Коши (5), (3) в релятивистском случае авторами разработана неявная схема 2-го порядка точности. Для ее

формулировки в терминах послойного перехода г0 ^ г1 = г0 + г необходимо ввести функцию

Г « ) =

Av) ф+1«\

и представить систему (5) в виде

f^=Г « (' )) « (' ).

= л 1 — \ v

(11)

^ = ✓ ( E ( x (t ), t ) + Г( « (t ))[ « ( t ) XB ( x (t ), t )]) ,

(12)

а также использовать вытекающие из нее следующие уравнения: = — ✓(Г( « (t )))3 ( « (t ) ; E( X (t ), t )) ,

й2 х (г)

йг£

+

е'г(и(г)) (е(х (г), г) +г(и (г)) [и (г) хв(х (г), г)] )-

йг[и (г)) йг

и

(г).

Схема основана на следующих конечноразностных соотношениях:

х1 - х0 _ 1 (Г 0 7 2

(

(Г0и0+Г1и1)-

и1 - и0 1

йи

у йг у

йи

у йг у

7

12

А

2 х^

и 2 х^

0

у Л 2 у

-О (72 ),

у йг2 у

О(74),

(15)

7 2

I. Стартовая итерация

1.1. Находим поля Е°= Е(х0,г0) , В°= В(х0,г°) , затем по ним и

формулам (12)-(14) находим

,1(0)

помощи экстраполяции:

1(0) 0 , 1-0 0, 7 х1 ;= х +7Г и +--

Г

у йг у

1.2. Находим х1^ с точностью О (т3) и Г1 ' с точностью

№ г Л0 йГ) ( й2

2

у

1(-)

йг

йг

и 2 х^0

2

у йг 2 у

+

О 73), Г1(-)=Г0 +7 ,.1(0)

Т0

у йг у

+

О 7) при

О7)-

находим

поля

Далее, используя

е1(0)= е(х1(0) , г), в1(0)= в (х1(0), г) с точностью О (г ).

1.3. Находим «1(0) с точностью О (73) следующим образом. Из соотношения (15) и второго уравнения в (12) можно получить следующее векторное уравнение относительно и1 :

Ш + О 7),

и1 +7 е Г1

2

^0 , 17/1

В1 хи1

где Ш = Ш и+ Ш1, Ш 0= и0 +7

2

с и

у йг у

Ш 1=7е 'Е1.

(16)

Это уравнение можно представить в форме системы линейных уравнений с матрицей +С)| = где I —единичная матрица и

сии = [шхи] .

0

Обратная

матрица

М"

имеет

вид

+ <Ь ^ = —-—^ - <Ь + су ® су ^ , где (о®(о— диада, образованная вектором

ю.

1+ |ю||

Поэтому решение векторного уравнения (16) дается формулой и = ^1 + ео) = ^ 2 (IV-[со х1Г] + (1Г:со)со)

Применяя эту формулу к уравнению (16) при ю = - Г1 Д1, получим

выражение для и через остальные слагаемые:

ч 2

W --в'г1

В1 XW] + [ -г I (W;В1)В1

и =

1 + [гг1

2

В1

-О(-3) .

Подставляя в эту формулу уже найденные г1(-), в1(0), Е1(0) , получим

следующие формулы:

W = W 0+ W1(0), W1(0) = - в'Е1(0).

и1(0) =

W - - в Г1(-) [В1(0) xW ] + Г1(-) 1 (W; В1(0)) В

?1(0)

+ [-Г1(-) |в1(0^2

о (-).

(17)

1.4. Используя полученные и1(0) и В1(0), Е1(0), с точностью О (г3)

находим

Г1(0) по формуле (10),

а также

^1(0)

V 1г J

по формуле (11) и

'<12 ^

1 J

по формуле (13). Присваиваем 50 = 1.

II. Обычная итерация к ^ к +1. Пусть уже известны х1К) , и1((к).

Г1«» в1К) Е1(к)

^ 1и 11(К) Г < 2^1(К)

у ёг J

V1г J

с точностью

О(-).

11.1. Находим л:1^к+1) с точностью О (г4) из соотношения (14):

х1^ = х0 +1(г0и0+г1КУК))-1-2^ ' 12

Гг ,,2\1(к) Г ,2 Л01

£х

vёг j

<Гх

vёг j

О(г4) .

2

с точностью

Далее находим поля £№1>= Е(х'1"+Ч,г1), В(х1"+11,г1)

0(Г4)'

11.2. Находим и"1» с точностью О (т3) по формуле (17) полностью аналогично (17):

и

Ж — те' ГЧК) _ 2_

Ж = Ж и+ Ж

ВК"+1) х Ж

1("+1) Ж 1("+1) = Т е'£1"+1)

+ '§ Г

1(")

( В1("+1);Ж ) В

1("+1)

(

1 + |ТГ1(")

ч2Л

В1

1("+1)

о т).

11.3. Используя полученные и

1"+1)

В

1("+1)

Е

1("+1)

с точностью

О (т3) находим Г1("+1) по формуле (10), а также

^ 1("+1)

у ёг J

по формуле (11) и

/ ,2 л1("+1)

ё х

у ёг j

по формуле (13).

х 1("+1) - х 1(")

11.4. Находим параметр относительной точности дг+1 = и параметр сходимости Ук+1 = дк+1 — дк. Итерации выполняются до достижения заданной относительной точности д^+1 < д. Если условие сходимости итерационного процесса Ук+\ = дк+х —дк< 0 на некоторой итерации нарушается или количество итераций превышает заданное максимально допустимое число ", то текущий шаг по времени т уменьшается

в два раза: т ^ т/2 и итерационный процесс начинается сначала. В конце выполняем присваивание х1 = х 1("+1), и1 = и1("+1}, Г1 = Г1("+1), Е1 = Е1^1»,

В1 = В1"1»

С ёи ^ | ёи

ёг

\1("+1)

ёг

(^Л1 Г

х

ёг2

х

у ёг2 J

Отметим, что изложенный алгоритм расчета решения задачи Коши (5), (3) в нерелятивистском случае широко используется нами с 2000-го года, был всесторонне протестирован и на практике доказал свою правильность и эффективность. Поэтому необходимо протестировать только новую неявную схему для релятивистского случая.

Для этого тестирования лучше всего подходит точное решение задачи Коши (1), (3) в случае постоянных параллельных друг другу магнитного и электрического полей, когда заряд испытывает ускорение электрическим полем. Отметим, что это решение обсуждается в монографии [7], но там

2

и

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

Тестирование метода расчета траекторий

Рассмотрим задачу Коши для безразмерной системы (11) в случае постоянных и параллельных полей В = Const, Е = Enb = Const (где

b = B/|B|). Ее можно представить в виде:

dx (t) u (t) du (t)

dt +1u(t)|2 ' dt

(t0 ) = x0, u (t0 ) =

= e

Enb +

[ u (t)

x (t ) = x

Обозначим

= u

(18)

y( u) = yj 1 +1 u|2

i ,> u = (u;b),

// (н;А)А. и и и . т I /".

Полученное нами точное решение задачи (18) удобно представить в следующей форме. Для фактора Лоренца получается формула:

/(г) = ^ 1 +1и° |2+ (и0 + е'Е г)' = ^(/°)2 + е'Е_г(2ии° + е'Е г) , где / = ^ 1 +1 н° |2.

Для фазы циклотронного вращения получаются формулы

dФ(r) IB / л ,|B|

-(-L = -Br, Ф(т) = e'B ln

dr у\т) E^

U0

+ e'Ej т + у(т)

л

b x u,

sin(Ф(г)) .

Для безразмерного импульса получается формула, аналогичная формуле для скорости в (8):

= +е'Е1тЬ + и°±со5(ф(т))-е'

Для координаты получается формула:

В 44 \в\

Были проведены расчеты решения задачи Коши (18) по изложенной неявной схеме 2-го порядка в широком диапазоне релятивистских начальных

условий и значений параметра Е\/В1 на время Т = 40@0, где

г о

x (t ) = x0 + ^ (г[т)-г0 )b + u- sin (ф(т))-e

L о b x u.

(l - cos (ф(т)) ) .

>

0O = 2щ-= —:—¡— . Эти расчеты показали высокую точность

/ dz Я

рассматриваемой схемы при достаточно малом шаге по времени. Заключение

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

Благодарности. Работа выполнена при поддержке подпрограммы V.15. «Динамика разреженной плазмы в космосе и лаборатории» Комплексной программы фундаментальных исследований ОФН РАН.

Литература

1. Мингалёв О. В., Мингалёв И. В., Мингалёв В. С. Двумерное численное моделирование динамики мелкомасштабной неоднородности в околоземной плазме // Космические исследования. 2006. Т. 44, № 5. С. 416-427.

2. Численное моделирование плазменного равновесия в одномерном токовом слое с ненулевой нормальной компонентой магнитного поля / О. В. Мингалёв [и др.] // Физика плазмы. 2007. Т. 33, № 11. С. 1028-1041.

3. Несимметричные конфигурации тонкого токового слоя с постоянной нормальной компонентой магнитного поля / О. В. Мингалёв [и др.] // Физика плазмы. 2009. Т. 35, № 1. С. 85-96.

4. Кинетические модели токовых слоев с широм магнитного поля / О. В. Мингалёв [и др.] // Физика плазмы. 2012. Т. 38, № 4. С. 329-344.

5. Formation of self-organized shear structures in thin current sheets / H. V. Malova [et al.] // J. Geophys. Res.: Space Physics. 2015. Vol. 120. DOI: 10.1002/2014JA020974.

6. Мингалёв И. В., Сецко П. В. Расчет траекторий заряженных частиц в магнитосфере Земли // Высокоширотные геофизические исследования: труды 4-й Школы молодых ученых (29-30 октября 2015 г., г. Мурманск) / Полярный геофизический ин-т. 2015. С. 157-162. (Труды Кольского научного центра РАН. Серия «Гелиогеофизика», вып. 1).

7. Ландау Л. Д., Лифшиц Е. М. Теория поля: 10 т. Т. 2: Теоретическая физика. Изд. 8-е. М.: Физматлит, 2012. 536 с.

Сведения об авторах Сецко Павел Владимирович

программист, Полярный геофизический институт, г. Апатиты E-mail: [email protected]

Мельник Михаил Николаевич

младший научный сотрудник, Полярный геофизический институт, г. Апатиты E-mail: [email protected]

Мингалёв Олег Викторович

к. ф.-м. н., заведующий сектором, Полярный геофизический институт,

г. Апатиты

E-mail: [email protected] Мингалёв Игорь Викторович

д. ф.-м. н., старший научный сотрудник, Полярный геофизический институт, г. Апатиты

E-mail: [email protected]