和其他流体力学中的数学方法类似,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}        好了,我们似乎推导出一大推复杂的公式,不用担心,你已经过了最困难的部分,下一节我们来点真的,让这些公式运行起来看看

20 thoughts on “SPH算法简介(三): 光滑核函数

  1. “不过不幸的是,这个公式是“不平衡”的,也就是说,位于不同压强区的两个粒子之间的作用力不等,所以计算中一般使用双方粒子压强的算术平均值代替单个粒子的压力”
    请问为什么一定要“平衡”?

    1. 这里的平衡是指根据牛顿的第三定律,粒子之间的作用力和反作用力应该是相等的

    1. 这些核函数背后的数学原理需要了解相应的流体力学知识才行,我只是玩票而已,这些公式都是从其他资料上抄来的,所以你要是有兴趣的话,还得自己学一下。这个系列是我几年前写的,当时忘了把参考资料记录下来了,不好意思

    1. 你好,我也是SPH方法的初学者,请问可以加个QQ好友吗,方便以后互相交流

  2. 你好,博主。
    我想问下在程序初始化时,粒子的尺寸 m_unitScale 起到什么作用?
    我看在计算粒子间距和碰撞时 都有考虑粒子尺寸。

  3. 博主,你好。感谢你的分享,写的太好了。我有个问题:
    公式3-16下面的,K(viscosity)在2D情况下是多少呢?

  4. 博主你好,请问计算密度的公式在向量ri-向量rj的模接近于h时会接近于0,而后又被作为除数,导致两个粒子在接近的一瞬间会被各自施加一个很大的速度,请问如何避免呢

    1. 实际运算中总是以某个实际的粒子为中心计算密度的,所以最极端的情况是这个粒子远离所有其他的粒子,那么就会只算上这个粒子自身,密度为315m/(64*pi*h^3),也不会是0,所以你担心的情况不会出现

  5. 十分感谢博主分享,是篇十分好的流体入门博客,源码对我在unity下重现给了很关键的帮助;
    再次谢谢。

  6. 请问累加公式中为什么m还要除以ρ,粒子的属性为什么还有ρ呀?粒子本身包含物质的多少的信息,那物质在空间中的分布信息不是也可以确认吗?就是说空间中ρ可以确认

发表评论

邮箱地址不会被公开。

This site uses Akismet to reduce spam. Learn how your comment data is processed.