MPSを用いた固液２相流モデルの計算

# 概要

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

https://qiita.com/fujitagodai4/items/2d86a632e4180195e128

https://qiita.com/fujitagodai4/items/65e186a845b9eeddc07b#mps-%E3%81%A8-cmps-hs-hl-%E3%81%AE%E6%AF%94%E8%BC%83

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

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

https://www.amazon.co.jp/%E6%95%B0%E5%80%A4%E6%B5%81%E7%A0%82%E6%B0%B4%E7%90%86%E5%AD%A6-POD%E7%89%88%E2%80%95%E7%B2%92%E5%AD%90%E6%B3%95%E3%81%AB%E3%82%88%E3%82%8B%E6%B7%B7%E7%9B%B8%E6%B5%81%E3%81%A8%E7%B2%92%E7%8A%B6%E4%BD%93%E3%81%AE%E8%A8%88%E7%AE%97%E5%8A%9B%E5%AD%A6%E2%80%95-%E5%BE%8C%E8%97%A4-%E4%BB%81%E5%BF%97/dp/4627495692/ref=pd_vtp_7/357-8805296-1691253?pd_rd_w=bQrZD&pf_rd_p=ae64b7f5-458b-4e9a-9b07-1feecb909091&pf_rd_r=TABR4Q08XTRS89HMEX9X&pd_rd_r=c8473a72-478a-4cbf-9aa7-e5bebf3f641c&pd_rd_wg=qe5Yt&pd_rd_i=4627495692&psc=1 


https://www.pwri.go.jp/jpn/results/report/report-program/2016/pdf/pro-3-6.pdf


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

```math
\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}$ は、流体-固体の相互作用項であり

```math
\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 のとき、方程式がナビエストークス方程式になると考えてプラスの符号にした。）

　質量が保存されるので、

```math
\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$ ステップの流速と圧力を使い、
```math
\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$ ステップの流速の更新を

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

```math
\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$ が液体の場合
```math
\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$ が固体の場合も同様に
```math
\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}$ は、

 ```math
\begin{align}
\nu_{mix} = \frac{\nu_l}{\sqrt{1+\frac{\rho_s}{\rho_l} c_i}}
\end{align}
```
と書ける。$c_i$ は濃度であり、

 ```math
\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}
```

```math
\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$ とする。

よって
```math
\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}
```
を
```math
\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}
```
と書き換える。

# 圧力項の計算

　仮速度・仮位置を
```math
\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$ が液体の場合

```math
\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$ が固体の場合も同様に

```math
\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$ の形にすると

```math
\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}
```
```math
\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$ を使い　圧力勾配を
 
```math
\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$ は、

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

# まとめ

 以下は、砂礫の沈降・堆積仮定のシミュレーションである。
![anime.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1521272/d897d13c-064e-087c-9316-e526a7d7ad86.gif)

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

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

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