非ニュートン流体とMPSを用いた粘性応力項の離散化

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

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

https://qiita.com/fujitagodai4/items/65e186a845b9eeddc07b

　土木分野では、水と土砂の混合物を１つの流体として扱う場合があり、その場合は、構成則（抵抗則）が必要となってきます。

　今回の記事では、非ニュートン流体の構成則について説明します。

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

https://www.amazon.co.jp/%E7%B2%92%E5%AD%90%E6%B3%95%E5%85%A5%E9%96%80-%E6%B5%81%E4%BD%93%E3%82%B7%E3%83%9F%E3%83%A5%E3%83%AC%E3%83%BC%E3%82%B7%E3%83%A7%E3%83%B3%E3%81%AE%E5%9F%BA%E7%A4%8E%E3%81%8B%E3%82%89%E4%B8%A6%E5%88%97%E8%A8%88%E7%AE%97%E3%81%A8%E5%8F%AF%E8%A6%96%E5%8C%96%E3%81%BE%E3%81%A7-C-%E3%82%BD%E3%83%BC%E3%82%B9%E3%82%B3%E3%83%BC%E3%83%89%E4%BB%98/dp/4621088343

https://www.amazon.co.jp/%E7%B2%92%E5%AD%90%E6%B3%95-%E9%80%A3%E7%B6%9A%E4%BD%93%E3%83%BB%E6%B7%B7%E7%9B%B8%E6%B5%81%E3%83%BB%E7%B2%92%E7%8A%B6%E4%BD%93%E3%81%AE%E3%81%9F%E3%82%81%E3%81%AE%E8%A8%88%E7%AE%97%E7%A7%91%E5%AD%A6-%E5%BE%8C%E8%97%A4-%E4%BB%81%E5%BF%97/dp/4627922310


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

# 粘性応力項
　ナビエストークス方程式は、粘性応力項を $\Pi$ として

```math
\begin{align}
&\frac{Du}{Dt} = -\frac{1}{\rho}\nabla P +  \frac{1}{\rho}\nabla \cdot \Pi  + g \\
&\Pi = 2\mu D \ \ , \ \ D=\frac{1}{2}\left(\nabla u + (\nabla u )^T \right)
\end{align}
```
と書ける。$D$ は、変形速度テンソルである。粘性係数 $\mu$ が定数ならば、ニュートン流体と呼ばれ、動粘性係数 $\nu=\mu/\rho$ を使い、

```math
\begin{align}
\frac{Du}{Dt} = -\frac{1}{\rho}\nabla P + \nu \nabla^2 u + g \\
\end{align}
```
とかける。（非圧縮条件も使っている。）
　非ニュートン流体は、粘性係数 $\mu$ がせん断速度 $\dot \gamma$ の関数であり、

```math
\begin{align}
\mu = \mu(|\dot \gamma|)
\end{align}
```
書ける。以下、非ニュートン流体について説明する。

# 代表的な非ニュートン流体
　粘性係数は、Ostwarld-de Wales のべき乗則を使うと

```math
\begin{align}
\mu = \mu_0 |\dot \gamma|^{n-1} 
\end{align}
```
と書ける。せん断速度 $\dot \gamma$ は、変形速度テンソル $D$ を用いて

```math
\begin{align}
|\dot \gamma| \equiv |2D| = \sqrt{2D_{\alpha\beta}D_{\alpha\beta}}
\end{align}
```
と定義される。

　$n=1$ のとき定数になるのでニュートン流体、$n<1$ のとき擬塑性流体、$n>1$ のときダイラタント流体と呼ばれる。

　以下の図は、Ostwarld-de Wales のべき乗則をプロットしたものである。
![Viscosity_coefficient.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1521272/2b80beff-ca22-bec5-17d7-c28f67cbb825.png)

Ostwarld-de Wales のべき乗則に関して、以下の特徴が挙げられる。
1. 　ニュートン流体（$n=1$）の場合は、一定である。
1. 　擬塑性流体（$n<1$）の場合、せん断速度が大きくなるにつれて粘性が小さくなる。つまり、せん断速度が大きくなれば、粘性が弱くなり流動しやすくなる。例として、トマトケチャップやマヨネーズ、タルタルソースなどが挙げられる。

1. 　ダイラタント流体（$n>1$）の場合、せん断速度が大きくなるにつれて粘性が大きくなる。つまり、せん断速度が大きくなれば、粘性が強くなり流動しにくくなる。例として、非常に濃い水溶き片栗粉が挙げられる。


　また、単純せん断応力とせん断速度の関係は、
```math
\begin{align}
\tau = \mu(|\dot \gamma|) \dot \gamma
\end{align}
```
と書ける。
　以下の図は、単純せん断応力とせん断速度の関係をプロットしたものである。

![Flow_curve.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1521272/adf8a543-af5e-8684-d03f-13094a6d7f65.png)

　図に登場するBingham 流体の構成則は、
```math
\begin{align}
\mu = \mu_B + \frac{\tau_B}{|\dot \gamma|}
\end{align}
```
とかける。$\tau_B$ は降伏応力であり、降伏応力$\tau_B$ まで流動しないモデルである。 $\mu_B$ は降伏後の粘性係数である。つまり、Bingham 流体の粘性応力項は、

```math
\begin{align}
\Pi_{\alpha\beta} &= 
\begin{cases}
2\left(\mu_B + \frac{\tau_B}{|\dot \gamma|} \right) D_{\alpha\beta}     & (|\Pi|> \tau_B)  \\
0           & (|\Pi|\leq \tau_B)
\end{cases}
\\
\\

|\Pi| &= \sqrt{\frac{1}{2}\Pi_{\alpha\beta}\Pi_{\alpha\beta} }
\end{align}
```
と計算される。



# MPSを用いた粘性項の離散化
　最初に、変形速度テンソルの離散化を考える。

```math
\begin{align}
D_{\alpha\beta} = \frac{1}{2} \left(\frac{\partial u_{\alpha}}{\partial x_{\beta}} + \frac{\partial u_{\beta}}{\partial x_{\alpha}} \right)
\end{align}
```
せん断速度は、対称テンソルなので

```math
\begin{align}
|\dot \gamma|^2 &= 2D_{\alpha\beta}D_{\alpha\beta}\\
& = D_{\alpha\beta}\left(\frac{\partial u_{\alpha}}{\partial x_{\beta}} + \frac{\partial u_{\beta}}{\partial x_{\alpha}} \right) \\
& =2D_{\alpha\beta} \frac{\partial u_{\alpha}}{\partial x_{\beta}}\\
& = \frac{\partial u_{\alpha}}{\partial x_{\beta}}\frac{\partial u_{\alpha}}{\partial x_{\beta}} + \frac{\partial u_{\alpha}}{\partial x_{\beta}}\frac{\partial u_{\beta}}{\partial x_{\alpha}}
\end{align}
```
と書ける。

　物理量を $\phi$ とすると、$\phi$ の微分は

```math
\begin{align}
\frac{\partial \phi}{\partial x_{\alpha}} = \frac{\partial \phi}{\partial |r|}\frac{\partial |r|}{\partial x_{\alpha}} = \frac{\partial \phi}{\partial |r|}\frac{x_{\alpha}}{|r|}
\end{align}
```
なので、MPSでは、

```math
\begin{align}
\left\langle \frac{\partial \phi}{\partial x_{\alpha}} \right\rangle_{i} 
&= \frac{1}{n_0} \sum_{j \neq i} \omega(|r_{ij}|) \frac{\phi_{ij}}{|r_{ij}|^2} {x_{\alpha}}_{ij}
\end{align}
```
と離散化され、この式を使いせん断速度が求まる。（方向に依るので次元数はつけない。）

　次に、粘性項の離散化について考える。非圧縮条件より
```math
\begin{align}
2\frac{\partial }{\partial x_{\alpha}} D_{\alpha\beta} &= \frac{\partial }{\partial x_{\alpha}}\left(\frac{\partial u_{\alpha}}{\partial x_{\beta}} + \frac{\partial u_{\beta}}{\partial x_{\alpha}} \right) \\
&=\frac{\partial^2 u_{\beta}}{\partial x_{\alpha}\partial x_{\alpha}}
\end{align}
```
となるので、粘性応力項の部分は
```math
\begin{align}
\frac{1}{\rho} \nabla \cdot \Pi &= \frac{2}{\rho}(\nabla\mu)\cdot D  + \frac{\mu}{\rho} \nabla^2 u \\
& = \frac{1}{\rho}(\nabla\mu)\cdot \left\{\nabla u + (\nabla u )^T \right\}  + \frac{\mu}{\rho} \nabla^2 u
\end{align}
```
となる。それぞれ、MPSで離散化すると、

```math
\begin{align}
\left\langle \frac{1}{\rho} \nabla\mu \right\rangle_{i} 
&= \frac{D}{n_0} \sum_{j \neq i} \omega(|r_{ij}|)\frac{1}{\bar \rho_{ij}} \frac{\mu_{j} -\mu_j}{|r_{ij}|^2} r_{ij} \\
\left\langle  \nabla u \right\rangle_{i} 
&= \frac{D}{n_0} \sum_{j \neq i} \omega(|r_{ij}|) \frac{u_{ij} \otimes r_{ij}}{|r_{ij}|^2} \\
\left\langle  (\nabla u)^T \right\rangle_{i} 
&= \frac{D}{n_0} \sum_{j \neq i} \omega(|r_{ij}|) \frac{r_{ij} \otimes u_{ij}}{|r_{ij}|^2} \\
\left\langle \frac{\mu}{\rho} \nabla^2 u \right\rangle_i &= \frac{5-D}{n_0}\sum_{j \neq i}  \frac{\bar \mu_{ij}}{\bar \rho_{ij}} \frac{r_e}{|r_{ij}|^3} u_{ij} \\
\bar \rho_{ij} &= \frac{\rho_i+\rho_j}{2} \ \ , \ \ 
\bar \mu_{ij}= \frac{2\mu_i\mu_j}{\mu_i+\mu_j}
\end{align}
```
となる。

# 擬塑性流体とダイラタント流体の比較
　単純に擬塑性流体（$n<1$）とダイラタント流体（$n>1$ ）の比較を行う。

　以下は、擬塑性流体の damBreak の結果である。
![plastic.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1521272/a6466c65-72e7-0377-f7ff-9c3110632686.gif)
以下は、ダイラタント流体の damBreak の結果である。
![dilatant.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1521272/01856675-25f0-dcee-4f97-f954531ab1ea.gif)
　
　擬塑性流体のシミュレーションは、シミュレーションの後半にかけて、せん断速度が大きくなるため粘性が小さくなり、流動しやすい流体になっていることが分かる。

　ダイラタント流体のシミュレーションは、シミュレーションの後半にかけて、せん断速度が大きくなるため粘性が大きくなり、流動しにくい流体になっていることが分かる。

# まとめ
　非ニュートン流体における、離散化されたナビエストークス方程式は
```math
\begin{align}
　u_i^{k+1} &= u_i^{k}-\frac{1}{\rho} \langle \nabla P \rangle_i^{k+1} \Delta t + g \Delta t  \\
&+ \left\langle \frac{1}{\rho} \nabla\mu \right\rangle_{i}^k \cdot \left\{\left\langle  \nabla u \right\rangle_{i}^k  +\left\langle  (\nabla u)^T \right\rangle_{i} ^k \right\}\Delta t
  +\left\langle \frac{\mu}{\rho} \nabla^2 u \right\rangle_i ^k \Delta t \\
\end{align}
```
である。粘性応力項の部分は、それぞれ
```math
\begin{align}
\left\langle \frac{1}{\rho} \nabla\mu \right\rangle_{i} ^k 
&= \frac{D}{n_0} \sum_{j \neq i} \omega(|r_{ij}|)\frac{1}{\bar \rho_{ij}} \frac{\mu_{j}^k  -\mu_j^k }{|r_{ij}|^2} r_{ij} \\
\left\langle  \nabla u \right\rangle_{i} ^k 
&= \frac{D}{n_0} \sum_{j \neq i} \omega(|r_{ij}|) \frac{(u_{j}^k -u_{i}^k ) \otimes r_{ij}}{|r_{ij}|^2} \\
\left\langle  (\nabla u)^T \right\rangle_{i}^k  
&= \frac{D}{n_0} \sum_{j \neq i} \omega(|r_{ij}|) \frac{r_{ij} \otimes (u_{j}^k -u_{i}^k )}{|r_{ij}|^2} \\
\left\langle \frac{\mu}{\rho} \nabla^2 u \right\rangle_i^k  &= \frac{5-D}{n_0}\sum_{j \neq i}  \frac{\bar \mu_{ij}^k }{\bar \rho_{ij}} \frac{r_e}{|r_{ij}|^3} (u_{j}^k -u_{i}^k ) \\
\bar \rho_{ij} &= \frac{\rho_i+\rho_j}{2} \ \ , \ \ 
\bar \mu_{ij}^k = \frac{2\mu_i^k \mu_j^k }{\mu_i^k +\mu_j^k }
\end{align}
```
と計算される。

　粘性係数 $\mu$ は、Ostwarld-de Wales のべき乗則では

```math
\begin{align}
\mu_i^k = \mu_0 |\dot \gamma_i^k|^{n-1} 
\end{align}
```
Bingham 流体の構成則では
```math
\begin{align}
\mu_i^k = \mu_B + \frac{\tau_B}{|\dot \gamma|_i^k}
\end{align}
```
となる。

　せん断速度 $|\dot \gamma|$ は、粘性応力項を計算する前に計算する。せん断速度 $|\dot \gamma|$ は、

```math
\begin{align}
|\dot \gamma_i^k|
& =\sqrt{ \left\langle\frac{\partial u_{\alpha}}{\partial x_{\beta}} \right\rangle_{i} ^k
\left\langle\frac{\partial u_{\alpha}}{\partial x_{\beta}}  \right\rangle_{i} ^k
+\left\langle \frac{\partial u_{\alpha}}{\partial x_{\beta}} \right\rangle_{i} ^k
\left\langle\frac{\partial u_{\beta}}{\partial x_{\alpha}} \right\rangle_{i} ^k }
\end{align}
```
と書ける。平方根の中の微分は、

```math
\begin{align}
\left\langle \frac{\partial u_{\beta}}{\partial x_{\alpha}} \right\rangle_{i} ^k
&= \frac{1}{n_0} \sum_{j \neq i} \omega(|r_{ij}|) \frac{{u_{\beta}}_{j}^k -{u_{\beta}}_{i}^k }{|r_{ij}|^2} {x_{\alpha}}_{ij}
\end{align}
```
と計算される。