КузминАМ Основы теории критичности 2008

(соответствующей	значению	i =1 )	и отражателя (соответствует
i = 2 ). В таком	реакторе	R1 , R2 –	радиусы активной зоны и

экстраполированной границы соответственно, а параметр ω принимает значения:

					(1)		1
	ω = β2	=		К∞		−1		– в активной зоне при r R	,
				К			L2
	1				эф			1
							1
	ω = −β 2		= −		1	– в отражателе лри r R ,			(2.28)
					L2
		2						2
					2
				R1 = R1 ,				R2 = R2 − R1 .

Общие решения уравнений (2.21) для сферической геометрии

можно записать в виде:
ψ	0,1	(r) = A sin(β1r)		+ B cos(β1r)			,	r R ,
		1	β r	1	β	r		1
			1		1
ψ		(r) = A	sh(β2r)	+ B ch(β2r)			,	r R ,
	0,2	2	β2r	2	β2r			2
где для определения коэффициентов					A1 , B1 ,			A2 , B2 имеем два

граничных условия (2.22) и два «условия сшивки» (2.23). Из равенств (2.22) получим:

B1 = 0, B2 = −A2 th(β2 R2 ) ,

а из условий (2.23) – однородную и линейную относительно неизвестных A1 , A2 систему уравнений. Раскрывая определитель

(2.26) этой системы и приравнивая его нулю, придём к условию критичности

1	− β	R ctg(β		R ) = D2		[1	+ β R cth(β		2	R )]	(2.29)
	1	1	1	1	D1		2	1		2

и следующему выражению для собственной функции ψ0 (r) (полагая A1 =1 ):

	ψ0 (r) =ψ0,1 (r) =	sin(β1r)	,	r R1 ,
		β1 r

31

ψ0 (r) =ψ0,2	(r) =	sin(β1R1 )		sh(β2 (R2	−r))	, r	R2 .
		sh(β2	R2 )	β1 r

Остаётся с помощью равенства (2.29) определить критические параметры реактора.

Рис.2.2. Определение корней уравнений (2.30) и (2.33)

Пусть задача о критичности реактора решается в постановке 2 (когда нужно найти Кэф при заданных размерах и свойствах зон). В этом случае равенство (2.29), рассматриваемое как уравнение относительно неизвестного Кэф , определяет собственные числа

задачи k0 , k1,…. Чтобы пояснить их получение,				сделаем замену
~	= β1R1	и запишем условие критичности в виде
α
		~	~	(2.30)
		1−α ctg(α) = h0 ,

где h0 известно и совпадает с правой частью равенства (2.29).

	Рис.2.2 иллюстрирует решение уравнения (2.30) и размещение
	~	~
его корней αn (n = 0, ±1, ± 2,…) на плоскости. Видно, что все		αn
действительны, располагаются симметрично относительно		оси
~	~
α	= 0 , а при значениях α ≥ 0 :

32

~	π, n = 0,1,2,…
nπ <αn < (n +1)

			~
Поскольку значения корней αn и собственных чисел kn связны
(как это следует из (2.28)) равенствами
(1)		R2
		1
kn = К∞	2	~2	2	, n = 0,1,…,	а	эффективный
	R1	+αn	L1

коэффициент размножения нейтронов Кэф совпадает с

максимальным собственным числом k0 , то получим

Кэф = К∞(1)

R2

= К∞(1)

1

1

,

(2.31)

2

~2

2

2

2

~2

R1

+α0

L1

1+α0

L1

~

где α

2

=

α0

наименьший по модулю корень уравнения

0

R2 , а

α0 –

1

~

(2.30),

удовлетворяющий неравенству

< π . Из равенства

0 <α0

(2.31) следует, что в критическом реакторе (когда Кэф =1 ) радиус активной зоны должен совпадать с величиной:

	~	L1
R1 =	α0
			.	(2.32)
	К∞(1) −1

Очевидно, к такому же значению радиуса придём, решая задачу о критичности реактора в постановке 1 и считая, что неизвестным

параметром является R1 . В этом случае вместо (2.30) придётся анализировать корни уравнения:

						~	~	~	(2.33)
						1−α ctg(α) = h(α) ,
~	D2			β2		~
где h(α) =		1	+			α cth(β2	R2 )	является линейной функцией
	D			β
	1				1

параметра α~ . При этом необходимо ограничиться лишь значениями α~ > 0 , чтобы не получить отрицательных значений для R1 , и выбирать, как и раньше, наименьшее число α~0 .

Сравнивая значения Кэф , R1 , рассчитанные по формулам (2.31), (2.32) для реактора с отражателем, со значениями Кэф(0) , R1(0) , полученными по формулам (2.3), (2.7) для реактора без отражателя,

33

полезно отметить следующее. В условно-критических реакторах с одинаковыми составами и радиусами R1 = R1(0) зон Кэф > Кэф(0) , а в

критических реакторах с одними и теми же составами активных зон R1 < R1(0) . Это связано с тем, что в реакторе с отражателем

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

Рис.2.3. Особенности функции f (x) = x cth x

Изменение свойств отражателя (прежде всего, толщины R2 и сечения поглощения Σ(a2) ) влияют на критический размер активной зоны. С ростом R2 (при сохранении остальных параметров)

уменьшается значение

производной

dh

~ cth(β

R )

и,

как

~

2

2

~

и R1

dα

R2

≥ 2L2

. Однако при

это

видно из рис.2.2, снижаются α0

снижение практически прекращается.

Для выяснения влияния

Σ(2)

на

значение R

правую

часть

a

1

равенства (2.33) преобразуем к виду

~

D2

~

R2

α x cth(x)

h(α)

=

1+

β

R

, x =

L

.

D

1

1

2

2

34

Учитывая зависимости, представленные на рис.2.2 и 2.3, нетрудно установить, что:

- при

Σ(2)

→ ∞ значения

L → 0, x → ∞,

x cth(x) → ∞ ,

a

2

~

(0)

π

а поэтому

α0 → π ,

R1

→ R1

=

,

β1

- при стремлении Σ(2) к нулю:

L

→ ∞,

x → 0,

x cth(x) →1, а

a

~

2

D2

~

α

+ β

R

функция h(α) → D

1

.

1

1

2

~

Если

в

последнем

случае

принять,

что

α

<<1 (а

это

β

R

1

2

выполняется,

когда

R >> R(0) )

и

D = D ,

то

~

~

~

2

1

2

1

корень

этого уравнения

1−α0 ctg(α0 ) = h(α0 ) ≈1. Наименьший

~

π

значит R1

≈

π

=

R1(0)

α0 ≈

,

а

.

Такое

сокращение радиуса

2

2β1

2

активной

зоны

R1 приведёт

почти

к

8-кратному

уменьшению

критической массы реактора. Для этого необходимо, чтобы отражатель не поглощал нейтроны и имел достаточно большую толщину.

В заключение отметим, что влияние отражателя на критический размер и распределение нейтронов в активной зоне можно оценить, пользуясь понятиями эффективной добавки δ и альбедо η

отражателя. В качестве эффективной добавки принимают разность характерных размеров (радиуса или полутолщины) двух критических реакторов: реактора без отражателя и активной зоны

реактора с отражателем (при одинаковых значениях β12 активных зон). Например, в случае сферического реактора:

δ = R(0)	− R .	(2.34)
1	1	[4] следующее
Используя определение (2.33),	можно получить

выражение для оценки δ активной зоны большого реактора (когда

δ << R1 и tg(β1δ) β1δ ):

35

1

(0)

D2

(0)

D2

2

(0)

δ =

+

δ0

−

+

δ0

,

2

R1

D1

R1

D1

−4R1 δ0

	D				R
где δ0 =		1	L2			2
	D				L
		2		th			.
					2

Отметим, что эффективная добавка δ в основном определяется свойствами отражателя и слабо зависит от свойств активной зоны. Если величина δ известна, то распределение нейтронов в активной зоне реактора с отражателем можно получить, рассматривая

реактор без отражателя с прежним значением параметра β12 , но увеличенным на δ размером Rэ = R1 +δ . Этим обстоятельством

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

Под альбедо

отражателя

η

подразумевают отношение

односторонних токов нейтронов J −

(из отражателя в активную

зону) и J + (из активной зоны в отражатель):

1 dφ

1

+ 2D

φ dr

J −

r=R

.

(2.35)

η = J + =

1 dφ

1

1−2D

φ dr

r=R

1

В частности, для сферического реактора с отражателем имеем

1−ζ

2D

R

R

η = 1+ζ ,

ζ =

2

1

2

1+

L

R

cth

L

.

1

2

2

Из определения (2.35) можно получить:
	1 dφ		1−η
D		= −		.	(2.36)
	φ dr r=R		2 (1+η)
		1

Выражениями (2.35) и (2.36) удобно воспользоваться при расчётах реакторов с несколькими отражателями и при выравнивании поля энерговыделения в активной зоне. В этих случаях достаточно, определив вначале η по формуле (2.35), использовать это значение

36

в последующих расчётах реакторов с разными свойствами активной зоны, формулируя условие на внешней границе активной зоны в виде (2.36).

2.4. Численный метод решения условно-критической задачи

Аналитический способ нахождения Кэф и асимптотического

потока нейтронов, рассмотренный в разделе 2.2, имеет ограниченную область применения. Его не удаётся распространить на случай многозонных реакторов в 2- и 3-мерных геометриях. Даже в одномерных геометриях, но с большим числом зон, получаем весьма сложного вида нелинейное уравнение (2.26), которое трудно решить аналитически. Приходится обращаться к численным итерационным методам нахождения его корней (например, методу Ньютона). Но в этом случае в реакторах с большим размером активной зоны среди всех собственных чисел

трудно выделить наибольшее k0 = Кэф .

Поясним это на примере сферического реактора без отражателя, собственные числа которого km рассчитываются по формулам (1.9)

с использованием	αm =	π (m +1) .	Нетрудно		видеть,	что
km →1 (m =1,2,…)		Rэ
	при	R → ∞ , т.е. происходит сближение
k0
собственных чисел с ростом радиуса R активной зоны. Поэтому в
процессе итерационного счёта вместо числа				k0 можно получить
какое-либо другое	число	km ≠ k0 .	Так	как	такому	числу

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

37

Напомним, что в методе итераций источников на каждой n-й

итерации решаются уравнения:
div (J (n) (r ))+Σa φ(n) (r ) = q(n−1) (r ) ,		φ(n) (r ) Ω ,	(2.37)
где ток	нейтронов J (n) (r ) связан с потоком		φ(n) (r )
соотношением:
	grad φ(n) (r ) +3Σtr J (n) (r ) = 0 .		(2.38)
Источник	q(0) (r ) на 1-й итерации	задаётся в виде	любой

неотрицательной функции, а на следующих итерациях рассчитывается по формуле: q(n−1) (r ) =ν f Σf φ(n−1) (r ) . В главе 1

показано, что такие итерации сходятся и в пределе дают Кэф и

асимптотический поток нейтронов φ(r ) .

На практике итерации прекращают, когда впервые с заданной погрешностью ε > 0 выполнится неравенство:

∫q(n) (r ) dV

К(n) − К(n−1)	≤ ε , где К(n) =	V		.	(2.39)
		∫q(n−1)
			(r ) dV

V

Если это произойдёт на итерации с номером n = p , то принимают:

Кэф = К( p) , φ(r ) = Cφ( p) (r ) .

(2.40)

При этом в задачах со значением Кэф , сильно отличающимся от

единицы, рекомендуется источник q(n) (r ) рассчитывать по формуле:

q(n) (r ) = К1(n) ν f Σf φ(n) (r ) .

Перейдём теперь к рассмотрению численного метода нахождения потоков нейтронов φ(n) (r ) . В основе метода лежит

преобразование дифференциальных уравнений (2.37) в систему линейных алгебраических уравнений относительно значений

φ(n) (rk ) в отдельных узлах rk заранее выбранной сетки разбиения.

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

38

пока случаем многозонного реактора в одномерной геометрии, в которой поток нейтронов зависит от одной переменной r и

выполняется условие симметрии		ν	dφ(n)		= 0 .
	r		dr
				r=0

Опуская в дальнейшем номер итерации n и учитывая, что в

рассматриваемой геометрии

divJ (r ) =

1

d

(rν J (r)),

gradφ

(r ) =

dφ(r)

r

,

rν

dr

dr

r

запишем уравнения

(2.37),

(2.38)

в виде системы

дифференциальных уравнений первого порядка:

−

d

(rν

J (r))−rν Σ(r)φ(r) + rν q(r) = 0 ,

(2.41)

dr

dφ(r)

~

dr

+3Σ(r) J (r) = 0 ,

(2.42)

где источник q(r) - кусочно-непрерывная функция, а транспортное

~		Σ(r) - кусочно-постоянные
сечение Σ(r) и сечение поглощения
функции, терпящие	разрывы на	границах раздела зон
r = Ri (i =1,2,…, n) .	Решения J (r) ,	φ(r) ищутся на множестве

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

J (0) = 0 , φ(Rn ) = 0 .

(2.43)

Здесь (как и раньше) принимается, что r = Rn находится на

экстраполированной границе.
Введём в	рассмотрение	две	сетки узловых			точек:
{rk , k =1,2,…, N} – основная		сетка	и	{rk−(1 / 2) ,	k =1,2,…, N} –
вспомогательная	сетка. Координаты			точек	этих	сеток
rk , rk−(1 / 2) чередуются и располагаются так, что:
	rk−(1 / 2) < rk < rk+(1 / 2) < rk+1 .

Проинтегрируем уравнение (2.41) по переменной r в пределах интервала rk−(1 / 2) ≤ r ≤ rk+(1 / 2) , а уравнение (2.42) – в пределах интервала rk ≤ r ≤ rk+1 . В результате получим равенства

39

	rk +(1 / 2 )
rkν+(1/ 2) J (rk+(1/ 2) ) − rkν−(1/ 2) J (rk−(1/ 2) ) +		∫(Σφ − q)rν dr = 0 ,	(2.44)
r k +1~	rk −(1/ 2)
φ(rk+1 ) −φ(rk ) +3 ∫Σ(r) J (r) dr	= 0 ,	k =1,2,…,	(2.45)

r k

которые являются точными аналогами уравнений (2.41) и (2.42). Причём, первое из них выражает условие баланса нейтронов в

rk +(1 / 2 )

объёме Vk = ∫rν dr .

rk −(1 / 2 )

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

располагаются равномерно с шагом r(i) . При этом возможны две

~

схемы размещения. По первой схеме на границы зон (где Σ, Σ, q

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

r(i) → 0 ) они эквивалентны, но в случае конечного шага

r(i) разностные уравнения, соответствующие первой схеме размещения точек, дают более высокий порядок локальной

аппроксимации [5]. Это означает, что при конечном шаге сетки первая схема позволяет получить значения потоков φ(rk ) с

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

Пусть Ni – число точек вспомогательной сетки в i -й зоне.

n

Тогда общее число узлов каждой сетки N = ∑Ni , а шаги

							i=1
r(i) будут равны:
r(1) =			R1	,	r(i) =	Ri , i = 2,3,…, n . (2.46)
	N		−0,5
		1				N	i

40