「Note」Basic Fluid Simulation
整理文件的时候翻出了之前学习流体模拟的一些笔记，先放一些在这里，后续的等有空的时候再整理一次(
微积分
梯度 (Gradient)
$$
\nabla f(x, y, z) = \left(\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}, \frac{\partial f}{\partial z}\right)
$$
其表示了函数 $f$ 上升最快的方向
多个函数的梯度则构成了一个矩阵：
$$ \nabla \boldsymbol F = \nabla(f, g, h) = \begin{pmatrix} \nabla f \\ \nabla g \\ \nabla h\end{pmatrix} $$散度 (Divergence)
对于向量 $\boldsymbol u = (u, v, w)$
$$
\mathrm {div} \boldsymbol u = \nabla \cdot \boldsymbol u = \frac {\partial u}{\partial x} + \frac {\partial v}{\partial y} + \frac {\partial w}{\partial z}
$$
其表示了某一点的发散情况。
旋度 (Curl)
$$
\nabla \times \boldsymbol u = \nabla \times (u, v, w) = \left(\frac {\partial w}{\partial y} - \frac {\partial v}{\partial z}, \frac {\partial u}{\partial z} - \frac{\partial w}{\partial x}, \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}\right)
$$
其表示了矢量场中某一点所包含的微元在场中的旋转程度。
拉普拉斯算子 (Laplacian)
$$
\nabla ^2 f = \nabla \cdot \nabla f = \frac {\partial ^ 2 f}{\partial x ^ 2} + \frac {\partial ^2 f}{\partial y ^ 2} + \frac {\partial ^ 2 f}{\partial z ^ 2}
$$
有向面积元
$$
\iint_S \boldsymbol v \cdot \boldsymbol n \mathrm d S
$$
为向量场 $\boldsymbol v$ 在定向曲面 $S$ 上的第二型曲面积分，令 $\mathrm d \boldsymbol S$ 表示有向量面积微元（这里方向均为外法向），即 $\mathrm d \boldsymbol S = \boldsymbol n \mathrm d S$，积分简化表示为
$$
\iint_S \boldsymbol v \cdot \mathrm d \boldsymbol S
$$
高斯公式
这里只列出了向量的形式，令 $V$ 表示一简单闭曲面 $S$ 为边界的体积，$\boldsymbol v$ 是定义在 $V$ 中和 $S$ 上连续可微的向量场，有
$$
\require{mediawiki-texvc}\oiint_S \boldsymbol v \cdot \mathrm d \boldsymbol S = \iiint_V \nabla \cdot \boldsymbol v \mathrm d V
$$
全微分
可微函数 $z = f(x, y)$，则有
$$
\mathrm d z = \frac {\partial f}{\partial x} \mathrm d x + \frac {\partial f}{\partial y}\mathrm d y
$$
拉格朗日描述和欧拉描述
研究流体通常有两种方法来描述：拉格朗日描述（Lagrangian description）和欧拉描述（Eulerian description）
拉格朗日描述是观察流体粒子的性质，而欧拉描述观察空间中一个固定点的性质。
物质导数 (Material Derivative)
许多定理是对于流体粒子成立的，而我们经常需要其在欧拉描述里的结果，所以我们需要将两种描述建立起联系，而这个联系就是物质导数。
每个粒子都有自己的位置 $\boldsymbol x$ 和速度 $\boldsymbol u$，令 $q$ 表示一个通用的物理量（如密度、速度等），每个粒子都有对应的 $q$ 值，$q(t, \boldsymbol x)$ 描述了在时间 $t$ 位置 $\boldsymbol x$ 的粒子对应的物理量的值 $q$。该物理量对于时间的变化率为（注意这里是全导数）
$$ \begin{aligned} \frac {\mathrm d }{\mathrm d t}q(t, \boldsymbol x) &= \frac {\partial q}{\partial t}\frac {\mathrm d t}{\mathrm d t} + \frac {\partial q}{\partial \boldsymbol x}\frac{\mathrm{d}\boldsymbol x} {\mathrm d t}\\ &= \frac {\partial q}{\partial t}+ \nabla q \cdot \boldsymbol u \end{aligned} $$于是就可以定义物质导数 $\frac{\mathrm D}{\mathrm D t}$
$$
\frac {\mathrm D}{\mathrm D t} \triangleq \frac {\mathrm d }{\mathrm d t} = \frac {\partial }{\partial t} + (\boldsymbol u \cdot \nabla)
$$
雷诺传输定理 (The Reynolds Transport Theorem)
用来定量描述流体场中流体性质的变化情况
令 $V$ 表示 control volume (CV)，$S$ 表示 control surface (CS)，令 $b$ 表示一物理性质，有
$$
\frac{\mathrm d B_{\text{sys}}}{\mathrm d t} = \iiint_{\text{V}}\frac{\partial}{\partial t}(\rho b) \mathrm d V + \oiint_{\text{S}} \rho b \boldsymbol u \cdot \boldsymbol n \mathrm d S
$$
其中
$$
B_{\text{sys}} = \iiint_V \rho b \mathrm d V
$$
连续体方程 (Continuity Equation)
对于 CV 和 CS，考虑质量守恒
$$
m_{\text{in}} - m_{\text{out}} = \Delta m_{\text{CV}}
$$
移向可得
$$
\Delta m_{\text{CV}} + m_{\text{net flow out}} = 0
$$
对时间求导，则
$$
\frac {\partial m_{\text{CV}}}{\partial t} + \dot{m}_{\text{net flow out}} = 0
$$
对于 $m_{\text{CV}}$，有
$$
m_{\text{CV}} = \iiint_V \rho \mathrm d V
$$
对于净流出质量，在 CS 上积分可得
$$
\dot m_{\text{net flow out}} = \oiint_{S}\rho\boldsymbol u \cdot \mathrm d \boldsymbol S
$$
这里是流出的净质量是由于 CS 的法向量是向外定义的。
于是可以得到
$$
\frac{\partial }{\partial t}\displaystyle \iiint_V \rho \mathrm d V + \oiint_{S} \rho \boldsymbol u \cdot \mathrm d \boldsymbol S = 0
$$
由高斯定理，有
$$
\begin{aligned}
& \frac{\partial}{\partial t} \iiint_V \rho \mathrm d V + \iiint_{V} \nabla \cdot (\rho \boldsymbol u) \mathrm d V = 0
\end{aligned}
$$
对于可微流体场（不可压缩流体），令 $V \rightarrow 0$，有
$$
\frac{\partial \rho}{\partial t} + \nabla \cdot \rho \boldsymbol u = 0
$$
动量方程
利用动量守恒推导，考虑 CV 的受力，包括作用于整个体积的 body force （如重力，电磁力等）和作用于 control surface 的 surface force（如压力，黏滞力等）
若令 $\boldsymbol F$ 表示 body force，则总 body force 为
$$
\iiint_V \boldsymbol F \mathrm d V
$$
对于 surface force，分成 normal stress 和 shear stress 来考虑，令 $\boldsymbol \sigma$ 表示 shear stress tensor，$p$ 表示压强，$\boldsymbol n$ 为外法向量，则 surface force 为
$$
-\oiint_S p\boldsymbol n \cdot \mathrm d \boldsymbol S + \oiint_S \boldsymbol \sigma \cdot \mathrm d \boldsymbol S
$$
类似连续体方程中的净流出质量，第一项为负号是由于外法向量是向外的
于是 total force 为
$$
\Sigma_F =-\oiint_S p\boldsymbol n \cdot \mathrm d \boldsymbol S + \oiint_S \boldsymbol \sigma \cdot \mathrm d \boldsymbol S + \iiint_V \boldsymbol F \mathrm d V
$$
根据牛顿第二定律有
$$
\Sigma_F = m \boldsymbol a = \frac{\partial}{\partial t}\iiint_{\text{sys}}\rho \boldsymbol u \mathrm d V
$$
根据 Reynolds transport theorem，有
$$
\frac{\mathrm d B_{\text{sys}}}{\mathrm d t} = \int_{\text{CV}}\frac{\partial}{\partial t}(\rho b) \mathrm d V + \int_{\text{CS}} \rho b \boldsymbol u \cdot \boldsymbol n \mathrm d S
$$
将速度场作为属性 $b$ 代入，即 $B_{\text{sys}} = m\boldsymbol u, b = \boldsymbol u$，有
其中 $\otimes$ 为张量积（Kronecker Product），令 $\boldsymbol u = (u_x, u_y, u_z), \boldsymbol n = (n_x, n_y, n_z)$，可以得到
$$ \begin{aligned} \boldsymbol u \otimes \boldsymbol u \cdot \boldsymbol n &= \begin{bmatrix}u_x u_x & u_xu_y & u_x u_z \\ u_yu_x & u_yu_y & u_yu_z \\ u_zu_x & u_zu_y & u_zu_z \end{bmatrix} \cdot \boldsymbol n \\ & = \begin{bmatrix}n_x & n_y & n_z \end{bmatrix}\begin{bmatrix}u_x u_x & u_xu_y & u_x u_z \\ u_yu_x & u_yu_y & u_yu_z \\ u_zu_x & u_zu_y & u_zu_z \end{bmatrix}\\ & = \begin{bmatrix}u_xu_xn_x + u_xu_y n_y + u_xu_zn_z\\ u_xu_yn_x + u_yu_yn_y + u_yu_zn_z\\ u_xu_zn_x + u_yu_zn_y + u_zu_z n_z \end{bmatrix}^T \end{aligned} $$而
$$ \begin{aligned} \boldsymbol u(\boldsymbol u \cdot \boldsymbol n) &= \boldsymbol u (u_xn_x + u_yn_y + u_zn_z) \\ & = \boldsymbol u \otimes \boldsymbol u \cdot \boldsymbol n \end{aligned} $$于是我们建立起了动量方程：
$$
\frac {\partial }{\partial t} \iiint_V \rho \boldsymbol u \mathrm d V + \oiint_S \rho \boldsymbol u \otimes \boldsymbol u \cdot \mathrm d \boldsymbol S = -\oiint_S p\boldsymbol n \cdot \mathrm d \boldsymbol S + \oiint_S \boldsymbol \sigma \cdot \mathrm d \boldsymbol S + \iiint_V \boldsymbol F \mathrm d V
$$
可以理解为 自身动量变化 + 周围传递动量变化 = surface force + body force
同样，有高斯定理可以推出
$$
\frac {\partial }{\partial t} \iiint_V \rho \boldsymbol u \mathrm d V + \iiint_V \nabla \cdot (\rho \boldsymbol u \otimes \boldsymbol u)\mathrm d V = -\iiint_V \nabla \cdot (p\boldsymbol n) \mathrm d V + \iiint_V \nabla \cdot \boldsymbol \sigma \mathrm d V + \iiint_V \boldsymbol F \mathrm d V
$$
令 $V \rightarrow 0$，可得
$$
\frac {\partial \rho \boldsymbol u}{\partial t} + \nabla \cdot (\rho \boldsymbol u \otimes \boldsymbol u) = -\nabla p + \nabla \cdot \boldsymbol \sigma + \boldsymbol F
$$
若 body force 只考虑重力，则 $\boldsymbol F = \rho \boldsymbol g$
牛顿流体 (Newtonian flows)
对于牛顿流体
$$
\boldsymbol \sigma = 2 \rho \nu \boldsymbol S - \rho \nu’ (\nabla \cdot \boldsymbol u)\boldsymbol I
$$
其中 $\nu$ 为 kinematic viscosity，$\nu’$ 为 bulk viscosity，$\boldsymbol S$ 为 strain rate
$$
\boldsymbol S = \frac 1 2 \left(\nabla \boldsymbol u + \nabla \boldsymbol u^T \right)
$$
不可压缩性 (Incompressibility)
考虑牛顿流体的动量方程，且流体不可压缩 $\nabla \cdot \boldsymbol u = 0$，我们可以化简动量方程为
$$
\rho\left(\frac {\partial \boldsymbol u}{\partial t} + \boldsymbol u \cdot \nabla \boldsymbol u\right) = -\nabla p + \rho \nu \nabla^2 \boldsymbol u + \boldsymbol F
$$
Navier−Stokes Equation
对于上式，令 body force 只考虑重力，即 $\boldsymbol F = \rho \boldsymbol g$，再将等式两边同时除以 $\rho$ 并整理即得到了 N-S 方程
$$ \begin{aligned} \frac {\partial \boldsymbol u}{\partial t} + \boldsymbol u \cdot \nabla \boldsymbol u + \frac {1}{\rho} \nabla p &= \boldsymbol g + \nu \nabla ^2 \boldsymbol u\\ \nabla \cdot \boldsymbol u & = 0 \end{aligned} $$压力方程 (Pressure Equation)
对等式两边取散度，有
$$
0 + \rho \nabla \cdot (\boldsymbol u \cdot \nabla \boldsymbol u) = -\nabla^2 p + \rho \nu \nabla \cdot \nabla ^2 \boldsymbol u + \nabla \cdot \boldsymbol F
$$
整理，得
$$
\nabla^2 p = \rho (\nu \nabla \cdot \nabla ^2 \boldsymbol u - \nabla \cdot (\boldsymbol u \cdot \nabla \boldsymbol u)) + \nabla \cdot \boldsymbol F
$$
至此，其实已经可以求解了，但阶数太高，会因为数值问题导致求解存在困难。
固体墙边界条件 (Solid wall boundary treatment)
No-slip condition: 若固体不移动，流体在固体墙法线上的速度分量为 $0$，即
$$
\boldsymbol u \cdot \boldsymbol n = 0
$$
固体移动则流体速度在法线方向上的分量与固体移动速度在法线方向上的分量保持一致
半拉格朗日对流算法 (Semi-Lagrangian Advection)
分步求解
对于微分方程通常采用分步求解的方法，这里只举一个例子，比如
$$
\frac {\mathrm d q}{\mathrm d t} = f(q) + g(q)
$$
假设我们已知 $n$ 时刻的解 $q^n$，我们要求 $n + 1$ 时刻的解 $q^{n + 1}$，我们可以先求 $f$ 的影响
我们可以用泰勒展开来分析其准确性
$$ \begin{aligned} q ^ {n + 1} &= (q ^ n + \Delta t f(q ^ n)) + \Delta t g (q ^ n + \Delta t f(q ^ n))\\ & = q ^ n + \Delta t f(q ^ n) + \Delta t (g(q ^ n) + O(\Delta t))\\ & = q ^ n + \Delta t(f (q ^ n) + g(q ^ n) + O(\Delta t^2))\\ & = q ^ n + \frac {\mathrm d q}{\mathrm d t}\Delta t + O(\Delta t ^ 2) \end{aligned} $$Helmholtz-Hodge decomposition
我们不会正向的去求解 N-S 方程中的压力项，而是利用流体的无压缩性质将速度场投影到散度为 $0$ 的空间上，间接地计算压力项。
对于一个向量场 $\boldsymbol w$，我们可以将其分解为无散场 $\boldsymbol u$ 和有散场 $\nabla q$
$$
\boldsymbol w = \boldsymbol u + \nabla q \qquad \nabla \cdot \boldsymbol u = 0
$$
定义投影算子 $\mathrm P$，$\mathrm P(\boldsymbol w)$ 将向量场 $\boldsymbol w$ 投影到无散向量场 $\boldsymbol u$ 上
对于这个算子可以通过泊松方程来求解，对于上面的等式两边取散度可以得到
$$
\nabla \cdot \boldsymbol w = \nabla^2 q
$$
求解该泊松方程，然后就能得到 $\boldsymbol u$
$$
\boldsymbol u = \mathrm P(\boldsymbol w) = \boldsymbol w - \nabla q
$$
可以发现一些性质，对于无散场，其投影后仍为自己本身，对于梯度场，其投影后为 $0$。
N-S 方程变形 (Modified Navier-Stokes Equation)
对于 N-S 方程：
$$
\frac {\partial \boldsymbol u}{\partial t} = - \boldsymbol u \cdot \nabla \boldsymbol u - \frac {1}{\rho} \nabla p + \nu \nabla^2 \boldsymbol u + \boldsymbol f
$$
通过投影算子 $\mathrm P$ 消除压力项，将 $\mathrm P$ 同时作用于 N-S 方程两边，由上面的性质可知，左边不变而压力项变为 $0$，
$$
\frac {\partial \boldsymbol u}{\partial t} = \mathrm P(-\boldsymbol u \cdot \nabla \boldsymbol u + \nu \nabla ^2 \boldsymbol u + \boldsymbol f)
$$
于是我们只用求解这个式子
拆分变形后的 N-S 方程
需要计算的项为体积力项 (body force) $\boldsymbol f$，对流项 (advection) $-\boldsymbol u \cdot \nabla \boldsymbol u$，粘度项 $\nu \nabla ^ 2\boldsymbol u$，以及投影项 $\mathrm P$
令当前时刻速度场为 $\boldsymbol w_0(\boldsymbol x)$，第一步我们求解 force $\boldsymbol f$ 项 (Add force)：
$$
\boldsymbol w_1(\boldsymbol x) = \boldsymbol w_0(\boldsymbol x) + \Delta t f(\boldsymbol x, t)
$$
第二步计算对流项 (Advection)：
$$
\boldsymbol w_2(\boldsymbol x) = \boldsymbol w_1(\mathrm p(\boldsymbol x, -\Delta t))
$$
暂时先忽略对流项的具体计算方法
第三步计算粘度项 (Diffusion):
$$
\boldsymbol w_3(\boldsymbol x) = \boldsymbol w_2(\boldsymbol x) + \nu \Delta t \nabla ^2 \boldsymbol u
$$
移向整理得
$$
(\mathrm{\boldsymbol I} - \nu \Delta t \nabla^2\boldsymbol u)\boldsymbol w_3(\boldsymbol x) = \boldsymbol w_2(\boldsymbol x)
$$
可以通过求解线性方程组得到
第四步做投影
$$
\nabla ^2 q = \nabla \cdot \boldsymbol w_3
$$
通过线性方程组求解，然后得到
$$
\boldsymbol w_4 = \boldsymbol w_3 - \nabla q
$$
即为下一时刻的速度
计算对流项
考虑对流方程
$$
\frac {\mathrm D q}{\mathrm D t} = \frac {\partial q}{\partial t} + \boldsymbol u \cdot \nabla q = 0
$$
移向可得
$$
\frac{\partial q}{\partial t} = -\boldsymbol u \cdot \nabla q
$$
半拉格朗日法的思想正如其名字，它用拉格朗日视角去解决欧拉视角下的对流方程
假设我们要求的是 $n+ 1$ 时刻下位置 $\boldsymbol x$ 处的物理量 $q$ 的新值。考虑拉格朗日视角，我们可以寻找 $n + 1$ 时刻的上一时刻（即 $n$ 时刻），哪一个粒子在速度场 $\boldsymbol u$ 的作用下流向了 $\boldsymbol x$，于是
$$
q^{n+1}(\boldsymbol x) = q^n(\mathrm p(\boldsymbol x, -\Delta t))
$$
其中 $\mathrm p(\boldsymbol x, -\Delta t)$ 即表示了上一时刻流向该位置粒子的位置
对于该位置的求法，最简单的方法就是直接用前向欧拉法进行倒推 $\mathrm p(\boldsymbol x, -\Delta t) = \boldsymbol x - \Delta t \boldsymbol u(\boldsymbol x)$
更高阶精度可以考虑龙格库塔法等
这里我们要求即是
$$
\frac {\partial \boldsymbol u}{\partial t} = -\boldsymbol u \cdot \nabla \boldsymbol u
$$
于是就有了第二步中的
$$
\boldsymbol w_2(\boldsymbol x) = \boldsymbol w_1(\mathrm p(\boldsymbol x, -\Delta t))
$$
参考
- Fluid Mechanics: Fundamentals and Applications, Yunus A. Cengel, D. and Cimbala, J.M.
「Note」Basic Fluid Simulation