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中，空心圓表示的是流體粒子。固體邊界通過自適應固體邊界粒子表示之后，固體粒子和流體粒子進行統一的連續性方程和動量方程計算，雖然解決了邊界截斷誤差問題，但不能保證固體的非穿透性。本文在固體粒子和流體粒子進行統一的連續性方程和動量方程計算的基礎上，在固壁邊界附近設置一個阻尼區；當流體粒子進入事先設置的阻尼區之后，邊界會給流體粒子一個彈簧阻尼力，防止粒子穿透邊界。傳統的邊界力是關于距離的非線性力，此種類型的邊界力會產生邊界擾動，邊界附近粒子不穩定，且一旦流體粒子穿透邊界，就不再會進入流體區域。本文通過對進入阻尼區的流體粒子施加一個類似于彈簧的彈性力；彈性力是一種關于距離的線性力，所以理論上比傳統的非線性力穩定?？拷吔绲牧黧w粒子速度可能很快，為了使得邊界處的流體粒子更加穩定，通過施加一個和流體粒子速度相反的力，使得靠近邊界具有較快速度的流體粒子速度降下來，使得邊界附近的粒子更加的穩定，即

圖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