авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

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

Pages:     | 1 || 3 | 4 |

«НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ ЦЕНТР КУРЧАТОВСКИЙ ИНСТИТУТ На правах рукописи ...»

-- [ Страница 2 ] --

( (2.54) Умножая выражение (2.51) на (k1), вычитая его из (2.52) и применяя обозначения (2.50) получим f k, j = (k1) jk, j + Ф k X (k3) cos(2 j ). (2.55) Умножая выражение (2.53) на (k), вычитая его из (2.54) и применяя обозначения (2.50) получим f k, j = (k1) jk, j +Ф k X (k3) cos(2 j ). (2.56) Подставим выражение (2.56) в равенство (2.42), получим ~ (0) (( ) ) S (k1) + (k1) j k, j (t ) + Ф k (t ) X (k3 ) (t ) cos(2 j ) = 4SФ k (t ) 2 (k1) C jk (t ) Fv(,0k) Ф k (t ).

t t j =1 j (2.57) Подставим выражение (2.56) в равенство (2.44), получим ~ (1) (( ) ) S (k1) + (k1) j k, j (t ) + Ф k (t ) X (k3 ) (t ) cos(2 j ) cos( j ) = Fv(,1k) X k (t ) 2 (k1) C (j1,)k (t ) t j t j = (2.58) Подставим выражение (2.56) в равенство (2.46), получим ~ ( 2) (( ) ) S (k1) + (k1) j k, j (t ) + Ф k (t ) X (k3 ) (t ) cos(2 j ) sin( j ) = Fv(,2k) X k (t ) 2 (k1) C (j2,k) (t ) t j t j = (2.59) Подставим выражение (2.56) в равенство (2.48), получим S (( (k1) + (k1) )jk, j + Ф k X (k3 ) cos(2 j ) )cos(2 j ) = j = (2.60) ~ = 4 SX ( t ) Fv(,3k) X (k3) (t ) 2 (k1) C(jk) (t ).

( 3) t j t k Предположим, что ток, по аналогии со стационарным случаем можно представить в виде ( ) ( Ф ) jk, j = (k1) + (k1) + X (k3) cos(2 j ) + Ф k X (k3) cos(2 j ) + k, j. (2.61) k Следует отметить, что в стационарном случае вектор k, j отвечает за влияние высших гармоник и в случае, когда число сторон равно числу пробных матриц, он равен нулевому вектору. Определим k, j в нестационарном случае. Для этого подставим выражение (2.61) в (2.57) – (2.60) и получим 1 ~ (0) (0) 1 (1) k,3 + k,1 + k, 4 + k, 2 = S 2 k t C jk (t ) S Fv,k t Ф k (t ), j 1 ~ (1) (1) (1) 1 (1) k,3 k,1 = S Fv,k t X k (t ) S 2 k t C jk (t ), (2.62) j = 1 F ( 2) X ( 2) (t ) 1 2 (1) C ( 2) (t ), ~ t k j t k,4 k,2 v,k k jk S S + = 1 F (3) X (3) (t ) 1 2 (1) ( 3) ~ k C jk (t ).

k,3 k, 4 t j t k,1 k,2 v,k k S S Видно, что в нестационарном случае ( k, j ) j = ( k, j )cos( j ) j = ( )sin( ) k, j j j = ( )cos(2 ) k, j j j = и в случае совпадения числа пробных матриц и сторон ячейки. Однако, можно сделать предположение о значительно меньшей степени влияния членов, стоящих в правой части равенств (2.62), по сравнению с другими членами в равенстве (2.61) и отбросить их при дальнейшем выводе уравнений. Таким образом, далее будем считать, что ( ) ( Ф ) jk, j = (k1) + (k1) + X (k3) cos(2 j ) + Ф k X (k3) cos(2 j ).





(2.63) k С учетом выражения для токов (2.63) и того, что для выделенных углов 0°, 90°, 180°, 270° (соответственно j = 1, 2, 3, 4) справедливо cos(2 j ) cos( j ) = cos( j ) (2.64) cos ( j ) = уравнения (2.31), (2.33), (2.35), (2.37) можно записать в виде ( 0 ),( 0 ) (0) (0) 4 SI k (t ) + 2Fv,k t I k (t ) + 2 t C jk (t ) (0) j ( )( ( ) ) (Ф k Ф k ) X k(3) X (k3 ) cos(2 j ) = (1) + S k + k (1) j = (1) (1),( 0 ) 4 SI k (t ) + 2Fv, k I (k1) (t ) + 2 C (j1,)k (t ) + t j t [( ] )( ( )) + S (k1) + (k1) (Ф k Ф k ) X (k3) X (k3) cos( j ) = j = 4 SI ( 2) (t ) + 2F ( 2),(0) I ( 2) (t ) + 2 C ( 2) (t ) + t j t k v,k k j,k [( ] )( ( ) ) + S (k1) + (k1) (Фk Фk ) X (k3) X (k3 ) cos(2 j ) sin( j ) = j = ( 3),( 0 ) ( 3) ( 3) 4 SI k (t ) + 2Fv,k t I k (t ) + 2 t C j,k (t ) + ( 3) j [( ))] )( ( (Ф k Ф k )cos(2 j ) X (k3) X (k3 ) = (1) + S k + k (1) (2.65) j = Представим матрицу ( (k1) + (k1) ) в виде ( ) ( ) () () D j (D j + D k ) D k = H j,k 1 2 1 1 1 + (j1) = (j1) (j1) (k1) + (j1) (k1) = (k1) (1) k h h где h (1) Dk = ( k ), (2.66) H j, k = 2D j (D j + D k ) D k.

(2.67) h – шаг решетки (для квадратной решетки совпадает с S). Тогда первое уравнение системы (2.65) имеет вид (0) I k (t ) + 2 C (jk) (t ) S H j, k (Ф k Ф k ) + S H j, k (X (k3) X (k3) )cos(2 j ) = 4 1 4SI (k0 ) (t ) + 2Fv(,0k),( 0 ) t t h j =1 h j = j Введем обозначения ( 0) 4 ( 0) (1) k = h ( k k ) (2.68) (1) = 4 ( (1) (3) ) k k k h и заметим, что величина ( ) h (0) I (k0) = (k0 ) (k1) Фk = k Фk.

Также введем обозначение ~ Fv(,lk),( 0 ) = Fv(,lk),( 0) kl),l = 0;

3 (если l=0, то l’=0, если l=3, то l’=1).

( (2.69) Посредством формул (2.68) и (2.69) преобразуем первые два члена в первом уравнении системы (2.65), получим 1 ~ ( 0 ), ( 0 ) Фk (t ) + C(j0, k (t ) 0Фk (t ) + (k0 )Фk (t ) + S(k0 ) (t ) = 0, 2 (2.70) Fv, k ) S t h j t где вектор S (k0) (t ) играет роль источника и имеет вид, S (k0 ) (t ) = 1 X k3) (t ) ( (2.71) 0 и 1 – конечно-разностные операторы, такие, что 0 A k = 2 H j,k ( A j A k ), h j = (2.72) A = H j,k (A j A k ) cos(2 j ), 1 k h 2 j = где A k – произвольный вектор.

Аналогичным образом преобразуем четвертое уравнение системы (2.65) и получим 1 ~ (3),( 0 ) (3) Fv, k X k (t ) + C (j3, k (t ) 0 X (k3) (t ) + (k1) X (k3) (t ) + S (k3) (t ) = 0, (2.73) ) t h j t S где S (k3) (t ) = 1Ф k (t ). (2.74) Уравнение для предшественников запаздывающих нейтронов.

Перейдем теперь к уравнению для предшественников запаздывающих нейтронов (1.4). Для наглядности запишем его еще раз G C j (r, t ) = j C j (r, t ) + j,g (r ) f, g (r, t ) g (r, t ), (1.4) t g = где g (r, t ) = g (r,, t )d.

Умножим уравнение (1.4) слева на и справа на ( k т (l ) ) 1 k т ( l ) (r ) + + d, j,k, проинтегрируем по объему ячейки (l ) G C j, k (t ) = j C (jl,)k (t ) + ( +,т0(l ) ) 1 +,т0(l ) (r ) j,g (r ) f, g (r, t ) g (r, t ) d, j,k dr (2.75) t k k g = Vk Учитывая формулу (1.31) L (kL ) (r,, t ) = (kl ) (r, )I (kl ) (t ) l = получим [ ] [ ] L L = (kl,)0(r )I (kl ) (t ) g.

g (r, t ) = g (r,, t )d = (l ) (r, )I (kl ) (t ) g d (2.76) k 4 l =0 l = Тогда, считая в (2.75), что gj (r) и gf (r, t ) симметричны относительно центра ячейки f (r ) k,0(r)dr = 0, l = 1, 2, т (l ) Vk где f (r ) – групповой вектор-строка, каждый элемент которого определяется как [ f (r)]g = j,g (r) f, g (r, t ).





и применяя обозначения из (2.50) для Ф k и (2.68) получим (l ) ) G h C j,k (t ) = j C (jl,)k (t ) + ( k,т0( l ) ) 1 k,т0( l ) (r ) j,g (r ) f, g (r, t ) (k0,0 (r ) (k0 ) k ( t ) d, j,k dr + + t g g = V k или в более краткой форме (l ) C j,k (t ) = j C (jl,)k (t ) + (j0k) Ф k (t ) l 0, (2.77) t, где [ ] h + т ( 0 ) 1 + (j0k) = ( k,0 ) dr k,т0( 0) (r ) d, j,k f (r ) (k0,0 (r ) (k0 ), т ) (2.78), 4 Vk Таким образом, имеем следующие системы уравнений для случая четырех пробных матриц 1 (0) Ф k (t ) + C (j0,k) (t ) 0 Ф k (t ) + (k0 ) Ф k (t ) + S k0 ) (t ) = v inv,k ( t h j t S 1 v ( 3) X ( 3) (t ) + 2 C ( 3) (t ) X ( 3) (t ) + (1) X ( 3) (t ) + S ( 3) (t ) = 0 k S inv,k t k h j t j,k k k k ( 0) (2.79) S k (t ) = 1 X k (t ) ( 3) S ( 3) (t ) = Ф (t ) 1 k k (l ) C j,k (t ) = j C j,k (t ) + j,k Ф k (t ) l 0, l = 0;

3.

(l ) ( 0) t ~ Отметим, что v (inv),k, v (inv),k – недиагональные матрицы v (inv,k = Fv(,lk),( 0 ). Для случая трех l) 0 пробных матриц уравнения имеют вид 1 (0) Ф k (t ) + C (j0k) (t ) 0 Ф k (t ) + (k0 ) Ф k (t ) = 0, v inv,k t h j t S, (2.80) ( 0) t C j,k (t ) = j C j,k (t ) + j,k Ф k (t ).

( 0) (0) (0) 2.2 Одномерные уравнения МПГ На первом этапе данной работы проводился вывод нестационарных уравнений МПГ для одномерной плоской геометрии. Данная геометрия интересна тем, что в плоской одномерной ячейке возможны только две пробные матрицы – симметричная (k0) (r, ) и антисимметричная (k1) (r, ). Таким образом, набор пробных групповых функций образует полную систему решений. Как результат, в этой геометрии проще исследовать основные особенности применения МПГ для нестационарного случая, не отвлекаясь на возможную недостаточность пробных решений.

В данном параграфе приведем лишь результаты, т.к. процедура получения уравнений аналогична двумерному случаю. Конечно-разностные уравнения МПГ выводятся в процессе минимизации невязки аналогичной вышеописанной.

Уравнения МПГ имеют вид ( 0) d 2d v inv, k dt k (t ) + h dt C j, k (t ) 0 k (t ) + k k (t ) = (0) (0) (2.81) j d C ( 0) (t ) = C ( 0) (t ) + (0), j (t ) ( 0) (t ) dt j, k j j,k k k (1) d (1) 2 d (1) v inv, k dt X k (t ) + h dt C j, k (t ) 1k (t ) + k X k (t ) = ( 0 ) (1) (2.82) j d C(1) (t ) = C(1) (t ) dt j, k j j,k Отметим, что в одномерных уравнениях 2 (0) ( k (k1) ) 1, (k0) = (2.83) h [ ] h + т ( 0 ) ( k,0 ) dr +,т0( 0) (r ) d, j, k f (r ) (k0, 0 (r ) (k0 ).

(j0k = т ) ) (2.84) k, 2 Vk Система уравнений (2.81) является основной системой конечно-разностных уравнений МПГ для пространственной кинетики в одномерной геометрии.

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

2.3 Итерационная схема Далее будем рассматривать двумерный случай для трех пробных матриц (2.80). Итерационную схему будем строить, используя полностью неявную схему Эйлера по времени. Для дискретизации временных производных используется аппроксимацию первого порядка:

(t ) k (ts 1 ) d k (t ) k s, (2.85) ts dt C (t ) C j, k (ts 1 ) (0) (0) d (0) C j, k (t ) j, k s. (2.86) ts dt В соответствии с (2.85), (2.86) и фактом использования полностью неявной схемы получим 1 C ( 0 ) (t ) C (j0k (t s 1 ) v ( 0 ) (t ) k (t s ) k (t s 1 ) + 2 j, k s ) 0 (t s )Ф k (t s ) + k (t s )Ф k (t s ) =, (0) s S t s t s inv, k hj (2.87) C j, k (t s ) C j, k (t s 1 ) (0) (0) = j C (j0, k (t s ) + (j0, k (t s )Ф (k0 ) (t s ) ) ) t s Сведем систему уравнений (2.87) к системе линейных алгебраических уравнений ~ M (t s )(t s ) = S (t s 1 ), (2.88) ~ где M (t s ) - матрица, S(t s 1 ) - вектор.

Преобразуем первое уравнение системы (2.87) в следующую форму C( 0 ) (t ) 1 v inv), k (ts ) ( 0 (t s ) + (k0 ) (t s )k (ts ) + 2 j, k s = S ts Sh j ts (2.89) C(j0k (ts 1 ) v inv), k (ts )k (ts 1 ) ) ( = +, Sh j ts S t s Получим выражение для C(j0,k) (t s ) из второго уравнения системы (2.87) C (j0,k (t s 1 ) + t s (j0,k (t s )Ф (k0) (t s ) ) ) C (j0,k (t s ) = ). (2.90) 1 + t s j Подставим (2.90) в (2.89) и получим 1 v inv),k (t s ) ( 0 ) (t ) ( 0 (t s ) + (k0) (t s ) + 2 j,k s k (t s )+ = Sh j 1 + t s j S t s j v ( 0 ) ( t ) ( t ) 2 C(j0k) (t s 1 ) + inv,k s k s = Sh j 1 + t s j St s, ~ Тогда матрица M (t s ) и вектор S(ts 1 ) из (2.88) имеют вид (j0k (ts ) ) (0) 1 v inv, k (ts ) 2,, M (t s ) = 0 (ts ) + (k0) (ts ) + (2.91) Sh j 1 + ts j S ts j v ( 0 ) (t ) (t ) ~ 2 C (j0,k (t s 1 ) + inv,k s k s 1.

S (t s 1 ) = ) (2.92) Sh j 1 + t s j St s ( 0) Для того чтобы определить начальное условие для C(j0,k) (t ) положим C j, k (t ) = 0, t тогда 1 (0) C (j0,k (t 0 ) = ) j, k (t 0 )Ф (k0) (t 0 ). (2.93) j Блок-схема решения нестационарных уравнений МПГ представлена на рисунке 2.1. Временная точка tS – последняя временная точка расчета.

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

Для решения системы уравнений (2.88) может быть использована традиционная схема – «группа за группой». В соответствии с этим подходом система (2.88) может быть переписана в следующем виде g 1 G ~ M gg (t s ) (gn ) (t s ) = S g (t s 1 ) M gg (t s ) (gn ) (t s ) gg (t s ) (gn 1) (t s ).

M (2.94) g =1 g = g + где (gn ) (t s ) – пространственный вектор, соответствующий энергетической группе g, M gg (t s ) – пространственная матрица с фиксированными энергетическими группами g и g, n – номер внешней итерации.

Для того чтобы решить уравнение (2.94) использовался метод верхней релаксации. Матрица M gg (ts ) раскладывается на диагональную матрицу D gg (ts ), Рисунок 2.1 – Блок-схема решения нестационарных уравнений МПГ нижнюю треугольную и верхнюю треугольную матрицы Lgg (t s ) и U gg (t s ) соответственно:

M gg (t s ) = D gg (t s ) + Lgg (t s ) + U gg (t s ). (2.95) В итоге (gn) (ts ) итерационная схема может быть записана в виде ) [ L ( (t s ) (gn,m1) (t s ) U gg (t s ) (gn,m ) (t s ) + (gn,m ) (t s ) = D gg (t s ) gg (2.96) g 1 G ~ + S g (t s 1 ) M gg (t s ) (gn, ) (t s ) gg (t s ) (gn 1, ) (t s ).

M g =1 g = g + где m – номер внутренней итерации.

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

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

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

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

Также в данной главе проведено построение итерационной схемы для решения нестационарных конечно-разностных уравнений МПГ для трех пробных матриц. Данное построение проведено с использованием полностью неявной схемы Эйлера по времени. Для дискретизации временных производных использовалась аппроксимацию первого порядка. Для решения полученной системы линейных алгебраических уравнений предложено использование традиционной схемы «группа за группой» с применением метода верхней релаксации.

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

Глава 3 Программный комплекс SUHAM-TD и его верификация 3.1 Программный комплекс SUHAM-TD Данная глава посвящена описанию и верификации программного комплекса SUHAM-TD на ряде классических пространственно-временных бенчмарках.

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

1) Решаемая задача – уравнения пространственной кинетики;

2) Возможность проведение расчетов без диффузионного приближения;

3) Возможность проведения расчетов без пространственной гомогенизации;

4) Моделируемая среда – полномасштабная активная зона;

5) Высокая производительность.

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

Программный комплекс проводит решение уравнений SUHAM-TD пространственной кинетики посредством МПГ. Все модули SUHAM-TD разработаны на языке FORTRAN (стандарты 77 и 90). На текущий момент в программном комплексе реализованы нестационарные конечно-разностные уравнения МПГ для трех пробных матриц в квадратной решетке блоков. В качестве составных частей следует отметить следующие SUHAM-TD диффузионные модули, реализующие метод поверхностных гармоник:

1) SUHAMD-PL-1D-SHM для одномерной плоской геометрии в стационарном случае;

для одномерной плоской геометрии в 2) SUHAMD-TD-PL-1D-SHM нестационарном случае (использует результаты работы SUHAMD-PL-1D SHM);

3) SUHAMD-SQ-2D-SHM для квадратной решетки блоков в двумерной геометрии в декартовых координатах в стационарном случае;

4) SUHAMD-TD-SQ-2D-SHM для квадратной решетки блоков в двумерной геометрии в декартовых координатах в нестационарном случае (использует результаты работы SUHAMD-SQ-2D-SHM).

Программы и прошли SUHAMD-SQ-2D-SHM SUHAMD-TD-SQ-2D-SHM государственную регистрацию. Копии свидетельств о регистрации представлены в приложении А.

Диффузионные модули, реализующие классический метод конечных разностей имеют аналогичное наименование, за исключением трех последних букв (они заменяются на FDM), например SUHAMD-TD-SQ-2D-FDM. Область применения этих модулей аналогична.

Модули проводящие расчеты без диффузионного SUHAM-TD, приближения методом поверхностных гармоник:

1) SUHAM-SQ-2D-SHM для квадратной решетки блоков в двумерной геометрии в декартовых координатах в стационарном случае;

2) SUHAM-TD-SQ-2D-SHM для квадратной решетки блоков в двумерной геометрии в декартовых координатах в нестационарном случае (использует результаты работы SUHAM-SQ-2D-SHM).

Как видно, SUHAM-TD позволяет проводить расчеты, как с использованием диффузионного приближения, так и без него. Здесь важной чертой является то, что уравнения МПГ инвариантны по отношению к исходным уравнениям, из которых они получены, а именно, к уравнению переноса нейтронов и диффузионному уравнению и факт использования диффузионного приближения определяется только использованием его в расчете пробных матриц. Для получения диффузионных пробных матриц SUHAM-TD использует собственные процедуры, а для расчета пробных матриц без диффузионного приближения – программы РАЦИЯ [25] и DICPN [26]. Программа РАЦИЯ используется для расчета нулевой пробной матрицы методом поверхностных псевдоисточников в G3-приближении [109]. Программа DICPN применяется для получения первой и второй пробных матриц методом сферических гармоник в P2-приближении [27].

На рисунке представлена общая схема нестационарного расчета, 3. реализованная в SUHAM-TD, которая позволяет проследить место отмеченных программ в расчетном процессе (блок «расчет пробных матриц»). Также отметим, что факт использования процедуры пространственной гомогенизации определяется только программами расчета пробных матриц. Программы РАЦИЯ и DICPN не используют данную процедуру.

Рисунок 3.1 – Общая блок-схема нестационарного расчета по SUHAM-TD Программа разрабатывается для произвольного числа SUHAM-TD энергетических групп, групп запаздывающих нейтронов и числа расчетных зон.

Требование, которое накладывает МПГ на рассчитываемую среду – периодичность системы присуще большинству современных ядерных – реакторов.

Достигнуть большей производительности, чем существующие коды, предлагается за счет использования МПГ, который демонстрирует это в стационарном случае. Далее в данной главе в процессе верификации SUHAM-TD проводится сравнение вычислительных затрат МПГ по сравнению с традиционным конечно-разностным методом и демонстрируется потенциал МПГ в вопросе сокращения вычислительных затрат.

3.2 Верификация программного комплекса SUHAM-TD Верификация программного комплекса SUHAM-TD проводилась на ряде классических бенчмарков, характеризующихся отсутствием неопределенностей в начальных данных и отсутствием обратных связей. Этот факт дает возможность выделить методическую погрешность МПГ, реализованного в SUHAM-TD. Это особенно важно, т.к. работа посвящена не разработке кода, реализующего хорошо известный и верифицированный метод, а развитию МПГ для нестационарных расчетов, исследованию его применения, а также разработке нового кода на его основе и его верификации.

Важным фактом при проведении верификации SUHAM-TD является соблюдение «чистоты эксперимента» в плане сравнения вычислительных затрат.

Для этого в рамках SUHAM-TD был специально разработан модуль, проводящий расчеты классическим конечно-разностным методом (МКР) по пространству и времени. Для сравнения вычислительных затрат МПГ и МКР вычисления по SUHAM-TD проводились с одинаковыми критериями точности и способами организации итераций (например, во всех программах матрицы обращаются с помощью метода верхней релаксации, везде используется полностью неявная схема). Также в программах не применяются какие-либо методы ускорения.

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

Функционалы для сравнения. Расчеты, проведенные программным комплексом SUHAM-TD с помощью традиционного конечно-разностного метода, будем обозначать а расчеты, SUHAM-FDM (Finite-Difference Method), проведенные с помощью МПГ – SUHAM-SHM (Surface Harmonics Method). При верификации рассматривались следующие функционалы:

• начальное keff;

• полная (P) и зонные (Pi) мощности в выбранных точках по времени;

• распределение групповых потоков в выбранных точках по времени;

• вычислительные затраты.

Таким образом, проводилось сравнение следующих функционалов:

1. Относительное отклонение эффективного коэффициента размножения в стационарном расчете k keff k keff SHM k eff =, (3.1) k keff где k keff – эффективный коэффициент размножения, полученный с помощью SHM МПГ, k keff – эффективный коэффициент размножения, полученный другим методом;

2. Относительное отклонение мощности в различных временных точках P SHM (t ) P (t ) P (t ) =, (3.2) P (t ) где P SHM (t ) – относительная мощность в момент времени t, полученная с помощью МПГ, P(t ) – относительная мощность в момент времени t, полученная другим методом;

3. Максимальное по модулю относительное отклонение по пространству потоков группы g в различных временных точках iSHM (t ) i, g (t ) max g (t ) = max,g, (3.3) i, g (t ) где i, g – нейтронный поток в зоне i группе g полученный с помощью МПГ, SHM i, g – нейтроны поток в зоне i группе g полученный другим методом;

4. Среднеквадратическое отклонение RMS (root mean square) групповых потоков iSHM i, g N RMS =,g, (3.4) N i =1 i,g где N – число зон сравнения;

5. Относительные вычислительные затраты t FDM t = SHM, (3.5) t где t FDM - время расчета по программе SUHAM-FDM, t SHM -время расчета по программе SUHAM-SHM.

В качестве варьируемых величин, влияющих на значения выше приведенных функционалов, были выбраны: h – пространственный шаг сетки, он же – шаг сетки при расчете пробных матриц МПГ в ячейках, H – размер ячейки в расчетах МПГ, для которых проводится расчет пробных матриц, шаг временной сетки t. Надо отметить, что все сравнения проводились при условии, что пространственный шаг сетки h в расчете SUHAM-FDM равен пространственному шагу сетки h, с которым в расчете SUHAM-SHM рассчитываются пробные матрицы.

В программном комплексе SUHAM-TD Критерии сходимости.

используются следующие критерии сходимости итераций:

1. Для эффективного коэффициента размножения n k keff k keff n ;

(3.6) n k keff 2. Для потока нейтронов (расчет по МКР) во всех расчетных ячейках i во всех группах g in, g in, ;

g (3.7) in, g 3. Для уровня нейтронов (расчет по МПГ) во всех расчетных ячейках i (ячейка МПГ) во всех группах g in, g in,, g (3.8) in, g где n – номер итерации, i, g – поток нейтронов в g-ой группе i-ой зоне, i, g – n n уровень МПГ в g-ой группе i-ой ячейке определенный по формуле (2.50).

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

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

Таким образом, при расчете исходного стационарного состояния был выбран критерий сходимости = 10-9 для внутренних итераций по потоку (расчет SUHAM-FDM) и уровню (расчет SUHAM-SHM). Этот же критерий выбран для внешних итераций по потоку, уровню и собственному значению и для значение = 10- сверхвнешних итераций МПГ. В нестационарном случае используется только для внутренних итераций, а для остальных итераций выбран критерий = 10-7.

При расчете двумерных тестов, как в стационарном, так и в нестационарном случаях критерий = 10-9 использовался только для внутренних итераций, а для остальных итераций был выбран критерий = 10-7.

3.2.1 Тест BSS- Описание задачи. Как было отмечено в главе 2, на первом этапе данной работы проводился вывод одномерных нестационарных уравнений МПГ. Для проведения верификации данных уравнений был выбран одномерный тест BSS- (Benchmark Source Situation) [20]. В этом тесте описана трехзонная размножающая среда (рисунок 3.2), при этом области 1 и 3 являются физически эквивалентными.

Граничные условия представлены в виде нулевых потоков на обеих границах.

Для гомогенных областей этого теста в двухгрупповом приближении задан полный набор макросечений взаимодействия, а также коэффициенты диффузии (таблица 3.1). Среднегрупповые скорости в среде для всех зон равны 7 v1 = 1·10 см/с и v 2 = 3·10 см/с для первой и второй энергетических групп соответственно. Спектр деления для мгновенных и запаздывающих нейтронов также одинаков для всей среды и равен 1 = 1 и 2 = 0. Запаздывающие нейтроны представлены в шести группах (таблица 3.2). Отметим, что суммарная доля запаздывающих нейтронов равна = 0,0075.

Рисунок 3.2 – Геометрия теста BSS- Таблица 3.1 – Начальные двухгрупповые константы задачи BSS- f 1, f 2, a1, a2, 12, 21, D1, D2, Зона см-1 см- см см см-1 см-1 см-1 см- 1, 3 1,500 0,500 0,011 0,180 0,010 0,200 0,015 2 1,000 0,500 0,010 0,080 0,005 0,099 0,010 Таблица 3.2 – Параметры запаздывающих нейтронов задачи BSS- № группы 1 2 3 4 5 i 0,00025 0,00164 0,00147 0,00296 0,00086 0, i 0,0124 0,0305 0,1110 0,3010 1,1400 3, Рассматривались следующие варианты бенчмарка с разными законами ввода реактивности:

• A1 - сечение поглощения нейтронов в первой зоне во второй энергетической группе линейно увеличивается на 3% в течение первой секунды;

• A2 - сечение поглощения нейтронов в первой зоне во второй энергетической группе линейно уменьшается на 1% в течение первой секунды;

• A3 - сечение поглощения нейтронов в первой зоне во второй энергетической группе линейно уменьшается на 5% в течение первой 0,01 секунды;

• A4 - сечение поглощения нейтронов в первой зоне во второй энергетической группе линейно уменьшается на 5% в течение первой 0,01 секунды, при этом скорости нейтронов на два порядка больше.

В дополнении к этим тестам был проведен расчет теста A0 (его нет в [20]) представляющего собой задачу аналогичную A1, но с возмущением в 1 и 3 зонах (симметричный тест). Динамическое поведение исследуется на интервале 0 t 4 с для всех рассматриваемых случаев.

По условию бенчмарка начальная Стационарное моделирование.

конфигурация приводилась к критическому состоянию путем деления сечений f на keff (квазистационарная задача). В таблице 3.3 приведены сравнительные результаты расчета keff начального состояния по двум программам по SUHAM-TD с помощью МПГ и МКР.

Таблица 3.3 – Эффективный коэффициент размножения в начальном состоянии SUHAM-FDM SUHAM-SHM h, keff, % t keff H, keff см см 1 0,901598785 5,9e-6 16, 2 0,901598786 6,0e-6 50, 0,25 0, 5 0,901598820 9,8e-6 7, 10 0,901598670 6,9e-6 3, 2 0,901605604 2,2e-6 8, 5 0,901605608 2,7e-6 1, 0,5 0, 10 0,901605752 1,9e-5* 0, 5 0,901632279 2,4e-6 0, 1,00 0, 10 0,901632293 8,9e-7 0, 8 0,901732256 5,5e-7 0, 2,00 0, 10 0,901732174 9,7e-6 0, 2,50 0,901801756 10 0,901801755 1,1e-7 0, *) Здесь и далее полужирным шрифтом выделено максимальное отклонение.

Из таблицы 3.3 видно:

• сходимость расчета keff в расчетах SUHAM-FDM и SUHAM-SHM с уменьшением пространственного шага сетки h;

• значения keff, в расчете SUHAM-SHM с разными размерами ячеек практически совпадают со значениями, рассчитанными в SUHAM-FDM (максимальное keff = 1,9·10-5 %);

• наличие практической независимости рассчитанного по МПГ значения keff от размера ячейки H;

• преимущество МПГ по вычислительным затратам проявляется для достаточно значительных пространственных размерностей (для шага сетки равного 0,25 см и 0,5 см), когда уменьшение размерности крупносеточного расчета не компенсируется увеличением размерности задач расчета пробных матриц.

Сравнение средних по локальным зонам групповых потоков показало, что для h = 0,25 см max не превышает 4,4·10-4 %, а RMS – не превышает 2,3·10-4 %. Также было проведено сравнение зонных мощностей (Pi, i – номер гомогенной зоны соответствующей рисунку 3.2) в начальной конфигурации.

Данное сравнение проводилось при h = 0,25 см. Относительные отклонения для первой, второй и третьей зон составили P1 = 1,6·10-4 %, P2 = 1,1·10-4 % и P3 = 1,4·10-4 % соответственно.

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

k eff k eff as = =, (3.9) k eff k eff где – реактивность и эффективный коэффициент размножения, k eff соответственно в начальный момент времени,, keff – реактивность и эффективный коэффициент размножения соответственно в момент окончания ввода реактивности. В таблице 3.4 приведены результаты расчета k eff, k eff и as (в относительных единицах и долях запаздывающих нейтронов) с помощью метода поверхностных гармоник на сетке с параметрами h = 0,25 см, H = 2 см для различных вариантов теста BSS-6.

Таблица 3.4 – Величина вводимой асимптотической реактивности Вариант as, отн. ед. as, k eff keff A0 0,891146917 -1,30E-2 -1, A1 0,900963249 -7,82E-4 -0, 0, A2 0,904949965 4,11E-3 0, A3, A4 0,922209964 2,48E-2 3, Первым этапом для проведения верификации Выбор эталона.

нестационарных уравнений МПГ является обоснование эталона для сравнительного анализа. За эталон будем принимать расчет классическим конечно-разностным методом. Однако возникает вопрос о размере шага разбиения по пространству и по времени для эталонного расчета по МКР.

Рассмотрим обоснование диффузионного эталона на примере варианта А0.

Как показали расчеты мощностей, величины t = 10-4 с достаточно для временной сетки, т.к. максимальное отклонение зонной мощности для расчета с данным шагом по времени от расчета с шагом по времени равным t = 10-3 с составило 8·10-3 %. Данные расчеты представлены в таблице Б.1.

В таблице Б.2 представлены расчеты мощности (t = 10-2 с) при различных расчетных шагах по пространству. Из данной таблицы видно, что шаг пространственной сетки h = 0,25 см является достаточным, т.к. максимальное относительное отклонение при этом шаге от расчета с шагом h = 0,5 см составляет 6,5·10-3 %. Следует отметить, что дальнейшее уменьшение шагов t и h приводит к тому, что отклонение становится меньше при достаточно значительном росте временных затрат. Зависимости P от t и от h представлены на рисунке 3.3.

Результаты, аналогичные представленным, были получены также для вариантов A1 и A2. Таким образом, эталонным будем считать в дальнейшем SUHAM-FDM расчет с шагами t = 10-4 с, h = 0,25 см.

Для тестов BSS-6-А3 и BSS-6-А4 проведено аналогичное исследование для выбора эталона. Следует отметить, что тесты BSS-6-А3 и BSS-6-A4 отличаются от выше рассмотренных, тем, что процесс постепенно переходит из процесса на запаздывающих нейтронах в процесс на мгновенных нейтронах (см. таблицу 3.4), а тест BSS-6-A4 имеет к тому же скорости нейтронов на два порядка больше ( v1 = 1·109 см/с и v 2 = 3·107 см/с), т.е. задача более жесткая.

При проведении расчета были выбраны следующие расчетные = 10-6 с, h = 0,25 см для задачи BSS-6-A3 и характеристики эталона: t t = 10-7 с, h = 0,25 см для задачи BSS-6-A4.

Рисунок 3.3 – Максимальное относительное отклонение мощности в расчете SUHAM-FDM: (а) – в зависимости от временного шага (при h = 0,25 см), (б) – в зависимости от пространственного шага (при t = 10-2 с).

Нестационарное моделирование. Перейдем к сравнению МПГ с прямым численным моделированием в нестационарных случаях. В таблице 3.5 приведены значения максимальных отклонений зонных мощностей (зоны представлены на рисунке 3.2), полученных в расчете SUHAM-SHM, от значений, полученных в расчете SUHAM-FDM при одинаковых значениях расчетного шага по времени t с пространственными шагами h = 0,25 см, H = 2 см.

На примере варианта A0 можно сделать вывод, что сравнение расчета SUHAM-SHM при t = 10-4 c, с эталонным дает значение максимального относительного отклонения 4,3·10-2 %. Данное значение практически не зависит от шага по времени и лежит для остальных значений t в интервале 4,2–4,3·10-2 %. Следует отметить тот факт, что значение максимального относительного отклонения всегда соответствует последней временной точке и зоне, где происходит возмущение (см. таблица Б.3). Из данных таблиц 3.5 и Б. можно сделать вывод, что расчет SUHAM-SHM имеет незначительные отклонения от прямого численного моделирования (при одинаковых шагах по времени), однако, они возрастают по мере движения к последним расчетным точкам исследуемого временного интервала.

Таблица 3.5 – Максимальные значения среди относительных отклонений зонных мощностей Pi, полученных в расчетах SUHAM-SHM от значений в расчетах SUHAM-FDM (t для SUHAM-SHM равен t для SUHAM-FDM), % t, c A0 A1 A2 A3 A 10-7 – – – – 2,4E+ 10-6 – – – 1,0E+00 2,4E+ 10-5 – – – 1,0E+00 – 10-4 4,3E-02 3,8E-02 7,6E-02 – – 10-3 4,2E-02 3,8E-02 7,5E-02 – – 10-2 4,2E-02 4,0E-02 7,5E-02 – – 10-1 4,2E-02 4,0E-02 7,6E-02 – – Проведем анализ отклонений мощностей в зависимости от размера ячейки.

Данный анализ необязательно проводить при мелком разбиении временной сетки, поэтому расчет проведен с шагом t = 10-2 с для вариантов А0–А2, и с шагами t = 10-5 с и t = 10-6 с для вариантов А3 и А4 соответственно. В таблице 3. приведены значения максимальных отклонений зонных мощностей в зависимости от шага крупносеточного разбиения Н при h = 0,25 см.

Из данной таблицы видно, что отклонение существенно возрастает с ростом шага крупносеточного разбиения в варианте от (например, A 1,3·10-2 % до 1%). Также необходимо отметить, что погрешность от увеличения размера ячейки наиболее заметна в тестах А3 и А4, что указывает на необходимость использования более мелких ячеек в случае быстрого ввода реактивности.

Таблица 3.6 – Максимальные отклонения зонных мощностей SUHAM-SHM от SUHAM-FDM Pi в зависимости от шага крупносеточного разбиения Н, % H, см A0 A1 A2 A3 A 1 1,3E-02 1,4E-02 2,9E-02 1,9E-01 1,8E- 2 4,2E-02 4,0E-02 7,5E-02 1,0E+00 2,4E+ 5 2,5E-01 2,2E-01 8,1E-01 7,5E+00 1,6E+ 10 1,0E+00 8,9E-01 3,4E+00 – – Таблица Б.4 иллюстрирует данное исследование более подробно для варианта A0. В ней приведены значения полной и зонных мощностей SUHAM FDM, а также отклонения значений, полученных в расчете SUHAM-SHM, от значений в расчете Из таблицы видно, что значение SUHAM-FDM.

максимального относительного отклонения практически не зависит от шага h (например, при размере ячейки равном H = 5 см, отклонение составило при различных значениях геометрического расчетного шага h величину примерно равную 2,4·10-1 %).

Проведем анализ отклонений мощностей в зависимости от расчетного шага по времени. В таблице 3.7 представлены максимальные отклонения мощностей, полученных МПГ при разных значениях t, от эталонного расчета при h = 0, см, H = 2 см. Из таблицы видно, что максимальное отклонение практически всегда соответствует максимальному шагу по времени (например, для варианта A0 оно максимально при t = 10-1 с и составляет 4,8·10-1 %).

Таблица Б.5 представляет данное исследование более детально для варианта A0. В ней приведены значения полной и зонных мощностей (SUHAM-SHM, h = 0,25 см) при различных шагах по времени, а также отклонения данных значений от эталонных значений. Из этой таблицы видно, что максимальное отклонение для определенного шага по времени в данном случае не обязательно относится к последней расчетной точки расчетного интервала и зоне, где происходит возмущение (при больших значениях t).

Таблица 3.7 – Максимальные относительные отклонения зонных мощностей Pi SUHAM-SHM от эталонного расчета при различных шагах по времени, % t, c A0 A1 A2 A3 A 10-7 – – – – 2,4E+ 10-6 – – – 1,0E+00 2,3E+ 10-5 – – – 2,6E+00 – 10-4 4,3E-02 3,8E-02 7,6E-02 – – 10-3 4,1E-02 3,7E-02 6,3E-02 – – 10-2 8,4E-02 3,7E-02 1,2E-01 – – 10-1 4,8E-01 4,2E-01 1,6E+00 – – Далее перейдем к рассмотрению распределения групповых потоков нейтронов. На рисунках 3.4 – 3.6 представлено распределение плотности потока нейтронов в различных временных точках в обеих энергетических группах для SUHAM-SHM при h = 0,25 см, H = 2 см, t = 10-4 c для вариантов А0, А1 и А соответственно.

В таблице Б.6 приведены максимальные среднеквадратические отклонения и относительные отклонения групповых потоков (SUHAM-SHM, h = 0,25 см, H = 2 см) от эталонного расчета. Из данной таблицы видно, что максимальные значения max и RMS соответствует случаю с наибольшим шагом по времени, однако не наблюдается постоянного роста отклонений в процессе укрупнения временного шага. Например, в варианте A1 во второй энергетическое группе для шага t = 10-3с max = 5,6·10-2 % что немного меньше значения для t = 10-4 c, составляющего max = 5,9·10-2 %.

Более детальные данные, посвященные этому исследованию, можно найти в таблице Б.7, представляющей значения только по варианту A0. Из нее можно видеть, что для малых временных шагов (t = 10-4 c, t = 10-3 c) значения max и RMS возрастают с ростом времени и их максимумы соответствуют последней Рисунок 3.4 – Тест BSS-6-A0. Пространственное распределение потоков в расчете SUHAM-SHM в первой (слева) и второй (справа) группах Рисунок 3.5 – Тест BSS-6-A1. Пространственное распределение потоков в расчете SUHAM-SHM в первой (слева) и второй (справа) группах Рисунок 3.6 – Тест BSS-6-A2. Пространственное распределение потоков в расчете SUHAM-SHM в первой (слева) и второй (справа) группах временной точке. Данная тенденция не наблюдается для временных шагов t = 10-2 c и t = 10-1 c. Тенденции RMS(t) для варианта A1 проиллюстрированы на рисунке 3.7.

Рисунок 3.7 – Тест BSS-6-A1. Среднеквадратические отличия пространственных распределений потоков во второй группе для SUHAM-SHM (t =10-1 c, 10-2 c, 10-3 c и 10-4 c) от SUHAM-FDM с шагом по времени равным t = 10-4 c Из данных графиков виден темп роста RMS с ростом времени для малых временных шагов. Следует обратить внимание на резкие скачки RMS для случаев t =10-2 c, t =10-1 c. Данные скачки обусловлены погрешностью, связанной с дискретностью ввода реактивности, и берут свое начало в момент этого ввода.

Скачок для случая t = 10-1 c более большой и его снижение продолжается до конца расчетного временного интервала. Скачок RMS для случая t = 10-2 c менее значителен и спадает практически сразу. Таким образом, можно сделать вывод, что ширина и амплитуда пиков RMS(t) обусловлена значительной погрешностью при небольших t, а их местоположение – началом ввода реактивности.

3.2.2 Тест PHWR Описание задачи. Тест PHWR предложен в 1977 в работе [110]. Данная проблема представляет собой идеализированное представление реактора CANDU, где переходный процесс моделирует аварию с потерей теплоносителя (LOCA) с последующим вводом поглощающих устройств. Тест имеет двумерную расчетную среду, представляющую собой четыре параллельно расположенные зоны одинакового размера (рисунок 3.8). Граничные условия представлены в виде нулевого потока на всех четырех границах.

Рисунок 3.8 – Геометрия теста PHWR Энергетическое распределение нейтронов представлено в данной задаче в виде двух энергетических групп. Нейтронно-физические константы представлены в таблице 3.8. Особенностью данного теста является то, что запаздывающие нейтроны в задаче представлены в двух вариантах: одна группа ( = 0,00572 и = 0,08 с-1) и шесть групп (таблица 3.9). Далее будем обозначать их как PHWR-1 и PHWR-6 соответственно. Среднегрупповые скорости в среде для всех зон равны 6 v1 = 7,16·10 см/с и v 2 = 2,86·10 см/с для первой и второй энергетических групп соответственно. Спектр деления для мгновенных и запаздывающих нейтронов также одинаков для всей среды и равен 1 = 1 и 2 = 0.

Таблица 3.8 – Свойства материалов для задачи PHWR f 1, f 2, D2, a1, a2, 12, D1, Зона -1 - см-1 - см- см см см см см 1и4 1,2 0,9 0,001 0,004 0,009 0 0, 2и3 1,2 0,9 0,001 0,004 0,009 0 0, Таблица 3.9 – Параметры шести групп запаздывающих нейтронов Группа 1 2 3 4 5 - i, 10 0,726620 0,819557 1,041405 2,335200 0,579031 0, - i, с 0,0137856 0,0305098 0,1336881 0,3155434 1,2250770 3, Переходный процесс развивается как следствие линейного роста f 2 на 0,5% в течение одной секунды для всей среды. На следующем этапе происходит линейное увеличение теплового сечения поглощения на 12,5% в течение 1,9 с в первой и второй зонах. Динамическое поведение исследуется на интервале 0 t 2,9 с.

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

В дополнение к приведенным значениям, отметим результаты расчетов с пространственными шагами 12,5 см и 100 см по осям абсцисс и ординат соответственно, приведенные в работе [110]: для программы ADEP – keff = 1,007489, для программы CERKIN – keff = 1,007484.

Исходя из данных расчетов, можно принять значение параметра h равное см для получения нестационарного эталонного расчета. Данное решение обусловлено незначительным отклонением keff для h = 4 см и h = 5 см, которое составило 3,7·10-5 %. В качестве оптимального значения параметра H для МПГ выбрана величина в 20 см, т.к. отклонение keff для H = 20 см составило 1,6·10-3 % при h = 4 см.

Таблица 3.10 – Эффективный коэффициент размножения keff в зависимости от параметров расчетной сетки SUHAM-FDM SUHAM-SHM h, см keff, % H, см keff keff 20 1,007429 1,62E- 4 1, 40 1,007480 6,72E- 20 1,007429 1,58E- 5 1,007413 40 1,007480 6,69E- 50 1,007519 1,06E- 8 1,007414 40 1,007480 6,53E- 40 1,007480 6,37E- 10 1, 50 1,007519 1,02E- 12,5 1,007418 50 1,007519 1,00E- Отметим, что среднеквадратичное отклонение среднего по ячейке потока в расчете SUHAM-SHM (h = 4 см и H = 20 см) от расчета SUHAM-FDM (h = 4 см) составило 0,15 % как для быстрой, так и для тепловой групп.

Максимальное относительное отклонение для средних по ячейке потоков для обеих энергетических групп составило 0,20 %.

Приведем значения асимптотической реактивности, рассчитанной по формуле (3.9). В таблице 3.11 приведены результаты расчета as с помощью МПГ на сетке с параметрами h = 4 см, H = 20 см на различных этапах ввода реактивности. Отметим, что приведенные значения справедливы для PHWR-1 и несмотря на незначительную разницу в относительной доли PHWR- запаздывающих нейтронов.

Определим шаг временной сетки для эталонного расчета. Для этого проведем расчет задачи посредством прямого конечно-разностного моделирования с различными временными шагами на выбранной пространственной сетке с шагом h = 4 см. Результаты данных расчетов приведены в таблицах В.1 и В.2. Из них можно видеть, что максимальные отклонения мощности между расчетами с временными шагами t = 10-5 с и t = 10-4 с составили -7,9·10-3 % и -7,2·10-3 % для PHWR-1 и PHWR-6 соответственно. Как результат, можно считать расчет с шагом временной сетки t = 10-5 с эталонным для дальнейшего анализа.

Таблица 3.11 – Величина вводимой асимптотической реактивности Временной as, отн. ед. as, k eff keff интервал 0–1c 1,007429 1,012466 4,9E-3 0, 1 – 2,9 c 1,012466 1,005713 6,6E-3 -1, 0 – 2,9 c 1,007429 1,005713 -1,7E-3 -0, Таким образом, в дальнейшем в качестве эталона принимаем расчет SUHAM-FDM, проведенный с помощью прямого численного моделирования на пространственной сетке с шагом h = 4 см и временной сетке c шагом t = 10-5 с.

Проанализируем отклонение мощности Нестационарный расчет.

полученной с помощью МПГ при шагах временной сетки 10-4 с, 10-3 с и 10-2 с по сравнению с результатами для t = 10-5 с. Максимальные относительные отклонения мощности в течение временного расчета SUHAM-SHM с различными шагами по времени от расчета с шагом t = 10-5 c при h = 4 см и H = 20 см представлены в таблице 3.12. Из данной таблицы можно видеть, что в обоих случаях погрешность падает с уменьшением шага временной сетки примерно на порядок, что свидетельствует о сходимости МПГ с уменьшением шага временной сетки. Детально по временным точкам отклонения представлены в таблицах В.3 и В.4.

Теперь проведем сравнение непосредственно МПГ и прямого конечно разностного моделирования. На рисунке 3.9 представлены графики изменения мощности с течением времени для обоих вариантов теста при t = 10-5 с, h = 4 см, H = 20 см. Из данных графиков видно, что визуально результаты по МПГ и прямому численному моделированию совпадают для обоих случаев.

Таблица 3.12 – Максимальное относительное отклонение мощности SUHAM SHM с различными шагами по времени от расчета SUHAM-SHM с шагом t=10-5 c, % t=10-4 c t=10-3 c t=10-2 c Вариант PHWR-1 -5,57E-03 -6,99E-02 -6,90E- PHWR-6 8,11E-03 7,17E-02 7,13E- Рисунок 3.9 – Изменения мощности в задаче PHWR Рассмотрим численные результаты. В таблице приведены 3. максимальные относительные отклонения мощности SUHAM-SHM (h = 4 см, H = 20 см) при различных временных шагах от расчетов SUHAM-FDM (h = 4 см) с различными временными шагами (сравнение проводится для одинаковых временных шагов).

Из данных таблиц видно, что от размера временного шага отклонение практически не изменяется и максимальная величина составляет 0,76 % и 0,84 % для одной и шести групп запаздывающих нейтронов соответственно. Эта величина характеризует степень приближения МПГ (в данном случае три пробных матрицы). Следует заметить, что данные отклонения относятся к точке с, которая соответствует смене знака вводимой реактивности. В остальных временных точках данное отклонение ниже и в последней временной точке составляет -0,35 % и -0,38 % для одной и шести групп запаздывающих нейтронов соответственно. В таблицах В.5 и В.6 относительные отклонения мощности приведены детально по временным точкам.

Таблица 3.13 – Максимальные относительные отклонения мощности SUHAM SHM от мощности по SUHAM-FDM при сохранении шага временной сетки, % t=10-5 c t=10-4 c t=10-3 c t=10-2 c Вариант PHWR-1 -7,58E-01 -7,54E-01 -7,48E-01 -7,45E- PHWR-6 -8,40E-01 -8,34E-01 -8,28E-01 -8,22E- Далее сопоставим данные SUHAM-SHM при различных временных шагах с эталонным расчетом (таблицы 3.14 и 3.15). Из данных таблиц видно, что максимальное относительное отклонение слабо меняется при изменении временного шага (за исключением грубой временной сетки с t = 10-2 с) и колеблется от 0,75 % до 0,77 % по абсолютной величине для PHWR-1 и от 0,83 % до 0,86 % по абсолютной величине для PHWR-6. В значительной мере временная дискретизация проявляет себя для временной сетки с шагом t = 10-2 с, когда относительное отклонение составляет -1,39% и -1,33% для PHWR-1 и PHWR- соответственно.

Таблица 3.14 – PHWR-1. Относительное отклонение мощности SUHAM-SHM при различных шагах временной сетке от значений в эталонном расчете Отклонение SUHAM-SHM от SUHAM-FDM SUHAM-FDM P, отн. ед. P, % t, c t = 10-5 c t = 10 c t = 10 c t = 10-3 c t = 10-2 c -5 - 0 1,00000E+00 - - - 0,10 1,02741E+00 -3,52E-01 -3,50E-01 -3,36E-01 -2,04E- 0,50 1,40688E+00 -7,06E-01 -7,03E-01 -6,86E-01 -5,09E- Продолжение таблицы 3. Отклонение SUHAM-SHM от SUHAM-FDM SUHAM-FDM P, отн. ед. P, % t, c t = 10-5 c t = 10 c t = 10 c t = 10-3 c t = 10-2 c -5 - 1,00 2,71885E+00 -7,58E-01 -7,53E-01 -7,07E-01 -2,45E- 1,25 1,83652E+00 -7,02E-01 -7,08E-01 -7,72E-01 -1,39E+ 1,50 9,92125E-01 -5,54E-01 -5,52E-01 -5,36E-01 -3,81E- 2,00 5,96599E-01 -4,54E-01 -4,52E-01 -4,40E-01 -3,18E- 2,50 5,03481E-01 -3,95E-01 -3,94E-01 -3,93E-01 -3,83E- 2,90 4,64385E-01 -3,51E-01 -3,50E-01 -3,50E-01 -3,54E- Таблица 3.15 – PHWR-6. Относительное отклонение мощности SUHAM-SHM при различных шагах временной сетки от значений в эталонном расчете Отклонение SUHAM-SHM от SUHAM-FDM SUHAM-FDM P, отн. ед. P, % t, c t = 10-5 c t = 10 c t = 10 c t = 10-3 c t = 10-2 c -5 - 0 1,00000E+00 – – – – 0,10 1,02746E+00 -3,53E-01 -3,51E-01 -3,38E-01 -2,04E- 0,50 1,41915E+00 -7,41E-01 -7,37E-01 -7,15E-01 -4,90E- 1,00 2,88514E+00 -8,40E-01 -8,32E-01 -7,69E-01 -1,33E- 1,25 2,05428E+00 -8,03E-01 -8,06E-01 -8,56E-01 -1,33E+ 1,50 1,16925E+00 -6,67E-01 -6,64E-01 -6,54E-01 -5,52E- 2,00 7,01766E-01 -5,29E-01 -5,26E-01 -5,18E-01 -4,39E- 2,50 5,68420E-01 -4,42E-01 -4,39E-01 -4,39E-01 -4,30E- 2,90 5,09428E-01 -3,81E-01 -3,78E-01 -3,78E-01 -3,79E- Проведем теперь анализ результатов для плотности потока нейтронов. В таблицах 3.16 и 3.17 представлены данные для средней плотности потока нейтронов во всей размножающей среде, полученные с помощью прямого конечно-разностного моделирования, МПГ, а также результаты расчетов по программам ADEP и CERBERUS/CERKIN (далее для краткости CERKIN) из [110]. Программа ADEP проводит расчеты прямым конечно-разностным моделированием используя ADE схему (alternating direction explicit method) с экспоненциальным разложением, а программа CERKIN по методу IQS (improved quasi-static method). Расчеты по обоим кодам (CERKIN и ADEP) проводились на пространственной сетке x = 12,5 см y = 100 см и временной сетке t = 2,5·10-3 с.

Также в таблицах 3.16 и 3.17 представлены относительные отклонения результатов по SUHAM-SHM от результатов от SUHAM-FDM. Относительное отклонение для среднего теплового потока имеет похожее поведение с относительным отклонением мощности.

Таблица 3.16 – PHWR-1.Средний тепловой поток, отн. ед.

, % t, c CERKIN ADEP SUHAM-FDM SUHAM-SHM 0 1,000 1,000 1,000 1,000 – 1,00 2,694 2,695 2,705 2,685 -0, 1,25 1,654 1,633 1,827 1,814 -0, 1,50 0,879 0,868 0,987 0,981 -0, 1,85 0,614 0,605 0,648 0,645 -0, 2,20 0,533 0,525 0,546 0,543 -0, 2,90 0,459 0,452 0,461 0,460 -0, Таблица 3.17 – PHWR-6. Средний тепловой поток, отн. ед.

, % t, c CERKIN ADEP SUHAM-FDM SUHAM-SHM 0 1,000 1,000 1,000 1,000 – 1,00 2,866 2,867 2,871 2,847 -0, 1,50 1,088 1,022 1,163 1,155 -0, 2,00 – 0,654 0,698 0,694 -0, 2,50 0,592 0,543 0,565 0,562 -0, 2,90 0,525 0,491 0,506 0,504 -0, Проведем сравнение средних потоков нейтронов по ячейке. В таблицах В. и В.8 представлены результаты сравнения средних по ячейке потоков нейтронов в виде модуля максимального отклонения и модуля среднеквадратичного отклонения. Из данных таблиц видно, что для среднеквадратичного отклонения наблюдается выход на плато со значением примерно 0,51 % и 0,68 % для PHWR- и PHWR-6 соответственно после скачка до величины 0,74 % и 0,85 % соответственно после смены знака вводимой реактивности. Данное поведение и характерные величины присущи обеим энергетическим группам. Что касается величины max, то она растет с величиной расчетного времени, также выходя на плато. Значение max в последней временной точке PHWR-1 составило 1,08 % (быстрая группа) и 1,19 % (тепловая группа). В случае PHWR-6 данная величина принимает значения 1,12 % и 1,23 % для быстрой и тепловой групп соответственно. Из таблиц В.7 и В.8 видно, что величина скачка max после смена знака вводимой реактивности незначительна. На рисунке 3. представлены графики плотности потока нейтронов в нижней полуплоскости во второй группе в моменты времени 0 с, 1,25 с и 2,9 с соответственно.

Рисунок 3.10 – Плотность потока тепловых нейтронов (МПГ, h = 4 см, H = 20 см, t = 10-5 с). Сверху вниз: 0 с, 1,25 с и 2,9 с В таблице 3.18 представлены относительные вычислительные затраты t для SUHAM-SHM (h = 4 см и H = 20 см) и SUHAM-FDM (h = 4 см). Видно, что вычислительные затраты для МПГ значительно меньше чем для прямого конечно разностного моделирования. Следует отметить, что показатели в таблице в значительной мере зависят от шага временной сетки.

Таблица 3.18 – Относительные вычислительные затраты t, отн. ед.

t, с PHWR-1 PHWR- 1,0E-03 32,3 22, 1,0E-04 10,3 9, 1,0E-05 4,5 1, 3.2.3 Тест TWIGL Описание задачи. В качестве другой тестовой проблемы был использован двумерный бенчмарк TWIGL [19]. Этот тест был предложен в 1969 году, но по прежнему активно используется для верификации расчетных кодов, что позволяет обеспечить сравнение с другими расчетными кодами и методами. Данная задача представляет собой идеализированное представление реактора типа LWR.

Геометрия задачи представлена в виде квадрата 80х80 см. Тест имеет регионы с тремя различными составами: возмущаемая зона с повышенным содержанием делящихся веществ (зона 1), невозмущаемая зона с повышенным содержанием делящихся веществ (зона 2) и зона с пониженным содержанием делящихся веществ (зона 3). Зоны 1 и 2 имеют идентичные константы в начальный момент времени. Расчетная зона с граничными условиями представлена на рисунке 3.11.

Нейтронно-физические характеристики среды представлены в виде двухгрупповых макроконстант (таблица 3.19). Среднегрупповые скорости нейтронов в тесте одинаковы для всей среды и составляют v1 = 107 см/с и v2 = 2·105 см/с. Спектр деления для мгновенных и запаздывающих нейтронов также одинаков для всей среды и равен 1 = 1 и 2 = 0.

Рисунок 3.11 – Геометрия теста TWIGL Таблица 3.19 – Нейтронно-физические константы для теста TWIGL f 1, f 2, D2, a1, a2, 12, D1, Зона s см см см-1 -1 - см см - см- см 1, 2 1,4 0,4 0,01 0,15 0,01 0,007 0, 3 1,3 0,5 0,008 0,05 0,01 0,003 0, В тесте представлена одна эффективная группа запаздывающих нейтронов.

Эффективная доля запаздывающих нейтронов составляет = 0,0075, а постоянная распада предшественников запаздывающих нейтронов = 0,08 с-1.

В задаче представлены два сценария развития переходного процесса:

возмущение скачком и линейное возмущение в зоне 1, далее TWIGL-S и TWIGL R соответственно. Возмущение скачком инициируется посредством уменьшения теплового сечения деления a2 на 0,0035 см-1 в нулевой момент времени. В случае линейного возмущения a2 линейно уменьшается на 0,0035 в течение 0,2 с.

Динамическое поведение исследуется на интервале 0 t 0,5 с для обоих случаев.

Стационарный расчет и выбор эталона. Для того чтобы получить оптимальную пространственную сетку стационарные расчеты были проведены для различных пространственных разбиений. Таблица 3.20 объединяет результаты для различных стационарных расчетов. Основываясь на результатах таблицы 3.20, далее для расчета этого бенчмарка будет использоваться пространственная сетки с параметрами h = 0,8 см и H = 4 см.

Таблица 3.20. Эффективный коэффициент размножения в зависимости от параметров расчетной сетки.

SUHAM-FDM SUHAM-SHM h, см k eff,% H, см k eff k eff 2 0,913184 -2,3E- 0,4 0,913204 4 0,913218 -1,5E- 8 0,913597 -4,3E- 4 0,913204 -1,8E- 0,8 0, 8 0,913581 -4,3E- 1,0 0,913176 8 0,913571 -4,3E- 2,0 0,913101 8 0,913495 -4,3E- На рисунке 3.12 – представлены плотности потока нейтронов, полученные с помощью МПГ на выбранной сетке в стационарном случае, для быстрой и тепловой групп. Среднеквадратичное отклонение среднего по ячейке потока от расчета на выбранной сетке составило 0,22 % для быстрой группы и 1,03 % для тепловой группы. Максимальное относительное отклонение для средних по ячейке потоков для быстрой группы составило 0,55 % и 3,29 % для тепловой группы.

Приведем значения асимптотической реактивности as, рассчитанной по формуле (3.9) с помощью МПГ на сетке с параметрами h = 0,8 см, H = 4 см.

Эффективный коэффициент размножения в возмущенной среде составил k eff = 0,916716. Таким образом, as = 4,2E-3 или as = 0,56.

При выборе шага временной сетки для получения эталонного решения были проведены расчеты с помощью прямого конечно-разностного моделирования. Эти расчеты были проведены с временными шагами от 10-2 с до 10-5 с. Расчет с временным шагом t = 10-5 с был выбран в качестве эталонного, т.к.

максимальное отклонение мощности между расчетом с шагом 10-5 с и расчетом с шагом 10-4 с составило 2,3·10-3 % для случая линейного возмущения. Будем считать, что данная точность является достаточной для дальнейших исследований. Относительные отклонения мощностей, полученных прямым конечно-разностным моделированием для различных временных шагов для обоих сценариев ввода реактивности, приведены в таблицах Г.1 и Г.2.

Таким образом, эталонным будем считать расчет, проведенный прямым конечно-разностным методом с параметрами h = 0,8 с, H = 4 см и t = 10-5 с.

Рисунок 3.12 – Стационарный быстрый (слева) и тепловой (справа) потоки для теста TWIGL Нестационарный расчет. Для данной задачи была проведена серия нестационарных расчетов без возмущения, как МПГ, так и прямым конечно разностным методом для того, чтобы проверить отсутствие изменения мощности при расчете без возмущения в рамках заданных критериев сходимости.

Рассматривался временной интервал от 0 с до 0,5 с. Расчеты были выполнены для следующих временных шагов: 10-3 с, 10-4 с и 10-5 с. Для всех случаев не получено никаких отклонений мощности от начального уровня для установленных критериев сходимости, что говорит об устойчивости алгоритмов, реализованных в обоих расчетах.

Проведем сравнение с эталонным расчетом. Расчет по МПГ был проведен с различными временными шагами и результаты были сопоставлены с эталонным решением. Максимальные относительные отклонения мощности SUHAM-SHM при различных временных шагах от эталонного расчета собраны в таблице 3.21.

Подробно по временным точкам эти данные представлены в таблицах Г.3 и Г. для TWIGL-R и TWIGL-S соответственно.

Таблица 3.21 – Максимальные относительные отклонения мощности SUHAM SHM от значений в эталонном расчете t = 10-5 c t = 10-4 c t = 10-3 c Вариант TWIGL-R 4,72E-02 4,78E-02 5,59E- TWIGL-S -9,54E-02 -1,41E-01 -9,81E- Проанализируем случай с линейным возмущением. В таблице Г.3 можно видеть, что в течение переходного процесса расчет SUHAM-SHM хорошо согласуется с расчетом SUHAM-FDM для временных шагов t = 10-5 с и t = 10-4 с. Для этих случаев относительное отклонение постепенно возрастает в течение переходного процесса. Однако, для более грубой временной сетки (t = 10-3 с) наблюдается два пика относительного отклонения после начала ввода реактивности и после окончания ввода реактивности. Следует заметить, что расчеты по МПГ дают большее значение мощности, чем конечно-разностное моделирование и максимальное относительное отклонение составляет 4,7·10-2 % и это отклонение соответствует последней временной точке t = 0,5 с.

Далее рассмотрим случай с мгновенным возмущением. Для этого случая наблюдаются пики относительного отклонения после мгновенного ввода реактивности для расчетов с любым шагом временной сетки (таблице Г.4).

Можно видеть, что увеличение временного шага приводит к значительному росту относительного отклонения для начальных временных точек. После этого наблюдается стабилизация относительного отклонения. Для случая с мгновенным возмущением максимальное относительное отклонение составило примерно -0,1 % для шага по времени t = 10-5 с.

В таблице Г.5 представлены отклонения для плотности потока нейтронов варианта TWIGL-R. Величина практически постоянна в течение max расчетного времени и равна 0,55 % и 3,3% для быстрой и тепловой групп соответственно. Среднеквадратичное отклонение RMS равно 0,22 % и 1,02 % для быстрой и тепловой групп соответственно. В таблице Г.6 представлены отклонения для плотности потока нейтронов варианта TWIGL-S. Величины max имеет пик 0,68 % для первой энергетической группы в начале расчета с последующим выходом на плато равное 0,54 %. В случае тепловой группы max не превышает 3,3 %. Среднеквадратичное отклонение RMS постоянно для тепловой группы и равно 1,02 %. В случае быстрой группы наблюдается пик до 0,29 % с последующим выходом на плато 0,22 %.

Как уже отмечалось, данный тест по-прежнему активно используется различными научными коллективами для верификации расчетных кодов, что позволяет обеспечить сравнение с другими расчетными кодами и методами. На рисунках 3.13 и 3.14 представлены данные из таблиц Г.7 и Г.8, в которых собраны результаты расчета теста TWIGL различными программами. На рисунках 3.13 и 3.14 имеются следующие обозначения: SUHAM-FDM – расчет SUHAM-FDM выполненный для временного шага t = 10-5 с и пространственного шага h = 0,8 см, SUHAM-SHM – расчет SUHAM-SHM с параметрами t = 10-5 с, h = 0,8 см и H = 4 см, MGNEM, NEM, AFEN, SPANDEX, JAR-IQS – результаты расчетов по одноименным кодам из источников [89], [111], [112] и [113] соответственно (результат по JAR-IQS получен автором этой программы).

Из вышеперечисленных работ доступна следующая информация о представленных кодах и результатах. Коды MGNEM, NEM и SPANDEX проводят расчеты с помощью метода NEM (Nodal Expansion Method). Расчет по MGNEM проводился с параметрами t = 10-2 с и с шагом пространственной сетки 8 см.

Расчет по коду NEM проводился на временной сетке t = 2,5·10-3 с. Код AFEN реализует и AFEN-метод (Analytic Function Expansion Nodal method) представленные по нему результаты получены на ячеечной сетке с размером нода 1х1 см с шагом по времени 10-3 с. Программа JAR-IQS использует в расчетах метод улучшенный квазистатический метод (IQS) для временной задачи. Данный расчет проведен на временной сетке t = 10-4 с.

Рисунок 3.13 – Изменение мощности с течение времени для TWIGL-R Рисунок 3.14 – Изменение мощности с течение времени для TWIGL-S Из представленных результатов по сторонним кодам можно заключить, что расчеты по SUHAM-TD хорошо согласуются с расчетами по другим кодам.

Особенно следует отметить хорошее согласие с результатами расчетов по кодам AFEN, SPANDEX и JAR-IQS.

В таблице 3.22 представлены относительные вычислительные затраты для обоих сценариев ввода реактивности, где tSHM – время на выполнение расчета SUHAM-SHM на пространственной сетке с параметрами h = 0,8 см и H = 4 см, tFDM – время расчета SUHAM-FDM с пространственной сеткой h = 0,8 см. Видно, что вычислительные затраты для МПГ значительно меньше чем для прямого конечно-разностного моделирования. Следует отметить, что показатели в таблице в значительной мере зависят от шага временной сетки.

Таблица 3.22 – Относительные вычислительные затраты t, отн. ед.

t, с TWIGL-R TWIGL-S 1,0E-03 40,1 35, 1,0E-04 11,2 9, 1,0E-05 3,5 4, Отметим, что показанная зависимость вычислительных затрат от шага временной сетки наблюдается во всех предыдущих задачах. Данный факт объясняется следующей причиной. С уменьшением временного шага обращение матрицы при решении реакторного уравнения производится быстрее как в случае прямого конечно-разностного моделирования, так и в случае МПГ по причине меньшего возмущения элементов в матрице. Однако в случае МПГ часть времени уходит на расчёт пробных матриц, продолжительность которого, как показывает расчёт этих тестов, слабо зависит от размера временного шага.

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

Например, в тесте TWIGL при текущей организации комплекса SUHAM-TD в течение первых 0,2 с при шаге t = 10-3 с отношение времени на расчёт пробных матриц к полному времени расчёта составляет 9 %, а при шагах t = 10-4 с и t = 10-5 с оно составляет 35 % и 49 %, соответственно.

3.2.4 Модифицированный тест 8-А Описание задачи. Тест 8-A1 [20], также получивший распространение в литературе под названием также был использован для «Dodds»-тест, верификации. Задача 8-А1 широко применяется по настоящее время для отладки кодов на начальном этапе разработки программ (например [85]). Здесь рассматривается модификация данного теста, заключающаяся в том, что координата r заменена на х, а координата z - на y. Интерес к расчету данного теста вызван тем фактом, что он является геометрически более масштабным.

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

Таблица 3.23 – Нейтронно-физические константы для теста 8-A s g, out, g *, g f, g, Dg, Зона Группа см см- см-1 см- 1, 15 1 1,0684 2,8E-2 0 2,6E- 2 3,2051E-1 3,3E-3 0 2, 14 1 1,3495 1,201E-2 0 1,2E- 2 8,7032E-1 1,9E-2 0 3, 4, 11 1 1,3052 1,0475E-2 1,1776E-2 8,0351E- 2 8,8857E-1 1,3063E-2 1,3268E-2 5, 12 1 1,3052 1,0475E-2 1,1776E-3 8,0351E- 2 8,8857E-1 1,2623E-2 1,3268E-2 6, 13 1 1,3052 1,0475E-2 1,1776E-3 8,0351E- 2 8,8857E-1 1,1183E-2 1,3268E-2 7, 8 1 1,3052 1,0475E-2 1,1776E-3 8,0351E- 2 8,8857E-1 1,3453E-2 1,3268E-2 Продолжение таблицы 3. f, g, s g, out, g *, g Dg, Зона Группа см см- см-1 см- 9 1 1,3052 1,0475E-2 1,1776E-3 8,0351E- 2 8,8857E-1 1,2973E-2 1,3268E-2 10 1 1,3052 1,0475E-2 1,1776E-3 8,0351E- 2 8,8857E-1 1,2933E-2 1,3268E-2 16 1 1,2997E+1 1,047E-2 1,2875E-3 7,9061E- 2 8,7951E-1 1,3065E-2 1,4246E-2 * out, g = с, g + f, g + sg g.

Рисунок 3.15 – Геометрия теста 8-A Запаздывающие нейтроны в тесте представлены в шести группах (таблица 3.24). Скорости нейтронов для всех зон одинаковые и равны v1 = 1·107 см/с и v2 = 2,2·105 см/с для первой и второй энергетических групп соответственно. Спектр мгновенных и запаздывающих нейтронов представлен как 1 = 1, 2 = 0.

Таблица 3.24 – Характеристики запаздывающих нейтронов задачи BSS-8-A Группа 1 2 3 4 5 i 2,470E–4 1,384E–3 1,222E–3 2,645E–3 8,320E–4 1,690E– - i, с 1,27E–2 3,17E–2 1,15E–1 3,11E–1 1,40 3, В тесте представлен следующий сценарий развития переходного процесса: в зонах 3 и 7 тепловое сечение поглощения линейно возрастает на 3%, а в зоне линейно уменьшается на 3% в течение 1 с. Переходный процесс исследуется в интервале 0 с t 4 с.

В данной задаче не проводилось выбора Результаты моделирования.

оптимальной пространственной сетки посредством сравнения расчетов различных вариантов. Прямой конечно-разностный расчет был проведен с пространственным шагом h = 1,25 см, а расчет по МПГ с шагами H = 2,5 см и h = 1,25 см. Отметим, что дальнейшее уменьшение пространственного шага ведет к значительному увеличению расчетного времени. Критерии точности, как для эффективного коэффициента размножения, так и для плотности потока нейтронов составили 10-7 %. Эффективный коэффициент размножения для SUHAM-FDM составил keff = 0,864095 и keff = 0,864101 для SUHAM-SHM. Таким образом, относительное отклонение составило keff = 7,4·10-4 %.

Приведем значения асимптотической реактивности as, рассчитанной по формуле (3.9) с помощью МПГ на сетке с параметрами h = 1,25 см, H = 2,50 см. Эффективный коэффициент размножения в возмущенной среде составил keff = 0,867136. Таким образом, as = 4,0E-3 или as = 0,62.

Рассмотрим поведение мощности во времени. В таблице 3.25 можно видеть, что в течение переходного процесса расчет SUHAM-SHM хорошо согласуется с расчетом SUHAM-FDM. Относительное отклонение постепенно возрастает в течение процесса возмущения и начинает падать по окончании процесса возмущения. Максимальное относительное отклонение мощности составляет -4,5·10-2 % и это отклонение соответствует временной точке окончания ввода реактивности t = 1 с. Отклонение в последней расчетной точке t = 4 с составляет -7,5·10-3 %.

Таблица 3.25 – Мощности в расчете SUHAM-FDM и относительное отклонение мощности SUHAM-SHM от SUHAM-FDM в различных временных точках (h = 1,25 см, H = 2,5 см, t = 10-3 c) P, отн. ед. P, % t, c 0,0E+00 1,000E+03 – 1,0E–01 9,928E+02 -7,30E- 2,0E–01 9,873E+02 -1,04E- 3,0E–01 9,918E+02 -1,13E- 4,0E–01 1,009E+03 -1,32E- 5,0E–01 1,040E+03 -1,66E- 6,0E–01 1,090E+03 -2,06E- 7,0E–01 1,165E+03 -2,51E- 8,0E–01 1,275E+03 -3,09E- 9,0E–01 1,439E+03 -3,71E- 1,0E+00 1,689E+03 -4,49E- 1,5E+00 2,313E+03 -1,23E- 2,0E+00 2,690E+03 -1,10E- 2,5E+00 3,081E+03 -1,09E- 3,0E+00 3,500E+03 -8,88E- 3,5E+00 3,957E+03 -8,47E- 4,0E+00 4,458E+03 -7,45E- Что касается вычислительных затрат, то отметим, что расчет SUHAM-FDM длился 78 часов 22 минуты, а расчет SUHAM-SHM – 3 часа 50 минут. Таким образом, расчет по МПГ проведен более чем в 20 раз быстрее, чем прямым конечно-разностным методом при незначительных отклонениях.

3.2.5 Транспортный тест TWIGL Помимо диффузионного варианта теста TWIGL, описанного в параграфе 3.2.3, существует транспортный вариант, описанный в работе [85]. Транспортный вариант отличается от диффузионного только набором констант, сохраняя геометрию, граничные условия, параметры запаздывающих нейтронов и законы возмущения. Нейтронно-физические характеристики для транспортного варианта представлены в таблице 3.26. Переход от транспортного варианта теста к ранее описанному диффузионному можно осуществить по формуле Dg =, (3.10) 3( t, g s, g ) P g = 0,01 см-1 [85]. Расчет данного варианта позволяет считая, что sg,P g верифицировать транспортный вариант SUHAM-TD.

Таблица 3.26 – Нейтронно-физические константы для транспортного варианта теста TWIGL f 1, f 2, a1, a2, 11, 12, s 1, 22, t1, t2, Зона s s s см-1 см-1 см-1 см-1 см-1 см-1 см-1 см-1 см-1 см- 1, 2 0,01 0,15 0,007 0,2 0,2281 0,01 0,0 0,8333 0,2481 0, 3 0,008 0,05 0,003 0,06 0,2464 0,01 0,0 0,6667 0,2644 0, На рисунках 3.16 и 3.17 представлены графики изменение мощности для транспортного варианта TWIGL для различных программ. На рисунках 3.16 и 3. имеются следующие обозначения: расчет SUHAM-TD – SUHAM-SHM выполненный для временного шага t = 10-5 с и размером ячейки H = 4 см, DORT-TD – расчет выполненный по программе DORT-TD методом дискретных ординат с шагом t = 10-3 с, MAF/MOC – расчет проведенный методом характеристик с шагом t = 10-4 с, MGSNM – расчет выполненный по программе MGSNM методом дискретных ординат c шагом t = 10-2 с. Важно отметить, что данные для DORT-TD и MAF/MOC были получены по запросу от авторов работ [3] и [114]. Данные для MGSNM взяты из [89]. Из полученных результатов можно отметить хорошее согласование со сторонними кодами.

Рисунок 3.16 – Изменение мощности с течением времени для транспортного TWIGL-R Рисунок 3.17 – Изменение мощности с течение времени для транспортного TWIGL-S Также необходимо отметить, что результаты расчета транспортного теста отличаются от результатов диффузионного теста. В таблице 3.27 приведена относительная мощность в разных временных точках, полученная по SUHAM-TD методом поверхностных гармоник, для диффузионного и транспортного вариантов TWIGL-R при параметрах t = 10-5 с, H = 4 см.

Таблица 3.27 – Относительная мощность по SUHAM-TD для диффузионного и транспортного TWIGL-R, отн. ед.

Диффузионный Транспортный t, c 0,0 1,0000 1, 0,1 1,3083 1, 0,2 1,9592 1, 0,3 2,0746 2, 0,4 2,0919 2, 0,5 2,1094 2, Заключение к главе В данной главе приведено описание разработанного программного комплекса В основу данного кода легли следующие SUHAM-TD.

основополагающие принципы:

1) решаемая задача – уравнения пространственной кинетики;

2) возможность проведение расчетов без диффузионного приближения;

3) возможность проведения расчетов без пространственной гомогенизации;

4) моделируемая среда – полномасштабная активная зона;

5) высокая производительность.

Программный комплекс проводит решение уравнений SUHAM-TD пространственной кинетики посредством метода поверхностных гармоник и метода конечных разностей по пространству и времени в 1D и 2D геометриях. На текущий момент в программном комплексе реализованы нестационарные конечно-разностные уравнения МПГ для трех пробных матриц в квадратной решетке блоков.

позволяет проводить расчеты, как с использованием SUHAM-TD диффузионного приближения, так и без него. Здесь важной чертой является то, что уравнения МПГ инвариантны по отношению к исходным уравнениям, из которых они получены, а именно, к уравнению переноса нейтронов и диффузионному уравнению и факт использования диффузионного приближения определяется только использованием его в расчете пробных матриц. Для получения диффузионных пробных матриц SUHAM-TD использует собственные процедуры, а для расчета пробных матриц без диффузионного приближения – программы РАЦИЯ (G3-приближение) и DICPN (P2-приближение).

Проведена верификация конечно-разностных нестационарных уравнений МПГ и программного комплекса на ряде классических SUHAM-TD пространственно-временных бенчмарках: BSS-6, PHWR, TWIGL (диффузионный и транспортный варианты) и на модифицированном бенчмарке 8-А1. Таким образом, рассматривались задачи в 1D и 2D геометриях с различными величинами вводимой реактивности по различным законам (мгновенно, линейно) и представлениями запаздывающих нейтронов (1 группа и 6 групп). Верификация МПГ-модулей SUHAM-TD проводилась посредствам сопоставления МПГ данных с данными, полученными средствами МКР-модулей, созданных также в рамках SUHAM-TD. Причиной создания собственных модулей, реализующих МКР, является обеспечение «чистоты эксперимента» при сравнении вычислительных затрат. Можно отметить, что результаты, полученные по МПГ, хорошо согласуются с эталонными данными, полученными прямым конечно-разностным моделированием.

Дополнительно проводилось сопоставление с результатами, полученными по сторонним программам. Тут следует отметить, что в большинстве публикуемых работ приводится только графическое представление результатов расчетов, что не позволило провести сопоставление результатов для всех рассмотренных тестов. Однако, следует заметить, что для диффузионного варианта задачи TWIGL представлены данные пяти сторонних кодов, а для транспортного варианта TWIGL – для трех сторонних кодов. Эти коды реализуют различные методы и разработаны различными коллективами. Данное сравнение позволило дополнительно убедиться в качестве результатов SUHAM-TD.

Вместе с хорошим согласованием с другими методами МПГ продемонстрировал значительно меньшие вычислительные затраты по сравнению с классическим МКР. Здесь важным фактом является соблюдение «чистоты эксперимента» в плане сравнения вычислительных затрат. Для этого данное сравнение проводилось только между МПГ-модулями и МКР-модулями комплекса SUHAM-TD. Здесь важным является то, что расчеты по МПГ и МКР проводились с одинаковыми критериями точности и способами организации итераций (например, во всех программах матрицы обращаются с помощью метода верхней релаксации, везде используется полностью неявная схема). Также в программах не применяются какие-либо методы ускорения.

Вместе с тем, показано уменьшение временного преимущества МПГ с уменьшение временного шага. Данный факт объясняется следующей причиной. С уменьшением временного шага обращение матрицы при решении реакторного уравнения производится быстрее как в случае прямого конечно-разностного моделирования, так и в случае МПГ по причине меньшего возмущения элементов в матрице. Однако в случае МПГ часть времени уходит на расчёт пробных матриц, продолжительность которого, как показывает расчёт этих тестов, слабо зависит от размера временного шага. Таким образом, с уменьшением временного шага, расчет пробных матриц занимает все большую долю времени.

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

Глава 4 Разработка и расчет пространственно-временного бенчмарка C5G7-TD для тестирования кинетических нейтронно-физических кодов 4.1 Обзор пространственно-временных бенчмарков Как уже отмечалось, современный нестационарный нейтронно-физический код должен базироваться на возможности решения уравнений пространственной кинетики без диффузионного приближения и пространственной гомогенизации в полномасштабной модели активной зоны. Важным вопросом в разработке такого кода является его верификация. Общей практикой для такой верификации является кросс-верификация, которая заключается в расчете определенных задач (бенчмарков) различными кодами и методами. Однако проведение такой процедуры осложняется следующей проблемой.

На сегодняшний день набор представленных в литературе пространственно временных бенчмарков с заданными групповыми нейтронно-физическими константами содержит своего рода пробел. С одной стороны, существует набор простых диффузионных бенчмарков. Расчетная область в таких задачах представляет собой несколько гомогенных зон, описанных, как правило, малогрупповыми (зачастую только две группы) диффузионными нейтронно физическими макроконстантами. К таким бенчмаркам можно отнести рассмотренные в предыдущей главе тесты: BSS-6, 8-A1, TWIGL, PHWR, а также, например, тест OBLONG [91] и др. Главная их особенность в том, что они просты, и при этом они сохраняют актуальность и в сегодняшние дни, находя свое применение, как правило, для отладки программ на начальном этапе разработки.

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

С другой стороны, имеется набор значительно более «сложных»

бенчмарков, которые описывают гетерогенную структуру среды, в качестве характеристик материалов содержат концентрации нуклидов, и, как правило, включают в себя характеристики обратных связей и т.д. К этой группе бенчмарков можно отнести PWR MOX/UO2 core transient benchmark [21], PBMR coupled neutronics/thermal-hydraulics transient benchmark the PBMR-400 core design [22], Prismatic coupled neutronics/thermal fluids transient benchmark of MHTGR- MW core design [115] и др. Результаты их расчетов содержат дополнительные неопределенности, связанные с неопределенностями ядерных данных, процедурой подготовки групповых сечений и погрешностью сторонних кодов (например, тепло-гидравлических кодов), применяемых для учета обратных связей и др. Этот факт не дает возможности выделить методическую погрешность метода, заложенного в расчетную программу.

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

существуют задачи обладающие признаками из обеих групп. В качестве примера можно привести тесты, которые можно, в целом, отнести к первой группе бенчмарков, но константы, которые содержатся в этих задачах, не являются диффузионными. К таким тестам можно отнести одномерный транспортный тест 16-A1 [20] и ранее рассмотренный транспортный вариант теста TWIGL.



Pages:     | 1 || 3 | 4 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.