朱晓临，何红虹，郭清伟，周韵若

(合肥工业大学数学学院，安徽 合肥 230009)

流体模拟是计算机图形学的一个重要分支；模拟固体入水是流体模拟的一个重要研究课题，有着广泛地应用背景；例如，运动员跳水，水上飞机降落、水下航行器的投放等。其中，在固体入水的模拟过程中，流体与固体的相互作用，流体与固体交互时交界面的压力处理、流体表面的大变形都是研究的重点和难点。

光滑粒子流体动力学(smoothed particle hydrodynamics，SPH)方法在处理大变形和快速移动的界面具有独特的优势。因此，SPH方法被广泛应用于各种流-固耦合问题。但传统SPH方法在模拟运动固体边界时，由于固体边界附近流体粒子支持域被截断，流体粒子密度不连续，使得计算精度较差，会导致边界处压力振荡等问题；而流体表面因为受到固体的冲击力，流体液面较粗糙。

传统SPH方法处理固体边界的方法有边界力法、虚粒子法和耦合力法。文献[1]最早使用边界力法处理固壁边界，但该方法存在多个未知参数，需要根据具体问题进行调整；若调整不当，则可能导致流体粒子穿透固壁边界和粒子混乱等问题。为了防止流体粒子穿透固体边界，文献[2]使用边界力法对流体粒子施加较大的排斥力，并提供有限的边界处理，给方程引入了刚度参数；排斥力只在流体粒子穿透固体边界时起作用，粒子到边界的距离随时间变化，且粒子可能在边界处出现异常加速，不符合事实。文献[3]提出一种双向的边界力法，防止流体粒子穿透固体边界，采用修正的SPH 压力公式，避免了负压；在流固发生碰撞时，对碰撞点处的相对速度施加一个边界条件，但容易发生粒子堆积现象。为提高计算精度，解决边界处流体粒子支持域被截断的问题。文献[4]采用虚粒子法处理固体边界，利用固定虚粒子表征固壁边界，除不发生位移外，虚粒子与流体粒子属性完全相同，可出现少许流体粒子穿透固壁边界的现象。文献[5]在模拟水下爆炸等问题时，提出了一种镜像生成虚粒子的方法，该方法守恒性较好、精度高，但动态生成虚粒子较复杂。文献[6]利用虚粒子法模拟固壁边界，直接将内部流场压力插值到虚粒子上，虚粒子通过压力梯度对流场施加边界力；但是，当流体中出现压力动荡时，插值得到的边界虚粒子压力值可能为负，虚粒子会对流体产生非正常的吸引力，导致流体粒子穿透固壁边界。

针对边界力法和虚粒子法的不足，文献[7]提出耦合力法，将边界力与虚粒子相结合施加边界条件，防止流体粒子穿透固体边界，提高了计算精度。该方法能够获得精确的光滑压力场，适用于复杂或移动的固体边界，但是单独使用边界粒子表征固壁边界，可能会导致流体粒子穿透固体边界。文献[8]在模拟固体入水和浮体出水时，采用了耦合力法，对非均匀分布粒子进行密度校正和核梯度校正；但其采用高阶SPH 格式，虽精度较高，但计算复杂，模拟速度较慢，实时性较差。

针对流体表面较粗糙的问题，文献[9]采用各向异性核函数重构更为光滑、平整的流体表面，真实地展现出流体表面的细节；但由于确定粒子各向异性核函数非常复杂，该方法重构流体表面的速度较慢。针对传统各向异性核函数构造流体表面速度慢的缺陷，文献[10]提出流体内部粒子采用各向同性核函数，边界粒子采用各向异性核函数，提高了流体表面绘制的速度，但模拟的逼真度不够。文献[11]根据流体粒子的特征值比值将粒子分为近表面粒子和内部粒子，表面重构时近表面粒子参与计算,内部粒子根据邻居粒子的数量直接对色函数赋值；该方法计算效率与流体粒子数有关，粒子数增多时，计算速度较慢。文献[12]使用δ+-SPH方法，通过对自由液面粒子位移的修正，对液面的压力和光滑度进行处理；采用自适应粒子细化算法，解决粒子无序的问题；但是靠近边界的粒子容易产生负压。

本文针对固体边界附近流体粒子支持域被截断，流体粒子密度不连续，导致流体粒子在固体边界处出现穿透、堆积以及流体表面变形难以处理等问题；通过改进压力计算方法，使流体粒子在固体边界附近正常流动，采用耦合力法，改进软斥力，对流体粒子进行排斥，可以很好的稳定交互界面压力场；再对流体表面粒子位置进行校正，使流体表面粒子分布均匀，液面流动自然，提升了模拟效果。

1 SPH方法简介

1.1 SPH方法

SPH方法是一种自适应、无网格的拉格朗日型粒子法。由LUCY[13]与GINGOLD 和MONAGHAN[14]在解决天体物理学时提出，后被广泛应用于流体力学和固体力学中。通过离散带有物理属性值的粒子，由核近似和粒子近似计算得到标量A在位置r处的对应值为

其中，mj为粒子j的质量；rj为位置；ρj为密度；Aj为在rj处的场量；函数W(r,h)是半径为h的光滑核函数。由式(1)可知，其梯度和拉普拉斯算子只影响核函数，即

其中，符号∇为梯度算子；∇2为拉普拉斯算子。

核函数是SPH方法的核心，本文关于不同作用力的计算，采用文献[15]中提出的核函数。

1.2 控制方程

将固体看作受刚体运动约束的流体，统一纳入Navier-Stokes方程中求解。基于弱可压缩SPH方法，动量方程和连续性方程为

其中，ρ为密度；v为速度；p为压强；μ为黏度系数；F为额外力。

密度计算为

压强由下列状态方程计算得到

针对固体边界附近流体粒子支持域被截断的问题，为提高模拟方法的稳定性和解决固体边界处流体粒子密度不连续，流体粒子i的密度由周围所有邻居粒子密度之和计算得到。

本文采用文献[16]的方法进行密度校正，将边界粒子密度引入流固相互作用力的计算中。边界粒子的质量统一为mb，从而使边界粒子的体积不依赖于质量，保证了流体密度的稳定性，边界粒子的体积为

2 本文工作

在固体入水实验中，需要对3 种边界进行处理，即运动的固体边界，模拟实验需设置水槽两侧的静止固壁边界、以及流体表面的自由流动液面边界。针对传统方法的不足，需对上述3 种不同的边界分别进行处理。

2.1 运动固体边界的处理

通过改进耦合力法，结合边界力易于计算以及虚粒子可以消除固体边界附近流体粒子密度不连续的优势，处理运动的固体边界。如图1 所示，固体边界粒子(红色)，有自己的质量和密度，在固体边界粒子上布置排斥力，用于给流体粒子施加排斥力防止渗透，运动固体边界内布置2 层虚粒子 (黑色)，虚粒子之间的距离设置为0.5h，用于弥补流体到达固体附近时，流体粒子的支持域h被截断的问题，虚粒子基本属性与流体粒子相同，速度与固体相同。

图1 粒子分布初始状态 Fig.1 Initial state of the particle distribution

流体粒子i的压力由周围邻居粒子j(j为流体粒子或固体粒子)的压力通过插值得到，如图2(a)所示。运动固体边界粒子b的压力只受到流体粒子的影响，通过邻居流体粒子的压力插值得到，如图2(b)所示。

图2 黑色是固体粒子，白色是流体粒子((a)流体粒子受力方式；(b)固体粒子受力方式) Fig.2 Black is solid particles and white is fluid particles ((a) Forced method of fluid particle;(b) Forced method of solid particle)

采用文献[17]中压力和黏性力的计算公式，则流体粒子i的压力和黏性力分别为

流体压力很大程度上受固体粒子的质量影响，因此本文在流体粒子压力计算中加入固体粒子质量及密度。因为固体质量与流体质量差异较大，直接加入固体质量的计算，容易产生压力振荡，所以本文在计算压力时，对于固体边界粒子质量的贡献选择使用体积ibV代替，使边界压力计算不依赖于固体质量，但又受到固体边界粒子的影响，可以有效提高数值稳定性，避免边界处压力振荡引起的流体粒子在固体边界附近分离或波动的现象。对于流体粒子i，改进的压力计算式为

相应地，边界附近流体粒子黏性力的计算同样考虑固体粒子的贡献，改进的流体粒子黏性力为

其中，mb为固体的质量；g为重力加速度。

本文对文献[9]中的耦合力法进行改进，原固体边界粒子对流体粒子施加的排斥力模型为

由于原始耦合力法仅通过边界粒子表征固体边界，边界粒子不具有任何物理属性，会使流体粒子穿透边界。而在运动固体与流体的交互过程中，流场的运动受固体的影响较大。本文在计算排斥力时，添加固体粒子属性的计算，可以更好地平衡排斥力的大小，减轻压力扰动。对接近固体边界的流体粒子进行有限距离的排斥，有效防止流体粒子穿透固体边界。针对运动固体入水问题，改进软斥力模型为

其中，cs为色函数；χ,f (η)均为可变参数系数，通过粒子之间的位置和核半径的关系，判断粒子之间距离的函数[9]，即

其中，xij为流体粒子i到固体粒子j之间的向量；rij为2 个粒子之间的距离，Δd为2 个粒子的初始距离。

2.2 静止固壁边界的处理

对于固体落入静水，水槽两边固壁边界防穿透的设置，采用虚粒子法处理。在此，由于固体落入水槽底部时，固体粒子推开流体粒子后，流体粒子支持域被截断，水槽的底部边界处会出现流体粒子与固体边界分离现象，故在水槽底部设置2 层虚粒子，用于固体到达水槽底部时，补充底部边界处流体粒子的支持域。

在固壁边界外2h(h为核半径)距离内铺设虚粒子，方便下一步对流体表面粒子位置进行修正时，搜索表面修正区域内的流体粒子，补充两侧流体粒子的支持域，避免搜索出两侧的流体粒子，节省计算。固壁边界外虚粒子设置如图3 所示，圆内是支持域半径为2h的流体粒子i的支持域。

图3 固壁边界虚粒子分布示意图 (蓝色为流体粒子，黑色为固壁边界粒子，白色为虚粒子) Fig.3 Schematic diagram of the distribution of virtual particles at solid wall boundary (Blue is fluid particles,black is solid boundary particles and white is virtual particles)

2.3 流体表面粒子的位置修正

在固体入水问题中，流体表面的流体粒子与运动固体的相互作用通常与流体表面的变化和破裂有关。因为固体粒子对流体粒子施加了排斥力，流体表面会比未施加排斥力的表面更加粗糙，所以需要对流体表面粒子位置进行修正。

当固体与流体交互时，搜索表面的流体粒子，同时也搜索出流-固交互界面的粒子。最终流体表面和流-固交互界面2h(h是核半径)距离内的所有流体粒子，形成一个修正区域[12](图4)。并通过修正流体表面粒子的法向速度，将凸出表面的流体粒子当前时刻的法向速度设为零，保留粒子切线方向的速度，再利用当前时刻的法向速度，更新流体粒子下一时刻的速度，修正表面流体粒子的位置，提升流体表面流动效果。

图4 修正区域示意图 Fig.4 Schematic diagram of correction area

本文采用文献[18]中提出的距离场法搜索表面粒子。水槽的固壁边界外设置的虚粒子刚好补充了左右边界的搜索范围，而运动固体边界的虚粒子设置不满足2h距离，不需要单独进行交界面粒子的搜索，故搜索出来的流体粒子均为表面流体粒子，减少了计算量。

计算修正区域内流体粒子相对表面的法线向量为

其中，ni是提取的自由表面流体粒子i的法向。若|ni(r)|＞2h，则粒子i为表面流体粒子；否则为内部流体粒子。搜索出自由表面流体粒子后，令粒子i当前时刻法向速度为零，保持切线方向速度，流体粒子i方向速度如图5 所示。

图5 流体粒子i方向速度示意图(蓝色为流体粒子，黑色为固体粒子) Fig.5 Schematic diagram of the velocity in the direction of fluid particle i(Blue is a fluid particle and black is a solid particle)

于是，在t时刻，表面流体粒子的速度为

修正位置只修正表面飞散的流体粒子当前时刻的位置，而t+1 时刻的速度计算中，利用t时刻的法向速度继续更新，即

更新下一步流体表面粒子的位置为

其中，vi是流体粒子i在t时刻的速度；vi+1是流体粒子i在t+1 时刻的速度。

3 实验结果与分析

实验在Windows 10 平台进行，开发环境为visual studio 2013，渲染实验在Unity2018 上完成。实验配置为Intel(R) Core(TM) i5-5200U，2.20 GHz CPU，4 G 内存，NVIDIA Geforce GT 620 显卡。

通过经典的二维圆柱入水实验的模拟，实现了边界力法、虚粒子法、耦合力法及本文方法以不同速度下落的入水场景，验证了本文方法的有效性。实验时间步长取0.05 s，粒子支持域半径h取0.04 m。流体粒子个数为5 200 个，二维圆柱粒子个数为180 个。流体粒子初始密度取1 000 kg/m3，固体粒子密度取2 600 kg/m3。根据不同帧数下的对比，表明本文方法优于其他3 种方法。

二维圆柱置于空中1 m 处，分别以6.5 m/s 和20 m/s 的速度下落。实验初始设置如图6 所示。

图6 二维圆柱入水初始位置示意图 Fig.6 Schematic diagram of the initial position of the two-dimensional cylinder entering the water

3.1 二维圆柱以6.5 m/s 和20 m/s 速度入水的实验结果

从图7～12 可以看出，本文方法与其他3 种方法相比，流体粒子与固体交互界面较光滑。边界力 法流体表面较粗糙且流体粒子与固体粒子之间存在较大空隙，即流体因受力较大，导致流体粒子与固体分离程度较大；虚粒子法在整个模拟过程中有少许流体粒子穿透固体边界的现象；耦合力法中，流体与固体的交互界面较粗糙，且边界附近出现压力振荡。本文方法在整个模拟过程中较稳定，固体边界处流体粒子流动较自然，界面分离清晰，无穿透或分离程度较大的现象。从表1和2 可以看出，本文方法的模拟速度较其他3 种方法快。

表1 二维圆柱以6.5 m/s 速度入水各方法在 不同帧数下的用时比较(s) Table 1 Comparison of the time consumption of the two-dimensional cylinder entering the water at a speed of 6.5 m/s under different frames (s)

图7 6.5 m/s 速度下各方法在第65 帧的模拟效果对比 ((a)边界力法；(b)虚粒子法；(c)耦合力法；(d)本文方法) Fig.7 Comparison of simulation effects of each method at frame 65 at a speed of 6.5 m/s ((a) Boundary force method; (b) Virtual particle method;(c) Coupling force method; (d) Method of this article)

图8 6.5 m/s 速度下各方法在第87 帧的模拟效果对比 ((a)边界力法；(b)虚粒子法；(c)耦合力法；(d)本文方法) Fig.8 Comparison of simulation effects of each method at frame 87 at a speed of 6.5 m/s ((a) Boundary force method; (b) Virtual particle method;(c) Coupling force method; (d) Method of this article)

图9 6.5 m/s 速度下各方法在第115 帧的模拟效果对比 ((a)边界力法；(b)虚粒子法；(c)耦合力法；(d)本文方法) Fig.9 Comparison of simulation effects of each method at frame 115 at a speed of 6.5 m/s ((a) Boundary force method;(b) Virtual particle method;(c) Coupling force method; (d) Method of this article)

图10 20 m/s 速度下各方法在第25 帧的模拟效果对比 ((a)边界力法；(b)虚粒子法；(c)耦合力法；(d)本文方法) Fig.10 Comparison of simulation effects of each method at frame 25 at a speed of 20 m/s ((a) Boundary force method; (b) Virtual particle method;(c) Coupling force method; (d) Ethod of this article)

图11 20 m/s 速度下各方法在第28 帧的模拟效果对比 ((a)边界力法；(b)虚粒子法；(c)耦合力法；(d)本文方法) Fig.11 Comparison of simulation effects of each method at frame 28 at a speed of 20 m/s ((a) Boundary force method; (b) Virtual particle method;(c) Coupling force method; (d) Method of this article)

图12 20 m/s 速度下各方法在第31 帧的模拟效果对比 ((a)边界力法；(b)虚粒子法；(c)耦合力法；(d)本文方法) Fig.12 Comparison of simulation effects of each method at frame 31 at a speed of 20 m/s ((a) Boundary force method; (b) Virtual particle method;(c) Coupling force method; (d) Method of this article)

表2 二维圆柱以20 m/s 速度入水各方法在 不同帧数下的用时比较 (s) Table 2 Comparison of the time consumption of the two-dimensional cylinder entering the water at a speed of 20 m/s under different frames (s)

3.2 小球入水渲染实验

本文分别渲染了以6.5 m/s 和20 m/s 速度下落的小球入水实验，小球半径为0.04 m，密度为2 600 kg/m3。并给出了不同速度下小球处于同一高度的实验效果图(图13～16)。

图13 小球位于0.45 m 高处((a) 6.5 m/s 速度第36 帧；(b) 20 m/s 速度第25 帧) Fig.13 The ball is at a height of 0.45 m ((a) Frame 36 at 6.5 m/s;(b) Frame 25 at 20 m/s)

图14 小球位于0.20 m 高处 ((a) 6.5 m/s 速度第44 帧；(b) 20 m/s 速度第28 帧) Fig.14 The ball is at a height of 0.20 m ((a) Frame 44 at 6.5 m/s;(b) Frame 28 at 20 m/s)

图15 小球落至底部((a) 6.5 m/s 速度第75 帧；(b) 20 m/s 速度第31 帧) Fig.15 The ball falls to the bottom ((a) Frame 75 at 6.5 m/s;(b) Frame 31 at 20 m/s)

图16 小球落至底部((a) 6.5 m/s 速度第80 帧；(b) 20 m/s 速度第65 帧) Fig.16 The ball falls to the bottom ((a) Frame 80 at 6.5 m/s;(b) Frame 65 at 20 m/s)

从渲染实验效果图可以看出，本文方法在模拟固体入水时，可以呈现出固体落水的完整过程，且模拟效果真实自然，细节展现较好。

4 结 论

本文通过对运动固体边界压力计算方法和耦合力法的改进，有效补充了边界处流体粒子的支持域、稳定压力场，防止流体粒子穿透固体边界等问题。进一步通过对流体表面粒子位置的修正，对流体液面进行优化，提升了模拟效果。通过模拟固体入水的二维实验和渲染实验，对提出的方法进行了验证，结果表明了本文方法的可行性与有效性。