TheCodeWay:我的博客

        和其他流体力学中的数学方法类似，SPH算法同样涉及到“光滑核”的概念，可以这样理解这个概念，粒子的属性都会“扩散”到周围，并且随着距离的增加影响逐渐变小，这种随着距离而衰减的函数被称为“光滑核”函数，最大影响半径为“光滑核半径”。

        反过来不难理解，尽管我们将流体视为一个个分散的粒子，但流体毕竟是连续充满整个空间的，流体中每个位置参与运算的值都是由周围一组粒子累加起来的。


        设想流体中某点\(\vec{r}\)（此处不一定有粒子）,在光滑核半径h范围内有数个粒子，位置分别是，\(\vec{r_0},\vec{r_1},\vec{r_2},\ldots\vec{r_j}\)，则该处某项属性A的累加公式为：\begin{equation}
A(\vec{r})=\sum_j{A_j\frac{m_j}{\rho_j}W(\vec{r}-\vec{r_j}, h)}\tag{3.1}
\end{equation}        其中\(A_j\)是要累加的某种属性，\(m_j\)和\(\rho_j\)是周围粒子的质量和密度，\(\vec{r}\)是该粒子的位置，\(h\)是光滑核半径。函数\(W\)就是光滑核函数。
        光滑核函数两个重要属性，首先一定是偶函数，也就是\(W(-r)=W(r)\)，第二，是“规整函数”，也就是\(\int{W(r)dr}=1\)

---

SPH推导过程

        我们假设流体中一个位置为\(\vec{r_i}\)的点，此处的密度为\(\rho(r_i)\)、压力为\(p(r_i)\)、速度为\(\vec{u}(r_i)\)，那么我们可以根据上一篇的公式2.8，可以推导出此处的加速度\(\vec{a}(r_i)\)为\begin{equation}
\vec{a}(r_i)=\vec{g}-\frac{\nabla{p(r_i)}}{\rho(r_i)}+\frac{\mu\nabla^2\vec{u}(r_i)}{\rho(r_i)}\tag{3.2}
\end{equation}        对于SPH算法来说，基本流程就是这样，根据光滑核函数逐个推出流体中某点的密度，压力，速度相关的累加函数，进而推导出此处的加速度，从而模拟流体的运动趋势，下面我们逐个来分析

---

密度

        根据公式3.1，用密度\(\rho\)代替\(A\)，可以得到\begin{equation}
\rho(r_i)=\sum_j{\rho_j\frac{m_j}{\rho_j}W(\vec{r_i}-\vec{r_j},h)}=\sum_j{m_jW(\vec{r_i}-\vec{r_j},h)}\tag{3.3}
\end{equation}        计算使用的光滑核函数称为Poly6函数，具体形式为:\begin{equation}
W_{Poly6}\ (\vec{r},h)=\begin{cases}
K_{Poly6}\ (h^2-r^2)^3 &, 0\le r\le h \\[2ex]
0 &, \text otherwise \\
\end{cases}\\ \text{with}\ r=|\vec{r}|\tag{3.4}
\end{equation}        其中\(K_{Poly6}\ \)是一个固定的系数，根据光滑核的规整属性，通过积分计算出这个系数的具体值，在2D情况下，在极坐标中计算积分：\begin{equation}
K_{Poly6}=1/{\int_0^{2\pi}\int_0^h{r(h^2-r^2)^3}dr\ d\theta}=\frac{4}{\pi h^8}\tag{3.5}
\end{equation}        3D情况下，在球坐标中计算:\begin{equation}
K_{Poly6}=1/{\int_0^{2\pi}\int_0^{\pi}\int_0^h{r^2sin(\varphi)(h^2-r^2)^3}dr\ d\varphi\ d\theta}=\frac{315}{64\pi h^9}\tag{3.6}
\end{equation}        由于所有粒子的质量相同都是\(m\)，所以在3D情况下，\(\vec{r_i}\)处的密度计算公式最终为：\begin{equation}
\rho(r_i)=m\frac{315}{64\pi h^9}\sum_j{\left(h^2-|\vec{r_i}-\vec{r_j}|^2\right)^3}\tag{3.7}
\end{equation}

---

压力

        根据上一节的结论，在位置\(r_i\)之处的由压力产生的作用力的计算公式为\begin{equation}
\vec{F_i}^{pressure}=-\nabla p(\vec{r_i})=-\sum_j{p_j\frac{m_j}{\rho_j}\nabla W(\vec{r_i}-\vec{r_j},h)}\tag{3.8}
\end{equation}        不过不幸的是，这个公式是“不平衡”的，也就是说，位于不同压强区的两个粒子之间的作用力不等，所以计算中一般使用双方粒子压强的算术平均值代替单个粒子的压力\(r_i\)之处的由压力产生的作用力的计算公式为\begin{equation}
\vec{F_i}^{pressure}=-\sum_j{\frac{m_j(p_i+p_j)}{2\rho_j}\nabla W(\vec{r_i}-\vec{r_j},h)}\tag{3.9}
\end{equation}        对于单个粒子产生的压力\(p\),可以用理想气体状态方程计算\begin{equation}
p=K(\rho-\rho_0)\tag{3.10}
\end{equation}        其中\(\rho_0\)是流体的静态密度，\(K\)是和流体相关的常数，只跟温度相关。
        压力计算中使用的光滑核函数称为Spiky函数\begin{equation}
W_{Spiky}\ (\vec{r},h)=\begin{cases}
K_{Spiky}\ (h-r)^3 &, 0\le r\le h \\[2ex]
0 &, \text otherwise \\
\end{cases}\\ \text{with}\ r=|\vec{r}|\tag{3.11}
\end{equation}        在3D情况下，\(K_{Spiky}=15/(\pi h^6)\)\begin{equation}
\nabla W_{Spiky}\ (\vec{r},h)=\frac{15}{\pi h^6}\nabla(h-r)^3=-\vec{r}\frac{45}{\pi h^6 r}(h-r)^2\tag{3.12}
\end{equation}        将公式3.12带入3.9，可以整理出公式3.2中压力产生的加速度部分\begin{equation}
\vec{a_i}^{pressure}=-\frac{\nabla p(\vec{r_i})}{\rho_i}=m\frac{45}{\pi h^6}\sum_j{\left(\frac{p_i+p_j}{2\rho_i\rho_j}(h-r)^2\frac{\vec{r_i}-\vec{r_j}}{r}\right)}\\ \text{with}\ r=|\vec{r_i}-\vec{r_j}|\tag{3.13}
\end{equation}

---

粘度

        现在把注意力集中到公式3.2中最后一部分，由粘度产生的作用力\begin{equation}
\vec{F_i}^{viscosity}=\mu\nabla^2\vec{u}(r_i)=\mu\sum_j{\vec{u_j}\frac{m_j}{\rho_j}\nabla^2W(\vec{r_i}-\vec{r_j},h)}\tag{3.14}
\end{equation}        这个公式同样有不平衡的问题，考虑到公式中的速度其实并不是绝对速度，而是粒子间的相对速度，所以这个公式的正确写法应该是：\begin{equation}
\vec{F_i}^{viscosity}=\mu\sum_j{m_j\frac{\vec{u_j}-\vec{u_i}}{\rho_j}\nabla^2W(\vec{r_i}-\vec{r_j},h)}\tag{3.15}
\end{equation}        其中的光滑核函数形式如下：\begin{equation}
W_{viscosity}\ (\vec{r},h)=\begin{cases}
K_{viscosity}\ \left(-\cfrac{r^3}{2h^3}+\cfrac{r^2}{h^2}+\cfrac{h}{2r}-1\right) &, 0\le r\le h \\[3ex]
0 &, \text otherwise \\
\end{cases}\\ \text{with}\ r=|\vec{r}|\tag{3.16}
\end{equation}        在3D情况下，\(K_{viscosity}=15/(2\pi h^3)\)\begin{equation}
\nabla^2 W_{viscosity}\ (\vec{r},h)=\nabla^2\frac{15}{2\pi h^3}\left(-\frac{r^3}{2h^3}+\frac{r^2}{h^2}+\frac{h}{2r}-1\right)=\frac{45}{\pi h^6}(h-r)\tag{3.17}
\end{equation}        由此可得到公式3.2的粘度部分\begin{equation}
\vec{a_i}^{viscosity}=\frac{\vec{F_i}^{viscosity}}{\rho_i}=m\mu\frac{45}{\pi h^6}\sum_j\frac{\vec{u_j}-\vec{u_i}}{\rho_i\rho_j}(h-\mid\vec{r_i}-\vec{r_j}\mid)\tag{3.18}
\end{equation}

---

        把公式3.13和3.17带入3.2，可以得到，对于粒子i，它的加速度可以由下面的公式计算\begin{equation}
\vec{a}(r_i)=\vec{g}+m\frac{45}{\pi h^6}\sum_j\left(\frac{p_i+p_j}{2\rho_i\rho_j}(h-r)^2\frac{\vec{r_i}-\vec{r_j}}{r}\right)+m\mu\frac{45}{\pi h^6}\sum_j\frac{\vec{u_j}-\vec{u_i}}{\rho_i\rho_j}(h-r)\\ \text{with}\ r=|\vec{r_i}-\vec{r_j}| \tag{3.19}
\end{equation}        好了，我们似乎推导出一大推复杂的公式，不用担心，你已经过了最困难的部分，下一节我们来点真的，让这些公式运行起来看看