CC BY
34
Вычислительные технологии
Том 8, № 5, 2003
ПРОНИКНОВЕНИЕ БЕГУЩЕЙ ВОЛНЫ В ПЛАЗМУ ПРИ НАЛИЧИИ РЕЗОНАНСНОЙ
ПОВЕРХНОСТИ*
В. П. Жуков, H. J. de Blank, И. В. Шваб Институт вычислительных технологий СО РАН, Новосибирск, Россия FOM-Instituut voor Plasmafysica, The Netherlands e-mail: zukov@ict.nsc.ru, shva@net.ict.nsc.ru
A penetration of a running magnetic field perturbation into plasma is investigated. It is significant that the perturbation is running and not simply oscillating. It is shown, that the penetration process has different forms, depending on the proportion of phase velocity and Alven velocity, and is accompanied by the plasma acceleration. The latter plays a crucial role in the process.
Введение
Динамический эргодический дивертор (DED) для токамака TEXTOR-94 представляет собой систему катушек переменного тока для создания возмущений магнитного поля, разрушающих магнитные поверхности вблизи резонансной поверхности q = 3, что приводит к стохастизации магнитного поля [1, 2]. В отсутствие токов плазмы, наведенных полем DED, стохастическое магнитное поле представляет собой суперпозицию равновесного поля и поля, создаваемого катушками эргодического дивертора в вакууме. Свойства этого стохастического поля описаны в работах [2-9]. Изучение влияния плазменных токов, индуцированных возмущением магнитного поля, и вычисление момента сил, действующего на плазму со стороны возмущения магнитного поля, производились в линейном приближении. Причем при описании распределения тока в окрестности резонансной поверхности использовались дополнительные предположения [10-13].
В настоящей работе описаны исследования влияния упомянутых выше токов на процессы проникновения поля DED в плазму и ускорения плазмы. Показано, что проникновение возмущения магнитного поля сопровождается ускорением плазмы, которое в свою очередь решающим образом влияет на проникновение магнитного поля. Этот процесс носит сугубо нелинейный характер. Поэтому изучения его в линейном приближении явно недостаточно.
Заметим, что стохастизация силовых линий магнитного поля имеет место только в трехмерном случае. Однако процесс проникновения магнитного поля и ускорения плазмы может быть рассмотрен в двумерном случае. В настоящей работе процесс проникновения изучается в двумерной плоской геометрии в модели как несжимаемой, так и сжимаемой магнитной гидродинамики.
* Работа выполнена при финансовой поддержке Euratom-FOM Association, NWO и Euratom.
( Институт вычислительных технологий Сибирского отделения Российской академии наук, 2003.
Необходимо также отметить, что используемые в настоящей работе модели имеют высокий уровень упрощений. Тем не менее в них учтены два главных для процесса проникновения момента. Во-первых, наличие резонансной поверхности в конфигурации магнитного поля. Во-вторых, то, что возмущение магнитного поля представляет собой бегущую, а не стоячую волну, что является причиной ускорения плазмы в целом.
Ниже дается описание основных предположений, которые позволяют получить двумерную модель из более общих трехмерных уравнений одножидкостной магнитной гидродинамики (МГД). Затем приводятся математическая постановка задачи и результаты расчетов. Далее обсуждается вопрос соответствия численных результатов и реальных параметров DED и сделано заключение.
1. Основные уравнения
Для полного описания DED в приближении МГД необходимо рассмотреть трехмерную задачу в тороидальной геометрии с учетом вакуумной области, где расположены катушки. Распределение тока jc(r, <,Z,t) в катушках можно приближенно описать как
jc ~ f (<) cos(m< — n(/Ro + ut) 8(r — rc),
где < — полоидальный угол; Z — тороидальная координата; R0 — большой радиус токамака; rc — радиус расположения катушек; m = 12 и n = 4 — номера основных полоидальных и тороидальных гармоник; u — частота тока катушках. Функция f определяется конструкцией DED. Она примерно равна единице на внутренней поверхности тора и близка к нулю на внешней [2, 13].
Трехмерная задача очень сложна. Ниже используется более простая двумерная постановка, которая хорошо описывает основной механизм ускорения плазмы и проникновения в нее магнитного поля. Основные этапы перехода от трехмерной к двумерной модели заключаются в следующем.
1. Прежде всего перейдем от тороидальной геометрии к цилиндрической, оставаясь в рамках трехмерного случая. Покажем, в каком случае и как это можно сделать. Для этого введем малый радиус токамака a, тороидальную компоненту магнитного поля Hz, скорость быстрого магнитного звука VZ = Hz/(4пр)1/2, плотность плазмы р, время tz = а/Vz, отклонение величины тороидального магнитного поля от его значения в вакууме hz = GHz — (Hz)/(G-1), где G = 1 + (r/R0)cos < — тороидальный метрический тензор,
a 2п
(• • •) = (па2)-1 / f • • • r dr d<.
0 0
Предположим, что аспектное отношение a/R0 ^ 1, и рассмотрим движение, которое определяется в основном полоидальным магнитным полем H ~ (a/R0)Hz и имеет следующие свойства:
а) характерный масштаб течения в полоидальном направлении < а, а в тороидальном направлении > R0;
б) характерная скорость плазмы ~ (a/R0)Vz, отклонение плотности плазмы от начального значения ~ (a/R0)p, hz = (a/R0)2Hz и т.д.
Предположим, что все функции состоят из двух частей. Одна — медленная с временным масштабом ~ (R0/a)tz, другая — быстрая с временным масштабом ~ tz. Предположим также, что быстрая часть всех функций, за исключением hz, в a/R0 раз меньше медленной. Обе части hz будут одного порядка.
Используя сделанные выше предположения, разложим исходные МГД-уравнения, записанные в тороидальных координатах, в ряд по a/R0, пренебрегая геометрическими факторами. В первом порядке разложения получим уравнения, которые совпадают с МГД-уравнениями в геометрии прямого цилиндра. Роль цилиндрических компонент r, p и z векторов играют компоненты r, p и ( векторов в тороидальной геометрии. Исключение составляет величина Hz, под которой подразумевается GHz.
Сделанные выше предположения действительно имеют место в случае DED, поскольку
i—1
частота возмущения и << t^ , магнитное поле катушек намного меньше тороидального магнитного поля, характерный пространственный масштаб возмущения магнитного поля катушек соответствует указанным выше требованиям и т.п.
2. Цилиндрическая геометрия, в отличие от тороидальной, позволяет перейти от трехмерной задачи к двумерной в случае винтовой симметрии, при которой все величины зависят от p — nz/(mR0), а не отдельно от p и z. Следовательно, можно записать
d/dz = —n/(mR0)d/dp.
На втором этапе упрощений положим f (p) = 1, что соответствует катушкам DED, расположенным вокруг всего токамака, а не только на его внутренней поверхности. Это позволяет совершить упомянутый переход к двумерной задаче. Уравнения МГД при наличии винтовой симметрии близки к обычным двумерным МГД-уравнениям в случае d/dz = 0. Приведем уравнения для магнитного поля, которые имеют наибольшее отличие от случая d/dz = 0:
dé тdé т 1 dé
— + Vr— + Vs--- =
dt dr r dp
= v (g ( 1 - ( r~—) + -d2é ^ 2H
r dr \g dr ) r2 dp2 J Rg
dHg = 1 d(g-1rSr ) 1 dS± dt g r dr r dp Rg dt '
Sr = — Vr Hg + Vg Hr + 'qdHg/dr,
V g Hs + vgr 1'
22
Ss = —VsHg + Vg Hs + Vgr-ldHg/dp,
д = 1 + г2/Я2.
Используемые здесь компоненты в и д векторов связаны с обычными цилиндрическими компонентами соотношениями
/а =^— ая!х'1д =+
где Я = тЯ0/(иа). В приведенных выше уравнениях ф и Нд — д-компоненты вектор-потенциала и магнитного поля; п — безразмерное сопротивление плазмы.
3. Упростим эти двумерные уравнения, используя разложения по параметру Я ^ 1 и тот же скейлинг, как и на этапе 1. Дополнительно предположим, что движение плазмы сосредоточено в тонком слое шириной I ^ а около границы г = а. Пренебрегая членами порядка Я-3 и (I/а)Я-2 в уравнениях для магнитного поля и скорости, получим следующие двумерные уравнения в плоской геометрии, записанные в безразмерном виде:
Р1Е = —V (р + — 3V ' v) + jzVé + vAv,
л -Ф
Р-Ж = (ну)я* + ^Аж* = -пь,
н = (Нх, Ну) = (дф/ду, -дф/дх), V = (ьх,Уу), Зх = -Аф + 2
Н я*
V
д2 д2 - д
+ 1Г2, -Л = Б7 + ^),
— —\ А
дХ ду / дх2 ' ду2' -Ь дЬ
(1.1)
дН 2 дф
+ V. №) = ^ №) + V. (п^) + я* %+^ (^=0,
1 + V ■ ы) = -р V ■ (V) + иЯ + п(з1 + ^Н*)2) + V ■ №(р/р)).
7 - дЬ
Здесь V — безразмерный коэффициент вязкости плазмы; Я — нагрев, связанный с работой вязких сил; х — безразмерный коэффициент теплопроводности. Существует следующее соответствие между величинами, входящими в (1.1), и величинами в цилиндрической системе координат. Координата х соответствует координате г, координата у — величине а<р, компоненты х, у, г векторов в (1.1) равны компонентам г, в и д в уравнениях, полученных на этапе 2.
Для уравнений (1.1) выполняется закон сохранения энергии с плотностью
V2 + ж2 ^ф)2 + Н2 р
Р—^^ + -—-£ +
2
2
7 - 1
В (1.1) была использована следующая нормировка. В качестве характерного магнитного поля используется величина Н*, которая равна Н3 на границе плазмы г = а, в качестве характерной длины — расстояние Хо = а - г3 между резонансной поверхностью (поверхность на которой Н3 = 0) и границей плазмы, в качестве скорости — жа = Н*/(4пр)1/2, в качестве времени — х0/жа и т.п. В этих единицах Нх приблизительно равно Я*/2, где Я* = Яа/хо 1.
Обычная процедура разложения в ряд (1.1) по Я* ^ 1 позволяет получить более простую модель несжимаемой плазмы:
дф/дЬ + V ■ Vф = п(Аф - 1),
дv/дt + V ■ Vv = VАv - АфVф - Vp,
(1.2)
V- V = 0.
4. В целях исключения области вакуума из рассмотрения зададим возмущение магнитного потока ф формулой
фр сов(2п(у - жрЬ)/уо). (1.3)
Это соответствует постоянной амплитуде возмущения магнитного потока, равной фр. В действительности подобная зависимость имеет место только для тока в катушках, а не для амплитуды магнитного потока на границе плазмы. Благодаря высокой проводимости
возмущенное магнитное поле не может немедленно проникнуть в плазму (даже на границе) . Поэтому возмущение магнитного поля на границе меняется самосогласованно от нуля до своего значения в вакууме.
5. Для сокращения времени счета будем предполагать, что задача симметрична относительно резонансной поверхности г3. В реальной геометрии токамака радиус г = г3 + х0 = а соответствует границе плазмы, а область г = г3 - х0 — ядру плазмы.
2. Математическая постановка задачи
Рассмотрим более подробно математическую постановку задачи в случае модели несжимаемой плазмы (1.2). Будем использовать координату у* = у - жрЬ вместо у. Подчеркнем, что преобразование Галилея при этом не производится. Скорость жу остается скоростью в лабораторной системе координат. Относительно координаты у* возмущение (1.3) является стационарным. Заметим, что в случае стоячей волны, когда возмущение задается формулой фр сов(2пу/у0) сов(^Ь), подобный переход к стационарному возмущению невозможен. Тогда уравнения (1.2) примут вид (далее по тексту * опускается)
дф/дЬ + V ■ Vф - Жрдф/ду = п(Аф - 1),
дv/дЬ + V ■ Vv - жрдуV = VАv - АфVф - Vp, V- V = 0.
Решение ищется в области 0 < х < 1, 0 < у < у0. Полагаются периодичность по у ^ у + у0 и симметрия относительно х = 0. При границе х =1 задаются условия на скорость жх = 0, джу/дх = 0 и магнитный поток
ф(х = 1) = 1/2 - фр(1 - ) со8(2пу/у0).
Характерное время включения возмущения Ь0 < 10. Расчеты показывают, что решение задачи практически не зависит от конкретного значения Ь0. Начальные условия имеют вид V = 0, ф = х2/2 (нейтральной поверхность х = 0). Коэффициент вязкости V полагается постоянным. Для определения сопротивления плазмы используется формула
п(х) = Пр + Пе ехр(-(х - 1)2/Ь2).
В расчетах полагалось Ь = 0.07 и пе = 0.1 ^ пР. Таким образом, п = пР =сопэ1 почти во всей области и резко возрастает до значения пе на границе. Резкое увеличение сопротивления на границе и постановка граничного условия джу/дх = 0 являются очень грубой попыткой учесть отрыв плазмы от стенки камеры.
Заметим, что вышеизложенная постановка задачи в случае жр = 0 очень близка постановке в работе [14], где исследовалось вынужденное присоединение, вызванное возбуждением альвеновских волн.
3. Конечно-разностная схема
Прежде всего отметим, что в настоящей работе не используются традиционные в динамике несжимаемой жидкости переменные функция потока — вихрь. Это связано с тем, что поток плазмы через сечение х = const меняется со временем и заранее не известен, соответственно, не известны и значения функции потока на границах. Для вычисления этих значений нужно прилагать дополнительные усилия. Использование переменных скорость — давление свободно от этого недостатка.
Поставленная выше задача решается конечно-разностным методом. Для этого введем сетку по x: Xi = ihx (i = 0, ...,io, hx = 1/(io + 1/2)) и y: yj = jhy (j = 1,...,jo, hy = yo/jo). Точки Xi, yj будем называть целыми. В них будем вычислять давление p и div v. Величины v, ф, у, V будем вычислять в полуцелых точках xi+1/2 = (i + 1/2)hx (i = -1,..., i0), yj+1/2 =
(j + l/2)hy.
Будем аппроксимировать Vp и div v с помощью разностных операторов Vh и divh
(VhP)i+1/2,j+1/2 =
pi+1,j+1 + pi+1,j — pi,j+1 — pi,j pi+1,j+1 + pi,j+1 — pi+1,j — pi,j
2hy
(divhv)i,j =
_ Vxi+1/2,j+1/2 + Vxi+1/2,j-1/2 — Vxi-1/2,j+1/2 — Vxi-1/2,j-1/2
2h + + Vyi+1/2,j+1/2 + Vyi-1/2,j+1 /2 — Vyi+1 /2 j-1 /2 — Vyi-1 /2 j-1 /2
2ку
Остальные члены будем аппроксимировать центральными разностями. Например,
(дф/дх)г+1/2,3+1/2 = (фг+3/2,3+1/2 — фг-1/2,3+1/2) / (2ЬХ) , (д2ф/дУ2)г+1,/2,3+1/2 = (фг+1/2,3+3/2 — 2фг+1/2,з+1/2 + ф%+1/2,з-1/2)/Щ и т.п. Производные по времени аппроксимировались как (рп+1 — )/г, где т = Ьп+1 — Ьп — шаг по времени. В дальнейшем для конечно-разностных операторов будем использовать те же обозначения, что и для дифференциальных.
Для повышения устойчивости схемы правую часть всех уравнений лучше брать неявно. Подставляя \п+1 = \п — тУрп+1 — тАфп+1Уф + ... (под "..." подразумеваются несущественные для изложения члены) в уравнения для дф/дЬ и V = 0, мы приходим к
необходимости решать уравнения фп+1 — фп
ф-^ = — ^У^1 + т(Уф)2Афп+1 + уАфп+1 + ...,
т
тУрп+1 = (Иу^п + ...).
Для решения этих уравнений использовался следующий итерационный процесс (ниже в будет означать номер итерации, а и — номер шага по времени):
ф8=0 = фп, vs=0 = Vп, р3=0 = рп, фз+1/2 — фп
т
lhs+1 — фп
ф-^ = хф ф^1/2 + щ фs+1 — т Ws)2^s, (3.1)
Хф г+1/2 + уф Ф — т (Vr)2^s
s 2 s
- = Хф ф ' + уф ф — т \ \ц
т
д д2 хФ = -vxir + (п + т (Ws)2)
дх дх2
д д2
Уф = -(vy _ vp) — + (n + т(Ws)2)
^ "Р> ду ' Vl ' v ^ ; ду2;
vs+1/2 — vn
vs+1 _ vn
Xvvs+1/2 + yvvs _ A^s+1V-0s+1 _ Vps
т
Xvvs+1/2 + yvvs+1 _ A^s+1V^s+1 _ Vps, (3.2)
д д2 , ч д д2
Xv = + v 7T2., yv = _(vy _ vP^ + v
дх дх2 ду ду2
V* = 0. (3.3)
Для ^+1/2 3+1/2 во внутренних точках (г = 0,..., г0 — 1) имеем
V* = V* — тУ^+\ V* = у^1 + тУ^. (3.4)
На границах х = 0 и х =1 для V, V* и V* имеем одинаковые граничные условия:
Ух -1/2,3+1/2 = —Ух +1/2,3+1/2, Уу -1/2,3+1/2 = Уу +1/2,3+1/2,
Ух го+1/2,3+1/2 = 0, Уу ¿0+1/2,3+1/2 1ЬпУу ¿0-1/2,3+1/2. (3.5)
Здесь 1ьп = 0 в случае условий прилипания (уу (х = 1) = 0) и 1ьп = 1 при условии проскальзывания (дУу/дх = 0).
Уравнения (3.3)-(3.5) позволяют найти р5+1. Для этого перейдем к Фурье-разложению по cos(2пmj/jo) (т = 0, ]0/2) и $,т(2пт] /]0) (т = 1,]0/2 — 1). Для Фурье-компонент давления p¿(m) после подстановки (3.4), (3.5) в (3.3) получаем
От^+1(т) — 2p¿(m) + p¿-l(m)) — 8т^+1(т) + 2p¿(m) + p¿-l(m)) =
= Q(i,m), г = 1,...,г0 — 1, (3.6)
2Ст^1(т) — po(m)) — 28т^1(т) + po(m)) = Q(i = 0,т), (3.7)
Ст (p¿о-l(m) — p¿о (т)) — (1ьп + 1)8т^о-1(т) + p¿о (т)) = Q(io,m). (3.8)
Здесь
С = сой2 (пт/]0) 8 = 81п2(пт/]0)
Ст = А2 , 8т = А2 ,
ху
Q¿(m) — Фурье-компоненты функции т-1 (см. уравнения (3.4), (3.5)).
Рассматривая уравнения (3.6)-(3.8), легко видеть, что случаи т = 0 и т = ]0/2 особые. Уравнения (3.6)-(3.8) в этих случаях либо не имеют решения, либо при наличии условий согласования с правой частью имеют бесконечно много решений. Учитывая то, что правая часть пропорциональна от некоторого вектора, можно показать, что условие согласования имеет место. Неоднозначность определения ps+1 при т = 0 связана с тем, что p определено с точностью до постоянной, а при т = ]0/2 — с особенностями схемы. Легко видеть, что прибавление к p функции f, пропорциональной (—1^+3, не влияет на решение, так как Уhf = 0. В программе при т = 0 и т = ]0/2 вместо (3.7) полагалось p¿=0(m) = 0.
Уравнения (3.6) и (3.8), дополненные условием рг=о(ш) = 0, однозначно разрешимы. При этом условие (3.7) будет выполняться автоматически благодаря наличию согласования между правой и левой частями (3.6)-(3.8).
Описанный итерационный процесс хорошо сходится к абсолютно неявной схеме. Можно показать, что в линейном приближении сходимость имеет место при любом т.
Отметим, что если рассматривать р как заданную функцию, то в качестве фп+1, \п+1 можно взять решение уравнений (3.1), (3.2) при любом в = 1, 2,... При этом уже при в = 1 получаемая схема будет обладать хорошим запасом устойчивости (в линейном приближении она будет абсолютно устойчивой).
В расчетах итерации прекращались при выполнении условия
тах < е тах [дьУ/ду1.
г,3 г,3
Расчеты показали, что е = 10-4 являются достаточно малыми. Дальнейшее уменьшение е практически не меняет решение.
Наличие итераций делает целесообразным использование координаты у* = у — урЬ. Как показывают расчеты, в рассматриваемых типах течений существуют длительные квазистационарные (в при использовании координаты у* ) фазы. На этих промежутках времени число итераций будет мало даже при большом т. В случае использования координаты у в соответствующей лабораторной системе координат зависимость от времени и, следовательно, необходимость большого числа итераций будут присутствовать всегда, так как граничное условие для ф будет зависеть от времени.
Количество операций, необходимое для продвижения на единицу времени, определяется величиной т и количеством итераций, необходимых для обеспечения условия несжимаемости. Оптимальный шаг по времени имеет величину порядка кх. Конкретное значение т зависит от стадии решения. В те моменты, когда решение сильно меняется со временем, т выгодно брать малым. На квазистационарных стадиях т можно значительно увеличить. При этом число итераций не превышает 1 на квазистационарных стадиях и 7 в моменты сильного изменения решения со временем.
Описанная схема имеет второй порядок аппроксимации по х и у. Для большинства расчетов было г0/]0 = 64/64. Такое число узлов обеспечивает погрешность получаемого решения заведомо меньше 5%, что для наших целей приемлемо. Точность расчетов контролировалась использованием последовательности узлов сетки. В отдельных случаях число узлов сетки достигало 256/256. Отметим, что пространственная аппроксимация является основным источником погрешности. Погрешность, связанная с т и е, значительно меньше.
Отметим, что уравнения (3.1), (3.2) можно переписать в более компактной форме. Например, вводя операторы X = тХф — 1, У = туф — 1 и функцию Г = —V ■ Уф + Урдгф/ду + п(Аф — 1), (3.1) можно записать как
X (ф3+1/2 — ф3) = ф3 — фп — тГ3,
У(ф3+1 — ф3) = —(ф3+1/2 — ф3).
4. Стационарное состояние
При Ь ^ ж решение задачи, поставленной в предыдущем разделе, стремится к устойчивому стационарному решению, соответствующему движущейся с фазовой скоростью плазме
Рис. 1. Стационарное распределение магнитного поля при фр = 0.05, yo = 10.
и распределению магнитного поля, совпадающему со своим вакуумным распределением. Это решение имеет вид гх = 0, гу = гр, ф = . Для имеем уравнение = 1 с
условиями на границе = 1/2 — фр со8(2пу/у0). Легко видеть, что
1 2 f 2ny\ cosh(2nx/yo)
Фж = _ x — фР cos - -—-;-—.
* 2 V Уo ) cosh(2n/yo)
Типичное распределение показано на рис. 1. Нас интересует процесс установления этого решения, для описания которого будут использованы следующие величины:
= vp 1 max vy — нормированная на фазовую скорость максимальная скорость
Р x,y
плазмы в направлении у;
— (гу) = (гру0)-1 / / гу ^ж^у — нормированная усредненная по всей области скорость
— (Vc) = (vpy0)-1 / vy(x = 0, y) dy — нормированная усредненная по резонансной поверхности x = 0 скорость vy;
— Jm и Jc — максимальные во всей области и на резонансной поверхности x = 0 значения величины | Jz — 1|;
— ь = myax^(x = 0,y)/dy ^ yo cosh(2n/yo) m d0(x = 0,y);
x max = 0,y)/dy 2пфр y dy '
y
— отношение максимального на резонансной поверхности значения компоненты x магнитного поля к той же величине в вакууме. Если магнитное поле внутри плазмы не возмущено, то bx = 0. Если возмущение полностью проникло в плазму (случай стационара), то bx = 1.
Из результатов расчетов следует, что развитие течения происходит по двум различным сценариям в зависимости от того, больше фазовая скорость vp альвеновской скорости на границе (равной 1 в используемых безразмерных переменных) или меньше.
vy;
5. Проникновение поля при ур < 1
В случае гр < 1 в первые моменты времени вблизи границы плазмы возникает скин-слой. Поскольку сопротивление плазмы резко возрастает около границы ж = 1, скин-слой образуется на некотором расстоянии от границы, а не непосредственно на ж = 1.
Благодаря движению плазмы возмущение тока (но не магнитного поля!) начинает проникать вглубь плазмы. В результате на резонансной поверхности х = 0 образуется токовый слой (рис. 2) подобно тому, как это происходит в случае ур = 0 [14], а скин-слой около х = 1 исчезает. Формирование токового слоя требует от нескольких десятков до нескольких сотен альвеновских времен.
В отличие от случая ур = 0, при ур = 0 токовый слой может существовать очень долго. При этом возмущение магнитного поля внутри плазмы, характеризуемое величиной Ьх, мало.
Токовый слой начинает распадаться (.с уменьшаться), а магнитное поле проникать (Ьх увеличиваться) только тогда, когда средняя по у скорость плазмы в токовом слое достигает величины, близкой к ур ((Ус) близко к 1).
Распределение скорости имеет ярко выраженный максимум в токовом слое (рис. 2,3). Это связано с тем, что ускорение плазмы происходит в областях наибольших сил, т.е. в окрестности токового слоя х = 0. В начальные моменты времени плазма ускоряется также в скин-слое.
Вне токового слоя плазма ускоряется в основном благодаря вязким силам и это требует большего времени, чем ускорение плазмы в токовом слое. Распределение магнитного поля, близкое к вакуумному распределению (Ьх близко к 1), достигается тогда, когда средняя по всей области скорость плазмы приближается к фазовой скорости ур ((уу) ^ 1).
Отметим интересную особенность поведения величины Ьх. Эта величина резко растет в
1.2
0.8
0.4
а
V /4 2 .............................3
1
1000 2000 3000
Рис. 2. Картина развития течения при V = 2 ■ 10
-4
= 0.1, фр = 0.025, Пр = 2 ■ 10-4, уо = 10; а
зависимость от времени величин 7с (кривая 1), Ьх (2), (г)у) (3), (гс) (4) и г)т (5); распределение гу при £ = 500 (б) и £ = 1900 (в); г — распределение плотности тока при £ = 1900.
г
р
Рис. 3. Картина развития течения при V = 2 ■ 10 5, гр = 0.1, фр = 0.025, пр = 2 ■ 10 4, уо = 10; а — зависимость от времени величин ,1С (кривая 1), Ьх (2), (г)у) (3), (гС) (4); б — распределение гу при г = 1500.
момент начала распада токового слоя, достигая некоторого максимума, а затем медленно стремится к 1. Этот максимум может существенно превышать 1. Таким образом, имеется эффект усиления возмущения магнитного поля плазмой.
Рассмотрим зависимость времени установления стационарного течения от параметров задачи.
1. Результаты расчетов показывают, что время проникновения магнитного поля тем меньше, чем больше амплитуда возмущения гфр. Причем эта зависимость очень сильная и нелинейная. Этого следует ожидать, поскольку средняя сила, действующая на плазму, (Дфдуф) ~ фр. Помимо амплитуды фр эта сила зависит от разности фаз по у между Дф и дуф. Например, для варианта изображенного на рис. 2, время пересоединения имеет порядок 2200 альвеновских времен. Для амплитуды фр, в два раза большей, это время уменьшится до 500 альвеновских времен.
2. Большее сопротивление пР и меньшая фазовая скорость ур приводят к более быстрому проникновению возмущения магнитного поля.
3. Вязкость V замедляет ускорение плазмы в токовом слое в результате трения с окружающими токовый слой неускоренными слоями плазмы. Поэтому увеличение V приводит к замедлению роста (уС) и увеличению времени начала проникновения магнитного поля
в плазму (начало роста Ьх). С другой стороны, окружающая токовый слой плазма при большем V ускоряется быстрее благодаря более быстрой передаче импульса из токового слоя вязкими силами. Поэтому релаксация распределения магнитного поля к своему вакуумному распределению (уменьшение Ьх от максимального значения до 1) требует меньше времени. В целом время проникновения слабо зависит от V. Картина течения при значениях V, отличающихся на порядок, представлена на рис. 2 и 3.
б. Проникновение поля при ур < 1. Устойчивый случай
В случае сверхальвеновской скорости (ур > 1) процесс проникновения можно разделить на две стадии. На первой стадии формируется скин-слой около границы плазмы х = 1 (рис. 4,
в, г), в котором происходит ускорение плазмы. Более глубокие слои плазмы ускоряются благодаря вязким силам. Токовый слой на х = 0 не образуется.
Дальнейшее развитие течения (стадия II) зависит от профиля уу. Если градиенты уу велики, то возникает неустойчивость, которая будет рассмотрена в следующем разделе. Если градиенты уу малы, как это имеет место в случае большой вязкости, то стадия II наступает после того, как значение ур — уу в расчетной области становится меньше 1 (£ = 1900) (рис. 4). Это связано с тем, что величины ур и уу входят в поставленную задачу только в виде разности ур — уу, поэтому конфигурация, в которой ур — уу < 1, а возмущение магнитного поля сосредоточено на границе, практически эквивалентна начальной (£ = 0) конфигурации в случае ур к 1.
Соответственно, в дальнейшем (£ > 1900) течение развивается, как в случае ур < 1: начинается проникновение возмущения тока вглубь плазмы, и на резонансной поверхности формируется токовый слой, в то время как скин-слой в окрестности х = 1 исчезает. После того как величина уу в токовом слое достигает значения ур (£ = 2500), токовый слой начинает разрушаться, а возмущение магнитного поля проникать вглубь плазмы (Ьх начинает возрастать) и т.д.
Рассмотрим, как зависит время установления конечной конфигурации от параметров задачи.
1. Поскольку ускорение плазмы меньше при меньшей амплитуде возмущения, при меньшем значении фр обе стадии длятся дольше. Например, в случае варианта, изображенного на рис. 4, длительность стадии I, как и стадии II, примерно равна 2000 альвеновским временам. При уменьшении фр в 1.5 раз длительность стадии I увеличивается до 20 000, а стадии II — до 9000 альвеновских времен.
2. Большее ур увеличивает длительность стадии I, но очень слабо влияет на стадию II, так как эта стадия начинается, когда ур — уу к 1 независимо от значения ур.
3. Большая вязкость ускоряет проникновение. Это связано со следующим. В случае большей вязкости увеличивается перенос уу от границы х = 1 вглубь плазмы на стадии I. Соответственно, в момент начала стадии II значение уу(х = 0) в случае большого V ближе к ур, чем в случае малого V. Последнее обстоятельство уменьшает длительность стадии II при большей вязкости. Этот эффект компенсирует дополнительное трение, связанное с вязкостью в скин-слое на стадии I.
4. Сопротивление пр слабо влияет на стадию I, так как ток концентрируется на краю плазмы, где п ^ Пр. Но на стадии II меньшее пр замедляет распад токового слоя х = 0.
Рис. 4. Устойчивый случай: vp = 5, фр = 0.05, = 2 ■ 10 , v = 10 , yo = 10; а — зависимость от времени величин Jm (кривая 1) и Jc (2); б — зависимость от времени величин bx (кривая 1), Jc (2), г)т (3), (г)у) (4), и (г)с) (5); распределение Jz (в) и vy (г) при t = 800 (стадия I); распределение Jz (д) и vy (е) при t = 2050 (начало стадии II).
7. Сильная неустойчивость при vp < 1
Если градиенты уу на стадии I велики, то возникает высокочастотная (с частотой порядка 1) мелкомасштабная неустойчивость. Заметим, что используемая в настоящей работе модель двумерная. Добавление третьей координаты может вызвать и другие неустойчивости.
В зависимости от конкретных параметров задачи неустойчивость может появиться как на стадии I (рис. 5), когда ур — уу > 1, так и в начале стадии II. Неустойчивость возмущает всю внутреннюю область и приводит к более быстрому ускорению плазмы и проникнове-
Рис. 5. Случай сильной неустойчивости гр = 5, фр = 0.05, пр = 2 ■ 10-3, V = 5 ■ 10-3, уо = 10; а — зависимость от времени .]т; б — зависимость от времени гт (кривая 1), (гу) (2), (гс) (3); в — зависимость от времени Ьх; распределение (г), гу (д) и ф (е) в момент времени Ь = 1400.
нию магнитного поля. В результате система достигает конечного стационарного состояния плазмы, минуя стадию токового слоя на х = 0. Неустойчивость исчезает, когда уу достигает значения ур, поскольку исчезает причина, ее порождающая.
Рассмотрим, как зависит от параметров задачи возможность возникновения неустойчивости.
1. Поскольку возникновение неустойчивости связано с градиентом уу, который уменьшается с увеличением V, при увеличении V неустойчивость исчезает.
2. Уменьшение фр уменьшает ускорение плазмы. Соответственно, стадия I продолжается дольше, давая вязкости больше времени для уменьшения градиентов уу. В результате
при достаточно малом фр неустойчивость исчезает.
3. Большее значение ур замедляет ускорение плазмы, и тем самым уменьшается значение . С другой стороны, при большем ур величина уу растет до больших значений, а следовательно, увеличивается и Ууу. Таким образом, большее значение ур задерживает наступление неустойчивости, но не предотвращает ее.
4. Величина пр слабо влияет на возможность развития неустойчивости по причине, указанной в конце предыдущего раздела.
8. Слабая неустойчивость при < 1
В области значений параметров, промежуточной между случаем сильной неустойчивости и случаем устойчивого течения, возникает неустойчивость, которую мы будем называть слабой. Она представляет собой быстрые колебания с периодом порядка одного альвеновского времени, которые имеют значительно меньшую амплитуду по сравнению с амплитудой возмущения фр (рис. 6). Как правило, она возникает в начале стадии II, когда ур — уу < 1. При этом течение развивается точно так же, как и в случае отсутствия неустойчивости, т.е. включает образование токового слоя на поверхности х = 0 и т.п.
1.2
0.8
0.4
■— 0
б
1 / j
/у/
\ 4
///з
2000
4000
6000
Рис. 6. Слабая неустойчивость ур = 5, фр = 0.05, пр = 2 ■ 10-3, V = 7.5 10-3, у0 = 10; а — зависимость от времени .]т (кривая 1), ,1С (2), б — зависимость от времени ут (кривая 1), {уу) (2), (г)С) (3) и Ьх (4); распределение (в) и уу (г) в момент времени Ь = 1709.
9. Влияние сжимаемости
Рассмотрим влияние сжимаемости на процесс проникновения. Для этого используем уравнения (1.1) с начальными условиями
ф = х2/2, V = уг = 0, р = в, Н = Я*/2.
Граничные условия такие же, как и в разд. 3.
Результаты расчетов показывают, что процесс проникновения в случае сжимаемой плазмы не отличается ни качественно, ни количественно от этого процесса в рассмотренном выше случае несжимаемой плазмы даже при Я* ~ 1. Это означает, в частности, что эволюция полоидального магнитного поля и полоидальной скорости плазмы практически не зависит от величины начального давления в и теплопроводности х. Отметим, что слабая зависимость от этих величин типична для процессов пересоединения [15].
Единственный новый эффект, связанный со сжимаемостью плазмы, — это появление ^-компоненты скорости, которая отсутствует в модели несжимаемой плазмы. Но значение уг достаточно мало и уменьшается с возрастанием Я*.
10. Соответствие безразмерных параметров и параметров DED TEXTOR
Выше при описании результатов использовались безразмерные переменные. В данном разделе рассмотрим соответствие между параметрами DED TEXTOR и этими безразмерными параметрами.
Оценим значение s-компоненты магнитного поля на границе плазмы Я,. В случае DED резонансная поверхность = 0 находится около границы плазмы, где плотность тока, текущего в плазме, мала. Поэтому в рассматриваемой области Я<Дг) ~ 1/r. Соответственно, « H^(r = a)a/r — r/(aR) = (r2 — r2)/(arR). Здесь rs = a(RH^(r = a)/#z)1/2 — радиус резонансной поверхности, R = mRo/(an). Принимая во внимание то, что х0 = а — rs ^ а, получим Я, « 2Hzx0/(aR) или
Я, « 2H^(r = а)хо/а. (9.1)
Значение x0 определяется тороидальным магнитным полем и полным током плазмы /р.
Оценим отношение Я, и амплитуды магнитного поля, создаваемого катушками DED. Радиус плазмы равен а = 46 см. Катушки размещены на радиусе rc = 53 см. Расстояние между катушками с противоположным током приблизительно 2nrc/(2m) ~ 12 см. Эта величина сопоставима с расстоянием между плазмой и катушками. Следовательно, в данном случае формула мультипольного разложения для магнитного поля |r — rc|-(m+1) не обладает большой точностью. Результат будет промежуточным между затуханием этого типа и затуханием в случае одной катушки вида 1 /1 r — rc |. В случае затухания типа 1/1 r — rc| магнитное поле катушек Hc на границе плазмы можно оценить как Яс/Я<Дг = а) ~ (/с/ТрХа/^с — а)) или, используя (9.1), — как Яс/Я, ~ (/¡//^(а/^ — а))(а/2х0). Тогда для тока плазмы /р = 600 кА, тока в катушке /с = 15 кА и х0 порядка нескольких сантиметров получим Яс/Я, ~ 1. Более точные оценки могут существенно уменьшить Яс, но в любом случае отношение Яс/Я, достаточно велико.
Оценим другие безразмерные параметры. Для х0 = 2.5 см, а = 46 см, R0/а = 4, m/n = 12/4, Hz = 25 кг, температуры 50 эВ, плотности 1019 м -3, частоты возмущения
w/(2n) = 10 кГц размер области вычисления в направлении y будет y0 = 2na/(mx0) = 10, безразмерное сопротивление (в случае спицеровской проводимости) — ~ 8 ' 10-4, безразмерный коэффициент ионной теплопроводности х ~ 0.009, vp = wa/(mva) = 1.5 ■ 10-2. В этом случае скорость va ~ 1.6 ■ 107 см/с.
Следовательно, DED TEXTOR соответствует случаю малой величины vp, достаточно высокого сопротивления и большой амплитуды возмущения. Это означает, что можно ожидать быстрого проникновения возмущения магнитного поля за время порядка сотен времен x0/va ~ 1.6 ■ 10-7 с.
Заключение
Ускорение плазмы играет важную роль в процессах проникновения возмущения магнитного поля, имеющего вид бегущей волны. Проникновение поля происходит только тогда, когда скорость плазмы достигает фазовой скорости возмущения.
Процесс проникновения протекает существенно различным образом в зависимости от того, меньше фазовая скорость возмущения альвеновской скорости или больше. В последнем случае течение может быть как устойчивым (ламинарным), так и турбулентным (неустойчивым). Заметим, что ламинарный режим имеет место только в случае большой вязкости плазмы и/или в случае очень маленькой амплитуды возмущения. Практически этот режим не может быть реализован в горячей плазме, но может представлять интерес в случае других проводящих жидкостей или газов.
Проникновение бегущего возмущения магнитного поля в плазму является сугубо нелинейным процессом. Учет сжимаемости плазмы не влияет на процесс проникновения поля DED катушек в плазму.
Результаты данной работы могут быть использованы для исследования не только плазмы, но и любых проводящих жидкостей или газов. В частности, бегущее магнитное поле можно использовать как возможный способ ускорения плазмы. При этом скорость плазмы будет определяться фазовой скоростью, а амплитуда возмущения магнитного поля может быть мала.
Список литературы
[1] Finken K.H., Wolf G.H. Background motivation, concept and scientific aims for building a dynamic ergodic divertor // Fusion Eng. and Design, 1997. Vol. 37. P. 337-340.
[2] Dynamic Ergodic Divertor (DED) for TEXTOR-94. IPP Euroatom-KFA, Juelich 1996 / Ed. by K.H. Finken.
[3] Finken K.H. Studies on a helical near field divertor// Nucl. Fusion. 1997. Vol. 37. P. 583593.
[4] Kaleck A., Hassler M., Evans T. Ergodization of the magnetic field at the plasma edge by the dynamic ergodic divertor of TEXTOR-94 // Fusion Eng. and Design. 1997. Vol. 37. P. 353-378.
[5] Abdullaev S.S., Eich Th., Finken K.H. Fractal structure of the magnetic field in the laminar zone of the dynamic ergodic divertor of the torus experiment for technology-oriented research (TEXTOR-94) // Phys. Plasmas. 2001. Vol. 8. P. 2739-2749.
[6] Abdullaev S.S., Finken K.H., Kaleck A., Spatschek К.Н. Twist mapping for the dynamics of magnetic field lines in a tokamak ergodic divertor // Phys. Plasmas. 1998. Vol. 5. P. 196-210.
[7] Abdullaev S.S., Finken K.H., Kaleck A., Spatschek К.Н. Asymptotic and mapping methods in study of ergodic divertor magnetic field in a toroidal system // Phys. Plasmas. 1999. Vol. 6. P. 153-174.
[8] Kobayashi M., Sewell G., Finken K.H. et al. Modeling analysis of the transport properties in TEXTOR-DED laminar zone with a finite element code // Contrib. Plasma Phys. 2002. Vol. 42. P. 163-168.
[9] Abdullaev S.S., Spatschek K.H. Rescaling invariance and anamalous transport in a stochastic layer // Phys. Rev. E. 1999. Vol. 60. P. R6287-R6290.
[10] Faulconer D.W., Koch R. Penetration of the rotating magnetic field into the plasma // Fusion Eng. and Design. 1997. Vol. 37. P. 399-409.
[11] Jensen T.H. Model for plasma response for a finite relative velocity between plasma and a rotating magnetic perturbation // Fusion Eng. and Design. 1997. Vol. 37. P. 437-441.
[12] Finken K.H., Abdullaev S.S., Eich T. et al. Plasma rotation induced by the Dynamic Ergodic Divertor // Nucl. Fusion. 2001. Vol. 41. P. 503-511.
[13] Finken K.H. Perturbation field penetration into the TEXTOR tokamak and the resulting torque // Nucl. Fusion. 1999. Vol. 39. P. 707-723.
[14] Rem J., Schep T.J. Nonlinear dynamics in forced reconnection // Plasma Phys. Controll. Fusion. 1998. Vol. 40. P. 139-145.
[15] Дудникова Г.И., Жуков В.П., Фукс Г. Поведение плотности и полного давления в процессе вынужденного и спонтанного пересоединения // ПМТФ. 1999. Т. 40, № 4. С. 11-15.
Поступила в редакцию 12 апреля 2003 г., в переработанном виде —18 июня 2003 г.
