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

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

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

引言

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

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

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

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

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

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

1 光滑粒子動力學基礎

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

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

圖1 SPH粒子算法

2 流體模擬方法

2.1 基于位置的流體

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

根據(jù)Macklin[8]的方法,為了實現(xiàn)密度計算的穩(wěn)定,在運動中不會發(fā)生較大的突變,對每一個粒子位置施加一個非線性約束系統(tǒng),約束的狀態(tài)方程為:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

A:吸附力光滑函數(shù)起始點.

B:吸附力光滑函數(shù)結(jié)束點.

ε,ζ:可調(diào)節(jié)吸附力,本文中分別為ε=3,ζ=2.

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

KDamp為阻尼系數(shù),0>KDamp>1,v0為工具球的速度,vj為與工具球碰撞的粒子j的速度,可以看出,碰撞力的大小除了與碰撞系數(shù)有關外,也和工具球與流體粒子的速度差有關,這也符合我們?nèi)粘Ｉ钪兴姷降那榫?

3 流體交互算法實現(xiàn)

3.1 流體粒子初始網(wǎng)格構(gòu)建

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

圖2 流體粒子網(wǎng)格構(gòu)建

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

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

3.2 算法流程

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

圖3 算法流程

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

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

4 實驗方案

4.1 實驗設備及環(huán)境

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

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

表2 實驗設備以及開發(fā)環(huán)境

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

4.2 實驗過程

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

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

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

4.2.2 交互力模擬實驗結(jié)果

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

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

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

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

4.2.3 交互力的連續(xù)性驗證及結(jié)果

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

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

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

5 總結(jié)

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

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.