Развитие метода неполной факторизации и его применение в практических задачах нейтронной кинетики
УДК 519.6:621.039.5
На правах рукописи
Троянова Надежда Михайловна
РАЗВИТИЕ МЕТОДА НЕПОЛНОЙ ФАКТОРИЗАЦИИ
И ЕГО ПРИМЕНЕНИЕ В ПРАКТИЧЕСКИХ ЗАДАЧАХ
НЕЙТРОННОЙ КИНЕТИКИ
Специальность
05.13.18 – Математическое моделирование,
численные методы и комплексы программ
АВТОРЕФЕРАТ
диссертации на соискание ученой степени
кандидата физико-математических наук
Обнинск 2013
Работа выполнена в Государственном научном центре Российской Федерации – Физико-энергетическом институте имени А. И. Лейпунского (ГНЦ РФ-ФЭИ).
Научный руководитель:
доктор физико-математических наук,
советник директора ИРМиТ ГНЦ РФ-ФЭИ
Гинкин Владимир Павлович
Официальные оппоненты:
доктор физико-математических наук, профессор,
главный научный сотрудник ГНЦ РФ-ФЭИ
Коробейников Валерий Васильевич
доктор физико-математических наук,
главный научный сотрудник ИПМ им. М.В.Келдыша РАН
Кулешов Андрей Александрович
Ведущая организация:
ИБРАЭ РАН, г.Москва
Защита диссертации состоится «____» __________ 2013 года в _____ часов на заседании диссертационного совета Д 201.003.01 при ГНЦ РФ-ФЭИ по адресу: 249033, г. Обнинск, Калужской обл., пл. Бондаренко, д. 1.
С диссертацией можно ознакомиться в библиотеке ГНЦ РФ-ФЭИ.
Автореферат разослан «____» ____________ 2013 г.
Ученый секретарь ГНЦ РФ-ФЭИ
доктор технических наук Т.Н.Верещагина
ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ
Решение задач нейтронной кинетики реакторов с учетом локальных изменений нейтронно-физических характеристик является одной из наиболее сложных проблем моделирования работы реакторных установок различного типа и назначения. При решении таких задач необходимо иметь эффективные методы решения систем линейных алгебраических уравнений большой размерности с разреженными матрицами коэффициентов.
Один из таких методов решения систем линейных разностных уравнений - метод неполной факторизации (МНФ) - был впервые предложен Н. И. Булеевым в 1959 году и положил начало серии итерационных схем, эффективных в определенных классах задач. В последующие годы МНФ интенсивно развивался и стал одним из основных инструментов решения пространственно-временных задач в различных приложениях, в том числе в задачах расчета атомных реакторов. В развитие этих методов автором диссертации были разработаны новые схемы и испытаны различные варианты МНФ в применении к аппроксимации групповых уравнений диффузии нейтронов системами конечно-разностных уравнений: сверхнеявная схема с периферийной компенсацией итерируемых членов, схема с мерцающим параметром, модификация схемы Шнайдера-Зедана, явная схема в качестве предобусловливателя для метода сопряженных градиентов; неявные и комбинированные схемы. Автором диссертации разработаны программы, реализующие предложенные схемы, и выполнены сравнительные численные исследования их сходимости.
Из числа приложений с использованием МНФ в диссертацию включены два: комплекс программ WIMS-ВОЛНА для расчета нейтронных полей реакторов на тепловых нейтронах типа ВВЭР в (hex, z) -геометрии и комплекс программ GVA (GRIF-SM, VOLNA, ARAMAKO) для расчета быстротекущих процессов в реакторах на быстрых нейтронах типа БН. Алгоритмы и программы нейтронно-физического расчета в этих комплексах разработаны автором диссертации. С помощью первого комплекса программ выполнялись реперные и прогнозные расчеты кампаний различных блоков реакторов ВВЭР-1000 применительно к проблеме деформаций тепловыделяющих сборок (ТВС) в активных зонах в процессе эксплуатации топлива и создания ТВС нового поколения ТВС-2 и ТВС-2М. С помощью второго комплекса были просчитаны быстрые переходные процессы в реакторе типа БН-800, вызванные несанкционированным движением стержней регулирования.
Актуальность работы. Работы по созданию новых эффективных численных алгоритмов для решения различных задач математической физики всегда были и остаются актуальными задачами вычислительной математики. Схема практического решения задачи пространственно-временного расчета реакторов с учетом локальных изменений нейтронно-физических характеристик активной зоны заключается в сведении исходной задачи к задачам, которые могут быть эффективно решены различными методами теории аппроксимации и линейной алгебры. В частности, задача, успешное и эффективное решение которой определяет в большой степени точность и скорость решения исходной задачи, - это решение системы разностных уравнений большого порядка с разреженной матрицей коэффициентов, часто обладающей плохой обусловленностью. Поэтому важным и актуальным является развитие методов решения плохо обусловленных систем линейных уравнений с разреженными матрицами.
Появляющиеся в практике реакторной эксплуатации новые сложные актуальные задачи требуют своего решения. К таким задачам относятся расчет нейтронно-физических характеристик реакторов типа ВВЭР с учетом деформации ТВС в активной зоне реактора в процессе его эксплуатации под воздействием неоднородных тепловых и нейтронных полей и расчет быстротекущих процессов в реакторах на быстрых нейтронах типа БН.
Среди работ команды специалистов, решавших техническую задачу термомеханической стабильности ТВС в составе активных зон ВВЭР-1000 и разработки новой геометрически стабильной при длительной эксплуатации конструкции ТВС нового поколения (ТВС-2 и ТВС-2М), расчеты локальных нейтронных полей и энерговыделения с учетом изменяющейся геометрии активной зоны и водо-урановых соотношений имели самую высокую актуальность.
Пространственно-временной расчет быстрых переходных процессов основан на решении нестационарного уравнения реактора квазистатическим методом, который позволяет выделить два временных масштаба: мелкий - для расчета быстроменяющегося амплитудного фактора (мощности реактора) и крупный - для расчета форм-функции (пространственных распределений групповых потоков нейтронов). Этот метод позволяет выполнить эффективный и экономичный пространственно-временной расчет быстрых переходных процессов в реакторе типа БН, особенно актуальный для анализа начальных стадий их развития.
Цели работы.
- Разработать новые эффективные варианты метода неполной факторизации для решения алгебраических систем уравнений, являющихся разностными аналогами уравнений эллиптического типа в (hex, z) –геометрии.
- Создать библиотеку подпрограмм для решения алгебраических систем уравнений на основе метода неполной факторизации.
- Разработать алгоритмы и программы нейтронно-физического расчета для расчета кампаний реакторов ВВЭР с использованием новых эффективных схем метода неполной факторизации; рассчитать пространственные нейтронно-физические характеристики в течение кампании реактора для различных блоков АЭС с различными загрузками и схемами перегрузок топлива с учетом формоизменения ТВС в процессе эксплуатации для обеспечения термомеханических расчетов активных зон ВВЭР и разработки топливных кассет нового типа.
- Разработать алгоритмы и программы нестационарного нейтронно-физического расчета переходных процессов в реакторах на примере БН-800 с использованием новых эффективных схем метода неполной факторизации; выполнить пространственно-временные расчетные исследования нейтронно-физических характеристик активной зоны реактора типа БН-800 при моделировании таких процессов, как несанкционированное извлечение стержней регулирования, в совместном теплогидравлическом и нейтронно-физическом расчёте.
Научная новизна результатов диссертации.
- Разработаны эффективные варианты МНФ с периферийной компенсацией итерируемых членов для решения уравнений диффузионного типа: сверхнеявная схема, схема с мерцающим параметром, модификация схемы Шнайдера-Зедана, явная схема в качестве предобусловливателя для метода сопряженных градиентов; неявные и комбинированные схемы.
- Предложен новый формализм описания МНФ для решения систем разностных уравнений, аппроксимирующих уравнения эллиптического типа, с произвольной нумерацией узлов разностной сетки.
- Создана библиотека подпрограмм IFML для решения систем линейных алгебраических уравнений на основе МНФ.
- Разработана трехмерная программа ВОЛНА расчета нейтронной кинетики реактора в многогрупповом диффузионном квазистационарном приближении с учетом формоизменения ТВС в составе комплекса программ WIMS-ВОЛНА и выполнены расчеты кампаний реакторов ВВЭР-1000 со смешанными загрузками и различными схемами перегрузок в обеспечение термомеханических расчетов активных зон ВВЭР.
- Разработана трехмерная программа VOLNA нестационарного расчета реактора в квазистатическом приближении с учетом теплового расширения активной зоны в составе комплекса программ GVA совместного нейтронно-физического и теплогидравлического расчета реакторов типа БН-800 и выполнены расчеты переходных процессов, обусловленных движением группы стержней.
Достоверность полученных результатов. Универсальность и высокая эффективность разработанных вариантов МНФ, особенно комбинированных схем, подтверждается численными сравнительными исследованиями и практическими приложениями. Достоверность результатов нейтронно-физических расчетов подтверждается хорошим согласием с расчетами тестовых задач по другим программам. Достоверность выводов диссертации подтверждается практикой проведения расчетов, результатами разработки и внедрения топливных кассет нового поколения на АЭС России, конструирование которых опиралось на результаты расчетов нейтронно-физических характеристик активных зон ВВЭР, результатами расчетов активных зон реакторов типа БН-800.
Научная ценность работы заключается в разработке новых эффективных схем МНФ и их программной реализации, что имеет существенное значение в области численных методов решения систем алгебраических уравнений с плохообусловленными матрицами. Эти методы нашли применение в решении практических задач нейтронной кинетики.
Практическая ценность разработанных эффективных методов решения плохо обусловленных систем линейных алгебраических уравнений большой размерности с разреженными матрицами коэффициентов и созданной на их основе библиотеки подпрограмм IFML состоит в применимости этих подпрограмм для решения широкого круга практических задач математической физики. Практическая ценность разработанных комплексов программ нейтронно-физического расчета заключается в большом количестве конкретных рекомендаций и выводов, сделанных по результатам расчетов по этим программам, и принятых проектно-конструкторских решений при разработке и внедрении нового топлива.
Автор выносит на защиту:
- Новые схемы, алгоритмы и программные реализации различных вариантов МНФ с периферийной компенсацией итерируемых членов и новые варианты комбинированных схем.
- Алгоритмы и программы решения систем трехмерных разностных уравнений диффузионного типа в (hex, z) – геометрии с использованием нового формализма описания МНФ для произвольной нумерации узлов разностной сетки.
- Библиотеку подпрограмм неполной факторизации IFML, предназначенных для решения разностных систем двумерных и трехмерных уравнений эллиптического типа.
- Алгоритм и программу ВОЛНА многогруппового нейтронно-физического расчета кампаний реакторов ВВЭР с учетом формоизменения ТВС в процессе эксплуатации в комплексе программ WIMS-ВОЛНА.
- Комплекс параметрических и поддерживающих расчетов кампаний с различными загрузками в обеспечение термомеханических расчетов активных зон ВВЭР-1000 и в обоснование внедрения на АЭС нового топлива.
- Алгоритм и программу VOLNA нестационарного многогруппового нейтронно-физического расчета быстрых переходных процессов в реакторе типа БН-800 в квазистатическом приближении с учетом теплового расширения активной зоны для комплекса программ GVA совместного нейтронно-физического и теплогидравлического расчета.
Личный вклад автора заключается в разработке новых явных и неявных комбинированных схем метода неполной факторизации с периферийной компенсацией итерируемых членов, обладающих более высокой эффективностью при решении систем линейных алгебраических уравнений с плохообусловленными матрицами по сравнению со схемами Булеева, Гинкина, Шнайдера-Зедана. Алгоритмы на основе этих схем реализованы в расчетных программах-решателях, которые в свою очередь были применены в алгоритмах и программах ВОЛНА и VOLNA, также написанных автором. Практические расчеты автор выполнял лично.
Апробация работы. Основные положения и результаты диссертационной работы докладывались на Всесоюзном совещании по динамике реакторов (г.Гатчина, 1990), на 13-ом Международном конгрессе по вычислительной и прикладной математике (г. Дублин, Ирландия, 1991), на VI Российской конференции по радиационной защите ядерных установок (г.Обнинск, 1994), на 9 международном совещании по безопасности ядерных реакторов (г.Москва, 1995), на семинарах «Нейтроника-97, -98», (г.Обнинск, 1997, 1998), на школе-семинаре МИФИ «Интегрированные математические модели и программные комплексы в ядерной энергетике» (г.Москва, 1998), на семинаре «Консультативная встреча специалистов IAEA по проблеме быстрых реакторов» (г.Обнинск, 1998), на международном семинаре “Деформация топливных сборок PWR и ВВЭР”, (г. Ржеж, Чехия, 1998), на международном семинаре «Анализ безопасности атомных станций с реакторами типа ВВЭР и РБМК» (г. Обнинск, 1998), на 4-ом международном конгрессе по прикладной математике (г. Эдинбург, Шотландия, 1999), на международной конференции «Математические идеи П.Л.Чебышева и их приложение к современным проблемам естествознания» (г. Обнинск, 2002), на Техническом комитете МАГАТЭ «Structural behaviour of fuel assemblies for water cooled reactors» (г.Кадараш, Франция, 2004).
Основное содержание диссертации изложено в шести публикациях в отечественных и международных журналах и сборниках трудов конференций, пяти препринтах, трудах отраслевых и международных семинаров, сборниках работ организаций и 58 научно-технических отчетах.
Структура и объем работы. Диссертация состоит из введения, 3 глав, заключения и выводов. Общий объем диссертации – 157 страниц, в том числе 51 рисунок и 31 таблица, список использованных источников содержит 47 наименований на 5 страницах.
СОДЕРЖАНИЕ РАБОТЫ
Во введении содержится обоснование актуальности темы диссертационной работы, сформулированы цели и основные положения, выносимые на защиту, приведена информация о структуре диссертации.
Первая глава содержит описание различных схем МНФ и предложения по их развитию. Описание дается на языке разностных схем, традиционно используемом в математическом моделировании при решении задач математической физики в различных приложениях, а также в предложенном автором формализме матриц связности в случае произвольной нумерации пространственных расчетных узлов.
Метод неполной факторизации для решения систем разностных уравнений эллиптического типа впервые был предложен Н.И. Булеевым.
Пусть требуется решить уравнение 
, где А – невырожденная, но труднообратимая матрица. Прибавим к обеим частям уравнения вектор Ву и выберем матрицу В такой, чтобы матрица (А+В) представлялась в виде произведения двух легкообратимых матриц М и N: 
. Тогда получим 
. Обозначая 
, получим систему двух уравнений, которую будем решать методом последовательных приближений:
,
.
Выбор различных матриц M, N и B определяет различные схемы метода неполной факторизации.
По типам факторизованных операторов схемы неполной факторизации условно делятся на явные и неявные схемы. К явным относятся схемы, в которых расчет в прямом и обратном направлениях ведется по обычным рекуррентным формулам, а к неявным – схемы, в которых по крайней мере в одном из направлений счета приходится использовать метод прогонок.
Практика расчетов показала, что для большинства задач, особенно с несимметричной исходной матрицей А, неявные схемы более устойчивы и обеспечивают более высокую скорость сходимости, чем явные схемы
Основное внимание в диссертации уделено неявным и комбинированным схемам, под руководством В.П. Гинкина и с его участием разработаны новые схемы МНФ и при численном решении тестовых задач исследована их эффективность.
Для симметричных матриц также могут быть построены эффективные схемы с использованием идей неполной факторизации. В работе описан эффективный вариант метода сопряженных градиентов с предобусловливателем по явной схеме неполной факторизации с диагональной и периферийной компенсацией итерируемых членов для решения симметричных задач в разных геометриях.
Пусть требуется решить уравнение Ay=f, где А – конечно-разностный оператор, аппроксимирующий уравнение эллиптического типа, и этот оператор действителен и симметричен, т.е. A = AT.
Метод сопряженных градиентов (CG) для решения этой системы уравнений имеет вид
, 
, 
, 
,
Выберем в качестве предобусловливающей матрицы К факторизованную матрицу из МНФ К=(A+B)=MN, где N=MT. При этом оператор B тоже симметричен, его конкретный вид определяет «качество» предобусловливающего оператора и существенно влияет на скорость сходимости итерационного процесса.
Построение предобусловливателя с выбором в качестве оператора M нижней треугольной матрицы для двумерного случая в (x-y)-геометрии на пятиточечном конечно-разностном шаблоне:
приводит к итерируемым выражениям By с неустранимыми членами, появляющимися от перемножения операторов M и N=MT и выходящими за рамки основного конечно-разностного шаблона:
Для ускорения сходимости МНФ применим периферийную компенсацию итерируемых членов по аналогии со схемой h-факторизации Гинкина для неявных схем, выбрав в качестве оператора В следующий оператор:
где 
- итерационный параметр. В отличие от схемы h-факторизации, приводящей к несимметричной матрице В, в данном случае матрица В остается симметричной, что позволяет использовать полученную схему МНФ в качестве предобусловливателя в методе сопряженных градиентов.
Аналогично были построены схемы МНФ с периферийной компенсацией итерируемых членов для трехмерного случая в (x-y-z)-геометрии на семиточечном шаблоне
и (hex-z)-геометрии на девятиточечном шаблоне
Здесь для каждого итерируемого члена в операторе В используется периферийная компенсация (PIF) по типу h-факторизации в виде тройки членов, взятых в соседних узлах разностной сетки, и введен параметр компенсации 
. При 
компенсация отсутствует; при 
- компенсация полная; при 
- компенсация частичная. Коэффициенты, задающие матрицу M, находятся из тождества By=(MMT-A)y и определяются по рекуррентным формулам. Использовалась также схема с диагональной компенсацией (DIF), предложенная Булеевым.
Численные исследования сходимости метода сопряженных градиентов с предложенным предобусловливателем при решении тестовых задач показали, что использование схем PIF и DIF в качестве предобусловливателя в методе сопряженных градиентов (CG) приводит к значительному ускорению сходимости этого метода, зависимость количества итераций от числа узлов n по одному направлению при оптимальном значении параметра 
приблизительно пропорциональна 
, оптимальное значение параметра 
лежит в области от 0.9 до 1.0 для различных типов задач. Наивысшей скоростью сходимости при 
из рассмотренных схем обладает схема CG-PIF.
Автором предложено новое описание схем неполной факторизации с произвольной нумерацией узлов разностной сетки, когда каждому узлу ставится в соответствие набор номеров узлов-соседей с использованием матрицы связности узлов сетки. Дана общая формулировка метода неполной факторизации с помощью предложенного формализма и построена новая сверхнеявная схема неполной факторизации с периферийной и диагональной компенсацией итерируемых членов для решения систем трехмерных уравнений диффузионного типа в (hex-z)-геометрии, в которой исходный девятиточечный оператор заменяется произведением двухточечного и восьмиточечного операторов, а для обращения восьмиточечного оператора используется разработанный автором модифицированный вариант явной схемы неполной факторизации Шнайдера-Зедана для решения двумерных систем разностных уравнений с семидиагональными матрицами.
Согласно сеточному шаблону факторизуемого оператора DA+B можно с использованием специальной матрицы связности MS(j,i) для каждого расчетного узла i, i=1,…N0, определить J0 номеров узлов-соседей, соответствующих остальным узлам шаблона, j=1,…J0. При этом трехмерной (hex-z)-геометрии соответствует девятиточечный шаблон с J0=9, j=9 – центральный узел, j=1,..6 – узлы гексагональной решетки в сечении z, j =7, 8 соответствуют предыдущему (j=7) и последующему (j=8) сечению по координате z.
Исходное разностное уравнение в (hex-z)-геометрии при использовании введенного формализма примет вид:
.
Для решения этого уравнения предложена новая сверхнеявная схема МНФ:
где оператор D выбирается диагональным, а операторы M и N – двухточечным и восьмиточечным, соответственно:
Сверхнеявность этой схемы заключается в том, что для решения последнего уравнения, определенного на восьмиточечном шаблоне, нельзя применить прямой экономичный метод с использованием только простых неявных схем типа одномерных прогонок. Поэтому для решения этого уравнения использована также итерационная схема неполной факторизации - новая модернизированная автором схема Шнайдера-Зедана, в которой с целью ускорения сходимости используется периферийная компенсация итерируемых членов по типу h-факторизации с итерационным параметром.
В оператор B входят дополнительные узлы с номерами MS(7,MS(j)), j=1,…6, появившиеся при факторизации трехмерных уравнений от перемножения M и N.
Вклад в оператор B дополнительных узлов можно компенсировать комбинацией членов, взятых в узлах основного шаблона. Пусть 
, оператор R определен на дополнительном шаблоне 
S компенсирующий оператор 
где 
и 
- параметр и матрица коэффициентов компенсации. Положим 
и определим компенсирующий оператор S по типу h-факторизации.
Тогда получим 
Двухточечный оператор M обращается простым пересчетом в сторону возрастания координаты z; обращение же оператора N ведется в сторону уменьшения координаты z, причем для каждого фиксированного z необходимо обратить двумерный семиточечный оператор вида:
размерностью M0, где M0 – количество узлов в hex-плоскости.
Для этого автором выбрана и модифицирована введением параметрической периферийной компенсации по типу h-факторизации схема Шнайдера-Зедана. Схема метода неполной факторизации для этого случая выглядит так:
где n - номер текущей итерации, 
с нечетными j– коэффициенты оператора 
, с четными j - коэффициенты оператора 
, причем 
и 
имеют единичные диагонали. Оператор 
будет содержать элементы в дополнительных узлах MS(1,MS(4)) и MS(5,MS(2)), вклад которых можно компенсировать комбинацией членов, взятых в узлах основного шаблона. Выберем 
, где R1 и R2 – операторы, определяющие итерируемые члены:
и 
, а S1 и S2 -компенсирующие операторы, определенные в узлах основного шаблона: 
и 
- параметр и матрица коэффициентов компенсации. Компенсирующие операторы выбираются по типу h-факторизации, и оператор B имеет вид:
Определив по рекуррентным формулам коэффициенты 
, при 
получаем схему Шнайдера-Зедана без компенсации, при 
получаем модифицированную автором схему Шнайдера-Зедана с периферийной компенсацией по типу h-факторизации, в которой первые члены разложения 
в ряд Тейлора в центральном узле шаблона равны нулю.
Уравнение для вспомогательной функции u решается простым пересчетом от узла к узлу в сторону возрастания координаты z, а уравнение для функции 
– в сторону убывания координаты z.
Кроме сверхнеявной схемы МНФ с постоянным параметром компенсации 
, автором сформулирован и исследован комбинированный МНФ, в котором значение параметра компенсации 
зависит от m - номера итерации.
Алгоритмы выбора параметра компенсации 
могут быть различными, например, простое чередование, циклический перебор и другие. Достаточно эффективной и экономной с точки зрения необходимости пересчета и (или) хранения коэффициентов операторов M и N для каждого 
оказалась комбинированная схема, получившая название схемы МНФ с мерцающим параметром - МНФ-МП, в которой для каждой нечетной итерации используется параметр компенсации, равный 0, а для каждой четной итерации – параметр, равный 1.
Эффект от чередования значений параметра компенсации (0 и 1) на каждой итерации оказался весьма высоким (зависимость числа итераций от числа узлов в схеме МНФ-МП при решении тестовых задач оказалась близкой к зависимости при оптимальном экспериментально подобранном параметре 
) и объясняется тем, что при параметре 
эффективно гасятся высокочастотные, а при 
– гладкие компоненты ошибки. Аналогичный эффект усиления сходимости наблюдался в двумерной комбинированной схеме HFPP В.П. Гинкина. Предложенный автором алгоритм чередования значений параметра не требует подбора оптимального значения параметра компенсации и особенно эффективен при решении плохообусловленных задач.
С целью систематизации и сохранения знаний, наличия большого количества различных схем и вариантов МНФ, реализованных разными разработчиками и предназначенных для решения различных классов задач, автором диссертации разработана структура и создана библиотека подпрограмм IFML для решения двумерных и трехмерных систем разностных уравнений методом неполной факторизации, в основном разработанных в математическом отделе ГНЦ РФ-ФЭИ. Библиотека IFML содержит 22 подпрограммы МНФ, 6 из которых разработаны автором, является открытой структурой и готова к дальнейшему наполнению по мере появления новых методов, разработки вариантов и схем неполной факторизации и их программной реализации.
Вторая глава посвящена алгоритмам и методам решения квазистационарных задач нейтронной физики. Основным содержанием этой главы является описание созданного трехмерного комплекса программ WIMS-ВОЛНА для расчета нейтронно-физических характеристик, локальных нейтронных полей и полей энерговыделений реакторов типа ВВЭР в течение кампании реактора с учетом перегрузок топлива. В этом комплексе расчет локальных нейтронных полей и энерговыделения производится по программе ВОЛНА, а константы рассчитываются в каждой точке реактора в каждый момент времени кампании по программе WIMSD4 с учетом изменения материального состава элементарных объемов реактора и средних температур топлива, оболочки твэла и теплоносителя. Результаты расчетов в виде данных по плотности потока нейтронов и энерговыделений в заданных пространственных точках в согласованных форматах передаются в программу расчета термомеханики активной зоны в заданные моменты кампании, а из программы расчета термомеханики в комплекс программ WIMS-ВОЛНА возвращаются данные о формоизменении ТВС. Предложен алгоритм учета влияния формоизменения ТВС на расчет нейтронно-физических характеристик реактора.
Автором разработаны алгоритмы и схемы конечно-разностной аппроксимации уравнения диффузии нейтронов в (hex, z) –геометрии и выполнены исследования по уточнению пространственной аппроксимации. В программе нейтронно-физического расчета реактора ВОЛНА в качестве решателей использовались разработанные автором эффективные схемы метода неполной факторизации.
Эти работы, а также работа по оценке влияния изменений величины водяных зазоров между ТВС во время эксплуатации реактора вследствие формоизменения ТВС на его нейтронно-физические характеристики, выполнялись в составе комплекса работ по созданию ТВС нового поколения (ТВС-2, ТВС-2М) и выработке методологии термомеханического анализа активных зон ВВЭР-1000, для чего выполнялись прогнозные расчеты термомеханических деформаций ТВС в активных зонах реакторов на различных блоках АЭС. Для основной цели комплексной работы – расчета термомеханического поведения активной зоны – нейтронно-физические расчеты являются обеспечивающими.
При решении уравнения реактора
для расчета кампании в квазистационарном приближении в программе ВОЛНА зададим временную сетку 
с переменным шагом 
, (i=1,...,m ), считая на шаге 
значения плотности потоков нейтронов не зависящими от времени. Трехмерное многогрупповое квазикритическое уравнение дискретизируется в (hex, z) -геометрии и решается методом итераций источника нейтронов деления с использованием трехслойной схемы с Чебышевским ускорением сходимости внешних итераций и быстросходящегося варианта МНФ на внутренних итерациях. В каждый момент времени в каждом элементарном расчетном объеме реактора производится прямой расчет по программе WIMSD4 изменения материального состава зоны вследствие выгорания топлива и других выгорающих материалов, перемещения органов регулирования, изменения температуры теплоносителя и соответствующих новым материальным композициям сечений взаимодействия. На каждом временном шаге кампании с помощью изменения концентрации бора в теплоносителе 
реактор выводится в критическое состояние при решении методом Ньютона нелинейного уравнения 
.
Для оценки ошибки, связанной с погрешностью пространственной аппроксимации в (hex-z) – геометрии при расчете распределений нейтронно-физических характеристик в программе ВОЛНА, были выполнены тестовые расчеты с заданными константами и известным точным решением модельной задачи. Наилучшие результаты получены в расчете с использованием семи точек на кассету в плане, когда каждая кассета аппроксимируется семью равными малыми шестиугольниками, эквивалентными по площади сечению ТВС.
Изменение величины водяного зазора между ТВС как в сторону увеличения, так и уменьшения приводит к изменению нейтронно-физических характеристик зоны: групповых макроконстант, радиальных распределений по сечению ТВС полей нейтронов и энерговыделений. Исходными предпосылками для необходимости разработки и исследования способа учета формоизменения ТВС послужили известные события по искривлению ТВС в активных зонах реакторов во время эксплуатации, связанные с этим трудности при перегрузках и перестановках ТВС во время ППР и база замеров искривлений (прогибов) ТВС разных блоков разных станций.
Необходимость учета формоизменения ТВС в активной зоне приводит к внесению изменений в методики и алгоритмы расчетов нейтронно-физических характеристик зоны от разработки специальных алгоритмов гомогенизации констант в рамках фиксированных сеток до построения следующих изгибам ТВС адаптивных пространственных сеток. Таким образом, при изучении нейтронно-физического поведения активной зоны с учетом обратных связей по термомеханике появляется еще один параметр, влияющий на константы, - водяной зазор.
Автором разработан алгоритм расчета макроконстант, позволяющий учесть неравномерный водяной зазор по периметру ТВС вследствие формоизменения ТВС в рамках фиксированной расчетной сетки, построенной в (hex-z) – геометрии, и исходной методики подготовки констант. Константы для каждой из семи малых гексагональных расчетных ячеек, аппроксимирующих ТВС в каждом слое по высоте, рассчитываются с учетом неравномерных зазоров между гранями ТВС, причем для центральной малой ТВС зазор полагается нулевым, для остальных шести малых ТВС - полусумме зазоров между соответствующими гранями ТВС.
Описание вариантных расчетов приведено в диссертации и иллюстрируется значениями длины кампании, максимального коэффициента неравномерности в энерговыделении ТВС Kq и по объему активной зоны Kv в зонах с формоизмененными ТВС по сравнению с неискривленной зоной. В ряде случаев эффект от воздействия смещения ТВС в активной зоне реактора на нейтронно-физические характеристики значителен и требует учета при проведении расчетов кампании.
Указанный комплекс программ и алгоритмы были использованы автором при решении задачи установления причин термомеханической нестабильности поведения активных зон реакторов ВВЭР-1000, наблюдавшейся в период с ~ 1990 по ~ 2006 гг. В составе группы специалистов, привлеченных к решению этой задачи, автор обеспечивал входными данными по пространственным распределениям нейтронных потоков и температур расчеты трехмерной термомеханики активных зон. В течение указанного периода в процессе проведения расчетно-параметрических исследований, верификации расчетных методов и программ, разработки обоснования новых конструктивных решений по ТВС и схемам загрузок, обоснования безопасности при повышении мощности энергоблоков до 104%, автором выполнены десятки параметрических и поддерживающих расчетов активных зон, в том числе нейтронно-физические характеристики активных зон моделировались в виде прогнозных расчетов, а также в виде тестовых расчетов по реально прошедшим кампаниям на конкретных энергоблоках. Расчеты нейтронно-физических характеристик реакторов, выполненные автором, использованы при разработке ТВС нового поколения (ТВС-2, ТВС-2М) и подготовке к их лицензированию. Эта работа была успешно выполнена и с 2003 года начата опытно-промышленная эксплуатация топливных кассет нового поколения, в разработку которых автор внёс свой небольшой вклад, а вскоре и промышленная их эксплуатация. Эти ТВС используются в промышленной эксплуатации на всех блоках Балаковской и Ростовской АЭС, предназначены для применения во всех строящихся блоках проекта АЭС-2006 и ВВЭР-ТОИ.
Третья глава посвящена трехмерной программе VOLNA нестационарного расчета нейтронно-физических характеристик активных зон во время быстрых переходных процессов в реакторах типа БН. Программа VOLNA интегрирована в комплекс программ GVA (GRIF-SM, VOLNA, ARAMAKO), позволяющий с константной поддержкой модуля АРАМАКО рассчитывать совместно с теплогидравлическим расчетом по программе GRIF-SM нейтронно-физические характеристики переходных процессов в реакторах на быстрых нейтронах. Комплекс программ GVA разрабатывался коллективом авторов под общим руководством Гинкина В.П. и Кузнецова И.А.; программа нейтронно-физического расчета VOLNA в составе этого комплекса разработана автором диссертации.
Программа расчета теплогидравлики GRIF-SM решает в (r, z) геометрии нестационарную систему уравнений сохранения массы, импульса и энергии в рамках модели пористого тела, причем пористость, проницаемость и другие коэффициенты системы уравнений модели зависят от пространственных координат. Пространственные распределения энерговыделения в каждый момент времени, рассчитываемые по программе VOLNA, являются для программы GRIF-SM входными данными.
Программа расчета пространственной кинетики VOLNA решает пространственно-временное уравнение переноса нейтронов в диффузионном квазистатическом приближении. Плотность потока нейтронов представляется в виде произведения пространственной форм-функции и амплитудного фактора, который предполагается не зависящим от пространственных координат. Исходное нестационарное уравнение тождественным преобразованием приводится к системе двух уравнений – системы уравнений типа точечной кинетики для амплитудного фактора P(t), интерпретируемого как мощность реактора:
;
и пространственно-временных уравнений для групповых форм-функций 
:
где 
; Q – источник нейтронов деления,
- скорость распада предшественников запаздывающих нейтронов в момент времени t.
Важной чертой квазистатического приближения является принципиальная возможность введения двух временных сеток: мелкой - для частого расчета амплитудного фактора и крупной - для относительно редкого пересчета трехмерной форм-функции, что приводит к экономии времени решения исходного нестационарного уравнения реактора.
Входными данными для программы VOLNA являются пространственные распределения материальных композиций и групповых макроконстант в каждый момент времени, скорости движения органов регулирования и другие параметры. Выходными данными программы VOLNA являются пространственные распределения энерговыделений по объему реактора, групповых плотностей потоков нейтронов и другие нейтронно-физические характеристики. Групповые макроконстанты для расчета по программе VOLNA вычисляются по программе AРAMAКO. Входными данными для этой программы являются пространственные распределения температур и плотностей топлива, конструкционных материалов и теплоносителя в каждый момент времени, получаемые по программе GRIF-SM.
Комплекс программ GVA состоит из двух блоков – стационарного расчета начального критического состояния реактора и блока нестационарного расчета реактора.
В расчете начального состояния в комплексе GVA производится итерационное согласование температурных и нейтронных полей по программам GRIF-SM и VOLNA. После достижения заданной точности реактор выводится в строго критическое состояние путем деления числа вторичных нейтронов на kэф. Затем решается сопряженная задача для определения функции ценности нейтронов.
Изменение макроконстант хотя бы в одной пространственной ячейке реактора вызывает изменение критического состояния и означает введение реактивности по той или иной физической причине. Реактивность вычисляется путем интегрирования по объему реактора изменений в распределениях макроконстант с весом форм-функции и функции ценности нейтронов деления для критического реактора.
На каждом крупном временном шаге, который выбирается на основе вычислительного опыта и требований программы GRIF-SM, в расчетной зоне производится полномасштабный пространственный пересчет температур по программе GRIF-SM, групповых макроконстант по программе АРАМАКО, нейтронных полей по программе VOLNA с итерациями по нелинейности для согласования полей энерговыделений и температур.
Система конечно-разностных уравнений с девятидиагональной матрицей коэффициентов для определения групповой пространственной форм-функции в программе VOLNA решается по схеме метода неполной факторизации МНФ-МП.
Алгоритм численного интегрирования пространственно-независимой системы уравнений типа уравнений точечной кинетики для амплитудного фактора P(t) построен на основе многошаговых методов Гира типа прогноза-коррекции и обеспечивает автоматический выбор шага мелкой временной сетки на каждом крупном временном шаге.
Важно заметить, что в комплексе GVA при расчетах используются две пространственные сетки: (r-z) - разбиение по кольцам в плоскости и высоте z расчетной зоны для программы GRIF-SM и (hex-z)- сетка для программы VOLNA, причем размеры слоев по высоте z в сетках могут не совпадать. Это обстоятельство требует использования алгоритмов сопряжения пространственных сеток и интерполяционных процедур при совместном расчете и обмене данными между программами GRIF-SM и VOLNA.
В программе GRIF-SM может быть учтено деформирование расчетной области реактора (РЗР) вследствие теплового расширения путем изменения размера расчетных кольцевых зон. В программе VOLNA тепловое деформирование РЗР моделируется изменением во времени эквивалентного размера под ключ расчетной гексагональной ячейки для каждого расчетного слоя по высоте РЗР, которое задается информацией, полученной из программы GRIF-SM.
В качестве примера в работе приведены результаты расчетов по комплексу GVA переходных процессов, обусловленных несанкционированным движением группы стержней в реакторе типа БН-800. Показано, что скорость извлечения стержней существенно влияет на характер поведения реактивности и мощности реактора. При движении стержней со скоростью 5 см/сек в реакторе на начальной стадии процесса не успевают отработать обратные связи, и значения введенной реактивности, вычисленные с учетом обратных связей, практически совпадают со значениями реактивности, рассчитанными без учета обратных связей. При скорости движения стержней 0.5 см/сек отрицательные обратные связи в значительной мере компенсируют вводимую реактивность. При моделировании процессов быстрого перемещения группы стержней наблюдается существенное формоизменение с течением времени пространственных нейтронно-физических полей и локальных энерговыделений в разных сечениях по высоте расчетной зоны реактора.
Как показали расчеты, применение в программе VOLNA разработанной автором эффективной сверхнеявной комбинированной схемы неполной факторизации с мерцающим параметром МНФ-МП на всех этапах, требующих решения систем конечно-разностных уравнений, существенно экономит время счета и делает возможным решать реальные пространственно-временные задачи на достаточно представительных интервалах времени.
Представленные результаты разработки алгоритмов и демонстрационных расчетов свидетельствуют, что комплекс программ GVA может служить составной частью создаваемого интегрального кода нового поколения для анализа безопасности РУ с реакторами на быстрых нейтронах. В этом случае применение новых эффективных схем метода неполной факторизации для решения сложных задач пространственной нейтронной кинетики становится особенно оправданным.
В заключении приведены основные результаты работы, сформулированы выводы.
- Новые эффективные варианты метода неполной факторизации, разработанные автором для решения алгебраических систем разностных уравнений эллиптического типа в (hex, z) –геометрии с использованием произвольной нумерации узлов разностной сетки, позволяют значительно экономить общее время счёта при решении задач нейтронной кинетики.
- Созданная библиотека подпрограмм IFML для решения систем линейных алгебраических уравнений на основе метода неполной факторизации дает возможность выбрать наиболее подходящие и эффективные в своем классе задач подпрограммы и использовать их в практических расчётах для различных физических приложений. Библиотека программ IFML имеет также общеобразовательное значение и служит цели сохранения знаний.
- Разработанные алгоритмы и программы нейтронно-физического расчета кампаний реакторов ВВЭР с использованием новых эффективных схем метода неполной факторизации позволили, в свою очередь, разработать алгоритм учета формоизменения ТВС в процессе эксплуатации реактора при проведении нейтронно-физических расчетов кампании с помощью комплекса программ WIMS-ВОЛНА. Результаты этих расчетов, выполненных во взаимодействии с группой разработчиков конструкции ТВС ВВЭР-1000 нового поколения, получивших название ТВС2 и ТВС2М, были использованы при расчетно-параметрических обоснованиях нового топлива и компоновок ТВС в активной зоне. Работы завершились созданием и широкомасштабным внедрением ТВС нового поколения на атомных станциях с реакторами ВВЭР-1000 и решением проблемы термомеханической нестабильности активных зон.
- Для нестационарного нейтронно-физического расчета переходных процессов в реакторах на быстрых нейтронах типа БН-800 с использованием новых схем метода неполной факторизации разработаны алгоритмы и программы, учитывающие обратные связи в случае заданного движения групп стержней в зоне реактора. Алгоритмы реализованы в рамках комплекса программ GVA нестационарного расчета реактора в квазистатическом приближении и апробированы на модельных задачах. Расчетные исследования и анализ, выполненные с помощью комплекса GVA, показали работоспособность комплекса и целесообразность совместного теплогидравлического и нейтронно-физического расчета при моделировании таких переходных процессов, как несанкционированное извлечение стержней регулирования.
основныЕ публикациИ по теме диссертации
Статьи в научных рецензируемых журналах:
- Гинкин В.П., Троянова Н.М., Новиков П.В. Новый вариант метода неполной факторизации для решения двумерных уравнений эллиптического типа с несимметричными матрицами // Математическое моделирование, т. 16, № 7, 2004, с.101-116.
Научные статьи в других изданиях:
- Гинкин В.П., Троянова Н.М. Использование метода неполной факторизации в трехмерной задаче нейтронно-физического расчета реактора ВВЭР-1000 // Всесоюзное совещание по динамике реакторов / Гатчина, 1990, 15 с.
- Гинкин В.П., Троянова Н.М. Вариант метода неполной факторизации для решения трехмерных девятиточечных уравнений эллиптического типа // 13-ый Международный конгресс по вычислительной и прикладной математике / Дублин, Ирландия. 1991. 6 с.
- Y.Likhachev, M.Meshkov, V.Troyanov, N.Troyanova. The results of thermo-mechanical calculation analysis of fuel assembly bowing in the WWER-1000 // International Conference «Safety analysis for NPPs of VVER and RBMK type», 26-30 October 1998, Obninsk, Russian Federation: Proceeding. – Vol.1, P.737-762.
- A.A.Bezborodov, A.V.Volkov, S.M.Ganina, V.P.Ginkin, I.A.Kuznetsov, N.M.Troyanova, Yu.E.Shvetsov. Spatial and Time-Dependent Calculation of Fast Reactor Transients // The Fourth International Congress on Industrial and Applied Mathematics, July 5-9, 1999, Edinburgh, Scotland.
- Troyanov V.M., Likhachev Y.I., Folomeev V.I., Demishonkov A.A., Troyanova N.M., Tutnov Al.A, Tutnov An.A, Kiselev A.S., Kiselev Al.S, Alekseev Е.Е., Ivanova O.I., Ulyanov A.I. Numerical and analytical investigation of WWER-1000 fuel assembly and reactor core thermal mechanics // Proceedings of a technical meeting held in Cadarache, France, 22- 26 November 2004, IAEA-TECDOC-1454, 2005, рр. 113-128.
- Гинкин В.П., Троянова Н.М. Использование метода неполной факторизации в трехмерной задаче нейтронно-физического расчета реакторов типа ВВЭР // 1990. Препринт ФЭИ-2104. 15с.
- Гинкин В.П., Троянова Н.М. Метод неполной факторизации для решения трехмерных девятиточечных разностных уравнений эллиптического типа // 1991. Препринт ФЭИ-2181. 17с.
- Гинкин В.П., Ваньков К.А., Троянова Н.М. Метод сопряженных градиентов в сочетании с неполной факторизацией и компенсацией итерируемых членов // 1993. Препринт ФЭИ-2324. 13с.
- Ваньков К.А., Гинкин В.П., Троянова Н.М. ВОЛНА – программа трехмерного нестационарного расчета реактора в квазистатическом приближении // Препринт ФЭИ-2360, 1994, 23 с.
- Гинкин В.П., Безбородов А.А, Волков А.В, Ганина С.М., Кузнецов И.А., Троянова Н.М., Швецов Ю.Е. Программа совместного решения уравнений пространственно-временного переноса нейтронов и теплогидравлических нестационарных и аварийных процессов в быстрых реакторах // Препринт ФЭИ-2637, 1997, 22 с.
- Лихачев Ю.И., Мешков М.Н., Троянов В.М., Троянова Н.М., Фоломеев В.И. Результаты расчетных исследований механизмов и причин искривления ТВС в активных зонах // Труды Международного Семинара «Деформация топливных сборок PWR и ВВЭР», Прага, Ржеж, Чехия, 17-20.02.1998.
- Гинкин В.П., Троянова Н.М. Алгоритм решения трехмерного нестационарного уравнения реактора методом неполной факторизации // Сб. научных трудов кафедры РКР ИАТЭ, вып.1 «Методы и средства моделирования физических процессов в ядерно-энергетических установках» под ред. Ю.А.Казанского / Обнинск: ИАТЭ, 1990, с.63-72.
- Гинкин В.П., Ваньков К.А., Троянова Н.М.Метод сопряженных градиентов и неполная факторизация с компенсацией // VI Российская конференция по радиационной защите ядерных установок // г.Обнинск, 1994.
- V.P.Ginkin, N.M. Troyanova, K.A.Vankov. A Quazystatic Approach to a Three-Dimensional Nonstationary Reactor Equipment Solution // Proceedings of the 9-th Topical Meeting on Problems of Nuclear Reactor Safety, Moscow, MEPhI, h/o «Volga», September 4-8, 1995, Vol.1, pp.63-65.
- Безбородов А.А., Волков А.В., Ганина С.М., Гинкин В.П., Кузнецов И.А., Троянова Н.М., Швецов Ю.Е Spatial and time-dependent calculation of fast reactor transients // «Consulting Meeting of IAEA specialists on fast reactor problem», 4 июня 1998 г., Обнинск, Россия.
- Гинкин В.П., Безбородов А.А., Волков А.В., Ганина С.М., Кузнецов И.А., Троянова Н.М., Швецов Ю.Е. Расчетные исследования нестационарных процессов в быстром реакторе, обусловленных несанкционированным извлечением из активной зоны компенсирующих стержней. // «Нейтроника-97», Обнинск, 28-30 октября 1997г.
- Безбородов А.А., Волков А.В., Ганина С.М., Гинкин В.П., Кузнецов И.А., Троянова Н.М., Швецов Ю.Е. Пространственно-временной расчет переходных процессов в быстрых реакторах // «Нейтроника-98», Обнинск, 27-29 октября 1998 г.
