基于位置的流體實時交互仿真①

王 坤,于 歌,梁 驥,郭麗麗

1(中國科學院 空間應用工程與技術中心 太空應用重點實驗室,北京 100094)

2(中國科學院大學,北京 100049)

引言

剛體與流體的交互是我們生活中常見到的,比如攪拌一杯水或者手直接與水接觸.在計算機仿真領域,與流體的力覺交互中,我們希望得到的是一種實時的逼真的交互體驗,這種交互手段在科學實驗仿真、醫療訓練等領域是迫切需要的.這對科學家提供更有效率、更接近真實結果的實驗或者訓練過程是非常有幫助的.

目前,模擬剛體與流體的實時交互在視覺與觸力覺的真實呈現仍然是一個頗具挑戰的問題.首先對于流體模擬,在傳統的方法中,通常分為兩類:基于網格模型的歐拉法[1]和基于粒子模型的拉格朗日法,目前采用最多的是基于粒子模型的拉格朗日法.由于常用的.其中,光滑粒子流體動力學(Smoothed Particle Hydrodynamic,SPH)是拉格朗日模型中最流行的粒子方法.Muller[2]表明SPH可用于模擬交互式流體的粘滯力和表面張力,他于2004年[3]將引力、斥力和粘滯力引入流體粒子與變形體網格之間的交互建模中實現了剛體-流體交互的實時仿真.

由于SPH算法需要較小的時間步長,對計算開銷較大,另外算法中采用的核函數的光滑半徑較大,實時性以及計算精度難以達到較高的要求.Muller[4]于2007年提出了一種基于位置的動力學方法,采用了較小的光滑半徑以及更大的時間步長,通過對粒子位置加以約束,直接更新粒子位置來模擬流體粒子,避免了需要先計算粒子受力然后才能計算粒子位置改變從而帶來的不穩定性,極大的提升了流體模擬的效率以及真實度.

一般來說,為了實現剛體與流體在混合環境下的交互,我們需要為每種類型的對象使用不同的算法,并在這些異構模型之間建立耦合接口,這將增加已經非常復雜和耗時的模擬的復雜性和計算成本.Keiser等[5]將這一算法過程進行改進,實現了固體、可變形和流體之間的相互作用.但依然難以達到實時性的要求.隨著硬件性能的不斷提升帶來的自然加速,更多的研究者將針對新硬件的體系結構來設計整個算法,以充分利用硬件新特性帶來的優勢;另外結合算法本身的不斷創新使整個流體模擬向實時性不斷邁進.而基于位置的流體算法的粒子模型也非常適合將算法實現并行化.

傳統的基于CPU的流體模擬,難以達到實時高效的流體仿真.GPU計算的核心思想是將可以被并行的指令盡可能多的移植到GPU上執行,充分利用GPU強大的并行計算能力批量執行指令.Hegeman[6]驗證了在GPU的并行架構上實現粒子模型的有效性表明GPU的運行速度對比純CPU算法的運行速度有超過一個數量級的提升.2013年,Domínguez[7]等人優化了現存基于CPU-GPU計算平臺的SPH計算模型,主要解決了CPU與GPU仿真數據的傳輸問題.

本文在各種模擬流體的算法中,綜合考慮了各種算法在交互過程中的利弊,采用一種SPH的改進算法,結合Müller的基于位置動力學(Position Based Dynamics)[4]得到的基于位置的流體算法(Position Based Fluid)[8]來實現流體的基本模擬,其次通過剛性體與流體粒子的碰撞檢測與響應計算交互力.本文著重探討了剛性體與流體的交互力覺渲染,基于CUDA平臺實現算法的并行化,提供高效穩定的計算結果,在視覺與力覺上實現真實與實時的流體交互.

1 光滑粒子動力學基礎

本文采用基于郎格朗日粒子法來模擬流體的運動,其中最常用的是光滑粒子流體動力學,該方法是一種用于流體動力學的插值方法,可以通過插值計算出流體的物理特性.用粒子代替流體來模擬流體運動,它把流體定義為在空間中離散分布并且位置可計算的粒子,這些粒子具有流體的質量、密度、速度等物理屬性.SPH方法通過積分近似和粒子近似將流體動力學方程進行離散,每個粒子上的屬性函數fi(p)通過支持域內相鄰粒子的疊加求和計算得到:

其中,j表示粒子i支持域內的第j個粒子,mj表示粒子j的質量,pj表示粒子j的位置,ρj表示粒子j的密度,fj表示在位置p處的某種屬性函數,如密度、壓力等,W表示核函數,h表示光滑核半徑,如圖1所示.

圖1 SPH粒子算法

2 流體模擬方法

2.1 基于位置的流體

在本文中,對于流體模擬,本文采用的是Macklin[8]提供的基于位置的流體模擬方法,在SPH方法的基礎上進行改進,采用了更大的時間步長以及較小的光滑半徑,可以直接更新粒子位置變化來模擬流體,避免了SPH流體模擬中的不穩定性.

根據Macklin[8]的方法,為了實現密度計算的穩定,在運動中不會發生較大的突變,對每一個粒子位置施加一個非線性約束系統,約束的狀態方程為:

其中,pn為i粒子的第n個鄰域粒子的位置,ρ0為算法中給定的初始密度,ρi為標準SPH算法的密度計算,根據公式(1)求得密度計算為:

其中所有粒子質量mj相同,核函數W采用Müller[3]給出的Poly6光滑核函數來求解密度值.

在該Müller[4]的基于位置動力學算法中,對粒子位置給定一個修正,所有粒子的位置滿足約束:

給定位置p的變化Δp在約束函數梯度的大小關系,其系數為λ.

結合公式(5)通過對約束的梯度計算得到:

結合Monaghan[9]給出的SPH方法,代入公式(2)和(3),推導出粒子i的位置約束梯度函數為:

結合公式(6)與(7)得到λ的解為:

因為公式(2)中的約束函數與密度計算有關,而密度計算是非線性的,當流體粒子太分散時容易導致密度計算結果的不穩定,根據Smith[10]提出方法中,可以通過混合一些約束來提高穩定性,公式(6)修改為:

其中ε為我們在模擬過程中定義的常量約束參數.λ的求解修改為:

在每一次循環計算過程中,都需要計算λ的值,然后結合公式(5),按照如下公式計算粒子i的位置更新,其計算為

如上所述,在基于位置的流體算法中,每次循環執行流體粒子位置的動態更新.當計算得到粒子的位置更新后,便可根據這一時間段內進一步計算得出粒子的速度以及加速度.每次循環中獨立解決每個位置約束,同時,在每次循環過程中,對剛性體進行碰撞檢測,每個步驟中我們重新計算每一個粒子的領域粒子,并且重新計算每個粒子位置的約束值.這樣可以避免對粒子密度的低估,導致最終結果出現較大偏差.

2.2 剛體與流體的交互力計算

根據Tse[11]給出的方法中,剛體與流體的交互過程中,工具球會受到流體的吸附力,當工具球與流體粒子的距離在支持域內,吸附力開始產生作用,流體粒子j與工具球的球心距為rj=p-pj,吸附力計算如下:

Ki,Kj為吸附力計算參數,WAdh為吸附力計算的光滑核函數,其計算為:

其中ε>ζ,Q=(B-A)..

A:吸附力光滑函數起始點.

B:吸附力光滑函數結束點.

ε,ζ:可調節吸附力,本文中分別為ε=3,ζ=2.

工具與流體在交互過程中除受到吸附力以外,還受到因碰撞產生的阻尼力,其大小計算方式如下:

KDamp為阻尼系數,0>KDamp>1,v0為工具球的速度,vj為與工具球碰撞的粒子j的速度,可以看出,碰撞力的大小除了與碰撞系數有關外,也和工具球與流體粒子的速度差有關,這也符合我們日常生活中所見到的情景.

3 流體交互算法實現

3.1 流體粒子初始網格構建

在流體模擬階段,計算每一個粒子的物理屬性需要改粒子鄰域內的粒子信息,因此,鄰域粒子搜索算法是整個算法中比較耗時的步驟,對流體粒子所在空間進行均勻網格劃分是目前被認為最為有效的方法.基于CUDA平臺我們可以快速的實現對流體粒子的網格構建,如圖2.

圖2 流體粒子網格構建

采用Nvidia[12]的粒子模擬方法來實現對流體粒子的初始網格劃分并記錄每一個流體粒子的領域粒子的索引值,記錄如表1,便于執行后續的鄰域粒子快速查找,同時能夠快速執行與交互工具的碰撞檢測和碰撞響應.在每一次循環計算中,我們都需要執行一個完整的CUDA程序.

表1 鄰域粒子的索引值記錄

3.2 算法流程

整個算法流程圖分為如圖3幾個階段.由于流體模擬算法的粒子性,我們可以并行的對每一個粒子的狀態并行的求解,在每一個階段通過CUDA平臺并行的對每一個粒子的狀態進行求解,包括對初始粒子網格的建立到鄰域粒子的搜索以及對于每一個流體粒子各項物理屬性的計算,在受力計算階段實時的通過設備輸出交互力,在最后階段計算粒子的位置、速度等狀態并更新繪制顯示.同時將這些狀態用于下一個循環中的位置預測計算.

圖3 算法流程

在每一個階段更為詳細的算法過程如下:

算法1.流體交互算法#初始化#1.設置粒子初始位置、速度等狀態量.2.根據粒子初始位置劃定網格.#位置預測#3.由外力以及設備反饋力更新速度.4.預測粒子位置.#鄰域粒子查找#5.根據網格數據查找粒子i的鄰域粒子并記錄索引.#迭代過程#6.在迭代次數內對粒子i執行.7.根據公式(3)計算粒子i的密度.8.根據公式(10)計算λi.9.根據公式(11)計算粒子i的位移Δpi.10.碰撞檢測.11.根據Δpi更新粒子位置.12.重新組織粒子和網格.#計算受力#13.執行碰撞響應.14.利用公式(12)計算附著力FAdh.15.利用公式(14)計算阻尼力FDamp.16.將工具所受合力反饋到力覺設備端.#粒子狀態更新#17.更新粒子i位置.18.更新粒子i速度.19.根據新的粒子位置渲染流體.

4 實驗方案

4.1 實驗設備及環境

開發環境如表2所示,算法通過CUDA并行化,運行在GPU端,采用的顯卡為NVIDIA Quadro M5000,包含2048個CUDA并行處理核心,所有并行計算的步驟都在GPU端執行,核心數越多,并行加速效率越高.

交互設備采用的是Geomagic公司開發的三自由度力覺交互設備Geomagic Touch,如圖4.該反饋設備具有六自由度操作位置檢測,三自由度力覺交互,通過RJ45以太網接口或USB接口連接終端.可以通過該設備實時控制工具球在流體中的運動并且反饋出工具球受到的合力,模擬出真實場景中在流體中交互的受力情況.

表2 實驗設備以及開發環境

圖4 Geomagic Touch三自由度力覺交互設備

4.2 實驗過程

4.2.1 剛體與流體的交互力模擬驗證

為了驗證基于位置的流體算法模擬流體與剛體的交互,我們建立了關于剛性工具在與流體的交互過程中的受力模擬實驗,實驗分為兩個部分,第一個部分模擬工具球在交互過程中的受力情況并進行分析,驗證交互過程的真實性;第二個部分通調整粒子數目來驗證力覺交互設備感知到的交互力是否穩定連續來驗證交互過程的穩定與實時性.

第一部分的實驗我們模擬的是工具球在流體交互過程中的受力情況,我們設置的初始粒子數目為20 k,從視覺角度來說,已經可以很好的滿足這一部分實驗的精細度.粒子的初始位置按照容器內部空間均勻排列,初始速度為零,初始受力為容器壁的約束受力以及重力.一個完整的模擬實驗過程包括從工具球向下方向勻速運動,從剛接觸流體表面到開始進入流體,然后完全沒入流體,在流體中勻速運動直到碰到杯壁返回到中間位置,最后再勻速向上離開流體,到工具球最后剛到達流體上表面因吸附力沾染流體,到最后完全離開流體.整個實驗過程的部分場景如圖5所示.

4.2.2 交互力模擬實驗結果

首先我們驗證的是工具球在流體交互中的受力情況,圖6模擬的是通過力覺交互設備控制工具球在流體中運動時的受力情況,如圖所示,X軸、Y軸和Z軸正方向的力分別用實線、點線和虛線來表示,其方向已經在圖6中標出.在整個過程中我們從圖中我們可以看到工具球受力是比較平滑的曲線,沒有發生大范圍的突變,說明在我們的算法中實現的力是穩定連續的.

圖5 工具在流體中的交互運動場景

圖6 工具球在流體中運動時的受力情況分析

實驗開始時,當工具球剛開始接觸流體表面時,接觸部分受到流體的吸附力,隨著工具球的逐漸進入,工具球還受到浮力的影響逐漸增大,到完全沒入流體,工具球受到的浮力趨于穩定.然后工具球在流體中沿-X水平方向做勻速運動,收到一個恒定的阻力.然后運動到邊界與杯壁發生碰撞,產生碰撞力,接著沿X方向運動到中間位置,期間受到-X方向的恒定阻力,最后工具球開始向上運動到開始離開流體表面,豎直方向上受到的流體浮力開始減少,到即將離開的一瞬間,工具球在此時沒有浸入流體已經不受浮力影響了,此時工具球還受到一個向下的吸附力,隨著工具球的進一步向上運動,最終工具球不受到流體交互力的影響.在整個運動過程中沒有Z方向的運動,因此Z方向的受力幾乎穩定為零,這也符合實驗預期結果.

4.2.3 交互力的連續性驗證及結果

對于流體模擬,需要一定數目的粒子,才能達到較好的模擬真實性,但是粒子數增加的同時會帶來計算開銷的增加,影響交互的實時性.在整個交互過程中,當力覺交互設備的輸出頻率在1 kHz左右時,我們通過設備能感知到連續的力覺,而不會有頓挫感.第二部分的模擬實驗中,在粒子數目能夠達到視覺仿真要求的前提下,通過改變粒子數目來驗證算法的效率,保證在每一次力覺交互循環中的時間步長在1 ms左右即可保證力輸出的連續性.

如下圖7所示,我們測試了不同粒子數目時,力覺交互設備的輸出頻率,粒子數目從10 k逐漸增加到80 k,每次都測試出整個實驗過程中力覺交互設備的刷新率,粒子數目維持在10 k的時候基本上可以滿足我們交互過程中的精細度,繼續增加粒子數目可以帶來更好的流體模擬效果,當然粒子數目越多可以模擬更大規模的流體交互,模擬的流體場景也更為復雜.當粒子數目達到60 k的時候已經遠遠達到流體模擬的精度要求,雖然頻率稍有下降,但仍然能夠很好的滿足交互的實時性,說明交互過程的實時性以及真實度達到了要求.

圖7 粒子數目與力覺交互輸出頻率的影響

5 總結

本文探討了剛體與流體交互過程中的力覺交互,實現了基于位置的流體模擬算法并且模擬了剛體與流體的交互過程,實時計算并渲染出剛體工具所受到的交互力,有效的避免了傳統SPH算法的精度較低、計算實時性較差等劣勢,同時剛性工具與流體的交互過程中很好的模擬了流體交互過程中的包括吸附力、浮力等影響比較重要的力覺.采用CUDA并行計算,極大的提高了流體算法的模擬效率以及精度.通過模擬實驗有效的驗證了交互的實時性以及真實度.對后續流體交互研究具有非常好的指導意義,我們未來將進一步挖掘流體交互中的諸如湍流等物理現象,提高交互的真實度.

1Foster N,Metaxas D.Realistic animation of liquids.Graphical Models and Image Processing,1996,58(5):471-483.[doi:10.1006/gmip.1996.0039]

2Müller M,Charypar D,Gross M.Particle-based fluid simulation for interactive applications.Proceedings of the 2003 ACM SIGGRAPH/Eurographics Symposium on Computer Animation.San Diego,CA,USA.2003.154-159.

3Müller M,Schirm S,Teschner M,et al.Interaction of fluids with deformable solids:Research articles.Computer Animation and Virtual Worlds,2004,15(3/4):159-171.

4Müller M,Heidelberger B,Hennix M,et al.Position based dynamics.Journal of Visual Communication and Image Representation,2007,18(2):109-118.[doi:10.1016/j.jvcir.2007.01.005]

5Keiser R,Adams B,Gasser D,et al.A unified lagrangian approach to solid-fluid animation.Proceedings of the Second Eurographics/IEEE VGTC Conference on Point-Based Graphics.New York,NY,USA.2005.125-133.

6Hegeman K,Carr NA,Miller GSP.Particle-based fluid simulation on the GPU.Proceedings of the 6th International Conference on Computational Science.Reading,UK.2006.228-235.

7Domínguez JM,Crespo AJC,Gómez-Gesteira M.Optimization strategies for CPU and GPU implementations of a smoothed particle hydrodynamics method.Computer Physics Communications,2013,184(3):617-627.[doi:10.1016/j.cpc.2012.10.015]

8Macklin M,Müller M.Position based fluids.ACM Transactions on Graphics (TOG),2013,32(4):104.

9Monaghan JJ.Smoothed particle hydrodynamics.Annual Review of Astronomy and Astrophysics,1992,30(1):543-574.[doi:10.1146/annurev.aa.30.090192.002551]

10Smith R.Open dynamics engine v0.5 user guide.Computer Graphics,2007,176(2):121-136.

11Tse B,Barrow A,Quinn B,et al.A smoothed particle hydrodynamics algorithm for haptic rendering of dental filling materials.Proceedings of 2015 IEEE World Haptics Conference (WHC).Evanston,IL,USA.2015.321-326.

12Green S.Particle simulation using CUDA.NVIDIA Whitepaper,2010.