Журнал вычислительной математики и математической физики, 2022, T. 62, № 7, стр. 1200-1208

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

Н. Я. Моисеев^1, *, В. М. Шмаков^1, **

¹ ФГУП РФЯЦ-ВНИИТФ им. академ. Е.И. Забабахина
456770 Снежинск, Челябинская обл., ул. Васильева, 13, а/я 245, Россия

^* E-mail: nik.moiseev.43@mail.ru
^** E-mail: v.m.shmakov@vniitf.ru

Поступила в редакцию 11.03.2021
После доработки 31.08.2021
Принята к публикации 11.02.2022

EDN: VVRUZS
DOI: 10.31857/S0044466922070079

Аннотация

Представлена дискретно-аналитическая разностная схема для решения нестационарного кинетического уравнения переноса частиц (нейтронов) в многогрупповом изотропном приближении методом расщепления. Особенность разностной схемы состоит в том, что решение уравнения переноса в многогрупповой модели сведено к решениям уравнений в одногрупповой модели. Эффективность схемы достигнута за счет вычисления интеграла столкновений из аналитических решений обыкновенных дифференциальных уравнений, которые описывают эволюцию нейтронов, пришедших в группу $g$ из всех групп $g{\kern 1pt} '$. Решения уравнений находятся без итераций по интегралу столкновений и без обращения матриц. Метод решения естественным образом обобщается на решение задач в многомерных пространствах и позволяет осуществить счет в параллельном режиме. Библ. 18. Фиг. 2.

Ключевые слова: кинетическое уравнение переноса нейтронов, интеграл столкновений, метод расщепления, аналитические решения.

1. ВВЕДЕНИЕ

Кинетическое уравнение переноса частиц (в частности, нейтронов) является интегро-дифференциальным уравнением “жесткого” типа. В общем случае правая часть уравнения содержит интеграл столкновений от искомой (неизвестной) функции. Одним из подходов к решению таких уравнений является подход, в котором решения находятся с применением методов дискретных ординат $Sn,\;Dsn$ (см. [1], [2]) по неявным разностным схемам итерациями по интегралу столкновений (см. [3]–[5]). Явные численные методы требуют неприемлемых для реальных расчетов ограничений на выбор шага интегрирования по времени. Основные трудности, которые возникают при решении уравнений итерационными методами по неявным схемам, связаны с точностью численных решений и медленной сходимостью итераций. Если систему уравнений решать методами Фотрие (см. [6]) или Райбики (см. [7]), то приходится иметь дело с решением систем линейных алгебраических уравнений с полными матрицами взаимодействия между группами (см. [8]). Как следствие, в случае решения реальных задач требуются большие объемы оперативной памяти и большие временные затраты.

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

Здесь представлена дискретно-аналитическая разностная схема (ДАРС) для решения интегро-дифференциального уравнения переноса нейтронов в многогрупповом изотропном приближении на основе модифицированного метода расщепления (см. [11]). Особенность разностной схемы состоит в том, что решение уравнения переноса в многогрупповом приближении сводится к решениям уравнений в одногрупповом приближении в каждой группе. Интеграл столкновений вычисляется из аналитических решений обыкновенных дифференциальных уравнений (ОДУ), которые описывают эволюцию нейтронов, пришедших в группу $g$ из всех групп $g{\kern 1pt} '$.

Представлены результаты расчетов модельной задачи из работы [12]. Задача имеет точные решения. Показано, что результаты расчетов по методике ДАРС удовлетворительно согласуются с точными решениями и с результатами расчетов по одномерному аналогу методики из [13]. Время решения модельной задачи по методике ДАРС по сравнению со временем решения по одномерному аналогу методики из [13] сократилось в 1.8 раза.

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

2.1. Постановка задачи

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

(1)

$\begin{gathered} \frac{1}{{{{{v}}_{g}}}}\frac{{\partial {{N}_{g}}}}{{\partial t}} + \mu \frac{{\partial {{N}_{g}}}}{{\partial x}} + {{\alpha }_{g}}{{N}_{g}} = 0.5{\kern 1pt} \mathop \sum \limits_{g' = 1}^G \,{{\beta }_{{g',g}}}{{S}_{{g'}}} + 0.5{{f}_{g}}, \\ {{S}_{g}} = \mathop \smallint \limits_{ - 1}^1 {{N}_{g}}\left( \mu \right)d\mu . \\ \end{gathered} $

Здесь $t,\;x$ – независимые переменные по времени и по пространству, $\mu $ – косинус угла между направлением полета нейтрона и осью $x$, $ - 1 \leqslant \mu \leqslant 1,$ $g$ – индекс группы, $g = 1,~2,...,G$, $G$ – число энергетических групп, ${{{v}}_{g}}$ – модуль вектора скорости нейтронов группы $g$, ${{N}_{g}}\left( {t,x,\mu } \right)$ – плотность потока нейтронов, ${{S}_{g}}\left( {t,x} \right)$ – полный поток нейтронов, ${{\alpha }_{g}}$ – коэффициент поглощения нейтронов группы $g$, ${{\beta }_{{g',g}}}$ – коэффициент размножения нейтронов с учетом переходов нейтронов между группами $g$ и $g{\kern 1pt} '$, ${{f}_{g}}$ – независимый источник нейтронов. Разностная сетка по переменной $\mu $ с центрами ${{\mu }_{m}} = 0.5\left( {{{\mu }_{{m + 1/2}}} + {{\mu }_{{m - 1/2}}}} \right)$ и шагами $\Delta \mu = {{\mu }_{{m + 1/2}}} - {{\mu }_{{m - 1/2}}} = 2{\text{/}}M$, $m = 1, \ldots ,M$, включает $M$ направлений движения частиц. Решается смешанная краевая задача с начальными и граничными условиями в области $D = \left\{ {{{x}_{L}} \leqslant x \leqslant {{x}_{R}},~~ - 1 \leqslant \mu \leqslant 1} \right\}$. Начальные условия: ${{N}_{g}}\left( {0,x,\mu } \right) = N_{g}^{0}\left( {x,\mu } \right)$. Граничные условия: на левой границе – ${{N}_{g}}\left( {t,{{x}_{L}},\mu } \right) = {{\varphi }_{L}}\left( {t,{{x}_{L}},\mu } \right)$ для $\mu > 0$, на правой – ${{N}_{g}}\left( {t,{{x}_{R}},\mu } \right) = {{\varphi }_{R}}\left( {t,{{x}_{R}},\mu } \right)$ для $\mu < 0$, где ${{\varphi }_{L}}\left( {t,{{x}_{L}},\mu } \right)$, ${{\varphi }_{R}}\left( {t,{{x}_{R}},\mu } \right)$ – известные функции. Предполагаем, что плотность вещества $\rho = 1$ . Требуется найти решение системы уравнений (1) для $t > 0$.

2.2. Полунеявная разностная схема

В пространстве $\left( {t,x} \right)$ построим разностную сетку с шагами интегрирования $\tau ,\;{{h}_{j}}$ вдоль координатных линий $t,\;x$ соответственно. Границы ячеек разностной сетки с центрами в точках ${{x}_{j}}$ обозначим как ${{x}_{{j - 1/2}}},~\;{{x}_{{j + 1/2}}}$, где $j$ – индекс ячейки. Введем обозначения: $N_{{g,j}}^{{n + 1}},\;N_{{g,j}}^{n}$ – плотности потока нейтронов, $S_{{g,j}}^{{n + 1}},\;S_{{g,j}}^{n}~$ – полные потоки нейтронов, $\alpha _{{g,j}}^{{n + 1}},\;\alpha _{{g,j}}^{n}~$ – коэффициенты поглощения, $\beta _{{g',g,j}}^{{n + 1}},\;\beta _{{g',g,j}}^{n}$ – коэффициенты размножения нейтронов, $f_{{g,j}}^{{n + 1}},\;f_{{g,j}}^{n}$ – источники. Величины с индексами $n + 1$, $n$ – это основные величины, которые относятся к центрам ячеек в моменты времени ${{t}^{{n + 1}}} = {{t}^{n}} + \tau $ и ${{t}^{n}}$ соответственно. Величины с дробными индексами ${{N}_{{g,j + 1/2}}},\;{{N}_{{g,j - 1/2}}}$, которые будем называть “большими”, как в схемах С.К. Годунова (см. [14]), относятся к граням ячеек и являются вспомогательными. В дальнейшем индекс $j$ у величин в центрах ячеек опущен.

Заменив интеграл в (1) квадратурной формулой, проинтегрировав первое уравнение в (1) по ячейке фазового пространства $\left( {t,x} \right)$ на интервале по времени $\left[ {{{t}^{n}},{{t}^{{n + 1}}}} \right]$ и применив теорему Гаусса–Остроградского, получим систему разностных уравнений, которые запишем в виде

(2)

$\begin{gathered} N_{g}^{{n + 1}} = \left\{ {N_{g}^{n} - \frac{{{{{v}}_{g}}\tau \mu }}{h}\left( {{{N}_{{g,j + 1/2}}} - {{N}_{{g,j - 1/2}}}} \right)} \right\} + {{{v}}_{g}}\tau \left[ { - \alpha _{g}^{*}N_{g}^{*} + 0.5~\mathop \sum \limits_{g' = 1}^G \,\beta _{{g',g}}^{*}S_{{g'}}^{*} + 0.5f_{g}^{*}} \right], \\ S_{g}^{{n + 1}} = \mathop \sum \limits_{m = 1}^M \,N_{{g,m}}^{{n + 1}}\Delta \mu . \\ \end{gathered} $

Здесь величины $\alpha _{g}^{*},\;\beta _{{g',g}}^{*},\;f_{g}^{*},\;N_{g}^{*},S_{{g'}}^{*}$ – это средние значения функций на интервале ${{t}^{n}},\;{{t}^{{n + 1}}}$, которые вычислены в момент времени ${{t}^{n}} + \tau {\kern 1pt} *$, где $0 \leqslant \tau {\kern 1pt} * \leqslant \tau $. Введем вспомогательные потоки нейтронов $N_{g}^{{n + 1/2}}$, которые вычисляются из выражения в фигурных скобках в (2) по разностным уравнениям

(3)

$N_{g}^{{n + 1/2}} = N_{g}^{n} - \frac{{{{{v}}_{g}}\tau \mu }}{{{{h}_{j}}}}\left( {{{N}_{{g,j + 1/2}}} - {{N}_{{g,j + 1/2}}}} \right).~$

Заменив выражение в фигурных скобках на вспомогательный поток $N_{g}^{{n + 1/2}}$ в первом уравнении в (2) и подставив выражение для $N_{g}^{{n + 1}}$ в квадратурную формулу для $S_{g}^{{n + 1}}$ в (2), получим для вычисления основных величин $N_{g}^{{n + 1}}$, $S_{g}^{{n + 1}}$ эквивалентную систему разностных уравнений в виде

(4)

$\begin{gathered} \frac{{N_{g}^{{n + 1}} - N_{g}^{{n + 1/2}}}}{{{{{v}}_{g}}\tau }} + \alpha _{g}^{*}N_{g}^{*} = 0.5{\kern 1pt} \mathop \sum \limits_{g' = 1}^G \,\beta _{{g',g}}^{*}S_{{g'}}^{*} + 0.5f_{g}^{*}, \\ \frac{{S_{g}^{{n + 1}} - S_{g}^{{n + 1/2}}}}{{{{{v}}_{g}}\tau }} + \alpha _{g}^{*}S_{g}^{*} = \mathop \sum \limits_{g' = 1}^G \,\beta _{{g',g}}^{*}S_{{g'}}^{*} + f_{g}^{*}. \\ \end{gathered} $

Вычтем первое уравнение в (4) из второго, умноженного на 0.5, получим для вычисления потоков нейтронов $N_{g}^{{n + 1}}$ дополнительное эквивалентное уравнение

$\frac{{(0.5S_{g}^{{n + 1}} - N_{g}^{{n + 1}}) - (0.5S_{g}^{{n + 1/2}} - N_{g}^{{n + 1/2}})}}{{{{{v}}_{g}}\tau }} = - \alpha _{g}^{*}(0.5S_{g}^{*} - N_{g}^{*}),$

в котором полные потоки $S_{g}^{{n + 1}}$ находятся из второго уравнения в (4). Дополнительное уравнение можно использовать, например, для контроля, либо как основное уравнение для вычисления потока нейтронов $N_{g}^{{n + 1}}$.

Особенность системы уравнений (4) в том, что полный поток нейтронов $S_{g}^{{n + 1}}$ вычисляется из разностного уравнения, а не по квадратурной формуле из (2). Уравнение (3) можно интерпретировать как разностное уравнение, описывающее перенос нейтронов по пространству в вакууме, уравнения (4) – взаимодействие нейтронов с веществом. Если “большие” величины ${{N}_{{g,j + 1/2}}},\;{{N}_{{g,j - 1/2}}}$ вычислять по величинам с нижнего временного слоя ${{t}^{n}}$ или с верхнего ${{t}^{{n + 1}}}$, то разностная схема (3) будет явной или неявной соответственно. В обоих случаях разностное уравнение (3) аппроксимирует дифференциальное уравнение переноса нейтронов по пространству в вакууме:

(5)

$\frac{1}{{{{{v}}_{g}}}}\frac{{\partial {{N}_{g}}}}{{\partial t}} + \mu \frac{{\partial {{N}_{g}}}}{{\partial x}} = 0.$

В дальнейшем предполагаем, что вспомогательные потоки $N_{g}^{{n + 1/2}}$ находятся по известной разностной схеме повышенной точности, которая аппроксимирует уравнение (5), например, по схемам из [15] или [16]. Предполагая, что вопрос с выбором разностной схемы (3) решен, сосредоточим внимание на подходах к решению уравнений (4).

Если потоки $N_{g}^{*}$ и $S_{g}^{*}$ относятся к верхнему временному слою, т.е. $N_{g}^{*} = N_{g}^{{n + 1}}$ и $S_{g}^{*} = S_{g}^{{n + 1}}$, то разностная схема (4) будет неявной схемой Эйлера первого порядка по времени с погрешностью аппроксимации $O\left( \tau \right)$. В случае одногрупповой модели решение по неявной разностной схеме (4) , в том числе и по дополнительному уравнению, не вызывает затруднений и находится без итераций по интегралу столкновений.

Если вспомогательные потоки нейтронов $N_{g}^{{n + 1/2}}$ вычисляются по явной разностной схеме (3) , а величины $N_{g}^{*}$, $S_{g}^{*}$ относятся к верхнему временному слою ${{t}^{{n + 1}}}$, то схему (4) будем называть полунеявной разностной схемой.

Первые дифференциальные приближения (см. [17]) разностных уравнений в (4) выпишем в виде

$\frac{1}{{{{{v}}_{g}}}}\frac{{d{{N}_{g}}}}{{dt}} + {{\alpha }_{g}}{{N}_{g}} = 0.5{\kern 1pt} \mathop \sum \limits_{g' = 1}^G \,\beta {{\,}_{{g',g}}}{{S}_{{g'}}} + 0.5{{f}_{g}} + \left( {\tau {\kern 1pt} {\text{*}} - 0.5\tau } \right)\frac{{d{{F}_{1}}}}{{dt}},$

(6)

$\frac{1}{{{{{v}}_{g}}}}\frac{{d{{S}_{g}}}}{{dt}} + {{\alpha }_{g}}{{S}_{g}} = \mathop \sum \limits_{g' = 1}^G \,{{\beta }_{{g',g}}}{{S}_{{g'}}} + {{f}_{g}} + \left( {\tau {\kern 1pt} {\text{*}} - 0.5\tau } \right)\frac{{d{{F}_{2}}}}{{dt}},$

${{F}_{1}} = 0.5{\kern 1pt} \mathop \sum \limits_{g' = 1}^G \,{{\beta }_{{g',g}}}{{S}_{{g'}}} + 0.5{{f}_{g}} - {{\alpha }_{g}}{{N}_{g}},\quad {{F}_{2}} = \mathop \sum \limits_{g' = 1}^G \,{{\beta }_{{g',g}}}{{S}_{{g'}}} + {{f}_{g}} - {{\alpha }_{g}}{{S}_{g}}.~~$

Из первых дифференциальных приближений (6) следует, что система разностных уравнений (4) с постоянными коэффициентами $\alpha _{g}^{*},\;\beta _{g}^{*},\;f_{g}^{*}$ аппроксимирует систему ОДУ

(7)

$\begin{gathered} \frac{1}{{{{{v}}_{g}}}}\frac{{d{{N}_{g}}}}{{dt}} + {{\alpha }_{g}}{{N}_{g}} = 0.5{\kern 1pt} \mathop \sum \limits_{g' = 1}^G \,{{\beta }_{{g',g}}}{{S}_{{g'}}} + 0.5{{f}_{g}}, \\ \frac{1}{{{{{v}}_{g}}}}\frac{{d{{S}_{g}}}}{{dt}} + {{\alpha }_{g}}{{S}_{g}} = \mathop \sum \limits_{g' = 1}^G \,{{\beta }_{{g',g}}}{{S}_{{g'}}} + {{f}_{g}}~ \\ \end{gathered} $

с первым порядком по времени. Начальные данные для ОДУ (7) – это решения уравнения (5). Очевидно, что дополнительное уравнение аппроксимирует ОДУ

$\frac{1}{{{{{v}}_{g}}}}\frac{{d{\text{ln}}\left( {0.5{{S}_{g}} - {{N}_{g}}} \right)}}{{dt}} = - {{\alpha }_{g}},$

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

2.3. Подход к вычислению интеграла столкновений

Рассмотрим подход к вычислению интеграла столкновений, основанный на анализе выражения

$\mathop \sum \limits_{g' = 1}^G \,{{\beta }_{{g',g}}}{{S}_{{g'}}}$

как источника нейтронов, которые появились в группе $g$ из всех групп $g{\kern 1pt} '$ с плотностью потока ${{\bar {N}}_{g}}$. Интеграл столкновений в правой части ОДУ (7) представим в виде

$\mathop \sum \limits_{g' = 1}^G \,{{\beta }_{{g',g}}}{{S}_{{g'}}}\left( {t,x} \right) = {{B}_{g}}\mathop \sum \limits_{g' = 1}^G \,{{\bar {\beta }}_{{g',g}}}{{S}_{{g'}}}\left( {t,x} \right) = {{B}_{g}}{{\bar {S}}_{g}}\left( {t,x} \right)~,$

где введены следующие обозначения:

${{B}_{g}} = \mathop \sum \limits_{g' = 1}^G \,{{\beta }_{{g',g}}},\quad {{\bar {\beta }}_{{g',g}}} = \frac{{{{\beta }_{{g',g}}}}}{{{{B}_{g}}}},\quad \mathop \sum \limits_{g' = 1}^G \,{{\bar {\beta }}_{{g',g}}} = 1,\quad {{\bar {S}}_{g}}\left( {t,x} \right) = \mathop \sum \limits_{g' = 1}^G \,{{\bar {\beta }}_{{g',g}}}{{S}_{{g'}}}\left( {t,x} \right).$

Выполнив преобразования

${{\bar {S}}_{g}} = \mathop \sum \limits_{g' = 1}^G \,{{\bar {\beta }}_{{g',g}}}{{S}_{{g'}}} = \mathop \sum \limits_{g' = 1}^G {\kern 1pt} \mathop \sum \limits_{m = 1}^M \,{{\bar {\beta }}_{{g',g}}}{{N}_{{g',m}}}\Delta \mu = \mathop \sum \limits_{m = 1}^M {\kern 1pt} \mathop \sum \limits_{g' = 1}^G \,{{\bar {\beta }}_{{g',g}}}{{N}_{{g',m}}}\Delta \mu = \mathop \sum \limits_{m = 1}^M \,{{\bar {N}}_{{g,m}}}\Delta \mu ,$

получим для плотности потока нейтронов ${{\bar {N}}_{g}}$ и полного потока ${{\bar {S}}_{g}}$ выражения в виде

(8)

${{\bar {N}}_{g}} = \mathop \sum \limits_{g' = 1}^G \,{{\bar {\beta }}_{{g',g}}}{{N}_{{g'}}},\quad {{\bar {S}}_{g}} = \mathop \sum \limits_{m = 1}^M \,{{\bar {N}}_{{g,m}}}\Delta \mu .$

Система ОДУ (7) с введенными обозначениями записывается в виде

(9)

$\begin{gathered} \frac{1}{{{{{v}}_{g}}}}\frac{{d{{N}_{g}}}}{{dt}} + {{\alpha }_{g}}{{N}_{g}} = 0.5{{B}_{g}}{{{\bar {S}}}_{g}} + 0.5{{f}_{g}}, \\ \frac{1}{{{{{v}}_{g}}}}\frac{{d{{S}_{g}}}}{{dt}} + {{\alpha }_{g}}{{S}_{g}} = {{B}_{g}}{{{\bar {S}}}_{g}} + {{f}_{g}} \\ \end{gathered} $

с начальными данными ${{N}_{g}}\left( {0,x,\mu } \right) = N_{g}^{{n + 1/2}}\left( {{{t}^{{n + 1}}},x,\mu } \right)$, ${{S}_{g}}\left( {0,x} \right) = S_{g}^{{n + 1/2}}\left( {{{t}^{{n + 1}}},x} \right)$.

Если полный поток ${{\bar {S}}_{g}}$ известен, то решение системы уравнений (9) не вызывает затруднений. В этом случае систему уравнений (9) можно рассматривать как одногрупповую модель для описания эволюции нейтронов ${{N}_{g}}$. Коэффициент ${{B}_{g}}$ – это суммарный коэффициент размножения нейтронов группы $g$. Квадратурную формулу (8) для потока ${{\bar {S}}_{g}}$ можно рассматривать как аппроксимацию соответствующего интеграла в системе уравнений (1) для потока нейтронов ${{\bar {N}}_{g}}$. Следовательно, поток ${{\bar {S}}_{g}}~$ есть полный поток нейтронов ${{\bar {N}}_{g}}$. Поскольку нейтроны ${{\bar {N}}_{g}}$ принадлежат группе $g$, то эволюция этих нейтронов должна описываться так же, как и нейтронов ${{N}_{g}}$ этой группы, а именно, уравнениями (9). Поэтому плотность потока нейтронов ${{\bar {N}}_{g}}$ и полный поток ${{\bar {S}}_{g}}$ должны удовлетворять уравнениям (9), которые можно записать в виде

(10)

$\begin{gathered} \frac{1}{{{{{v}}_{g}}}}\frac{{d{{{\bar {N}}}_{g}}}}{{dt}} + {{\alpha }_{g}}{{{\bar {N}}}_{g}} = 0.5{{B}_{g}}{{{\bar {S}}}_{g}} + 0.5{{f}_{g}}, \\ \frac{1}{{{{{v}}_{g}}}}\frac{{d{{{\bar {S}}}_{g}}}}{{dt}} = - \left( {{{\alpha }_{g}} - {{B}_{g}}} \right){{{\bar {S}}}_{g}} + {{f}_{g}} \\ \end{gathered} $

с начальными данными

${{\bar {N}}_{g}}\left( {0,x,\mu } \right) = \bar {N}_{g}^{n} = \mathop \sum \limits_{g' = 1}^G \,{{\bar {\beta }}_{{g',g}}}N_{{g'}}^{{n + 1/2}}\left( {{{t}^{{n + 1}}},x,\mu } \right),$

${{\bar {S}}_{g}}\left( {0,x} \right) = \bar {S}_{g}^{n} = \mathop \sum \limits_{g' = 1}^G \,{{\bar {\beta }}_{{g',g}}}S_{{g'}}^{{n + 1/2}}\left( {{{t}^{{n + 1}}},x} \right).$

Система уравнений (10) – это аналог одногрупповой системы уравнений для потока нейтронов ${{\bar {N}}_{g}}$. Для вычисления полного потока ${{\bar {S}}_{g}}$ достаточно выписать решение второго уравнения в (10). Уравнение для полного потока ${{\bar {S}}_{g}}$ в (10) линейное и интегрируется в квадратурах. Поэтому решение системы уравнений (10) не вызывает каких-либо затруднений. Если коэффициенты постоянные, то полный поток ${{\bar {S}}_{g}}$ вычисляется из аналитического решения по формулам

(11)

$\begin{gathered} \bar {S}_{g}^{{n + 1}} = \left\{ \begin{gathered} {{\gamma }_{3}}\bar {S}_{g}^{n} + \left( {1 - {{\gamma }_{3}}} \right)\frac{1}{{{{\alpha }_{g}} - {{B}_{g}}}}{{f}_{g}},\quad {{\alpha }_{g}} \ne {{B}_{g}},~ \hfill \\ \bar {S}_{g}^{n} + {{{v}}_{g}}\tau {{f}_{g}},\quad {{\alpha }_{g}} = {{B}_{g}}, \hfill \\ \end{gathered} \right. \\ {{\gamma }_{3}} = {\text{exp}}\left( { - {{{v}}_{g}}\tau \left( {{{\alpha }_{g}} - {{B}_{g}}} \right)} \right). \\ \end{gathered} $

Следовательно, интеграл столкновений также вычисляется аналитически. Если данный подход к вычислению интеграла столкновений применить к разностным уравнениям (4), то систему разностных уравнений (4) запишем в виде

$\frac{{N_{g}^{{n + 1}} - N_{g}^{{n + 1/2}}}}{{{{{v}}_{g}}\tau }} + \alpha _{g}^{*}N_{g}^{*} = 0.5B_{g}^{*}\bar {S}_{g}^{*} + 0.5f_{g}^{*},$

$\frac{{S_{g}^{{n + 1}} - S_{g}^{{n + 1/2}}}}{{{{{v}}_{g}}\tau }} + \alpha _{g}^{*}S_{g}^{*} = B_{g}^{*}\bar {S}_{g}^{*} + f_{g}^{*}.$

Здесь

$B_{g}^{*} = \mathop \sum \limits_{g' = 1}^G \,\beta _{{g',g}}^{*},\quad \bar {\beta }_{{g',g}}^{*} = \frac{{\beta _{{g',g}}^{*}}}{{B_{g}^{*}}},\quad \bar {S}_{g}^{*} = \mathop \sum \limits_{g' = 1}^G \,\bar {\beta }_{{g',g}}^{*}S_{{g'}}^{*}.$

Уравнение для вычисления полного потока $\bar {S}_{g}^{*}$ запишем в виде

$\frac{{\bar {S}_{g}^{{n + 1}} - S_{g}^{{n + 1/2}}}}{{{{{v}}_{g}}\tau }} = - (\alpha _{g}^{*} - B_{g}^{*})\bar {S}_{g}^{*} + f_{g}^{*}.$

Следовательно, если величины $N_{g}^{*} = N_{g}^{{n + 1}}$ и $S_{g}^{*} = S_{g}^{{n + 1}}$, $\bar {S}_{g}^{*} = \bar {S}_{g}^{{n + 1}}$, то решение преобразованной неявной системы уравнений (4) не вызывает каких-либо затруднений, находится без итераций по интегралу столкновений и без обращения матриц.

2.4. Дискретно-аналитическая разностная схема

Рассмотрим подход к построению разностной схемы ДАРС, в которой для вычисления основных величин $N_{g}^{{n + 1}},\;S_{g}^{{n + 1}}$ используются аналитические решения ОДУ (9), (11) вместо разностных уравнений (4). Подставив выражение для ${{\bar {S}}_{g}}$ из уравнения (11) в уравнения (9), получим для вычисления потока нейтронов ${{N}_{g}}$ и полного потока ${{S}_{g}}$ линейную систему ОДУ:

(12)

$\begin{gathered} \frac{1}{{{{{v}}_{g}}}}\frac{{d{{N}_{g}}}}{{dt}} + {{\alpha }_{g}}{{N}_{g}} = 0.5{{B}_{g}}{{\gamma }_{3}}\bar {S}_{g}^{n} + 0.5\left( {1 + \frac{{1 - {{\gamma }_{3}}}}{{{{\alpha }_{g}} - {{B}_{g}}}}{{B}_{g}}} \right){{f}_{g}}, \\ \frac{1}{{{{{v}}_{g}}}}\frac{{d{{S}_{g}}}}{{dt}} + {{\alpha }_{g}}{{S}_{g}} = {{B}_{g}}{{\gamma }_{3}}\bar {S}_{g}^{n} + \left( {1 + \frac{{1 - {{\gamma }_{3}}}}{{{{\alpha }_{g}} - {{B}_{g}}}}{{B}_{g}}} \right){{f}_{g}}. \\ \end{gathered} $

Поскольку при построении разностных схем предполагается, что коэффициенты поглощения, размножения и функция источника постоянны на интервале интегрирования по времени $\left[ {{{t}^{n}},{{t}^{{n + 1}}}} \right]$, а разностные уравнения (4) аппроксимируют ОДУ (7), то аналитические решения ОДУ (7) могут быть взяты как численные решения уравнений (1) вместо решений по разностной схеме (4) . Аналитические решения ОДУ (7) сводятся к решениям ОДУ (12) и записываются в виде

$N_{g}^{{n + 1}} = {{\gamma }_{2}}N_{g}^{{n + 1/2}} + 0.5\left( {{{\gamma }_{3}} - {{\gamma }_{2}}} \right)\bar {S}_{g}^{n} + 0.5\frac{{1 - {{\gamma }_{3}}}}{{{{\alpha }_{g}} - {{B}_{g}}}}{{f}_{g}},$

(13)

$S_{g}^{{n + 1}} = {{\gamma }_{2}}S_{g}^{{n + 1/2}} + \left( {{{\gamma }_{3}} - {{\gamma }_{2}}} \right)\bar {S}_{g}^{n} + \frac{{1 - {{\gamma }_{3}}}}{{{{\alpha }_{g}} - {{B}_{g}}}}{{f}_{g}},$

$\bar {S}_{g}^{n} = \mathop \sum \limits_{g' = 1}^G \,\bar {\beta }_{{g',g}}^{n}S_{{g'}}^{{n + 1/2}},\quad {{\gamma }_{2}} = \exp \left( { - {{{v}}_{g}}\tau {{\alpha }_{g}}} \right).$

Разностная схема (13) – это схема с положительными коэффициентами. Поэтому решения всегда положительные и нет ограничения на выбор шага интегрирования по времени. Однако этот плюс может стать минусом в ситуациях, когда шаг интегрирования по времени большой, а искомые функции и коэффициенты имеют большой градиент. В таком случае усреднение величин будет “грубым” на интервале интегрирования по времени $\left[ {{{t}^{n}},{{t}^{{n + 1}}}} \right]$, а погрешность решения – большой. Поэтому ограничения на выбор шага интегрирования по времени могут быть связаны с требованиями к точности численных решений. Простейший подход к устранению таких ситуаций – это уменьшение шага интегрирования по времени, что приводит к увеличению времени счета. Второй подход – это построение схемы предиктор–корректор повышенной точности типа Рунге–Кутты. В этом случае коэффициенты и полный поток ${{\bar {S}}_{g}}$ предварительно рассчитываются на этапе предиктор в момент времени $t{\kern 1pt} * = {{t}^{n}} + \tau {\kern 1pt} *,\;\tau {\kern 1pt} * \geqslant 0.5\tau $. Основные величины рассчитываются на этапе корректор в момент времени ${{t}^{{n + 1}}}$. Если $\tau {\kern 1pt} * = 0.5\tau $, то погрешность аппроксимации будет порядка $O({{\tau }^{2}})$. Если на этапе предиктор применить двухточечную формулу Гаусса (см. [18]), то получим разностную схему с погрешностью аппроксимации порядка $O({{\tau }^{3}})$.

Третий кардинальный подход – это вычисление основных величин с контролем точности в каждой точке по критерию точности, либо проведение двух расчетов с шагами интегрирования $\tau $ и $0.5\tau $ и сравнение результатов расчетов. Если результаты не удовлетворяют заданной точности, то шаг интегрирования уменьшается в два раза, и расчет в точке повторяется. Если результаты удовлетворяют заданной точности, то в качестве решения принимается последний результат. Такая технология счета является стандартной для решения ОДУ, обеспечивает гарантированную точность и эффективность расчетов (см. [18]).

3. РЕЗУЛЬТАТЫ ЧИСЛЕННЫХ РЕШЕНИЙ МОДЕЛЬНОЙ ЗАДАЧИ

Работоспособность методики ДАРС проверялась сравнением результатов расчетов с результатами расчетов по одномерному аналогу методики из [13] и с точными решениями задачи из [12]. Одномерный аналог методики [13] назовем методикой ГРАНАТ.

Задача. Дана сферически-симметричная область $\left| r \right| \leqslant 1$. В начальный момент времени задано распределение нейтронов, которое рассчитывается по формуле

${{N}_{g}} = {{N}_{{0,g}}}\left( {1 + \mu r} \right),\quad {{N}_{{0,g}}} = 1,\quad g = 1,\;2,\;3,\;4.$

Источники ${{f}_{g}} = 0$. Скорости в группах: ${{{v}}_{1}} = 10,$ ${{{v}}_{2}} = 5$, ${{{v}}_{3}} = 2,$ ${{{v}}_{4}} = 1$. Коэффициенты поглощения в группах: ${{\alpha }_{1}} = 0.1$, ${{\alpha }_{2}} = 0.2,$ ${{\alpha }_{3}} = 0.5$, ${{\alpha }_{4}} = 1$. На правой границе $r = 1$ задан поток нейтронов $N = \left( {1 + \mu } \right)\exp \left( { - t} \right)$. Задана матрица коэффициентов ${{\beta }_{{g',g}}}$:

${{M}_{4}} = \left( {\begin{array}{*{20}{c}} {0.8~~~~~0.08~~0.07~~0.05} \\ {0.1~~~~~0.7~~~~0.12~~0.08} \\ {0.08~~0.15~~0.6~~~~0.17} \\ {0.05~~0.15~~0.3~~~~0.5} \end{array}} \right).$

Решения с матрицей ${{M}_{4}}$ соответствуют общему случаю, когда все группы взаимодействуют между собой. Точное решение для потока нейтронов ${{N}_{g}}$ записывается в виде (см. [13])

${{N}_{g}} = {{N}_{{0,g}}}\left( {1 + \mu r} \right){\text{exp}}\left( { - \lambda t} \right),\quad \lambda = {{\alpha }_{g}}{{{v}}_{g}},$

для полного потока нейтронов – $S = 2{{N}_{0}}{\text{exp}}\left( { - \lambda t} \right)$. Константы ${{\alpha }_{g}},\;{{{v}}_{g}}$ подобраны так, что $\lambda = 1$. Поэтому решения в группах совпадают. Требуется рассчитать полные потоки в момент времени $t = 1$.

Задачи считались на последовательности сгущающихся разностных сеток с числом интервалов по пространству и по углам 12 × 12, 25 × 25, 50 × 50, 100 × 100, 200 × 200 с шагами интегрирования по времени 0.002, 0.001, 0.0005, 0.00025, 0.000125 соответственно.

На фиг. 1 представлены зависимости от шагов интегрирования по пространству погрешностей аппроксимации по четырем группам, которые получены в расчетах на сходимость по методикам ДАРС и ГРАНАТ в момент времени $t = 1$.

Фиг. 1.

Зависимости от $h{\text{/}}{{h}_{0}}$ погрешностей аппроксимации по четырем группам 1, 2, 3, 4: квадраты – методика ГРАНАТ, кружочки – методика ДАРС.

На фиг. 2 представлены распределения полных потоков в области по четырем группам в расчетах по методике ДАРС на трех различных разностных сетках по пространству и по углу в момент времени $t = 1$.

Фиг. 2.

Зависимости от $x$ полных потоков $S$ по четырем группам: кружочки – разностная сетка 12 × 12, сплошные линии – 50 × 50, пунктирные – 200 × 200, квадраты – точное решение.

Из анализа графиков на фиг. 1 следует, что погрешности аппроксимации по методикам ДАРС и ГРАНАТ удовлетворительно согласуются между собой, стремятся к нулю при уменьшении шагов интегрирования по времени и по пространству, а численные решения сходятся к точному решению.

Из анализа графиков на фиг. 2 следует, что результаты расчетов по методике ДАРС сходятся к точному решению при уменьшении шагов интегрирования по времени и по пространству и удовлетворительно согласуются между собой на разностных сетках 50 × 50 и далее. На “грубых” разностных сетках полные потоки в группах далеки не только от аналитического решения, но и различаются между собой. Различие исчезает на подробных разностных сетках.

Время счета задачи на разностной сетке 200 × 200 по методике ДАРС в 1.8 раза меньше, чем по методике ГРАНАТ.

Удовлетворительное совпадение результатов расчетов по методике ДАРС с точными решениями задачи и с решениями по методике ГРАНАТ позволяет сделать вывод, что уравнения (10) адекватно описывают эволюцию потока нейтронов ${{\bar {N}}_{g}}$ и полного потока ${{\bar {S}}_{g}}$. Поэтому предложенный подход к вычислению интеграла столкновений можно признать оправданным и эффективным.

ЗАКЛЮЧЕНИЕ

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

Решение уравнения переноса нейтронов в многогрупповой модели сведено к решениям в одногрупповой модели.

Эффективность достигнута за счет вычисления интеграла столкновений из аналитических решений ОДУ, которые моделируют эволюцию нейтронов, появившихся в группе $g$ из всех групп $g{\kern 1pt} '$.

Результаты расчетов модельной задачи согласуются с точными решениями и с решениями, полученными по методике ГРАНАТ.

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

Список литературы

Карлсон Б., Белл Дж. Решение транспортного уравнения Sn-методом. Сб. “Физика ядерных реакторов”. М.: Атомиздат, 1959. С. 408–432.
Карлсон Б., Белл Дж. Численное решение задач кинетической теории нейтронов. Сб. “Теория ядерных реакторов”. М.: Госатомиздат, 1963. С. 243–258.
Марчук Г.И., Лебедев В.И. Численные методы в теории переноса нейтронов. М.: Атомиздат, 1981.
Шагалиев Р.М., Шумилин В.А., Алексеев А.В. и др. Математические модели и методики решения многомерных задач переноса частиц и энергии, реализованные в комплексе “САТУРН-3” // ВАНТ. Сер. Матем. моделирование физ. процессов. 1999. Вып. 4. С. 20–26.
Гаджиев Ф.Д., Кондаков И.А., Писарев В.Н., Стародумов О.И., Шестаков А.А. Метод дискретных ординат с искусственной диссипацией (DDAD-схема) для численного решения уравнения переноса нейтронов // ВАНТ. Сер. Матем. моделирование физ. процессов. 2003. Вып. 4. С. 13–24.
Feautrier P.C.R. Sur la resolution numerique de I’equation de transfert // Acad. Sci. Paris. 1964. V. 258. P. 3198–3210.
Rybicki G. A modified Feautrier method // J. of Quantitative Spectroscopy and Radiative Transfer. 1971. V. 11. P. 589–596.
Михалс Д. Звездные атмосферы. М.: Мир, 1982.
Яненко Н.Н. Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, 1967.
Марчук Г.И. Методы расщепления. М.: Наука, 1988.
Моисеев Н.Я., Шмаков В.М. Модифицированный метод расщепления для решения кинетического уравнения переноса частиц // Ж. вычисл. матем. и матем. физ. 2016. Т. 56. № 8. С. 1480–1490.
Кондаков И.А., Селезнев В.Н., Стародумов О.И., Шестаков А.А. Аналитические тесты для решения задач переноса частиц численными методами // ВАНТ. Сер. Матем. моделирование физ. процессов. 2003. Вып. 2. С. 28–43.
Арсентьев А.П., Писарев В.Н. Особенности применения $TVD$-подхода к $D{{S}_{n}}$ методу решения трехмерного уравнения переноса нейтронов в криволинейной системе координат // ВАНТ. Сер. Матем. моделирование физ. процессов. 2011. Вып. 1. С. 13–39.
Годунов С.К., Забродин А.В., Прокопов Г.П., Иванов М.Я. Численное решение многомерных задач газовой динамики. Под ред. С.К. Годунова. М.: Наука, 1976.
Моисеев Н.Я., Силантьева И.Ю. Разностные схемы произвольного порядка аппроксимации для решения линейных уравнений переноса с постоянными коэффициентами методом Годунова с антидиффузией // Ж. вычисл. матем. и матем. физ. 2008. Т. 48. № 7. С. 1282–1293.
Моисеев Н.Я. Неявные разностные схемы бегущего счета повышенной точности // Ж. вычисл. матем. и матем. физ. 2011. Т. 51. № 5. С. 920–935.
Шокин Ю.И. Метод дифференциального приближения. Новосибирск: Наука, 1979.
Каханер Д., Моулер К., Неш С. Численные методы и программное обеспечение. М.: Мир, 2001.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал