Цветков Применение численных методов 2007

	D		c2		δ 2
где		=		=	sc , δ	sc	− так называемая глубина скин-слоя
	m		4πσ		τsc
(глубина,		на которую проникает магнитное поле в плазму), τsc −

характерное «скиновое» время (время проникновения).

6.2. Применение МГД приближения для расчета ускорения плазмы в коаксиальном плазменном ускорителе

С учетом соотношения divρv = v gradρ + ρ divv , с учетом

того, что в цилиндрической системе divv =

1

∂

(rv

) + 1 ∂vθ

+ ∂vz ,

r ∂r

r

r ∂θ

∂z

для аксиально-симметричной системы

B = Bθ (рис. 21) уравнение

непрерывности примет вид:

∂vr

vr

∂vz ) = 0 .

∂ρ

+v ∂ρ +v

∂ρ + ρ(

+

+

∂t

r

r ∂r

z ∂z

∂r

∂z

С учетом того, что в уравнении движения полная производная по времени имеет вид dtd = ∂∂t +v ∂∂r , уравнение движения для ради-

альной составляющей имеет вид:

		∂v		∂v	∂v	∂		B2		B2
	ρ(	r +v		r +v	r ) = −		( p +	θ	) −	θ	,
						∂r		8π		4πr
		∂t	r ∂r		z ∂z

для осевой составляющей:

ρ(∂vz +v		∂vz +v			∂vz ) = −	∂
						∂z
∂t	r	∂r		z	∂z
С учетом того, что			Bθ = div(gradBθ
симметричного поля примет вид:

( p + Bθ2 ) . 8π

) , уравнение для аксиально-

∂Bθ

=

c2

(

1 ∂

(r

∂Bθ ) +

∂2 Bθ ) −

1

(v B ) −

1

(v B ) .

∂t

4πσ

r ∂r

∂r

∂z

∂r

∂z2

r θ

z θ

61

Таким образом, имеем систему уравнений:

∂ρ

+v

∂ρ

+v

∂ρ

+ ρ(

∂vr

+ vr

+

∂vz ) = 0;

∂t

r ∂r

z ∂z

∂r

r

∂z

∂v

∂v

∂v

∂

B2

B2

ρ(

+v

+v

= −

( p

+

) −

r

r

r )

θ

θ

;

∂r

8π

4πr

∂t

r ∂r

z ∂z

∂v

∂v

∂v

∂

B2

ρ(

+vr

+vz

z ) = −

( p

+

);

z

z

θ

∂z

8π

∂t

∂r

∂z

ρ

p =

kT;

mi

∂Bθ

=

c2

(

1

∂

(r

∂Bθ

) +

∂2 Bθ

)

−

1

(v B

)

−

1

(v B

).

∂z2

∂t

4πσ

r

∂r

∂r

∂r

∂z

r

θ

z θ

Расчет каждого из уравнений системы является задачей Коши, тре-
L0	R0		бующей зада-
			ния		началь-
			ных		условий.
			Считаем,				что
	Bθ	r	в	начальный
c0			момент				вре-
		z	мени		плазма
				заполняет
		θ	слой		толщи-
			ны		l,	тогда
				начальные
			условия				для
	r2		0 ≤ z ≤ l				и
				r1 ≤ r ≤ r2
	r1	z	запишем				в
			виде:
	z0	Zк
	l
Рис. 21. Схема коаксиального ускорителя

62

vz = v0 ;vr = 0;Bθ = 0;T =T0 ;

ρ = ρ ,0

для l ≤ z ≤ zk и r1 ≤ r ≤ r2 в виде:

vz = 0;vr = 0;Bθ = 0;T = 0;

ρ = 0.

Численно эти уравнения можно рассчитывать в конечных разностях, введя сетку (см. рис. 21). Зададим граничные условия для

r = r1;			r = r2 ;
	< z	≤ z0	и для	< z	≤ zk
0			0

для z = 0 в виде:

vr = 0;

∂vz = 0;∂r

в виде:

∂∂ρr = 0;Bθ = 0,

vr			= 0;
v		z	= 0;
						B	(r	, 0,t) r
Bθ (r, 0,t) =						θ	1	1	;
								r
	ρ		= ρ		,
				0

63

cr1

2I (t)

где значение магнитного поля в точке (r1, 0) задается соотношени-

ем Bθ (r1, 0,t) = . Зависимость тока от времени I (t) опреде-

ляется из уравнения электрической цепи

		d
	u =			(L0 I ) + R0 I;
		dt
				1
	u = −				I,
			c
				0

где c0, L0, R0 – емкость, индуктивность и сопротивление цепи, u – напряжение.

6.3.Ускорение плазмы в рельсотроне (одномерный случай)

Водномерном случае система уравнений МГД приближения примет вид:

∂ρ

= −v

∂ρ

− ρ

∂v

;

∂t

∂x

∂x

∂v

∂v

1

∂p

B ∂B

= −v

−

−

;

∂t

∂x

ρ

∂x

4πρ ∂x

∂B

= −v ∂B − B ∂v ;

∂t

∂x

∂x

ρ

p = m kT.

i

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

вых значений в момент времени n+1 на пространственной сетке

1 ≤ j ≤ k :

64

n +1

n

n

t

n

n

n

t

n

n

ρ j

= ρ j

− vi

(ρ j +1 − ρ j −1 ) − ρ j

(v j +1

− v j −1 );

2

x

2 x

n +1

n

n

t

n

n

1

t

n

n

v j

= v j

− v j

(v j +1

− v j −1 ) −

(ρ j +1

− ρ j −1 ) −

2 x

n

2

x

ρ j

n

t

n

n +1

B j

−

(B +

− B

−

);

4πρnj 2

x

j

1

j

1

n +1

n

n

t

n

n

n

t

n

n

B

=

B

− v

(B

− B

) − B

(v

− v

);

j

j

j

2 x

j +1

j −1

j

2 x

j +1

j

−1

n +1

n +1

ρ j

kT

n +1

pi

=

.

m

i

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

Можно пойти на дальнейшие упрощения, приняв, что плазма движется как единый плазменный поршень, внутри которого параметры плазмы однородны по длине поршня. Подводимая в канал

плотность мощности pε = ∂∂εt = −divS , где S = 4cπ [E ×H ] − век-

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

∂

H 2

∂

H 2

pε =

(

y

) +v

(

y

) − E j =

∂t

8π

∂x

8π

∂

H y2

v c

∂H y

H

j2

∂

H y2

v

jH

j2

=

(

) +

y

+

=

(

) +

y

+

.

∂t

8π

c 4π

∂x

σ p

∂t

8π

c

σ p

Интегрируя по объему, получим:

∂

H 2

∂

H

2

P =UI =

y

dV

+

y

vdV + Rp I 2

,

∂t ∫ 8π

∂x ∫

8π

где I − ток через плазменный поршень Rp − сопротивление плаз-

менного поршня. Первое слагаемое соответствует мощности, расходуемой на заполнение объема магнитным полем, второе слагае-

65

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

c U 2

c U 2

L I 2

L xI 2

mv2

t

W0

=

0 0

=

0

+

0

+

x

+

+ ∫RI 2dt = const ,

2

2

2

2

2

0

где

c0

− емкость батареи питания,

L0 − индуктивность внешней

электрической цепи, Lx − погонная индуктивность канала (рис.22), R = R0 + Rp − сумма сопротивления внешней цепи и сопротивления плазмы, U0 ,U − начальное напряжение на батареи питания и

напряжение в момент времени t . Мощность, подводимая к плазменному поршню в терминах электротехнических макропарамет-

ров P =UI =	1	L v I 2	+ L x I I + E		+ E		. Первое слагаемое
	2	x	x	k		R

соответствует мощности по перемещению плазменного поршня.

Учитывая, что мощ-

L0

R0

ность есть сила, умно-

женная

на

скорость,

I

ускоряющую

силу

y

можно

выразить

через

C0

x

a

макропараметры

элек-

трической

цепи

F = 1

L I 2 . Тогда сис-

b

2

x

lp

тему

уравнений,

опи-

L0

R0

I

Lx

сывающая данную мо-

дель,

можно

предста-

v

H y

вить в виде:

C0

Rp

Рис. 22. Схема рельсотрона

66

x = v;

1

v =

L I 2

;

2mp

x

1

U = −

I;

c0

= (L0 + Lx x)I +

U

1 Lx

I

2

p

=

;

2

a b

a

R

=

;

σ p

b lp

p

p

n

=

;

kT

(1+γ

)

p

p

R I 2 =σ T 4 2(a

p

St

p

I (R + Lx v);

b +b lp + a lp )

где a,b,lp − высота, длина, ширина поршня (см. рис. 22), σp −

проводимость плазмы, определяемая, например, по формуле Спитцера, γ p − степень ионизации плазмы, определяемая уравнением

Саха. Последнее уравнение системы задает тепловой баланс в предположении, что все омически выделяемое в плазме тепло идет на излучение плазмы как излучение черного тела. Система замкнута. Точность расчета проверяется сохранением энергии:

		c U 2		c U 2		L I 2		L xI 2		mv2	t
	W0 =	0 0	=	0	+	0	+	x	+		+ ∫RI 2dt = const .
		2		2		2		2		2
											0

67

6.4. Применение МГД приближения для расчета пристеночного падения потенциала

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

поля, вполне можно свестиϕ к электростатической, то есть описы-

0

ваемой уравнением Пуассона: ϕ = −4πe(ni −ne ) . Будем считать,

	плазма	je	стенка
		ji
			ϕ
		l
	x		0

ϕW

Рис. 23. Пристеночное падение потенциала

что протяженность пристеночного слоя равна l, и на его границе плазма квазинейтральна, то есть

ni0 = ne0 = n0 , а потенциал на границе ϕ(l) =ϕ0 (рис. 23). Задачу

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

ji = const = ni (x)vi (x) = ni0 (x)vi0 (x) ,

где vi0 (x) – скорость ионов на гра-

нице пристеночного слоя. Из этого соотношения с учетом

68

m v 2	= e		ϕ		получим n (x) = n		ϕ	0
i i									. Распределение электронов
2					i	0	ϕ(x)

в присутствие потенциального барьера электрического поля опи-

сывается распределением Больцмана ne (x) = n0 exp(e(ϕ(x) −ϕ0 )) , kTe

где Te – температура электронов плазмы. Тогда уравнение Пуассона перепишем в виде:

d 2ϕ	= −4πen (	ϕ0	−exp(	e(ϕ(x) −ϕ0 )	)) .
dx2		ϕ(x)
	0			kT
				e

Помножив это уравнение на 2			dϕ	dx и проинтегрировав, получим
			dx
( dϕ)2	= −4πen (4ϕ1/ 2ϕ1/ 2		(x) −2		kTe	exp(	e(ϕ(x) −ϕ0 )	)) +C .
dx	0	0			e		kTe

Константу интегрирования определяем из граничного условия для потенциала ϕ(l) =ϕ0 , и с учетом того, что электрическое поле в

плазме

E = − dϕ

= 0 , то есть

0 = −4πen (4ϕ

0

−2

kTe

) +C , сле-

dx

l

0

e

kTe

довательно, C = 8πen

(2ϕ

0

−

) . Тогда, получим уравнение для

0

e

расчета потенциала:

( dϕ)2

= −8πen (2ϕ

(

ϕ(x) −1) −

kTe

(exp(

e(ϕ(x) −ϕ0 )

) −1)) .

dx

0

0

ϕ

0

e

kT

e

Это задача Коши неявного вида (правая часть зависит от ϕ(x) , численно ее можно решить, например, разобранным ранее методом итераций, необходимо только знать значение ϕ0 и задать гранич-

ное значение потенциала на стенке ϕ(0) =ϕw . Граничное условие со стороны плазмы задается критерием Бома: ϕ0 ≥ kT2ee , который физически означает необходимость «доускорения» ионов в пред-

69

пристеночном слое до скорости звука в плазме v	= c	S	=	kTe	(для
0i				mi
Ti = 0 ). Математически это условие вытекает			из	условия

( dϕ)2 ≥ 0 . Действительно, если разложить по малому параметру

dx

ϕ =ϕ −ϕ0 , сохранив второй порядок в разложении, то

dϕ

)2

≈ −8πen

(2ϕ

ϕ

( ϕ)2

kT

e ϕ

(

dx

(1+

−

8ϕ2

−1) −

e

(1+

kT

+

2ϕ

0

0

0

e

0

e

+ 1 (e ϕ)2 −1)) = 2πen ( ϕ)2

kTe + 2eϕ0

.

2

kTe

0

kTeϕ0

kTe

Так как ϕ

0

< 0 , то kT + 2eϕ

0

≥ 0 ,

то есть eϕ

0

≤ −

. Граничное

e

2

условие на стенке определяется из равенства потоков ионов и электронов, при условии, что стенка «не заземлена»:

1 n v exp(

e(ϕW −ϕ0 )

)

= n v

,

4

0 e

kT

0 i0

e

где v =

8kTe

, v

=

k(Te +Ti )

= c

S

(по критерию Бома). С уче-

e

πme

0i

mi

том

ϕ

0

= −

kTe

ϕ

=

kTe

(ln(

2πme (Te

+Ti )

) −1) .

Например, при

2e

W

2e

miTe

Ti =Te

kTe

пристеночное

падение

потенциала

ϕ

≈ −3.8

= −3.8T

[эВ] .

W

e

e

Учет электронной эмиссии с поверхности первой стенки

Поток электронов со стенки может существенно изменять распределение потенциала в пристеночной области. Это можно учесть в МГД приближении, если ввести коэффициент вторичной

электронной эмиссии γe , связывающий плотность потока вторич-

70