姚　亮，陳國(guó)棟

(福州大學(xué)物理與信息工程學(xué)院，福建福州　350116)

基于格子渦(VIC)方法的煙霧仿真研究

姚亮，陳國(guó)棟*

(福州大學(xué)物理與信息工程學(xué)院，福建福州350116)

煙霧仿真一直是計(jì)算機(jī)圖形學(xué)研究的難點(diǎn)，基于歐拉法的煙霧仿真存在嚴(yán)重的數(shù)值耗散的問(wèn)題，導(dǎo)致仿真的細(xì)節(jié)表現(xiàn)不盡如人意，基于純拉格朗日的方法計(jì)算量大實(shí)時(shí)性難以保證。為了達(dá)到較好的效果并且控制計(jì)算量提高計(jì)算速度，本文使用基于格子渦（VIC)的方法對(duì)煙霧進(jìn)行仿真。在歐拉框架下計(jì)算速度場(chǎng)，然后在拉格朗日框架下追蹤渦粒子沿跡線的運(yùn)動(dòng)。在歐拉框架下使用FFT（快速Fourier變換)快速求解Poisson方程，得到渦量輸運(yùn)方程中從渦量到速度的轉(zhuǎn)換，加快計(jì)算速度。對(duì)于粘性擴(kuò)散項(xiàng)中的粒子渦量交換，本文使用粒子強(qiáng)度交換法（PSE)來(lái)求解方程。最終實(shí)現(xiàn)了對(duì)煙霧的仿真，細(xì)節(jié)效果保留較好，實(shí)時(shí)性可以接受，驗(yàn)證了該方法的有效性。

渦量輸運(yùn)方程；格子渦（VIC)方法；快速Fourier變換（FFT)；粒子強(qiáng)度交換（PSE)

隨著計(jì)算機(jī)硬件水平的不斷提高，計(jì)算機(jī)圖形學(xué)近年來(lái)發(fā)展迅速，現(xiàn)實(shí)世界中的眾多現(xiàn)象和景物都可以在計(jì)算機(jī)中進(jìn)行仿真。然而流體一直是計(jì)算機(jī)圖形學(xué)中較為困難的仿真對(duì)象。主要是流體動(dòng)力學(xué)在概念上和計(jì)算復(fù)雜度方面的挑戰(zhàn)使得流體仿真存在困難。煙霧是流體的一部分，在電影特效，視頻游戲，虛擬手術(shù)系統(tǒng)等方面都需要產(chǎn)生煙霧效果。因此對(duì)于煙霧這樣的流體仿真仿真是計(jì)算機(jī)圖形學(xué)研究的重要方向。

經(jīng)過(guò)近幾十年的研究，中外學(xué)者提出了許多煙霧仿真的方法。歸納起來(lái)煙霧仿真的方法主要分為三類:基于粒子系統(tǒng)的煙霧仿真、基于物理模型的煙霧仿真以及基于紋理視覺(jué)效果的煙霧仿真[1]。

Reeves[2]于1983年提出了粒子系統(tǒng)的概念，利用形狀簡(jiǎn)單的微小粒子作為基本元素來(lái)構(gòu)建復(fù)雜的物體，尤其是云、水波、煙霧等自由形變物體?；谝曈X(jué)效果的方法主要是通過(guò)紋理片[3]，紋理球[4]的方法對(duì)煙霧進(jìn)行仿真。

基于物理模型的方法，通過(guò)建立煙霧的數(shù)理模型，從而對(duì)煙霧進(jìn)行仿真，精確度高，細(xì)節(jié)效果好，是近年煙霧仿真的熱點(diǎn)方法。但是由于計(jì)算復(fù)雜度較高，必須進(jìn)行適當(dāng)?shù)暮?jiǎn)化，提高仿真的實(shí)時(shí)性。本文就是基于物理模型方法來(lái)對(duì)煙霧進(jìn)行仿真?；谖锢砟Ｐ偷姆椒ǖ睦碚摶A(chǔ)是流體動(dòng)力學(xué)基本公式，著名的Navier-Stokes方程。基于Navier-Stokes方程三維流體仿真方法最早是1999年Stam[5]提出的“穩(wěn)定流體”的方法。該方法利用半拉格朗日方法大大降低了計(jì)算量，但是會(huì)帶來(lái)不可避免的數(shù)值耗散，小尺度渦旋消失太快，細(xì)節(jié)保持不是很好。Fedkiw[6]等對(duì)煙霧進(jìn)行仿真時(shí)，使用“渦旋約束”的方法把數(shù)值耗散的能量損失補(bǔ)回到流場(chǎng)中，提高了仿真效果。2004年X Wei，W li，K Mueller[7]等人將LBM的方法運(yùn)用到煙霧等氣體的仿真中，該方法不直接求解Navier-Stokes方程，而是從微觀上的遷移和碰撞中恢復(fù)宏觀的流體方程求解。近年來(lái)基于“渦絲”的方法開(kāi)始出現(xiàn)在煙霧仿真中，把煙霧描述為空間中不斷生成、運(yùn)動(dòng)、纏繞的細(xì)絲，取得了非常好的效果。

1　物理模型

煙霧是一種流體，要建立煙霧的物理模型就需要從牛頓流體的動(dòng)力學(xué)方程開(kāi)始。我們知道對(duì)于粘性可壓縮牛頓流體的動(dòng)力學(xué)方程，即可壓縮N-S方程為:

然后對(duì)方程(2)的兩邊都做旋度運(yùn)算，即可得到不可壓縮流體的渦量動(dòng)力學(xué)方程，又稱為渦量輸運(yùn)方程:

求解這個(gè)方程就可以得到對(duì)煙霧的完整描述。下面說(shuō)明方程右邊各項(xiàng)的物理意義:表示由于流場(chǎng)的速度梯度引起渦線的伸縮和彎曲，從而使渦量的大小和方向發(fā)生變化。表示渦量的粘性擴(kuò)散效應(yīng)，渦量因粘性而產(chǎn)生，因粘性而擴(kuò)散，最終因粘性而耗散。對(duì)于是比容，對(duì)于正壓流體，壓強(qiáng)是比容的單值函數(shù)，則該項(xiàng)為零，反之為斜壓流體，即有兩個(gè)獨(dú)立的熱力學(xué)變量，則將發(fā)生熱對(duì)流，引起渦量的變化。本文主要任務(wù)是對(duì)煙霧進(jìn)行仿真而不是對(duì)計(jì)算流體力學(xué)的研究。對(duì)于熱對(duì)流方面的問(wèn)題，本文不再做進(jìn)一步的討論，因此在文中將煙霧簡(jiǎn)化為正壓流體，暫時(shí)忽略對(duì)該項(xiàng)的求解。

2　物理模型的求解

2.1格子渦(VIC)方法

求解上述方程的關(guān)鍵在于確定渦元的速度，在拉格朗日觀點(diǎn)中渦量被離散為渦粒子，確定渦粒子速度的方法是使用Biot-Savart公式直接求和，若渦粒子的數(shù)目為N則有:

因此每求解一遍渦粒子的速度就要計(jì)算O(N2)次，當(dāng)N較大時(shí)計(jì)算量就變得異常巨大，因此要實(shí)現(xiàn)煙霧的實(shí)時(shí)仿真，使用快速算法勢(shì)在必行。

格子渦[8]的方法是在歐拉框架內(nèi)用流函數(shù)的Poisson方程計(jì)算速度場(chǎng)，再在拉格朗日框架內(nèi)追蹤渦粒子沿著跡線的運(yùn)動(dòng)。

首先在計(jì)算平面內(nèi)劃分網(wǎng)格，為了求解Poisson方程，先要將渦粒子攜帶的渦量分配到相應(yīng)的網(wǎng)格節(jié)點(diǎn)上，通常使用面積加權(quán)法，如圖1所示。

圖1　面積加權(quán)法示意圖

圖中1，2，3，4為網(wǎng)格節(jié)點(diǎn)，對(duì)于網(wǎng)格中的每一個(gè)節(jié)點(diǎn)k有:

其中S( k)是渦粒子在一個(gè)網(wǎng)格內(nèi)的第k個(gè)節(jié)點(diǎn)的面積權(quán)重，T是該渦元的強(qiáng)度。把所有渦粒子的渦量按如上方法分配到相應(yīng)的網(wǎng)格的各個(gè)節(jié)點(diǎn)上以后，就可以知道每個(gè)網(wǎng)格節(jié)點(diǎn)的合渦量。

2.2Helmholtz分解

對(duì)于該方程的求解有多種方法，任何一種歐拉框架的數(shù)值方法都可以用來(lái)求解Poisson方程。在渦方法中可以借助FFT(快速傅里葉變換)[9]來(lái)快速求解Poisson方程求得Ai的值，然后利用簡(jiǎn)單的差分形式就可以得到節(jié)點(diǎn)的速度:

得到每個(gè)節(jié)點(diǎn)的速度之后就要回到拉格朗日的框架中追蹤渦粒子的運(yùn)動(dòng)，根據(jù)已求得的每個(gè)網(wǎng)格節(jié)點(diǎn)的速度，再次運(yùn)用面積加權(quán)法求得各渦粒子的速度:

至此就完成了渦量到速度的轉(zhuǎn)換，由此可以計(jì)算出渦粒子在下一時(shí)刻的位置。VIC的方法最大優(yōu)點(diǎn)在于速度場(chǎng)的計(jì)算效率高。每個(gè)時(shí)間步長(zhǎng)內(nèi)只要進(jìn)行O(N)+O(M log2M)次操作，其中M是網(wǎng)格數(shù)，比直接求和法的O(N2)要小一個(gè)數(shù)量級(jí)，這就可以引入更多的粒子數(shù)目，提高仿真的分辨率，尤其是在渦量集中的區(qū)域，格子渦的方法充分發(fā)揮了拉格朗日和歐拉方法各自的優(yōu)點(diǎn)。

2.3渦流的伸縮和彎曲

圖2　有限差分法示意圖

對(duì)其離散化后為:

由此公式(4)中的伸縮和彎曲項(xiàng)求解完畢。

2.4擴(kuò)散:粒子強(qiáng)度交換

圖3　渦粒子與臨近粒子的強(qiáng)度交換

運(yùn)動(dòng)粘性系數(shù)v，每個(gè)渦粒子P與它相鄰的N個(gè)臨近渦粒子交換渦量，對(duì)于第i個(gè)渦粒子其形式為:

2.5增加煙霧效果

在流體仿真中為了達(dá)到較好的效果加入一定量的普通粒子，這被稱為示蹤劑，這些粒子不攜帶渦量，通過(guò)其在均勻網(wǎng)格中的位置，使用速度插值的方法賦予這些粒子速度，使其跟隨流體一起運(yùn)動(dòng)，從而達(dá)到較好的視覺(jué)效果。

2.6仿真流程圖(圖4)

圖4　煙霧仿真流程圖

3　仿真結(jié)果

仿真環(huán)境:

英特爾Core i5-2300 CPU＠2.8GHz 6G RAM，NVidia Quadro 600獨(dú)立顯卡1GB顯存。

Windows 8.1 64bit操作系統(tǒng)，Visual studio 2013程序采用使用C++編寫，OpenGL開(kāi)發(fā)庫(kù)GLUT 3.7版本。

3.1渦旋細(xì)節(jié)效果(圖5)

圖5　煙霧中的渦旋效果

3.2煙霧的渦環(huán)效果(圖6)

圖6　煙霧中的渦環(huán)效果，截取第10，20，30，40，50，60，70，80，90幀的效果

3.3噴射煙霧的效果(圖7)

圖7　噴射煙霧的效果截取第20，40，60，80，100，120幀的效果圖

4　總結(jié)

為了達(dá)到較好的效果并且控制計(jì)算量提高計(jì)算速度，本文使用基于格子渦(VIC)的方法對(duì)煙霧進(jìn)行仿真。在歐拉框架下計(jì)算速度場(chǎng)，然后在拉格朗日框架下追蹤渦粒子沿跡線的運(yùn)動(dòng)。在歐拉框架下使用FFT(快速Fourier變換)快速求解Poisson方程，得到渦量輸運(yùn)方程中從渦量到速度的轉(zhuǎn)換，提高了速度場(chǎng)的計(jì)算效率。對(duì)于粘性擴(kuò)散項(xiàng)中的粒子渦量交換，本文使用粒子強(qiáng)度交換法(PSE)來(lái)求解方程。本文基本實(shí)現(xiàn)了細(xì)節(jié)效果較好的煙霧仿真，但是沒(méi)有考慮燃燒場(chǎng)景。對(duì)于熱對(duì)流以及外力(如風(fēng)等)作用的效果會(huì)在今后的研究中進(jìn)一步探討，從而實(shí)現(xiàn)對(duì)灼燒煙霧的逼真仿真。

[1]夏輝麗，王繼州.煙霧模擬方法的研究綜述[J].電腦知識(shí)與技術(shù)，2014，10(21):5083-5086.

[2]William T.Reeves.Particle systems A technique for modeling a class of fuzzy objects[J].ACM Transactions on Graphics，1983.2(2): 91-108.

[3]袁雪霞，尹新富.煙霧的快速模擬.計(jì)算機(jī)工程與設(shè)計(jì)，2008，29 (9):2392-2393，2396.

[4]王繼州，袁雪霞.紋理球方法的煙霧模擬[J].小型微型計(jì)算機(jī)系統(tǒng)，2013，34(7):1680-1685.

[5]Stam J.Stable fluids[C]//Proceedings of ACM SIGGRAPH 1999.USA:ACM Press，Addison-Wesley Publishing Co.，1999:121 -128.

[6]Ronald Fedkiw，Jos Stam，Henrik Wann Jensen.Visual simulation of smoke[C]//Computer Graphics Proceedings，Annual Conference Series.Los Angeles:ACM SIGGAPH，2001:15-22.

[7]X Wei，W li，K Mueller，et al.The lattice-boltzmann method for simulating gaseous phenomena[J].IEEE Transactions on Visualization and Computer Graphics，2004，10(2):164-176.

[8]童秉綱，尹協(xié)遠(yuǎn)，朱克勤.渦運(yùn)動(dòng)理論[M].合肥:中國(guó)科學(xué)技術(shù)大學(xué)出版社，2009.

[9]Hockney R W.The Potential calculation and some applications [M].New York:Method of Comput.Phys.Academic press，1970: 135-212.

[10]Degond P.，Mas Gallic S.The Weighted Particle Method For Convection、diffusion Equations.Part I:the Case of all Isotropic Viscosity[J].Mathematics of Computation，1989，53:485-507.

[11]Degond P.，Mas-Gallic S.The Weighted Particle Method for Convection-diffusion Equations.Part II:the Antisotroic Case[J]，Mathematics of Computation，1989，53:509-525.

(責(zé)任編輯:曾晶)

Smoke Simulation Based on VIC Method

YAO Liang，CHEN Guodong*
(College of Physics and Information Engineering，F(xiàn)uzhou University，F(xiàn)uzhou 350116，China)

Fluid simulation is always a key problem for computer graphics research.Simulations using Eulerian view discretization can suffer from unwanted numerical diffusion，with fewer visual details.In order to achieve better effect and decrease computational complexity and improve computing speed，the method used in the study was based on VIC simulation of smoke.In the Eulerian framework velocity field was calculated，and then trace the route of vortex particles was traced along the path line.In the Eulerian framework，F(xiàn)FT was used to speed up the process of solving Poisson equation，convert vorticity to velocity in the vorticity transport equation，speed up the calculation.For the viscous term，considering the particle vorticity exchange，the PSE was used to solve the equation.Finally the smoke simulation was realized，details of the effect was better retained.Real-time performance is satisfactory，verifying the effectiveness of the method.

vorticity transport equation；VIC method；FFT；PSE

TP391.41

A

1000-5269(2016)01-0070-05DOI:10.15958/j.cnki.gdxbzrb.2016.01.17

2015-06-24

國(guó)家自然科學(xué)基金(61471124)；福建省自然科學(xué)基金(2013J05090)

姚亮(1989-)，男，在讀碩士，研究方向:數(shù)字媒體技術(shù)，Email:969220344＠qq.com.

陳國(guó)棟，Email:390693515＠qq.com.