CN116933674A - 基于lbm的快速流体模拟方法 - Google Patents

Abstract

基于LBM的快速流体模拟方法属于计算机图形学领域，具体是模拟了龙卷风的动画效果。首先针对目前大多数龙卷风模拟方法仅局限于粉尘颗粒本身，在以往龙卷风物理模型的基础上，利用双流体系统思想构建空气与粉尘颗粒双流体龙卷风模型；然后针对标准LBM无法求解双流体问题，基于二元流体碰撞思想，构建两相格子玻尔兹曼方程，使用D3Q19模型，并根据边界单元在网格模型中的空间位置，将边界单元分为两类，分别对两类边界单元的物理量进行求解；最后，设计和实现龙卷风模拟系统，设计总体框架和各个功能模块，实现模拟。

Description

基于LBM的快速流体模拟方法

技术领域

本发明基于LBM格子玻尔兹曼方法，研究了一种针对快速流动的流体模拟的方法，并最终实现对龙卷风的模拟。首先针对目前大多数龙卷风模拟方法仅局限于粉尘颗粒本身，在以往龙卷风物理模型的基础上，利用双流体系统思想构建空气与粉尘颗粒双流体龙卷风模型；然后针对标准LBM无法求解双流体问题，基于二元流体碰撞思想，构建两相格子玻尔兹曼方程，使用D3Q19模型，并根据边界单元在网格模型中的空间位置，将边界单元分为两类，分别对两类边界单元的物理量进行求解；最后，设计和实现龙卷风模拟系统，设计总体框架和各个功能模块，实现模拟。本发明属于计算机图形学领域，具体涉及LBM格子玻尔兹曼方法，流体模拟等技术。

背景技术

对快速流动且不可预测的洪水、火焰火灾等快速流体进行研究和模拟，能够使研究者和救援队伍迅速掌握灾情现场的情况，精准开展救援，能够最大程度减少自然灾难带来的损害。另一方面，随着科技的发展，我们可以看到通过计算机模拟动画在虚拟仿真、电影特效、广告制作、游戏场景设计等各个方面都有广泛的应用。例如电影中的大洪水，游戏中的瀑布、火焰场景等，使得电影和游戏的效果更加逼真，用户体验感更好。

流体中的烟雾、水体、沙尘暴、龙卷风等自然场景一直是我们生活生存必不可少的一部分，龙卷风更是如此。

在计算机图形学领域，像龙卷风等高速旋转流动的流体的模拟一直是一个巨大的挑战。随着虚拟现实的爆炸式发展，流体模拟技术被广泛应用于虚拟现实的内容开发，像2014年好莱坞大片《星际穿越》能在票房以及口碑上名利双收，除了穿越式的剧情以及脑洞大开的天文学元素外，计算机图像学模拟的特效比如龙卷风、海洋风暴以及黑洞效果等着实为影片加分不少，栩栩如生的虚拟海洋场景、逼真的龙卷风场景的加入使我们亲身感受到环境的变化，提升我们对环境的爱护意识。

龙卷风作为流体的一种，其内部结构和形状都是不规则的，所以对龙卷风的建模或者控制非常具有挑战性。首先，龙卷风本身是由无数的小粒子构成，且这些粒子在随机运动的过程中会受到气压、温度、重力等各种外界因素影响。而龙卷风是由于粒子受到漩涡力的作用绕着中心轴做高速运动形成，由于粒子数量众多，大量粒子的数值计算对计算机的开支特别大，龙卷风的实时性是一个重要挑战。此外，龙卷风在运动的过程中，大量粉尘颗粒粒子在做随机无规则运动的时候，如何保持相对的龙卷风的形状极具挑战，而且其形状影响龙卷风的真实性。最后，龙卷风的形态非常复杂，我们很难用传统的单相流模型模拟出来，但是多相流的模型的引入又会增加流体之间的耦合，增加了计算的开支，影响模拟的实时性。但是不考虑流体之间的耦合效果，将会导致龙卷风的很多细节的丢失，影响其真实性。综上所述，龙卷风的模拟是一项非常有挑战性的课题。

科学家们对于自然现象的模拟开展了广泛的研究，比如具有丰富细节的烟雾模拟，大尺度的海洋模拟、水面流动模拟，具有真实感的各类火焰模拟等。在追求物理准确性的同时，更要追求高效率，要用尽可能少的资源，绘制出更好的视觉效果，这也是科学研究所追求的目标。

在传统的流体计算方法中，一般是先基于N-S方程将宏观控制方程离散，然后使用某种数值方法来求解离散方程，最终得到宏观的各个物理量。传统的流体模拟方法已经发展得比较成熟，模拟效果也比较好，但对于快速流动的流体，由于其流体粒子运动速度较快，然而基于N-S方程的计算模拟速度较慢，从而在计算机模拟速度与真实的流体运动状态之间产生偏差，影响最终效果。LBM是近十年来新发展起来的一种流体模拟方法，它用简单的微观动力学模型来模拟复杂的宏观现象，具有高并行性、演化过程简单、微观碰撞容易实现等优点。作为解决流体力学数值模拟的方法之一，LBM所提出的建模手段和建模思想，为流体模拟开创了新的领域。

对于任何一个流体模型，都是需要大量的计算维持的，当实际情况越来越复杂时，计算量更成指数型增长，因此流体模拟长期以来都受着计算机计算性能和资源的限制。但随着高性能计算机的发展，具有相同计算能力的GPU成本远小于CPU，因此可以将强大的图形硬件用于处理流体模拟的复杂问题，由于GPU非常适合处理那些能够表示为并行计算的问题，因此GPU也从一开始的图像视觉处理逐渐发展到越来越广泛的科学研究领域。充分利用GPU的优势，会给流体模拟的研究带来巨大的便利，研究GPU技术在LBM流体模拟中有着重要的意义。

另外，目前使用LBM对于快速流动的流体，比如龙卷风的模拟依然比较少，有必要基于LBM对快速流体的模拟进行探索和研究，基于这一目标，本发明提出了一种以模拟龙卷风为例，基于LBM模拟快速流体的方法。

发明内容

针对快速流体的模拟，本发明选择龙卷风作为模拟对象。目前对龙卷风的研究主要集中在实验室观测、数学建模和数值模拟等方面。这些方法都不能生成逼真的可视龙卷风。本发明的目的是利用计算机图形学技术，在真实物理模型的基础上探究一个模拟真实龙卷风流动的算法。主要包括两个模块：一是设计一个流体物理模型，能够较真实的描述龙卷风的流动。龙卷风实际上是高速旋转的空气，与周围的粉尘颗粒相互作用产生的效果。目前比较在模拟龙卷风效果的研究中，采用的比较多的物理模型，仍然是以N-S方程为主的单流体模型，不能很好的模拟逼真的龙卷风效果。故本发明提出了一种以N-S方程为主的双流体模型，以高速运动的空气作为第一种流体，周边的粉尘颗粒作为第二种流体，以此建立龙卷风的物理模型。二是基于LBM的思想求解，由于需要求解的是双流体模型，本发明在标准LBM的基础上做了拓展，构建了两相LBM模型。该方法的提出，主要解决了LBM在模拟快速流动的流体(如龙卷风)方面的技术问题。本发明过程如附图1所示。可分为以下三个步骤：基于空气与粉尘颗粒的双流体龙卷风模型构建、两相流LBM方程构建、系统的设计与实现。

(1)基于空气与粉尘颗粒的双流体龙卷风模型构建

本发明研究对象为龙卷风，为了有效模拟龙卷风旋转运动的物理状态，首先构建一种基于双流体系统的龙卷风物理模型，该模型结构清晰，遵循物理规则，能有效保证流体流动真实性。针对龙卷风的两大构成—空气、粉尘颗粒，构建了双流体N-S模型，模型表达了空气与粉尘颗粒之间的相互作用力，以及除此之外粉尘颗粒受到的外力约束。

(2)两相流LBM方程构建

对于龙卷风物理模型，使用两相LBM求解，对于碰撞模型，在标准LBM碰撞项的基础上，拓展为自碰撞项和互碰撞项。对于外力项，则以物理模型中粉尘颗粒这一层流体为主体，引入这一层流体所受到的外力作为LBM方程的外力项。并按照格子在空间中的位置，将边界格子分成两类，一类是空间顶层和底层的边界格子，一类是空间四个侧面的边界格子，对两类不同的边界格子作出不同的边界约束。

(3)系统的设计与实现

根据前面建立的基于边界约束的龙卷风模型以及LBM求解过程，设计并实现龙卷风模拟系统，同时对该系统进行实验，并对实验结果进行详细分析。整个系统分为用户层、数据处理层和数据存储层。用户层主要负责对用户操作的接收和传递，数据处理层主要负责整体计算流程，数据存储层主要负责存储计算流程中的关键物理量，便于数据处理层开始下一层循环，以提高整体计算效率。

本发明与现有技术相比，具有以下明显的优势和有益效果：

本发明提出了一种基于LBM的快速流体模拟方法，以高速旋转的龙卷风为例，扩展了单流体模型来模拟龙卷风的流动及其与周围碎片的相互作用。然后，构建了两相LBM模型来求解龙卷风模型。这个方法相比于图形学领域的传统的进行界面捕捉并求解N-S方程的方法，具有全新的特点和优势，不需要设计复杂的界面捕捉算法。其次，本发明提出的方法的所有步骤都很容易并行化，并在实验中通过在GPU上实现了相关系统，成功绘制出了龙卷风效果。本发明提出的方法，有助于了解龙卷风的气流运动，也有助于了解龙卷风的形成环境，

附图说明：

图1基于LBM的快速流体模拟技术流程图

图2系统整体框架图

图3系统各功能模块图

具体实施方式

根据上述描述，以下是一个具体的实施流程，但本专利所保护的范围并不限于该实施流程。

步骤1：基于空气与粉尘颗粒的双流体龙卷风模型构建

步骤1.1：在传统模型的基础上，构建一个空气与粉尘颗粒的双流体模型

本发明使用的龙卷风模型，基于N-S方程改造得来，具体如下：

其中为哈密顿算符，其计算公式为：/>对任意一个量执行哈密顿算符，等于对该量的所有分量求偏导之和。

把空气视为第一个流体，空气流动速度为u₁,式(1)为空气这一层流体的连续性方程；把粉尘尘埃等颗粒物视为第二个流体，流动速度为u₂，式(2)为粉尘颗粒物这一层流体的连续性方程。对于式(3)，p₁为空气这一层流体的压力场。对于式(4)，p₂为粉尘颗粒这一层流体的压力场；在式(3)和式(4)中，f₁是空气与粉尘颗粒之间的相互作用力，f₂是粉尘颗粒所受到的涡度约束。

对于式(3)、式(4)和式(5)中这一项，为空气与粉尘颗粒间相互作用力，ρ₂为粉尘颗粒这一层流体的密度,m₂为粉尘颗粒的质量，v₂为粉尘颗粒这一层流体的运动粘度。

对于式(4)和式(6)中的为粉尘颗粒在运动中受到的涡度约束，v_x,v_y,v_z为粉尘颗粒速度矢量的三个分量，r为粉尘颗粒运动时绕所在平面做旋转运动的半径。

步骤2：两相流LBM方程构建

步骤2.1：确定碰撞模型

本发明使用的格子模型为D3Q19模型(3维空间，19个离散速度)，以下为使用的两相LBM模型：

f_i(r+e_iΔt，t+Δt)-f_i(r，t)＝Ω_二元碰撞 (8)

Ω_二元碰撞＝Ω_1,i+Ω_2，i+F (9)

τ＝3v₂+0.5 (11)

ρ＝ρ₁+ρ₂ (14)

式(11)是两相格子Boltzmann方程，f_i就是代表第i个离散方向上的分布函数，与空间坐标和时间有关，没有具体的函数表达式，其实质意义为记录每一时间步内的计算状态，在计算过程中，每一个时间步计算完成后，都会根据Ω_二元碰撞和f_i ^eq去更新f_i的值。其中i取值为0，1，2…，18，对应D3Q19中的19个离散速度。

r为网格的空间坐标，t为时间，Δt为时间步长，e_i是一组离散速度方向，每一个e_i有3个分量，其取值用矩阵表示为：

其中每一列分别对应一个e_i，第一列对应e₀＝(0，0，0)，第二列对应e₁＝(1，0，0)，以此类推e₂＝(-1，0，0)，e₃＝(0，1，0)，e₄＝(0，-1，0)，e₅＝(0，0，1)，e₆＝(0，0，-1)，e₇＝(1，1，0)，e₈＝(-1，-1，0)，e₉＝(1，-1，0)，e₁₀＝(-1，1，0)，e₁₁＝(1，0，1)，e₁₂＝(-1，0，-1)，e₁₃＝(1，0，-1)，e₁₄＝(-1，0，1)，e₁₅＝(0，1，1)，e₁₆＝(0，-1，-1)，e₁₇＝(0，1，-1)，e₁₈＝(0，-1，1)。

式(12)为f_i最终趋近的平衡分布函数f_i ^eq，w_j是权系数，j的取值与空间维数有关，在D3Q19模型中，j的取值为0，1，2，对应 ρ₂是粉尘颗粒的密度，u是整个流场的宏观速度。

Ω_二元碰撞为碰撞模型，分为自碰撞项Ω_1，i和互碰撞项Ω_2，i。对于自碰撞项Ω_1，i，采用LBGK模型τ是松弛时间，取值与粉尘颗粒的运动粘度v₂有关，τ＝3v₂+0.5。

对于互碰撞项，本发明引入二元碰撞模型：

上式代表了两种流体间的相互作用。D是空间维度，本发明是在3维空间下模拟，故D取值为3。C_σ被定义为分子的特异速度，在本发明中，视为LBM网格模型规定好的离散速度相对于流体宏观速度的速度e_i-u，其中u是整个流场的宏观速度，取值为ρ是整个流场的宏观密度，为空气的密度和粉尘颗粒的密度之和。Δf_i可以理解为分布函数f_i在这一碰撞过程中发生的变化，碰撞的效应是改变f_i使其趋于麦克斯韦平衡态分布f_i ^eq，设改变率的大小和f_i与f_i ^eq的差值成正比，可得本发明使用的互碰撞项如下：

其中ρ₂是粉尘颗粒的密度，e_i为LBM计算过程中预设的离散速度，u是整个流场的宏观速度，u₁和u₂分别为空气的速度和粉尘颗粒的速度。

步骤2.2：确定外力项

式(9)中的F为外力作用项，表示为：

其中w_j是LBM计算过程中每个离散速度方向的权重系数，前文已经提到，是已知量。为前面提到的空气与粉尘颗粒的相互作用力以及粉尘颗粒受到的涡度约束。

整体计算流程如附图1所示。

步骤3：系统的设计与实现

步骤3.1：确定系统开发平台

本发明所实现的室内烟雾扩散模拟系统是使用C++语言和开放性图形接口OpenGL在VScode中编程实现的，通过在VScode中配置OpenGL图形函数库，充分发挥OpenGL较高的图像渲染能力和可移植性的特点，使得本发明所实现的龙卷风模拟系统能够较好地满足真实性的要求。

步骤3.2：流场数据的组织

由于计算着色器的计算结构与流场的三维网格划分非常相似，在进行方程求解并行算法设计时，本文将流场划分的三维网格与计算着色器结构进行统一，每个网格对应一个工作项，这样在方程求解时，将龙卷风属性数据以三维纹理格式进行存储，使得每个网格区域每阶段的计算作为一个独立的任务分配给各个工作项独立执行，从而使得方程求解各步计算以并行的方式执行。在计算流程中需要表述龙卷风状态的属性数据包括速度、密度以及每个网格点的位置坐标。其中密度、温度等只包含一种数据，速度和位置数据则分为x，y，z三个方向的分量，这样所有的龙卷风属性共需要七个分量。为将龙卷风属性数据传递到GPU中，结合GL_RGBA四分量纹理数据格式，算法将所有数据组合成两个四个分量三维纹理：将速度三个分量数据组织到一个纹理数据中，第四个分量设置为0。将位置坐标和密度数据组织到一个纹理数据中，直接传入GPU与帧缓存进行绑定供后续计算使用。假设网格规模为4³，在CPU中分配64个规模大小的数据空间，用data1,data2...data64表达，与之对应的三维纹理的大小则为4ⅹ4ⅹ4。CPU到GPU的数据映射方式为：CPU中的数据将从data1开始按x方向、y方向、z方向的顺序依次映射到三维纹理缓存中，最后形成每一个三维纹理数据点代表一个网格点的效果。

步骤3.3：系统总体设计

龙卷风模拟系统总体上分为三层，分别是用户层，数据处理层和数据存储层，附图2为模拟系统的总体框架图。

下面分别对系统总体框架图中每一层的功能进行详细说明。

(1)用户层

用户通过执行程序，可以对系统界面进行操作，并通过操作鼠标或者键盘向系统发出命令，例如旋转或者移动观察视角，系统将收到的命令传递给处理层。

(2)数据处理层

数据处理层是模拟系统的关键部分，当处理层接收到由用户传给系统的命令后，首先会对得到的命令进行解析，并进行相应的计算，包括求解龙卷风的速度场，密度场等，计算空气与粉尘颗粒间的相互作用力，粉尘颗粒的涡度约束和处理边界，同时和数据存储层进行数据传输，最后将计算得到的密度场数据绘制到屏幕上，并反馈给用户。

(3)数据存储层

数据存储层主要与处理层进行数据交换和数据传输，存储流场速度场和密度场等数据。

步骤3.4：各分功能设计

根据附图2模拟系统的总体框架图，在满足龙卷风模拟真实性的要求下，对该系统的功能进行设计，如附图3所示，该系统的功能主要分为五个模块，分别是参数设置模块，LBM碰撞模块，合外力计算模块，LBM流动模块及边界约束模块以及龙卷风绘制模块。

(1)参数设置模块

在求解过程中，为了计算龙卷风的速度，压力和密度在任意时间步长下的值，需要设置一些初始条件，包括所有网格单元物理量的初始值，龙卷风起始位置，通过传入的参数计算初始平衡分布。另外，通过对鼠标进行监听，通过GLFW内置函数mouse_callback获取鼠标移动的俯仰角pitch和偏航角yaw，其中俯仰角是描述摄像机镜头往上看和往下看的角，偏航角是描述摄像机镜头往左看和往右看的角。计算过程中取龙卷风的底部中心为整个渲染界面的中心点，具体取为(0，0，0)。通过下式计算摄像机视角与渲染界面中心点的方向向量direction：

将该方向向量传入GLM内置函数lookAt中，即可实时根据鼠标移动转换摄像机视角。上述提到的GLFW和GLM，均为openGL图形库内置的图形框架。

(2)LBM碰撞模块

根据式(9)，执行碰撞过程，这一过程中，格子的物理量会发生变化，将计算结果及时储存，以供当前时间步内的流动过程以及下一个时间步的碰撞过程使用。

(3)合外力计算模块

为了增强模拟的细节，将两层流体之间的相互作用力以及涡度约束做为外力的一部分，分别进行求解，将它们合并后加入到LBM方程中。

(4)LBM流动模块及边界约束模块

执行流动过程，将网格上的分布函数移动到相邻网格，并计算流体宏观密度、速度，以及两个流体的密度、速度。再根据不同类型的边界单元，计算龙卷风遇到边界后的速度和密度的变化，从而对龙卷风的运动进行约束。

(5)龙卷风绘制模块

绘制模块根据龙卷风的密度场将龙卷风的旋转状态渲染到屏幕上。本发明使用了基于多色RGB值的纹理生成算法。

首先根据渲染场景的背景颜色选定一组烟雾的颜色的RGB值，为了能够突出龙卷风运动细节，选择与背景颜色对比明显的颜色来作为渲染色。之后根据选定龙卷风颜色数量对龙卷风密度进行分组，其具体方法为：假设选定龙卷风RGB值个数为N，初始密度为M，则将0到M之间数值分成N组，每一组对应一种RGB值。在本发明中，渲染场景背景色为暗蓝色，RGB值为(115,140,153)，选择4个颜色作为龙卷风的渲染色，分别为银白色，RGB(192,192,192)；深灰色，RGB(169,169,169)；灰色，RGB(128,128,128)；暗灰色，RGB(105,105,105)，是由浅到深的4种颜色。每个网格点在纹理计算时根据当前密度值所在密度分组确定该点的RGB值，该点RGBA值中透明度计算方式如下所示：

式中a为最终的透明度值，d是该网格龙卷风密度值，d_max为当前网格所处密度数据范围内最大值，d_min为当前网格所处密度数据范围内最小值。

Claims (2)

1.基于LBM的快速流体模拟方法，其特征在于：

步骤1：基于空气与粉尘颗粒的双流体龙卷风模型构建

步骤1.1：在传统模型的基础上，构建一个空气与粉尘颗粒的双流体模型

使用的龙卷风模型，基于N-S方程改造得来，具体如下：

其中为哈密顿算符，其计算公式为：/>对任意一个量执行哈密顿算符，等于对该量的所有分量求偏导之和；

把空气视为第一个流体，空气流动速度为u₁,式(1)为空气这一层流体的连续性方程；把粉尘尘埃等颗粒物视为第二个流体，流动速度为u₂，式(2)为粉尘颗粒物这一层流体的连续性方程；对于式(3)，p₁为空气这一层流体的压力场；对于式(4)，p₂为粉尘颗粒这一层流体的压力场；在式(3)和式(4)中，f₁是空气与粉尘颗粒之间的相互作用力，f₂是粉尘颗粒所受到的涡度约束；

对于式(3)、式(4)和式(5)中这一项，为空气与粉尘颗粒间相互作用力，ρ₂为粉尘颗粒这一层流体的密度,m₂为粉尘颗粒的质量，v₂为粉尘颗粒这一层流体的运动粘度；

对于式(4)和式(6)中的为粉尘颗粒在运动中受到的涡度约束，v_x,v_y,v_z为粉尘颗粒速度矢量的三个分量，r为粉尘颗粒运动时绕所在平面做旋转运动的半径；

步骤2：两相流LBM方程构建

步骤2.1：确定碰撞模型

使用的格子模型为D3Q19模型，包括3维空间，19个离散速度，以下为使用的两相LBM模型：

f_i(r+e_iΔt，t+Δt)-f_i(r，t)＝Ω_二元碰撞 (8)

Ω_二元碰撞＝Ω_1，i+Ω_2，i+F (9)

τ＝3v₂+0.5 (11)

f_i ^eq＝w_jρ₂(1+3e_iu+4.5(eu)²-1.5u²) (12)

ρ＝ρ₁+ρ₂ (14)

式(11)是两相格子Boltzmann方程，f_i就是代表第i个离散方向上的分布函数，与空间坐标和时间有关，没有具体的函数表达式，其实质意义为记录每一时间步内的计算状态，在计算过程中，每一个时间步计算完成后，都会根据Ω_二元碰撞和f_i ^eq去更新f_i的值；其中i取值为0，1，2…，18，对应D3Q19中的19个离散速度；

r为网格的空间坐标，t为时间，Δt为时间步长，e_i是一组离散速度方向，每一个e_i有3个分量，其取值用矩阵表示为：

其中每一列分别对应一个e_i，第一列对应e₀＝(0，0，0)，第二列对应e₁＝(1，0，0)，以此类推e₂＝(-1，0，0)，e₃＝(0，1，0)，e₄＝(0，-1，0)，e₅＝(0，0，1)，e₆＝(0，0，-1)，e₇＝(1，1，0)，e₈＝(一1，-1，0)，e₉＝(1，-1，0)，e₁₀＝(-1，1,0)，e₁₁＝(1，0，1)，e₁₂＝(-1,0，-1)，e₁₃＝(1，0，-1)，e₁₄＝(-1，0，1)，e₁₅＝(0，1，1)，e₁₆＝(0，-1，-1)，e₁₇＝(0，1，-1)，e₁₈＝(0，-1,1)；

式(12)为f_i最终趋近的平衡分布函数f_i ^eq,w_j是权系数，j的取值与空间维数有关，在D3Q19模型中，j的取值为0,1,2，对应 ρ₂是粉尘颗粒的密度，u是整个流场的宏观速度；

Ω_二元碰撞为碰撞模型，分为自碰撞项Ω_1,i和互碰撞项Ω_2，i；对于自碰撞项Ω_1，i，采用LBGK模型τ是松弛时间，取值与粉尘颗粒的运动粘度v₂有关，τ＝3v₂+0.5；

对于互碰撞项，引入二元碰撞模型：

上式代表了两种流体间的相互作用；D是空间维度，是在3维空间下模拟，故D取值为3；C_σ被定义为分子的特异速度，在中，视为LBM网格模型规定好的离散速度相对于流体宏观速度的速度e_i-u，其中u是整个流场的宏观速度，取值为ρ是整个流场的宏观密度，为空气的密度和粉尘颗粒的密度之和；Δf_i可以理解为分布函数f_i在这一碰撞过程中发生的变化，碰撞的效应是改变f_i使其趋于麦克斯韦平衡态分布f_i ^eq，设改变率的大小和f_i与f_i ^eq的差值成正比，可得使用的互碰撞项如下：

其中ρ₂是粉尘颗粒的密度，e_i为LBM计算过程中预设的离散速度，u是整个流场的宏观速度，u₁和u₂分别为空气的速度和粉尘颗粒的速度；

步骤2.2：确定外力项

式(9)中的F为外力作用项，表示为：

其中w_j是LBM计算过程中每个离散速度方向的权重系数，前文已经提到，是已知量；为前面提到的空气与粉尘颗粒的相互作用力以及粉尘颗粒受到的涡度约束。

2.根据权利要求1所述的基于LBM的快速流体模拟方法，其特征在于还包括：

步骤3：系统的设计与实现

步骤3.1：确定系统开发平台

所实现的室内烟雾扩散模拟系统使用C++语言和开放性图形接口OpenGL在VScode中编程实现的；

步骤3.2：流场数据的组织

将流场划分的三维网格与计算着色器结构进行统一，每个网格对应一个工作项，这样在方程求解时，将龙卷风属性数据以三维纹理格式进行存储，使得每个网格区域每阶段的计算作为一个独立的任务分配给各个工作项独立执行，从而使得方程求解各步计算以并行的方式执行；在计算流程中需要表述龙卷风状态的属性数据包括速度、密度以及每个网格点的位置坐标；其中密度、温度只包含一种数据，速度和位置数据则分为x，y，z三个方向的分量，这样所有的龙卷风属性共需要七个分量；为将龙卷风属性数据传递到GPU中，结合GL_RGBA四分量纹理数据格式，算法将所有数据组合成两个四个分量三维纹理：将速度三个分量数据组织到一个纹理数据中，第四个分量设置为0；将位置坐标和密度数据组织到一个纹理数据中，直接传入GPU与帧缓存进行绑定供后续计算使用；假设网格规模为4³，在CPU中分配64个规模大小的数据空间，用data1,data2...data64表达，与之对应的三维纹理的大小则为4x4x4；CPU到GPU的数据映射方式为：CPU中的数据将从data1开始按x方向、y方向、z方向的顺序依次映射到三维纹理缓存中，最后形成每一个三维纹理数据点代表一个网格点的效果；

步骤3.3：系统总体设计

龙卷风模拟系统总体上分为三层，分别是用户层，数据处理层和数据存储层；

下面分别对系统总体框架图中每一层的功能进行详细说明；

(1)用户层

用户通过执行程序，可以对系统界面进行操作，并通过操作鼠标或者键盘向系统发出命令，系统将收到的命令传递给处理层；

(2)数据处理层

数据处理层是模拟系统的关键部分，当处理层接收到由用户传给系统的命令后，首先会对得到的命令进行解析，并进行相应的计算，包括求解龙卷风的速度场，密度场，计算空气与粉尘颗粒间的相互作用力，粉尘颗粒的涡度约束和处理边界，同时和数据存储层进行数据传输，最后将计算得到的密度场数据绘制到屏幕上，并反馈给用户；

(3)数据存储层

数据存储层与处理层进行数据交换和数据传输，存储流场速度场和密度场等数据；

步骤3.4：各分功能设计

该系统的功能主要分为五个模块，分别是参数设置模块，LBM碰撞模块，合外力计算模块，LBM流动模块及边界约束模块以及龙卷风绘制模块；

(1)参数设置模块

在求解过程中，为了计算龙卷风的速度，压力和密度在任意时间步长下的值，需要设置一些初始条件，包括所有网格单元物理量的初始值，龙卷风起始位置，通过传入的参数计算初始平衡分布；另外，通过对鼠标进行监听，通过GLFW内置函数mouse_callback获取鼠标移动的俯仰角pitch和偏航角yaw，其中俯仰角是描述摄像机镜头往上看和往下看的角，偏航角是描述摄像机镜头往左看和往右看的角；计算过程中取龙卷风的底部中心为整个渲染界面的中心点，具体取为(0，0，0)；通过下式计算摄像机视角与渲染界面中心点的方向向量direction：

将该方向向量传入GLM内置函数lookAt中，即可实时根据鼠标移动转换摄像机视角；上述提到的GLFW和GLM，均为openGL图形库内置的图形框架；

(2)LBM碰撞模块

根据式(9)，执行碰撞过程，这一过程中，格子的物理量会发生变化，将计算结果及时储存，以供当前时间步内的流动过程以及下一个时间步的碰撞过程使用；

(3)合外力计算模块

为了增强模拟的细节，将两层流体之间的相互作用力以及涡度约束做为外力的一部分，分别进行求解，将它们合并后加入到LBM方程中；

(4)LBM流动模块及边界约束模块

执行流动过程，将网格上的分布函数移动到相邻网格，并计算流体宏观密度、速度，以及两个流体的密度、速度；再根据不同类型的边界单元，计算龙卷风遇到边界后的速度和密度的变化，从而对龙卷风的运动进行约束；

(5)龙卷风绘制模块

绘制模块根据龙卷风的密度场将龙卷风的旋转状态渲染到屏幕上；

首先根据渲染场景的背景颜色选定一组烟雾的颜色的RGB值，为了能够突出龙卷风运动细节，选择与背景颜色对比明显的颜色来作为渲染色；之后根据选定龙卷风颜色数量对龙卷风密度进行分组，其具体方法为：假设选定龙卷风RGB值个数为N，初始密度为M，则将0到M之间数值分成N组，每一组对应一种RGB值；在中，渲染场景背景色为暗蓝色，RGB值为(115，140,153)，选择4个颜色作为龙卷风的渲染色，分别为银白色，RGB(192,192,192)；深灰色，RGB(169,169,169)；灰色，RGB(128,128，128)；暗灰色，RGB(105，105，105)，是由浅到深的4种颜色；每个网格点在纹理计算时根据当前密度值所在密度分组确定该点的RGB值，该点RGBA值中透明度计算方式如下所示：

式中a为最终的透明度值，d是该网格龙卷风密度值，d_max为当前网格所处密度数据范围内最大值，d_min为当前网格所处密度数据范围内最小值。