概要

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

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

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

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

粘性応力項

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

\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$ を使い、

\begin{align}
\frac{Du}{Dt} = -\frac{1}{\rho}\nabla P + \nu \nabla^2 u + g \\
\end{align}

とかける。（非圧縮条件も使っている。）
　非ニュートン流体は、粘性係数 $\mu$ がせん断速度 $\dot \gamma$ の関数であり、

\begin{align}
\mu = \mu(|\dot \gamma|)
\end{align}

書ける。以下、非ニュートン流体について説明する。

代表的な非ニュートン流体

　粘性係数は、Ostwarld-de Wales のべき乗則を使うと

\begin{align}
\mu = \mu_0 |\dot \gamma|^{n-1} 
\end{align}

と書ける。せん断速度 $\dot \gamma$ は、変形速度テンソル $D$ を用いて

\begin{align}
|\dot \gamma| \equiv |2D| = \sqrt{2D_{\alpha\beta}D_{\alpha\beta}}
\end{align}

と定義される。

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

　以下の図は、Ostwarld-de Wales のべき乗則をプロットしたものである。

Ostwarld-de Wales のべき乗則に関して、以下の特徴が挙げられる。

　ニュートン流体（$n=1$）の場合は、一定である。
　擬塑性流体（$n<1$）の場合、せん断速度が大きくなるにつれて粘性が小さくなる。つまり、せん断速度が大きくなれば、粘性が弱くなり流動しやすくなる。例として、トマトケチャップやマヨネーズ、タルタルソースなどが挙げられる。
　ダイラタント流体（$n>1$）の場合、せん断速度が大きくなるにつれて粘性が大きくなる。つまり、せん断速度が大きくなれば、粘性が強くなり流動しにくくなる。例として、非常に濃い水溶き片栗粉が挙げられる。

　また、単純せん断応力とせん断速度の関係は、

\begin{align}
\tau = \mu(|\dot \gamma|) \dot \gamma
\end{align}

と書ける。
　以下の図は、単純せん断応力とせん断速度の関係をプロットしたものである。

　図に登場するBingham 流体の構成則は、

\begin{align}
\mu = \mu_B + \frac{\tau_B}{|\dot \gamma|}
\end{align}

とかける。$\tau_B$ は降伏応力であり、降伏応力$\tau_B$ まで流動しないモデルである。 $\mu_B$ は降伏後の粘性係数である。つまり、Bingham 流体の粘性応力項は、

\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を用いた粘性項の離散化

　最初に、変形速度テンソルの離散化を考える。

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

せん断速度は、対称テンソルなので

\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$ の微分は

\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では、

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

と離散化され、この式を使いせん断速度が求まる。（方向に依るので次元数はつけない。）

　次に、粘性項の離散化について考える。非圧縮条件より

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

となるので、粘性応力項の部分は

\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で離散化すると、

\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 の結果である。

以下は、ダイラタント流体の damBreak の結果である。

　
　擬塑性流体のシミュレーションは、シミュレーションの後半にかけて、せん断速度が大きくなるため粘性が小さくなり、流動しやすい流体になっていることが分かる。

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

まとめ

　非ニュートン流体における、離散化されたナビエストークス方程式は

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

である。粘性応力項の部分は、それぞれ

\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 のべき乗則では

\begin{align}
\mu_i^k = \mu_0 |\dot \gamma_i^k|^{n-1} 
\end{align}

Bingham 流体の構成則では

\begin{align}
\mu_i^k = \mu_B + \frac{\tau_B}{|\dot \gamma|_i^k}
\end{align}

となる。

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

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

と書ける。平方根の中の微分は、

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

と計算される。