概要

土木の分野において、粒子法の研究が盛んにおこなわれております。前回の記事でMPSについて説明しました。

　今回の記事では、水と土砂を扱うような、2相流モデルについて説明したいと思います。このモデルは、固相の支配方程式を流体のように扱うことで、2相流の計算を行います。

以下の書籍・論文の内容を参考にしています。

固液2相流の支配方程式

　固液2相流の支配方程式は、抵抗則などがない場合

\begin{align}
\frac{D\rho_l}{Dt}& + \rho_l\nabla \cdot u_l = 0 \\
\frac{D\rho_s}{Dt}& + \rho_s\nabla \cdot u_s = 0 \\

\frac{Du_l}{Dt} &= \left(-\frac{1}{\rho_l}\nabla P_l + \nu_l \nabla^2 u_l \right)_l + F_{ls-l} + g  \\
\frac{Du_s}{Dt} &= \left(-\frac{1}{\rho_s}\nabla P_s + \nu_s \nabla^2 u_s \right)_s + F_{ls-s} + g  \\
\end{align}

と書ける。$u$ は流速、$P$ は圧力、$ρ$ は密度、$g$ は重力加速度である。下付き添え字 $l,s$ は、液相と固相を意味する。括弧 $()$ の下付き添え字は、添え字の相の粒子のみで計算することを意味する。
　$F_{ls-l},F_{ls-s}$ は、流体-固体の相互作用項であり

\begin{align}
F_{ls-l} &= \left(-\frac{1}{\rho_l}\nabla P_l + \nu_l \nabla^2 u_l \right)_s  \\
F_{ls-s} &= \left(-\frac{1}{\rho_s}\nabla P_s + \nu_s \nabla^2 u_s \right)_l  \\
\end{align}

と書ける。（論文によっては、マイナス符号が付いている。l=s のとき、方程式がナビエストークス方程式になると考えてプラスの符号にした。）

　質量が保存されるので、

\begin{align}
\frac{D\rho_l}{Dt}  = 0 & \Rightarrow  \ \ \nabla \cdot u_l = 0 \\
\frac{D\rho_s}{Dt}  = 0 & \Rightarrow  \ \ \nabla \cdot u_s = 0 \\
\end{align}

となり非圧縮条件が得られる。

MPSの離散化

　固液2相流の支配方程式は、$k+1$ ステップの流速は、$k$ ステップの流速と圧力を使い、

\begin{align}
u_l^{k+1} &= u_l^{k}-\frac{1}{\rho_l} \left\{ \langle \nabla P_l \rangle_l + \langle \nabla P_l \rangle_s \right\} ^{k+1}\Delta t + \left\{ \langle \nu_l\nabla^2 u_l  \rangle_l + \langle \nu_l\nabla^2 u_l \rangle_s \right\}^k \Delta t + g \Delta t \\
u_s^{k+1} &= u_s^{k}-\frac{1}{\rho_s} \left\{ \langle \nabla P_s \rangle_l + \langle \nabla P_s \rangle_s \right\} ^{k+1}\Delta t + \left\{ \langle \nu_s\nabla^2 u_s  \rangle_l + \langle \nu_s\nabla^2 u_s \rangle_s \right\}^k \Delta t + g \Delta t \\
\end{align}

と離散化される。$\langle \rangle$ の下付き添え字は、添え字の相の粒子のみで計算することを意味する。 $k+1$ ステップの流速の更新を

\begin{align}
u_m^{k+1} & = u_m^{k} + \delta u_m^p + \delta u_m^c \ \ , \ \ m=l,s
\end{align}

\begin{align}
\delta u_m^p & = \left\{ \langle \nu_m\nabla^2 u_m  \rangle_l + \langle \nu_m\nabla^2 u_m \rangle_s \right\}^k \Delta t + g \Delta t  \\
\delta u_m^c & = -\frac{1}{\rho_m} \left\{ \langle \nabla P_m \rangle_l + \langle \nabla P_m \rangle_s \right\} ^{k+1}\Delta t
\end{align}

と書くことにする。

粘性項の計算

　ラプラシアンの計算は、粒子 $i$ が液体の場合

\begin{align}
{\langle \nabla^2 u_l \rangle_m^k}_{i}  & = \frac{5-D}{n_0}\sum_{j \neq i}  \frac{r_e}{|r_{ij}|^3}({u_m}_j^k  - {u_l}_i^k ) \\
r_{ij} &= {r_m}_j - {r_l}_i 
\end{align}

と書ける。$m$ に関しては、$m=l$ ならば粒子 $j$ は液体、 $m=s$ ならば粒子 $j$ は固体である。粒子 $i$ が固体の場合も同様に

\begin{align}
{\langle \nabla^2 u_s \rangle_m ^k}_{i} & = \frac{5-D}{n_0}\sum_{j \neq i}  \frac{r_e}{|r_{ij}|^3}({u_m}_j^k  - {u_s}_i^k ) \\
r_{ij} &= {r_m}_j - {r_s}_i 
\end{align}

と書ける。

　動粘性係数$\nu$ は、液相・固相のどちらの値もする必要があるが、混合体の粘性 $\nu_{mix}$ を使うことで液相の動粘性係数だけで計算を行う。
　混合体の粘性 $\nu_{mix}$ は、

\begin{align}
\nu_{mix} = \frac{\nu_l}{\sqrt{1+\frac{\rho_s}{\rho_l} c_i}}
\end{align}

と書ける。$c_i$ は濃度であり、

\begin{align}
c_i = \frac{\sum_{i\neq j} \delta_{sj} \omega_{c}(|r_{ij}|)}{\sum_{i\neq j} \omega_{c}(|r_{ij}|)}
\end{align}

\begin{align}
\omega_c(|r_{ij}|) = 
\begin{cases}
1 & (0 \leq |r_{ij}| \leq r_e)  \\
0             & (r_e \leq |r_{ij}|)
\end{cases}
\end{align}

と書ける。$\delta_{sj}$ は、$j$ の粒子が固体なら $1$ 、$j$ の粒子が液体なら $0$ とする。

よって

\begin{align}
\delta u_m^p & = \left\{ \langle \nu_m\nabla^2 u_m  \rangle_l + \langle \nu_m\nabla^2 u_m \rangle_s \right\}^k \Delta t + g \Delta t  \\
\end{align}

を

\begin{align}
\delta u_m^p & = \nu_{mix}\left\{ \langle \nabla^2 u_m  \rangle_l + \langle \nabla^2 u_m \rangle_s \right\}^k \Delta t + g \Delta t  \\
\end{align}

と書き換える。

圧力項の計算

　仮速度・仮位置を

\begin{align}
u_m^* &= u_m^k + \delta u_m^p \\
r_m^* &= r_m^k + u_m^* \Delta t 
\ \ , \  m=l,s
\end{align}

とする。

ポアソン方程式は、粒子 $i$ が液体の場合

\begin{align}
\frac{1}{n_0}\sum_{j \neq i}  \left(\frac{6-D}{\rho_l} -\frac{1}{\rho_m} \right) \frac{r_e}{|r_{ij}^*|^3}({P_m}_j^{k+1} - {P_l}_i^{k+1}) &= \frac{1}{n_0 \Delta t} \sum_{j \neq i} \frac{r_e}{|r_{ij}^*|^3} r_{ij}^*\cdot u_{ij}^* \\
r_{ij}^* = {r_m}_j^* - {r_l}_i^* \\
u_{ij}^* = {u_m}_j^* - {u_l}_i^* \\
\end{align}

と書ける。$m$ は、粒子 $j$ が液体ならば $m=l$ であり、粒子 $j$ が固体ならば $m=s$ である。粒子 $i$ が固体の場合も同様に

\begin{align}
\frac{1}{n_0}\sum_{j \neq i}  \left(\frac{6-D}{\rho_s} -\frac{1}{\rho_m} \right) \frac{r_e}{|r_{ij}^*|^3}({P_m}_j^{k+1} - {P_s}_i^{k+1}) &= \frac{1}{n_0 \Delta t} \sum_{j \neq i} \frac{r_e}{|r_{ij}^*|^3} r_{ij}^*\cdot u_{ij}^* \\
r_{ij}^* = {r_m}_j^* - {r_s}_i^* \\
u_{ij}^* = {u_m}_j^* - {u_s}_i^* \\
\end{align}

と書ける。連立方程式$Ax=b$ の形にすると

\begin{align}
&Ax=b　\\
&x_i = P_i^{k+1} , \ \  b_i = \frac{1}{n_0 \Delta t} \sum_{j \neq i} \frac{r_e}{|r_{ij}^*|^3} r_{ij}^*\cdot u_{ij}^*
\end{align}

\begin{align}
A_{ij} = 
\begin{cases}
 -\frac{1}{n_0}\sum_{j \neq i}  \left(\frac{6-D}{\rho_i} -\frac{1}{\rho_j} \right) \frac{r_e}{|r_{ij}^*|^3}& (i=j)  \\
 \frac{1}{n_0}  \left(\frac{6-D}{\rho_i} -\frac{1}{\rho_j} \right) \frac{r_e}{|r_{ij}^*|^3} & (i\neq j)
\end{cases}
\end{align}

の形になる。$\rho_i$ または $\rho_j$ は、粒子 $i$ または $j$ が液体のとき $\rho_l$、粒子 $i$ または $j$ が固体のとき $\rho_s$ とする。上記の連立方程式の係数行列 $A$ は、非対称行列なのでBi-CGSTAB 法などの解法で解く必要がある。

　求まった圧力 $P$ を使い　圧力勾配を

\begin{align}
\langle \nabla P \rangle_{i}
&= \frac{D}{n_0} \sum_{j \neq i} \omega(|r_{ij}|) \frac{\left(P_i+P_j \right)-\left(\hat{P}_i +\hat{P}_j\right) }{|r_{ij}|}\frac{r_{ij}}{|r_{ij}|} \\
\end{align}

を用いて計算する。$\hat P$ は、

\begin{align}
\hat{P}_i &= \min_{j\in \Omega} (P_i,P_j) \\
\Omega &= \{j:\omega(|r_{ij}|) \neq 0\}
\end{align}

である。

まとめ

以下は、砂礫の沈降・堆積仮定のシミュレーションである。

　シミュレーションの左側の粒子を砂と仮定し、シミュレーションの右側の粒子を水と仮定して計算を行った。砂の方が重いので、砂が沈降していることが分かる。しかし、恐らく、実際の現象で砂の流動は、右側の壁に当たって上昇することはないと考えられる。固相の振る舞いも、液相に近い振る舞いになっている。

　この記事で引用した論文では、抵抗則を支配方程式に入れることで、実際の現象に近い結果が得られている。その抵抗則は、濃度に依存する。

　今後は、抵抗則を勉強したい。