[ «РАЗВИТИЕ МЕТОДА ПОВЕРХНОСТНЫХ ГАРМОНИК ДЛЯ РЕШЕНИЯ ЗАДАЧ НЕЙТРОННОЙ ПРОСТРАНСТВЕННОЙ КИНЕТИКИ В ЯДЕРНЫХ РЕАКТОРАХ ...»
НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ ЦЕНТР
КУРЧАТОВСКИЙ ИНСТИТУТ
На правах рукописи
УДК 621.039.514
Кондрушин Антон Евгеньевич
РАЗВИТИЕ МЕТОДА ПОВЕРХНОСТНЫХ ГАРМОНИК ДЛЯ РЕШЕНИЯ
ЗАДАЧ НЕЙТРОННОЙ ПРОСТРАНСТВЕННОЙ КИНЕТИКИ
В ЯДЕРНЫХ РЕАКТОРАХ
Специальность: 05.13.18 «Математическое моделирование, численные методы и комплексы программ»Диссертация на соискание ученой степени кандидата физико-математических наук
Научный руководитель:
д.т.н. Бояринов В.Ф.
Москва – Оглавление Введение…………………………………………………………………………... ….. Глава 1 Обзор литературы……………………………………………………....... 1.1 Общие уравнения………………………………………………………......... 1.2 Обзор проекционно-сеточных методов решения уравнения переноса нейтронов…………………………………………………………………...... 1.2.1 Метод конечных элементов…………………………………….......... 1.2.2 Метод граничных элементов………………………………………...... 1.2.3 Метод конечных суперэлементов Федоренко……………………..... 1.2.4 Метод матриц отклика……………………………………………....... 1.3 Метод поверхностных гармоник………………………………………….... 1.3.1 Основы метода поверхностных гармоник и определения………..... 1.3.2 Обзор программных комплексов SUHAM и SVS………………....… 1.3.3 Обзор работ по нестационарным уравнениям МПГ……………....… 1.4 Обзор методов решения нестационарного уравнения переноса нейтронов………………………………………………………………….. …. 1.4.1 Полностью явный метод…………………………………………… …. 1.4.2 Полностью неявный метод………………………………………….… 1.4.3 -метод……………………………………………………………….… 1.4.4 Метод переменных направлений………………………………….. …. 1.4.5 Улучшенный квазистатический метод………………………….….… 1.4.6 SCM метод………………………………………………………….. …. 1.5 Обзор программ для решения нестационарного уравнения переноса нейтронов………………………………………………………………….. …. Заключение к главе 1…………………………………………………………. …. Глава 2 Нестационарные уравнения метода поверхностных гармоник…. …. 2.1 Двумерные нестационарные уравнения МПГ…………………………...… 2.1.1 Поверхностная невязка……………………………………………...… 2.1.2 Объемная невязка…………………………………………………....… 2.1.3 Вывод уравнений МПГ………..…………………………………….… 2.2 Одномерные уравнения МПГ……………………………………………..… 2.3 Итерационная схема………………………………………………………..… Заключение к главе 2…………………………………………………………..… Глава 3 Программный комплекс SUHAM-TD и его верификация………...… 3.1 Программный комплекс SUHAM-TD…………………………………… …. 3.2 Верификация программного комплекса SUHAM-TD………………….. …. 3.2.1 Тест BSS-6………………………………………………………….. …. 3.2.2 Тест PHWR…………………………………………………………. …. 3.2.3 Тест TWIGL………………………………………………………… …. 3.2.4 Модифицированный тест 8-A1……………………………………. …. 3.2.5 Транспортный тест TWIGL………………………………………...... Заключение к главе 3………………………………………………………….... Глава 4 Разработка и расчет пространственно-временного бенчмарка C5G7-TD для тестирования кинетических нейтронно-физических кодов.... 4.1 Обзор пространственно-временных бенчмарков…………
4.2 Описание бенчмарка C5G7………………………………………………..... 4.3 Расчет кинетических характеристик для теста C5G7-TD………………... 4.4 Законы ввода реактивности для теста C5G7-TD………………………...... 4.5 Результаты моделирования теста C5G7-TD……………………………..... Заключение к главе 4………………………………………………………….... Заключение………………………………………………………………………..... Обозначения………………………………………………………………………... Список литературы………………………………………………………………... Приложение А Копии свидетельств о регистрации модулей SUHAM-TD.... Приложение Б Результаты расчета теста BSS-6………………………….….... Приложение В Результаты расчета теста PHWR……………………………... Приложение Г Результаты расчета теста TWIGL…………………………..... Приложение Д Результаты расчета теста C5G7-TD………………………...... Развитие современной атомной энергетики требует повышенного внимания к характеристикам надежности и безопасности ядерных реакторов. Это внимание обуславливается наличием потенциальной возможности возникновения аварии и значительными затратами при строительстве атомных станций, по причине проектирования с запасом с целью предотвращения потенциально возможных аварий. Важнейшую роль в проектировании надежных, безопасных и вместе с тем исследовательских и проектных расчетов, позволяющих приблизится к оптимальному соотношению этих показателей. При этом важнейшее значение имеет проведение качественного нейтронно-физического расчета.
В последние годы все большее внимание уделяется развитию кодов, позволяющих проводить нестационарный расчет ядерного реактора. Это связано с нестационарных кодов содержит ряд серьезных приближений в своей нейтроннофизической составляющей на фоне высоких требований к безопасности реакторов. К таким приближениям в первую очередь относятся пространственная гомогенизация, расчет в малом числе энергетических групп, диффузионное приближение, а также применение разного рода поправок, полученных расчетноэкспериментальным путем и призванных уменьшить ошибки в моделировании по заложенной в расчет математической модели (например [1, 2]). Вместе с тем, следует также отметить, что ряд исследователей (например, [3–5]) указывает на возможность получения улучшенных результатов путем решения более универсального газокинетического уравнения, особенно для анализа безопасности. Результаты, полученные таким путем, позволили бы обрести большую уверенность в качестве получаемых результатов. Таким образом, можно заключить, что существует необходимость в создании нестационарных вышеуказанные приближения.
Как известно, оперативность моделирования процесса переноса нейтронов главным образом зависит от сложности рассматриваемой модели (геометрические производительности вычислительной машины и метода, заложенного в основу моделирующего кода. Отказ от вышеизложенных приближений приводят к значительному усложнению математической модели и увеличению размера решаемой системы уравнений, что напрямую влияет на вычислительные затраты.
Вместе с этим, несмотря на очень интенсивное развитие компьютерной техники в последние десятилетия, моделирование нестационарного процесса переноса вышеперечисленных приближений является трудоемкой задачей даже для современных компьютеров. Что касается методов моделирования переноса нейтронов, то наиболее широко используемые на сегодняшний день методы, такие как метод дискретных ординат, метод характеристик, метод вероятности вычислительных затрат. При этом, применение кода, основанного на методе, позволяющем быстро получать основные нейтронно-физические функционалы в нестационарном расчете с достаточной для практики точностью, могло бы повысить оперативность получения результатов и позволило бы применять его для сложных масштабных моделей с умеренными вычислительными затратами.
эффективных и вместе с тем экономичных алгоритмов.
В качестве метода, который мог бы лежать в основе нейтронно-физической вышеизложенными характеристиками, предлагается использовать предложенный Н.И. Лалетиным [6–10] метод поверхностных гармоник (МПГ). МПГ является процедурой построения математической модели для описания нейтроннофизических процессов в ядерном реакторе, учитывающей особенности нейтронно-физических реакторных задач [9].
Главной причиной для обращения к МПГ стал тот факт, что данный метод, реализованный в стационарных программных комплексах SUHAM [11, 12] и SVS [13], позволяет получать основные нейтронно-физические функционалы с точностью сравнимой с точностью детерминистических методов при вычислительных затратах сравнимых с затратами инженерных подходов. Данные свойства МПГ проявил как при расчете классических сильно идеализированных бенчмарков, так и при расчете моделей как проектируемых (БРЕСТ-ОД-300, БН-1200 [14], ГТ-МГР [15]), так и реальных реакторов (РБМК-1000 [14], ВВЭР-1000 [16, 17]). Данный метод реализован в одно-, двух- и трехмерной геометриях [14, 17], в системах с квадратной и треугольной решеткой. Все это позволяет рассматривать МПГ как хорошо апробированный и состоявшийся метод, проявивший свои положительные черты на широком классе реакторных задач. Кроме этого, МПГ очень удобен для применения алгоритмов параллельных вычислений [18].
Для применения МПГ в нестационарных задачах необходимо провести детальный математический вывод нестационарных уравнений МПГ и их верификацию посредством программной реализации и расчета ряда бенчмарков.
Вышеописанная ситуация делает актуальной работу по развитию МПГ на задачи пространственно-временной кинетики, созданию алгоритмов и кода для их программной реализации, которые можно было бы рассматривать как базовое ядро для создания современного вычислительного инструмента.
Здесь следует заметить, что верификация такого кода на этапе разработки осложняется следующей проблемой. Существующий на сегодняшний день набор пространственно-временных бенчмарков содержит своего рода пробел. С одной стороны существует ряд классических диффузионных бенчмарков (например, TWIGL [19], 8-A1 [20]) представляющих собой несколько гомогенных зон и содержащих диффузионные нейтронно-физические константы в виде малогрупповых макроконстант и прочих величин. Такие тесты не позволяют провести верификацию программы, проводящей расчеты без пространственной гомогенизации и без диффузионного приближения. С другой стороны, имеется набор бенчмарков, которые описывают гетерогенную структуру среды, содержат концентрации нуклидов в качестве характеристик материалов, включают в себя характеристики обратных связей и т.д. Например, к таким задачам можно отнести PWR MOX/UO2 core transient benchmark [21], PBMR coupled neutronics/thermalhydraulics transient benchmark the PBMR-400 core design [22]. Результаты их расчетов содержат дополнительные погрешности (погрешность ядерных данных, погрешности тепло-гидравлических кодов и другие) что не позволяет выделить погрешность метода, заложенного в основу нейтронного кода.
Тест, описывающий гетерогенную среду, содержащий заданные нейтроннофизические макроконстанты материалов для решения уравнения переноса нейтронов без диффузионного приближения и обратных связей, позволил бы избежать излишних сложностей при верификации и мог бы дать больше возможностей для верификации непосредственно метода. Задача с такого рода характеристиками позволяет выделить методическую составляющую погрешности.
Исходя из вышеописанных проблем, формируется цель диссертационной работы – разработка алгоритмов и расчетных программ для решения нейтроннофизических пространственно-временных задач в ядерных реакторах на основе метода поверхностных гармоник для повышения надежности, точности и оперативности предсказания важнейших нейтронно-физических нестационарных характеристик ядерного реактора. Для достижения этой цели решены следующие задачи:
1. проведение детального вывода пространственно-временных уравнений метода поверхностных гармоник посредством классического метода минимизации невязки на основе нестационарного уравнения переноса нейтронов;
2. построение численного алгоритма решения нестационарного уравнения переноса нейтронов в одномерном и двумерном ядерном реакторе на основе полученных уравнений МПГ;
3. программная реализация разработанных алгоритмов в рамках программного комплекса SUHAM-TD [23];
4. верификация созданного кода и проведение с его помощью исследований применения метода поверхностных гармоник в пространственно-временном расчете;
5. разработка нестационарного бенчмарка C5G7-TD для решения уравнения переноса на основе стационарного бенчмарка C5G7 (benchmark on deterministic transport calculations without spatial homogenization) [24] посредством получения нестационарных характеристик материалов.
Проведение расчета предложенного бенчмарка средствами программного комплекса SUHAM-TD.
Научная новизна результатов, представленных в работе, состоит:
• в разработке алгоритмов и их реализации в расчетной программе SUHAM-TD для решения одномерных и двумерных конечно-разностных нестационарных уравнений МПГ в ядерных реакторах с квадратной решеткой;
• в проведении верификации разработанного кода SUHAM-TD на ряде бенчмарков с демонстрацией эффективности уравнений и разработанного кода;
• в создании пространственно-временного бенчмарка C5G7-TD для решения уравнения переноса на основе международного стационарного бенчмарка C5G7 с приведением результатов расчета.
Достоверность полученных результатов. Разработанные алгоритмы реализованы в программном комплексе SUHAM-TD и верифицированы на ряде пространственно-временных бенчмарков в одномерной и двумерной геометриях с различными сценариями ввода реактивности в систему. Достоверность полученных результатов подтверждается сравнением с результатами, полученными другими методами и программами, как автором, так и другими научными коллективами.
Практическая ценность полученных результатов определяется тем, что:
• полученные нестационарные уравнения МПГ могут быть применены для любого типа реакторов с квадратной регулярной решеткой блоков;
• разработанные алгоритмы и модули программного комплекса SUHAM-TD предлагаются в качестве базового ядра для создания современного вычислительного инструмента для анализа нейтронных переходных процессов;
• созданный пространственно-временной бенчмарк C5G7-TD предоставляет возможность для кросс-верификации любых пространственно-временных кодов с возможностью исследования эффекта гомогенизации и применения диффузионного приближения в нестационарном случае.
Личный вклад автора. Диссертант является соавтором научных работ по теме исследования и все основные результаты диссертации получены при непосредственном участии автора, а именно:
• нестационарные уравнения МПГ и расчетные алгоритмы, реализованные в созданных кодах;
• разработанный программный комплекс SUHAM-TD за исключением программ РАЦИЯ [25], DICPN [26];
• верификационные исследования разработанного алгоритма и программ;
• созданный нестационарный тест C5G7-TD и его расчет по программе SUHAM-TD.
На защиту выносятся • алгоритмы решения нестационарных уравнений МПГ с тремя пробными матрицами, реализованные в созданных кодах, для случая реактора с квадратной решеткой блоков;
• разработанные программы комплекса SUHAM-TD;
• результаты верификации разработанного программного комплекса;
• созданный пространственно-временной бенчмарк C5G7-TD.
Апробация работы. Основные положения диссертации докладывались на следующих российских и международных конференциях, школах и семинарах:
• межведомственный ежегодный семинар по нейтронно-физическим проблемам атомной энергетики «Нейтроника-2010» (г. Обнинск, 26 – 28 октября 2010 г.);
• 8-ая Курчатовская молодежная научная школа (г. Москва, 22 – 25 ноября 2010 г.);
• международная конференция по физики ядерных реакторов «PHYSOR»
(г. Ноксвилл, Теннесси, США, 15 – 20 апреля 2012 г.);
• научная сессия НИЯУ МИФИ (г. Москва, 1–6 февраля 2013 г.);
• международная конференция по математическому моделированию и расчету ядерных реакторов «M&C» (г. Сан-Валли, Айдахо, США, 5 – 9 мая 2013 г.);
• межведомственный ежегодный семинар по нейтронно-физическим проблемам атомной энергетики «Нейтроника-2013» (г. Обнинск, 6 – 8 ноября 2013 г.);
• 11-ая Курчатовская молодежная научная школа (г. Москва, 12 – 15 ноября 2013 г.).
рецензируемом научном журнале «Вопросы атомной науки и техники. Серия:
Физика ядерных реакторов».
Структура и объем диссертации. Диссертация состоит из введения, четырех глав, заключения, списка литературы из 122 наименований и пяти приложений, содержит 171 страницу, 66 таблиц и 34 рисунка.
Первая глава содержит краткий обзор проекционно-сеточных методов, нашедших применение для решения нейтронно-физических задач, к которым относится и метод поверхностных гармоник. Дано краткое описание МПГ.
Описаны программные комплексы SUHAM и SVS с рассмотрением круга основных решенных ими задач. Приведен краткий обзор методов для решения нестационарного уравнения переноса нейтронов и работ по нестационарным уравнениям МПГ. Дан краткий обзор основных современных недиффузионных нестационарных нейтронно-физических кодов для расчета реакторной кинетики.
Во второй главе приведен детальный вывод двумерных нестационарных конечно-разностных уравнений МПГ для квадратной решетки блоков посредством применения процедуры минимизации невязки. Приведен результат вывода одномерных уравнений. Проведено построение итерационной схемы.
Третья глава посвящена описанию созданного программного комплекса SUHAM-TD и его верификации на классических пространственно-временных бенчмарках. Проводится исследование применения МПГ в нестационарном случае и демонстрация его эффективности.
В четвертой главе приводится описание созданного пространственновременного бенчмарка C5G7-TD. Приведено описание алгоритма расчета кинетических параметров. Представлены результаты расчета предложенного бенчмарка посредством программного комплекса SUHAM-TD.
Данная глава посвящена обзору методов и программ для решения одной из основных задач анализа нейтронных переходных процессов в ядерном реакторе, а запаздывающими нейтронами [27]:
где v ( E ) – скорость нейтронов с энергией E, (r,, E, t ) – функция распределения нейтронов, t (r, E, t ) – полное макроскопическое сечение взаимодействия, внешний источник нейтронов, p ( E ) – спектр мгновенных нейтронов деления, (r, E ) – доля запаздывающих нейтронов, f (r, E, t ) – макроскопическое сечение деления нейтронов, – число вторичных нейтронов на акт деления, d, j ( E ) – спектр запаздывающих нейтронов в j-ой группе запаздывающих нейтронов, j – постоянная распада предшественников запаздывающих нейтронов группы j, C j (r, t ) – концентрация предшественников запаздывающих нейтронов группы j.
Запишем данные уравнения с использованием группового энергетического приближения и транспортного приближения при учете анизотропии рассеяния:
где v inv – диагональная матрица, элементами которой являются обратные скорости, операторы L и K имеют следующий вид:
K(r,, t ) = Для обеспечения единственности решения системы уравнений (1.3) и (1.4) дополнительно ставится граничное условие на границе расчетной среды VS и начальное условие Далее в работе будет рассматриваться система уравнений (1.3–1.6) описывающая пространственно-временную кинетику в ядерном реакторе.
Как уже было отмечено, данная глава посвящена анализу методов и программ, используемых для решения задачи (1.3–1.6). Так как решению проводится обзор методов решения стационарной формы уравнения переноса нейтронов (1.3). При этом внимание уделяется только методам, принадлежащих, как и МПГ, к классу проекционно-сеточных методов [28, 29], что позволяет рассмотреть метод поверхностных гармоник в совокупности с методами из этого класса. Приводится описание основ МПГ и его текущей реализации. Приводится рассматриваются методы решения нестационарного уравнения переноса и наиболее известные современные нейтронно-физические программы, проводящие пространственно-временные расчеты переходных процессов в ядерном реакторе.
1.2 Обзор проекционно-сеточных методов решения Одна из основных идей МПГ, а именно поиск решения краевой задачи в виде линейной комбинации пробных функций и некоторых коэффициентов, является широко распространенным подходом и реализуется в классе методов, называемых проекционными (projection methods) [30, 31]. Для описания сути этих методов будим их рассматривать на примере односкоростного диффузионного уравнения с граничным условием где D – коэффициент диффузии нейтронов, Vs – граница расчетной области V.
Идея проекционных методов заключается в том, что точное решение (r ) краевой задачи (1.7), (1.8) ищется в виде функции (r) удовлетворяющей граничным условиям (1.8) и являющейся линейной комбинацией пробных (базисных) функций:
где ci – неизвестные коэффициенты, i (r ) (i = 1,2,..., I ) – пробные функции или независимыми [32]. Данный подход поиска решения краевых задач является альтернативным подходом [32] по отношению к конечно-разностному методу (МКР) [33]. Необходимо отметить, что характер получаемого решения в определенной степени определяется характером пробных функций.
В качестве преимущества проекционных методов можно отметить, что в проекционных методах решение находится сразу во всей области изменения независимых переменных, а не в отдельных точках.
Проекционные методы отличаются по способу поиска коэффициентов cn в (1.9). Их подразделяют на вариационные методы (например, метод Релея-Ритца (Rayleigh-Ritz method) [32]) и методы взвешенных невязок (method of weighted residuals, MWR). К последним относят, например, метод коллокаций (collocation method) [32, 34] и метод Галеркина (Galerkin weighted residual method) [32, 34].
Проекционно-сеточные методы, о которых далее пойдет речь, являются своего рода синтезом проекционных и разностных методов. Их преимущество заключается в том, что к плюсам проекционных методов добавляется то, что эти алгоритмы приводят к «системам уравнений, подобным возникающим в разностных методах (т.е. незначительное число элементов матриц этих систем были бы ненулевыми)» [28, с.17]. Для получения данных методов, достаточно в проекционных методах в качестве базисных функций i (r ) (i = 1,2,..., I ) брать финитные функции, т.е. функции, которые отличны от нуля лишь на небольшой части расчетной области решаемой задачи, что приводит к системе некоторых разностных уравнений, близкой к системе, возникающей в разностном методе [28].
Интерес к рассмотрению этих методов в контексте данной работы вызван тем, что вышеописанное свойство проекционно-сеточных методов, является одним из базовых моментов МПГ. Здесь кратко рассмотрим проекционносеточные методы, более полно нашедшие применение в нейтронно-физических расчетах ядерных реакторов, а именно метод конечных элементов [32], метод граничных элементов [35], метод конечных суперэлементов Федоренко [36] и метод матриц отклика [37].
В методах Релея-Ритца и Галеркина пробные функции i (r ) обычно являются полиномами и чтобы повысить точность расчета необходимо увеличивать степень этих полиномов, что в свою очередь ведет к повышению сложности алгоритма. Как альтернативу этому можно рассматривать возможность разделения всей расчетной области на подобласти (конечные элементы), что позволяет повысить точность расчета при сохранении порядка пробных функций i (r ), благодаря учету свойств элементов по отдельности. Т.е. «сплошная среда разбивается на ряд элементов, которые можно рассматривать как конкретные ее части» [35, с.9]. Это является фундаментальной идеей метода конечных элементов (МКЭ, finite element method, FEM) [32].
Рассмотрим основные моменты данного метода. В МКЭ расчетная область V разбивается на элементы Vn. В результате глобальное решение (r) представляет собой сумму решений по всем элементам:
Внутри n-го элемента функция n (r) ищется в виде следующей линейной комбинации В литературе функции n,i (r ) часто называют функциями формы (shape functions) [32]. Функции формы равны нулю за пределами рассматриваемого элемента и должны быть линейно независимыми. Точки пересечения границ элементов называют узлами.
Отметим, что МКЭ «может основываться как на вариационных принципах, так и на более общих выражениях метода взвешенных невязок» [35, с.9]. Таким образом, для получения уравнений МКЭ используется как метод Релея-Ритца, так и метод Галеркина, основанные на идее наличия некого функционала J который зависит от коэффициентов cn,i [32]. В случае МКЭ вместо всей расчетной области c функционалом J(ci) рассматриваются элементы с функционалами Jn(сn,i) которые в сумме равны J. Первым способом для получения уравнений МКЭ является метод Релея-Ритца при котором, для каждого m-го узла получают узловое уравнение (nodal equation), которое имеет вид:
где n – номер элемента, содержащего m-ый узел. Все узловые уравнения собираются в одно матричное уравнение, которое называется уравнением системы (system equation).
Другой способ для получения уравнений МКЭ – метод Галеркина. В этом случае общий интеграл системы где W j (r ) – j-ая весовая функция, R(r ) – невязка, получаемая при подстановке в решаемое уравнение приближенного решения (1.10). В качестве весовых функций w j (x) выбираются функции формы n,i (r ) [32] и получают систему уравнений для n-го элемента (element equations):
Уравнения (1.14) называются уравнениями элементов (element equations). В качестве альтернативы также возможно получение узловых уравнений. Можно соответствующих уравнений элементов. Необходимо отметить, что в двух- и трехмерных случаях использование уравнений элемента существенно проще, чем узловых уравнений [32].
Важно отметить, что МКЭ позволяет работать с областями произвольной формы. Особенностью МКЭ является то, что в нем можно использовать неравномерные сетки. Такие сетки могут быть статичными, т.е. неизменными в течение всего расчета, а могут быть динамичными, т.е. изменяться по мере расчета, отслеживая особенности результатов. Это делает МКЭ очень гибкой технологией для систем с произвольной геометрией.
Если рассматривать применение МКЭ для уравнения (1.7) с применением метода Галеркина, то отправной точкой является ослабленная формулировка (weak form) [38] уравнения (1.7) которую можно получить, применяя первую формулу Грина [39]:
где V – некоторый объем, ограниченный достаточно гладкой поверхностью S, производными в внутри V + S и имеющие непрерывные вторые производные внутри S, n – внешняя нормаль к поверхности S. Таким образом, имеем следующую ослабленную формулировку уравнения (1.7):
Следует также отметить, что, как правило, в роли функций формы выступают полиномы [40].
Применение метода взвешенных невязок Галеркина накладывает различные требования на непрерывность функции (r ) и на ее производные. Применение данного метода «напрямую» к уравнению (1.7) требует, чтобы функция (r ) была интегрируема с квадратом второй производной, т.е.
При этом весовой функции достаточно быть лишь интегрируемой с квадратом, т.е.
непрерывности искомой функции, что в данном случае делается посредством применения первой формулы Грина. В результате имеем, что функции w j (r ) и (r ) должны быть непрерывными вместе со своими первыми производными.
Одинаковое требование непрерывности к функциям w j (r ), (r ) позволяет использовать одни и те же базисные функции при их разложении, что приводит к симметричным матрицам [35].
Первое применение метода конечных элементов к теории диффузии нейтронов, как отмечается в [38], датируется 1970-и годами. Эти разработки в области применения МКЭ к диффузии нейтронов и в частности к уравнениям переноса были обобщены и описаны в [41–43]. МКЭ реализован в ряде программ нейтронно-физического расчета и здесь отметим лишь некоторые из них:
NEDPCM [38], FTRAN, FTRAN2 [44, 45], EVENT [46].
В отличие от МКЭ, который использует дискретное представление, как самой области, так и ее границы, метод граничных элементов (МГЭ, boundary element method, BEM) основывается на дискретном представлении лишь внешней границы [35]. Если в случае МКЭ используется ослабленная формулировка, снижающая порядок производной на единицу по отношению к исходной задаче, то в МГЭ используется обратная формулировка, которая снижает порядок производной искомой функции еще на единицу, но при этом имеется повышение порядка производной весовой функции. Для получения обратной формулировки уравнения (1.7) применим вторую формулу Грина [38] где V, S, u = u (r ), v = v (r ), n – те же что и для (1.15). Зная, что ток нейтронов на поверхности равен получим В качестве весовой функции применяются функция Грина G (r, ) исходного уравнения где (r ) – дельта-функция Дирака.
В МГЭ проводится дискретизация границы S посредством множества узлов, разделяющих ее на отдельные отрезки (граничные элементы). Заменяя весовую функцию w(r ) на фундаментальное решение для узла n-го граничного элемента G (r, n ) в (1.21) и используя (1.22) и свойства функции Дирака, получаем [47, 48] где Таким образом, использование функций Грина в качестве весовых функций позволяет обеспечить обращение в нуль интегралов по области и сводит задачу только к поиску неизвестных на границе области [35]. Так, в интегральном уравнении (1.23) содержатся неизвестные только на границе, т.е. (1.23) является Поскольку в результате BIE включает неизвестные только на границе системы, размерность задачи практически сокращается на единицу. Поток на границе S представляется как где n (r ) – пробная функция n-го граничного элемента, I n – значение потока в узле n-го граничного элемента. Аналогичным образом представляется функция тока j (r ) и после ее подстановки вместе с (1.24) в BIE получают систему относительно значений потоков и токов в узлах [35]. Основным преимуществом МГЭ по сравнению с МКЭ является снижение размерности задачи на единицу.
Следствием этого является сокращение объема расчетов. Первое применение МГЭ к уравнению диффузии нейтронов, как отмечается в [35], датируется серединой 1980-х [49]. В качестве программы, реализующей МГЭ в нейтроннофизических расчетах, можно отметить NEDPCM [35].
1.2.3 Метод конечных суперэлементов Федоренко Метода конечных суперэлементов Федоренко (МКСЭ) был предложен [36, 50] для решения задач теории диффузии, ядерных реакторов, кинетики и т.д.
Метод нацелен для работы с областями, имеющими резкие неоднородности в малых зонах по сравнению со всей рассматриваемой областью. Такие проблемы, как правило, решаются с использованием адаптированных сеток, либо с использованием очень мелких сеток. При этом в первом случае требуются специально разработанные алгоритмы, а во втором случае – большие вычислительные ресурсы. Однако, беря за базовый принцип идею, что резкие неоднородности, как правило, расположены в малых подобластях, становится возможным использование крупных сеток, грани и узлы которых проходят через относительно гладкие части решения [51].
Другой особенностью метода является то, что в методе Федоренко функции n,i (r ) являются решениями исходного уравнения в элементе с определенными граничными условиями на границе этого элемента, в отличие от МКЭ, где базис строится из общих, не связанных с решаемой задачей, соображениях, и функции n,i (r ) записываются в явном виде. Из условий сшивки на координатных линиях сетки, состоящих из требований непрерывности искомой функции и непрерывности нормальных к линии сетки производных, строится система «балансных соотношений». Следует заметить, что функция (r ) непрерывна и почти всюду, за исключением координатной сетки, удовлетворяет решаемому уравнению.
Таким образом, в МКСЭ базис состоит из решений исходной задачи и это определяет его преимущество. МКСЭ претендует на расчет с большим шагом, что вместе с тем определяет его ограниченность. Также следует заметить, что МКСЭ требует предварительного расчета базисов [52].
Первое рассмотрение метода матриц отклика (response matrix method) как отмечается в [53], относится еще к 1963 г. [37]. Метод имеет сходные с методом Федоренко подходы. Рассмотрим основные моменты данного метода на примере уравнения переноса нейтронов в соответствии с работой [53].
Проводится пространственная дискретизация грубой сеткой. На локальном уровне расчета проводится решение уравнения переноса нейтронов в каждом элементе для определенного набора граничных условий для каждого элемента.
Набор граничных условий представляет собой втекающие односторонние токи нейтронов равные некоторому ортогональному базису m, m = 0,1... в фазовом пространстве, определенному для границы n-го элемента. Таким образом, решается уравнение переноса нейтронов в элементе с граничным условием где J ( ns ) – ток падающих нейтронов на поверхность S n-го элемента, – вектор, представляющий фазовое пространство. Результатом локального этапа является набор функций m ( ns ) ( m ( ns ) – функция распределения нейтронов в n-ом элементе при граничном условии m заданном на поверхности элемента S).
Граничное условие (1.25) накладывается на все элементы среды независимо от того являются они граничными или нет.
Рассмотрим уравнение метода матриц отклика, решаемое на глобальном уровне расчета. Пусть J nt + ( nt ) – ток вылетающих частиц с поверхности t, инициированный налетающим током m ( ns ) на поверхность S. Его можно представить в виде разложения по базису m :
где rimt = m ( nt ) J nt + ( nt )d nt. Отличительной чертой метода матриц отклика является то, что связь между током втекающих частиц и током вытекающих частиц можно определить следующим уравнением где R n – матрица отклика (response matrix). Вектор J n представляет собой вектор коэффициентов разложения для всех поверхностей n-го элемента Поскольку вытекающие токи для ячейки являются втекающими токами для смежных ячеек, вводится матрица связи M (connectivity matrix) такая что где J + и J – блочные вектора, состоящие из векторов J n+ и J n соответственно для всех ячеек. Матрица M имеет не более чем один ненулевой элемент на строку и столбец. Таким образом, уравнение метода матриц отклика имеет вид где R – это блочная диагональная матрица, состоящая из блоков R n. Данный метод реализован, например, в программе HELIOS [54].
Б.П. Кочуровым развивался схожий с вышеперечисленными методами подход – метод -матриц [55–57]. В качестве отправной точки для этого метода служит метод источников–стоков [58]. Итоговые уравнения данного метода по своей сути близки к уравнениям МПГ. Эти уравнения были получены только для трех поперечных пробных матриц, что является основным ограничением этого метода.
После рассмотренных выше проекционно-сеточных методов рассмотрим метод поверхностных гармоник, который также можно отнести к данному классу методов и которому посвящена данная работа. Данный метод был предложен профессором Н.И. Лалетиным в 1976 году в работе [6].
1.3.1 Основы метода поверхностных гармоник и определения Приведем здесь кратко описание базовых моментов МПГ и определения, использующиеся в работах по МПГ.
Пространственная декомпозиция. Так же как и в представленных в обзоре проекционно-сеточных методах, в МПГ расчетная область V разбивается на пространственные подобласти Vk (ячейки). В результате глобальное решение (r, ) (групповой вектор) представляет собой сумму решений по всем ячейкам:
где k (r, ) – решение в k-ой ячейке. Тут следует заметить, что МПГ опирается на идею, что активная зона, большинства современных ядерных реакторов может быть представлена в виде набора ячеек одинаковой формы на разных уровнях вложенности, а именно ячеек тепловыделяющих сборок и ячеек твэлов. Такая периодическая структура активной зоны позволяет проводить пространственную декомпозицию естественны образом.
Пробные функции. Групповое решение в каждой ячейке представляется в виде линейной комбинации пробных матриц и амплитуд:
где (kl )(r, ) – l-ая пробная матрица для k-ой ячейки, I (l ) – l-ая амплитуда k-ой ячейки в виде группового энергетического вектора. Отметим, что в дальнейшем будем использовать символ ( L ) (r, ) для решения во всей активной зоне полученной при использовании L + 1 количества пробных матриц.
Остановимся подробнее на пробных матрицах МПГ. Пробные матрицы МПГ состоят из G пробных векторов (kl,)g (r, ) Каждый (kl,)g (r, ) пробный вектор представляет собой решение уравнения переноса нейтронов в k-ой ячейке с граничным условием в виде тока, заданного в виде l-ой координатной функции Wl (rs ) на поверхности ячейки rs в g’-ой энергетической группе (в остальных энергетических группах ток равен нулю).
Определим следующие связанные с пробными функциями матрицы (l ) элементы которых имеют вид [7]:
где где n – внутренняя нормаль к границе ячейки. Таким образом, ( (kl ) ) gg представляет собой деленный на норму N l интеграл по границе области ячейки и по от решения уравнения переноса в g-ой энергетической группе с весом координатной функции Wl (rs ) и множителя 3( n) 2 при граничном условии, распределенном по l-ой координатной функции в g’-ой энергетической группе.
Координатные функции. Приведем описание координатных функций для ячеек с квадратной внешней границей, т.к. дальнейшее получение уравнений проводится непосредственно для квадратной решетки блоков. Координатные функции МПГ Wl (rs ) описываются следующей формулой (например [9]):
где l – номер координатной функции l=0, 1,…, L, Pp ( j ) – полиномы Лежандра p-го порядка, M – число боковых сторон (граней) ячейки (в случае квадратной решетки блоков равно M = 4), a – длина боковой стороны ячейки, j – угол между нормалью, построенной из центра ячейки на j-ю боковую сторону ячейки, и осью абсцисс (см. рисунок 1.1).
Рисунок 1.1 – Параметры координатных функций В таблице 1.1 представлены первые восемь координатных функций для квадратной ячейки [8].
Таблица 1.1 – Первые восемь координатных функций для квадратной ячейки На рисунке 1.2 представлена схема втекания токов для первых восьми координатных функций для ячейки с квадратной границей.
Рисунок 1.2 – Схема втекания токов, соответствующая первым восьми координатным функциям для ячейки с квадратной границей 1.3.2 Обзор программных комплексов SUHAM и SVS В данном параграфе кратко рассмотрим программные комплексы SUHAM и SVS, реализующие стационарные уравнения метода поверхностных гармоник.
SUHAM. «Комплекс программ SUHAM предназначен для реализации конечно-разностных уравнений метода поверхностных гармоник, для расчета нейтронно-физических процессов в ядерных реакторах с треугольной и квадратной решетками блоков (ТВС)» [14, c.10]. «Модули комплекса SUHAM, реализующие конечно-разностные уравнения МПГ, применяются для решения многогруппового уравнения переноса нейтронов во всем объеме активной зоны реактора методом поверхностных гармоник» [14, c.10]. Комплекс состоит из двух независимых частей SUHAM-W [11] и SUHAM-U [12].
SUHAM-W взаимодействует с компонентами программы WIMS-SH (раннее название WIMS-SU) [61] для подготовки групповых сечений изотопов и материалов и использует в качестве библиотеки ядерных данных библиотеку программы WIMS-D [62]. К основным возможностям SUHAM-W по решению двумерного конечно-разностного уравнения МПГ можно отнести:
1. Решение группового уравнения диффузии (прямого и сопряженного) с различными граничными условиями (традиционный метод гомогенизации);
2. Решение уравнений МПГ с разным числом пробных матриц (3, 4, 7 и 8 для квадратной решетки и 3, 4, 5 и 6 для треугольной решетки);
3. Решение сопряженных уравнений МПГ с тремя пробными матрицами.
SUHAM-U взаимодействует с компонентами комплекса UNK [63] для подготовки групповых сечений изотопов и материалов и решения уравнений изотопной кинетики, использует в качестве библиотеки ядерных данных библиотеку программы UNK. Комплекс SUHAM-U позволяет решать двумерные и трехмерные нейтронно-физические задачи. Двумерный модуль имеет следующие возможности:
1. Решение группового уравнения переноса в двухэтапном расчете топливных сборок с помощью МПГ с разным числом пробных матриц: от 3-х до 8-и для квадратной решетки и от 3 до 6 для треугольной решетки.
2. Трехэтапный расчет активной зоны ядерного реактора с шестигранными 3. Расчет выгорания ТВС ВВЭР-1000.
Трехмерный модуль решает трехмерные конечно-разностные уравнения МПГ с тремя поперечными и двумя продольными пробными матрицами на каждую ячейку для реакторов с квадратной решеткой блоков.
Верификация программного комплекса SUHAM проводилась на следующих задачах: ТВС PWR [64], ТВС ВВЭР-1000 с урановым и MOX топливом [65], бенчмарк C5G7 [66], бенчмарк VENUS-2 [67], полномасштабная двумерная зона ВВЭР [16], полиячейки и модельные сборки РБМК [14], двумерная зона БРЕСТ-ОД-300 [14], реактор ГТ-МГР [15].
SVS. Комплекс программ SVS (Surface Values System) и его составные части в виде самостоятельных модулей SVL [13, 68] и SVC [13] «используются в исследованиях по физике ВВЭР и РБМК и позволяют оценивать вклад от различных приближений, используемых в проектных исследованиях и эксплуатационных расчетах» [1, c.8]. SVS позволяет проводить исследования как локального нейтронного поля, так и всей активной зоны целиком.
SVL (Surface Values Lattice) представляет собой решеточный код который «предназначен для решения задач теории переноса нейтронов в ячейках и полиячейках для основных типов решеток ядерных реакторов» [1, c.12]. Данная программа разрабатывалась в основном для реакторов легководного типа, но также применима для реакторов графитового и тяжеловодного типов. SVL использует МППИ для расчета характеристик ячеек и МПГ для расчета кассет (полиячеек). Программа SVL позволяет проводить расчеты на критичность и расчеты на выгорание ячеек в цилиндрической и кластерной геометриях, 2-D и 3-D сборок для гексагональной, квадратной и треугольной решеток блоков. Также следует отметить что программа позволяет работать со сборками, имеющими нарушения регулярности решетки.
SVC (Surface Values Core) «предназначена для решения полномасштабных нейтронно-физических расчетов активных зон реакторов ВВЭР» [1, c.17] используя МПГ. Входные данные для SVC готовит программа SVL. В обеих программах используется последовательный подход, т.е. «используются методы, в которых начальное приближение удовлетворяет требованиям проектных расчетов, а результаты расчетов получаются по точности сравнимые с реперными, что можно проверить, переходя к следующим приближениям» [1, с.13].
Общая схема расчета по SVS выглядит следующим образом: расчет характеристик кассет, расчет активной зоны, восстановление потвэльных полей, учет обратных связей [1]. На рисунке 1.3 представлена блок схема комплекса SVS [13].
Комплекс прошел множественную верификацию, в том числе на ячейках реактора BBЭР-1000 [69], сборках реакторов типа ВВЭР (в частности ZR6) [13, 70]. Также следует отметить, что проводился расчет кампании реактора ВВЭР-1000 первого блока Волгодонской АЭС [13].
В дополнение к рассмотренным программам, следует отметить, что имеется также гетерогенная версия программы STEPAN, использующая МПГ [2, 71].
1.3.3 Обзор работ по нестационарным уравнениям МПГ В данном параграфе рассмотрим основные работы [72–77] в которых проводились попытки получения конечно-разностных уравнений метода поверхностных гармоник для нестационарного уравнения переноса нейтронов.
Первая такая попытка была сделана в 1990 г. Н.В. Султановым и В.Ф.
Бояриновым. Суть этой работы заключается в том, что члены уравнения переноса, связанные с запаздывающими нейтронами и с производной по времени были перенесены в правую часть и представлены в виде источников. Введены пробные функции от объемных источников (производная по времени и запаздывающие нейтроны). По аналогии со стационарным случаем, невязка функционала возникала только на границах ячеек. Как итог были получены конечноразностные уравнения МПГ для пространственной кинетики в 1-D (плоской), 2-D (квадратной и треугольной) и 3-D геометриях. Дальнейшего развития этот подход не получил и опубликован не был.
В 1993 г. в работе [72, 73] был предложен другой подход. Главной особенностью этого подхода является то, что при построении конечноразностных уравнений МПГ используются только пробные функции стационарного уравнения переноса нейтронов. Как итог, помимо поверхностной невязки на границах ячеек возникает еще и объемная невязка, которая связанна с производной функции распределения нейтронов по времени и с запаздывающими нейтронами. Конечно-разностные уравнения получаются из минимизации суммарной невязки.
пространственной кинетики метода поверхностных гармоник с помощью совместного рассмотрения прямой и сопряженной задач в одинаковом приближении МПГ.
В работе [77] представлены полученные еще в 90-х годах конечноразностные уравнения пространственной кинетики метода поверхностных гармоник для плоской одномерной решетки.
Следует отметить, что ни в одной из этих работ полученные уравнения не были программно реализованы и верифицированы. Первой работой, в которой проведена программная реализация нестационарных уравнений МПГ, является [78]. Программная реализация проведена для конечно-разностных уравнений пространственной кинетики в одномерной плоской геометрии.
Следует также отметить, что в [2] без вывода приведены двухгрупповые 11-точечные по пространству нестационарные уравнения МПГ. Реализованные в программе STEPAN эти уравнения позволили повысить надежность расчета ряда нейтронно-физических функционалов.
1.4 Обзор методов решения нестационарного представлено в виде где (t ) - первая производная по времени функции распределения нейтронов, H(t ) - полный оператор переноса в матричной форме. В уравнение (1.35) все независимые аргументы, кроме времени, для простоты опущены. Имеется ряд подходов для работы с временной зависимостью в уравнение (1.35), которые подразделяются на две категории.
Первая категория – это прямые методы, использующие конечно-разностную дискретизацию по времени, в которых временная ось представляется в виде набора небольших временных промежутков, на каждом из которых константы постоянны. Прямые конечно-разностные методы являются наиболее простыми методами для решения пространственно-временных задач. Эти методы делают единственное предположение, что дифференциальные операторы адекватно описываются конечными разностями [79]. Как следствие отсутствия ряда серьезных приближений, эти методы характеризуются наличием большого числа неизвестных в системе уравнений. По этой причине они требуют большого объема вычислительных затрат. К таким методам относят, например, полностью явный и полностью неявный методы.
Вторая категория – это непрямые методы использующие процедуру разделения переменных. Среди этих методов можно назвать улучшенный квазистатический метод [27] и SCM [80] методы. Далее кратко рассмотрим нестационарные методы решения переноса нейтронов.
При использовании полностью явного метода Эйлера (fully explicit method) уравнение (1.35) принимает вид и решение представляется в виде Данный подход требует минимума вычислительных затрат на каждый временной шаг по сравнению с другими конечно-разностными методами, однако он неустойчив если шаг t недостаточно мал. Данный факт требует такого уменьшения временного шага, что преимущество, даваемое простотой этого метода, как правило, компенсируется большим количеством временных шагов [80].
Дополнительно отметим, что В.И. Лебедевым был разработан вариант устойчивой разностной явной схемы с переменными шагами по времени реализованный в программе DUMKA [81].
Проблема устойчивости явного метода может быть решена посредством использования неявного метода. При использовании неявного метода Эйлера (fully implicit method) уравнение (1.35) принимает вид и решение представляется в виде Данный подход является устойчивым, однако требует обращения матрицы I tH (t n+1 ) на каждом временном шаге, что влечет большие вычислительные затраты.
Наиболее известным среди полунеявных методов является -метод (-method) [80]. В -методе уравнение (1.35) представляется в виде где – диагональная матрица, элементы которой ii удовлетворяют неравенству 0 ii 1. При = 0 метод переходит в полностью явный метод, при = I – полностью неявный метод, а при ii = 0,5 (для всех i) переходит в схему КранкаНиколсон (Crank-Nicholson). Решение в -методе имеет вид Матрица выбирается так чтобы повысить точность конечно-разностной аппроксимации и может изменяться от шага к шагу.
Рассмотрим суть метода переменных направлений (МПН) [33]. Метод был предложен D.W. Peaceman, H.H. Rachford [82] и J. Douglas [83]. Идея МПН заключается в том, чтобы в одном из направлений использовать полностью неявную схему (в остальных направлениях полностью явную схему) и изменять это направление от временной точки к временной точке [80]. Для простоты будем рассматривать двумерный случай. Матрица H разделяется по следующему принципу где X – симметричная трехдиагональная матрица, отвечающая за перенос в одном направлении, Y – симметричная трехдиагональная матрица, отвечающая за перенос в перпендикулярном первому направлении [84].
В методе переменных направлений конечно-разностная схема имеет вид а, значит, решение представляется в виде Разложение матрицы H на две матрицы, отвечающие за перенос в различных направлениях, производится с целью попеременного применения явной и неявной схем для различных направлений в соответствии со схемой Писмена-Рекфорда [82]. Матрицы R i выбираются таким образом, чтобы сделать процесс обращения матриц в системе (1.44) наименее трудоемким [84].
Данный метод является довольно эффективным при решении временных задач, однако, он требует введения конечно-разностной аппроксимации расчетной пространственной области, что делает его несовместимым, например, с методом характеристик [4].
Данный метод вместе с -методом был применен в программе TWIGL для решения двумерного группового уравнения диффузии [19].
Основная идея улучшенного квазистатического метода (improved quasistatic method, IQS) заключается в том, что выделяется две функции, имеющие различную скорость изменения во временном масштабе: быстро меняющаяся амплитуда (амплитудный фактор) и форм-функция, которая обычно изменяется медленнее. В этом методе нестационарная функция распределения (r, E,, t ) представляется в виде произведения форм-функции (r, E,, t ), которая слабо изменяется с течением времени, и амплитуды T (t ), которая зависит только от времени и с течением времени, обычно, изменяется быстрее [85]:
Это представление функции распределения дополняется следующим условием нормировки:
где + (r, E, ) – функция, сопряженная стационарной функции (r, E, ).
Уравнение для расчета форм-функции (уравнение формы) получается после подстановки (1.45) в нестационарное уравнение переноса нейтронов (1.1):
После ряда преобразований нестационарного уравнения переноса нейтронов с использованием формулы (1.45) и сопряженного стационарного уравнения переноса нейтронов получают уравнение для расчета амплитуды:
где функционалы (t ), (t ), (t ) и ci (t ) определяются через функции (r, E,, t ) и + (r, E,, t ) [27, 85]. Уравнение для ci (t ) имеет вид Обычно форм-функция изменяется более медленно во времени, чем амплитуда и, следовательно, форм-функция рассчитывается с более крупными временными шагами, чем амплитуда. Результат решения уравнения формы используется для расчета параметров точечной кинетики, которые используются при решении уравнения (1.48).
(r,, E, t ) в уравнении (1.47), то придем к квазистатическому приближению, использующему факт слабой зависимости форм-функции (r, E,, t ) от времени.
Жесткость системы пространственной кинетики определяется временем жизни мгновенных нейтронов, которое на несколько порядков меньше времени жизни запаздывающих нейтронов [80]. Метод SCM (stiffness confinement method) был разработан [86] для уменьшения ограничений на временной шаг, и, следовательно, повышения производительности вычислений. В данном методе рассматриваются отдельно посредством введения динамических частот (dynamic frequencies) Используя эти определения, производится замена производных в исходной поглощения и деления, которые включают динамические частоты. Для того чтобы получить решение во времени, делается оценка динамических частот, решается модифицированное уравнение переноса, пересчитываются концентрации предшественников. Такие итерации проводятся до сходимости [86]. Данный метод реализован в таких кодах как SPNOVA [87] и PANTHER [88].
В дополнении к рассмотренным методам необходимо отметить класс линейных многошаговых методов (linear multistep methods) [32, 34]. В таких методах значение решения в n-ой временной точке определяется значениями в ряде предыдущих точек. В качестве примеров методов из этого класса можно назвать методы Адамса-Башфорта (Adams-Bashforth), Адамса-Мултона (AdamsMoulton) и др. Недостатком этих методов является возрастание величины необходимой памяти по сравнению с одношаговыми методами для хранения моделировании полномасштабных моделей. Также отметим наличие класса методов типа прогноз-коррекция (predictor-corrector methods) [32] обеспечивающие контроль ошибки на каждом шаге расчета, что также является важной чертой современного кода.
1.5 Обзор программ для решения нестационарного Как уже отмечалось ключевой задачей для анализа нейтронных переходных процессов в ядерном реакторе является решение нестационарного уравнения переноса нейтронов с запаздывающими нейтронами. Исторически, из-за ограниченных возможностей вычислительных технологий на ранних этапах атомной науки, нейтронные переходные процессы описывались в значительной степени упрощенно, с использованием 1D подходов, либо с использованием уравнений точечной кинетики [3, 89].
Со значительным ростом вычислительных возможностей компьютерной техники в течение второй половины XX века и развитием вычислительных пространственно-временного группового уравнения переноса нейтронов с запаздывающими нейтронами в диффузионном приближении. Среди таких кодов наиболее известными иностранными кодами являются TWIGL [19], LUMAC [90], MATIN [91], SADI [84], DIF3D-K [92] и множество других. Среди российских кодов, проводящих расчеты ядерных реакторов и использующих данные приближения можно отметить STEPAN [2], БИПР-8КН [93], JAR-IQS [94], ShIPR [95], ГЕФЕСТ [96] и др. Коды данного поколения содержат, в том или ином пространственная гомогенизация и расчет в малом числе энергетических групп.
Несмотря на то, что качество расчетов, проводимых этими кодами, часто является достаточным для анализа многих сценариев поведения различных типов реакторов, ряд исследователей (например, [3–5]) указывает на возможность получения улучшенных результатов путем решения более универсального газокинетического уравнения, особенно для анализа безопасности.
пространственно-временного уравнения переноса нейтронов, были программы TIMEX [97] и TRANZIT [98] для 1-D и 2-D геометрий в начале 70-х годов. На текущий момент уравнение переноса нейтронов, как правило, применяется в основном в стационарном расчете, по причине высоких вычислительных затрат для нестационарных задач даже на современных компьютерах. Вместе с тем, растущая необходимость применения уравнения переноса нейтронов в расчете переходных процессов, особенно ввиду строгих стандартов безопасности, предъявляемых к новым разрабатываемым реакторам, привела к созданию ряда современных программных комплексов, позволяющих решать пространственновременное уравнение переноса нейтронов с запаздывающими нейтронами. Далее рассмотрим кратко наиболее известные из них.
DORT-TD. Код DORT-TD [3, 4] представляет собой нестационарную расширенную версию хорошо известного кода DORT [99]. DORT является двухмерным кодом (также есть одномерный плоский вариант), разработанным для x-y, r-z или r- геометрий. Он может быть использован для решения прямой или сопряженной форм уравнения переноса нейтронов. Решается либо диффузионное уравнение, либо уравнение переноса нейтронов методом дискретных ординат.
TORT-TD. Трехмерный нестационарный код TORT-TD [100, 101] проводит расчеты с помощью метода дискретных ординат. Данный код основывается на стационарном коде TORT [102, 103]. TORT-TD решает нестационарное групповое уравнение переноса нейтронов в декартовых координатах и цилиндрической геометрии. Стабильность временной схемы достигается посредством использования полностью неявной схемы.
нестационарного уравнения переноса нейтронов в декартовых координатах методом дискретных ординат. PARTISN использует полунеявную временную схему.
EVENT. Код EVENT [45] использует методы конечных элементов и сферических гармоник для решения нестационарного уравнения переноса нейтронов.
BARS. Программа BARS [105] использует метод -матриц Кочурова и улучшенный квазистационарный метод в нестационарных расчетах.
В дополнение к этим детерминистическим кодам, следует отметить, наиболее известные программы по решению нестационарного уравнения переноса нейтронов методом Монте-Карло, а именно TDKENO [106] и TRIPOLI [107].
В результате выполненного в данной главе обзора литературы по проекционно-сеточным методам (методы конечных элементов, граничных элементов, конечных суперэлементов Федоренко, матриц отклика), применяемых при проведении нейтронно-физических расчетов ядерных реакторов, можно сделать вывод, что проекционно-сеточные методы нашли широкое применение в нейтронно-физических расчетах ядерных реакторов. Из всех рассмотренных проекционно-сеточных методов особое внимание уделено методу поверхностных гармоник, выбранному в качестве метода для разработки нового кода для проведения пространственно-временных расчетов.
Выбор в пользу МПГ на фоне остальных проекционно-сеточных методов сделан на основе факта, что метод поверхностных гармоник является одним из наиболее разработанных методов для решения стационарного уравнения переноса нейтронов. МПГ занимает промежуточное место между детерменистическими и инженерными методами и позволяет получать основные нейтронно-физические функционалы с точностью первых при вычислительных затратах сравнимых со вторыми. Метод поверхностных гармоник предоставляет возможность последовательного уточнения решения искомой задачи с увеличением номера МПГ-приближения. Вместе с тем, наиболее важным фактом для выбора МПГ является то, что метод поверхностных гармоник разрабатывался специально для решения нейтронных реакторных задач и, как следствие, учитывает физические особенности реакторов при построении уравнений. В качестве примера такого учета следует отметить следующие особенности МПГ:
• упорядоченность координатных функций по степени их влияния на нейтронно-физические функционалы, что позволяет достигать приемлемой приближениях метода;
• учет гладкости функции потока нейтронов на границе ячеек при выводе конечно-разностных уравнений.
В дополнение к вышеперечисленному имеется многолетний опыт использования МПГ в НИЦ «Курчатовский институт» с реализацией в программных комплексах идеализированных бенчмарках (например C5G7), так на реакторах ВВЭР-1000, РБМК-1000, БРЕСТ-ОД-300, ГТ-МГР и БН-1200. При этом МПГ демонстрировал детерминистическими методами при сохранении качества их результатов.
нестационарном случае.
В данной главе отмечается, что начало работы над нестационарными уравнениями МПГ было положено еще в начале 90-х годов с публикациями ряда работ на эту тему различными авторами. Вместе с тем, следует отметить, что ни в одной из этих работ полученные уравнения не были программно реализованы и верифицированы. Также следует отметить, что сокращенное изложение выводов данных уравнений в этих работах, не позволяет воспользоваться этими результатами напрямую.
Также в данной главе рассмотрены наиболее известные методы решения нестационарного уравнения переноса нейтронов, а именно: полностью явная и неявная схемы, -метод, метод переменных направлений, улучшенный квазистационарный метод и SCM-метод. Здесь необходимо отметить, что безусловно устойчивым и при этом одним из наиболее простых в реализации подходов для решения нестационарного уравнения переноса нейтронов является полностью неявная схема Эйлера. Данный подход выбран для дальнейшего использования при построении нестационарных алгоритмов на основе метода поверхностных гармоник.
При анализе работ по современным нестационарным нейтроннофизическим кодам, видно, что одним из важных требований, предъявляемых к кинетической части таких кодов, является возможность проведения расчета без использования диффузионного приближения. Также важно отметить, что количество таких программ довольно ограничено. В этом свете в качестве положительной черты метода поверхностных гармоник можно отметить инвариантность конечно-разностных уравнений МПГ по отношению к исходным уравнениям, из которых они получены, а именно, к уравнению переноса нейтронов и диффузионному уравнению, а также отсутствие процедуры пространственной гомогенизации.
Как итог данной главы, можно заключить, что метод поверхностных гармоник является подходящим и перспективным методом для применения его в разработке современного нейтронно-физического пространственно-временного кода.
Глава 2 Нестационарные уравнения метода поверхностных гармоник В данной главе проводится детальный вывод двумерных конечноразностных нестационарных уравнения метода поверхностных гармоник для трех и четырех пробных матриц в квадратной решетке блоков. Получение данных уравнений проводится традиционным для МПГ методом минимизации невязки.
Также основываясь на полностью неявной схеме Эйлера, представлено получение итерационной схемы решения уравнений МПГ для трех пробных матриц.
Решение (r,, t ) уравнений (1.3) и (1.4) будем представлять в виде функции ( L ) (r,, t ) т.е. ( L ) (r,, t ) удовлетворяет уравнению (1.3) с некоторой невязкой которая является точным решением уравнения, сопряженного стационарной форме уравнения (1.3) во всем объеме V рассчитываемого объекта, а именно где Умножим скалярно уравнение (2.2) на + (r, ), уравнение (2.3) на ( L ) (r,, t ) и вычтем второе из первого. Результат проинтегрируем по всему ( L ) (r,, t ) точно не удовлетворяет условиям непрерывности на границах ячеек и граничным условиям на внешней границе реактора, то равенство бесконечность на внутренних поверхностях, поскольку на этих поверхностях решение ( L ) (r,, t ) имеет разрывы [27].
где Vk – объем k-ой ячейки, (a; b ) – скалярное произведение векторов a и b. Далее [6–10], как N. Отметим, что величина N представляет собой сумму двух компонентов, связанных с поверхностной невязкой и объемной невязкой Далее остановимся подробнее на этих величинах.
Величина N s (2.6) связана с поверхностной невязкой, т.е. невязкой, обусловленной тем фактом, что ( L ) (r,, t ) точно не удовлетворяет уравнению переноса нейтронов на границах ячеек из-за ограниченного числа пробных матриц. Значения данной величины приведено, например, в работе [9] и в данном параграфе приводится вывод данного представления только для полноты описания уравнений метода поверхностных гармоник. Преобразуем (2.6) используя теорему Остроградского-Гаусса [108] В последнем члене в нижних индексах номер k’ – номер ячейки соседней с k-ой ячейкой, для которой грань j’ совпадает с гранью j ячейки k, при этом j’ = j + 2, для j = 1, 2 и j’ = j – 2, для j = 3, 4. Знак минус перед вторым членом появляется из-за перегруппировки членов невязки.
Разложим прямую (kL ) (r,, t ) и сопряженную + ( L ) (r, ) функции по формулам Здесь следует отметить, что в формуле (2.9) опущена временная зависимость у пробных функций (kl ) (r, ), т.к. на каждом временном интервале ts в данном подходе они являются решениями стационарного уравнения переноса нейтронов, соответствующего некоторой временной точке на интервале ts (см. далее уравнение (2.87).
Разложение функций (kl ) (r j, ), k т ( l ) (r j, ) по полиномам Лежандра на границе ячеек имеет вид где Pn ( µkj ) – полином Лежандра n-го порядка, µ kj = ( ; n kj ), разложением (2.12), оставляя только первые два члена [7], проинтегрируем по и учтем (2.13), получим (2.15) Далее введем обозначения для средних по грани ячейки тока и уровня нейтронов [9] где S – длина грани ячейки. Также отметим, что для первых четырех сопряженных пробных матриц можно записать Подставим выражения (2.16) и (2.17) в (2.15) и учтем (2.18). Отметим, что знак минус при токе jk j появляется из-за обратной нормали, а в уровне fk j знак сохраняется из-за наличия квадрата (n kj ; ) 2. Таким образом, выражение (2.15) принимает вид:
Здесь использован факт, что Wl (r j ) постоянен на каждой грани ячейки для l = 0, …, 3. Далее отметим, что = Wl (r j ) Wl (r j )I (kl ) (t ) Тогда (2.19) примет вид где Величина N v связана с объемной невязкой. Данная невязка возникает в нестационарном случае и обусловлена тем фактом, что пробные матрицы рассчитываются посредством решения стационарного уравнения переноса нейтронов. Как уже отмечалось, величина N v имеет вид:
Разложим ( L ) (r,, t ) в соответствии с формулой (2.9), а функцию k ( L ) (r,, t ) по формуле (2.10).
Отметим, что в случае объемной невязки разложение матриц (kl ) (r, ), + (r, ) в ряд по полиномам Лежандра производится в объеме всей ячейке:
где Pn (µ) - полином Лежандра n-го порядка, Ограничим сумму ряда (2.23) до двух членов [7] и подставим формулы (2.9), (2.10) и (2.23) в (2.7) Далее отметим, что полиномы Лежандра ортогональны [108], т.е.
где nm - символ Кронекера. Также отметим, что (2.26), заменяя (kl ) (r, ) по формуле (2.22).
Проинтегрируем второй член Таким образом, при L = L объемная невязка равна Далее, при условии получим Как было показано в предыдущих двух параграфах величину (2.5) можно записать как Получение уравнений МПГ будем проводить традиционным для МПГ способом – методом минимизации невязки. Приравняем нулю коэффициенты при I k ( l ) +,(nl ). Для этого в объемной части невязки вставим выражения +,т0( l ) +,т0(l ) +,т1( l ) ( k,т ( l ) ) после вектора I k ( l ) (t ). В результате получим при l = 0 получим Для l = Для l = Для l = В уравнениях (2.31) – (2.38) использованы следующие обозначения где Умножим уравнение (2.31) на (k1) и вычтем его из уравнения (2.32), получим где Умножим уравнение (2.34) на (k1) и вычтем его из уравнения (2.35), получим где Умножим уравнение (2.35) на (k1), вычтем его из уравнения (2.36) и сделаем предположение о симметрии ячеек относительно поворота на угол 90°, в результате чего (k2) = (k1). Тогда где Умножим уравнение (2.37) на (k3) и вычтем его из уравнения (2.38), получим где В формулах (2.42), (2.44), (2.46) и (2.48) использованы следующие обозначения Запишем средние по грани ток jkj и уровень (поток) f kj и для граничной ячейки Умножая выражение (2.51) на (k1), вычитая его из (2.52) и применяя обозначения (2.50) получим Умножая выражение (2.53) на (k), вычитая его из (2.54) и применяя обозначения (2.50) получим Подставим выражение (2.56) в равенство (2.42), получим Подставим выражение (2.56) в равенство (2.44), получим Подставим выражение (2.56) в равенство (2.46), получим Подставим выражение (2.56) в равенство (2.48), получим представить в виде Следует отметить, что в стационарном случае вектор k, j отвечает за влияние высших гармоник и в случае, когда число сторон равно числу пробных матриц, он равен нулевому вектору. Определим k, j в нестационарном случае. Для этого подставим выражение (2.61) в (2.57) – (2.60) и получим Видно, что в нестационарном случае и в случае совпадения числа пробных матриц и сторон ячейки. Однако, можно сделать предположение о значительно меньшей степени влияния членов, стоящих в правой части равенств (2.62), по сравнению с другими членами в равенстве (2.61) и отбросить их при дальнейшем выводе уравнений. Таким образом, далее будем считать, что С учетом выражения для токов (2.63) и того, что для выделенных углов 0°, 90°, 180°, 270° (соответственно j = 1, 2, 3, 4) справедливо уравнения (2.31), (2.33), (2.35), (2.37) можно записать в виде Представим матрицу ( (k1) + (k1) ) в виде где h – шаг решетки (для квадратной решетки совпадает с S). Тогда первое уравнение системы (2.65) имеет вид Введем обозначения и заметим, что величина Также введем обозначение Fv(,lk),( 0 ) = Fv(,lk),( 0) kl),l = 0;3 (если l=0, то l’=0, если l=3, то l’=1).
Посредством формул (2.68) и (2.69) преобразуем первые два члена в первом уравнении системы (2.65), получим где вектор S (k0) (t ) играет роль источника и имеет вид, 0 и 1 – конечно-разностные операторы, такие, что где A k – произвольный вектор.
Аналогичным образом преобразуем четвертое уравнение системы (2.65) и получим где Перейдем теперь к уравнению для предшественников запаздывающих нейтронов (1.4). Для наглядности запишем его еще раз где проинтегрируем по объему ячейки Учитывая формулу (1.31) получим Тогда, считая в (2.75), что gj (r) и gf (r, t ) симметричны относительно центра ячейки где f (r ) – групповой вектор-строка, каждый элемент которого определяется как и применяя обозначения из (2.50) для Ф k и (2.68) получим или в более краткой форме где Таким образом, имеем следующие системы уравнений для случая четырех пробных матриц Отметим, что v (inv),k, v (inv),k – недиагональные матрицы v (inv,k = Fv(,lk),( 0 ). Для случая трех пробных матриц уравнения имеют вид На первом этапе данной работы проводился вывод нестационарных уравнений МПГ для одномерной плоской геометрии. Данная геометрия интересна тем, что в плоской одномерной ячейке возможны только две пробные матрицы – симметричная (k0) (r, ) и антисимметричная (k1) (r, ). Таким образом, набор пробных групповых функций образует полную систему решений. Как результат, в этой геометрии проще исследовать основные особенности применения МПГ для нестационарного случая, не отвлекаясь на возможную недостаточность пробных решений.
В данном параграфе приведем лишь результаты, т.к. процедура получения уравнений аналогична двумерному случаю. Конечно-разностные уравнения МПГ выводятся в процессе минимизации невязки аналогичной вышеописанной.
Уравнения МПГ имеют вид Отметим, что в одномерных уравнениях Система уравнений (2.81) является основной системой конечно-разностных уравнений МПГ для пространственной кинетики в одномерной геометрии.
Система уравнений (2.82) решается только в том случае, если необходимо восстановить ход функций распределения нейтронов и предшественников запаздывающих нейтронов внутри ячеек.
Далее будем рассматривать двумерный случай для трех пробных матриц (2.80). Итерационную схему будем строить, используя полностью неявную схему Эйлера по времени. Для дискретизации временных производных используется аппроксимацию первого порядка:
В соответствии с (2.85), (2.86) и фактом использования полностью неявной схемы получим Сведем систему уравнений (2.87) к системе линейных алгебраических уравнений где M (t s ) - матрица, S(t s 1 ) - вектор.
Преобразуем первое уравнение системы (2.87) в следующую форму Получим выражение для C(j0,k) (t s ) из второго уравнения системы (2.87) Подставим (2.90) в (2.89) и получим Тогда матрица M (t s ) и вектор S(ts 1 ) из (2.88) имеют вид тогда Блок-схема решения нестационарных уравнений МПГ представлена на рисунке 2.1. Временная точка tS – последняя временная точка расчета.
Далее рассмотрим решение системы уравнений (2.88). Отметим, что для ее решения может быть применен любой способ решения систем линейных алгебраических уравнений.
Для решения системы уравнений (2.88) может быть использована традиционная схема – «группа за группой». В соответствии с этим подходом система (2.88) может быть переписана в следующем виде где (gn ) (t s ) – пространственный вектор, соответствующий энергетической группе g, M gg (t s ) – пространственная матрица с фиксированными энергетическими группами g и g, n – номер внешней итерации.
Для того чтобы решить уравнение (2.94) использовался метод верхней релаксации. Матрица M gg (ts ) раскладывается на диагональную матрицу D gg (ts ), Рисунок 2.1 – Блок-схема решения нестационарных уравнений МПГ нижнюю треугольную и верхнюю треугольную матрицы Lgg (t s ) и U gg (t s ) соответственно:
В итоге (gn) (ts ) итерационная схема может быть записана в виде где m – номер внутренней итерации.
В данной главе представлен детальный вывод двумерных конечноразностных нестационарных уравнения метода поверхностных гармоник для трех и четырех пробных матриц в квадратной решетке блоков. Отметим, что при выводе данных уравнений использован подход, заключающийся в том, что для получения пробных функций решается стационарная форма уравнения переноса нейтронов с возмущенными коэффициентами, соответствующими определенной временной точке.
Представленные уравнения метода поверхностных гармоник получены традиционным для МПГ методом минимизации невязки. В рассмотренном подходе невязка состоит из поверхностной и объемной невязок. Поверхностная невязка обусловлена тем фактом, что искомая функция точно не удовлетворяет уравнению переноса нейтронов на границах ячеек из-за ограниченного числа пробных матриц. Значения данной величины получено создателями метода и в данной главе приведено только для полноты описания процесса получения уравнений метода поверхностных гармоник. Объемная невязка обусловлена тем фактом, что пробные матрицы рассчитываются посредством решения стационарного уравнения переноса нейтронов. Также отметим, что данная невязка возникает только в нестационарном случае.
Отметим, что получению двумерных нестационарных уравнений МПГ предшествовал этап получения одномерных уравнений для одномерной плоской геометрии. Данная геометрия интересна тем, что в плоской одномерной ячейке возможны только две пробные матрицы – симметричная и антисимметричная.
Таким образом, набор пробных групповых функций образует полную систему решений. Как результат, в этой геометрии проще исследовать основные особенности применения метода поверхностных гармоник для нестационарного случая, не отвлекаясь на возможную недостаточность пробных решений. В данной главе приведены только результаты, т.к. процедура получения уравнений аналогична двумерному случаю.
Также в данной главе проведено построение итерационной схемы для решения нестационарных конечно-разностных уравнений МПГ для трех пробных матриц. Данное построение проведено с использованием полностью неявной схемы Эйлера по времени. Для дискретизации временных производных использовалась аппроксимацию первого порядка. Для решения полученной системы линейных алгебраических уравнений предложено использование традиционной схемы «группа за группой» с применением метода верхней релаксации.
Таким образом, в данной главе представлена математическая модель для решения нестационарного уравнения переноса нейтронов с запаздывающими нейтронами на основе метода поверхностных гармоник.
Глава 3 Программный комплекс SUHAM-TD и его верификация Данная глава посвящена описанию и верификации программного комплекса SUHAM-TD на ряде классических пространственно-временных бенчмарках.
Как было отмечено в главе 1, требования, предъявляемые современным нейтронно-физическим кинетическим кодам для моделирования переходных процессов в основном состоят в следующем:
1) Решаемая задача – уравнения пространственной кинетики;
2) Возможность проведение расчетов без диффузионного приближения;
3) Возможность проведения расчетов без пространственной гомогенизации;
4) Моделируемая среда – полномасштабная активная зона;
5) Высокая производительность.
Именно эти требования заложены в разрабатываемый программный комплекс SUHAM-TD. Опишем его основные характеристики, последовательно описывая пути удовлетворения этим требованиям.
пространственной кинетики посредством МПГ. Все модули SUHAM-TD разработаны на языке FORTRAN (стандарты 77 и 90). На текущий момент в программном комплексе реализованы нестационарные конечно-разностные уравнения МПГ для трех пробных матриц в квадратной решетке блоков. В диффузионные модули, реализующие метод поверхностных гармоник:
1) SUHAMD-PL-1D-SHM для одномерной плоской геометрии в стационарном 2) SUHAMD-TD-PL-1D-SHM нестационарном случае (использует результаты работы SUHAMD-PL-1DSHM);
3) SUHAMD-SQ-2D-SHM для квадратной решетки блоков в двумерной геометрии в декартовых координатах в стационарном случае;
4) SUHAMD-TD-SQ-2D-SHM для квадратной решетки блоков в двумерной геометрии в декартовых координатах в нестационарном случае (использует результаты работы SUHAMD-SQ-2D-SHM).
государственную регистрацию. Копии свидетельств о регистрации представлены в приложении А.
Диффузионные модули, реализующие классический метод конечных разностей имеют аналогичное наименование, за исключением трех последних букв (они заменяются на FDM), например SUHAMD-TD-SQ-2D-FDM. Область применения этих модулей аналогична.
приближения методом поверхностных гармоник:
1) SUHAM-SQ-2D-SHM для квадратной решетки блоков в двумерной геометрии в декартовых координатах в стационарном случае;
2) SUHAM-TD-SQ-2D-SHM для квадратной решетки блоков в двумерной геометрии в декартовых координатах в нестационарном случае (использует результаты работы SUHAM-SQ-2D-SHM).
Как видно, SUHAM-TD позволяет проводить расчеты, как с использованием диффузионного приближения, так и без него. Здесь важной чертой является то, что уравнения МПГ инвариантны по отношению к исходным уравнениям, из которых они получены, а именно, к уравнению переноса нейтронов и диффузионному уравнению и факт использования диффузионного приближения определяется только использованием его в расчете пробных матриц. Для получения диффузионных пробных матриц SUHAM-TD использует собственные процедуры, а для расчета пробных матриц без диффузионного приближения – программы РАЦИЯ [25] и DICPN [26]. Программа РАЦИЯ используется для расчета нулевой пробной матрицы методом поверхностных псевдоисточников в G3-приближении [109]. Программа DICPN применяется для получения первой и второй пробных матриц методом сферических гармоник в P2-приближении [27].
реализованная в SUHAM-TD, которая позволяет проследить место отмеченных программ в расчетном процессе (блок «расчет пробных матриц»). Также отметим, что факт использования процедуры пространственной гомогенизации определяется только программами расчета пробных матриц. Программы РАЦИЯ и DICPN не используют данную процедуру.
Рисунок 3.1 – Общая блок-схема нестационарного расчета по SUHAM-TD энергетических групп, групп запаздывающих нейтронов и числа расчетных зон.
периодичность системы присуще большинству современных ядерных реакторов.
Достигнуть большей производительности, чем существующие коды, предлагается за счет использования МПГ, который демонстрирует это в стационарном случае. Далее в данной главе в процессе верификации SUHAM-TD традиционным конечно-разностным методом и демонстрируется потенциал МПГ в вопросе сокращения вычислительных затрат.
3.2 Верификация программного комплекса SUHAM-TD Верификация программного комплекса SUHAM-TD проводилась на ряде классических бенчмарков, характеризующихся отсутствием неопределенностей в начальных данных и отсутствием обратных связей. Этот факт дает возможность выделить методическую погрешность МПГ, реализованного в SUHAM-TD. Это особенно важно, т.к. работа посвящена не разработке кода, реализующего хорошо известный и верифицированный метод, а развитию МПГ для нестационарных расчетов, исследованию его применения, а также разработке нового кода на его основе и его верификации.
Важным фактом при проведении верификации SUHAM-TD является соблюдение «чистоты эксперимента» в плане сравнения вычислительных затрат.
Для этого в рамках SUHAM-TD был специально разработан модуль, проводящий расчеты классическим конечно-разностным методом (МКР) по пространству и времени. Для сравнения вычислительных затрат МПГ и МКР вычисления по SUHAM-TD проводились с одинаковыми критериями точности и способами организации итераций (например, во всех программах матрицы обращаются с помощью метода верхней релаксации, везде используется полностью неявная схема). Также в программах не применяются какие-либо методы ускорения.
Отметим, что по условию всех рассматриваемых здесь тестов, для того чтобы рассчитать нестационарный процесс начальное состояние реактора делается критическим посредством деления сечения f во всех энергетических группах на keff.
Функционалы для сравнения. Расчеты, проведенные программным комплексом SUHAM-TD с помощью традиционного конечно-разностного метода, проведенные с помощью МПГ – SUHAM-SHM (Surface Harmonics Method). При верификации рассматривались следующие функционалы:
• начальное keff;
• полная (P) и зонные (Pi) мощности в выбранных точках по времени;
• распределение групповых потоков в выбранных точках по времени;
• вычислительные затраты.
Таким образом, проводилось сравнение следующих функционалов:
1. Относительное отклонение эффективного коэффициента размножения в стационарном расчете где k keff – эффективный коэффициент размножения, полученный с помощью МПГ, k keff – эффективный коэффициент размножения, полученный другим методом;
2. Относительное отклонение мощности в различных временных точках где P SHM (t ) – относительная мощность в момент времени t, полученная с помощью МПГ, P(t ) – относительная мощность в момент времени t, полученная другим методом;
3. Максимальное по модулю относительное отклонение по пространству потоков группы g в различных временных точках где i, g – нейтронный поток в зоне i группе g полученный с помощью МПГ, i, g – нейтроны поток в зоне i группе g полученный другим методом;
4. Среднеквадратическое отклонение RMS (root mean square) групповых потоков где N – число зон сравнения;
5. Относительные вычислительные затраты где t FDM - время расчета по программе SUHAM-FDM, t SHM -время расчета по программе SUHAM-SHM.
приведенных функционалов, были выбраны: h – пространственный шаг сетки, он же – шаг сетки при расчете пробных матриц МПГ в ячейках, H – размер ячейки в расчетах МПГ, для которых проводится расчет пробных матриц, шаг временной сетки t. Надо отметить, что все сравнения проводились при условии, что пространственный шаг сетки h в расчете SUHAM-FDM равен пространственному шагу сетки h, с которым в расчете SUHAM-SHM рассчитываются пробные матрицы.
используются следующие критерии сходимости итераций:
1. Для эффективного коэффициента размножения 2. Для потока нейтронов (расчет по МКР) во всех расчетных ячейках i во всех группах g 3. Для уровня нейтронов (расчет по МПГ) во всех расчетных ячейках i (ячейка МПГ) во всех группах g где n – номер итерации, i, g – поток нейтронов в g-ой группе i-ой зоне, i, g – уровень МПГ в g-ой группе i-ой ячейке определенный по формуле (2.50).
Как уже отмечалось в главе 2, в одномерном случае возможны только две координатные функции и в этом смысле система пробных решений полна. На основании этого факта расчет по МПГ и расчет классическим конечноразностным методом должны давать очень близкие результаты и возможно получение решения задачи практически с любой, заранее заданной точностью.
Чтобы увидеть разницу в этих расчетах были выбраны заведомо избыточные, чем в обычных случаях, критерии точности, а в таблицах числа приведены с большим, чем принято в практических расчетах, числом значащих цифр.
Таким образом, при расчете исходного стационарного состояния был выбран критерий сходимости = 10-9 для внутренних итераций по потоку (расчет SUHAM-FDM) и уровню (расчет SUHAM-SHM). Этот же критерий выбран для внешних итераций по потоку, уровню и собственному значению и для сверхвнешних итераций МПГ. В нестационарном случае используется только для внутренних итераций, а для остальных итераций выбран критерий = 10-7.
При расчете двумерных тестов, как в стационарном, так и в нестационарном случаях критерий = 10-9 использовался только для внутренних итераций, а для остальных итераций был выбран критерий = 10-7.
Описание задачи. Как было отмечено в главе 2, на первом этапе данной работы проводился вывод одномерных нестационарных уравнений МПГ. Для проведения верификации данных уравнений был выбран одномерный тест BSS- (Benchmark Source Situation) [20]. В этом тесте описана трехзонная размножающая среда (рисунок 3.2), при этом области 1 и 3 являются физически эквивалентными.
Граничные условия представлены в виде нулевых потоков на обеих границах.
Для гомогенных областей этого теста в двухгрупповом приближении задан полный набор макросечений взаимодействия, а также коэффициенты диффузии (таблица 3.1). Среднегрупповые скорости в среде для всех зон равны v1 = 1·10 см/с и v 2 = 3·10 см/с для первой и второй энергетических групп соответственно. Спектр деления для мгновенных и запаздывающих нейтронов также одинаков для всей среды и равен 1 = 1 и 2 = 0. Запаздывающие нейтроны представлены в шести группах (таблица 3.2). Отметим, что суммарная доля запаздывающих нейтронов равна = 0,0075.
Таблица 3.1 – Начальные двухгрупповые константы задачи BSS- Таблица 3.2 – Параметры запаздывающих нейтронов задачи BSS- Рассматривались следующие варианты бенчмарка с разными законами ввода реактивности:
• A1 - сечение поглощения нейтронов в первой зоне во второй энергетической группе линейно увеличивается на 3% в течение первой секунды;
• A2 - сечение поглощения нейтронов в первой зоне во второй энергетической группе линейно уменьшается на 1% в течение первой секунды;
• A3 - сечение поглощения нейтронов в первой зоне во второй энергетической группе линейно уменьшается на 5% в течение первой 0,01 секунды;
• A4 - сечение поглощения нейтронов в первой зоне во второй энергетической группе линейно уменьшается на 5% в течение первой 0,01 секунды, при этом скорости нейтронов на два порядка больше.
В дополнении к этим тестам был проведен расчет теста A0 (его нет в [20]) представляющего собой задачу аналогичную A1, но с возмущением в 1 и 3 зонах (симметричный тест). Динамическое поведение исследуется на интервале 0 t 4 с для всех рассматриваемых случаев.
Стационарное моделирование.
конфигурация приводилась к критическому состоянию путем деления сечений f на keff (квазистационарная задача). В таблице 3.3 приведены сравнительные результаты расчета keff начального состояния по двум программам по SUHAM-TD с помощью МПГ и МКР.
Таблица 3.3 – Эффективный коэффициент размножения в начальном состоянии
SUHAM-FDM SUHAM-SHM
2,50 0,901801756 10 0,901801755 1,1e-7 0, Здесь и далее полужирным шрифтом выделено максимальное отклонение.Из таблицы 3.3 видно:
• сходимость расчета keff в расчетах SUHAM-FDM и SUHAM-SHM с уменьшением пространственного шага сетки h;
• значения keff, в расчете SUHAM-SHM с разными размерами ячеек практически совпадают со значениями, рассчитанными в SUHAM-FDM (максимальное keff = 1,9·10-5 %);
• наличие практической независимости рассчитанного по МПГ значения keff от размера ячейки H;
• преимущество МПГ по вычислительным затратам проявляется для достаточно значительных пространственных размерностей (для шага сетки равного 0,25 см и 0,5 см), когда уменьшение размерности крупносеточного расчета не компенсируется увеличением размерности задач расчета пробных матриц.
Сравнение средних по локальным зонам групповых потоков показало, что для h = 0,25 см max не превышает 4,4·10-4 %, а RMS – не превышает 2,3·10-4 %. Также было проведено сравнение зонных мощностей (Pi, i – номер гомогенной зоны соответствующей рисунку 3.2) в начальной конфигурации.
Данное сравнение проводилось при h = 0,25 см. Относительные отклонения для первой, второй и третьей зон составили P1 = 1,6·10-4 %, P2 = 1,1·10-4 % и P3 = 1,4·10-4 % соответственно.
Кроме этого были рассчитаны значения асимптотической реактивности as, которая определялась по формуле:
где – реактивность и эффективный коэффициент размножения соответственно в начальный момент времени,, keff – реактивность и эффективный коэффициент размножения соответственно в момент окончания ввода реактивности. В таблице 3.4 приведены результаты расчета k eff, k eff и as (в относительных единицах и долях запаздывающих нейтронов) с помощью метода поверхностных гармоник на сетке с параметрами h = 0,25 см, H = 2 см для различных вариантов теста BSS-6.
Таблица 3.4 – Величина вводимой асимптотической реактивности сравнительного анализа. За эталон будем принимать расчет классическим конечно-разностным методом. Однако возникает вопрос о размере шага разбиения по пространству и по времени для эталонного расчета по МКР.
Рассмотрим обоснование диффузионного эталона на примере варианта А0.
Как показали расчеты мощностей, величины t = 10-4 с достаточно для временной сетки, т.к. максимальное отклонение зонной мощности для расчета с данным шагом по времени от расчета с шагом по времени равным t = 10-3 с составило 8·10-3 %. Данные расчеты представлены в таблице Б.1.
В таблице Б.2 представлены расчеты мощности (t = 10-2 с) при различных расчетных шагах по пространству. Из данной таблицы видно, что шаг пространственной сетки h = 0,25 см является достаточным, т.к. максимальное относительное отклонение при этом шаге от расчета с шагом h = 0,5 см составляет 6,5·10-3 %. Следует отметить, что дальнейшее уменьшение шагов t и h приводит к тому, что отклонение становится меньше при достаточно значительном росте временных затрат. Зависимости P от t и от h представлены на рисунке 3.3.
Результаты, аналогичные представленным, были получены также для вариантов A1 и A2. Таким образом, эталонным будем считать в дальнейшем SUHAM-FDM расчет с шагами t = 10-4 с, h = 0,25 см.
Для тестов BSS-6-А3 и BSS-6-А4 проведено аналогичное исследование для выбора эталона. Следует отметить, что тесты BSS-6-А3 и BSS-6-A4 отличаются от выше рассмотренных, тем, что процесс постепенно переходит из процесса на запаздывающих нейтронах в процесс на мгновенных нейтронах (см. таблицу 3.4), а тест BSS-6-A4 имеет к тому же скорости нейтронов на два порядка больше ( v1 = 1·109 см/с и v 2 = 3·107 см/с), т.е. задача более жесткая.
характеристики эталона: t t = 10-7 с, h = 0,25 см для задачи BSS-6-A4.
Рисунок 3.3 – Максимальное относительное отклонение мощности в расчете SUHAM-FDM: (а) – в зависимости от временного шага (при h = 0,25 см), (б) – в зависимости от пространственного шага (при t = 10-2 с).
Нестационарное моделирование. Перейдем к сравнению МПГ с прямым численным моделированием в нестационарных случаях. В таблице 3.5 приведены значения максимальных отклонений зонных мощностей (зоны представлены на рисунке 3.2), полученных в расчете SUHAM-SHM, от значений, полученных в расчете SUHAM-FDM при одинаковых значениях расчетного шага по времени t с пространственными шагами h = 0,25 см, H = 2 см.
На примере варианта A0 можно сделать вывод, что сравнение расчета SUHAM-SHM при t = 10-4 c, с эталонным дает значение максимального относительного отклонения 4,3·10-2 %. Данное значение практически не зависит 4,2–4,3·10-2 %. Следует отметить тот факт, что значение максимального относительного отклонения всегда соответствует последней временной точке и зоне, где происходит возмущение (см. таблица Б.3). Из данных таблиц 3.5 и Б. можно сделать вывод, что расчет SUHAM-SHM имеет незначительные отклонения от прямого численного моделирования (при одинаковых шагах по времени), однако, они возрастают по мере движения к последним расчетным точкам исследуемого временного интервала.
Таблица 3.5 – Максимальные значения среди относительных отклонений зонных мощностей Pi, полученных в расчетах SUHAM-SHM от значений в расчетах SUHAM-FDM (t для SUHAM-SHM равен t для SUHAM-FDM), % Проведем анализ отклонений мощностей в зависимости от размера ячейки.
Данный анализ необязательно проводить при мелком разбиении временной сетки, поэтому расчет проведен с шагом t = 10-2 с для вариантов А0–А2, и с шагами t = 10-5 с и t = 10-6 с для вариантов А3 и А4 соответственно. В таблице 3. приведены значения максимальных отклонений зонных мощностей в зависимости от шага крупносеточного разбиения Н при h = 0,25 см.
Из данной таблицы видно, что отклонение существенно возрастает с ростом 1,3·10-2 % до 1%). Также необходимо отметить, что погрешность от увеличения размера ячейки наиболее заметна в тестах А3 и А4, что указывает на необходимость использования более мелких ячеек в случае быстрого ввода реактивности.
Таблица 3.6 – Максимальные отклонения зонных мощностей SUHAM-SHM от SUHAM-FDM Pi в зависимости от шага крупносеточного разбиения Н, % Таблица Б.4 иллюстрирует данное исследование более подробно для варианта A0. В ней приведены значения полной и зонных мощностей SUHAMFDM, а также отклонения значений, полученных в расчете SUHAM-SHM, от максимального относительного отклонения практически не зависит от шага h (например, при размере ячейки равном H = 5 см, отклонение составило при различных значениях геометрического расчетного шага h величину примерно равную 2,4·10-1 %).
Проведем анализ отклонений мощностей в зависимости от расчетного шага по времени. В таблице 3.7 представлены максимальные отклонения мощностей, полученных МПГ при разных значениях t, от эталонного расчета при h = 0, см, H = 2 см. Из таблицы видно, что максимальное отклонение практически всегда соответствует максимальному шагу по времени (например, для варианта A0 оно максимально при t = 10-1 с и составляет 4,8·10-1 %).
Таблица Б.5 представляет данное исследование более детально для варианта A0. В ней приведены значения полной и зонных мощностей (SUHAM-SHM, h = 0,25 см) при различных шагах по времени, а также отклонения данных значений от эталонных значений. Из этой таблицы видно, что максимальное отклонение для определенного шага по времени в данном случае не обязательно относится к последней расчетной точки расчетного интервала и зоне, где происходит возмущение (при больших значениях t).
«