CC BY
40
Вычислительные технологии
Том 11, № 3, 2006
ГИБРИДНАЯ МОДЕЛЬ РАСПРОСТРАНЕНИЯ АЛЬФВЕНОВСКОЙ ВОЛНЫ СДВИГА В БЕССТОЛКНОВИТЕЛЬНОЙ ПЛАЗМЕ*
Г. И. ДудниковА Институт вычислительных технологий СО РАН, Новосибирск, Россия
e-mail: [email protected]
Л. В. ВшивковА, Р. Рэнкин Университет Альберты, Эдмонтон, Канада
A numerical model of propagation of shear Alfven waves in collisionless plasma is presented. The model is based on the kinetic description for electrons in plasma and the hydrodynamical approximation for ions (hybrid model). Vlasov kinetic equation is solved by the modified particle in cell method that allows reducing numerical noise inherent to PIC-method. Algorithms for updating the particle position (that ensures plasma quasi-neutrality) and an algorithm for calculation of the electric field are proposed.
Введение
Изучение процессов генерации и распространения альфвеновских волн сдвига в бесстолк-новительной плазме имеет большое значение для прикладных и фундаментальных исследований, поскольку данные процессы связаны с переносом импульса и энергии, развитием неустойчивостей и ускорением заряженных частиц. Альфвеновская волна сдвига представляет собой поперечную магнитогидродинамическую волну, которая имеет частоту ниже циклотронной частоты ионов. Данный тип волн существует повсеместно, начиная от Солнца, магнитосферы Земли (и других планет) и кончая лабораторными установками термоядерного синтеза [1-4]. Особое значение имеют волновые процессы в магнитосфере Земли, где распространение волн сопровождается магнитными суббурями, полярными сияниями и нарушением радиосвязи. Результаты измерений в авроральной зоне, выполненных на спутниках, показали, что наблюдаемые низкочастотные электромагнитные волны представляют собой альфвеновские волны сдвига [5]. Эти волны являются источником мелкомасштабных флуктуаций тока, который играет существенную роль в комплексной динамике ионосферы Земли. Исследование того, как распространяются альфвеновские волны и как происходит взаимодействие волна — электроны, необходимо для понимания динамики резонансов в магнитосфере Земли. При этом наиболее важно изучение процессов, когда
* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 04-01-00850).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.
поперечный размер альфвеновской волны существенно меньше, чем плазменный скин-слой или ларморовский радиус ионов. В каждом из этих случаев (называемых соответственно инерциальным и кинетическим) для описания эффектов необходимо привлекать кинетические модели, поскольку гидродинамические модели становятся неприменимыми. Существенная разница в пространственно-временных масштабах для электронов и ионов приводит к трудностям использования полностью кинетического описания и применения метода частиц (PIC) при численной реализации модели. Иногда при решении задач на основе такого подхода используется модельное отношение масс электронов и ионов, как, например, в работе [6], где m¿/me = 16. Альтернативой полному PIC-моделированию является использование гибридных (комбинированных) моделей, которые существенно сокращают время счета. Традиционные гибридные модели объединяют кинетические уравнения для ионов и магнитогидродинамические уравнения для электронов и имеют длинную историю успешного применения при решении ряда проблем бесстолкновительной плазмы [7-9]. Однако в задачах распространения альфвеновских волн сдвига важную роль играют эффекты их взаимодействия с электронами окружающего фона, которые приводят к необходимости кинетического описания электронной компоненты плазмы и создания другого типа комбинированных моделей.
В настоящей работе представлена гибридная модель, в которой ионный компонент плазмы описывается магнитогидродинамическими уравнениями, а электроны — кинетическим уравнением Власова. Данная модель может быть использована для численного моделирования распространения альфвеновских волн и динамики заряженных частиц в магнитосфере Земли.
1. Постановка задачи
Рассмотрим задачу распространения альфвеновской волны в полярной зоне магнитосферы Земли в следующей постановке. В начальный момент времени £ = 0 в цилиндрическую область 0 < г < Я, 0 ^тах, заполненную плазмой с постоянной плотностью п = п0
и магнитным полем дипольного типа, входит альфвеновская волна с амплитудой А и частотой ш. При своем движении она возмущает максвелловскую функцию распределения электронов, формируя группу быстрых частиц, которые движутся в направлении к Земле. Исследуем низкочастотные колебания квазинейтральной плазмы (п = пе = п), состоящей из электронов и ионов водорода, поэтому током смещения пренебрегаем. В данной постановке предполагается, что в направлении распространения альфвеновской волны имеется достаточно сильное магнитное поле, которое делает движение электронов практически одномерным. В цилиндрических координатах ось г направлена вдоль силовых линий магнитного поля. Таким образом, движение электронов происходит только в направлении оси г.
Выпишем исходную систему уравнений предлагаемой гибридной модели. Для описания электронов используется кинетическое уравнение Власова, которое с учетом сделанных предположений имеет следующий вид:
д/ + дf Е д/ 0 (1)
т +дГ - ^ Е дШ = 0' (1)
где / = /(£, г, Ш) — функция распределения электронов; уг — г-компонента скорости частиц; е, те — заряд и масса электрона; Ег, Вг — г-компоненты напряженности элек-
трического и магнитного полей; Ш = ^т^2 + МВГ — кинетическая энергия частицы;
2
„ те^2_ ,
М = ——— — магнитный момент (задается для каждой частицы в начальный момент 2ВГ
времени и сохраняется в течение всего расчета); — скорость электронов в поперечном направлении.
Движение ионной компоненты рассматривается в приближении магнитной гидродинамики:
¿П + ДУ [пЩ = 0; (2)
dVi у7 + miUi— = -vpi + en dt
E + 1V x B
c
(3)
где V, П и р — скорость, плотность и давление ионной компоненты плазмы, р = р0п/; 7 — показатель адиабаты; тг — масса иона; ^ = {У, У^, У}, В = {Вг, 0, Вг}, Е = {Ег, Е^, 0}, V = , 0, 0}. К этим уравнениям добавляются уравнения Максвелла
^ = -Ц- (4)
гоШ = 4Ппе (У - К) , (5)
Здесь У =< ve > — средняя скорость электронной компоненты, У = {Уег, 0, 0}.
Кинетическое уравнение Власова для электронов решается методом частиц в ячейках, поэтому вместо данного уравнения используются уравнения его характеристик. С учетом сделанного замечания запишем исходную систему уравнений (1)-(5) в цилиндрических координатах (г, и безразмерных переменных:
¿г ¿Ш 1 ^ . .
= ^, ~Ж = -Ер; (6)
"БГ + ГТИ (rnVr) + -— (nV,) = 0; (7)
dn + 1 d ( у ) + 1 d
dt rdr r r
E, = — V Br + Vr Bz; (8)
"TTT + Vr^--1---= — V,Br; (9)
dt dr r
dV dV V dV V2 1 dP
+ У dV^ + ^ dV" — — = —1 dP + Er + V,Bz; (10)
dt dr r r ndr
dBz = 1 dE — 1A (rE),
dt r r dr ^ '
V, = —1 fz, Ver = Vr - ^. (11)
n dr nr
В качестве нормировочных выбраны следующие характерные величины: no, Во, ро, VA = Bo/V4nmino, Wpi = v/4nnoe2/mi, L = c/wPi, to = L/Va, Eo = (Va/C)Bо, в = me/mi. В соответствии с постановкой задачи начальные данные имеют следующий вид:
n = no, Er = E, = 0, Vr = V, = Vz = 0, Br = Br (r) = ^, Bz = 0. (12)
При решении задачи для всех функций выбираются периодические граничные условия по < и невозмущенные граничные условия при г = гт;п. На границе г = гтах задаются условия входа альфвеновской волны
к
Е^ = А вт(ш£), Вх =--А вт(ш£),
ш
к2 к2 В
К =--А ео8(ш£), К = —^ А в1п(ш£), (13)
шп0 ш2п0
где А — амплитуда волны, а ш и к связаны соотношением [12]
ш
2_ к2^2
1 + к2Л2'
ку — компонента к в направлении магнитного поля; Ле — длина волны. Граничные условия для остальных функций при г = гтах имеют вид
дп = 0 дЕ = 0 — = 0
дг дг дг
В^ = 0, Ех = 0. (14)
2. Описание алгоритма
В расчетной области гт;п < г < гтах, 0 < < < <тах введена равномерная сетка с шагами Л,1, Л,2 по осям < и г соответственно. Функция Вх определяется в узлах сетки (< = ¿Л,ь гк = к^2,), функции , Е^, Вг, ^ — в точках (^¿, гк-1/2 = (к - 0.5)^^, функции Ут, Ет, В^]г — в точках (^¿-1/2 = (г — 0.5)^1, гк), а функции п, У, Ег, ^, р — в серединах ячеек (<£¿-1/2 = (г — 0.5)Л,Ь гк-1/2).
На первом этапе для каждой частицы решается уравнение
^ = — 1 Ег — МдВт, ^ в Г дг
которое получается из уравнений (6) в предположении, что магнитный момент каждой частицы определяется перед началом расчета и остается постоянным. Для решения уравнения используется конечно-разностная схема 0(т, Л|):
„,т+1 „,т 1 Вт , — Вт
Т = — вЕт'*-1/2'к — м н2 ,
гт+1 _ гт
--— = ^т+1.
Т
Средние значения скоростей и плотности плазмы определяются по формулам
пе,г-1/2,к = ^ Я^ (<¿-1/2 — <') Ят (гк — О) ,
3
]С3 (<¿-1/2 — <з) Ят (гк — гз)
Ует,г-1/2,к
ne,¿—1/2,к
где j — номер частицы. Ядро частицы имеет вид
Rr(r) н i 0 — Э при |r| -"2'
0 при |r| > h2.
Остальные уравнения решаются конечно-разностным методом. Для удобства записи последующих формул введем обозначения:
—1,k—1/2) ,
K-.i—1/2,k—1/2 = ^ (Vr,i-1/2,fc + Vr,i-1/2,fc-1) j — 1,fc—1 /2
,i,k+1/2
—1,k+1/2) j
Vr,i,fc—1/2 = 4 (Vr,i—1/2,fc + Vr,i—1/2,fc—1 + Vr,i+1/2,fc + Vr,i+1/2,fc—1) ,
а r \ fi— 1/2,k—1/2 — fi—3/2,k—1/2j где 1/2,k—1/2 > 0j
A/i—1/2,k—1/2 = S r- f f> , n
[ /i+1/2,k—1/2 — /i—1/2,k—1/2, где l/^,i_1/2,fc—1/2 < 0,
а (• J fi— 1/2,k—1/2 — fi—1/2,k—3/2j где K,i— 1/2,k—1/2 > 0j
Ar /i—1/2,k—1/2 = S f t f> ' n
[ Ji—1/2,k+1/2 — Ji—1/2,k—1/2, где Vr,i—1/2,k—1/2 <
Далее с помощью конечно-разностных схем определяются: 1) плотность ионов плазмы
т
m+1 _ m
nm+1 = nm __v/m A (rn)m
ni—1/2,k—1/2 = ni—1/2,k—1/2 h2rfc 1/2 Vr,i—1/2,k—1/2Ar (r n)i—1/2,fc—1/2
T V/m A nm
h^'k 1/2 i—1/2'fc— 1/2A^' i— 1/2,fc—1/2 Tnm / Vm _ Vm Vm _ Vm
' ni— 1/2,k—1/2 I V,i,k—1/2 V,i—1,k—1/2 Vr,i—1/2,k Vr,i—1/2,k—1
+ rk—1/2
rk-1/2 h1 h2 2) r-компонента средней скорости ионов
Vm+1 _ V m _ T т/m a vm
Vr,i—1/2,k = Vr,i—1/2,k h1rk 1/2,kA Vr,i—1/2,k
т T
2
__!_т/m A V m 1 _!_ / V/m
7 Vr,i—1/2,k Ar Vr,i—1/2,k + r l V,i—1/2,k h2 rk
T /^pm+1 _ pm+1 + T f Em + Т/m Rm^ •
h^r+z1^ Vpi—1/2>k+1/2 pi—1/2'k—1/V + T rr + > 1/2,k •
3) компонента электрического поля
Em+1 _ _т/m r 1 т/m+1 Rm ;
E^,i,k—1/2 = Vz,i,k—1/2Br,i,k—1/2 + Vr,i,k—1/2Bz,i,k—1/2;
4) z-компонента средней скорости ионов
Vm+1 _ V m __T_V/m A V m
Vz,i—1/2,k—1/2 = Vz,i—1/2,k—1/2 h1rk 1/2 1/2,k—1/2Vz,i—1/2,k—1/2
— — Vт+1 Л Vт
^ Уг,г-1/2,к-1/2Л: ^,¿—1/2^—1/2
-1/2 - т(СВг)г_
5) г-компонента напряженности магнитного поля
Вт+1 _ Вт +
В.г,г,к = В* г к +
т
,г,г,к
Гк
Е
тт г,г+1/2,к - Е:,г—1/2,к
Гк Л-2
г—1/2,к-1/2
(гЕ^) - (гЕ„)
г,к—1/2
т+1
6) г-компонента средней скорости электронов
вт+1 _ вт+1
1 В.г,г,к В.г,г—1,к
кт+\/0 „ = кт+1о,-
е:,г— 1 /2,к г,г— 1/2,к -т+1 Г
"г- 1/2,к 'к
7) компонента У^ средней скорости ионов
V т+1
1/2
Вт+1 _ вт+1
1 В.г,г, к В* г к—
.г, г, к— 1
" т+1 "г , к—1/2
^2
(15)
Система конечно-разностных уравнений имеет первый порядок аппроксимации по времени и пространству. Алгоритм вычисления г-компоненты электрического поля Е: приведен ниже.
2.1. Вычисление Ег
Основную сложность в данном алгоритме представляет вычисление электрического поля Е:, поскольку оно определяется из уравнения движения для электронов с использованием метода частиц в ячейках. Для каждой частицы выполняется уравнение
¿V: = _1 ег — мдВ
в дг
Если в каждой ячейке просуммировать уравнения для всех частиц и результат разделить на количество частиц, находящихся в ячейке, получим уравнение
^ < ^ > 1 „, дВ:
: = — Е: - < М > :
в дг
где < IV >= — средняя скорость частиц, находящихся в ячейке, а < М > — средний магнитный момент. Так как поле Е:, действующее на частицу, интерполируется в местоположение частицы, средние значения скорости в узлах сетки зависят от значений Е: в окружающих узлах. Поэтому для определения Е: получается система линейных алгебраических уравнений, в каждом уравнении которой между собой связаны 15 значений. Для каждого узла получаем свое уравнение вида
15
У: ^ ^ арЕ:
Р=1
Здесь ар — коэффициенты, возникающие при интерполяции электрического поля в местоположении данной частицы.
Для решения полученной системы уравнений члены со значениями в трех центральных точках пятнадцатиточечного шаблона переносятся в левую часть системы, а остальные — в правую, и данная система решается методом прогонки для каждого значения индекса г.
2.2. Корректировка положения частиц
В рассматриваемой модели плотность ионов в ионной компоненте плазмы определяется из уравнения движения. Из-за квазинейтральности плазмы плотности ионов и электронов должны совпадать (пг = пе), однако восстановление плотности электронов в соответствии со стандартным методом частиц по формуле
Пг-1/2,к-1/2 = ^ тЯ (<£¿-1/2 - я (гк-1/2 - Г.) ,
где
Я(<)= / 1 — ^, когда М< Н, [ 0, когда |<| > Н,
не позволяет удовлетворить требуемому условию квазинейтральности. Это связано с дискретностью временного шага и ограниченным количеством частиц. Поэтому необходимо дополнительно корректировать положение частиц, что и было осуществлено в нашей модели.
Рассмотрим соответствующий алгоритм корректировки в одномерном случае. Пусть частицы в какой-то момент времени имеют координату х. из интервала [хг-1/2, £¿+1/2]. Тогда масса частиц в ячейке с номером г — 1/2 определяется по формуле
т .
Мг-1/2 = (Х*+1/2 — £3^
.
где т — масса одной частицы. Обозначим координату частиц после сдвига, корректирующего положение частиц через х. = х. + ^, где ^ — соответствующий сдвиг. Масса частиц перепишется следующим образом:
Мг-1/2 = (хг+1/2 — £ — .
.
Для частиц из левого интервала, т.е. х. € [хг-3/2,хг-1/^, имеем
т
Мг-1/2 = у (х. — хг-3/2),
Н
М-1/2 = Н (х. + ^-1 — х-3/2).
.
Здесь, координата после сдвига будет иметь вид х. = х. + ^-1. Обозначим через х] все х. € [хг-1/2, хг+1/2], а через х2 — все х. € [х»-3/2,х»-1/2]. Далее, суммируя Мг-1/2 и Мг'-1/2 для данных интервалов, получаем
т
= X
— хг - 3/2) + г+1/2 — х. )
2 - 1 7,Ж2 7,Ж ■
М -1/2
т
У> — х»-3/2 + ^-1) + ^(х,+ 1/2 — х. — ^)
И, составляя их разницу, подсчитываем ДМг-1/2:
ДМ,_1/2 = М,_1/2 - М_1/2 = ^
+ ^Ы,)
'
- (Жг_1^г_1 - N п
т. е.
п
- Жг_1^г_1 =--ДМ,_1/2)
'
где N1 — число частиц в интервале Рг = (хг_1/2) Х+1/2).
п
Обозначим — ДМ,_1/2 = Е,_1/2,
'
г, и перепишем последнее выражение
гг_ 1 +
г—1/2
О, г = 1)...) г^
Поскольку количество ячеек равно гт, а интервалов Рг — гт +1, получим систему, состоящую из гт уравнений для гт + 1 неизвестного (г = 0,...,гт). Для однозначного определения необходимо дополнительное условие, которое может быть выбрано из следующих:
0
это означает, что средние смещения всех частиц равны нулю;
г=0
б) г0 находится из минимума функционала Ф ^ ^ г2;
г=0
в) смещения в граничных ячейках определяются пропорционально дисбалансу плотности в граничных ячейках.
В двумерном случае корректировка ведется с помощью итераций по £ и г, причем по координате £ используется первый способ (а) определения г0, а по координате г — третий — (в). Необходимость итераций связана с тем, что при определении масс в ячейках используются значения масс в соседних ячейках в перпендикулярном направлении.
2.3. Устранение счетного шума
В связи с вынужденно ограниченным числом модельных частиц приходится сталкиваться с явлениями, не имеющими отношения к физике моделируемых процессов. Рассмотрим алгоритм, который позволяет существенно уменьшить счетные шумы.
Средняя скорость электронов в модели вычисляется по формуле Р1С-метода
Уг,г_1/2Л = —-- У^Ш?VЯ (£г_1/2 - £,) Я (гЛ - Г,), (16)
пг_1/2>к ?
Пг_1/2,к_1/2 = ^ 'Я (£г_1/2 - £,) Я (г&_1/2 - Г?) . 3
Флуктуации в скоростях и положениях частиц приводят к тому, что средние скорости частиц, вычисленные по формуле (16), являются негладкими функциями и не совпадают со скоростью электронов, вычисленной по формуле (15) (рис. 1, кривая 1). Хотя амплитуда наблюдаемых колебаний уменьшается с увеличением числа частиц в ячейке N как 1 /N2, но эти колебания всегда приводят к аналогичным колебаниям электрического поля ЕГ, которое вычисляется по скоростям электронов. Проблема шума всегда существовала в методе
п пп?
г
Рис. 1. Скорость электронов Ver вдоль радиуса т.
частиц в ячейках. Различные способы уменьшения уровня счетного шума предлагаются в работах [9-11] и, как правило, определяются спецификой решаемой задачи. Предложенный нами метод подавления шума состоит в следующем.
1. Предположим, что функция распределения средних скоростей частиц Ver, вычисленная по формуле (16), является на шаге m гладкой.
2. За шаг т эти частицы, имея различные скорости, переместятся в новое положение. Если вновь подсчитать средние скорости, то Vr будет негладкой функцией, даже если на частицы не действуют никакие силы. Разница средних скоростей в последовательные моменты времени и будет шумовой добавкой в средние скорости
AVer = к; - Ver.
3. Поэтому, вычитая значение шумовой добавки из скорости отдельной частицы в данной ячейке, получим скорректированное значение скорости частиц
Vr j Vr j /^v ver (rj ) .
4. Теперь, если подсчитаем среднюю скорость в начальный момент времени по формуле (16), получим Vjm с шумом.
5. Но тогда после прохождения шага т получаем VJT+1, но уже без шума.
Данный алгоритм реализуется на каждом временном шаге перед вычислением электрического поля Er по скоростям электронов. Проведенные расчеты показали, что этот алгоритм позволяет существенно уменьшить счетный шум (рис. 1). Кривая 2 соответствует пространственному распределению скорости, полученной с использованием описанного алгоритма.
Заключение
Приведем некоторые результаты тестовых расчетов, позволяющих судить о работоспособности созданного алгоритма для решения задачи распространения альфвеновской волны в бесстолкновительной плазме. На рис. 2 приведено пространственное распределение поперечной компоненты магнитного поля альфвеновской волны сдвига в последовательные моменты времени. Видно, что передний фронт волны распространяется с альфвеновской скоростью Va и имеет в соответствии с законом дисперсии [12] осцилляторную структуру. Характерный размер осцилляций S ~ c/opi. Заметим, что представленные моменты
Рис. 2.
времени соответствуют начальной стадии генерации альфвеновской волны. Исследование дальнейшей эволюции альвеновской волны требует проведения расчетов на существенно больших размерах расчетной области и не является целью данной статьи.
В настоящей статье описана новая гибридная модель, основанная на кинетическом описании электронной компоненты плазмы и гидродинамическом приближении для ионов, и представлены алгоритмы ее численной реализации. Модифицирован Р1С-метод для решения кинетического уравнения Власова, что позволило существенно уменьшить присущий Р1С-методу счетный шум. Разработаны алгоритмы корректировки положения частиц для выполнения условия квазинейтральности плазмы и вычисления электрического поля из уравнений движения отдельных частиц. Проведенная серия тестовых расчетов задачи распространения альфвеновской волны сдвига в неоднородном магнитном поле позволяет сделать вывод о возможности использования созданной модели при решении задач распространения альфвеновской волны в бесстолкновительной плазме. В дальнейшем на основе созданной модели предполагается изучение спектральных характеристик функции распре-
деления электронов и динамики прохождения альфвеновской волны в полярной области магнитосферы Земли.
В заключение авторы благодарят В.Т. Тихончука и В.А. Вшивкова за полезные обсуждения.
Список литературы
[1] Паркер Е. Динамические процессы в межпланетной среде. М.: Мир, 1965.
[2] Rezania V., Samson J.C. Quasi-periodic oscillations: Resonant shear Alfven waves in neutron star magnetospheres // Astronomy and Astrophysics. 2005. Vol. 436. P. 999-1008.
[3] Tikhonohuk V.T., Rankin R. Electron kinetic effects in standing shear Alfven waves in the dipolar magnetosphere // Phys. of Plasmas. 2000. Vol. 7, N 6. P. 2630-2645.
[4] VanZeeland M., Gekelman W., Vincena S., Maggs J. Current and shear Alfven wave radiation generated by exploding laser-produced plasma: Perpendicular incidence // Phys. of Plasmas. 2003. Vol. 10, N 5. P. 1243-1252.
[5] Louarn P., Wahlund J.E., Chust T., de Feraudy H. et.al. Observation of kinetic Alfven waves by the Freja spacecraft // Geophys. Res. Lett. 1994. Vol. 21, N 17. P. 1847-1850.
[6] Sakai J.I., Tsiklauri D., Yamamura W. et.al. Simulation of collisionless damping of shear wave — applications to coronal heating and solar wind acceleration // Proc. of ISSS-7, Kyoto, Japan, 2005.
[7] Chodura R. A hybrid fluid-particle model of ion heating in high-Mach-number shock waves // Nuclear Fusion. 1975. N 15. P. 55-61.
[8] Березин Ю.А., ДудниковА Г.И. Численные модели плазмы и процессы пересоединения. М.: Наука, 1985.
[9] Хокни Р., Иствуд Дж. Численное моделирование методом частиц. М.: Мир, 1987.
[10] Березин Ю.А., Вшивков В.А. Метод частиц в динамике разреженной плазмы. Новосибирск: Наука, 1980.
[11] Вшивков В.А., Снытников В.Н. О методе частиц для решения кинетического уравнения Власова // Журн. вычисл. матетатики и мат. физики. 1998. Т. 38, № 11. С. 1877-1883.
[12] Шафранов В.Д. Электромагнитные волны в плазме // Вопр. теории плазмы. 1963. Вып. 3. C. 3-140.
Поступила в редакцию 1 февраля 2006 г., в переработанном виде — 21 февраля 2006 г.
