Кинетика переноса нейтронов в групповом методе Монте-Карло тема диссертации и автореферата по ВАК РФ 00.00.00, кандидат наук Дьячков Иван Игоревич

• Дьячков Иван Игоревич
• кандидат наук
• 2023

Оглавление диссертации кандидат наук Дьячков Иван Игоревич

Введение

Глава 1. Обзор литературы

1.1 История развития и современное состояние

1.2 Используемое программное обеспечение

Вывод к главе

Глава 2. Групповой физический модуль

2.1 Описание модуля MAGMA

2.2 Алгоритмы взаимодействия нейтронов с веществом по групповым сечениям и макропараметрам

2.4 Верификация группового физического модуля

2.4.1 Ра]счет бесконечной решетки ТВС

2.4.2 Численный бенчмарк-тест C5G7

2.4.2.1 Описание сборки C5G7

2.4.2.2 Описание тестов C5G7

2.4.2.3 Сравнение результатов с эталонными данными

2.4.2.4 Использование результатов программы КИР в качестве реперных

Вывод к главе

Глава 3.Алгоритм моделирование кинетических процессов прямым методом Монте-Карло

3.1 Описание алгоритма моделирования кинетики прямым методом Монте-Карло

3.2 Верификация алгоритмов моделирования переходных процессов программой КИР на основании численного бенчмарк-теста RPIGS-KIN

3.2.1 Описание теста RPIGS-KIN

3.2.2 Результаты расчёта

3.3 Верификация алгоритмов моделирования переходных процессов программой

КИР на основании численного бенчмарк-теста С507-ТЭ

3.3.1. Краткое описание тестов С507-ТО

3.3.2. Вариант 0 (ТОО)

3.3.3. Вариант 3 (ТО3)

Вывод к главе

Глава 4. Алгоритмы моделирования кинетического источника нейтронов методом Монте-Карло

4.1 Описание метода

4.2. Программная реализация кинетического источника нейтронов

4.3. Верификация алгоритма кинетического источника нейтронов

4.3.1. Предварительный вариант

4.3.2. Результаты двумерных тестов С507-ТЭ

4.3.3. Вариант ТОО

4.3.4. Вариант ТО1

4.3.5. Вариант ТО2

4.3.6. Вариант ТО3

Вывод к главе

Заключение

Список сокращенийи условных обозначений

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

ПРИЛОЖЕНИЕ А

ПРИЛОЖЕНИЕ Б

ПРИЛОЖЕНИЕ В

ПРИЛОЖЕНИЕ Г

ВВЕДЕНИЕ

С момента 1932 года, когда Джеймс Чедвик открыл такую частицу как нейтрон, до настоящего времени, физика взаимодействия частиц стремительно развивается. В настоящее время постоянно разрабатываются, производятся, вводятся в эксплуатацию и эксплуатируются новые ядерные энергетические установки как исследовательские, так и промышленные АЭС. В таких странах как Россия, США, Швеция, Франция и в ряде других стран атомная энергетика вносит значительный вклад в объемы производства электроэнергии и является неотъемлемой частью энергетики данных стран.

Но с развитием таких сложных и потенциально опасных систем как ядерный реактор возрастают требования по защите и обеспечению безопасности, которые необходимо обосновывать для всех этапов жизненного цикла ЯЭУ. В связи с развитием вычислительных технологий стало возможным проводить высокоточные расчеты динамики ЯЭУ [1-2] при пуске, выходе на мощность, работе на мощности и проектных и запроектных авариях в целях обеспечения безопасности установки на ранних этапах ее разработки и сопровождать ее до вывода из эксплуатации.

Основной частью динамики ЯР является кинетика [3] - нейтронно-физическая составляющая расчета. В данный момент существует множество методов решения нестационарного уравнения переноса нейтронов. Их можно разделить на две основные группы: детерминистические (включают в себя ряд приближений: диффузионные, точечной и распределенной кинетики, метод дискретных ординат и поверхностных гармоник и др., и высокоточный метод-статистический (прямой метод Монте-Карло)).

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

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

В настоящее время работы по созданию и усовершенствованию программ, основанных на решении нестационарного уравнения переноса нейтронов методом Монте-Карло ведутся по всему миру (Россия - КИР [4] (НИЦ КИ) и TDMCC (ВНИИЭФ-) [5], Иран - TDCM [6], Финляндия - SERPENT 2 [7], Нидерланды -TRIPOLI [8-9], Китай - TMCC [10]).

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

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

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

Цель диссертационной работы

Разработка алгоритмов и программных средств решения пространственно-

временных задач переноса нейтронов методом Монте-Карло в групповом

приближении для повышения точности и надежности расчета нейтронно-

физических характеристик ядерных реакторов.

Для достижения данной цели были решены следующие задачи:

• Разработка и программная реализация алгоритмов группового физического модуля с учетом запаздывающих нейтронов и изотропного и анизотропного рассеяния для решения стационарного и нестационарного уравнения переноса нейтронов методом Монте-Карло в групповом приближении.

• Разработка методики и алгоритмов формирования кинетического источника нейтронов, и их программная реализация.

• Верификация алгоритмов на ряде тестовых задач.

Научная новизна результатов работы

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

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

Практическая значимость

• Показана практическая применимость данных алгоритмов и методов при моделировании переходных процессов ЯЭУ, которые позволяют сократить время счета и ресурсозатраты, и направить их на увеличение точности счета.

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

Основные положения, выносимые на защиту

• Алгоритмы группового физического модуля.

• Методика и алгоритмы формирования кинетического источника нейтронов.

• Программная реализация алгоритмов в модулях MAGMA и KS.

• Результаты верификации.

Личный вклад автора

• Разработка алгоритмов и их реализация в виде группового физического модуля MAGMA.

• Разработка методики формирования кинетического источника и реализация в виде модуля KS.

• Создание моделей и проведение этапов верификации разработанных алгоритмов.

Введение диссертации (часть автореферата) на тему «Кинетика переноса нейтронов в групповом методе Монте-Карло»

Апробация работы

Основные положения диссертации доложены на следующих российских и

международных конференциях:

• «Международная конференция молодых специалистов, ученых и аспирантов по физике ядерных реакторов "Волга-2018"» (3-7 сентября 2018г.), спортивно-оздоровительный лагерь «Волга», Тверская область, Россия (1 доклад)

• «Международные научно-технические конференции "обеспечение безопасности АЭС с ВВЭР"» (21 - 24 мая 2019г.) пр. Берёзовой Рощи 12, Москва, Россия (1 доклад)

• Всероссийская научно-техническая конференция «Нейтронно-физические проблемы атомной энергетики» («Нейтронника-2019»), (27 - 29 ноября 2019г.), пл. Бандаренко, г. Обнинск, Калужская обл., Россия (1 доклад)

• Всероссийская научно-техническая конференция «Нейтронно-физические проблемы атомной энергетики» («Нейтронника-2022»), (31 мая - 5 июня 2022г.), пл. Бандаренко, г. Обнинск, Калужская обл., Россия (1 доклад)

• «Международная конференция молодых специалистов, ученых и аспирантов по физике ядерных реакторов "Волга-2022"» (5-9 сентября 2022 г.), спортивно-оздоровительный лагерь «Волга», Тверская область, Россия (1 доклад)

Достоверность результатов работы

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

Публикации

По теме диссертации опубликовано шесть статей, пять из которых в рецензируемых научных изданиях ВАК:

• I.I. Dyachkov, M.V. Ioannisian. Calculation of the numerical benchmark NURISP based on the complex MRNK+KEDRD // Journal of Physics Conference Series, 1133(1):012031. 2018

• М.В. Иоаннисиан, И.И. Дьячков, В.П. Быков, С.Ю. Закиров. Моделирование динамических процессов с использованием метода многозонной кинетики на примере численных бенчмарк-тестов // Атомная энергия, Том 126 (№2), 2019

• И.И. Дьячков, М.В. Иоаннисиан. Верификация группового физического модуля в программе расчета нейтронной кинетики КИР2 на основе бенчмарк-теста C5G7 // ВАНТ: Серия: Ядерно-реакторные константы, № 3. 2020

• М. В. Иоаннисиан, В. Д. Давиденко и И. И. Дьячков. О моделировании начального источника мгновенных нейтронов для решения нестационарных задач методом Монте-Карло // ВАНТ: Серия: Ядерно-реакторные константы, №3. 2020.

• В. В. Белоусов, М. И. Гуревич, В. Д. Давиденко, И. И. Дьячков, М. В. Иоаннисиан, А. А. Ковалишин, М. Р. Малков, К. Ф. Расчкач , К. Г. Чернов и Р. В. Широков. Программный комплекс КИР2 для моделирования стационарного и нестационарного переноса частиц методом Монте-Карло. // ВАНТ: Серия: Физика ядерных реакторов, № 1. 2022.

• В.И. Белоусов, В.Д. Давиденко, И.И. Дьячков, Н.Н. Елкин, В.Г. Зборовский, М.В. Иоаннисиан, В.В. Лиханский, М.Р. Малков, О.В. Хоружий и М.Г.

Чернецкий. Связанный нейтронно-физический и теплогидравлический расчёт тестовой модели реактора с водой сверхкритического давления с применением метода Монте-Карло // ВАНТ: Ядерно-реакторные константы, № 1, рр. 24-40. 2023.

Структура и объем диссертации

Диссертация состоит из введения, обзора литературы, 4 глав и заключения. Работа содержит 126 страниц печатного текста, 45 рисунков и 51 таблиц. Библиография насчитывает 58 наименования.

ГЛАВА 1. ОБЗОР ЛИТЕРАТУРЫ 1.1 История развития и современное состояние

До недавних пор метод Монте-Карло использовался только для решения стационарного уравнения переноса нейтронов, а решение нестационарного уравнения полностью основывалось на детерминистических методах [11]. Использование метода Монте-Карло было ограничено в связи с его требовательностью к ресурсам и времени счета, а именно, с недостаточными вычислительными мощностями того времени. А детерминистические коды, которые обычно быстры при расчете стационарных состояний, имели огромное преимущество.

При решении нестационарного уравнения переноса нейтронов программами, основанными на детерминистических методах, в связи с большим числом неизвестных (1014) обычно используется двухэтапный подход [12-13] . Моделируется два состояния для рассматриваемой модели в различные моменты времени, после чего происходит вычисление изменения потока нейтронов. Также, использование различных методов решения уравнения переноса детерминистическими методами и используемые в них приближения, которые имеют определенные ограничения, сильно зависят от рассматриваемой задачи и геометрической модели, для которых решается данная задача. Именно по данной причине требуются программные средства, которые будут служить эталонными, сводя ограничения данных методов и приближений к минимуму и позволив использовать данный код, как «золотой стандарт», для кросс-верификации программ, основанных на детерминистических методах.

Метод Монте-Карло хорошо зарекомендовал себя как код прецизионного (высокоточного) класса при решении стационарного уравнения переноса нейтронов. И в связи с вычислительными возможностями настоящего времени работы по созданию и усовершенствованию программ, в которых решается нестационарное уравнение переноса нейтронов по методу Монте-Карло, ведутся по всему миру (Россия -TDMCC [5], Иран - TDCM [6], Финляндия - SERPENT 2 [7], Нидерланды - TRIPOLI [8-9], Китай - TMCC [10] и т.д.).

Нередко в программах, основанных на методе Монте-Карло, также вводятся различные приближения, но, в отличие от детерминистических методов, данные приближения применяются для ускорения расчета и уменьшения затрат ресурсов. Так, например, в программах TDMC [5], TRIPOLY[8-9], TMCC [8] и TDMCC [6] и SERPENT 2 [7] используются весовые методы, которые моделируют взаимодействия нейтрона с веществом, присваивая ему вес. В программе TDKENO (T-REX - усовершенствованная версия [14]) используется квазистатическое приближение. Также в большинстве программ существует контроль популяции нейтронов, т.е. временная область делится на интервалы, при пересечении которых число нейтронов нормируется на необходимое число, что позволяет экономить вычислительные ресурсы. В некоторых программах используют гибридные методы, в которых происходит взаимодействие метода Монте-Карло и детерминистических программ. По первому методу происходит расчет распределения нейтронов и получение функционалов в стационарном состоянии, по второму - расчет эволюции нейтронов во времени. При создании программ, основанных на методе Монте-Карло с использованием каких-либо приближений, снова появляется потребность в эталоне, на который требуется ориентироваться при их создании и кросс-верификации.

Но не смотря на введение определенных приближений и на современную вычислительную мощность, на сегодняшний день, выполнение расчетов кинетических процессов методом Монте-Карло является ограниченным в применении. Расчеты, проводимые по данному методу, ограничены либо небольшими моделями ЯЭУ, либо бесконечной решеткой ТВС [15-16]. В данный момент не существует программы, которая смогла бы провести высокоточный расчет с достаточной точностью кинетического процесса для реактора типов ВВЭР, РБМК, PWR и др. в разумные временные сроки.

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

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

Данное соотношение возможно получить несколькими способами. В первом случае приходится проводить предварительный стационарный расчет с получением распределения нейтронов по модели с корректными соотношениями запаздывающих нейтронов, как сделано в программе SERPENT 2 [17], после чего полученные результаты используются при решении нестационарного уравнения переноса нейтронов. Как правило, такое моделирование источника носит двумерный характер и ограничивается зонным представлением, а не точечным.

Другой способ, это моделирование процесса выхода в критическое состояние от внешнего источника нейтронов. Данный метод в настоящее время реализован в программе КИР [4]. Этот способ является наиболее временно и ресурсозатратным, т.к. при моделировании выхода в критическое состояние требуется моделирования нестационарного процесса длительностью около 100 секунд, чтобы получить корректное соотношение мгновенных и запаздывающих нейтронов и равновесное состояние ядер-предшественников запаздывающих нейтронов. Так, в программе RMS [18] был реализован данный метод с использованием метода ускорения сходимости источника деления за счет экстраполяции источника [19]. Существуют и другие детерминистические методы ускорения сходимости начального состояния от источника нейтронов [20-22].

Например, метод предложенный Sjenitzer и Hoogenboom и реализованный в Tripoli [23]. В данной программе реализован алгоритм неактивных итераций мощности, при которых происходит накопление запаздывающих нейтронов и изъятие их в конце каждой итерации. Tripoli усредняет каждую выборку. После

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

В данной диссертационной работе предложена методика формирования кинетического источника нейтронов, которая позволяет:

• получить в начальный момент кинетического процесса критическое состояние реактора

• избежать длительных предварительных расчетов

• уйти от зонного представления распределения нейтронов, приблизив расчет «максимально» к реальному.

К числу программных средств (ПС), основанных на Монте-Карло, относится программа КИР [24-25], которая предназначена для получения решения на одноядерных и многоядерных компьютерах не только стационарного, но и нестационарного уравнения переноса нейтронов с учётом запаздывающих нейтронов [26]. В программе реализован аналоговый метод Монте-Карло [27], который моделирует историю нейтрона на основе файлов оценённых ядерных данных в системах с трёхмерной геометрией в однородной и неоднородной среде.

В программе КИР при моделировании кинетики ядерных реакторов решается стандартное уравнение с интегрированием по лучу, которое является достаточно известным и представлено, например, в работах [28-29]. Приведём это уравнение из работы [29], в котором расписаны члены, обеспечивающие генерацию мгновенных и запаздывающих нейтронов:

1.2 Используемое программное обеспечение

Ф (г, П, Е, X)

(1)

где

q(r ', Q, E, t ') = J dE ' J dQ '

, t ') =

E' Q '

(EW Д^ (r ',E',t')

Ф(г ', Q ', E ', t') +

+J dE ' J dQ ' 2 Zx (r ', Q ', E ', t ') f (r ', Q ', E ' ^ Q, Е)Ф(г ', Q ', E ', t ') + , (2)

E ' Q ' fission

Nd y ( E )

+2. ^T^Wj (r ', t ') + Q(r ' Q, E, t ')

j=1

gC J<f' *) Cj (r, t) = J dE J dQPj (r, £>(r, £')S fissiœi (r, E', t)0(r, Q', E, t) (3)

^ E ' Q'

Здесь все обозначения общеприняты.

(r,^,E,t) - фазовые координаты частицы (r - радиус-вектор местоположения, О - вектор направления, E - энергия, t - время); v - скорость нейтрона; r ' = r - s Q ' ; s'

t ' = t---

v ;

Xprompt(E) - нормированный спектр мгновенных нейтронов; v(r, E) - ожидаемое полное число нейтронов, испускаемых на одно деление в точке r, вызванное нейтроном с энергией E;

G(r, t) - плотность предшественников запаздывающих нейтронов в /-ой группе;

Nd - число групп предшественников запаздывающих нейтронов, обычно рассматривается 6 или 8 групп;

Ai - постоянная распада в группе /;

Xi(E) - энергетический спектр запаздывающих нейтронов группы i; Д(г, E) - доля величины v(r, E), отнесённая к i-й группе предшественников запаздывающих нейтронов;

ei(r,E)^v(r,E) - ожидаемое число предшественников i-ой группы, образующихся в результате деления ядра в точке r нейтроном с энергией E;

Д( r,E) = 2 Д. (r, e) , где i = l,...Nd;

i

[1 - e(r, E)]v(r, E) - число мгновенных нейтронов, появляющихся в результате деления ядра в точке r нейтроном с энергией E;

N - число групп запаздывающих нейтронов, при решении задач реакторной кинетики обычно рассматривается 6 или 8 групп запаздывающих нейтронов;

в и в] относятся к энергии Е', а Хр и /у - к энергии Е:

£0«/(г, Е, ?) - полное макроскопическое сечение взаимодействия нейтронов;

Хх(г;П',Е'^П,Е;?) - макроскопическое сечение взаимодействия типа х нейтронов с веществом;

,Е' ^ &,Е) - вероятность того, что при столкновении в точке г нейтрона, имеющего направление полёта Л' и энергию Е', происходит реакция типа х и в той же точке появится нейтрон с направлением движения Л и энергией Е (вероятность перехода).

При решении данных уравнений в программе КИР находится полное число взаимодействий в заданной области фазового пространства - интервале по времени |/0, ^ + аг] , энергетическом интервале [Ео, Ео + АЕ] и выделенном объёме

^ е А

которое описывается интегралом по потоку:

Е0+Ае ^+М

; ее[еое+Ае] = |Лг | ЛП | ЛЕ \ Лт ап(г,Е,т)-Ф(г,П,Е,т), (4)

У2 |п|=1 е0 г0 '

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

тУег; те[/0 ,/0 +А/]; Ее[Е0 Е +АЕ] гй = -2--(5)

Аt (5)

Для получения усреднённого интегрального потока нужно положить <г„ (г, Е,т) = 1 в выражении (4).

Также программа КИР используется для расчета динамики в программном комплексе ДАРЕУС [30-32] (связка КИР и ГАРД), предназначенного для моделирования процессов в экспериментальных растворных реакторах. Кроме того, проводились и другие связанные расчеты - связка нейтронно-физической программы КИР и теплогидравлической ТК-СКД [33].

ПС КИР написано на языках программирования FORTRAN стандарта 2008 и С++ стандарта 2011 года (физический модуль PHM). В программе используются и усовершенствованы некоторые особенности известного российского кода MCU-5, и в этом смысле КИР может считаться одним из продолжений развития данного кода, хотя КИР в своём развитии уже очень отдалился от MCU-5: общим между этими кодами на данный момент является лишь базовый геометрический модуль и базовые принципы организации расчёта [4].

Отличительная особенность кода КИР связана с применением методов объектно-ориентированного программирования (ООП) для разработки кода. Наличие классов и объектов позволяет упростить понимание кода и его составляющих, сделать код гибким и более открытым для совершенствования, что облегчает расширение и упрощает поддержку кода. Кроме того, использование ООП упрощает переход с одного языка программирования на другой и позволяет воспользоваться объектами из другого языка через программный интерфейс (Application Program Interface, API). Также надо отметить, что, с другой стороны, использование ООП может ускорить разработку и поддержку программы.

КИР позволяет выполнять расчёты на многопроцессорных компьютерах. Распараллеливание программ выполнено на базе программного интерфейса MPI (Message Passing Interface). Он является наиболее распространённым стандартом интерфейса обмена данными в параллельном программировании, а его реализации существуют для большого числа компьютерных платформ.

С учётом возможности решения, наряду со стандартными статическими задачами, полноценной нестационарной задачи комплекс существенно расширяет область применения метода Монте-Карло в проектной практике. Ниже перечислены основные прикладные задачи комплекса [4]:

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

• расчёт эффектов и коэффициентов реактивности;

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

• создание математических бенчмарков для верификации программ инженерного класса.

Вывод к главе 1

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

Одной из основных проблем программ, основанных на методе Монте-Карло, является получение начального состояния системы (условной критики). Существует различные методики получения условной критичности, но, как правило, они являются трудоёмкими.

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

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

ГЛАВА 2. ГРУППОВОЙ ФИЗИЧЕСКИЙ МОДУЛЬ 2.1 Описание модуля MAGMA

Для программы КИР были разработаны алгоритмы взаимодействия нейтрона с веществом по групповым сечениям и макропараметрам с учетом запаздывающих нейтронов для расчета кинетики ЯЭУ [34]. Алгоритмы реализованы в виде физического модуля MAGMA (Many-Group Module with anisotropic scattering -МногоГрупповой Модуль с учетом Анизотропного рассеяния).

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

Ввод данных возможен как из общего входного файла КИР, так и отдельными файлами для каждого материала. Входными параметрами модуля являются многогрупповые константы материалов, которые готовятся как на основе модулей подготовки констант комплекса UNK, так и непосредственно с использованием физического модуля PHM программы КИР или другой сторонней программой. Предусмотрен ввод всех основных макроскопических сечений (полное, поглощения, деления, захвата, упругое и неупругое рассеяния, и реакций (n, 2n) и (n, 3n)), а также спектры запаздывающих и мгновенных нейтронов, матрицы рассеяния и косинусов рассеяния. Так как модуль предназначен и для расчета кинетических процессов, то присутствует ввод и долей запаздывающих нейтронов,

и времен запаздывания. Также во входном файле присутствует возможность использовать функции зависимости сечений от времени.

2.2 Алгоритмы взаимодействия нейтронов с веществом по групповым

сечениям и макропараметрам

На рисунке 1 представлена схема алгоритма группового модуля. В процессе счета из банка извлекается произвольный нейтрон и определяется в каком материале произошло его взаимодействие со средой по средней длине пробега в транспортном модуле программы КИР. Далее групповой модуль определяет тип реакции взаимодействия - деление, рассеяние (упругое или неупругое) или захват (также присутствуют реакции (п, 2п) и (п, 3п)).

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

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

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

Рисунок 1 - Схема алгоритма розыгрыша реакции нейтрона в групповом модуле

MAGMA

Далее возможно два варианта моделирования направления полета нейтрона. Если в исходных данных группового модуля MAGMA отсутствует ввод средних косинусов рассеяния, которые определяют процесс моделирования анизотропного рассеяния, то направление нейтрона (u, v, w) после взаимодействия считается изотропным и разыгрывается по стандартному алгоритму:

Разыгрывается пара случайных чисел ^ и и определяются коэффициенты:

A = 2£ -1 (6)

в = £ (7)

C = A2 + B2 (8)

Если C > 1, то разыгрывается новая пара ^ , ^2 до тех пор, пока не будет C<1, и далее:

cos^= (A2 - B2)/C (9)

sin^ = 2AB/C (10)

cosé= 2^3 - 1 (11)

D = (cosé)2 (12)

sin#= Vi - D (13)

u= cosé (14)

v= sinécos^ (15)

w= sin é sin ^ (16)

- случайное число, равномерно распределённое на [0,1]). Если в исходных данных (ПРИЛОЖЕНИЕ А) MAGMA присутствует ввод карты COSMATR с направляющими косинусами, то направление нейтрона (u, v, w) после взаимодействия разыгрывается по следующему алгоритму.

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

2.4 Верификация группового физического модуля

2.4.1 Расчет бесконечной решетки ТВС На первом этапе верификации алгоритмов модуля MAGMA была выбрана бесконечная решетка тепловыделяющих сборок. Геометрическая конфигурация кассеты была задана на основе одного из вариантов водо-водяного реактора со сверхкритическим давлением [33] (рисунок 22). ТВС представляет собой правильную шестигранную призму высотой 390 см с размером «под ключ» равным 14.46 см. Сборка ТВС состоит из 199 тепловыделяющих элементов,

расположенных в гексагональной решётке с шагом 9.42 мм, расстояние между которыми 1.02 мм. Сборка окружена чехлом толщиной 2 мм.

Рисунок 2 - Поперечное сечение аксиального слоя ТВС в области твэлов.

В данном тесте было проведено сравнение результатов программы КИР с использованием PHM модуля (реперные результаты), который использует при решении уравнения переноса нейтронов оцененные ядерные данные, с результатами программы КИР с модулем MAGMA. В данном тесте рассматривалась зависимость отклонений результатов программы КИР с модулем MAGMA от числа групп, на которые делится энергетическая область. Также было проведено исследование влияния изотропии и анизотропии рассеяния на отклонение результатов модуля MAGMA при изменениях числа энергетических групп.

Для расчётов по групповому модулю MAGMA был выбран ряд разбиений энергетической области на 2, 4, 6, 12, 25, 32, 45, 59 и 299 групп (границы энергетических групп представлены в Приложении Б) для оценки зависимости эффективного коэффициента размножения от числа энергетических групп и

характера разбиения. Разбиения сформированы следующим образом. За основу взяты энергетические границы 25 групп БНАБ [35]. Для формирования разбиений с большим числом групп (32, 45, 59) уточнялись границы в области резонансов и энергий порогового деления чётных изотопов Pu и U-238. Разбиение на 299 групп выбрано для оценки точности расчёта по модулю MAGMA.

Для каждого энергетического разбиения был сделан расчёт на основе физического модуля PHM (прямой расчёт), в ходе которого по модулю регистрации MAGMAFunctionals были получены макроконстанты для физического модуля MAGMA. После этого по модулю MAGMA были выполнены две группы расчётов:

• расчёт с анизотропным рассеянием, в котором использовались матрицы средних косинусов;

• расчёт с изотропным рассеянием.

В результате расчётов получена оценка влияния анизотропного и изотропного рассеяния при расчёте на групповых константах и при прямом расчёте на эффективный коэффициент размножения и спектр плотности потока нейтронов. Все результаты получены со статистикой 300 млн. историй нейтронов.

На рисунке 3 представлены результаты расчёта спектра плотности потока нейтронов, полученные при расчётах с применением модулей PHM и MAGMA с 299 - групповым разбиением для плотности потока в топливе.

При сравнении результатов проведённых расчётов наблюдаются расхождения плотности потока в топливе с использованием анизотропного и изотропного рассеяния с 10 по 70 границ энергетических групп. Данным границам соответствуют тепловая область энергий с 0эВ до 3,83эВ. Максимальное отклонение в плотности потока для расчёта с изотропией наблюдается в 9 группе (9,66%.).

Список литературы диссертационного исследования кандидат наук Дьячков Иван Игоревич, 2023 год

// // //

О.----

// // ж

о" о-__

■■■■■■■■

¿г

СХ

С

4х

КИР* -о-КеШСХР

Рисунок 9 - Зависимость значений ЯМБ от программы

0.35 0.30 0.25

0)

ас' 0.20

н

о

! 0.15

к

3

° 0.10 0.05 0.00

SSS *// */

_,_КИР* -о- Ref.MCNP Рисунок 10 - Зависимость значений MRE от программ

Таким образом, результаты расчета по программе КИР хорошо согласуются с результатами Reference MCNP по двум выбранным методикам сравнения и указывают на корректную работу разработанного модуля MAGMA.

2.4.2.4 Использование результатов программы КИР в качестве реперных

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

В таблице 5 показан пример, как зависимости процентной доли твэлов в конфигурации Rodded A по всей области активной зоны, энерговыделения которых

вошли в доверительные интервалы (1а - 68%) программы КИР от статистики проводимого расчета.

Таблица 5 - Зависимость процентной доли твэлов, энерговыделения которых вошли в доверительные интервалы, для конфигурации Коёёеё А от статистики расчета КИР *

Статистика расчета КИР

Программа 1*108 5*108 2,5*1010

ATTILA 86% 53% 26%

CRX 71% 39% 37%

PARTISN 93% 86% 73%

TORT-PSU 81% 44% 23%

UNKGRO 21% 9% 6%

Исходя из этого, более корректно использовать "взвешенные" отклонения AVG, RMS и MRE, не зависящие от статистики расчета реперной программы.

В таблицах 6, 7, 8 представлены результаты отклонений AVG, RMS, MRE программ детерминистического класса ATTILA, MCCG3D, CRX, PARTISN, UNKGRO и программ, основанных на методе Монте-Карло, VIM, MCNP*, UNKMK, от реперных результатов КИР * и Reference MCNP для трех сборок Unrodded, Rodded A, Rodded B. Данные представлены для трех конфигураций по трем слоям.

Таблица 6 - Ошибки AVG, RMS, MRE для варианта UNRODDED

Программа AVG RMS MRE

КИР* Ref.MCNP КИР* Ref.MCNP КИР* Ref.MCNP

СЛОЙ 1

UNKMK 0,122 0,173 0,156 0,224 0,052 0,065

PARTISN 0,153 0,180 0,184 0,230 0,083 0,080

MCNP* 0,159 0,380 0,205 0,400 0,064 0,140

VIM 0,175 0,224 0,228 0,290 0,065 0,087

ATTILA 0,196 0,210 0,218 0,270 0,072 0,070

CRX 0,216 0,240 0,277 0,320 0,084 0,090

MCCG3D 0,226 0,260 0,292 0,330 0,129 0,120

UNKGRO 0,905 0,960 1,230 1,290 0,367 0,400

Программа AVG RMS MRE

КИР* Ref.MCNP КИР* Ref.MCNP КИР* Ref.MCNP

СЛОЙ 2

UNKMK 0,125 0,192 0,171 0,249 0,036 0,058

MCNP* 0,170 0,430 0,224 0,440 0,051 0,130

VIM 0,195 0,192 0,254 0,249 0,059 0,058

PARTISN 0,188 0,190 0,206 0,240 0,064 0,060

ATTILA 0,211 0,200 0,236 0,260 0,062 0,060

CRX 0,226 0,260 0,282 0,350 0,066 0,070

MCCG3D 0,271 0,280 0,319 0,500 0,089 0,090

UNKGRO 0,985 1,060 1,433 1,530 0,277 0,300

СЛОЙ 3

UNKMK 0,181 0,260 0,242 0,344 0,030 0,044

MCNP* 0,257 0,570 0,332 0,590 0,043 0,100

VIM 0,269 0,335 0,353 0,435 0,043 0,057

ATTILA 0,203 0,400 0,238 0,460 0,050 0,080

CRX 0,205 0,280 0,260 0,370 0,034 0,050

PARTISN 0,333 0,270 0,373 0,360 0,049 0,040

MCCG3D 0,551 0,440 0,635 0,560 0,090 0,070

UNKGRO 1,283 1,400 1,728 1,860 0,189 0,200

Таблица 7 - Ошибки AVG, RMS, MRE для варианта RODDED A

Программа AVG RMS MRE

КИР* Ref.MCNP КИР* Ref.MCNP КИР* Ref.MCNP

СЛОЙ 1

UNKMK 0,114 0,168 0,144 0,218 0,052 0,069

MCNP* 0,147 0,370 0,185 0,380 0,064 0,150

VIM 0,182 0,221 0,234 0,287 0,080 0,093

ATTILA 0,187 0,180 0,208 0,230 0,091 0,080

PARTISN 0,291 0,310 0,323 0,350 0,178 0,170

MCCG3D 0,260 0,290 0,359 0,400 0,123 0,120

CRX 0,349 0,370 0,390 0,420 0,159 0,180

UNKGRO 0,938 0,970 1,259 1,290 0,428 0,450

Программа AVG RMS MRE

КИР* Ref.MCNP КИР* Ref.MCNP КИР* Ref.MCNP

СЛОЙ 2

UNKMK 0,127 0,185 0,167 0,241 0,038 0,059

MCNP* 0,170 0,420 0,223 0,430 0,052 0,130

VIM 0,193 0,240 0,251 0,316 0,061 0,074

PARTISN 0,071 0,180 0,088 0,230 0,021 0,050

ATTILA 0,197 0,200 0,226 0,250 0,054 0,060

MCCG3D 0,204 0,280 0,282 0,390 0,060 0,080

CRX 0,246 0,320 0,321 0,430 0,085 0,110

UNKGRO 1,036 1,130 1,491 1,600 0,294 0,320

СЛОЙ 3

UNKMK 0,199 0,265 0,255 0,338 0,028 0,038

MCNP* 0,231 0,600 0,299 0,620 0,032 0,090

VIM 0,285 0,351 0,368 0,439 0,040 0,051

ATTILA 0,185 0,320 0,236 0,390 0,034 0,050

MCCG3D 0,602 0,600 0,671 0,690 0,103 0,100

PARTISN 0,927 0,870 1,059 1,050 0,165 0,160

CRX 1,172 1,260 1,212 1,320 0,178 0,190

UNKGRO 1,429 1,510 1,850 1,950 0,188 0,200

Таблица 8 - Ошибки AVG, RMS, MRE для варианта RODDED B

Программа AVG RMS MRE

КИР* Ref.MCNP КИР* Ref.MCNP КИР* Ref.MCNP

СЛОЙ 1

UNKMK 0,116 0,145 0,146 0,186 0,065 0,072

MCNP* 0,137 0,340 0,175 0,350 0,073 0,170

VIM 0,167 0,191 0,213 0,247 0,085 0,096

ATTILA 0,075 0,170 0,089 0,210 0,044 0,100

MCCG3D 0,315 0,340 0,399 0,420 0,164 0,190

CRX 0,700 0,740 0,778 0,850 0,468 0,500

PARTISN 0,833 0,820 0,865 0,890 0,417 0,390

UNKGRO 1,135 1,200 1,415 1,490 0,674 0,710

Программа AVG RMS MRE

КИР* Ref.MCNP КИР* Ref.MCNP КИР* Ref.MCNP

СЛОЙ 2

UNKMK 0,134 0,180 0,177 0,230 0,038 0,054

MCNP* 0,157 0,420 0,205 0,430 0,047 0,130

VIM 0,214 0,251 0,278 0,318 0,063 0,076

ATTILA 0,137 0,200 0,161 0,250 0,045 0,060

MCCG3D 0,387 0,430 0,490 0,530 0,135 0,140

CRX 0,732 0,770 0,813 0,890 0,214 0,220

PARTISN 0,628 0,720 0,819 0,900 0,221 0,250

UNKGRO 1,148 1,200 1,650 1,740 0,306 0,320

СЛОЙ 3

UNKMK 0,237 0,287 0,297 0,370 0,028 0,034

MCNP* 0,277 0,660 0,353 0,670 0,033 0,080

VIM 0,325 0,370 0,408 0,464 0,038 0,044

ATTILA 0,285 0,470 0,314 0,570 0,033 0,050

MCCG3D 0,573 0,710 0,691 0,870 0,067 0,080

UNKGRO 1,753 1,910 2,241 2,430 0,189 0,210

CRX 1,998 2,190 2,129 2,350 0,230 0,050

PARTISN 2,152 0,710 2,466 0,870 0,286 0,080

Из таблиц видно, что при использовании результатов КИР* в качестве реперных, значения отклонений AVG, RMS и MRE в большинстве случаев уменьшаются по сравнению с Reference MCNP. Для более наглядного представления на рисунках 11, 12 и 13 некоторые данные из таблиц 6, 7, 8 приведены в графическом виде.

1.4 1.2 1.0

0.8

га Ю

а

g о.б

0.4 0.2

0.0

# ^ /

4 С

КИР* --0--Ref.YCNP

riF

Рисунок 11 - Зависимость значений ЛУО от программы для конфигурации

Коёёеё Л, 2 слой

Рисунок 12 - Зависимость значений RMS от программы для конфигурации

Rodded A, 2 слой

0.35

^ 0.20 -сз

ю

I 0.15 -О

0.05 -

0.30 -

0.25 -

0.10 -

0.00

• КИР* - о - Ы.е£МСХР

Рисунок 13 - Зависимость значений МКЕ от программы для конфигурации

Яоёёеё А, 2 слой

Результаты этого анализа свидетельствуют о целесообразности использования расчетных данных программы КИР со статистикой 2,5-1010 в качестве реперных. Более подробные результаты расчета программы КИР представлены в ПРИЛОЖЕНИИ В.

Верификация программы КИР с разработанным групповым модулем (МАОМА) продемонстрировала корректность работы описанных алгоритмов как при расчете решетки ТВС водо-водяного реактора, так и при расчете международного численного бенчмарк-теста С5О7.

Результаты эффективного коэффициента размножения согласуются с результатами как при сравнении с расчетом по программе КИР с модулем РНМ, так и с другими прецизионными программами, основанных на методе Монте-Карло.

Вывод к главе 2

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

Результаты распределений энерговыделения, полученных по программе КИР, хорошо согласуются с результатами реперной программы Ref.MCNP как в пределах доверительных интервалов, так и по взвешенным отклонениям AVG, RMS и MRE.

При анализе результатов было замечено, что при использовании в качестве реперных расчетные данные программы КИР со статистикой 25 млрд. наблюдается снижение оценок отклонений AVG, RMS, MRE для большинства результатов расчета по детерминистическим кодам.

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

Анализ результатов показал возможность использования результатов программы КИР (25 млрд.) в качестве реперных, а для сопоставления результатов желательнее применять оценки AVG, RMS и RME.

ГЛАВА 3. АЛГОРИТМ МОДЕЛИРОВАНИЕ КИНЕТИЧЕСКИХ ПРОЦЕССОВ ПРЯМЫМ МЕТОДОМ МОНТЕ -КАРЛО 3.1 Описание алгоритма моделирования кинетики прямым методом Монте-

Карло

Для моделирования кинетических процессов протекающих в ЯЭУ в программе КИР были модернизированы алгоритмы для учета временных зависимостей движения нейтрона. Были реализованы алгоритмы линейного изменения свойств материалов во времени и замены материалов во времени, для моделирования изменения положения органов регулирования через физический модуль, изменения плотности материала и мгновенного погружения органов регулирования.

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

При моделировании кинетических процессов, в начальный момент которых рассматриваемая модель находится в условно-критическом состоянии, в программе КИР используется метод, который моделирует выход в критическое состояние от внешнего источника нейтронов. Для данного источника нейтронов были реализованы алгоритмы розыгрыша нейтронов.

Опишем процесс переноса нейтронов в программе КИР от первичного прародителя (внешнего источника) и его потомков (рисунок 14-15) [38]. В общем случае историей нейтрона (запаздывающего, мгновенного и самого внешнего источника) называется его путь от рождения до гибели в результате поглощения (захват и деление) или утечки из области А (рисунок 15). В этом процессе, кроме координат, энергии и направления, отслеживается его время. Если время нейтрона

превысило заданное рассматриваемое время нестационарного процесса ^епа , то

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

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

Рисунок 14 - Схематическое представление истории нейтрона

Рисунок 15 - Пример цепочки нейтронов от внешнего источника

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

3.2 Верификация алгоритмов моделирования переходных процессов программой КИР на основании численного бенчмарк-теста RPIGS-KIN

3.2.1 Описание теста RPIGS-KIN Для проверки алгоритмов расчёта нестационарных функционалов и их среднеквадратичных отклонений использовалась тестовая модель «Rectangular Prismwith 1-Group Constants Critical» (RPIGS), описанная в [39] и представленная в

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

Таблица 9 - Макроскопические параметры теста КРЮБ-КШ

Физические параметры теста Значения

Макроскопическое сечение поглощения 0,4932 см-1

Макроскопическое сечение деления 0,18 см-1

Макроскопическое сечение рассеяние 0,5 см-1

Среднее число вторичных нейтронов 2,74

Средняя скорость нейтронов 22000 см/с

Время жизни мгновенных нейтронов 9,22 10-5 с

Таблица 10 - Параметры запаздывающих нейтронов

Р. 2,40212-10-4 1,23990-10-3 1,18372-10-3 2,65400-10-3 1,08811-10-3 4,55805-10-4

1,33512-10-2 3,26637-10-2 1,20919-10-1 3,03934-10-1 8,52665-10-1 2,86399

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

нейтронов е] ( = 0) = 0 в момент времени ? = 0 запускаются нейтроны внешнего

источника, далее рассматривается процесс до ? = 100 с, в котором физические параметры среды не меняются.

На момент времени ? = 0 система является подкритической (отсутствуют запаздывающие нейтроны), что приводит к уменьшению потока нейтронов на начальном временном интервале моделирования. В процессе взаимодействия нейтронов со средой образуются предшественники запаздывающих нейтронов,

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

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

Более того, данному процессу можно сопоставить численное решение модели одноточечной кинетики (21) с начальными условиями (22).

Р) _ Рд) -р(г) йг Л(г)

йс} (г) р (г)

+ ХЛ С (г) + )

Р(г)-Хс с 3 (г), ] = 1,2,...6.

^ йг Л(г) = 0) = 1; > 0) = 0; Р(г = 0) = 0; с(г = 0) = 0.

(21)

(22)

<

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

3.2.2 Результаты расчёта Для получения результатов расчёта по программе КИР использовалась временная сетка:

• 0 - 10-5 с (100 временных интервалов),

• 10-5 - 10-4 с (90 временных интервалов),

• 10-4 - 10-3 с (90 временных интервалов),

• 10-3 - 10-2 с (90 временных интервалов),

• 10-2 - 10-1 с (90 временных интервалов),

• 10-1 - 1 с (90 временных интервалов),

• 1 - 10 с (90 временных интервалов),

• 10 - 100 с (90 временных интервалов).

По программе КИР было проведено три расчёта с разными статистиками (106, 107 и 108 цепочек прародителя нейтронов) на основе оценки по столкновениям для исследования сходимости результатов при увеличении статистики.

На рисунке 16 представлены результаты расчёта изменения интегрального потока по программе КИР и точечной модели на основе численной схемы M3K2 [40-41] с относительной точностью 10-6.

10 поток, o.e. 1_ __

101 \

- KIR 10б историй нейтронов \

--- KIR 107 историй нейтронов \

--- KIR 10s историй нейтронов \

............- Уравнение точечной кинетики \

ю-2 ю-3

ю-7

Рисунок 16 - Результаты расчёта интегральной мощности потока теста RPIGS-KIN

В начальный момент времени нейтроны генерируются внешним источником [38]. Так как среднее время жизни мгновенных нейтронов составляет около 10-4 секунд, и система имеет нулевую мощность, происходит поглощение

нейтронов и резкий спад плотности потока нейтронов после 10-4 секунд. Происходит накопление ядер-предшественников запаздывающих нейтронов, и на начальном интервале времени, около 0,1 секунды, они начинают появляться и порождать новые мгновенные и запаздывающие нейтроны, что приводит, в конечном итоге, к стабилизации плотности потока, и позволяет остановить резкое уменьшение плотности потока. В результате процесса после 60 секунд система входит в критическое состояние, соответствующее физическим параметрам среды.

Наибольшее среднеквадратичное отклонение, полученное по программе КИР со статистикой расчёта 106 историй нейтронов, составляет 22,5 % во временной точке 0,645 с, а для статистик 107 и 10 8 историй нейтронов - 5,72% при ? = 0,715 ^ и 1,73 % при ? = 0,25 ^ соответственно.

На рисунке 17 представлены графики изменения плотности интегрального потока нейтронов на разных временных интервалах, на которых отмечены среднеквадратичные отклонения (полученные по программе КИР) в положительную и отрицательную сторону, а также результаты решения точечного уравнения кинетики. Результаты, полученные по программе КИР, разделены на временные интервалы в логарифмическом масштабе (10-6, 10-5, 10-4, 10-3, 10-2, 10-1, 1, 10, 100). Для рисунков со шкалами 10-3, 10-2 была добавлена вспомогательная ось для более наглядного понимания, так как в данной временной области среднеквадратичные отклонения графически совпадают с графиком точечной кинетики.

Рисунок 17 (а) - Результаты расчёта по программе КИР для разных временны интервалов: (а)- уравнение точечной кинетики, (б) - КИР со статистикой 106 историй нейтронов, (в) - КИР со статистикой 107 историй нейтронов, (г) - КИР со статистикой 108 историй нейтронов

Поскольку для моделирования процесса и получения результатов программа использует оценку по столкновениям, их число в интервалах порядка 10-6 мало, отклонения в данных интервалах принимают значения большие, чем в интервалах 10-5 и т.д. При возрастании временных интервалов число столкновений в регистрируемых интервалах увеличивается, и, следовательно, происходит уменьшение среднеквадратичного отклонения. С каждым следующим интервалом времени значение отклонения уменьшается. На временных интервалах с 0,1 по 100 секунд масштаб был увеличен для более наглядного понимания. Как уже было описано выше, начиная с 0,1 секунды, начинают появляться запаздывающие нейтроны, которые в свою очередь из-за малого количества вносят возмущения в среднеквадратичное отклонение потока нейтронов.

Для сопоставления результатов программы КИР и решения уравнений точечной кинетики использовалась методика с использованием доверительных интервалов (таблица 11). Данная методика заключается в подсчёте процентной доли временных точек «точного» решения, значения которых вошли в доверительные интервалы программы КИР (а ^ 68,3 %, 2а ^ 95,4 % и 3а ^ 99,72 %).

Таблица 11 - Доля точек, которые попали в решение точечной кинетики

Статистика расчёта по программе КИР а Кол-во точек 2а Кол-во точек 3а Кол-во точек

68,30% 499 95,40% 696 99,72% 728

106 историй нейтронов 72,33% 528 97,81% 714 99,59% 727

107 историй нейтронов 62,19% 454 93,42% 682 99,73% 728

108 историй нейтронов 71,92% 525 94,52% 690 99,18% 724

Полученные значения а хорошо согласуются с установленными процентными соотношениями, что в свою очередь указывают на корректный расчёт в программе КИР. Также в таблице 11 представлено число точек,

соответствующих процентной доли. Ста процентам соответствует 730 временных точек.

На рисунке 18 показано изменение среднеквадратичного отклонения результатов, полученных по программе КИР, при увеличении статистики. Каждый расчёт число историй нейтронов увеличивалось в 10 раз, а значение отклонения

уменьшалось в 710 раз, что соответствует определению средней квадратичной погрешности согласно теории вероятностей.

¡5 Отношение отклонений, отн.ед. 5

(а)

время, с

Ю-7 Ю-6 Ю-5 Ю-4 Ю-3 Ю-2 Ю-1 1 101 102

Рисунок 18 - Отношение среднеквадратичных отклонений результатов расчётов по программе КИР с разными статистиками: (а) - отношение 106 к 107 историй нейтронов, (б) - отношение 107 к 108 историй нейтронов.)

3.3 Верификация алгоритмов моделирования переходных процессов программой КИР на основании численного бенчмарк-теста С5С7-ТБ

3.3.1. Краткое описание тестов C5G7-TD Важным этапом при создании и усовершенствовании программ для нейтронно-физических расчетов для обеспечения безопасности ЯЭУ является их верификация на различных бенчмарк-тестах. В случае программ, направленных на решение стационарного уравнения переноса нейтронов, библиотека бенчмарков весьма обширна и давно сформирована, но в случае решения нестационарного уравнения переноса нейтронов, моделирования быстрых кинетических процессов, проблема малого количества бенчмарков является одной из актуальных. Одним из немногих бенчмарк-тестов, существующих на сегодняшний день, является численный бенчмарк тест C5G7-TD

Бенчмарк C5G7-TD представляет собой сборник тестовых задач пространственно-временной нейтронной кинетики с гетерогенным описанием пространственных областей для решения нестационарного уравнения переноса нейтронов без учета обратных связей [42]. Бенчмарк-тест C5G7-TD основан на хорошо известной спецификации, изданной в 2016 г. [36], описывающей стационарные транспортные тесты (C5G7). В настоящее время официальных результатов бенчмарка C5G7-TD не представлено. Но расчеты по данному бенчмарк-тесту введутся по всему миру [43-44].

Бенчмарк поделен на две области - тесты для двумерной и трехмерной модификации реактора. Для двумерной модели реактора представлено 4 варианта, которые включают в себя 19 тестовых задач. Для трехмерной модели представлено 2 варианта, которые включают в себя 9 тестовых задач [45]. При верификации алгоритмов моделирования кинетических процессов использовались двумерные тесты C5G7-TD. Геометрическая модель двумерной сборки, используемой в тесте C5G7-TD, полностью совпадает с двумерной моделью бенчмарк-теста C5G7, описанного во второй главе диссертации.

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

В данной главе представлены результаты расчета по программе КИР для двух рядов двумерных тестов: для тестов с мгновенным изменением положения органов регулирования и изменением плотности теплоносителя. Для оценки корректности работы алгоритмов моделирования кинетических процессов прямым методом Монте-Карло в групповом приближении используются результаты программ PANDAS-MOC (метод характеристик) и MPACT [46] (метод дискретных ординат). Результаты данных программ были оцифрованы с графиков приведенных в [47] и [43]. Следовательно, сравнение результатов данных программ с результатами программы КИР носят строго демонстрационный характер в связи с недостаточной точностью полученных значений. Детальный анализ данных, полученных по программе КИР, будет приведен в главе 4.

Для проведения нестационарных расчетов тестов C5G7-TD рассматриваемая сборка должна находиться в критическом состоянии. Для получения данного состояния был проведен предварительный стационарный расчет с целью получения эффективного коэффициента размножения системы со статистикой 100 млрд. историй нейтронов (1.186582 ± 0.000005).

При моделировании кинетических процессов программа КИР моделирует внешний источник мгновенных нейтронов, после чего происходит моделирование процесса выхода в равновесное состояние системы, в процессе которой происходит накопление ядер-предшественников запаздывающих нейтронов и образование корректного соотношения мгновенных и запаздывающих нейтронов для получения критичности. На рисунке 19 продемонстрированы результаты расчета выхода в критическое состояние рассматриваемой сборки, полученные по программе КИР в логарифмическом масштабе.

Для каждого теста рассчитываемое время кинетического процесса составило 10 секунд. Среднее время расчета варианта составляет 2150 минут с использованием распараллеливания расчета на 480 процессорах. Источником в момент времени I = 0 с было смоделировано 250 млн. мгновенных нейтронов. Выходными данными программы КИР служат интегральные значения скорости реакции деления на рассматриваемых интервалах времени.

Рисунок 19 - Изменение скорости реакции деления сборки во времени при выходе

в критическое состояние

3.3.2. Вариант 0 (ТОО) Данный вариант представляет собой тесты для двумерной задачи, моделирующее мгновенное погружение поглощающих стержней на глубину 10% активной зоны, и которые находятся в нем в течении 1-ой секунды. После чего стержни извлекаются на 5% высоты и в таком положении находятся еще 1 секунду. К концу 2-ой секунды переходного процесса стержни полностью извлекаются. На рисунке 20 представлено изменение доли вставленной части контрольных стержней во времени для данных тестов.

1.0 время, с

Рисунок 20 - Изменение доли вставленной части контрольных стержней во

времени для варианта TD0

Данные тесты возможно смоделировать физическим модулем путем мгновенного изменения материалов в зонах нахождения СУЗ-ов и описать формулами (23).

2 х (г) = 20т, г = о,г > 25,

2х (г) = ГГ + 0.1(2Х -2°хт ),0 < г < 15, 2х (г) = + 0.05(2хХ - 20т),15 < г < 25,

(23)

где

2

от

макроскопические сечения канала, при извлеченных поглощающих

стержнях -х

х - макроскопические сечения стержней СУЗ.

Тесты варианта TD0 различаются между собой группами вводимых органов регулирования. Нумерация групп органов регулирования была представлена в

х

бенчмарке C5G7 (рисунок 7). Конфигурации тестовых задач для варианта TD0 представлены в таблице 12.

Таблица 12 - Тестовые задачи для вариантов TD0 - TD2

Номер тестовой задачи Описание тестовой задачи

1 Ввод/вывод контрольных стержней группы 1

2 Ввод/вывод контрольных стержней группы 3

3 Ввод/вывод контрольных стержней группы 4

4 Ввод/вывод контрольных стержней групп 1, 3, и 4

5 Ввод/вывод контрольных стержней групп 1 - 4

На рисунках 21 и 22 представлены результаты тестов варианта TD0 программы КИР. Также на рисунках приведены результаты программ PANDAS-MOC и MPACT. Результаты по программе MPACT приведены только для 6-ти первых секунд.

При анализе результатов можно заметить, что скорости реакции деления, полученные для теста №2 и №3, имеют значительный разброс значений, что объясняется малым вносимым возмущением из-за изменения положения органов регулирования в периферийных ТВС. Следовательно, для расчетов изменения положения органов регулирования в периферийных ТВС требуется увеличить статистическую точность расчета, чтобы избежать значительных флуктуаций результатов.

При сопоставлении результатов среднее отклонение программы PANDAS-MOC от программы КИР не превышает 6%, программы MPACT от КИР не превышает 7.5% для всех тестов варианта TD0.

1.0

0.8

ш

о «

В

п ч 8 8 8

ш а л н о о а о и

и

0.6 -

0.4 -

0.2

0.0

КИР

РЛММБ МРАСТ

• ТБ0-1

□ □

ТБО-2 ТБ0-5

9

10

время, с

Рисунок 21 - Результаты расчета варианта TD0 (1-го, 2-го и 5-го тестов)

1.0

0.8 -

ш о

ш ч 8 8 8 И й

Л

н о о а о и

и

0.2

0.0

алд.——в-

КНР РАМ) ДБ МРЛС'1 - □ ♦ ТС 0-3

□

ТВ 0-4

4 5 6

время, с

10

0

1

2

3

4

5

6

7

8

й 0.4 -

0

1

2

3

7

8

9

Рисунок 22 - Результаты расчета варианта TD0 (3-го и 4-го тестов)

3.3.3. Вариант 3 (ГО3) В данном двумерном варианте моделируется процесс изменения плотности теплоносителя во времени. В начальный момент времени плотность теплоносителя во всех тепловыделяющих элементах равняется своему номинальному значению - р . В первую секунду плотность теплоносителя линейно уменьшается до значения сор , после чего в течении 2-ой секунды линейно возрастает до начального значения. Вариант состоит из 4-х тестов, с различным значением с . На рисунке 23 приведена зависимость изменения плотности теплоносителя от времени для всех тестов варианта TD3.

время, с

Рисунок 23 - Изменение плотности теплоносителя во времени для тестов

варианта TD3

Как и в предыдущем варианте, программа КИР сравнивается с результатами программ PANDAS-MOC и MPACT. Результаты данных программ и программы КИР приведены на рисунке 24.

Результаты программы КИР согласуются с результатами программ. Для всех четырех тестов характерен разброс значений во временном интервале [2;3] секунд

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

Среднее отклонение программ МРАСТ и РАМОАБ-МОС не превышает 5%. Проведение более детального и подробного анализа сравнения результатов данных программ не возможен в связи с их отсутствием, а данные скоростей реакции деления, полученные для программ МРАСТ и РАМОАБ-МОС, были оцифрованы, как уже упоминалось ранее. Более детальное сравнение и верификация методик расчета будут представлены в четвертой главе диссертации.

1.0 ш-

]

0.0 -I-,-,-,-,-,-,-,-,-,-

0123456789 10

Время, с

Рисунок 24 - Результаты расчета варианта ТОО (3-го и 4-го тестов)

Вывод к главе 3

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

Результаты теста RPIGS-KIN показали, что программа КИР с групповым модулем MAGMA корректно моделирует процесс выхода в критическое состояние от внешнего источника мгновенных нейтронов. Показано, что при увеличении статистики расчета по программе КИР, результаты программы сходятся к решению уравнения точечной кинетики.

Анализ результатов моделирования нестационарных процессов на тесте C5G7-TD и сравнение результатов с результатами программы MPACT и PANDAS-MOC показали корректность и адекватность построенных методик расчета. Результаты программ согласуются друг с другом.

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

ГЛАВА 4. АЛГОРИТМЫ МОДЕЛИРОВАНИЯ КИНЕТИЧЕСКОГО ИСТОЧНИКА НЕЙТРОНОВ МЕТОДОМ МОНТЕ-КАРЛО.

4.1 Описание метода

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

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

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

Рассмотрим тест BSS-6 (Benchmark Source Situation) [48-50], который представляет собой трехзонную размножающую пластину (рисунок 25), где зоны 1 и 3 эквивалентны между собой. Физические свойства материалов представлены в виде двух-групповых констант, а характеристики ядер-предшественников

Рисунок 25 - Геометрия теста ВББ-б

На рисунке 26 представлен выход данной модели в критическое состояние от внешнего источника мгновенных нейтронов в течении 100 секунд. При расчете использовался временной шаг 0.05 секунд. Для более детального рассмотрения и фиксации временного момента достижения критического состояние был проведен анализ изменения значений скоростей реакции деления, а именно отклонение скорости реакции деления на временном интервале I от ¡-1 (рисунок 27). Для более наглядного представления было проведено усреднение результатов на данных временных интервалах: 0.5 секунд (10 точек), 5 секунд (100 точек) и 10 секунд (200 точек). Исходя из представленных результатов, критическое состояние данной модели достигается около 60 секунд, а дальнейшие расхождения носят статистический характер, свойственный методу Монте-Карло.

Следовательно, для небольшой модели, такой как В8Б-6, процесс накопления ядер-предшественников запаздывающих нейтронов и выхода в равновесное состояние требует предварительного расчета 60 секунд кинетического процесса, что сильно удлиняет время счета или увеличивает затраты ресурсов. При рассмотрении более сложных моделей, таких как C5G7-TD, требуемое время составляет около 90 секунд (рисунок 19).

1.000

<и о