2a.tex

\section{Алгоритм моделирования пространственно-неоднородного уравнения
Больцмана весовой схемой Бернулли.}

Рассмотрим область $F$ заполненную смесью из $N$ сортов газов. Основными характеристиками частиц газа будут их масса $m_i$, диаметр $d_i$, и функция
распределения  по скоростям $ f_i (\vec x, \vec v, t)$ в точке $\vec x$ в момент времени $t$. Эволюция частиц газа описывается системой 
уравнений Больцмана:
\begin{align}
\frac{\partial f_i(\vec v)}{\partial t} + \left( \vec v, \frac{\partial f_i(\vec v)}{\partial x}\right) = \sum_{k=1}^{N} I[f_i, f_k]\\
I[f_i, f_k] = \int (f_i(\vec v)^*f_k(\vec w)^* - f_i(\vec v)f_k(\vec w)) Bb\left( |\vec u|, \frac{(\vec u, \vec n)}{|\vec u|} \right) d\vec n \, d\vec w
\end{align}

Здесь $\vec n \in \mathbb{S}^2$ - вектор единичной сферы, указывающий направление относительной скорости частиц после столкновения, 
$\vec u = \vec v - \vec w$ - относительная скорость частиц.

Для уравнения Больцмана, по-видимому, впервые Грэдом в (3) было предложено искать решение на 
каждом временном шаге $\tau$ чередованием свободномолекулярного движения и пространственно-однородной релаксации. Эта же идея и используется в схеме испытаний Бернулли.


Алгоритм расщепления заключается в следующем: временная переменная разбивается на интервалы $\Delta t$, и на каждом интервале вместо уравнения Больцмана (1) рассматривается последовательность более простых задач: вначале решается уравнение свободномолекулярного течения:
\begin{equation}
\frac{\partial f_i}{\partial t} + \left( v, \frac{\partial f}{\partial x}\right) = 0
\end{equation}

а затем на том же временном интервале - уравнение пространственно-однородной релаксации:
\begin{equation}
\frac{\partial f_i}{\partial t} = \sum_{i = 1}^{N} I[f_i, f_k]
%\key
\end{equation}

При этом начальным условием для (3) служит решение уравнения (2), а в свою очередь решение
уравнения (3) дает начальное условие (2) на следующем шаге по времени.


Для упрощения дальнейшего обсуждения введем следующие определения.
Модельная частица - набор из 6 переменных, вес $w_i \in \mathbb{R}_+$ - количество реальных частиц, обладающих одинаковыми параметрами, которые 
представляет модельная, скорость $v \in \mathbb{R}^3$, масса $m_i \in \mathbb{R}_+$ и диаметр $d_i \in \mathbb{R}_+$ - характеристики реальных частиц.\\
Ячейка - набор из $n$ частиц.\\
Расчетная область - множество ячеек, задающие пространственную область, в которой ищется решение.
Реализация - серия состояний расчетной области, проиндексированных параметром $\{t_i\}$.


Численное решение задачи Коши методом испытаний Бернулли с весами строится следующим образом. Фиксируется $K$ независимых реализаций случайного процесса, 
расчетная область разбивается на ряд ячеек, в которых содержится постоянное число $n$ модельных частиц. Функцией распределения скоростей частиц
сорта $l$ в ячейке будет $\phi_l(v) = \sum_{i=1}^{N}w_i[\vec v_i = v]$. Алгоритм моделирования можно представить в виде следующей последовательности
действий на каждом шаге $\Delta t$ для всех реализаций:
\begin{description}
\item[Шаг 1.] Столкновительная релаксация частиц.
\item[Шаг 2.] Перемещение частиц в соседние ячейки.
\item[Шаг 3.] Перемешивание частиц между независимыми реализациями случайного процесса.
\end{description}


Случайный процесс строится как последовательность столкновений, разыгрываемых по Монте-Карло. Отличительной чертой схемы Бернулли является возможность использования в расчетах предельно малого числа представителей каждого компонента с сохранением больцмановской частоты столкновений. Для каждой пары частиц
в ячейке с вероятностью\cite{Korolev85}:
\begin{equation}
p_{ij}^{lk} = K \frac{n_k \sigma_{lk} g_{ij}^{lk} \Delta t}{N_k} = K \frac{n_l \sigma_{lk} g_{ij}^{kl} \Delta t}{N_l}, \quad K = \left\{\begin{array}{cc}
1&\text{, при }l \neq k \\
\\
\dfrac{N_l}{N_l - 1}&\text{, при }l = k \\
\end{array}
 \right.
\end{equation}
\begin{description}
 \item $N_l, N_k$ – количество модельных представителей компонент $l$ и $k$
\end{description}
происходит столкновение, а скорости определяются в соответствии с механикой упругого удара. Для частиц с разными весами при достоверно произошедшем  
столкновении скорости частиц не обязательно изменятся. Модельные частицы, обладающие меньшим весом,$w_{min}$ изменят скорости с вероятностью 1, а 
частицы с большим весом $w_{max}$ изменят скорость с вероятностью $w_{min}/w_{max}$. Такая особенность моделирования связана с парностью столкновений.
Общее число столкнувшихся пар частиц будет $w_{min}$, а $w_{max} - w_{min}$ частиц не будут участвовать в столкновении. Таким образом образуется
3 множества частиц, и для сохранения постоянного количества модельных частиц используется метод Монте-Карло, чтобы определить новые параметры модельных
представителей.


Введение весовых множителей позволяет представлять каждый из $S$ компонентов $A_l$ концентраций $n_l, l = 1, 2, \ldots,S$ рассматриваемой газовой смеси одним и тем же числом модельных частиц $N$ независимо от значений $n_l$. В результате все многообразие реальных молекул каждого компонента заменяется $N$ подмножествами частиц, обладающих одним и тем же вектором скорости внутри каждого подмножества (см. рис. \ref{img1}). Любое подмножество частиц может быть подвергнуто внешнему воздействию различной природы и интенсивности.


Таким образом, каждый модельный представитель фактически характеризует собой поведение значительного числа реальных частиц. Мощности каждого подмножества соответствует "вес"\ $w$  его модельного представителя. 
\begin{figure}[ht]
%   \centering
 \begin{center}
  \includegraphics[width=0.7\textwidth]{img1.pdf}
  % img1.eps: 0x0 pixel, 300dpi, 0.00x0.00 cm, bb=0 0 595 842
 \end{center}
 \caption{Весовая схема моделирования}
 \label{img1}
\end{figure}
В начальный момент моделирования все "веса"\ для простоты принимаются равными $w_l(0) = n_l / N$. При химических превращениях в процессе расчета значения "весов"\ могут изменяться, однако неизменно для каждого компонента смеси выполняется условие: $\sum w_l(t) = n(t)$.


Существенно подчеркнуть, что при взаимодействии модельных молекул скорость частицы, обладающей большим "весом"\ изменяется не обязательно, а лишь с вероятностью, равной отношению меньшего "веса"\ к большему. При этом всегда строго соблюдается принцип парности взаимодействия реальных молекул. На рис. \ref{img2} показана схема реализации весового подхода при моделировании упругого взаимодействия двух модельных представителей различной мощности.

\begin{figure}[ht]
 \centering
 \includegraphics[width=0.7\textwidth]{img2.pdf}
 \caption{Весовой подход при моделировании столкновений}
 \label{img2}
\end{figure}


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


Была разработана вычислительная процедура, при реализации которой в течение всего процесса моделирования в отдельной пространственной ячейке каждый компонент смеси представляется одним и тем же постоянным числом частиц. Алгоритм расчета пространственного сдвига на малом временном шаге $\Delta t$ сводится к моделированию методом Монте-Карло акта обмена частицами между соседними ячейками. Далее для упрощения символики индексы компонентов газовой смеси опускаются.


Для иллюстрации алгоритма ниже рассмотрено двумерное плоское течение. Пусть все пространственные ячейки имеют прямоугольную форму и одинаковые размеры $\Delta x$ и $\Delta y$, число ячеек по осям $OX$ и $OY$ одинаково и равно $K$. Без потери общности рассуждений будем считать, что частица имеет положительные составляющие вектора скорости $cx>0$, $c_y>0$. Тогда схема реализации пространственного сдвига частицы $A_{kn}$  весом $w_{kn}$, находящейся в пространственной ячейке с номером $k$ по оси $OX$ и номером $n$ по оси $OY$ будет такова:
\begin{description}
 \item[Шаг 1.] Вычисляется "доля"\ частиц в группе $A_{kn}$, перешедших в ячейки
с номерами $(k+1, n)$, $(k+1, n+1)$, и $(k, n+1)$:
$$\delta_{(k+1)n} = w_{kn} (\frac{c_x \Delta t}{\Delta x}) (1 - \frac{c_y \Delta t}{\Delta y});$$
$$\delta_{(k+1)(n+1)} = w_{kn} (\frac{c_x \Delta t}{\Delta x}) (\frac{c_y \Delta t}{\Delta y});$$
$$\delta_{k(n+1)} = w_{kn} (1-\frac{c_x \Delta t}{\Delta x}) (\frac{c_y \Delta t}{\Delta y})$$.
 \item[Шаг 2.] "Доля"\ частиц, оставшихся в ячейке $(k, n)$, уменьшается на величину: 
$$\delta_s = \delta_{k(n+1)} + \delta_{(k+1)(n+1)} + \delta_{(k+1)n},$$
т.е. их вес становится равным: 
$w^{*}_{kn} = w_{kn} - \delta_s$.
 \item[Шаг 3.] Веса частиц $A_{(k+1)n}$, $A_{k(n+1)}$, $A_{(k+1)(n+1)}$  увеличиваются соответственно на величины $\delta_{(k+1)n}$, $\delta_{k(n+1)}$, $\delta_{(k+1)(n+1)}$.
 \item[Шаг 4.] Скорости частиц $A_{(k+1)n}$, $A_{k(n+1)}$, $A_{(k+1)(n+1)}$  изменяются с вероятностями, соответствующими весовому подходу: 
$$P_{(k+1)n} (c_{kn} \rightarrow c_{(k+1)n} ) = \frac{\delta_{(k+1)n}}{w_{(k+1)n}}$$
$$P_{k(n+1)} (c_{kn} \rightarrow c_{k(n+1)} ) = \frac{\delta_{k(n+1)}}{w_{k(n+1)}}$$
$$P_{(k+1)(n+1)} (c_{kn} \rightarrow c_{(k+1)(n+1)} ) = \frac{\delta_{(k+1)(n+1)}}{w_{(k+1)(n+1)}}$$
\end{description}

\begin{figure}[t]
 \centering
 \includegraphics[width=0.95\textwidth]{move2.png}
 \caption{Схема пространственного сдвига.}
 \label{img_move}
\end{figure}

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


Алгоритм «поддерживает» в ячейках постоянное число частиц независимо от интенсивности переходов как на равномерных сетках так и на переменных. При этом, алгоритм естественным образом сохраняет это условие при работе с адаптивными сетками. Без каких-либо принципиальных изменений в алгоритме возможно рассмотрение осесимметричных и трехмерных течений.