Сравнение аналитической модели бесстолкновительного магнитного пересоединения, построенной на базе решения уравнения Грэда-Шафранова, с кинетической, рассчитанной методом particle-in-cell Текст научной статьи по специальности «Физика»

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

Сер. 4. 2009. Вып. 1

УДК 537.84, 533.9.072

Д. Б. Коровинский, А. В. Дивин, В. С. Семёнов СРАВНЕНИЕ АНАЛИТИЧЕСКОЙ МОДЕЛИ

БЕССТОЛКНОВИТЕЛЬНОГО МАГНИТНОГО ПЕРЕСОЕДИНЕНИЯ, ПОСТРОЕННОЙ НА БАЗЕ РЕШЕНИЯ УРАВНЕНИЯ ГРЭДА-ШАФРАНОВА, С КИНЕТИЧЕСКОЙ, РАССЧИТАННОЙ МЕТОДОМ PARTICLE-IN-CELL *>

Введение. Согласно результатам исследований, проведённых в рамках проекта Geospace Environment Modelling (GEM) Magnetic Reconnection Challenge [1], учёт вклада эффекта Холла является минимальным необходимым требованием к моделям, описывающим быстрое магнитное пересоединение в бесстолкновительной плазме. Действие этого эффекта создает характерную квадрупольную структуру магнитного поля, предсказанную Соннерапом [2], и обеспечивает высокую скорость пересоединения (порядка 0,2) вне зависимости от физического механизма процессов, протекающих в ближайшей окрестности X-линии. Эта окрестность, называемая электронной диффузионной областью (EDR), характеризуется повышенной ролью эффектов, связанных с инерцией электронов, а также с анизотропией тензора электронного давления, что приводит к непригодности жидкостных приближений внутри EDR и необходимости использовать кинетическое описание плазмы. Существенно, однако, что эта область весьма сжата; в направлении поперёк слоя её размер характеризуется инерционной длиной электрона le = c/rae, где с есть скорость света, Юе = (4nne2/me)1/2 - плазменная частота электрона, а me - его масса. В направлении ускорения плазмы EDR демонстрирует составную структуру с размерами порядка lp и 10lp, где lp = с/(4кпе2/mp)1/2 - инерционная длина протона, а mp - его масса. На масштабе lp наблюдается сильный электронный ток в направлении X-линии, а на масштабе 10lp - электронные джеты в обе стороны от неё [3].

Упрощающие предположения. Другим существенным фактом является то, что вблизи X-линии, на масштабах порядка lp, протоны из-за инертности движутся много медленнее электронов, что позволяет считать полный ток примерно равным электронному [4]

j = -neVe, (1)

где j обозначает плотность тока, e - абсолютную величину заряда электрона, n - концентрацию электронов, а Ve - их массовую скорость. Магнитогидродинамическое описание плазмы при таком допущении и с учётом эффекта Холла составляют приближение электронной холловской магнитогидродинамики (eHMHD). Для рассматриваемого случая стационарной квазинейтральной бесстолкновительной плазмы закон Ома в приближении eHMHD имеет вид

Е+ -Ve х В = — — VPe, (2)

с ne

Работа выполнена при поддержке Российского фонда фундаментальных исследований, грант № 07-05-00776-а. PIC-моделирование было проведено А. В. Дивиным во время стажировки в университете Мэриленда, США.

© Д. Б. Коровинский, А. В. Дивин, В. С. Семёнов, 2009

где Е и В есть векторы напряжённости электрического и индукции магнитного полей, соответственно, а Ре - тензор электронного давления. Из выражения (1) и уравнения (2) следует, что в уравнении движения протонов

е

пшр (Vр ■ \7)Л^р — —Рр + еЕ И—Vр х В,

(3)

где Ур есть массовая скорость протонов, а Рр - тензор протонного давления, можно пренебречь третьим слагаемым в правой части, по сравнению со вторым. Согласно данным численного моделирования [1, 3], дивергенция тензора электронного давления в уравнении (2) вносит весомый вклад только внутри ЕБИ, где этот тензор становится существенно анизотропным; вне ЕБИ электронное давление можно считать скалярной величиной, чем мы и воспользовались в нашем анализе. То же справедливо и в отношении тензора протонного давления в уравнении (3). Остальными уравнениями еНМНБ являются уравнения Максвелла и уравнение неразрывности.

Аналитическая модель. Несущественность физического механизма процессов внутри ЕБИ для общей картины и описание плазмы в еНМНБ-приближении позволяют получить самосогласованное аналитическое решение для задачи стационарного бесстолкновительного магнитного пересоединения в квазинейтральной несжимаемой плазме. Здесь мы приводим схему решения, обсуждаемого детально в работе [5]. Выберем оси координат следующим образом: за направление оси X примем направление вектора индукции магнитного поля на бесконечности (в верхней полуплоскости), за ось У примем направление электрического тока вдоль X-линии, и ось Z направим перпендикулярно им обеим. Будем считать слой однородным вдоль X-линии, т. е. положим, что для всех плазменных величин производная д/ду = 0. Для плоскости XZ и векторов, принадлежащих этой плоскости, введём обозначение ^. Далее, приведём величины к безразмерному виду посредством нормировки индукции магнитного поля В' = В/Во, протонных и электронных скоростей Ур е = V' ре/Уа, напряжённости электрического поля Е' = Е/Еа, газовых давлений Рр е = Рре/Ро, и расстояний г' = г/1р. Здесь Во - значение индукции магнитного поля на верхней границе рассматриваемой области над X-линией, Уа = Во/(4ппшр)1/2 - соответствующая протонная альфве-новская скорость, Еа = (\/с)УаВ0 - альфвеновское электрическое поле, Р0 = В2/4п, и г = (х,у, г). Ниже мы будем использовать нормированные величины, опуская штрихи. Наконец, введём потенциал магнитного поля А и потенциал электрического поля ф:

где еу есть орт оси У. Заметим, что из однородности слоя в направлении оси У и закона Фарадея в стационарном случае следует, что компонента электрического поля Еу, играющая роль скорости пересоединения магнитного потока, равна константе, которую мы обозначаем е. Компонента магнитного поля Ву является также потенциалом - функцией тока для скорости электронов в плоскости XZ [4],

С использованием введённых обозначений решение задачи описывается [5] системой уравнений (5)-(10):

В± = (Вх, Вг) = Ух (Аеу);

Е± = (Ех ,Ег) = -У±ф(х,г),

Уе± = (Уех ,Уег) = -Ух (Ву )еу.

(4)

(5)

r

ву(г) = (- i)k+h J ^L-^+ВуЫ; (6)

ro

Ф eft = т^Ву + G(A)] (7)

1-V^+Il-1-\V±A\2 + G(A)=CtT] (8)

(V± ■ Vp±)=0; (9)

r

VPy(r)=eJ^- + Vpy(r0)- (10)

ro

Здесь магнитный потенциал A и компонента скорости электронов Vey находятся из решения уравнения Грэда-Шафранова (7), где G(A) - неизвестная моделирующая функция, о которой будет сказано ниже. Компонента магнитного поля By находится интегрированием вдоль проекций силовых линий на плоскость XZ согласно уравнению (8), в котором ro обозначает начальную точку интегрирования, а k - номер квадранта. Эффективный электрический потенциал фед = ф — (1/ne)Pe находится из уравнения (9). Движение протонов в плоскости XZ подчиняется закону Бернулли (10), где

П = Pp + (1/2)B2 есть полное давление, а Ctr обозначает константу вдоль траекто-

рии; также их движение подчинено условию неразрывности (11). Наконец, компонента массовой скорости протонов Vpy находится интегрированием вдоль проекций их траекторий на плоскость XZ согласно уравнению (12).

Таким образом, ключом к построенной модели является неизвестная функция G(A), стоящая в правой части уравнения Грэда-Шафранова (7). Заметим, что эта функция имеет ясный физический смысл, согласно уравнению (9) она даёт основной вклад в эффективный электрический потенциал, а её производная есть электрический ток в направлении X-линии. Это обстоятельство позволяет на основе простых физических соображений задать функцию dG/dA и получить решение задачи. Пример такого решения и его анализ приведены в [5], здесь же мы только суммируем основные черты этой модели.

Решение построено без учёта диссипации, поэтому, строго говоря, EDR не входит в область его применимости. В частности, это видно из того, что в начале координат знаменатель подынтегрального выражения в уравнении (8) обращается в ноль, а величина By, соответственно, в бесконечность. Поэтому при вычислении By по формуле (8) интегрирование должно обрываться на границе EDR, а её вклад должен учитываться отдельно. В то же время, в начале координат значение By равно нулю по условию симметрии, поэтому, учитывая гладкость функции By (x,z) и малые размеры EDR, вклад EDR можно считать пренебрежимо малым. Согласно той же формуле (8), экстремальные значения функция By (x, z) должна принимать на сепаратрисах магнитного поля, а в ближайшей их окрестности достигает экстремальных значений градиент By, т. е., согласно (6), компоненты электронной скорости Vex,ez. Кроме того, на сепаратрисах принимает экстремальное значение градиент функции G[A(x, z)], а также функция Vey (x, z), для которой закон Ампера позволяет получить следующую оценку: max\Vey\ = VA(mp/me)1/2/5, где 5 есть полуширина (размер вдоль оси Z) EDR, измеренная в инерционных длинах электрона le. В частности, при 5 =1, max \ Vey \ достигает электронной альфвеновской скорости. Из вышесказанного следует, что вблизи сепаратрис магнитного поля локализованы сильные электрические поля с напряжённостью

порядка 10Ea. Заметим также, что вне EDR вклад градиента Pe пренебрежимо мал, поэтому Е ~— grad фе^.

Уравнениями (7)-(9) полностью описывается холловская структура магнитного поля и соответствующих токов. Движение протонов находится из подсистемы (10)—(12), для решения которой необходимо привлечь информацию о распределении полного давления П. Конфигурация задачи позволяет использовать приближение погранслоя, поэтому полное давление в рассматриваемой области является функцией лишь одной переменной, П = n(x), и исполняет роль граничного условия, которое должно определяться из решения для объемлющего пространства. При простейшем предположении П = const формула (10) описывает движение протонов с ускорением вплоть до альф-веновской скорости [5]. Приближение погранслоя позволяет записать решение уравнения (7) в виде:

где функция Вг(х, 0) в уравнении (14) является граничным условием.

Резюмируя, система уравнений (7-14) описывает самосогласованное решение задачи стационарного бесстолкновительного пересоединения в квазинейтральной несжимаемой плазме, получаемое на основе решения уравнения Грэда-Шафранова для магнитного потенциала в приближении погранслоя при граничных условиях, задаваемых функциями Вг (х, 0) и П(х, гтах).

Для проверки эффективности представленной модели мы воспользовались результатами численного моделирования, проведённого методом Рагйс1е-ш-Се11 (Р1С^1ти1а-^оп), реализованном в коде Р3Б [6], который является кинетическим электромагнитным кодом для расчёта динамики протонов и электронов методом частиц. Код адаптирован для расчётов на многопроцессорной вычислительной машине с использованием библиотеки параллельных вычислений МР1.

Р1С-моделирование. Начальное состояние задавалось в виде нейтрального токового слоя типа слоя Харриса [7] с добавлением фоновой плазмы однородной плотности. Для возбуждения процесса пересоединения в центре вычислительной области на исходную конфигурацию магнитного поля Вх = Во tanh(z/^) накладывалось возмущение типа X-точки:

где параметры (Ъх,Ъг), определяет размер вычислительной области. Все физические величины были нормированы на их соответствующие значения в начальном невозмущённом состоянии, пространственный масштаб - на инерционную длину протона, а время - на обратную гирочастоту И-1 = (трс)/(еВо). При заданной интенсивности возмущения Фо = 0, 3 за время £ « 15 процесс входил в квазистационарный режим, параметры которого мы сравнили с аналитической моделью. Отношение масс составляло тр/те = 64, отношение температур Тр/Те = 3/2. Шаг по времени для электромагнитного поля осуществлялся методом с перешагиванием, шаг по времени для частиц

A

(11)

x

(12)

o

(13)

0,8 0,7 0,6 0,5 а 0,4 0,3 0,2 0,1 0

4

2

0

2

4

Рис. 1. Данные Р1С-моделирования: а) ]еу (0, г) - жирная линия, ]ру (0, г) - тонкая линия; б) \Уеу (х, 0)| - жирная линия, Уру (х, 0) - тонкая линия

рассчитывался методом Бориса [8]. Вычисления проводились в квадрате со стороной Ьх = Ь = Ь = 38,4 при открытых граничных условиях [9, 10] для частиц:

dne

0.

д Ve

дТр

дх ’ дх и электрического, и магнитного полей:

dt

дВ

*,у

дх

0,

дЕ.

дх

= 0. Ex.z = 0. Bz = 0.

(14)

(15)

На верхней и нижней границах z = ±19. 2 были поставлены граничные условия типа проводящей поверхности (perfect electric conductor). При таких граничных условиях за рассматриваемое вычислительное время из системы утекает не более 15 % всего магнитного потока и частиц.

Сравнение моделей. На рис. 1а представлено полученное распределение токов в направлении X-линии; отчётливо выделяется EDR как область преобладания электронного тока над протонным; её полуширина 8 « (3/4)lp, что при данном соотношении масс соответствует примерно 6le. Согласно приведённой оценке, при таком размере EDR величина max \Vey \ должна достигать значения 7VA; результат PIC-моделирования даёт 6VA (рис. 1б). Видно, что отношение электронного тока к протонному монотонно уменьшается при удалении от X-линии от примерно 11 в начале координат до 5 при х = 4. Поскольку условие применимости eHMHD формулируется как \Ve\ ^ \Vp\, мы приняли значение xmax =4 за границу моделируемой области. Из тех же соображений верхняя граница области eHMHD по данным PIC-моделирования определена нами как прямая zmax = 4, что соответствует примерно 30le при mp/me = 64, т. е. примерно 0. 7lp при реальном соотношении масс. Что касается продольного размера EDR, то её полудлина составила приблизительно 2lp. В качестве исходных данных для аналитической модели были также привлечены полученные в PIC-моделировании функции Vey (А) и Bz(х. 0) (рис. 2), а также данные по полному давлению П(х. zmax), которое

- 1

- 2

а - 3

- 4

- 5

- 6

- 1 - 0,5 0 0 1 2 3 4

Рис. 2. Данные Р1С-моделирования - а) производная йО/йЛ = Уеу(Л); б) магнитное поле Бг (х, 0):

точка А = 0 на рис. 2а соответствует сепаратрисам магнитного поля; тонкой линией показаны результаты моделирования, жирной линией — аппроксимирующие кривые

явилось линейно нарастающей, но мало меняющейся величиной: П(0,гтах) = 0, 61, П(жтах,2:тах) = 0, 66. Наконец, полученное в кинетическом моделировании значение скорости пересоединения е = Еу составило 0,2.

На рис. 3 представлены траектории электронов по результатам обеих моделей. В обоих случаях просматривается сепаратриса магнитного поля, вдоль которой выстилаются траектории, что являет собой классическую структуру холловских токов [2], зафиксированную лабораторными экспериментами [11] и спутниковыми наблюдениями [12]. Также видны электронные джеты вдоль оси X в согласии с результатами аналогичных исследований, проведённых другими авторами [3, 13]. Профили скоростей в этих джетах представлены на рис. 4. При этом значения электронной скорости, полученные из теоретической модели, оказываются занижены примерно вдвое по сравнению с данными Р1С-моделирования, что объясняется, очевидно, действием дивергенции тензора электронного давления, не учитываемого в еНМНБ-приближе-нии. Как видно из рисунка, в обеих моделях происходит монотонное увеличение скорости протонов до одинаковых, примерно, скоростей; упрощения, принятые в еНМНБ-модели, приводят к сокращению пространственных масштабов этого процесса, примерно втрое. Несмотря на количественное расхождение, качественно картины движения протонов, представленные на рис. 5, демонстрируют хорошее согласие. Также хорошо согласуются структуры магнитного поля, представленные на том же рисунке.

Резюме. Графики, представленные на рис. 3 и 5 показывают, что допущенное пренебрежение тензором электронного давления приводит к искажению (сжатию) масштабов процесса и снижению точности результата, особенно вблизи оси X - области

Рис. 3. Траектории электронов в плоскости XZ по результатам аналитической модели (а) и Р1С-моделирования (б)

2,5

1,5

о, 5

б

о

1

3

4

Рис. 4- Скорости электронов Vex (x, 0) - жирная линия и протонов Vpx (x, 0) - тонкая линия по данным аналитической модели (а), по данным PIC-моделирования (б):

внутри EDR, т. е. при x < 2, электронная скорость в аналитической модели не определена

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

2

1

о

Рис. 5. Траектории протонов в плоскости XZ - жирные линии и магнитные силовые линии - тонкие линии, по данным аналитической модели (а) и по данным Р1С-моделирования (б)

2

3

4

сепаратрисы магнитного поля, и соответствующей квадрупольной структуры этого поля; формирование узкого электронного джета вдоль X-линии, джетов в направлении ускорения протонов и само их ускорение.

Литература

1. Birn J., Drake J. F., Shay M. A., Rogers B. N., Denton R. E. et al. Geospace Environmental Modeling (GEM) Magnetic Reconnection Challenge // J. Geophys. Res. 2001. Vol. 106. N A3. P. 3715-3719.

2. Sonnerup B. U. O. Magnetic field reconnection // Solar System Plasma Physics. 1979. Vol. 3. P. 45-108.

3. Shay M.A., Drake J.F., Swisdak M. Two-scale structure of the electron dissipation region during collisionless reconnection // Phys. Rev. Lett. 2007. Vol. 99. P. 155002-(1)-155002-(4).

4. Biskamp D. Magnetic reconnection in Plasmas. Cambridge: University Press, 2000. P. 207-209.

5. Korovinskiy D. B., Semenov V. S., Erkaev N. V., Divin A. V., Biernat H. K. 2.5D analytical model of steady-state Hall magnetic reconnection // J. Geophys. Res. 2008. Vol. 113. N A0. P. 4205-(1)-4205-(13).

6. Zeiler A., Biskamp D., Drake J. F., Rogers B. N., Shay M. A. et al. Three-dimensional particle simulations of collisionless magnetic reconnection // Ibid. 2002. Vol. 107. N A9. P. 1230-1236.

7. Harris E. G. On a plasma sheet separating regions of oppositely directed magnetic field // Nu-ovo Cimento. 1962. Vol. 23. P. 115-121.

8. Birdsall C. K., Langdon A. B. Plasma Physics via Computer Simulation. Bristol: Adam Hilger, 1991. P. 123-145.

9. Divin A. V., Sitnov M.I., Swisdak M., Drake J. F. Reconnection onset in the magneto-tail: Particle simulations with open boundary conditions // Geophys. Res. Lett. 2007. Vol. 34. P. L09109-(1)-L09109-(14).

10. Pritchett P. L. Collisionless magnetic reconnection in a three-dimensional open system // J. Geophys. Res. 2001. Vol. 106. N A11. P. 25961-25977.

11. Cothran C. D., Landreman M., Brown M.R., Matthaeus W.H. Generalized Ohm’s law in a 3D reconnection experiment // Geophys. Res. Lett. 2005. Vol. 32. P. L03105-(1)-L03105-(4).

12. Alexeev I. V., Owen C. J., Fazakerley A. N., Runov A., Dewhurst J. P. et al. Cluster observations of currents in the plasma sheet during reconnection // Ibid. 2005. Vol. 32. P. L03101-(1)-L03101-(4).

13. Daughton W., Scudder J. and Karimabadi H. Fully kinetic simulations of undriven magnetic reconnection with open boundary conditions // Phys. Plasmas. 2006. Vol. 13. P. 072101-(1)-072101-(15).

npHHSTO k пу6aнкацнн 2 ceHTa6pa 2008 r.