朱曉臨，殷競(jìng)存，汪歡歡

?

改進(jìn)的流體模擬固體邊界處理算法

朱曉臨，殷競(jìng)存，汪歡歡

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

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

光滑粒子流體動(dòng)力學(xué)方法；流體模擬；固體邊界處理；鏡像粒子法；復(fù)雜場(chǎng)景

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

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

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

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

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

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

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

1 流固耦合算法

1.1 SPH流體

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

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

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

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

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

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

通過式(2)求得粒子密度后，需要對(duì)粒子壓強(qiáng)進(jìn)行計(jì)算，以求得粒子壓力。

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

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

1.2 鏡像粒子法改進(jìn)

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

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

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

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

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

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

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

對(duì)于非滑移邊界，有

對(duì)于滑移邊界，有

利用傳統(tǒng)SPH方法對(duì)插值點(diǎn)的場(chǎng)量進(jìn)行求解，即

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

1.3 算法描述

本文算法如圖3所示。

圖3 本文算法流程圖

仿真算法步驟如下：

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

步驟3. 通過式(2)計(jì)算插值點(diǎn)密度inter，式(5)計(jì)算壓強(qiáng)inter。

步驟4. 通過式(10)計(jì)算插值點(diǎn)速度inter。

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

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

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

步驟10. 重復(fù)步驟3。

2 實(shí)驗(yàn)

2.1 二維潰壩

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

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

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

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

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

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

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

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

2.2 三維潰壩

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

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

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

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

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

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

表1 本文方法與鏡像粒子法耗時(shí)比較1

表2 本文方法與鏡像粒子法耗時(shí)比較2

表3 本文方法與鏡像粒子法耗時(shí)比較3

表4 本文方法與鏡像粒子法耗時(shí)比較4

3 結(jié) 論

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

本文方法渲染效果未能達(dá)到滿意程度，是受當(dāng)前硬件限制和軟件條件，后期將進(jìn)行相關(guān)改進(jìn)。此外，在容器較大、流體運(yùn)動(dòng)比較平緩的情況下，預(yù)先固定鏡像粒子將造成一定的浪費(fèi)。因此，接下來將對(duì)這兩方面進(jìn)行改進(jì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] 邵緒強(qiáng), 周忠, 張勁松. 微可壓縮SPH流體的穩(wěn)定性固體邊界處理算法[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 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

國(guó)家自然科學(xué)基金項(xiàng)目(61272024)

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