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_{ab} + \beta \mu_{ab}^2}{\bar{\rho}_{ab}} & \mathbf{v}_{ab} \cdot \mathbf{r}_{ab} < 0 \\ 0 & \text{otherwise} \end{cases}\]

其中 \(\mu_{ab} = h \frac{\mathbf{v}_{ab} \cdot \mathbf{r}_{ab}}{|\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_{ab}\]

动量方程¶

\[\frac{d\mathbf{v}_a}{dt} = -\sum_b m_b \left( \frac{P_a}{\rho_a^2} + \frac{P_b}{\rho_b^2} + \Pi_{ab} \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)\]