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

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

Том 9, № 6, 2004

ОБ АЛГОРИТМЕ РАСЧЕТА НАГРЕВА ПЛАЗМЫ ЭЛЕКТРОННЫМ ПУЧКОМ*

В. М. КОВЕНЯ

Институт вычислительных технологий СО РАН, Новосибирск, Россия

e-mail: kovenya@ict.nsc.ru

Т. В. КозлинскАЯ Новосибирский государственный университет, Россия

Е^исш^! difference scheme of predictor-corrector type is constructed for a numerical solution of plasma physics equations. The numerical calculations of plasma heating in a modelling facility are carried out.

Введение

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

Во второй части работы построена схема типа предиктор-корректор для численного решения уравнения физики плазмы в двухтемпературном приближении [5, 6]. Как ив [2], она реализуется на дробных шагах скалярными прогонками, а на этапе корректор явно и аппроксимирует исходные уравнения со вторым порядком. На основе проведенных расчетов нагрева и движения плазмы в установке с единой магнитной ячейкой получены оценки нагрева плазмы и распределения полей, подобные результатам экспериментов на установке ГОЛ-3 [5, 6].

* Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований (грант № 02-01-01029) и Министерства образования РФ (грант № Е 02-1.0-25).

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

1. Физико-математическая модель

Нагрев и движение плазмы в магнитной системе экспериментально исследовались на установке ГОЛ-3 [5, 6]. Быстрый нагрев плазмы осуществляется релятивистским электронным пучком высокой мощности. В гофрированном магнитном поле это приводит к появлению периодического продольного изменения плазменного давления. Градиент давления инициирует движение плазмы по направлению к середине каждой магнитной ячейки. Эксперименты показали, что происходит быстрая термализация направленной энергии движения плазмы. Поскольку на начальной стадии нагрева температура ионов плазмы мала и длина свободного пробега ионов много меньше длины одной ячейки, численное моделирование динамики двухкомпонентной плазмы может быть проведено в гидродинамическом приближении. Динамика плазмы описывается уравнениями неразрывности и движения в приближении замагниченной плазмы и уравнениями энергии для электронов и ионов, которые в рамках одномерного движения могут быть представлены в безразмерном виде в следующей форме [5, 6]:

др др д (V

т + ^ + Врд~г Кб дv дv А0 др о дЬ дг р дг

д'Рг + v др- + а д

дЬ дг г дг

др др д (V

т + ^ + адг Ь

о,

д_

дг

~ дТ-

6д- ) + А1 (Р - 2Р-)

(1)

д (? дТ\ ^

дг 6 + ^

где р- = рТ, ре = рТе, р = р- + ре; 7 = 5/3; I = 7 - 1; а = 7Вр, а- = 6 = К^Т^"2/1, £е = Ке^еТе/2/(¡() — коэффициенты теплопроводности как функции температур Т и Те; А0 = с0Мр/(^М), А1 = 1.3014 ■ 10-3^е2нМр(Те3/2М1) — безразмерные

параметры; с0 — ионно-звуковая скорость волн в плазме; К0 = ((0/1 +

_ дг{(6 6е) д!

описывает изменение энергии в компонентах плазмы, а 6(о — нагрев плазмы электронами пучка. Вывод уравнений и механизм взаимодействия плазмы приведены в [6]. Представим систему уравнений (1) в векторном виде:

I

К

(2)

где

П

д_

дг

0 0 о

Вр д

I

д. 1

дгВ

р

V р-

р

К

/0 \

0 0

\Ко/

а

а

дг

д. 1

- дгВ

д. 1

дг В

0 0

д

V---

0

0

Ао д

р дг

д 6 д 1

+ 2А1 -А1 д

дг дг дг р

^ —6 — 1

дг дг е дг р /

Решение задачи (2) будем искать в области 0 < г < Ь, г > 0. Граничные условия на концах системы определяются свободным вытеканием плазмы через торцевые магнитные пробки и задаются в виде [5, 6]

-РТ ■ -и

-РТ*и=о,ь = 0, —и=о,ь = 0,Те,г|^=с,ь = шах(Тое;г, 0.05Те"Г). (3)

-г -г

Начальные условия соответствуют заполнению системы однородным потоком с параметрами Т0е,г = 1 эВ, п0 = 1015 см-3, -0 = 0.

2. Разностная схема для уравнений газовой динамики

Численный метод решения системы уравнений (2) опишем вначале для уравнений газовой динамики в квазиодномерном приближении

-г

+ П/ = б,

(4)

где

/

Р -

Р

( V-

-г

П

V

в -1

-г

в -1

0

1 -

Р -г

-г

/

1

Система уравнений (4) — частный случай уравнений (2) при рг = р, А0 = 1, Б = 0, ; £е = 0, а под В понимается площадь поперечного сечения канала. Очевидно, при В = уравнения (4) есть система одномерных нестационарных уравнений газовой динамики.

Численное решение задачи (4) будем отыскивать в области П= {г0 < г < г1; 0 < г < Т0}. Введем в П расчетную сетку с постоянными шагами к и т. Аппроксимируем оператор П разностным оператором П:

П

( -пЛ ВрпЛ— 0

В РП Л

Рп

1

0 7 ВрпЛ— -пЛ Б '

0

-пл

Здесь Л

Л-, Л

+ ,

если - > 0, если - < 0,

Л

Л

Л-,

+,

в

если - > 0, если - < 0,

а Л±/з = ±(Л'±1 - /з)/к — ап-

проксимация первой производной с учетом знака скорости -п с порядком О (к). Следуя работам [2, 3], представим оператор П в виде суммы операторов П = П1^ + П2^, где

П

1к =

( -пЛ ВрпЛ1 0 \

в

0 -пЛ 0

0 0 0

П

2Н

00 00

0

1

\

V

0 7ВрпЛ

1

в

Л

Р

-пЛ

(5)

Разностная схема типа предиктор-корректор

т+1/4 _ /п /п+1/2 _ £п+1/4

£п+1/4 _ £п £п+1/2 _ £п+1/4

/-^ + Пш/п+1/4 = 0, /-^-+ ^2^/п+1/2 = 0; (6)

та та

/п+1 _ /п

/-^ + П /п+1/2 = о (7)

т

,к\

аппроксимирует исходные уравнения с порядком 0(тт + тк + кк), где т =2 при а = 0.5 + 0(т) и т =1 при а = 0.5, а к — порядок аппроксимации оператора П на этапе корректор. Для повышения точности расчета оператор П на этапе корректор аппроксимируем симметричным оператором со вторым порядком. Если уравнения газовой динамики на этом этапе аппроксимировать в дивергентной форме по схеме

ип+1 _ ттп

--— + ЛWn+1/2 = ^п+1/2, (8)

т

то разностная схема (6), (8) аппроксимирует исходные уравнения с порядком 0(т2 + к2) при а = 0.5 + 0(т). Здесь

/ Вру \ ( 0 \

и = , W = В(ру2 + р) , Г = рЛВ

\ Ву(Е + р)

В 0

В силу симметричной аппроксимации вектора W схема немонотонна и для устранения осцилляций в нее на этапе корректор вводится сглаживающий оператор второго порядка малости подобно [4]:

. г (а/)?+1 _ (а/)?-1 к , ,

Ла/ = -- _ 2 (Ь^+1/2Л+/^ _ ^-1/2Л-Л) , (9)

где Ьу±1/2 = (Ьу±1 + Ьу) /2, Ьу = е2 | /, d = \aj+1 _ ау | + _ ау-1| , а е = 0, если д = 0, иначе е = (ау+1 _ 2ау + ау-1)/д.

Разностные схемы (6), (7) или (6), (8) на дробных шагах реализуются трехточечными скалярными прогонками, а на этапе корректор явно [2] и в линейном приближении безусловно устойчивы при а > 0.5. Тестирование схем проведено на решении двух задач: задаче о распаде произвольного разрыва и течении в канале переменного сечения. Для первой задачи начальные и краевые условия задавались в виде

р(х, 0) = 1, и(х, 0) = 0, р(х, 0) = 1 при 4 < х < 0,

р(х, 0) = 0.125, и(х, 0) = 0, р(х, 0) = 0.1 при 0 <х < 6,

р(_4, г) = 1, и(_4, г) = 0, р(_4, г) = 1, р(б, г) = 0.125, и(б, г) = 0, р(б, г) = 0.1

На рис. 1 приведены результаты расчетов по схеме второго порядка со сглаживанием на момент времени г = 2.5 на сетке, содержащей 200 узлов по г, при а = 0.505 со сглаживанием решения в правой части схемы (7) или (8).

Видно, что ударная волна размазывается на 2-3 узла, а контактный разрыв — на ~ 6 узлов. Схема правильно передает скорость движения ударной волны и волны разрежения. Расчеты проведены при числе Куранта к =1.

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

■3 -2-10 1 2 3 '4

Рис. 1.

Density

1.000

0.374 4.3S7

1.023 3.000

2.020

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.G 0.Э

1 Velocity

'0.1 0.2 '0.3 0.4 0.5 '0.6 '0.7 '0.8 '0.9

Pressure

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.Э

Рис. 2.

при £ = 0 в канале переменного сечения, задаваемого по формуле В (х) = 0.5 + 0.5(1 — 2х)2, 0 < х < 1, задавались значения функций на входе (х = 0) и выходе (х =1) из канала, удовлетворяющие точному решению

Р u

Р

x=0

1

1.0237

8

Р

u Р

1 ' x=l

0.8835 1.1586 7.0315

которое содержит разрыв газодинамических величин. На рис. 2 приведены точное (сплошная линия) и численное (значки (о)) решения после установления по схеме второго порядка (6), (8) со сглаживанием на расчетной сетке, содержащей 200 узлов.

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

3. Разностная схема для уравнений физики плазмы

Подобно (3) представим разностный оператор П в виде суммы П^ и П2^, где

П

ш

( уЛ ВрЛ1 0 0 ^ В

0 уЛ 0 0

0 0 0 0

V 0 0 0 0 /

п

0 0 0

0 0 0

0 А1 а*ЛВ аЛ В уЛ

0 0

1

к.

р

0

Ао

Р

(10)

Л

уЛ - Л&Л1 р

а ЛаЛ — аппроксимация вторых производных симметричным оператором. С учетом введенных обозначений разностная схема

/п+1/4 _ f п /п+1/2 _ f п+1/4

та

уп уп+1/2 _ уп+1/4

у + Пш/п+1/4 = 0, у-у-+ П2^/п+1/2 = 0,

та

/ п+1 _ /п

(11)

т

+ п /п+1/2 = ^п+1/2

аппроксимирует систему уравнений (2) с порядком 0(тт + тЛ + Л2), где т = 2 при а = 0.5 + 0(т). Здесь П^ — аппроксимация уравнений симметричными разностями. Для устранения осцилляций в него подобно (9) добавлены сглаживающие члены. В линейном приближении при ^ = 0 схема (11) безусловно устойчива. На дробных шагах она реализуется скалярными трехточечными прогонками, а на этапе корректор — явно. Действительно, на первом дробном шаге система разностных уравнений

рп+1/4 _ рп 1

Р-^ + упЛрп+1/4 + ВрпЛ-уп+1/4 = 0,

та В

п+1/4 _ _п

+ упЛуп+1/4 = 0,

та

рп+1/4 _ рп рп+1/4 _ рп "г_1 г _ 0 __" _ 0

та

та

может быть решена независимо для каждой компоненты вектора /п+1/4. На втором дробном шаге разностные уравнения

рп+1/2 _ рп+1/4

та

0,

vn+1/2 _ vn+1/4 A

^ ^ - + ^ Äpn+1/2 = 0,

та

n+1/2 n+1/4

p _ Рг - + (vn+1/2/B) + vnAp™+1/2 = легл ( рГ'7рп ) + A1 (2pi

та

+ 1/2 _ pn+1/2

Р

pn+1/2 _ рП+1/4

та

+ anA (vn+1/2/B) + vnApn+1/2 = леегл (pn+1/2/pn)

решаются в такой последовательности: исключая ^га+1/2 из уравнения для полного давления, получим разностное уравнение относительно рп+1/2 в виде

I + та^гл _ т2а2агл—— л _ тал&л

Bpn

1

pn

Р

n+1/2 = pn+1/4 _ тааnЛ■

n+1/2

B '

Его решение может быть получено скалярной трехточечной прогонкой, после чего явно вычисляется <уга+1/2. Новое значение рп+1/2 определяется скалярной прогонкой из уравнения энергии для р при уже известных значениях рп+1/2 и ^га+1/2. Значение плотности на слое п +1/2 переносится с предыдущего шага.

4. Моделирование нагрева плазмы

С целью проверки механизма нагрева плазмы в рамках предложенной модели (2) проведены расчеты в модельной установке ГОЛ-3 с одиночной магнитной ячейкой. В начальный момент времени при t = 0 в канале L = 12.3 задавались невозмущенные значения искомых величин р = 1, v = 0, Ti = Te = 0.001. Внешнее магнитное поле задавалось по формуле B = 2 + cos(2(x _ 1.5п)) для 1.5п < x < 2.5п, иначе B = 3. Система предполагалась открытой, допускающей свободное вытекание плазмы на границах канала и удовлетворяющей краевым условиям (3). Численное решение уравнений (2) находилось по схеме предиктор-корректор (11) с порядком 0(т2 + h2) при а = 0.505 и со сглаживанием решения на этапе корректор по формуле (9). Расчетная сетка содержала 400 узлов по координате z. Измельчение шагов сетки по пространству в 2-4 раза не приводило к изменению результатов расчетов. Выше отмечалось, что разностная схема (11) безусловно устойчива при отсутствии внешних сил. Однако для ненулевой правой части F = 0 (при наличии внешних источников Q), значение которой вычисляется явно на этапе корректор, разностная схема может терять свойство безусловной устойчивости. В этом случае временной параметр т определяется из устойчивости схемы по правой части. За счет влияния F электронная температура резко возрастает (~ 1000) и источниковый член в уравнении энергии становится определяющим, значительно превосходящим конвективные и другие члены уравнения. Численные расчеты проведены при значении k ~ (0.01 _ 0.1), допускаемом устойчивостью схемы. На рис. 3, 4 представлены распределения скорости, плотности, электронной и полной температуры в моменты времени t = 2 и t = 4 соответственно.

В соответствии с теорией и экспериментальными данными, полученными в установке ГОЛ-3 (см. [5, 6]), нагрев плазмы электронным пучком происходит неравномерно из-за неравномерности распределения магнитного поля: температура значительно повышается по краям пробкотрона и в меньшей степени в центре ямы. Неоднородность давления приводит к ускорению плазмы к центру ячейки. При t > 4 два потока сталкиваются и, как следует из эксперимента и теории, они проходят друг через друга.

Рис. 3.

Рис. 4.

В рамках односкоростной модели (1) проведенные численные расчеты дали следующие результаты. При решении задачи за счет нагрева плазмы электронным пучком происходит существенный рост электронной температуры (~ 1000 раз), причем ее возрастание максимально вне магнитной ямы (см. рис. 3). Ионная же температура возрастает почти в 10 раз и достигает максимального значения на границах магнитной ямы. За счет ускорения плаз-

мы к центру магнитной ямы плотность возрастает в зоне ямы до значений порядка 2.8 от начального и уменьшается вне ямы. Полное давление целиком определяется электронным давлением (температурой), оно возрастает на концах системы больше, чем в центре ямы. Вплоть до времени t = 4 численное решение качественно и количественно близко к полученному в эксперименте [5, 6]. Однако для t > 4 из расчетов следует, что два потока плазмы сталкиваются в центре ямы, что сопровождается резким возрастанием температур и давлений. Для таких времен расчета односкоростная модель перестает быть справедливой, и для расчета таких течений необходимо использовать многопотоковые модели.

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

[1] ковеня в.м. Разностные методы решения многомерных задач: Курс лекций. Новосибирск: Изд-во НГУ, 2004. 146 с.

[2] Ковеня в.м., Лебедев а.с. Модификация метода расщепления для построения экономичных разностных схем // Журн. вычисл. математики и мат. физики. 1994. Т. 34, № 6. С. 886-897.

[3] КозлинскАЯ Т.В. Схема предиктор-корректор для численного решения уравнений газовой динамики // Прогр. и тез. докл. IV Всерос. конф. молодых ученых по мат. моделированию и информ. технологиям, 30.10 — 5.11 2003 г., Красноярск. С. 29.

[4] Ковеня в.м., Лебедев а.с. Численное исследование отрывного течения в ближнем следе. Новосибирск, 1987 (Препр. / АН СССР. Сиб. отд-ние; № 14-87. 48 с.)

[5] Астрелин В.Т., БурдАков А.В. Численное моделирование коллективного ускорения ионов плазмы в ячейках многопробочной ловушки ГОЛ-3 // Тез. докл. XXXI конф. по физике плазмы и УТС. Звенигород, 2004. С. 102.

[6] Arzhannikov a.v., ästrelin a.v., bürdakov a.v. et al. Recent results on the GOL-3-II facility // Trans. of Fusion Technology. 1999. Vol. 35, N 1T. P. 112-120.

Поступила в редакцию 29 сентября 2004 г-