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

□ОЗ466444

На правах рукописи

КСс/Яу— '

СЫЧУГОВА Елена Павловна

АЛГОРИТМЫ РЕШЕНИЯ УРАВНЕНИЯ ПЕРЕНОСА НЕЙТРОНОВ И ГАММА-КВАНТОВ В ЗАДАЧАХ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ ЯДЕРНЫХ РЕАКТОРОВ И ИХ ЗАЩИТЫ

05.13.18 Математическое моделирование, численные методы и комплексы программ

Автореферат диссертации па соискание ученой степени кандидата физико-математических наук

О 9 АПР 2009

Москва - 2009 год

003466444

Работа выполнена в Институте прикладной математики им. М.В. Келдыша Российской академии наук

Научный руководитель: кандидат физико-математических наук,

старший научный сотрудник Воронков Александр Васильевич

Официальные оппоненты: доктор физико-математических наук,

Зизин Михаил Николаевич.

кандидат физико-математических наук, старший научный сотрудник Аристова Елена Николаевна.

Ведущая организация: Физико-энергетический институт

им. А.И. Лейпунского, ГНЦ РФ ФЭИ.

Защита состоится « Сс ¿у^а^_2009 г. в и часов

па заседании диссертационного совета Д 002.024.02 при Институте прикладной математики им. М.В. Келдыша РАН по адресу: 125047, Москва, Миусская пл., 4.

С диссертацией можно ознакомиться в библиотеке Института прикладной математики им. М.В.Келдыша РАН.

Автореферат разослан » -/¿¿.Л/уУ»^ 2009 г.

Ученый секретарь диссертационного совета кандидат физико-математических паук

/¿^У ЩЕРИЦА О.В.

Общая характеристика работы

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

Актуальность темы.

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

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

Цель работы.

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

Достоверность результатов.

Достоверность полученных результатов подтверждается сравнительными численными исследованиями радиационных полей в защите реактора СВБР 75/100 в X-Y-Z и R-ç-Z геометрии, выполненными по разработанным программам в пакете «РЕАКТОР» и по известной программе TORT (США) с использованием одной и той же системы констант, а также

сопоставлением с результатами исследований, проведенных ранее в одномерной геометрии по другим методикам.

Новизна работы.

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

Практическая значимость работы.

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

Реализация и внедрение результатов работы.

Созданные программные модули К[N31), КШЗОб, К1ИКТ7. пакета «РЕАКТОР» переданы в ФГУП ОКБ «ГИДРОПРЕСС» РОСАТОМ'а для проведения массовых расчетов переноса нейтронов и гамма-квантов в задачах математического моделирования ядерных реакторов и их защиты в 3„Рт приближении. Массовые расчеты стали возможными благодаря тому, что в этих программах используются эффективные методы ускорения внешних и внутренних итераций.

Апробация.

Основные результаты диссертации докладывались и обсуждались на следующих семинарах и научных конференциях:

1. Семинар им. К.И. Бабенко ИПМ РАН (рук. В.К. Брушлинский);

2. Семинар Института математического моделирования РАН (рук. Е.И. Леванов);

3. 15-й семинар «Нейтроника-2004» - «Нейтроипо-физические проблемы атомной энергетики» г. Обнинск, 27-30 октября 2004 г.

4. IX Российская научная конференция «Радиационная защита и радиационная безопасность в ядерных технологиях» 24-26 октября 2006 г., Федеральное Агентство по Атомной Энергии, ГНЦ РФ Физико-Энергетический Институт им. А.И. Лейпунского, г. Обнинск.

5. 18-й семинар «Нсйтроника-2007» - «Нейтроипо-физические проблемы атомной энергетики» г. Обнинск, 30 октября - 2 ноября 2007 г.

6. 19-й семинар «11ейтроника-2008» - «Нейтронно-физичсские проблемы атомной энергетики» г. Обнинск, 28-31 октября 2008 г.

Структура и объем работы.

Диссертация состоит из введения, четырех глав, заключения, грех приложений и списка литературы. Материал диссертации изложен на 120 страницах, включает 41 рисунок, 11 таблиц и список литературы из 75 наименований.

Краткое содержание работы

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

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

Другой проблемой является ускорение сходимости внутренних итераций при решении неоднородной задачи расчета потока частиц в одной энергетической группе методом итераций по столкновениям [2]. Одним из распространенных методов является метод ребаланса. Первоначальная версия этого метода на мелкой сетке является устойчивой в узком диапазоне изменения пространственных шагов [4]. Когда этот метод был усовершенствован, и численные расчеты показали его эффективность,

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

В первой главе описаны две основные задачи математического моделирования ядерных реакторов. Математической моделью для описания переноса нейтронов в теории ядерных реакторов является линеаризованное уравнение Больцмана, которое записывается относительно функции потока нейтронов и гамма-квантов. В задаче без внешнего источника временная зависимость поля нейтронов описывается линейным уравнением вида 8Ч'/д1 = Лх¥. Предположим, что оператор А не зависит от времени и поэтому общее решение может быть записано в экспоненциальной форме Ч/ = е/". Можно предположить существование единственной функции Ч,та!1(?,у)>0 и такого вещественного числа а, что е'"Ч/тах(?,у) = ев'Ч'1пах(?,у), и для произвольной функции ч*0 (г, V) > 0 функция Ч'тах(7, V) является асимптотическим пределом решения с'"'1/(|(г,:у)/|е'"1Н0(»:/у)| -» Ч'тах(л,у) при

/->оо. В конечномерном пространстве эволюционный оператор ем аппроксимируется квадратной матрицей еш с неотрицательными элементами. По теореме Фробениуса [2] такая матрица имеет хотя бы одно неотрицательное собственное значение. Более того, из теоремы Перрона -Фробениуса [5] следует, что если неотрицательная матрица неразложима, то существует одно положительное собственное число Лт„=еаЫ, равное ее спектральному радиусу, и ему соответствует положительный собственный вектор 4*т„. (Матрица является неразложимой, если у нее все элементы отличны от нуля). В стационарном состоянии реактора а = О, Л„„ = 1 и им соответствует поток Ч,тм(?,у)> 0, который называется критическим.

Главной проблемой является задача нахождения критического потока, соответствующего стационарному состоянию реактора. После дифференцирования по времени выражения е'"Ч'т„(?,у) = Ч'1п„(г,у) для критического потока, соответствующего Лтю = 1, оно превращается в систему линейных однородных уравнений = для решения которой

используется метод итерации источника [2].

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

Л,У +—Л2Ч" = 0 или = ,

кф

где введен постоянный множитель и оператор А представлен в виде

суммы двух операторов Л = А, + Л2 так, что у оператора Л, существует обратный оператор и А, легко обращается, а оператор А=-А~'А2 обладает наибольшим по модулю простым собственным значением. Число называется эффективным коэффициентом размножения реактора. Величина Кг1/ показывает, как следует изменить оператор А2, чтобы получить стационарное решение при заданных параметрах реактора. В случае Кг]Г = 1 полученное решение совпадает с критическим потоком, являющемся нетривиальным решением ^„„(^.у) системы линейных однородных уравнений АЧ'тт(?,у) = 0.

Оператор Т = -А2А~' обладает лучшими свойствами, по сравнению с оператором А=-А~'А2. Он является интегральным, неотрицательным и действует па функции, заданные в трехмерном пространстве. Оба оператора имеют одинаковое наибольшее по модулю простое собственное значение Кг// = А, > |Я,|, 1 = 1,2,..., которому соответствует единственная (с точностью до положительного множителя) положительная собственная функция (> = А2Ч'. Поэтому лучше перейти от задачи поиска функции ^„„(Я.у) распределения потока нейтронов в шестимерном фазовом пространстве к задаче поиска функции <2(?) скорости генерации нейтронов деления путем решения эквивалентной системы уравнений переноса нейтронов: тд = Ке/Гд.

В диссертации приведена полная постановка однородной и неоднородной задач переноса нейтронов и гамма-квантов в мпогогрупповом приближении. Уравнение переноса частиц в многогрупповом приближении в некоторой пространственной области б запишем для энергетической группы £ = 1,2,...,МЗ в следующем виде:

(О ■ УЧ' * (?, П)) + 2 £ (?) ■ Ч"' (?, П) = Я1,у (?,С1). (1) Правая часть уравнения (1) для нейтронов группы # = 1,2,имеет следующий вид:

(г,О) = X |Р(г,Й • П')Ч'" (?,П')с/П' +Хр £ VI.1}(г)Ф* (?) + 5"(г), (2)

А=1 4т А-1

где введен скалярный поток нейтронов:

ф ц?) = — Гч,'(?,й)Л2. (3)

4гс

Для групп гамма-квантов Л', +1, Ы, + 2,+ Ыг = правая часть уравнения (1) имеет следующий вид:

к

4» 4»

На границе Г пространственной области с заданы нулевые значения углового потока для направлений ¿1 внутрь этой области:

где п - внешняя нормаль к границе области о. В (1) - (5) использованы следующие обозначения:

Й - единичный вектор в направлении полета частиц £2 = {0,а>) в трехмерной геометрии, где 0 - полярный угол между вектором и осью Z, со - азимутальный угол между его проекцией на плоскость Х-У и осью X; - число групп нейтронов;

Nу - число групп гамма-квантов;

N0 - полное число энергетических групп;

У-'^П) - плотность потока частиц в точке г в направлении П в группе g со скоростью Vй;

Ф*(?) - скалярный поток частиц в точке г в группе g со скоростью Vе ;

Ц(г) - полное макроскопическое сечение взаимодействия частиц;

- макроскопическое сечение рассеяния нейтронов из группы И в группу я;

Р*">8(г,пй') - макроскопическое сечение рассеяния гамма-квантов из группы Л в группу g;

Р*г''(г,пп') - макроскопическое сечение образования гамма-квантов в

группе g при столкновении нейтронов из группы и с ядром;

гХ'(7) - число нейтронов деления, возникающих при одном акте деления;

х1 - спектр деления нейтронов;

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

Система уравнений (1) - (5) описывает распределение нейтронов и гамма-квантов с учетом заданных внутренних источников, а также с учетом процесса деления. Если в правой части (2) член с делением отсутствует (у.=0 для всех g), система уравнений (1) - (5) описывает распределение нейтронов и гамма-квантов в зависимости от заданных источников (неоднородная задача). Задача на КеЛ (однородная задача) описывается системой уравнений (1) - (3) с граничными условиями (5) для нейтронов g = \,2,...,N„ с нулевыми внутренними источниками 5е(г) и множителем 1/АГе# перед вторым слагаемым в правой части (2).

Индикатриса рассеяния задана в виде ряда по полиномам Лежандра Рт до степени т, т.е. в Рт — приближении. Угловая зависимость решения описана набором дискретных направлений О, с соответствующими весами и>, для вычисления интегралов, стоящих в правой части уравнения переноса,

(5)

в Л'„ приближении метода дискретных ординат, где п - порядок угловой квадратуры (и = 2,4,...).

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

Многогрупновая система уравнений переноса нейтронов (1) для расчета Кгд может быть записана в операторном виде:

/л'= /*!'+ х1-мА'

к*

где Ь - оператор переноса, Р - оператор рассеяния нейтронов из верхних энергетических групп в нижние группы и внутри группы, х - спектральный оператор деления мгновенных нейтронов, /*" - оператор деления, м„ -оператор расчета нулевого момента от Т по угловой переменной (3), 'И -функция плотности потока нейтронов. Преобразуем эту систему к виду:

тй = кс1й, (6)

где <2 - функция плотности источника деления, () = /-"Л^Т, Т = 1<'М0(Ь-Р)

В конечномерном евклидовом пространстве Т - матрица, {) -собственный вектор, соответствующий собственному значению Кс,г. Компонентами вектора <2 являются значения источника деления в / - той пространственной ячейке 6, / = 1,...,/ где V, - ее элементарный

Л 1

объем. Матрица Т в левой части (6) удовлетворяет условиям теоремы Перроиа-Фробсниуса [5], т.е. она неотрицательна и неразложима. Тогда она имеет единственное простое положительное наибольшее собственное значение Я, = (А, > для всех к + \ ), которому соответствует собственная

функция <2 с неотрицательными компонентами, () ^ сЕ'ф' , где Ф* = Л/01Р*.

А-1

Степенной метод [6] для нахождения наибольшего собственного значения Кс/Г и соответствующего ему собственного вектора 0 имеет вид:

р = 1,2,...; *">=!, (7)

где начальное приближение дт рассчитывается по единичному начальному приближению скалярного потока ФЛ<0) = (1,...,1) во всех группах, т.е.

вг = ^^/¡^К, и параметр К(р) вычисляется по формуле:

/1=1

В методе простой итерации параметр К{р) вычисляется по формуле:

а"\ (9)

I-1

а начальное приближение источника деления @<0) и соответствующие ему скалярные потоки нейтронов Ф'(0) во всех группах задаются так, чтобы

1

сумма нейтронов по объему реактора равнялась единице, т.е. =1. Для

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

Для любого произвольного положительного начального приближения дт метод (7) - (8) сходится [6], [7] к решению:

/>-»00 p—i да

где с0 - произвольное положительное число, а Q - собственная функция. Сходимость гараитирована тем фактом, что Krff по модулю больше всех других собственных значений матрицы Т [6].

Справедливо следующее свойство включения [6], позволяющее оценить на каждой внешней итерации верхнюю и нижнюю границы величины Kef. Если определить значения и по следующим формулам:

и AT« = К^ттЦ^г, (10)

шах ( . q(p-|) > V )

то справедливы неравенства:

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

В случае «больших» задач, когда число пространственных точек 1 больше миллиона, итерационный процесс (7) и (8) или (9) сходится очень медленно. Часто возникают случаи, когда в ходе итераций относительная ошибка значения К<р) па двух соседних итерациях на один - два порядка меньше относительной ошибки расчета Kcg по его верхней и нижней границам:

еъа =К<")«ев =(К{Р1 и процесс медленно сходится. В этом случае можно говорить о том, что в ходе итерационного процесса получено приближение (псевдорешение),

соответствующее собственному значению К1''1, близкому к КсП. Возникает необходимость в использовании метода ускорения.

В диссертации описан новый метод - процесс) ускорения сходимости внешних итераций в задаче расчета и источника деления, используемый после того, как вычислены значения по (7) и К{г) по (8) или (9). Алгоритм ускорения заключается в следующем.

На итерации с номером р = 1,2,... после того, как вычислены К{р) и приближение , полагается Я(/1> = 0 и к{,,)=0. Вычисляется значение параметра Яр) по формуле:

| п^) _ о{р ^ II

-%-¡, (11)

где используется окгаэдрическая или евклидова норма. Если Я''"' <0 или Л<р1 > 1, то делается переход к расчету следующей внешней итерации. В противном случае для значений:

0<Я,"'<1 (12)

вычисляется параметр экстраполяции к>р):

я<" (13)

\-л{р)'

и проверяется выполнение условия: , тЩк^У" у-2))

<6, (14)

тах(кг(/') ,к{р~1]) где значения *-„,„„ и б являются параметрами метода. Если условие (14) выполняется, то делается экстраполяция с использованием приближений 2</>-о и Qi.pt по формуле:

е=е(')+^)(е<')-ё<"-,)) (15)

и за приближенное значение <2{р) принимается вектор:

(16)

Одновременно уточняются скалярные потоки Фе = М0Ч'г во всех группах с номерами g = \,2,...,N„ по формулам:

Ф^Ф^+^ЧФ^'-Ф^"). (17)

Угловые моменты группового потока с индексами V > 1 и и>0, ы<V<тк переопределяются пропорционально изменению скалярных потоков:

(к) (*) 1р> фЯ (Я) (Я) (Р> ф*

Ф1"'=Ф1". -^Г' ф2".=ф2". -^рг-

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

групповых скалярных потоков Фк в каждой точке пространства но формулам

(11) - (17) и угловых моментов потока по формулам (18).

Алгоритм ускорения сходимости внешних итераций эффективен в случае удачно подобранных критериев. Из (13) видно, что параметр экстраполяции к^ может быть очень большим положительным числом, поэтому значение необходимо ограничивать сверху. Обычно гта1 =3000. Оптимальное значение параметра 5 зависит от решаемой задачи. Из неравенства (14) следует, что значение параметра 5 находится в интервале 3 6 (0,1). В задачах физики реакторов 3 = 0.1 ^ 0.8.

Когда Я{р> « 1, т.е. вблизи точного решения, 3 - процесс не позволяет делать экстраполяцию (15) - (18). В этом случае значения Я{р) на трех последовательных итерациях не удовлетворяют какому-нибудь из неравенств

(12) или (14). Поэтому итерационный процесс заканчивается без ускорения.

Численные расчеты показали, что 3 - процесс сходится быстрее

степенного метода (7), (8) или метода простых итераций (7), (9) к наибольшему собственному значению и соответствующему

собственному вектору 2 неотрицательной неразложимой матрицы Т.

В диссертации приведены результаты тестовых расчетов Ке1 критической сборки 001)1УА [8] в одномерной сферической геометрии. Экспериментальное значение =1.000±0.001. В задаче заданы 30-ти групповое приближение по энергии нейтронов, угловая квадратура Гаусса-Лежандра I] приближение индикатрисы рассеяния и сетка из 800 интервалов. В качестве начального приближения заданы единичные значения скалярного потока. Для окончания внутренних итераций задана величина максимальной относительной ошибки скалярных потоков £„,=106 и максимальное число итераций в группе 20. Критерием окончания внешних итераций служит одновременное выполнение условий ек1Г<\0~6 и £й<10~6. Для ускорения сходимости внешних итераций использовался 3 - процесс с вычислением параметра Яи'' по октаэдрической норме.

Результаты расчетов без ускорения и с ускорением для различных значений параметра в совпали. Получено значение Ке1 = 0.99990. Расчеты показали, что наилучшее ускорение происходит при 5 = 0.9 в 1.6 раза по общему числу итераций (11 итераций с ускорением и 18 без него). На рис. 1 и 2 показаны в логарифмическом масштабе значения К(р' и относительных ошибок сИ] и Ев, полученные на внешних итерациях без ускорения и с ускорением с параметром 5 = 0.9, соответственно. Из этих рисунков видно, что при расчете без ускорения порядок ошибки убывает, как линейная функция, а при расчете с ускорением ближе к параболической. На рис. 3 и 4 показано соответствующее поведение параметра Л(р}. Из рис. 3 видно, что при расчете без ускорения с 3-й по 16-ю итерацию Л'р) практически не меняется, т.к. «мешает» псевдорешение, соответствующее значению

Я(р) к 0.44. На рис. 4 символом «*» показаны 5-я, 6-я и 10-я итерации, на которых использовалась линейная экстраполяция пссвдорсшения.

Расчет без ускорения показал, что пссвдорсшение, соответствующее собственному значению Л00 ю 0.44, препятствует быстрой сходимости внешних итераций.

Рис. 1. Изменение Klp), cícff и cQ без ускорения.

2 3 4 5 в 7 8 9 10 11 12 13 14 15 1в 17 1 Iteration р

Рис. 2. Изменение Kip), ricjr и r,Q с ускорением, 3 = 0.9.

-Л

-1—■ ' ■ ■_I_

3 4 5 в 7 в 9 10 11 12 13 14 15 1в 17 10 Iteration р

Рис.3. Изменение Л1р) без ускорения.

экстраполяция

]' ' > ■_i_i_■ у ■ ■_i_i_i_i_i_i_i_i_i

0 1 2 3 4 5 б 7 в 9 10 11 12 13 14 15 1в 17 1В Iteration р

Рис.4. Изменение Áíp> с ускорением, 3 = 0.9.

В диссертации приведены результаты расчетов и источника деления исходного состояния критической сборки BZD/1 в экспериментах «ZEBRA» [9] в X-Y-Z геометрии в 30-ти групповом приближении по энергии нейтронов для S4P, и SsP, приближений метода дискретных ординат на сетке, состоящей из 24304 точек, которые показывают эффективность 5 - процесса. Приведены соответствующие таблицы с числом внешних итераций при 5 = 0.2 в зависимости от формулы, используемой для расчета параметра Л<р): с октаэдрической нормой или евклидовой нормой. Из расчетов следует, что в StP¡ приближении эффективнее расчет с евклидовой нормой, а в SsP,

приближении - расчет с октаэдрической нормой. Выигрыш примерно в 2.2 раза по числу итераций. Сделан вывод о том, что любая из двух норм для расчета параметра Яы предпочтительнее формулы З2 - процесса [10].

В таблицах 1 и 2 приведены результаты расчетов с различными значениями 3 от 0.1 до 0.9, вычисленными по октаэдрической норме.

Таблица 1. Число внешних итераций в S4P, приближении при расчете Kt, задачи ZEBRA с ускорением для различных значений 5 и без него.

3 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Без ускорения

Число внешних итераций 73 33 35 32 29 34 26 28 37 77

Коэффициент ускорения 1.1 2.1 2.0 2.2 2.4 2.2 2.7 2.6 2.1 1.

Таблица 2. Число внешних итераций в SJ\ приближении при расчете Kt

задачи ZEBRA с ускорением для различных значений 3 и без него.

5 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Без ускорения

Число внешних итераций 34 30 34 32 26 31 27 23 38 79

Коэффициент ускорения 2.0 2.2 2.0 2.1 2.5 2.3 2.7 3.0 2.1 1.

Из таблиц видно, что наилучшее ускорение в 3 раза по общему числу итераций получено в £8Р3 приближении при 3 = 0.8 (23 итерации с ускорением и 79 без него). При 3 = 0.7 получено одинаковое ускорение в 2.7 раза в обоих приближениях. Это означает, что параметр 3 можно выбрать, проведя расчеты в более низком угловом приближении ив приближении индикатрисы рассеяния Г]. Поэтому для эффективного расчета однородной задачи на КгИ в более высоком приближении нужно сначала выбрать Зщл оптимальное значение параметра 3 путем проведения расчетов в более низком приближении, например и затем использовать его для более точных расчетов.

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

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

В диссертации описан алгоритм решения уравнения переноса в одной группе методом итерации по столкновениям [2] с использованием «взвешенной» схемы [2] в трехмерной Х-У-Х геометрии. Соотношение баланса частиц в ячейке [*,_|/2.*,ч|/2] х [уНП'У^мг] х К Для

направления £2т =(/'„, ) получается интегрированием уравнения переноса (1) с правой частью Я" (г, £2) = дх (г, О) I Л'" (Ч'х (г, С})) по объему ячейки и по направлениям из углового диапазона у/„ вокруг направления йт:

кИте. - -Фгл+^ч» =щ,м.т. (19)

В правой части (19) член с внутригрупповым рассеянием отделен от остальных слагаемых, а символом «~» для краткости обозначен промежуточный результат на итерации д. Направление полета частиц учтено в знаках величин Для записи результата использовано

приращение индекса, определенное с помощью функции sign(x), которая принимает значения ±1 в соответствии со знаком ее аргумента: а = 1/2л7£я(//„); Ь = 1 /); с = \12хщп(1;1п). Коэффициенты Л = Ьу1Ьг1, В = Лх Лг,, С = Л*¡AyJ - площади боковых поверхностей ячейки, V = Лг,ЛууЛг, -объем ячейки.

Система уравнений для поправок получается следующим образом. Умножим (19) на угловой вес и просуммируем по всем направлениям £2п. В результате получим следующее уравнение:

7+? 7-ч ,7-ч 7+« , 7+« 7-« .7-« 7+«

и /и/ш —и /и/ад +•' ¡-чал н/ад —и у+| /у +и ц-чхк /л/у +

~ , (20)

+гЬми 2 -гЬм +гЬл-ч2 -^цмп+уъД^ =удм

где

слагаемое J*шlг,j,^-J'шl2,j,^ представляет собой ток из рассматриваемой ячейки через грань х11112. Поправки вводятся так, чтобы восстановить баланс для промежуточного результата итерации с номером q, предполагая, что токи, выходящие из рассматриваемой ячейки через ее грани, пропорциональны скалярному потоку, который образуется в этой ячейке, и коэффициент пропорциональности равен :

Ум/и* и*%чгил = и*1и*чи =^ЬмигМ*Ьми1 (21)

В результате получим уравнение для поправок /л,, справедливое в каждой внутренней ячейке:

+7~1иш (22)

+Г1н,и +7'Ь+Ш + УЪмп +ТЬмиг -Ф&)] = Щ,Д

Для решения системы диффузионно-подобных уравнений (22) относительно поправок /,», используется итерационный метод верхней релаксации [2]. Точным решением является единичная поправка , если

на итерации по столкновениям (19) в каждой внутренней точке потоки до и после итерации совпадают ф* ,=фИтерация ц решения уравнения (19) заканчивается перенормировкой (21) скалярных потоков:

(23)

Угловые моменты потока также нормируются умножением на /,»_,, что позволяет сохранять угловое распределение промежуточного результата итерации ц. Для повышения эффективности ускорения используется алгоритм [12] увеличения фиктивных токов на границе каждой ячейки.

Исследована устойчивость метода пространственного ребаланса с помощью Фурье-анализа его линеаризованного приближения совместно с «алмазной» схемой [2] в одномерной плоской геометрии на равномерной сстке. Устойчивость этого метода ранее не исследовалась. Рассмотрен случай изотропного рассеяния частиц на специальном классе задач с постоянными сечениями и постоянным изотропным источником в бесконечной области. Точным решением этой задачи является постоянный скалярный поток

Полная система конечно-разностных уравнений в одномерной геометрии для решения этой задачи с ускорением методом пространственного ребаланса в Л*„ приближении метода дискретных ординат на равномерной сетке по переменной х с шагом А = хУ(„2-хнп может быть записана на

внутренней итерации с номером q в следующем виде:

кК^.,,» -Ф:., = Аб+Иад®Г'. (24)

= 1 + ), (25)

м

(26)

т=1

7*%п = ЕЧ'+ии.М.». , = X *'-./*»- , (27)

т т

//„> 0 //я<0

Т?+|/2=^-т1п(7+у+|/2,/>.1/2), (28)

Т]-\12 =J*Чj-U2+JЧj_^/2, У+1 /2 у+1/2+^+)/2, (29)

Ф; =/;•$}, (31)

гДе - направляющие косинусы и и>ш - веса квадратуры, т = 1,2,..., М . В правой части (24) слагаемое с впутригрупповым рассеянием определяется по приближению скалярного потока на предыдущей итерации Ф* 1. В эту систему входит уравнение (25) «алмазной» схемы. По квадратурной формуле (26) вычисляется промежуточный скалярный поток Ф*, который в случае отсутствия ускорения подставляется в правую часть (24) для расчета на следующей итерации. При расчете с ускорением вычисляются односторонние токи 3*)±м1 по формулам (27), которые увеличиваются но формулам (28) и (29) с множителем <1 на каждой внутренней итерации. Увеличенные токи являются коэффициентами системы конечно-разностных уравнений (30) относительно поправок /'. Итерация </ заканчивается перенормировкой (26) скалярных потоков.

Для исследования устойчивости метода пространственного ребаланса (24) - (31) предполагается следующее приближение скалярного потока вблизи точного решения:

Ф*"'=(1 + (32)

где £« 1. Будем считать, что промежуточные значения углового и скалярного потока на итерации q имеют аналогичный вид, а поправки имеют линейный вид:

/;=!+*•*;. (зз)

Тогда можно оцепить устойчивость линеаризованного приближения, определяя, куда стремится линейная добавка к решению Ф^ = (?/£_,: к пулю или к бесконечности. Подставим (32) и (33) в уравнения (24) - (31) и приравняем члены при одинаковых степенях с. Члены порядка 0(1) автоматически удовлетворятся. Пренебрегая членами порядка 0(с2), получим систему линейных уравнений с постоянными коэффициентами для линейных добавок к решению, к которой можно применить Фурье-анализ устойчивости:

+ I.,- = h■LйJ^^'\ (34)

(35)

— м

*)= (36)

~ №%, О + с!) + 8Ч1 + 2/0 + <!)) ~ Г ■ 8]-, 0 + <0 = ^ (?/ - ф]''), (37)

(38>

где введено значение у = ^А.*'« -1/4.

11т>а

Для анализа устойчивости системы (34) - (38) применяем дискретное преобразование Фурье по L гармоникам (L<J, J - число пространственных точек):

= = I,

J=о

где Л, - амплитуды гармоник, которые суммируются в точках xj = nY.Tj/J, j = 0,1,..., J-\. Введем следующие гармоники рассматриваемых функций:

ф]~' = ХГе"' , Vh = ХГате*', ^ = X?bme"» , ф< = X^F^', g) = Я^ве"'. (39) Подставляя (39) в (34) - (38) и решая полученную систему уравнений, найдем ф] =Л, -ф] ', где амплитуда Л, определяется по формуле:

ФУ w" -1)

Л,-С±_^_+ . (40)

ti\+i-2fimtgT,III \-c+Ay(\+d)s\x\ т,1 II v

В (40) введен параметр т,=1Н/2, где Н = я2.ти. Амплитуда Л, определяет скорость, с которой гармоника решения с номером / усиливается при |Л,|>1 или затухает при |Л,|<1. Нас интересует спектральный радиус метода пространственного рсбаланса для заданной квадратуры Sn. Из формулы (40) следует, что амплитуда нулевой гармоники равна нулю, а остальные амплитуды являются периодическими функциями / с периодом л/II. Поэтому достаточно рассмотреть их изменение только на интервале 0<г, <я72 для I = 1,2,...,L-\, чтобы оценить значение спектрального радиуса:

р=шах supU,(r)|.

0<r<W2 '

В таблице 3 показаны результаты исследования поведения модуля амплитуд (40) при d = 0 и с = 105/ Хг =1 для различных значений L от 1 до 6

в области максимальных значений шага // = я7£. При вычислении \Л,\ значения //„ и wm соответствовали квадратуре S4 Гаусса - Лежандра. В этой таблице приведены: общее число гармоник L, шаг Я, значения первых пяти амплитуд и номер амплитуды /mix, па которой достигается максимальное значение.

Таблица 3. Значения модуля амплитуд |д,л""'(Я)| для L = 1,...,6 и Н = л! L.

L // = 71 i L \л?приЦ |дГ(я)| /„„

1 3.1415 л л/2 л/3 л/4 л/5 1

2 1.5708 0.17951 л/2 л/3 л/4 л/5 2

3 1.0472 0.25079 0.43796 л/3 л/4 л/5 3

4 0.78540 0.45181 0.31867 0.52379 л/4 л/5 4

5 0.62832 0.68696 0.35530 0.35657 0.52039 л/5 1

6 0.52360 0.93611 0.43889 0.32583 0.38069 0.47759 1

Из таблицы 3 видно, что за исключением трех случаев, максимальное значение у первой амплитуды. Результаты численного исследования показали, что начиная с ¿>5 (при II </г/5) максимальное значение всегда достигается на первой амплитуде. Поэтому достаточно рассматривать первую амплитуду в качестве спектрального радиуса метода, т.к. при решении уравнения переноса грубые сетки из менее пяти интервалов не используются. Исследования показали, что при близких к 1 значениях с для любого фиксированного значения II увеличением параметра (1 можно добиться того, чтобы радиус сходимости метода стал меньше единицы.

На Рис. 5 показано в логарифмическом масштабе изменение спектрального радиуса р в зависимости от (I для шага Я =/г/1000. Из рисунка видно, чем меньше пространственный шаг, тем больше значение множителя (1, при котором р< 1. Это означает, что увеличение коэффициентов уравнения (30) для поправок приводит к уменьшению спектрального радиуса сходимости линеаризованного метода пространственного ребаланса до значений р< 1, при которых этот метод устойчив.

Рис. 5. Зависимость спектрального радиуса р от параметра с! для II = /¡-/1000,

Таким образом, с помощью Фурье-анализа доказана устойчивость линеаризованного метода пространственного ребаланса совместно с «алмазной» схемой на равномерной сетке в одномерной плоской задаче в приближении метода дискретных ординат для бесконечной области при изотропном рассеянии частиц с постоянными сечениями и постоянным изотропным источником.

В диссертации приведены результаты тестовых расчетов железоводной борированной защиты в одномерной геометрии с использованием в^ГШ схемы [15] совместно с методом пространственного ребаланса. Для сравнения эффективности ускорения приведены результаты расчетов этой задачи ЛиТЮ схемой совместно с Р^А методом ускорения, выполненные ранее в пакете «РЕАКТОР» [16]. Результаты расчетов показали, что метод пространственного ребаланса эффективнее Р^А метода и обеспечивает

р

а

Р™„№ = 0,1545.

существенный вычислительный выигрыш примерно в 5 раз по времени, в то время как P,SA метод дает меньший выигрыш по времени в 3 - 4 раза.

В третьей главе описаны трехмерные программы метода дискретных ординат с ускорением сходимости внутренних и внешних итераций, представляющие собой модули пакета «РЕАКТОР» и предназначенные для расчета значений Kejr и источника деления, а также задач защиты. Приведено описание каждого модуля. Дана общая блок-схема модулей и их основной подпрограммы DSNPN. Модуль KIN3D предназначен для расчета потока частиц в X-Y-Z геометрии в выпуклых областях, горизонтальное сечение которых задается на равномерной квадратной сетке. Модуль KINRTZ предназначен для расчета потока частиц в R-(p-Z геометрии в цилиндрических областях с неравномерной сеткой. Модуль KIN3D6 предназначен для расчета уравнения переноса нейтронов в HEX-Z геометрии в выпуклых областях, горизонтальное сечение которых задается на сетке, состоящей из правильных шестиугольников.

В четвертой главе приведены результаты решения задач математического моделирования проектируемого в настоящее время ядерного энергетического реактора СВБР 75/100 (Свинцово-Висмутовый Быстрый Реактор) [17] и его защиты в трехмерной геометрии. Продемонстрирована эффективность предложенных алгоритмов и созданных программ.

В первой части главы приведены результаты расчетов Keïï и источника деления реактора СВБР 75/100 с использованием схемы с «нулевой» коррекцией в HEX-Z геометрии в 30-ти групповом SsP, приближении по энергии нейтронов для различных значений параметра S от 0.2 до 0.9. Пространственная сетка состояла из 54 плоскостей по вертикальной оси Z и из 22789 точек (гексагональных ячеек) в каждой плоскости OXY. Более высокое приближение индикатрисы рассеяния SK /' рассчитано с параметром S = 0.7, при котором результаты расчетов в S^P, приближении были в 2 раза эффективнее по времени счета, чем расчет без ускорения. По результатам расчетов разработаны рекомендации о выборе значения параметра 5.

Во второй части главы приведены результаты расчета радиационных полей в защите реакторной установки СВБР 75/100 (30 нейтронных групп и 19 гамма групп) с использованием ©-WDD схемы [15] в трехмерной геометрии. В X-Y-Z геометрии пространственная область представляет собой четвертую часть реакторной установки в виде прямоугольного параллелепипеда 252.5 см по оси X и по оси Y и 752 см по оси Z. Равномерная сетка состояла из 101 интервала по X и по Y и 188 интервалов по Z, что составляло 1917788 ячеек. Использована полностью симметричная квадратура Ss и Ръ - приближение индикатрисы рассеяния. Константы рассчитаны по программе TRANSX 2.0 (США) [18]. Получено хорошее совпадение результатов расчетов плотности полного потока нейтронов с аналогичными результатами, полученными по программе TORT [19].

На основании проведенной серии расчетов задачи защиты реактора СВБР 75/100 сделан вывод о том, что на эффективность ускорения внутренних итераций методом пространственного ребаланса влияют в основном два фактора: выбор номера итерации, с которой начинается ускорение, и выбор параметра 0 из интервала (0,1) схемы 0-\УОО. Проведена серия расчетов для оптимального выбора этих параметров. Для анализа результатов расчетов и оценки точности алгоритма выбраны три характерные точки области расчета, в разной степени удаленные от источника. Результаты расчетов показали, что разброс средних значений полного потока находится в пределах 20%. При ослаблении потока в 1015 раз такой уровень разброса значений считается приемлемым. Таким образом, при расчете задачи защиты можно использовать любые значения параметра 0 для получения приемлемого решения. Это важно, т.к. чем больше значение 0, тем больше влияние осцилляций, мешающих сходимости итераций. С точки зрения величины времени счета эффективной оказалась схема с параметром 0 = 0.3 и с ускорением во всех группах, начиная с первой итерации.

Аналогичные результаты были получены в Я-ф-/ геометрии. Рассчитано ЗД - приближение. Использована сетка из 110 интервалов по радиусу, 120 интервалов по углу ф и 137 интервалов по оси Ъ, что составляет 1808400 ячеек.

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

ЗАКЛЮЧЕНИЕ

Предложен новый метод (6 - процесс) ускорения сходимости внешних итераций при решении задач расчета эффективного коэффициента размножения и источника деления. Реализованы эффективные и адаптивные методики ускорения сходимости итераций для решения уравнения переноса нейтронов и гамма-квантов в задачах математического моделирования ядерных реакторов и их защиты. Эти методики универсальны, для их применения не требуется предварительных сведений о поведении и свойствах решения. Они могут быть использованы для решения задач математического моделирования ядерных реакторов различного типа и назначения и для расчета их защиты в трехмерной геометрии. Методы ускорения требуют на каждой итерации относительно небольшого количества дополнительных арифметических действий. 3 - процесс может быть использован при решении больших задач поиска наибольшего собственного значения неотрицательной неразложимой матрицы.

Выполнены расчеты проектируемого в настоящее время энергетического ядерного реактора СВБР 75/100 в трехмерной геометрии по разработанным программам решения уравнения переноса методом дискретных ординат.

ОСНОВНЫЕ РЕЗУЛЬТАТЫ

1. Предложен новый метод {8 - процесс) ускорения сходимости внешних итераций при решении задач расчета эффективного коэффициента размножения и источника деления, возникающих при математическом моделировании ядерных реакторов методом дискретных ординат в трехмерной геометрии.

2. Исследована устойчивость линеаризованного приближения метода пространственного ребаланса па одногрупповой задаче с изотропным рассеянием в бесконечной однородной среде с помощью Фурье-анализа. Доказана его устойчивость.

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

4. Проведены численные исследования па тестовых и реальных задачах эффективности реализованных методов ускорения сходимости внутренних и внешних итераций. Разработаны рекомендации по выбору параметра 8 для эффективного ускорения внешних итераций при расчетах в SnPm приближении. Разработаны рекомендации по оптимальному выбору параметра в взвешенной 0-WDD схемы.

5. Выполнены трехмерные расчеты проектируемого в настоящее время энергетического ядерного реактора СВБР 75/100, которые были использованы в проектных работах.

Публикации в реферируемых журналах.

1. Сычугова Е.П. Исследование устойчивости и эффективности метода пространственного ребаланса для ускорения сходимости итераций в задачах переноса частиц // Математическое моделирование, 2008, Т. 20, №9, С. 75-93.

2. A.V. Voronkov, Е.Р. Sychugova. А second-order finite volumc discretization of the time-dependent transport equation on arbitrary quadrilaterals in R-Z geometry // Nuclear Science and Engineering, American Nuclear Society, 2004, Vol. 148, №1,P. 186-194.

3. A.V. Voronkov, E.P. Sychugova. CDSN-Method for Solving the Transport Equation // Transport Theory and Statistical Physics, 1993, Vol. 22, Iss. 2-3, P. 221-245.

Другие публикации но теме диссертации.

4. Е.П. Сычугова. Численные методы решения уравнения переноса в многогрупповом приближении в трехмерной геометрии в пакете «РЕАКТОР» // Препринт ИПМ РАН, 2007, №78, С. 1-24.

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

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

Цитируемая литература

1. Marvin L. Adams, Edward W. Larsen. Fast iterative methods for discrete-ordinates particle transport calculations // Progress in Nuclear Energy, 2002, V. 40, Issue l,p. 3-159.

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

3. J1.A. Люстерник. Замечания к численному решению краевых задач уравнения Лапласа и вычислению собственных значений методом сеток // Труды математического института им. Стеклова, 1947, Т. 20, С.49-64.

4. G. Cefus, E.W. Larsen. Stability Analysis of Fine-Mesh Rebalance // Trans. Am. Nucl. Soc., 1988, V. 56, P. 309-310.

5. В.В. Воеводин, Ю.А. Кузнецов. Матрицы и вычисления // Паука, М., 1984, с. 130.

6. Варга Р.С. Численные методы решения многомерных миогогрупповых диффузионных уравнений // В кн.: Теория ядерных реакторов. Пер. с англ. под ред. Г.А. Батя, Госатомиздат, М., 1963.

7. Шишков Л.К. Методы решения диффузионных уравнений двумерного ядерного реактора // Атомиздат, М., 1976.

8. International Handbook of Evaluated Criticality Safety Benchmark Experiments // NEA/NSC/DOC(95)03.

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

10.H.C. Бахвалов. Численные методы // Изд-во «Наука», М., 1973.

11.W. W. Engle Jr., F. R. Mynatt. A Comparison of Two Methods of Inner Iteration Convergence Acceleration in Discrete Ordinates Codes // J. Trans. Am. Nucl. Soc., 1968, V. 11, P. 193-194.

12.W. H. Reed. The Effectiveness of Acceleration Techniques for Iterative Methods in Transport Theory // J. Nucl. Sci. Eng., 1971, V. 45, P. 245-254.

13.W. A. Rhoades, R. L. Childs, W. W. Engle Jr. Comparison of Rebalance Stabilization Methods for Two-Dimensional Transport Calculations // J. Trans. Am. Nucl. Soc., v. 30, 1978, p. 583.

14.W.A. Rhoades. Improvements in Discrete Ordinates Acceleration // J. Trans. Am. Nucl. Soc., 1981, V. 39, P. 753.

15.W. A. Rhoades, W. W. Engle. A New Weighted-Difference Formulation for Discrete Ordinates Calculations // J. Trans. Am. Nucl. Soc., 1977, V. 27, P. 776-777.

16.A. Voronkov, V. Arzhanov. "REACTOR" - program System for Neutron -Physical Calculation" // Proc. Int. Top. Meting: Advances in Mathematics, Computations and Reactor Physics, 1991, V. 1, April 28 - May 2, Pittsburg, USA.

17.Yu. G. Dragunov, V. S. Stepanov, N. N. Klitnov, A. V. Dedul, S. N. Bolvanchikov, A. V. Zrodnikov, G. I. Toshinsky, O. G. Komlev. Project of SVBR-75/100 reactor plant improved safety for nuclear sources of small and medium power // 5lh International Conference on Nuclear Option in Countries with Small and Medium Electricity Grids, Dubrovnik, Croatia, 2004, May 1620.

18.R. E. MacFarlane. "TRANSX-2: Code for Interfacing MATXS Cross-Section Libraries to Nuclear Transport Codes" // LA-12312-MS, 1992.

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

Подписано в печать 23.03.2009 г. Печать лазерная цифровая Тираж 100 экз.

Типография Aegis-Print 115230, Москва, Варшавское шоссе, д. 42 Тел.: (495) 785-00-38 www.autoref.webstolica.ru

ВВЕДЕНИЕ

На защиту выносятся результаты.

ГЛАВА I. Задачи математического моделирования ядерных реакторов и их защиты.

1.1. Математическое моделирование ядерных реакторов.

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

1.3. Приближение индикатрисы рассеяния полиномами Лежандра

Рт — приближение).

1.4. Угловые квадратуры Sn и ESn метода дискретных ординат.

ГЛАВА II. Итерационные методы решения.

2.1. Ускорение сходимости внешних итераций.

2.1.1. 8 - процесс ускорения сходимости итераций в задачах расчета значений Keff и источника деления.

2.1.2. Исследование скорости сходимости 8 - процесса в задаче поиска наибольшего собственного значения неотрицательной неразложимой квадратной матрицы А

2.1.3. Результаты тестовых расчетов критической сборки GODIVA в одномерной сферической геометрии.

2.1.4. Результаты расчета значений KcjJ и источника деления исходного состояния критической сборки BZD/1 в экспериментах «ZEBRA» в Х-Y-Z геометрии.

2.2.Ускорение сходимости внутренних итераций.

2.2.1. Дискретизация уравнения переноса и метод решения.

2.2.2. Метод пространственного ребаланса.

2.2.3. Фурье анализ устойчивости метода пространственного ребаланса совместно с «алмазной» схемой на примере одногрупповой плоской задачи для изотропного рассеяния в бесконечной среде.

2.2.4. Результаты тестовых расчетов железоводной борированной защиты в одномерной геометрии.

ГЛАВА III. Реализация метода дискретных ординат и эффективных методов ускорения сходимости в трехмерных модулях KIN3D, KINRTZ и KIN3D6 пакета

РЕАКТОР» для моделирования ядерных реакторов и их защиты.

3.1. Описание модуля KIN3D.

3.2. Описание модуля KINRTZ.

3.3. Описание модуля KIN3D6.

3.4. Входная информация к модулям.

3.5. Выходная информация модулей.

3.6. Информация для продолжения работы модулей.

3.7. Общая блок-схема управляющей программы KIN модулей KIN3D,

KINRTZ и KIN3D6 и основной подпрограммы DSNPN.

ГЛАВА IV. Проверка эффективности предлагаемых алгоритмов на математических моделях реальной ядерно-энергетической установки.

4.1. Расчет значений Keff и источника деления реактора СВБР 75/100.

4.2. Расчет радиационных полей в защите реактора СВБР 75/100.

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

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

В течение более 20-ти лет в Институте Прикладной Математики им. М.В. Келдыша РАН под руководством A.B. Воронкова проводились работы по созданию численных методов и программ решения стационарных и нестационарных систем многогрупповых уравнений переноса частиц в приближении метода дискретных ординат [1] с учетом структуры расчетной области в различных геометриях [2] - [37]. Использование таких программ позволило в 1994 году впервые в нашей стране провести трехмерные расчеты уравнения переноса [8] для получения натриевого пустотного коэффициента реактивности критической сборки BZD/1 [9] - [10] в экспериментах «ZEBRA» [38] в SA приближении метода дискретных ординат, в Р{ приближении индикатрисы рассеяния и в 26 групповом разбиении по энергии. Численное решение получено методом итерации источника [1] без использования алгоритмов ускорения сходимости итераций в рамках пакета прикладных программ «РЕАКТОР» [39]. Пакет «РЕАКТОР» разработан по модульному принципу и использовался ранее для математического моделирования ядерных реакторов на быстрых нейтронах в диффузионном приближении. Поэтому все программы решения многогрупповых систем уравнений переноса оформлены в виде модулей. В состав пакета «РЕАКТОР» входят модули различного типа и назначения. Обмен информацией между модулями осуществляется через общий Архив. Системное обеспечение пакета служит для создания Архива данных, ввода информации и сопровождения работы каждого модуля. Имеются геометрические модули, модуль ввода ядерно-физических данных, модули расчета констант, модуль задания параметров рассчитываемой задачи, модуль обработки результатов расчета и др. Большинство модулей используются до расчета решения систем многогрупповых уравнений переноса частиц, подготавливая и записывая в Архив исходные данные.

Для решения уравнения переноса нейтронов и гамма-квантов в одномерной геометрии в пакете «РЕАКТОР» использовалась положительная AWDD схема метода дискретных ординат совместно с согласованной PXSA [И] схемой ускорения сходимости внутренних итераций. Пакет программ KINRZ дополнен двумя методами коррекции «алмазной» схемы [1] для получения положительного решения: схемой с «нулевой» коррекцией потока (DDO схема) и методом коррекции по шаговой схеме (DDST схема) [1], программами для расчета индикатрисы рассеяния в Рт приближении и DSA — методом ускорения внутренних итераций [16]. Разработан алгоритм [17] и создана программа решения многогруппового стационарного уравнения переноса нейтронов и гамма-квантов на сетках, согласованных со структурой расчетной области в двумерной R-Z геометрии. Разработана нестационарная методика решения уравнения переноса в R-Z геометрии на сетках, состоящих из произвольных выпуклых четырехугольников [18], [20] и получены результаты тестовых расчетов. Для проведения расчетов активной зоны ядерного реактора на быстрых нейтронах с гексагональной решеткой ТВЭЛов разработан метод и создана программа решения уравнения переноса в HEX-Z геометрии [21] с использованием DDO схемы на сетках, максимально согласованных с конфигурацией расчетной области. Первоначально написанные в однопроцессорном виде трехмерные программы решения уравнения переноса перенесены на параллельные ЭВМ кластерного типа [22] - [26] в пакете «REACTOR-P».

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

Наиболее полный обзор современных методов ускорения сходимости итераций для решения уравнений переноса частиц методом дискретных ординат дан в работе [40]. Плохую сходимость внешних итераций обычно связывают с близостью двух наибольших по модулю собственных значений матрицы Т системы уравнений ТО = АО, когда доминантное отношение модуля второго собственного значения к наибольшему первому близко к единице. Одним из известных методов ускорения сходимости внешних итераций является метод полиномов Чебышева [41], [42]. Этот метод организует внешние итерации вида = Сп(Т)^(р) с помощью полиномов Чебышева Сп(Т) степени п, построенных на основе информации о двух наибольших по модулю собственных значениях матрицы Т. Для его использования надо предварительно оценить величину доминантного отношения. Получение оценки сравнимо по сложности с решением исходной задачи. Кроме того, метод полиномов Чебышева слишком чувствителен к вычислительной погрешности [1], [43].

В методе сдвига [40] для ускорения сходимости внешних итераций предлагается подобрать положительную константу к: > для того, чтобы доминантное отношение у новой задачи стало меньше, чем у исходной. Новая задача получается путем добавления к обеим частям системы уравнений отрицательного слагаемого, записанного па разных итерациях в виде источника деления, деленного на к. В этом случае приходится переопределять источник деления на каждой внутренней итерации, что приводит к усложнению процесса сходимости внутренних итераций в каждой группе. Кроме того, в задачах с доминантным отношением близким к единице трудно заранее угадать наиболее подходящую константу к.

В методе обратных итераций Виланда [41] предлагается выбрать константу «/с» близкую к наибольшему собственному значению и затем организовать итерационный процесс, как в методе сдвига, чтобы получить соответствующий собственный вектор. Полученную задачу нельзя свести к последовательному решению одногрупповых уравнений [41].

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

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

В методе квазидиффузии [45], [46] в основном решается многогрупповая система квазидиффузионных уравнений, а решение системы уравнений переноса используется для поправки диффузионного дробно-линейного функционала. Для расчета Кс/Г используется система кназидиффузионных уравнений, полученная суммированием исходной системы по энергетическим группам. Расчет Ке(Г происходит быстрее за счет того, что более трудоемкое решение системы уравнений переноса осуществляется редко. Метод квазидиффузии хорошо зарекомендовал себя в одномерной геометрии и в двумерной К-7, геометрии. Неизвестно, насколько этот метод эффективен в трехмерной геометрии.

На этом список известных методов ускорения сходимости внешних итераций практически заканчивается. Среди других методов ускорения выделяются методы, основанные на идее Л.А. Люстерника [47], предложенной для решения систем неоднородных линейных уравнений

Установлено [47], что при большом числе итераций р векторы х^р) - х(р+1) и х(р+,) — х(я+2) практически пропорциональны, а отношение компонент этих векторов приблизительно равно наибольшему собственному значению \ матрицы В. Если известны последовательные приближения х(/,),х(р+1),х(р+2), то можно оценить Л, по какой-либо /той компоненте: х = Вх + с в.1) методом простой итерации: х(р) = Вх{р-Х) +с.

Р+1) „(Р+2)

7 ">» в.2) а затем найти х, используя приближенное равенство: х(/,+1)-Зс(/,)). (в.З)

1-Я,

Метод Люстерника (в.2) - (в.З) был модифицирован и использован для ускорения сходимости итераций при решении неоднородной одногрупповой задачи (в.1) переноса нейтронов [48]. Величина параметра Л(р) < 1 вычислялась по формуле:

Л(р) = — j]^ -(p(p-X)\dV dV где (р(р) - приближение скалярного потока на итерации р , G - пространственная область, в которой рассматривается процесс переноса нейтронов. Линейная экстраполяция выполнялась через каждые четыре итерации по формуле:

V = (P(p)+~^{<P(p)-(P(p-X)) (в.4) с последующей нормировкой полученного результата для того, чтобы сохранить полный баланс частиц. Сделан вывод о том, что метод Люстерника может быть более эффективным, если его использовать в ходе итерационного процесса [48].

За рубежом метод Люстерника известен, как метод асимптотической экстраполяции источника ASE (Asymptotic Source Extrapolation) [40]. Этот метод применяется после выполнения нескольких простых итераций по источнику рассеяния. Оценка собственного значения Я с наибольшей амплитудой гармоники ошибки делается по формуле отношения двух норм: Я х(р+1)-х(р) в.5) а линейная экстраполяция по формуле: х = х(р+1/2) ч—^—(х{р+{12) —х(р)). 1-Я

Сделан вывод о том, что метод ASE не препятствует неприемлемо медленной сходимости внутренних итераций [40], если имеется много гармоник.

При использовании метода Люстерника для ускорения сходимости итераций возникает проблема разработки критериев его эффективного применения. В [49] описан близкий к методу Люстерника алгоритм ускорения сходимости итераций при решении неоднородной системы линейных уравнений (в.1), который называется «S2 - процесс». Этот метод рассмотрен для случая, когда матрица В имеет простую структуру. В 82 процессе компоненты приближения х(р) уточняются в ходе итерационного процесса после выполнения условия пропорциональности двух векторов [49]: х{р~1) -х(р~2\х{р) -х(р^) м(р) = не ii

1, (в.6) где в числителе стоит скалярное произведение векторов, а в знаменателе их евклидовы нормы. Если условие (в.6) выполнено, то вычисляется собственное значение Я\р) матрицы В задачи (в.1) по формуле:

П -Гп-2) -(п) ~Лп- 1К • (В-7)

Это собственное значение используется для уточнения приближения х к решению х : в.8)

В [49] получены следующие оценки ¿\р) = \ + / А, \р), х = х(р) + 0(|Я,|/') и х = у(р) +0(\А2\р), где х - точное решение (в.1). Сделан вывод о том, что у(р) является лучшим начальным условием для последующих итераций по сравнению с х(р). Производя время от времени уточнения (в.7) и (в.8), иногда удается существенно уменьшить общее число итераций. Для того чтобы делать такие уточнения, необходимо чтобы в разложении в ряд по собственным векторам п -мерной матрицы В невязки начального приближения:

1=1 одно из слагаемых преобладало над другими [49]. Условие (в.6) называется условием практической применимости линейной экстраполяции (в.7), (в.8). Кроме того, в [49] предложена возможная схема метода простой итерации с применением д2 - процесса ускорения сходимости. Отмечено наличие ряда моментов, затрудняющих реализацию метода и требующих серьезной математической подготовки и проведения большой серии численных экспериментов [49].

Ь(Р+1) -(Р+2)

В [50] предложено использовать метод Люстерника (в.З), где Я, = -ух<р) -х(/,+1) для ускорения сходимости итераций при решении задачи Ах = х, чередуя выполнение нескольких обычных итераций х(р) = Ах(р'1) с расчетом уточненного приближения, которое каждый раз принимается за новое начальное приближение.

В данной работе предложен разработанный автором новый метод ускорения сходимости итераций, который называется 5 - процесс, для решения задачи нахождения наибольшего по модулю собственного значения матрицы системы линейных уравнений. В задачах расчета Kcff этот метод служит для ускорения сходимости внешних итераций в случае их замедления и заключается в использовании линейной экстраполяции полученного приближения источника деления и соответствующих ему скалярных потоков вблизи положительных корней характеристического многочлена итерируемой матрицы, а также пересчета моментов угловых потоков. Этот метод является дальнейшим развитием метода Люстерника [47] и д2 - процесса Бахвалова [49] для решения однородных задач.

Проблема ускорения сходимости внутренних итераций возникает при решении уравнений переноса в одной энергетической группе методом итераций по столкновениям. Этот метод сходится быстро в задачах с оптически тонкой средой или с высокими скоростями реакции захвата. Однако он сходится медленно в задачах с оптически толстыми средами и в тех случаях, когда отношения с - s /Т,г для разных веществ близки к единице, где Eos - нулевой момент сечения рассеяния и Еу. - полное сечение.

Одним из способов ускорения сходимости внутренних итераций является метод ребаланса [51], [52], дополнительно требующий на каждой итерации относительно небольшого количества арифметических действий. В тех случаях, когда этот метод устойчив, он значительно ускоряет расчеты. Это нелинейный метод, основанный на дополнительном расчете на каждой итерации по столкновениям решения системы диффузионно-подобных разностных уравнений для мультипликативных пространственных поправок к нулевому и первому угловым моментам решения в Р] приближении индикатрисы рассеяния. Этот метод привлекателен еще и тем, что если решается задача в Рт приближении индикатрисы рассеяния, то для сохранения баланса частиц при ускорении все угловые моменты решения умножаются на соответствующие поправки в каждой пространственной точке. Такой возможности лишены DSA метод [53] и PXSA метод [54] ускорения сходимости внутренних итераций, т.к. в них используются аддитивные пространственные поправки. В DSA методе поправки ищутся к нулевому моменту решения, а в PXSA методе - к нулевому и первому угловым моментам решения.

Различают два способа ускорения методом ребаланса: на мелкой сетке (fine-mesh rebalance) и на грубой сетке (coarse-mesh rebalance). В методе ребаланса на мелкой сетке [55] уравнение для поправок записывается на сетке, заданной для решения уравнения переноса. В методе ребаланса на грубой сетке для расчета поправок используется более крупная сетка, объединяющая две или несколько ячеек мелкой сетки. Автором проведено теоретическое исследование устойчивости метода ребаланса в обоих случаях [56], [57] с помощью Фурье-анализа линеаризованного метода пространственного ребаланса совместно с «алмазной» схемой ^ метода дискретных ординат в одномерной плоской геометрии на равномерной сетке для случая изотропного рассеяния частиц на специальном классе задач: с постоянными сечениями и постоянным изотропным источником в бесконечной области. Результаты анализа показали, что для любого фиксированного 0 < с < 1 область устойчивости метода зависит от величины параметра Ег/г, где И - шаг сетки. Для достаточно малых значений с метод ребаланса устойчив и эффективно ускоряет в обоих случаях. При приближении значений с к 1 оба метода могут быть неустойчивыми. Скорость сходимости метода ребаланса становится медленной и метод даже расходится на сетках, образованных очень большим числом мелких ячеек [57]. Метод ребаланса на мелкой сетке устойчив в очень узкой области изменения параметра «1 для значений с « 1 [56], [40].

В [58] предложен способ повышения устойчивости метода ребаланса на мелкой сетке в одномерной геометрии путем введения фиктивных граничных токов на сторонах ячейки. Подтверждением устойчивости этого метода служат результаты численных расчетов [58]. Этот усовершенствованный метод называется методом пространственного ребаланса или методом частичных токов. Этот алгоритм обобщен на случай двумерной геометрии и предложена другая формула для фиктивных токов [59]. На основе результатов расчетов сделан вывод о том, что новая формула предпочтительнее для значений «с», близких к единице. Дальнейшее усовершенствование методики ускорения состоит из двух модификаций [60]: получены граничные условия для нахождения мультипликативных поправок, согласованные с граничными условиями уравнений переноса и способствующие ускорению и разработан алгоритм увеличения фиктивных токов с целью достижения эффективности ускорения в наиболее трудных задачах. Устойчивость окончательного варианта метода пространственного ребаланса исследована автором [61] с помощью Фурье-анализа линеаризованного метода в одномерной плоской геометрии.

Усовершенствованы имеющиеся программы в Х-У-^ геометрии (модуль К1№0) и в НЕХ^ геометрии (модуль КГЫЗБб) [26], [28] - [31] и разработаны новые программы (модуль КГЫЮ^) в трехмерной геометрии [27] для проведения трехмерных расчетов эффективного коэффициента размножения КС]Т ядерных реакторов на быстрых нейтронах и источника деления в однородных задачах, а также потоков нейтронов и гамма-квантов в задачах их защиты. В модулях КШЗБ и KINRTZ реализована взвешенная 0 - схема [62] (0-\¥ОО схема) метода дискретных ординат для расчета задач защиты. Во все трехмерные модули добавлены программы расчета индикатрисы рассеяния в Рт приближении. Для ускорения сходимости внутренних итераций реализован метод пространственного ребаланса, а для ускорения сходимости внешних итераций разработан и реализован новый метод 8 - процесс [36]. Результаты трехмерных расчетов реакторной установки СВБР 75/100 [63] (Свинцово-Висмутовый Быстрый Реактор) эквивалентной энергетической мощностью 75/100 МВт и его радиационной защиты, полученные по этим программам, доложены на конференциях [32] - [35], [37]. Данная работа построена следующим образом.

В первой главе описаны две основные задачи математического моделирования ядерных реакторов и их защиты. Раздел 1 посвящен проблеме математического моделирования ядерных реакторов на основе использования линеаризованного уравнения Больцмана для описания переноса нейтронов и гамма-квантов. В разделе 2 приведена постановка этих задач в многогрупповом приближении. В разделе 3 описано приближение индикатрисы рассеяния в виде ряда по полиномам Лежандра Рт до степени т, т.е. в Рт приближении. В разделе 4 описаны угловые квадратуры метода дискретных ординат, используемые в программах.

Во второй главе рассмотрены итерационные методы решения основных задач. Первый раздел посвящен новому методу 8 - процессу ускорения сходимости внешних итераций. В параграфе 1 описан степенной метод решения задачи расчета значений Ке]Г и источника деления, приведено свойство включения, позволяющее оценить на каждой внешней итерации верхнюю и нижнюю границы величины перечислены критерии окончания итерационного процесса. Обоснована возможность рассмотрения эквивалентной задачи, имеющей наибольшее собственное значение, равное единице, путем разложения начального приближения в ряд по собственным и присоединенным векторам итерируемой матрицы. Затем описан предлагаемый алгоритм ускорения сходимости внешних итераций и приведены его формулы. В параграфе 2 проведено исследование скорости сходимости 8 - процесса на примере решения вспомогательной задачи Ах = х поиска наибольшего собственного значения неотрицательной неразложимой квадратной матрицы А. Показано, что формулу линейной экстраполяции 8 - процесса можно интерпретировать, как шаг метода Ньютона [49], сделанный вблизи положительного корня характеристического уравнения итерируемой матрицы А. Представлена блок-схема 5 - процесса. В параграфе 3 приведены результаты тестовых расчетов Kejr критической сборки GODIVA [64] в одномерной сферической геометрии в SSP, приближении с различными значениями параметра S и без ускорения. Получена оценка эффективности ускорения. Показаны графики поведения значений параметров К(р) и Л(р) на внешних итерациях без ускорения и с использованием S - процесса ускорения для S = 0.9. В параграфе 4 рассмотрены результаты расчетов Kcjf и источника деления исходного состояния критической сборки BZD/1 в экспериментах «ZEBRA» [38] в X-Y-Z геометрии в 30-ти групповом приближении по энергии для SAPX и SsP3 приближений методом дискретных ординат. Сделан вывод о возможности выбора эффективного значения параметра д на основе анализа результатов предварительных расчетов в более низком SAP{ приближении.

Второй раздел посвящен методу пространственного ребаланса. В параграфе 1 описан алгоритм решения уравнения переноса в одной группе на одной внутренней итерации с использованием «взвешенной» схемы [1] в X-Y-Z геометрии. В параграфе 2 приведены формулы метода пространственного ребаланса и описана эффективная методика его использования. В параграфе 3 проводится исследование устойчивости метода пространственного ребаланса с «алмазной» схемой [61] с помощью Фурье-анализа линеаризованного метода в одномерной плоской геометрии на равномерной сетке. Рассмотрен случай изотропного рассеяния частиц на специальном классе задач: с постоянными сечениями и постоянным изотропным источником в бесконечной области. Устойчивость этого метода ранее не исследовалась. Результаты расчетов показали, что Фурье-анализ корректно предсказывает свойства устойчивости и сходимости нелинейного метода пространственного ребаланса [61]. В параграфе 4 приведены результаты тестовых расчетов железоводной борированной защиты [11] в одномерной геометрии 0-WDD схемой совместно с методом пространственного ребаланса. Для сравнения эффективности ускорения приведены результаты расчетов этой задачи AWDD схемой совместно с Р] SA методом, выполненные ранее в пакете «РЕАКТОР».

Третья глава посвящена реализации метода дискретных ординат и эффективных методик ускорения сходимости внутренних и внешних итераций [35] - [37], [61] в трехмерных программных модулях пакета «РЕАКТОР», предназначенных для расчета значений Keff и источника деления и задач защиты. Приведено описание каждого модуля.

В первом разделе описан модуль KIN3D расчета потока частиц в X-Y-Z геометрии в выпуклых областях, горизонтальное сечение которых задается на равномерной квадратной сетке. Во втором разделе описан модуль K1NRTZ, предназначенный для расчета потока частиц в R-(p-Z геометрии в цилиндрических областях с неравномерной сеткой. В третьем разделе описан модуль KIN3D6 расчета уравнения переноса в HEX-Z геометрии в выпуклых областях, горизонтальное сечение которых задается па сетке, состоящей из правильных шестиугольников. В четвертом разделе приведена общая блок-схема модулей и их основной подпрограммы DSNPN.

В четвертой главе приведены результаты решения задач математического моделирования проектируемого в настоящее время энергетического ядерного реактора СВБР 75/100 [63] и его защиты в трехмерной геометрии. Продемонстрирована эффективность предложенных алгоритмов и созданных программ. В первом разделе приведены результаты расчетов Ке(( и источника деления реактора СВБР 75/100 в HEX-Z геометрии в 30-ти групповом приближении для различных значений параметра 8 от 0.2 до 0.9 по DDO схеме. Эта же задача рассчитана в более высоком SsP\ приближении индикатрисы рассеяния с параметром S = 0.7, при котором результаты расчетов в SH Р] приближении были в 2 раза эффективнее по времени. По результатам расчетов разработаны рекомендации о выборе значения параметра 5. Во втором разделе рассмотрены результаты расчета радиационных полей в защите реакторной установки СВБР 75/100 (30 нейтронных групп и 19 гамма групп) в X-Y-Z и R-(p-Z геометрии с использованием 0-WDD схемы, в Ss приближении метода дискретных ординат и в Ръ — приближении индикатрисы рассеяния. Результаты расчетов сопоставлены с аналогичными результатами, полученными по программе TORT [65] при использовании одной и той же системы констант. Сделан вывод о том, что на эффективность ускорения внутренних итераций методом пространственного ребаланса влияют в основном два фактора: выбор номера итерации, с которой начинается ускорение итераций, и выбор 0 параметра схемы ©-WDD. Приведены рекомендации по оптимальному выбору этих параметров.

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

В Приложении II приведены основные формулы DDO схемы в HEX-Z геометрии.

В Приложении III описана 0-WDD схема и приведены ее основные формулы в X-Y-Z и R-cp-Z геометрии.

На защиту выносятся следующие результаты:

1. Предложен новый метод (3 - процесс) ускорения сходимости внешних итераций при решении задач расчета эффективного коэффициента размножения и источника деления, возникающих при математическом моделировании ядерных реакторов методом дискретных ординат в трехмерной геометрии.

2. Исследована устойчивость линеаризованного приближения метода пространственного ребаланса на одногрупповой задаче с изотропным рассеянием в бесконечной однородной среде с помощью Фурье-анализа. Доказана его устойчивость.

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

4. Проведены численные исследования на тестовых и реальных задачах эффективности реализованных методов ускорения сходимости внутренних и внешних итераций. Разработаны рекомендации по выбору параметра 8 для эффективного ускорения внешних итераций при расчетах в 5\Рт приближении. Разработаны рекомендации по оптимальному выбору параметра в взвешенной схемы.

5. Выполнены трехмерные расчеты проектируемого в настоящее время энергетического ядерного реактора СВБР 75/100, которые были использованы в проектных работах.

ЗАКЛЮЧЕНИЕ

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

Выполнены расчеты проектируемого в настоящее время энергетического ядерного реактора СВБР 75/100 в трехмерной геометрии по разработанным программам решения уравнения переноса методом дискретных ординат.

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

2. А.В. Воронков, Е.П. Сычугова. Характеристический метод дискретных ординат решения уравнения переноса // Препринт ИПМ АН СССР, 1987, №50, С. 1-20.

3. А.В. Воронков, Е.П. Сычугова. Комбинированные методы решения уравнения переноса//Препринт ИПМ АН СССР, 1988, №18, С. 1-31.

4. A.V. Voronkov, Е.Р. Sychugova. Combined Characteristic Methods for Solving the Transport Equation in (X,Y) Geometry // Тезисы докладов на Международном симпозиуме "Численные методы решения уравнения переноса", Москва, 26-28 мая 1992, С. 260-263.

5. Е.Р. Sychugova. Improvement of the Nodal CDSN Method of Discrete Ordinates for the Transport Equation in (X,Y) Geometry // Тезисы докладов на Международном симпозиуме "Численные методы решения уравнения переноса", Москва, 26-28 мая 1992, С. 244-247.

6. Е.П. Сычугова. Композиционный линейный характеристический CLC метод дискретных ординат для уравнения переноса в (X,Y) — геометрии // Препринт ИПМ РАН, 1992, №48, С.1-.15.

7. A.V. Voronkov, Е.Р. Sychugova. CDSN-Method for Solving the Transport Equation // Transport Theory and Statistical Physics, 1993, V. 22, Iss. 2-3, P. 221-245.

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

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

10. A.M. Волощенко, А.В. Воронков, Е.П. Сычугова. Согласованная P,SA схема ускорения внутренних и внешних итераций для уравнения переноса нейтронов и фотонов в одномерных геометриях в пакете РЕАКТОР // Препринт ИПМ РАН, 1996, №2, С. 1-29.

11. Л.В. Воронков, Е.П. Сычугова. CDSN метод решения уравнения переноса в трехмерной X-Y-Z геометрии // Препринт ИПМ РАН, 1996, №66, С. 1-15.

12. А.В. Воронков, Е.П. Сычугова. Линейный характеристический метод дискретных ординат для решения уравнения переноса в X-Y геометрии // Препринт ИПМ РАН, 1996, №91, С. 1-16.

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

14. А.В. Воронков, Е.П. Сычугова. Нестационарная методика для расчета многогрупповых двумерных кинетических уравнений на сетках, состоящих из произвольных четырехугольников // Препринт ИПМ РАН, 2000, №80, С. 1-28.

15. А.В. Воронков, Ю.Н. Орлов, Е.П. Сычугова. Уравнение переноса нейтронов в движущейся среде // Препринт ИПМ РАН, 2001, №91, С. 1-16.

16. А.В. Воронков, Е.П. Сычугова. Решение уравнений переноса нейтронов в трехмерной HEX-Z геометрии // Тезисы доклада на 15-м семинаре «НЕЙТРОНИКА 2004» Нейтронно-физические проблемы атомной энергетики, 27-30 октября 2004, г. Обнинск.

17. V.I. Arzhanov, A.V. Voronkov, A.S. Golubev, N.A. Konovalov, V.A. Krukov, E.P. Sychugova. Development of portable program system for modeling neutron-physical processes in nuclear reactor on parallel computers // Proc. of Joint Int. Conf. on

18. Mathematical Methods and Supercomputing for Nuclear Applications, 5-9 October 1997, Saratoga, P. 1454-1463.

19. A.B. Воронков, Е.П. Сычугова, B.K. Смирнов, С.JT. Головков, А.Г. Рубин. Верификация пакета REACTOR-P на взаимосогласованных тестовых задачах // Отчет ИПМ РАН, 2003, №7-16-2003, С. 1-11.

20. A.B. Воронков, A.C. Голубев, Е.П. Сычугова, и др. Решение уравнения переноса в SnPm приближении в X-Y-Z и в HEX-Z геометриях па скалярных и параллельных ЭВМ // Отчет ИПМ PAII, 2005, № 7-09-2005, С. 1-65.

21. A.B. Воронков, Е.П. Сычугова. Решение уравнения переноса в SnPm приближении в R-FI-Z геометрии // Отчет ИПМ РАН, 2005, № 7-19-2005, С. 1-40.

22. A.B. Воронков, В.И. Журавлев, Е.А. Земсков, Е.П. Сычугова, Е.В. Аверченков, П.Б. Афанасьев, A.B. Дедуль, В.В. Кальченко. Новые возможности GNDL системы константного обеспечения расчета реакторов и защиты // Тезисы доклада на IX

23. Российской научной конференции «Радиационная защита и радиационная безопасность в ядерных технологиях» 24-26 октября 2006 г., Федеральное Агентство по Атомной Энергии, ГНЦ РФ Физико-Энергетический Институт им. А.И. Лейпунского, Обнинск, С. 78-80.

24. Е.П. Сычугова. Численные методы решения уравнения переноса в трехмерной геометрии в пакете «РЕАКТОР» // Тезисы докладов на 18-м семинаре НЕЙТРОНИКА-2007, Обнинск, 30 октября 2 ноября, 2007.

25. Е.П. Сычугова. Численные методы решения уравнения переноса в многогрупповом приближении в трехмерной геометрии в пакете «РЕАКТОР» // Препринт ИПМ РАН, 2007, №78, С. 1-24.

26. Е.П. Сычугова. Методы ускорения сходимости итераций в задачах математического моделирования ядерных реакторов и защиты в пакете «РЕАКТОР» // Тезисы докладов на 19-м семинаре НЕЙТРОНИКА-2008, Обнинск, 28-31 октября, 2008.

27. A.D. Knipe. Status of UKAEA Low Power Reactors // Proc. of Int. Top. Meeting on Safety, Status and Future of Non-Commercial Reactors and Irradiation Facilities in Boise, Idaho, 30 Sept.-4 Oct. 1990.

28. A. Voronkov, V. Arzhanov. "REACTOR" program System for Neutron - Physical Calculation // Proc. Int. Top. Meting: Advances in Mathematics, Computations and Reactor Physics, 1991, V. 1, April 28 - May 2, Pittsburg, USA.

29. Marvin L. Adams, Edward W. Larsen. Fast iterative methods for discrete-ordinates particle transport calculations // Progress in Nuclear Energy, v. 40, Issue 1, 2002, p. 3-159.

30. Варга P.C. Численные методы решения многомерных многогрупповых диффузионных уравнений // В кн.: Теория ядерных реакторов. Пер. с англ. под ред. Г.А. Батя, Госатомиздат, М., 1963.

31. Г.И. Марчук, В.И. Лебедев. Численные методы в теории переноса нейтронов // Атомиздат, М., 1971.

32. Шишков Л.К. Методы решения диффузионных уравнений двумерного ядерного реактора// Атомиздат, М., 1976.

33. J. S. Warsa, Т.А. Wareing, J.E. Morel, J.M. McGhee, R.B. Lehoucq. Krylov subspace iterations for deterministic k-eigenvalue calculations // J. Nucl. Sci. Eng., 2004, V. 147, P. 26-42.

34. В.Я. Гольдин, Д.Ю. Анистратов. Реактор на быстрых нейтронах в саморегулируемом нейтронно-ядерном режиме // Математическое моделирование, 1995, Т. 7, №10, С. 1232.

35. Е.Н. Аристова, В.Я. Гольдин. Экономичный расчет многогруппового уравнения переноса нейтронов для пересчета усредненных по спектру сечений // Математическое моделирование, 2008, Т. 20, №11, С. 41-54.

36. Л.А. Люстерник. Замечания к численному решению краевых задач уравнения Лапласа и вычислению собственных значений методом сеток // Труды математического института им. Стеклова, 1947, Т. 20, С.49-64.

37. В.Н. Морозов. О решении кинетических уравнений с помощью Sn -метода // В сб. «Теория и методы расчета ядерных реакторов», Госатомиздат, М., 1962, С. 91.

38. Н.С. Бахвалов. Численные методы // Наука, М., 1973.

39. М.А. Красносельский, Г.М. Вайникко, П.П. Забрейко, Я.Б. Рутицкий, В.Я. Стеценко. Приближенное решение операторных уравнений // Наука, М., 1969, С.36.

40. E.L. Wachspress. Iterative Solution of Elliptic Systems and Applications to the Neutron Diffusion Equations of Reactor Physics // Prentice-Hall, Englewood Cliffs, New Jersey, 1966.

41. E.E. Lewis, W.F. Miller. Jr. Computational Methods of Neutron Transport // John Wiley & Sons, New York, 1984.

42. R.E. Alcouffe. Diffusion Synthetic Acceleration Methods for the Diamond-Differenced Discrete-Ordinates Equations //J. Nucl. Sci. Eng., 1977, V. 64, P. 344-355.

43. A.M. Voloschenko. Completely consistent P} synthetic acceleration scheme for chargedparticle transport calculations // Proc. 1996 ANS Topical Meeting Radiation Protection and Shielding, April 21-25, 1996, No. Falmouth, USA. V. 1, P. 408-415.

44. W. W. Engle Jr., F. R. Mynatt. A Comparison of Two Methods of Inner Iteration Convergence Acceleration in Discrete Ordinates Codes // J. Trans. Am. Nucl. Soc., 1968, V. 11, P. 193-194.

45. G. Cefiis, E.W. Larscn. Stability Analysis of Fine-Mesh Rebalance // J. Trans. Am. Nucl. Soc., 1988, V. 56, P. 309-310.

46. G. R. Cefus, E. W. Larsen. Stability Analysis of Coarse-Mesh Rebalance // J. Nucl. Sci. Eng., 1990, V. 105, P. 31-39.

47. W. H. Reed. The Effectiveness of Acceleration Techniques for Iterative Methods in Transport Theory// J. Nucl. Sci. Eng., 1971, V. 45, P. 245-254.

48. W. A. Rhoades, R. L. Childs, W. W. Engle Jr. Comparison of Rebalance Stabilization Methods for Two-Dimensional Transport Calculations // J. Trans. Am. Nucl. Soc., 1978, V. 30, P. 583.

49. W.A. Rhoades. Improvements in Discrete Ordinates Acceleration // J. Trans. Am. Nucl. Soc., 1981, V. 39, P. 753.

50. Сычугова Е.П. Исследование устойчивости и эффективности метода пространственного ребаланса для ускорения сходимости итераций в задачах переноса частиц // Математическое моделирование, 2008, Т. 20, №9, С. 75-93.

51. W. A. Rhoades, W. W. Engle. A New Weighted-Difference Formulation for Discrete Ordinates Calculations // J. Trans. Am. Nucl. Soc., 1977, V. 27, P. 776-777.

52. International Handbook of Evaluated Criticality Safety Benchmark Experiments // NEA/NSC/DOC(95)03.

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

54. B.B. Воеводин, Ю.А. Кузнецов. Матрицы и вычисления // Наука, М., 1984.

55. Lathrop K.D., Carlson B.G. Discrete Ordinates Angular Quadrature of the Neutron Transport Equation // LASL Report LA-3186 (February 1965).

56. B.G. Carlson. A Method of Characteristics and Other Improvements in Solution Methods for the Transport Equation // J. Nucl. Sci. Eng., 1976, V. 61, P. 408-425.

57. B.A. Треногин. Функциональный анализ // Наука, М., 1980, С.386.

58. А. Исимару. Распространение и рассеяние волн в случайно-неоднородных средах // Мир, М., 1981, Т. 1.

59. R. Е. MacFarlane. TRANSX-2: Code for Interfacing MATXS Cross-Section Libraries to Nuclear Transport Codes // LA-12312-MS, 1992.72. http://www.new.fT/dbforms/data/eva/evatapes/endfb6 re!8 Библиотека оцененных ядерных данных.

60. E.W. Larsen. Diffusion-Synthetic Acceleration Methods for Discrete-Ordinates Problems // J. Transport Theory and Statistical Physics, 1984, V. 13, Iss. 1-2, P. 107-126.

61. Г. H. Мантуров, M. H. Николаев, A. M. Цибуля. Программа подготовки констант CONSYST. Описание применения // Препринт ФЭИ-2828, 2000, Обнинск.