SPH求解器¶

基本原理¶

光滑粒子流体动力学（Smoothed Particle Hydrodynamics, SPH）是一种无网格Lagrangian方法，通过核函数插值来求解偏微分方程。

核函数近似¶

对于任意函数 $A(\mathbf{r})$，SPH插值定义为：

$$A(\mathbf{r}) = \sum_{b} A_b W(\mathbf{r} - \mathbf{r}_b, h)$$

其中： - $W$ 是光滑核函数 - $h$ 是光滑长度 - $b$ 遍历邻域粒子

核函数¶

PhysEngine支持Wendland C2核函数：

$$W(q, h) = \alpha_d \begin{cases} (1 - q/2)^2 (1 + 2q) & 0 \leq q \leq 2 \ 0 & q > 2 \end{cases}$$

归一化系数： - 1D: $\alpha_1 = 4/3$ - 2D: $\alpha_2 = 15/7\pi$ - 3D: $\alpha_3 = 21/16\pi$

算法流程¶

flowchart TD
    A[初始化粒子位置] --> B[时间步长估计]
    B --> C[八叉树邻居搜索]
    C --> D[计算右端项 RHS]
    D --> E[状态方程更新]
    E --> F[边界条件处理]
    F --> G[位置/速度更新]
    G --> H{收敛?}
    H -->|否| B
    H -->|是| I[输出结果]

状态方程¶

采用Tait方程描述流体状态关系：

$$P = B \left[ \left( \frac{\rho}{\rho_0} \right)^\gamma - 1 \right]$$

参数说明： - $B = \rho_0 c_0^2 / \gamma$ - $c_0$ 为参考声速 - $\gamma$ 通常取7

人工粘度¶

为处理激波，引入Lagrangian人工粘度：

$$\Pi_{ab} = \begin{cases} \frac{-\alpha \bar{c}{ab} \mu} + \beta \mu_{ab}^2}{\bar{\rho{ab}} & \mathbf{v} < 0 \ 0 & \text{otherwise} \end{cases}$$} \cdot \mathbf{r}_{ab

其中 $\mu_{ab} = h \frac{\mathbf{v}{ab} \cdot \mathbf{r}$}}{|\mathbf{r}_{ab}|^2 + \epsilon h^2

核心方程¶

密度更新（连续方程）¶

$$\frac{d\rho_a}{dt} = \rho_a \sum_b \frac{m_b}{\rho_b} (\mathbf{v}a - \mathbf{v}_b) \cdot \nabla_a W$$

动量方程¶

$$\frac{d\mathbf{v}a}{dt} = -\sum_b m_b \left( \frac{P_a}{\rho_a^2} + \frac{P_b}{\rho_b^2} + \Pi$$} \right) \nabla_a W_{ab} + \mathbf{g

时间步长¶

基于CFL条件和声波传播条件：

$$\Delta t = \min\left( \text{CFL} \cdot \frac{h}{c_0 + v_{\max}}, 0.25 \cdot \frac{h}{\sqrt{\nu_{\max}}} \right)$$