Интегро дифференциальное уравнение переноса нейтронов

СОДЕЖАНИЕ ВВЕДЕНИЕ. 3 Глава 1 Обзор литературы Уравнение переноса нейтронов в интегро-дифференциальной форме Глава 2 Алгоритмы

Юлия Ширкова 5 лет назад Просмотров:

2 СОДЕЖАНИЕ ВВЕДЕНИЕ. 3 Глава 1 Обзор литературы Уравнение переноса нейтронов в интегро-дифференциальной форме Глава 2 Алгоритмы расчёта кинетики ядерных реакторов методом Монте Карло Начальные условия Вывод интегральных, нестационарных уравнений переноса нейтронов Алгоритмы моделирования кинетических процессов прямым методом Алгоритмы моделирования кинетических процессов в приближениях Алгоритмы модуля источников Глава 3 Описание программы Назначение и структура программы Описание геометрии расчетной системы Моделирование взаимодействий нейтронов с веществом Условия применения программы Изменение параметров реактора при моделировании динамического процесса Глава 4 Тестирование программы Тест RP1GC критический прямоугольный параллелепипед Моделирование кинетического процесса в критической системе Одногрупповой кинетический тест Тест RP1SC подкритический прямоугольный параллелепипед Тест RPCEU235 прямоугольный параллелепипед с изотопом U Тест BSS-6 одномерная модель реактора Расчет бесконечной решетки твэлов ВВЭР прямым и приближенными методами ЗАКЛЮЧЕНИЕ СПИСОК СОКРАЩЕНИЙ И УСЛОВНЫХ ОБОЗНАЧЕНИЙ СПИСОК ЛИТЕРАТУРЫ

3 ВВЕДЕНИЕ Актуальность работы. Перспективы развития ядерной энергетики напрямую зависят от надежности проектируемых реакторов. Постоянно повышающиеся требования к безопасности ЯУ влекут за собой и требования повышения точности расчетного обоснования характеристик ЯУ не только новых, но и существующих. Авария на АЭС «Фукусима-Дайичи» показала, что обоснование ядерной безопасности ЯЭУ требует повышенного внимания. При анализе ЯБ реакторов моделируются различные проектные и запроектные аварийные ситуации, т.е. фактически расчетным путем моделируются динамические процессы, протекающие в реакторной установке. Основной частью расчета динамики ЯР является его нейтронно-физическая часть — кинетика. Как правило, используются программы инженерного класса, реализующие приближенные методы решения нестационарного уравнения переноса нейтронов, к ним относятся: диффузионные приближения, приближения точечной или распределенной кинетики и др. В последнее время, в связи с бурным развитием вычислительной техники, стали появляться программы, реализующие транспортное приближение, метод дискретных ординат, метод поверхностных гармоник. Решение нестационарных уравнений переноса нейтронов методом Монте-Карло реализовано в ряде программ: TDMC (Иран) [1], SERPENT 2 (Финляндия) [2], TRIPOLI (Нидерланды) [3 8], TMCC (Китай) [9], TDMCC (Россия) [10, 11], а также корейская программа [13, 14]. В условиях недостаточного количества экспериментальных данных по динамике ЯР необходимо создание прецизионных программ, расчеты которых использовались бы как бенчмарки при кросс-верификации широкого круга инженерных программ, необходимость в которых обусловлена требованиями по оперативности получения расчетных характеристик ЯР. На основании всего вышесказанного становится очевидным актуальность и перспективность работы по созданию прецизионных программ, реализующих метод Монте-Карло для расчета динамики ЯР.

4 Цель диссертационной работы разработка алгоритмов и расчетных программ для решения нейтронно-физических пространственно-временных задач на основе метода Монте-Карло для повышения надежности, точности расчета нейтроннофизических нестационарных характеристик ядерных реакторов. Для достижения этой цели решены следующие задачи: 1) Разработан алгоритм прямого метода решения нестационарных уравнений, описывающих кинетику ЯР с учетом запаздывающих нейтронов, используя аналоговый метод Монте-Карло; 2) Получены интегральные уравнения переноса нейтронов в квазистатическом и усовершенствованном квазистатическом приближениях; 3) Разработан алгоритм расчета форм-функции при использовании квазистатического приближения; 4) Разработано программное обеспечение для расчета прямым и приближенными методами кинетики ядерных реакторов с использованием метода Монте-Карло; 5) Проведена апробация алгоритмов и программы на некоторых реперных расчетных моделях. Научная новизна. 1) Выведены интегральные уравнения переноса нейтронов в квазистатическом и усовершенствованном квазистатическом приближениях; 2) Разработаны оригинальные алгоритмы и программная реализация прямого моделирования кинетики ЯР методом Монте-Карло без использования весовых коэффициентов; 3) Разработаны алгоритмы и программное обеспечение приближенных методов моделирования кинетики ЯР. Практическая значимость. 1) Разработаны программы КИР и КИР-П, реализующие решение нестационарных уравнений переноса нейтронов методом Монте-Карло, которые используются в программном комплексе ДАРИЙ (Динамика Атомных Реакторов), предназначенном для моделирования динамических процессов, протекающих в активных зонах ядерных реакторов с жидкометаллическим теплоносителем. Для проведения теплогидравлических расчётов в этом комплексе используется программа ТЕИСП, разработанная в АО «НИКИЭТ». 4

5 2) Программный комплекс находится в опытной эксплуатации в НИЦ «Курчатовский институт». 3) Практическая значимость результатов работы подтверждена актом о внедрении в АО «НИКИЭТ». Основные положения, выносимые на защиту. 1) Интегральные уравнения переноса нейтронов в квазистатическом и усовершенствованном квазистатическом приближениях; 2) Алгоритмы расчета форм-функции в квазистатическом приближении; 3) Алгоритмы прямого моделирования кинетики ЯР методом Монте-Карло; 4) Результаты тестирования программ КИР и КИР-П. Апробация работы. Основные результаты диссертационной работы докладывались и обсуждались на следующих научных семинарах и конференциях: Школа-семинар по проблемам физики реакторов («Волга-2014»), 2014 г.; 9-я Международная научно-техническая конференция «Обеспечение безопасности АЭС с ВВЭР», Научно-техническая конференция «Нейтроника-2016» «Нейтронно-физические проблемы ядерной энергетики». Публикации. По теме диссертации опубликовано 4 статьи в реферируемых научных журналах из перечня ВАК РФ. Личный вклад автора. 1) Выведены уравнения переноса нейтронов в квазистатическом и усовершенствованном квазистатическом приближениях; 2) При непосредственном участии автора разработаны все представленные алгоритмы моделирования кинетики реакторов; 3) Программная реализация алгоритмов кинетики в адиабатическом и квазистатическом приближении; 4) Разработан алгоритм расчета форм-функции в квазистатическом приближении; 5) Разработка расчетных моделей и расчет тестовых задач, анализ полученных результатов; 5

6 6) Разработка новой тестовой задачи по расчету реакторов типа ВВЭР, которая может служить реперной при верификации других программ; 7) Разработка модуля источников и общего модуля управления; 8) Распараллеливание программы и оптимизация счета. Структура и объем диссертации. Диссертация состоит из введения, обзора литературы, трех глав и заключения. Работа содержит 93 страницы печатного текста, 32 рисунка и 11 таблиц. Библиография насчитывает 57 наименований. 6

7 Глава 1 Обзор литературы В скором времени после открытия явления [15], что при делении ядра урана в результате захвата им нейтрона часть нейтронов деления вылетает не мгновенно, а с некоторым запаздыванием, пришло осознание решающей роли запаздывающих нейтронов для создания установок с управляемой цепной самоподдерживающейся реакцией деления [16]. Последующие измерения показали, что мгновенные нейтроны испускаются за время порядка с, в то время как периоды полураспада предшественников запаздывающих нейтронов могут достигать значения более 50 с [17]. За прошедшие с тех пор десятилетия опубликовано большое количество работ, посвящённых исследованиям кинетики и динамики ядерных реакторов. Результаты данных исследований обобщены в нескольких специально этим вопросам посвящённых монографиях [17 20]. В книгах описана общая теория ядерных реакторов [21 23], а в учебных пособиях [24 26] по физике ЯР обязательно присутствуют главы по кинетике и динамике ЯР. Обсуждение вопросов кинетики и динамики реакторов в приложении к импульсным исследовательским ЯР содержится в работах [27 31]. «Область собственно кинетики реакторов объединяет явления, при которых плотность нейтронов и интенсивность делений в реакторе достаточно высоки, чтобы можно было не учитывать статические флуктуации, и в то же время не настолько велики, чтобы нужно было учитывать влияние реактивностной обратной связи, а область собственно динамики реакторов объединяет все случаи нестационарных процессов, при которых нейтронные переходные явления сопровождаются изменениями температуры среды реактора» [28]. До настоящего времени исследования кинетики и динамики реакторов проводятся с использованием тех или иных алгоритмов решения приближённых уравнений точечной или распределённой кинетики, ниже указаны некоторые из них: JAR-IQS [32], ГЕФЕСТ [33], PRISET-MBIR [34] программы использующие диффузионное приближение; LYCKY_TD [35], DORT-TD [36], TORT-TD [37], PARTISN [38]. программы использующие метод дискретных ординат; SUHAM [40], SVS [41] программы использующие метод поверхностных гармоник. В программном комплексе SVS также используется метод поверхностных псевдоисточников; 7

8 EVENT [39] программа использующая методы конечных элементов и сферических гармоник. В последнее время стали появляться работы с описанием программ, в которых для решения уравнения переноса нейтронов с временной зависимостью плотности нейтронов используется метод Монте-Карло. На данный момент в мире реализованы несколько программ по расчету кинетики реакторов методом Монте-Карло. К таким программам относятся коды: TDMC (Иран) [1]. Программа для расчета кинетики реакторов прямым методом (без приближений). Алгоритмы программы не описаны; SERPENT 2 (Финляндия) [2]. Программа для расчета кинетики и динамики реакторов прямым методом. Программа была расширена для расчета кинетики, добавив переменную времени для нейтронов. Как указано в работе, возможен расчет только быстро протекающих процессов на быстрых нейтронах, т.к. отсутствует учет запаздывающих нейтронов; TRIPOLI (Нидерланды) [3 8]. Программа для расчета кинетики и динамики реакторов прямым методом. Алгоритмы программы описаны достаточно подробно. TMCC (Китай) [9]. Программа для расчета кинетики реакторов прямым методом. Алгоритмы программы частично описаны; TDMCC (Россия) [10, 11] Программа для расчета кинетики реакторов прямым методом. Алгоритмы программы не описаны; TDKENO (Россия) [12] Программа для расчета кинетики реакторов в квазистатическом приближении, используя метод Монте-Карло. Алгоритмы программы не описаны; Корейская программа [13, 14]. В программе решаются уравнения усовершенствованного квазистатического приближения, а метод Монте-Карло используется для определения параметров точечной кинетики и вычисления формфункции. Алгоритмы программы описаны достаточно подробно; В большинстве из работ реализованные по перечисленным программам алгоритмы описываются очень кратко или совсем не описываются, и приводятся результаты некоторых тестовых расчётов. Исключением являются работы [3 8], в которых реализованные в программе TRIPOLI (Нидерланды) алгоритмы описаны достаточно подробно и [13, 14] с описанием алгоритмов корейской программы. 8

9 В программах TDMC, TRIPOLI, TMCC и TDMCC методом Монте-Карло решается временное уравнение переноса нейтронов с учётом запаздывающих нейтронов. Все эти программы используют весовые методы, т.е. вес нейтрона или вклад их в оцениваемые функционалы в течение кинетического процесса может меняться. Создание таких программ стало возможным благодаря развитию в последнее время многопроцессорной вычислительной техники [42]. Самый мощный в мире компьютер (информация на декабрь 2015 г. [42]) китайский Tianhe-2 (MilkyWay-2) имеет более 3 млн. 100 тыс. вычислительных ядер с тактовой частотой 2,2 ГГц и обладает производительностью

33,9 пфлоп оп./с (где 1 петафлоп (пфлоп) = 1 квадрильону = ; 1 эксафлоп (эфлоп) = 1 квинтильону = ), на втором месте американский Titan (более 560 тыс. ядер с тактовой частотой 2,2 ГГц, производительность

17,6 пфлоп оп./с), на третьем месте американский Sequoia (более 1,5 млн. ядер с тактовой частотой 1,6 Ггц, производительность

17,2 пфлоп оп./с), на четвёртом японский К (более 700 тыс. ядер с тактовой частотой 2,0 Ггц, производительность более 10 пфлоп оп./с). К 2020 г. в Японии планируют построить компьютер с производительностью 1 эфлоп оп./с, т.е. он будет иметь

ядер. В России (декабрь 2015 г.) самым мощным компьютером является «Ломоносов» -, который функционирует в МГУ им. М.В. Ломоносова и изначально был рассчитан на производительность 1 пфлоп оп./с. Сейчас он имеет более 37 тыс. вычислительных ядер с тактовой частотой 2,6 ГГц, его производительность

1,85 пфлоп оп./с и он занимает 22-ое место в TOP500 [42]. Исходя из того, что для оценки энерговыделения в регистрационной зоне приемлема статистическая погрешность 1%, тогда для покассетного расчёта распределения по объёму активной зоны реакторов слоями по высоте достаточно промоделировать 10 8 историй. Ниже даются некоторые приблизительные расчётные характеристик реакторов с жидкометаллическим теплоносителем. Оценка необходимого расчётного времени для моделирования процесса длительностью 1 с даётся в предположении, что для каждого временного отрезка длительностью, равной времени жизни нейтрона в реакторе, необходимо промоделировать 10 8 историй нейтронов. 9

10 Таблица Расчётные характеристики реакторов с жидкометаллическим теплоносителем. Расчётная характеристика Тип реактора БН-600 БРЕСТ-300-ОД Время жизни нейтрона, мкс

0,6 Средняя энергия нейтрона, вызвавшего деление, кэв

200 Произведение времени жизни нейтрона на его среднюю 2 3,5 скорость, м Число столкновений на историю

200 Время моделирования одной истории на вычислительном ядре с тактовой частотой 2,2 ГГц по

0,005 программе MCU, с Время моделирования 100 млн. историй на одном вычислительном ядре с тактовой частотой 2,2 ГГц, с Время моделирования 100 млн. историй на компьютере Titan при коэффициенте 0,5 1 распараллеливания близком к 1, с Время моделирования процесса длительностью 1 с на компьютере Titan 9 суток 1, с = 400 ч = 18 суток Приведённые в Таблице 1 оценки расчётного времени, необходимого для моделирования протекания динамического процесса в реакторах с жидкометаллическим теплоносителем, являются явно завышенными, т.к. получены при очень жёстком предположении о необходимости прослеживания 100 млн. историй нейтронов для каждого жизненного цикла поколения нейтронов, равного времени жизни нейтрона в реакторе. На практике может оказаться, что набирать такую статистику необходимо на временных интервалах, в раз бόльших, чем время жизни нейтрона, подобные ожидания основываются на следующих рассуждениях: в реакторе ВВЭР-1000 нейтрон испытывает в среднем

30 столкновений, его среднее время жизни около 4, с, и время моделирования одной истории по программе MCU составляет

1, с; в реакторе РБМК-1000 нейтрон испытывает в среднем

150 столкновений, среднее время жизни нейтрона деления

7, с, и время моделирования одной истории по программе MCU приблизительно такое же, как для БРЕСТ-300-ОД. Если сопоставить 10

11 эти данные с данными Таблицы В.1, то получится, что для моделирования динамики больших энергетических реакторов потребуется на 3 порядка меньше расчётного времени, чем для моделирования динамики с жидкометаллическим теплоносителем при условии рассмотрения процесса одинаковой длительности. Даже при том, что для набора достаточной статистики требуется промоделировать судьбу не 10 8 нейтронов, а 10 9 всё равно остаётся разница на 2 порядка, что выглядит довольно парадоксально. Объяснение такого парадокса состоит в том, что характерным «квантом» времени для реактора является не время жизни мгновенного нейтрона деления, а более длительный промежуток времени. Изменения в активной зоне происходят либо под действием извне (движение стержней, изменение расходов и т.д.), либо являются реакцией системы на вводимые изменения (обратные связи по температуре и плотности, нарушение геометрии конструкционных элементов и т.п.). Масштаб времени таких изменений давно определён детерминистскими динамическими расчётами и составляет

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

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

12 процессов, протекающих в активной зоне реактора МБИР, пересчёт форм-функции осуществляется через каждые 0,01 с. Если такая же технология расчётов будет применяться и при использовании программы КИР-П, то приведённые в последней строке Таблицы 1 оценки времени моделирования процесса длительностью 1 с на компьютере Titan необходимо уменьшить на 4-5 порядков. 1.1 Уравнение переноса нейтронов в интегро-дифференциальной форме Пусть: Φ(r,Ω,E,t) плотность потока нейтронов; C i (r,t) плотность предшественников запаздывающих нейтронов в группе i, где i = 1,2 N d ; λ i постоянная распада в группе i; χ i (E) энергетический спектр запаздывающих нейтронов группы i; χ p (E) нормированный спектр мгновенных нейтронов; ν(r,e) ожидаемое полное число нейтронов, испускаемых на одно деление в точке r, вызванное нейтроном с энергией E; β i (r,e) доля этого числа, отнесённая к i-й группе предшественников запаздывающих нейтронов; N d количество групп запаздывающих нейтронов (обычно N d = 6 или 8). Для упрощения записи приведённых ниже в этом разделе формул допустим, что λ i, β i, χ i не зависят от типа деления. Тогда система уравнений переноса нейтронов в интегро-дифференциальной форме имеет следующий вид [17,22,20]: где 12

13 (r,ω,e,t) фазовые координаты частицы; r радиус-вектор местоположения; Ω вектор направления полёта; E энергия; t время; Σ = Σ(r,E,t) полное макроскопическое сечение взаимодействия нейтронов; Σ x (r;ωʹ,eʹ Ω,E;t) = Σ x (r;ωʹ,eʹ)f x (r,ω,e Ω,E;t) дифференциальное макроскопическое сечение взаимодействия нейтронов с веществом; вероятность того, что при столкновении в точке нейтрона, имеющего направление полёта и энергию, в той же точке появится нейтрон с направлением движения и энергией (вероятность перехода); скорость нейтрона; Q(r,Ω,E,t) внешние источники;. Если чисто формально представить поток нейтронов в виде произведения амплитудного фактора n(t) и форм-функции, т.е. 414]: то можно получить уравнение переноса нейтронов для форм-функции [28 стр. 13

14 Решив неоднородное дифференциальное уравнение первого порядка (1.3), можно получить выражение для концентраций эмиттеров запаздывающих нейтронов: При записи потока нейтронов в виде произведения двух сомножителей [19,22], предполагается, что амплитудный фактор n(t) должен описывать основную временную зависимость, в то время как форм-функция должна слабо меняться со временем. Впервые метод был предложенный Генри [43] в 1958 г. Сопряжённая функция (или функция ценности) удовлетворяет следующему уравнению для критической системы где нулевой индекс означает, что данная величина относится к критическому состоянию 14

15 Форм-функция нормируется так, что выполняется соотношение где скорость нейтронов. Объёмный интеграл берётся внутри выпуклой поверхности, на которой определены граничные условия. Цель нормировки удовлетворение требованию Такую нормировку всегда можно осуществить [22]. Хотя функция отнормирована, уравнение (1.9) не определяет никакой нормировки для функций и n и, следовательно, амплитудный фактор n(t) может быть нормирован независимо любым подходящим для этого способом. Амплитудный фактор n(t) можно нормировать, положив его в момент времени t 0 равным полному числу нейтронов, присутствующих в это время в системе. Это в свою очередь определит нормировку форм-функции в момент t = t 0, и из уравнения (1.4) будет определена нормировка во все другие времена. Величина n(t) будет представлять полное число нейтронов в момент времени t, если форм-функция меняется незначительно. Можно отождествить n(t 0 ) с мощностью реактора в некоторый момент времени t 0. Функция n(t) остаётся почти равной мощности реактора, пока форм-функция не сильно отличается от начальной. Во многих практических случаях нормировка на мощность более удобна, особенно когда учитываются обратные связи, так как они определяются уровнем мощности реактора. При выводе уравнений, описывающих поведение во времени точечного реактора (кинетика), используется следующая процедура [44 46, 22]. Сначала уравнение (1.1) умножается на, а уравнение (1.7) — на. Результаты затем вычитаются и интегрируются по объёму, углам и энергии с учётом уравнения (1.9), которое 15

16 используется в члене, содержащем. Члены с градиентом затем уничтожаются с использованием теоремы Гаусса-Остроградского и граничных условий. Окончательный результат включает члены, описывающие источники мгновенных и запаздывающих нейтронов, и некоторые разности, например, между Σ и Σ0. Он может быть записан в виде где определены ниже. Кроме того, если уравнение (1.3) умножить на и проинтегрировать по всем переменным, то получим Уравнения (1.10) и (1.11) описывают кинетику реактора. Параметр ρ(t) называется реактивностью и определяется следующим соотношением: где Δ означает разности между соответствующими величинами в данный момент времени и в стационарном состоянии, т.е. ΔΣ = Σ Σ 0. Другие параметры в уравнениях (1.10) и (1.11) определяются выражениями: 16

17 Множитель в предыдущих выражениях более или менее произволен. Он не влияет на решения уравнений (1.10) и (1.11), так как всегда сокращается. Однако на практике значения F выбираются так, чтобы различные параметры допускали физическую интерпретацию в простейших ситуациях. В этом отношении разумным выбором для F является следующий [22]: Можно убедиться, что при таком определении F входящие в уравнение (1.10) величины имеют ясную физическую интерпретацию. Уравнения кинетики (1.10), (1.11) можно записать в несколько ином виде, используя вместо величин ρ и Λ непосредственно K эф и время жизни мгновенных нейтронов l [22]. Под временем жизни l мгновенных нейтронов в реакторе понимается величина, обратно пропорциональная скорости генерации ценности в эквивалентном критическом реакторе с числом вторичных нейтронов на деление νʹ(e)=ν(e)/k эф : 17

18 Время жизни мгновенных нейтронов и время генерации мгновенных нейтронов связаны соотношением: l = ΛK эф. (1.20) Параметры, введённые выше, являются в некоторой степени произвольными, вопервых, вследствие того, что функция не полностью определена, так как соответствующая критическая система отчасти произвольна, и, во-вторых, из-за выбора величины [22]. Различные параметры между собой связаны, поэтому необходима осторожность для того, чтобы обеспечить их согласование, но если согласованный выбор сделан, то приведённые выше параметры полностью определены. 18

19 Глава 2 Монте Карло Алгоритмы расчёта кинетики ядерных реакторов методом 2.1 Начальные условия Любой динамический процесс начинается с некоего начального состояния. Для определения параметров начального состояния решается однородное уравнение переноса нейтронов методом Монте-Карло и получается пространственное потвэльное или покассетное, с учётом высотной зависимости, распределение скорости реакции рождения нейтронов, с которого и начинается моделирование динамики. Предполагается, что перед началом процесса мощность реактора поддерживалась постоянной достаточно долгое время, чтобы нейтроны и предшественники запаздывающих нейтронов достигли равновесия, соответствующего этому уровню мощности, т.е. достигли стационарного состояния. Соотношения между потоком нейтронов и концентрациями предшественников в стационарном состоянии получаются из уравнения (1.3) в предположении, что производные по времени равны нулю: нейтронов. где суммарное количество предшественников запаздывающих Плотность нейтронов выражается формулой: где v скорость нейтрона. 19

20 Соотношения между стационарными средними величинами плотности нейтронов и концентрациями предшественников можно получить также из уравнения (1.11) в предположении нулевой производной по времени предшественников: Мощность реактора и число нейтронов в нём связаны соотношением [53]: n 0 = 3, νλp, где P мощность реактора в ваттах; ν среднее число нейтронов деления; Λ время генерации. Здесь предполагается, что при одном делении выделяется энергия 200 МэВ, тогда для скорости выделения энергии (мощности) 1 Вт требуется 3, дел/с. При этих условиях количество предшественников запаздывающих нейтронов в реакторе равно C = ΣC i = 3, P/Σλ i /β i. При мощности 1 Вт их количество порядка и не зависит от среднего времени жизни мгновенных нейтронов. Количество нейтронов в реакторе намного меньше и при одной и той же мощности зависит от времени генерации Λ. Для мощности 1 Вт в реакторах на быстрых нейтронах n составляет величину

Вывод интегральных, нестационарных уравнений переноса нейтронов Проведенный анализ работ [13 10] показывает, что в них рассматривается лишь нестационарное уравнение переноса нейтронов в интегро-дифференциальной форме. Используя метод Монте-Карло набирается статистика процессов в фазовом 20

21 пространстве, т.о. решается интегральное уравнение переноса нейтронов. Доказательство возможности перехода от интегро-дифференциального уравнения к интегральному уравнению Пайерса [48] приводится в монографии Владимирова [49]. Методом Монте-Карло естественным образом моделируются случайные процессы, происходящие в реакторе, которые в свою очередь описываются уравнением переноса нейтронов. Для реализации приближенных методов, необходимо рассмотреть, каким образом зависящие от времени параметры отражаются в кинетическом уравнении, записанном в интегральном виде. В книге Белла и Глесстона [22 стр. 21] приводится вывод уравнения (2.4) в интегральной форме стандартным методом характеристик. Вводится полная производная ППН следующим образом: s это расстояние вдоль направления перемещения Из интегро-дифференциального уравнения (1.1) нетрудно увидеть, что В итоге уравнение для ППН приводится к интегральному виду: где определяется формулами (1.2) и (1.3). Если расписать каждый член для источника нейтронов q и поделить уравнение (2.4) на (согласно произведению (1.4)), то можно получить уравнение для форм- 21

22 функции без каких-либо приближений: Выражение в показателе экспоненты в (2.5) — коэффициент ослабления. В адиабатическом и квазистатических приближениях весь рассчитываемый интервал времени делится на подинтервалы на которых форм-функция не зависит от времени, поэтому параметр s(r) на каждом из подинтервалов не зависит от времени. Производная форм-функции по параметру s вводится следующим образом: В таком случае производится замена переменных (2.3), а замена переменных (2.2) исключается. В усовершенствованном квазистатическом приближении обычно используется линейная апроксимация производной форм-функции по времени, которая рассчитывается следующим образом: 22

23 где значение времени, для которого последний раз пересчитывалась форм-функция. Возможны и другие аппроксимации производной форм функции [50 стр. 61], такие как экспоненциальная и квадратичная, которые также несложно использовать в расчетах. После произведения замены переменных (2.3) и (2.7) из уравнения (1.5) получается следующее уравнение: Уравнение (2.8) можно записать следующим образом, убирая полученную на предыдущем шаге функцию в источник, а член при внутрь квадратных скобок. 23

24 Уравнение представляет собой линейное дифференциальное уравнение первого порядка, которое может быть проинтегрировано с помощью интегрирующего множителя. Был выбран следующий множитель: домножив на него с обеих сторон равенства получаем: Полученное выражение далее интегрируется по от — до текущего, в предположении что 0 при. Это верно, т.к. когда в реакторе не был пущен не один нейтрон (т.е. не один из них не прошел хотя бы какой-то путь) форм функция равна 0. Тогда получаем: 24

25 Домножив обе стороны равенства на получаем: Если положить и заменить знак переменной интегрирования так, чтобы пределы интегрирования были (0, ) и (0, ), получаем: Если в уравнение подставить выражения для и (уравнения (2.9) и (2.10)), умножив и разделив правую часть уравнения на, можно получить интегральное уравнение для форм-функции в усовершенствованном квазистатическом приближении: 25

26 В квазистатическом приближении производная форм-функции по времени в уравнении (14) полагается равной нулю. Интегральное уравнение для форм-функции, полученное в квазистатическом приближении, выглядит следующим образом: Адиабатическое приближение вводит три упрощения [28 стр ]. В нем не различаются формы распределения источников запаздывающих и мгновенных нейтронов, следовательно, пренебрегается временным сдвигом в распределении предшественников запаздывающих нейтронов. Тогда слагаемое 26

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

28 описываются членом, который в квазистатическом приближении должен быть исключен, что сделать сравнительно просто в процессе расчета количества нейтронов следующего поколения. В то же время должен быть распределен источник запаздывающих нейтронов, который в уравнении (2.12) описывается членом. Помимо тех операций, которые должны быть выполнены для квазистатического приближения, из уравнений (2.11) и (2.12) видно, что для численного решения уравнения методом Монте-Карло (2.11) следует учесть при расчете коэффициента ослабления слагаемое. В то же время должны быть распределены источники нейтронов по распределению на предыдущем временном шаге, которые в уравнении (2.11) выражаются членом. Данный функционал оценивается на предыдущем временном интервале. 2.3 Алгоритмы моделирования кинетических процессов прямым методом Для моделирования динамических процессов, протекающих в ядерных реакторах с использованием метода Монте-Карло, необходимо многократно решить неоднородное уравнение переноса нейтронов. Более того, можно сказать, что расчёт кинетики методом Монте-Карло это прямое моделирование физических процессов, протекающих в реакторах нулевой мощности при заданном внешнем источнике нейтронов. Внешний источник нейтронов формируется из двух составляющих: ИНН (источник начальных нейтронов) нейтроны находящиеся в реакторе на начало моделируемого процесса. ИЗНН (источник запаздывающих начальных нейтронов) нейтроны появляющиеся от предшественников запаздывающих нейтронов, накопленных в реакторе на начало моделируемого процесса. После того как источники сформированы одна за другой прослеживаются истории каждого нейтрона и его потомков до момента исчезновения из системы (захват, утечка или выход за временной интервал моделирования динамического процесса [0, T] с). В случае возникновения в процессе деления от этих нейтронов или от их потомков 28

29 вторичных нейтронов, все они записываются в банк нейтронов и их истории также будут прослежены. Если в процессе расчета появляется запаздывающий нейтрон, то по известным выходам запаздывающих нейтронов β i разыгрывается группа эмиттера и сразу разыгрывается время появления запаздывающего нейтрона, по известной постоянной распада λ i. После чего он записывается в банк данных и в случае если время его появления лежит в рассчитываемом временном интервале его история также будет прослежена. Вес каждого нейтрона равен 1, нейтрон либо появился, либо нет. Полный поток нейтронов в реакторе представляет собой сумму потоков, образующихся от ИНН и ИЗНН. Таким образом, процесс моделирования кинетического процесса сводится к расчёту историй нейтронов от ИНН и ИЗНН, а так же всех их потомков. При этом, если поток и другие функционалы от нейтронов ИНН и ИЗНН и их потомков регистрируются отдельно, то нет необходимости в процессе расчёта поддерживать соотношение между числом нейтронов ИНН и ИЗНН. В этом случае полный поток в реакторе можно определить как линейную комбинацию потоков от ИНН и ИЗНН: Ф(r,Ω,E,t)=Ф ИНН (r,ω,e,t)+αф ИЗНН( r,ω,e,t) Коэффициент α представляет собой следующее соотношение: Числитель данного соотношения можно выразить следующим образом: где T длительность рассчитываемого временного интервала. Множитель учитывает то, что только часть эммитеров распадется на рассчитываемом временном интервале. 29

30 Т.к. расчет любого процесса начинается с критического состояния, то для С i, можно воспользоваться формулой (2.1). После подставления ее в выражение для С i, можно получить следующую формулу для коэффициента α Данная формула получена для случая, если в системе находится только один делящийся нуклид. В случае если в системе существует n делящихся нуклидов, то формула для коэффициента α будет выглядеть следующим образом: где P j средняя по реактору вероятность деления на j-м нуклиде, эффективная доля запаздывающих нейтронов i-й группы, постоянная распада i-й группы, j-го нуклида. Полученная формула для α справедлива, если рассчитанные Ф ИНН и Ф ИЗНН нормированы на один нейтрон источников ИНН и ИЗНН соответственно. В противном случае в формуле (1) нужно учитывать количество промоделированных нейтронов источников ИНН и ИЗНН. Т.к параметры P j, и необходимые для расчета коэффициента α определяются со статистической погрешностью, в программе КИР также реализован другой алгоритм расчета α, согласно которому, сначала проводится расчёт временного интервала [0, t 1 ] c, в течение которого реактор находится в критическом состоянии, а непосредственное моделирование динамического процесса начинается с момента времени t 1. В интервале [0, t 1 ] c поток нейтронов не зависит от времени, учитывая, что этот поток равняется сумме потоков от ИНН и ИЗНН, можно определить коэффициент α. Разбивая временной интервал [0, t 1 ] на N равновеликих подинтервалов, в каждом из которых будем регистрировать поток от ИНН и ИЗНН. После чего подбирается коэффициент α исходя из независимости от времени результирующего потока. Очевидно, что 30

31 где Ф 0, t ] — полный поток в интервале [0, t 1 ], [ i ИНН Фi и ИЗМН Ф i потоки в i-ом временном интервале от ИЗНН и ИНН соответственно. На рисунке 2.1 проиллюстрировано распределение по времени ИНН и ИЗНН (по экспонентам), на начало моделируемого процесса. Рисунок 2.1 носит только иллюстративный характер и не показывает численные соотношения между ИНН и ИЗНН. Рисунок Распределение ИЗНН и ИНН по времени Такой подход позволяет использовать любое число нейтронов как в ИНН, так и в ИЗНН, в зависимости от характеристик моделируемого кинетического процесса. При использовании прямого метода расчета точность моделирования и допустимая длительность рассматриваемого временного интервала протекания процесса определяется только точностью оценённых ядерных данных и имеющимися вычислительными ресурсами, т.е. точностью описания источника. 31

32 Реализованный в программе КИР алгоритм решения нестационарного уравнения переноса состоит из трех частей: 1. Расчёт стационарного критического состояния реактора для определения пространственного распределения скорости реакции рождения нейтронов; 2. Формирование ИНН и ИЗНН с учётом пространственного распределения скорости реакции рождения нейтронов; 3. Расчёт нестационарного кинетического процесса. Расчёт стационарного критического состояния (однородной задачи) реактора осуществляется стандартным методом Монте-Карло по поколениям и не нуждается в пояснениях. Для моделирования кинетического процесса и удобства анализа результатов временной интервал [0, T] разбивается на ряд интервалов i tmax, не обязательно равновеликих, в каждом их которых регистрируются заданные пользователем функционалы. В течение каждого временного интервала свойства реактора считаются постоянными. Введем следующие обозначения: банк первой очереди пакет нейтронов, находящихся в текущем расчётном временном интервале, банк второй очереди пакет потомков от ИНН находящихся за пределами текущего временного интервала, банк третьей очереди пакет потомков от ИЗНН. Для расчёта нестационарного кинетического процесса во временном интервале [0,T] с устанавливаются следующие значения: — NHISTP число нейтронов в ИНН; — NHISTD число нейтронов в ИЗНН; — NPROC число доступных процессоров; — IPROMPT параметр указывающий что нейтрон или мгновенный IPROMPT=1 или запаздывающий IPROMPT=0; — i tek = 1 текущий временной интервал; — i tmax полное число рассчитываемых интервалов. Исходя из данных, полученных при расчёте стационарного состояния, для каждого нейтрона источника модуль источников программы КИР генерирует: координаты точки рождения из пространственного распределения энерговыделения; 32

33 энергию из распределения скоростей генерации нейтронов деления по делящимся нуклидам и соответствующего спектра мгновенных нейтронов, если нейтрон мгновенный, и запаздывающего, если нейтрон запаздывающий; направляющие косинусы полёта из изотропного распределения; время рождения запаздывающего нейтрона исходя из равновесных концентраций предшественников запаздывающих нейтронов в соответствующей координатам точки рождения регистрационной зоне и нуклида во временном интервале [0, T] c моделирования кинетического процесса. Время рождения мгновенного нейтрона устанавливается равным нулю, это время рождения принимается за начальный возраст нейтрона. Весь алгоритм расчета кинетики реактора по программе КИР описывается следующими пунктами: 1) Генерируется NHISTP мгновенных нейтронов, которые помещаются в банк первой очереди. IPROMT устанавливается равным 1. Переход на 10). 2) Генерируется NHISTD запаздывающих нейтронов. Если время рождения запаздывающего нейтрона превышает первый временной интервал, то нейтрон помещается в банк третьей очереди, иначе в банк первой очереди. IPROMPT устанавливается равным нулю; переход на 10); 3) Разыгрывается длина свободного пробега нейтрона. 4) По длине свободного пробега и скорости нейтрона определяется время его жизни до точки столкновения, которое добавляется к суммарному возрасту данного нейтрона. 5) Если суммарный возраст нейтрона превышает верхнюю границу временного интервала моделирования кинетического процесса T, то переход на 10); 6) Делаются вклады в оценки по столкновениям по регистрационной временной и пространственной сетке. Если IPROMPT=1, вклады делаются для функционалов от ИНН, или от ИЗНН. 7) Стандартным способом моделируется тип столкновения; 8) В случае если происходит рассеяние продолжается прослеживание траектории переход на 3), если утечка или захват без деления нейтрона переход на 10), иначе переход на 9) 33

34 9) Разыгрывается число N f родившихся в результате деления нейтронов; определяется, сколько из них мгновенных и запаздывающих. Мгновенным нейтронам присваивается полный возраст нейтрона, вызвавшего деление. Для запаздывающих нейтронов разыгрывается время рождения (появление из предшественника). Запаздывающим нейтронам присваивается полный возраст нейтрона, вызвавшего деление, к которому добавляется время рождения запаздывающего нейтрона. Нейтроны, суммарный возраст которых превышает верхнюю границу временного интервала i tek при IPROMPT=1, помещаются в банк второй очереди, при IPROMPT=0 в банк третьей очереди. Нейтроны, суммарный возраст которых превышает верхнюю границу временного интервала моделирования кинетического процесса T, отбрасываются, остальные нейтроны помещаются в банк первой очереди; переход на 10); 10) Проверяется, пуст ли банк первой очереди. Если нет выбирается нейтрон из банка, переход на 3), иначе переход на 11); 11) Если IPROMPT=1 и i tek =1 переход на 2), иначе переход на 12); 12) Если IPROMPT=0, присваивается i tek = i tek + 1. Если i tek >i tmax, то переход на 15); 13) Если IPROMPT=1 присваивается IPROMPT=0, иначе присваивается IPROMPT=1; переход на 14); 14) Если IPROMPT=0, нейтроны из банка третьей очереди с суммарным временем жизни, лежащим в пределах временного интервала i tek, пересылается в банк первой очереди. Если IPROMPT=1 пересылается соответствующие нейтроны из банка второй очереди; переход на 10); 15) Окончание расчёта. Следует отметить, что реализованный в программе КИР алгоритм регистрации рассчитываемых функционалов позволяет использовать несколько временных сеток вплоть до максимально подробных, с временным шагом сравнимым со средним временем жизни нейтронов деления. Это позволяет переопределять, например, положение органов регулирования или температур и плотностей материалов не в определенные моменты времени, заданные как исходные данные, а непосредственно при моделировании кинетического процесса, используя некоторые критерии, такие как ограничения по изменению мощности реактора, допустимые максимальные изменения энерговыделения в твэлах или их частях, и т.п. 34

35 При наличии соответствующих вычислительных ресурсов можно организовать параллельный расчёт ИНН и ИЗНН, задав NHISTD равным значению С (см. формулу(**)), т.е. все предшественники распадаются с вероятностью 1), проследив эти запаздывающие нейтроны на своих процессорах. В этом случае можно точно сымитировать работу реактора-суперкомпьютера, к котором предшественники распадаются по своим законам вне зависимости от летящих в реакторе нейтронов. 2.4 Алгоритмы моделирования кинетических процессов в приближениях Моделирование кинетического процесса начинается с вычисления начального состояния, которое определяется из решения однородного уравнения переноса нейтронов. Результатами расчётов являются величина K эф, распределение энерговыделения, скорости реакции генерации нейтронов и параметры точечной кинетики. В программе КИР-П метод Монте-Карло применяется для расчёта параметров точечной кинетики и вычисления форм-функции. Для моделирования кинетических процессов в программе КИР-П используются адиабатическое, квазистатическое и усовершенствованное квазистатическое приближения. Для этих приближений зависимости от времени амплитудного множителя и концентраций предшественников запаздывающих нейтронов описываются уравнениями ( ), а форм-функции уравнениями ( ). Разные приближения могут использоваться как при решении разных задач кинетики быстрых реакторов, так и при моделировании одного кинетического процесса на разных временных отрезках его протекания. На выбор приближения должно влиять то, как сильно изменяются на данном временном отрезке мощность реактора, его геометрия, распределение энерговыделения и др. Плотность потока нейтронов представляется в виде произведения амплитудного множителя n(t) и форм-функции, т.е. 35

36 Схема решения приближённых уравнений кинетики традиционна и заключается в том, что для каждого очередного временного интервала вычисляются амплитудный множитель и концентрации предшественников, после чего пересчитывается формфункция. При расчёте форм-функции вычисляются и параметры точечной кинетики. В программе КИР-П имеется возможность рассчитывать, с учетом функции ценности, следующие параметры точечной кинетики: время жизни мгновенных нейтронов [29] и эффективная доля запаздывающих нейтронов [30]. Но на данном этапе развития программы, при расчете указанных параметров, на этапах расчета кинетики, функция ценности не использовалась. В КИР-П используется то же самое константное обеспечение, что и в программе КИР. Преимущество программы КИР-П перед программой КИР в более быстром, на несколько порядков, моделировании кинетического процесса. А преимущество перед инженерными программами, в которых реализованы те же приближения, заключается в прецизионном расчёте форм-функции и параметров точечной кинетики методом Монте- Карло. Общий принцип расчета форм-функции в квазистатическом приближении Необходимо решить уравнение: В виду принятого отсутствия взаимодействия нейтрон с нейтроном уравнение (2.14) можно решать в два этапа. На первом этапе рассчитать форм-функцию обусловленную источниками мгновенных нейтронов существующих в момент времени t, а на втором этапе обусловленную источниками запаздывающих нейтронов появившихся в системе на рассчитываемый момент времени. 36

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

38 рисунке 2.2. Алгоритм и схема расчета Схема расчёта кинетического процесса в приближениях представлена на Рисунок Схема расчёта, используя квазистатическое приближение 1) Для момента времени t 0 =0 рассчитывается стационарное состояние реактора, для которого вычисляются: коэффициент размножения К эф0, скорости генерации нейтронов по делящимся нуклидам и зонам, средние энергии нейтронов, вызывающих деление, для каждого нуклида; 2) Рассчитывается возмущенное состояние реактора с измененной геометрией или измененными физическими свойствами материалов, из которого получаем: Коэффициент размножения К эф1 (реактивность), β эфф, время жизни мгновенных нейтронов, распределение энерговыделения и другие необходимые функционалы. 3) Расчет амплитудного множителя n(t) и слагаемого необходимого для расчета в квазистатическом приближении коэффициента ослабления; 4) Если приближение адиабатическое то переход на 5), если квазистатическое, то производится расчет форм функции: 38

39 4.1 Производится расчет однородного уравнения (2.16); 4.2 Распределяются нейтроны по распределению эмиттеров запаздывающих нейтронов (по распределению генерации запаздывающих нейтронов); 4.3 Производится расчет неоднородного уравнения (2.17); 4.4 Оба расчета складываются с посчитанным коэффициентом. 5) Расчет результирующей ППН на интервале: ; 6) Далее, если конец рассматриваемого процесса не достигнут, переходим к пункту 2, в противном случае окончание расчета. Учет распределения запаздывающих нейтронов. Распределение эмиттеров запаздывающих нейтронов, для расчета уравнения (2.17) формируется следующим образом: Входящие в эту сумму распределения концентраций эмиттеров рассчитываются следующим образом: где рассчитываемый момент времени; — какой-либо из предыдущих моментов времени; — распределение эмиттеров запаздывающих нейтронов, посчитанное для какого-либо из предыдущих моментов времени; постоянная распада для i-1 группы запаздывания; выход запаздывающих нейтронов для i-1 группы запаздывания. 39

40 2.5 Алгоритмы модуля источников Модуль источников ИСТ используется при решении двух разных задач. Первая задача это расчёт критичности при решении однородного уравнения переноса нейтронов, а вторая задача это расчёт кинетики, который является составной частью моделирования динамического процесса, при котором многократно решается неоднородное уравнение. При моделировании динамики однородное уравнение решается на начальном этапе. Для решения однородного уравнения методом Монте- Карло необходимо определить фазовые координаты нейтронов первого поколения. При решении обеих задач принимается, что угловое распределение нейтронов источника всегда изотропно в лабораторной системе координат. При решении первой задачи энергия нейтрона E определяется из спектра деления Уатта с константами из библиотеки BNAB/MCU. При моделировании динамики определяются вероятности генерации нейтронов при делении ядер изотопов U 235, U 238, Pu 239, Pu 240, U 233, U 236, U 232, U 234, Pu 241, Pu 242, Pu 238, Pu 243 и средние энергии нейтронов, вызывающих деление этих ядер. Энергия мгновенного нейтрона разыгрывается из спектра Уатта с вероятностями: для деления ядра U 235 нейтроном с энергией, для деления ядра U 238 нейтроном с энергией, для деления ядра Pu 239 нейтроном с энергией, для деления ядра Pu 240 нейтроном с энергией, для деления ядра Pu 243 нейтроном с энергией =1. Параметры и пересчитываются на каждом временном шаге. Для розыгрыша энергии нейтрона сначала определяется, какой изотоп претерпел деление. Далее для этого изотопа из спектра Уатта определяются энергии мгновенных нейтронов. Для определения энергии запаздывающего нейтрона используются спектры из библиотеки BNAB/MCU. 40

41 При решении однородного уравнения принимается, что время рождения мгновенного нейтрона t=0. Алгоритм розыгрыша времени рождения запаздывающих нейтронов описан ниже. Для вычисления пространственных координат точки рождения нейтрона при решении однородного уравнения имеются три возможности: — можно задать «точечный источник», и тогда все нейтроны рождаются в одной точке. По умолчанию считается, что все нейтроны источника рождаются в точке с координатами (0.,0.,0.); — имеется опция, позволяющая разыгрывать координаты нейтронов источника равномерно по объёму тел из следующего списка: прямой круговой цилиндр с образующей, параллельной оси OZ, параллелепипед с рёбрами, параллельными осям координат; правильная шестиугольная призма с образующей, параллельной оси OZ, сфера; — имеется опция, позволяющая разыгрывать координаты нейтронов источника по вероятностям генерации нейтронов в выбранных зонах, вероятности должны быть посчитаны ранее, при каждом расчете они записываются в файл. При моделировании кинетического процесса, вычисляя пространственные координаты точки рождения нейтрона, предполагается, что элементы топлива (ТВС, твэлы) располагаются по правильной треугольной или квадратной решётке и нейтроны источника равномерно распределены по объёмам этих элементов (зон). Вероятности генерации нейтрона в зонах определяются для каждого временного шага. При моделировании кинетического процесса отдельно учитываются вклады в рассчитываемые функционалы от мгновенных нейтронов, имеющихся в реакторе перед началом кинетического процесса, и запаздывающих нейтронов, предшественники которых накопились к этому моменту. По определённым вероятностям генерации нейтронов при делении ядер определяется изотоп, в результате деления ядра которого образуются предшественники запаздывающих нейтронов. При работе реактора на постоянной мощности на рассматриваемом временном интервале от 0 до Т в каждый момент t * при делении ядра выбранного изотопа рождается β i запаздывающих нейтронов группы i, т.к. плотность вероятности распада предшественников запаздывающих нейтронов группы i есть, то 41

42 вероятность того, что эти предшественники доживут до конца интервала Т равна В результате всех делений на интервале от 0 до Т доля эмиттеров запаздывающих нейтронов группы i будет равна Обычно предполагается, что перед началом кинетического процесса реактор работает на мощности достаточно долгое время (T ), тогда. Для определения группы запаздывающих нейтронов, образующихся на интервале времени от t 1 до t 2 после начала динамического процесса, доля умножается на вероятность того, что эти предшественники доживут до момента времени t 1 и на вероятность их распада на этом интервале. В результате чего итоговая вероятность P i распада предшественников запаздывающих нейтронов группы i на интервале от t 1 до t 2 будет равна: Для выбранного изотопа и полученной группы i определяется энергия запаздывающего нейтрона и время его рождения. Формулы для получения выборочного значения времени распада t на интервале времени от t 1 до t 2 следующие: где ξ случайное число, равномерно распределённое на [0, 1]. 42,

43 Глава 3 Описание программы 3.1 Назначение и структура программы Программа КИР предназначена для решения неоднородного стационарного и нестационарного и однородного уравнений переноса нейтронов аналоговыми методами Монте-Карло на основе оценённых ядерных данных в системах с трёхмерной геометрией на одноядерных и многоядерных компьютерах. В программе КИР реализовано решение нестационарного уравнение переноса нейтронов методом Монте- Карло с временными зависимостями положения органов регулирования, плотностей и температур материалов конструкционных элементов. Учитываются запаздывающие нейтроны как от предшественников, которые накопились к началу моделируемого временного процесса, так и предшественников, которые образуются в течение процесса. Программа КИР-П предназначена для расчётов кинетики ядерных реакторов с использованием метода факторизации, когда поток нейтронов представляется в виде произведения амплитудного множителя и форм-функции, в рамках классических адиабатического, квазистатического и усовершенствованного квазистатического приближений [22, 28]. Форм-функция рассчитывается методом Монте-Карло по алгоритмам, аналогичным тем, которые реализованы в программе КИР для расчётов потока нейтронов, входящие в уравнения перечисленных приближений эффективная доля запаздывающих нейтронов и время жизни мгновенных нейтронов также рассчитываются методом Монте-Карло с учётом функции ценности нейтронов деления. За основу разработки программ КИР и КИР-П были взяты модули ППП MCU-4 [51] и MCU-5 [52], из которых собираются рабочие программы, дающие возможность решать только однородное и стационарное неоднородное уравнения переноса нейтронов. Нестационарное уравнение эти программы решать не позволяют. Так же программа КИР имеет модульную структуру; номенклатурный состав модулей совпадает с составом, используемым в программах, собираемых из модулей ППП MCU-5. — под модулем [51, 52] имеется в виду одна или несколько программных единиц с предписанным архитектурой пакета предназначением и межмодульными связями, которые регламентируются интерфейсами с полностью описанными точками входа и выхода. 43

44 — под специализированной рабочей программой исполняемая программа, собранная из вполне определенных модулей и предназначенная для решения конкретных задач математического моделирования физических процессов. В этих определениях КИР это специализированная рабочая программа, которая собираются из модулей следующих типов. — Управляющий модуль (УМ) организует совместную работу всех модулей. — Модуль траектории (МТ) осуществляет организацию моделирования отдельных траекторий нейтронов в пространстве и времени. — Модуль источников (МИ) моделирует фазовые координаты частиц источника или начального пакета нейтронов при решении однородной и неоднородной задачи. — Геометрический модуль (ГМ) определяет местонахождения частицы, отслеживает её перемещение по пространству и вычисляет все необходимые функции, аргументами которых служат пространственные координаты. — Физический модуль (ФМ) вырабатывает необходимые сечения взаимодействия нейтронов со средой и моделирует их столкновения с ядрами среды. — Регистрационный модуль (РМ) накапливает статистику событий, необходимую для оценки функционалов, в том числе и зависящих от времени протекания кинетического процесса, а потом обрабатывает эту статистику. — Модуль оборудования (МО) включает программы, которые объединяются в подмодули, зависящие от типа компьютера и применяемой операционной системы: — датчик псевдослучайных чисел; — программы открытия файлов, необходимых для работы, и назначения их логических номеров; — программы датчика времени и даты; — программы ввода/вывода; — программы бесформатного ввода. Можно записать, что КИР = УМ+ МТ + МИ + ГМ + ФМ + РМ + МО. 44

45 Так же как и в программах семейства MCU, организацией счёта занимается модуль траекторий, последовательно вызывая соответствующие подпрограммы модулей источников, геометрического, физического и регистраций. Для программы КИР использовались тексты геометрического модуля NCG и модуля оборудования ENV из пакета MCU-5 без внесения изменений. Были модифицированы тексты модуля управления УМ и подмодули физического модуля СОФИЗМ, FARION, ФИМБРОЭН и ФИМТОЭН [52]. Выбранная компоновка физического модуля обеспечивает возможность непрерывного слежения за энергией нейтрона и определения изотопа, с которым взаимодействует нейтрон. Модули траекторий ТРАКТ и источников ИСТ написаны специально для программы КИР. Кроме того, модуль регистрации RGS дополнен подмодулем пользователя РУГА. Таким образом, компоновка программы КИР определяется следующей формулой: КИР = УМ + ТРАКТ + ИСТ + NCG + +(СОФИЗМ+FARION+ФИМБРОЭН+ФИМТОЭН+RAPAN+FIMUL) + + (RGS+РУГА) + ENV. Программу КИР-П составляют две программы КИР-АМ и КИР-ФФ: КИР-П = КИР-АМ + КИР-ФФ. Первая из них предназначена для расчёта амплитудного множителя, а вторая для расчёта форм-функции. Программа КИР-ФФ и по структуре, и по модульному наполнению, и по константному обеспечению полностью совпадает с программой КИР: КИР-ФФ КИР. В текущей версии программы КИР-П в качестве КИР-АМ используется подпрограмма TFUNK из программы PRISET-MBIR [34]. 45

46 3.2 Описание геометрии расчетной системы Программа позволяет моделировать системы, состоящие из объёмных элементов практически произвольной формы. Геометрический модуль универсального типа позволяет моделировать трёхмерные системы с произвольной геометрией, используя комбинаторный подход, основанный на описании сложных пространственных форм (зоны) комбинациями простых тел или поверхностей, с помощью теоретикомножественных операций пересечения, дополнения и объединения. Существует некоторый набор типов тел-примитивов, для каждого типа тела имеется система параметров, полностью описывающая форму конкретного тела и его положение в пространстве. Каждому типу тела соответствует программа, определяющая точки пересечения любой прямой с любым конкретным телом этого типа. Для обеспечения эффективности счёта все тела ограничены плоскостями или поверхностями второго порядка. Возможно также задание решёток, получаемых размножением некоторых исходных элементов, заданных с помощью комбинаторики. 3.3 Моделирование взаимодействий нейтронов с веществом Константное обеспечение программы составляет банк данных MCUDB50. Программа позволяет рассчитывать нейтронно-физические процессы с учетом непрерывного изменения энергии частицы при столкновениях, а также имеется возможность использовать ступенчатую зависимость сечений от энергии. При моделировании переноса нейтронов учитываются следующие эффекты: — при генерации нейтронов деления допускается использование спектра деления мгновенных и запаздывающих нейтронов; — в быстрой энергетической области учитывается анизотропия упругого рассеяния в системе центра масс, имеется возможность проводить моделирование неупругих столкновений с учётом законов, содержащихся в файлах оценённых ядерных данных. — в области неразрешённых резонансов, сечения вычисляются по подгрупповым параметрам или с использованием f — факторов Бондаренко, в обоих случаях с учётом температурной зависимости используемых параметров. — в области разрешённых резонансов допускается как подгрупповое, так и поточечное описание сечений. Для наиболее важных изотопов используются 46

47 непрерывные зависимости сечений, так как при моделировании в каждой энергетической точке они вычисляются по резонансным параметрам. Такая схема позволяет проводить расчёты непосредственно с использованием данных по резонансным параметрам без предварительной подготовки таблиц сечений и оценивать температурные эффекты через аналитические зависимости сечений от температуры. Моделирование столкновений в области термализации проводится по модели непрерывного изменения энергии с учётом корреляций между изменением энергии и угла при рассеянии, при этом учитываются химические связи, тепловое движение ядер и когерентные эффекты для упругого рассеяния. 3.4 Условия применения программы Программа написана на языке Фортран 90/95. Для компиляции итогового текста необходимо использовать компилятор Intel (Intel Visual Fortran Compiler). Для многопроцессорного режима счета программа поставляется с настройками библиотеки MPI, в которой используется пакет MPICH2. Распараллеливание программы осуществлено на базе программного интерфейса MPI (Message Passing Interface), являющимся наиболее распространённым стандартом интерфейса обмена данными в параллельном программировании, и его реализации существуют для большого числа компьютерных платформ. При работе в режиме многопроцессорных вычислений программа задействует для расчёта все доступные ей процессоры (точнее ядра, далее в тексте процессоры предполагаются одноядерными). Коэффициент распараллеливания программы при отсутствии промежуточных записей приблизительно равен 1. Общая схема расчёта при этом остается такой же, как и при расчёте на одном процессоре. Основным процессором является нулевой процессор. Помимо собственно счёта он контролирует прохождение задачи. Таким образом, однопроцессорный расчёт с точки зрения программы — это расчёт, выполняемый только на нулевом процессоре. Сбор информации, накопленной процессорами, осуществляется в моменты промежуточной и финальной записи на диск. В момент промежуточной записи происходит наибольшая потеря скорости вычислений, поскольку возникают затраты времени на синхронизацию процессоров (ожидание самого медленного) и на собственно 47

48 запись на диск. Для расчетов необходимо выбирать максимально возможный, но приемлемый с точки зрения аварийного прерывания счёта интервал для сохранения промежуточных результатов на диск. 3.5 Изменение параметров реактора при моделировании динамического процесса Для моделирования динамического процесса был разработан алгоритм автоматического изменения характеристик ЯР. Каждый из параметров, описывающих динамический процесс (положение стержней, плотность теплоносителя, температуры материалов и т.п.), могут быть изменены перед началом расчета следующего временного шага. Реализованный алгоритм использует возможность задания части данных, зависящих от времени протекания динамического процесса, в специальных файлах, которые подключаются к основному файлу с описанием расчетной модели с помощью оператора INCLUDE. Перед каждым шагом INPUT программа автоматически выполняет в директории с файлом исходных данных поиск файла с именем вида DYNinnnn и при наличии такого файла преобразует его в файл с именем вида DYNi. Здесь i число от 0 до 9, а nnnn четырёхзначное число от 0 до 9999 с нулями в пустых позициях слева, равное номеру шага динамического процесса. Как написано в разделе 3.4 расчет одного временного интервала происходит в два шага, в связи с чем для одного временного интервала необходимо задать два одинаковых файла с именами DYNinnnn и DYNi[nnnn+1]. Необходимо отметить, что нумерация шагов начинается с нуля, т.о. перед самым первым шагом input программа выполняет поиск всех файлов вида DYNi0000 и, если такие файлы имеются, переименовывает их в файлы DYNi. Такого вида автоматизация в совокупности с использованием стандартной для исходных данных строки #INCNLUDE позволяет изменять параметры расчётной модели в процессе расчёта динамического процесса в зависимости от номера шага. 48

49 Описанные возможности могут быть использованы, например, для того, чтобы автоматизировать процесс изменения положения органов регулирования (поглощающих стержней) системы управления и защиты реактора и других параметров модели в процессе собственно расчёта, не изменяя данные исходного базового варианта. Для этого в файл исходных данных с помощью строки #INLCUDE в необходимые места добавляются ссылки на файлы вида DYNi, в соответствующих им файлах вида DYNinnnn определяются изменяемые параметры для соответствующего шага по времени nnnn, при этом файлы для нулевого шага должны быть обязательно. Файлы для остальных шагов могут появляться только в момент действительного изменения параметров. 49

50 Глава 4 Тестирование программы В настоящем разделе представлены результаты тестирования программного комплекса КИР. Рассмотрены международные расчётные бенчмарки, которые используются для верификации как диффузионных программ, так и основанных на методе Монте-Карло прецизионных программ, предназначенных для моделирования кинетики ядерных реакторов. 4.1 Тест RP1GC критический прямоугольный параллелепипед В работах [3 8] рассматривается модельная почти критическая с учётом источников запаздывающих нейтронов задача (К эф =1,00012±0,00005 или 0,99990±0,00001) с заданными одногрупповыми константами, в этих работах тест никак не идентифицирован. Для удобства дальнейшего изложения будем называть этот вариант тестом RP1G (Rectangular Prism with 1-Group Constants). Расчётная область представляет собой прямоугольный параллелепипед размерами 10х20х24 см 3, помещённый в вакуум. Одногрупповые сечения расчётной области представлены в Таблице 4.1. Скорость нейтронов v=22000 см/с, число вторичных нейтронов на один акт деления ν=2,5. Рассматривается 6 групп эмиттеров запаздывающих нейтронов, характеристики которых представлены в Таблице 4.2. Доля запаздывающих нейтронов составляет β=0, Таблица 4.1 Одногрупповые сечения расчётной области Сечение см -1 Σ tot 1,0 Σ fis 0,25 Σ s 0,4118 Σ abs =(Σ tot — Σ s ) 0,

51 Таблица 4.2 Характеристики групп запаздывающих нейтронов группы Постоянная распада, с -1 Относительный выход 1 0,0127 0, ,0317 0, ,1156 0, ,311 0, ,4 0, ,87 0,00178 Среднее/суммарное 0,0784 0, Моделирование кинетического процесса в критической системе В критической системе поток нейтронов не зависит от времени. Также не зависят от времени концентрации предшественников запаздывающих нейтронов, а скорость их распада равна скорости образования новых. Соотношение между мгновенными и запаздывающими нейтронами, образующимися в расчётной области, также постоянно. Таким образом, расчётное моделирование кинетического процесса в критической системе является достаточно хорошим тестом для первичной верификации программы. Все расчётные функционалы в пределах своих статистических погрешностей не должны зависеть от времени при моделировании кинетического процесса. В общем случае, для моделирования кинетического процесса на временном отрезке [0 — T] необходимо сформировать начальный пакет нейтронов, «влетающих» в него в момент времени 0 с и начальный пакет нейтронов (ИНН), образующихся в этом интервале, за счёт накопленных на момент времени 0 с предшественников запаздывающих нейтронов (ИНЗН). Также следует учесть, что часть нейтронов, «влетающих» во временной интервал в момент времени 0 с родились в результате распада предшественников, а характеристики запаздывающих нейтронов (например, их энергия) отличаются от мгновенных. Таким образом, если характеристики нейтронов обоих пакетов и соотношение между ними выбрано правильно, то при расчёте критической системы должен быть получен постоянный нейтронный поток на любой момент времени интервала [0 — T]. 51

52 Следует отметить, что доля мгновенных нейтронов вычисляется в процессе стационарного расчёта, в частности, зависит от таких характеристик, как средняя скорость нейтрона, среднего числа вторичных нейтронов деления и сечения деления. Как и любой функционал, рассчитанный методом Монте-Карло, эти величины имеют расчётную погрешность, следовательно, и соотношение между пакетами также будет определено с некоторой погрешностью. За счёт этой погрешности, т.е. приближённого определения соотношения между начальными пакетами, при моделировании нейтронный поток будет меняться в зависимости от времени протекания кинетического процесса. Но в конечном итоге через некоторое время нейтронный поток стабилизируется и становится не зависящим от времени — это следует из нейтроннофизических характеристик расчётной области, которая является критической. Другими словами, в расчётной области устанавливается равновесное соотношение между мгновенными и запаздывающими нейтронами, именно в этот момент и можно начинать непосредственное моделирование кинетического процесса. Время стабилизации нейтронного потока должно зависеть от статистической погрешности рассчитанных в стационарной задаче функционалов. Естественно эту погрешность можно уменьшить, увеличив статистику в расчёте стационарного состояния, тем не менее, она всегда будет присутствовать. Таким образом, перед началом моделирования необходим расчёт некоторого начального временного интервала, в течении которого нейтронный поток и другие функционалы становятся независящими от времени. Такой подход к моделированию кинетического процесса в целом аналогичен моделированию стационарного состояния методом Монте-Карло, при котором, как известно, отбрасывается некоторое количество начальных поколений (пакетов) нейтронов. Величина временного интервала, в течение которого нейтронный поток стабилизируется, зависит от точности определения соотношения между стартовыми пакетами. Интерес представляют крайние случаи, когда в качестве начального пакета используются либо ИНН, либо ИЗНН. 52

53 Без учёта ИЗНН расчётная область является подкритической и для такой системы нейтронный поток должен падать по экспоненциальному закону. На рисунке 4.1 приведён поток мгновенных нейтронов в зависимости от времени, рассчитанный по программе КИР, без учёта образования предшественников запаздывающих нейтронов. Регистрация потока осуществлялась во временных интервалах размером 0,001 с, нейтронный поток нормирован на единицу в первом временном регистрационном интервале. Как следует из представленных результатов, за 0,1 с нейтронный поток снижается примерно на четыре порядка. Рисунок 4.1 ППН без учета запаздывающих нейтронов и без учета ИЗНН 53

54 Учёт образующихся запаздывающих нейтронов, но без учета ИЗНН приводит к тому, что нейтронный поток на начальном этапе моделирования кинетического процесса снижается более плавно (см. рисунок 4.2), и по мере выхода предшественников запаздывающих нейтронов на равновесные концентрации в критической расчётной области становится постоянным, т.е. не зависящим от времени. При этом время выхода потока на не зависящий от времени уровень должно быть сравнимо с периодом полураспада предшественников запаздывающих нейтронов. Начиная примерно с 0,1 с основной вклад в поток нейтронов приходится на образующиеся запаздывающие нейтроны и их потомки (как запаздывающие, так и мгновенные), этот вклад в поток в зависимости от времени приведен на рисунке 4.3. Результаты расчёта нейтронного потока с учётом запаздывающих нейтронов от предшественников, образующихся в процессе расчёта на временном интервале с, приведены на рисунке 4.4. Регистрация осуществлялась во временных интервалах, равных одной секунде, результаты нормированы на единицу в первом временном регистрационном интервале. Для большей наглядности значение потока в этой временной зоне на рисунке 4.4 не приведено. Рисунок 4.2 ППН без учета ИЗНН, но с учётом образующихся запаздывающих нейтронов 54

55 Доля, % сек. Рисунок 4.3 Доля нейтронов и их потомков, образовавшихся от предшественников, которые образовались во время протекания кинетического процесса (t>0) в зависимости от времени, % Рисунок 4.4 ППН в зависимости от времени с учётом образования запаздывающих нейтронов, без учета ИЗНН Поскольку в расчётной схеме не учитывался ИЗНН, система является подкритической, поэтому на интервале от 0 до приблизительно 50 с абсолютное значение потока нейтронов падает. После t = 50 с источник нейтронов становится не 55

56 зависящим от времени (за счёт накопления и выхода на равновесные концентрации предшественников запаздывающих нейтронов), т.е. система выходит в критическое состояние. Источники (число рождений) 6-ти групп запаздывающих нейтронов в зависимости от времени представлены на рисунке 4.5. Как следует из полученных результатов, предшественники запаздывающих нейтронов, за исключением 1-ой группы (постоянная распада равна 0,0133 с -1 ) через 50 с выходят на равновесные концентрации. Доля нейтронов 1-ой группы составляет 3,6 % в общем источнике запаздывающих нейтронов. На момент времени 50 с источник запаздывающих нейтронов этой группы сформирован более чем на 90 %, поэтому влияние его роста на общий источник нейтронов становится незначительным. 1.E-02 Источник нейтронов, отн.ед. 1.E-03 1.E-04 1.E-05 1.E сек. Группа 1 Группа 2 Группа 3 Группа 4 Группа 5 Группа 6 Сумма Рисунок 4.5 Источник 6-ти групп запаздывающих нейтронов в зависимости от времени Если в качестве стартового пакета нейтронов использовать только ИНН, то требуется расчёт достаточно большого начального временного интервала (порядка с) для формирования равновесных концентраций предшественников запаздывающих нейтронов. Представленные выше результаты фактически показывают вклад ИНН в полный источник нейтронов во всем временном интервале. В критической системе источник 56

57 нейтронов постоянный и не зависит от времени, поэтому спад значения источника нейтронов, показанный на рисунках 4.2, 4.4, связан с отсутствием в расчётной схеме ИЗНН. Как следует из рисунка 4.4, этот вклад составляет около 5 %. В действительности он существенно меньше, поскольку представленные выше результаты не учитывают спад источника нейтронов во временных регистрационных интервалах, которые имели размер 1 с. Более детальная регистрация по времени показывает, что спад источника в интервале от 0 до 1 с превышает 2 порядка. Результат расчёта суммарного источника нейтронов (ИНН + ИЗНН) в начальном временном интервале в регистрационных зонах размером 0,01 с показан на рисунке 4.6. Рисунок 4.6 Суммарный источник нейтронов (ИНН + ИЗНН) в зависимости от времени (регистрационные зоны размером 0,01 с) Результат расчета интегральной ППН в зависимости от времени представлен на рисунке 4.7. Составляющие источников запаздывающих нейтронов представлены на рисунке 4.8. Временной интервал для регистрации источников равен одной секунде. 57

58 Рисунок 4.7 ППН в зависимости от времени (регистрационный интервал 1 с) Рассмотрим другой крайний случай, когда стартовый пакет нейтронов формируется исключительно из ИЗНН. Стартовый пакет нейтронов в этом случае формируется следующим образом. Сначала разыгрывается номер группы — i запаздывающих нейтронов, затем время t появления запаздывающего нейтрона этой группы в системе: где случайное число, равномерно распределённое на отрезке 1; постоянная распада i-ой группы запаздывающих нейтронов. Координаты рождения запаздывающего нейтрона в данном расчёте разыгрывались равномерно по содержащим делящиеся нуклиды зонам. Так же, как и при формировании стартового пакета, состоящего только из ИНН, требуется некоторое время, чтобы мгновенные нейтроны вышли на равновесные концентрации, и установилось не зависящее от времени соотношение между мгновенными и запаздывающими нейтронами, однако время этого процесса незначительно и составляет примерно 0,1 с. 58

59 Рисунок 4.8 ИЗНН (Источник t 0) и полный источник (сумма) Как следует из представленных выше результатов, в расчётной области сразу (в пределах статистической погрешности) формируется постоянный источник мгновенных и запаздывающих нейтронов, а, следовательно, и нейтронный поток. Регистрационные временные интервалы размером в одну секунду не позволяют отследить его изменение. Поэтому на рисунке 4.9 представлен полный источник нейтронов на временном отрезке [0-0,1с] с более детальной регистрационной сеткой. Ширина регистрационного временного интервала 10-5 с, источник нормирован на единицу на момент времени 0,1 с. 59

60 Рисунок 4.9 ППН в зависимости от времени (регистрационный интервал 10-5 с) Как следует из рисунка 4.9., начиная с момента времени 0,05 с, полный источник нейтронов в расчётной области становится постоянным и не зависящим от времени. Это означает, что система вышла в критическое состояние, и что в ней установилось не зависящее от времени соотношение между рождающимися мгновенными и запаздывающими нейтронами. Отклонение ППН от единицы в интервале 0 0,05 с связано с отсутствием в расчётной схеме ИНН. Вклад этих нейтронов вместе с их потомками в общий источник к моменту времени t = 0,05 с становится практически равным нулю. Для дальнейшего сокращения начального временного интервала, в течение которого нейтронный поток стабилизируется необходимо, как было отмечено выше, формировать в качестве начальных данных два пакета нейтронов с соответствующим соотношением между ними. Вопрос, имеет ли это смысл, на сегодняшний момент является открытым, и требуется некоторый объём практических расчётных исследований, чтобы на него ответить однозначно. В любом случае, представленные в настоящем разделе результаты расчётов по программе КИР подтверждают правильность работы реализованных алгоритмов моделирования кинетических процессов. 60

61 4.1.2 Одногрупповой кинетический тест 1 Для модельной призмы RP1G рассматривается следующий кинетический процесс. В интервале [0 10] c расчётная область находится в критическом состоянии. На 10-ой секунде сечение поглощения (Σ abs ) мгновенно уменьшается с 0,5882 см -1 до 0,5870 см -1 с соответствующим уменьшением Σ tot. На 40-ой секунде кинетического процесса оно также мгновенно становится равным первоначальному значению. Данная задача фактически моделирует процесс перевода реактора с одного уровня мощности на другой. Результаты расчёта по комплексу КИР данного процесса приведены на рисунке Регистрация потока осуществлялась во временных интервалах равных 0,1 с. Для сравнения также представлены результаты расчёта, приведённые в работе [4]. Следует отметить, что в работе [4] результаты расчёта данной задачи представлены в виде графика, поэтому более детальное сравнение не представляется возможным. Тем не менее, как следует из рисунка 4.10 данные, полученные по программе КИР, достаточно хорошо согласуется с результатами работы [4]. На рисунках 4.11 и 4.12 приведены изменения мощности в районе времени изменения нейтронно-физических свойств расчётной области. Скачки мощности после 10 и 40 с кинетического процесса связаны соответственно с увеличением и уменьшением потока мгновенных нейтронов, дальнейшие изменения мощности связаны соответственно с увеличением (после 10-ой секунды) и уменьшением (после 40-ой) концентраций предшественников запаздывающих нейтронов. Вклад в мощность запаздывающих нейтронов показан на рисунке

62 Рисунок 4.10 Зависимость мощности от времени кинетического процесса Рисунок 4.11 Зависимость мощности от времени кинетического процесса при уменьшении Σ abs 62

63 Рисунок 4.12 Зависимость мощности от времени кинетического процесса при увеличении Σ abs Рисунок 4.13 Мощность, выделяемая в расчётной области за счёт запаздывающих нейтронов в зависимости от времени 63

64 4.2 Тест RP1SC подкритический прямоугольный параллелепипед В этом тесте, предложенном в работе [3], рассматривается помещённый в вакуум подкритический прямоугольный параллелепипед, заполненный тем же материалом, что и в предыдущем тесте. Размеры призмы составляют 8х18х22 см 3. Моделирование начинается с решения условно-критической задачи, после чего процесс прослеживается в течение 1000 с. По условиям тестовой задачи предполагается, что источники запаздывающих нейтронов находятся в равновесных концентрациях на начало моделирования кинетического процесса. В работе [3] приведено значение эффективного коэффициента размножения нейтронов, который составляет К эф =0,97823±0, В расчёте по программе КИР эффективный коэффициент размножения составил К эф =0,9790±0,0005. Оценка нейтронного потока при моделировании кинетического процесса осуществлялась на временных интервалах, которые экспоненциально возрастают, начиная с интервала длительностью 10-4 с, что примерно соответствует времени жизни нейтронов в системе. В расчёте по программе КИР нейтронный поток регистрировался на двух временных сетках. Первая, детальная сетка, с шириной временного интервала 10-4 с использовалась в начале кинетического процесса (от 0 до 2 с). Вторая сетка с шириной временного интервала 0,1 с применялась на всем промежутке времени протекания процесса. Расчёты проводились на персональном компьютере. Прослеживались истории 1 млн. мгновенных нейтронов и запаздывающих, предшественники которых накопились к началу кинетического процесса. Статистическая погрешность интегральной плотности потока нейтронов без учёта корреляции последовательных поколений нейтронов

0,1 %. Учёт корреляции должен несколько увеличить статистическую погрешность. Результаты расчётов представлены на рисунке 4.14 в сравнении с данными, полученными по программе TRIPOLI-4.7 [3]. 64

65 Рисунок 4.14 Плотность потока нейтронов в подкритическом параллелепипеде в зависимости от времени протекания кинетического процесса Результаты, полученные по программам КИР и TRIPOLI, хорошо согласуются между собой, как в начале процесса, где преобладающий вклад в поток дают мгновенные нейтроны, так и в конце процесса, когда всё определяют запаздывающие нейтроны. Однако при более внимательном рассмотрении самого начала кинетического процесса (примерно до 0,001 с) можно заметить, что полученная в расчётах по программе КИР плотность потока нейтронов возрастает, что, очевидно, противоречит физике решаемой задачи (см. рисунок 4.15). Такое поведение плотности потока нейтронов объясняется тем, что в расчёте по программе КИР при формировании источника мгновенных нейтронов всем нейтронам присваивалось стартовое время 0 с. В действительности они имеют некоторое распределение по времени, в процессе моделирования движения мгновенных нейтронов и их потомков происходит их перераспределение по временной шкале. Примерно через 0,005 с они выходят на «реальные», соответствующие физическим свойствам расчётной области временные значения. Временной интервал, в котором происходит перераспределение, зависит от времени жизни нейтрона в расчётной области и составляет примерно несколько десятков времён жизни. В данном тесте этот интервал примерно равен [0 0,005] с. 65

66 .Рисунок 4.15 Плотность потока нейтронов в подкритическом параллелепипеде в интервале [0 0,01] c кинетического процесса Очевидно, что эта некорректность в начальном распределении нейтронов по временной шкале никак не сказывается на нейтронном потоке при моделировании длительных динамических процессов. При рассмотрении регистрационных временных интервалов размером 0,1 с начальный скачок потока просто нельзя обнаружить (см. рисунок 4.14). Однако, при моделировании динамических процессов, физические свойства в которых заметно меняются во временных интервалах, сравнимых со временем жизни нейтрона, неточности в определении начального распределения во времени мгновенных нейтронов может внести заметные погрешности в окончательный результат. Для корректного решения этой проблемы требуется проведение дополнительных исследований, выходящих за рамки данной работы. В частности, необходимо разработать и провести расчёты специальных математических бенчмарков с детальным учётом времени жизни мгновенных нейтронов на начало кинетического процесса. 66

67 4.3 Тест RPCEU235 прямоугольный параллелепипед с изотопом U 235 Тестовая задача RPCEU235 [3] представляет собой прямоугольный параллелепипед размерами 10х20х24 см 3. Параллелепипед заполнен изотопом U 235. В работе [3] приведена критическая концентрация U 235 равная 0,044925х /см 3. Критическая концентрация была получена по программе TRIPOLI-4.7 с использованием констант, подготовленных из файлов оценённых ядерных данных JEFF [54]. Расчёт по программе КИР при данной концентрации даёт значение коэффициента размножения 1,003±0,001. Для расчёта использовался подмодуль ФИМБРОЭН составного физического модуля с библиотекой BNAB/MCU в качестве константного обеспечения. Для приведения расчётной области к критическому состоянию концентрация U 235 была уменьшена до 0,044630х /см 3. Коэффициент размножения составил при этом 1,0001±0,0001. По условиям тестовой задачи во временном интервале с концентрация U 235 увеличивается на 0,000075х /см 3, а в интервале с снова принимает первоначальное значение. При моделировании кинетического процесса по программе КИР увеличение концентрации U 235 было скорректировано пропорционально уменьшению концентрации U 235 до критического состояния, т.е. концентрация U 235 в интервале с в расчёте по программе КИР была увеличена на 0, х /см 3. Расчёты проводились на суперкомпьютере на 240 ядрах. Статистическая погрешность без учёта корреляции нейтронов между поколениями меньше 0,01 %, что примерно соответствует 1 млрд. делений во временном интервале шириной 0,1 с. На рисунке представлена зависимость мощности от времени кинетического процесса. Для сравнения представлены результаты, полученные по программе TRIPOLI. Регистрация потока и мощности при расчётах по программе КИР осуществлялась во временных интервалах шириной 0,1 с. 67

68 1.E+03 Мощность, отн. ед. Мощность, отн.ед. 1.E+02 1.E+01 1.E сек. Время, с КИР TRIPOLI Рисунок 4.16 Зависимость мощности от времени кинетического процесса Результаты, полученные по программе КИР, достаточно хорошо согласуются с данными, представленными в [3]. Следует отметить, что в работе [3] результаты расчёта приведены в виде графика, поэтому более детальное сравнение не представляется возможным. Расчёты по программам КИР и TRIPOLI проводились с использованием библиотек, полученных на основе разных файлов оценённых ядерных данных, которые, в частности, отличаются характеристиками запаздывающих нейтронов. Кроме того, в работе [3] не указано с какой точностью было получено критическое состояние, а именно отсутствует значение коэффициента размножения. Время жизни нейтрона в данной задаче составляет порядка 10-8 с, поэтому в течение 70 секунд кинетического процесса сменяется порядка поколений. Небольшая неточность в определении критической концентрации U 235 может привести к значительным изменениям мощности в рассматриваемом временном интервале. На рисунке 4.17 приведены сравнительные результаты расчётов выделяемой мощности в параллелепипеде при концентрации урана 0,044630х /см 3 и 0,044640х /см 3. Эти расчёты носили оценочный характер и поэтому были проведены на небольшой статистике. 68

69 Мощность, отн.ед с А Б Рисунок 4.17 Зависимость мощности от времени кинетического процесса при концентрации U 235 равной 0,044630х /см 3 (А) и 0,044640х /см 3 Фактически на этом рисунке приведён дополнительный вклад в мощность при увеличенной на 0,02 % концентрации U 235 в расчётной области. На временном интервале в 60 с такое незначительное изменение концентрации приводит к увеличению мощности на

%. Сравнивать результаты, полученные по программам КИР и TRIPOLI, несмотря на общее хорошее согласие, можно только качественно, поскольку существует ряд неопределённостей в интерпретации исходных данных, которые касаются использования разных библиотек ядерных данных и точности определения критической концентрации U Тест BSS-6 одномерная модель реактора Тестовая задача BSS-6 (BenchmarkSourceSituation) [55] представляет собой одномерный двухгрупповой диффузионный тест. Расчётная область состоит из трёх размножающих зон, размеры которых представлены на рисунке Рисунок 4.18 Геометрия тестовой задачи BSS-6 В стационарном случае физические свойства зон 1 и 3 эквивалентны, на обеих границах по условиям тестовой задачи заданы нулевые потоки в обеих группах. 69

70 Спектр деления для мгновенных и запаздывающих нейтронов одинаков для всей расчётной области и равен χ 1 =1 и χ 2 =0. Скорость движения нейтронов первой группы равна см/c, второй см/c. В Таблице 4.3 приведены диффузионные двухгрупповые макроконстанты тестовой задачи, здесь под Σ а понимается сечение увода, т. е. Σ а1 =Σ с1 +Σ f1 +Σ s 1->2 Σ а2 =Σ с2 +Σ f2. Таблица 4.3 Начальные диффузионные двухгрупповые константы задачи BSS-6 Зона Σ а1, см -1 Σ а2, см -1 νσ f1, см -1 νσ f2, см -1 Σ s, см -1 Σ s, см -1 D 1, см 1 и 3 0,011 0,180 0,010 0,200 0,015 0,0 1,5 0,5 2 0,010 0,080 0,005 0,099 0,010 0,0 1,0 0,5 Учитывая, что рассеяние в диффузионном приближении принимается изотропным, полное сечение Σ tot можно определить через коэффициент диффузии D как Σ tot =1/(3*D), а внутригрупповые сечения рассеяния в первой и второй группах как Σ s1 =Σ tot1 -Σ a1 -Σ s 1->2 и Σ s2 =Σ tot2 -Σ a2, соответственно. Приняв ν=2.5 во всех группах и зонах, можно определить сечения захвата (Σ s ) и деления (Σ f ). Полученные транспортные макросечения тестовой задачи приведены в Таблице 4.4. Таблица 4.4 Начальные транспортные двухгрупповые константы задачи BSS-6 Зона Σ с1, см -1 Σ с2, см -1 Σ f1, см -1 Σ f2, см -1 Σ s, см -1 Σ s1, см -1 Σ s2, см -1 1 и 3 0,007 0,100 0,004 0,0800 0,015 0,1962 0, ,5 2 0,008 0,0404 0,002 0,0396 0,010 0,3133 0, ,5 Значение коэффициента размножения начального состояния при использовании диффузионного приближения равно 0,90155 [55]. В работе [56], в которой представлены результаты, полученные по программе SUHAM-МПГ с использованием метода поверхностных гармоник, приводится значение К эф =0, По программе КИР с использованием макросечений, представленных в Таблице 4.4, расчётное значение коэффициента размножения составило 0,9095±0,0005. Завышенное значение коэффициента размножения связано с ненулевыми потоками на границе рассматриваемой области при расчётах методом Монте-Карло, в которых ставятся условия вылета нейтронов на границе системы. Поэтому для получения с использованием транспортной системы констант коэффициента размножения, соответствующего диффузионной задаче, было скорректировано сечение 70 D 2, см ν и

71 внутригруппового рассеяния. Во всех группах и всех зонах оно было уменьшено на 0,04 см -1. Такой подход к сведению транспортных констант к диффузионным (и наоборот), хорошо известен [56, 57]. Значение коэффициента размножения с откорректированным внутригрупповым сечением рассеяния, полученное по программе КИР, составило 0,9019±0,0005. Пространственные распределения плотности потока нейтронов первой и второй группы, полученные по программе КИР, представлены на рисунках 4.19 и В расчёте по программе КИР размеры регистрационных зон равны 1 см, для сравнения также приведены результаты расчётов по программе SUHAM-МПГ, взятые из работы [56]. Расчёты проводились на суперкомпьютере на 240 ядрах. Статистическая погрешность без учёта корреляции нейтронов между поколениями

0,01 %, что примерно соответствует 100 млн. делений во временном интервале шириной 0,1 с. Как следует из рисунков, результаты с приемлемой точностью совпадают. Различия связаны, главным образом, с «краевыми эффектами» на границе расчётной области и на границах физических зон, что достаточно очевидно при сравнении диффузионных и транспортных расчётов. 71

72 Рисунок 4.19 Пространственное распределение плотности потока нейтронов первой группы в стационарном состоянии в сравнении с результатами работы [56] (Ref) В работе [55] рассматривалось несколько вариантов ввода реактивности, из которых были выбраны задачи A1 и A2, отличающиеся разными знаками вводимой реактивности. Задача А1 сечение поглощения нейтронов в первой зоне во второй энергетической группе линейно увеличивается на 3 % в течение первой секунды. Задача А2 сечение поглощения нейтронов в первой зоне во второй энергетической группе линейно уменьшается на 1 % в течение первой секунды. 72

73 Рисунок 4.20 Пространственное распределение плотности потока нейтронов второй группы в стационарном состоянии в сравнении с результатами работы [56] (Ref) По условиям теста начальное состояние приводилось к критическому путём деления величины νσ f на K эф. В расчётах по программе КИР критическое состояние было получено делением ν на K эф. Зависимость полной мощности в расчётной области от времени кинетического процесса для задачи А1 приведена на рисунке На рисунке приведены изменения мощности в каждой физической зоне. Представленные на рисунках зависимости мощностей нормированы на 1 в начале кинетического процесса, абсолютные значения мощности в сравнении с результатами работы [55] приведены в Таблицах 4.5 и 4.6. Регистрация функционалов в расчёте по программе КИР осуществлялась на временных интервалах длительностью 0,1 с. 73

74 Рисунок 4.21 Относительное изменение мощности в тестовой задаче А1 в зависимости от времени кинетического процесса в сравнении с результатами работы [11] (Ref) Мощность, отн.ед с Зона 1 Зона 2 Зона 3 Зона 1 (Ref) Зона 2 (Ref) Зона 3 (Ref) Рисунок 4.22 Относительное изменение мощности в физических зонах тестовой задачи А1 в зависимости от времени кинетического процесса в сравнении с результатами работы [55] (Ref) 74

75 Таблица 4.5 Относительная мощность в зависимости от времени кинетического процесса задачи А1 Время, с Р, КИР, отн.ед. Р, Ref, отн.ед. (1- Ref/КИР),% 0, ,1 0,9638 0,9299 3,52 0,2 0,8939 0,8733 2,30 0,5 0,7685 0,7597 1,15 1,0 0,6571 0,6588-0,25 1,5 0,6364 0,6433-1,08 2,0 0,6224 0,6307-1,33 Таблица 4.6 Относительная мощность в физических зонах расчётной области в зависимости от времени кинетического процесса задачи А1 Зона 1 Время, с Р, КИР, отн.ед. Р, Ref, отн.ед. (1- Ref/КИР), % 0, ,1 0,9230 0,8621 6,60 0,2 0,7902 0,7521 4,82 0,5 0,5515 0,5336 3,25 1,0 0,3505 0,3452 1,52 1,5 0,3186 0,3235-1,53 2,0 0,3004 0,3066-2,05 Зона 2 Время, с Р, КИР, отн.ед. Р, Ref, отн.ед. (1- Ref/КИР),% 0, ,1 0,9666 0,9340 3,37 0,2 0,9000 0,8805 2,17 0,5 0,7816 0,7724 1,18 1,0 0,6740 0,6753-0,19 1,5 0,6527 0,6588-0,93 2,0 0,6380 0,6455-1,18 75

76 Продолжение таблицы 4.6 Зона 3 Время, с Р, КИР, отн.ед. Р, Ref, отн.ед. (1- Ref/КИР),% 0, ,1 0,9986 0,9910 0,76 0,2 0,9995 0,9831 1,64 0,5 0,9652 0,9655-0,03 1,0 0,9383 0,9463-0,85 1,5 0,9298 0,9383-0,92 2,0 0,9211 0,9312-1,10 Аналогичные результаты для задачи А2 представлены на рисунках 4.23 и В Таблицах 4.7 и 4.8 приведены отклонения результатов расчёта по программе КИР от данных теста [55]. Отметим, что в тестовых результатах для задачи А1 значения мощности представлены в интервале 0 2 с, а для задачи А2 этот интервал составляет 0 4 с. Рисунок 4.23 Относительное изменение мощности в тестовой задаче А2 в зависимости от времени кинетического процесса в сравнении с результатами работы [55] (Ref) 76

77 6 5 Мощность, отн.ед ,5 1 1,5 2 2,5 3 3,5 4 с Зона 1 Зона 2 Зона 3 Зона 1 (Ref) Зона 2 (Ref) Зона 3 (Ref) Рисунок 4.24 Относительное изменение мощности в физических зонах тестовой задачи А2 в зависимости от времени кинетического процесса в сравнении с результатами работы [55] (Ref) Таблица 4.7 Относительная мощность расчётной области в зависимости от времени кинетического процесса задачи А2 Время, с Р, КИР, отн.ед. Р, Ref, отн.ед. (1- Ref/КИР),% 0, ,1 1,0184 1,028-0,94 0,2 1,0552 1,063-0,74 0,5 1,1856 1,205-1,64 1,0 1,7038 1,74-2,13 1,5 1,9970 1,959 1,90 2,0 2,2300 2,166 2,87 3,0 2,6890 2,606 3,09 4,0 3,1947 3,108 2,71 77

78 Таблица 4.8 Относительная мощность в физических зонах расчётной области в зависимости от времени кинетического процесса задачи А2 Зона 1 Время, с Р, КИР, отн.ед. Р, Ref, отн.ед. (1- Ref/КИР),% 0, ,1 1,0750 1,0560 1,77 0,5 1,3550 1,3980-3,17 1 2,3450 2,4350-3,84 2 3,3100 3,2160 2,84 3 4,1650 4,0170 3,55 4 5,1000 4,9280 3,37 Зона 2 Время, с Р, КИР, отн.ед. Р, Ref, отн.ед. (1- Ref/ КИР),% 0, ,1 1,0184 1,0270-0,85 0,5 1,1766 1,1930-1,39 1 1,6753 1,7010-1,53 2 2,1849 2,1130 3,29 3 2,6248 2,5400 3,23 4 3,1275 3,0270 3,21 Зона 3 Время, с Р, КИР, отн.ед. Р, Ref, отн.ед. (1- Ref/КИР), % 0, ,1 1,0022 1,0040-0,18 0,5 1,0222 1,0290-0,67 1 1,1124 1,1070 0,48 2 1,2378 1,1990 3,14 3 1,3443 1,2980 3,45 4 1,4818 1,4160 4,44 78

79 Полученные результаты по программе КИР достаточно хорошо согласуются с данными, полученными по другим программам, которые представлены в международном сборнике математических бенчмарков [55]. Максимальное отклонение в оценке интегральной мощности не превышает 3,5 %, в физических зонах 6 %. Наибольшие различия наблюдается в зоне 1 для временного интервала [0 1] c в задаче А1, т. е. в процессе линейного увеличения сечения поглощения. В зонах 2 и 3 различия в мощности не превышают 4,5 %. Следует также отметить, что процесс увеличения мощности (задача А2) идёт несколько быстрее, чем в тестовых результатах [55]. 4.5 Расчет бесконечной решетки твэлов ВВЭР прямым и приближенными методами Был разработан тест представляющий собой бесконечную решетку полномасштабных твэлов ВВЭР, который рассчитан, используя прямой метод решения уравнения переноса нейтронов, адиабатическое приближение и квазистатическое приближение. Данный тест представляется интересным в связи с возможностью получения по высоте твэла значительных перекосов профиля энерговыделения. Расчетная модель представляет собой шестигранную призму с полной утечкой на торцах ячейки и трансляционной симметрией на боковых гранях ячейки (см. рисунок 4.25). Концентрации нуклидов во всех материалах представлены в таблице Рисунок 4.25 Поперечное сечение ячейки с твэлом реактора ВВЭР 79