朱曉臨，殷競存，汪歡歡

?

改進的流體模擬固體邊界處理算法

朱曉臨，殷競存，汪歡歡

(合肥工業(yè)大學數(shù)學學院，安徽 合肥 230009)

流體模擬是計算機圖形學的一個重要研究分支，流體的固體邊界處理一直是流體模擬的研究重點，光滑粒子流體動力學(SPH)方法中的鏡像粒子法是處理固體邊界的一個重要方法。鏡像粒子法通過靠近邊界的流體粒子在邊界外動態(tài)生成對應的鏡像粒子來處理固體邊界問題，但隨著邊界復雜程度的提高，傳統(tǒng)的鏡像粒子法生成鏡像粒子的復雜度也隨之提高，模擬效率隨之降低。為此，文章對鏡像粒子法進行改進，提出一種新的鏡像粒子場量求值方法，有效地降低了復雜邊界情況下生成鏡像粒子的復雜度，且使靠近邊界的流體粒子場量更加均勻。仿真實驗結(jié)果表明，隨著流體模擬粒子數(shù)的增加以及邊界復雜程度的提高，該方法比傳統(tǒng)鏡像粒子法效率高的優(yōu)勢也更加明顯。

光滑粒子流體動力學方法；流體模擬；固體邊界處理；鏡像粒子法；復雜場景

基于物理的流體模擬方法一直是計算機動畫領(lǐng)域的一個研究熱點，在游戲、廣告和影視特效中得到了廣泛應用。流體的固體邊界處理更是計算機動畫中流體模擬技術(shù)的重要研究內(nèi)容。

光滑粒子流體動力學(smoothed particle hydrodynaimics，SPH)方法的流體模擬中存在3類流固耦合方法：懲罰力方法、鏡像粒子法和靜態(tài)粒子法。

懲罰力方法是通過構(gòu)造一種基于距離和速度的懲罰力函數(shù)，對靠近邊界的流體粒子施加反作用力以防止流體穿透固體邊界。但是，當懲罰力較大時，會導致邊界處流體粒子的壓強出現(xiàn)明顯變化，懲罰力較小時，流體粒子就會非物理性穿透邊界。此外，由于邊界截斷，流固交界面上流體粒子數(shù)量會相對不足，流體粒子就會附著在固體邊界上。鏡像粒子法，即靠近邊界的流體粒子在邊界外動態(tài)生成對應的鏡像粒子，一定程度上能夠提高流體邊界粒子精度以及流固交界面處流體的壓強分布均勻，但存在著對于復雜和可變形固體邊界難以生成鏡像粒子的問題。靜態(tài)粒子法在邊界外分布幾層固定的邊界粒子，并將其看作位置固定的特殊流體粒子參與流體的物理計算過程。該方法實施簡單，但是在運動劇烈時流體粒子容易穿透。

MONAGHAN[1]使用L-J形式的基于粒子距離計算粒子間的相互作用力，但是忽略了邊界處流體粒子運動方向與邊界平行，易造成流體紊動現(xiàn)象。MONAGHAN和PRICE[2]提出新的構(gòu)造函數(shù)，解決了上述問題。MüLLER等[3]基于懲罰力方法實現(xiàn)了微可壓縮SPH流體與可變形固體之間的交互模擬。使用較大的懲罰力能夠保證流體不會穿透固體邊界，但同時必須使用較小的時間步以保證模擬的穩(wěn)定性。BECKER等[4]提出了直接力方法，通過預測-校正技術(shù)，對粒子速度和位置進行分布計算。上述方法都存在一個問題，由于流固交界面上流體粒子數(shù)量不足而造成流體粒子附著在固體邊界上。HARADA等[5]針對這類問題提出了一種邊界權(quán)重函數(shù)，通過計算流體粒子到邊界的垂直距離，可以在流體密度計算中加入固體邊界的貢獻以彌補邊界粒子的不足。該方法解決了流體粒子的附著問題，但只適合規(guī)則固體邊界的情況，且會導致固體邊界處流體密度分布不均勻。

HU和ADAMS[6]使用鏡像粒子法模擬了規(guī)則固體邊界條件下流體運動。SOLENTHALER和GROSS[7]重新定義接近自由表面的粒子采樣方法。SCHECHTER和BRIDSON[8]提出一種ghost方法用于模擬流體的自由表面和邊界條件。BECKER和TESCHNER[9]提出泰特方程去降低不可壓縮性，并且介紹利用高阻尼方程初始化粒子使鏡像粒子達到穩(wěn)定的密度。但是，在復雜或可變形固體邊界處難以生成鏡像粒子是以上鏡像粒子法均存在的問題。

SOLENTHALER等[10]基于靜態(tài)粒子法的思想，對固體邊界進行粒子采樣，并將這些粒子看作位置固定的特殊流體粒子，參與流體的物理計算過程，以保證流固交界面附近流體粒子壓強和速度分布均勻。但是，在流體粒子速度較大時會出現(xiàn)粒子穿透固體邊界的現(xiàn)象，因此需要使用較小的時間步長來解決此問題。DALRYMPLE和KNIO[11]通過對邊界進行多層靜態(tài)粒子采樣，以避免流體粒子穿透固體邊界。IHMSEN等[12]在靜態(tài)粒子法的基礎上，對穿透固體邊界的粒子實施位置校正，同時在物理計算過程中加入了固體邊界對流體粒子的貢獻，有效避免了流固交界面附近流體粒子密度變化過大的問題。邵緒強等[13]提出一種靜態(tài)粒子與位置-速度校正技術(shù)相結(jié)合的方法，解決了流體與固體邊界交互時流體粒子的穿透與衰減問題；但靜態(tài)粒子作為特殊的流體粒子，仍然存在邊界截斷所造成的誤差。

在流固耦合模擬中，鏡像粒子法由于其本身特性越來越受到廣泛應用，但在復雜邊界下卻存在計算量過大、算法復雜等缺點。針對鏡像粒子法的上述缺點，本文提出一種改進方法：改變鏡像粒子分布、采用新的鏡像粒子場量計算方法。大量實驗結(jié)果表明，隨著模擬粒子數(shù)和邊界復雜度的提高，本文方法比傳統(tǒng)鏡像粒子法效率高的優(yōu)勢更加明顯。

1 流固耦合算法

1.1 SPH流體

SPH方法是一種無網(wǎng)格拉格朗日型的粒子方法，在1977年分別由LUCY[14]與GINGOLD和MONAGHAN[15]提出，從開始用于解決三維開放空間的天體物理學問題，到目前被廣泛應用于流體動力學、固體力學及其他工程學科各種問題的數(shù)值仿真。

DESBRUN和GASCUEL[16]首次將SPH方法引進圖形學領(lǐng)域，模擬可變形固體；MüLLER等[17]則采用SPH方法來模擬水，為該方法在模擬水等流體方面奠定了算法基礎。

SPH是一種無網(wǎng)格插值方法，通過對離散粒子的物理屬性進行插值計算，而得到連續(xù)區(qū)域內(nèi)任一點處的對應值

其中，()是位置處的物理屬性值；為鄰居粒子總數(shù)；m、、A、分別為位置處支持域范圍內(nèi)粒子的質(zhì)量、密度、物理屬性值和位置；(–,)為光滑核函數(shù)；為光滑核半徑。

將式(1)應用于流體密度可得

在每個時間步內(nèi)，流體粒子的速度、位置等信息根據(jù)其所受合力進行計算。合力包括壓力、粘性力和外力，外力為重力，即

通過式(2)求得粒子密度后，需要對粒子壓強進行計算，以求得粒子壓力。

文獻[14]使用氣體狀態(tài)方程式來計算粒子壓強，即

其中，為氣體常數(shù)；0為常量密度。

1.2 鏡像粒子法改進

傳統(tǒng)鏡像粒子法中，在每一時間步，靠近邊界的流體粒子在邊界外鏡像生成虛擬粒子(圖1)，且將靠近邊界的流體粒子場量賦給相應的鏡像粒子場量，即。

其中，j為靠近邊界的流體粒子；i為對應的鏡像粒子；A(i)ghost為鏡像粒子i的場量；A(j)fluid為流體粒子j的場量。

因此，傳統(tǒng)鏡像粒子法在每一時間步需動態(tài)生成鏡像粒子，并重新確定鏡像粒子的位置，計算量有所增加。而且在復雜邊界情況下，難以確定鏡像粒子的位置，使得計算量和算法復雜度大大提高。

與之相反，本文提出一種固定的鏡像粒子方法。首先，將提出的鏡像粒子固定分布在固體邊界表面的外側(cè)，保持位置不變，分布層數(shù)與粒子影響域和粒子初始分布有關(guān)其次，為計算每個固定鏡像粒子的場量，本文引入插值點這一概念，使每個固定鏡像粒子都有其對應的插值點，通過計算插值點處的場量而得出固定鏡像粒子的場量。每個固定鏡像粒子對應插值點的位置與固定鏡像粒子的位置關(guān)于固體邊界相互對稱(圖2)。

圖2 本文方法鏡像粒子分布

每個固定插值點的變量由傳統(tǒng)的SPH方法通過相鄰的流體粒子獲得，見式(10)。然后將插值點的場量賦值給相應的ghost粒子，見式(7)。

對于不同要求的滑移條件，對速度可采用式(8)或式(9)。

對于非滑移邊界，有

對于滑移邊界，有

利用傳統(tǒng)SPH方法對插值點的場量進行求解，即

其中，為插值粒子；為領(lǐng)域內(nèi)的流體粒子。

1.3 算法描述

本文算法如圖3所示。

圖3 本文算法流程圖

仿真算法步驟如下：

步驟1. 初始化流體粒子和鏡像粒子。

步驟3. 通過式(2)計算插值點密度inter，式(5)計算壓強inter。

步驟4. 通過式(10)計算插值點速度inter。

步驟5. 通過式(7)~(9)得到鏡像粒子密度ghost、壓強ghost和速度ghost。

步驟6. 通過式(2)計算流體粒子密度，式(5)計算流體粒子壓強P。

步驟7.對流體粒子施加力，即=pressure+viscosity+other。

步驟10. 重復步驟3。

2 實驗

2.1 二維潰壩

潰壩是經(jīng)典的自由表面流體運動模型，常用來驗證壁面邊界施加方法的有效性。

計算參數(shù)有：初始密度為1 000 kg/m3，邊界粒子數(shù)目為292，流體粒子數(shù)目1 681，流體粒子初始間距為0.018 m，時間步長為1.0′10–3s，粒子的影響域=0.04 m。

2.1.1 本文方法與不施加邊界處理相比較

由圖4可知，圖4(a)流體粒子不論是在固體邊界還是在流體表面，粒子分布都比圖4(b)更加均勻。因此，本文方法能夠使粒子分布更加均勻，使場量更加穩(wěn)定。

圖4 本文方法與不施加邊界力法比較

2.1.2 本文方法與靜態(tài)粒子法相比較

由圖5可知，圖5(a)流體粒子在邊界處沒有穿透現(xiàn)象；圖5(b)流體粒子在邊界處有穿透現(xiàn)象。因此，本文方法比靜態(tài)粒子法更有效的防止粒子滲透。

圖5 本文方法與靜態(tài)粒子法比較

2.2 三維潰壩

實驗在Window 7平臺上進行，采用visual studio 2013，cuda8.0和POV。實驗環(huán)境為Intel(R) Core(TM) i5-2400 3.10 GHz CPU，4 G內(nèi)存，NVIDIA Geforce GT 620顯卡。

由圖6~9可以看出，本文方法能夠有效模擬各種復雜邊界下的流體運動；由表1~4可知，隨著邊界復雜程度的提高和模擬粒子數(shù)的增加，本文方法的耗時比鏡像粒子法耗時短的效果更加明顯。

圖6 本文方法與鏡像粒子法比較1

圖7 本文方法與鏡像粒子法比較2

圖8 本文方法與鏡像粒子法耗時比較3

圖9 本文方法與鏡像粒子法耗時比較4

表1 本文方法與鏡像粒子法耗時比較1

表2 本文方法與鏡像粒子法耗時比較2

表3 本文方法與鏡像粒子法耗時比較3

表4 本文方法與鏡像粒子法耗時比較4

3 結(jié) 論

本文基于SPH方法提出一種新的流固邊界處理算法，給出一種虛擬粒子場量新的求值方法，大大降低了處理邊界的復雜程度；隨著模擬粒子數(shù)和邊界復雜度的提高，本文方法比傳統(tǒng)鏡像粒子法效率高的優(yōu)勢更加明顯。

本文方法渲染效果未能達到滿意程度，是受當前硬件限制和軟件條件，后期將進行相關(guān)改進。此外，在容器較大、流體運動比較平緩的情況下，預先固定鏡像粒子將造成一定的浪費。因此，接下來將對這兩方面進行改進。

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

[2] MONAGHAN J J, PRICE D J. Toy stars in one dimension [J]. Monthly Notices of the Royal Astronomical Society, 2004, 350(4): 1449-1456.

[3] MüLLER M, SCHIRM S, TESCHNER M. Interaction of fluids with deformable solids [J]. Journal of Visualization and Computer Animation, 2004, 15(3-4): 159-171.

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

[5] HARADA T, KOSHIZUKA S, KAWAGUCHI Y. Real-time fluid simulation coupled with cloth [C]//Theory and Practice of Computer Graphics. Coslar: Eurographics Press, 2007: 13-18.

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

[7] SOLENTHALER B, GROSS M. Two-scale particle simulation [J]. ACM Transactions on Graphics (TOG), 2011: 30(4): 1-8.

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

[9] BECKER M, TESCHNER M. Weakly compressible SPH for free surface flows [C]//ACM SIGGRAPH/Eurographics Symposium on Compute Animation. New York: ACM Press, 2007: 209-217.

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

[11] DALRYMPLE R A, KNIO O. SPH modelling of water waves [C]//Proceedings of the Fourth Conference on Coastal Dynamics. New York: ACM Press, 2000: 779-787.

[12] IHMSEN M, AKINCI N, GISSLER M, et al. Boundary handling and adaptive time-stepping for PCISPH [C]//Proceedings of the Seventh Workshop on Virtual Reality Interactionsand Physical Simulations. Switzerland: Eurographics Association Press, 2011: 79-88.

[13] 邵緒強, 周忠, 張勁松. 微可壓縮SPH流體的穩(wěn)定性固體邊界處理算法[J]. 計算機輔助設計與圖形學學報, 2014, 26(11): 1915-1922.

[14] LUCY L B. A numerical approach to the testing of the fission hypothesis [J]. The Astronomical Journal, 1977, 82(82): 1013-1024.

[15] GINGOLD R A, MONAGHAN J J. Smoothed particle hydrodynamics: theory and application to non-spherical stars [J]. Monthly Notices of the Royal Astronomical Society, 1977, 181(3): 375-389.

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

[17] MüLLER M, CHARYPAR D, GROSS M. Particle-based fluid simulation for interactive applications [C]// Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Compute Animation. Switzerland: Eurographics Association Press, 2003: 154-159.

An Improved Solid Boundary Treatment Algorithm of Fluid Simulation

ZHU Xiaolin, YIN Jingcun, WANG Huanhuan

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

Fluid simulation is an important research branch of computer graphics, solid boundary treatment of fluid has always been the focus of fluid simulation, image particle method in smooth particle hydrodynamics is an important method to deal with solid boundary. Image particle method deals with solid boundary problems by dynamically generating corresponding mirror particles outside the boundary by fluid particles near the boundary. However, with the improvement of the complexity of the boundary, the complexity of the traditional mirror particle method to generate the mirror particles is also increased, and the simulation efficiency is reduced. In order to solve the above problem, the mirror particle method is improved in this paper by giving a new technique for calculating values of mirror particle field, which effectively reduces the complexity of the generation of traditional mirror particles with complex boundary conditions, and the field of fluid particles near the boundary is more uniform. A lot of simulation results show that this method becomes more efficient compared with the mirror particle method with the increasement of the number of particles and the complexity of the scene.

smoothed particle hydrodynaimics method; fluid simulation; solid boundary treatment; ghost particle field; complex scene

TP 391.9

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

A

2095-302X(2018)02-0263-06

2017-06-22；

2017-11-08

國家自然科學基金項目(61272024)

朱曉臨(1964–)，男，安徽池州人，教授，博士。主要研究方向為數(shù)值逼近、CAGD、圖形圖像處理。E-mail：zxl_hfut@126.com