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

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

Том 12, № 2, 2007

ТРЕХМЕРНАЯ МОДЕЛЬ ТИРИНГ-НЕУСТОЙЧИВОСТИ В ОТКРЫТЫХ ЛОВУШКАХ С ЭЛЕКТРОННЫМ ПУЧКОМ*

А.В. Бурдлков Институт ядерной физики СО РАН, Новосибирск, Россия e-mail: [email protected]

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

Pioneering results on the numerical modeling of the tearing instability in an open trap with the electron beam (GOL-3, Budker Institute of Nuclear Physics SB RAS) are presented. The reduced magnetohydrodynamics approach in 3-dimensional cylindrical geometry had been used. To solve the problem a finite-difference code has been created.

Введение

Хорошо известная тиринг-неустойчивость изучалась многими исследователями с использованием разнообразных физических моделей (кинетических, магнитогидродинамических) в различных геометриях (плоской, цилиндрической, двух- и трехмерной) [1-11]. Постановка исследуемых ранее задач была связана с процессами, происходящими в космической плазме либо в установках типа токамак, компактный тор, "Токовый слой" и т. д. В последнее время выяснилось, что тиринг-неустойчивость может возникать также в открытых ловушках с электронным пучком. Более того, эта неустойчивость играет существенную роль в динамике плазмы в таких установках. В частности, по-видимому, она ответственна за внезапное разрушение токовой структуры, сопровождающееся иногда значительным разрушением внутренней поверхности стенок камеры установки ГОЛ-3 в Институте ядерной физики СО РАН [12].

Подчеркнем, что тиринг-неустойчивость на установке ГОЛ-3 развивается в условиях, имеющих принципиальные отличия от изучавшихся ранее случаев. Эти отличия состоят в следующем.

1. На установке ГОЛ-3 необходимая для развития тиринг-неустойчивости конфигурация магнитного поля с токовым слоем возникает в результате действия двух факторов: наличие пучка (внешнего тока) и наличие сильного градиента проводимости плазмы на

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 04-01-00244-А).

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

краю пучка. Градиент проводимости вызывается пучковой турбулентностью, возникающей в пространстве, занятом пучком [13].

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

3. Ранее изучались периодические по одной из координат задачи (токамак) либо процессы, протекающие в неограниченной плазме. В случае ГОЛ-3 плазма ограничена с торцов. Соответственно возникает необходимость изучения влияния торцов на течение. Заметим, что особенности задачи, отмеченные выше в пп. 1 и 2, могут быть изучены в двумерной постановке. Однако ответы на вопрос о том, какая мода будет доминирующей, и вопрос о взаимодействии мод могут быть получены только при решении трехмерной задачи.

4. Магнитное поле на установке ГОЛ-3 является гофрированным.

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

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

Так как в открытых ловушках продольное магнитное поле велико, будем описывать течение плазмы уравнениями редуцированной магнитной гидродинамики [4], которые в цилиндрической геометрии в безразмерных переменных имеют вид

дФ д Ф

— + (УУ)Ф = — + п(АФ + л,);

ди + (УУ)и = ¿МН*) + Ц +

V = (К, Уф) = (г-1дФ/др, -дФ/дг); Н = (И, Нф) = (г-1дФ/др, -дФ/дг); ДФ = -и;

л = -ДФ;

. 1 д ( д\ 1 д2

Д = -тг- г— +

1.1)

1.2)

1.3)

1.4)

1.5)

г дг \ дг) г2 др2

Здесь V — скорость плазмы, которая в данной модели является несжимаемой и имеет только полоидальные (г и ф) компоненты. Соответственно, вводятся функция тока Ф и завихренность ш; H — полоидальное магнитное поле. В этой модели тороидальная компонента (z-компонента) магнитного поля постоянна. В использованной в настоящей статье нормировке Hz = 1; Ф — вектор-потенциал магнитного поля; jz — плотность тороидального тока, которая включает в себя как ток плазмы, так и ток пучка; j — плотность тока пучка; п и v — безразмерные сопротивление и вязкость плазмы.

Для нормировки использовались параметры установки ГОЛ-3. Как уже упоминалось, магнитное поле нормировано на величину его тороидальной компоненты, равной 50 кГс, скорость — на скорость Альвена Va = Hz/(4пр)1/2 = 3.5 • 108 см/с, где р = 1015 см-3 — плотность плазмы, которая в данной модели постоянна. В качестве единиц длины использовался радиус плазмы a = 5 см, времени — величина a/Va = 1.5 • 10-8 с и т. д.

Решение рассматривается в области 0 < r < 1, 0 < z < z0. В размерных величинах z0 = 12 м. На границах z = 0 и z = zo скорость плазмы полагалась равной нулю (Ф = 0). В используемой модели при z = 0 и z = z0 для вектор-потенциала Ф граничных условий не требуется. Граница r = 1 полагалась идеально проводящей стенкой (Ф = const) с условиями прилипания для скоростей Vr = V^ = 0. Последнее дает дФ/dr = дФ/дф = 0 или Ф = 0 и граничные условия для ш (условия Тома, см. далее).

В начальный момент времени плазма покоится (V = 0), а распределение магнитного поля задается с помощью вектор-потенциала вида

Ф = Ф0(г) + Ф1 rm(1 - r2) cos(m^ - 2nnz/z0). (1.6)

Здесь Ф0(г) — распределение вектор-потенциала, соответствующее равномерному распределению по площади r < 4 см фонового тока jz = —5 кА, а Ф1 — амплитуда возмущения. Само возмущение, как легко видеть, соответствует винтовой моде с полоидальным номером m и тороидальным номером n. Зависимость возмущения от радиуса выбиралась из следующих соображений. На границе r = 1 величина Ф должна быть постоянной. При r ^ 0 возмущение должно быть регулярным, т. е. быть представлено в виде степенного ряда по декартовым координатам. Для этого необходимо, чтобы оно содержало степень r не ниже m. Заметим, что при малых амплитудах зависимость возмущения от радиуса принципиального значения не имеет. Кроме того, использовалось также возмущение, в котором Ф отличалось от Ф0(г) только в одном узле расчетной сетки, в котором Ф = Ф0(г) + Ф1.

Плотность тока в пучке jb постоянна при r < Гь = 2.5 см и соответствует полному току пучка 30 кА. При r > Гь пучок отсутствует. Такое распределение jb имеет место при t < tb (tb = 600 или 9 • 10-6 с в размерной форме). При t > tb пучок отсутствует.

Сопротивление п при t < tb имеет следующее распределение:

П = п0 при r > rb + ô,

П = Пь П0 при r < гь, П = Пь + (П0 — Пь)(г — гь)/0 при гь < г < гь + ô.

Ширина переходной области ô связана с размытием пучка и в представленных ниже расчетах равнялась 0.25 см. В безразмерном виде п0 = 4.5 • 10-6, что соответствует спитце-ровской проводимости при температуре 100 эВ. Величина пь в 10... 1000 раз превышала П0. Увеличение сопротивления в области пучка соответствует турбулентному значению проводимости, точное значение которого неизвестно.

При t > tb сопротивление п = По во всей области, так как после выключения пучка турбулентность исчезает. Следует отметить, что температура плазмы в области пространства, занятой пучком, за время существования пучка может существенно повышаться, превышая 1 кэВ. После выключения пучка температура в области, где был пучок, уменьшается до 100 эВ за счет теплопроводности на торцы камеры за время примерно 50 мкс [13]. Это время намного больше характерного времени (< 10 мкс) рассматриваемых в настоящей работе процессов. Высокая температура в отсутствие турбулентности означает высокую проводимость плазмы. Однако развитие тиринг-неустойчивости определяется проводимостью плазмы в окрестности резонансной поверхности, которая находится вне пучка. Здесь температура не сильно отличается от 100 эВ, а сопротивление — от по. Поэтому в расчетах для простоты положено n(t > tb) = по.

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

Вязкость плазмы в рассматриваемом случае мала [14]. В расчетах полагалось v = const = 0. 1 п0. Конкретное ее значение при v ^ п не влияет на результат.

2. Разложение по винтовым модам

Одна из основных проблем, рассматриваемых в настоящей работе, — проблема состава мод возникающего течения. В случае винтовой симметрии, когда d/dz = —R-1d/dp (R = 2nn/(mz0)), МГД-уравнения практически совпадают со случаем d/dz = 0, но в роли Яф выступает величина Hs = H^ — (r/R)Hz [6]. Соответственно, если конфигурация магнитного поля такова, что Hs имеет разные знаки по разную сторону поверхности rmn, называемой резонансной, то развивается процесс пересоединения (аннигиляция магнитных потоков разного знака, тиринг-неустойчивость). Условие существования поверхности Hs = 0 можно переписать в виде q = Hz/(НфД) = 1. При q > 1 конфигурация неустойчива по отношению к тиринг-модам. Расчеты показывают, что при параметрах плазмы, соответствующих установке ГОЛ-3, неустойчивыми могут быть возмущения с различными m и n. Например, m = n =1; m =1, n = 0; m =1, n = 2 и т. д. В случае периодических по z граничных условий (токамак) возможно развитие (в том числе нелинейное) одной моды. В случае открытой ловушки различные моды связаны друг с другом в силу граничных условий. Поэтому встает вопрос о составе мод возникающего течения.

Для того чтобы характеризовать количественно значимость различных мод, заметим, что в случае одной моды какая-либо величина f (r, р, z) является постоянной вдоль линий p(z) = 2nnz/(mz0)+ const и периодичной по р с периодом 2n/m. Поэтому можно говорить о том, что функция f содержит в себе моду m, n вида

m— i zo

fmn (r,p) = (mzo)-1V f I r, p +---1--,z)dz. (2.1)

t^i \ m mzo J

t=0 о

В качестве амплитуды моды можно взять величину Amn = max | fmn |. Ниже в качестве

f используется w, поскольку эта функция отлична от нуля только при наличии неустойчивости.

В примененном для решения поставленной выше задачи численном методе различные функции f задавались в узлах равномерной сетки pj, . Вычислить fmn в точках pj

можно следующим образом. Находим /р(^к) — дискретный фурье-образ по ], где р — номер

т- 1 ко

фурье-гармоники. Затем вычисляем величину / = (тк0)-1 ^ /(гк)е-гра1к. Здесь аш =

1=0 к=1

2п//т + /(т^0), к0 — число точек сетки по г. Искомая функция /тп получается

обратным фурье-преобразованием по р функции /*.

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

Рассмотрим конечно-разностный алгоритм, разработанный для решения поставленной задачи.

Расчетная сетка была равномерной по < и г с шагами к^ и кг и неравномерной по г. Поскольку наибольшие градиенты возникают при г > гь, сетка по радиусу сгущалась именно в этой области. Функции Ф и ш вычислялись в целых узлах гг, <, гк, а Ф и — в полуцелых узлах Гг+1/2, <£¿+1/2, ^+1/2. При этом го = 0, гко = го, Гг=о = 0, (гго+1 + Гго)/2 = 1, гг+1/2 = (гг+1 + гг)/2, г—1/2 = —г1/2. Здесь г0, и к0 соответствуют числу узлов в направлениях г, < и г.

Пространственные производные, входящие в уравнения (1.1)—(1.5), аппроксимировались со вторым порядком точности естественным образом. Например, член д^/дг в уравнении (1.2) вычислялся в целых точках. Это делалось следующим образом:

_ < ^'.г >г,^,к+1/2 — < .¿г >г,^,к—1/2

где

< ¿г >г,^,к+1/2 =

_ (^гг+1/2,^+1/2,к+1/2 + ^'гг—1/2,^+1 /2,к+1/2 + ^гг+1 /2,^ — 1 /2,к+1 /2 + ,Ы-1/2,^ — 1/2,к+1/2) = 4 '

Исключение составляли члены (VV)Ф и Они аппроксимировались тремя различ-

ными способами.

1. Член (VV)Ф вычислялся с первым порядком точности по ближайшим узлам с учетом знака скорости. Член (VV)ш записывался в виде ^у ^ш) и аппроксимировался так:

г'3'к гг (гг+1/2 — гг-1/2) 2гг /2) ' '

В центре координат при г = 0, где функции не зависят от индекса ], имеем

2 ¿о

а1у ^ш)г=о,к =-- У^ ?н=1/2л. (3.2)

г1/2^0

Здесь 5гг+1/2,^',к = Кг+1/2,.?,кгде ш* = шг,^,к при Кг,+1/2,.?,к > 0 и ш* = шг+1,^,к при Кг+1/2,.?,к < 0. Аналогично аппроксимировалась величина ^ = У^ш.

Подробно об особенностях аппроксимации различных дифференциальных операторов в цилиндрической системе координат см. в работе [11].

2. При аппроксимации второго порядка с помощью центральных разностей для завихренности использовались формулы (3.1) и (3.2), но в качестве (и, аналогично, использовалось выражение 5гг+1/2,^,к = УТг+1/2,^,к(шг+1,^-,к + шг,^-,к)/2.

3. При аппроксимации третьего порядка с учетом знака скорости для завихренности член, связанный с переносом, записывался как (VV)ш. Производная дш/д< вычислялась в зависимости от знака скорости следующим образом:

2шг,^'+1,к + 3шг,^',к — 6шг,^-1,к + шг,^-2,к п -—- при ^¿к > 0;

—шг,?+2,к + 6шг,?+1,к — 3шг,?,к — 2шг,?-1,к ^ п -^^-----—- при У™ к < 0'

Аналогично вычислялись остальные производные, входящие в интересующие нас выражения. Эта аппроксимация обеспечивает третий порядок точности, причем остаточные члены, например для У^д/д<, имеют вид — 2|У^|д4/д<4. Остаточные члены такого вида слабо влияют на решение в области малых градиентов и эффективно сглаживают мелкомасштабные возмущения.

Заметим, что при У < 0 и вычислении Ф в точке г0 — 1/2 и ш в точке г0 необходимо знать Ф в точке г0 + 3/2 и ш в точке г0 + 2, где они не заданы. В этом случае использовалась аппроксимация первого порядка. Поскольку скорость около границы г = 1 стремится к нулю, это мало влияет на точность решения.

Заметим также, что при вычислении производных в центре координат полагалось, что

/г=-а,^',к /г=а,^'+^о/2,к' (3'3)

Кроме того, в центре координат г = 0 выражение (VV)ш аппроксимировалось как

2 А дш,

2 \—^ дш

(VVш)г=0,k = ~ Угг=0,7,к"д" |г=0,7,к7

/ , ' /г=0,7,к г\

^о дг

где

гг=о,^',к

Фг=1,^-+^о/4,к — Фг=1,^+3^о/4,к

2гг=1

Здесь дш/дг|г=0,^-,к записывалась так же, как и в узлах г = 0, с учетом (3.3).

Остановимся на граничных условиях для ш на стенке. Использовалось условие, аналогичное условию Тома, которое можно получить, разлагая Ф в точке г0 в ряд Тэйлора и учитывая граничные условия на скорость и связь (1.3) и (1.4):

, 4Фго ¿,к . , . _ „

шго+1,^',к = —шго,^-,к + 7.-^' фго+1,^,к + фго,^,к = 0'

(1 — гго)

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

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

Рассмотрим алгоритм перехода с временного слоя ¿га на слой ¿""+1 = ¿га + т. Он состоит из следующих шагов и использует итерации (в — номер итерации, в + а — этапы внутри одной итерации, причем /= /5).

1. При ^ = 0 полагаем ^ = Ф5 = Фп.

2. Находим Ф5 по ^ с помощью конечно-разностного аналога уравнения (1.4). При этом используем разложение в дискретный ряд Фурье по ] (координата и прогонкой по г (координата г). Напомним, что в данной модели оператор Лапласа двумерный.

3. Вычисляем

/ дФ5 \

ф.+° = фп + г / _(узу)ф. + __ + п(Дф5 + ^-пИ .

4. Вычисляем ф5+1/3 посредством уравнения

(т^ _ /)(Ф5+1/3 _ Ф5) = _(Ф5+° _ ф5). Здесь I — единичная матрица, а — конечно-разностный аналог оператора п*Д, где

П*(г, г, ¿) = ш1п(п(г) + т(УФ)2).

Член с т(УФ)2 добавлен для улучшения устойчивости и ускорения сходимости итераций по отношению к альвеновским возмущениям.

Данный этап алгоритма осуществлялся с помощью разложения в дискретный ряд Фурье по индексу ] (координате и с последующей прогонкой по индексу г (координате г). Индекс к (координата г) здесь выступает в роли параметра.

5. Вычисляем ф5+2/3 посредством уравнения

(тЬ2 _ I)(Ф5+2/3 _ Ф5) = _(ф5+1/3 _ ф5),

¿2 = _К£ + (п + т(УФ)2 _ П*) 1 д (г л •

дг г дг \ дг/

При реализации данного этапа использовалась пятиточечная (в случае аппроксимации членов уравнений, связанных с переносом с третьим порядком точности) прогонка по индексу г.

6. Вычисляем Ф5+1 посредством уравнения

(тЬ3 _ I)(Ф5+1 _ Ф5) = _(Ф5+2/3 _ Ф5),

1 д 1 д2 ¿3 = — ^ — + (п + т(УФ)2 _ п*)

г ф д^ г2 д^2

При реализации данного этапа использовалась пятиточечная циклическая прогонка по индексу .

7. Вычисляем посредством уравнения

= + т ^_(уУ V + а1у(н8+1^+1) + + Vд^^ .

8. Вычисляем ф5+1/3 посредством уравнения

(тЬ4 _ I)(^+1/3 _ = _(Ф5+° _ Ф5),

¿4 = V Д.

На этапах 7 и 8 в формуле Тома в граничных условиях для ш использовалась Ф*.

9. Находим Ф*+1 по ш*+1/3 с помощью конечно-разностного аналога уравнения (1.4).

10. Вычисляем ш*+2/3 посредством уравнения

(тЬ5 — I )(ш*+2/3 — ш*) = — (ш*+1/3 — ш"),

д

¿2 = — К —' дг

11. Вычисляем ш*+1 посредством уравнения

(тЬ6 — I )(ш'+1 — ш*) = —(ш*+2/3 — ш*),

¿2 = — Ч «

г д^

12. Вычисляем ш в граничных точках с помощью формулы Тома по Ф*+1.

13. Если максимальное по пространственным узлам значение функций |ш*+1 — ш*| и |Ф*+1 — ф*| меньше некоторого заданного числа, то полагаем шга+1 = ш*+1, Фга+1 = Ф*+1. В противном случае возвращаемся на этап 2.

Каждый из приведенных этапов полностью аппроксимирует исходные уравнения. В качестве мга+1 можно выбрать любое и", в > 1. Уже после первой итерации получается схема, обладающая достаточным запасом устойчивости. В расчетах число итераций редко превышало 3. При этом условие устойчивости можно приблизительно записать как т < СЛ,Г/Улр. Здесь С — число порядка единицы; Л,г — минимальный шаг сетки по радиусу; Уар — характерная скорость Альвена, вычисленная по полоидальному магнитному полю. Это ограничение на т носит физический характер и не является обременительным.

Заметим, что выделение этапов 4, 8 и 9 связано с историей создания кода. По-видимому, этап 9 можно опустить, а члены, связанные с вязкостью и сопротивлением плазмы, расщепить только по пространству, добавив их в этапы 10, 11 и 5, 6.

Описанный алгоритм может быть достаточно просто "распараллелен".

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

Ниже для иллюстрации результатов расчетов будет использован вариант с возмущением (1.6) при параметрах т =1, п = 0, Ф1 = 10-4. Здесь сопротивление плазмы вне пучка П0 = 4'5• 10-6, а в области, занятой пучком, Пь = 1000^о. Расчеты показали, что происходит следующее.

Благодаря наличию пучка и градиентов сопротивления на временах £ > 100 в отсутствие возмущения (Ф1 = 0), нарушающего аксиальную симметрию, возникает следующая конфигурация (рис. 1). Плотность тока ^ в области пучка близка к ^ь. В переходной области на краю пучка находится токовый слой — узкая область пространства с большой плотностью тока. Причем ток в слое направлен в противоположном направлении по отно-

1 2п

шению к току пучка (^ < 0). На периферии ^ мал, а изменения полного тока JJ ^г^г^

00

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

Рис. 1. Типичное распределение плотности тока к моменту выключения пучка (Ь = 600) в случае отсутствия возмущения (одномерная задача).

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

Благодаря наличию возмущения (1.6) возникает неустойчивость, приводящая к появлению сильной асимметрии в распределении (рис. 2).

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

Пусть дано уравнение А1 = пАхх. Тогда для величины ] = —Ахх имеем уравнение + с]х = П3хх + 'Цхх3, где с = — 2^х. Таким образом, наличие градиентов сопротивления П эквивалентно сносу тока ] со скоростью с. В случае нашей задачи эта величина очень большая. Она влияет так, что распределение прижимается к точке г = Гь + 5. При г > Гь величина п постоянна и описываемый эффект отсутствует. Наличие членов, связанных с Уп, не только является причиной формирования конфигурации с резонансной поверхностью, но и приводит к сжатию поверхности Гь + 5 возмущения тока , возникающего благодаря неустойчивости. В результате плотность тока достигает больших величин, а толщина по радиусу возмущения становится чрезвычайно малой. На рис. 3 показана зависимость величины ]т = тах | от времени при различных параметрах конечно-разностной сетки.

Всплеск ]т в начальный момент времени, а затем уменьшение этой величины связаны с формированием аксиально симметричной конфигурации с обратным током на границе пучка и постепенным его затуханием благодаря сопротивлению плазмы. Реально ]т достигается в области токового слоя (обратного тока). При £ > 100 величина ]т опять растет, что связано уже с развитием неустойчивости и отклонением от аксиальной симметрии. При £ > 200 токовый слой становится настолько узким, что дальнейший расчет не представляется возможным. Поэтому мы поступаем следующим образом.

Заметим, что возмущение концентрируется в окрестности нейтральной поверхности, вблизи границы пучка. Другие области пространства возмущены слабо. Можно ожидать, что дальнейшее развитие неустойчивости при наличии Уп будет приводить только к увеличению плотности тока в этой окрестности. Магнитные потоки в результате пересоединения уменьшатся не намного. Поэтому в момент времени £ = 200 выделяем составляющие

Рис. 2. Распределение плотности тока в сечениях г = го/4 (вверху), г = го/2 (в середине), г = 3го/4 (внизу) в моменты времени Ь = 200 (слева) и Ь = 700 (справа).

функций Ф и ш, зависящие от

2п 2п

Ф = Ф — (2п)-^ У Ф^, ш = ш — (2п)-^ У ш^

о о

Рис. 3. Зависимость ]т от времени: в случае отсутствия возмущения (одномерная задача) — кривая 0; в случае амплитуды возмущения Ф1 = 10-4 для различных сеток (кривая 1 — число узлов равномерной сетки по г, р, г 200 х 64 х 32, кривая 2 — сетка в 2 раза реже, кривая 3 — сетка в 2 раза чаще).

Рис. 4. Силовые линии поля Нг, при т — п — 1 (изолинии функции Ф + (2п/г°)г2/2) в момент времени Ь — 760 в сечении г — го/2.

Рис. 5. Линии тока (изолинии Ф) в момент времени t = 700 в сечении z = Zo/2: знаками "+" и " —" обозначены области, соответствующие Ф > 0 и Ф < 0.

Находим решение одномерной задачи (возмущение Ф1 = 0) в момент времени £ = 600, которое обозначим как Ф0 и и°. Функции Ф = Ф0 + ф и и = и° + и используем в качестве начальных данных для расчета течения после выключения пучка (£ > ¿ь = 600).

При таком подходе к задаче в момент времени £ > 600 наблюдается достаточно резкое падение максимума тока (рис. 3), что связано с исчезновением Уп и, соответственно,

Рис. 6. Распределение плотности тока в сечении z = Zo/2 в момент времени t = 860.

уширением токового слоя.

Дальнейшее развитие тиринг-неустойчивости близко к классическому: магнитные потоки разных знаков взаимно уничтожаются, на картине силовых линий Hr, Hs четко выделяется характерный остров (рис. 4), плазма втекает в токовый слой, где происходит пересоединение, и вытекает вдоль него (рис. 5). В образующейся конфигурации внешний поток меньше внутреннего. Поэтому в определенный момент времени возмущение выходит на границу r =1 (рис. 6), где затухает после полного уничтожения внешнего потока.

Остановимся на составе мод. В случае начального возмущения (1.6) с m = n = 1 эта мода остается доминирующей и в последующие моменты времени, хотя появляются и другие гармоники. Только в момент времени t > 900, когда выходит на стенку камеры, начинает доминировать мода m =1, n = 0, но в этом случае влияние стенки велико и обычные представления о процессе пересоединения не применимы. Заметим, что вблизи стенки температура плазмы падает, а ее сопротивление растет. В настоящей работе увеличение п при r ^ 1 не учитывается. Изучение пристеночных процессов выходит за рамки настоящей работы.

В случае возмущения (1.6) с m = 1, n = 0 в начальный момент времени отлична от нуля только эта гармоника. Однако к моменту t = 200 амплитуды гармоник m =1, n = 0 и m = n =1 близки друг другу (рис. 7). При t > 620 гармоника m = n =1 становится доминирующей и остается таковой до момента выхода возмущения на стенку (t ~ 900). После этого начинает преобладать гармоника m = 1, n = 0.

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

Доминирующее положение моды m = n =1 объясняется следующим. Анализ одномерной конфигурации показывает, что мода m = n = 1 в сравнении с другими неустойчивыми

rmn 0

модами имеет наибольший как внутренний J Hsdr, так и внешний J Hsdr магнитные по-

0 rmn

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

Рис. 7. Распределение итп (состав винтовых мод для завихренности) в момент времени Ь = 200 и Ь = 670 в случае начального распределения с т = 1, п = 0.

мода т = п = 1 доминирует.

В заключение остановимся на точности получаемых расчетов. Представление о ней можно составить по рис. 3, где приведена зависимость ]т от времени в случае различных пространственных сеток. Из этого рисунка можно сделать вывод о сходимости и о неплохой точности расчетов. Заметим, что ]т является локальной величиной и содержит вторые производные по пространству. Поэтому зависимость от сетки именно этой величины достаточно показательна. При проведении расчетов точность контролировалась также по другим величинам.

Заключение

Расчеты позволяют сделать следующее выводы.

1. Рразрушение токовой структуры, наблюдаемое в экспериментах, может быть объяснено развитием тиринг-неустойчивости. Время развития неустойчивости, полученное в расчетах (порядка нескольких микросекунд), совпадает с наблюдаемым в экспериментах [12] временем разрушения тока.

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

2. Доминирует мода m = n = 1.

3. Поскольку внутренний магнитный поток меньше внешнего, в результате развития тиринг-неустойчивости плазма может выбрасываться на стенку камеры. Возможно, этим объясняются эффекты повреждения стенок, иногда наблюдаемые в экспериментах [12]. Дисбаланс магнитных потоков невелик. Поэтому тепловая энергия теряется далеко не полностью. Тем не менее возникает задача оптимизации параметров установки для избежания этих выбросов.

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

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

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

[2] Березин Ю.А., Федорук М.П. Моделирование нестационарных плазменных процессов. Новосибирск: Наука, 1993.

[3] Porcelli F., Borgogno D., Califani F., Grasso D., Ottayiani M., Pegoraro F. Recent advances in collisionless magnetic reconnection // Plasma Phys. Control. Fusion. 2002. Vol. 44. P. B389-B405.

[4] Днестровский Ю.Н., Костомаров Д.П. Математическое моделирование плазмы. М.: Наука, 1982.

[5] Жуков В.П. Слияние магнитных ячеек в модели электронной магнитной гидродинамики // Физика плазмы. 2002. Т. 28, № 5. С. 451-461.

[6] Жуков В.П. Моделирование тиринг-неустойчивости в нередуцированной двухжидкостной магнитной гидродинамике // Физика плазмы. 2005. Т. 31, № 8. С. 721-732.

[7] Lazzaro E., Ferrero M., Gianoli L., Valdettaro L. Four-field description of fast reconnection of m = 1, n = 1 modes // Physica Scripta. 2000. Vol. 61. P. 624-627.

[8] Магнитное пересоединение в двумерных и трехмерных конфигурациях // Тр. ИОФАН. Т. 51. М.: Наука, 1996.

[9] Lerbinger R., Lucian J.F. A New Semi-implicit method for MHD Computations // J. Comput. Phys. 1991. Vol. 97. P. 444-459.

[10] Жуков В.П., ДудниковА Г.И., Франк А.Г. Численное моделирование формирования и эволюции токового слоя в ограниченной плазме // Вычисл. технологии. 2004. Т. 9, № 3. С. 39-49.

[11] Жуков В.П. Конечно-разностная схема для решения двухжидкостных МГД-уравнений в цилиндрической системе координат // Журн. вычисл. математики и мат физики. 2005. Т. 45, № 1. С. 156-169.

[12] PoSTUPAEY V.V., ÄRZHANNIKOY A.V., ÄSTRELIN V.T. et AL. Rol of q profile for plasma confinement in the multimirror trap GOL3 // Transaction of Fusion Science and Technology. 2005. Vol. 47, N 1T. P. 84-96.

[13] Астрелин В.Т., БурдАков А.В., Поступаев В.В. Подавление теплопроводности и генерация ионно-звуковых волн при нагреве плазмы электронным пучком // Физика плазмы. 1998. Т. 24, № 5. С. 450-462.

[14] Брагинский С.И. Явления переноса в плазме // Вопр. теории плазмы. М.: Атомиздат, 1963. Вып. 1. С. 183-272.

Поступила в редакцию 1 декабря 2006 г., в переработанном виде — 21 декабря 2006 г.