概要
土木の分野において、粒子法の研究が盛んにおこなわれております。前回の記事で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}
である。
まとめ
以下は、砂礫の沈降・堆積仮定のシミュレーションである。
シミュレーションの左側の粒子を砂と仮定し、シミュレーションの右側の粒子を水と仮定して計算を行った。砂の方が重いので、砂が沈降していることが分かる。しかし、恐らく、実際の現象で砂の流動は、右側の壁に当たって上昇することはないと考えられる。固相の振る舞いも、液相に近い振る舞いになっている。
この記事で引用した論文では、抵抗則を支配方程式に入れることで、実際の現象に近い結果が得られている。その抵抗則は、濃度に依存する。
今後は、抵抗則を勉強したい。