边界截断问题的解决

边界截断是SPH模拟的一个重要问题。如果支撑域被边界截断，就会使得粒子属性近似值偏低。SPHinXsys是如何解决的呢？本文将分壁面附近、自由表面、开放边界、压力边界、流壳耦合几个场景分别阐述。

壁面附近的边界截断

注意：这里的壁面不是壳体。对于这种情况，应对方法是最简单的，用户保证壁面至少要有三层粒子即可。因为默认使用的Wendland Quintic（C2）核的默认截断半径是$2.6\Delta L$，所以壁面三层粒子一般都能够让壁面附近流体粒子有完整的支撑域。

自由表面的边界截断

液体与空气（未模拟）接触的界面通常视为自由液面。

density summation

对于这种液面，在做density summation的时候要使用FreeSurface类型的，例如DensitySummationComplexFreeSurface。它会强制令流体的密度大于参考密度：

$$\langle\rho_i\rangle=\max\left\{\rho_0\frac{\sum_j W_{ij}}{\sum_j W_{ij}^{0}}, \rho_0\right\}$$

transport velocity correction

如果要做transport velocity correction，ParticleScope须为BulkParticles。

参考案例

test_2d_dambreak，test_2d_water_entry_exit。

开放边界的边界截断

开放边界包括自由流（free stream）和非周期性的入口出口（不包括压力边界）。

density summation

在做density summation的时候，要使用FreeStream类型的，例如DensitySummationFreeStreamComplex。它会把

$$\rho_i^\mathrm{sum}=\rho_0\frac{\sum_j W_{ij}}{\sum_j W_{ij}^{0}}$$

$$\langle\rho_i\rangle=\begin{cases} \rho_i^\mathrm{sum}+\max\{0,\rho_i-\rho_i^\mathrm{sum}\}\rho_0/\rho_i,& i\text{ has free surface neighbor(s)};\\ \rho_i^\mathrm{sum},& \text{else}. \end{cases}$$

等价于

$$\langle\rho_i\rangle=\begin{cases} \rho_i^\mathrm{sum},& i\text{ has no free surface neighbor or }\rho_i < \rho_i^\mathrm{sum};\\ \rho_i^\mathrm{sum}(1-\rho_0/\rho_i)+\rho_0,& \text{else}.\\ \end{cases}$$

对于表面粒子或者近表面粒子，一般会有

1. $\rho_0>\rho_i$，于是$\langle\rho_i\rangle<\rho_0$。

2. $\rho_i > \rho_i^\mathrm{sum}$，于是$\langle\rho_i\rangle > \rho_i^\mathrm{sum}$。

也即$\rho_i^\mathrm{sum} < \langle\rho_i\rangle < \rho_0$。所以我们可以认为FreeStream类型的density summation将边界粒子的$\langle\rho_i\rangle$拉高了一些，但没有像FreeSurface类型那样直接拉高到参考密度，因此这是一种更软的处理。

transport velocity correction

如果要做transport velocity correction，ParticleScope须为BulkParticles。

参考案例

test_2d_free_stream_around_cylinder（自由流）、test_3d_poiseuille_flow_shell（速度入口）。

压力边界的边界截断

压力边界也是开放边界，但是SPHinXsys有专门的处理。

压力梯度计算

使用PressureBoundaryCondition<TargetPressure>。它会在动量方程中添加一项$2 p_\mathrm{b} \sum_j m_j/(\rho_i \rho_j) \nabla W_{i j}$，以弥补边界截断给压力梯度带来的损失。

具体理论见Zhang et al（2025）。假想buffer区域之外有一堆粒子，如果这些假想的粒子存在，那么buffer区域的粒子的支撑域就完整了。可这些假想粒子实际上是缺失的。作者通过一些巧妙的数学变换，使得在这些粒子缺失的情况下也能正确考虑它们的影响。

首先，根据梯度的反对称性，有以下近似（前提是粒子均匀分布）：

$$\sum_j \frac{m_j}{\rho_j}\nabla W_{ij}\approx \boldsymbol{0}.$$

对buffer区域的粒子来说，如果算上边界外缺失的粒子，那么上式依然成立。我们记buffer区域的粒子下标为$j$，缺失粒子的下标为$k$，于是

$$\sum_j \frac{m_j}{\rho_j} \nabla W_{i j}+\sum_k \frac{m_k}{\rho_k} \nabla W_{i k} \approx \boldsymbol{0}.$$

$$\sum_k \frac{m_k}{\rho_k} \nabla W_{i k}\approx -\sum_j \frac{m_j}{\rho_j} \nabla W_{i j}.\tag{1}$$

SPHinXsys使用Riemann solver进行物理场的稳定。对于buffer区域粒子，带有Riemann solver的压力梯度离散形式为

$$\nabla p_i=2 \sum_j P_{i j}^* \frac{m_j}{\rho_j} \nabla W_{i j}+2 \sum_k P_{i k}^* \frac{m_k}{\rho_k} \nabla W_{i k},$$

其中$P^*$是任意两个粒子对的压力Riemann solution。

对于边界外的粒子，我们指定它具有恒定的压力$p_\mathrm{b}$（第一类边界条件），也即$P_{i k}^*=p_\mathrm{b}$。结合式(1)，可以得出

$$\nabla p_i=2 \sum_j P_{i j}^* \frac{m_j}{\rho_j} \nabla W_{i j}-2 p_\mathrm{b} \sum_j \frac{m_j}{\rho_j} \nabla W_{i j} .$$

将其代入动量方程的离散形式：

$$\begin{aligned} \frac{d \mathbf{v}_i}{d t}= & -2 \sum_j m_j\left(\frac{P_{i j}^*}{\rho_i \rho_j}\right) \nabla W_{i j}+2 p_\mathrm{b} \sum_j\left(\frac{m_j}{\rho_i \rho_j}\right) \nabla W_{i j} \\ & +2 \sum_j m_j \frac{\eta \mathbf{v}_{i j}}{\rho_i \rho_j r_{i j}} \frac{\partial W_{i j}}{\partial r_{i j}}+\mathbf{g} . \end{aligned}\tag{2}$$

注意，动量方程中多出的一项是PressureBoundaryCondition干的第一件事（也是最重要的）：

$$2 p_\mathrm{b} \sum_j\left(\frac{m_j}{\rho_i \rho_j}\right) \nabla W_{i j}$$

它干的第二件事，是令速度只有边界的法向分量，没有切向分量，也就是做了个投影：

$$\mathbf{v}_i=\left(\mathbf{v}_i \cdot \hat{u}\right) \hat{u},$$

其中$\hat{u}$是buffer区域边界的单位法向量。

density summation

使用Pressure版本的density summation，例如DensitySummationPressureComplex。它只对内部流动区域的做density summation。buffer区域的流体使用target pressure计算密度（已在BidirectionalBuffer中实现）。

transport velocity correction

如果要做transport velocity correction，ParticleScope须为BulkParticles。

参考案例

test_2d_pulsatile_poiseuille_flow。

流壳耦合中的边界截断

壳是单层没错，但每个壳粒子携带了厚度 Thickness，并在邻域/核积分时用若干法向 dummy 层把“壳背后那一侧缺失的邻域”补成“体积分”。所以对流体粒子而言，靠近壳面时并不会出现“核支撑域空一半”的问题——支撑域是通过这种等效体积贡献被补齐的，而不是靠真实的多层离散。

具体可以看neighborhood.cpp跟shell相关的代码，这里就不展开了。

参考案例

test_2d_channel_flow_fluid_shell。

参考文献

Zhang S., Fan Y., Wu D., Zhang C., Hu X.; Dynamical pressure boundary condition for weakly compressible smoothed particle hydrodynamics. Physics of Fluids 2025; 37 (2): 027193. https://doi.org/10.1063/5.0254575