一種基于物理特性的雨滴飛濺模擬算法

，，

(浙江工業(yè)大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，浙江 杭州 310023)

近年來(lái),有關(guān)自然景物的模擬[1-2]已成為計(jì)算機(jī)圖形學(xué)領(lǐng)域的一個(gè)研究熱點(diǎn).降雨過(guò)程及雨滴與固體表面的交互也日益受到該領(lǐng)域研究者的重視[3-4].下落的雨滴跟固體表面接觸后,因沖擊力而向四周迸射的物理現(xiàn)象叫做飛濺.對(duì)雨滴飛濺過(guò)程的模擬，有著非常重要的意義.主要表現(xiàn)在：1) 逼真地模擬雨滴飛濺過(guò)程，是計(jì)算機(jī)圖形學(xué)領(lǐng)域中重要的研究課題，它在數(shù)字娛樂(lè)、虛擬仿真、影視及廣告中具有重要的應(yīng)用價(jià)值和廣闊的市場(chǎng)潛力；2) 模擬雨滴在植物冠層內(nèi)滴落、飛濺與沉積的動(dòng)態(tài)過(guò)程，對(duì)研究病害在植物冠層中的傳播途徑，進(jìn)而確定最佳的防治策略，具有重要的指導(dǎo)意義.

目前，已有不少學(xué)者致力于雨滴飛濺現(xiàn)象的研究.100多年前，Worthington第一次通過(guò)頻閃觀測(cè)儀拍攝飛濺水滴的照片對(duì)飛濺現(xiàn)象進(jìn)行分析[5-6].Stow C等對(duì)飛濺水滴的數(shù)量、大小和速度的關(guān)系函數(shù)進(jìn)行了數(shù)值研究[7].S.Saint-Jean等對(duì)雨滴與植物碰撞后飛濺的傳輸路徑進(jìn)行二維的數(shù)值分析[8].但以上的研究并未涉及雨滴的三維可視化模擬.K.Garg等建立一種經(jīng)驗(yàn)概率模型來(lái)模擬基于傾斜角和表面材質(zhì)的水滴飛濺現(xiàn)象[9]，并采用Kshitiz[10]中的算法來(lái)繪制雨景，但未考慮空氣阻力的影響.汪繼文等[11]實(shí)現(xiàn)了單個(gè)水滴的濺落過(guò)程，但沒(méi)對(duì)雨的濺落效果進(jìn)行繪制.Thornton等編寫(xiě)了一款可在Maya平臺(tái)上使用的插件程序，通過(guò)手工繪制，允許用戶構(gòu)建卡通風(fēng)格的液體潑濺場(chǎng)景[12].由于他們主要是針對(duì)非真實(shí)感動(dòng)畫(huà)的繪制和卡通化模擬，所以并沒(méi)有考慮雨滴行為的物理過(guò)程及雨滴的受力情況.針對(duì)目前雨水飛濺研究的現(xiàn)狀，提出一種基于物理特性的雨滴飛濺跟蹤算法.該算法考慮到空氣阻力及風(fēng)力對(duì)雨滴的影響，同時(shí)基于粒子系統(tǒng)，利用GPU加速雨場(chǎng)景的繪制.這種方法不僅能完善雨場(chǎng)繪制[13]的模擬方法，也對(duì)雨滴飛濺路徑的研究[14-16]及潑濺效果的模擬提供方法上的借鑒.

1 雨滴的模擬及雨滴受力分析

1.1 雨滴的模擬

粒子系統(tǒng)被認(rèn)為是模擬不規(guī)則模糊物體最為成功的一種圖形生成算法[17-19].它由大量不規(guī)則的微小粒子圖元組成，而每個(gè)粒子均具有生命周期和屬性值.在生命周期內(nèi)，經(jīng)歷產(chǎn)生、活動(dòng)、消亡三個(gè)階段，并按照一定的規(guī)則改變形狀，不斷運(yùn)動(dòng)[11].

該算法采用粒子系統(tǒng)對(duì)雨場(chǎng)進(jìn)行建模.在粒子系統(tǒng)中構(gòu)建2種不同的粒子：落地前的雨滴粒子和落地后的飛濺粒子.其中一個(gè)雨滴粒子包含N個(gè)飛濺粒子.雨滴粒子在空間的一定高度產(chǎn)生，由隨機(jī)函數(shù)控制，呈均勻分布，粒子數(shù)目由全局變量NRain控制.飛濺粒子在雨滴跟地面碰撞的瞬間產(chǎn)生，繼承其所屬雨滴的當(dāng)前落地的位置，粒子數(shù)目由全局變量N控制.在無(wú)風(fēng)的情況下，雨滴的運(yùn)動(dòng)方向默認(rèn)為豎直向下，而飛濺粒子的方向由其高度角和方位角兩個(gè)參數(shù)決定.兩種粒子的數(shù)據(jù)結(jié)構(gòu)可分別描述為：

Typedef struct{

float pos[3];//雨滴粒子的空間坐標(biāo)

float vel[3];//雨滴粒子的空間速度

float gra[3];//雨滴粒子的空間加速度

float alpha;//雨滴粒子的透明值

float radius;//雨滴粒子的半徑大小

int type;//雨滴粒子的類型

int life;//雨滴粒子的生命周期

PARTICLES pa[Num];//飛濺粒子

}Rain;//雨滴粒子

Typedef struct{

float p[3];//飛濺粒子的空間坐標(biāo)

float v[3];//飛濺粒子的空間速度

float g[3];//飛濺粒子的空間加速度

float alpha;//飛濺粒子的透明值

float r;//飛濺粒子的半徑大小

float elevationAngle;//粒子的高度角

float azimuthalAngle;//粒子的方位角

int life;//飛濺粒子的生命周期

} PARTICLES；//飛濺粒子

為了提高運(yùn)行效率并節(jié)省計(jì)算時(shí)間，做出如下假設(shè)：1) 雨滴在下落過(guò)程中只考慮空氣阻力、水平風(fēng)力以及重力的影響；2) 粒子在飛濺過(guò)程中只考慮重力對(duì)粒子的影響；3) 雨滴粒子和飛濺粒子在運(yùn)動(dòng)過(guò)程中無(wú)能量損失，無(wú)水分的蒸發(fā)，不考慮形變.

1.2 雨滴受力分析

為簡(jiǎn)化雨滴在空中的下落過(guò)程，假設(shè)其初始時(shí)刻的速度為零.根據(jù)受力情況，分析得到雨滴開(kāi)始做加速運(yùn)動(dòng)，當(dāng)其受到的空氣阻力與重力平衡時(shí)，轉(zhuǎn)為勻速下落.圖1為雨滴下落過(guò)程的受力簡(jiǎn)化圖.

圖1 作用在雨滴上的外力

圖1中，F(xiàn)w為雨滴在水平方向上所受的風(fēng)力，F(xiàn)a為雨滴在水平方向上所受的空氣阻力，F(xiàn)v為雨滴在豎直方向上遇到的空氣阻力，G為雨滴的重力.假設(shè)雨滴的形狀為圓球形，半徑為r，密度為ρ，質(zhì)量為m，重力系數(shù)為g，雨滴速度為v，空氣密度ρa(bǔ)=1.21 kg/m3.根據(jù)牛頓動(dòng)力學(xué)公式，可以計(jì)算出雨滴在下落過(guò)程中任意時(shí)刻的位移和速度，即

(1)

G=mg

(2)

(3)

(4)

其中：Fd為空氣阻力是Fv和Fa的和值；CD為空氣阻力系數(shù)跟雷諾系數(shù)Re成一定比例，計(jì)算公式[20]為

(5)

其中：Re=ρa(bǔ)vd/μa；d為飛濺粒子的直徑；μa為運(yùn)動(dòng)粘滯系數(shù).

在雨場(chǎng)模擬中，水平風(fēng)力的大小由隨機(jī)函數(shù)控制.由于考慮到水平風(fēng)力和空氣阻力的影響，該算法采用二階龍格庫(kù)塔法來(lái)求解雨滴行為的微分方程，計(jì)算公式為

(6)

其中參數(shù)h為時(shí)間間隔.

二階龍格庫(kù)塔法是顯示單步式，每前進(jìn)一步需要計(jì)算兩個(gè)函數(shù)值k1和k2.故適當(dāng)?shù)倪x取參數(shù)h，可使每步計(jì)算得到的結(jié)果達(dá)到所需的精度.這里時(shí)間間隔h選取為0.02 s.

2 雨滴飛濺算法

2.1 飛濺粒子屬性的初始化和更新

雨滴的飛濺過(guò)程與雨滴的物理屬性和接觸物的表面屬性有關(guān).研究表明[21]：雨滴的下落速度越大，飛濺粒子的數(shù)目越多；雨滴的半徑越大，飛濺粒子的平均半徑也越大.而物體表面的粗糙度、硬度、濕潤(rùn)度、傾斜度和親水性等物理屬性也會(huì)影響雨滴的飛濺過(guò)程.

假設(shè)半徑為r的雨滴從y軸方向自由落體跟表面S發(fā)生碰撞后，破碎過(guò)程中產(chǎn)生的飛濺粒子的數(shù)量為n.然后n個(gè)飛濺粒子將以各自的初速度v朝著不同的方向運(yùn)動(dòng).運(yùn)動(dòng)方向由高度角θ和方位角φ組成.圖2所示為雨滴粒子跟表面S碰撞后，其中一個(gè)飛濺粒子的傳播路徑.表面S的法向量為N，表面傾斜角為α.

圖2 粒子飛濺的示意圖

Gary等通過(guò)實(shí)驗(yàn)得到飛濺現(xiàn)象的數(shù)據(jù)，發(fā)現(xiàn)飛濺粒子的屬性值很大程度上符合高斯分布[9].因此，為了描述飛濺粒子的屬性，該算法采用高斯分布函數(shù)對(duì)飛濺粒子進(jìn)行初始化[22].

1) 飛濺粒子的數(shù)目

雨滴碰撞破碎后，產(chǎn)生的飛濺粒子的數(shù)目n符合高斯分布函數(shù)G(x;u,δ)，u,δ為已知變量，則

n=n0+nvar·rand()

(7)

其中：n0為平均值；nvar為方差.

2) 飛濺粒子的方向

粒子的高度角θ符合高斯分布函數(shù)G(x;u,δ)，u,δ為已知變量.方位角φ每隔x°統(tǒng)一分布，其中x=360/n.

θ=θ0+θvar·rand()

(8)

其中：θ0為平均值；θvar為方差.

3) 飛濺粒子的初速度

粒子的初速度大小v符合高斯分布函數(shù)G(x;u,δ)，其中u,δ為已知變量，則

v=v0+vvar·rand()

(9)

其中：v0為平均值；vvar為方差.

根據(jù)式(9)求得的v值，可分別得到粒子在x,y,z軸上的分速度，即

vx=v·cosθ·cosφ

(10)

vy=v·sinθ

(11)

vz=v·cosθ·sinφ

(12)

為了真實(shí)的反映飛濺的物理過(guò)程，參數(shù)n0，v0，θ0，nvar，vvar和θvar的取值參照文獻(xiàn)[9]中的實(shí)驗(yàn)數(shù)據(jù).從其他學(xué)者的物理實(shí)驗(yàn)結(jié)果[8]可知：雨滴的最小直徑為0.5 mm，最大直徑不超過(guò)6 mm，而直徑在0.5 mm和5 mm之間的雨滴能產(chǎn)生飛濺現(xiàn)象，飛濺粒子的半徑r不大于1.5 mm，飛濺粒子的速度不大于20 m/s.由此可以限定飛濺粒子屬性的取值范圍.

4) 飛濺粒子的速度更新

飛濺過(guò)程中，粒子保持x和z方向上的速度不變.由于雨滴在飛濺過(guò)程中經(jīng)歷的時(shí)間較短，所以風(fēng)力和空氣阻力對(duì)它的影響很小，可以忽略不計(jì).為增強(qiáng)粒子濺落的效果，使用變量k來(lái)控制垂直方向上的速度.故y軸方向上的速度表示為

particle[i]·v[2]=particle[i]·v[2]+

particle[i]·g[2]·t·k

(13)

其中：0≤i≤n-1；k為速度調(diào)節(jié)因子；t為時(shí)間變量.

5) 飛濺粒子的位置更新

飛濺過(guò)程中只考慮重力對(duì)粒子的影響，根據(jù)牛頓第二定律，求得飛濺粒子的運(yùn)動(dòng)軌跡為

particle[i]·p[0]=particle[i]·p[0]+

particle[i]·v[0]·t

(14)

particle[i]·p[1]=particle[i]·p[1]+

particle[i]·g[1]·t·t

(15)

particle[i]·p[2]=particle[i]·p[2]+

particle[i]·v[2]·t

(16)

式中：0≤i≤n-1;t為時(shí)間變量.

2.2 雨滴路徑跟蹤算法

跟蹤雨滴粒子飛濺路徑的算法步驟如下：

Step1初始化雨場(chǎng)的雨量、半徑、空間位置、空間加速度、空間速度和生命周期等屬性值.

Step2根據(jù)牛頓運(yùn)動(dòng)定律，計(jì)算得到雨滴下落過(guò)程的位置和速度.

Step3檢測(cè)雨滴是否跟表面發(fā)生碰撞.如發(fā)生碰撞，則將雨滴粒子的生命周期設(shè)為0，轉(zhuǎn)至Step 4；若沒(méi)發(fā)生碰撞，則轉(zhuǎn)至Step 2.

Step4雨滴粒子碰撞后破碎成飛濺粒子，修改雨滴的類型type，初始化飛濺粒子各個(gè)屬性值(坐標(biāo)、加速度、速度、高度角和方位角、生命周期等).

Step5根據(jù)飛濺粒子的初速度和受力情況，更新計(jì)算得到飛濺粒子的運(yùn)動(dòng)軌跡.

Step6判斷飛濺粒子的生命周期.若為0，粒子消亡.

Step7顯示所有現(xiàn)存的粒子，即繪制粒子.

3 基于GPU加速技術(shù)的場(chǎng)景繪制

3.1 粒子屬性的存儲(chǔ)方式

在模擬雨場(chǎng)飛濺效果時(shí)，該算法將粒子的所有屬性都存儲(chǔ)到相應(yīng)的紋理像素中.將速度和位置屬性的分量x，y，z分別存儲(chǔ)到相應(yīng)浮點(diǎn)紋理的R，G，B三個(gè)分量中[23].為避免GPU對(duì)顯示存儲(chǔ)器訪問(wèn)的讀寫(xiě)沖突，該算法采用雙緩沖區(qū)的辦法.即創(chuàng)建一對(duì)紋理貼圖，其中一塊用于輸入，以讀取信息供cg代碼計(jì)算，另一塊則作為輸出，以保存計(jì)算結(jié)果.粒子的速度和位置信息都存儲(chǔ)在一對(duì)紋理貼圖中.由于粒子的其他屬性(大小、顏色、生命周期、透明度等)不需很復(fù)雜的計(jì)算.因此，可將它們存儲(chǔ)在單一紋理內(nèi)即可.粒子屬性的存儲(chǔ)方式如圖3所示.

圖3 粒子屬性的存儲(chǔ)方式

3.2 粒子的繪制過(guò)程

粒子繪制算法分為CPU運(yùn)行階段和GPU運(yùn)行階段.粒子的產(chǎn)生和消亡需要借助CPU完成.首先，在內(nèi)存中建立一個(gè)空閑堆棧當(dāng)作索引存儲(chǔ)器.然后當(dāng)新粒子產(chǎn)生時(shí)，從索引分配器中取出一個(gè)空閑的索引號(hào)，并將粒子的屬性值寫(xiě)入到索引號(hào)指向的粒子堆棧中，這就完成了粒子的產(chǎn)生.當(dāng)粒子的生命周期結(jié)束時(shí)，將粒子的索引號(hào)壓入空閑堆棧中，即完成粒子的消亡.而粒子屬性的更新和粒子的繪制需要借助GPU完成.具體實(shí)現(xiàn)步驟如下：

Step1以紋理的形式，將粒子屬性讀入GPU，記為紋理X.

Step2根據(jù)前一幀的速度和位置值，運(yùn)用相應(yīng)的運(yùn)動(dòng)學(xué)公式計(jì)算出下一幀的速度和位置.

Step3將Step 2中最后得到的結(jié)果保存到另一個(gè)紋理中，記為紋理Y.

Step4交換紋理X和紋理Y，即將保存結(jié)果的紋理Y作為下次計(jì)算的輸入，而紋理X則作為輸出，保存下次計(jì)算的結(jié)果.

以上粒子屬性的更新是以頂點(diǎn)的表示方式在Vertex Shader中完成.最后通過(guò)Fragment Shader來(lái)完成繪制.

4 算法實(shí)現(xiàn)與結(jié)果

依據(jù)上述建模思想，在Visual Studio 2005的集成開(kāi)發(fā)環(huán)境下，該算法使用C++編程語(yǔ)言并采用OpenGL三維圖形渲染機(jī)制，對(duì)雨滴飛濺現(xiàn)象進(jìn)行可視化模擬.采用的計(jì)算機(jī)硬件配置：CPU為Inter i5-750；主頻為2.67 G；內(nèi)存為4 GB；顯卡為NVIDIA GeForce3.

圖4 單個(gè)粒子隨時(shí)間變化的潑濺效果

圖5 不同雨量的飛濺效果

圖4為單個(gè)粒子隨時(shí)間變化的潑濺效果.圖5為繪制不同雨量時(shí)的渲染效果.從圖5中可以看出：該算法逼真的模擬了不同雨量時(shí)下雨的場(chǎng)景和雨滴落地的飛濺效果，具有較好的視覺(jué)和實(shí)時(shí)效果.在雨量和場(chǎng)景幾何復(fù)雜度不同的情況下，系統(tǒng)運(yùn)行的性能也有所改變.在圖5中，場(chǎng)景的粒子系統(tǒng)共有20 000雨滴粒子，整個(gè)模擬過(guò)程繪制的速度在10～25幀/s的范圍內(nèi)，能滿足實(shí)時(shí)性的需求.

與以往的工作相比，在繪制效果上，該方法更加逼真地展現(xiàn)了雨滴飛濺水花的大小和形狀.在算法上，該方法考慮空氣阻力和水平風(fēng)力對(duì)雨滴的影響，構(gòu)建了兩種不同的粒子來(lái)描述其運(yùn)動(dòng)的物理過(guò)程，并對(duì)雨滴落地的飛濺路徑進(jìn)行跟蹤模擬.在繪制效率上，該算法利用GPU可編程管線的數(shù)據(jù)并行處理方式以及OPENGL的內(nèi)存對(duì)象，采用雙緩存的數(shù)據(jù)存儲(chǔ)結(jié)構(gòu).跟基于CPU的方法相比，大大提高了繪制效率.

5 結(jié)束語(yǔ)

通過(guò)粒子系統(tǒng)構(gòu)建雨場(chǎng)，結(jié)合統(tǒng)計(jì)學(xué)知識(shí)和前人的實(shí)驗(yàn)數(shù)據(jù)模擬雨滴的行為軌跡，提出一種基于物理特性的雨滴飛濺算法.實(shí)驗(yàn)結(jié)果表明：該算法得到了良好的視覺(jué)效果，利用GPU的并行性和可編程性，提高了實(shí)驗(yàn)的運(yùn)行速度.在未來(lái)的工作中，將進(jìn)一步完善目前的模型，考慮雨水下落和飛濺過(guò)程中雨滴自身形狀的變化.此外，表面濕潤(rùn)度、粗糙度、親水性等屬性對(duì)飛濺現(xiàn)象的影響也將應(yīng)用到雨場(chǎng)的繪制中.同時(shí)加入更多的隨機(jī)函數(shù)處理方法來(lái)提高效果的渲染，這些都有待后續(xù)進(jìn)一步的工作和完善.

參考文獻(xiàn)：

[1] 丁維龍，陳敏智，程志君．植物器官可視化建模系統(tǒng)的設(shè)計(jì)與實(shí)現(xiàn)[J]．浙江工業(yè)大學(xué)學(xué)報(bào)，2009，37(3)：295-310．

[2] 馬培良，丁維龍，古輝．基于OpenGL和雙三次貝塞爾曲面的稻葉可視化建模[J]．浙江工業(yè)大學(xué)學(xué)報(bào)，2010，38(1)：36-40．

[3] TATARCHUK N, ISIDORO J. Artist-directable real-time rain rendering in city environments[C]//In Proceedings of SIGGRAPH 2006．New York：ACM Press，2006：23-64．

[4] HUBER L, FITT B D L, MCCARTNEY H A. The incorporation of pathogen spores into rain-splash droplets[J]. Plant Pathol，1996，45(3)：506-517．

[5] WORTHINGTON A M. On the forms assumed by drops falling vertically on a horizontal plate[J]．Proc of Royal Society of London，1876，25：261-271．

[6] WORTHINGTON A M. A second paper on the forms assumed by drops of liquids falling vertically on a horizontal plate[J]. Proc of Royal Society of London，1877，25：498-503．

[7] STOW C D, STAINER R D. The physical products of a splashing water drop[J]. Meteorological Society of Japan Journal，1977，55(5)：518-532．

[8] SAINT-JEAN S, CHELLE M, HUBER L. Modelling water transfer by rain-splash in a 3D canopy using Monte Carlo integration[J]．Agricultural and Forest Meteorology，2004，121(1/2)：183-196．

[9] GRAG K, KRISHNAN G, NAYAR S K. Material based splashing of water drops[C]//In Proceedings of Eurographics Symposium on Rendering. Switzerland: Eurographics Association，2007：171-182．

[10] KSHITIZ G, SHREE K N. Photorealistic rendering of rain streaks[J]. ACM Trans on Graphics，2006，25：996-1002．

[11] 汪繼文，陸和軍．基于粒子系統(tǒng)的水滴濺落模擬[J]．科學(xué)技術(shù)與工程，2010，10(6)：1547-1550．

[12] THORNTON J. Directable simulation of stylized water splash effects in 3D space[C]//In Proceedings of SIGGRAPH 2006-Material Presented at the ACM SIGGRAPH 2006 Conference. Boston: ACM Press，2006：94-94．

[13] PIERRE R, VINCENT J. Realistic real-time rain rendering[J]. Computer & Graphics，2006，30(4)：507-518．

[14] YANG Meng, HUANG Meng-cheng, YANG Gang, et al. Physically-Based animation for realistic interactions between tree branches and raindrops[C]//In Proc of the 17thACM Symposium on Virtual Reality Software and Technology (ACM VRST 2010)．New York：ACM Press，2010：83-86．

[15] YANG Meng, JIANG Long-sheng, LI Xiao- sheng, et al. Interactive coupling between a tree and raindrops[J]. Computer Animation and Virtual Worlds，2012，23(3/4)：267-277．

[16] WANG H, MUCHA J P, TURK G. Water drops on surfaces[J]. ACM Trans on Graphics，2005，24(3)：921-929．

[17] XU Y Z, LEE K, KIM E, et al. Simulation of smoke to improve unity 3D games engine particle system based on FDS[J]． Soft Technology，2013，19：183-186．

[18] KIM M, CHEN G N. Dynamic particle system for mesh extraction on the GPU[C]//In Proceedings of the 5th Annual Workshop． New York：ACM Press，2012：38-46．

[19] DING Wei, ZHU Zhong-bin. Real-time rain and snow rendering[C]//In 2013 Second International Conference on Agro-Geoinformatics．Fairfax：Agro-Geoinformatics Press，2013：32-35．

[20] MERCER G, SWEATMAN W L, ELVIN A, et al. Process driven models for spray retention by plants[C]//In Proceedings of the 2006 Mathematics in Industry Study Group. New Zealand: Mathematics in Industry Study Group Press，2006：57-85．

[21] STOW C, HADFIELD M. An experimental investigation of fluid flow resulting from the impact of a water drop with an unyielding dry surface[J]. Proc of Royal Society of London 1981,373：419-441．

[22] MUNDO C, SOMMERFELD M, TROPEA C. Droplet-wall collisions: experimental studies of the deformation and breakup process[J]. International Journal Multiphase Flow，1995，21(2)：151-173．

[23] 許楠，郝愛(ài)民，王莉莉．一種基于GPU的粒子系統(tǒng)[J]．計(jì)算機(jī)工程與應(yīng)用，2006，19(3)：77-79．