第4回「粒子の動かし方」と「加速度の求め方」について | 粒子法入門

今回のコラムでは、粒子法で基本となる粒子の動かし方についてお話します。粒子の移動量（距離）は「距離=速度×時間※1」で表されますので速度の情報が必要です。その速度を求める（更新する）のに必要な加速度についてもお話します。

1. 粒子の動かし方

簡単に書くと、次のような手順で粒子を動かします。

[0] 初期時刻における全ての粒子の座標※2と速度$ \boldsymbol{u}$※3を定める。

[1] ナビエ・ストークス方程式により粒子の加速度$ \boldsymbol{a}$を求める。

[2] 求めた加速度に従って粒子の速度を更新（加速、減速、等速）させる。

[3] 更新した速度で粒子を少し移動させる。

[4] 手順[1]に戻って[1]～[3]を繰り返す※４。

いかがでしょうか。そんなに難しくないですよね。

2. 粒子の加速度の求め方

前節の手順[1]で、「ナビエ・ストークス方程式により粒子の加速度を求める」と書きました。ナビエ・ストークス方程式は、流体の運動方程式です。まずは運動方程式から説明します。運動方程式は、ある質量の物体に力$ \boldsymbol{F}$が加わったときの物体の加速度$ \boldsymbol{a}$を表す式であり、次式のように表されます。

$$ \boldsymbol{m} \boldsymbol{a}=\boldsymbol{F} $$

この式の両辺を$\boldsymbol{m}$で割ると

$$ \boldsymbol{a}=\frac{\boldsymbol{F}}{\boldsymbol{m}} $$

となります。この式は、物体に働く加速度は物体に働く力$ \boldsymbol{F}$に比例し、質量$ \boldsymbol{m}$に反比例するということを表します。式(2)の右辺（物体に働く力$ \boldsymbol{F}$を質量$ \boldsymbol{m}$で割った値）を計算すれば、加速度$ \boldsymbol{a}$が求まるということです。流体の場合も同様です。ある流体の塊に働く力をその流体の塊の質量で割ることで、次式のように加速度が求まります。

$$ \boldsymbol{a}=-\frac{1}{\rho} \nabla P+\nu \nabla^{2} \boldsymbol{u}+\boldsymbol{g} $$

この式(3)がナビエ・ストークス方程式です。様々な記号が使われていますが、拒絶反応を起さないでください。この式を簡単に説明すると、「流体の場合にはいくつかの力※５（この場合３つの力）が働くため、右辺のように複数の加速度成分の和として加速度は書ける」ということを表しています。粒子法では、式(3)の右辺を計算することで流体粒子の加速度$ \boldsymbol{a}$を求めているのです。式(3)の右辺の各項の意味については補足※６を参照ください。各項の具体的な計算方法については参考文献1)に詳しく書かれていますので、興味がある人はぜひご参照ください。

参考文献：

越塚誠一、柴田和也、室谷浩平、”粒子法入門”、丸善出版、2014年6月25日

補足：

※1：ここで時間は、時間刻みの幅を意味します。流速や粒子の大きさにもよりますが、例えば $ \Delta\boldsymbol{t}=\boldsymbol{0.001}$など十分小さな値に設定します。

※2：粒子の初期座標は、通常は粒子が等間隔に並ぶように定めます。この間隔は、初期粒子間距離と呼ばれ、粒子の大きさ（空間解像度）を表します。

※3：粒子の初速度は、解く問題に合わせて設定します。例えば、初期時刻に流体が静止していれば、全ての粒子の初期速度を$\boldsymbol{0m/s}$に設定します。

※4：圧力を陰的に求める場合は、若干手順が増えます。

※5：表面張力や他の外力が加わる場合は、さらに項の数が増えます。

※6：式(3)の右辺の第１項の$-\frac{1}{\rho} \nabla P$は圧力による加速度（圧力$ P$の傾きに比例する力。マイナスがついているのは、圧力が高い方向ではなく低い方向に力が働くためである。）を表します。右辺第２項の$\nu \nabla^{2} \boldsymbol{u}$は粘性による加速度の項（流体の粘り気により運動量が拡散される項。$\nu$は動粘性係数。）を表す。右辺第３項の $\boldsymbol{g}$は重力項（鉛直下向きに約$\boldsymbol{9.81m/s}$の大きさを持つベクトル）です。