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

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

Том 9, № 5, 2004

ПОЛНОСТЬЮ КОНСЕРВАТИВНАЯ РАЗНОСТНАЯ СХЕМА ДЛЯ РАСЧЕТОВ ТЕРМОЯДЕРНОГО РЕАКТОРА НА ОСНОВЕ АСИММЕТРИЧНОЙ ЦЕНТРОБЕЖНОЙ

ЛОВУШКИ

В. И. Волосов, С. С. Клепов Институт ядерной физики им. Г.И. Будкера СО РАН,

Новосибирск, Россия e-mail: V.I.Volosov@inp.nsk.su, kljopov@ngs.ru

A numerical method for solution of the two-dimensional problem for the Fokker —

Plank equation describing the behavior of collisional plasma is considered. The absolutely conservative difference scheme is realized, which allows one to find solutions at a rather large ratio of the plasma confinement time and the Coulomb scattering time. Numerical computations for the modernized open magnetic trap of the ACT type have been carried out.

Использование полностью консервативных разностных схем (ПКР-схем) для решения уравнения Фоккера — Планка необходимо при решении задачи о нахождении параметров плазмы в модифицированной открытой магнитной термоядерной ловушке типа асимметричной центробежной (АЦЛ) [1-5], поскольку время удержания в этих ловушках много больше характерного времени жизни частиц. В этих задачах требуется решать двумерное в пространстве скоростей уравнение Фоккера — Планка (УФП). Как было показано в [6], при достаточно большом отношении величины потенциального барьера, удерживающего плазму в ловушке, к температуре плазмы невозможно получить корректные решения этих задач, используя разностные схемы, в которых не выполняется принцип полной консервативности, поскольку аппроксимационные ошибки приводят к появлению фиктивных источников энергии и, как следствие, охлаждению плазмы в расчетах.

В данной работе произведен расчет АЦЛ [2, 3] для безнейтронной реакции puB с использованием ПКР-схемы аналогично тому, как это было проведено в работе М. Пеккера

[7]. В отличие от ранних расчетов, выполненных в [7, 8], в данном варианте ловушки время жизни отдельных компонентов плазмы на два порядка (по оценке) превышает характерные (кулоновские) времена жизни этих частиц в такой системе, что требует в расчетах тщательного соблюдения консервативности схемы. Другая особенность рассматриваемой задачи — необходимость детального анализа двумерной функции распределения каждого компонента плазмы. Эти расчеты ранее не проводились, так как в [7] использовались только интегральные характеристики плазмы. Плазма, состоящая из DT ионов, рассматривалась как среда, состоящая из ионов с одинаковой массой (2.47Mp) и энергией, определяемой одним параметром T (максвелловское распределение ионов по энергиям). В данных

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

расчетах каждый компонент плазмы (Р,В,е) рассматривался отдельно. Находилось распределение плотности ионов в фазовом пространстве Г(Е,в). Для каждого компонента использовалась сетка 50 х 30 ячеек с разбиением по амплитуде Г не менее 50 точек. Задача решалась до установления устойчивого стационарного состояния, т. е. использовалось несколько тысяч итераций по всем параметрам задачи. Число итераций зависит от е, поэтому потребовалась некоторая переработка схемы расчетов. В отличие от ранних работ, в расчете были использованы дополнительные условия для определения потерь на тормозное излучение. Также необходимым является и изменение функции источника исходя из требований по налагаемым условиям удержания частиц в ловушке, уходу их из системы и потерям по энергии.

Таким образом, постановка задачи отличается от постановки задачи в [7, 8]. Вместо однокомпонентной рассчитывается двухкомпонентная плазма. Появились новые члены, описывающие специальные условия потерь энергии, изменились члены, описывающие приход энергии в систему вследствие обмена энергией между разными типами частиц. В данной работе реализована полностью консервативная разностная схема с учетом перечисленных выше специальных требований, вытекающих из условий данной задачи. Приведены результаты численных расчетов. Показана применимость ПКР-схем для расчета АЦЛ.

Двумерное уравнение Фоккера — Планка, описывающее временную эволюцию функции распределения частиц в пространстве скоростей в безразмерных координатах, имеет вид [8, 9]

f = I--1-.-; ч ~ д

дt m2 дv

j,k

k

дНа +1 д f 02Ga f

дvk f +2 дvj \дvkдvjfa

+ Sa. (1)

Здесь О и Н — потенциалы Розенблюта; а — тип частицы; та, га — масса и заряд частицы; 1а — функция распределения частиц по скоростям; Ба — функция источника:

2

На ^ ^ 1пЛа^ 1 + ^;

X-^ -

&а / у 2 9 в 1п Ла/3;

в га

1пЛад — кулоновский логарифм; кр = 2 / |^ — ^'І-1^3^'; 9р = / |^ — ,У/|^3'У/.

В магнитных ловушках функция распределения частиц не зависит от азимутального угла и симметрична относительно плоскости в = 0. Граничные условия имеют вид

ди(0у) 0 0

= 0 при V = 0,

дv

ди{Ov)

О при 0 = О и п/2,

дв

fa = 0 на (RH — 1)v2 sin2 в — v2 cos2 в + Ua = 0,

где Rh — пробочное отношение [2, 3]; Ua — потенциальный барьер, удерживающий частицы типа а. Явный вид этого потенциала и функции источника определяются типом магнитной ловушки и методом инжекции.

Как было показано в [7, 8], построенная система разностных уравнений на равномерной сетке Vj = ]5У,ви = к8в, аппроксимирующих УФП (1), неполностью консервативна,

так как не имеет разностного аналога закона сохранения энергии. Энергия сохраняется с точностью 6У2. Нетрудно видеть, что аппроксимационные ошибки по V являются фиктивными источниками энергии, которые приводят к охлаждению плазмы.

Для построения ПКР проанализируем УФП в прямоугольных координатах, как это было сделано в [7]. Выпишем в явном виде производные от функций О и Н:

д 29а дviдvj

[ , Г Л (у — у')2 — ('V — ^)(^' - ^Ь /

!а(У)------------,--------

I \у — У/\3

£ - —2 /

Разностное уравнение, прямо соответствующее форме (1), имеет вид

i,j

дНа £ , 1 ( д2Са £

д Iа + о I а_. а_. !о

2 \ дvi дvj

где (')г означает разностную производную по скорости vx, например:

(п)х — (п(Ух + Н,Уу) — п(Ух — Н,Уу))/(2к),

(')t — разностную производную по времени: (')t — ('(£ + т) — '(^))/т. Будет ли эта схема консервативной, зависит от того, каким образом определить производные дНа/дУг и д2Са/дь1дь^. В самом деле, выражения для изменения энергии и импульса по времени:

дНа

дvi

+ - АСа fак3

fviк3

Е

дН

дv.

а fаh3

Равенство нулю правых частей можно обеспечить, если вычислять производные дHa/дvi и д2Ga/дviдvj как аппроксимацию интегралов.

Выпишем разностную схему, которая использовалась для расчетов параметров плазмы в АЦЛ. Так как задача является аксиально-симметричной, рассмотрение ведется в сферических переменных V, в, р (зависимость от р отсутствует). Двумерное УФП в этих переменных для каждого из компонентов плазмы перепишется в виде

ди

дЬ

I д_

УМ дУ

~ 1 д — УНа]^а + 2у ду (Аа/а)

+

+

1д

яд! + г f + п дд£а

Яа дв + Гаіа + дУ

+ Ба-

(2)

вт в дв

Функции определяются двумерными интегралами по V и в. Аппроксимируя эти интегралы двумерными суммами, несложно построить разностную схему, соответствующую (2), которая сохраняет энергию и число частиц. Функции На, На, да (а также производные да по V) можно вычислять, используя разложение в ряды по полиномам Лежандра:

ка(У,в) = ’£

4п

1=0

21 + 1

(У-(1+1)Ыа1 + У 1МаІ)Рі(еС8 в),

j

г

ha(V,0) = Y,4n{VlMai - V-(1+l)Nal)Pi(cos

l=0

,(v.o) = '£

4n

l=0

2l + 1

1

2l + 3

(V-(1+l)Eal + V l+2Mal)-

1

2l - 1

(V1-1 Nal + V1 Ral)

Pl (cos 9),

d9a

£■

4n

dV ^ 2l + 1

l=0

1

2l + 3

((l + 2)Vl+1)Mal - (l + 1)V-(l+2)Eal)-

1

2l - 1

(lVl-1Ral - (l - 1)V-lNal)

Pl (cos 9),

где

d2ga

L . 4n

dV2 ^ 2l + 1

l=0

(l +22 + 3+ 1) (VlMal + V-(l+3)Nal)-

(l - 1)l l-

2l - 1

(Vl-2Ral + V-(l+1)Nal)

Pl (cos 9),

V

Mal(V) = fal(V’)V’1-ldV’; Nal(V) = fal(V')V,2+ldV’;

V

V

Ral(V) = fal(V’) V’3—ldV’; Eal(V) = fal(V’)V’4+l dV’;

V

(3)

2l + 1

fal (V) = -^ У fa(V, 9) sin 9d9.

0

Верхний предел L в суммах, вообще говоря, равен то. При подстановке функций (3) с конечным L в (2) получится приближенное дифференциальное уравнение. Существенным для нас является то обстоятельство, что это приближенное уравнение сохраняет энергию при любом L (сохранение числа частиц для этого уравнения очевидно и обеспечивается его дивергентной формой). В самом деле,

СЮ (/ 1

2n^2aap[ 1 + ma

a,e

4n

mp J -f—' 2l + 1

fa(V,9)Pl(cos 9)d(cos 9) I x

l=0

0 -1

x

l+

m„

ma + me

V2+lMet - f 1 +—me—) V1-N

ma + me

*@l

dV

g

^ 4ne4z2zp(ma + me) I ^ 16п2

ав mame x I to(2/ +1)2

_ ( 11 ~---- ) falNei

ma + mp/

I + m ) fei Nal (V )_ ma + me/

V1-ldV = О. (4)

Итак, исходным для построения разностной схемы берется уравнение (2) с коэффициентами переноса, вычисляемыми с помощью (3). Граничные условия для уравнения определяются из соотношений

д!а{0у) 0 0

= 0 при V = 0,

дv

dfa(0v)

О при в = О и п/2,

дв

fa = 0 на (RH — 1)v2 sin2 в — v2 cos2 в + Ua = 0.

Для решения этой системы вводится равномерная разностная сетка по V, в ,t (область решения 0 < V < V0, 0 < в < п, t > 0): Vi+1/2 = (i + 1/2)hV, i = 0,1, 2,..., I — 1; в^+\/2 = (j + 1/2)hg, j = 0,1, 2,..., 2J — 1, hv ,h$ — шаг сетки по модулю скорости и углу соответственно; tn = пт, п = 0,1, 2... — шаг по времени, который является итерационным параметром. Метод решения (метод продольно-поперечной прогонки) выбран как один из наиболее простых в реализации методов, обеспечивающих приемлемые скорость нахождения и устойчивость решения. Введем вспомогательные сеточные функции:

Ra,i,j+1/2 = — 2 (Vi+1/2H'n,i+1/2,j+1/2 + ^Vi-1/2Ha,i-1/2,j+1/2) +

І

'/ , \7 (Aa,i+1/2,j+1/2fa,i+1/2,j+1/2

(Xi+1/2 + Xi-1/2)hv

_Aa,i-1/2,j+1/2fa,i-1/2,j+1/2) при * > 1, Ra,0,j+1 /2 — 0,

Pa,i+1/2,j = 2h (Ba,i+1/2,j+1/2 + Ba,i+1/2,j-1/2)x (5)

I'd

x(fZ

і+1/2,7+1/2 fа л+т^-т)+

I 1 ( c n fn I C'n f n )

+ 2 (Ca,i+1/2,j+1/2 Ja,i+1/2,j+1/2 + Ca,i+1/2,j — 1/2Ja,i+1/2,j—1/2) ,

8h (Da,i+1/2,j+1/2 + Dn,i+1/2,j-1/2) x

iV

i+3/2,j+1/2 + J a,i+3/2,j-1/2 _ J a,i+1/2,j+1/2 _ f a,i-1/2,j-1/2) при 2J _ 1 > j > 1,

x (fan ,i+3/2,j+1/2 I fa,i+3/2,j —1/2 fо:,і+1|2,7+1|2 fa,i—1|2,j—1|2)

nn Pa ,i+1/2,0 = Pa ,i+1/2,2J — О.

Искомая разностная схема в обозначениях (5) имеет вид

fn+1 f n

J a,i+1/2,j+1/2 J a,i+1/2,j+1/2 _

V 2

Уг+1/2

тэи _ ТЭП Г)П _ Г)П

Ла,г+1,.7+1/2 Ла,і,.7+1/2 + Ра,г+1/2,і+1 Ра,г+1/2,і

Ну вІП в^+1/2Не

Закон сохранения числа частиц выполняется автоматически:

(6)

д I -1 2 J -1 г П+1 _ г П т V 4

^Па _ ^ ^ га,і+1/2,і+1/2 г а,і+1/2,і+1/2 ,ПаУі+1/2 „ а и и _ п

------ — 2_^ 2-^ 2---------------------------- Х вІП "з+1/2НУ НЄ — 0Т-----------------------------------------------------------------Т-2

і=0 J=0

Проверим закон сохранения энергии:

АЕ 1 -1 2J-1

Т

а,і+1/2,7+1/2 Г®-,і+1/2,з+1/2 Х

а і=0 і=0

Ша^+1/2 . п , ,

Х------2--- вІП ^+1/2 Ну Не —

I-1 2 J-1 т V 2

та ^і+1/2

Г) Hа,г+1/2,j+1/2fа,г+!/2л+1/2 вj+1/2hV .

а )=0 j=0

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

Расчеты с использованием вышеприведенной схемы позволяют получить в значительной степени более детальную информацию о функции распределения, чем проведенные вычисления в [7]. Для демонстрации работоспособности метода приведем результаты некоторых расчетов, выполненных на основе разностных уравнений (6), для р11В реакции в АЦЛ системе для некоторых параметров установки. Для проверки и отладки программы расчета проведен ряд некоторых простых математических тестов с измельчением сетки, показавших корректность и верность решения с математической точки зрения.

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

Рис. 1. Вид функции распределения протонов в скоростном пространстве для разных углов 9.

о 0.5 1 1.5 V/V0

Рис. 2. Вид функции распределения протонов в энергетическом пространстве в зависимости от скорости ионов.

/

4

3

2

1

0

Л0 0.5 1 1.5 V/V0

Рис. 3. Функция распределения ионов 11B в зависимости от поперечной скорости при в = 0.

решения уравнения Фоккера — Планка по вышеописанной схеме выбран метод продольнопоперечной прогонки как один из наиболее простых в реализации и оптимальных по скорости установления решения.

Авторы благодарят Е.В. Бобуха за существенный вклад и помощь в работе над программой расчетов.

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

[1] Будкер Г.И. Термоядерные реакции в системе с магнитными пробками. К вопросу о непосредственном преобразовании ядерной энергии в электрическую // Физика плазмы и проблемы управляемых термоядерных реакций. Т. 3. М.: Изд-во АН СССР, 1958.

[2] Волосов В.И. Асимметричная центробежная ловушка // Физика плазмы. 1997. Т. 23, № 9. С. 811-815.

[3] VOLOSOV V.I. Some problems of magnetic fusion reactor and asymmetrical centrifugal magnetic confinement device // Transaction of Fusion Technology. Jan. 1999. Vol. 35. P. 258-262.

[4] Димов Г.И., ЗАКАйДАКОВ В.В., Кишиневский М.Е. Термоядерная ловушка с двумя пробками // Физика плазмы. 1976. Т. 2, № 4. С. 597-610.

[5] Lehnert B. Rotating plasmas // Nucl. Phys. 1971. Vol. 11. P. 485-533.

[6] Волосов В.И., Пеккер М.С. О точности численных расчетов потерь плазмы из открытых магнитных ловушек // Числ. методы механики сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ВЦ; ИТПМ. 1979. Т. 10, вып. 1. С. 45-56.

[7] Пеккер М.С. Продольное удержание в магнитной ловушке с вращающейся плазмой (численный расчет). Дис. ... канд. физ.-мат. наук. Новосибирск, 1983.

[8] Волосов В.И., Пеккер М.С. О численных методах решения двумерной задачи для уравнения Фоккера — Планка // Журн. вычисл. математики и мат. физики. 1980. Т. 20. С. 1341-1346.

[9] Rosenbluth M.N., MacDonald W.M., Judd D.L. Fokker — Planck equation for an inverse-square force // Phys. Rev. 1957. Vol. 107, N 1. P. 1-6.

[10] Пеккер М.С. Полностью консервативная разностная схема для двумерного уравнения Фоккера — Планка. Новосибирск, 1980 (Препр. ИЯФ; № 80-38). 11 с.

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