Развитие прецизионных и инженерных методов и программ расчета ядерных реакторов с использованием алгоритмов Монте-Карло

Как правило, для источника используется изотропное приближение, т.е. считается, что нейтроны деления рождаются изотропно в лабораторной системе координат. Разложение в ряд потока нейтронов характеризует анизотропию потока, в то время, как разложение сечения рассеяния - анизотропию рассеяния. На практике суммирование по l обрывают при достижении заданной величины L, поэтому в транспортном уравнении суммирование проводится: по индексу l _ от 0 до L, а по индексу m _ от -l до l. При этом говорят, что транспортное уравнение решается в P_L или в P_N приближении, где N = L.

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

, v=0,3.(2)

Здесь приняты следующие обозначения:

i,j - номера геометрических зон, i=1,NR;

g - номер энергетической группы, g=1,NG;

_ угловые гармоники плотности столкновений, ;

_ угловые гармоники потока нейтронов;

_ полное сечение в зоне i в группе g;

V_i - объем зоны i.

Угловые гармоники потока нейтронов при l=0,1 имеют ясный физический смысл:

- скалярный поток нейтронов Ф;

, - x- и y- компоненты тока нейтронов J_x , J_y ;

- z- компонента тока нейтронов J_z .

_ источники нейтронов в группе g в зоне i , появившиеся как в процессе деления (S_f), так и в процессе рассеяния из других групп:

.(3)

Т.к. принято, что распределение нейтронов деления изотропно в лабораторной системе координат, все высшее угловые гармоники источника деления равны нулю. При v=l=m=0 источники деления определяются по формуле:

(4)

здесь - спектр деления изотопа k , содержащегося в материале зоны i.

Анизотропные вероятности определены формулой:

,(5)

где u=0,1,2,3; v=0,1,2,3; l=0,L; m=0,l; l'=0,L ;m'=0,l'.

, v=0,1,2,3; l=1,L; m=1,l - функция двух переменных, зависящая от азимутального и полярного углов :

,(6)

где:

_ полиномы Лежандра;

_ присоединенные функции Лежандра;

Из формулы (4) видно, что в P₁ приближении при l=1, m=1:

.

Коэффициенты , определены по формулам:

.(7)

Количество неизвестных в уравнении (2) равно N_вектNRNG, где N_вект - число угловых гармоник потока нейтронов, определенное порядком P_L приближения. Для трехмерной (3D) и двумерной (2D) геометрий N_вект вычисляется по формулам:

N_вект(3D) = (1+L)²

N_вект(2D) = (1+L) (1+L)/2

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

Тензор обобщенных вероятностей симметричен относительно главной диагонали:

.

Имеет место следующее соотношение взаимности:

.

Для изотропных вероятностей первых столкновений P_ij, сумма строки матрицы равна единице: . Для обобщенных вероятностей произведение коэффициента и суммы обобщенных вероятностей равно единице только для диагональных элементов тензора:

.

Во всех остальных случаях это произведение равно 0.

В разделе 2.2 описаны алгоритмы расчета вероятностей первых столкновений. В настоящее время в программе MCU-FCP изотропные вероятности первых столкновений вычисляются для двумерной и трехмерной геометрий, обобщенные вероятности первых столкновений - только для двумерной геометрии. Интеграл по объемам V_i и V_j в формуле (5) сводится к двойному интегралу:

,(8)

где

. (9)

В формуле (8) предполагается, что проекция системы на плоскость, которая перпендикулярна оси Z, окружена прямоугольником. Основание этого прямоугольника перпендикулярно направлению траектории, которое определяется азимутальным углом .

Используются следующие обозначения:

H() _ длина основания прямоугольника, заключающего в себе проекцию системы на плоскость XOY;

h _ координата, характеризующая начала траектории в одной из точек этого основания; S_i - площадь зоны i;

_ij(h,) _ оптическая длина пути частицы от выхода из зоны i до входа в зону j для траектории, начинающейся в точке h отрезка H() в направлении ; _i, _j _ оптические длины отрезков той же траектории в зонах i, j;

Ki_n _ функции Бикли, которые определяются по формуле:

.

Двойной интеграл в формуле (8) вычисляется методом Монте-Карло со специальным выбором узлов сетки интегрирования. С использованием геометрического модуля из пакета программ MCU проводятся N лучей (в западной литературе - ray tracing).

Для каждого из N лучей, направление траектории и координата начала траектории h выбираются с использованием специальных сеток, например ЛП _ последовательностей Соболя. Это позволяет существенно ускорить вычисление интеграла. Так, для двумерной ячейки РБМК необходимо провести всего около 1000 лучей.

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

Для всех зон i, j, которые пересекает луч n, по формуле (9) вычисляются величины T_ij. По окончании счета обобщенные вероятности определяются по формуле:

.

В формуле (9) используются следующие обозначения:

С, b_n - константы, зависящие от порядка P_L приближения;

_ функция, зависящая от азимутального угла :

.

Для каждой пары i,j, можно оценить размерность тензора анизотропных вероятностей N_тенз в двумерной системе. В Таблице 1 приведена размерность тензора с учетом свойства симметрии относительно главной диагонали, т.е. количество различных элементов в тензоре.

Таблица 1. Размерность тензора обобщенных вероятностей для двумерной геометрии

L	Nвект	Nтенз
1	3	6
2	6	21
3	10	55
4	15	120
5	21	231

В разделе 2.3 описаны алгоритмы решения транспортного уравнения методом вероятностей первых столкновений с учетом анизотропии рассеяния.

Уравнение (2) решается стандартным методом итерации источника. На каждой итерации решение ищется отдельно для тепловой и эпитепловой области энергий, после чего источник нейтронов деления пересчитывается.

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

Метод GMRES (generalized minimal residual algorithm) используется для решения системы линейных уравнений с разреженной несимметричной матрицей вида: Ax=b, и является одним из методов подпространств Крылова. Следует отметить основные особенности алгоритма GMRES.

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

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

- Достигается существенная экономия оперативной памяти компьютера: для обращения матрицы nn требуется всего k векторов размерностью n. Так, для задачи с n=10⁵ используется значение k=25.

В пакете MCUFCP используется подпрограмма GMRES из свободно-распространяемой библиотеки LAPACK (Linear Algebra PACKage).

В разделе 2.4 описан статус пакета прикладных программ MCUFCP, в котором для решения транспортного уравнения переноса нейтронов используется метод вероятностей первых столкновений.

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

Приводятся результаты верификации применительно к расчету ячеек и полиячеек РБМК.

Константное обеспечение пакета MCUFCP составляет банк ядерных данных FCPDAT, который содержит информацию для 282-х изотопов.

В него входят только те библиотеки банка DLC/MCUDAT-2.2 пакета программ MCU, которые содержат сечения в групповом представлении. Это библиотеки:

БНАБ/MCU - расширенная и модифицированная версия 26-групповой системы констант БНАБ-93;

ТЕПКОН - многогрупповые сечения в области термализации;

ABBNL - 40 групповые сечения, которые используется для получения сечений “суммарного изотопа” при решении задач выгорания;

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

Кроме того, потребовалось создать новые специализированные библиотеки, содержащие микросечения в многогрупповом представлении: ELSM, ELSMT, FSP и GRC.

Методики и программы генерации этих библиотек были разработаны автором диссертации. Библиотеки содержат:

ELSM угловые моменты (l=0,1,...,5) сечения рассеяния упругих переходов в эпитепловой области энергий, нормированные на полное сечение рассеяния в группе;

ELSMT угловые моменты (l=0,1,...,5) сечения рассеяния упругих переходов в тепловой области энергий, нормированные на полное сечение рассеяния в группе;

FSP групповые спектры деления актиноидов во всей области энергии;

GRC - групповые микросечения в резонансной области энергии, рассчитанные методом Монте-Карло.

Матрицы упругого рассеяния библиотек ELSM и ELSMT рассчитывались методом Монте-Карло с использованием автономного запуска сегмента PHEAN из составного физического модуля пакета программ MCU.

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

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

В эпитепловой области начальная энергия нейтрона разыгрывалась по стандартному спектру, принятому в системе констант БНАБ, а именно: в первых трех энергетических группах - спектр деления, в остальных группах вплоть до 1 эВ - спектр Ферми. В тепловых группах, как это принято в программе TERMAC, использовался спектр, определяемый функцией Ферми, которая обеспечивает непрерывный переход от спектра Максвелла к спектру Ферми без разрыва первой производной по энергии.

Далее определялся номер энергетической группы, в которую попадает энергия нейтрона после рассеяния - в эпитепловой области используется групповое разбиение системы констант БНАБ, в тепловой области используется групповое разбиение библиотеки ТЕПКОН.

Для каждой начальной энергетической группы сегмент PHEAN вызывался NTOT раз. Параметр NTOT выбирался таким, чтобы статистическая ошибка расчета не превышала 0.01%. Для всех изотопов NTOT=10⁷.

Библиотека FSP содержит спектры деления делящихся изотопов. Энергетическая зависимость спектра описывается аналитической формулой (спектр Уатта), параметры которой для каждого делящегося изотопа выбирались из библиотеки БНАБ/MCU. Расчет спектров деления проводился с использованием стандартной подпрограммы численного интегрирования с автоматическим выбором шагов.

Библиотека GRC является единственным проблемно-ориентированным разделом в банке данных MCUDB-5.0. Групповые сечения изотопов рассчитываются методом Монте-Карло по программе MCU-REA/1 в системе, характерной для конкретного канала РБМК.

Для изотопов, содержащихся в топливных каналах (далее ТК) РБМК, рассматривается двумерная ячейка с ТК. Для изотопов, содержащихся в нетопливных каналах (далее НК), рассматривается двумерная полиячейка 3х3, в центре которой располагается НК, окруженный восемью ТК.

Групповые сечения зависят от типа канала. Например, в ТК с урановым топливом с обогащением 2.0% сечения изотопов в топливе будут отличаться от сечений в ТК с уран-эрбиевым топливом с обогащением 2.6%.

Кроме того, сечения одного и того же изотопа могут отличаться в различных элементах конструкции канала (сечения U-235 в топливе 1-го и 2-го ряда твэлов, сечения Zr в трубе ТК и в оболочках твэлов, и т.д.).

Поэтому при использовании библиотеки GRC в константном модуле MCUFCP действуют следующие соглашения:

· Библиотека GRC состоит из файлов с расширением;

· Имя файла библиотеки GRC совпадает с именем изотопа в формате MCU, расширение файла определяет тип канала, например, файл U235.U24 содержит сечения U235 в ТК с обогащением топлива 2.4%, файл U235.E26 _ сечения U235 в ТК с обогащением эрбиевого топлива 2.6%;

· Для каждого материала (физической зоны в терминологии MCU) указывается тип канала (расширение файла из библиотеки GRC, откуда будут читаться сечения резонансных изотопов, находящихся в данном материале). По умолчанию расширение файла сечений `GRC';

· Для каждого резонансного изотопа в материале имеется возможность указать номер секции в файле GRC. Различные секции соответствуют разным элементам конструкции канала. По умолчанию, номер секции совпадает с номером материала в исходных данных. Например, сечения изотопов в топливе 1-го ряда твэлов располагаются в секции IMAT=1, 2-го ряда - в секции IMAT=2.

С использованием данных соглашений каждому резонансному изотопу в каждом материале однозначно ставится в соответствие секция в файле GRC. Таким образом, секция файла GRC содержит групповые сечения изотопа из определенного элемента конструкции канала (номер секции) в определенном типе канала (расширение файла GRC). Данные сечения зависят от трех параметров активной зоны РБМК: концентрации изотопа (выгорание), температуры топлива и физической плотности теплоносителя.

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

В настоящее время создана промышленная версия библиотеки GRC, которая содержит сечения резонансных нуклидов, входящих в материалы практически всех существующих на сегодняшний день типов ячеек РБМК, а также пилотная версия библиотеки для ВВЭР.

Использованные подходы подготовки данных в библиотеках ELSM, ELSMT, FSP и GRC обеспечивают модели описания взаимодействия нейтронов с ядрами среды сходные с теми, которые используются в программе MCU-REA/1, что позволяет исключить константную составляющую погрешности по сравнению с MCU-REA/1. Это оказывается полезным при оценке вклада приближения плоских потоков в погрешность расчета, который зависит от детализации разбиения системы на регистрационные зоны.

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

модуль управления организует совместную работу всех модулей;

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

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

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

транспортный модуль решает транспортное уравнение переноса нейтронов методом ВПС;

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

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

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

Необходимо отметить, что в пакете программ MCUFCP используется ряд модулей из пакета MCU-4, например, управляющий модуль; групповые подмодули физического модуля; геометрический модуль; модули финальной обработки, обеспечивающие интерфейс с модулем расчета выгорания; модули расчета выгорания.

Из модулей пакета MCUFCP собираются специализированные рабочие программы MCU-FCP/1 и MCU-FCP/2.

MCU-FCP/1 программа расчета нейтронно-физических характеристик (далее НФХ) в процессе кампании в ячейках и полиячейках РБМК в дву- и в трехмерной геометрии.

Анизотропия рассеяния учитывается в транспортном приближении. Программа прошла процедуру верификации: для различных рабочих каналов по программам MCU-REA/1 и MCU-FCP/1 были посчитаны около семисот вариантов.

Варианты отличались выгоранием топлива, температурой топлива и плотностью теплоносителя.

Максимальное отклонение от значений эффективного коэффициента размножения (далее _ K_eff), полученных по MCU-REA, составило 0.3%.

Отклонения в эффекте Доплера и эффекте обезвоживания находятся в пределах 5% и 10% соответственно.

Сравнение энергетического и пространственного распределения нейтронов в ячейке с РК также показало хорошее согласие с результатами расчетов по MCU-REA.

В настоящее время программа MCU-FCP/1 используется для подготовки и уточнения библиотек малогрупповых констант для реакторных программ нейтронно-физического расчета РБМК SADCO и ТРОЙКА в НИКИЭТ и во ВНИИАЭС, соответственно.

Приводятся некоторые численные результаты, полученные по новой программе MCU-FCP/1. Рассматривается расчетная модель ячейки РБМК-1000 (см. Рис. 4), угол симметрии 30°. Особенность реализованных алгоритмов метода ВПС - приближение плоских потоков. Была проделана работа по выбору разбиения ячейки РБМК на регистрационные зоны.

Рассматривались несколько вариантов. На Рис. 5 показан вариант с 41-й регистрационной зоной.

Отклонение теплового потока в топливе, полученного методом ВПС от реперного значения, полученного методом Монте-Карло, составляет 0.3%. Расхождение в коэффициенте размножения - 0.3%. Время счета ~ 5 сек.



Рис. 4. Ячейка РБМК-1000. Угол симметрии 30

1,2 - топливо 3 - оболочка твэл 4 - центральный стержень 5 - теплоноситель 6 - труба ТК 7 - замедлитель

Рис. 5. Расхождение в потоках тепловых нейтронов, полученных методом ВПС и методом Монте-Карло, %

MCU-FCP/2 программа расчета НФХ в процессе кампании в ячейках и кассетах ВВЭР в дву- и в трехмерной геометрии. Для двумерной геометрии анизотропия рассеяния учитывается в P_N приближении до P₅ включительно.

Рабочие программы MCU-FCP/1 и MCU-FCP/2 используются вместе с банком данных FCPDAT, который содержит библиотеки сечений в многогрупповом представлении.

Модуль расчета вероятностей ВЕПС и транспортный модуль ПЕРСТ были включены в программу ТВС-М в качестве альтернативных модулей для расчета пространственного распределения нейтронов в кассетах ВВЭР. Эти модули также используются для разработки программы ТВС-КВАДРО спектрального расчета реакторов PWR и BWR. В настоящее время ведутся работы по тестированию и верификации программы MCU-FCP/2, модулей ВЕПС и ПЕРСТ применительно к расчетам НФХ в кассетах ВВЭР. Планируется, что программа MCU-FCP/2, модули ВЕПС и ПЕРСТ в составе программ ТВС-М и ТВС-КВАДРО будут применяться для подготовки и уточнения малогрупповых констант, которые используются в инженерных программах расчета активных зон реакторов типа ВВЭР, PWR и BWR.

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

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

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

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

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

Метод Бенуа с утечкой нейтронов заданной геометрическим параметром, используется для вычисления коэффициентов диффузии по программе MCU-REA/1.

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

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

, , (10)

где ,, - средний квадрат проекции смещения нейтрона по координатным осям.

Этот способ получения коэффициентов диффузии просто распространяется на малогрупповой случай, но вместо сечения поглощения в формуле (10) следует использовать сечение увода из группы.

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

Алгоритм L² заключается в следующем. Рассмотрим сначала бесконечную решетку одинаковых ячеек. Проекция смещения нейтрона на одну из осей, например, Х претерпевшего от момента рождения до момента поглощения N столкновений с ядрами среды , и квадрат этой проекции смещения , тождественно равны:

; ,(11)

где x_i _ проекция на ось X пробега нейтрона от точки столкновения с номером (i _1) до точки столкновения с номером i; n - число столкновений между i-м и j-м пробегом нейтрона, n N - параметр, задаваемый в исходных данных.

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

.(12)

Коэффициенты диффузии, полученные методом М ², в дальнейшем будем обозначать D(М ²), методом среднеквадратичных пробегов - D(L²). Из формулы (11) следует, что для бесконечной решетки одинаковых ячеек и при учете всех перекрестных членов в сумме методы М ² и L² эквивалентны, т.е. D(L²) D(М ²).

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

,(13)

где J - номер ячейки. Будем рассматривать суммы S₁ и S₂ отдельно. Для их вычисления предлагается следующий алгоритм.

Сумма S₁(J) формируется следующим образом: имеем пробег l_i, с двумя точками - начальной и конечной. Номера ячеек, в которых они находятся, известны. Тогда в сумму S₁(J) добавляется , если начальная точка пробега принадлежит ячейке J. Если конечная точка пробега принадлежит ячейке K, в сумму S₁(K) добавляется также . Если обе точки лежат в одной ячейке, то в соответствующую сумму добавляется .

В сумме S₂(J) фигурируют два пробега l_i и l_j, определяемые четырьмя точками - двумя начальными и двумя конечными. Номера ячеек для этих точек известны. В сумму S₂ каждой из этих четырех ячеек добавляется . При таком способе формирования сумм S₁(J) и S₂(J) для бесконечной решетки одинаковых ячеек любой сложности коэффициенты диффузии D(L²) и D(М ²) будут в точности равны.

Предложенный алгоритм реализован автором в программном модуле USER_L2 в составе программы MCU-REA/1 и может быть применен для таких задач, как:

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

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

- оценка влияния угловых корреляций между пробегами нейтрона на анизотропию диффузии.

Рассмотрим применение предложенного алгоритма для решения третьей задачи применительно к уран-графитовым реакторам. Для численного моделирования были выбраны характерные ячейки, вычисление тензора диффузии в которых представляет трудности. Для разделения корреляций, возникающих из-за гетерогенности среды, от корреляций, вызванных анизотропией рассеяния, применяли опцию программы MCU-REA/1 со сферически-симметричным рассеянием.

Статистическая ошибка рассчитываемых коэффициентов диффузии не превышала 0,1%.

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

Сравнивали D_r(L²), вычисленные при разных значениях числа столкновений между пробегами нейтрона (n) c точными коэффициентами диффузии D_r(М ²). Разница между ними и определяет вклад в анизотропию диффузии корреляций, вызванных гетерогенностью системы. Радиальный коэффициент диффузии определялся по формуле:

D_r=(D_x+D_y)/2.

Рассматривается бесконечная решетка квадратных ячеек уран-графитового газоохлаждаемого реактора. Для такой ячейки рассчитывались радиальные коэффициенты диффузии: D_r(М ²) и зависимость D_r(L²) от числа столкновений n между пробегами нейтрона. Ясно, что при достаточно большом n, отношение D_r(L²)/D_r(М ²) стремится к единице. Рассматривалось два типа ячеек а) и б), схемы которых показаны на Рис. 6.

В ячейке а) блок, состоящий из металлического урана, окружен кольцевым воздушным зазором. Радиус блока 2,2 см, воздушного зазора 5,4 см, шаг решетки 21 см. Материал замедлителя - графит.

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

Ячейка а)	Ячейка б)

Рис. 6 Схема ячеек а) и б): 1 - топливо; 2 - воздушный зазор; 3 - графит

Результаты расчетов представлены на Рис. 7. Для ячейки а) вклад угловых корреляций в радиальный коэффициент диффузии для тепловой группы немного превышает вклад для быстрой группы и составляет -7% (здесь и далее граница групп 1эВ).

Ячейка а)	Ячейка б)

Рис. 7. Зависимость отношения D_r(L²)/D_r(М ²) от числа столкновений между пробегами нейтрона для ячеек а) и б): 1, 2 - быстрая, тепловая группа соответственно

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

Алгоритм вычисления тензора диффузии методом среднеквадратичных пробегов, реализованный в программе MCU-REA/1, применяется для подготовки библиотек малогрупповых констант программы SADCO, предназначенной для проектных и эксплутационных расчетов реакторов РБМК. В частности, данный алгоритм применялся для расчетов коэффициентов диффузии ячеек пятого энергоблока Курской АЭС.

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

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

В Главе 1 диссертационной работы описана программа MCU-PR первая программа из семейства программ MCU с возможностью решать задачи выгорания. На основе программы MCU-PR была разработана, верифицирована и аттестована в Ростехнадзоре программа MCU-REA применительно к расчету НФХ реакторов типа ВВЭР с учетом выгорания топлива в процессе кампании. Верификация программы MCU-REA для решения задач выгорания проводилась непосредственно автором диссертационной работы.

В программах MCU-PR и MCU-REA для решения уравнений изотопной кинетики используется модуль BURNUP. Автором диссертации был разработан альтернативный модуль расчета выгорания ORIMCU. Возможность проведения расчетов изотопного состава с использованием двух альтернативных модулей выгорания в рамках одной программы позволяет повысить ее надежность и оценить методическую погрешность каждого из модулей, так как спектральная задача решается с использованием одного и того же алгоритма решения транспортного уравнения и одной библиотеки ядерных данных.

В данной главе дано описание применения программы MCU-REA с модулем ORIMCU для решения ряда прикладных задач, в том числе для определения радиационных характеристик облученного ядерного топлива (далее _ ОЯТ). Радиационными характеристиками ОЯТ называются такие величины, как активность, остаточное тепловыделение, интенсивность источников нейтронного и гамма-излучения, радиотоксичность.

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

Зоны могут включать 689 легких изотопов (например, входящих в состав конструкционных материалов), 129 актиноидов, 879 продуктов деления. Для каждого изотопа может учитываться до 7 типов реакций. В качестве расчетного блока изотопной кинетики на одном временном шаге в модуле ORIMCU используется известная программа ORIGEN-S из системы программ SCALE-4.3 (США). Программа ORIGEN-S в составе MCU-REA используется в РНЦ "Курчатовский Институт" начиная с 1999 года, когда было достигнуто соглашение с Окриджской лабораторией (США) об объединении компьютерных кодов MCU и ORIGEN-S.

В 2005 году программа MCU-REA/1, в которой используются модули расчета выгорания BURNUP и ORIMCU, была аттестована в Ростехнадзоре в качестве реперной для расчета НФХ реакторов типа ВВЭР в процессе кампании. Вопросы верификации программы MCU-REA/1 применительно к расчету НФХ реакторов ВВЭР в процессе кампании рассматриваются в Главе 5. Работы по верификации и аттестации программы проводились при непосредственном участии автора диссертации.

Раздел 4.2 посвящен использованию модуля ORIMCU в программах MCU_REA/1 и ТВС_РАД для расчета выгорания и радиационных характеристик ОЯТ реакторов ВВЭР.

Одной из задач обоснования безопасности хранения и транспортировки отработавшего топлива реакторов типа ВВЭР является расчет радиационных характеристик топлива выгружаемых кассет в зависимости от времени выдержки. Начиная с 1999 года, в РНЦ «Курчатовский институт» для этих целей широко используется программный комплекс MCU-REA/1 с модулем ORIMCU. Однако наличие реперной программы не отменяло необходимости иметь инженерную программу для выполнения расчетов радиационных характеристик отработавшего топлива реакторов ВВЭР. Поэтому вскоре был разработан программный комплекс ТВС-РАД, в котором в качестве спектральной программы, рассчитывающей распределение мощности и одногрупповые сечения нуклидов, использовалась программа ТВС-М, а в качестве модуля расчета выгорания _ модуль ORIMCU. Использование уже существующей программы в качестве модуля расчета изменения нуклидного состава топлива и его радиационных характеристик позволило создать необходимый расчетный инструмент в короткие сроки. Программа ТВС-РАД включена в инженерный программный комплекс КАСКАД, предназначенный для эксплутационных и проектных расчетов реакторов ВВЭР. В настоящее время программа ТВС-РАД проходит процедуру аттестации в Ростехнадзоре.

Верификация реперной программы MCU_REA/1 основана на сравнении с результатами двенадцати измерений остаточного энерговыделения десяти тепловыделяющих сборок (ТВС) американских реакторов PWR Turkey Point Unit 3 и Point Beach Unit 2. Выгорание топлива находилось в интервале от 25 до 39 МВтсут/кг. Суммарная величина отклонения расчета от эксперимента по всем двенадцати измерениям составила 1.2 1.6 %, что является очень хорошим согласием с экспериментом (1.6% _ удвоенное значение стандартного отклонения). Верификация инженерной программы ТВС-РАД основывалась на сравнении с результатами расчетов по другим программам, в том числе MCU_REA/1. На основании рассмотренных материалов сопоставления, можно сделать вывод о том, что погрешность расчетного прогнозирования радиационных характеристик ОЯТ не превышает 10% для полной активности и 15% для остаточного тепловыделения. Если же ограничиться рассмотрением долговременного хранения топлива (более 200 суток), то погрешность расчета остаточного тепловыделения ОЯТ не превышает 8-10%. В настоящее время программа ТВС-РАД проходит процедуру аттестации в Ростехнадзоре.

В разделе 4.3 описан предложенный и программно реализованный автором алгоритм численного моделирования выгорания оксидного смешанного топлива в ВВЭР-1000 с учетом двойной гетерогенности. Топливо состоит из двух различных материалов: обедненной урановой матрицы с малым количеством плутония, и агломератов - частиц с повышенной концентрацией плутония. Важной задачей является корректный учет двойной гетерогенности, позволяющий оценить отношение энерговыделения в агломератах и матрице к среднему энерговыделению по таблетке твэла. Приводятся численные результаты распределения выгорания, энерговыделения, а также концентрации гелия, криптона и ксенона по радиальным зонам твэла в зависимости от среднего выгорания топлива. Расчеты проводились по программе MCU-REA/1 с модулем расчета выгорания ORIMCU.

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

Микроструктура смешанного топлива с плутониевыми агломератами, которое предполагается использовать в ВВЭР-1000, на сегодняшний день точно неизвестна. Поэтому в предварительных исследованиях при учете двойной гетерогенности предполагалось отсутствие спектральной составляющей эффекта. Это означает, что спектр нейтронов в двухкомпонентной и в гомогенизированной среде совпадает. В этом случае нет необходимости в знании размеров агломератов и единственным параметром микроструктуры топлива, который используется в расчетах, является объемная доля агломератов в нем. Таким образом, задавшись минимальным набором характеристик топлива в качестве исходных данных, можно определить долю агломератов F_V и смоделировать выгорание смешанного топлива с плутониевыми агломератами по упрощенному алгоритму, который учитывает только геометрическую составляющую эффекта. Такое приближение позволяет использовать в расчетах стандартную аттестованную версию программы MCU_REA/1 с модернизированными интерфейсами модуля выгорания.

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

Cмешанное оксидное топливо изготавливают в два этапа. На первом этапе готовят смесь из диоксида обедненного урана обогащением по ²³⁵U 0,2 % по массе с добавлением некоторого количество диоксида плутония (около 20 % по массе). На втором этапе ее смешивают с диоксидом обедненного урана в пропорции, позволяющей получить требуемое среднее содержание плутония в таблетке. В результате образуется смешанное топливо из двух или более различных материалов, различающихся содержанием плутония: матрица со случайным расположением частиц с повышенным содержанием плутония - агломератов. Матрица состоит из смеси диоксида обедненного урана и малого количества плутония - около 1 % по массе.

В качестве исходных данных приняты следующие характеристики топлива:

_Т = 4,5 % по массе - среднее содержание плутония в таблетке,

_Т = 10,3 г/см³ - физическая плотность таблетки,

_А = 20 % по массе - содержание плутония в агломератах,

_М = 1 % по массе - содержание плутония в матрице,

k = 0,96 - отношение плотностей диоксида урана и плутония,

Х₅ = 0,2 % по массе - обогащение двуокиси урана по ²³⁵U.

На основании этих данных определяли объемную долю, плотность агломератов, плотность матрицы и среднее содержание ²³⁵U в таблетке по следующим формулам (14):

, , ,(14)

,

где - массовая доля плутония в агломератах и матрице соответственно.

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

Для расчета выгорания выбрана модель бесконечной решетки твэлов ВВЭР-1000, расположенных в узлах треугольной решетки с шагом 1,275 см, средняя плотность топлива принята равной 10,3 г/см³. При такой плотности топлива объемная доля агломератов в топливе составляет 18,4 %. Геометрические характеристики твэлов, изотопный состав материалов ячейки, температура и плотность теплоносителя, содержание бора в теплоносителе и средняя удельная мощность топлива взяты из описания международного расчетного теста “VVER LEU and MOX Computational Benchmark”, NEA/NSC/DOC(2002)10. Топливо в твэле разбивали на восемь радиальных зон. Выгорание рассчитывали с постоянной удельной мощностью топлива до среднего выгорания 60 МВтсут/кг.

Рассчитывали выгорание, относительное энерговыделение, концентрацию He, Kr, Xe в агломератах, матрице и в среднем по таблетке. Кроме того, вычислялись коэффициенты гетерогенности выгорания, которые в данной радиальной зоне равны отношению выгорания в агломерате (матрице) к среднему выгоранию по зоне. Аналогично определяли коэффициенты гетерогенности энерговыделения F_P, криптона, ксенона и гелия: F_Kr, F_Xe, F_He.

Коэффициент гетерогенности F является характеристикой эффекта двойной гетерогенности, который обусловлен наличием агломератов в топливе. Если свойства агломератов и матрицы одинаковы, то коэффициенты гетерогенности в агломератах и матрице равны единице. Чем больше коэффициент гетерогенности в агломератах, тем больше эффект двойной гетерогенности. Так, коэффициенты гетерогенности выгорания, энерговыделения, криптона и ксенона в агломератах в начале кампании достигают 4. По мере выгорания топлива свойства агломератов и матрицы сближаются, и коэффициент гетерогенности в агломератах уменьшается, в матрице _ увеличивается. При бесконечно большом выгорании они стремятся к 1 (см. Рис. 8а).

Характер зависимости от среднего выгорания в таблетке коэффициентов F_Kr , F_Xe такой же, как и коэффициентов F_B, F_P. Иначе изменяется F_He. Это объясняется тем, что криптон и ксенон образуются по одному каналу - как осколки деления, поэтому их концентрация в топливе практически пропорциональна мощности или выгоранию топлива. Гелий образуется по трем каналам: как осколок деления (1), в результате -распада трансурановых элементов (2) и в результате реакции (n,) на кислороде (3). Поэтому его накопление (4) носит нелинейный характер (см. Рис. 8б).

a)	б)

Рис. 8. а) Коэффициенты гетерогенности энерговыделения F_P в матрице (1) и в агломератах (2), б) накопления гелия в зависимости от выгорания топлива

С использованием описанного алгоритма по программе MCU-REA/1 с модулем ORIMCU были подготовлены константы для новой инженерной программы МОХТАВ расчета распределения выгорания, энерговыделения, а также концентрации гелия, криптона и ксенона по радиусу твэлов с оксидным смешанным топливом ВВЭР-1000. Указанные физические характеристики рассчитываются в зависимости от обогащения топлива, плотности топлива, среднего выгорания, числа радиальных зон в топливе, плотности теплоносителя и содержания бора в теплоносителе. Программа МОХТАВ передана во ВНИИНМ им. А.А. Бочвара для использования в составе программного комплекса СТАРТ-3, предназначенного для термомеханических расчетов при обосновании работоспособности проектируемых твэлов с оксидным смешанным топливом.

В разделе 4.4 приводится описание разработанного автором комплекса программ ORIPKD, который используется для численного моделирования радиационных полей от энергетических источников внутри защитной оболочки АЭС при аварии. Излагается методика расчета мощности дозы гамма-излучения от распределенного по пространству и энергии источника. Расчет мощности дозы реализован инженерным методом точечного ядра с использованием факторов накопления. Знание величины мощности дозы гамма-излучения в боксах системы локализации аварий необходимо для моделирования радиационно-химических реакций и поведения изотопов йода. Для расчета нуклидного состава отработавшего ядерного топлива используется программа MCU-REA с модулем ORIMCU.

Задача ставится следующим образом. При авариях на АЭС с ВВЭР из активной зоны в систему локализации аварий могут выходить радионуклиды, которые создают в боксах системы локализации поля ионизирующих излучений. Под влиянием этих полей в пространстве бокса происходят различные радиационно-химические реакции, которые влияют на поведение изотопов йода. Для расчета величины мощности дозы гамма-излучения (фотонного излучения) в боксах, необходимой для моделирования радиационно-химических реакций и поведения изотопов йода, была разработана автором инженерная программа PKD.

Программа PKD входит в программный комплекс ORIPKD, который позволяет проводить расчеты: накопления радионуклидов в топливе ядерных реакторов в зависимости от времени облучения; радиационных характеристик ОЯТ в зависимости от времени выдержки; мощности дозы гамма-излучения внутри защитной оболочки АЭС при аварии.

Программный комплекс ORIPKD является развитием программы MCU-REA с модулем ORIMCU. Он состоит из следующих программ:

- программа MCU-REA с модулем ORIMCU, предназначенная для расчета изменения нуклидного состава топлива в зависимости от времени облучения и определения радиационных характеристик ОЯТ в зависимости от времени выдержки;

- программа PKD, предназначенная для математического моделирования переноса фотонов методом точечного ядра и вычисления дозовых характеристик.

Задача определения дозовых характеристик ставится следующим образом. Из активной зоны и первого контура реактора ВВЭР в заданный момент времени после начала аварии в защитную оболочку АЭС выходит некоторое количество радионуклидов (актиноидов и продуктов деления). В результате этого в газовой фазе, на поверхностях защитной оболочки АЭС и в воде на полу помещений происходит накопление радионуклидов, являющихся источником гамма- излучения. Необходимо оценить радиационную обстановку в помещениях, находящихся внутри защитной оболочки АЭС. Программа PKD позволяет по известным значениям масс каждого из вышедших из первого контура радионуклидов определить мощность дозы гамма-излучения в каждом из помещений защитной оболочки АЭС.

Программа PKD была передана в Отдел Радиационной Безопасности Отделения ВВЭР для использования в работах по оценке радиационных последствий при проектных и запроектных авариях с учетом проектных решений по системе безопасности и локализации проекта АЭС с ВВЭР-1500.

Глава 5 посвящена верификации прецизионных и инженерных программ расчета реакторов.