基于粒子的水沸騰模擬

劉洋洋, 朱曉臨, 梁欣鑫, 范承凱(合肥工業(yè)大學(xué) 數(shù)學(xué)學(xué)院,安徽 合肥 230009)

基于粒子的水沸騰模擬

劉洋洋, 朱曉臨, 梁欣鑫, 范承凱
(合肥工業(yè)大學(xué) 數(shù)學(xué)學(xué)院,安徽 合肥 230009)

近年來(lái),定位流體(position based fluids,PBF)方法以其高效的計(jì)算速率和真實(shí)的模擬效果成為模擬大場(chǎng)景流體運(yùn)動(dòng)的主流方法。文章基于PBF方法使用少量的流體粒子模擬水的沸騰場(chǎng)景。首先采用流體粒子作為熱量傳輸?shù)拿浇?通過(guò)SPH方法中的核函數(shù)理論處理模型中的熱傳導(dǎo),使得模擬區(qū)域的熱量分布更連續(xù),相比于傳統(tǒng)的以網(wǎng)格作為熱量傳輸媒介的加熱方式,模擬效果更好;其次在流體粒子模型的屬性計(jì)算中引入粒子效果精度,加快模型計(jì)算速度。實(shí)驗(yàn)結(jié)果驗(yàn)證了該模型具有更好的模擬效果及穩(wěn)定性。

流體模擬;沸騰；定位流體方法;氣泡

自然現(xiàn)象的流體中,水流的場(chǎng)景模擬以其豐富的視覺(jué)感在電影特效和游戲制作里具有重要的應(yīng)用價(jià)值。目前流體模擬中較流行的2種方法分別是歐拉網(wǎng)格法和拉格朗日粒子法。歐拉網(wǎng)格法是固定在模擬對(duì)象所處的空間上的模擬對(duì)象在固定網(wǎng)格單元上運(yùn)動(dòng)。文獻(xiàn)[1]使用交錯(cuò)網(wǎng)格開(kāi)創(chuàng)了一種穩(wěn)定性較差的三維水流模擬方法; 文獻(xiàn)[2]模擬三維流體運(yùn)動(dòng)時(shí)介紹了一種無(wú)條件穩(wěn)定的方法來(lái)近似求解流體運(yùn)動(dòng)方程,該方法結(jié)合了半拉格朗日方法和隱式求解方法,因此可以得到實(shí)時(shí)穩(wěn)定的模擬結(jié)果,但是半拉格朗日方法卻會(huì)造成大量數(shù)值耗散而導(dǎo)致一些典型的“漩渦”效果消失太快。這一缺點(diǎn)后來(lái)被文獻(xiàn)[3]采用高階插值和渦流限制解決。盡管如此,作為基于網(wǎng)格的模擬方法主要的局限性仍是普遍的計(jì)算復(fù)雜度太高。拉格朗日粒子法是通過(guò)追蹤一系列任意分布的“節(jié)點(diǎn)”來(lái)求解具有各種邊界條件的積分方程或偏微分方程組,從而得到精確穩(wěn)定的數(shù)值解。光滑粒子動(dòng)力學(xué)方法(smoothed particle hydrodynamics,SPH)是最早提出的一種粒子方法,文獻(xiàn)[4-5]將該方法應(yīng)用到連續(xù)固體力學(xué)和流體力學(xué)中,并將SPH方法解釋為核函數(shù)法,成功模擬了流場(chǎng)中的激波強(qiáng)間斷現(xiàn)象。文獻(xiàn)[6]構(gòu)建了一種可以對(duì)流體的自由表面交互模擬的SPH模型;此后,SPH方法才正式廣泛應(yīng)用于模擬不可壓縮流體以及多相流體模擬。然而由于SPH本身精度和穩(wěn)定性等原因,流體模擬的效果很一般。文獻(xiàn)[7]提出了一種新穎的基于定位動(dòng)態(tài)學(xué)的流體模擬方法(position based dynamics，PBD),通過(guò)物體運(yùn)動(dòng)的限制方程來(lái)修正物體的位置,這種方式的運(yùn)算效率很高,模擬流體的運(yùn)動(dòng)也具有很真實(shí)的動(dòng)態(tài)感。文獻(xiàn)[8]結(jié)合SPH方法和定位動(dòng)態(tài)學(xué)方法提出了定位流體方法(position based fluids，PBF),通過(guò)粒子密度限制方程來(lái)校正位置,模型中人工添加的壓力項(xiàng)很好地避免了粒子小范圍飛濺或聚集的問(wèn)題,模擬效果很好,適用于大規(guī)模的場(chǎng)景模擬。

模擬流體內(nèi)部的氣泡運(yùn)動(dòng)同樣是流體模擬中的一個(gè)重要組成部分,而氣泡模擬中常用的方法有volume of fluids(VOF)[9]和水平集方法(level set method)[10-11],它們同屬于基于歐拉網(wǎng)格下構(gòu)造的界面捕獲類方法,并且都用一個(gè)標(biāo)量函數(shù)來(lái)描述界面的幾何特征,可以很好地處理拓?fù)浣Y(jié)構(gòu)。本文采用粒子方法來(lái)模擬流體內(nèi)部氣泡的運(yùn)動(dòng)過(guò)程,氣泡粒子在上升中逐漸變大并融合附近的氣泡粒子。這種氣泡模擬方式無(wú)須考慮氣泡形狀的變化,可以加速模型的計(jì)算效率。

1 定位流體簡(jiǎn)介

模擬不可壓縮流體的關(guān)鍵是如何在快速的模擬計(jì)算條件下保持流體對(duì)象的某一物理量(壓強(qiáng)、密度或者能量等)守恒,例如N-S(Navier-Stokes)方程是通過(guò)控制流體運(yùn)動(dòng)時(shí)的動(dòng)量和能量守恒來(lái)精確描述流體運(yùn)動(dòng),同時(shí)它也是流體模擬的物理模型中典型的控制方程,但是其偏微分方程組的數(shù)值解法復(fù)雜而且耗時(shí),難以保證實(shí)時(shí)模擬的要求。文獻(xiàn)[8]提出的定位動(dòng)態(tài)學(xué)方法通過(guò)控制流體粒子運(yùn)動(dòng)時(shí)的密度保持不變,先預(yù)測(cè)模擬對(duì)象的位置,再通過(guò)數(shù)次位置的迭代校正,最終滿足給定的密度限制方程。

在PBF方法的粒子模型中,密度限制要求模擬中的粒子密度值保持不變。首先令Ci(p)=ρi/ρ0-1,其中，ρ0為粒子的靜止密度；粒子i的密度計(jì)算公式為：

其中,h為支集半徑，粒子i的支集半徑表示以粒子i的坐標(biāo)為中心,h為半徑的一個(gè)球形空間范圍;p=[p1,p2,…,pn]T,pj(j=1,2,…,n)表示粒子i支集半徑內(nèi)其他粒子的位置;mj為粒子j的質(zhì)量;W1(pi-pj,h)為核函數(shù)。在下面的計(jì)算過(guò)程中用到的核函數(shù)W1、W2、W3參考文獻(xiàn)[7]。然后確定粒子的位置修正值Δp,新位置滿足如下密度限制方程:

Ci(p+Δp)=0

(1)

顯然Δp變化最快的方向和函數(shù)梯度Ci(p)的方向一致,因此可以確定一個(gè)標(biāo)量λi滿足：

(2)

(1)式的泰勒一階近似為:

Ci(p+Δp)≈Ci(p)+Ci(p)Δp=0

(3)

將(2)式代入(3)式,解得:

(4)

其中，M=diag(m1,m2,…,mn)表示鄰近粒子的質(zhì)量構(gòu)成的對(duì)角矩陣。為了避免計(jì)算(4)式時(shí)分母為0,在分母上增加一個(gè)調(diào)節(jié)參數(shù)ε(ε>0),(4)式改為:

(5)

將(5)式代入(2)式得:

為了使粒子的新位置更快地滿足密度限制方程,可以把支集半徑內(nèi)的粒子密度限制系數(shù)λj加入到粒子的位置更新中,同時(shí)引入一個(gè)非負(fù)的人工壓力項(xiàng)s以避免粒子的“聚簇”現(xiàn)象,其數(shù)值參數(shù)同文獻(xiàn)[8]。

綜上,Δp的最終計(jì)算如下:

(6)

對(duì)(6)式數(shù)次迭代后,新位置的粒子在給定的誤差范圍內(nèi)滿足密度限制方程,同時(shí)加入的人工壓力項(xiàng)使得表面的模擬效果更自然。

為了減小計(jì)算中數(shù)值耗散對(duì)模擬的影響,給每個(gè)粒子加入渦流限制(vorticity confinement),粒子添加的旋轉(zhuǎn)力形式如下:

ωi=W3(pi-pj,h),

vij=vj-vi,η=|ωi|,N=,

(7)

最后給每個(gè)粒子添加適當(dāng)?shù)酿ば粤?以保證連續(xù)運(yùn)動(dòng)的實(shí)時(shí)效果:

(8)

2 模型框架

因?yàn)楸疚牡哪M場(chǎng)景是水逐漸加熱至沸騰產(chǎn)生氣泡,所以本文模型包括熱傳導(dǎo)模型和粒子模型,其中粒子模型有流體粒子模型和氣泡粒子模型。流體粒子額外加入一個(gè)溫度屬性T,當(dāng)流體粒子溫度達(dá)到極大值時(shí),流體粒子分解形成氣體粒子,在壓力作用下氣體粒子快速上升并在液體表面處“破裂消失”。

2.1 熱傳導(dǎo)模型

本文模型將水槽流體粒子模擬區(qū)域網(wǎng)格化,每個(gè)網(wǎng)格只有一個(gè)溫度屬性T(x,y,z),初始值為20,最大值為100;然后虛擬一個(gè)在水槽底面網(wǎng)格上的熱源區(qū)域Ω(X±φ，0，Z±φ),其中(X,0,Z)為熱源的中心位置；φ和φ分別控制熱源在X和Z方向上的加熱范圍。水槽內(nèi)部網(wǎng)格的熱量傳輸遵從熱傳遞原理:

=αΔ2T

(9)

其中,α為一個(gè)和熱量傳輸速率相關(guān)的常量,本文取α=0.1。

熱傳導(dǎo)粒子模型中的受熱對(duì)象是粒子,熱量也是通過(guò)粒子的運(yùn)動(dòng)進(jìn)行擴(kuò)散。首先假定與受熱區(qū)域直接接觸的流體粒子溫度與熱源溫度相同,然后通過(guò)SPH方法中的光滑核函數(shù)可以近似得到區(qū)域內(nèi)任意位置的熱量大小,將傳統(tǒng)的小立方體表示的離散型溫度分布轉(zhuǎn)化為空間連續(xù)型的溫度分布。

使用SPH方法近似后,(9)式變?yōu)椋?/p>

(10)

其中，α′為一個(gè)和熱傳導(dǎo)相關(guān)的調(diào)節(jié)參數(shù)；xi、xj分別為粒子i和粒子j的空間位置?？紤]到SPH方法自身存在的數(shù)值耗散以及在邊界處計(jì)算精度不足等引起的熱量丟失,將(10)式最終改為如下形式:

(11)

假設(shè)模擬區(qū)域內(nèi)粒子的總數(shù)量為M,由(11)式可得粒子模型熱傳導(dǎo)方法的計(jì)算復(fù)雜度為O(M)。

在真實(shí)的水沸騰場(chǎng)景中,水分子的熱量直接影響其內(nèi)能,進(jìn)而使得水分子劇烈運(yùn)動(dòng);而在流體沸騰的模型中,同樣也要設(shè)置流體粒子的溫度對(duì)其速度值的增幅。本文使用的速度增幅為:

(12)

2.2 氣泡粒子模型

本文的氣泡模擬包括氣泡生成、氣泡上升、氣泡破碎的模擬。在氣泡上升的過(guò)程中,其體積逐漸增加,碰撞多個(gè)氣泡融合生成一個(gè)新的大氣泡粒子,氣泡在液體表面“爆炸”后推動(dòng)附近的流體粒子,產(chǎn)生液面沸騰效果。

2.2.1 氣泡生成

在加熱過(guò)程中,水由于溶解度的降低會(huì)在受熱區(qū)域逐漸析出空氣形成氣泡,因此模型中受熱網(wǎng)格的溫度在達(dá)到最大溫度后,與之相接觸的流體粒子會(huì)在下一時(shí)刻“汽化”為氣泡粒子。初始生成的氣泡粒子半徑固定,在上升的過(guò)程中隨著壓力的減小逐漸變大,半徑變化公式為:

r=αh+β(y-h)tb,

其中,h為氣泡粒子初始生成高度;y為當(dāng)前高度;tb為氣泡的生命;α、β為2個(gè)調(diào)節(jié)參數(shù),本文取值分別為0.75、0.25。

2.2.2 氣泡上升

氣泡粒子的運(yùn)動(dòng)過(guò)程同樣遵循PBF方法中的密度限制方程,但其質(zhì)量相對(duì)較小,因此在周圍流體粒子的“壓力”作用下會(huì)緩慢上升?，F(xiàn)實(shí)中,除了流體本身的限制外,流體中的氣泡在生成后的上升過(guò)程中,由于氣泡區(qū)域的壓力作用導(dǎo)致氣泡粒子會(huì)靠近周圍氣泡數(shù)量更密集的區(qū)域,因此在氣泡粒子運(yùn)動(dòng)的水平面上添加適當(dāng)?shù)奈恢眯Ｕ蠈?shí)際的模擬情形。本文氣泡粒子在水平面上的運(yùn)動(dòng)方程如下:

(13)

其中,X+為目標(biāo)氣泡粒子支集半徑內(nèi)X坐標(biāo)大于目標(biāo)氣泡粒子的其他粒子的數(shù)量;X-為目標(biāo)氣泡粒子支集半徑內(nèi)X坐標(biāo)小于目標(biāo)氣泡粒子的其他粒子的數(shù)量;Z+為目標(biāo)氣泡粒子支集半徑內(nèi)Z坐標(biāo)大于目標(biāo)氣泡粒子的其他粒子的數(shù)量;Z-為目標(biāo)氣泡粒子支集半徑內(nèi)Z坐標(biāo)小于目標(biāo)氣泡粒子的其他粒子的數(shù)量;κ為一個(gè)調(diào)節(jié)系數(shù),本文取κ=0.1。

2.2.3 氣泡融合

2個(gè)氣泡粒子融和的觸發(fā)條件為粒子間距小于2個(gè)粒子半徑之和,即d≤r1+r2。當(dāng)一個(gè)氣泡粒子需要與2個(gè)以上的粒子融合時(shí),則出現(xiàn)多個(gè)粒子融合情形。粒子融合的過(guò)程嚴(yán)格遵從體積守恒、動(dòng)量守恒。新氣泡粒子的位置、半徑及上升速度如下:

s=(s1+s2+…+sn)/n,

其中,融合后的氣體粒子初始速度方向向上。

2.2.4 氣泡破碎

氣泡粒子在離開(kāi)網(wǎng)格模擬氣泡區(qū)域進(jìn)入液體表面模擬區(qū)域時(shí)即開(kāi)始執(zhí)行氣泡破碎程序。若氣泡半徑r≤rs(rs為人工設(shè)置的一個(gè)常量)時(shí),則氣泡直接“消失”;若氣泡半徑r>rs,則氣泡粒子發(fā)生“爆炸”,記錄爆炸粒子位置(x,ymax,z),然后以該位置為中點(diǎn),kR(k為人工控制量)為半徑,搜索范圍內(nèi)的流體粒子,對(duì)其產(chǎn)生“推動(dòng)”效果。設(shè)流體粒子坐標(biāo)為(x′,y′,z′),則其對(duì)應(yīng)的分量速度增量為:

(14)

其中,λ1、λ2為調(diào)節(jié)系數(shù),本文分別為0.2、0.3。

3 粒子效果精度

文獻(xiàn)[12]在SPH方法的基礎(chǔ)上引入了粒子精細(xì)度的概念,通過(guò)粒子在模擬中的位置、黏性力和壓力3個(gè)因素確定粒子的精細(xì)度,并根據(jù)粒子精細(xì)度來(lái)分配粒子下一時(shí)刻的支集半徑。類似地,本文加入了粒子效果精度F∈(0,1],根據(jù)粒子當(dāng)前時(shí)刻的效果精度值確定下一時(shí)刻目標(biāo)粒子的計(jì)算屬性,其數(shù)值大小反映了粒子在模擬效果中的重要程度。

本文基于PBF粒子方法的液體粒子模型,對(duì)渲染效果有影響的因素如下: ① 液體粒子是否處于液體表面;② 液體粒子的動(dòng)量。

對(duì)于因素①,當(dāng)粒子處于流體表面時(shí),粒子效果精度取值為1;其他情形和粒子的高度值成正比。然而在動(dòng)態(tài)的實(shí)時(shí)模擬中,每一時(shí)間間隔內(nèi)確定模擬流體的表面函數(shù)是一個(gè)很耗時(shí)的過(guò)程,因此本文直接采用高度值來(lái)判斷粒子是否處于流體表面,估計(jì)流體表面略下方的粒子平均高度值為H。這種判斷方式簡(jiǎn)單且計(jì)算快,適用于流體區(qū)域在高度場(chǎng)的運(yùn)動(dòng)變化幅度較小的場(chǎng)景,如文中模擬的流水沸騰。對(duì)于因素②,根據(jù)液體粒子的動(dòng)量來(lái)計(jì)算粒子效果精度,而不直接采用粒子的速度值作為參考因素,這是考慮到以后的粒子模型的粒子種類會(huì)更豐富;沒(méi)有將粒子的溫度作為參考因素是因?yàn)榱Ｗ訙囟纫呀?jīng)對(duì)粒子速度進(jìn)行了增幅,所以選擇最合適的液體粒子的動(dòng)量作為計(jì)算粒子效果精度的因素。綜上所述,液體粒子效果精度的計(jì)算公式為:

(15)

其中,yi為目標(biāo)粒子的高度值;Ei為目標(biāo)粒子的動(dòng)量;Emax為目標(biāo)粒子的支集半徑內(nèi)所有粒子動(dòng)量的最大值;α+β=1,α>0，β>0。

在PBF方法模擬流體當(dāng)前幀的計(jì)算全部結(jié)束后,通過(guò)最終計(jì)算出的每個(gè)粒子的效果精度值,可以動(dòng)態(tài)調(diào)節(jié)粒子的支集半徑和迭代次數(shù),進(jìn)而影響最終模擬的效果和計(jì)算效率。具體如下:

Ri=FPBFR,Ji=FPBFJ,

其中,Ri、Ji分別為目標(biāo)粒子下一時(shí)刻計(jì)算的支集半徑和迭代次數(shù);R、J為2個(gè)定值,分別表示支集半徑和迭代次數(shù)的最大值。

同樣,本文的水蒸氣粒子模型也對(duì)SPH方法的模擬粒子加入了粒子效果精度的屬性,但是相關(guān)的影響因素不同于文獻(xiàn)[12]。一般在煙霧模擬場(chǎng)景中可以影響最終效果的因素是粒子的密度和壓力2個(gè)屬性;密度決定了粒子對(duì)周圍其他粒子的作用大小;壓力決定了粒子的運(yùn)動(dòng)軌跡。因此,本文的水蒸氣粒子效果精度計(jì)算如下:

(16)

其中,max{ρi} 、max{Pi}分別為目標(biāo)粒子支集半徑內(nèi)所有粒子的最大密度值和壓力值。目標(biāo)粒子的壓力值計(jì)算公式如下:

其中，p(xj)為目標(biāo)粒子的壓強(qiáng)。

得到水蒸氣粒子的效果精度值后就可以確定下一時(shí)刻每個(gè)粒子的支集半徑，即

Ri′=FSPHR′,

其中,Ri′為目標(biāo)粒子下一時(shí)刻屬性計(jì)算的支集半徑;R′為定值,表示支集半徑的最大值。

文獻(xiàn)[12]中的實(shí)驗(yàn)?zāi)M在SPH模型中加入了粒子精細(xì)度后雖然可以提速86%，但是模擬的細(xì)節(jié)卻沒(méi)有得到提升。本文模型對(duì)PBF模型加入了改進(jìn)的粒子精細(xì)度控制方法，在多次的調(diào)試數(shù)據(jù)模擬實(shí)驗(yàn)后，最終的結(jié)果是在同等模擬場(chǎng)景下，本文方法可以將PBF模型的計(jì)算速度提升16%～24%，同時(shí)還提高了模擬的細(xì)節(jié)效果。

4 程序步驟

本文模型的具體步驟如下：

(1) 初始化粒子模型后,在水槽底部中心區(qū)域網(wǎng)格進(jìn)行加熱。

(2) 判斷流體粒子的溫度T。若T≥100,則執(zhí)行步驟(3);否則執(zhí)行步驟(4)。

(3) 將溫度達(dá)到極限的流體粒子在周圍空間隨機(jī)生成4個(gè)半徑固定的氣體粒子。

(4) 為每個(gè)粒子添加外力:

vi+Δtfext(pi)/mi→vi,

其中,fext(pi)為粒子i受到的重力。

(5) 預(yù)測(cè)下一時(shí)刻粒子的位置:

(6) 搜尋每個(gè)粒子支集半徑內(nèi)的所有粒子。

(7) 根據(jù)(5)式計(jì)算支集半徑內(nèi)每個(gè)粒子的λi。

(8) 根據(jù)(6)式計(jì)算支集半徑內(nèi)每個(gè)粒子的Δpi。

(10) 判斷是否達(dá)到迭代步數(shù),否則回到步驟(4)。

(11) 更新每個(gè)粒子的速度:

特別地，對(duì)氣體粒子還要加入(13)式。

(13) 判斷氣泡粒子是否需要“破碎”。若是,則根據(jù)(14)式添加相關(guān)流體粒子的增量;否則執(zhí)行步驟(14)。

(14) 由(12)式更新氣體粒子的半徑大小。

(15) 更新水槽網(wǎng)格的溫度。

(16) 由步驟(15)的溫度值及(12)式最后更新所有流體粒子的速度值。

(17) 根據(jù)流體粒子和氣體粒子的位置用OpenGL進(jìn)行渲染。

(18) 重復(fù)步驟(2)。

5 實(shí)驗(yàn)結(jié)果

本文仿真模型實(shí)驗(yàn)平臺(tái)為Windows 7系統(tǒng),程序開(kāi)發(fā)的軟件選用Microsoft Visual Studio 2008,三維圖形的編程接口和最終的粒子渲染采用OpenGL。算法實(shí)現(xiàn)的硬件環(huán)境為1臺(tái)Intel(R) Core(TM) i3-4130的CPU、4 G內(nèi)存的PC機(jī)。為了驗(yàn)證本文模型的優(yōu)點(diǎn),本文將PBF的粒子運(yùn)動(dòng)程序換為傳統(tǒng)的SPH方法,具體的模擬結(jié)果如圖1所示。

圖1 沸騰模型基于不同方法幀數(shù)為67的表面效果

在程序運(yùn)行中通過(guò)截取圖片測(cè)得圖1中67幀次前后若干圖的幀率，見(jiàn)表1所列。

由表1可見(jiàn),本文的模型在基于SPH的粒子方法時(shí)模擬的平均幀率約為1.31,而PBF方法約為1.96,相比于傳統(tǒng)的SPH方法,計(jì)算效率提高了49.6%。

表1 SPH方法和PBF方法的幀率比較

流體沸騰模擬中,從氣泡生成到氣泡消失過(guò)程的氣泡效果如圖2所示。

圖2 流體沸騰模型中的氣泡模擬效果

6 結(jié) 論

本文在PBF方法的粒子框架基礎(chǔ)上加入了氣泡粒子,并通過(guò)2種不同粒子間的相互運(yùn)動(dòng)模擬了流體加熱沸騰的場(chǎng)景。相比于傳統(tǒng)的SPH方法,本文模型的計(jì)算速度更快,表面的沸騰效果更真實(shí),模擬中加入的氣泡模擬也很好地提高了模擬場(chǎng)景的細(xì)節(jié)。然而由于粒子模型中另一種粒子的加入使得穩(wěn)定性有所降低,本文模型模擬的時(shí)間步長(zhǎng)受到了限制,表面模擬的效果變化劇烈以致氣泡粒子在即將出水時(shí)的“破裂”效果不明顯,后續(xù)的研究工作將主要放在氣泡模擬方面,采用新的粒子運(yùn)動(dòng)方法或者新的氣泡粒子框架來(lái)加強(qiáng)流體沸騰表面的模擬。

[1] FOSTER N，METXAS D.Realistic animation of liquids.[J].Graphical Models & Image Processing，1996，58(5):204-212.

[2] STAM J.Stable fluids[J].Acm transactions on graphics，2001，1999:121-128.

[3] FEDKIW R，Stam J，Jensen H W.Visual simulation of smoke[C]//Proceedings of the 28th Annual Conference on Computer Graphics and Interactive Techniques.[S.l.]:ACM，2001:15-22.

[4] MONAGHAN J J.An introduction to SPH[J].Computer Physics Communications，1988，48(1):89-96.

[5] MONAGHAN J J.Why particle methods work[J].Siam Journal on Scientific & Statistical Computing，1982，3(4):422-433.

[9] HONG J M，KIM C H.Animation of bubbles in liquid [J].Computer Graphics Forum，2003，22(3):253-262.

[10] OSHER S，SETHIAN J A.Fronts propagating with curvature-dependent speed:algorithms based on hamilton-jacobi formulations[J].Journal of Computational Physics，1988，79(1):12-49.

[11] ZHENG W，YONG J H，PAUL J C.Simulation of bubbles[C]//Proceedings of the 2006 ACM SIGGRAPH/Eurographics Symposium on Computer Animation.Aire-la-Ville,Switzerland:Eurographics Association,2006:325-333.

[12] 譚詩(shī)瀚,段茗,楊紅雨.非均勻粒子流體模擬[J].計(jì)算機(jī)工程與設(shè)計(jì)，2011，32(8):2760-2763.

(責(zé)任編輯 張 镅)

Boiling water simulation based on particles

LIU Yangyang, ZHU Xiaolin, LIANG Xinxin, FAN Chengkai
(School of Mathematics， Hefei University of Technology， Hefei 230009， China)

In recent years， the method of position based fluids(PBF) has become the mainstream in simulating huge scene of fluid motion for its efficient calculation and perfect effect. In this paper， the model is established to simulate the boiling water by using few of fluid particles based on PBF. In the model， the fluid particles are taken as the medium of heat conduction to handle heat exchange with the theory of kernel function in smoothed particle hydrodynamics(SPH)， which makes the distribution of heat in the simulating area more continuous and therefore obtains better simulation result compared to the traditional method using the grid as the medium of heat transference. Then the fluid particles are given a concept of precision of effect, which accelerates the computation and enhances the details on the surface. The experimental result verifies better simulation effect and stability of the model.

fluid simulation; boiling; position based fluids(PBF) method; bubble

2016-03-22；

2016-06-09

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

劉洋洋(1989-),男,安徽蚌埠人,合肥工業(yè)大學(xué)碩士生; 朱曉臨(1964-),男,安徽池州人,博士,合肥工業(yè)大學(xué)教授,碩士生導(dǎo)師.

10.3969/j.issn.1003-5060.2017.06.025

TP301

A

1003-5060(2017)06-0854-06