Участник:S s serov/DSn-метод решения уравнения переноса

Автор текущей версии статьи: Серов Сергей

Свойства и структура алгоритма

Общее описание алгоритма

Уравнение переноса нейтронов описывает процессы поглощения, рассеяния и размножения, происходящие в ядерном реакторе. Одним из популярных детерминистических методов его решения, как известно, является DSn-метод . Далее описывается DSn-метод решения стационарного одногруппового уравнения переноса нейтронов в 1D плоско-параллельной геометрии с изотропным рассеянием.

Геометрическая модель решаемой задачи имеет следующий вид: перпендикулярно единственной участвующей в рассмотрении оси [math]\displaystyle{ x }[/math] (будем считать, что она направлена слева направо) расположено несколько бесконечных по осям [math]\displaystyle{ y }[/math] и [math]\displaystyle{ z }[/math] плоско-параллельных слоев некоторого вещества (например, урана, свинца, воды, воздуха и т.д.) заданной толщины, прилегающих друг к другу.

Само уравнение переноса (УП) нейтронов записывается в следующем виде:

[math]\displaystyle{ \mu \frac{dN(x,\mu)}{dx} + \sigma(x)N(x,\mu) = Q(x,\mu), }[/math]

где

[math]\displaystyle{ Q(x,\mu) = \frac{\sigma_S^0(x)}{2}\int_{-1}^1 N(x,\mu')d\mu' + \frac{1}{2}q(x). }[/math]

Используемые обозначения:

• [math]\displaystyle{ \mu = cos(\theta) }[/math] --- косинус угла [math]\displaystyle{ \theta }[/math] между направлением полета нейтрона и положительным направлением оси [math]\displaystyle{ x }[/math];
• [math]\displaystyle{ N(x,\mu) }[/math] --- плотность потока нейтронов, летящих в заданном направлении [math]\displaystyle{ \mu }[/math] в слое с координатой [math]\displaystyle{ x }[/math];
• [math]\displaystyle{ SN(x) = \int_{-1}^1 N(x,\mu')d\mu' }[/math] --- скалярный поток нейтронов в точке [math]\displaystyle{ x }[/math];
• [math]\displaystyle{ \sigma(x) }[/math] --- макроскопическое сечение взаимодействия нейтронов со средой в точке [math]\displaystyle{ x }[/math], т.е. коэффициент обобщенного поглощения нейтронов физическим материалом, расположенным в этой точке среды;
• [math]\displaystyle{ \sigma_S^0(x) }[/math] --- изотропное макроскопическое сечение обобщенного рассеяния нейтронов в точке [math]\displaystyle{ x }[/math] среды, т.е. коэффициент обобщенного рассеяния нейтронов в этой точке среды;
• [math]\displaystyle{ q(x) }[/math] --- независимый изотропный источник нейтронов в точке [math]\displaystyle{ x }[/math].

Одногрупповое УП является интегро-дифференциальным. Левая его часть представляет собой обыкновенный дифференциальный оператор, который отвечает за пролет нейтронов через физическую среду. Правая часть УП содержит интеграл столкновений, который отвечает за рождение вторичных частиц в реакциях упругого и неупругого рассеяния, деления нейтронов, реакций типа [math]\displaystyle{ (n,2n) }[/math] и других, а также независимый источник нейтронов. Целью решаемой задачи является нахождение плотности потока нейтронов [math]\displaystyle{ N(x,\mu) }[/math] или скалярного потока [math]\displaystyle{ SN(x) }[/math]. Для решения задачи достаточно вычислить только скалярный поток [math]\displaystyle{ SN(x) }[/math], так как дифференциальная плотность потока нейтронов [math]\displaystyle{ N(x,\mu) }[/math] может быть получена из решения ОДУ с известной правой частью, в которую входит скалярный поток [math]\displaystyle{ SN(x) }[/math]. Таким образом, целью расчета является нахождение стационарного (на бесконечности по времени) распределения плотности потока нейтронов [math]\displaystyle{ N(x,\mu) }[/math] и скалярного потока [math]\displaystyle{ SN(x) }[/math].

Для каждого слоя вещества необходимый набор параметров ([math]\displaystyle{ \sigma }[/math] и [math]\displaystyle{ \sigma_S^0 }[/math]), задающих его физические свойства, известен, а также известны краевые условия на левом и правом концах рассматриваемой системы. В рассматриваемой постановке задачи мы считаем, что правее [math]\displaystyle{ x_{max} }[/math] находится вакуум, то есть плотность потока нейтронов равна нулю по всем направлениям полета нейтронов, влетающих справа налево через правую границу системы: [math]\displaystyle{ N(x_{max},\mu) = 0, \forall \mu \in [-1, 0] }[/math]; а на левой границе [math]\displaystyle{ x_{min} }[/math] заданы краевые условия зеркального отражения, при котором нейтроны, прилетевшие справа на левую границу, отражаются от нее, изменяя направление своего полета с [math]\displaystyle{ \mu }[/math] на [math]\displaystyle{ -\mu }[/math], то есть краевое условие имеет вид: [math]\displaystyle{ N(x_{min},-\mu) = N(x_{min},\mu), \forall \mu \in [-1, 0] }[/math]. Также значения параметров [math]\displaystyle{ \sigma(x), \sigma_S^0(x) }[/math] и [math]\displaystyle{ q(x) }[/math] будем считать постоянными в каждом из слоев. Слой с одинаковыми макросечениями и независимым источником нейтронов будем также называть областью.

Для численного решения задачи вводятся 2 одномерные дискретные сетки: по оси [math]\displaystyle{ x }[/math] и по косинусу угла [math]\displaystyle{ \mu }[/math].

Поэтому данный метод решения УП нейтронов и получил название метода дискретных ординат (по [math]\displaystyle{ \mu }[/math] на плоскости [math]\displaystyle{ (x,\mu) }[/math]), или DSn – метода.

Альтернативными методами решения этого уравнения являются метод характеристик, метод Монте-Карло и некоторые другие.

Математическое описание алгоритма

Входные данные:

• геометрическая модель задачи (количество слоев, их расположение на оси [math]\displaystyle{ x }[/math] и толщина);
• [math]\displaystyle{ \sigma(x) }[/math] для каждого слоя;
• [math]\displaystyle{ \sigma_S^0(x) }[/math] для каждого слоя;
• [math]\displaystyle{ q(x) }[/math] для каждого слоя.

Параметры:

• [math]\displaystyle{ I }[/math] --- количество ячеек дискретной сетки по оси [math]\displaystyle{ x }[/math];
• [math]\displaystyle{ M }[/math] --- количество ячеек дискретной сетки по параметру [math]\displaystyle{ \mu }[/math];
• [math]\displaystyle{ \nu_0 }[/math] --- количество итераций по правой части УП.

Выходные данные:

• [math]\displaystyle{ SN(x) }[/math] для каждой ячейки введенной сетки.

Введем дискретную сетку по [math]\displaystyle{ x }[/math] и по [math]\displaystyle{ \mu }[/math]:

[math]\displaystyle{ \mathbb{X} = \{x_i, i \in \{1, \dots, I + 1\} | x_{i+1} = x_i + \Delta x, x_1 = x_{min}, \Delta x = \frac{x_{max} - x_{min}}{I}\}, }[/math]

[math]\displaystyle{ \mathbb{M} = \{\mu_m, m \in \{1, \dots, M + 1\} | \mu_{m+1} = \mu_m + \Delta \mu, \mu_1 = -1, \Delta \mu = \frac{2}{M}\}. }[/math]

Середины ячеек соответствующих сеток будем обозначать как [math]\displaystyle{ x_{i+1/2} }[/math] и [math]\displaystyle{ \mu_{m+1/2} }[/math].

Алгоритм решения задачи выглядит так:

1. Задается некоторое начальное неотрицательное распределение нейтронов в системе [math]\displaystyle{ SN_{i+1/2} }[/math].
2. Для решения интегро-дифференциального УП нейтронов используем метод итераций источника, то есть организуем итерации [math]\displaystyle{ \nu }[/math] по правой части УП, которая называется источником (он является суммой интеграла столкновений и независимого источника нейтронов). Независимый источник [math]\displaystyle{ q_{i+1/2} }[/math] от номера итерации не зависит и является постоянным на всех итерациях. Интеграл столкновений [math]\displaystyle{ SN_{i+1/2} }[/math] нужно перевычислять на каждой итерации. Более конкретно, по уже известному с предыдущей итерации или из краевого условия распределению [math]\displaystyle{ SN_{i+1/2} }[/math] вычисляем правую часть УП, то есть [math]\displaystyle{ Q_{i+1/2} }[/math] для всех [math]\displaystyle{ i \in \{1, \dots, I\} }[/math]. Далее задача сводится к решению ОДУ на каждой итерации по правой части УП.
3. Алгоритм решения ОДУ на каждой итерации по правой части является неоднородным. Сначала рассчитываются все направления полета нейтронов с отрицательными [math]\displaystyle{ \mu }[/math]. При этом на правой границе системы задаются вакуумные краевые условия, то есть [math]\displaystyle{ N(x_{max},\mu)=0 }[/math] для всех [math]\displaystyle{ \mu \lt 0 }[/math]. Расчет по [math]\displaystyle{ x }[/math] ведется справа налево, то есть по убыванию индекса [math]\displaystyle{ i+1/2 }[/math]. На левой границе системы [math]\displaystyle{ x_{min} }[/math] формируются краевые условия отражения [math]\displaystyle{ N_{Left}(\mu_{m+1/2}) }[/math] (для всех [math]\displaystyle{ \mu \lt 0 }[/math]). Далее они будут использоваться как краевые условия на левой границе системы [math]\displaystyle{ x_{min} }[/math] при расчете всех направлений [math]\displaystyle{ \mu \gt 0 }[/math].
4. Рассчитываем все направления [math]\displaystyle{ \mu_{m+1/2} \gt 0 }[/math]. При этом порядок расчета по координате [math]\displaystyle{ x }[/math] будет следующим: от левой границы системы [math]\displaystyle{ x_{min} }[/math] к правой [math]\displaystyle{ x_{max} }[/math], то есть по возрастанию индекса [math]\displaystyle{ i+1/2 }[/math]. Чтобы "оттолкнуться" от левой границы [math]\displaystyle{ x_{min} }[/math] используем уже сформированные при расчете случая [math]\displaystyle{ \mu \lt 0 }[/math] краевые условия [math]\displaystyle{ N_{Left}(\mu_{m+1/2}) }[/math] на левой границе системы. Продолжаем для случая [math]\displaystyle{ \mu \gt 0 }[/math] накапливать массив [math]\displaystyle{ SN_{i+1/2}^{new} }[/math]. В конце расчета всех [math]\displaystyle{ \mu_{m+1/2} }[/math] в массиве [math]\displaystyle{ SN_{i+1/2}^{new} }[/math] уже насчитаны новые значения [math]\displaystyle{ SN }[/math] для следующей итерации [math]\displaystyle{ \nu }[/math]. Их нужно только присвоить в нужном месте программы снова в массив [math]\displaystyle{ SN_{i+1/2} }[/math].
5. Теперь опишем метод решения обыкновенного дифференциального уравнения с известной правой частью и используемые для этого разностные схемы. Также учитываем, что в указанной постановке задачи параметры [math]\displaystyle{ \sigma(x), \sigma_S^0(x) }[/math] и [math]\displaystyle{ q(x) }[/math] заданы постоянными для каждого слоя.

   Уравнение имеет вид:

   :[math]\displaystyle{ \mu_{m+1/2} \frac{dN(x,\mu_{m+1/2})}{dx} + \sigma(x)N(x,\mu_{m+1/2}) = Q(x,\mu_{m+1/2}). }[/math]

   Проинтегрируем его по ячейке [math]\displaystyle{ T_{i+1/2} = [x_i, x_{i+1}] }[/math] и получим уравнение баланса нейтронов (сечение взаимодействия нейтронов со средой [math]\displaystyle{ \sigma(x) = Const }[/math] в каждой области вещества и может быть вынесено за знак интеграла):

   :[math]\displaystyle{ \mu_{m+1/2}(N_{i+1} - N_i) + \sigma(x) \int_{x_i}^{x_{i+1}} N(x,\mu_{m+1/2})dx = \int_{x_i}^{x_{i+1}} Q(x,\mu_{m+1/2})dx. }[/math]

   Далее по теореме о среднем:

   :[math]\displaystyle{ \int_{x_i}^{x_{i+1}} N(x,\mu_{m+1/2})dx = \bar{N}\cdot \Delta x, \int_{x_i}^{x_{i+1}} Q(x,\mu_{m+1/2})dx = \bar{Q}\cdot \Delta x, \Delta x = x_{i+1} - x_i. }[/math]

   Тогда уравнение (называемое уравнением баланса) примет вид:

   :[math]\displaystyle{ \mu_{m+1/2}(N_{i+1} - N_i) + \sigma \bar{N} \Delta x = \bar{Q} \Delta x. }[/math]

6. При решении УП нейтронов на основе традиционного DSn-метода обычно используются 2 основные конечно-разностные схемы: DD (diamond difference) – алмазная и St (step) – шаговая. Приведем далее формулы, соответствующие случаю расчета всех направлений полета нейтронов с положительными [math]\displaystyle{ \mu }[/math]. Напомним, что в этом случае расчет по координате [math]\displaystyle{ x }[/math] ведется по возрастанию индекса [math]\displaystyle{ i+1/2 }[/math]. Для случая [math]\displaystyle{ \mu \lt 0 }[/math] формулы для DD– и St– схем выводятся аналогично.

   St-схема (шаговая) ([math]\displaystyle{ \bar{N} = N_{i+1}) }[/math]:

   :[math]\displaystyle{ \bar{N} = \frac{\bar{Q} \Delta x + \mu_{m+1/2} N_i}{\sigma \Delta x + \mu_{m+1/2}}, N_{i+1} = \frac{\bar{Q} \Delta x + \mu_{m+1/2} N_i}{\sigma \Delta x + \mu_{m+1/2}}. }[/math]

   DD-схема (алмазная) ([math]\displaystyle{ \bar{N} = \frac{1}{2}(N_i + N_{i+1}), N_{i+1} = 2\bar{N} - N_i }[/math]):

   :[math]\displaystyle{ \bar{N} = \frac{\bar{Q} \Delta x + 2\mu_{m+1/2} N_i}{\sigma \Delta x + 2\mu_{m+1/2}}, N_{i+1} = \frac{2\bar{Q} \Delta x + (2\mu_{m+1/2} - \sigma \Delta x) N_i}{\sigma \Delta x + 2\mu_{m+1/2}}. }[/math]

   Алмазная DD-схема имеет 2-ой порядок аппроксимации и поэтому является более точной, чем St-схема 1-го порядка аппроксимации, однако не обладает свойствами положительности и монотонности, в отличие от последней, вследствие чего в некоторых случаях может давать отрицательные значения плотности потока, что недопустимо с точки зрения физики явления (невозможно в реальной физической ситуации). В связи с этим традиционный DSn–метод использует обе конечно-разностные схемы для обеспечения, с одной стороны, максимальной точности разностного решения, а с другой – гарантированной его неотрицательности, то есть физичности.

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

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

Вычислительное ядро алгоритма

Вычислительное ядро алгоритма состоит из вычисления значений плотности потока нейтронов в каждой ячейке сетки по оси [math]\displaystyle{ x }[/math] для каждого значения косинуса направления полета нейтронов в соответствии с выбранными разностными схемами.

Макроструктура алгоритма

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

Схема реализации последовательного алгоритма

1. Задать неотрицательное распределение скалярного потока нейтронов [math]\displaystyle{ SN }[/math] в системе.
2. На каждой итерации:
   1. Вычислить правую часть УП для каждой ячейки двумерной сетки по [math]\displaystyle{ x }[/math] и [math]\displaystyle{ \mu }[/math]: [math]\displaystyle{ Q(x,\mu) = \frac{\sigma_S^0(x)}{2}\int_{-1}^1 N(x,\mu')d\mu' + \frac{1}{2}q(x) }[/math].
   2. Провести вычисление плотностей потока нейтронов в каждой ячейке двумерной сетки справа налево по оси [math]\displaystyle{ x }[/math] с использованием алмазной и шаговой (при необходимости) разностных схем:

      St-схема (шаговая) ([math]\displaystyle{ \bar{N} = N_{i+1}) }[/math]:

      :[math]\displaystyle{ \bar{N} = \frac{\bar{Q} \Delta x + \mu_{m+1/2} N_i}{\sigma \Delta x + \mu_{m+1/2}}, N_{i+1} = \frac{\bar{Q} \Delta x + \mu_{m+1/2} N_i}{\sigma \Delta x + \mu_{m+1/2}}. }[/math]

      DD-схема (алмазная) ([math]\displaystyle{ \bar{N} = \frac{1}{2}(N_i + N_{i+1}), N_{i+1} = 2\bar{N} - N_i }[/math]):

      :[math]\displaystyle{ \bar{N} = \frac{\bar{Q} \Delta x + 2\mu_{m+1/2} N_i}{\sigma \Delta x + 2\mu_{m+1/2}}, N_{i+1} = \frac{2\bar{Q} \Delta x + (2\mu_{m+1/2} - \sigma \Delta x) N_i}{\sigma \Delta x + 2\mu_{m+1/2}}. }[/math]

   3. Вычислить скалярный поток нейтронов в каждой ячейке сетки по [math]\displaystyle{ x }[/math] путем интегрирования плотностей потока нейтронов по всем значениям косинуса угла полета [math]\displaystyle{ \mu }[/math]: [math]\displaystyle{ SN(x) = \int_{-1}^1 N(x,\mu')d\mu' }[/math].

Входные и выходные данные алгоритма

Входные данные: массив [math]\displaystyle{ \sigma(x) }[/math] для каждого слоя, массив [math]\displaystyle{ \sigma_S^0(x) }[/math] для каждого слоя, массив [math]\displaystyle{ q(x) }[/math] для каждого слоя. Объем входных данных: [math]\displaystyle{ 3k }[/math], где [math]\displaystyle{ k }[/math] --- количество слоев в геометрии задачи.

Выходные данные: массив [math]\displaystyle{ SN(x) }[/math] для каждой ячейки введенной сетки. Объем выходных данных: [math]\displaystyle{ N }[/math], где [math]\displaystyle{ N }[/math] --- количество ячеек введенной сетки по [math]\displaystyle{ x }[/math].

Литература

<references \>