CN112347708B - 一种基于物理过程的实时光学烟幕的渲染方法 - Google Patents

Abstract

本发明实施例提供的一种基于物理过程的实时光学烟幕的渲染方法，利用光照烟雾模型以及红外烟雾模型渲染烟雾的过程中，根据背景辐射，外部光源，烟雾自身透过率等因素对烟雾渲染的影响，计算获得密度网格、温度网格以及亮度网格，同时，红外模型引入了自发辐射模型，基于密度网格、温度网格以及亮度网格计分别计算烟雾的透明度以及亮度，基于射线的亮度和透明度对前时间帧的当前位置视线方向的像素进行混合渲染，在兼顾效率的同时提高了烟雾渲染的视觉真实性和光照准确性。

Description

一种基于物理过程的实时光学烟幕的渲染方法

技术领域

本发明属于光学烟幕领域，具体涉及一种基于物理过程的实时光学烟 幕的渲染方法。

背景技术

烟雾模拟是计算流体力学中的一个关于不可压缩低速流动的子问题。 这类低速流体在现实的自然环境中广泛存在，为了提高场景的渲染真实感， 如何高效、真实地渲染这类低速流体成了计算机图形学中不可绕过的热点 问题。生活中最常见的低速流体就是烟雾和流水，两者都是不可压缩的NS 方程的解，主要的区别是水体多了一个液面的边界条件，这导致烟雾更注 重整体的外形，而流水更注重表面的变化，这导致两者在网格划分和粒子 处理上具有不同的优化手段和策略。烟雾的类型有很多，生活中燃烧有机 物生成的黑烟，核电站冷却塔蒸法形成的白色雾气，军事上为了遮蔽目标 使用的各种可见光、红外烟幕。现实中存在于有大量的烟雾，如何在实时 的虚拟世界中更好地将烟雾真实地还原并渲染出来是一个值得不断探索和 研究的问题。

现有技术建立了基于NS方程建立了火焰的流体模型，但该模型没有考 虑压力项，且求解算泊松方程的算法比较原始，结果的真实性和运行效率 都有所欠缺。后面研究人员将流体力学模型应用于红外烟幕实时仿真中， 其流体力学模型采用了Stable Fluid模型，使用CUDA作为开发工具加速流 体的仿真，达到了实时性的要求，但是其红外光照模型相对简单，渲染效 果较差，与现实的烟雾差距较大，仿真结果不够真实。

发明内容

为了解决现有技术中存在的上述问题，本发明提供了一种基于物理过 程的实时光学烟幕的渲染方法。本发明要解决的技术问题通过以下技术方 案实现：

本发明实施例提供的一种基于物理过程的实时光学烟幕的渲染方法， 包括：

S1，获取初始时间帧的前一时间帧网格结构中烟雾发生器的位置、浓 度和温度以及碰撞器的边界条件，初始化当前时间帧的网格结构；

其中，所述网格结构是直角坐标规则的矢量面心网格，所述网格结构 中每个单元网格的6个面心位置存储该位置的速度，单元网格中心位置存 储密度，压强和温度；

S2，基于已获取的当前时间帧的前一时间帧的速度场，使用半拉格朗 日法回溯计算获得当前时间帧的速度场；

S3，基于网格结构的旋量场的涡旋力、烟雾的浮力、重力以及当前时 间帧的速度场，计算得到当前时间帧添加外力的速度场；

S4，基于当前时间帧添加外力的速度场，使用隐式的后向欧拉法迭代 计算获得添加粘滞力的速度场；

S5，基于添加粘滞力的速度场和第一边界条件，获得满足不可压缩条 件的速度场；

S6，基于当前时间帧的前一时间帧的浓度、温度以及满足不可压缩条 件的速度场，使用半拉格朗日法回溯计算获得当前时间帧的密度场和温度 场；

S7，确定添加第二边界条件约束的密度场和温度场，获得密度网格以 及温度网格；所述密度网格是在所述密度场作用下的网格结构，所述温度 网格是在所述温度场作用下的网格结构；

S8，针对当前时间帧，基于红外光线建立的亮度网格，将从每个像素 出发的射线携带的透过率初始化为1，将每个像素携带的亮度初始化为0；

S9，在射线进入亮度网格并按照固定步长移动时，根据温度网格采样 获得的温度值以及密度网格采样获得的密度值，计算亮度网格的自发辐射 值；

S10，计算太阳光射线进入所述亮度网格的太阳光辐射；

S11，基于自发辐射、太阳光辐射、消光系数以及透过率，计算射线在 当前时间帧的当前位置的透明度和亮度；

S12，重复步骤S9至S11，直到射线离开亮度网格范围；

S13，基于射线的亮度和透明度对前时间帧的当前位置视线方向的像素 进行Alpha混合渲染。

本发明实施例提供的一种基于物理过程的实时光学烟幕的渲染方法， 利用光照烟雾模型以及红外烟雾模型渲染烟雾的过程中，根据背景辐射， 外部光源，烟雾自身透过率等因素对烟雾渲染的影响，计算获得密度网格、 温度网格以及亮度网格，同时，红外模型引入了自发辐射模型，基于密度 网格、温度网格以及亮度网格计分别计算烟雾的透明度以及亮度，基于射 线的亮度和透明度对前时间帧的当前位置视线方向的像素进行混合渲染，在兼顾效率的同时提高了烟雾渲染的视觉真实性和光照准确性。

以下将结合附图及实施例对本发明做进一步详细说明。

附图说明

图1是本发明实施例提供的一种一种基于物理过程的实时光学烟幕的 渲染方法的流程图；

图2是本发明实施例提供的基于光线投射的体绘制渲染结构示意图；

图3是本发明实施例提供的基于物理过程的实时光学烟雾的渲染过程 的实现流程图。

具体实施方式

下面结合具体实施例对本发明做进一步详细的描述，但本发明的实施 方式不限于此。

实施例一

如图1所示，本发明实施例提供的一种基于物理过程的实时光学烟幕 的渲染方法，包括：

S1，获取初始时间帧前一时间帧网格结构中烟雾发生器的位置、浓度 和温度以及碰撞器的边界条件，初始化当前时间帧的网格结构；

其中，所述网格结构是直角坐标规则的矢量面心网格，所述网格结构 中每个单元网格的6个面心位置存储该位置的速度，单元网格中心位置存 储密度，压强和温度；

S2，基于已获取的当前时间帧的前一时间帧的速度场，使用半拉格朗 日法回溯计算获得当前时间帧的速度场；

S3，基于所述网格结构的旋量场的涡旋力、烟雾的浮力、重力以及当 前时间帧的速度场，计算得到当前时间帧添加外力的速度场；

S4，基于当前时间帧添加外力的速度场，使用隐式的后向欧拉法迭代 计算获得添加粘滞力的速度场；

S5，基于添加粘滞力的速度场和第一边界条件，获得满足不可压缩条 件的速度场；

其中，第一边界条件是速度场的边界条件。

S6，基于当前时间帧的前一时间帧的浓度、温度以及满足不可压缩条 件的速度场，使用半拉格朗日法回溯计算获得当前时间帧的密度场和温度 场；

S7，确定添加第二边界条件约束的密度场和温度场，获得密度网格以 及温度网格；所述密度网格是在所述密度场作用下的网格结构，所述温度 网格是在所述温度场作用下的网格结构；

其中第二边界条件是密度场以及温度场的边界条件。

S8，针对当前时间帧，基于红外光线建立的亮度网格，将从每个像素 出发的射线携带的透过率初始化为1，将每个像素携带的亮度初始化为0；

S9，在射线进入所述亮度网格并按照固定步长移动时，根据温度网格 采样获得的温度值以及密度网格采样获得的密度值，计算所述亮度网格的 自发辐射值；

可以理解，亮度网格的建立过程在后续流场的差分格式和网格结构优 化中具体介绍，针对亮度网格中每一个单元格中的粒子，从相机视口的方 向求取亮度以及透明度。

S10，计算太阳光射线进入所述亮度网格的太阳光辐射；

S11，基于所述自发辐射、太阳光辐射、消光系数以及透过率，计算射 线在当前时间帧的当前位置的透明度和亮度；

S12，重复步骤S9至S11，直到射线离开亮度网格范围；

S13，基于射线的亮度和亮度对前时间帧的当前位置视线方向的像素进 行Alpha混合渲染。

本发明实施例提供的一种基于物理过程的实时光学烟幕的渲染方法， 利用光照烟雾模型以及红外烟雾模型的渲染烟雾的过程中，根据背景辐射， 外部光源，烟雾自身透过率等因素对烟雾渲染影响，计算获得密度网格、 温度网格以及亮度网格，同时，红外模型引入了自发辐射模型，基于密度 网格、温度网格以及亮度网格计分别计算烟雾的透明度以及亮度，基于射 线的亮度和亮度对前时间帧的当前位置视线方向的像素进行混合渲染，在兼顾效率的同时提高了烟雾渲染的视觉真实性和光照准确性。

实施例二

作为本发明可选的一种实施方式，在所述S3的步骤之前，所述渲染方 法还包括：

S31，对当前时间帧的速度场的每一点求涡矢量；

S32，从所述涡矢量中确定旋度梯度矢量；

S33，将所述涡矢量和所述旋度梯度矢量的叉积，与网格尺寸、涡旋系 数的乘积，确定为所述网格结构的旋量场的涡旋力。

实施例三

作为本发明可选的一种实施方式，所述S5的步骤包括：

S51，基于添加粘滞力的速度场和第一边界条件建立压强迭代矩阵；

S52，使用并行超松弛高斯塞德尔迭代法对所述压强迭代矩阵，获得压 强场；

S53，将所述压强场叠加在所述速度场上，获得满足不可压缩条件的速 度场。

实施例四

作为本发明可选的一种实施方式，在所述S51的步骤之前，所述渲染 方法还包括：

S51a，基于无流入流出边界条件的约束表达式以及投影原理，确定运 动边界的相对速度；

S51b，对所述相对速度进行修正；

S51c，使用库伦摩擦去除滑动摩擦加速度对所述修正后的相对速度的 影响，获得第一边界条件。

实施例五

作为本发明可选的一种实施方式，所述S6的步骤包括：

基于当前时间帧的前一时间帧的密度、温度以及当前时间帧满足不可 压缩条件的速度场，在使用半拉格朗日法每次回溯计算过程中，利用中点 法对速度场对进行采样后，基于采样后速度场、前时间帧满足不可压缩条 件的速度场、当前时间帧的前一时间帧的密度以及温度，计算获得当前时 间帧的密度场和温度场。

实施例六

作为本发明可选的一种实施方式，所述S9的步骤包括：

S91，在射线进入所述亮度网格并按照固定步长每移动一次时，在所述 温度网格采样温度值，并在所述密度网格中采样密度值；

S92，基于密度值，计算消光系数和发射截面；

S93，基于所述消光系数、发射截面以及温度值，计算所述亮度网格的 自发辐射值。

实施例七

作为本发明可选的一种实施方式，所述S93的步骤包括：

在所述温度值下，对不同红外波段的单位面积的黑体在单位立体角和 单位波长间隔内的辐射功率进行积分，获得积分结果；

将所述消光系数、所述发射截面以及所述积分结果乘积确定为所述亮 度网格的自发辐射值。

实施例八

作为本发明可选的一种实施方式，所述S10的步骤包括：

S101，针对目标亮度网格，将每个到达亮度网格的射线的透过率初始 化为1，追踪每条射线的到达位置；

S102，在射线到达目标亮度网格之前，当所述射线每移动一个固定步 长时，在所述密度网格中采样密度值，并更新透过率；

S103，在射线到达目标亮度网格之后，使用已经获得的射线透过和光 源的亮度值，计算出每个目标亮度网格中亮度单元格的由散射造成的辐射 出射度，并将该辐射出射度确定为太阳光辐射。

实施例九

作为本发明可选的一种实施方式，所述S11的步骤包括：

将当前时间帧的消光系数、透过率、散射相函数以及当前时间帧的前 一时间帧的透明度输入透明度计算公式，获得当前时间帧的当前位置的透 明度；

将所述自发辐射、太阳光辐射、散射相函数以及当前时间帧的前一时 间帧亮度，计算射线在当前时间帧的当前位置的亮度。

结合实施例一到实施例九，先对本发明实施例提供的基于物理过程的 实时光学烟雾的渲染方法的原理以及实施过程作详细介绍。

(a)流体动力学的控制方程组

将烟雾的运动看作空间中的一个个大小有限不可分割的微团在运动， 这是拉格朗日法的视角，即空间是连续的，而微团是离散的，在空间中运 动的微团携带着速度、质量、温度、体积等等我们所关心的物理量。为了 方便分析和后续的计算和推导，用密度替换质量。现在只研究一个微团， 它的运动导致其物理量和空间坐标都在随着时间变化，微团开始时位于 (x₀，y₀，z₀)，其速度、温度、密度、体积大小分别为

T₀，ρ₀，V₀，经过Δt时 间后，微团运动到坐标(x₁，y₁，z₁)，其速度、温度、密度、体积大小分别为

T₁，ρ₁，V₁。任意一个微团内部的物理量都与空间和时间4个维度相关， 可以将任一时间内的物理量表示为f(x，y，z，t)。为了描述流体微团各物理量 随时间的变化，需要对物理量f计算其关于时间的全导数：

注意到右边三项与速度的关系，上式可写成：

上式的左边是全导数，也被称作物质导数，代表拉格朗日视角下的微 团物理量随时间的变化率，右边是一个新的等价微分式，这个微分式的内 涵就是欧拉表示。可以看出，在欧拉表示中，一个流体微团的物理量随时 间的变化率可以看作两个部分：1.空间上固定点物理量随时间的变化率(局 部导数)2.流体微团运动时，由空间分布不均匀导致的物理量随时间的变化 率(迁移导数)。将上式变形为：

上式就是流动中最基本的对流项的数学形式，这是流体力学中最常用 的公式，它使得拉格朗日表示和欧拉表示可以很简单地进行转换。

接下来考虑流场的质量守恒性，由于微团不可分割，质量随着微团本 身运动，所以质量的物质导数为0，但微团的密度和体积会随着时间变化， 于是有如下形式：

上式方括号内的项是单位体积微团在运动过程中体积的相对变化率， 即速度散度

的物理意义，由此，就能得到(拉格朗日)连续性方程：

带入(1-3)中，将拉格朗日表示转化为欧拉表示：

考虑流体微团表面受到的压力和粘滞力(只考虑正应力)，由牛顿第二 定律，得到动量方程，这里直接给出拉格朗日和欧拉表示下的动量方程的 数学形式：

动量方程即Navier-Stokes方程，其中，p为压强，v为粘性系数，a为 体积力的加速度(如重力，惯性力等)。N-S方程也是Cauchy动量方程排 除了偏应力后的特殊形式。

最后，设单位体积微团内能量的变化率P，单位时间流入单位体积微团的净 能量为P₀，体积力和表面力对单位体积微团的做功的功率为P₁，根据能量 守恒定律：

其中e为单位质量的内能。对于烟雾等大部分常见的流体，流入功率P₀为热功率，根据傅立叶热传导方程：

其中

是单位质量微团的加热率，k为热扩散系数。继续考虑体积力和 表面力对单位体积微团的功率，有：

结合以上两个式子，得到能量方程的拉格朗日表示：

引入对流项，得到能量方程的欧拉表示：

通过以上讨论，我们一共得到六个流体力学控制方程，组成两组具有 等价的欧拉方程组和拉格朗日方程组。不论是哪个方程组，其中的连续性 方程，动量方程，能量方程和都具有4个未知流场物理量：速度、密度、 压强。

在n维空间中，由于动量为矢量，能量和密度为标量，所以三个控制 方程联立能够得到包含2+n个方程的方程组。而未知流场变量共有3+n个 (n个方向的速度、密度、压强和内能)，所以需要引入额外的方程使得解 封闭。对于满足克拉珀龙状态方程且比热为常数的完全气体，可以引入理 想气体状态方程和量热状态方程：

使用这两个状态方程后，压强和内能建立了热力学关系，使得原始方 程组封闭。

(b)烟雾的流体力学模型

对于烟雾类型的流场计算，通常采用不可压粘性(或无粘性)流体模 型。由于流体的不可压条件导致流体微团的密度是恒定的，这里的密度不 变是指拉格朗日视角下的密度不变，这也导致了微团的体积和质量不变， 但拉格朗日微团能够运动，这导致欧拉视角下的密度仍然是变化的。用数 学上描述：

由于dρ/dt＝0，其对流项变为：

这时候欧拉视角下的密度变化仅仅是因为对流导致的质量搬运，同时， 使用拉格朗日方法中对烟幕进行模拟的时候，往往假设粒子大小和密度是 不变的(这就是为什么大部分粒子系统都不需要记录密度分布，只需要使 用一个密度常数计算外力即可)。同样的，动量方程和连续性方程变为如 下形式。

拉格朗日表示：

欧拉表示：

注意上述两个方程组是封闭的，对于拉格朗日表示，只有三个速度变 量和三个方程(密度是常数)，对于欧拉表示，有三个速度变量，一个压强 变量和四个方程。通过假设微团的密度为常数，能量方程已从方程组中完 全解耦，这表明求解不可压流体的速度场和压力场，只需用到连续性方程 和动量方程。对于涉及传热的问题，只需要在得到速度场和压强场后，使 用能量方程和量热状态方程直接得到温度场。

接下来需要计算欧拉表示下的密度分布。虽然流体微团的密度是不变 的，但由于微团在空间上是运动的，空间上任意坐标的流场密度会随流体 运动而改变，由(2-14)可以完成的密度对流的计算，从而得到空间任意位置 的密度。

(c)基于网格的数据结构

基于网格的仿真与基于粒子的仿真在架构上截然不同，网格是对空间 的离散化，将每一时刻空间的状态记录下来，网格的分辨率直接影响着仿 真的精细程度。网格可以是动态的，也可以是静态的，动态网格通常是由 三角形(四面体)组成的框架结构，每个时刻的参数会记录在顶点上，同 时顶点也会随着时间运动，这种网格通常被使用在一些欧拉-拉格朗日的混 合架构或者自适应网格算法中。本发明实施例主要使用的是静态网格，静态网格又可以分为规则网格(regular grid)和不规则网格(irregular grid)，规则 网格的每个点在对应的坐标轴对齐且均匀分布，在三维空间中常见的规则 网格有直角坐标规则网格、极坐标规则网格和圆柱坐标规则网格。与此相 反，不规则网格是不与坐标轴对齐或不均匀分布的网格。规则网格的数据 结构比较简单直观：一个三维空间的规则网格的数据结构就是一个三维数 组，只需要记录数组的大小，通过下标就能完成对数组元素的访问。本发 明实施例的烟雾模型采用直角坐标规则的矢量面心网格结构，网格在其面 心位置保存矢量数据，在其体心位置保存标量数据。单元的6个面心位置 存储该位置的速度矢量，单元中心位置存储密度，压强和温度三个标量。 网格的速度矢量的三个分量是错开保存的，这种面心网格结构在CFD领域 被称为MAC(marker-and-cell)网格，这样做有两个好处：1.能够避免“棋盘 格”式的压力分布，导致后续求解泊松方程的迭代效果不佳甚至失效。2.计 算压力项的时候存在压力与速度的耦合，错开之后能够提高差分计算的数 值精度。

由于使用了面心网格，矢量网格的比标量网格在其矢量方向上多出额 外一个单元格，假设标量网格的大小为(n_x，n_y，n_z)，则x方向的速度网格大 小为(n_x+1，n_y，n_z)，y方向的速度网格大小为(n_x，n_y+1，n_z)，z方向的速 度网格大小为(n_x，n_y，n_z+1)。所以，本发明实施例在具体实现时，对速度， 密度，压强和温度网格采用了6个独立的体素贴图。

在使用网格和粒子系统时，为了保证单位上的统一，本发明实施例使 用粒子数×10⁸/m³作为密度单位来衡量微团区域或网格区域的烟雾浓度， 使用开尔文(K)作为温度单位来描述微团区域或网格区域的平均分子热运 动，使用m/s作为速度单位来描述微团区域或网格区域的平均速度。由于这 些量都与体积无关，所以可以在粒子系统的分布模型和网格系统的分布模 型中自由切换，保证后续计算过程的准确性。

(d)基于网格的烟雾流体仿真方法

在上述描述中通过理论得到不可压粘性流动的NS方程，并且阐述了数 学上计算速度场，压力场和温度场的方法。在CFD计算中，需要将数学上 连续的偏微分方程进行某种方式的离散化，欧拉表示的NS方程采用有限差 分的方法，将方程的偏微分项用某种差分格式表示，差分格式依赖于变量 采用的离散化网格结构，最终得到用于时间推进的迭代方程。如果该迭代 方程是显式的，则可以直接由前一时刻的变量计算得到下一时刻的变量，不断推进求解。如果该迭代方程是隐式且线性(或能够线性化)的，需要 解一个线性方程组，这时候通常采用雅克比迭代法或共轭梯度下降法，经 过多次迭代得到下一时刻的状态变量。对于烟雾这类不可压粘性(非粘性) 流体，其密度场，速度场，温度场的计算采用的是显式方法，其压力场的 计算采用的是隐式方法。

1、流场的差分格式和网格结构优化

上述提到面心网格结构具有的两个好处，下面将应用差分格式进行对 其进行分析。求解有关不可压流动的问题，其关键在于由连续性方程导出 的对速度场的约束，其数学表述为：

对该方程写成其对应的中心差分格式：

在这种差分格式下，类似“棋盘格”分布的速度场将无法被更新，这 类速度场任何一个分量在其方向上以2个网格长度呈周期分布，网格上任 意一个顶点的中心差分都为0，导致其散度为0，但对于实际的物理流动， 这种“棋盘格”分布的流场是没有意义的。

同样，在求解压力项时，若对压力网格采用中心差分的方法，也会导 致计算出来的压力梯度无法修正“棋盘格”式的压力分布，造成压力梯度 为0的假象。所以，若将中心差分格式应用在不可压NS方程时，一旦出现 “棋盘格”分布的速度场这种无意义的速度压强分布，差分方程将无法改 变它。解决该问题的方法一般有两个，一种方法是使用迎风差分格式代替 中心差分格式，另一种方法是使用交错网格。迎风差分格式一般适用于已 知流场信息传递方向且流场存在非常陡的间断，但由于烟雾的弥散过程较 为缓慢且不可压流体的声速(压力传播的速度)为无穷大，NS方程对空间 变量为椭圆型方程，其空间上任意位置的流场信息来源于整个网格空间， 不能够使用迎风格式。所以本发明实施例采用交错网格来解决这个问题。

交错网格法将速度网格和压力网格的空间坐标错开排布，对速度场采 用面心网格，压力场采用以格点或者单元为中心的网格(体心网格)，其关 键在于将压力网格上相邻格点的差分用于计算速度网格所需的压力梯度， 速度场采用面心网格结构。压力场采用体心网格结构，压力点恰好位于三 个错开的速度网格单元的中心。采用交错网格后，压力网格和速度网格间 的耦合可以采用近邻差分替代中心差分，不仅仅避免了“棋盘格”分布的出现，也使得差分间隔变小，精度提高。

2、压强修正与实时解算

在上述过程中，本发明实施例介绍和分析了PCISPH，其使用隐式的迭 代法提高了时间步长，大幅加速了SPH的解算速度，在欧拉网格中，使用 隐式的方法计算压力项更加直接且有效。在本节中，将首先根据NS方程推 导出与压强相关的泊松方程的数学形式，然后根据该数学形式建立线性方 程组，最后将提出一种用于求解泊松方程的并行迭代算法，使用该算法能 够大大加速压力项的求解速度。

在不可压流场的连续性条件(式2-17)求解压力场和速度场，本发明 实施例采用的方法是压力(速度)修正法：现在有一个不满足(不可压缩 的)连续性条件的速度场u^*，称其为原速度场，现在需要寻找一个速度修 正项，使得原速度场加上一个速度修正项后，得到满足(不可压缩的)连 续性条件的修正速度。即：

速度修正项Δu由压力产生，由NS方程可知压力导致的加速度为：

上述加速度在一个时间步长Δt中产生了修正速度Δu：

因为密度得变化远小于压力得变化，所以密度ρ在局部可以看作常数， 并且时间步长Δt也是常数，可设

于是上式就变成

将上式代入到(2-3)中，得到：

在实际的迭代计算中，一般把u^*看作是在第n个时间步中，对流计算、 重力作用计算、扩散计算等步骤后产生的中间值

这些步骤使得速度 网格各点散度不为0，不满足不可压缩的连续性条件。压力修正项-Δp^*创 造了大小与

相反的散度，从而使下一时间步的速度场uⁿ⁺¹的散度 为0，满足了连续性条件，于是可以使用该速度场去对流标量场。公式(2-6) 中的u^*’为已知项，该方程只有压力修正项为未知，所以该方程为一个泊 松方程，泊松方程为椭圆型方程，这也从数学上验证了由于不可压流体的 声速为无穷大而导致其压力信息能够扩散到整个网格中。

方程右边是已知项，其差分方法采用交错网格差分，得到隐式的数值 迭代方程(注意该压力修正项中已经包含时间步常数和密度常数)：

解该线性方程组可以使用雅克比迭代法，或共轭梯度下降法(PCG) 计算近似解。上式是偏微分方程离散化得到的差分方程，将差分方程应用 在空间网格中，实际上就得到了一个线性方程组：设p为整个空间的离散 压力值展开的向量，b为由速度的散度展开的向量，该线性方程组可以写成：

Mp＝b (2-26)

方程中p和b的值很容易得到，关键在于建立矩阵M，对上述线性方程 组分析可知，系数矩阵每一行分布代表网格中某一个点的压力方程，若该 点的相邻点的速度场有效(碰撞体和计算区域的边界的速度场为无效速度 场)，则该相邻点视为有效相邻点。矩阵上的对角元素的值与有效相邻点数 目呈正比，同样的，非对角元素的值与其对应的相邻点是否有效相关，有：

上式中，n_x，n_y，n_z对应x，y，z方向上的有效相邻点数，i，j为网格三维坐 标展开的一维索引值，每个值代表网格上一个点。对于一个三维的大小为 32*32*32的压力网格，会生成一个大小为32768*32768的7对角矩阵矩阵 (6个相邻面对应的系数和1个中心系数)，由于矩阵是稀疏的对称矩阵， 能够将大小压缩为32768*4，在GPU上恰好可以使用一个32*32*32的RGBA体素纹理存储。为了将该线性方程用于共轭梯度下降法，需要令其 系数矩阵M为正定矩阵，在这里只需要改变矩阵M上每个元素的符号，并 使线性方程右边同时发生变化即可。压强矩阵变为：

依据上述方程，可以很容易建立矩阵M，并将其应用在雅克比迭代法 或共轭梯度下降法中。

对于CFL数较小的流场，一般采用雅可比迭代法或高斯-塞德尔迭代法， 雅可比迭代法收敛速度比较慢，高斯-塞德尔迭代法收敛速度较快但由于每 个分量存在耦合性，无法进行并行加速。本发明实施例采用高斯-塞德尔迭 代法的思想，按下标将压强网格分为互相耦合的两组，组内的每个元素在 单次迭代过程中不存在耦合性，所以可以对每一组分别使用并行加速算法， 并将迭代结果留给下一组使用，满足并行算法要求的同时也加快了迭代速 度。两组下标分别满足：

在实现中，本发明实施例对下标x和下标y取网格大小范围内的任意值， 以z值区分两组下标，这样下标就可以写成如下形式：

其中，l的取值范围为[0，z_max/2)，也可以代表z方向的线程数。接着， 对每组下标对应的压强向量和压强矩阵应用超松弛迭代法，其数学表述如 下：

x_n+1＝(1-sor)*x_n+sor*D^-1(b-(L+U)x_n) (2-31)

上式中的sor是松弛因子，L和U分别是压强的下三角矩阵与上三角矩 阵，它们的和是压强矩阵除去对角的部分。在三维空间中，压强矩阵的每 一行除去对角部分只有两两对称的六个分量，也就是说在迭代过程的每一 步中，耦合只发生在相邻的压强位置上，这也证明了上述的分组方法避免 了传统高斯塞德尔迭代法带来的耦合性。在本发明实施例的实践过程中， 将松弛因子设置为1.15，当迭代进行到第12次时，平均相对误差就能降低到1％以下，而一般的并行雅可比迭代法需要30次以上的迭代才能有同样 的效果。

3、流场边界条件处理

在上述介绍了流场的控制方程，控制方程描述了流场信息在时空中传 输的规律，但并不决定流场的形态。流场的形态由边界条件和控制方程共 同决定，形成了流场绕不同形状和材质的物体时表现出的不同形态。描述 物理表面的边界条件，在算法上一般将其抽象为碰撞体(Collider)，描述场 源的边界条件，在算法上一般将其抽象为发生器(Emitter)。碰撞体会在其 空间边界上保持约束条件，但它可能是运动且能够离开流场计算区域。发 射器有连续和非连续两种形态，连续的发射器和碰撞体类似，非连续的发 射器只会在初始化阶段对它占据区域的每个网格设定每个物理量的初始状 态，且不在后续的时间推进上保持它们。

对于无粘性流体，流体在碰撞体表面切向自由流动，若运动碰撞体表 面为非渗透壁，则碰撞体法向不存在质量的流入和流出，该边界条件被称 为无流入流出边界条件(No-flux Condition)，约束的数学形式如下：

(u-u_c)·n＝u_rel·n＝0 (2-32)

式中，u为流场处于边界的速度，u_c为运动边界的速度，u_rel为碰撞体 和流场的相对速度，n为碰撞处的表面法向，上述约束条件说明了流场在边 界处的流速中垂直于边界的分量为0，根据投影定理：

得到平行于运动边界的相对速度，假设u_rel是约束前的速度，其垂直于 边界方向的相对速度不为零，需要对其进行修正：

上述操作类似于一个空间滤波器，仅对碰撞边界与流场相接触部分进 行处理和更新。需注意的是，这种边界的接触在数值计算中容易被忽视， 需要对流场于碰撞边界的重合部分执行几次外延(extrapolation)迭代操作。

对于粘性流体，采用滑动边界条件处理碰撞，这时候需要额外考虑表 面摩擦力对切向速度的影响。使用库伦摩擦作为滑动摩擦模型，对速度的 切向分量减去由速度的垂直分量导致的滑动摩擦加速度对速度的影响：

式中u_t代表速度于边界表面的切向分量。其物理意义为：法向速度是 由压力造成的，若垂直表面速度与表面法向相反(-u.n＞0)，则短时间 正压力造成的法向速度“增量”(在压力作用前，法向速度假定为0)应该 呈比例转化为单位时间切向速度的减少量。在实际计算中，边界上速度的 垂直分量是由之前的迭代步骤(如压力计算)产生的，在碰撞处理步骤后 边界上的速度应只剩下切向分量。

前述就讨论了烟雾处理中常用的速度场边界条件，标量场的边界条件 一般采用Neumann条件和Dirichlet条件：

Neumann条件描述了垂直于边界方向上f变化情况，若c＝0，则在跨 越边界的法向上f不变，将无旋速度场用速度势表示，受Neumann条件约 束的速度势与受无流入流场条件(No-flux Condition)约束的无旋速度场等 价。Dirichlet条件则将标量场f约束在固定值，当c＝0时，该条件与无滑动 边界条件类似(摩擦力系数为无穷大的滑动边界条件)。对于c≠0的 Neumann条件和Dirichlet条件，通常被使用在流场的发射源上。

4、对流项计算

从物质导数公式(2-2)可知，守恒形式的流场计算中为了描述流体微 团运动产生的固定点变量随时间的变化，必须引入对流项(漂移导数)，其 偏微分方程为：

传统的欧拉方法直接在物理量网格上计算中心差分，然后和速度场相 乘并用于欧拉时间积分上，其迭代形式为：

这种方法只有在Δt很小时才能满足稳定性，在实时渲染中，Δt通常为 1/60s，对于一般网格大小和流速的流体，达不到上述方法的稳定性要求， 所以需要另寻其他方法。

对于基于粒子的流场仿真，其物理参数依赖于运动粒子的空间坐标， 研究对象是运动粒子携带的物理量，不需要考虑对流项，这种基于粒子的 方法被称为拉格朗日方法。介于拉格朗日法和欧拉法之间的半拉格朗日法 (Semi-Lagrangian Method)常常被使用在欧拉网格对流计算中，其迭代格 式如下：

f(x)ⁿ⁺¹＝f_sample(x-Δtu)ⁿ (2-38)

这里x为网格点坐标，f(x)为网格点处的变量值，两者都是离散的。 f_sample是采样后的网格函数，x-Δtu为回溯项。其物理意义为：流场流 动造成的随时间变化是由流场的空间分布不均匀造成的，回溯的目的在于 寻找由于流动影响当前位置的流场变量的空间坐标。为了提高回溯的精确 度，可以采用中点法或BFECC方法，中点法对速度进行了一次额外的采样：

u＝u_sample(x-0.5Δtu) (2-39)

由于半拉格朗日法不断利用前一时刻的结果进行插值，这种操作包含 了平滑效应，可以防止在迭代过程中数值的震荡和发散，所以它对时间步 长的要求没有其他方法严格(在数值计算上具备无条件稳定性，这是该方 法的主要优势)。但也造成了流场细节信息在多次迭代后被“抹平”，即便 使用了中点法和BFECC提高回溯的精度，这种平滑效应仍然存在。这种平 滑效应是由采样算法中使用的线性插值带来的，在网格精度不高的情况下会导致场的扩散(平滑效应)，这种扩散被称为数值扩散，导致其细节特征 消失。对于速度场，这种扩散会带来额外的黏滞性，导致流场失去涡流。 有多种方法解决这个问题，采用粒子和网格混合的方法就是其中之一，这 种混合方法使用基于粒子(有限元)的方法解算对流项，并使用基于网格 (有限差分)的方法解算剩余项目。另一种完全基于网格的方法是使用立 方插值(Cubic Interpolation)代替线性插值。基本思想是通过一系列 Catmull-Rom插值算法替代三次线性插值。但由于立方插值破坏了线性插值 的单调性，需要对其结果进行额外的截断处理。

另外，在半拉格朗日法的回溯过程中，需要对回溯得到的坐标进行一 次额外的边界检测。

除了使用半拉格朗日法外，另一种典型的对流项计算方法将迎风格式 作为对流方程的差分格式，这种方法被称为迎风法。将式(2-14)写成一维 形式：

迎风格式的思想在于根据流场信息的流动方向(这里是速度的方向) 选择合适的单侧差分，将上式写成迎风差分格式：

根据速度的方向得到相应的迎风格式后，将其使用在欧拉法中，得到 迎风法的迭代公式：

迎风法采用基于有限差分的欧拉法进行迭代，需要考虑其稳定性，稳 定性要求满足CFL条件(Courant-Friedrichs-Lewy Condition)：

其中uΔt/Δx被称为库朗数(Courant Number)。所以迎风法的稳定性条 件为库朗数小于1。迎风法相较半拉格朗日法在同样时间步长具备更高的精 度，但它的时间步长选择受到速度和网格密度的限制，采用自适应网格和 自适应步长对其进行优化可以取得很好的效果，如果需要更高的精度，可 以采用二阶迎风格式代替一阶迎风格式。迎风方法一般用于高速流场的仿 真中，由于本发明实施例仿真的烟雾是低速流场，所以在流场网格中使用 基于半拉格朗日的对流方法。

5、外力项计算

体积力和浮力等外力的计算相对简单，一般都是整合在外力项的步骤 中完成。由于体积力是覆盖整个网格范围的非接触力，计算重力惯性力等 体积力，只需要在每个时间推进步上对每个网格添加体积力造成的加速度， 然后再执行一次边界检查即可：

u^*＝u+gΔt (2-44)

由于爆炸或者燃烧造成的的密度分布不均匀，压力修正法在计算时需 要引入不均匀的密度分布，表示压力修正项的泊松方程(2-23)需要修改为：

采用有限差分法离散化后，上式是个复杂的线性方程组，对于密度分 布极度不均匀的流场，解这个方程组可能是必要的，但带来的计算成本也 较为巨大。为了简化该问题，将浮力近似为与温度相关的体积力，有经验 公式：

f_b＝-αρy+β(T-T_amb)y (2-46)

式中，α和β分别表示密度和温度对浮力的影响系数，T_amb为环境温度， 可以由计算平均温度得到，也可以设定一个固定值。还有其他浮力计算的 经验公式：

6、降低流场的涡量消失现象

半拉格朗日法带来了一种无条件稳定的网格对流方法，但也带来了不 可避免的数值粘性耗散，该数值粘性耗散是由于半拉格朗日方法的采样过 程会将原本剧烈的流场波动平滑化，使得流场的熵增大，剧烈波动的区域 趋于平缓。这种耗散对高雷诺数(低粘性)的流动的影响尤为明显，特别 是烟雾，失去涡旋特征的烟雾会显得单调，让观看者感到细节丢失。一般 有两种办法可以加强烟雾的涡旋效果，这里将介绍第一种方法，第二种方法将在混合粒子和网格的章节进行讨论。

在这里介绍一种基于网格的涡量恢复方法，该方法被称为涡量约束法(Vorticity confinement)。其思想是将流动的旋量从对流采样前的网格中提 取出来，在对流过后将旋量按一定比例转化为涡旋力并作为外力项的一部 分施加在网格上，达成涡量的恢复。从流场上提取旋量非常简单，只需要 对速度网格每一点求其涡矢量即可：

然后从涡矢量中获得其大小的梯度矢量：

旋度梯度矢量指向旋度增大的方向(即旋转中心的方向)，所以涡矢量 和旋度梯度矢量的叉积就是当前位置涡旋力的方向。该涡旋力需要乘上网 格尺寸h才能保持对不同的网格大小的不变性，最后再用涡旋系数∈控制涡 旋放大的倍率，即：

f_v＝∈h(N×ω) (2-49)

涡量约束法本质是对旋量的提取放大后的重引入，旋量本身不会对流， 只会不断生成和消失。

7、粘滞力的解算方法

由推导的动量方程(2-15)中分离出粘滞力加速度，该部分在NS方程 中也被称为扩散项：

其中，μ为粘性系数，使用欧拉法可以得到速度场迭代方程：

粘滞力项可以使用显式方法或隐式方法求解，显式方法一般用于粘滞 力较小的情况，隐式方法一般用于粘滞力较大的情况。下面将对两种方法 进行讨论。

如果直接在u_n处计算粘滞力，该迭代法成为显式方法，又被称为前向 欧拉法。显式方法在每个时间步对速度场可以直接计算出其改变量。当然， 显式方法的缺陷是受到的步长限制比较严格，若拉普拉斯算子的差分格式 为中心差分：

则它的稳定性条件如下：

这种计算扩散项的方法被称为前向欧拉法(Forward Euler)，常用于粘 性系数较低的情况。当粘性系数较大时，其稳定性条件就成了一个问题， 步长精度不够的情况下，会造成不稳定的扩散解。

不稳定是显式欧拉法的常见问题，通常的解决手段是找到合适的隐式 方法替代显式方法。将方程(2-25)右边第二项移动到左边，并且在n+1 的网格处计算拉普拉斯算子的中心差分，得到后向欧拉迭代方程：

uⁿ⁺¹-ΔtvL(uⁿ⁺¹)＝uⁿ (2-54)

其中，L(·)为拉普拉斯算子的中心差分表示，该方程为一个线性方程组， 和压力修正项计算一样，求解该线性方程组可以使用共轭梯度下降法(当

较大时)，或雅克比迭代法、高斯-赛德尔迭代法(当

较小时)。

体绘制通常用在查看三维物体内部的结构，而烟雾呈现出来的亮度和 透过率就是将在相机视线方向上不同深度的烟雾“切片”通过数学手段叠 加在一起形成的。本发明实施例采用了体绘制种最常用的光线投射算法， 其过程非常简单，首先需要一个表示烟雾占据区域的长方体，当光线命中 长方体前表面时开始体绘制：光线不断地在烟雾体积内步进，不断累计亮 度和透过率，当光线离开长方体后表面，体绘制结束，其过程中累计的亮度和透过率就是烟雾的亮度透过率。

图2是基于光线投射的体绘制渲染结构示意图。首先，在现代GPU渲 染管线中，物体的透明度通道会被使用在透明度混合中，所以体绘制烟雾 时不需要考虑背景的颜色，只需要保证烟雾的渲染顺序在背景之后，并根 据体数据(譬如密度网格)计算出消光系数和前向散射系数，通过在路径 上的积分，计算出透过率，该透过率包含了烟雾的遮挡和前向散射效果。 其次，输入像素着色器的顶点仅仅是代表网格占据空间的长方体的局部坐 标，通过将相机的世界坐标变换到网格空间中，并将像素位置的网格坐标 减去相机坐标，就能得到从该像素出发沿着视线方向的矢量。最后，数值 积分沿着该像素视线方向进行，以网格间隔为步长采样体纹理，当采样点 离开网格区域后，积分停止，将积分计算的结果(亮度和透过率)作为该 像素的结果，从像素着色器输出到渲染管线的颜色混合阶段，与背景进行 混合。

(e)光学烟雾的渲染

1、可见光烟雾的光照渲染

渲染可见光烟雾需要主要考虑的因素是透过率和对光源的散射。在本 章的开头介绍过，烟雾粒子对光子的相互作用概率决定了其透过率，从视 线方向穿透烟雾体的总透过率就是最终被渲染出来的Alpha通道。光子在 单位路径上与烟雾粒子的相互作用的平均概率也被称为消光系数，粒子的 消光系数与烟雾的消光截面和密度相关，有：

σ_t＝C_extρ (3-1)

上式中σ_t为消光系数，C_ext为消光截面。当光子与粒子进行相互作用时， 有一部分被吸收，一部分被散射，散射的概率被称为反射率Ω。密度大且厚 的烟雾的透过率一般比较低，所以决定烟雾颜色的是反射率。燃烧有机物 产生的浓烟含碳量高，反射率低，所以呈现黑色，云层和水蒸气包含的主 要是水汽，反射率高，所以呈现亮白色。

对于不同的烟雾成分，其散射特效也不同。一般来说，散射主要集中 在光束传播的方向上。粒子的散射概率在不同角度上的分布在数学上用相 函数(phase function)来描述，本发明实施例采用Henyey-Greenstein函数 作为烟雾的散射相函数，该函数的特点是比较平滑：

上式中，g是一个无量纲的控制散射相函数分布的常数。当g＜1时， 散射是各项异性的，g越大，分布越集中，前向散射占据的比率越大。当g＝1 时，散射方向完全集中在前向。由于背景的前向散射是发生在体绘制的积 分路径上的，所以前向散射在实际应用中可以当作透过率来处理。

Henyey-Greenstein相函数在4π空间内的积分是归一化的，即：

可见光烟雾发出的光子大多来自外界光源，比如太阳。所以渲染可见 光烟雾需要进行两个渲染通道。第一个通道对烟雾的每个亮度网格单元， 在它到光源的路径上进行路径积分，亮度网格的大小不一定要与密度网格 的大小相同，一般会使用较小的亮度网格以减少计算的时间复杂度，本发 明实施例将该通道称为亮度网格通道。亮度网格通道的具体计算过程和体 绘制类似，区别在于体绘制是针对相机视口的每个像素积分它的亮度和透过率，该通道是对网格的每个单元积分它的透过率。积分过程也是只考虑 前向散射，其具体计算过程如下：

第一步，将每个到达目标亮度网格的射线的透过率被初始化为1，然后 开始追踪每条射线。

第二步，在射线到达目标亮度网格之前，每移动一个固定步长，就在 密度网格中采样并更新透过率，该透过率包括前向散射的部分，假设前向 散射的夹角很小，前向散射率为：

亮度网格的透过 率更新公式如下：

在上式中，

代表的是前向散射贡献的透过率。

第三步，在射线到达目标亮度网格之后，使用已经获得的射线透过和 光源的亮度值，计算出每个亮度单元格的由散射造成的辐射出射度M_s，其 数学表述如下：

需要注意的是，该辐射出射度M_s是包括整个4π球面空间的。上式中， 每个单元网格除了自己吸收一部分光子，其余的光子都以散射的形式向4π 立体角空间发射，这些散射的光子最后部分进入人眼和相机，形成烟雾外 形的亮度。第二个渲染通道将使用体绘制算法对烟雾进行绘制，在积分路 径上将分别对密度网格和亮度网格进行采样，采样到的密度被用来计算透 明度，采样到的散射辐出度被用来计算亮度，其路径积分过程中的透明度和亮度更新公式如下：

上式中的ψ_θ，Δθ是以θ为散射角，以Δθ为(像素对网格点的)散射张角， 积分计算相函数得到的散射率。为了保证计算的实时性，采用预计算的方 法获得一张表示不同方向和张角的散射率的贴图，该帖图通过对不同散射 角度和散射张角计算散射相函数的数值积分值得到，在实时渲染中通过采 样这张二维贴图得到所需的散射率。最终，累加的亮度和累乘的透明度将 作为像素着色器的输出，被用于Alpha混合阶段。

2、红外烟雾的光照与自发辐射渲染

红外烟雾特别是军事上用于遮挡军事目标而使用的红外烟幕，发烟器 或者其他烟雾发射源都会具有比环境温度高的初始温度，烟雾的红外自发 辐射不可以忽略。所以，相较于可见光烟雾，红外烟雾需要额外考虑一个 自发辐射项。烟雾内部的自发辐射一部分会被外围的烟雾吸收和散射，另 外一部分会透射出来并进入相机。烟雾内部的红外散射同样遵循着和可见 光一样的原理，在这里，本发明实施例同样只考虑在视线积分路径上的前向散射。体绘制过程中，透明度和亮度的路径积分公式为：

上式与可见光的区别在于，亮度更新公式中多了一项L_e，该项代表路 径积分经过的位置处的烟雾粒子团自发辐射的亮度，该亮度与该粒子团的 等效辐射面积有关，而粒子团的等效辐射面积又与单个粒子的大小，粒子 数密度，粒子团的大小有关。考虑将一个密度单元格作为粒子团，单元格 中随机分布着具有相同粒径的大量粒子，粒子被抽象为球形，其辐射截面 在任意方向上都是面积相同的圆。粒子群在视线方向的正投影是一个由许多圆球相互重叠而形成的阴影，该阴影的面积可以看作单元格内粒子群的 等效辐射面积，也就是说，可以把粒子群的辐射等效为一个位于网格中心， 法向方向与视线方向相同，面积与粒子群等效面积相等的正方形黑体的辐 射。设粒子的半径为r，粒子的辐射截面为：

A_p＝2πr² (3-8)

当视线从单元格某个方向上观察时，粒子的辐射截面会发生重叠，重 叠是一个随机的事件，在这里本发明实施例使用条件概率的思想，假设当 前有n个粒子，网格内粒子的等效辐射截面为S_n，当加入一个新的粒子时， 网格截面上任意一点落在新粒子截面上的概率为A_p/D²，D²是网格截面的 大小。原等效辐射截面S_n在发生遮挡后的数学期望为：

加入新的粒子并考虑该粒子的遮蔽效应后，第n+1个粒子的等效辐射 面积为：

将该数列展开，就得到如下形式：

很显然，上式是一个等差数列的前N项和，求和之后有：

假设粒子的截面远小于网格大小，这时候就有

上式可以进行泰 勒展开，得到等效辐射面积的计算公式：

接下来，由于本发明实施例在流体计算部分使用的是粒子密度作为密 度网格的参量，所以这里需要将粒子数参量转化为粒子密度，根据N＝ ρD³/m_p，这里m_p为单个粒子的质量。转化后的等效辐射面积函数S(ρ)为：

单位网格面积的黑体等效辐射面积又被称为发射截面，由上式可以得 到发射截面R_e(ρ)的数学表述为：

接着假设粒子的等效发射截面是一个发射率与粒子相同，大小为σ的朗 伯体，根据黑体辐射公式，可以得到网格在视线方向的辐射出射度M_r：

上式中的M_b(λ，T)是普朗克公式，它的物理意义为：在一定温度下，单 位面积的黑体在单位立体角和单位波长间隔内的辐射功率。一般的红外相 机都是针对特定波长段的红外光，本发明实施例主要考虑3-5um和8-12um 波长的下的黑体辐射公式，对普朗克公式在对应的波长范围内进行数值积 分就能得到。

与散射相函数一样，该黑体在不同波段下的辐射曲线可以存储在格式 为单通道、32bit浮点数精度的DDS纹理中，在体绘制的渲染过程中，通过 采样特定的纹理贴图，获得不同波段范围的黑体辐射出射度函数

最后，考虑发射截面和发射率和朗伯体的2π发射空间，得到网格上由自发 辐射产生的亮度L_r：

参考图3所示，红外烟雾的整体渲染过程如下：

步骤一：将从每个像素出发的射线携带的透过率初始化为1，将其携带 的亮度初始化为0。

步骤二：在射线到离开网格范围之前，每移动一个固定步长，就在密 度网格和温度网格中采样，密度值用来计算消光截面和发射截面，温度值 用来采样黑体辐射贴图。然后根据(3-17)计算出网格的自发辐射亮度。如果 考虑太阳光辐射的影响，还需要额外采样亮度网格(亮度网格的计算过程 与可见光相同)，并计算太阳光的散射亮度。

步骤三：将自发辐射、太阳光辐射、消光系数以及上个积分位置的透 过率综合考虑，应用公式(3-7)，计算射线在当前位置的透过率和亮度。

步骤四：不断地重复步骤2和步骤3，直到射线离开网格范围。离开网 格范围后，将射线的亮度和透过率作为该像素的输出，用于Alpha混合阶 段。

3、本发明实施例提供的基于物理过程的实时光学烟雾的渲染过程的实 现。

首先将处于等待初始化状态(刚刚加入场景或者刚刚完成一帧计算) 的烟雾进入第一个计算通道即初始化通道，该通道的计算着色器根据烟雾 发生器(Emitter)的位置、浓度和温度，以及其他碰撞器(Collider)的边 界条件初始化前一时刻的网格，每个发生器或碰撞器只负责初始化其覆盖 的网格区域，其余区域保持不变，仿真流程如图3所示：

步骤一：初始化后的网格将进入第二个计算通道——速度场对流通道， 在这里使用的对流方法是半拉格朗日法(Semi-Lagrange Method)，采样前 一时刻的速度场，然后回溯计算即对流项计算，得到当前时刻的速度场。

步骤二：完成该步骤的烟雾会接着进入外力计算通道，在该通道内， 涡旋力从旋量场中被提取出来，与烟雾的浮力、重力一起作为合外力施加 在速度网格的每个单元上。获得外力加速度的烟雾速度网格，即先通过外 力项计算，更新速度场，使得原速度场上加上外力，并将速度场存储在每 个网格中；

步骤三：继续进入粘滞力计算通道，使用隐式的后向欧拉法进行2-5 步的迭代(因为烟雾大多是低粘性流体，粘性系数很低)，就能得到加入粘 性后的速度场；该粘性计算方法也会将外力加速度带来的影响传递到周边 网格；

步骤四：计算完成所有外力和粘滞力后，这时的速度场不满足不可压 缩条件，所以烟雾数据将进入第4个通道——“压强修正通道”，这个通道 主要分三步进行：1.从当前速度场(外力+粘性)和边界条件建立压强迭代 矩阵，2.使用并行超松弛高斯塞德尔迭代法解算压强场(见压强修正与实时 解算)，3.将压强场梯度作用在速度网格上，得到满足不可压缩条件的速度 场，得到压强修正的速度场；

步骤五：标量场对流通道开始处理烟雾数据，采样(满足不可压缩条 件)的压强修正的速度场，上一帧的标量场，然后使用对流项计算回溯得 到当前时刻的标量场；同时回溯步骤使用了中点法来提高回溯精度，同样 的，在有碰撞体的情况下，需要额外考虑碰撞边界条件，得到温度网格以 及密度网格。

上述过程就是基于网格的烟雾流场模型的计算过程，但在进行体绘制 渲染之前，还需要根据太阳光的方向和强度生成亮度网格，生成亮度网格 的方法在红外烟雾的光照与自发辐射渲染中由详细介绍，在这里就不赘述 了。获得当前的亮度网格，密度网格和温度网格后，这些资源将绑定在像 素着色器上，在图形通道进行体绘制渲染。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明， 不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域 的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简 单推演或替换，都应当视为属于本发明的保护范围。

Claims (9)

1.一种基于物理过程的实时光学烟幕的渲染方法，其特征在于，包括：

S1，获取初始时间帧前一时间帧网格结构中烟雾发生器的位置、浓度和温度以及碰撞器的边界条件，初始化当前时间帧的网格结构；

其中，所述网格结构是直角坐标规则的矢量面心网格，所述网格结构中每个单元网格的6个面心位置存储该位置的速度，单元网格中心位置存储密度，压强和温度；

S2，基于已获取的当前时间帧的前一时间帧的速度场，使用半拉格朗日法回溯计算获得当前时间帧的速度场；

S3，基于所述网格结构的旋量场的涡旋力、烟雾的浮力、重力以及当前时间帧的速度场，计算得到当前时间帧添加外力的速度场；

S4，基于当前时间帧添加外力的速度场，使用隐式的后向欧拉法迭代计算获得添加粘滞力的速度场；

S5，基于添加粘滞力的速度场和第一边界条件，获得满足不可压缩条件的速度场；

S6，基于当前时间帧的前一时间帧的浓度、温度以及满足不可压缩条件的速度场，使用半拉格朗日法回溯计算获得当前时间帧的密度场和温度场；

S7，确定添加第二边界条件约束的密度场和温度场，获得密度网格以及温度网格；所述密度网格是在所述密度场作用下的网格结构，所述温度网格是在所述温度场作用下的网格结构；

S8，针对当前时间帧，基于红外光线建立的亮度网格，将从每个像素出发的射线携带的透过率初始化为1，将每个像素携带的亮度初始化为0；

S9，在射线进入所述亮度网格并按照固定步长移动时，根据温度网格采样获得的温度值以及密度网格采样获得的密度值，计算所述亮度网格的自发辐射值；

S10，计算太阳光射线进入所述亮度网格的太阳光辐射；

S11，基于所述自发辐射、太阳光辐射、消光系数以及透过率，计算射线在当前时间帧的当前位置的透明度和亮度；

S12，重复步骤S9至S11，直到射线离开亮度网格范围；

S13，基于射线的亮度和透过率对前时间帧的当前位置视线方向的像素进行Alpha混合渲染。

2.根据权利要求1所述的渲染方法，其特征在于，在所述S3的步骤之前，所述渲染方法还包括：

S31，对当前时间帧的速度场的每一点求涡矢量；

S32，从所述涡矢量中确定旋度梯度矢量；

S33，将所述涡矢量和所述旋度梯度矢量的叉积，与网格尺寸、涡旋系数的乘积，确定为所述网格结构的旋量场的涡旋力。

3.根据权利要求1所述的渲染方法，其特征在于，所述S5的步骤包括：

S51，基于添加粘滞力的速度场和第一边界条件建立压强迭代矩阵；

S52，使用并行超松弛高斯塞德尔迭代法对所述压强迭代矩阵，获得压强场；

S53，将所述压强场叠加在所述速度场上，获得满足不可压缩条件的速度场。

4.根据权利要求3所述的渲染方法，其特征在于，在所述S51的步骤之前，所述渲染方法还包括：

S51a，基于无流入流出边界条件的约束表达式以及投影原理，确定运动边界的相对速度；

S51b，对所述相对速度进行修正；

S51c，使用库伦摩擦去除滑动摩擦加速度对所述修正后的相对速度的影响，获得第一边界条件。

5.根据权利要求1所述的渲染方法，其特征在于，所述S6的步骤包括：

基于当前时间帧的前一时间帧的密度、温度以及当前时间帧满足不可压缩条件的速度场，在使用半拉格朗日法每次回溯计算过程中，利用中点法对速度场对进行采样后，基于采样后速度场、前时间帧满足不可压缩条件的速度场、当前时间帧的前一时间帧的密度以及温度，计算获得当前时间帧的密度场和温度场。

6.根据权利要求1所述的渲染方法，其特征在于，所述S9的步骤包括：

S91，在射线进入所述亮度网格并按照固定步长每移动一次时，在所述温度网格采样温度值，并在所述密度网格中采样密度值；

S92，基于密度值，计算消光系数和发射截面；

S93，基于所述消光系数、发射截面以及温度值，计算所述亮度网格的自发辐射值。

7.根据权利要求6所述的渲染方法，其特征在于，所述S93的步骤包括：

在所述温度值下，对不同红外波段的单位面积的黑体在单位立体角和单位波长间隔内的辐射功率进行积分，获得积分结果；

将所述消光系数、所述发射截面以及所述积分结果乘积确定为所述亮度网格的自发辐射值。

8.根据权利要求1所述的渲染方法，其特征在于，所述S10的步骤包括：

S101，针对目标亮度网格，将每个到达亮度网格的射线的透过率初始化为1，追踪每条射线的到达位置；

S102，在射线到达目标亮度网格之前，当所述射线每移动一个固定步长时，在所述密度网格中采样密度值，并更新透过率；

S103，在射线到达目标亮度网格之后，使用已经获得的射线透过和光源的亮度值，计算出每个目标亮度网格中亮度单元格的由散射造成的辐射出射度，并将该辐射出射度确定为太阳光辐射。

9.根据权利要求1所述的渲染方法，其特征在于，所述S11的步骤包括：

将当前时间帧的消光系数、透过率、散射相函数以及当前时间帧的前一时间帧的透明度输入透明度计算公式，获得当前时间帧的当前位置的透明度；

将所述自发辐射、太阳光辐射、散射相函数以及当前时间帧的前一时间帧亮度，计算射线在当前时间帧的当前位置的亮度。

Images