Решение уравнений переноса нейтронов в трехмерной HEX-Z геометрии

Аннотация

В данной работе описаны численные методы, используемые в пакете «РЕАКТОР» для решения уравнения переноса в многогрупповом приближении в прямоугольной X-Y-Z геометрии, гексагональной HEX-Z геометрии и цилиндрической R-φ-Z геометрии. Это наиболее трудоемкая часть пакета, предназначенная для расчета нейтронно-физических характеристик активной зоны ядерного реактора (значений Κ_eff) и радиационных полей S_n - методом дискретных ординат с P_m - приближением сечений рассеяния. Отличие от предыдущей версии пакета состоит в использовании наиболее эффективных взвешенных схем для решения уравнения переноса и методов ускорения внутренних и внешних итераций. В качестве начального приближения потоков для расчета Κ_eff используется численное решение, полученное в диффузионном приближении. Для расчета защиты задается распределенный по активной зоне источник делений. В качестве примера представлены результаты трехмерных расчетов и радиационных полей реактора СВБР 75/100.

Abstract

The numerical methods used in program package “REACTOR” for solving the transport equation in multi-group approximation in rectangular X-Y-Z geometry, hexagonal HEX-Z geometry and cylindrical R-φ-Z geometry are described in this work. It is the most working time component of program package, intended for calculation of neutron - physical characteristics of reactor core (eigenvalue Κ_eff) and radiation fields by the S_n - method of Discrete Ordinates with P_m - approximation of scattering cross sections. The distinction from the previous version of program package consists in using the most effective weighted schemes for solving the transport equation and inner and outer iteration acceleration methods. To calculate the eigenvalue the numerical solution obtained in diffusion approximation is used for initial approximation of fluxes. For the shielding problems the distributed fission source in reactor core is used. For example the three-dimensional results of calculation the eigenvalue Κ_eff and radiation fields are given for SVBR-75/100 reactor plant.

3. S_n приближение метода дискретных ординат. 7

6. Ускорение внешних итераций по источнику деления. 15

7. Ускорение внешних итераций для подкритических систем. 18

Введение

В течение многих лет задачи переноса нейтронов по реакторной тематике в ИПМ РАН решались численно методом Дискретных ординат с использованием пакета «РЕАКТОР» в двумерной RZ геометрии [1,2] и в трехмерной прямоугольной X-Y-Z геометрии [3,4]. В настоящее время в связи с бурным развитием вычислительной техники численное решение уравнения переноса в трехмерной геометрии находит все большее применение при решении реальных задач. Наряду с классическими прямоугольными сетками в X-Y-Z геометрии стали использоваться криволинейные сетки в R-φ-Z геометрии и гексагональные сетки в HEX-Z геометрии [5], поскольку они максимально согласованы с конфигурацией расчетной области.

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

§ KIN1D – расчет во всех одномерных геометриях: плоской, цилиндрической, сферической;

§ KIN2D – расчет в двумерной R-Z геометрии для цилиндрической области;

§ KIN3D – расчет в трехмерных геометриях. В X-Y-Z геометрии расчетная область представляет собой прямой параллелепипед с неравномерной сеткой по переменным X, Y и Z. В HEX-Z геометрии расчетная область представляет собой прямую призму, в основании которой задана равномерная гексагональная сетка, с неравномерной сеткой по вертикальной оси Z. Предполагается, что в плоскости O-X-Y область имеет произвольную выпуклую форму. В R-φ-Z геометрии расчетная область представляет собой прямой цилиндр или его часть с неравномерной сеткой по переменным R, φ и Z.

Данная работа посвящена описанию численных методов, реализованных в последнее время в пакете «РЕАКТОР» с целью повышения эффективности численного решения реальных трехмерных задач. Не вдаваясь в детали геометрии каждого модуля в отдельности, опишем наиболее общую их часть на примере трехмерной геометрии. Имеется в виду наиболее трудоемкая задача решения многогрупповой системы уравнений переноса нейтронов и g-квантов S_n методом дискретных ординат с P_m приближением индикатрисы рассеяния.

Приводятся результаты реальных трехмерных расчетов для реактора СВБР 75/100 [6-8].

1. Постановка задачи.

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

Правая часть уравнения (1) для нейтронов группы () имеет следующий вид:

Правая часть уравнения (1) для g-квантов групп имеет следующий вид:

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

В (1) - (4) использованы следующие обозначения:

- единичный вектор в направлении полета нейтронов,

- плотность потока в точке в направлении в группе ,

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

- функция распределения внешних источников.

Система уравнений (1) – (5) описывает стационарное распределение нейтронов и g-квантов с учетом процесса деления, а также заданных внутренних источников. Если в правой части (2) член с делением отсутствует ( для всех ), система уравнений (1) – (5) описывает стационарное распределение нейтронов и g-квантов в зависимости от заданных источников. Задача на (однородная задача) описывается системой уравнений (1) – (2) с нулевыми внутренними источниками и множителем перед членом с делением.

2. P_m приближение источника рассеяния.

Предполагается, что макроскопическое сечение рассеяния задается в виде ряда по полиномам Лежандра степени (в приближении):

где - присоединенные функции Лежандра первого рода с нормировкой:

Предполагается, что искомая функция разлагается в ряд по присоединенным функциям Лежандра:

где угловые моменты по косинусам и синусам решения уравнения переноса в группе , определяются по формулам:

3. S_n приближение метода дискретных ординат.

В качестве угловой сетки в X-Y-Z и HEX-Z геометриях используется полностью симметричная квадратура [9] для небольших значений . Общее число направлений в расчете составляет: , где - геометрическая размерность.

В криволинейной R-φ-Z геометрии используется квадратура Карлсона с равными весами [10].

Угловая сетка разбивает сферу направлений на сектора одинаковой площади. На интервале задаются узлы , . Для каждого интервал изменения угла , делится на подинтервалы , , где , . Полученное разбиение симметрично отображается на часть полусферы соответствующую интервалу . Диаграмма расположения точек для угловой квадратуры на октанте , представлена на Рис.1.

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

Рис. 1. Диаграмма расположения точек угловой квадратуры на октанте.

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

В пакете «РЕАКТОР» используются эффективные конечно-разностные схемы. Численные методы различаются набором дополнительных соотношений, замыкающих систему уравнений. Для решения задач на используется известная схема fix-up с обнулением отрицательных потоков на выходе и последующим учетом баланса в ячейке. Для решения задач защиты используется схема WDIF, которая хорошо зарекомендовала себя при расчетах по американской программе TORT [11]. Практические расчеты показали, что схема WDIF наиболее приспособлена для расчета этих задач, особенно в криволинейной геометрии. Хотя эта схема не является положительной для отдельно взятого направления, но в скалярных потоках дает практически всюду положительное решение при оптимальном выборе параметра . Рассмотрим эту схему более подробно.

Балансное уравнение метода дискретных ординат в X-Y-Z геометрии имеет вид:

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

Подставляя соотношения (2) в (1), получим формулу для расчета плотности потока в ячейке. Существует проблема выбора параметров (), которые принадлежат интервалу . В схеме WDIF выбор весовых параметров осуществляется на основе решения вспомогательных задач. Например, для направления по оси X значение берется из решения следующей системы уравнений:

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

Выбор параметра существенно влияет на результаты расчета. Параметр эмпирически подбирается так, чтобы получить более точное решение. В пакете «РЕАКТОР» по умолчанию полагается .

Для ускорения внутренних итераций в пакете «РЕАКТОР» первоначально [2] использовался диффузионный синтетический метод (DSA) [12]. В некоторых задачах использование этого линейного метода приводило к сокращению времени счета. К сожалению, существуют задачи, для решения которых использование этого метода может привести к неустойчивому алгоритму. К таким задачам относятся задачи с оптически толстыми областями, а также с областями с высокой анизотропией рассеяния. Именно для этих задач наиболее необходимо использовать ускорение.

В настоящее время в пакете «РЕАКТОР» реализован нелинейный метод ребаланса [13]. Существенное различие между этими двумя методами состоит в том, что в методе ребаланса используется поправка, полученная из решения диффузионной системы уравнений, коэффициенты которой находятся из угловых потоков, полученных на этой же внутренней итерации. Усовершенствование метода ребаланса, предложенное в работе [13] в одномерной плоской геометрии, позволило повысить его эффективность. Численные расчеты показали, что метод ребаланса эффективнее DSA метода. Ниже приводится описание этого метода.

На первых итерациях эффективно работает Метод группового ребаланса. Разобьем источник в правой части уравнения переноса на две части, отделяя внутригрупповое рассеяние () от остальной части источника ():

Просуммируем уравнение переноса для группы по пространственной области и по направлениям. Результат может быть записан в операторном виде:

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

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

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

Исключая из (4) и (5) и , получим формулу для вычисления множителя :

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

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

Результатом итерации с номером является следующее уравнение, где волной обозначен неребалансный результат, который не удовлетворяет уравнению (8):

Если предположить, что скалярный ток пропорционален потоку в ячейке:

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

получим линейную систему уравнений относительно для внутренних точек пространственной области:

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

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

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

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

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

Как уже отмечалось выше, при решении уравнения (11) может быть получено отрицательное, либо расходящееся решение , если значения сечения рассеяния велики. Благодаря усовершенствованию метода ребаланса, предложенному в работе [13], эта неустойчивость должна исчезнуть. Как известно, необходимым и достаточным условием сходимости итерационного метода является условие, при котором все собственные значения матрицы перехода по модулю меньше 1. Применительно к диффузионному уравнению (11) это условие выполняется только при достаточном увеличении его коэффициентов. Тогда амплитуда ребалансных поправок уменьшится, что позволит приблизить поправки к значению 1 в тех областях, где они велики.

Процедура усовершенствования заключается в следующем. Сначала как обычно вычисляются скалярные токи. Затем каждый из них увеличивается на фиктивный ток следующим образом:

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

С учетом этих изменений уравнение ребаланса (11) примет следующий вид:

Полученная система уравнений (12) решается итерационным методом верхней релаксации с автоматическим выбором параметра релаксации в каждой группе.

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

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

1) максимальная относительная ошибка решения уравнения переноса до ускорения:

2) относительная ошибка в этой же точке после ускорения методом ребаланса:

а) , т.е. ребаланс действует в противоположную сторону, чем основная итерация,

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

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

6. Ускорение внешних итераций по источнику деления.

Рассмотрим степенной метод вычисления и источника деления :

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

В пакете «РЕАКТОР» для ускорения внешних итераций до недавнего времени использовался метод полиномов Чебышева. Для использования этого метода ускорения надо предварительно оценить отношение двух наибольших собственных чисел матрицы перехода , что является затруднительным. С другой стороны, этот алгоритм слишком чувствителен к вычислительной погрешности [15].

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

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

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

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

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

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

где символом обозначено значение параметра экстраполяции:

Алгоритм ускорения внешних итераций основан на использовании формул (2) и (3) для линейной экстраполяции источника деления в каждой точке пространства. Одновременно экстраполируются скалярные потоки по такой же формуле. Для вычисления доминантного соотношения более удобно использовать глобальную величину ошибки :

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

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

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

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

В заключение заметим, что в [15] этот метод ускорения описан для решения систем линейных уравнений и называется - процесс ускорения сходимости.

7. Ускорение внешних итераций для подкритических систем.

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

то уравнение переноса может быть записано в операторном виде:

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

Исключая из (3) и (4) , получим формулу для вычисления множителя :

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

Результаты расчетов по пакету «РЕАКТОР» показали хорошее совпадение с расчетами реактора СВБР 75/100 по программе TORT [7] на одинаковых константах.

Расчет эффективного коэффициента размножения выполнялся в HEX-Z геометрии в диффузионном приближении и в кинетическом S4P1 приближении для 30 групп нейтронов. В задаче использовалось 18 различных горизонтальных картограмм. Общее число ячеек сетки 1791582. На Рис. 1 и 2 показаны результаты расчета энерговыделения в диффузионном и кинетическом S4P1 приближениях соответственно.

На Рис. 3 показано поведение значений на внешних итерациях при решении этой задачи в диффузионном приближении. Задавалось ограничение на число внутренних итераций не больше 9. На рисунке видно влияние экстраполяции источника деления при ускорении внешних итераций (кривая немонотонно изгибается и поднимается вверх к точному решению).

Расчет защиты для реакторной установки СВБР 75/100 выполнялся в X-Y-Z и R-FI-Z геометриях в приближении S8P3, 30 групп нейтронов и 19 групп g-квантов. В X-Y-Z геометрии 2417637 пространственных точек: 237 слоев по Z, 101 точек по X и 101 точек по Y. В R-FI-Z геометрии рассчитывались меньшие по объему области с более подробной сеткой. На Рис. 4. показано поведение плотности полного потока нейтронов в R-FI-Z геометрии.

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

1. A.V. Voronkov, A.N. Chebeskov, E.P. Sychugova, I.Y. Krivitsky, E.V. Matveeva, Y.N. Mironovich, A.D. Knipe. Low Reactivity Sodium-Void Benchmark Study in an Annular Heterogeneous Assembly. Proceeding of the International Topical Meeting on Sodium Cooled Fast Reactor Safety, IPPE, Obninsk, 1994.

2. А.В. Воронков, Л.П. Басс, Е.П. Сычугова. Пакет программ KINRZ для решения многогруппового уравнения переноса нейтронов и γ – квантов в двумерной R-Z геометрии методом дискретных ординат. Препринт ИПМ РАН, 1999, №83.

3. А.В. Воронков, Е.П. Сычугова, Е.В. Матвеева, А.Н. Чебесков. Расчетные исследования натриевого пустотного коэффициента реактивности (НПЭР), измеренного в экспериментах «ZEBRA». Препринт ИПМ РАН, 1996, №65.

4. А.В. Воронков, Е.П. Сычугова. Решение уравнения переноса нейтронов в двумерной R-Z и трехмерной X-Y-Z геометриях методом дискретных ординат. Препринт ИПМ РАН, 1995, №6.

5. А.В. Воронков, Е.П. Сычугова, А.С. Голубев. Решение стационарного уравнения переноса нейтронов в трехмерной Hex-Z геометрии. Отчет ИПМ РАН, 2004, № 7-03-2004 (37c.)

6. А.В. Воронков, Е.П. Сычугова, П.Б. Афанасьев, А.В. Дедуль, В.В.Кальченко. Трехмерные расчеты радиационной защиты в пакете программ РЕАКТОР-ГП. Тезисы доклада на IX Российской научной конференции «Радиационная защита и радиационная безопасность в ядерных технологиях» 24-26 октября 2006 г., Федеральное Агентство по Атомной Энергии, ГНЦ РФ Физико-Энергетический Институт им. А.И. Лейпунского, г. Обнинск, с. 52-54.

7. Е.В. Аверченков, П.Б. Афанасьев, А.В. Дедуль, В.В. Кальченко, А.В.Воронков, Е.П. Сычугова. Трехмерные расчеты радиационной защиты реакторной установки СВБР 75/100. Тезисы доклада на IX Российской научной конференции «Радиационная защита и радиационная безопасность в ядерных технологиях» 24-26 октября 2006 г., Федеральное Агентство по Атомной Энергии, ГНЦ РФ Физико-Энергетический Институт им. А.И. Лейпунского, г. Обнинск, с. 255-257.

8. Аверченков Е.В., Афанасьев П.Б., Дедуль А.В., Кальченко В.В., ВоронковА.В., Сычугова Е.П. Совместные расчеты «РЕАКТОР» + «ЗАЩИТА» реакторной установки СВБР 75/100 в трехмерной геометрии. Тезисы докладов на 17-м семинаре «НЕЙТРОНИКА – 2006» Нейтронно-физические проблемы атомной энергетики, 30 октября -03 ноября 2006, г. Обнинск.

9. «Вычислительные методы в физике реакторов», Сб. статей под ред. Х. Гринспена, К. Келбера, Д. Окрента, Москва, Атомиздат, 1972.

10. B.G. Carlson «A Method of Characteristics and Other Improvements in Solution Methods for the Transport Equation», Nucl. Sci. Eng., v. 61, 408 (1976).

11. DOORS3.2 One, Two- and Tree Dimensional Discrete Ordinates Neutron/Photon Transport Code System. RSIC Computer Code Collection CCC-650, 1998.

12. E.M. Gelbard and L.A. Hegeman. The Synthetic Method as Applied the S_n Equations. Nucl. Sci. Eng., v. 37, 288, (1969).

13. W. H. Reed. The Effectiveness of acceleration Techniques for Iterative Methods in Transport Theory. Nucl. Sci. Eng., v. 45, 245, (1971).

14. W.A. Rhoades. Improvements in Discrete Ordinates Acceleration. Trans. Am. Nucl. Soc., v. 39, 753 (November 1981).

15. Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков. Численные методы. Изд-во МГУ им. М.В. Ломоносова, 3-е издание, дополненное и переработанное, Москва, 2004.