Closed24

粒子法

ピン留めされたアイテム
ksttrksttr

参考文献

  1. 後藤仁志 ，粒子法 連続体・混相流・粒状体のための計算科学，森北出版（2018）
  2. 越塚誠一，柴田和也，室谷浩平，粒子法入門 流体シミュレーションの基礎から並列計算と可視化まで，丸善出版（2014）
ksttrksttr

粒子法とは

流体を粒子の軍団として離散化して数値計算したもの．

移流の効果を粒子の移動として表現するので，非線形項を陽に離散化して計算する必要がない．
そのため，非線形項を離散化したときに現れる数値誤差（高波数の振動とか）が抑えられる．

また，格子点によって離散化する方法と比べて，自由境界が激しく変形する場合にも簡単に適応できる．
例えば非線形波動の突っ立ちのように，格子点の方法では計算が破綻してしまう状況でも，粒子法では特に問題なく計算が続けられる．

ksttrksttr

流体の粒子の集団としての離散化

粒子方では，流体を仮想的な粒子の集団として離散化する．
この点は，Lagrange的な描像を採用していると表現できる．Lagrange的な描像での粒子（流体粒子またはLagrange粒子）は互いに相互作用しないが，粒子法では，粒子同士が相互作用する．

相互作用の具体系は，粒子法のバリエーション(SPH, MPS, ...)によって異なる．

SPH

SPH法では，まず，任意の物理量 \phi を重み w との畳み込み積分をすることで，一定の大きさを持つ量になます．
このようになますと，離散化ができる．さらに，この物理量は粒子によって運ばれていると表現する．
重み w がある範囲に広がっていることで，圧力項と拡散項の離散化によって，粒子同士に相互作用が生まれる．

MPS

MPS法での相互作用，および，離散化は，なかなか発見的なようである．
MPS法では，まず初めに粒子間の相互作用を定義する．

最近では，高精度化の試みに伴って，より論理的に導出されてきているらしい．

ksttrksttr

MPS法

まずは，MPS法を調べる．
とりあえず，お気持ち的な導出だけをして，とっとと実装までしてしまいたい．

MPS法では，物理量は粒子の上で定義される．
したがって，ある物理量 \phi のある位置 \boldsymbol{r} の時刻 t での値 \phi(\boldsymbol{r},t) を求めるには，その場所・時刻に粒子がいなければならない．

また，物理量の空間微分に関する量を計算するには，他の地点の粒子との差分を取る必要があるが，近傍粒子との差分だけでなく，影響範囲を定義しておいて，この影響範囲内にある粒子全てとの重み付きの差分をとる．

このような意図で，MPS法で登場する粒子には，カーネルによって定義される他粒子への影響範囲が存在する．

重み関数

この影響範囲は重み関数 w(r) で表される．
MPS法では，重み関数 w は以下のように与えれるのが普通らしい．

w(r) = \begin{cases} r_e/r - 1 & \mathrm{for} \quad 0 \leq r \leq r_e \\ 0 & \mathrm{otherwise} \end{cases}

r は粒子間の距離を表す．
ここで，r_e が影響範囲を決めていて，これより近い粒子が影響し，これよりも遠い粒子は影響しない．

粒子数密度

MPS法では，粒子数密度 n を定義する．
前にも述べた通り，全ての物理量は，粒子上でのみ定義される．

粒子を適当に番号づけて，i = 1, \dots, N とする．また，粒子 i の位置を \boldsymbol{r}_i とする．
粒子 i での粒子数密度 n_i は，

n_i = n(\boldsymbol{r}_i) \equiv \sum_{i\not=j} w(|\boldsymbol{r}_i - \boldsymbol{r}_j|)

と定義される．

流体の質量密度 \rho はこの粒子数密度に比例すると仮定する．

\rho(\boldsymbol{r}) \propto n(\boldsymbol{r})

（この気持ちはまだよくわかっていない．．．）

つまり，非圧縮流体では，この粒子数密度は時間によらず一定であるとする．

粒子数密度の基準値 n_0

密度が均一な非圧縮性流体の場合，粒子数密度 n は一様であり，時間に対して一定である．
したがって実用上では，初期時刻において，どこか一点だけ計算しておけば十分である．
このような粒子数密度を基準値と呼び，n_0 と表す．粒子の番号 0 に対する粒子数密度と同じ記号であるが，以下で混同することはない．

n_0 は，初期時刻において，影響半径内に均一に粒子があるような粒子を選び，その粒子に対して，

n_0 = \sum_{j\not=i^\prime} w(|\boldsymbol{r}_j - \boldsymbol{r}_{i^\prime}|)

を計算する．ここで，i^\prime は上記のように選ばれた粒子の番号である．

ksttrksttr

MPS法での微分演算子

Navier-Stokes方程式を解くために，微分演算子を離散化する必要がある．

格子がないので，粒子同士の差分として，微分を離散化する．

MPS法の離散化には2ステップある．

  1. 粒子 i と粒子 j との間での微分演算子の離散化を導入する．
  2. 1の離散化を全粒子にわたって重み付き平均をとる．

1による離散化を (\cdot)_{ij} と書く．2による重みつき平均を \langle \cdot \rangle_i と書き，以下で定義する．

\langle f \rangle_i \equiv \frac{\sum_{j\not=i} (f)_{ij} w(|\boldsymbol{r}_j - \boldsymbol{r}_i|)}{\sum_{j\not=i} w(|\boldsymbol{r}_j - \boldsymbol{r}_i|)} = \frac{1}{n_i}\sum_{j\not=i} (f)_{ij} w(|\boldsymbol{r}_j - \boldsymbol{r}_i|)
ksttrksttr

勾配

勾配 \nabla \phi の計算には，まず，物理量 \phi が粒子iの影響範囲内で一様かつ等方であると仮定する．
影響範囲とは，wの台，つまり，w > 0となる範囲であった．

これは，\phiの激しく変動せず，影響範囲を決める量r_eが十分小さい時に妥当な条件だと思われるが，厳密な条件は不明．

この仮定を満たす関数\phiに対して，勾配の離散化は以下のように行われる．
まず，ステップ1では，

(\nabla \phi)_{ij} = d\frac{\phi_j - \phi_i}{|\boldsymbol{r}_j - \boldsymbol{r}_i|} \frac{\boldsymbol{r}_j - \boldsymbol{r}_i}{|\boldsymbol{r}_j - \boldsymbol{r}_i|}

とする．ただし，d は空間次元である．d がつく理由はよくわからない．
\phi(r)として勾配を計算すると，\nabla \phi(r) = \phi^\prime(r) \boldsymbol{r}/rであって，dは出てこない．．．

ステップ2では，これを重み付き平均する．

\langle \nabla \phi \rangle_i = \frac{d}{n_i} \sum_{j\not=i} \frac{(\phi_j - \phi_i) (\boldsymbol{r}_j - \boldsymbol{r}_i) }{|\boldsymbol{r}_j - \boldsymbol{r}_i|^2} w(|\boldsymbol{r}_j - \boldsymbol{r}_i|)

発散

同じようにして，発散は，

\langle \nabla \cdot \boldsymbol{u} \rangle_i = \frac{d}{n_i} \sum_{j\not=i} \frac{(\boldsymbol{u}_j - \boldsymbol{u}_i) \cdot (\boldsymbol{r}_j - \boldsymbol{r}_i) }{|\boldsymbol{r}_j - \boldsymbol{r}_i|^2} w(|\boldsymbol{r}_j - \boldsymbol{r}_i|)

となる．

ksttrksttr

ラプラシアン

ラプラシアン \Delta \phi の計算も，これまでと同様に，\phi が影響範囲内で一様かつ等方であると仮定する．

この仮定のもとでは，粒子 i の位置を原点とした極座標表示で，ラプラシアンは，

\Delta \phi = \frac{1}{r^2} \frac{\partial}{\partial r} r^2 \frac{\partial \phi}{\partial r} = \frac{1}{r^2}\left(2r\frac{\partial \phi}{\partial r} + r^2 \frac{\partial^2 \phi}{\partial r^2} \right)

となる．これを離散化する．

これまでと同様にステップ1で，粒子 i と 粒子 j との間で，ラプラシアンを表現する．すなわち，

(\Delta \phi)_{ij} = \frac{1}{r_{ij}^2} \left(2r_{ij}\left(\frac{\partial \phi}{\partial r}\right)_{ij} + r_{ij}^2 \left(\frac{\partial^2 \phi}{\partial r^2}\right)_{ij} \right)

を構成する．

r_{ij}については，r_{ij} = |\boldsymbol{r}_j - \boldsymbol{r}_i| とする．

rによる微分を含む部分に関しては，以下のようにテイラー展開から構成する．

粒子 i の位置 \boldsymbol{r}_i 周りの \phi のテイラー展開は，

\phi(\boldsymbol{r}) = \phi(\boldsymbol{r}_i) + (r-r_i) \left.\frac{\partial \phi}{\partial r}\right|_{r=r_i} + \frac{1}{2} (r-r_i)^2 \left.\frac{\partial^2 \phi}{\partial r^2}\right|_{r=r_i} + \cdots

となる．ただし，ここでも \phi の一様等方性をつかって，\phi は r=|\boldsymbol{r}| のみの関数として計算している．

このテイラー展開を2次までで打ち切って，r=r_j，すなわち，粒子 j での \phi を評価すると，

\phi_j = \phi_i + (r_j - r_i) \left. \frac{\partial \phi}{\partial r}\right|_{r=r_i} + \frac{1}{2} (r_j-r_i)^2 \left.\frac{\partial^2 \phi}{\partial r^2}\right|_{r=r_i}

すなわち，

2(r_j - r_i) \left. \frac{\partial \phi}{\partial r}\right|_{r=r_i} + (r_j-r_i)^2 \left.\frac{\partial^2 \phi}{\partial r^2}\right|_{r=r_i} = 2(\phi_j - \phi_i)

ここで，r_j - r_i \rightarrow |\boldsymbol{r}_j - \boldsymbol{r}_i| として辺々 |\boldsymbol{r}_j - \boldsymbol{r}_i|^2 で割ると，ラプラシアンの粒子i,jでの表現を得る．

(\Delta \phi)_{ij} = \frac{2(\phi_j - \phi_i)}{|\boldsymbol{r}_j - \boldsymbol{r}_i|^2}

ステップ2によって，粒子i でのラプラシアンの表現は，

\langle \Delta \phi \rangle_i = \frac{2d}{n_i} \sum_{j\not=i} \frac{(\phi_j - \phi_i)}{|\boldsymbol{r}_j - \boldsymbol{r}_i|^2} w(|\boldsymbol{r}_j - \boldsymbol{r}_i|)

ここでも d をつける．

MPS法では，さらに，

\frac{1}{n_i} \sum_{j\not=i} \frac{(\phi_j - \phi_i)}{|\boldsymbol{r}_j - \boldsymbol{r}_i|^2} w(|\boldsymbol{r}_j - \boldsymbol{r}_i|) \rightarrow \frac{\sum_{j\not=i}(\phi_j - \phi_i)w(|\boldsymbol{r}_j - \boldsymbol{r}_i|)}{\sum_{j\not=i} |\boldsymbol{r}_j - \boldsymbol{r}_i|^2 w(|\boldsymbol{r}_j - \boldsymbol{r}_i|)}

と置き換える． 非圧縮流体であれば，n_i \simeq n_0 であるから，結局，

\begin{align*} \langle \Delta \phi \rangle_i &= \frac{2d}{\lambda n_0} \sum_{j\not=i} (\phi_j - \phi_i)w(|\boldsymbol{r}_j - \boldsymbol{r}_i|) \\ \lambda &\equiv \frac{1}{n_0} \sum_{j\not=i} |\boldsymbol{r}_j - \boldsymbol{r}_i|^2 w(|\boldsymbol{r}_j - \boldsymbol{r}_i|) \end{align*}

がMPS法でのラプラシアンの離散化である．

ksttrksttr

Navier-Stokes方程式の時間発展

Navier-Stokes方程式の時間発展は，射影法を用いて，2段階で行われる．
ここで，Navier-Stokes方程式は，

\begin{align*} \frac{D \boldsymbol{u}}{D t} &= -\nabla \left(\frac{p}{\rho}\right) + \nu \Delta \boldsymbol{u} + \boldsymbol{g} \\ \nabla \cdot \boldsymbol{u} &= 0 \end{align*}

である．ここで，\rhoは流体の質量密度，\nuは粘性係数，\boldsymbol{g}は重力加速度ベクトルである．

ここで，粒子法における時間発展する変数は，速度場 \boldsymbol{u} と 圧力場 p に加え，粒子数密度 n があり，非圧縮性の要請から n を時間に関して一定に保つように計算する．

ksttrksttr

Navier-Stokes方程式の物質微分を，

\frac{D\boldsymbol{u}}{Dt} = \frac{\boldsymbol{u}^{k+1} - \boldsymbol{u}^{k}}{\Delta t}

と時間方向に離散化すると，粒子法での時間発展は，

\boldsymbol{u}^{k+1} = \boldsymbol{u}^k + (\nu\Delta \boldsymbol{u})^k\Delta t + \boldsymbol{g} \Delta t - \left(\nabla\left(\frac{p}{\rho}\right)\right)^{k+1} \Delta t

となる．
粒子法の特徴として，物質微分をそのまま差分で表し，移流の効果を粒子位置の更新で表す．

非圧縮性を保つために，射影法によって時間発展を2段階で行う．
ここで，上記の離散化したNavier-Stokes方程式の圧力勾配項がk+1となっているのは，射影法による時間発展の第2段階で陰的に解くことによる．

ksttrksttr

第1ステップ

まず，圧力勾配項を除いた項による時間発展を行う．すなわち，

\boldsymbol{u}^* = \boldsymbol{u}^k + (\nu\Delta \boldsymbol{u})^k\Delta t + \boldsymbol{g} \Delta t

として，さらに，粒子位置の更新も行う．

\boldsymbol{r}^* = \boldsymbol{r}^k + \boldsymbol{u}^* \Delta t

第2ステップ

第2ステップでは，残りの圧力勾配項による時間発展を行うのであるが，このとき，非圧縮性 \nabla\cdot \boldsymbol{u}^{k+1} = 0 を保つために，圧力場を陰的に求める．

今，離散化したNavier-Stokes方程式は，

\boldsymbol{u}^{k+1} = \boldsymbol{u}^* - \left(\nabla\left(\frac{p}{\rho}\right)\right)^{k+1} \Delta t

と書ける．非圧縮性を保つためには，この両辺の発散をとって，

0 = \nabla\cdot\boldsymbol{u}^* - \Delta \left( \frac{p}{\rho} \right)^{k+1} \Delta t

を満たす必要がある．これは，圧力場 p に関するPoisson方程式になっている．

ちなみに，このような性質を持つ圧力勾配 \nabla p^{k+1} が存在することは，Helmholtzの分解定理が保証している．すなわち，任意のベクトル場 \boldsymbol{u} は，

\boldsymbol{u} = \boldsymbol{u}^s + \nabla \phi, \quad \mathrm{for} \quad \nabla\cdot\boldsymbol{u}^s = 0

と表すことができる．ここで，\boldsymbol{u} \rightarrow \boldsymbol{u}^*, \boldsymbol{u}^s \rightarrow \boldsymbol{u}^{k+1}，\phi \rightarrow p^{k+1} と置き換える．

ksttrksttr

圧力場の計算

圧力場 p^{k+1} の計算を行うために，上記のPoisson方程式を解くのだが，そのまま数値的にPoisson方程式を解くのではなく，もう少し計算をしておく．

\boldsymbol{u}^*は時間発展の中間で生成された仮想的な速度場であって，非圧縮性を保っていない．これによって，流体の質量密度 \rho の値が変化している．この時の質量密度を \rho^* としよう．質量密度は連続の式

\frac{D \rho}{Dt} + \rho\nabla\cdot\boldsymbol{u} = 0

を満たす．したがって，\rho^* は，

\frac{D \rho^*}{Dt} + \rho^*\nabla\cdot\boldsymbol{u^*} = 0

を満たす．これを同じように離散化すると，

\frac{\rho^* - \rho^0}{\Delta t} + \rho^* \nabla\cdot\boldsymbol{u}^* = 0

となって，結局，

\nabla \cdot \boldsymbol{u}^* = -\frac{1}{\Delta t} \frac{\rho^* - \rho^0}{\rho^*} \simeq -\frac{1}{\Delta t} \frac{\rho^* - \rho^0}{\rho^0}

となる．今，質量密度 \rho は粒子数密度 n に比例しているので，

\nabla \cdot \boldsymbol{u}^* = -\frac{1}{\Delta t} \frac{n^* - n^0}{n^0}

となる．したがって，圧力場 p^{k+1} の満たすPoisson方程式は，

\Delta p^{k+1} = -\frac{\rho^0}{\Delta t^2} \frac{n^* - n^0}{n^0}

となる．ここで，n^* は，

n^*_i = \sum_{j\not=i} w(|\boldsymbol{r}^*_j - \boldsymbol{r}^*_i|)

と計算する．

ksttrksttr

Poisson方程式の数値解法

以上で，MPS法の計算方法は全て説明したが，最後に，圧力場の計算で必要となるPoisson方程式の数値計算方法を書いておく．MPS法では，ラプラシアンを独自に離散化しているので，格子法によるものとは違うことに注意．

MPS法でのラプラシアンの離散化によって，圧力場のPoisson方程式は，

\frac{2d}{\lambda n_0} \sum_{j\not=i} (p^{k+1}_j - p^{k+1}_i) w(|\boldsymbol{r}^*_j - \boldsymbol{r}^*_i|) = -\frac{\rho^0}{\Delta t^2} \frac{n^*_i - n^0}{n^0}, \quad \mathrm{for} \quad i = 1,\dots, N

となる． \boldsymbol{p}^{k+1} \equiv \left( p^{k+1}_1, \dots , p^{k+1}_N \right)^\mathrm{T}，\boldsymbol{n}^* \equiv (n_1^*, \dots, n_N^*)^\mathrm{T}，(w_{ij})_{ij} = w(|\boldsymbol{r}^*_{j}-\boldsymbol{r}^*_i|)，

A \equiv \begin{pmatrix} -\frac{2dn_1^*}{n_0\lambda} & \frac{2d}{n_0\lambda} w_{12} & \cdots & \frac{2d}{n_0\lambda} w_{1 N-1} & \frac{2d}{n_0\lambda} w_{1N} \\ \frac{2d}{n_0\lambda} w_{21} & -\frac{2dn_2^*}{n_0\lambda} & \cdots & \frac{2d}{n_0\lambda} w_{2 N-1} & \frac{2d}{n_0\lambda} w_{2N} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \frac{2d}{n_0\lambda} w_{N-1 1} & \frac{2d}{n_0\lambda} w_{N-1 2} & \cdots & -\frac{2dn_{N-1}^*}{n_0\lambda} & \frac{2d}{n_0\lambda} w_{N-1 N} \\ \frac{2d}{n_0\lambda} w_{N1} & \frac{2d}{n_0\lambda} w_{N2} & \cdots & \frac{2d}{n_0\lambda} w_{N N-1} & -\frac{2dn_N^*}{n_0\lambda} \end{pmatrix}

とすると，A^{\mathrm{T}} = A であって，

A \boldsymbol{p}^{k+1} = -\frac{\rho^0}{\Delta t^2} \left(\frac{\boldsymbol{n}^*}{n_0} - \boldsymbol{1} _N \right) \equiv \boldsymbol{b}

と表すことができる．

したがって，k+1 ステップ目の圧力場を求めるには，上記の線形連立方程式を解けば良い．
線形連立方程式の数値解法には，不完全Cholesky分解前処理付き共役勾配法を用いるのが一般的のようである．

ksttrksttr

k+1 ステップ目の速度場

以上で圧力場 p^{k+1} が求められるので，これに粒子法における勾配の計算方法を用いて，\nabla p^{k+1} を求めれば，k+1 ステップ目の速度場 \boldsymbol{u}^{k+1} を求めることができる．すなわち，

\boldsymbol{u}^{k+1} = \boldsymbol{u}^* - \frac{1}{\rho} \nabla p^{k+1} \Delta t

これを用いて粒子位置を更新して，MPS法の1サイクルが完了する．

\boldsymbol{r}^{k+1} = \boldsymbol{r}^k + \boldsymbol{u}^{k+1} \Delta t
ksttrksttr

MPS法アルゴリズム

以上をまとめると，MPS法のアルゴリズムは以下のようになる

for k = 1 to T do

  1. \Delta \boldsymbol{u}^k_i = \frac{2d}{\lambda n_0} \sum_{j\not=i} (\boldsymbol{u}^{k+1}_j - \boldsymbol{u}^{k+1}_i) w(|\boldsymbol{r}^*_j - \boldsymbol{r}^*_i|) for i=1 to N
  2. \boldsymbol{u}^*_i = \boldsymbol{u}^k_i + (\nu\Delta \boldsymbol{u}^k_i )\Delta t + \boldsymbol{g} \Delta t, \quad \boldsymbol{r}^*_i = \boldsymbol{r}^k_i + \boldsymbol{u}^*_i \Delta t\quad for i = 1 to N
  3. n^* = \sum_{i=1}^N w(|\boldsymbol{r}^*_j - \boldsymbol{r}^*_i|)\quad for i = 1 to N
  4. Solve \frac{2d}{\lambda n_0} \sum_{j\not=i} (p^{k+1}_j - p^{k+1}_i) w(|\boldsymbol{r}^*_j - \boldsymbol{r}^*_i|) = -\frac{\rho^0}{\Delta t^2} \frac{n^*_i - n^0}{n^0} \quad in terms of p^{k+1}_i for i=1 to N
  5. \nabla p^{k+1}_i = \frac{d}{n_i} \sum_{j\not=i} \frac{(p^{k+1}_j - p^{k+1}_i) (\boldsymbol{r}_j - \boldsymbol{r}_i) }{|\boldsymbol{r}_j - \boldsymbol{r}_i|^2} w(|\boldsymbol{r}_j - \boldsymbol{r}_i|)
  6. \boldsymbol{u}^{k+1} = \boldsymbol{u}^* - \frac{1}{\rho} \nabla p^{k+1} \Delta t, \quad \boldsymbol{r}^{k+1} = \boldsymbol{r}^k + \boldsymbol{u}^{k+1} \Delta t\quad for i = 1 to N

end for

実際の計算には，これに加えて境界条件が必要である．

ksttrksttr

境界条件

以上に加えて，境界条件を粒子ごとに与えれば，MPS法の実装ができるようになる．
ここでは，基本的な境界条件である，固定壁境界条件・自由表面境界条件・流入・流出について書いておく．

境界条件の設定はなんだかよくわからないことが多い

固定壁

  • 粒子法では壁も粒子で構成する．
  • 流体粒子と接する面を構成する粒子を壁面表層粒子と呼び，
  • 壁面表層粒子の内側で，流体粒子とは接することのない粒子のことをダミー粒子と呼ぶ．
  • 2つを合わせて壁粒子と呼ぶ．

壁面表層粒子とダミー粒子の位置の更新は行われない．

ダミー粒子の必要性について

ダミー粒子が必要な理由は以下．

  • ダミー粒子がないと壁近傍の粒子の粒子数密度が小さくなってしまい，非圧縮性を保つことができない．
  • 自由表面との違いが判定できなくなる．

壁粒子の圧力

  • 壁表面粒子の圧力は，流体粒子と同様に更新する
  • ダミー粒子の圧力は更新しない．ダミー粒子の圧力は常に 0 としておく．
  • 圧力を決めるPoisson方程式はダミー粒子を含まないものとする．詳細は以下．

圧力を決める連立方程式の行列 A = (a_{ij}) の非対角成分のうち，ダミー粒子との関係を表す部分は全て0にする．

圧力を決める連立方程式の行列 A = (a_{ij}) の対角成分のうち，カーネルの影響領域内にダミー粒子を含むものに対して，

a_{ii} = \frac{2d}{\lambda n_0} \sum_{j\not=i, p_j \not= 0} w(|\boldsymbol{r}_j^* - \boldsymbol{r}_i^*|)

と修正する．これによって，非対角成分と体格成分の関係

\sum_{j\not=i} a_{ij} = -a_{ii}

を保つ．

以上は結局，ダミー粒子に関する量は圧力には関与しないことを意味する．Poisson方程式から完全にダミー粒子の影響は削除しておくのである．

これによって，壁面での圧力勾配は0になる（？？？）．したがって，固定壁において，圧力場にNeumann境界条件を課したことになる．

壁粒子の速度

  • no-slip条件を課すときは，壁面表層粒子の速度をすべて0にしておく．
  • free-slip条件を課すときは，粘性項の計算のときに，壁粒子と流体粒子の相互作用の部分を0に設定する．これは，流体粒子と壁粒子の速度が等しいとしていることと同値である．

自由表面

まだ

流入・流出

まだ

ksttrksttr

とりあえずMPS法を実装してみる

  • 言語はC++
  • dam break問題を解く

dam break問題

とは？

初期に設定する定数 (思いついた順)

物理名 記号 プログラム変数名
動粘性係数 \nu VISCOSITY
重力加速度 g GRAVITATIONAL_ACCELERATION
流体の質量密度 \rho FLUID_DENSITY
時間刻み \Delta t dt
粒子数 N N
計算するタイムステップ END_OF_TIME
初期の粒子間距離 l_0 l_0

このうち，粒子数 N と初期の粒子間距離 l_0 は，粒子の初期配置を決めることによって計算する．

また，粒子の影響半径 r_e は，l_0 の定数倍とし，演算子ごとにその値を変える．つまり，

r_e = \alpha l_0

としたとき， \alpha の値は，演算子ごとに変える．

ksttrksttr

無次元化?

まずは，方程式を無次元化するとともに，特徴的なスケールを決める．

計算に関与する次元を持つ量は，

  1. 粒子の初期速度 U
  2. 重力加速度の大きさg
  3. 粘性係数 \nu
  4. 質量密度\rho
  5. 境界の大きさ L
  6. 初期の粒子間距離 l_0
  7. 影響範囲 r_e
  8. 時間刻み幅 \Delta t

また，粒子数 N と粒子数密度 n もある．

ダム決壊問題に絞る

ダム決壊問題に絞って考えよう

この場合，U=0である．

この場合，重力加速度 g と粘性係数 \nu および，質量密度\rho で無次元化することになるだろうが，
これらは，ほとんど決まった値である．
つまり，実質，変化するパラメータは，境界の大きさ L くらいである．
言い換えると，無次元量

\left(\frac{\nu^2}{L^3g}\right)^{1/3}

を定義するのと，L をそのまま使うのは同じことであろう．

そこで今回は，水に対する値，

\begin{align*} \nu &= 0.000001 & [\mathrm{m^2/s}] \\ \rho &= 1 & [\mathrm{kg/m^3}] \end{align*}

を用いる．

また，g = 9.8 \, [\mathrm{m/s^2}] とする．

ksttrksttr

計算領域

計算領域を設定し，その中に（壁粒子・ダミー粒子もふくめて）粒子を置く．

計算領域は，2次元の場合は長方形，3次元の場合は直方体にする．

ダクトのような周囲が壁で囲まれている場合は，ダミー粒子を全て覆えるくらいの計算領域があればよい．

壁で覆われていない場合は，計算時間の間，粒子が領域からはみ出ない程度の計算領域が必要である．
ここで，初期条件や重力定数，粘性係数から見積もってみよう．

-> 見積もれる？？？

  • 流体が破砕する場合の水飛沫の飛距離はかなり自明でなさそう．
  • 次元解析で見積もる方法を探す

当面は，十分に大きな領域をとっておく．
例えば，

  • ダム崩壊問題だと，ダムの壁の２倍程度
  • 球体流体が落ちてくる問題だと，球体流体の初期位置の２倍程度の高さ，
  • あるいは，固定壁粒子（およびダミー粒子）で計算領域を囲む

など．

ksttrksttr

近傍粒子探索の方法

探索範囲を均等な格子に区切る方法

  1. 計算領域を均等な格子に分割し番号を振る．格子の大きさは粒子の影響範囲程度にする．
  2. ある格子の近傍格子の番号を計算する関数を用意する．
  3. 粒子が属する格子の番号を求め，記憶しておく．
  4. 粒子の属する格子の近傍格子を求める．
  5. 自分の格子及び近傍格子の中にある粒子についてだけ近傍粒子探索をする．
ksttrksttr

初期条件

計算領域は初期の状態を決めてから，決定するほうがよさそうなので，まずは初期条件を決める．

近傍粒子探索の方法

探索範囲を均等な格子に区切る方法

  1. 計算領域を均等な格子に分割し番号を振る．格子の大きさは粒子の影響範囲程度にする．
  2. ある格子の近傍格子の番号を計算する関数を用意する．
  3. 粒子が属する格子の番号を求め，記憶しておく．
  4. 粒子の属する格子の近傍格子を求める．
  5. 自分の格子及び近傍格子の中にある粒子についてだけ近傍粒子探索をする．
ksttrksttr

粒子間の距離が近すぎるときの対処方法

粒子間の距離が近くなりすぎると，カーネルの値が大きくなったり，圧力勾配の値が大きくなったりして，計算が不安定になりやすい．

そこで，粒子間の距離が近すぎるときは，人工的に何らかの処理をおこなって計算を安定化させるようである．

今回は，粒子同士がある一定の距離より近づいた時は，非弾性衝突したと考える．
衝突の反発係数を e としよう．

粒子 i に粒子 j が衝突した場合を考える．この時，粒子 i から見れば，粒子 j が相対速度速度
\boldsymbol{u}_j - \boldsymbol{u}_i でやってくるように見える．したがって，進行方向の速度

V \equiv -(\boldsymbol{u}_j - \boldsymbol{u}_i) \cdot \frac{\boldsymbol{r}_j - \boldsymbol{r}_i}{|\boldsymbol{r}_j - \boldsymbol{r}_i|}

で 静止した粒子 i に粒子 j がぶつかったと考えれば良い．

ksttrksttr

非弾性衝突時の速度変化

運動量変化と反発係数の定義から，粒子 i および j の衝突方向の相対速度を v_i, v_j とすると，

\begin{align*} m V &= - mv_i + mv_j \\ eV &= v_i + vj \end{align*}

であるから，衝突後の粒子 i の（衝突前に対する）相対速度は，

v_i = - \frac{1+e}{2} V

これは，衝突方向，(\boldsymbol{r_j} - \boldsymbol{r}_i)/|\boldsymbol{r_j} - \boldsymbol{r}_i| 方向の速度変化であるから，結局，衝突後の粒子 i の速度 \boldsymbol{u}_i^\prime は，

\boldsymbol{u}_i^\prime = \boldsymbol{u}_i - \frac{1+e}{2} V \frac{\boldsymbol{r}_j - \boldsymbol{r}_i}{|\boldsymbol{r}_j - \boldsymbol{r}_i|}

となる．

したがって，粒子 i と粒子 j がある距離よりも近づいた場合は，粒子 i について，速度を上記のように変更する．（粒子 j についても同様に行う．重複しないようにここでは，i だけ更新する．）

衝突した粒子全部による寄与は，

\boldsymbol{u}_i^\prime = \boldsymbol{u}_i - \sum_{j=1}^{N} \theta(R_c - |\boldsymbol{r}_j - \boldsymbol{r}_i|) \frac{1+e}{2} V \frac{\boldsymbol{r}_j - \boldsymbol{r}_i}{|\boldsymbol{r}_j - \boldsymbol{r}_i|}

ここで，R_c は衝突と判定する距離であり，\theta(x) はHeaviside関数である．

ksttrksttr

衝突判定距離

衝突と判定する距離はどのように選ぶか．

これはよくわからない．初期粒子間距離 l_0 の半分くらいが相場ようだ．

反発係数の値

これもよくわからないけど，あまり大きくないほうが良さそうである． ~ 0.2 ?

このスクラップは2022/09/03にクローズされました