基于格子涡(VIC)方法的烟雾仿真研究

姚　亮，陈国栋

(福州大学物理与信息工程学院，福建福州　350116)

基于格子涡(VIC)方法的烟雾仿真研究

姚亮，陈国栋*

(福州大学物理与信息工程学院，福建福州350116)

烟雾仿真一直是计算机图形学研究的难点，基于欧拉法的烟雾仿真存在严重的数值耗散的问题，导致仿真的细节表现不尽如人意，基于纯拉格朗日的方法计算量大实时性难以保证。为了达到较好的效果并且控制计算量提高计算速度，本文使用基于格子涡（VIC)的方法对烟雾进行仿真。在欧拉框架下计算速度场，然后在拉格朗日框架下追踪涡粒子沿迹线的运动。在欧拉框架下使用FFT（快速Fourier变换)快速求解Poisson方程，得到涡量输运方程中从涡量到速度的转换，加快计算速度。对于粘性扩散项中的粒子涡量交换，本文使用粒子强度交换法（PSE)来求解方程。最终实现了对烟雾的仿真，细节效果保留较好，实时性可以接受，验证了该方法的有效性。

涡量输运方程；格子涡（VIC)方法；快速Fourier变换（FFT)；粒子强度交换（PSE)

随着计算机硬件水平的不断提高，计算机图形学近年来发展迅速，现实世界中的众多现象和景物都可以在计算机中进行仿真。然而流体一直是计算机图形学中较为困难的仿真对象。主要是流体动力学在概念上和计算复杂度方面的挑战使得流体仿真存在困难。烟雾是流体的一部分，在电影特效，视频游戏，虚拟手术系统等方面都需要产生烟雾效果。因此对于烟雾这样的流体仿真仿真是计算机图形学研究的重要方向。

经过近几十年的研究，中外学者提出了许多烟雾仿真的方法。归纳起来烟雾仿真的方法主要分为三类:基于粒子系统的烟雾仿真、基于物理模型的烟雾仿真以及基于纹理视觉效果的烟雾仿真[1]。

Reeves[2]于1983年提出了粒子系统的概念，利用形状简单的微小粒子作为基本元素来构建复杂的物体，尤其是云、水波、烟雾等自由形变物体。基于视觉效果的方法主要是通过纹理片[3]，纹理球[4]的方法对烟雾进行仿真。

基于物理模型的方法，通过建立烟雾的数理模型，从而对烟雾进行仿真，精确度高，细节效果好，是近年烟雾仿真的热点方法。但是由于计算复杂度较高，必须进行适当的简化，提高仿真的实时性。本文就是基于物理模型方法来对烟雾进行仿真。基于物理模型的方法的理论基础是流体动力学基本公式，著名的Navier-Stokes方程。基于Navier-Stokes方程三维流体仿真方法最早是1999年Stam[5]提出的“稳定流体”的方法。该方法利用半拉格朗日方法大大降低了计算量，但是会带来不可避免的数值耗散，小尺度涡旋消失太快，细节保持不是很好。Fedkiw[6]等对烟雾进行仿真时，使用“涡旋约束”的方法把数值耗散的能量损失补回到流场中，提高了仿真效果。2004年X Wei，W li，K Mueller[7]等人将LBM的方法运用到烟雾等气体的仿真中，该方法不直接求解Navier-Stokes方程，而是从微观上的迁移和碰撞中恢复宏观的流体方程求解。近年来基于“涡丝”的方法开始出现在烟雾仿真中，把烟雾描述为空间中不断生成、运动、缠绕的细丝，取得了非常好的效果。

1　物理模型

烟雾是一种流体，要建立烟雾的物理模型就需要从牛顿流体的动力学方程开始。我们知道对于粘性可压缩牛顿流体的动力学方程，即可压缩N-S方程为:

然后对方程(2)的两边都做旋度运算，即可得到不可压缩流体的涡量动力学方程，又称为涡量输运方程:

求解这个方程就可以得到对烟雾的完整描述。下面说明方程右边各项的物理意义:表示由于流场的速度梯度引起涡线的伸缩和弯曲，从而使涡量的大小和方向发生变化。表示涡量的粘性扩散效应，涡量因粘性而产生，因粘性而扩散，最终因粘性而耗散。对于是比容，对于正压流体，压强是比容的单值函数，则该项为零，反之为斜压流体，即有两个独立的热力学变量，则将发生热对流，引起涡量的变化。本文主要任务是对烟雾进行仿真而不是对计算流体力学的研究。对于热对流方面的问题，本文不再做进一步的讨论，因此在文中将烟雾简化为正压流体，暂时忽略对该项的求解。

2　物理模型的求解

2.1格子涡(VIC)方法

求解上述方程的关键在于确定涡元的速度，在拉格朗日观点中涡量被离散为涡粒子，确定涡粒子速度的方法是使用Biot-Savart公式直接求和，若涡粒子的数目为N则有:

因此每求解一遍涡粒子的速度就要计算O(N2)次，当N较大时计算量就变得异常巨大，因此要实现烟雾的实时仿真，使用快速算法势在必行。

格子涡[8]的方法是在欧拉框架内用流函数的Poisson方程计算速度场，再在拉格朗日框架内追踪涡粒子沿着迹线的运动。

首先在计算平面内划分网格，为了求解Poisson方程，先要将涡粒子携带的涡量分配到相应的网格节点上，通常使用面积加权法，如图1所示。

图1　面积加权法示意图

图中1，2，3，4为网格节点，对于网格中的每一个节点k有:

其中S( k)是涡粒子在一个网格内的第k个节点的面积权重，T是该涡元的强度。把所有涡粒子的涡量按如上方法分配到相应的网格的各个节点上以后，就可以知道每个网格节点的合涡量。

2.2Helmholtz分解

对于该方程的求解有多种方法，任何一种欧拉框架的数值方法都可以用来求解Poisson方程。在涡方法中可以借助FFT(快速傅里叶变换)[9]来快速求解Poisson方程求得Ai的值，然后利用简单的差分形式就可以得到节点的速度:

得到每个节点的速度之后就要回到拉格朗日的框架中追踪涡粒子的运动，根据已求得的每个网格节点的速度，再次运用面积加权法求得各涡粒子的速度:

至此就完成了涡量到速度的转换，由此可以计算出涡粒子在下一时刻的位置。VIC的方法最大优点在于速度场的计算效率高。每个时间步长内只要进行O(N)+O(M log2M)次操作，其中M是网格数，比直接求和法的O(N2)要小一个数量级，这就可以引入更多的粒子数目，提高仿真的分辨率，尤其是在涡量集中的区域，格子涡的方法充分发挥了拉格朗日和欧拉方法各自的优点。

2.3涡流的伸缩和弯曲

图2　有限差分法示意图

对其离散化后为:

由此公式(4)中的伸缩和弯曲项求解完毕。

2.4扩散:粒子强度交换

图3　涡粒子与临近粒子的强度交换

运动粘性系数v，每个涡粒子P与它相邻的N个临近涡粒子交换涡量，对于第i个涡粒子其形式为:

2.5增加烟雾效果

在流体仿真中为了达到较好的效果加入一定量的普通粒子，这被称为示踪剂，这些粒子不携带涡量，通过其在均匀网格中的位置，使用速度插值的方法赋予这些粒子速度，使其跟随流体一起运动，从而达到较好的视觉效果。

2.6仿真流程图(图4)

图4　烟雾仿真流程图

3　仿真结果

仿真环境:

英特尔Core i5-2300 CPU＠2.8GHz 6G RAM，NVidia Quadro 600独立显卡1GB显存。

Windows 8.1 64bit操作系统，Visual studio 2013程序采用使用C++编写，OpenGL开发库GLUT 3.7版本。

3.1涡旋细节效果(图5)

图5　烟雾中的涡旋效果

3.2烟雾的涡环效果(图6)

图6　烟雾中的涡环效果，截取第10，20，30，40，50，60，70，80，90帧的效果

3.3喷射烟雾的效果(图7)

图7　喷射烟雾的效果截取第20，40，60，80，100，120帧的效果图

4　总结

为了达到较好的效果并且控制计算量提高计算速度，本文使用基于格子涡(VIC)的方法对烟雾进行仿真。在欧拉框架下计算速度场，然后在拉格朗日框架下追踪涡粒子沿迹线的运动。在欧拉框架下使用FFT(快速Fourier变换)快速求解Poisson方程，得到涡量输运方程中从涡量到速度的转换，提高了速度场的计算效率。对于粘性扩散项中的粒子涡量交换，本文使用粒子强度交换法(PSE)来求解方程。本文基本实现了细节效果较好的烟雾仿真，但是没有考虑燃烧场景。对于热对流以及外力(如风等)作用的效果会在今后的研究中进一步探讨，从而实现对灼烧烟雾的逼真仿真。

[1]夏辉丽，王继州.烟雾模拟方法的研究综述[J].电脑知识与技术，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]袁雪霞，尹新富.烟雾的快速模拟.计算机工程与设计，2008，29 (9):2392-2393，2396.

[4]王继州，袁雪霞.纹理球方法的烟雾模拟[J].小型微型计算机系统，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]童秉纲，尹协远，朱克勤.涡运动理论[M].合肥:中国科学技术大学出版社，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.

(责任编辑:曾晶)

Smoke Simulation Based on VIC Method

YAO Liang，CHEN Guodong*
(College of Physics and Information Engineering，Fuzhou University，Fuzhou 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，FFT 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

国家自然科学基金(61471124)；福建省自然科学基金(2013J05090)

姚亮(1989-)，男，在读硕士，研究方向:数字媒体技术，Email:969220344＠qq.com.

陈国栋，Email:390693515＠qq.com.