Апсе Исползование программы ТИМЕ26 2008

−

1 ∂

r D (r, z)

∂

φ

(r, z) −

∂

D (r, z)

∂

φ

(r, z) +

r ∂r

K

∂r

K

1

∂z

K

∂z

K

+Σcfd , K (r, z)φK (r, z)=

χK ∑νΣm (r, z)φm (r, z)+

(5)

Kэф

m

+K∑−1Σm→r (r, z)φm (r, z).

m=1

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

	φK (r, z) =ϕK (r) cos	πz	,
		HАЗ +2δэф, Z

где HАЗ – высота активной зоны; δэф,Z – эффективная добавка в аксиальном направлении.

Подставляя φK (r, z) в уравнение (5), получим:

−

1 d

[rDK (r)

dϕK (r)

2

DK (r)]ϕK (r) =

]+[Σcfd , K + ωZ

r dr

dr

(6)

1

=

χK ∑νΣf ,m (r)ϕm (r)+K∑−1Σm→K (r)ϕm (r),

Kэф

m=1

где ωz

2 =(

π

)2 −

баклинг (или лапласиан), описываю-

H АЗ +2δэф, Z

щий аксиальную утечку нейтронов из активной зоны.

Можно показать, что уравнение (6) может быть обобщено на случай произвольной одномерной геометрии (плоскость, цилиндр, сфера) и записано следующим образом:

11

−

1

d

[r

α

DK (r)

dϕK (r)

]+[Σcfd , K (r)+ωZ

2

DK (r)]ϕK (r)=

rα

dr

dr

(7)

=

1

χK ∑νΣf ,m (r)ϕm

(r)+K∑−1Σm→K (r)ϕm (r),

Kэф

m=1

где α – показатель геометрии, равный 0 для плоскости, 1 для цилиндра и 2 для сферы; ωz2 – баклинг, описывающий утечку нейтро-

нов. Для плоскости и цилиндра ωz 2 = ( πHАЗ + 2δэф, Z )2 , а для сферы ωz2 = 0. В более общем случае значения аксиального баклинга ωz2 могут быть взяты из других, например, двумерных расчетов и введены в программу как исходные данные.

Физический смысл членов, входящих в уравнение (7), вполне объясним. Первый член левой части описывает перенос нейтронов в радиальном направлении, второй член – поглощение, замедление и утечку нейтронов в аксиальном направлении. Два члена правой части уравнения отвечают за генерацию нейтронов деления и «приход» нейтронов в рассматриваемую энергетическую группу из всех предшествующих за счет замедления.

Решение уравнения (7) осуществляется методом итераций источников, в соответствии с которым член, описывающий генерацию нейтронов деления, рассчитывается на нейтронных потоках предыдущей итерации:

−	1		[	d	rα D (r)	dϕK (n)	] +
	rα			dr	K	dr
+[Σcfd , K +ωZ		2 DK (r)]ϕK (n) =qK (n−1) +K∑−1						Σm→K ϕm(n) ,
							m=1

где n – номер текущей итерации источников;

qK	(n−1) (r) =	1	χK ∑νΣ f ,m ϕm(n−1) (r);
		(n−1)
		Kэф

12

Kэф(n−1) =∑ ∫ νΣ f ,m (r) ϕm(n−1) (r) rα dr. m

Процесс итераций продолжается до тех пор, пока значение Kэф не стабилизируется в рамках заданной точности ε, т.е.

Kэф(n) −Kэф(n−1) ≤ε.

Kэф(n)

Таким образом, в рамках каждой итерации источников необходимо решить следующее уравнение:

−	1	[	d	rα DK (r)	dϕK	] +ΣK (r)ϕK =qK +K∑−1Σm→K ϕm =qK ', (8)
	α		dr
	r				dr		m=1
где ΣK (r) =Σcfd, K (r) +ωZ							2 DK (r) , т.е. найти ϕK (r) при известной

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

Первым шагом к решению уравнения (8) является получение его конечно-разностного аналога, т.е. замена непрерывной пространственной зависимости ϕK (r) значениями этой функции в дискрет-

ных пространственных точках. Для этого в каждой зоне системы вводится равномерная (т.е. с постоянным шагом) сетка пространственных точек. На границах раздела зон размещаются две точки с одинаковой координатой: одна слева от границы, другая – справа от нее, т.е. эти точки принадлежат разным зонам системы. Затем уравнение (8) интегрируется в окрестности произвольной точки l, т.е. в пределах интервала [rl −0,5∆ri , rl +0,5∆ri ], где rl – координа-

та точки, а ∆ri −шаг между точками в зоне i. Такое интегрирование приводит к следующим результатам:

13

∫

1 d

α

dϕK

α

α

dϕK

α

dϕK

r

DK

r

dr =(r

DK

)l+0,5 −(r

DK

)l−0,5

=

rα

dr

dr

dr

dr

=(r +0,5∆r )α D

ϕK ,l+1 −ϕK ,l

−(r −0,5∆r )α D

ϕK ,l −ϕK ,l−1

;

K ,i

K ,i

l

i

∆ri

l

i

∆ri

∫ΣK (r)ϕK (r) rαdr =ΣK ,i ϕK ,l rl α ∆ri ;

∫qK (r) rα dr =qK ,l rl α ∆ri .

Окончательно конечно-разностный аналог уравнения (8) может быть записан в следующем 3-точечном виде (опуская индекс энергетической группы K):

Al ϕl−1 +Bl ϕl

+Cl ϕl+1 =Ql ;

(9)

где

A =−

Di

(r −0,5∆r )α;

C =−

Di

(r +0,5∆r )α;

∆r

∆r

l

l

i

l

l

i

i

i

B =Σ

r α

∆r +

Di

[(r +0,5∆r )α +(r −0,5∆r )α];

∆r

l

i l

i

l

i

l

i

i

Ql =ql rl α ∆ri .

Итак, дифференциальное уравнение (8) преобразовалось в систему алгебраических уравнений (9), число которых равно количеству точек конечно-разностной сетки. Решение этой системы находится так называемым методом «двойной прогонки», суть которого заключается в следующем. Предположим, что значение плотности потока нейтронов в точке l может быть выражено через ее значение в следующей точке l+1 следующим образом:

ϕl =αl ϕl+1 +βl .

(10)

Подставив это выражение в уравнение (9), получим систему рекуррентных уравнений для нахождения коэффициентов αl и βl

через их значения в предыдущей точке αl −1 и βl−1 :

14

Al (αl−1 ϕl +βl−1) +Bl ϕl +Cl ϕl+1 =Ql ;

(Al αl−1 + Bl )ϕl = −Cl ϕl+1 +(Ql − Al βl−1);

	αl =−	Cl		;
		Al αl−1	+Bl
				(11)

βl =Ql − Al βl−1 .

Al αl−1 +Bl

Для того чтобы начать вычисления коэффициентов αl ,βl во всех пространственных точках системы, необходимо определить их значения в первой точке (значения α1,β1 ), а затем использовать рекуррентные соотношения (11).

Для вычисления значений α1 и β1 применяется краевое условие

уравнения (1) в центре системы. В большинстве случаев расчета БР и ЭЛЯУ – это условие симметрии нейтронного потока, т.е.

rϕK (r = 0) = 0.

Интегрируя уравнение (8) в окрестности первой точки (она размещена на расстоянии полшага от центра системы), можно получить соотношение:

[Σ (0,5)	α ∆r 2	− D ]ϕ =q (0,5)α ∆r 2											− D ϕ	2	;
1	1				1			1	1	1			1
то есть									D1
	α1 =											;
				D			−Σ		(0,5)α ∆r 2
		1						1		1
						q		(0,5)α ∆r 2
	β =						1			1	.
							(0,5)α ∆r 2
	1			Σ						− D
		1							1	1

Теперь, используя рекуррентные соотношения (11), можно определить значения коэффициентов αl ,βl во всех точках первой зо-

ны, кроме последней точки. Последняя точка первой зоны n и первая точка второй зоны n+1 имеют одну и ту же координату

15

rn = rn+1 , но расположены по разные стороны границы раздела зон. Для определения значений αl ,βl в последней точке первой зоны и

в первой точке второй зоны используются условия непрерывности плотности потока и тока нейтронов на границе раздела зон:

ϕ(rn ) = ϕ(rn+1);

D1 r ϕ(rn ) = D2 r ϕ(rn+1).

Из условия непрерывности плотности потока нейтронов следует, что αn =1;βn = 0.

Условие непрерывности тока нейтронов позволяет определить значения αn+1 и βn+1 следующим образом. Уравнение (8) необходимо проинтегрировать сначала слева от границы зоны, т.е. в интервале [rn −0,5∆r1 , rn ] , а потом справа от границы зоны, т.е. в интервале [rn+1 ,rn+1 +0,5∆r2 ] , и результаты сложить.

Первое интегрирование дает:

−r

α D (

dϕ

)

+ r

α D ϕn −ϕn−1

+ 0,5Σ ϕ

r α ∆r = 0,5q

r

α ∆r .

n

1 dr

n

n−0,5

1

∆r

1

n

n

1

n

n

1

1

Второе интегрирование дает:

−r

α D

ϕn+2 −ϕn+1

+r

α D (

dϕ

)

+

n+1.5

2

∆r

n+1

2

dr

n+1

2

+0,5Σ

2

ϕ

r

α

∆r

= 0,5 q

r

α ∆r .

n+1 n+1

2

n+1 n+1

2

После сложения результатов интегрирования получим:

2(

rn−0,5

)α D ϕn −ϕn−1

+2(

rn+1,5

)α

D

ϕn+1 −ϕn+2

+

r

r

1

∆r

2

∆r

n

1

n

2

+(Σ1 ∆r1 +Σ2 ∆r2 )ϕn =qn ∆r1 + qn+1 ∆r2.

16

Учитывая, что

ϕn+1 = ϕN , а ϕт−1 = αт−1 ϕт +βт−1 = αт−1 ϕт+1 +βт−1,

можно связать значения ϕn+1 и ϕn+2 следующим уравнением:

rn−0,5

α

D

rn+1,5

α D

ϕ

[2(

)

1

(1−α

n−1

)

+ 2(

)

2

+Σ ∆r

+Σ

2

∆r ] =

n+1

rn

∆r1

rn

∆r2

1 1

2

rn+1,5

α

D

rn−0,5

α

D

=ϕ

[2(

)

2

]+(q

∆r

+ q

∆r

) + 2(

)

2

β

.

r

∆r

r

∆r

n+2

n

1

n+1

2

n−1

n

2

n

2

Из последнего уравнения легко найти значения коэффициентов αl ,βl в первой точке второй зоны:

rn+1,5 α

D

2(

)

2

r

∆r

αn+1 =

n

2

;

(12)

rn−0,5

α

D

rn+1,5

α

D

2(

)

1

(1−α

n−1

) +

2(

)

2

+ Σ ∆r +Σ

2

∆r

rn

∆r1

rn

∆r2

1 1

2

rn−0,5

α D

qn∆r1 + qn+1

∆r2 + 2(

)

1

βn+1 =

rn

∆r1

.

rn−0,5

α

D

rn+1,5

α D

2(

)

1

(1−α

n−1

)

+ 2(

)

2

+ Σ ∆r

+Σ

2

∆r

rn

∆r1

rn

∆r2

Затем, можно снова воспользоваться рекуррентными формулами (11) для нахождения значений коэффициентов αl ,βl во всех

точках второй зоны, кроме последней ее точки. В последней точке этой зоны коэффициентов α =1,β = 0 , а в первой точке следующей

зоны надо воспользоваться соотношениями (12). В результате использования такого алгоритма в каждой зоне можно будет найти значения коэффициентов αl ,βl для всех пространственных точек

системы.

17

Таким образом, так называемой «прямой» прогонкой, т.е. в направлении от центра к периферии, вычисляются коэффициенты αl ,βl для уравнения (10). Но для того чтобы воспользоваться этим

уравнением, необходимо «зацепиться» за значение плотности потока нейтронов в последней точке системы и рассчитать значения ϕl во всех предыдущих точках уже в направлении от периферии к

центру, т.е. осуществить «обратную» прогонку.

Значение плотности потока нейтронов в последней точке может быть найдено из краевого условия задачи. Например, в большинстве расчетов используется нулевое граничное условие, т.е. ϕN = 0 , где

N – номер последней пространственной точки, находящейся на экстраполированной границе системы. Тогда, используя уравнение (10), можно найти значения ϕl во всех предшествующих точках. Однако

программа TIME26 делает возможным использование и другого краевого условия, а именно условия симметрии на границе реактора, т.е. r ϕ(rN ) = 0. Для нахождения значения ϕN необходимо проин-

тегрировать уравнение (8) в окрестности последней точки. В этом случае она отстоит на полшага от границы системы, то есть интегрирование производится в интервале [rN −0,5∆ri , rN +0,5∆ri ]. Интег-

рирование дает следующий результат:

(rαDi ddrϕ)N −0,5 +Σi ϕN rN α ∆ri =qN rN α ∆ri .

Поскольку

dϕ

(r

) = ϕN −ϕN −1 = ϕN (1−αN −1) −βN −1 , зна-

dr

N −0,5

∆ri

∆ri

чение ϕN может быть вычислено по следующей формуле:

rN −0,5

α

D

q

N

+

(

)

i

β

N −1

rN

∆r

2

ϕN =

i

.

rN −0,5

α

D

(

)

i

(1−α

N −1

) +Σ

i

rN

∆r 2

i

18

Итак, «прямой» прогонкой, т.е. в направлении от центра к периферии, находятся значения коэффициентов αl и βl во всех про-

странственных точках системы. Затем, из краевых условий определяется значение плотности потока нейтронов в последней точке системы, и «обратной» прогонкой, т.е. в направлении от периферии к центру, вычисляются значения ϕl во всех точках. Этот алгоритм

и получил название метода «двойной» прогонки.

1.2. Решение уравнений кинетики изотопного состава

Ключевым моментом в решении уравнения (3), описывающего временное поведение концентраций изотопов ρ(t) при работе системы на постоянной мощности, является вычисление экспоненци-

ˆ ∆

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

ˆ

A также не зависят от времени.

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

ˆ

сто при работе БР или ЭЛЯУ на мощности, матрица A может быть записана в нижнетреугольном виде, то есть ее элементы ai, j = 0 приi < j. Хотя программа TIME26 способна решать урав-

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

ˆ

A .

			ˆ
Известно, что матрицу A можно представить в следующем виде:
		ˆ	ˆ	ˆ	−1	;
		A=D diag(λ) D
где	ˆ
	D – матрица, столбцами которой являются собственные векто-
	ˆ	– диагональная матрица, элементами кото-
ры матрицы A ; diag(λ)

ˆ

рой являются собственные значения матрицы A .

19

Нетрудно показать, что, если известны собственные векторы и

собственные значения λ матрицы		ˆ
		A , то экспонента от нее может
быть вычислена следующим образом:
ˆ	ˆ	ˆ	−1	(13)
exp(A) = D diag[exp(λ)] D			.

Действительно:

ˆ = ˆ + ˆ + exp(A) E A

ˆ 2		ˆ n	ˆ ˆ ˆ	−1
A		A
	+... +		= DED		+
2!		n!

ˆ	ˆ	−1	+
Ddiag(λ) D

ˆ	λ2 ˆ −1	ˆ	λn ˆ	−1
+ Ddiag(	2! ) D	+... + Ddiag( n! ) D			=
	ˆ	ˆ	−1
	=Ddiag[exp(λ)]D .

В случае нижнетреугольной матрицы определение ее собственных значений и собственных векторов существенно упрощается. Легко показать, что собственными значениями нижнетреугольной матрицы являются ее диагональные элементы ai,i , а элементы мат-

ˆ	ˆ −1	можно вычислить по следующим формулам:
риц D и	D
				i−1
		ˆ		∑ai,l dl, j
				l= j
		(D)i, j = di, j =					; di,i =1;
					a j, j −ai,i
					i
		ˆ −1			∑al, j di,l '
					l= j+1
		(D	)i, j = di, j ' =					; di,i ' =1.
					ai,i −a j, j

ˆ ∆

После определения экспоненты от матрицы A t новый изотопный состав может быть определен по формуле:

r	ˆ	r
ρ(t0	+∆t) = exp(A∆t) ρ(t0 ).

20