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

1(中國(guó)科學(xué)院 空間應(yīng)用工程與技術(shù)中心 太空應(yīng)用重點(diǎn)實(shí)驗(yàn)室,北京 100094)

2(中國(guó)科學(xué)院大學(xué),北京 100049)

引言

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

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

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

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

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

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

1 光滑粒子動(dòng)力學(xué)基礎(chǔ)

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

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

圖1 SPH粒子算法

2 流體模擬方法

2.1 基于位置的流體

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

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

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

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

在該Müller[4]的基于位置動(dòng)力學(xué)算法中,對(duì)粒子位置給定一個(gè)修正,所有粒子的位置滿足約束:

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

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

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

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

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

其中ε為我們?cè)谀M過程中定義的常量約束參數(shù).λ的求解修改為:

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

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

2.2 剛體與流體的交互力計(jì)算

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2 算法流程

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

圖3 算法流程

在每一個(gè)階段更為詳細(xì)的算法過程如下:

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

4 實(shí)驗(yàn)方案

4.1 實(shí)驗(yàn)設(shè)備及環(huán)境

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

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

表2 實(shí)驗(yàn)設(shè)備以及開發(fā)環(huán)境

圖4 Geomagic Touch三自由度力覺交互設(shè)備

4.2 實(shí)驗(yàn)過程

4.2.1 剛體與流體的交互力模擬驗(yàn)證

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

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

4.2.2 交互力模擬實(shí)驗(yàn)結(jié)果

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

圖5 工具在流體中的交互運(yùn)動(dòng)場(chǎng)景

圖6 工具球在流體中運(yùn)動(dòng)時(shí)的受力情況分析

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

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

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

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

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

5 總結(jié)

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

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.