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

Автор текущей версии статьи: Серов Сергей, гр. 417, 2017-ый год

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

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

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

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

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

где

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

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

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

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

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

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

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

Входные данные: \begin{itemize}

\item геометрическая модель задачи (количество слоев, их расположение на оси [math]x[/math] и толщина);

\item [math]\sigma(x)[/math] для каждого слоя;

\item [math]\sigma_S^0(x)[/math] для каждого слоя;

\item [math]q(x)[/math] для каждого слоя.

\end{itemize}

Параметры: \begin{itemize}

\item [math]IR[/math] --- количество ячеек дискретной сетки по оси [math]x[/math];

\item [math]M_0[/math] --- количество ячеек дискретной сетки по параметру [math]\mu[/math].

\end{itemize}

Выходные данные: \begin{itemize}

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

\end{itemize}

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

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

[math]\mathbb{M} = \{\mu_m, m \in \{1, \dots, M_0 + 1\} | \mu_{m+1} = \mu_m + \Delta_{\mu}, \mu_1 = -1, \Delta_{\mu} = \frac{1 - (-1)}{M_0}\}.[/math]

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

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

\begin{enumerate}

\item Рассмотреть значения [math]\mu \lt 0[/math]. Для каждого [math]m \in \{1, \dots, M_0 + 1\}:\mu_m \lt 0[/math]:

\begin{enumerate}

\item Использовать метод итерации источника. Учесть краевые условия в точке [math]x = x_{IR+1}[/math]. Для каждого [math]i \in \{IR + 1, \dots, 2\}[/math] последовательно от [math]IR+1[/math] в порядке убывания:

\begin{enumerate}

\item Рассчитать правую часть уравнения переноса по формуле: :[math]Q(x_i,\mu_m) = \frac{\sigma_S^0(x_i)}{2}SN(x_i) + \frac{1}{2}q(x_i),[/math] где :[math]SN(x_i) = \int_{-1}^1 N(x_i,\mu_m)d\mu_m = \sum_{m=1}^{\frac{M_0}{2}} N(x_i,\mu_m) \cdot \Delta_{\mu} = \Delta_{\mu} \cdot \sum_{m=1}^{\frac{M_0}{2}} N(x_i,\mu_m).[/math]

\item Решить обыкновенное дифференциальное уравнение с известной правой частью относительно [math]N(x,\mu_m)[/math] и вычислить [math]N(x_{i-1},\mu_m)[/math], [math]\bar{N}(x,\mu_m) = N(x_{i-1/2},\mu_m)[/math] и [math]SN_{Left}(x_{i-1/2})[/math] по выбранной разностной схеме. При необходимости использовать правило переключения схем.

\end{enumerate}

\end{enumerate}

\item Рассмотреть значения [math]\mu \geq 0[/math]. Для каждого [math]m \in \{1, \dots, M_0 + 1\}:\mu_m \geq 0[/math]:

\begin{enumerate}

\item Использовать метод итерации источника. Учесть полученные на предыдущем шаге значения [math]N(x_1,\mu_m), \mu_m \lt 0,[/math] при вычислении краевых условий в точке [math]x = x_{1}[/math]. Для каждого [math]i \in \{1, \dots, IR\}[/math] последовательно от [math]1[/math] в порядке возрастания:

\begin{enumerate}

\item Рассчитать правую часть уравнения переноса по формуле: :[math]Q(x_i,\mu_m) = \frac{\sigma_S^0(x_i)}{2}SN(x_i) + \frac{1}{2}q(x_i),[/math] где :[math]SN(x_i) = \int_{-1}^1 N(x_i,\mu_m)d\mu_m = \sum_{m=1}^{\frac{M_0}{2}} N(x_i,\mu_m) \cdot \Delta_{\mu} = \Delta_{\mu} \cdot \sum_{m=1}^{\frac{M_0}{2}} N(x_i,\mu_m).[/math]

\item Решить обыкновенное дифференциальное уравнение с известной правой частью относительно [math]N(x,\mu_m)[/math] и вычислить [math]N(x_{i+1},\mu_m)[/math], [math]\bar{N}(x,\mu_m) = N(x_{i+1/2},\mu_m)[/math] и [math]SN_{Right}(x_{i+1/2})[/math] по выбранной разностной схеме. При необходимости использовать правило переключения схем.

\end{enumerate}

\end{enumerate}

\item В качестве ответа вычислить [math]SN(x_{i+1/2}) = SN_{Left}(x_{i+1/2}) + SN_{Right}(x_{i+1/2})[/math] для каждого [math]i \in \{1, \dots, IR\}[/math].

\end{enumerate}

Теперь опишем метод решения обыкновенного дифференциального уравнения с известной правой частью и используемые для этого разностные схемы. Также учитываем, что в указанной постановке задачи параметры [math]\sigma(x), \sigma_S^0(x)[/math] и [math]q(x)[/math] заданы константными для каждого слоя.

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

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

Проинтегрируем его по отрезку [math][x_i, x_{i+1}][/math]:

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

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

[math]\int_{x_i}^{x_{i+1}} N(x,\mu_m)dx = \bar{N}\cdot \Delta_x, \int_{x_i}^{x_{i+1}} Q(x,\mu_m)dx = \bar{Q}\cdot \Delta_x, \Delta_x = x_{i+1} - x_i.:[/math]

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

[math]\mu_m(N_{i+1} - N_i) + \sigma \bar{N} \Delta_x = \bar{Q} \Delta_x.[/math]

Далее приведем формулы для шага 2 алгоритма (то есть для [math]\mu_m \geq 0[/math]). Найдем решения этого разностного уравнения по двум разностным схемам. \\ St-схема (шаговая) ([math]\bar{N} = N_{i+1})[/math]:

[math]\bar{N} = \frac{\bar{Q} \Delta_x + \mu_m N_i}{\sigma \Delta_x + \mu_m}, N_{i+1} = \frac{\bar{Q} \Delta_x + \mu_m N_i}{\sigma \Delta_x + \mu_m}.[/math]

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

[math]\bar{N} = \frac{\bar{Q} \Delta_x + 2\mu_m N_i}{\sigma \Delta_x + 2\mu_m}, N_{i+1} = \frac{2\bar{Q} \Delta_x + (2\mu_m - \sigma \Delta_x) N_i}{\sigma \Delta_x + 2\mu_m}.[/math]

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

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