DSMC法の特徴の一つは，物理空間をセルと呼ばれる微細な空間要素に分割し，同一のセルに存在する分子同士で分子間衝突の計算を実行するというものである．そのため，セルを非常に小さくする必要があり，理想的には，セル一辺は平均自由行程以下にせねばならない（計算対象によっては平均自由行程の1/3以下が必要との報告もある）．密度が非常に低い宇宙空間なら，この条件は大して問題にはならないであろうが，それが例えば大気圧状態であるなら，平均自由行程は約0.066 $\mu$mとなり，1mmの長さには15000個のセルが必要で，1mm$^{2}$の面積なら2.3×10$^{8}$個となる．分子数としては，さらにその数倍が必要になるので，小型のPCでは扱えないし，計算時間も膨大なものになる．もちろん，セルを平均自由行程より大きくとって計算してもその結果が全く使い物にならないというわけではないが，厳密さに欠けるものになるのは避けられない．
私が大学の講義で平衡状態の分子速度の分布関数（マクスウェル分布）を解説しているとき，このマクスウェル分布をうまく利用してセルを大きくできないだろうかと考えたのが，このU_system を考案するきっかけである．DSMC法の通常の分子間衝突では，同じセルにある分子は，セル内で別々の位置にあるにもかかわらず，全て，同じ位置（例えばセル中心）にあると考えて衝突計算にかける．セルが平均自由行程を超えて大きくなるほど，セル内で遠く離れた位置にあって本来衝突するはずのない分子が衝突計算され，結果にひずみを作ることになる．そこで，もし，二つの分子を衝突させる際，元に位置に応じて分子の速度成分に修正を加えたらどうだろうか．その際，平衡状態の分布関数を利用できないだろうか．これが，U_systemの根底にある考え方です．
以下ではU_systemの説明をするにあたって，U_system用の例題プログラムを使用するが，これは，5.1節で解説したDSMC法の基礎プログラムにU_system部分を追加したものです．したがって，あらかじめ5.1節とそれに付属する基礎プログラムを十分に理解している必要がある．U_systemの説明中，例題プログラムの「〇行～〇行参照」と表記するが，5.1節で行ったような丁寧な説明でなく，大雑把なものに留めており，また，必ずしも解説上の変数記号とFortranプログラムの変数記号が一致しているとは限らないので，細心の注意を払って読んでほしい．
　先に述べたように，従来の衝突計算法では，セルが大きくなるとそれに比例して二つの衝突分子の位置は離れるため，シミュレーションが忠実に行われなくなる．本来，衝突する2つの分子は同じ位置で選ばれるのが最も好ましいが，有限個の分子しか用いないシミュレーションでは，たとえセルを小さくしたとしても全く同じ位置で2つの分子を見つけるのは困難である．
　U_systemでは，衝突計算直前に一方の分子Pを，あたかも，もう一方の分子Qの位置に一致させたかのように分子速度の補正を行う．いま，流れ場中の分子は，全て，ある温度と流れ速度を持つ局所平衡状態の速度分布関数（マクスウェル分布）に従うものと仮定する．なお，以下では，2つの分子位置における速度分布関数をマクスウェル分布で仮定して補正式を導いているが，得られた補正式自体は，2つの位置の分布関数が必ずしもマクスウェル分布でなくとも有効である（例えば，マクスウェル分布式に現れる指数関数が別の関数に変わってもよい）．唯一必要な条件は，「分布中心が，流速に応じて左右に平行移動し，分子の熱運動速度が，温度の平方根に比例するような分布形状」というものである．大きいセルといっても同じセル内にある2つの分子であるから，マクスウェル分布でないとしても，上記の条件を満たすのはそれほど困難なことではない．
　さて，分子Pの位置の分布関数を$f_1$，分子Qの位置の分布関数を$f_2$とし，分子Pの速度が$f_1$上の一点Aにあるものとしよう（Fig.1）．衝突計算に際して，このA点を分布関数$f_2$上のどこかに移せばよいことになるが，私が最も自然であると考えた方法は，$f_1$上のA点位置と$f_2$上の位置が相対的に同じになるよう移動させることである（B点）．このB点位置を求める計算は簡単である．いま，$x$方向速度成分$u$の分布関数について考えると，$f_1$と$f_2$は，次に示す式(1)および式(2)のように書ける．ここで，$R$は気体定数，$T_{x1}$, $T_{x2}$は，分子の並進運動に基づく温度3成分のうちの$x$方向速度成分により定義される温度（以下，$x$方向温度と呼ぶ）で，$U_1$，$U_2$は，流速の$x$成分である．

\begin{equation}
f_{1}=\frac{1}{\sqrt{2\pi R T_{x1}}}\exp\{-\frac{(u_{A}-U_{1})^2}{2 R T_{x1}}\}
\end{equation}

\begin{equation}
f_{2}=\frac{1}{\sqrt{2\pi R T_{x2}}}\exp\{-\frac{(u_{B}-U_{2})^2}{2 R T_{x2}}\}
\end{equation}

$f_1$および$f_2$をそれぞれ最大値で無次元化し，それらを等しいとおくと，次式となる．

\begin{equation}
\frac{(u_{A}-U_{1})^2}{T_{x1}}=\frac{(u_{B}-U_{2})^2}{T_{x2}}
\end{equation}

したがって，次式が得られる．

\begin{equation}
u_{B}=U_{2}+(u_{A}-U_{1})\sqrt{\frac{T_{x2}}{T_{x1}}}
\end{equation}

この修正された速度$u_B$（なお$u_A$は修正前の速度）を用いて，これ以降の計算は一般的な衝突計算法（Birdの計算法）で行う．結果として得られた衝突後の速度$u_C$ (C点) は，分子Pを元の位置に戻すのに伴い，次式を用いて再び修正される(D点)．

\begin{equation}
u_{D}=U_{1}+(u_{C}-U_{2})\sqrt{\frac{T_{x1}}{T_{x2}}}
\end{equation}

なお，U_systemという名称であるが，上記のように見掛け上，衝突ペアのうちの片方の分子は，一旦，衝突相手の分子位置に移動し，衝突後に元の位置に戻ることになるので，Uターンするという意味でUを用い，また，この操作をすべての衝突ペアに対して実行する仕組み全体をsystemと表現して，このように名付けている．なお，セル内で流速と温度に変化がなければ，一般的な衝突計算法と同じものになる．
　ここで，式(4)(5)を用いるためには，分子PとQそれぞれの位置における流速$U$と温度$T$のデータが必要である．セルの端のデータが分かっているとして，一次元問題では単なる線形補間を使い，二次元問題では双線形補間を，三次元問題では三線形補間を使うことになる．例題プログラムは軸対称問題であり，セルの４隅の値を用いた双線形補間となるので下記にその計算式を示す．
＜双線形補間＞ 二次元問題におけるセル内の任意位置の流速と温度の求め方
長方形セルの四隅の座標を ($x$, $y$), ($x+A$, $y$), ($x$, $y+B$), ($x+A$, $y+B$) とし，その点の温度を $T_1$, $T_2$, $T_3$, $T_4$ とすると， セル内の任意位置 ($x+a$, $y+b$) における温度は，

\begin{equation}
T=(1－b/B)\{(1－a/A) T_{1}+ (a/A) T_{2}\} + (b/B)\{(1－a/A) T_{3} + (a/A) T_{4}\}
\end{equation}

で求められる．流速についても同様である．
なお，三次元問題に適用する８隅データを用いた三線形補間の計算式は少し複雑になるが，練習として読者諸君が導くとよい．さて，最初の計算だけは，従来のDSMC法（Bird法）によって各セルの流速と温度を求める．DSMC法の長い繰返し計算の中で，従来の方法による計算は最初の１回だけでよく，以降は，U_systemによって流速と温度のデータが逐次更新される．なお，U_systemの計算１回に要する時間は，従来の計算法より増加するが，セルを大きくできる点を考慮すれば，計算全体での時間比較ではU_systemが有利となる．
　さて，分子速度の補正を行うと，各衝突分子ペアの持つ運動量とエネルギが衝突前後で厳密には保存されない．もちろん流れ場全体としてあるいは長時間にわたって平均をとれば，統計的ばらつきの範囲内で保存される．しかし，多くの計算例を検討した結果，一回の衝突において，衝突前後のエネルギに大きな差が生じると，特に，数個のセル範囲で衝撃波が生じるような変化の激しい流れ場では，結果にひずみを成長させる原因となることが判明した．そこで，衝突前後でエネルギが保存されるような改良を加えたところ，U_systemの計算安定化は大きく改善された(U_system-3)．しかし，運動量に関する保存が保証されてないとの指摘があり，また，計算条件によっては温度上昇が予想外に生じるといった事例も発生して，更なる改善が望まれていた．Sunら（参考文献7）は，分子間衝突を無次元速度空間で行うことを考え，速度を無次元化するためのスケール変数に流速と温度を用いるという概念を導入してDSMC解析を行っている（LCDSMC法）．その計算法は，基本的には，著者らのU_systemと同じ算法ではあるが，特筆すべきは，分子の運動量と並進エネルギを各セルごとに厳密保存するために，Pareschiら（参考文献8）が提案した補正方法を採用したことである．次に，その方法を説明する．
ある流速と温度が与えられた局所平衡状態の速度分布関数から，乱数を用いて，複数の分子の速度を再現した場合，再現された分子の速度から，逆に流速と温度を求めても，厳密には元に戻らない．これを厳密に元に戻すために，Pareschiらは，運動量とエネルギの厳密保存の方法を考案した．いま，この方法をU-systemに適用することを考える．衝突前と衝突後の分子速度成分を$u$, $v$, $w$および$u^{*}$, $v^{*}$, $w^{*}$ とし，これを，Pareschiらの補正法を使って運動量とエネルギを厳密保存して $u^{+}$, $v^{+}$, $w^{+}$ となるものとしよう．この場合，補正は式(6)のように，4つのパラメータ $\lambda_X$, $\lambda_Y$, $\lambda_Z$, $\tau$ によって行うことになる．

ここで， $\lambda_X$, $\lambda_Y$, $\lambda_Z$, $\tau$ を求めるために，次式(7)(8)に示す運動量保存式とエネルギ保存式を用いる．なお，mは分子質量，添え字iは，全ての分子について和を求めることを意味する．

                                         (7)
        　　　　　　　　　　　　　　　　　　　　　                        

　　　　　　(8)
　　　　　　　　　　　　　　　　　　　　　
結果として，$\tau$, $\lambda_X$, $\lambda_Y$, $\lambda_Z$ は下記の一連の式で求められる．

        　　　　　　　　　　　　　
                                             (9)

ここで，

本ブログに収められたU_systemは，上記に示した運動量とエネルギの厳密保存を施したものである．ちなみに，この厳密保存の方法は，南部法（修正南部法）に適用すれば南部法の欠点を補えるのではないかとも考えられる．なお，$\tau$ を求める際に，平方根の前の符号には必ず正を用いる必要がある．負であっても保存則は満足するが，負であると，補正前後の速度方向を一致させるために，修正用パラメータ $\lambda_X$, $\lambda_Y$, $\lambda_Z$ が非常に大きな値になってしまう（補正値 $\lambda_X$, $\lambda_Y$, $\lambda_Z$ と流速の絶対値がほぼ同じ大きさになる）．一方，正の場合，$\lambda_X$, $\lambda_Y$, $\lambda_Z$ は流速の1/100以下である．もともと式(6)は修正に使うものであり，修正前後で差異が小さいことを想定しているので，差異の小さい方を採用する．実際に適用した結果でも，平方根の前の符号に負を用いると，保存則は成立するがU_systemの効果は殆どなくなってしまう．
最後に，U_systemという名称であるが，そのUが，先に述べたようなUターンという意味ではなく，私UsamiのUから取ったのではないかと勘繰りを入れる人がいるかも知れない．下賤な私をよく知る人からすればもっともな話である．しかし，これに関しては下記と同様のことを言っておきたいと思う．日本の文豪「森鴎外」の小説「雁」に登場する絶世の美女「お玉」のその後について，僕（鴎外）の情人になったのではないかと訝しがる人に対して，小説の最後で鴎外は次のように述べている．「読者は無用の憶測をせぬが好い．」