一种基于涡方法的烟雾扩散动画模型

卞玉龙　王保云　周文　袁青涛

摘要：烟雾的扩散现象的模拟是一个令人感兴趣的领域。本文在分析了真实的烟雾在室内扩散的特点的基础上，提出了一种基于涡粒子法的烟雾模拟方法。充分考虑了烟雾粒子的边界条件，利用八叉树网格来划分空间，计算了涡对速度场的影响。用粒子系统实时模拟了烟雾在室内无风条件下从室内扩散到室外的效果。实验结果表明，用本文中的方法模拟的烟雾具有丰富的流体细节，和较高的计算效率。

关键词： 流体动画；烟雾；涡；粒子

中图分类号：TP311 文献标识码：A 文章编号：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 引言

在计算机图形学中，烟雾的模拟一直是一个既有意义又很有挑战的一个课题。它在各个领域都有着广泛的应用。如电影和广告特效，火灾的模拟，材料的耐火性实验等等。

由于纳维-斯托克方程组（NSES）的复杂性，使得实时模拟流体变得困难。目前，烟雾的模拟方法主要分为三大类：基于网格的方法、基于粒子的方法和网格和粒子相结合的方法。而基于涡方法的烟雾模拟有着其他方法不可比拟的优势，它可以非常逼真的模拟出具有丰富细节的流体效果。

为了观察烟雾从屋内扩散到屋外的真实效果，我们做做了真实的烟雾扩散场景，如图1所示。从图中可以看出，当烟雾扩散到室外的时候，具有明显的丝状和漩涡效果。本文用基于涡粒子法来模拟这一效果。

2 相关工作

2.1 基于网格的方法

基于网格的离散化是一种欧拉方法。最简单的网格就是均匀网格[1， 2]，就是把空间分割成长短相等的小空间。这有利于快速检索，因为可以直接根据空间中的位置来确定存储位置。但是均匀网格会造成浪费，所以有人提出用用自适应网格，或者八叉树网格[3]等等。也就是根据需要调节网格的分辨率。

基于网格的方法的优点是简单，容易实现。缺点是在计算平移项的时候，会出现数值耗散，所以，如果网格的分辨率太低，误差会被放大，从而影响到视觉效果。

2.2 基于粒子的方法

与基于网格的离散化方法相比，基于粒子的方法[4， 5， 6]是一种拉格朗日方法，大体上分为两类：SPH方法[7， 8， 9]和离散涡方法（DVM）[10]。SPH方法用流体粒子来表示流场。涡方法的特点是流体细节特别丰富。

基于粒子法的优点就是不存在数值耗散。但是粒子法不能非常精确的跟踪和重构流体的表面，导致流体表面不够平滑，于是有人引入表面跟踪算法，如level-set方法[11]。

2.3基于网格和粒子相结合的方法

为了克服了网格法和粒子法的缺点，Stam[12]首先引入半拉格朗日法，它的基本思想是，用无网格法计算平移，用网格来进行差值。这样即避免数值耗散，又提高计算效率。

3 模拟所用的基本方法

3.1 涡方法的基本数学表示

给定一个三维空间中的一个速度场[u（x，t）=（u，v，w）]，该速度场内的任一点的角速度[Ω]可以表示为以下形式：

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

涡度[ω]表示为角速度的两倍，即

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

目前，只有一个方程能很好地描述各種流体，那就是著名的Navier-Stokes方程，该方程的拉格朗日形式可以写为：

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

[??u=0] （4）

其中，[u=（u，v，w）]表示速度，[p]表示压力，[ρ]表示加速度，[f]表示外力，[?2]为拉普拉斯项。将方程（1）的两边分别求旋度得到涡方法的控制方程：

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

其中，[μ?2ω]为涡的粘性扩散，[（ω??）u]为涡的拉伸和变形。[?×?pρ2]为粒子受到的浮力，[?×f]表示涡受到的外力作用。

3.2 速度场的计算

基于涡方法的流体仿真中，需要特殊计算平移过程，因为涡度方程不会直接提供速度，但是平移的过程需要速度值，如公式（5）所示。所以在计算涡的平移的时候，需要通过速度来得到涡度。

每一个涡粒子内都有一个速度场来推动涡内的粒子绕着涡做旋转运动。并且这个涡的速度场也会影响其他区域的速度。计算单个涡对周围速度的影响的时候，我们用毕奥-萨法尔定律，表示如下：

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

如果有多个涡对同一区域都有作用，则将这些涡的作用累加。

3.3 热量的扩散

把烟雾近似地看成是一种特殊的气体。气体的密度受温度影响，密度决定气体上升或下沉。因此，我们根据这一原理建立一个温度和密度之间的联系：

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

其中，[ρ0]是在温度为[T0]的时候流体的密度。解这个方程得：[ρ=ρ0[1-α（T-T0）]]。这里，[α]是指热膨胀系数，通常定义为：[α≡1V?V?Tp=1ρ?p?Tp]。

热膨胀描述了一种材料的体积怎样随着温度的变化而变化。为了得到理想气体的热膨胀系数，把上式带入[pV=nRT]得[α=1T]。其中，[n]是指物质的量，[R]是指普通的理想气体的常数。气体的密度和温度之间的关系近似为：

[ρ=ρ0T0T] （8）

有热量的扩散就会牵扯到气体内能的变化，在空间中，内能随着时间的变化表示为：

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

其中，[e]表示内能，[q]表示热通量。当气体分子碰撞的时候，它们会把一部分动能转换为内能，[?]是指由动能转换为内能的量。

4 离散化求解方法

为了避免数值耗散，粒子的平移用拉格朗日无网格法计算。速度场的变化受到涡的影响，并且速度场中的一个位置可能受到多个涡的影响，当涡的数量很多的时候，计算量是巨大的，为了降低计算量，我们采用了一种八叉树网格，如图2所示：

在图3的树形结构中，叶子节点表示最下层网格上的涡度，上层节点表示一簇涡。先设一个关于距离的阈值[r0]，如果涡和查询点之间的距离小于[r0]，则逐一计算每个涡对该点的影响，否则用涡簇表示。根据公式（6），单个涡对周围环境的影响表示为：

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

如果把每个涡看成是点，就可以简单的表示为：

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

如果涡和所查询的影响点的距离大于阈值[r0]，则将涡簇看成是单个的涡，这个涡的涡度为这一簇涡的每个涡度的矢量和，表示如下：

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

这个涡的位置表示为：

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

5 粒子系统

这里的粒子分两种，一种是涡粒子，另一种是普通的粒子。一个涡粒子由许多个普通的粒子组成。涡粒子也有和普通粒子相同的属性，如位置、速度、大小等。涡的平移也遵循相同的规则，因此也可以把涡看成是普通的粒子。

图4 粒子系统的计算流程图

6 边界条件

6.1 涡与固体的碰撞

涡与固体碰撞的时候要满足边界无滑移和无渗透条件。假设碰撞后反射向量为[b]，在接触点固体的表面法线为[n]，在接触点的流体的运动方向为[u]。涡[ω]要同时垂直于[b]和[u]，即：

[ω=u×n] （14）

也就是说[b=ω×u=u×n×u]，另涡的半径[r=b]。最后，分配涡度[ω]大小，使涡在接触点的速度与[u]相反。所以：

[ω=4πr2uV] （15）

6.2 粒子与固体的碰撞

粒子与固体的碰撞也要满足边界无滑移和无渗透条件。烟雾粒子碰到障碍物的时候方向发生偏移，如果碰撞前烟雾粒子的运动方向为[u]，障碍物面的法线为[n]，那么烟雾粒子的因为碰撞而偏移后的速度变为为：

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

由于粒子和固体进行碰撞的时候，粒子和固体之间会有摩擦力，这个摩擦力会使粒子的反射角[b]大于入射角[a]其中，所以引入表面摩擦系数[β]，当表面为自由滑移边界的时候[β=1]，当表面为无滑移条件的时候[β=0]。

7 实验结果

本实验的硬件环境如下，处理器：英特尔Core i7 950 @ 3.07GHz 四核。内存：8GB。显卡：Nvidia Quadro FX 5800。不同涡粒子和粒子的生成速度所对应的帧率如表1所示：

得到的烟雾扩散效果图如图6所示：

图5 涡粒子法生成的烟雾扩散效果

从图5中可以看出，与真实的烟雾扩散特性相似，当烟雾从门窗冒出的时候，具有明显的丝状和漩涡效果。

8 結论及展望

本文用主要用涡粒子法实时模拟了烟雾从室内扩散到室外的效果，充分考虑了烟雾粒子和固体障碍物的碰撞作用。引入八叉树网格来计算速度场。实验结果表明，用本文的方法能够实时模拟具有丰富的流体细节的烟雾效果。

下一步的工作包括，将火焰加入到仿真中，并进一步考虑用GPU进行加速；采用更好的空间划分方法，如自适应树结构，提高划分效率。

参考文献：

[1]唐勇，李晓艳，吕梦雅，等.MacCormack方法优化烟雾模拟中NavierStokes方程对流项的求解[J].计算机辅助设计与图形学学报，2010，22（04）： 724-728.

[2]吕梦雅，翟松青，唐勇，等.风环境下烟雾的实时模拟[J].小型微型计算机系统，2012，33（03）：165-168.

[3]F. Losasso，F.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]常元章，柳有權，鲍凯，等.一种基于粒子的牛顿流体与粘弹性流体统一模拟方法[J]. 计算机学报， 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]延诃，王章野，廖斌斌，等.基于物理的海浪场景的真实感建模与绘制[J].计算机辅助设计与图形学学报， 2008， 20（09）：1117-1125.

[9]陈志宏，张正兰.动画烟雾动态模拟算法研究[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的火焰动画生成方法研究[J].计算机研究与发展， 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.