吕梦雅，孙旭东，甄志华，唐 勇,*，郭慧玲,3

(1. 燕山大学 信息科学与工程学院，河北 秦皇岛066004； 2. 河北省计算机虚拟技术与系统集成重点实验室，河北 秦皇岛066004； 3. 河北环境工程学院 信息工程系，河北 秦皇岛066102)

0 引言

流体仿真一直是计算机图形学领域研究的热点之一，烟雾作为流体的一种重要表现形式，具有多变性以及易受外界环境影响等特点，因此对于烟雾的真实感模拟非常困难，而关于烟雾融合的实时绘制则更具挑战性。

近年来，国内外的研究者陆续投身到对烟雾的模拟当中，他们利用粒子系统、元胞自动机、扩散过程模型、光照模型等方法进行烟雾的实时性绘制。2005年，Feldman等人[1]改进了物理网格，使用不规则四面体网格搭建烟雾模型，在很大程度上提高了仿真效果的真实性，同时，该方法可以减少计算量，仿真结果也比较稳定。2007年，秦培煜等人[2]提出了一种基于物理模型的烟雾的实时数值模拟方法，利用破开算子法分别对外力项、对流项和投影项进行求解，保证了模拟的实时性，但是丢失了流体表面的部分细节。2010年，Lentine等人[3]通过粗化欧拉流体网格求解压力项，降低了求解泊松方程所需的时间，同时仍可保持流体模拟的细节，但该算法在精细网格上不能保证实时性。2010年，Kang等人[4]提出了一种将距离函数和体积分数结合起来混合多流体模拟的方法，使混溶流体之间可以平滑过渡，但是忽视了不同流体间的相互作用力。2012年，吕梦雅等人[5]提出了一种基于流体力学方程的烟雾模拟方法，实现了烟雾与风场的交互，但是对于烟雾在风场中的描述不准确。2013年，Ando等人[6]使用离散化与空间自适应流体隐式粒子相结合的方法来有效地模拟具有极端空间适应性的液体，极大地加快了大规模流体的模拟。2014年,任博等人[7]提出了一种基于光滑的多流混合粒子流体动力学方法，使得流体融合模拟的效果较为逼真，但是对于大量的粒子模拟，实时性能将会降低。2016年，唐勇等人[8]提出了一种基于欧拉模型的多烟雾自由融合实时仿真算法，能够快速逼真地表现烟雾自由融合的过程，但烟雾融合的边界之间过渡不平滑，缺乏细腻感。2018年，Sato等人[9]提出了一种通过利用不可压缩感知的现有速度场来合成湍流运动的方法，可以直观地创建具有湍流运动的高分辨率流体动画，但是该方法的模拟形态空间较为单一。2019年，Ganter等人[10]提出了一种直接体绘制方法，使用BVH构建应用聚类算法以降低叶节点复杂性，对于图像的渲染有了显著改善。2019年，Kim等人[11]提出了一种基于感知的神经风格传输算法，能够很好地将烟雾特征从自然图像传输到烟雾模拟，实现单相流体的一般内容感知操作，但是对于多相流体的混合操作，不能生成准确的烟雾形态。

因此，本文提出一种改进的多烟雾融合实时绘制算法，解决多烟雾融合时边界过渡不自然的问题，并且通过刻画小尺度可变风场，实现可变风场中多烟雾的融合。

1 引入可变风场的烟雾模型

使用基于物理的方法来模拟烟雾运动，为了准确描述可变风场中多烟雾融合的过程，通过刻画小尺度可变风场，并将描述风力的时变函数引入到N-S方程的外力项中，从而实现可变风场中多烟雾的融合。

1.1 N-S方程组的求解

在流体力学中，N-S方程是描述不可压缩流体的最为全面的方程。本文对N-S方程进行简化，采用其不可压缩非粘性形式的方程(即欧拉方程)来模拟烟雾运动的物理模型，其公式为

∇·u=0

，

(1)

，

(2)

式(1)为烟雾运动质量守恒方程，式(2)为动量守恒方程，u为速度场，p为压强，ρ为流体密度，∇为微分算子，“·”是矢量点积，f为外力项。

简化后的N-S方程包含外力项、对流项、压力项。首先计算外力项对烟雾建模的影响，然后分别计算对流项与压力项。外力项包括热浮力以及风力，当流体受到外力作用时，速度场也会动态改变，因此，在每一时间步的开始，计算外力项对速度场的影响：

ut1(x)=ut(x)+(Δt)f，

(3)

t为当前时刻，t1为经过外力项计算后的下一时刻，根据当前时刻的速度场对下一时刻的速度场实时更新。其中外力由热浮力和风力组成，即

f=fbuoy+fwind，

(4)

其中，fbuoy为热浮力，fwind为风力。热浮力的计算需考虑浮力、烟雾重力等一些因素的影响，由文献[1]得出热浮力计算式为

fbuoy=-αρz+β(T-Tamb)z，

(5)

其中，α、β都为常数，ρ为烟雾密度，z为垂直向上的分量，T为烟雾浓度，Tamb为环境温度。

对流项可使流体通过对流，在传输过程中实现相邻区域间的能量传输，对流项的公式为

(6)

其中，λ为取值范围在0到1之间的常系数，针对不同的区域，λ的取值不同，经过外力项计算，得到速度场u1。

计算流体压力项时，通过压力梯度来修正速度散度，从而使散度为0。压力泊松方程的半离散形式表示公式为

(7)

其中，t2为经过对流项计算后的时刻，通过修正散度进而计算该时刻的压力梯度，最后完成一个时间步长的计算。

1.2 可变风场的引入

为了解决多烟雾融合速率低的问题，通过引入小尺度可变风场，并将描述风力的时变函数加入到N-S方程的外力项中，从而实现可变风场中多烟雾的融合。求得可变风场的随机过程，并将其代入到式(4)中的fwind项中，从而根据fwind项进一步计算得到外力项。本文使用一种基于谱分解的非平稳过程风场仿真算法，利用非负矩阵分解的方法，充分利用了快速傅立叶变换的优点。另外，在进行谱分解后，利用快速傅立叶变换技术提高了风场仿真的效率。其功率密度谱矩阵为

(8)

其中，s(ω,t)为Hermitian矩阵，其非对角线元素为ω的复函数，如下式所示：

(9)

s(ω,t)使用Cholesky方法可分解为H(ω,t)与HT*(ω,t)的乘积，HT*(ω,t)为H(ω,t)的共轭转置矩阵。由此可得出当前时刻模拟风场的随机过程：

(10)

通过上式，即可求解基于当前时刻的随机风速，从而计算式(4)中的fwind项，进而求得每一时间步的速度场，通过调整不同时刻的速度场，改变烟雾的流速，使得烟雾融合的速率得以改善。

2 烟雾的体绘制算法

为改善不同烟雾融合时的边界过渡不自然的问题，对传统的体绘制算法进行改进，引入错切变形法对烟雾融合进行绘制，其主要思想是把三维离散数据场的直接投影分解为三维数据场的错切变换和二维图像的变形。该方法是利用空间数据结构加快空区域的跳跃速度，减少相邻颜色间空区域的间隔，保证相邻颜色之间的平滑过渡，并且利用垂直投影的方法来加速对流体的绘制。

2.1 错切变形算法

传统的体绘制算法一般分为两类。其中一类以足迹表法为例，按照对象空间排序，此方法利用空间数据结构加快了空区域的跳跃速度，进而提升了绘制有效区域的效率。另外一类以光线投射法为例，按图像空间排序，图像空间中的每个像素沿法线方向通过三维体积发射光，然后沿光线的路径进行采样，通过合成各采样点的颜色值，即可得到最终的图像。

本文对光线投射法进行改进，同时结合了足迹表法的优点进行烟雾融合的绘制，在很大程度上可以改善烟雾融合边界过渡不平滑的问题。错切变形法结合了足迹表法与光线投射法的优点，既可以利用空间数据结构加快空区域的跳跃速度，又可以利用提前终止光的方法来加速对流体的绘制，同时采用基于剪切扭曲因子分解的方法进行烟雾渲染，将体数据转换为剪切的对象空间，对每个切边进行平移与重采样。将重新采样的切片按从前到后的顺序组合到一起，然后对中间图像进行扭曲，将中间图像转换为图像空间，第二个重采样步骤生成正确的最终图像，如图1所示。

图1 错切变形算法Fig.1 The method of shearing deformation

2.2 多烟雾融合颜色的叠加

不同烟雾渲染不同的颜色及其不透明度的计算式为

Ci=C0i(α′+Fi)，

(11)

α′=f(d,|∇d|)，

(12)

其中，Ci表示第i种烟雾最终渲染的颜色，C0i代表第i种烟雾的初始颜色，Fi表示第i种颜色烟雾的贡献度，α′及其函数f代表不透明度及其标量场和梯度大小的函数，d为浓度场，|∇d|为d的梯度函数，浓度场的平流带动烟雾的颜色变化，浓度的取值为0～1，浓度为0时表示无烟，不需要绘制；浓度为1时表示此时烟雾浓度达到峰值，绘制颜色最深。使用α′及其函数来快速划分烟雾区域中的透明、半透明、不透明区域。由此可得出最终烟雾混合的计算式：

(13)

采取以上策略，实现不同颜色烟雾的融合状态。

3 烟雾融合的形态空间

为了丰富多烟雾融合时的形态空间，引入包围盒技术实现多烟雾融合区域的网格划分，增加多烟雾融合空间的表现形态，有效地减少网格模型对形态空间的限制。由于AABB包围盒的碰撞检测算法简单实用，所以采用其进行形态空间的划分。AABB包围盒的每一平面都与坐标平面平行，其表达式为

D={(x,y,z)|xmin≤x≤xmax,

ymin≤y≤ymax,zmin≤z≤zmax}，

(14)

其中，x、y、z为烟雾的坐标，xmin、xmax、ymin、ymax、zmin、zmax分别为包围盒在x、y、z轴的最小与最大坐标。

检测烟雾粒子某个点是否属于AABB包围盒范围：

boolean contains(float[]p){

return

(p[0]>=xmin[0])&&(p[0]<=xmax[0])&&

(p[1]>=ymin[1])&&(p[1]<=ymax[1])&&

(p[2]>=zmin[2])&&(p[2]<=zmax[2])

}

检测值返回一个布尔类型，若是烟雾粒子属于包围盒范围，则对该区域烟雾粒子进行物理建模，通过修改包围盒的不同形状，从而展现不同形状空间的烟雾形态。

4 实验结果与分析

本文实验的硬件：处理器为Intel Core i7 4790 3.60GHz，内存为4G，显卡为AMD Radeon R7 250，软件环境为Windows 10操作系统，unity 5.6.1，Visual Studio 2017。

图2分别展示了文献[9]与本文方法的效果图，图2(a)为文献[9]采用基于感知的神经风格传输算法，可以实现单相流体的精确模拟，但是针对多相流体的融合模拟，可以观察到不同相的烟雾在融合过程中会出现明显的分隔，图2(b)为本文的方法，针对不同相的烟雾，应用错切变形绘制算法，着重处理烟雾融合的边缘部分，可以使不同相的烟雾融合更加平滑。

图2 双股烟雾融合效果Fig.2 Effect of double smoke fusion

图3引入可变风场，有效提升局部烟雾的流速，加快烟雾融合的速率。图3(a)为未加入风场的效果图，图3(b)为加入可变风场后的效果图。当未加入风场时，三股烟雾呈自然蔓延的趋势，对比加入小尺度可变风场后，在相同的时间下，加入风场后三股烟雾的流速显著提升，从而很大程度上提高了烟雾的融合速率。

图3 加入风场与未加风场烟雾融合效果对比图Fig.3 Comparison of smoke fusion effect in a field with or without wind

图4为文献[12]与本文方法烟雾融合的细节对比，可以看到文献[12]中烟雾融合的边界过渡不平滑，因此本文对其改进。其中图4(a)为文献[12]的模拟，图4(b)为本文的烟雾融合绘制方法，针对文献[12]烟雾融合的边界过渡不平滑的问题，采用错切变形绘制算法进行烟雾的绘制，着重处理烟雾融合边缘的过渡，减少相邻颜色间空区域的间隔，可以得到很好的效果，并且用本文方法可使实验的帧率提升1倍。

图4 七彩烟雾融合的对比图Fig.4 Contrast pictures of colorful smoke fusion

图5是本文方法模拟不同条件下的七彩烟雾融合效果图，图5(a)是未添加风场的效果图，图5(b)是在图5(a)的基础上添加了可变风场，且风速系数为2的效果图，其中风力加快了烟雾的流速，提升了融合速率。

图6为文献[7]、文献[8]与本文方法对比，图6(a)为文献[7]通过布朗扩散来模拟流体的混合，但是模拟大量粒子的混合，将不能保证实时性。图6(b)是文献[8]中通过欧拉方法模拟烟雾融合，可以看到烟雾融合边界的过渡不够平滑。图6(c)为本文使用错切变形法来进行烟雾融合的绘制，着重处理烟雾融合边缘的过渡，减少相邻颜色间空区域的间隔，同时对于大量粒子的模拟，使不同烟雾融合处更平滑，过渡更自然的同时仍可保证较高的帧率。

图5 加入风场与未加风场对比Fig.5 Comparison of wind field with and without wind field

图6 流体、烟雾融合的细节对比Fig.6 Comparison of the details of fluid and smoke fusion

图7为添加了镜面以及不同数量烟雾的融合效果，在相同时间下，随着烟雾数量增多，相邻烟雾之间平滑地融合，且烟雾的融合速率并没有降低。

图8应用AABB包围盒技术，使得烟雾模拟的空间展示不再单调。本文模拟了多种空间的烟雾融合，其中图8(a)的模拟空间为立方体,图8(b)的模拟空间为空心球体，图8(c)的模拟空间为外八角方管内圆管环状体。通过引入AABB包围盒技术，实现多烟雾融合区域的网格划分，增加多烟雾融合空间的表现形态。

图7 不同条件下的烟雾融合效果Fig.7 Smoke fusion effect under different conditions

图8 不同空间下的烟雾融合效果Fig.8 Smoke fusion effect in different spaces

为检测改进的烟雾融合效率，表1为部分实验的帧率，从表1可知，使用相同网格或粒子数的情况下，本文方法在实时性上更具有优势。

表1 场景帧率统计表Tab.1 Scene frame rate statistics

5 结论

本文改进了可变风场中多烟雾融合实时绘制算法。通过加入的小尺度可变风场与改进的错切变形绘制算法，不仅使得烟雾融合的速率得到了提升，还改善了烟雾融合的效果，使得多烟雾融合时更加平滑自然。另外本文引入了AABB包围盒技术，增加了多烟雾融合时空间的表现形状。在加入风场与未加风场的对比试验中帧率增长了42%，在七彩烟雾融合场景中，本文方法使得七股烟雾之间的融合更加自然。从所设计的与前人相关研究成果的对比实验可以看出，可变风场中多烟雾融合的绘制速度得到了提升，融合边界更加自然，烟雾表现空间形态得到了丰富。