CC BY
154
ТЕХНИЧЕСКИЕ НАУКИ
УДК 519.63
А. В. Аввакумов, А. О. Васильев, П. Е. Захаров
ПРОГРАММНАЯ РЕАЛИЗАЦИЯ МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ ДЛЯ УРАВНЕНИЯ ДИФФУЗИИ НЕЙТРОНОВ
Для инженерных расчетов ядерных реакторов широко используется уравнение переноса нейтронов в диффузионном приближении. Обычным методом решения является конечно-разностная аппроксимация. Для получения результатов с требуемой точностью необходимо использовать либо очень подробную сетку, либо дополнительные подходы, например, нодальные методы. Мы предлагаем другой подход к решению [1] уравнения диффузии нейтронов, который базируется на методе конечных элементов (МКЭ) и использовании современных методов расчета на базе свободного программного обеспечения. Здесь для повышения точности есть 2 способа: дробить сетку и/или повышать порядок конечного элемента. Разработан нейтронно-физический код (программа), который решает стационарное уравнение диффузии в 2D и 3D гексагональной геометрии в приближении двух энергетических групп, имеет графическую оболочку (интерфейс). Для численного решения используется метод конечных элементов, реализованный в вычислительном пакете FEniCS [2], библиотека для решения спектральных задач SLEPc [3], а для построения и генерации сетки - программа Gmsh [4]. Для тестирования вычислительного алгоритма было рассмотрено несколько численных тестов (получены очень хорошие результаты) и проведено распараллеливание кода в тестовом режиме.
Ключевые слова: ядерный реактор, активная зона, водо-водяной энергетический реактор (ВВЭР), математическое моделирование, уравнение диффузии нейтронов, двухгрупповое приближение, задача на собственные значения, эффективный коэффициент размножения, метод конечных элементов, параллельное программирование, нейтронно-физический код.
АВВАКУМОВ Александр Владимирович - с. н. с. научно-исследовательского центра «Курчатовский Институт».
E-mail: avvakumov2009@rambler.ru
AVVAKUMOV Alexander Vladimirovich - Senior Researcher, National Research Center Kurchatov Institute.
E-mail: avvakumov2009@rambler.ru
ВАСИЛЬЕВ Александр Олегович - ведущий инженер научно-исследовательской кафедры вычислительных технологий ИМИ СВФУ им. М. К. Аммосова.
E-mail: haska87@gmail.com
VASILYEVAlexander Olegovich - Leading Engineer, Research Department of Computing Technologies, Institute of Mathematics and Information Science, North-Eastern Federal University named after M. K. Ammosov.
E-mail: haska87@gmail.com
ЗАХАРОВ Петр Егорович - к. ф.-м. н., доцент-исследователь научно-исследовательской кафедры вычислительных технологий СВФУ им. М. К. Аммосова.
E-mail: zapetch@gmail.com
ZAKHAROVPetr Egorovich - Candidate of Physical and Mathematical Sciences, Associate Professor of Research Department of Computing Technologies, Institute of Mathematics and Information Science, NorthEastern Federal University named after M. K. Ammosov.
E-mail: zapetch@gmail.com
A. A. Avvakumov, A. O. Vasilyev, P. E. Zakharov
Software Implementation of the Finite Element Method for Neutron Diffusion Equation
The most of engineering codes, used in nuclear reactor calculations, are based on diffusion approximation of the neutron transport equation. The finite-difference method is widely used to solve the diffusion equation. Acceptable accuracy of the results reached only by using fine mesh or complementary approaches, such as nodal methods. We propose another approach based on finite-element method (FEM) and the use of modern methods of calculation based on free software. There are two ways to improve the accuracy: split up the grid and / or increase the order of finite element. There was created the neutronic code (software), which solves the stationary diffusion equation in 2D and 3D hexagonal geometry in the approximation of two energy groups, and has a GUI (interface). For the numerical solution it is used the finite element method implemented in the computer package FeniCS, a library for solving spectral problems SLEPc, and for the construction and mesh generation it is used the program Gmsh. For the testing of computational algorithm authors considered several numerical tests that produced good results. There was also conducted parallelization of code (in test mode).
Keywords: nuclear reactor, active zone, water-water energetic reactor (VVER), mathematical modeling, neutron diffusion equation, two-group approximation, eigenvalue problem, effective multiplication factor, finite element method, parallel programming, neutronic code.
Введение
Стремительное развитие атомной энергетики во второй половине прошлого века стимулировало разработку эффективных методов математического моделирования уравнения переноса нейтронов. Вопрос безопасности ядерных реакторов стал важнейшим после нескольких аварий на АЭС в мире. Новые стандарты безопасности поставили перед инженерами и учеными, которые занимаются эксплуатацией и проектированием ядерных реакторов, важные цели повышения качества моделирования физических процессов в ядерном реакторе. В связи с этим разработка новых методов и алгоритмов расчета реакторов получила дополнительное ускорение.
Другим важным фактором, влияющим на развитие методов моделирования ядерных реакторов, является энергетика. Каждый день простоя АЭС из-за нештатных ситуаций ведет к финансовым потерям. Большая доля нештатных ситуаций связана с человеческим фактором, ошибками операторного персонала, численность которого может быть снижено повышением квалификации. Для этого следует создавать новые и совершенствовать существующие тренажеры, для чего необходимы более совершенные модели и современные про -граммные продукты. Требование работы тренажера в режиме реального времени зачастую приводит к значительному упрощению использующихся в тренажере моделей основных физических процессов. Постоянное ускорение вычислительной техники и использование параллельных вычислений позволяет применять более сложные модели.
Основу современной российской атомной отрасли составляют реакторы типа ВВЭР. В ближайшие десятилетия этот тип реакторов останется доминирующим, поэтому разработки эффективных методов и современных программных продуктов нейтронно-физических процессов в активной зоне реактора типа ВВЭР остается актуальной.
Физические процессы, происходящие в активной зоне ядерного реактора [5], зависят от распределения нейтронного потока, математическое описание которого основывается на уравнении переноса нейтронов [6]. В общем виде это уравнение имеет интегро-дифферен-циальную форму, а искомое распределение потока нейтронов зависит от времени, энергии, пространственных и угловых переменных. Для практических расчетов ядерных реакторов, как правило, используют упрощенные формы уравнения переноса нейтронов. Наибольшее распространение для анализа реакторов получило уравнение, известное как групповое диффузионное приближение, [7-8] которое используется в подавляющем большинстве
инженерных нейтронно-физических кодов. В практике реакторных расчетов особое место занимает решение условно-критической задачи [9] (которая в математической формулировке сводится к задаче на собственные значения), характеризующее эффективный коэффициент размножения нейтронов. Собственной функцией этой задачи является стационарное распределение потока нейтронов.
Постановка задачи
Текущая версия разработанного программного продукта решает уравнение диффузии нейтронов в двухгрупповом приближении:
^(ДУФ,) + ( + Ег ) = К- (1Ф1 +у2 Е, 2 Ф 2),
(1)
-У(2УФ 2 ) + Еа2Ф 2 = Ег Ф1.
С граничным условиям альбедного типа:
дФ 8
= ^Ф 8'8 =1>2' (2)
Здесь групповые параметры О (г), Хвд(г), ^(г), Х4(г) - коэффициент диффузии, макросечение поглощения, число вторичных нейтронов и макросечение деления соответственно, а X (г) - макросечение рассеяния, п - внешняя нормаль границы дО, у^ - групповой альбедный параметр (логарифмическая производная). Решением уравнений (1)-(2) является эффективный коэффициент размножения нейтронов Ке// и стационарное распределение нейтронного потока Ф(г).
Для численного решения задачи методом конечных элементов уравнения (1)-(2) необходимо привести к вариационной постановке [10]. Стандартным способом перевода дифференциальной задачи в вариационную является умножение уравнения на некую функцию V, интегрирование полученного уравнения по области, замена производных второго порядка через интегрирование по частям. Функция V называется тестовой функцией, а искомая функция - пробной функцией. В нашем случае мы каждое уравнение умножаем на тестовую функцию, первое - на v1, второе - на v2 и интегрируем полученные уравнения по области О. Тогда получаем:
-[ V (ЦУФ!) У^г +[ ( + 2) ф^г = — [ (+ у2 2 г 2Ф 2) У2 ^г а а К а (3)
-^(^Ф 2 ) dr +[2а 2Ф 2 У2 dr = [ 2 Фу dr.
Далее заменяем первые интегралы с помощью интегрирования по частям и используем формулу Гаусса-Остроградского для перехода к поверхностному интегралу:
(Д УФ1 )г = | ((УФ1, Vv1) ^г - | D1v1 ds
дп
дФ (4)
2 )v2dr = 2, Vv2 )dr - |D2^-2ds.
а а да дп
Тогда из (2) и (3) получаем следующую систему уравнений: |а (ЦУ^УУ^Г + £ (Ха1 + £ )ф^г + = £ £{1 ф + У2 £{2 ф
" (5)
£ (Ц2УФ + а2 Ф 2 У2 ^ + £^2 Ф2 У2 ^ = £ £г Ф^2 ^.
Полученная вариационная задача формулируется следующим образом: нужно найти такие
функции Фg е V, которые удовлетворяют системе уравнений (5) для любых vg е V, где V - пространство пробных функций, а V - пространство тестовых функций. Здесь V=H1(Q), V = Н1 , где Н1(О) - пространство Соболева, состоящее из функций vg таких, что и | | имеют конечный интеграл в О.
Далее мы должны перейти от непрерывной вариационной задачи (5) к дискретной задаче. Введем конечномерные пространства Ук с V, V с V и определим в них дискретную задачу: найти Ф^ е Vh такие, что
/п (АУФ1к У1к +(Х а1 + Хг )Ф1 Л^г + \еаУ1ф1„У1кЖ = 1 с
= (У1 Х/1 Ф1к + ^ Х/2 Ф2к (6) I (о2 УФ 2 к ^ к + |аЕ
а 2 Ф2кУ2к
Перепишем уравнение (6) как задачу на собственные значения в следующей операторной форме:
Ш = XF Ф, (7)
где М - оператор, описывающий убыль (потерю) нейтронов за счет процессов переноса (утечки), поглощения и рассеяния, а F - оператор, описывающий образование (генерацию) нейтронов за счет процессов деления и рассеяния из верхней области энергий. Операторы М и F являются блочными:
м.(М 0 1, г.(г г-1. (8)
IМл М„1 { 0 0 1 ()
Ищем наибольшее собственное число Х1, характеризирующий эффективный коэффициент размножения нейтронов К =Ш а соответствующая ему собственная функция Ф(г) описы-
еЛ 1
вает стационарное распределение нейтронного потока в данной системе (г е
Одной из основных задач в физике ядерных реакторов является оценка различных функционалов нейтронного потока. Определим нейтронную мощность Р(г) как следующий функционал:
Р = А ( / 1Ф1 / 2 Ф 2 ), (9)
где А - коэффициент нормировки на заданное значение интегральной мощности.
Нейтронно-физический код (программа)
Существующие инженерные нейтронно-физические коды разработаны, как правило, под один тип реакторных установок (РУ): есть коды, предназначенные только для расчета реакторов типа ВВЭР (например, БИПР-7). Подобная тенденция разработки нейтронно-фи-зических кодов применительно к конкретному типу РУ сохраняется до настоящего времени. Используемые в этих кодах подходы и приближения, применяемые для определенного типа РУ, как правило, заложены в самом структуре алгоритмов. Это может быть определенный тип геометрии, фиксированная расчетная сетка, ограничение на число энергетических групп и т. д. В результате имеется набор разнородных кодов, обладающих ограниченным набором возможностей. Это делает практически невозможным использование этих кодов не только для другого типа РУ, обладающего своими специфическими особенностями, но также в случае модернизации проекта данного типа РУ (например, при переходе от одного типа геометрии к другой).
В этой связи возникает вопрос: можно ли разработать некую универсальную расчетную платформу нейтронно-физического кода, на базе которой можно было бы проводить
инженерные расчеты РУ различных типов, при этом не видоизменяя структуры кода, а лишь используя соответствующие расчетные модули. Из всех известных методов, пожалуй, лишь метод Монте-Карло соответствует этой цели. Нейтронно-физические коды, использующие метод Монте-Карло, такие как MCNP или MCU, без каких-либо ограничений применяются для расчета реакторных систем любого типа. Однако использование подобных кодов в инженерных расчетах РУ по различным причинам очень проблематично. Рассмотрим подходы и модели, используемые в современных инженерных нейтронно-физических кодах.
Для большинства инженерных нейтронно-физических кодов геометрическая модель активной зоны имеет схожую структуру. Это, как правило, набор призм (гексагональных или прямоугольных), которые моделируют тепловыделяющую сборку (ТВС), причем поперечное сечение призмы повторяет форму ТВС. Торцевой и радиальный отражатели, экраны моделируются призмами такой же формы. По высоте геометрическая модель активной зоны разбивается на несколько аксиальных слоев или зон, в пределах которых физические и теплогидравлические свойства, определяющие коэффициенты уравнений диффузии, считаются одинаковыми. Отдельный вопрос - формирование расчетной сетки для моделирования нейтронно-физических процессов в РУ Большинство кодов имеют дело с собственной сеткой, сгенерированной самим кодом, использование которой другими кодами практически исключено. Обычно сетки в плане имеют регулярную структуру (например, шестиугольная ячейка ТВС ВВЭР может разбиваться на несколько равносторонних треугольников). По высоте сетка может быть нерегулярной, однако ее структура остается постоянной по всему расчетному объему, т. е. возможность локального сгущения сетки по высоте, например в пределах отдельной ТВС, не допускается.
Таким образом, учитывая схожесть и различия в подходах и моделях, используемых в современных инженерных нейтронно-физических кодах, сформулируем принципиальные положения, которые мы заложили в расчетную платформу нейтронно-физического кода.
Во-первых, на данном этапе рассматриваем уравнение переноса нейтронов в диффузионном групповом приближении. В перспективе возможно расширение функциональных возможностей кода за счет использования недиффузионных моделей, например, 8Р3-приближения. Учитываются последние достижения в области численных методов и алгоритмов, позволяющие обеспечить необходимую точность и быстродействие нейтрон-но-физических расчетов, включая применение современных технологий программирования с использованием объектно-ориентированного подхода. Используются свободные программные обеспечения, которые дают ощутимые преимущества, обусловленные мультиязычным (Qt, Python) интерфейсом (рис. 1), большой надежностью и модификацией программ.
а— 1ЧИ" Neutron Diffusion
Рис. 1. Стартовая страница программы
Во-вторых, применяем традиционный подход при создании геометрической модели активной зоны как набора призм, имитирующих сборку ТВС, с использованием гомогенизированного представления физических зон. В перспективе возможно использование более подробных моделей (например, более корректное описание бокового отражателя как негомогенной среды).
В-третьих, нейтронно-физический код работает с произвольными расчетными сетками (сгенерированными в стандартном формате). Расчетная сетка может формироваться независимо (т. е. вне нейтронно-физического кода) с помощью внешних генераторов сеток. Данный подход позволяет построить и протестировать расчетную сетку любой сложности. Базовая сетка формируется самим кодом с помощью встроенного генератора сеток.
В-четвертых, переходим от конечно-разностной аппроксимации к методу конечных элементов. В этом случае отпадает необходимость в дополнительных методах повышения точности конечно-разностного расчета (нодальные методы), поскольку повысить точность расчета можно как сгущением сетки, так и повышением порядка конечного элемента.
В нейтронно-физическом коде использовалось свободное программное обеспечение, включающее следующие компоненты:
• генератор сеточных объектов — Gmsh;
• вычислительный пакет FEniCS как инструмент для решения краевых задач, описываемых дифференциальными уравнениями в частных производных, методом конечных элементов;
• решатель задачи на собственные значения на основе пакета высокопроизводительной линейной алгебры, поддерживающей параллельные вычисления — SLEPc;
• система пост-процессорной обработки и визуализации — библиотека VTK.
Тестовые расчеты
Для тестирования вычислительного алгоритма рассмотрим несколько численных тестов, имитирующих различные конфигурации загрузок типа ВВЭР с гексагональной кассетной структурой. В расчетах варьировались следующие параметры:
• npa - параметр, характеризующий детализацию расчетной сетки - число расчетных ячеек (конечных элементов) на кассету;
• p - порядок конечных элементов (диапазон изменения p: от 1 до 3);
Вычислялись следующие параметры:
• эффективный коэффициент размножения Kf
• распределение нейтронной мощности P(9) по кассетам с нормировкой на среднее значение по активной зоне.
С целью анализа сходимости и эффективности разработанного алгоритма все тестовые расчеты были выполнены с фиксированной точностью отгонки собственного числа, равной 10-15. Сравнение полученных результатов проводилось с результатами расчетов по диффузионной мелкосеточной программе DIF3D-FD [11] (эталонное решение было получено путем экстраполяции результатов на бесконечно малый размер элементарной ячейки расчетной сетки). Будем рассматривать следующие отклонения в расчетных параметрах:
• для эффективного коэффициента размножения абсолютное отклонение от «эталонного» значения Keff:AK=\Keff-KrJ выражается вpcm (percent-milli, т. е. 10-5);
• для распределения покассетных мощностей P вычисляются относительные отклонения е. (выражаются в %):
P — Pref
£i =——i—,
i рЩ
pref
где ' - «эталонное» значение мощности в кассете /(/-1,..., N).
• по отклонениям е. рассчитываются интегральные отклонения:
■ среднеквадратическое отклонение RMS:
Рис. 2. Геометрическая модель теста IAEA-2D с отражателем
RMS=J„E
Nt?
среднее по модулю отклонение AVR:
1 N
AVR = —Ykl, Nti1 '''
■ максимальное по модулю отклонение MAX:
MAX = max .
Модифицированный тест IAEA-2D с отражателем
Тестовая задача является модификацией на случай гексагональной геометрии известной тестовой задачи IAEA-2D [11]. Геометрическая модель активной зоны реактора моделируется набором кассет гексагональной формы (рис. 2), где различными цветами показаны кассеты различных сортов. Активная зона имеет 13 стрежней СУЗ (устройства систем управления и защиты реактора) и 1/12 зеркальную симметрию. Размер кассеты «под ключ» равен 20 см.
Диффузионные нейтронно-физические константы заданы в табл. 1. Граничные условия задаются в виде логарифмической производной.
На рис. 3 показана расчетная сетка при npa=6. Результат расчета плотности нейтронного потока для быстрых нейтронов при 7=0,5, p=1 показан на рис. 4, а покассетная мощность при таких же параметрах - на рис. 5.
Таблица 1
Диффузионные константы для модифицированного теста IAEA-2D
2
Материал 1 2 3 4
D1 1,50 1,50 1,50 1,50
D2 0,40 0,40 0,40 0,40
Xa1 0,01 0,01 0,01 0,00
Xa2 0,08 0,085 0,13 0,01
Xr 0,02 0,02 0,02 0,04
Xf1 0,00 0,00 0,00 0,00
Xf2 0,056 0,056 0,056 0,00
V1Xf1 0,00 0,00 0,00 0,00
V2Xf2 0,135 0,135 0,135 0,00
П1Ч- т- Ч*.'-■ Й . ._.. ■ ,- — -... V -
<5 ■ ■ ■ ф.
Ш $
■-,]„ ' Ф ©. ■ ж
:;к Й ^ Ф §
Рис. 3. Расчетная сетка теста IAEA-2D с отражателем
Рис. 4. Плотность нейтронного потока теста IAEA-2D с отражателем
Рис. 5. Покассетная мощность теста IAEA-2D с отражателем
Таблица 2
Результаты расчета модифицированного теста IAEA-2D с отражателем при у=0,5
npa Р Keff(pCm) AK(pcm) RMS(%) AVR(%) MAX(%) t(sec)
1 1,0104126 490,56 13,29 11,13 23,73 0,03
6 2 1,0062265 71,95 1,88 1,59 3,40 0,11
3 1,0055754 6,84 0,22 0,18 0,41 0,27
1 1,0069873 148,03 4,54 3,77 8,45 0,10
24 2 1,0056090 10,20 0,30 0,25 0,57 0,54
3 1,0055135 0,65 0,02 0,02 0,04 1,46
1 1,0059079 40,90 1,28 1,06 2,41 0,50
96 2 1,0055186 1,16 0,04 0,03 0,07 3,00
3 1,0055097 0,27 0,01 0,01 0,02 8,30
1 1,0056119 10,49 0,34 0,28 0,64 2,95
384 2 1,0055102 0,32 0,01 0,01 0,02 18,70
3 1,0055096 0,26 0,01 0,01 0,02 52,00
1 1,0055354 2,84 0,09 0,08 0,17 18,65
1536 2 1,0055096 0,26 0,01 0,01 0,02 132,30
3 1,0055096 0,26 0,01 0,01 0,02 370,50
ref. 1,0055070
В табл. 2 показана сводная таблица результатов расчета модифицированного теста IAEA-2D с отражателем при разных параметрах. Здесь приняты следующие обозначения: npa - число ячеек на кассету; p - порядок конечного элемента; Keff- эффективный коэффициент размножения; AK - абсолютное отклонение от «эталонного» значения; RMS - среднеквадратичное отклонение; AVR - среднее отклонение; MAX - максимальное отклонение; t - время счета. На рис. 6 показано ускорение расчетов в зависимости от числа процессоров.
Из табл. 2 можно сделать следующие выводы:
• наблюдается устойчивая сходимость решения тестовых задач при увеличении числа конечных элементов на кассету npa и порядка конечных элементов p;
• с точки зрения экономичности расчета, увеличение порядка конечных элементов p намного эффективнее увеличения числа конечных элементов на кассету npa.
Рис. 6. Эффективность распараллеливания теста IAEA-2D c отражателем
Заключение
Разработана программа, которая решает стационарное уравнение диффузии в 2D и 3D гексагональной геометрии в приближении двух энергетических групп:
• используются современные вычислительные средства на базе свободного программного обеспечения;
• имеется графическая оболочка — интерфейс;
• просчитаны численные тесты, получены очень хорошие результаты.
Л и т е р а т у р а
1. Аввакумов А. А., Вабищевич П. Н., Васильев А. О. Метод конечных элементов для уравнения диффузии нейтронов в гексагональной геометрии // Вестник СВФУ. - Т. 11, № 5. - 2014.
2. Logg, K. A. Mardal, G. Wells. Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book. Lecture Notes in Computational Science and Engineering. Springer, 2012.
3. Campos, J. E. Roman, E. Romero, A. Tomas. SLEPc Users Manual, 2013.
4. C. Geuzaine, J. F. Remacle. Gmsh Reference Manual, 2014.
5. Климов А. Н. Ядерная физика и ядерные реакторы. - М.: Атомиздат, 1971.
6. Stacey W. M. Nuclear Reactor Physics. John Wiley & Sons, 2007.
7. Ганев И. Х., Доллежаль Н. А. Физика и расчет реактора. - М., Энергоиздат, 1981.
8. Марчук Г И., Лебедев В. И. Численные методы в теории переноса нейтронов. - М., Атомиздат, 1981.
9. S. Gonzalez-Pintor, G. Verdu, D. Ginestar. Approximation of the neutron diffusion equation on hexagonal geometries. - New York, 2009.
10. A. Hebert. Application of a dual variational formulation to finite element reactor calculations. Annals of Nuclear Energy, 20:823-845, 1993.
11. Y. A. Chao, Y. A. Shatilla. Conformal mapping and hexagonal nodal methods-ii: Implementation in the anc-h code. Nuclear Science and Engineering, 121:210-225, 1995.
R e f e r e n c e s
1. Avvakumov A. A., Vabishchevich P. N., Vasil'ev A. O. Metod konechnykh elementov dlia uravneniia diffuzii neitronov v geksagonal'noi geometrii // Vestnik SVFU. - T. 11, № 5. - 2014.
2. Logg, K. A. Mardal, G. Wells. Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book. Lecture Notes in Computational Science and Engineering. Springer, 2012.
3. Campos, J. E. Roman, E. Romero, A. Tomas. SLEPc Users Manual, 2013.
4. C. Geuzaine, J. F. Remacle. Gmsh Reference Manual, 2014.
5. Klimov A. N. Iadernaia fizika i iadernye reaktory. - M.: Atomizdat, 1971.
6. Stacey W. M. Nuclear Reactor Physics. John Wiley & Sons, 2007.
7. Ganev I. Kh., Dollezhal' N. A. Fizika i raschet reaktora. - M., Energoizdat, 1981.
8. Marchuk G. I., Lebedev V. I. Chislennye metody v teorii perenosa neitronov. - M., Atomizdat, 1981.
9. S. Gonzalez-Pintor, G. Verdu, D. Ginestar. Approximation of the neutron diffusion equation on hexagonal geometries. - New York, 2009.
10. A. Hebert. Application of a dual variational formulation to finite element reactor calculations. Annals of Nuclear Energy, 20: 823-845, 1993.
11. Y. A. Chao, Y. A. Shatilla. Conformal mapping and hexagonal nodal methods-ii: Implementation in the anc-h code. Nuclear Science and Engineering, 121: 210-225, 1995.
^SHir^ir
