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

Размещено на http://www.allbest.ru/

2

Размещено на http://www.allbest.ru/

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

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

кандидата физико-математических наук

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

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

Сычугова Елена Павловна

Москва - 2009 год

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

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

старший научный сотрудник

Воронков Александр Васильевич

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

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

кандидат физико-математических наук,

старший научный сотрудник

Аристова Елена Николаевна.

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

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

Защита состоится «___» _______________2009 г. в _____ часов

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

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

Автореферат разослан «___» _______________2009 г.

Ученый секретарь диссертационного совета

кандидат физико-математических наук ЩЕРИЦА О.В.

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

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

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

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

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

Цель работы.

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

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

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

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

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

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

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

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

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

Апробация.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

или ,

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

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

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

. (1)

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

, (2)

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

. (3)

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

. (4)

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

, (5)

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

- единичный вектор в направлении полета частиц в трехмерной геометрии, где - полярный угол между вектором и осью Z, - азимутальный угол между его проекцией на плоскость X-Y и осью X;

- число групп нейтронов;

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

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

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

- скалярный поток частиц в точке в группе со скоростью ;

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

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

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

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

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

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

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

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

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

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

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

,

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

, (6)

где - функция плотности источника деления, , .

В конечномерном евклидовом пространстве - матрица, - собственный вектор, соответствующий собственному значению . Компонентами вектора являются значения источника деления в - той пространственной ячейке , где - ее элементарный объем. Матрица в левой части (6) удовлетворяет условиям теоремы Перрона-Фробениуса [5], т.е. она неотрицательна и неразложима. Тогда она имеет единственное простое положительное наибольшее собственное значение ( для всех ), которому соответствует собственная функция с неотрицательными компонентами, , где .

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

, (7)

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

. (8)

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

, (9)

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

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

, ,

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

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

и , (10)

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

.

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

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

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

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

На итерации с номером после того, как вычислены и приближение , полагается и . Вычисляется значение параметра по формуле:

, (11)

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

(12)

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

, (13)

и проверяется выполнение условия:

, (14)

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

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

(15)

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

. (16)

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

. (17)

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

, . (18)

Таким образом, предлагаемый алгоритм ускорения сходимости внешних итераций основан на линейной экстраполяции источника деления , групповых скалярных потоков в каждой точке пространства по формулам (11) - (17) и угловых моментов потока по формулам (18).

Алгоритм ускорения сходимости внешних итераций эффективен в случае удачно подобранных критериев. Из (13) видно, что параметр экстраполяции может быть очень большим положительным числом, поэтому значение необходимо ограничивать сверху. Обычно . Оптимальное значение параметра зависит от решаемой задачи. Из неравенства (14) следует, что значение параметра находится в интервале . В задачах физики реакторов .

Когда , т.е. вблизи точного решения, - процесс не позволяет делать экстраполяцию (15) - (18). В этом случае значения на трех последовательных итерациях не удовлетворяют какому-нибудь из неравенств (12) или (14). Поэтому итерационный процесс заканчивается без ускорения.

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

В диссертации приведены результаты тестовых расчетов критической сборки GODIVA [8] в одномерной сферической геометрии. Экспериментальное значение . В задаче заданы 30-ти групповое приближение по энергии нейтронов, угловая квадратура Гаусса-Лежандра , приближение индикатрисы рассеяния и сетка из 800 интервалов. В качестве начального приближения заданы единичные значения скалярного потока. Для окончания внутренних итераций задана величина максимальной относительной ошибки скалярных потоков и максимальное число итераций в группе 20. Критерием окончания внешних итераций служит одновременное выполнение условий и . Для ускорения сходимости внешних итераций использовался - процесс с вычислением параметра по октаэдрической норме.

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

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

Рис. 1. Изменение , и Рис. 2. Изменение , и

без ускорения. с ускорением, .

Рис.3. Изменение без ускорения. Рис.4. Изменение с ускорением,.

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

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

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

	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. Число внешних итераций в приближении при расчете задачи ZEBRA с ускорением для различных значений и без него.

	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 раза по общему числу итераций получено в приближении при (23 итерации с ускорением и 79 без него). При получено одинаковое ускорение в 2.7 раза в обоих приближениях. Это означает, что параметр можно выбрать, проведя расчеты в более низком угловом приближении и в приближении индикатрисы рассеяния . Поэтому для эффективного расчета однородной задачи на в более высоком приближении нужно сначала выбрать оптимальное значение параметра путем проведения расчетов в более низком приближении, например , и затем использовать его для более точных расчетов.

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

В диссертации описан алгоритм решения уравнения переноса в одной группе методом итерации по столкновениям [2] с использованием «взвешенной» схемы [2] в трехмерной X-Y-Z геометрии. Соотношение баланса частиц в ячейке х х для направления получается интегрированием уравнения переноса (1) с правой частью по объему ячейки и по направлениям из углового диапазона вокруг направления :

. (19)

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

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

, (20)

где

, , ,

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

. (21)

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

(22)

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

. (23)

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

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

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

, (24)

, (25)

, (26)

, , (27)

, , (28)

, , (29)

, (30)

, (31)

где - направляющие косинусы и - веса квадратуры, . В правой части (24) слагаемое с внутригрупповым рассеянием определяется по приближению скалярного потока на предыдущей итерации . В эту систему входит уравнение (25) «алмазной» схемы. По квадратурной формуле (26) вычисляется промежуточный скалярный поток , который в случае отсутствия ускорения подставляется в правую часть (24) для расчета на следующей итерации. При расчете с ускорением вычисляются односторонние токи по формулам (27), которые увеличиваются по формулам (28) и (29) с множителем на каждой внутренней итерации. Увеличенные токи являются коэффициентами системы конечно-разностных уравнений (30) относительно поправок . Итерация заканчивается перенормировкой (26) скалярных потоков.

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

, (32)

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

. (33)

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

, (34)

, (35)

, (36)

, (37)

, (38)

где введено значение .

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

, ,

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

, , , , . (39)

Подставляя (39) в (34) - (38) и решая полученную систему уравнений, найдем , где амплитуда определяется по формуле:

. (40)

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

.

В таблице 3 показаны результаты исследования поведения модуля амплитуд (40) при и для различных значений от 1 до 6 в области максимальных значений шага . При вычислении значения и соответствовали квадратуре Гаусса - Лежандра. В этой таблице приведены: общее число гармоник , шаг , значения первых пяти амплитуд и номер амплитуды , на которой достигается максимальное значение.

Таблица 3. Значения модуля амплитуд для и .

1	3.1415						1
2	1.5708	0.17951					2
3	1.0472	0.25079	0.43796				3
4	0.78540	0.45181	0.31867	0.52379			4
5	0.62832	0.68696	0.35530	0.35657	0.52039		1
6	0.52360	0.93611	0.43889	0.32583	0.38069	0.47759	1

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

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

Рис. 5. Зависимость спектрального радиуса от параметра для , .

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

В диссертации приведены результаты тестовых расчетов железоводной борированной защиты в одномерной геометрии с использованием И-WDD схемы [15] совместно с методом пространственного ребаланса. Для сравнения эффективности ускорения приведены результаты расчетов этой задачи AWDD схемой совместно с методом ускорения, выполненные ранее в пакете «РЕАКТОР» [16]. Результаты расчетов показали, что метод пространственного ребаланса эффективнее метода и обеспечивает существенный вычислительный выигрыш примерно в 5 раз по времени, в то время как метод дает меньший выигрыш по времени в 3 - 4 раза.

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

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

В первой части главы приведены результаты расчетов и источника деления реактора СВБР 75/100 с использованием схемы с «нулевой» коррекцией в HEX-Z геометрии в 30-ти групповом приближении по энергии нейтронов для различных значений параметра от 0.2 до 0.9. Пространственная сетка состояла из 54 плоскостей по вертикальной оси Z и из 22789 точек (гексагональных ячеек) в каждой плоскости ОXY. Более высокое приближение индикатрисы рассеяния рассчитано с параметром , при котором результаты расчетов в приближении были в 2 раза эффективнее по времени счета, чем расчет без ускорения. По результатам расчетов разработаны рекомендации о выборе значения параметра .

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

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

Аналогичные результаты были получены в R-ц-Z геометрии. Рассчитано - приближение. Использована сетка из 110 интервалов по радиусу, 120 интервалов по углу ц и 137 интервалов по оси Z, что составляет 1808400 ячеек.

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

ЗАКЛЮЧЕНИЕ

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

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

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

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

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

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

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

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

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

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

2. A.V. Voronkov, E.P. Sychugova. A second-order finite volume 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 1, p. 3-159.

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

3. Л.А. Люстерник. Замечания к численному решению краевых задач уравнения Лапласа и вычислению собственных значений методом сеток // Труды математического института им. Стеклова, 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. Н.С. Бахвалов. Численные методы // Изд-во «Наука», М., 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. А. 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. Klimov, 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 // 5^th International Conference on Nuclear Option in Countries with Small and Medium Electricity Grids, Dubrovnik, Croatia, 2004, May 16-20.

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.

Размещено на Allbest.ru

...