基于速度修正的固壁边界处理方法

朱晓临，陈 薇

基于速度修正的固壁边界处理方法

朱晓临，陈 薇

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

固壁边界处理方法的研究一直是流体模拟中的难点问题，常见的固壁边界处理方法有边界力法和虚粒子法。边界力法通过对靠近边界的流体粒子施加作用力防止其穿透边界，但模型参数较多，力的大小难以调控，且在计算中会产生边界截断误差问题；虚粒子法通过在边界外生成虚粒子解决了边界截断误差问题，但在处理复杂边界问题时，由于外部的虚粒子的生成较困难，且分布不均，计算精度受到影响，出现粒子飞散的情况。针对2种方法存在的问题，提出一种基于速度修正的固壁边界处理方法，无需求解边界力或在边界外生成虚粒子，直接利用动量方程和计算速度耗损量求出流体粒子碰到边界后的反弹速度，大大降低了处理边界的复杂程度，也克服了2种方法在边界拐角处粒子不均匀采样而导致的算法不稳定的问题。模拟仿真验证了该方法在稳定性、计算效率方面均较传统边界力法和虚粒子法更好；随着粒子数的增加，该方法耗时更少、计算效率更高，对复杂场景的模拟效果更好。

流体模拟；固壁边界；边界力法；虚粒子法；速度修正

光滑粒子流体动力学(smooth particle hydrodynamics, SPH)方法是流体模拟，尤其是液体模拟的主要方法。流体模拟的固壁边界处理方法主要包括边界力法和虚粒子法2种。

边界力法应用一层边界粒子来表示固壁边界，通过该层边界粒子对邻近边界的流体粒子施加作用力，防止流体粒子穿透固壁边界。其关键是壁面作用力的施加模型。

1994年，MONAGHAN[1]最先提出通过施加边界力的方式进行固壁边界处理，并提出兰纳-琼斯(Lennard-Jones，L-J)势函数作为作用力的施加模型，但该模型中存在许多未知参数，需要根据具体问题进行调整，调整不当可能会导致流体粒子穿透边界或粒子飞乱等问题。2004年，MÜLLER等[2]基于L-J势函数实现了在微可压缩情况下，将SPH方法用于模拟流体与可变形固体的相互作用，并且在作用力较大时，流体粒子不会穿透固体边界，但需要更小的时间步长来确保模拟的稳定性。2007年，HARADA等[3]针对边界粒子不足的问题，提出一种边界粒子密度权重函数，将边界粒子的密度加入到流体密度的计算中，从而减小了因边界粒子的不足造成的计算误差。2009年，BECKER等[4]以直接力的方法处理固壁边界，同时采用预测-校正技术计算粒子的速度和位置。上述方法仅适合规则固体边界的情况，且在计算邻近边界粒子的属性时由于边界外测粒子的缺失导致计算会产生截断误差。

2009年，MONAGHAN和KAJTAR[5]采用一种简化的边界力计算公式克服了L-J势函数方法计算边界力参数较多的不足，但力的大小和作用范围仍难以控制。2012年，LIU等[6]提出了一种新型边界力模型，并将边界粒子与虚粒子结合施加边界条件，但若单独使用一层边界粒子表示固壁边界，仍会出现流体粒子穿透边界的问题。边界力法可以方便地表征复杂的边界，且实现方法相对简单，易于编程。然而，主要缺点为边界力的范围和强度难以控制，若边界力太大，则会导致边界附近粒子属性的数值计算波动较大造成粒子飞乱；反之则将导致流体粒子穿透边界。

虚粒子法应用多层虚粒子表征固壁边界，通过插值或镜像方式将流体粒子的部分属性(质量、密度、压强等)赋予虚粒子，虚粒子与流体粒子通过流体力学控制方程相互作用。虚粒子法的关键是其属性值的获取及与流体粒子的相互作用。

1996年，RANDLES和LIBERSKY[7]提出通过镜像方式生成虚粒子进行边界处理，如果流体粒子到边界的距离小于其支持域半径，那么在固壁外就会自动生成镜像虚粒子，流体粒子和镜像虚粒子的速度大小相等但方向相反，其他物理量相同。1997年，MORRIS和MONAGHAN[8]使用插值虚粒子法模拟SPH流体与固壁边界的相互作用。2002年，LIU等[9]提出在固壁边界处布置2种虚粒子，第一种是固定在边界上，且对靠近边界的流体粒子施加一个排斥力，从而防止流体粒子穿透边界；第二种是由邻近边界的流体粒子对称产生的，并赋予其相应物理量，从而使得边界附近的流体粒子在进行积分计算时不会产生截断误差。2006年，HU和ADAMS[10]使用镜像虚粒子法模拟流体在规则固体边界条件下的运动，虚粒子设置在流体粒子作用域与边界重叠的部分，并且虚粒子与流体粒子具有相同的密度、质量、压力和黏度，将虚粒子看作是流体粒子进行处理。上述方法只适用于规则固壁直边界，对于复杂边界情况，虚粒子的生成过于复杂，计算量较大。

2012年，SCHECHTER和BRIDSON[11]提出了一种新镜像粒子法来模拟流体的自由表面和进行边界处理。该方法在液体和固体之间以及液体和空气之间产生镜像粒子，并使用液体粒子通过插值运算推导出镜像粒子的属性值。该方法解决了存在于流体自由表面附近的非自然张力视觉形变问题和在固壁边界处无法真实模拟的流体的黏附现象。2015年，刘虎等[12]提出了应用一系列的虚粒子来表征固壁边界的处理方法，当求解方程式时，让虚粒子具有流体粒子的部分属性(密度、质量、压力等)，求解其连续方程和状态方程，并将边界问题耦合到方程的求解过程中，无需提出显式的边界条件。该方法的优点是无需对边界进行特殊的处理，缺点是通过求解连续方程和状态方程得到的虚拟粒子的压力，与真实压力之间存在误差，导致模拟的不稳定，出现粒子飞散的情况。虽然虚粒子法的守恒性较好，精度也较高，但在生成自由表面的虚粒子时过于复杂，计算时间大大增加，固壁边界形状复杂时更加如此。

总体而言，边界力法适合处理复杂的固体边界，然而模型的参数较多，力的大小难以调控，并且在计算过程中会产生边界截断误差问题；虚粒子法能够很好地解决边界截断误差，但在处理复杂边界问题时，边界外部的虚粒子的生成较困难，计算精度和稳定性难以保证，会出现粒子飞乱甚至程序运行崩溃的情况。

针对上述问题，本文提出一种基于速度修正的固壁边界处理方法，无需在边界外生成虚粒子或求解边界力，直接利用动量方程和计算速度耗损量求出流体粒子遇到边界后的反弹速度，大大降低了处理边界的复杂程度，也克服了虚粒子法和边界力法在拐角处粒子不均匀采样而导致的算法不稳定的问题。模拟仿真验证了本文方法在稳定性、计算效率方面均比传统方法更好；随着粒子数的增加，本文方法的耗时增长慢于传统方法，计算效率高的优势更明显；对复杂场景的模拟效果更好。

1 SPH方法简介

1.1 SPH基本公式

在SPH方法中，函数()定义为

其中，为光滑函数的影响和支持范围即光滑长度；则称为光滑核函数(smoothing kernel function)或光滑函数(smoothing function)，简称为核函数(kernel)。

本文采用三次样条的核函数

其中，=/；α为函数在维空间中的归一化系数，其在一维、二维和三维空间中的取值分别为1/，15/(7π2)，3/(2π3)。

通过将式(3)离散化，用光滑函数的支持域内所有粒子叠加求和的形式表示任一点处的粒子近似式为

其中，为粒子的支持域内的粒子总数；m，，分别为支持域范围内粒子的质量、密度和位置；(–,)为粒子对粒子产生影响的光滑函数，且该函数与光滑长度有关。

1.2 控制方程

本文用到的是拉格朗日形式的Navier-Stokes方程(简称N-S方程)

本文采用SPH方法的计算过程流程如图1所示。

图1 SPH方法的流程图

使用SPH方法求解N-S方程，对于任一粒子，其密度为

求得粒子密度后，还需要计算粒子压强，进而求得粒子的压力。

利用泰特方程求压强，即

其中，为常数，一般取=7；为参数，用于限制密度的最大改变量；0为参照密度。

SPH方法中求解粒子所受的压力和黏性力的表达式分别为

流体粒子的速度、位置等信息是根据其所受合力计算所得，利用式(10)和式(11)，将N-S方程中的动量守恒方程式(7)转化为

2 基于速度修正的固壁边界处理方法

文献[1]最早提出应用边界力的方式施加流固边界条件，并提出一种基于距离的L-J势函数的方法施加作用力

其中，系数1一般取12，2一般取4；的取值一般为流体最大速度的平方量级；0为光滑核半径；r为流体粒子到固体边界上位置点的距离；为流体粒子位置到固体边界上位置点之间的向量。

文献[7]在处理直壁边界问题时，能够很好地解决边界截断误差，但在处理圆弧表面、倾斜壁面等边界问题时，边界外部的镜像粒子的生成较困难，并且镜像粒子分布不均，计算精度也会受影响，出现粒子跃出甚至程序崩溃的情况。

目前，已有的边界处理方法都很难兼顾计算精度、边界复杂度以及计算效率。本文提出一种新的基于速度修正的固壁边界处理方法：先在边界上设置一个阻尼区，当流体粒子运动到阻尼区内时，且流体粒子速度与边界法向夹角大于π/2时，将该区域分为4个小区域，并分别对流体速度作不同处理：将流体粒子速度从，坐标系转化为，坐标系表示(其中，分别为边界的法向与切向)；然后用流体粒子与相关的边界粒子的相对位置计算速度耗损量；之后将不同区域的流体粒子反弹到相应区域；最后将处理后的流体粒子速度转回，坐标系。具体步骤如下：

(1) 添加阻尼区。以边界斜坡为例(图2)：在斜坡边界(黑色线条)上方添加阻尼区边界(蓝色线条)，阻尼区的宽度为，其数值与粒子的支持域长度(即核函数的光滑长度)相同。

图2 边界阻尼区示意图

(2) 流体粒子到边界的距离。在边界上等距分布粒子，边界粒子间距与初始流体粒子间距有关。流体粒子到边界的距离定义为：流体粒子到边界粒子的2个最短距离的平均值作为该流体粒子到边界的距离，如图3所示。

图3 流体粒子到边界粒子的距离

(3) 划分区域。当流体粒子在阻尼区内，且其速度方向与边界法向夹角大于π/2时，将该区域分为①、②、③、④4个小区域，如图4所示，其中，①~④为第1~4个小区域。

依本文方法，区域①中的粒子将反弹到区域⑧中；区域②中的粒子反弹到区域⑦中；区域③的粒子反弹到区域⑥中；区域④的粒子反弹到区域⑤中。在流体粒子碰到边界反弹的过程中，其反弹速度会受到周围流体粒子的影响，故需要计算速度的耗损量。

(4) 计算速度耗损量。在流体粒子周边的边界粒子数为，记为1,2,···,，如图5所示。

图5 流体粒子与边界粒子的夹角

由图5得

其中，为流体粒子到边界的距离(式(14))，r为粒子到粒子α的距离(=1,2,···)。

计算粒子速度耗损量为

当粒子到达边界时会发生反弹，各种复杂边界均可用类似的方法进行处理。

将本文方法与传统边界力法和虚粒子法进行对比，可以发现：①在稳定性方面，边界力法是通过基于距离的L-J函数来计算邻近边界的流体粒子受到的排斥力，当流体粒子运动到拐角处时，因受到来自相同距离边界粒子施加的力的影响导致粒子的反弹速度的计算不稳定；虚粒子法在处理复杂边界问题时，边界外部的镜像虚粒子的生成较困难，且分布不均匀，计算精度受影响。本文方法可对边界的形状做检测，通过计算速度耗损量对速度进行修正，可以更好地控制粒子在拐角处运动的稳定性，克服了传统方法在边界拐角处粒子不均匀采样而导致的算法不稳定的问题。②在计算效率方面，与边界力法相比，本文方法无需使用较多参数且通过较复杂的边界力模型来计算力，而是直接对速度进行修正，通过计算速度耗损量求出流体粒子遇到边界后的反弹速度，减小了计算量；与虚粒子法相比，本文方法无需在边界外部生成虚粒子，降低了处理边界的复杂程度，尤其在处理复杂边界问题时，虚粒子的生成较困难，且在计算流体粒子的速度及位置信息时，虚粒子也参与其中，使得计算量大大增加。因此，本文方法在稳定性和计算效率方面均要优于传统的边界力法和虚粒子法。

3 实 验

3.1 二维溃坝模拟

以二维溃坝实验模拟斜坡边界和弧形边界场景，验证本文方法在稳定性和计算效率方面的有效性。相关参数：流体粒子的初始密度为1 000 kg/m3，粒子初始间距为0.002 92 m，时间步长为10×10–3s，粒子的支持域长度=0.006 m。

(a) 边界力法(b) 虚粒子法(c) 本文方法

由图6可知，图6(a)为边界力法，稳定性略差，粒子撞到边界后，上爬升过程中出现轻微的粒子飞散现象；图6(b)为虚粒子法，稳定性差，粒子撞到边界后，上爬升过程中粒子飞散情况明显；图6(c)为本文方法，在粒子撞到右侧边界后，爬升时的稳定性好，没有出现明显的粒子飞散情况，与实际中水撞到固壁边界后爬升的效果更为接近。

上述3种方法模拟实验进行流体粒子总数分别为2 500个、5 000个和7 500个粒子，每组实验进行5次取运行到2 000帧时的平均值作为运行时间，模拟二维溃坝斜坡边界运行时间的结果见表1。

表1 3种方法模拟二维溃坝斜坡边界的运行时间

由表1和图7可以看出，本文方法明显比传统的边界力法和虚粒子法的计算效率高，并且随着粒子数的增加，耗时增长也比其他2种方法慢，计算效率高的优势更明显。

图7 随着粒子数的增加，边界力法、虚粒子法和本文方法耗时增长对比

同表1模拟实验条件，3种方法模拟二维溃坝弧形边界的运行时间见表2。

由图8和图9及表2可知，本文方法在稳定性和计算效率方面要优于传统边界力法及虚粒子法，并且随着粒子数的增加，本文方法的耗时增长慢于其他2种方法，计算效率更高。

(a) 边界力法(b) 虚粒子法(c) 本文方法

表2 3种方法模拟二维溃坝弧形边界的运行时间

3.2 复杂场景模拟

图10和图11为3种方法模拟液体在漏斗和弓形水管的运行场景以及边界处理效果。

由图10和图11可知，图10(c)为本文方法模拟漏斗的场景，上壁未出现明显的粒子飞散情况；图11(c)为模拟弓形水箱场景，流体流动的整体稳定性较好，且粒子分布更均匀。说明本文方法在模拟复杂场景时的稳定性更好，适用范围更广。

图9 随着粒子数的增加，边界力法、虚粒子法和本文方法耗时增长对比

(a) 边界力法(b) 虚粒子法(c) 本文方法

(a) 边界力法(b) 虚粒子法(c) 本文方法

4 结 论

针对传统的边界力法和虚粒子法的不足，本文提出一种基于速度修正的固壁边界处理方法，既不需要求解边界力也无需在边界外生成虚粒子，直接利用动量方程和计算速度耗损量求出流体粒子遇到边界后的反弹速度，大大降低了处理边界的复杂程度，也克服了边界力法和虚粒子法在边界拐角处粒子不均匀采样而导致的算法不稳定的问题。模拟仿真结果表明本文方法的稳定性、计算效率均比传统的2种方法更好，且随着粒子数的增加，本文方法的耗时增长也比上述2种方法慢，计算效率高的优势更明显；通过模拟一些复杂场景，本文方法未出现粒子飞散情况，稳定性好，适用范围更广。

[1] MONAGHAN J J. Simulating free surface flows with SPH [J]. Journal of Computational Physics, 1994, 110(2): 399-406.

[2] MÜLLER M, SCHIRM S, TESCHNER M, et al. Interaction of fluids with deformable solids [J]. Computer Animation and Virtual Worlds, 2004, 15(34): 159-171.

[3] HARADA T, KOSHIZUKA S, KAWAGUCHI Y. Smoothed particle hydrodynamics in complex shapes [J/OL].[2018-10-11].https://dl.acm.org/citation.cfm?id=2614375.

[4] BECKER M, TESSENDORF H, TESCHNER M. Direct forcing for Lagrangian rigid-fluid coupling [J]. IEEE Transactions on Visualization and Computer Graphics, 2009, 15(3): 493-503.

[5] MONAGHAN J J, KAJTAR J B. SPH particle boundary forces for arbitrary boundaries [J]. Computer Physics Communications, 2009, 180(10): 1811-1820.

[6] LIU M B, SHAO J R, CHANG J Z. On the treatment of solid boundary in smoothed particle hydrodynamics [J]. Science China (Technological Sciences), 2012, 55(1): 244-254.

[7] RANDLES P W, LIBERSKY L D. Smoothed particle hydrodynamics: Some recent improvements and applications [J]. Computer Methods in Applied Mechanics and Engineering, 1996, 139(1-4): 375-408.

[8] MORRIS J P, MONAGHAN J J. A switch to reduce SPH viscosity [J]. Journal of Computational Physics, 1997, 136(1): 41-50.

[9] LIU M B, LIU G R, LAM K Y. Investigations into water mitigation using a meshless particle method [J]. Shock Waves, 2002, 12(3): 181-195.

[10] HU X Y, ADAMS N A. A multi-phase SPH method for macroscopic and mesoscopic flows [J]. Journal of Computational Physics, 2006, 213(2): 844-861.

[11] SCHECHTER H, BRIDSON R. Ghost SPH for animating water [J]. ACM Transactions on Graphics, 2012, 31(4): 1-8.

[12] 刘虎, 强洪夫, 陈福振, 等. 一种新型光滑粒子动力学固壁边界施加模型[J]. 物理学报, 2015, 64(9): 384-397.

Treatment of Solid Boundary Based on Velocity Correction

ZHU Xiao-lin, CHEN Wei

(School of Mathematics, Hefei University of Technology, Hefei Anhui 230009, China)

The study of solid boundary treatment method has been a difficult problem in fluid simulation. The common methods of solid boundary treatment are boundary force method and virtual particle method. The boundary force method prevents fluid particles from penetrating the boundary by exerting force on the fluid particles near the boundary, but this method has too many parameters and the force is difficult to adjust, and the boundary truncation error will occur in the calculation. The virtual particle method solves the problem of boundary truncation error by generating virtual particles outside the boundary. However, when dealing with complex boundary problems, the generation of virtual particles is difficult, and the calculation accuracy will be affected due to the uneven distribution of virtual particles, which leads to particles to disperse. To solve these problems, the paper presents a new method for the treatment of solid boundary based on velocity correction, which does not need to solve the boundary force or generate virtual particles outside the boundary. The momentum equation and the velocity consumption are used directly to calculate the rebound velocity of the fluid particles when they hit the boundary, which greatly reduces the complexity of the boundary treatment. It also overcomes the problem of the instability of the boundary force method and the virtual particle method caused by the uneven sampling of particles at the corner of the boundary. Simulation results show that the proposed method is more stable and more efficient than traditional methods above, and with the increase of the number of particles, the time consuming of this method is also slower than that of the two methods, and the advantage of high computational efficiency is more obvious and the simulation effect of complex scene is also better.

fluid simulation; solid boundary; boundary force method; virtual particle method; velocity correction

TP 391.9

10.11996/JG.j.2095-302X.2019040637

A

2095-302X(2019)04-0637-07

2019-01-30；

定稿日期：2019-03-20

国家自然科学基金项目(61272024)

朱晓临(1964-)，男，安徽池州人，教授，博士，硕士生导师。主要研究方向为数值逼近、计算机图形学、CAGD、图形图像处理等。 E-mail：zxl_hfut@126.com