広告
計算の流れ
半陰解法
計算は、タイムステップごとに半陰解法で計算されます。 現時刻のタイムステップをk,次の時刻のタイムステップをk+1とすると支配方程式の各項は、次式となります。
\[
\left[\frac{D\rho}{Dt}\right]^{k+1}=0
\]
\[
\frac{D\vec v}{Dt}
=
-\left[\frac{1}{\rho}\nabla P\right]^{k+1}
+\left[\upsilon\nabla^2\vec v\right]^k
+[g]^k
\]
上記の様に各項を陽的部分、陰的部分に分けて計算を行います。 計算の流れとしては、まず、陽的部分の計算を行い、その値を用いて陰的部分の計算を行います。 ここで、各粒子の速度、位置を次の様に定義します。
タイムステップkの速度・位置
\[
\vec v^k,\qquad \vec r^k
\]
陽的計算が終わった時点での速度・位置
\[
\vec v^*,\qquad \vec r^*
\]
タイムステップk+1の速度・位置
\[
\vec v^{k+1},\qquad \vec r^{k+1}
\]
陽的部分の計算
陽解法で解く項には、粘性項、重力項があります。
\[
\vec v^*
=
\vec v^k+\Delta t\,(\upsilon\nabla^2\vec v+g)^k
\]
\[
\vec r^*
=
\vec r^k+\Delta t\,\vec v^*
\]
\[
\left[\nabla^2\vec v\right]_i^k
=
\frac{2d}{\lambda n^0}
\sum_{j\ne i}
(\vec v_j^k-\vec v_i^k)
w(|\vec v_j^k-\vec v_i^k|)
\]
陰的部分の計算
次に陰解法で解く部分を説明します。
\[
\begin{aligned}
n^0&=n^{k+1}=n^*+n'\\[6pt]
\vec v^{k+1}&=\vec v^*+\vec v'\\[6pt]
\vec r^{k+1}&=\vec r^*+\vec r'
\end{aligned}
\]
陽解法で解いた項は、
\[
\vec v^*,\qquad \vec r^*
\]
陰解法で解く項は、
\[
\vec v',\qquad \vec r'
\]
になります。陰解法で解く項には、圧力項があります。圧縮性流れの質量収支式より
\[
\frac{D\rho}{Dt}+\rho\nabla\cdot\vec v=0
\]
左辺第2項の密度r[kg/m3]を一定値で近似すると
\[
\frac{D\rho}{Dt}+\rho_0\nabla\cdot\vec v=0
\]
流体の密度と粒子数密度は比例していることから
\[
\frac{1}{n^0}\frac{Dn}{Dt}+\nabla\cdot\vec v=0
\]
従って、修正量(')について解き、時間について離散化すると
\[
\frac{1}{n^0}\frac{n'}{\Delta t}+\nabla\cdot\vec v'=0
\]
となります。ここで、運動量収支式から速度の修正量は
\[
\vec v'
=
-\frac{\Delta t}{\rho_0}\nabla P^{k+1}
\]
となるので、勾配をとると
\[
(\nabla\cdot\vec v')
=
-\frac{\Delta t}{\rho_0}\nabla^2P^{k+1}
\]
となります。従って、
\[
\begin{aligned}
\frac{\Delta t}{\rho_0}\nabla^2P^{k+1}
&=
\frac{1}{n^0}\frac{n'}{\Delta t}\\[6pt]
&=
-\frac{1}{n^0}\frac{n^*-n^0}{\Delta t}
\end{aligned}
\]
となり
\[
\nabla^2P^{k+1}
=
-\frac{\rho_0}{\Delta t^2}
\frac{n^*-n^0}{n^0}
\]
が得られます。右辺は既知です。左辺は次式で表されます。
\[
\begin{aligned}
\nabla^2P_j^{k+1}
&=
\frac{2d}{\lambda n^0}
\sum_{j\ne i}(P_j^{k+1}-P_i^{k+1})
w(|\vec r_j-\vec r_i|)\\[6pt]
&=
\frac{2d}{\lambda n^0}
\sum_{j\ne i}w(|\vec r_j-\vec r_i|)P_j^{k+1}
-
\frac{2d}{\lambda n^0}
\sum_{j\ne i}w(|\vec r_j-\vec r_i|)P_i^{k+1}
\end{aligned}
\]
従って、連立方程式を解くことにより、k+1における圧力が計算されます。 計算された圧力から、次式より速度の修正量が計算されます。
\[
\vec v'
=
-\frac{\Delta t}{\rho_0}\nabla P^{k+1}
\]
ここで、圧力の勾配モデルは数値安定性のため、次式を使用します。
\[
\left[\nabla P_i^{k+1}\right]
=
\frac{d}{n^0}
\sum_{j\ne i}
\frac{P_j^{k+1}-\hat P_i^{k+1}}
{|\vec r_j-\vec r_i|^2}
(\vec r_j-\vec r_i)
w(|\vec r_j-\vec r_i|)
\]
これで、タイムステップk+1の速度が計算されました。
| prev | | | up | | | next |
