SPH流固耦合模拟自适应采样固壁虚粒子边界处理方法

朱晓临，汪欢欢，殷竞存



SPH流固耦合模拟自适应采样固壁虚粒子边界处理方法

朱晓临，汪欢欢，殷竞存

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

固壁虚粒子边界处理方法是流体模拟中一种主要边界处理方法，但其不能确保流体粒子不穿透固体边界，并且计算量较大。为防止流体粒子穿透边界，在边界附近设置一个阻尼区，阻尼区内的流体粒子被边界施加一个弹性力和一个和流体粒子运动速度方向相反的阻尼力，使得边界附近流体粒子更加稳定。为减少计算量，提出两种边界粒子自适应采样法：一种是依据边界周围粒子数目的不同，边界粒子自适应地采样质量不同的大小粒子；另一种是依据边界周围粒子数目的不同，边界粒子自适应的采样不同层数的相同质量粒子。与传统的固体边界粒子采样方法相比，该方法减少了边界粒子数目，加快了模拟速度，节省了计算机内存，基于GPU加速技术实现的三维流体模拟，能够进行实时交互。

固壁虚粒子；阻尼区；弹性力；阻尼力；边界粒子自适应采样；GPU加速技术

流体和固体的交互在日常生活中是普遍存在的，例如石头掉进水里、衣服浮在水面上、冰块的融化、海绵吸水等。流固耦合动力学涉及到流体动力学和固体动力学，对于流固耦合问题中的固体研究总体上可以分为刚体和可变形的物体；流体主要分为气体与液体。从总体来看，流固耦合问题按照耦合机理可以分为两类：第一类是指流体与固体的相互作用只发生在流体与固体的交界面上，即固体不具有渗水性，如石头和水的作用；第二类是指流体区域与固体区域部分或全部重叠在一起，如渗流问题。在计算机图形学领域，流固耦合问题的模拟主要是基于物理的偏微分方程求解方法的数值模拟。主要的流固耦合数值模拟方法有有限元法、边界元法、光滑粒子流体动力学(smoothed particle hydrodynamics，SPH)方法等。然而，有限元法需要全局离散，导致问题的自由度和原始信息量大；边界元法只在边界上离散问题域，降低了问题的复杂度；有限元法和边界元法都是欧拉网格方法，难以模拟具有较大变形的流固耦合自然现象。SPH方法是一种拉格朗日无网格插值方法，能方便地处理具有大变形流固耦合问题。对于SPH方法进行流固耦合模拟，边界处理是研究的重点和难点。对于SPH流固耦合的边界处理，目前主要有4类方法：基于距离的排斥力方法、固壁虚粒子法、虚粒子法、预测矫正法前3类方法是基于固体分别给予流体不同的力而实现耦合效应；第4类方法是计算机图形学领域中常用的正法。基于距离的排斥力方法计算的固体流体之间的力不是物理力，导致模拟时粒子堆积问题；固壁虚粒子法由于将固体看成由流体粒子组成，能避免了粒子堆积问题，但计算量较大，并且不能确保流体粒子不穿透固体边界；虚粒子法由于需要计算固体边界处的法向量，并且动态的生成边界粒子，计算量很大，对于没有确定形状的边界很难计算边界法向量，也不能保证粒子不穿透边界；预测矫正法是对每个流体粒子进行预测位置，对于即将发生穿透固体的流体粒子进行位置以及速度矫正，流体与固体之间没有相互作用力，没有物理规律支持。目前，对于SPH流固耦合的边界处理主要使用固壁虚粒子法。

1 相关工作

早期的流体模拟主要通过参数建模来模拟，文献[1]采用对波浪进行函数建模方法模拟波浪效果，但并非基于物理规律计算，效果不真实。FOSTER和METAXAS[2]首次通过对流体运动控制方程数值求解进行流体模拟。STAM和FIUME[3]首次将SPH方法引入到计算机图形学领域，模拟气体与火焰效果。随后DESBRUN和GASCUEL[4]应用SPH方法模拟可变形物体；MULLER等[5]运用SPH方法模拟流体。SPH方法在流体模拟领域得到快速发展，并进一步发展到流固耦合问题中。在SPH方法进行流固耦合模拟时，研究的难点和重点在于固体边界的处理。MONAGHAN等[6]最先提出通过固体边界上的一系列位置点对靠近固体边界的流体粒子施加一个排斥力以防止流体粒子穿透发生。后来很多学者对此进行研究，并提出了一系列的排斥力公式。HU和ADAMS[7]基于虚粒子方法实现了固体边界是平面条件下SPH流固耦合的稳定模拟；SOLENTHALER等[8]对固体边界进行粒子采样，通过提出一种统一化的粒子模型，即固壁虚粒子法，在流体粒子的物理属性值计算过程中，将固体边界粒子看作流体粒子，参与流体计算，以保证流固交界面附近流体粒子的密度和压强分布均匀，模拟出液体固化和固体融化的流固耦合效果。在固壁虚粒子方法的基础上，SUN和HAN[9]对每个流体粒子实行侦察机制，将即将发生穿透刚体的流体粒子速度改为为刚体速度，流体粒子位置改为边界点，但是局部效果不真实。AKINCI等[10]通过考虑固体边界粒子对流体粒子属性计算的相对贡献，有效解决了边界粒子分布不均匀导致的边界粒子堆积，并模拟了流体与刚性体状、面状及线状物体的耦合。以上固壁虚粒子法都不能保证流体粒子不穿透固壁边界，且计算量较大。

本文针对SPH流固耦合边界处理最常用的固壁虚粒子方法进行改进，提出SPH流固耦合模拟自适应采样固壁虚粒子边界处理方法，即为阻尼区结合自适应边界粒子采样边界处理方法。

2 SPH流固耦合建模

2.1 SPH方法

SPH方法是一种积分插值的方法，将模拟区域离散成一系列的带有物理属性值的粒子。这些物理属性值包括粒子的速度、位置、所受的力、压强和温度等。模拟区域中位置处的物理属性值()可以通过包含位置的支持域中的粒子属性值进行插值计算，即

其中，m为第粒子的质量；A为第个粒子的属性值；为第个粒子的密度；函数(-x,)为光滑核函数；x为第个粒子的位置；为光滑核函数的支持域半径。

2.2 SPH流体固体建模

因为流体的运动规律遵守Navier-Stokes方程，所以流体的运动模拟就是求解Navier-Stokes方程，即

因为本文的流体模拟中含有自由面，而自由面附近处流体粒子周围的流体粒子数目偏少，若采用式(5)进行密度计算，将导致自由面附近处流体粒子密度偏小，所以本文采用修正的密度计算公式，即

压强计算公式采用不可压缩性条件，即

因为流体粒子对之间的压力和粘性力满足牛顿第二定律，所以本文采用对称化的方式进行压力和粘性力计算，即

(1) 密度核函数

(2) 压强核函数

(3) 粘性核函数

2.3 边界粒子自适应采样

传统固壁虚粒子方法的固体边界粒子采样通常是在模拟的初始时刻就给出所有边界粒子的分布，粒子的层数通常是两层粒子，并且边界粒子数目以及分布不会随着模拟的进行而改变。在每一个模拟时间步，所有边界粒子都要进行全部物理属性值的计算，计算量较大。随着模拟的进行，流体并不总是和固体有接触，因此，没有和流体区域接触的固体边界的计算是多余的。本文提出两种边界粒子自适应采样的方法：一种是依据边界周围粒子数目的不同，边界粒子自适应地采样质量不同的大小粒子，称之为方法一；另一种是依据边界周围粒子数目的不同，边界粒子自适应的采样不同层数的相同质量粒子，称之为方法二。

2.3.1 传统的固壁虚粒子法

图1 固壁虚粒子法边界粒子采样

2.3.2 边界粒子自适应采样方法一

图2 边界粒子采样一

2.3.3 边界粒子自适应采样方法二

图3 边界粒子采样二

2.4 边界弹簧阻尼力模型

图4中，空心圆表示的是流体粒子。固体边界通过自适应固体边界粒子表示之后，固体粒子和流体粒子进行统一的连续性方程和动量方程计算，虽然解决了边界截断误差问题，但不能保证固体的非穿透性。本文在固体粒子和流体粒子进行统一的连续性方程和动量方程计算的基础上，在固壁边界附近设置一个阻尼区；当流体粒子进入事先设置的阻尼区之后，边界会给流体粒子一个弹簧阻尼力，防止粒子穿透边界。传统的边界力是关于距离的非线性力，此种类型的边界力会产生边界扰动，边界附近粒子不稳定，且一旦流体粒子穿透边界，就不再会进入流体区域。本文通过对进入阻尼区的流体粒子施加一个类似于弹簧的弹性力；弹性力是一种关于距离的线性力，所以理论上比传统的非线性力稳定。靠近边界的流体粒子速度可能很快，为了使得边界处的流体粒子更加稳定，通过施加一个和流体粒子速度相反的力，使得靠近边界具有较快速度的流体粒子速度降下来，使得边界附近的粒子更加的稳定，即

图4 弹簧阻尼模型

2.5 本文自适应流固耦合模拟算法步骤

步骤1.流体粒子初始布局。

步骤2.自适应边界粒子采样。

3 实验结果

本节通过三维溃坝模拟来验证本文提出的方法。为了更清晰、直观地演示流体粒子穿透固体边界的情形，采用三维溃坝模拟的二维剖面进行演示，并与传统的固壁虚粒子法进行比较，包括模拟时间的比较。最后给出三维溃坝模拟效果的渲染图，并与固壁虚粒子法进行比较。

3.1 二维溃坝模拟

二维溃坝实验是在windows 7 平台上进行模拟，实验中的流体粒子数是1 200，固壁虚粒子的粒子数为200，实验开发环境为visual studio 2013，实验配置为Intel (R) Core (TM) i5-2400，3.10 GHz CPU，4 G内存，NVIDIA Geforce GT 620显卡。

图5中红色粒子表示的是流体粒子，蓝色粒子表示的是固体边界粒子。从图5可以看出，固壁虚粒子边界方法中的流体粒子穿透了固体边界，而本文的两种自适应边界粒子采样方法没有粒子穿透边界。

从图6中可以看出，本文的两种方法所得到的压力更加稳定，过渡更加平滑(表1)，本文方法模拟的精度要高，模拟速度更快。

图5 二维溃坝不同方法模拟图

3.2 三维溃坝模拟

三维溃坝实验是在windows 7平台上进行模拟，实验中的流体粒子数是16 384，固壁虚粒子的粒子数为2 048，实验开发环境为visual studio 2013，以及cuda8.0进行加速计算，最后的渲染部分本实验采用了povray进行渲染，实验配置为Intel (R) Core (TM) i5-2400，3.10 GHz CPU、4 G内存，NVIDIA Geforce GT 620显卡。

从图7可以看出，本文的两种自适应边界粒子采样方法的模拟结果和固壁虚粒子边界法的模拟效果相当，其比固壁虚粒子方法模拟速度更快(表2)。

图6 二维溃坝不同方法模拟的压力云图

表1 二维溃坝不同方法模拟时间比较

图7 三维溃坝不同方法模拟渲染图

表2 三维溃坝不同方法模拟时间比较

4 结 论

本文提出一种改进的固壁虚粒子边界处理方法，通过在边界附近设置一个阻尼区，分别采用两种边界粒子自适应采样方法，不但有效地防止了粒子穿透固体边界，而且减少了计算量和内存。利用GPU加速技术，通过统一化的粒子模型并运用SPH方法实现了固体和流体的耦合模拟。

[1] PEACHEY D R. Modeling waves and surf [J]. ACM Computer Graphics, 1986, 20(4): 65-74.

[2] FOSTER N, METAXAS D. Relastic animation of liquids [J]. Graphical Models and Image Processing, 1996, 58(5): 471-483.

[3] STAM J, FIUME E. Depicting fire and other gaseous phenomena using diffusion processes [C]//Proceedings of 22ndAnnual Conference on Computer Graphics and Interactive Techiniques. New York: ACM Press, 1995: 129-136.

[4] DESBRUN M, GASCUEL M P. Smoothed particles: a new paradigm for animating highly deformable bodies [C]//Proceedings of Eurographics Workshop on Computer Animation and Simulation’ 96, Berlin: Springer, 1996: 61-76.

[5] MULLER M, CHARYPAR D, GROSS M. Particle-based fluid simulation for interactive applications [C]//Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation. Goslar: Eurographics Association, 2003: 154-159.

[6] MONAGHAN J J, THOMPSON M C, HOURIGAN K. Simulating free surface flow with SPH [J]. Journal of Computational Physics, 1994, 110(2): 399-406.

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

[8] SOLENTHALER B, SCHLAFLI J, PAJAROLA R. A unified particle model for fluid-solid interactions [J]. Computer Animation and Virtual Worlds, 2007, 18(1): 69-82.

[9] SUN H Q, HAN J Q. Particle-based realistic simulation of fluid-solid interaction [J]. Computer Animation and Virture World, 2010: 21(6): 589-595.

[10] AKINCI N, IHMSEN M, AKINCI G, et al. Versatile rigid-fluid coupling for incompressible SPH [J]. ACM Transactions on Graphics, 2012, 31(4): 62: 1-62: 8.

Adaptive Sampling Method for Solid Wall Virtual Particle Boundary Handling in SPH Fluid Solid Coupling Simulation

ZHU Xiaolin, WANG Huanhuan, YIN Jingcun

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

The method for solid wall virtual particle boundary handling is a main boundary handling method in fluid simulation, while it cannot ensure that no fluid penetrates the solid boundary, and requires a huge amount of computation. To prevent fluid particles from penetrating the boundary, in this paper a damp zone is set up near the boundary, and each fluid particle in this zone will be applied by the boundary with an elastic force and a damping force that is opposite to the moving velocity of the fluid particle. Therefore, the particles near the boundary become more stabilized. To reduce the amount of computation, this paper proposes two adaptive sampling methods for boundary particles based on different numbers of particles around the boundary: one is to adaptively sample particles of different masses; the other is to adaptively sample particles of the same mass in different layers. Compared with the traditional method for solid boundary particle sampling, these two adaptive sampling methods reduce the number of boundary particles, accelerate the simulation speed and save the computer memory. In addition, all of the 3D fluid simulations in this paper are based on the GPU acceleration technique, and can interact in real time.

solid wall virtual particle; damping zone; elastic force; damping force; adaptive sampling of boundary particles; GPU acceleration technique

TP 391

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

A

2095-302X(2018)03-0381-08

2017-01-14；

2017-05-11

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

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