唐 勇，張利衛(wèi)，呂夢(mèng)雅，毛菊珍

（1.燕山大學(xué)信息科學(xué)與工程學(xué)院，河北秦皇島066004；2.河北省計(jì)算機(jī)虛擬技術(shù)與系統(tǒng)集成重點(diǎn)實(shí)驗(yàn)室，河北秦皇島066004）

自由場(chǎng)景下煙霧動(dòng)態(tài)追蹤的實(shí)時(shí)仿真

唐 勇1，2，*，張利衛(wèi)1，2，呂夢(mèng)雅1，2，毛菊珍1，2

（1.燕山大學(xué)信息科學(xué)與工程學(xué)院，河北秦皇島066004；2.河北省計(jì)算機(jī)虛擬技術(shù)與系統(tǒng)集成重點(diǎn)實(shí)驗(yàn)室，河北秦皇島066004）

針對(duì)煙霧在自由場(chǎng)景下跟隨物體的運(yùn)動(dòng)形態(tài)刻畫困難的問題，本文提出一種動(dòng)態(tài)追蹤方法。首先，針對(duì)追蹤煙霧因空氣阻力產(chǎn)生運(yùn)動(dòng)態(tài)勢(shì)變化問題，建立空氣阻力模型，實(shí)現(xiàn)自然逼真的煙霧運(yùn)動(dòng)；其次，構(gòu)建自適應(yīng)網(wǎng)格場(chǎng)，引入邊界盒子覆蓋網(wǎng)格區(qū)域，實(shí)時(shí)捕捉物體中心與自適應(yīng)網(wǎng)格場(chǎng)中心在可接受偏差內(nèi)轉(zhuǎn)換，實(shí)現(xiàn)煙霧的位置追蹤；再次，采用旋轉(zhuǎn)矢量法，通過影響熱浮力方向向量有效解決煙霧方向偏轉(zhuǎn)；然后，利用有限差分法快速求解N-S方程，對(duì)流項(xiàng)采用MacCormack方法實(shí)現(xiàn)無插值，減少數(shù)值耗散，提高實(shí)時(shí)性；最后，架設(shè)多種煙霧動(dòng)態(tài)追蹤實(shí)例。實(shí)驗(yàn)表明，本文方法能夠真實(shí)實(shí)時(shí)模擬出煙霧動(dòng)態(tài)追蹤運(yùn)動(dòng)物體的自然逼真效果。

煙霧動(dòng)態(tài)追蹤；空氣阻力模型；自適應(yīng)網(wǎng)格場(chǎng)；旋轉(zhuǎn)矢量；N-S方程

0　引言

煙霧等流體在自然場(chǎng)景中扮演著舉足輕重的角色，基于物理的煙霧仿真也一直是計(jì)算機(jī)圖形學(xué)與虛擬現(xiàn)實(shí)領(lǐng)域研究的熱點(diǎn)與難點(diǎn)。近年來，煙霧模擬的真實(shí)性與實(shí)時(shí)性得到不斷地提升與突破。1999年Stam結(jié)合半拉格朗日和隱式方法來近似求解流體物理方程［1］，獲得了穩(wěn)定實(shí)時(shí)的模擬結(jié)果，但因?yàn)閿?shù)值耗散導(dǎo)致小尺度漩渦消失太快。2008年，F(xiàn)rank Losasso等利用網(wǎng)格細(xì)化技術(shù)在自由的八叉樹上捕捉小尺度的視覺細(xì)節(jié)［2］，減少壓力項(xiàng)求解時(shí)間，然而不能實(shí)時(shí)確定需要細(xì)化的計(jì)算域。同年，康奈爾大學(xué)Theodore Kim等人提出一種小波湍流方法［3］，使大、小尺度細(xì)節(jié)單獨(dú)處理。2011年，斯坦福大學(xué)Michael Lentine等人改進(jìn)半拉格朗日方法［4］，將質(zhì)量守恒與動(dòng)量守恒用于流體模擬，有效改善了大時(shí)間步長(zhǎng)流體細(xì)節(jié)缺失的情況。2013年，斯坦福大學(xué)Zhu Bo等人通過動(dòng)態(tài)延伸均勻網(wǎng)格有效提升數(shù)據(jù)存儲(chǔ)與計(jì)算速度并動(dòng)態(tài)捕捉煙霧細(xì)節(jié)［5］，但延伸的網(wǎng)格會(huì)使煙霧產(chǎn)生虛影。2014年，美國(guó)威斯康星大學(xué)Rajsekhar Setaluri等人用自適應(yīng)的稀疏均勻網(wǎng)格結(jié)構(gòu)來優(yōu)化數(shù)據(jù)存儲(chǔ)與處理［6］，從而加快煙霧模擬速度，但仍然存在數(shù)值耗散的問題。

前人的研究中大部分只考慮了密閉空間中自然騰升煙霧的真實(shí)感與實(shí)時(shí)模擬，而對(duì)其跟隨運(yùn)動(dòng)物體實(shí)時(shí)產(chǎn)生的煙霧變化仍然甚少涉及。煙霧的動(dòng)態(tài)追蹤中，需要實(shí)時(shí)獲取運(yùn)動(dòng)物體的位置，其算法涉及每幀數(shù)據(jù)的獲取與合理轉(zhuǎn)換，復(fù)雜度高；運(yùn)動(dòng)物體方向的隨機(jī)變化，使得煙霧速度如果不能及時(shí)更新，會(huì)導(dǎo)致視覺上的延遲等，從而將使得煙霧的跟隨運(yùn)動(dòng)復(fù)雜多樣。

針對(duì)上述問題，本文提出一種煙霧動(dòng)態(tài)追蹤方法，以期實(shí)現(xiàn)煙霧自主跟隨運(yùn)動(dòng)物體的繪制仿真。首先，針對(duì)煙霧因空氣阻力產(chǎn)生運(yùn)動(dòng)態(tài)勢(shì)變化問題，根據(jù)尼古拉茲實(shí)驗(yàn)曲線，建立空氣阻力模型；其次，引入邊界盒子覆蓋網(wǎng)格區(qū)域，重新定義三維歐拉網(wǎng)格結(jié)構(gòu)中心，通過與運(yùn)動(dòng)物體中心數(shù)據(jù)的算法傳遞，解決煙霧的位置追蹤；然后，采用旋轉(zhuǎn)矢量法，由運(yùn)動(dòng)物體的運(yùn)動(dòng)方向作用于煙霧，建立力場(chǎng)分析方程更新其熱浮力，實(shí)現(xiàn)煙霧的方向追蹤；最后，優(yōu)化改進(jìn)的N-S流體力學(xué)方程，保證模擬的實(shí)時(shí)性。

1　煙霧物理模型的構(gòu)建

模型的構(gòu)建對(duì)煙霧的模擬舉足輕重。在流體力學(xué)領(lǐng)域中，N-S方程組是最常用的物理模型，并且是表達(dá)不可壓縮流體最全面的微分方程。合理的建立與求解煙霧物理方程至關(guān)重要。本文采用不可壓縮非粘性形式的N-S方程：

式中，u為速度，p為壓力，ρ為流體的密度，f為合外力，“·”為矢量點(diǎn)積，“?”為梯度算子。方程（1）為連續(xù)性方程保證質(zhì)量守恒，方程（2）為動(dòng)量方程保證動(dòng)量守恒。方程（1）和（2）是N-S基本方程組。

考慮到在物體帶動(dòng)下的煙霧受到空氣阻力影響較大，在此引入空氣阻力fresis以便真實(shí)刻畫煙霧的運(yùn)動(dòng)態(tài)勢(shì)。其中fbuoy為熱浮力，fconf為漩渦限制力，合外力f計(jì)算方程如下：

其中，漩渦限制力［7］表述為

為簡(jiǎn)單計(jì)，視物體為質(zhì)點(diǎn)，空氣密度均勻，空氣阻力方程為

其中，λ為阻力系數(shù)，ρ為空氣密度，l/d為物體迎風(fēng)面積，v為煙源速度。根據(jù)尼古拉茲試驗(yàn)曲線，流體阻力與速度多少次方成正比取決于流態(tài)，而流態(tài)主要取決于流速。當(dāng)空氣速度較低時(shí)，流態(tài)處于層流區(qū)λ=64/Re，其中Re為雷諾數(shù)，表達(dá)式為Re=ρvl/μ，將其帶入式（5），得到阻力fresis和速度v的一次方成正比，即

至此，方程（1）～（6）構(gòu)成了本文基本的煙霧物理方程。

2　煙霧的動(dòng)態(tài)追蹤

由于物體運(yùn)動(dòng)位置的不定性、發(fā)煙方向偏轉(zhuǎn)的未知性，使得正確模擬動(dòng)態(tài)追蹤的煙霧極其復(fù)雜。為克服未知因素而引起的系列問題，引入自適應(yīng)網(wǎng)格場(chǎng)與旋轉(zhuǎn)矢量分塊解決上述問題。

2.1煙霧的位置追蹤

每個(gè)模擬開始，在物體C的附近選擇一個(gè)邊界盒子，預(yù)先計(jì)算出盒子的中心與物體C中心的轉(zhuǎn)化，包含有一定的偏移量。設(shè)定三維坐標(biāo)系（x，y，z），即有

式中，S（x，y，z）為邊界盒子的中心，Sc（xc，yc，zc）為物體C中心，σ為二者中心轉(zhuǎn)化的偏移量。對(duì)邊界盒子所覆蓋區(qū)域自適應(yīng)離散為均勻網(wǎng)格，三維網(wǎng)格空間大小為Size.x*Size.y*Size.z。轉(zhuǎn)化圖如圖1所示，其中，立方體C即為物體C的模擬，邊緣網(wǎng)格線范圍為盒子邊界。物體C可自由移動(dòng)，其運(yùn)動(dòng)方程為

圖1　邊界盒子的形成Fig.1 Formation of the bounding box

追蹤的煙霧每幀由物體C確定其中心位置，即在每時(shí)刻t的速度v等于物體C的速度vc，則每幀邊界盒子始終與物體C保持σ的位移差，故邊界盒子的運(yùn)動(dòng)方程為

式中，Xt為t時(shí)刻邊界盒子的位置，Xt-1為t-1時(shí)刻邊界盒子的位置，Xct為t時(shí)刻物體C的位置，Xct-1為t-1時(shí)刻物體C的位置，至此構(gòu)建出實(shí)時(shí)追蹤物體C的煙霧自適應(yīng)網(wǎng)格場(chǎng)。

圖2分別為靜止火炬與運(yùn)動(dòng)中火炬煙霧不同的運(yùn)動(dòng)形態(tài)，運(yùn)動(dòng)火炬上的煙霧因受空氣阻力而發(fā)生偏轉(zhuǎn)。煙霧的自適應(yīng)網(wǎng)格跟隨火炬而運(yùn)動(dòng)，實(shí)現(xiàn)其位置的動(dòng)態(tài)追蹤。

圖2　靜動(dòng)態(tài)火炬下煙霧不同的運(yùn)動(dòng)形態(tài)Fig.2 Smoke′s movements under static and dynamic flare

在自適應(yīng)網(wǎng)格場(chǎng)中由于對(duì)一組偏微分方程構(gòu)成的N-S方程組進(jìn)行精確求解十分困難，由于有限差分法的計(jì)算區(qū)域被劃分為具有多個(gè)子區(qū)域的規(guī)則網(wǎng)格，而在規(guī)則網(wǎng)格上又具有求解快速簡(jiǎn)單等特點(diǎn)，且本文構(gòu)建的自適應(yīng)網(wǎng)格場(chǎng)為均勻網(wǎng)格，故采用有限差分法求解N-S方程實(shí)現(xiàn)煙霧的模擬，邊界采用諾伊曼邊界條件處理網(wǎng)格。

2.2煙霧的動(dòng)態(tài)方向偏轉(zhuǎn)

物體的旋轉(zhuǎn)偏移量，決定了煙霧的走向。引入旋轉(zhuǎn)矢量，通過動(dòng)態(tài)影響熱浮力方向因子，處理煙霧的方向偏轉(zhuǎn)?？紤]到煙霧粒子溫度、煙霧密度及煙霧重力作用下，熱浮力與重力的計(jì)算公式為

式中右邊依次為煙霧的熱浮力與煙霧的重力。Δt為時(shí)間步長(zhǎng)，ΔT為煙霧溫度與外界溫度差，ρ為煙霧的密度，y=（0，1，0，0）表示垂直向上的方向。α、β為常數(shù)，表示溫度、密度對(duì)煙霧速度的影響程度。

以初始靜態(tài)物體C為參考系，計(jì)算出物體中心到物體出煙口的向量ac=（xc，yc，zc），物體以角θ發(fā)生旋轉(zhuǎn)，向量ac實(shí)時(shí)隨之偏轉(zhuǎn)。歸一化ac為煙霧熱浮力的方向向量RUp，

方向向量RUp作為煙霧受控旋轉(zhuǎn)的旋轉(zhuǎn)因子，令

式中，ω為旋轉(zhuǎn)因子，控制煙霧的旋轉(zhuǎn)強(qiáng)度，RUpx，RUpy，RUpz分別表示方向向量的x、y、z分量。故煙霧的動(dòng)態(tài)方向偏轉(zhuǎn)公式為

3　復(fù)雜場(chǎng)景下煙霧動(dòng)態(tài)追蹤運(yùn)動(dòng)方程的計(jì)算過程

為保證煙霧模擬的實(shí)時(shí)性，需對(duì)N-S方程進(jìn)行優(yōu)化。對(duì)于方程（2），采用破開算子法將其拆分為對(duì)流項(xiàng)、外力項(xiàng)和壓力投射項(xiàng)并依次求解。不同階段煙霧的速度場(chǎng)用u0、u1、u2表示，則方程求解過程如下：

3）壓力項(xiàng)，即求解泊松方程?·u2=?2p［10］，方程左右分別用一階二階導(dǎo)數(shù)中心差分格式展開，后采用雅各比迭代法求解式，第l次迭代時(shí)得到，將其應(yīng)用到壓力修正方程u3=u2+?p中，得到投射作用后的速度場(chǎng)u3。

至此，不斷地進(jìn)行循環(huán)計(jì)算，直到煙霧模擬結(jié)束。

4　仿真實(shí)例

實(shí)驗(yàn)環(huán)境為Windows系統(tǒng)，在Unity3D平臺(tái)上［11］實(shí)現(xiàn)煙霧動(dòng)態(tài)追蹤的效果展示。硬件系統(tǒng)為：Intel Core i5-2320 CPU，3.0 GHz，4G內(nèi)存，顯卡為NVIDIA GeForce GTX 550 Ti。實(shí)驗(yàn)幀率在41.1幀/秒左右，能夠滿足實(shí)時(shí)性要求。

圖3為煙霧動(dòng)態(tài)追蹤受力的炸彈，受力方向?yàn)橛疑戏?5度，（a）為未受力前靜止的炸彈與煙，（b）為炸彈受力后煙霧的形態(tài)變化，（c）為兩股煙霧在同質(zhì)量的炸彈受到不同方向力的不同運(yùn)動(dòng)狀態(tài)。左側(cè)炸彈受力為F（1，2，0），右側(cè)炸彈受力為F（1，1，0），其中水平向右為x正方向，又由于漩渦力的添加，煙霧細(xì)節(jié)很好的保留下來。

圖3　煙霧動(dòng)態(tài)追蹤受力運(yùn)動(dòng)的炸彈Fig.3 Smoke dynamic tracking force movement of the bomb

圖4為采用不同shader、不同擴(kuò)散因子使得逆時(shí)針轉(zhuǎn)動(dòng)風(fēng)車上動(dòng)態(tài)追蹤的煙霧呈現(xiàn)不同的顏色與累積。紅色煙霧擴(kuò)散因子為1，藍(lán)色煙霧擴(kuò)散因子為0.998。對(duì)比可看出，擴(kuò)散強(qiáng)度影響煙霧尾跡的旋轉(zhuǎn)細(xì)節(jié)，但同時(shí)，不同的擴(kuò)散因子可營(yíng)造不同的累積效果。

圖4　不同顏色與擴(kuò)散因子的煙霧風(fēng)車Fig.4 Smoke of different colors and diffusion factors

圖5為火柴棒上動(dòng)態(tài)追蹤的火焰。其中圖（e）為火柴運(yùn)行軌跡圖，圖（a）、（b）、（c）、（d）分別展示了在不同狀態(tài)點(diǎn)，火焰隨火柴的不同運(yùn)動(dòng)形態(tài)。本文效果在分辨率為64×128×64下每幀所用時(shí)間為0.016 s，充分滿足實(shí)時(shí)性需求，且真實(shí)感效果符合自然場(chǎng)景下的火柴燃燒。

5　結(jié)束語

本文提出一種煙霧的動(dòng)態(tài)追蹤實(shí)現(xiàn)方案。構(gòu)建空氣阻力模型，使得動(dòng)態(tài)的煙霧模擬逼真自然；引入自適應(yīng)網(wǎng)格場(chǎng)與旋轉(zhuǎn)矢量法，成功地解決了煙霧動(dòng)態(tài)追蹤的位置與方向偏轉(zhuǎn)問題；利用有限差分法快速簡(jiǎn)單求解N-S方程，采用MacCormack方法計(jì)算對(duì)流項(xiàng)，實(shí)現(xiàn)無插值，減少數(shù)值耗散，很好地提高了實(shí)時(shí)性；下一步著重處理煙霧的擴(kuò)散以及高精度下提高煙霧的動(dòng)態(tài)追蹤模擬實(shí)時(shí)性問題。

圖5　動(dòng)態(tài)追蹤的火焰細(xì)節(jié)展示Fig.5 Dynamic tracing flame details

［1］Stam J.Stable fluids［C］//Proceedings of Computer Graphics，Annual Conference Series，ACM Transactions on Graphics，Los Angeles，1999，121-128.

［2］Losasso F，Gibou F，F(xiàn)edkiw R.Simulating water and smoke with an octree data structure［J］.ACM Transactions on Graphics.2004，23（3）：457-462.

［3］Kim T，Thurey N，James D.Wavelet turbulence for fluid simulation［J］.ACM Transactions on Graphics，2008，27（3）：1-6.

［4］Lentine M，Aanjaneya M，F(xiàn)edkiw R，et al.Mass and momentum conservation forfluidsimulation［C］//ProceedingsofACM Siggraph/EurographicsSymposiumonComputerAnimation，Vancouver，2011：91-100.

［5］Zhu B，Lu W，Cong M，et al.A new grid structure for domain extension［J］.ACM Transactions on Graphics，2013，32（4）：1-12.

［6］Setaluri R，Aanjaneya M，Bauer S，et al.SP grid：a sparse paged grid structure applied to adaptive smoke simulation［J］.ACM Transactions on Graphics，2014，33（6）：1-12.

［7］唐勇，吳娛，呂夢(mèng)雅，等.一種改進(jìn)的自適應(yīng)漩渦限制實(shí)時(shí)煙霧模擬［J］.小型微型計(jì)算機(jī)系統(tǒng)，2012，33（12）：2676-2679.

［8］Selle A，F(xiàn)edkiw R，Kim B，et al.An unconditionally stable MacCormack method［J］.Journal of Scientific Computing，2008，35（2/3）：350-371.

［9］Fedkiw R，Stam J，Jensen H.Visual simulation of smoke［C］//Proceedings of Computer Graphics，Annual Conference Series，ACM Siggraph，Los Angeles，2001：15-22.

［10］柳有權(quán)，王章野，朱鑒，等.基于物理的流體動(dòng)畫加速技術(shù)的研究進(jìn)展［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)報(bào)，2013，15（3）：312-321.

［11］Dai Qing，Yang Xubo.Interactive smoke simulation and rendering on the GPU［C］//Proceedings of ACM Siggraph International Conference on Virtual-Reality Continuum and Its Applications in Industry，New York，2013：177-182.

Real-time simulation of smoke dynamic tracking in free scene

TANG Yong1，2，ZHANG Li-wei1，2，Lü Meng-ya1，2，MAO Ju-zhen1，2
（1.School of Information Science and Engineering，Yanshan University，Qinhuangdao，Hebei 066004，China；2.The Key Laboratory for Computer Virtual Technology and System Integration of Hebei Province，Qinhuangdao，Hebei 066004，China）

The real-time simulation of smoke is a great challenge in computer graphics all the time.In order to realize moving ways of the smoke dynamic tracking in a free scene，some steps are presented as follows：firstly，the air resistance model is built on account of the changes of the smoke′s tracing moving ways due to air resistance；secondly，a self-adaptive grid field is established covered by the boundary box，then in order to successfully position smoke tracking，it is captured real-timely of the conversion relationship between the center of object and the center of self-adaptive field in the acceptable deviation；thirdly，rotating vector method is adopted to solve the direction-deflexion of smoke by influencing the buoyancy vector；Fourthly，finite difference method can solve the N-S equation quickly and simply，and MacCormack method is adopted to solve the convective term with no interpolation，reduction of numerical dissipation and improvement of real-time performance；Finally，a variety of instances of smoke dynamic tracking are simulated；Experiments show that the above method can realize the vividly effect of the smoke dynamic tracking.

smoke dynamic tracking；air resistance model；self-adaptive grid field；rotating vector；N-S equation

TP391.9

A DOI：10.3969/j.issn.1007-791X.2015.06.012

1007-791X（2015）06-0541-05

2015-04-15 基金項(xiàng)目：國(guó)防重大專項(xiàng)資助項(xiàng)目（2014AT09）

*唐勇（1964-），男，四川遂寧人，博士，教授，博士生導(dǎo)師，主要研究方向?yàn)橛?jì)算機(jī)圖形學(xué)、虛擬現(xiàn)實(shí)技術(shù)及應(yīng)用，Email：tangyong@ysu.edu.cn。