概要
土木の分野において、粒子法の研究が盛んにおこなわれています。前回の記事で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}
と計算される。