基于位置的流體實(shí)時(shí)交互仿真①

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

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)練過(guò)程是非常有幫助的.

目前,模擬剛體與流體的實(shí)時(shí)交互在視覺與觸力覺的真實(shí)呈現(xiàn)仍然是一個(gè)頗具挑戰(zhàn)的問(wè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),通過(guò)對(duì)粒子位置加以約束,直接更新粒子位置來(lái)模擬流體粒子,避免了需要先計(jì)算粒子受力然后才能計(jì)算粒子位置改變從而帶來(lái)的不穩(wěn)定性,極大的提升了流體模擬的效率以及真實(shí)度.

一般來(lái)說(shuō),為了實(shí)現(xiàn)剛體與流體在混合環(huán)境下的交互,我們需要為每種類型的對(duì)象使用不同的算法,并在這些異構(gòu)模型之間建立耦合接口,這將增加已經(jīng)非常復(fù)雜和耗時(shí)的模擬的復(fù)雜性和計(jì)算成本.Keiser等[5]將這一算法過(guò)程進(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)行速度有超過(guò)一個(gè)數(shù)量級(jí)的提升.2013年,Domínguez[7]等人優(yōu)化了現(xiàn)存基于CPU-GPU計(jì)算平臺(tái)的SPH計(jì)算模型,主要解決了CPU與GPU仿真數(shù)據(jù)的傳輸問(wèn)題.

本文在各種模擬流體的算法中,綜合考慮了各種算法在交互過(guò)程中的利弊,采用一種SPH的改進(jìn)算法,結(jié)合Müller的基于位置動(dòng)力學(xué)(Position Based Dynamics)[4]得到的基于位置的流體算法(Position Based Fluid)[8]來(lái)實(shí)現(xiàn)流體的基本模擬,其次通過(guò)剛性體與流體粒子的碰撞檢測(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é)的插值方法,可以通過(guò)插值計(jì)算出流體的物理特性.用粒子代替流體來(lái)模擬流體運(yùn)動(dòng),它把流體定義為在空間中離散分布并且位置可計(jì)算的粒子,這些粒子具有流體的質(zhì)量、密度、速度等物理屬性.SPH方法通過(guò)積分近似和粒子近似將流體動(dòng)力學(xué)方程進(jìn)行離散,每個(gè)粒子上的屬性函數(shù)fi(p)通過(guò)支持域內(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)通過(guò)對(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]提出方法中,可以通過(guò)混合一些約束來(lái)提高穩(wěn)定性,公式(6)修改為:

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

在每一次循環(huán)計(jì)算過(guò)程中,都需要計(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)過(guò)程中,對(duì)剛性體進(jìn)行碰撞檢測(cè),每個(gè)步驟中我們重新計(jì)算每一個(gè)粒子的領(lǐng)域粒子,并且重新計(jì)算每個(gè)粒子位置的約束值.這樣可以避免對(duì)粒子密度的低估,導(dǎo)致最終結(jié)果出現(xiàn)較大偏差.

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

根據(jù)Tse[11]給出的方法中,剛體與流體的交互過(guò)程中,工具球會(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.

工具與流體在交互過(guò)程中除受到吸附力以外,還受到因碰撞產(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è)階段通過(guò)CUDA平臺(tái)并行的對(duì)每一個(gè)粒子的狀態(tài)進(jìn)行求解,包括對(duì)初始粒子網(wǎng)格的建立到鄰域粒子的搜索以及對(duì)于每一個(gè)流體粒子各項(xiàng)物理屬性的計(jì)算,在受力計(jì)算階段實(shí)時(shí)的通過(guò)設(shè)備輸出交互力,在最后階段計(jì)算粒子的位置、速度等狀態(tài)并更新繪制顯示.同時(shí)將這些狀態(tài)用于下一個(gè)循環(huán)中的位置預(yù)測(cè)計(jì)算.

圖3 算法流程

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

算法1.流體交互算法#初始化#1.設(shè)置粒子初始位置、速度等狀態(tài)量.2.根據(jù)粒子初始位置劃定網(wǎng)格.#位置預(yù)測(cè)#3.由外力以及設(shè)備反饋力更新速度.4.預(yù)測(cè)粒子位置.#鄰域粒子查找#5.根據(jù)網(wǎng)格數(shù)據(jù)查找粒子i的鄰域粒子并記錄索引.#迭代過(guò)程#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所示,算法通過(guò)CUDA并行化,運(yùn)行在GPU端,采用的顯卡為NVIDIA Quadro M5000,包含2048個(gè)CUDA并行處理核心,所有并行計(jì)算的步驟都在GPU端執(zhí)行,核心數(shù)越多,并行加速效率越高.

交互設(shè)備采用的是Geomagic公司開發(fā)的三自由度力覺交互設(shè)備Geomagic Touch,如圖4.該反饋設(shè)備具有六自由度操作位置檢測(cè),三自由度力覺交互,通過(guò)RJ45以太網(wǎng)接口或USB接口連接終端.可以通過(guò)該設(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)過(guò)程

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

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

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

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

首先我們驗(yàn)證的是工具球在流體交互中的受力情況,圖6模擬的是通過(guò)力覺交互設(shè)備控制工具球在流體中運(yùn)動(dòng)時(shí)的受力情況,如圖所示,X軸、Y軸和Z軸正方向的力分別用實(shí)線、點(diǎn)線和虛線來(lái)表示,其方向已經(jīng)在圖6中標(biāo)出.在整個(gè)過(guò)程中我們從圖中我們可以看到工具球受力是比較平滑的曲線,沒(méi)有發(fā)生大范圍的突變,說(shuō)明在我們的算法中實(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)入,工具球還受到浮力的影響逐漸增大,到完全沒(méi)入流體,工具球受到的浮力趨于穩(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í)沒(méi)有浸入流體已經(jīng)不受浮力影響了,此時(shí)工具球還受到一個(gè)向下的吸附力,隨著工具球的進(jìn)一步向上運(yùn)動(dòng),最終工具球不受到流體交互力的影響.在整個(gè)運(yùn)動(dòng)過(guò)程中沒(méi)有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è)交互過(guò)程中,當(dāng)力覺交互設(shè)備的輸出頻率在1 kHz左右時(shí),我們通過(guò)設(shè)備能感知到連續(xù)的力覺,而不會(huì)有頓挫感.第二部分的模擬實(shí)驗(yàn)中,在粒子數(shù)目能夠達(dá)到視覺仿真要求的前提下,通過(guò)改變粒子數(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)過(guò)程中力覺交互設(shè)備的刷新率,粒子數(shù)目維持在10 k的時(shí)候基本上可以滿足我們交互過(guò)程中的精細(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í)性,說(shuō)明交互過(guò)程的實(shí)時(shí)性以及真實(shí)度達(dá)到了要求.

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

5 總結(jié)

本文探討了剛體與流體交互過(guò)程中的力覺交互,實(shí)現(xiàn)了基于位置的流體模擬算法并且模擬了剛體與流體的交互過(guò)程,實(shí)時(shí)計(jì)算并渲染出剛體工具所受到的交互力,有效的避免了傳統(tǒng)SPH算法的精度較低、計(jì)算實(shí)時(shí)性較差等劣勢(shì),同時(shí)剛性工具與流體的交互過(guò)程中很好的模擬了流體交互過(guò)程中的包括吸附力、浮力等影響比較重要的力覺.采用CUDA并行計(jì)算,極大的提高了流體算法的模擬效率以及精度.通過(guò)模擬實(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.