Численное моделирование образования и эволюции токового слоя в ограниченной плазме Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 9, № 3, 2004

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ФОРМИРОВАНИЯ И ЭВОЛЮЦИИ ТОКОВОГО СЛОЯ В ОГРАНИЧЕННОЙ ПЛАЗМЕ*

Г. И. ДудниковА, В. П. Жуков Институт вычислительных технологий СО РАН, Новосибирск, Россия e-mail: [email protected], [email protected]

А. Г. ФРАНК

Институт общей физики им. А.М. Прохорова РАН, Москва, Россия

e-mail: [email protected]

We consider peculiarities of the evolution of a current sheet formed in a vicinity of the X-type null-line in the plasma, which is confined by the rigid walls. It is shown that such boundary conditions result in an abrupt and significant change of the current sheet magnetic structure and transformation of the current sheet into the elliptical Z-pinch.

Введение

Известно, что распространение магнитозвуковых возмущений в окрестности нулевой линии X-типа двумерного (2D) магнитного поля обычно завершается формированием токового слоя и значительной концентрацией магнитной энергии (см. [1-4] и приведенную там библиографию). Эта энергия может трансформироваться в тепловую и кинетическую энергию плазмы, в потоки излучений и ускоренных частиц в случае быстрого разрушения токового слоя, и тогда реализуются явления вспышечного типа, подобные вспышкам на Солнце [1-7]. Из процессов, которые могут инициировать разрушение слоя, в литературе обсуждаются тиринг-неустойчивость, тепловая неустойчивость, а также аномальное сопротивление плазмы [1, 6]. Как было установлено экспериментально, во многих случаях причинами разрушения слоя являются локальный сверхбыстрый нагрев плазмы и нарушение поперечного равновесия токового слоя, вызванное таким нагревом [8]. Указанные процессы обусловлены образованием в центральной области слоя магнитной структуры с замкнутыми силовыми линиями, т. е. локального Z-пинча, или магнитного острова, и подавлением теплопроводности плазмы вдоль поверхности слоя [8-11].

В настоящей работе на основе численного моделирования показано, что образование магнитного острова может быть связано при определенных условиях с ограниченным поперечным размером плазмы, который соизмерим с шириной токового слоя. Постановка задачи в настоящей работе близка к постановке в [12-15], где было проведено численное

* Работа выполнена при частичной поддержке Российского фонда фундаментальных исследований (грант № 03-02-17282) и Федеральной программы ведущих научных школ (НШ-1965.2003.2).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2004.

моделирование эволюции токовых слоев, развивающихся в окрестности нулевой линии X-типа 2Б магнитного поля под действием цилиндрической магнитозвуковой волны, приходящей из бесконечности. В этих работах моделировалось распространение возмущений в неограниченной области, что характерно главным образом для астрофизических приложений. В результате моделирования было показано, что образующийся токовый слой является устойчивым относительно конечных возмущений на протяжении интервалов времени порядка нескольких десятков альвеновских времен, т. е. слой находится в квазистационарном состоянии. В настоящей работе в отличие от [12-15] предполагается, что плазма занимает ограниченный объем, что характерно в основном для лабораторной плазмы. При определенных граничных условиях в такой ограниченной системе возможны быстрое перераспределение тока и существенное изменение магнитной структуры слоя, связанное с формированием в центральной области крупномасштабного магнитного острова. Отметим, что в принципе астрофизическая плазма также является ограниченной, однако по сравнению с лабораторными условиями ширина токового слоя обычно существенно меньше, чем расстояние до границы. Например, в качестве этого расстояния для солнечной короны в расчетах обычно принимается длина магнитной петли. Параметры, использованные в настоящей работе при численном моделировании, типичны для лабораторных экспериментов, которые проводятся в ИОФ РАН на установке "Токовый Слой" в рамках изучения процессов магнитного пересоединения, динамики и эволюции токовых слоев [10, 11].

1. Постановка задачи

Математическая постановка задачи такова. Плазма описывалась двумерными (д/дг = 0) уравнениями одножидкостной магнитной гидродинамики:

дЛ/дг + (УУ)Л = пДА; (1)

р(дУ/дг + (УУ)У) = -Ур - АЛУЛ + V(АУ + 1 Уа1уУ); (2)

др/дг + а1у(рУ) = 0; (3)

1 (др/дг + &у(Ур)) = -р &у(У) + ¿\у(хУТ) + п(АЛ)2 + (4)

Y - 1

2

Q = 2 ((dVx/dx)2 + (dVy/ду)2) - 3(divV)2+

+ (dVx/dy + dVy/дх)2 + (ÖVZ/дх)2 + (öVz/ду)2, T = p/p, H = (Hx, Hy) = (дА/ду, -дА/дх), V = (Vx,Vy), Y = 5/3.

Здесь А — z-компонента вектор-потенциала магнитного поля; остальные обозначения общепринятые. Использованы следующие безразмерные переменные. В качестве единицы длины выбран радиус камеры R, единицы плотности — начальная плотность p0, магнитного поля — магнитное поле на границе H0, скорости — скорость Альвена Va = Ho /^/4npö, времени — ta = R/Va, давления — магнитное давление H^/(4п), электрического поля H0Va/c. Безразмерные коэффициенты переноса имеют следующий вид: ц = с2/(4поRVa), v = v*/(p0RVa), X = X*Mi/(p0RVa). Здесь о — проводимость плазмы, v* — размерный

коэффициент динамической вязкости, х* — размерный коэффициент теплопроводности, Mi — масса иона.

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

A = A0 = (x2 - y2)/2 (X-точка), р = ß, р = 1, V = 0.

Расчетная область представляла собой квадрат —1 < x < 1, —1 < y < 1. Течение инициировалось возмущением вектор-потенциала A на границе этой области, имевшем вид

A = Ao — Eit2/(t + Д), (5)

что соответствует электрическому полю Ez = —dA/dt, растущему за характерное время Д (ниже везде Д = 0.5) от нуля до Ei.

Границы расчетной области соответствовали жесткой стенке с заданной температурой, т. е.

V = 0, р/р = ß. (6)

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

Поставленную задачу достаточно решать в 1/4 области, полагая оси координат осями симметрии.

2. Конечно-разностная схема

Для численного решения поставленной задачи использовалась следующая явная конечно-разностная схема.

Функции А и V вычислялись в целых, а р и р — в полуцелых узлах сетки. Это позволяло более компактно аппроксимировать производные по пространству и таким образом несколько улучшить точность расчетов и более естественно реализовать граничные условия в конечно-разностном случае. Целые узлы сетки имели координаты хг = гЛ,х, Уз = , а полуцелые узлы — координаты £¿+1/2 = (г + 1/2)Л,х, уз+1/2 = (З + 1/2)^у. Здесь г = —1, 0,1,...,г0, З = —1, 0,1,...,Зо, г0 Ьх = З0Л,у = 1. Функции Агз- Vij при г = —1 или З = —1 вычислялись с учетом симметрии задачи по значениям этих функций в узлах г =1 или З = 1. Функции рг+1/2)3+1/2 и рг+1/2)3+1/2 при г = —1 или З = —1 вычислялись с учетом симметрии задачи по значениям этих функций в узлах г = 0 или З = 0. Постановка граничных условий (5), (6) в конечно-разностном случае также не вызывает труда.

Во внутренних точках для перехода от момента времени ¿га к моменту ¿""+1 = ¿га + т использовались следующие формулы:

n

An+1 _ A

-j = LA(Vn, An),

Vn+1 _ V'n ij = LV(Vn, An+\pn, pn),

т

pn+1 _ nn Pi+1/2,j+1/2 pi+1/2j+1/2 _ fc, +1

т

Lh(Vn+1, pn),

р+1/2,]+1/2 Р'г+1/2,3+1/2 = ^/1(уп+1 Лп+1 рп рп+1)

т V ' ' '

Здесь ¿А, ¿V, Ь^, Ь^ — разностные аналоги операторов Ьд, ЬV, Ьр, Ьр для уравнений (1)-(4), записанных в виде дЛ/дг = ЬА(У, Л), дУ/дг = ЬV(V, Л, р, р) и т.п.

Оператор Лапласа аппроксимировался на пятиточечном шаблоне обычным образом. Градиенты от различных функций вычислялись с помощью центральных разностей, использующих значение функции в ближайших четырех узлах. Исключение составляли члены вида (УУ)/, где / = (Л, V), и ^у(/У), где / = (р,р). В первом случае выражение Ухд//дх и аналогично Ууд//ду аппроксимировались разностью вперед или назад в зависимости от знака скорости. Во втором случае &у(/V) аппроксимировалось выражением

!А' ({ЛТ\\ Ях г +1,3+1/2 — Яхг,3+1/2 Яу г+1/2,3+1 — Яу г+1/2,3

(Шу(/У ))г +1/2,3+1/2 = -7--1--7-•

Пх Пу

Здесь Ях = У*/*, где У* = (Ухг,з + УХг,з+1)/2, а /* = /-1/2,3, если У* > 0, или /+3/2,3, если У * < 0. Для Яу имеем аналогичные выражения.

Приведенная схема имеет первый порядок аппроксимации по т, Т х и Т у и является условно устойчивой. Ограничения на шаг т имеют вид типа т < (|Ух| + с)Пх, (\Уу| + с)Пу, тах(п, V, х) тах(Пх, Пу). Здесь с — максимальная скорость звука (альвеновская скорость).

Точность расчетов контролировалась решением задачи на сетках с различными шагами Пх и Пу. Большинство расчетов проводилось на сетке 80 х 80, что обеспечивало погрешность не более 5 %.

3. Результаты расчетов

Рассмотрим результаты расчетов для варианта Е1 = 0.16, п = 0.015, х = V = 10-4 и в = 0.05. Такие значения Е1 и V соответствуют, например, следующим размерным параметрам, характерным для установки "Токовый Слой" [10, 11]: плотность электронов Ме = 1015см-3, Н0 = 10 кГс, радиус установки 5 см, плазма создавалась в аргоне (Мг/Мр = 40), степень ионизации Z = 2 и проводимость а = 0.6 • 1014 с-1. Электрическое поле равно 250 В/см.

В начальные моменты времени получена стандартная картина течения: к центру координат сходится магнитозвуковая волна. Примерно к моменту г =1.5 она достигает центра координат и на оси х формируется токовый слой с плотностью тока (рис. 1-3, г = 4.7). Плазма втекает в этот слой поперек и вытекает вдоль него. Это течение мало меняется в течение нескольких альвеновских времен, т. е. является метастабильным. Оно мало отличается от течения, наблюдаемого в случае неограниченной плазмы.

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

В токовом слое происходит процесс пересоединения, благодаря которому на оси х постоянно генерируется у-компонента магнитного поля из х-компоненты. Изменение магнитного потока Л на границе (5) таково, что оно стремится уменьшить этот поток (дЛ/дг < 0). В силу граничных условий непротекания (6) поток магнитного поля не может выноситься за пределы расчетной области. Поэтому уменьшение Л на границе означает уменьшение Ну в окрестности точек х = ±1, у = 0. Более того, в окрестности этих точек происходит генерация Ну обратного знака. Этот процесс проиллюстрирован на рис. 4, где изображена

функция А(х, у = 0) — А(х = 0, у = 0) = — / Ну(£, у = в различные моменты времени.

о

Заметим, что в случае неограниченной плазмы скорость плазмы и изменение А согласованы таким образом, что уменьшение А на границе равно уменьшению магнитного потока вследствие его выноса из расчетной области. Соответственно, поля Ну обратного знака не возникает.

Описанную выше тенденцию к образованию в окрестности точек х = ±1, у = 0 поля Ну противоположного знака можно также проиллюстрировать на примере стационарного (дА/д£ = —Е1) уравнения для Ну = —дА/дх на оси х. Оно имеет вид

пдНу/дх = Ух Ну + / (х), (7)

где / = Е + пд2А/ду2, Ну(х = 0) = 0.

Внутри токового слоя / < 0, что обеспечивает Ну < 0. Вблизи границы х = 1 в силу (5) д2А/ду2 = д2А0/ду2 = —1. Поскольку интерес представляет случай Е1 ^ п, то /(х ^ 1) ~ Е1 > 0. Скорость плазмы при х ^ 1 стремится к нулю. Рассмотрим окрестность 8* точки х = 1, у = 0, в которой УцНу ^ Е1. Согласно (7) в этой окрестности Ну увеличивается на 8*Е1/п. При достаточно большом Е1/п имеем |Ну(х = 1 — 8*)| < 8*Е1/п, т.е. знак Ну вблизи границы изменяется. В случае неограниченной плазмы УХ (х = 1) > 0 достаточно велико, поэтому правая часть (7) отрицательна или по крайней мере положительна, но не очень велика при всех х. Соответственно, изменения знака Ну не происходит.

Появление Ну противоположного знака связано с увеличением тока в окрестности х = ±1, у = 0 и образованием О-точек магнитного поля, которые со временем становятся все более ярко выраженными (рис. 3, £ = 5.1). В результате токовый слой оказывается окружен замкнутыми силовыми линиями. Сила натяжения такого магнитного поля стремится сжать токовый слой вдоль оси х к центру.

Рассмотрим баланс магнитных сил в окрестности точек х = ±1, у = 0 до образования магнитного поля обратного знака (£ < 5) и после (£ > 5). Поскольку V около границы стремится к нулю, в силу (1) и (5) З^(х = ±1) = Е1/п. Соответственно, сила ^Х(х = ±1, у = 0) = — , действующая на плазму со стороны магнитного поля в окрестности этих точек при £ < 5, направлена к границам расчетной области. При £ > 5 из-за изменения знака Ну (х = ±1, у = 0) эта сила направлена к центру области.

В связи с изменением направления действия сил при £ > 6.5 возникает волна, которая распространяется к центру расчетной области (см. рис. 1-3, 6.5 < £ < 8.0), и в результате происходит быстрое увеличение размера токового слоя по оси у. Фронт волны движется со скоростью порядка половины альвеновской скорости. Скорость плазмы за фронтом направлена от границы х =1 к центру координат. На рис. 5 приведена зависимость З^ в центре от времени: из рис. 1 и 5 видно, что плотность тока за фронтом волны уменьшается. В момент достижения волной центра координат ток может принимать отрицательные значения. Отметим, что на рассматриваемых этапах эволюции полный ток, протекающий по всей расчетной области, со временем увеличивается.

После достижения волной центра координат и некоторого переходного процесса возникает магнитная конфигурация, которая в дальнейшем не претерпевает принципиальных изменений: ток (см. рис. 1, £ = 25) и плазма сосредоточены в центре и занимают значительную часть расчетной области. В результате происходит полная перестройка магнитного поля. Подчеркнем, что начальное магнитное поле и поле на этапе существования токового слоя топологически эквивалентны Х-точке, а на больших временах возникает О-точка. При этом Ну на оси х на больших временах изменяет знак по сравнению с начальным

х

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

Рис. 2. Поле скорости (в 1/4 области) в различные моменты времени.

Рис. 3. Силовые линии (изолинии А) в различные моменты времени.

Рис. 4. Распределение А(х,у = 0) — А(х = 0, у = 0) в различные моменты времени.

Рис. 5. Зависимость в центре области (х у = 0) от времени.

(см. рис. 1 и 4, £ = 25). Конфигурация магнитного поля, устанавливающаяся на больших временах, близка к конфигурации, получающейся при решении задачи о распределении магнитного поля в случае неподвижной плазмы с граничными условиями (5).

Остановимся на количественных закономерностях. В определенном диапазоне параметров с погрешностью 10% справедливы следующие соотношения (в = 0.05, х = 0):

и = 2.06£Г°-94п0Л2 ~ 1/Е, Зтах = 1.06£Гп"0-63,

5 = 1.47ЕО'05П0'58 - /п, А = 1.2Е0'4. (8)

Здесь — момент прихода волны от границ плазмы к центру координат, а Зтах — максимальное значение (х = у = 0) (см. рис. 5); 5 — толщина слоя, определяемая как расстояние вдоль оси у, на котором значение равно 1/2 от его значения в центре области; А — ширина токового слоя вдоль оси х, определяемая аналогично его ширине. Величина А во время метастабильной фазы (2 < £ < 6) практически постоянна. Отметим слабую зависимость рассматриваемых величин от в.

Остановимся на поведении полной массы плазмы. При условии гибели плазмы на стенках камеры эта величина на этапах существования и разрушения токового слоя уменьшается со временем. На больших же временах (£ > 25) она постоянна. Это связано с тем, что плазма оказывается сосредоточенной в центре области и поток плазмы на стенки мал.

При некотором соотношении между Е1 и п возможен сценарий развития течения, при котором плазма, имеющаяся в расчетной области, сосредоточивается в токовом слое, а затем выбрасывается из него в окрестность точек х = ±1, у = 0, где гибнет. Отражение не успевает происходить. Характеристикой этого процесса может служить время при котором плотность плазмы в центре слоя становится меньше 0.1. При не слишком маленьком > 5.5) с погрешностью 10% справедливо соотношение

1Р = 1.8Е-0'43п-0'1 5 - Е-1 /2.

Заключение

В настоящей работе на основе численного моделирования исследовались особенности поведения ограниченной плазмы в 2D магнитном поле c нулевой линией X-типа. Показано, что в условиях гибели плазмы на границе, т. е. когда плотность плазмы на границе не может превышать начальную, на определенном этапе эволюции слоя происходят быстрое перераспределение тока и существенное изменение магнитной структуры. Этот процесс начинается с увеличения плотности тока вблизи боковых концов слоя, так что слой оказывается окруженным замкнутыми магнитными силовыми линиями. В результате направление сил, воздействующих на плазму, изменяется на противоположное, что приводит к сжатию плазмы по оси x (вдоль большего размера слоя) при одновременном существенном увеличении его меньшего размера по оси y. Перестройка магнитной структуры слоя приводит к формированию в центральной области крупномасштабного магнитного острова, где сосредоточены большая часть плазмы и тока, и такая конфигурация в дальнейшем не претерпевает принципиальных изменений в течение по крайней мере 20 альвеновских времен.

Список литературы

[1] SyrüVATSKII S.I. Pinch sheets and reconnection in astrophysics // Annu. Rev. Astron. Astrophys. 1981. Vol. 19. P. 1б3-229.

[2] BISCAMP D. Nonlinear Magnetohydrodynamics. Cambridge Univ. Press, 1993.

[3] Priest E.R., FüRBES T. Magnetic reconnection. MHD theory and applications. Cambridge Univ. Press, 2000.

[4] КАДОМЦЕВ Б.Б. Перезамыкание магнитных силовых линий // Успехи физ. наук. 1987. Т. 151. С. 3-29.

[5] SVESTKA Z. Solar flares. Dordrecht: Reidel, 197б.

[6] СыРОВАТСКИй С.И. Ключевые вопросы теории вспышек // Изв. АН СССР. Сер. физ. 1979. Т. 43. С. б95-707.

[7] ФРАНК А.Г. Формирование, эволюция и взрывное разрушение токовых слоев в плазме // Вопросы физики плазмы и плазменной электроники: Тр. ФИАН. М.: Наука, 1985. Т. 1б0. С. 93-121.

[8] КИРИй Н.П., МАРКОВ B.C., ФРАНК А.Г. Локальный импульсный нагрев плазмы и разрушение токового слоя // Письма в ЖЭТФ. 1992. Т. 5б(2). С. 82-8б.

[9] БОГДАНОВ С.Ю., МАРКОВ B.C., ФРАНК А.Г. Изменение топологии магнитного поля в процессе взрывного разрушения токового слоя // Письма в ЖЭТФ. 1982. Т. 35. С. 232-235.

[10] БОГДАНОВ С.Ю., ДРЕйДЕН Г.В., КИРИй Н.П. И ДР. Динамика плазмы в токовых слоях. II. Нелинейные режимы формирования токовых слоев // Физика плазмы. 1992. Т. 18. С. 1283-1295.

[11] Богданов С.Ю., Кирий Н.П., Франк А.Г. Эволюция двумерных токовых слоев в линейных и нелинейных режимах // Магнитное пересоединение в двумерных и трехмерных конфигурациях: Тр. ИОФАН. М.: Наука, 1996. Т. 51, вып. 2. С. 5-76.

[12] Брушлинский К.В., Заборов А.М., СыровАтский С.И. Численный анализ токового слоя в окрестности магнитной нулевой линии // Физика плазмы. 1980. Т. 6, вып. 2. C. 297-311.

[13] Буланов С.В., ДудниковА Г.И., Жуков В.П. и др. Токовые слои в окрестности критических точек магнитного поля / / Магнитное пересоединение в двумерных и трехмерных конфигурациях: Тр. ИОФАН. М.: Наука, 1996. Т. 51, вып. 2. C. 101-123.

[14] Буланов С.В., ДудниковА Г.И., Жуков В.П. и др. Численное моделирование токового слоя в окрестности нулевой линии магнитного поля / / Краткие сообщения по физике. 1994. № 5, 6. C. 28-33.

[15] Буланов С.В., ДудниковА Г.И., Жуков В.П. и др. Пересоединение магнитных силовых линий в окрестности критических точек // Физика плазмы. 1996. Т. 22, № 10. C. 867-895.

[16] Березин Ю.А., ДудниковА Г.И. Численные модели плазмы и процессы пересоединения. М.: Наука, 1985.

Поступила в редакцию 10 декабря 2003 г., в переработанном виде — 27 марта 2004 г.