流体計算における支配方程式
連続の式
$$
\nabla\cdot \mathbf{u}=0\tag{1}
$$
湧き出し量（単位体積当たりの流れの増減）が$\mathrm{div}(\mathbf{u})$であるため、非圧縮性流体であることを指す。
空気等は圧縮性流体だが、音速を超えて移動しない限り大きな密度変化は起こらないので、湧き出し量はゼロとして考える。
参考
運動量保存則
Navier-Stokes方程式
$$
\begin{aligned}\rho\frac{d\mathbf{u}}{d t}&=\mathbf{F}_p+\mathbf{F}_{\mu}+\mathbf{F}_g+\mathbf{F}_s\\\frac{\partial\mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}&=-\frac{\nabla p}{\rho}+\nu\cdot\nabla^2\mathbf{u}+\mathbf{g}+\frac{\sigma}{\rho}(\kappa_1+\kappa_2)\mathbf{n}_s\end{aligned}\tag{2}$$
$F_p$：圧力項
$F_{\mu}$：粘性項
$F_g$：重力項
$F_s$：表面張力項
左辺
$\partial\mathbf{u}/\partial t$は時間項、$(\mathbf{u}\cdot\nabla)\mathbf{u}$は移流項と呼ばれる
2つ合わせて加速度を表すもので、速度$\mathbf{u}$を時間微分したら出てくる式
以下導出
$\mathbf{u}$は$t,x(t),y(t),z(t)$の4変数を持ち、$x,y,z$が$t$従属変数なので連鎖律によって次式のように変化量を表せる
$$
\begin{aligned}\frac{d\mathbf{u}}{dt}&=\frac{\partial\mathbf{u}}{\partial t}+\frac{\partial\mathbf{u}}{\partial x}\frac{dx}{dt}+\frac{\partial\mathbf{u}}{\partial y}\frac{dy}{dt}+\frac{\partial\mathbf{u}}{\partial z}\frac{dz}{dt}\\\frac{d\mathbf{u}}{dt}&=\frac{\partial\mathbf{u}}{\partial t}+\left(u_x\frac{\partial}{\partial x}+u_y\frac{\partial}{\partial y}+u_z\frac{\partial}{\partial z}\right)\mathbf{u}\\\frac{d\mathbf{u}}{dt}&=\frac{\partial\mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}\end{aligned}\tag{3}
$$
$\nabla\cdot\mathbf{u}$が各要素の偏微分であることに対し、$\mathbf{u}\cdot\nabla$はナブラ演算子の各要素に係数を追加する操作である
非圧縮性流体における連続の式$\nabla\cdot \mathbf{u}=0$から、
$(\mathbf{u}\cdot\nabla)\mathbf{u}=(\nabla\cdot\mathbf{u})\mathbf{u}=0$
と計算できるように見えるが、順番によって操作が異なるため成立しない
また非圧縮性流体にのみ、以下の変形が成り立つ
$$
\begin{aligned}\nabla\cdot(\mathbf{u}\mathbf{u})&=(\mathbf{u}\cdot\nabla)\mathbf{u}+\mathbf{u}(\nabla\cdot \mathbf{u})\\\nabla\cdot(\mathbf{u}\mathbf{u})&=(\mathbf{u}\cdot\nabla)\mathbf{u}\end{aligned}
$$全微分を使った導出方法
ナビエ-ストークス方程式の導出
右辺
流体にかかる力の和
運動方程式立てるときと同じ
圧力項
圧力の向きは面の外側が正で、単位体積にかかる圧力を求めたいのでマイナスになる
単位体積の1面積は単位面積なので、$-\nabla p\cdot S=-\nabla p\cdot1=-\nabla p$
左辺が加速度のため、密度で割ることで加速度としている
$$\begin{aligned}ma=\rho Va=\rho a&=F\\a&=\frac{F}{\rho}\end{aligned}\tag{4}$$
$(単位体積V=1)$
粘性項
圧力項が正面からの力だとすると、粘性項は横からの力
ニュートンの粘性法則より
$$\tau=\mu(\nabla\mathbf{u}+\nabla\mathbf{u}^{\mathrm{T}})\tag{5}$$
分解すると
$$\tau=\mu\begin{bmatrix}2\partial_xu&\partial_yu+\partial_xv&\partial_zu+\partial_xw\\\partial_yu+\partial_xv&2\partial_yv&\partial_zv+\partial_yw\\\partial_zu+\partial_xw&\partial_zv+\partial_yw&2\partial_zw\\\end{bmatrix}\tag{6}$$
$(\partial_n=\frac{\partial}{\partial n})$
これらを簡易化して
$$\tau_{ij}=\mu(\partial_{j} u_i+\partial_{i} u_j)\tag{7}$$
（$i,j$にはそれぞれ$x,y,z$が入る）
単位位体積あたりにかかる粘性力は
$$
\mathbf{F}_\mu=\begin{bmatrix}\partial_x\tau_{xx}+\partial_y\tau_{yx}+\partial_z\tau_{zx}\\\partial_x\tau_{xy}+\partial_y\tau_{yy}+\partial_z\tau_{zy}\\\partial_x\tau_{xz}+\partial_y\tau_{yz}+\partial_z\tau_{zz}\end{bmatrix}^{\mathrm{T}}\tag{8}
$$
簡易化すると
$$
(F_\mu)_i=\partial_x\tau_{xi}+\partial_y\tau_{yi}+\partial_z\tau_{zi}\tag{9}
$$式(8)に式(6)を代入すると、
$$
\begin{aligned}(F_\mu)_i&=\mu(\partial_x(\partial_{i} u_x+\partial_{x} u_i)+\partial_y(\partial_{i} u_y+\partial_{y} u_i)+\partial_z(\partial_{i} u_z+\partial_{z} u_i))\\(F_\mu)_i&=\mu(\partial_x(\partial_{i} u_x)+\partial_x(\partial_{x} u_i)+\partial_y(\partial_{i} u_y)+\partial_y(\partial_{y} u_i)+\partial_z(\partial_{i} u_z)+\partial_z(\partial_{z} u_i))\\(F_\mu)_i&=\mu(\partial_x(\partial_{i} u_x)+\partial_{x}^2 u_i+\partial_y(\partial_{i} u_y)+\partial_{y}^2 u_i+\partial_z(\partial_{i} u_z)+\partial_{z}^2 u_i)\\(F_\mu)_i&=\mu(\partial_{i}(\partial_x u_x+\partial_y u_y+\partial_z u_z)+(\partial_{x}^2+\partial_{y}^2+\partial_{z}^2 )u_i)\\(F_\mu)_i&=\mu(\partial_{i}(\nabla\cdot\mathbf{u})+(\partial_{x}+\partial_{y}+\partial_{z} )^2u_i)\\(F_\mu)_i&=\mu\partial_{i}(\nabla\cdot\mathbf{u})+\mu\nabla^2 u_i\\
\end{aligned}\tag{10}
$$
非圧縮性流体と仮定すると、式(1)より、
$$
(F_\mu)_i=\mu\nabla^2 u_i\tag{11}
$$
即ち、
$$
\mathbf{F}_\mu=\mu\nabla^2 \mathbf{u}\tag{12}
$$
圧力項と同じように密度で割り、加速度としてやると、
$$
\begin{aligned}\frac{\mathbf{F}_\mu}{\rho}&=\frac{\mu\nabla^2 \mathbf{u}}{\rho}\\\frac{\mathbf{F}_\mu}{\rho}&=\frac{\mu}{\rho}\nabla^2 \mathbf{u}\\\frac{\mathbf{F}_\mu}{\rho}&=\nu \nabla^2 \mathbf{u}\end{aligned}\tag{13}
$$重力項
$$
\begin{aligned}\mathbf{F}_g=\rho\mathbf{g}\\\frac{\mathbf{F}_g}{\rho}=\mathbf{g}\end{aligned}
$$
特に言うことなし
表面張力項
ヤングラプラスの式
$$
\mathbf{F}_s=\Delta p_s=\sigma\left(\frac{1}{R_1}+\frac{1}{R_2}\right)
$$
$\sigma$：表面張力係数
$R$：曲率半径
即ち、
$$
\frac{\mathbf{F}_s}{\rho}=\frac{\sigma}{\rho}(\kappa_1+\kappa_2)\mathbf{n}_s
$$
$1/R=\kappa$：曲率
$\mathbf{n}_s$：単位法線ベクトル
VOF法における体積分率方程式
ナビエストークス方程式の左辺と同じように解く
$$
\frac{d\alpha}{dt}=0
$$
$\alpha$が$t,x(t),y(t),z(t)$の関数だとして、連鎖律から、
$$
\begin{aligned}\frac{\partial\alpha}{\partial t}\frac{dt}{dt}+\frac{\partial\alpha}{\partial x}\frac{dx}{dt}+\frac{\partial\alpha}{\partial y}\frac{dy}{dt}+\frac{\partial\alpha}{\partial z}\frac{dz}{dt}&=0\\\frac{\partial\alpha}{\partial t}+(\mathbf{u}\cdot\nabla)\alpha&=0\end{aligned}
$$
変形して、
$$
\frac{\partial \alpha}{\partial t}+\nabla\cdot(\alpha\mathbf{u})=0
$$https://arxiv.org/pdf/1703.00937
https://www.sciencedirect.com/science/article/pii/S0021999112001830?via%3Dihub