авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ

Pages:   || 2 | 3 | 4 |
-- [ Страница 1 ] --

НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ ЦЕНТР

КУРЧАТОВСКИЙ ИНСТИТУТ

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

УДК 621.039.514

Кондрушин Антон Евгеньевич

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

ЗАДАЧ НЕЙТРОННОЙ ПРОСТРАНСТВЕННОЙ КИНЕТИКИ

В ЯДЕРНЫХ РЕАКТОРАХ

Специальность: 05.13.18

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

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

Научный руководитель:

д.т.н. Бояринов В.Ф.

Москва – 2014 Оглавление Введение…………………………………………………………………………... ….. Глава 1 Обзор литературы……………………………………………………....... 1.1 Общие уравнения………………………………………………………......... 1.2 Обзор проекционно-сеточных методов решения уравнения переноса нейтронов…………………………………………………………………...... 1.2.1 Метод конечных элементов…………………………………….......... 1.2.2 Метод граничных элементов………………………………………...... 1.2.3 Метод конечных суперэлементов Федоренко……………………..... 1.2.4 Метод матриц отклика……………………………………………....... 1.3 Метод поверхностных гармоник………………………………………….... 1.3.1 Основы метода поверхностных гармоник и определения………..... 1.3.2 Обзор программных комплексов SUHAM и SVS………………....… 1.3.3 Обзор работ по нестационарным уравнениям МПГ……………....… 1.4 Обзор методов решения нестационарного уравнения переноса нейтронов………………………………………………………………….. …. 1.4.1 Полностью явный метод…………………………………………… …. 1.4.2 Полностью неявный метод………………………………………….… 1.4.3 -метод……………………………………………………………….… 1.4.4 Метод переменных направлений………………………………….. …. 1.4.5 Улучшенный квазистатический метод………………………….….… 1.4.6 SCM метод………………………………………………………….. …. 1.5 Обзор программ для решения нестационарного уравнения переноса нейтронов………………………………………………………………….. …. Заключение к главе 1…………………………………………………………. …. Глава 2 Нестационарные уравнения метода поверхностных гармоник…. …. 2.1 Двумерные нестационарные уравнения МПГ…………………………...… 2.1.1 Поверхностная невязка……………………………………………...… 2.1.2 Объемная невязка…………………………………………………....… 2.1.3 Вывод уравнений МПГ………..…………………………………….… 2.2 Одномерные уравнения МПГ……………………………………………..… 2.3 Итерационная схема………………………………………………………..… Заключение к главе 2………………………………………………………….





.… Глава 3 Программный комплекс SUHAM-TD и его верификация………...… 3.1 Программный комплекс SUHAM-TD…………………………………… …. 3.2 Верификация программного комплекса SUHAM-TD………………….. …. 3.2.1 Тест BSS-6………………………………………………………….. …. 3.2.2 Тест PHWR…………………………………………………………. …. 3.2.3 Тест TWIGL………………………………………………………… …. 3.2.4 Модифицированный тест 8-A1……………………………………. …. 3.2.5 Транспортный тест TWIGL………………………………………...... Заключение к главе 3………………………………………………………….... Глава 4 Разработка и расчет пространственно-временного бенчмарка C5G7-TD для тестирования кинетических нейтронно-физических кодов.... 4.1 Обзор пространственно-временных бенчмарков…………......................... 4.2 Описание бенчмарка C5G7………………………………………………..... 4.3 Расчет кинетических характеристик для теста C5G7-TD………………... 4.4 Законы ввода реактивности для теста C5G7-TD………………………...... 4.5 Результаты моделирования теста C5G7-TD……………………………..... Заключение к главе 4………………………………………………………….... Заключение………………………………………………………………………..... Обозначения………………………………………………………………………... Список литературы………………………………………………………………... Приложение А Копии свидетельств о регистрации модулей SUHAM-TD.... Приложение Б Результаты расчета теста BSS-6………………………….….... Приложение В Результаты расчета теста PHWR……………………………... Приложение Г Результаты расчета теста TWIGL…………………………..... Приложение Д Результаты расчета теста C5G7-TD………………………...... Введение Развитие современной атомной энергетики требует повышенного внимания к характеристикам надежности и безопасности ядерных реакторов. Это внимание обуславливается наличием потенциальной возможности возникновения аварии и значительными затратами при строительстве атомных станций, по причине проектирования с запасом с целью предотвращения потенциально возможных аварий. Важнейшую роль в проектировании надежных, безопасных и вместе с тем экономически выгодных ядерных реакторов играет проведение исследовательских и проектных расчетов, позволяющих приблизится к оптимальному соотношению этих показателей. При этом важнейшее значение имеет проведение качественного нейтронно-физического расчета.

В последние годы все большее внимание уделяется развитию кодов, позволяющих проводить нестационарный расчет ядерного реактора. Это связано с наличием факта, что большинство существующих на данный момент нестационарных кодов содержит ряд серьезных приближений в своей нейтронно физической составляющей на фоне высоких требований к безопасности реакторов. К таким приближениям в первую очередь относятся пространственная гомогенизация, расчет в малом числе энергетических групп, диффузионное приближение, а также применение разного рода поправок, полученных расчетно экспериментальным путем и призванных уменьшить ошибки в моделировании по заложенной в расчет математической модели (например [1, 2]). Вместе с тем, следует также отметить, что ряд исследователей (например, [3–5]) указывает на возможность получения улучшенных результатов путем решения более универсального газокинетического уравнения, особенно для анализа безопасности. Результаты, полученные таким путем, позволили бы обрести большую уверенность в качестве получаемых результатов. Таким образом, можно заключить, что существует необходимость в создании нестационарных нейтронно-физических программ, в алгоритмы которых не заложены вышеуказанные приближения.





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

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

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

В качестве метода, который мог бы лежать в основе нейтронно-физической составляющей современного нестационарного кода, обладающего вышеизложенными характеристиками, предлагается использовать предложенный Н.И. Лалетиным [6–10] метод поверхностных гармоник (МПГ). МПГ является процедурой построения математической модели для описания нейтронно физических процессов в ядерном реакторе, учитывающей особенности нейтронно-физических реакторных задач [9].

Главной причиной для обращения к МПГ стал тот факт, что данный метод, реализованный в стационарных программных комплексах SUHAM [11, 12] и SVS [13], позволяет получать основные нейтронно-физические функционалы с точностью сравнимой с точностью детерминистических методов при вычислительных затратах сравнимых с затратами инженерных подходов. Данные свойства МПГ проявил как при расчете классических сильно идеализированных бенчмарков, так и при расчете моделей как проектируемых (БРЕСТ-ОД-300, БН-1200 [14], ГТ-МГР [15]), так и реальных реакторов (РБМК-1000 [14], ВВЭР-1000 [16, 17]). Данный метод реализован в одно-, двух- и трехмерной геометриях [14, 17], в системах с квадратной и треугольной решеткой. Все это позволяет рассматривать МПГ как хорошо апробированный и состоявшийся метод, проявивший свои положительные черты на широком классе реакторных задач. Кроме этого, МПГ очень удобен для применения алгоритмов параллельных вычислений [18].

Для применения МПГ в нестационарных задачах необходимо провести детальный математический вывод нестационарных уравнений МПГ и их верификацию посредством программной реализации и расчета ряда бенчмарков.

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

Здесь следует заметить, что верификация такого кода на этапе разработки осложняется следующей проблемой. Существующий на сегодняшний день набор пространственно-временных бенчмарков содержит своего рода пробел. С одной стороны существует ряд классических диффузионных бенчмарков (например, TWIGL [19], 8-A1 [20]) представляющих собой несколько гомогенных зон и содержащих диффузионные нейтронно-физические константы в виде малогрупповых макроконстант и прочих величин. Такие тесты не позволяют провести верификацию программы, проводящей расчеты без пространственной гомогенизации и без диффузионного приближения. С другой стороны, имеется набор бенчмарков, которые описывают гетерогенную структуру среды, содержат концентрации нуклидов в качестве характеристик материалов, включают в себя характеристики обратных связей и т.д. Например, к таким задачам можно отнести PWR MOX/UO2 core transient benchmark [21], PBMR coupled neutronics/thermal hydraulics transient benchmark the PBMR-400 core design [22]. Результаты их расчетов содержат дополнительные погрешности (погрешность ядерных данных, погрешности тепло-гидравлических кодов и другие) что не позволяет выделить погрешность метода, заложенного в основу нейтронного кода.

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

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

1. проведение детального вывода пространственно-временных уравнений метода поверхностных гармоник посредством классического метода минимизации невязки на основе нестационарного уравнения переноса нейтронов;

2. построение численного алгоритма решения нестационарного уравнения переноса нейтронов в одномерном и двумерном ядерном реакторе на основе полученных уравнений МПГ;

3. программная реализация разработанных алгоритмов в рамках программного комплекса SUHAM-TD [23];

4. верификация созданного кода и проведение с его помощью исследований применения метода поверхностных гармоник в пространственно-временном расчете;

5. разработка нестационарного бенчмарка C5G7-TD для решения уравнения переноса на основе стационарного бенчмарка C5G7 (benchmark on deterministic transport calculations without spatial homogenization) [24] посредством получения нестационарных характеристик материалов.

Проведение расчета предложенного бенчмарка средствами программного комплекса SUHAM-TD.

Научная новизна результатов, представленных в работе, состоит:

• в разработке алгоритмов и их реализации в расчетной программе SUHAM-TD для решения одномерных и двумерных конечно-разностных нестационарных уравнений МПГ в ядерных реакторах с квадратной решеткой;

• в проведении верификации разработанного кода SUHAM-TD на ряде бенчмарков с демонстрацией эффективности уравнений и разработанного кода;

• в создании пространственно-временного бенчмарка C5G7-TD для решения уравнения переноса на основе международного стационарного бенчмарка C5G7 с приведением результатов расчета.

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

Практическая ценность полученных результатов определяется тем, что:

• полученные нестационарные уравнения МПГ могут быть применены для любого типа реакторов с квадратной регулярной решеткой блоков;

• разработанные алгоритмы и модули программного комплекса SUHAM-TD предлагаются в качестве базового ядра для создания современного вычислительного инструмента для анализа нейтронных переходных процессов;

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

Личный вклад автора. Диссертант является соавтором научных работ по теме исследования и все основные результаты диссертации получены при непосредственном участии автора, а именно:

• нестационарные уравнения МПГ и расчетные алгоритмы, реализованные в созданных кодах;

• разработанный программный комплекс SUHAM-TD за исключением программ РАЦИЯ [25], DICPN [26];

• верификационные исследования разработанного алгоритма и программ;

• созданный нестационарный тест C5G7-TD и его расчет по программе SUHAM-TD.

На защиту выносятся • алгоритмы решения нестационарных уравнений МПГ с тремя пробными матрицами, реализованные в созданных кодах, для случая реактора с квадратной решеткой блоков;

• разработанные программы комплекса SUHAM-TD;

• результаты верификации разработанного программного комплекса;

• созданный пространственно-временной бенчмарк C5G7-TD.

Апробация работы. Основные положения диссертации докладывались на следующих российских и международных конференциях, школах и семинарах:

• межведомственный ежегодный семинар по нейтронно-физическим проблемам атомной энергетики «Нейтроника-2010» (г. Обнинск, 26 – 28 октября 2010 г.);

• 8-ая Курчатовская молодежная научная школа (г. Москва, 22 – 25 ноября 2010 г.);

• международная конференция по физики ядерных реакторов «PHYSOR»

(г. Ноксвилл, Теннесси, США, 15 – 20 апреля 2012 г.);

• научная сессия НИЯУ МИФИ (г. Москва, 1–6 февраля 2013 г.);

• международная конференция по математическому моделированию и расчету ядерных реакторов «M&C» (г. Сан-Валли, Айдахо, США, 5 – 9 мая 2013 г.);

• межведомственный ежегодный семинар по нейтронно-физическим проблемам атомной энергетики «Нейтроника-2013» (г. Обнинск, 6 – 8 ноября 2013 г.);

• 11-ая Курчатовская молодежная научная школа (г. Москва, 12 – 15 ноября 2013 г.).

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

Физика ядерных реакторов».

Структура и объем диссертации. Диссертация состоит из введения, четырех глав, заключения, списка литературы из 122 наименований и пяти приложений, содержит 171 страницу, 66 таблиц и 34 рисунка.

Первая глава содержит краткий обзор проекционно-сеточных методов, нашедших применение для решения нейтронно-физических задач, к которым относится и метод поверхностных гармоник. Дано краткое описание МПГ.

Описаны программные комплексы SUHAM и SVS с рассмотрением круга основных решенных ими задач. Приведен краткий обзор методов для решения нестационарного уравнения переноса нейтронов и работ по нестационарным уравнениям МПГ. Дан краткий обзор основных современных недиффузионных нестационарных нейтронно-физических кодов для расчета реакторной кинетики.

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

Третья глава посвящена описанию созданного программного комплекса SUHAM-TD и его верификации на классических пространственно-временных бенчмарках. Проводится исследование применения МПГ в нестационарном случае и демонстрация его эффективности.

В четвертой главе приводится описание созданного пространственно временного бенчмарка C5G7-TD. Приведено описание алгоритма расчета кинетических параметров. Представлены результаты расчета предложенного бенчмарка посредством программного комплекса SUHAM-TD.

Глава 1 Обзор литературы 1.1 Общие уравнения Данная глава посвящена обзору методов и программ для решения одной из основных задач анализа нейтронных переходных процессов в ядерном реакторе, а именно решению нестационарного уравнения переноса нейтронов с запаздывающими нейтронами [27]:

(r,, E, t ) + (r,, E, t ) + t (r, E, t )(r,, E, t ) = v( E ) t = dE d s (r, E E,, t )(r,, E, t ) + qextern (r,, E, t ) + (1.1) p (E) J dE d (1 (r, E ) ) f (r, E, t )(r,, E, t ) + + ( E ) j C j (r, t ) d,j 4 j = C (r, t ) = j C j (r, t ) + dE d (r, E ) f (r, E, t ) (r,, E, t ). (1.2) t j где v ( E ) – скорость нейтронов с энергией E, (r,, E, t ) – функция распределения нейтронов, t (r, E, t ) – полное макроскопическое сечение взаимодействия, s (r, E E,, t ) – макроскопическое сечение рассеяния, – qextern (r,, E, t ) внешний источник нейтронов, p ( E ) – спектр мгновенных нейтронов деления, (r, E ) – доля запаздывающих нейтронов, f (r, E, t ) – макроскопическое сечение деления нейтронов, – число вторичных нейтронов на акт деления, d, j ( E ) – спектр запаздывающих нейтронов в j-ой группе запаздывающих нейтронов, j – постоянная распада предшественников запаздывающих нейтронов группы j, C j (r, t ) – концентрация предшественников запаздывающих нейтронов группы j.

Запишем данные уравнения с использованием группового энергетического приближения и транспортного приближения при учете анизотропии рассеяния:

(r,, t ) = K(r,, t ) + 1 [ C (r, t ) ] J v inv (r,, t ) + L (1.3) t jj d, j 4 j = G C j (r, t ) = j C j (r, t ) + j,g (r ) d f, g (r, t ) g (r,, t ). (1.4) t g =1 где v inv – диагональная матрица, элементами которой являются обратные скорости, операторы L и K имеют следующий вид:

L(r,, t ) = (r,, t ) + tr (r, t )(r,, t ) s,tr (r, t )(r, t ), p J G (1 j,g (r)) d f,g (r, t )g (r,, t ), K(r,, t ) = 4 j =1 g =1 (r, t ) = (r,, t )d.

Для обеспечения единственности решения системы уравнений (1.3) и (1.4) дополнительно ставится граничное условие на границе расчетной среды VS (r,, t ) = s (r,, t ), r VS, (1.5) и начальное условие (r,,0) = 0 (r, ). (1.6) Далее в работе будет рассматриваться система уравнений (1.3–1.6) описывающая пространственно-временную кинетику в ядерном реакторе.

Как уже было отмечено, данная глава посвящена анализу методов и программ, используемых для решения задачи (1.3–1.6). Так как решению нестационарной задачи предшествует решение стационарной задачи, то проводится обзор методов решения стационарной формы уравнения переноса нейтронов (1.3). При этом внимание уделяется только методам, принадлежащих, как и МПГ, к классу проекционно-сеточных методов [28, 29], что позволяет рассмотреть метод поверхностных гармоник в совокупности с методами из этого класса. Приводится описание основ МПГ и его текущей реализации. Приводится краткий анализ работ по нестационарным уравнениям МПГ. Также рассматриваются методы решения нестационарного уравнения переноса и наиболее известные современные нейтронно-физические программы, проводящие пространственно-временные расчеты переходных процессов в ядерном реакторе.

1.2 Обзор проекционно-сеточных методов решения уравнения переноса нейтронов Одна из основных идей МПГ, а именно поиск решения краевой задачи в виде линейной комбинации пробных функций и некоторых коэффициентов, является широко распространенным подходом и реализуется в классе методов, называемых проекционными (projection methods) [30, 31]. Для описания сути этих методов будим их рассматривать на примере односкоростного диффузионного уравнения D 2 (r ) a (r ) + q (r ) = 0 (1.7) с граничным условием (r ) = s (r ), r VS, (1.8) где D – коэффициент диффузии нейтронов, Vs – граница расчетной области V.

Идея проекционных методов заключается в том, что точное решение (r ) краевой задачи (1.7), (1.8) ищется в виде функции (r) удовлетворяющей граничным условиям (1.8) и являющейся линейной комбинацией пробных (базисных) функций:

I (r ) (r ) = ci i (r ), (1.9) i = где ci – неизвестные коэффициенты, i (r ) (i = 1,2,..., I ) – пробные функции или базисные функции Пробные функции должны быть линейно (базис).

независимыми [32]. Данный подход поиска решения краевых задач является альтернативным подходом [32] по отношению к конечно-разностному методу (МКР) [33]. Необходимо отметить, что характер получаемого решения в определенной степени определяется характером пробных функций.

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

Проекционные методы отличаются по способу поиска коэффициентов cn в (1.9). Их подразделяют на вариационные методы (например, метод Релея-Ритца (Rayleigh-Ritz method) [32]) и методы взвешенных невязок (method of weighted residuals, MWR). К последним относят, например, метод коллокаций (collocation method) [32, 34] и метод Галеркина (Galerkin weighted residual method) [32, 34].

Проекционно-сеточные методы, о которых далее пойдет речь, являются своего рода синтезом проекционных и разностных методов. Их преимущество заключается в том, что к плюсам проекционных методов добавляется то, что эти алгоритмы приводят к «системам уравнений, подобным возникающим в разностных методах (т.е. незначительное число элементов матриц этих систем были бы ненулевыми)» [28, с.17]. Для получения данных методов, достаточно в проекционных методах в качестве базисных функций i (r ) (i = 1,2,..., I ) брать финитные функции, т.е. функции, которые отличны от нуля лишь на небольшой части расчетной области решаемой задачи, что приводит к системе некоторых разностных уравнений, близкой к системе, возникающей в разностном методе [28].

Интерес к рассмотрению этих методов в контексте данной работы вызван тем, что вышеописанное свойство проекционно-сеточных методов, является одним из базовых моментов МПГ. Здесь кратко рассмотрим проекционно сеточные методы, более полно нашедшие применение в нейтронно-физических расчетах ядерных реакторов, а именно метод конечных элементов [32], метод граничных элементов [35], метод конечных суперэлементов Федоренко [36] и метод матриц отклика [37].

1.2.1 Метод конечных элементов В методах Релея-Ритца и Галеркина пробные функции i (r ) обычно являются полиномами и чтобы повысить точность расчета необходимо увеличивать степень этих полиномов, что в свою очередь ведет к повышению сложности алгоритма. Как альтернативу этому можно рассматривать возможность разделения всей расчетной области на подобласти (конечные элементы), что позволяет повысить точность расчета при сохранении порядка пробных функций i (r ), благодаря учету свойств элементов по отдельности. Т.е. «сплошная среда разбивается на ряд элементов, которые можно рассматривать как конкретные ее части» [35, с.9]. Это является фундаментальной идеей метода конечных элементов (МКЭ, finite element method, FEM) [32].

Рассмотрим основные моменты данного метода. В МКЭ расчетная область V разбивается на элементы Vn. В результате глобальное решение (r) представляет собой сумму решений по всем элементам:

(r ) = n (r ). (1.10) n Внутри n-го элемента функция n (r) ищется в виде следующей линейной комбинации I n (r ) = c n,i n,i (r ). (1.11) i = В литературе функции n,i (r ) часто называют функциями формы (shape functions) [32]. Функции формы равны нулю за пределами рассматриваемого элемента и должны быть линейно независимыми. Точки пересечения границ элементов называют узлами.

Отметим, что МКЭ «может основываться как на вариационных принципах, так и на более общих выражениях метода взвешенных невязок» [35, с.9]. Таким образом, для получения уравнений МКЭ используется как метод Релея-Ритца, так и метод Галеркина, основанные на идее наличия некого функционала J который зависит от коэффициентов cn,i [32]. В случае МКЭ вместо всей расчетной области c функционалом J(ci) рассматриваются элементы с функционалами Jn(сn,i) которые в сумме равны J. Первым способом для получения уравнений МКЭ является метод Релея-Ритца при котором, для каждого m-го узла получают узловое уравнение (nodal equation), которое имеет вид:

J J = n = 0, (1.12) cn,i n cn,i где n – номер элемента, содержащего m-ый узел. Все узловые уравнения собираются в одно матричное уравнение, которое называется уравнением системы (system equation).

Другой способ для получения уравнений МКЭ – метод Галеркина. В этом случае общий интеграл системы J = w j (r) R (r )dr, (1.13) V где W j (r ) – j-ая весовая функция, R(r ) – невязка, получаемая при подстановке в решаемое уравнение приближенного решения (1.10). В качестве весовых функций w j (x) выбираются функции формы n,i (r ) [32] и получают систему уравнений для n-го элемента (element equations):

(r )R(r )dr = 0, i=1,..,I. (1.14) n,i Vn Уравнения (1.14) называются уравнениями элементов (element equations). В качестве альтернативы также возможно получение узловых уравнений. Можно отметить, что уравнение узла может быть получено комбинацией соответствующих уравнений элементов. Необходимо отметить, что в двух- и трехмерных случаях использование уравнений элемента существенно проще, чем узловых уравнений [32].

Важно отметить, что МКЭ позволяет работать с областями произвольной формы. Особенностью МКЭ является то, что в нем можно использовать неравномерные сетки. Такие сетки могут быть статичными, т.е. неизменными в течение всего расчета, а могут быть динамичными, т.е. изменяться по мере расчета, отслеживая особенности результатов. Это делает МКЭ очень гибкой технологией для систем с произвольной геометрией.

Если рассматривать применение МКЭ для уравнения (1.7) с применением метода Галеркина, то отправной точкой является ослабленная формулировка (weak form) [38] уравнения (1.7) которую можно получить, применяя первую формулу Грина [39]:

v uvdV = u n dS (u;

v )dV, (1.15) V S V где V – некоторый объем, ограниченный достаточно гладкой поверхностью S, u = u (r ) и – функции, непрерывные вместе со своими первыми v = v (r ) производными в внутри V + S и имеющие непрерывные вторые производные внутри S, n – внешняя нормаль к поверхности S. Таким образом, имеем следующую ослабленную формулировку уравнения (1.7):

(w j (r );

(r )) + a w j (r )(r ) 1 w j (r)q(r ) dV w j (r) (r ) dS = 0.

(1.16) n V D D S Следует также отметить, что, как правило, в роли функций формы выступают полиномы [40].

Поясним необходимость применения ослабленной формулировки.

Применение метода взвешенных невязок Галеркина накладывает различные требования на непрерывность функции (r ) и на ее производные. Применение данного метода «напрямую» к уравнению (1.7) требует, чтобы функция (r ) была интегрируема с квадратом второй производной, т.е.

d d 2 + dr + dr 2 dV.

(1.17) V При этом весовой функции достаточно быть лишь интегрируемой с квадратом, т.е.

dV.

w (1.18) j V Отметим, что во многих случаях предпочтительнее снижать степень непрерывности искомой функции, что в данном случае делается посредством применения первой формулы Грина. В результате имеем, что функции w j (r ) и (r ) должны быть непрерывными вместе со своими первыми производными.

Одинаковое требование непрерывности к функциям w j (r ), (r ) позволяет использовать одни и те же базисные функции при их разложении, что приводит к симметричным матрицам [35].

Первое применение метода конечных элементов к теории диффузии нейтронов, как отмечается в [38], датируется 1970-и годами. Эти разработки в области применения МКЭ к диффузии нейтронов и в частности к уравнениям переноса были обобщены и описаны в [41–43]. МКЭ реализован в ряде программ нейтронно-физического расчета и здесь отметим лишь некоторые из них:

NEDPCM [38], FTRAN, FTRAN2 [44, 45], EVENT [46].

1.2.2 Метод граничных элементов В отличие от МКЭ, который использует дискретное представление, как самой области, так и ее границы, метод граничных элементов (МГЭ, boundary element method, BEM) основывается на дискретном представлении лишь внешней границы [35]. Если в случае МКЭ используется ослабленная формулировка, снижающая порядок производной на единицу по отношению к исходной задаче, то в МГЭ используется обратная формулировка, которая снижает порядок производной искомой функции еще на единицу, но при этом имеется повышение порядка производной весовой функции. Для получения обратной формулировки уравнения (1.7) применим вторую формулу Грина [38] v u (uv vu )dV = u n v n dS, (1.19) V S где V, S, u = u (r ), v = v (r ), n – те же что и для (1.15). Зная, что ток нейтронов на поверхности равен (r ) j (r ) = D, (1.20) n получим w(r )[D (r ) (r ) + q(r )]dV = 0, a V w(r ) [(r ) w(r ) w(r ) (r ) + w(r )q(r )]dV + D j (r) (r ) n dS = 0.

w(r ) (1.21) a V S В качестве весовой функции применяются функция Грина G (r, ) исходного уравнения D 2 G (r, ) a G (r, ) + (r ) = 0, (1.22) где (r ) – дельта-функция Дирака.

В МГЭ проводится дискретизация границы S посредством множества узлов, разделяющих ее на отдельные отрезки (граничные элементы). Заменяя весовую функцию w(r ) на фундаментальное решение для узла n-го граничного элемента G (r, n ) в (1.21) и используя (1.22) и свойства функции Дирака, получаем [47, 48] q(r ) G (r, n ) G (r, n ) с( n )( n ) + G (r, n ) j (r ) (r ) dV + dS = 0, (1.23) n S D D V где 0, r V + S, c(r ) = 0,5, r S 1, r V, r S Таким образом, использование функций Грина в качестве весовых функций позволяет обеспечить обращение в нуль интегралов по области и сводит задачу только к поиску неизвестных на границе области [35]. Так, в интегральном уравнении (1.23) содержатся неизвестные только на границе, т.е. (1.23) является граничным интегральным уравнением (boundary integral equation, BIE).

Поскольку в результате BIE включает неизвестные только на границе системы, размерность задачи практически сокращается на единицу. Поток на границе S представляется как (r ) = n (r )I n, (1.24) n где n (r ) – пробная функция n-го граничного элемента, I n – значение потока в узле n-го граничного элемента. Аналогичным образом представляется функция тока j (r ) и после ее подстановки вместе с (1.24) в BIE получают систему относительно значений потоков и токов в узлах [35]. Основным преимуществом МГЭ по сравнению с МКЭ является снижение размерности задачи на единицу.

Следствием этого является сокращение объема расчетов. Первое применение МГЭ к уравнению диффузии нейтронов, как отмечается в [35], датируется серединой 1980-х [49]. В качестве программы, реализующей МГЭ в нейтронно физических расчетах, можно отметить NEDPCM [35].

1.2.3 Метод конечных суперэлементов Федоренко Метода конечных суперэлементов Федоренко (МКСЭ) был предложен [36, 50] для решения задач теории диффузии, ядерных реакторов, кинетики и т.д.

Метод нацелен для работы с областями, имеющими резкие неоднородности в малых зонах по сравнению со всей рассматриваемой областью. Такие проблемы, как правило, решаются с использованием адаптированных сеток, либо с использованием очень мелких сеток. При этом в первом случае требуются специально разработанные алгоритмы, а во втором случае – большие вычислительные ресурсы. Однако, беря за базовый принцип идею, что резкие неоднородности, как правило, расположены в малых подобластях, становится возможным использование крупных сеток, грани и узлы которых проходят через относительно гладкие части решения [51].

Другой особенностью метода является то, что в методе Федоренко функции n,i (r ) являются решениями исходного уравнения в элементе с определенными граничными условиями на границе этого элемента, в отличие от МКЭ, где базис строится из общих, не связанных с решаемой задачей, соображениях, и функции n,i (r ) записываются в явном виде. Из условий сшивки на координатных линиях сетки, состоящих из требований непрерывности искомой функции и непрерывности нормальных к линии сетки производных, строится система «балансных соотношений». Следует заметить, что функция (r ) непрерывна и почти всюду, за исключением координатной сетки, удовлетворяет решаемому уравнению.

Таким образом, в МКСЭ базис состоит из решений исходной задачи и это определяет его преимущество. МКСЭ претендует на расчет с большим шагом, что вместе с тем определяет его ограниченность. Также следует заметить, что МКСЭ требует предварительного расчета базисов [52].

1.2.4 Метод матриц отклика Первое рассмотрение метода матриц отклика (response matrix method) как отмечается в [53], относится еще к 1963 г. [37]. Метод имеет сходные с методом Федоренко подходы. Рассмотрим основные моменты данного метода на примере уравнения переноса нейтронов в соответствии с работой [53].

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

Набор граничных условий представляет собой втекающие односторонние токи нейтронов равные некоторому ортогональному базису m, m = 0,1... в фазовом пространстве, определенному для границы n-го элемента. Таким образом, решается уравнение переноса нейтронов в элементе с граничным условием J ( is ) = m ( ns ), (1.25) где J ( ns ) – ток падающих нейтронов на поверхность S n-го элемента, – вектор, представляющий фазовое пространство. Результатом локального этапа является набор функций m ( ns ) ( m ( ns ) – функция распределения нейтронов в n-ом ns ns элементе при граничном условии m заданном на поверхности элемента S).

Граничное условие (1.25) накладывается на все элементы среды независимо от того являются они граничными или нет.

Рассмотрим уравнение метода матриц отклика, решаемое на глобальном уровне расчета. Пусть J nt + ( nt ) – ток вылетающих частиц с поверхности t, ms инициированный налетающим током m ( ns ) на поверхность S. Его можно представить в виде разложения по базису m :

J nt + ( nt ) = m ( it ), r ms ms (1.26) imt m = где rimt = m ( nt ) J nt + ( nt )d nt. Отличительной чертой метода матриц отклика ms ms является то, что связь между током втекающих частиц и током вытекающих частиц можно определить следующим уравнением J n+ = R n J n (1.27) где R n – матрица отклика (response matrix). Вектор J n представляет собой вектор коэффициентов разложения для всех поверхностей n-го элемента J n = ( jn 0 1 0 jnS )T.

M jn1... jn 2 jn 2...

Поскольку вытекающие токи для ячейки являются втекающими токами для смежных ячеек, вводится матрица связи M (connectivity matrix) такая что J = MJ +, (1.28) где J + и J – блочные вектора, состоящие из векторов J n+ и J n соответственно для всех ячеек. Матрица M имеет не более чем один ненулевой элемент на строку и столбец. Таким образом, уравнение метода матриц отклика имеет вид J = MRJ, (1.29) где R – это блочная диагональная матрица, состоящая из блоков R n. Данный метод реализован, например, в программе HELIOS [54].

В дополнение к рассмотренным методам следует отметить, что Б.П. Кочуровым развивался схожий с вышеперечисленными методами подход – метод -матриц [55–57]. В качестве отправной точки для этого метода служит метод источников–стоков [58]. Итоговые уравнения данного метода по своей сути близки к уравнениям МПГ. Эти уравнения были получены только для трех поперечных пробных матриц, что является основным ограничением этого метода.

1.3 Метод поверхностных гармоник После рассмотренных выше проекционно-сеточных методов рассмотрим метод поверхностных гармоник, который также можно отнести к данному классу методов и которому посвящена данная работа. Данный метод был предложен профессором Н.И. Лалетиным в 1976 году в работе [6].

1.3.1 Основы метода поверхностных гармоник и определения Приведем здесь кратко описание базовых моментов МПГ и определения, использующиеся в работах по МПГ.

Пространственная декомпозиция. Так же как и в представленных в обзоре проекционно-сеточных методах, в МПГ расчетная область V разбивается на пространственные подобласти Vk (ячейки). В результате глобальное решение (r, ) (групповой вектор) представляет собой сумму решений по всем ячейкам:

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

Пробные функции. Групповое решение в каждой ячейке представляется в виде линейной комбинации пробных матриц и амплитуд:

L (kL ) (r, ) = (kl ) (r, )I (kl ), (1.31) l = где (kl )(r, ) – l-ая пробная матрица для k-ой ячейки, I (l ) – l-ая амплитуда k-ой k ячейки в виде группового энергетического вектора. Отметим, что в дальнейшем будем использовать символ ( L ) (r, ) для решения во всей активной зоне полученной при использовании L + 1 количества пробных матриц.

Остановимся подробнее на пробных матрицах МПГ. Пробные матрицы МПГ состоят из G пробных векторов (kl,)g (r, ) ( ) (kl ) (r, ) = (kl,)1(r, ) (kl,)2 (r, ) (kl,)G (r, ). (1.32) Каждый (kl,)g (r, ) пробный вектор представляет собой решение уравнения переноса нейтронов в k-ой ячейке с граничным условием в виде тока, заданного в виде l-ой координатной функции Wl (rs ) на поверхности ячейки rs в g’-ой энергетической группе (в остальных энергетических группах ток равен нулю).

Определим следующие связанные с пробными функциями матрицы (l ) элементы k которых имеют вид [7]:

( (kl ) ) gg = drs Wl (rs )( k (rs )) gg, (l ) (1.33) r N l rs где ((kl ) (r )) gg = 3 d ( n)2 ((kl ) (r, )) gg, N l = Wl (rs )Wl (rs )drs, где n – внутренняя нормаль к границе ячейки. Таким образом, ( (kl ) ) gg представляет собой деленный на норму N l интеграл по границе области ячейки и по от решения уравнения переноса в g-ой энергетической группе с весом координатной функции Wl (rs ) и множителя 3( n) 2 при граничном условии, распределенном по l-ой координатной функции в g’-ой энергетической группе.

Координатные функции. Приведем описание координатных функций для ячеек с квадратной внешней границей, т.к. дальнейшее получение уравнений проводится непосредственно для квадратной решетки блоков. Координатные функции МПГ Wl (rs ) описываются следующей формулой (например [9]):

cos ( l j ) Wl (rs ) = Pp ( j ), p=0, 1, 2,…, P;

l=0, 1, 2,…L (1.34) sin ( l j ) где j 1 2 a j = (1) a [rs 2 (2 j 1)], 1 j 1, 0 rs Ma = 2 ( j 1) jM l – номер координатной функции l=0, 1,…, L, Pp ( j ) – полиномы Лежандра p-го порядка, M – число боковых сторон (граней) ячейки (в случае квадратной решетки блоков равно M = 4), a – длина боковой стороны ячейки, j – угол между нормалью, построенной из центра ячейки на j-ю боковую сторону ячейки, и осью абсцисс (см. рисунок 1.1).

Рисунок 1.1 – Параметры координатных функций В таблице 1.1 представлены первые восемь координатных функций для квадратной ячейки [8].

Таблица 1.1 – Первые восемь координатных функций для квадратной ячейки l l Wl (rs ) Wl (rs ) P1 ( j ) cos(2 j ) 0 1 cos( j ) P1 ( j ) cos( j ) 1 sin( j ) P1 ( j ) sin( j ), 2 cos(2 j ) P1 ( j ) 3 На рисунке 1.2 представлена схема втекания токов для первых восьми координатных функций для ячейки с квадратной границей.

Рисунок 1.2 – Схема втекания токов, соответствующая первым восьми координатным функциям для ячейки с квадратной границей 1.3.2 Обзор программных комплексов SUHAM и SVS В данном параграфе кратко рассмотрим программные комплексы SUHAM и SVS, реализующие стационарные уравнения метода поверхностных гармоник.

SUHAM. «Комплекс программ SUHAM предназначен для реализации конечно-разностных уравнений метода поверхностных гармоник, для расчета нейтронно-физических процессов в ядерных реакторах с треугольной и квадратной решетками блоков (ТВС)» [14, c.10]. «Модули комплекса SUHAM, реализующие конечно-разностные уравнения МПГ, применяются для решения многогруппового уравнения переноса нейтронов во всем объеме активной зоны реактора методом поверхностных гармоник» [14, c.10]. Комплекс состоит из двух независимых частей SUHAM-W [11] и SUHAM-U [12].

SUHAM-W взаимодействует с компонентами программы WIMS-SH (раннее название WIMS-SU) [61] для подготовки групповых сечений изотопов и материалов и использует в качестве библиотеки ядерных данных библиотеку программы WIMS-D [62]. К основным возможностям SUHAM-W по решению двумерного конечно-разностного уравнения МПГ можно отнести:

1. Решение группового уравнения диффузии (прямого и сопряженного) с различными граничными условиями (традиционный метод гомогенизации);

2. Решение уравнений МПГ с разным числом пробных матриц (3, 4, 7 и 8 для квадратной решетки и 3, 4, 5 и 6 для треугольной решетки);

3. Решение сопряженных уравнений МПГ с тремя пробными матрицами.

SUHAM-U взаимодействует с компонентами комплекса UNK [63] для подготовки групповых сечений изотопов и материалов и решения уравнений изотопной кинетики, использует в качестве библиотеки ядерных данных библиотеку программы UNK. Комплекс SUHAM-U позволяет решать двумерные и трехмерные нейтронно-физические задачи. Двумерный модуль имеет следующие возможности:

1. Решение группового уравнения переноса в двухэтапном расчете топливных сборок с помощью МПГ с разным числом пробных матриц: от 3-х до 8-и для квадратной решетки и от 3 до 6 для треугольной решетки.

2. Трехэтапный расчет активной зоны ядерного реактора с шестигранными ТВС.

3. Расчет выгорания ТВС ВВЭР-1000.

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

Верификация программного комплекса SUHAM проводилась на следующих задачах: ТВС PWR [64], ТВС ВВЭР-1000 с урановым и MOX топливом [65], бенчмарк C5G7 [66], бенчмарк VENUS-2 [67], полномасштабная двумерная зона ВВЭР [16], полиячейки и модельные сборки РБМК [14], двумерная зона БРЕСТ-ОД-300 [14], реактор ГТ-МГР [15].

SVS. Комплекс программ SVS (Surface Values System) и его составные части в виде самостоятельных модулей SVL [13, 68] и SVC [13] «используются в исследованиях по физике ВВЭР и РБМК и позволяют оценивать вклад от различных приближений, используемых в проектных исследованиях и эксплуатационных расчетах» [1, c.8]. SVS позволяет проводить исследования как локального нейтронного поля, так и всей активной зоны целиком.

SVL (Surface Values Lattice) представляет собой решеточный код который «предназначен для решения задач теории переноса нейтронов в ячейках и полиячейках для основных типов решеток ядерных реакторов» [1, c.12]. Данная программа разрабатывалась в основном для реакторов легководного типа, но также применима для реакторов графитового и тяжеловодного типов. SVL использует МППИ для расчета характеристик ячеек и МПГ для расчета кассет (полиячеек). Программа SVL позволяет проводить расчеты на критичность и расчеты на выгорание ячеек в цилиндрической и кластерной геометриях, 2-D и 3-D сборок для гексагональной, квадратной и треугольной решеток блоков. Также следует отметить что программа позволяет работать со сборками, имеющими нарушения регулярности решетки.

SVC (Surface Values Core) «предназначена для решения полномасштабных нейтронно-физических расчетов активных зон реакторов ВВЭР» [1, c.17] используя МПГ. Входные данные для SVC готовит программа SVL. В обеих программах используется последовательный подход, т.е. «используются методы, в которых начальное приближение удовлетворяет требованиям проектных расчетов, а результаты расчетов получаются по точности сравнимые с реперными, что можно проверить, переходя к следующим приближениям» [1, с.13].

Общая схема расчета по SVS выглядит следующим образом: расчет характеристик кассет, расчет активной зоны, восстановление потвэльных полей, учет обратных связей [1]. На рисунке 1.3 представлена блок схема комплекса SVS [13].

SVL SVC Feedbacks Рисунок 1.3 – Блок-схема комплекса SVS Комплекс прошел множественную верификацию, в том числе на ячейках реактора BBЭР-1000 [69], сборках реакторов типа ВВЭР (в частности ZR6) [13, 70]. Также следует отметить, что проводился расчет кампании реактора ВВЭР-1000 первого блока Волгодонской АЭС [13].

В дополнение к рассмотренным программам, следует отметить, что имеется также гетерогенная версия программы STEPAN, использующая МПГ [2, 71].

1.3.3 Обзор работ по нестационарным уравнениям МПГ В данном параграфе рассмотрим основные работы [72–77] в которых проводились попытки получения конечно-разностных уравнений метода поверхностных гармоник для нестационарного уравнения переноса нейтронов.

Первая такая попытка была сделана в 1990 г. Н.В. Султановым и В.Ф.

Бояриновым. Суть этой работы заключается в том, что члены уравнения переноса, связанные с запаздывающими нейтронами и с производной по времени были перенесены в правую часть и представлены в виде источников. Введены пробные функции от объемных источников (производная по времени и запаздывающие нейтроны). По аналогии со стационарным случаем, невязка функционала возникала только на границах ячеек. Как итог были получены конечно разностные уравнения МПГ для пространственной кинетики в 1-D (плоской), 2-D (квадратной и треугольной) и 3-D геометриях. Дальнейшего развития этот подход не получил и опубликован не был.

В 1993 г. в работе [72, 73] был предложен другой подход. Главной особенностью этого подхода является то, что при построении конечно разностных уравнений МПГ используются только пробные функции стационарного уравнения переноса нейтронов. Как итог, помимо поверхностной невязки на границах ячеек возникает еще и объемная невязка, которая связанна с производной функции распределения нейтронов по времени и с запаздывающими нейтронами. Конечно-разностные уравнения получаются из минимизации суммарной невязки.

В работах были получены конечно-разностные уравнения [74–76] пространственной кинетики метода поверхностных гармоник с помощью совместного рассмотрения прямой и сопряженной задач в одинаковом приближении МПГ.

В работе [77] представлены полученные еще в 90-х годах конечно разностные уравнения пространственной кинетики метода поверхностных гармоник для плоской одномерной решетки.

Следует отметить, что ни в одной из этих работ полученные уравнения не были программно реализованы и верифицированы. Первой работой, в которой проведена программная реализация нестационарных уравнений МПГ, является [78]. Программная реализация проведена для конечно-разностных уравнений пространственной кинетики в одномерной плоской геометрии.

Следует также отметить, что в [2] без вывода приведены двухгрупповые 11-точечные по пространству нестационарные уравнения МПГ. Реализованные в программе STEPAN эти уравнения позволили повысить надежность расчета ряда нейтронно-физических функционалов.

1.4 Обзор методов решения нестационарного уравнения переноса нейтронов Нестационарное уравнение переноса нейтронов может быть (1.1) представлено в виде (t ) = H(t )(t ), (1.35) & где (t ) - первая производная по времени функции распределения нейтронов, & H(t ) - полный оператор переноса в матричной форме. В уравнение (1.35) все независимые аргументы, кроме времени, для простоты опущены. Имеется ряд подходов для работы с временной зависимостью в уравнение (1.35), которые подразделяются на две категории.

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

Вторая категория – это непрямые методы использующие процедуру разделения переменных. Среди этих методов можно назвать улучшенный квазистатический метод [27] и SCM [80] методы. Далее кратко рассмотрим нестационарные методы решения переноса нейтронов.

1.4.1 Полностью явный метод При использовании полностью явного метода Эйлера (fully explicit method) уравнение (1.35) принимает вид (t n+1 ) (t n ) = H(t n ) (t n ) (1.36) t и решение представляется в виде (t n+1 ) = (I + tH (t n )) (t n ). (1.37) Данный подход требует минимума вычислительных затрат на каждый временной шаг по сравнению с другими конечно-разностными методами, однако он неустойчив если шаг t недостаточно мал. Данный факт требует такого уменьшения временного шага, что преимущество, даваемое простотой этого метода, как правило, компенсируется большим количеством временных шагов [80].

Дополнительно отметим, что В.И. Лебедевым был разработан вариант устойчивой разностной явной схемы с переменными шагами по времени реализованный в программе DUMKA [81].

1.4.2 Полностью неявный метод Проблема устойчивости явного метода может быть решена посредством использования неявного метода. При использовании неявного метода Эйлера (fully implicit method) уравнение (1.35) принимает вид (t n+1 ) (t n ) = H (t n+1 ) (t n+1 ) (1.38) t и решение представляется в виде (tn+1 ) = (I tH(t n+1 ))1 (t n ). (1.39) Данный подход является устойчивым, однако требует обращения матрицы I tH (t n+1 ) на каждом временном шаге, что влечет большие вычислительные затраты.

1.4.3 -метод Наиболее известным среди полунеявных методов является -метод (-method) [80]. В -методе уравнение (1.35) представляется в виде (t n+1 ) (t n ) = H(t n+1 ) (t n+1 ) + (I )H(t n ) (t n ), (1.40) t где – диагональная матрица, элементы которой ii удовлетворяют неравенству 0 ii 1. При = 0 метод переходит в полностью явный метод, при = I – полностью неявный метод, а при ii = 0,5 (для всех i) переходит в схему Кранка Николсон (Crank-Nicholson). Решение в -методе имеет вид (t n+1 ) = [I tH(t n+1 )] [I + t (I )H(t n )](t n ).

(1.41) Матрица выбирается так чтобы повысить точность конечно-разностной аппроксимации и может изменяться от шага к шагу.

1.4.4 Метод переменных направлений Рассмотрим суть метода переменных направлений (МПН) [33]. Метод был предложен D.W. Peaceman, H.H. Rachford [82] и J. Douglas [83]. Идея МПН заключается в том, чтобы в одном из направлений использовать полностью неявную схему (в остальных направлениях полностью явную схему) и изменять это направление от временной точки к временной точке [80]. Для простоты будем рассматривать двумерный случай. Матрица H разделяется по следующему принципу H = H1 + H 2 = H 3 + H 4 (1.42) где H1 = X + R 1, H 2 = Y + R 2, H 3 = Y + R 3, H 4 = X + R 4, X – симметричная трехдиагональная матрица, отвечающая за перенос в одном направлении, Y – симметричная трехдиагональная матрица, отвечающая за перенос в перпендикулярном первому направлении [84].

В методе переменных направлений конечно-разностная схема имеет вид (t n +1 / 2 ) (t n ) = H 1 (t n +1 / 2 ) + H 2 (t n ) t / (1.43) (t n +1 ) (t n +1 / 2 ) = H (t ) + H (t n +1 / 2 ) n + t / 3 а, значит, решение представляется в виде t t (t n+1/ 2 ) = I H1 I + H 2 (t n ) 2 (1.44) t t (t n+1 ) = I 2 H 3 I + 2 H 4 (t n+1/ 2 ) Разложение матрицы H на две матрицы, отвечающие за перенос в различных направлениях, производится с целью попеременного применения явной и неявной схем для различных направлений в соответствии со схемой Писмена-Рекфорда [82]. Матрицы R i выбираются таким образом, чтобы сделать процесс обращения матриц в системе (1.44) наименее трудоемким [84].

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

Данный метод вместе с -методом был применен в программе TWIGL для решения двумерного группового уравнения диффузии [19].

1.4.5 Улучшенный квазистатический метод Основная идея улучшенного квазистатического метода (improved quasi static method, IQS) заключается в том, что выделяется две функции, имеющие различную скорость изменения во временном масштабе: быстро меняющаяся амплитуда (амплитудный фактор) и форм-функция, которая обычно изменяется медленнее. В этом методе нестационарная функция распределения (r, E,, t ) ~ представляется в виде произведения форм-функции (r, E,, t ), которая слабо изменяется с течением времени, и амплитуды T (t ), которая зависит только от времени и с течением времени, обычно, изменяется быстрее [85]:

~ (r, E,, t ) = (r, E,, t )T (t ). (1.45) Это представление функции распределения дополняется следующим условием нормировки:

~ (r, E,, t ) + (r, E, ) dVdEd = const, (1.46) v( E ) где + (r, E, ) – функция, сопряженная стационарной функции (r, E, ).

Уравнение для расчета форм-функции (уравнение формы) получается после подстановки (1.45) в нестационарное уравнение переноса нейтронов (1.1):

~ 1 ~ T (t ) (r,, E, t ) ~ ~ (r,, E, t ) + + (r,, E, t ) + t (r, E, t )(r,, E, t ) = v( E ) t t T (t ) = dE d s (r, E E,, t )(r,, E, t ) + ~ (1.47) p (E) 1J dE d (1 (r, E ) ) f (r, E, t )(r,, E, t ) + d, j ( E ) j C j (r, t ) ~ + 4 4 j = После ряда преобразований нестационарного уравнения переноса нейтронов с использованием формулы (1.45) и сопряженного стационарного уравнения переноса нейтронов получают уравнение для расчета амплитуды:

dT (t ) (t ) (t ) = (t ) T (t ) + i ci (t ). (1.48) dt i ~ где функционалы (t ), (t ), (t ) и ci (t ) определяются через функции (r, E,, t ) и + (r, E,, t ) [27, 85]. Уравнение для ci (t ) имеет вид dci (t ) i (t ) = ci (t ) i ci (t ). (1.49) (t ) dt Обычно форм-функция изменяется более медленно во времени, чем амплитуда и, следовательно, форм-функция рассчитывается с более крупными временными шагами, чем амплитуда. Результат решения уравнения формы используется для расчета параметров точечной кинетики, которые используются при решении уравнения (1.48).

Дополнительно отметим, что если положить нулю производную ~ (r,, E, t ) в уравнении (1.47), то придем к квазистатическому приближению, t ~ использующему факт слабой зависимости форм-функции (r, E,, t ) от времени.

1.4.6 SCM метод Жесткость системы пространственной кинетики определяется временем жизни мгновенных нейтронов, которое на несколько порядков меньше времени жизни запаздывающих нейтронов [80]. Метод SCM (stiffness confinement method) был разработан [86] для уменьшения ограничений на временной шаг, и, следовательно, повышения производительности вычислений. В данном методе уравнение переноса нейтронов и уравнение для предшественников рассматриваются отдельно посредством введения динамических частот (dynamic frequencies) 1 g 1 C m (r, t ), C (r, t ) m g. (1.50) g t C m t Используя эти определения, производится замена производных в исходной системе уравнений. В результате получается уравнение, аналогичное стационарному уравнению переноса, но с модифицированными членам поглощения и деления, которые включают динамические частоты. Для того чтобы получить решение во времени, делается оценка динамических частот, решается модифицированное уравнение переноса, пересчитываются концентрации предшественников, далее улучшенная оценка динамических частот рассчитывается с использованием новых потоков и концентраций предшественников. Такие итерации проводятся до сходимости [86]. Данный метод реализован в таких кодах как SPNOVA [87] и PANTHER [88].

В дополнении к рассмотренным методам необходимо отметить класс линейных многошаговых методов (linear multistep methods) [32, 34]. В таких методах значение решения в n-ой временной точке определяется значениями в ряде предыдущих точек. В качестве примеров методов из этого класса можно назвать методы Адамса-Башфорта (Adams-Bashforth), Адамса-Мултона (Adams Moulton) и др. Недостатком этих методов является возрастание величины необходимой памяти по сравнению с одношаговыми методами для хранения решений, полученных на предыдущих шагах, что особо важно при моделировании полномасштабных моделей. Также отметим наличие класса методов типа прогноз-коррекция (predictor-corrector methods) [32] обеспечивающие контроль ошибки на каждом шаге расчета, что также является важной чертой современного кода.

1.5 Обзор программ для решения нестационарного уравнения переноса нейтронов Как уже отмечалось ключевой задачей для анализа нейтронных переходных процессов в ядерном реакторе является решение нестационарного уравнения переноса нейтронов с запаздывающими нейтронами. Исторически, из-за ограниченных возможностей вычислительных технологий на ранних этапах атомной науки, нейтронные переходные процессы описывались в значительной степени упрощенно, с использованием 1D подходов, либо с использованием уравнений точечной кинетики [3, 89].

Со значительным ростом вычислительных возможностей компьютерной техники в течение второй половины XX века и развитием вычислительных алгоритмов были разработаны диффузионные коды для решения пространственно-временного группового уравнения переноса нейтронов с запаздывающими нейтронами в диффузионном приближении. Среди таких кодов наиболее известными иностранными кодами являются TWIGL [19], LUMAC [90], MATIN [91], SADI [84], DIF3D-K [92] и множество других. Среди российских кодов, проводящих расчеты ядерных реакторов и использующих данные приближения можно отметить STEPAN [2], БИПР-8КН [93], JAR-IQS [94], ShIPR [95], ГЕФЕСТ [96] и др. Коды данного поколения содержат, в том или ином количестве такие приближения как диффузионное приближение, пространственная гомогенизация и расчет в малом числе энергетических групп.

Несмотря на то, что качество расчетов, проводимых этими кодами, часто является достаточным для анализа многих сценариев поведения различных типов реакторов, ряд исследователей (например, [3–5]) указывает на возможность получения улучшенных результатов путем решения более универсального газокинетического уравнения, особенно для анализа безопасности.

К одним из первых программ, разработанных для решения пространственно-временного уравнения переноса нейтронов, были программы TIMEX [97] и TRANZIT [98] для 1-D и 2-D геометрий в начале 70-х годов. На текущий момент уравнение переноса нейтронов, как правило, применяется в основном в стационарном расчете, по причине высоких вычислительных затрат для нестационарных задач даже на современных компьютерах. Вместе с тем, растущая необходимость применения уравнения переноса нейтронов в расчете переходных процессов, особенно ввиду строгих стандартов безопасности, предъявляемых к новым разрабатываемым реакторам, привела к созданию ряда современных программных комплексов, позволяющих решать пространственно временное уравнение переноса нейтронов с запаздывающими нейтронами. Далее рассмотрим кратко наиболее известные из них.

DORT-TD. Код DORT-TD [3, 4] представляет собой нестационарную расширенную версию хорошо известного кода DORT [99]. DORT является двухмерным кодом (также есть одномерный плоский вариант), разработанным для x-y, r-z или r- геометрий. Он может быть использован для решения прямой или сопряженной форм уравнения переноса нейтронов. Решается либо диффузионное уравнение, либо уравнение переноса нейтронов методом дискретных ординат.

TORT-TD. Трехмерный нестационарный код TORT-TD [100, 101] проводит расчеты с помощью метода дискретных ординат. Данный код основывается на стационарном коде TORT [102, 103]. TORT-TD решает нестационарное групповое уравнение переноса нейтронов в декартовых координатах и цилиндрической геометрии. Стабильность временной схемы достигается посредством использования полностью неявной схемы.

Код проводит решение трехмерного PARTISN. PARTISN [104] нестационарного уравнения переноса нейтронов в декартовых координатах методом дискретных ординат. PARTISN использует полунеявную временную схему.

EVENT. Код EVENT [45] использует методы конечных элементов и сферических гармоник для решения нестационарного уравнения переноса нейтронов.

BARS. Программа BARS [105] использует метод -матриц Кочурова и улучшенный квазистационарный метод в нестационарных расчетах.

В дополнение к этим детерминистическим кодам, следует отметить, наиболее известные программы по решению нестационарного уравнения переноса нейтронов методом Монте-Карло, а именно TDKENO [106] и TRIPOLI [107].

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

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

• упорядоченность координатных функций по степени их влияния на нейтронно-физические функционалы, что позволяет достигать приемлемой для практики точности основных функционалов уже в низших приближениях метода;

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

В дополнение к вышеперечисленному имеется многолетний опыт использования МПГ в НИЦ «Курчатовский институт» с реализацией в программных комплексах и Метод показал качественные результаты как на SUHAM SVS.

идеализированных бенчмарках (например C5G7), так на реакторах ВВЭР-1000, РБМК-1000, БРЕСТ-ОД-300, ГТ-МГР и БН-1200. При этом МПГ демонстрировал преимущество по вычислительным затратам по сравнению с детерминистическими методами при сохранении качества их результатов.

Вышеперечисленные факты делают перспективным применение МПГ в нестационарном случае.

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

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

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

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

Глава 2 Нестационарные уравнения метода поверхностных гармоник 2.1 Двумерные нестационарные уравнения МПГ В данной главе проводится детальный вывод двумерных конечно разностных нестационарных уравнения метода поверхностных гармоник для трех и четырех пробных матриц в квадратной решетке блоков. Получение данных уравнений проводится традиционным для МПГ методом минимизации невязки.

Также основываясь на полностью неявной схеме Эйлера, представлено получение итерационной схемы решения уравнений МПГ для трех пробных матриц.

Решение (r,, t ) уравнений (1.3) и (1.4) будем представлять в виде функции ( L ) (r,, t ) ( L ) (r,, t ) (r,, t ), (2.1) т.е. ( L ) (r,, t ) удовлетворяет уравнению (1.3) с некоторой невязкой ( L) [ j C j (r, t ) d, j ] = R.

1J (r,, t ) + L ( L ) (r,, t ) K ( L ) (r,, t ) (2.2) v inv t 4 j = Введем сопряженную групповую функцию + (r, ) (вектор-столбец), которая является точным решением уравнения, сопряженного стационарной форме уравнения (1.3) во всем объеме V рассчитываемого объекта, а именно L+ + (r, ) = K + + (r, ), (2.3) где L+ + (r, ) = + (r, ) + tr (r ) + (r, ) т,tr (r ) + (r ), s d f (r) (r, ).

K + + (r, ) = + т 4 Умножим скалярно уравнение (2.2) на + (r, ), уравнение (2.3) на ( L ) (r,, t ) и вычтем второе из первого. Результат проинтегрируем по всему объему V и по всем направлениям, при этом отметим, что, т.к. функция ( L ) (r,, t ) точно не удовлетворяет условиям непрерывности на границах ячеек и граничным условиям на внешней границе реактора, то равенство + (r,, t ) L+ + (r, )d dV (r, ) L( L ) (r,, t )d dV = ( L) (2.4) 4V 4V несправедливо, т.к. члены, содержащие производные, обращаются в бесконечность на внутренних поверхностях, поскольку на этих поверхностях решение ( L ) (r,, t ) имеет разрывы [27].

[( )]) ( ) ( + ;

grad + (r, );

(kL ) (r,, t ) drd + = (r, );

R drd k k Vk 4V (2.5) + + (r, );

v inv (kL ) (r,, t ) + t C j,k (r, t ) d, j,k drd t k 4 j 4 k Vk где Vk – объем k-ой ячейки, (a;

b ) – скалярное произведение векторов a и b. Далее ( ) + величину будем обозначать, придерживаясь классических работ (r, );

R drd 4V [6–10], как N. Отметим, что величина N представляет собой сумму двух компонентов, связанных с поверхностной невязкой [( )]) ( ;

grad + (r, );

(kL ) (r,, t ) drd N s = (2.6) k k Vk и объемной невязкой + k (r, );

v inv (kL ) (r,, t ) + 1 C j,k (r, t ) d, j,k drd.

N v = (2.7) t 4 j t 4 k Vk Далее остановимся подробнее на этих величинах.

2.1.1 Поверхностная невязка Величина N s (2.6) связана с поверхностной невязкой, т.е. невязкой, обусловленной тем фактом, что ( L ) (r,, t ) точно не удовлетворяет уравнению переноса нейтронов на границах ячеек из-за ограниченного числа пробных матриц. Значения данной величины приведено, например, в работе [9] и в данном параграфе приводится вывод данного представления только для полноты описания уравнений метода поверхностных гармоник. Преобразуем (2.6) используя теорему Остроградского-Гаусса [108] [ ] dr (n ;

) + k (r j, ) ( L ) k (r j,, t ) ( L ) k (r j,, t ).

N s = 2 (2.8) d j kj j =1 j k В последнем члене в нижних индексах номер k’ – номер ячейки соседней с k-ой ячейкой, для которой грань j’ совпадает с гранью j ячейки k, при этом j’ = j + 2, для j = 1, 2 и j’ = j – 2, для j = 3, 4. Знак минус перед вторым членом появляется из-за того что нормали n k j заменены на нормали n kj. Коэффициент появляется из-за перегруппировки членов невязки.

Разложим прямую (kL ) (r,, t ) и сопряженную + ( L ) (r, ) функции по k формулам L (kL) (r,, t ) = (kl ) (r, )I (kl ) (t ), (2.9) l = L + + + k ( L ) (r, ) = I k (l ) k т (l ) (r, ). (2.10) l = Здесь следует отметить, что в формуле (2.9) опущена временная зависимость у пробных функций (kl ) (r, ), т.к. на каждом временном интервале ts в данном подходе они являются решениями стационарного уравнения переноса нейтронов, соответствующего некоторой временной точке на интервале ts (см. далее уравнение (2.87).

+ Разложение функций (kl ) (r j, ), k т ( l ) (r j, ) по полиномам Лежандра на границе ячеек имеет вид 1 2n + 2 Pn (µkj )(kl,)n (rj ), (kl ) (r j, ) = (2.11) 2 n = 1 2n + Pn ( µ kj ) +,тn( l ) (r j ), + k т (l ) (r j, ) = (2.12) k 2 n=0 где Pn ( µkj ) – полином Лежандра n-го порядка, µ kj = ( ;

n kj ), P (µ ) (kl,)n (r j ) = (l ) (r j, ) d, (2.13) n k k P (µ ) + +(l ) k,(nl ) (r j ) = (r j, ) d. (2.14) n k k Подставим разложения (2.9), (2.10) в Также воспользуемся (2.8).

разложением (2.12), оставляя только первые два члена [7], проинтегрируем по и учтем (2.13), получим L 1 L L dr j I k+(l) 2 k+,т0(l) (r j )(n kj ;

) (kl ) (r j, )I (kl ) (t ) (kl) (r j, )I (kl) (t ) + N s = k j =1 l=0 l = 0 l =0 4 4 j L L + (n kj ;

) 2 k,т1(l ) (r j ) (kl ) (r j, )I (kl ) (t ) (kl) (r j, )I (kl) (t ).

+ l =0 l = 2 (2.15) Далее введем обозначения для средних по грани ячейки тока и уровня нейтронов [9] L dr j (n kj ;

) (kl ) (r j, )I (kl ) (t ) d j kj = 4 (2.16) l =0 S j L dr j (n kj ;

) 2 (kl ) (r j, )I (kl ) (t ) d f kj = 4, (2.17) l =0 S j где S – длина грани ячейки. Также отметим, что для первых четырех сопряженных пробных матриц можно записать + + k,тn(l ) (r j ) = Wl (r j ) k,тn( l ), n = 1, 2. (2.18) Подставим выражения (2.16) и (2.17) в (2.15) и учтем (2.18). Отметим, что знак минус при токе jk j появляется из-за обратной нормали, а в уровне fk j знак сохраняется из-за наличия квадрата (n kj ;

) 2. Таким образом, выражение (2.15) принимает вид:

[ ] S I k (l ) (t ) Wl (r j ) k,т0(l ) (j kj + j k j ) + Wl (r j ) k,т1(l ) (f kj f k j ).

4 L 1 + + + N s = (2.19) 4 4 k j =1 l = 0 Здесь использован факт, что Wl (r j ) постоянен на каждой грани ячейки для l = 0, …, 3. Далее отметим, что L (l ) L Wl (r j ) j kj = Wl (r j ) (n kj ;

) (kl) (r j, )I (kl) (t ) d = Wl (r j ) (n kj ;

) (kl ) (r j, )d I k (t ) = l = 0 l = 0 4 L = Wl (r j ) Wl (r j )I (kl ) (t ) l = Тогда (2.19) примет вид 1 + L 4 I k ( l ) (t ) k,т0( l ) N l I (kl ) (t ) + S Wl (r j ) jk j + k,т1( l ) N l (kl ) I (kl ) (t ) S Wl (r j )fk j + + N s = k l =0 2 4 4 j =1 j = (2.20) где 4 S, l = 0, Nl = (2.21) 2 S, l = 1, 2.

2.1.2 Объемная невязка Величина N v связана с объемной невязкой. Данная невязка возникает в нестационарном случае и обусловлена тем фактом, что пробные матрицы рассчитываются посредством решения стационарного уравнения переноса нейтронов. Как уже отмечалось, величина N v имеет вид:

t C j,k (r, t ) d, j,k drd + N v = (r, );

v inv ( L ) k (r,, t ) +. (2.7) k t 4 j k Vk Разложим ( L ) (r,, t ) в соответствии с формулой (2.9), а функцию k ( L ) (r,, t ) по + формуле (2.10).

Отметим, что в случае объемной невязки разложение матриц (kl ) (r, ), (l ) + (r, ) в ряд по полиномам Лежандра производится в объеме всей ячейке:

k 1 2n + Pn (µ ) (kl,)n (r ), (kl ) (r, ) = (2.22) 2 n=0 1 2n + + т (l ) (r, ) = 2 Pn (µ) +,тn(l ) (r), (2.23) k k 2 n = где Pn (µ) - полином Лежандра n-го порядка, + + (kl,)n (r ) = (kl ) (r, )d, k,(nl ) (r ) = k ( l ) (r, )d.

4 Ограничим сумму ряда (2.23) до двух членов [7] и подставим формулы (2.9), (2.10) и (2.23) в (2.7) L ( ) L l =0 k dr I + (l ) (t ) P0 (µ ) +,т0(l ) (r ) + 3P1 (µ) +,т1(l ) (r ) v inv (kl ) (r, ) I (kl ) (t ) + 4 N v = d t k k l =0 k Vk L + ( l ) ( ( )) + I k (t ) P0 (µ) k,т0(l ) (r ) + 3P1 (µ) k,т1(l ) (r ) t C (r, t ) d, j,k + + j,k l =0 j (2.24) Далее отметим, что полиномы Лежандра ортогональны [108], т.е.

P (µ)P (µ)dµ = 2n + 1, (2.25) n m nm где nm - символ Кронекера. Также отметим, что 2 1 d = d dµ = 2 dµ. (2.26) 4 1 Распишем отдельно интеграл по d в первом члене суммы с учетом (2.25) и (2.26), заменяя (kl ) (r, ) по формуле (2.22).

L + ( l ) L 1 2n + ( ) Pn (µ kj ) (kl,)n (r ) I (kl ) (t ) d I k (t ) P0 (µ kj ) k,т0( l) (r ) + 3P1 (µ kj ) k,т1( l) (r ) v inv + + = l=0 2 t l =0 n = L + L 1 2n + 1 Pn (µ kj ) (kl,)n (r ) I (kl ) (t ) d + = I k (l ) (t ) P0 (µ kj ) k,т0(l ) (r ) v inv + 2 t 4 l = 0 l =0 n = L + ( l ) L 1 2n + 1 Pn (µ kj ) (kl,)n (r ) I (kl ) (t ) d I k (t )3P1 (µ kj ) k,т1( l ) (r ) v inv + + = l=0 2 t l =0 n= L + L 1 = 4 I k ( l ) (t ) k,т0( l ) (r ) v inv (kl,)0 (r ) I (kl ) (t ) + + 4 t l = 0 l =0 L + L 3 + 4 I k ( l ) (t ) k,т1( l ) (r ) v inv (kl,1 (r ) I (kl ) (t ) + ) 4 t l = 0 l =0 Проинтегрируем второй член L + ( l ) ( ) 41 t C I k (t ) P0 (µ kj ) k,т0( l) (r ) + 3P1 (µ kj ) k,т1( l ) (r ) + + = (r, t ) d, j,k d j,k l =0 j L + ( l ) I k (t ) P0 (µ kj ) +,т0(l ) (r ) C j,k (r, t ) d, j,k d = = 4 j t k l = L = I + (l ) (t ) +,т0(l ) (r ) C j,k (r, t ) d, j, k j t k k l = Таким образом, при L = L объемная невязка равна 1 L L dr (2n + 1) I +(l ) (t ) +,тn(l ) (r ) v inv (kl,n (r ) I (kl ) (t ) + n = N v = ) t k k l =0 l =0 4 k Vk (2.27) L + + I k (l ) (t ) k,т0( l ) (r ) C j,k (r, t ) d, j,k + j t l =0 Далее, при условии (r ) (kl,n (r )dr = 0 если l l.

+ т (l ) ) (2.28) k,n Vk получим 1 L I +(l ) (t ) (2n + 1) +,тn( l ) (r ) v inv (kl,)n (r )dr I (kl ) (t ) + k n = N v = t k 4 k l =0 Vk (2.29) + +,т0(l ) (r ) C j,k (r, t ) d, j,k dr j t k Vk 2.1.3 Вывод уравнений МПГ Как было показано в предыдущих двух параграфах величину (2.5) можно записать как L +(l ) 1 + т (l ) I k (t ) k, 0 N l I (kl ) (t ) + S Wl (r j ) j k j + N = l = 4 k j = + k,т1(l ) N l (kl ) I (kl ) (t ) S Wl (r j )f k j + j = (2.30) 1 L + I k (l ) (t ) (2n + 1) k,тn(l ) (r ) v inv (kl,)n (r )dr I (kl ) (t ) + + + n =0 t l =0 Vk C j,k (r, t ) d, j,k dr + I +( l ) (t ) +,т0( l ) (r ) t k k j Vk Получение уравнений МПГ будем проводить традиционным для МПГ способом – методом минимизации невязки. Приравняем нулю коэффициенты при ( ) I k ( l ) +,(nl ). Для этого в объемной части невязки вставим выражения +,т0( l ) +,т0(l ) + и k k k +,т1( l ) ( k,т ( l ) ) после вектора I k ( l ) (t ). В результате получим при l = 0 получим 1 + + k C (t ) (0) (0) 4 SI (k0 ) (t ) + S jk j (t ) + 2Fv(,0k),( 0) I k (t ) + 2 j, k = 0, (2.31) t t j =1 j (0) 4 S (k0 ) I (k0) (t ) S f k j (t ) + 6Fv(,0k),(1) I k (t ) = 0. (2.32) t j = Для l = C (t ) (1) (1) 2 SI (t ) + S cos( j ) jk j + 2F I k (t ) + 2 j, k = 0, (1) (1),( 0 ) (2.33) t t k v,k j =1 j (1) 2S (k1) I (k1) (t ) S f k j cos( j ) + 6Fv(,1k),(1) I k (t ) = 0. (2.34) t j = Для l = C (t ) ( 2) ( 2) 2 SI (k2 ) (t ) + S jk j sin( j ) + 2Fv(,2k),( 0 ) I k (t ) + 2 j, k =0, (2.35) t t j =1 j ( 2) 2S (k2) I (k2 ) (t ) S f k j sin( j ) + 6Fv(,2k),(1) I k (t ) = 0. (2.36) t j = Для l = C (t ) ( 3) ( 3) 4 SI (k3) (t ) + S jk j cos(2 j ) + 2Fv(,3k),( 0) I k (t ) + 2 j, k =0, (2.37) t t j =1 j ( 3) 4S (k3) I (k3) (t ) S f k j cos(2 j ) + 6Fv(,3k),(1) I k (t ) = 0. (2.38) t j = В уравнениях (2.31) – (2.38) использованы следующие обозначения Fv(,lk), ( n ) = ( k,тn( l ) ) + + т(l ) (r ) v inv (kl,)n (r )dr, (2.39) k,n Vk C(jl,)k (t ) = FСl,)j, k d, j, k, ( (2.40) где FСl,)j, k = ( k,т0( l ) ) + + т(l ) ( (r )C j, k (r, t )dr. (2.41) k, Vk Умножим уравнение (2.31) на (k1) и вычтем его из уравнения (2.32), получим ~ ( 0) ( ) S f k, j (t ) + (k1) jk j (t ) = 4SФ k (t ) 2 (k1) C j,k (t ) Fv(,0k) Ф k (t ), (2.42) t t j =1 j где ( )( ) ~ Fv(,0k) = 2 (k1) Fv(,0k),( 0 ) 3Fv(,0k),(1) (k0) (k1). (2.43) Умножим уравнение (2.34) на (k1) и вычтем его из уравнения (2.35), получим ~ (1) ( ) S f k, j + (k1) j k j cos( j ) = Fv(1k) X k (t ) 2 (k1) C (j1,)k (t ), (2.44) t j t, j = где ( )( ) ~ Fv(,1k) = 2 (k1) Fv(,1k),( 0) 3Fv(,1k),(1) (k0 ) (k1). (2.45) Умножим уравнение (2.35) на (k1), вычтем его из уравнения (2.36) и сделаем предположение о симметрии ячеек относительно поворота на угол 90°, в результате чего (k2) = (k1). Тогда ~ ( 2) ( ) S f k, j + (k1) j k j sin( j ) = Fv(,2k) X k (t ) 2 (k1) C (j2,k) (t ), (2.46) t j t j = где ( )( ) ~ Fv(,2k) = 2 (k1) Fv(,2),( 0) 3Fv(,2),(1) (k0) (k1). (2.47) k k Умножим уравнение (2.37) на (k3) и вычтем его из уравнения (2.38), получим ~ ( 3) ( ) S f k, j + (k1) j k j cos(2 j ) = 4 SX (k3) (t ) Fv(,3k) X k (t ) 2 (k1) C (j3, k (t ), ) (2.48) t j t j = где ( )( ) ~ Fv(,3k) = 2 (k1) Fv(,3k),( 0 ) 3Fv(,3k),(1) (k1) (k3). (2.49) В формулах (2.42), (2.44), (2.46) и (2.48) использованы следующие обозначения ( ) Ф k = (k0 ) (k1) I (k0 ) (1) ( ) X k = k k I k (0 ) (1) (1) (2.50) ( 2) ( ) X k = k k I k ( 0) (1) ( 2 ) ( ) X ( 3) = (1) ( 3) I ( 3) k k k k Запишем средние по грани ток jkj и уровень (поток) f kj jk, j =I (k0 ) + I (k1) cos j + I (k2 ) sin j + I (k3) cos(2 j ), (2.51) f k, j = (k0) I (k0) + (k1) I (k1) cos j + (k2 )I (k2) sin j + (k3) I (k3) cos(2 j ), (2.52) и для граничной ячейки jk, j =I (k0 ) (I (k1) cos j + I (k2 ) sin j ) + I (k3 ) cos(2 j ), (2.53) f k, j = (k0 ) I k0) (k1) ( I (k1) cos j + I (k2 ) sin j ) + (k3) I (k3 ) cos(2 j ).



Pages:   || 2 | 3 | 4 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.