test_merthods.tex

\newpage

\section{Методы исследования}

\subsection{Общая постановка задачи}

\subsubsection{Взаимодействие газа и оптической решетки}

\hspace{0.5cm} Газ нейтральных частиц, находящийся в поле оптической решетки, испытывает нерезонансное воздействие со стороны электрического поля и поляризуется, преобретая ненулевой дипольный момент. Потенциал взаимодействия пространственно неоднородного электрического поля и диполя в квазиэлектростатическом приближении \cite{Boyd},\cite{Takekoshi} имеет вид \eqref{eq:fieldU}, где $\alpha$ - статическая поляризуемость частиц газа.

\begin{equation} \label{eq:fieldU}
U(x,t) = - \int_0^{\vec{E}} \alpha \vec{ \mathcal{E}} d\vec{ \mathcal{E}} = - \alpha \vec{E}(x,t)^2 /2 
\end{equation}

Поляризуемость в общем виде имеет симметричный тензорный вид, однако его можно выбором осей описать тремя числами $\alpha_{xx},\alpha_{yy},\alpha_{zz}$ \cite{Parcell}. Она может быть записана \cite{Larsen1999} в виде \eqref{eq:polariz}, где $\alpha_{\parallel}$ и $\alpha_{\perp}$ описывают компоненты параллельную и перпендикулярную молекулярной оси, а $\theta$ задает угол между осью молекулы и направлением поляризации. Можно упростить вид величины поляризуемости, принебрегая вкладом ориентации молекулы и усреднив по всем возможным направлениям оси молекулы.

\begin{equation} \label{eq:polariz}
\alpha = (\alpha_{\parallel} - \alpha_{\perp})cos^2 \theta + \alpha_{\perp} =_{(aver)} (\alpha_{\parallel} + 2 \alpha_{\perp})/3
\end{equation}

Распространение электромагнитного излучения в среде описывается системой уравнений Максвелла \eqref{eq:MaxwellFull} с отсутствующими токами и свободными зарядами (частицы имеют нейтральный заряд и ионизации нет).

\begin{equation} \label{eq:MaxwellFull}
\left\{ \begin{array}{rl}
&
div(\varepsilon E) = 4\pi\rho, \rho \equiv 0, \\&
div B = 0,  \mu = 1, \\&
rot E = -\frac{1}{c} \frac{\partial B}{\partial t}, \\&
rot H = 4\pi \frac{j}{c} + \frac{1}{c} \frac{\partial (\varepsilon E)}{\partial t}, j \equiv 0
       \end{array} \right.
\end{equation}

Если нас интересует лишь электрическая компонента, то волна описывается волновым уравнением. В среде с заданной диэлектрической проницаемостью $\varepsilon(x,t)$ оно имеет вид \eqref{eq:waveeq}. 

\begin{equation} \label{eq:waveeq}
\Delta E - \frac{1}{c^2}\frac{\partial^2 (\varepsilon E)}{\partial t^2} = 0
\end{equation}

В рассматриваемой задаче диэлектрическая проницаемость газа как оптической среды определяется локальным значением плотости газа. Так, диэлектрическая проницаемость суть квадрат показателя преломления, можно записать в форме \eqref{eq:epsilonN}. 

\begin{equation} \label{eq:epsilonN}
\varepsilon (x,t) = n^2 (x,t) = 1 + \frac{\alpha N(x,t)}{\varepsilon_0}
\end{equation}

Влияние неоднородного распределения плотности газа на решение волнового уравнения или уравнений Максвелла часто не учитывается ввиду того, что отклонение $\Delta n$ в $n =  n_0 + \Delta n(x,t)$ мало, т.е. $|\Delta n| / n \ll 1$, и считается, что волна распространяется через среду с постоянной проницаемостью $\varepsilon = 1 + \frac{\alpha N_0}{\varepsilon_0}$, где $N_0$ - средняя плотность газа в исследуемом объеме. Несмотря на то, что отклонение показателя преломления от средней величины достаточно малы для подобного упрощения (например, порядка $10^{-5}$), амплитуда возмущений плотности газа могут достигать $30-40\%$ от средней величины плотности \cite{MyRGD}. Решение уравнений, описывающих распространение электромагнитных волн в неоднородной среде, представляет особенный интерес в данной задаче.

Оптическая решетка создается двумя источниками лазерного излучения, поэтому полное поле $E(x,t)$ является суперпозицией двух электрических полей от каждого источника. В качестве основы для дальнейшего анализа рассмотрим случай интерференции двух противонаправленных монохроматических источников одинаковой интенсивности в приближении однородной среды. Тогда две плоские волны с одинаковыми частотами образуют полное поле $E = E_0 sin(kx-\omega t) + E_0 sin(-kx-\omega t)$, а квадрат поля можно записать в виде \eqref{eq:fieldSimple}, осреднив быстро осциллирующую временную компоненту.

\begin{equation} \label{eq:fieldSimple}
<E^2 (x,t)>_{\tau} \ = \ <4 E_0^2 cos^2(kx) sin^2(\omega t)>_{\tau} \ = \ E_0^2 (1 + cos(2kx))
\end{equation}

Амплитуда поля напрямую связана с интенсивностью излучения - чаще всего, именно эта величина задается в эксперименте. Это соотношение может быть записано как \eqref{eq:Etointens}.

\begin{equation} \label{eq:Etointens}
E_0^2 = \frac{2 I}{c \varepsilon_0}
\end{equation}

Частицы газа испытывают ускорение, находясь во внешнем неоднородном потенциале оптической решетки; сила, действующая на молекулы может быть записана в виде \eqref{eq:forceX}. Изменение скорости частицы в каждый момент времени кроме градиента электрического поля определяется параметром $\alpha/m$, где $m$ - масса молекулы. Очевидно, что частицы с разным отношением поляризуемости к массе будут вести себя по разному - этот факт может быть использован, например, при разделении газовых смесей с помощью оптических решеток \cite{ShevIvanSepar}.

\begin{equation} \label{eq:forceX}
F(x) = -\triangledown_x U(x,t) = \alpha \triangledown_x E^2(x) /2 = - \alpha k E_0^2 sin(2kx)
\end{equation}

В более общем виде, конечно, частицы могут испытывать силу, зависящую не только от координаты $x$ вдоль оси пересечения лучей (рис. \ref{fig:simpleField}) , но направленную по радиусу в случае осесимметричной задачи. Другими словами,

\begin{equation} \label{eq:forceFull}
\vec{F}(\vec{r},t) = -\vec{ \triangledown} U(\vec{r},t) = - (\vec{e_r}\triangledown_r + \vec{e_x}\triangledown_x) U(\vec{r},t)
\end{equation}

\begin{figure}[h]
\centering
\includegraphics[width=250px]{simpleField}
\caption{Пространственная зависимость максимальной силы, действующей на молекулы газов с разными поляризуемостями для лазерных пучков мощностью $0.13J$.}
\label{fig:simpleField}
\end{figure}

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


Радиальная компонента силы возникает в случае неоднородности интенсивности лазерного пучка в плоскости, перпендикулярной направлению распространения электромагнитной волны. При гауссовом профиле радиальной зависимости интенсивности излучения (см. \S 3.1.2) градиентная сила направлена в центральную область пучка и собирает молекулы в область наибольшей интенсивности.


\subsubsection{Характеристики оптических решеток}\label{sec:laserQ}

\hspace{0.5cm} Рассмотрим два источника лазерного излучения с частотами $\omega_{1}$ и $\omega_{2}$, направленных навстречу друг другу. Соответствующие волновые вектора $\textbf{k}_{1}$ и $\textbf{k}_{2}$ связаны дисперсионным соотношением с частотами лазеров, и скалярное произведение $(\textbf{k}_{1},\textbf{k}_{2}) < 0$. В общем случае, конечно, направления распространения волн не должны быть параллельными. На практике, оптическая решетка создается двумя лазерами, направленными под углом $\varphi < 180^\circ$ друг к другу. При этом оптическая решетка не является одномерной, но в дальнейшем мы всегда будет говорить о центральной области пересечения лазерных лучей вдоль оси $x$. В этом предположении задача становится пространственно одномерной, и все интересующие нас величины будут иметь зависимость только от координаты $x$. Для сложных сред, например, с меняющимся показателем преломления, будем считать $\textbf{k}_{1,2} = \textbf{k}_{1,2}(x,t)$. Тогда электрическую компоненту каждого из электромагнитных полей лазеров можно в общем виде записать как (\ref{eq:fullE}). Здесь $\mathbb{E}_{1,2}$ - полное поле, распространяющееся в направлении $\textbf{k}_{1,2}$, а $E_{1,2}$ - соответствующие комплексные амплитуды.

\begin{equation} \label{eq:fullE}
\left\{ \begin{array}{rl}
\mathbb{E}_{1} = E_{1}(x,t)\ e^{ik_1(x,t)x-i \omega_1 t} + c.c, \\
\mathbb{E}_{2} = E_{2}(x,t)\ e^{ik_2(x,t)x-i \omega_2 t} + c.c
       \end{array} \right.
\end{equation}

Как было показано ранее, нас интересует не мгновенное распределение полей, а осредненный за некоторое время $\tau \gg \omega_{1,2}^{-1}$ (и меньшее, чем характерное время процессов в газе) квадрат полного поля ${E^2\sim <(\mathbb{E}_1+\mathbb{E}_2)^2>_\tau}$, связанный с дипольным взаимодействием газа и излучения. На рис.\ref{fig:fullE1} и \ref{fig:fullE2} представлены примеры вида интенсивности поля, создаваемый оптической решеткой при пересечении лучей лазерного излучения под углом $\varphi$. 

\begin{figure}[!tbp]
  \centering
  \subfloat[$\varphi = 135^\circ$]{\includegraphics[width=0.5\textwidth]{figure_1}\label{fig:fullE1}}
  \hfill
  \subfloat[$\varphi = 178^\circ$]{\includegraphics[width=0.5\textwidth]{figure_2}\label{fig:fullE2}}
  \caption{Мгновенное распределение интенсивности интерференции двух плоских монохроматических волн, распространяющихся под углом $\varphi$ с радиальным гауссовым профилем интенсивности.}
\end{figure}

В случае, если частоты излучения отличаются, то оптическая решетка приобретает ненулевую скорость $\xi = \Omega \ /\ q$, где $\Omega = \omega_1 - \omega_2$ и $q = |\textbf{k}_1-\textbf{k}_2| = k_1+k_2$ (в случае $\varphi = 2\pi$). С другой стороны, оптическая решетка может быть создана не двумя лазерами, а единственным лазером и движущимся со скоростью  $  \textbf{v}_m \ || \ \textbf{k}$ зеркалом. Тогда, за счет эффекта Доплера, отраженный луч будет иметь частоту $\omega_2 \approx \omega_1(1-2v_m/c)$ при $|v_m|/c \ll 1$. Несложно показать, что и в этом случае скорость оптической решетки определяется соотношением (\ref{eq:latticeSpeed}).

\begin{equation}\label{eq:latticeSpeed}
\xi = \frac{\omega_1-\omega_2}{k_1+k_2} = v_m
\end{equation}

Вид поля оптической решетки определяется множителем ${cos(q(x-\xi_0t)-\beta t^2)}$, где ${(x-\xi_0t)}$, как уже было сказано ранее, задает ненулевую скорость оптической решетки, и, при условии $\Omega \ll \omega_{1,2}$, коэффициент $\beta = \frac{d\Omega}{dt}$ возникает при наличии линейной модуляции частоты лазерного излучения (\textit{chirped laser pulse}). Это означает, что разность частот $\Omega = \Omega(t) = \omega_1(t)-\omega_2(t)$ может зависеть от времени. Скорость решетки тогда определяется как $\xi = \xi_0 + 2\beta t/q$. Изменение скорости оптической решетки может быть использовано для увеличения эффектов, связанных с передачей энергии и импульса газу \cite{GraulChirp}.

Анализ встречных лазерных пучков строится на предположении, что амплитуда волн отдельно зависит и от временной, и от пространственной радиальной координаты. С точки зрения практического применения, во многих случаях амплитуды поперечного электрического и магнитного полей импульсного излучения описываются гауссовым профилем ($TEM_{00}$ мода), и электромагнитное поле не имеет компонеты, направленной вдоль направление распространения волны. В численных исследованиях обычно используется либо прямоугольная форма сигнала, либо гауссов профиль, более обоснованный с точки зрения экспериментального применения \cite{Graul2013},\cite{Coppendale2011}. В первом случае зависимость интенсивности от времени и радиальной координаты может быть описана константой - амплитуды полей неизменны за все время взаимодействия, во-втором профиль интенсивности описывается в виде (\ref{eq:IntensityGaus}) \cite{Diels2006}. Энергия лазерного импульса обозначается как $\Sigma$, $D$ - диаметр луча, $\tau$ - FWHM ширина импульса.

\begin{equation}\label{eq:IntensityGaus}
I_i(r,t) = \Sigma_i \ \left( \frac{4 \ ln2}{\pi}\right) ^{3/2} \frac{1}{D_i^2\tau_i} exp(-4 \ ln2\frac{r^2}{D_i^2}) \ exp(-4 \ ln2 \frac{t^2}{\tau_i^2})
\end{equation}

\subsection{Метод прямого статистического моделирования}

\hspace{0.5cm} Моделирование сложных аэродинамических процессов является крайне актуальной задачей нашего времени. Экспериментальные методы изучения разреженных течений обладают существенными техническими и технологическими недостатками и крайне затратны. Разработка современных аэрокосмических аппаратов, исследование неравновесных гиперзвуковых течений и взаимодействие струй двигателей, изучение конструкций спутников и космических станций, исследование явления оптического захвата газа - это далеко не полный список приложений, требующих новых, эффективных подходов, обладающих высокой точностью и широкой областью применимости. Метод \textit{прямого статистического моделирования} (ПСМ) является многообещающим кандидатом и универсальным методом подобных расчетов, основанный на численном решение уравнения Больцмана.

\subsubsection{Уравнение Больцмана и схема мажорантной частоты}

\hspace{0.5cm} Газ, для которого справедливо условие малости занимаемого молекулами объема относительного полного объема газа, называется разреженным. Если $n$ - число молекул газа в единице объема и $d$ - эффективный диаметр молекулы, $L$ - характерный размер течения, то условие разреженности можно записать как $nd^3 \rightarrow 0$ и $nLd^2 \leqslant const$. Воздух остается разреженным вплоть до сотни атмосфер, в соответствии с этими критериями. В таких газах осуществляются в основном лишь парные столкновения молекул, в то время как столкновения трех и более частиц маловероятны и ими пренебрегают \cite{Kogan}.

Одночастичная функция распределения $f_1(\vec{r},\vec{p},t)$ описывает плотность вероятности для одной выбранной частицы газа, для двух частиц эта плотность определяется двухчастичной функцией распределения $f_2(\vec{r_1},\vec{p_1},\vec{r_2},\vec{p_2},t)$. Для разреженного газа выполняется гипотеза о молекулярном хаосе \eqref{eq:MolecChaos}.

\begin{equation}
\label{eq:MolecChaos}
f_2(\vec{r_1},\vec{p_1},\vec{r_2},\vec{p_2},t) = f_1(\vec{r_1},\vec{p_1},t) f_1(\vec{r_2},\vec{p_2},t)
\end{equation}

Другими словами, гипотеза о молекулярном хаосе говорит о независимости нахождения двух частиц в некоторых точках $6N$-мерного фазового пространства.

Одночастичная функция распределения должна удовлетворять условию нормировки на плотность частиц $<f_1(\vec{r},\vec{p},t)>_{\vec{p}} \  = n(\vec{r},t)$ и описывается нелинейным интегро-дифференциальным уравнением Больцмана \eqref{eq:BoltzmanEq}. Правая часть уравнения Больцмана описывает источниковую часть уравнения переноса и называется интегралом столкновений.

\begin{equation}
\label{eq:BoltzmanEq}
\frac{\partial f}{\partial t} + \frac{\vec{p}}{m} \nabla_{\vec{r}} f + \frac{\vec{F}}{m} \nabla_{\vec{p}} f = \int (f' f'_1 - f f_1) \vec{g} b \ db d\varepsilon d\vec{v}_1
\end{equation}

Здесь $\vec{g} = \vec{v_1} - \vec{v}$ - относительная скорость сталкивающихся молекул, $b$ и $\varepsilon$ - параметры столкновения (прицельное расстояние и угол направления относительно $\vec{g}$), штрихи обозначают функции распределения после столкновения. На газ может действовать внешняя сила $\vec{F}$, в задаче о взаимодействии газа и оптической решетки равная $\vec{F} = -\vec{\nabla} U(\vec{r},t) = \frac{\alpha}{2} \vec{\nabla}  E^2(\vec{r},t)$.

Столкновение молекул не меняет модуля относительной скорости $\vec{g}$, а изменение её направления определяется величиной дифференциального сечения рассеяния $\sigma$. Полное эффективное сечение рассеяние определяется как \eqref{eq:ScatterFull}.

\begin{equation}
\label{eq:ScatterFull}
\sigma_{tot} = \int_{0}^{\Omega_{max}} \sigma d\Omega = \int b \ db d\varepsilon
\end{equation}

Частота столкновения молекул определяется как $\nu = n \overline{\sigma_{tot} g}$, а средняя длина свободного пробега $\lambda = \overline{v} / \nu$ определяется средней тепловой скоростью молекул газа. Число столкновений пропорционально квадрату плотности и равно $N_c = n^2 \overline{\sigma_{tot}g}/2$.

Метод прямого статистического моделирования является приближенным стохастическим подходом к решению кинетического уравнения. Трудоемкость алгоритма моделирования состоит в том, что для вычисления вероятностей и проведения этапа столкновений количество операций пропорционально квадрату числа частиц. Для построения линейного по сложности алгоритма применяется схема мажорантной частоты, впервые предложенный в \cite{IvanRogas}. В настоящей работе используется приближенный по времени метод подсчета величины мажорантной частоты.

В схеме с расщеплением перенос частиц происходит не после каждого столкновения, а в дискретные моменты времени с шагом в $\Delta t$. Тогда, так как внутри одной вычислительной ячейки число частиц не меняется за время шага, столкновения могут выполняться в каждой ячейке независимо от других. Мажоратная частота столкновений для $j$-ой ячейки записывается как

\begin{equation}
\nu_j = \frac{N_j (N_j - 1)}{2} \frac{(g \sigma_{tot})_j^{max}}{d_k}
\end{equation}

Для определения количества столкновений$N_c$ за временной шаг случайным образом разыгрывается случайная величина $\tau_{col}$, имеющая экспоненциальное распределение $\nu_j exp(-\nu_j \tau_{col})$, а затем определяется число столкновений $N_c$ так, что $\sum\limits_{k\leq N_p} \tau_{col,k} \leq \Delta t$ \ и \ $\sum\limits_{k \leq N_p+1} \tau_{col,k} > \Delta t$.


\subsubsection{Описание метода ПСМ}

\hspace{0.5cm} Применимость методов, использующихся для решения задач газодинамики, определяется безразмерным параметром - числом Кнудсена $Kn$. Оно определяется как соотношение между длиной свободного пробега $\lambda$ к характерному размеру задачи $L$. При числах Кнудсена, соответствующих континуальному режиму ($Kn < 0.01$) можно принебречь внутренней стуктурой газа и использовать уравнение Навье-Стокса для описания течения. Уравнение Навье-Стокса, основное уравнение гидродинамики, выводится из нелинейного кинетического уравнения Больцмана методом Чепмена-Энскога в приближении малого отклонения функции распределения скоростей от равновесной. Данных подход, с учетом дополнительных условий, можно расширить до $Kn < 0.1$. При числе Кнудсена $Kn \geqslant 0.1$ наступает переходный режим и уравнение Навье-Стокса заведомо не применимо, так как распределение скоростей частиц становится существенно неравновесным. Поэтому при больших числах Кнудсена необходимо рассматривать уравнение Больцмана.

Одним из методов численного решения уравнения Больцмана является прямое численное интегрирование \cite{Aristov},\cite{Tcheremissine}, обладающий, однако, рядом существенных недостатков: существенное повышение трудоемкости метода при увеличении размерности задачи, сложность обобщения метода при наличии внутренних степеней свободы и химических реакций течений и др.

В работе \cite{Kogan},\cite{Grad} описан метод Греда решения задач газодинамики, построенный на полученной из кинетического уравнения системе 13 моментных уравнений. Было показано \cite{Holway}, что метод Греда имеет важные недостатки. В дальнейшем были разработаны подходы на большем числе уравнений.

Метод прямого статистического моделирования \cite{Bird} является основным инструментом исследования сложных, многомерных течений разреженного газа. Это стохастический метод численного решения уравнения Больцмана при конечном числе Кнудсена. Основными преимуществами метода являются:

\begin{itemize}
\item Простота перехода от одномерной постановки задачи к двух- и трехмерным течениям,
\item эффективность использования метода на современных параллельных архитектурах CPU и GPU,
\item возможности использования различных подходов моделирования взаимодействия частиц газа, структуры их внутренних степеней свободы, химических реакций.
\end{itemize}

В задачах, решаемых методом ПСМ, обычно огромное число реальных молекул ($10^{10} - 10^{15}$) представляется меньшим числом ($\sim 10^5 - 10^9$) пробных частиц. Расчетная область разбивается на ячейки с размерами, много меньшими локальной длины свободного пробега молекул. Моделирование течения за дискретный шаг времени $\Delta t$ можно разбить на два этапа:

\begin{itemize}
\item Перемещение частиц, определяемое их текущими скоростями, на расстояние, которое бы они пролетели за шаг $\Delta t$. На этом этапе также проводится определение, взаимодействует ли частица с твердой поверхностью - если да, то частица изменяет направление скорости и отражается. За очередной временной шаг скорость частицы может меняется, например, из-за движения частицы во внешнем потенциале.
\item Межмолекулярные столкновения частиц, находящихся в одной вычислительной ячейке. Пара сталкивающихся частиц выбирается случайно, затем проверяется возможность столкновения этих частиц. Вероятность столкновения определяется относительной скоростью частиц, но не их положением внутри ячейки. Число выбираемых пар может определяться различными методами, однако в любом случае оно пропорционально квадрату числа частиц в ячейке.
\end{itemize}

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

Уменьшение числа Кнудсена приводит к существенному увеличению требуемых вычислительных ресурсов расчета. Поэтому исследование течений в континуальном и около-континуальном режимах в прошлом являлось непосильной задачей. Однако, с развитием произодительности современных вычислительных систем, в особенности параллельных методов с использованием многоядерных CPU и графических карт \cite{Vashen1},\cite{Vashen2}, область применения метода ПСМ расширилась даже на континуальные течения.

\subsubsection{Алгоритм метода ПСМ}

\hspace{0.5cm} Алгоритм метода ПСМ, используемый в частности в настоящей работе в программном комплексе \textit{SMILE++}, принципиально состоит из двух шагов - перемещение частиц и моделирование их столкновений. На практике, однако, данные два основные этапа состоят из большего числа подэтапов, полностью описывающих все необходимые процедуры расчета. На рис. \ref{fig:AlgoSmile} представлена диаграмма с базовыми этапами алгоритма метода ПСМ. Некоторые из подэтапов могут пропускаться и не использоваться в ряде задач, например, если область расчетов замкнута и внешнее течение, традиционно используемое для расчетов задачи обтекания тел, отсутствует.

\begin{figure}[h]
\centering
\includegraphics[width=400px]{AlgoSmile}
\caption{Диаграмма этапов алгоритма метода ПСМ.}
\label{fig:AlgoSmile}
\end{figure}

Начальная геометрия и разбиение расчетной области, характеристики газового течения, свойства внешних полей и другие стартовые данные задаются до начала расчета на этапе \textbf{Reader}. Затем на этапе \textbf{Gen} происходит заселение расчетной области, а также генерация частиц в случае наличия внешнего потока в расчетную область. Этап перемещения частиц \textbf{Move} описывает сдвиг пробных молекул за заданный временной шаг, в том числе удаление вылетевших за пределы расчетной области частиц. Если расчет производится в параллельном режиме, с использованием нескольких вычислительных потоков, то методы этапа \textbf{Paralleler} производят обмен частицами между соответствующими потоками. В целом, \textbf{Paraller} ответственнен за все процедуры, связанные с обменом данными между потоками. Далее, этапы \textbf{Collide} и \textbf{Chem} совершают вычисление столкновений между частицами и производят химические реакции, если такие подразумеваеются задачей. Сбор статистических данных и обработка полей течения проходит на этапе \textbf{ADC}. Расчетный цикл повторяется либо по достижению стационарного распределения полей течения в исследуемой области, либо при необходимом количестве итераций и требуемом физическом времени в эксперименте.

Задачи, решаемые методом ПСМ, обычно формулируются как равномерное течение, взаимодействующее с заданной геометрией и меняющее свою структуру. На рис. \ref{fig:SmileExample} продемонстрирована эволюция течения со скоростью 7500 м/с на высоте 90 км при обтекании цилиндра.

\begin{figure}[h]
\centering
\includegraphics[width=400px]{SmileExample}
\caption{Эволюция течения к стационарному состоянию при обтекании цилиндра.}
\label{fig:SmileExample}
\end{figure}

В подавляющем большинстве случаев использование однопоточного режима расчетов невозможно ввиду высокой требовательности метода к величине вычислительных ресурсов и памяти. В прикладных задачах метод ПСМ дополняется алгоритмами параллелизации, например, разбиением расчетной области на отдельные подзадачи \cite{IvanMark} или разделением данных между процессорами \cite{Grishin}.

\subsubsection{Параметры экспериментов и расчетной области}

\hspace{0.5cm} В каждом численном эксперименте задается набор параметров, определяющих геометрию и разбиение моделируемой области, свойства газа и излучения, вычислительные параметры расчета, например, количество используемых параллельных потоков. Изменение каждого из свойств приводит к изменению расчета, и поэтому описание задаваемых характеристик является неотъемлемым пунктом постановки задачи.

\textsl{Расчетная область.} Моделируемая область взаимодействия представляет собой двумерную прямоугольную область размером $d \times (d\cdot NCell)$ м$^2$, где $d$ - размер одной квадратной вычислительной ячейки, длина области больше ширины в $NCell$ раз - количество моделируемых ячеек вдоль оси распространения излучения. Размер $d$ определяется длиной волны излучения $\lambda$ и составляет $\lambda / 100$.

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

Число ячеек $NCell$ обычно выбирается кратным 100, таким образом моделируется целое число волн. Размер коротких областей составляет до 200$\lambda$, длиных достигает нескольких тысяч, например 2000$\lambda$. Для длины волны 532 нм размер такой области чуть более 1 мм.

\textsl{Характеристики излучения.} Каждый из двух источников излучения характеризуется следующими параметрами:

\begin{itemize} 
\itemsep0em 
\item Длина волны излучения [м] - также определяет производные величину (частоту, волновой вектор)
\item Интенсивность [Вт/см$^2$] - эквивалентно энергии или мощности пучка при заданной ширине импульса
\item Ширина импульса [с] - ширина гауссова профиля сигнала
\item Центр импульса [с] - момент времени, при котором интенсивность излучения максимальна
\item Диаметр пучка [м] - размер пучка, совместно с энергией и шириной импульса определяет интенсивность излучения
\end{itemize}

Длина волны излучения во всех экспериментах со стоящей решеткой устанавливалась равной 532 нм. Интенсивность варьировалась в диапазоне от $2.5 \cdot 10^{15}$ до $10^{17}$ Вт/см$^2$. Центр импульса и ширина импульса использовались в задачах с медленным (адиабатическим) включением поля, обычно ширина импульса $\tau$ равнялась 5 нс. Диаметр пучка был равен 45 $\mu$м.

\textsl{Моделируемые частицы.} В методах прямого статистического моделирования большое число реальных частиц заменяется на одну моделируемую, соотношение $MolRatio$ между числом реальных и моделируемых может составлять $10^{-13} - 10^{-15}$, то есть одна пробная частица представляет собой 1/$MolRatio$ реальных. В настоящих экспериментах это число равнялось $10^{-5}$, область размером $L = 2000\lambda$ заселялась более чем $1.4 \cdot 10^8$ моделируемыми частицами. Такое число модельных частиц соответствует уровню статистических флуктуаций порядка $0.1\%$. Для работы с таким количеством объектов использовалось от $1\times 8$ до $15 \times 8$ параллельных потоков, один шаг по времени $\Delta t = 10^{-12}\div10^{-13}$ с моделировался не более 6 секунд (для большой моделируемой области размером около 1 см).

Свойства частиц газа, используемых при моделировании, представлены в таблице. Модель твердых сфер, приближенно описывающая молекулы газа как сферические частицы с диаметром $d = d_r (c_{ref} / c_r)^\nu$, зависящем от относительной скорости при соударении $c_r$ (показатель вязкости газа равен 0.81 и 0.84 для агрона и метана соответственно), используется для моделирования межмолекулярных столкновений. Модель Ларсена-Борнакке \cite{Borgnakke} применялась для моделирования обмена энергией между трансляционными и внутренними степенями свободы.

\begin{center}
    \begin{tabular}{ | l | l | l | p{4cm} |}
    \hline
    Газ & $m$, Масса [$10^{-26}$ кг] & $\alpha$, Поляризуемость [$10^{-40}$ см$^2$ V$^{-1}$] & $d_r$, Диаметр [$10^{-10}$ м]\\ \hline
    Ar & 6.63 & 1.83 \cite{BirdMFP} & 4.17 \\ \hline
    CH$_4$ & 2.656  & 2.88 \cite{BirdMFP} & 4.83 \\ 
    \hline
    \end{tabular}
\end{center}

\subsection{Макропараметры течения}

\hspace{0.5cm} Получение информации о пространственном распределении макропараметров течения - среднего вектора скорости, температуры, давления и многих других - является существенным этапом обработки и анализа результатов, полученных методом ПСМ. В то время, как традиционно накопление таких данных осуществляется после установления стационарного течения в некоторых задачах, иногда необходимо знать значения осредненных по пространственной ячейке моментов функции распределения на каждом временном шаге.

Интересующие нас макропараметры могут быть получены осреднением некоторой функции скорости по набору частиц \cite{Havyland}. Мы будем говорить о выборке макропараметров по ячейкам - вычисляется значение некоторой величины по всем частицам в ячейке. Если $\phi(v)$ - интересующий нас функционал скорости, то среднее значение $\overline{\phi}$ может быть получено как:

\begin{equation}
n \overline{\phi} = \int \phi f(r,v) dv
\end{equation}

Для набора $N_c$ частиц оценка момента скорости $\phi$ очевидно записывается как $n \overline{\phi} = \sum ^{N_c} \phi F_n / V_{c}$, где $F_n$ - число реальных молекул, приходящихся на одну модельную, а $V_c$ - объем ячейки.

Средняя скорость молекул - первый момент скорости, может быть получен простым осреднением по частицам в ячейке $v_0 = \left\langle v \right\rangle$. Выражение также справедливо для $x,y,z$ компонент вектора скорости. Для конечного набора частиц средняя скорость является средним арифметическим скоростей частиц в ячейке:

\begin{equation}
v_{0,x} = \left\langle v_x \right\rangle = \frac{\sum\limits_j v_x^j}{\sum\limits_j 1} = \frac{\sum\limits_j v_x^j}{N_c}
\end{equation}

Если необходимо получить значение величины средней скорости газа в некотором наборе ячеек, достаточно знать данные о числе частиц и средней скорости в каждой $i$-ой отдельной ячейке. Тогда:

\begin{equation}
\left\langle v_{0,x} \right\rangle = \frac{\sum\limits_i v_{0,x}^i N_c^i}{\sum\limits_i N_c^i}
\end{equation}

Температура поступательного движения молекул может быть записана функция первого и второго моментов скорости. Определив термальную скорости $v'$ как разность между полной величиной скорости и средней скоростью движения газа как целого $v' = v - v_0$, температура может быть записана как

\begin{equation}
\frac{3}{2} k_B T_{tr} = \frac{m}{2} \left\langle v'^2 \right\rangle 
\end{equation}

Можно также определить $x$-компоненту поступательной температуры (вдоль выделенного направления действия градиентной силы) как ${T_{x,tr} = \frac{m}{k} \left\langle v'^2_x \right\rangle}$. Средний квадрат термальной скорости можно удобно получить, зная среднее значение квадрата полной скорости и среднюю скорость газа: $\left\langle v'^2_x \right\rangle = \left\langle v_x^2 \right\rangle - \left\langle v_{x} \right\rangle ^2$.

Аналогично вопросу о средней скорости течения в наборе ячеек, порой необходимо знать среднее значение температуры газа не в одной ячейке, а в нескольких. В этом случае осреднение проводится по всему объему, либо может быть выражено через средние значения скорости, температуры и плотности в отдельных ячейках \eqref{eq:AvgTemp}.

\begin{equation}
\label{eq:AvgTemp}
T_{x,tr} = \frac{m}{k_B} \left( \frac{\sum\limits_{cells} \frac{k_B}{m} T_x^j + N_j (v_{0,x}^j)^2}{\sum\limits_{cells} N_j} - \left\langle v_{0,x} \right\rangle^2 \right) 
\end{equation}

Вычисление температур $T_{rot}$, $T_{vib}$, определяющих вращательные и колебательные состояния молекул, вычисляются аналогично средней скорости частиц в объеме.