卞玉龍　王保云　周文　袁青濤

摘要：煙霧的擴(kuò)散現(xiàn)象的模擬是一個(gè)令人感興趣的領(lǐng)域。本文在分析了真實(shí)的煙霧在室內(nèi)擴(kuò)散的特點(diǎn)的基礎(chǔ)上，提出了一種基于渦粒子法的煙霧模擬方法。充分考慮了煙霧粒子的邊界條件，利用八叉樹網(wǎng)格來(lái)劃分空間，計(jì)算了渦對(duì)速度場(chǎng)的影響。用粒子系統(tǒng)實(shí)時(shí)模擬了煙霧在室內(nèi)無(wú)風(fēng)條件下從室內(nèi)擴(kuò)散到室外的效果。實(shí)驗(yàn)結(jié)果表明，用本文中的方法模擬的煙霧具有豐富的流體細(xì)節(jié)，和較高的計(jì)算效率。

關(guān)鍵詞： 流體動(dòng)畫；煙霧；渦；粒子

中圖分類號(hào)：TP311 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1009-3044（2016）36-0274-03

A Model Of Smoke Diffusion Animation Based On Vertex Method

BIAN Yu-long，WANG Bao-yun，ZHOU Wen，YUAN Qing-tao

（ZhengYuan Geographic Information Co.， LTD， Beijing 101300，China）

Abstract：Smoke diffuse simulation is an interesting area in the field of computer graphic. After analyzing the appearance of smoke diffuse out of the room， a new method which is based on the vortex method is presented in this paper. Fully taking into account the boundary condition of particle， a Lagrangian method is used to advect particles. In order to compute the influence of vortex on the velocity of other field， an Octree grid is used to divide the space. We simulate the appearance of smoke spreading out of a room in real-time. The result shows that our result has a good visual effect， and is efficient.

Key words：fluid animation； smoke； vortex； particle

1 引言

在計(jì)算機(jī)圖形學(xué)中，煙霧的模擬一直是一個(gè)既有意義又很有挑戰(zhàn)的一個(gè)課題。它在各個(gè)領(lǐng)域都有著廣泛的應(yīng)用。如電影和廣告特效，火災(zāi)的模擬，材料的耐火性實(shí)驗(yàn)等等。

由于納維-斯托克方程組（NSES）的復(fù)雜性，使得實(shí)時(shí)模擬流體變得困難。目前，煙霧的模擬方法主要分為三大類：基于網(wǎng)格的方法、基于粒子的方法和網(wǎng)格和粒子相結(jié)合的方法。而基于渦方法的煙霧模擬有著其他方法不可比擬的優(yōu)勢(shì)，它可以非常逼真的模擬出具有豐富細(xì)節(jié)的流體效果。

為了觀察煙霧從屋內(nèi)擴(kuò)散到屋外的真實(shí)效果，我們做做了真實(shí)的煙霧擴(kuò)散場(chǎng)景，如圖1所示。從圖中可以看出，當(dāng)煙霧擴(kuò)散到室外的時(shí)候，具有明顯的絲狀和漩渦效果。本文用基于渦粒子法來(lái)模擬這一效果。

2 相關(guān)工作

2.1 基于網(wǎng)格的方法

基于網(wǎng)格的離散化是一種歐拉方法。最簡(jiǎn)單的網(wǎng)格就是均勻網(wǎng)格[1， 2]，就是把空間分割成長(zhǎng)短相等的小空間。這有利于快速檢索，因?yàn)榭梢灾苯痈鶕?jù)空間中的位置來(lái)確定存儲(chǔ)位置。但是均勻網(wǎng)格會(huì)造成浪費(fèi)，所以有人提出用用自適應(yīng)網(wǎng)格，或者八叉樹網(wǎng)格[3]等等。也就是根據(jù)需要調(diào)節(jié)網(wǎng)格的分辨率。

基于網(wǎng)格的方法的優(yōu)點(diǎn)是簡(jiǎn)單，容易實(shí)現(xiàn)。缺點(diǎn)是在計(jì)算平移項(xiàng)的時(shí)候，會(huì)出現(xiàn)數(shù)值耗散，所以，如果網(wǎng)格的分辨率太低，誤差會(huì)被放大，從而影響到視覺(jué)效果。

2.2 基于粒子的方法

與基于網(wǎng)格的離散化方法相比，基于粒子的方法[4， 5， 6]是一種拉格朗日方法，大體上分為兩類：SPH方法[7， 8， 9]和離散渦方法（DVM）[10]。SPH方法用流體粒子來(lái)表示流場(chǎng)。渦方法的特點(diǎn)是流體細(xì)節(jié)特別豐富。

基于粒子法的優(yōu)點(diǎn)就是不存在數(shù)值耗散。但是粒子法不能非常精確的跟蹤和重構(gòu)流體的表面，導(dǎo)致流體表面不夠平滑，于是有人引入表面跟蹤算法，如level-set方法[11]。

2.3基于網(wǎng)格和粒子相結(jié)合的方法

為了克服了網(wǎng)格法和粒子法的缺點(diǎn)，Stam[12]首先引入半拉格朗日法，它的基本思想是，用無(wú)網(wǎng)格法計(jì)算平移，用網(wǎng)格來(lái)進(jìn)行差值。這樣即避免數(shù)值耗散，又提高計(jì)算效率。

3 模擬所用的基本方法

3.1 渦方法的基本數(shù)學(xué)表示

給定一個(gè)三維空間中的一個(gè)速度場(chǎng)[u（x，t）=（u，v，w）]，該速度場(chǎng)內(nèi)的任一點(diǎn)的角速度[Ω]可以表示為以下形式：

[Ω（x，t）=12?×u（x，t）=12ijk??x??y??zuvw] （1）

渦度[ω]表示為角速度的兩倍，即

[ω（x，t）=2Ω（x，t）=?×u（x，t）] （2）

目前，只有一個(gè)方程能很好地描述各種流體，那就是著名的Navier-Stokes方程，該方程的拉格朗日形式可以寫為：

[DuDt=μ?2u-?pρ+f] （3）

[??u=0] （4）

其中，[u=（u，v，w）]表示速度，[p]表示壓力，[ρ]表示加速度，[f]表示外力，[?2]為拉普拉斯項(xiàng)。將方程（1）的兩邊分別求旋度得到渦方法的控制方程：

[DωDt=（ω??）u+μ?2ω+?ρ×?pρ2+?×f] （5）

其中，[μ?2ω]為渦的粘性擴(kuò)散，[（ω??）u]為渦的拉伸和變形。[?×?pρ2]為粒子受到的浮力，[?×f]表示渦受到的外力作用。

3.2 速度場(chǎng)的計(jì)算

基于渦方法的流體仿真中，需要特殊計(jì)算平移過(guò)程，因?yàn)闇u度方程不會(huì)直接提供速度，但是平移的過(guò)程需要速度值，如公式（5）所示。所以在計(jì)算渦的平移的時(shí)候，需要通過(guò)速度來(lái)得到渦度。

每一個(gè)渦粒子內(nèi)都有一個(gè)速度場(chǎng)來(lái)推動(dòng)渦內(nèi)的粒子繞著渦做旋轉(zhuǎn)運(yùn)動(dòng)。并且這個(gè)渦的速度場(chǎng)也會(huì)影響其他區(qū)域的速度。計(jì)算單個(gè)渦對(duì)周圍速度的影響的時(shí)候，我們用畢奧-薩法爾定律，表示如下：

[du=ω×r4πr3] （6）

如果有多個(gè)渦對(duì)同一區(qū)域都有作用，則將這些渦的作用累加。

3.3 熱量的擴(kuò)散

把煙霧近似地看成是一種特殊的氣體。氣體的密度受溫度影響，密度決定氣體上升或下沉。因此，我們根據(jù)這一原理建立一個(gè)溫度和密度之間的聯(lián)系：

[ρ0-ρρ0=α（T-T0）] （7）

其中，[ρ0]是在溫度為[T0]的時(shí)候流體的密度。解這個(gè)方程得：[ρ=ρ0[1-α（T-T0）]]。這里，[α]是指熱膨脹系數(shù)，通常定義為：[α≡1V?V?Tp=1ρ?p?Tp]。

熱膨脹描述了一種材料的體積怎樣隨著溫度的變化而變化。為了得到理想氣體的熱膨脹系數(shù)，把上式帶入[pV=nRT]得[α=1T]。其中，[n]是指物質(zhì)的量，[R]是指普通的理想氣體的常數(shù)。氣體的密度和溫度之間的關(guān)系近似為：

[ρ=ρ0T0T] （8）

有熱量的擴(kuò)散就會(huì)牽扯到氣體內(nèi)能的變化，在空間中，內(nèi)能隨著時(shí)間的變化表示為：

[ρDeDt=-??q-p（??u）+?] （9）

其中，[e]表示內(nèi)能，[q]表示熱通量。當(dāng)氣體分子碰撞的時(shí)候，它們會(huì)把一部分動(dòng)能轉(zhuǎn)換為內(nèi)能，[?]是指由動(dòng)能轉(zhuǎn)換為內(nèi)能的量。

4 離散化求解方法

為了避免數(shù)值耗散，粒子的平移用拉格朗日無(wú)網(wǎng)格法計(jì)算。速度場(chǎng)的變化受到渦的影響，并且速度場(chǎng)中的一個(gè)位置可能受到多個(gè)渦的影響，當(dāng)渦的數(shù)量很多的時(shí)候，計(jì)算量是巨大的，為了降低計(jì)算量，我們采用了一種八叉樹網(wǎng)格，如圖2所示：

在圖3的樹形結(jié)構(gòu)中，葉子節(jié)點(diǎn)表示最下層網(wǎng)格上的渦度，上層節(jié)點(diǎn)表示一簇渦。先設(shè)一個(gè)關(guān)于距離的閾值[r0]，如果渦和查詢點(diǎn)之間的距離小于[r0]，則逐一計(jì)算每個(gè)渦對(duì)該點(diǎn)的影響，否則用渦簇表示。根據(jù)公式（6），單個(gè)渦對(duì)周圍環(huán)境的影響表示為：

[ux14πω（x'）×rr3dx'] （10）

如果把每個(gè)渦看成是點(diǎn)，就可以簡(jiǎn)單的表示為：

[ux14πi=1Nωi×riri3] （11）

如果渦和所查詢的影響點(diǎn)的距離大于閾值[r0]，則將渦簇看成是單個(gè)的渦，這個(gè)渦的渦度為這一簇渦的每個(gè)渦度的矢量和，表示如下：

[ω=i=1Nωi] （12）

這個(gè)渦的位置表示為：

[P=i=1NPiωi] （13）

5 粒子系統(tǒng)

這里的粒子分兩種，一種是渦粒子，另一種是普通的粒子。一個(gè)渦粒子由許多個(gè)普通的粒子組成。渦粒子也有和普通粒子相同的屬性，如位置、速度、大小等。渦的平移也遵循相同的規(guī)則，因此也可以把渦看成是普通的粒子。

圖4 粒子系統(tǒng)的計(jì)算流程圖

6 邊界條件

6.1 渦與固體的碰撞

渦與固體碰撞的時(shí)候要滿足邊界無(wú)滑移和無(wú)滲透條件。假設(shè)碰撞后反射向量為[b]，在接觸點(diǎn)固體的表面法線為[n]，在接觸點(diǎn)的流體的運(yùn)動(dòng)方向?yàn)閇u]。渦[ω]要同時(shí)垂直于[b]和[u]，即：

[ω=u×n] （14）

也就是說(shuō)[b=ω×u=u×n×u]，另渦的半徑[r=b]。最后，分配渦度[ω]大小，使渦在接觸點(diǎn)的速度與[u]相反。所以：

[ω=4πr2uV] （15）

6.2 粒子與固體的碰撞

粒子與固體的碰撞也要滿足邊界無(wú)滑移和無(wú)滲透條件。煙霧粒子碰到障礙物的時(shí)候方向發(fā)生偏移，如果碰撞前煙霧粒子的運(yùn)動(dòng)方向?yàn)閇u]，障礙物面的法線為[n]，那么煙霧粒子的因?yàn)榕鲎捕坪蟮乃俣茸優(yōu)闉椋?/p>

[Dir=u-β?2?u?n?nn×n 0≤β≤1] （16）

由于粒子和固體進(jìn)行碰撞的時(shí)候，粒子和固體之間會(huì)有摩擦力，這個(gè)摩擦力會(huì)使粒子的反射角[b]大于入射角[a]其中，所以引入表面摩擦系數(shù)[β]，當(dāng)表面為自由滑移邊界的時(shí)候[β=1]，當(dāng)表面為無(wú)滑移條件的時(shí)候[β=0]。

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

本實(shí)驗(yàn)的硬件環(huán)境如下，處理器：英特爾Core i7 950 @ 3.07GHz 四核。內(nèi)存：8GB。顯卡：Nvidia Quadro FX 5800。不同渦粒子和粒子的生成速度所對(duì)應(yīng)的幀率如表1所示：

得到的煙霧擴(kuò)散效果圖如圖6所示：

圖5 渦粒子法生成的煙霧擴(kuò)散效果

從圖5中可以看出，與真實(shí)的煙霧擴(kuò)散特性相似，當(dāng)煙霧從門窗冒出的時(shí)候，具有明顯的絲狀和漩渦效果。

8 結(jié)論及展望

本文用主要用渦粒子法實(shí)時(shí)模擬了煙霧從室內(nèi)擴(kuò)散到室外的效果，充分考慮了煙霧粒子和固體障礙物的碰撞作用。引入八叉樹網(wǎng)格來(lái)計(jì)算速度場(chǎng)。實(shí)驗(yàn)結(jié)果表明，用本文的方法能夠?qū)崟r(shí)模擬具有豐富的流體細(xì)節(jié)的煙霧效果。

下一步的工作包括，將火焰加入到仿真中，并進(jìn)一步考慮用GPU進(jìn)行加速；采用更好的空間劃分方法，如自適應(yīng)樹結(jié)構(gòu)，提高劃分效率。

參考文獻(xiàn)：

[1]唐勇，李曉艷，呂夢(mèng)雅，等.MacCormack方法優(yōu)化煙霧模擬中NavierStokes方程對(duì)流項(xiàng)的求解[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2010，22（04）： 724-728.

[2]呂夢(mèng)雅，翟松青，唐勇，等.風(fēng)環(huán)境下煙霧的實(shí)時(shí)模擬[J].小型微型計(jì)算機(jī)系統(tǒng)，2012，33（03）：165-168.

[3]F. Losasso，F(xiàn).Gibou and R.Fedkiw， Simulating Water and Smoke with an Octree Data Structure[J]. SIGGRAPH， 2004（23）：457-462.

[4]Reeves W T.Particle systems-a technique for modeling a class of fuzzy objects[J].Computer Graph， 1983， 17（3）： 359-376.

[5]常元章，柳有權(quán)，鮑凱，等.一種基于粒子的牛頓流體與粘彈性流體統(tǒng)一模擬方法[J]. 計(jì)算機(jī)學(xué)報(bào)， 2010， 33（07）： 1286-1295.

[6]Yuanzhang Chang，Kai Bao，Youquan Liu，Jian Zhu，and Enhua Wu.A particle-based method for viscoelastic fluids animation[J].VRST' 09： Proceedings of the 16th ACM Symposium on Virtual Reality Software and Technology，2009， 111-117.

[7]Bo Zhu，Xubo Yang and Ye Fan.Creating and Preserving Vortical Details in SPH Fluid[J].Pacific Graphics 2010， 29（07）.

[8]延訶，王章野，廖斌斌，等.基于物理的海浪場(chǎng)景的真實(shí)感建模與繪制[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)， 2008， 20（09）：1117-1125.

[9]陳志宏，張正蘭.動(dòng)畫煙霧動(dòng)態(tài)模擬算法研究[J].計(jì)算機(jī)仿真， 2009， 26（02）： 217-220.

[10]A. Angelidis.Simulation of Smoke Based on Vortex Filament Primitives[J].SCA' Proc.of the Symposium on Computer Animation， 2005：87-96.

[11]洪義，王兆其，朱登明，等.基于Level Set的火焰動(dòng)畫生成方法研究[J].計(jì)算機(jī)研究與發(fā)展， 2010， 47（11）：1849-1856.

[12]J.Stam.Stable fluids[J]. Proceedings of the 26th Annual Conference on Computer Graphics and Interactive Techniques， SIGGRAPH 99. New York： ACM Press/Addison Wesley Publishing Co， 1999. 121-128.