Smoothed-particle hydrodynamics

1. SPH Kernel

常用的sph kernel函数如下式所示

$$W(r)=c \left\{\begin{matrix} 1-6(\frac{r}{h})^2 + 6(\frac{r}{h})^3 &(0 \le \frac{r}{h} < \frac{1}{2})\\ 2(1-\frac{r}{h})^3 &(\frac{1}{2} \le \frac{r}{h} < 1) \\ 0 &else \end{matrix}\right.$$

c为归一化系数，对于1d sph 有:

$$I_1 =\int _0^{\frac{h}{2}} 1 - \frac{6}{h^2}r^2 + \frac{6}{h^3} r^3 dr + \int_\frac{h}{2}^{h} 2(1-\frac{r}{h})^3 dr\\ = (r-\frac{2}{h^2} r^3 + \frac{3}{2h^3} r^4)|_{0}^{\frac{h}{2}} - h\frac{(1-\frac{r}{h})^4}{2}|^{h}_{\frac{h}{2}}\\ = \frac{h}{2} - \frac{h}{4} + \frac{3h}{32} + \frac{h}{32} \\ =\frac{3}{8}h\\ c_1 = \frac{1}{2I_1}=\frac{4}{3h}$$

对于2d sph有:

$$I_2= \int_0^{\frac{h}{2}} (1 - \frac{6}{h^2}r^2 + \frac{6}{h^3}r^3)rdr + \int_{\frac{h}{2}}^h (1-\frac{1}{h}r)^3rdr=A+B \\ A = \int_0^\frac{h}{2} r-\frac{6}{h^2}r^3+\frac{6}{h^3} r^4 dr\\ = \frac{1}{2}r^2-\frac{3}{2h^2}r^4+\frac{6}{5h^3}r^5 |_0^{\frac{h}{2}}\\ =\frac{1}{8}h^2-\frac{3}{32} h^2+\frac{6}{160}h^2=\frac{11}{160}h^2\\ B=\int_{\frac{h}{2}}^{h}2(1-\frac{3}{h}r+\frac{3}{h^2}r^2 - \frac{1}{h^3}r^3)rdr\\ =\int_{\frac{h}{2}}^{h}2r-\frac{6}{h}r^2 +\frac{6}{h^2}r^3-\frac{2}{h^3}r^4dr\\ =r^2 -\frac{2}{h}r^3+\frac{3}{2h^2}r^4-\frac{2}{5h^3}r^4|_{\frac{h}{2}}^{h}\\ =h^2-2h^2+\frac{3}{2}h^2-\frac{2}{5}h^2-\frac{1}{4}h^2+\frac{1}{4}h^2-\frac{3}{32}h^2+\frac{6}{160}h^2\\ =\frac{3}{160}h^2\\ \\ c_2=\frac{1}{\int_{0}^{2\pi} I_2d\theta}=\frac{40}{7\pi h^2}$$

2d sph kernel的代码为:

@ti.func
def sph_kernel(x : vec2, x0 : vec2):
    r = (x - x0).norm()
    r_div_h = r / h
    r_div_h2 = r_div_h * r_div_h
    r_div_h3 = r_div_h2 * r_div_h

    rv = 0.0
    
    if r_div_h < 1 / 2:
        rv = 1 - 6 * r_div_h2 + 6 * r_div_h3
    elif r_div_h < 1:
        one_m_r_div_h = (1 - r_div_h)
        rv = 2 * one_m_r_div_h * one_m_r_div_h * one_m_r_div_h
    else:
        rv = 0.0

    return rv * c

对2d sph kernel求梯度有

$$\nabla W=\frac{\partial W(r,h)}{\partial x} + \frac{\partial W(r,h)}{\partial y}\\ =\frac{\partial W(r,h)}{\partial r} (\frac{\partial r}{\partial x} + \frac{\partial r}{\partial y})\\ \frac{\partial W}{\partial r} = c \cdot \left\{\begin{matrix} -\frac{12}{h}(\frac{r}{h}) + \frac{18}{h}(\frac{r}{h})^2 &(0 \le \frac{r}{h} < \frac{1}{2})\\ -\frac{6}{h}(1-\frac{r}{h})^2 &(\frac{1}{2} \le \frac{r}{h} < 1) \\ 0 &else \end{matrix}\right.\\ r^2=x^2+y^2 \to \frac{\partial r}{\partial x}=\frac{x}{r},\frac{\partial r}{\partial y} = \frac{y}{r}$$

python代码为:

@ti.func
def gradient_sph_kernel(x: vec2, x0: vec2):
    x_m_x0 = x - x0
    r = x_m_x0.norm()

    drdx = x_m_x0[0] / (r + 4e-3)
    drdy = x_m_x0[1] / (r + 4e-3)

    r_div_h = r / h

    dSdr = 0.0

    if r_div_h < 1 / 2:
        dSdr = 18 / h3 * r * r - 12 / h2 * r
    elif r_div_h < 1:
        one_m_r_div_h = (1 - r_div_h)
        dSdr = -6 * one_m_r_div_h * one_m_r_div_h / h 
    else:
        dSdr = 0.0


    return vec2(dSdr * drdx, dSdr * drdy) * c

Sph用采样点的加权平均作为函数估计值

$$f(x)\approx \sum_i f(x_i) V_i W(r_i)\\ r_i = ||x-x_i||_2$$

对微分算子 ,由于Sph中只有 与 有关，因此

$$\nabla f\approx \sum_i f(x_i) V_i \cdot \nabla W(r_i)=\sum_i \frac{m_i}{\rho_i}f(x_i)\nabla W(r_i)$$

这对任意微分运算成立

通常，在估计梯度时，我们常常会用sph的拓展非对称形式和对称形式。由密度公式

非对称:

$$\nabla (f\cdot 1) =\nabla f \cdot 1 + f\cdot \nabla 1\\ \sum f(x_i)V_i\nabla W(r_i)=\nabla f +f(x)\sum V_i \nabla W(r_i)\\ \nabla f=\sum V_i\nabla W \cdot (f(x_i)-f(x))\\ =\sum m_i\nabla W \frac{f(x_i)-f(x)}{\rho_i}$$

对称:

$$\nabla (f\cdot q^n)=\nabla f \cdot q^n + f\cdot nq^{n-1} \nabla q\\ \nabla f = q^{-n}\nabla (f\cdot q^n)-f\cdot nq^{-1} \nabla q\\ =q^{-n} \sum f(x_i)q_i^n V_i \nabla W-nfq^{-1} \sum q_i V_i\nabla W$$

令

$$\nabla f = \rho \sum f(x_i) \rho_i^{-1} V_i \nabla W + f\rho^{-1} \sum \rho_iV_i\nabla W\\ =\sum m_i\frac{f(x_i)\rho}{\rho_i^2} \nabla W + \sum m_i\frac{f(x)}{\rho}\nabla W\\ =\rho \sum m_i(\frac{f(x_i)}{\rho_i^2} + \frac{f(x)}{\rho^2}) \nabla W$$

这两种形式不存在绝对的优劣。

一般在计算压力，作用力相关的物理量的时候会用对称形保证系统动量守恒。

2. wcsph weakly compressible smooth particle hydrodynamics

​ 2.1 理论

​ 在模拟流体模拟时，往往需要求解带不可压约束的navier-stokes方程

$$\frac{Dv}{Dt}=g-\frac{1}{\rho}\nabla p+\mu \nabla^2 v (1)\\ \nabla \cdot v = 0 (2)$$

​ 在传统方法中，为了满足不可压约束往往会同过联立方程(1),(2)同时求解 同时更新压力和速度。具体做法为通过spliting把微分方程(1)分为两个部分 advection和projection，逐步求解方程。这些方程被称为advection-projection方法。其中projection涉及把空间离散化后建立稀疏矩阵并求解。

$$\frac{Dv}{Dt} = g+\mu \nabla ^2 v (advection)\\ \left\{\begin{matrix} \frac{\partial v}{\partial t} = -\frac{1}{\rho}\nabla p \\ \nabla \cdot v = 0 \end{matrix}\right. (projection)$$

​ 此时，压强可以看成流体为了满足约束(2)产生的内力。这个过程中 为常量。

​ 可以看到，由于不可压约束,advection-projection多了求解线性系统这一复杂步骤。wcsph将不可压条件改为弱可压条件，避免了复杂计算。它将 删去，将压强建模为 ,其中 为rest状态下密度。方程变为，其中 为变量。

$$\frac{dv}{dt}=g-\frac{1}{\rho} \nabla p+ \mu \nabla^2 v$$

​ wcsph将空间离散化为一系列带属性的小粒子(定值粒子质量 ,粒子体积 变量：密度,速度,位置) 。该方法利用拉格朗日视角追踪粒子，通过sph方法采样空间中各种属性用以更新速度。算法为(其中):

1) 遍历所有粒子估计密度

$$\rho_i=\sum_{neighbor} \frac{m_j}{\rho_j} \rho_j W_{ij}=\sum_{neighbor}m_jW_{ij}$$

1. 遍历所有粒子估计压强

$$p_i = max(k_1(\frac{\rho_i}{\rho_0})^{k_2} - k_1, 0)$$

为了解决free surface边界条件，将 限制为大于0可以防止边缘的粒子密度太小导致压强为负

1. 遍历所有粒子计算方程其余项:

$$\frac{1}{\rho_i}\nabla p_i=\frac{1}{\rho_i} \cdot \rho_i \sum_{neighbor}m_i (\frac{p_i}{\rho_i^2}+\frac{p_j}{\rho_j^2})\nabla W \\ =\sum_{neighbor} m_i(\frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2})\nabla W \\ \mu \nabla^2 v=\mu\cdot 2(d+2) \sum_{j=0}^n \frac{m_j}{\rho_j}\left(\frac{v_{i j} \cdot x_{i j}}{\left\|x_{i j}\right\|^2+0.01 h^2}\right) \nabla W_{i j}$$

这里用了 的拟合化形式 。

​ 4. 遍历所有粒子，带入方程，更新速度

$$\frac{dv}{dt} = g-\frac{1}{\rho} \nabla p + \mu \nabla^2 v$$

2.2 踩坑

模拟过程中参数的一组合理值

$$\begin{matrix} \rho_0 = 1 & radius=0.25 & support\_radius(h)=4\cdot radius=1 &volume=\pi\cdot radius^2 =\frac{16}{\pi}\\ mass=volume\cdot \rho_0 =\frac{\pi}{16} & g=1 & dt\in(10^{-4},2\times10^{-4}) & \mu \in (0.05, 1)\\ k_1=10 & k_2=7 \end{matrix}$$

时间步长对wcsph模拟稳定性至关重要。当 过大时会出现数值爆炸现象：

3 DFSPH (Divergence-Free Smooth Particle Hydrodynamics)

DFSPH提出了一种更稳定的流体模拟方法，它在wcsph给出的压力模型上修改压力计算方式，保证了系统的密度恒定以及速度散度恒定，从而使系统满足不可压约束。

*.1 Rigid-Fluid Coupling

Versatile Rigid-Fluid Coupling for Incompressible SPH提出了一种能同时模拟刚体和sph流体的方法，其方法核心为给刚体生成模拟粒子。在模拟过程中同时计算流体粒子和刚体粒子相互的作用力。

对于任意形状的刚体，作者用 Particle-Based Simulation of Granular Materials 中的方法生成粒子。具体做法为:

1. 对于刚体表面（如三角网格），计算其sdf用于生成不同偏移参数下的粒子。
2. 利用marching tetrahedra(2d 下为marching squre)，为sdf表面生成三角形网格(2d下为线段)
3. 对每个三角形(2d下为线段),在其内部等概率随机生成 个粒子其中 为密度参数 为粒子半径, 为三角形面积 (2d下为 其中 为线段长度)
4. 调整粒子位置让其在网格表面分布均匀。做法为以速度 更新粒子位置。其中 ， 为刚体的sdf场， 为sdf在 位置的法线。 为核函数(可选用sph函数)。 把更新后的粒子约束在sdf表面， 使粒子之间的距离尽可能的远。

利用该方法生成的圆形为:

在模拟刚体过程中，作者采用如下形式计算两个粒子之间的压强力:

$$F^p = -m_i\frac{1}{\rho_i} \sum_j m_j \frac{p_j-p_i}{\rho_j}\nabla W_{ij}\\ =-m_i\frac{1}{\rho_i}\sum_j m_j\frac{p_j}{\rho_j}\nabla W_{ij} +\frac{1}{\rho_i}\sum_jm_j\frac{p_i}{\rho_j}\nabla W_{ij}\\ 则 F^p_{i\leftarrow j}=m_i \frac{1}{\rho_i} \cdot m_j \frac{p_j}{\rho_j} \nabla W_{ij}=m_im_j \frac{p_j}{\rho_i\rho_j}\nabla W_{ij}$$

对于刚体粒子j和流体粒子i。由于连续性，不妨令 和 则

$$F^p_{i\leftarrow j} = m_i m_j \frac{p_i}{\rho_i ^2} \nabla W$$

对于 我们采用 计算。其中 为流体粒子在reference状态下的密度。 在由预计算提前计算。

虽然我们在生成粒子的过程中尽可能的保证了粒子的均匀分布。由于刚体表面的几何形式，仍会导致粒子在空间分布的不均匀。因此，在计算粒子体积时，可对密集区域的粒子分配较小体积，对稀疏区域粒子分配较大体积。体积的形式为:

$$V_j = \frac{m_j}{m_j\sum_k W_{jk}}=\frac{1}{\sum_k W_{jk}}$$

其中 为刚体粒子k的邻域刚体粒子的sph核函数的加权平均。

对于处于液体刚体交界处的液体粒子i的压强 。其密度 需改写为

$$\rho_{f_i} = m_{f_i} \sum_j W_{ij}+\sum_k M_k (\rho_0^i)W_{ik}$$

其中j为粒子i的液体粒子邻居，k为粒子i的刚体粒子邻居,

在实现时，若采用symetric form sph则其形式为:

$$F^p = -m_i\sum_j m_j (\frac{p_j}{\rho_j^2}+\frac{p_i}{\rho_i^2})\nabla W_{ij}\\ F^p_{i\leftarrow j}=-m_im_j\frac{p_j}{\rho_j^2}\nabla W_{ij}\\ 令p_j = p_i, \rho_j = \rho_i\\ F_{i\leftarrow j} ^p=-m_iM_j(\rho_0^i)\frac{p_i^2}{\rho_i^2}\nabla W_{ij}\\ F_p=-m_i\sum_j M_j(\rho_0^i)(\frac{p_i}{\rho_i^2} + \frac{p_i}{\rho_i^2})\nabla W_{ij}$$

实现后效果如下图所示，可以看到除了少部分高速粒子漏入刚体内部，大部分粒子都与刚体表面交互并被阻挡。

TODO: 摩擦力

对于刚体受力，可以累加各个刚体粒子受力,并将这些粒子的力和力矩累加作为刚体受力结果。

$$F_j = \sum_i F_{j\leftarrow i}^p+\sum_i F_{j\leftarrow i} ^v=\sum_i F_{i\leftarrow j}^p+\sum_i F_{i\leftarrow j} ^p\\ F = \sum_j F_j\\ T = \sum (x-x_j) \times F_j$$