CN111428423A - 实现总能量守恒的晶格玻尔兹曼求解器 - Google Patents

Abstract

本公开涉及实现总能量守恒的晶格玻尔兹曼求解器。描述了用于使用用于求解标量传输方程以及解决总能量的晶格玻尔兹曼(LB)方法模拟流体流的技术。除了用于流体流的晶格玻尔兹曼(LB)函数之外，该技术还包括通过将特定总能量添加到要平流的粒子的状态并且从在时间间隔将不平流的粒子的状态减去特定总体能量来修正粒子的状态向量集，并且根据修正的状态集执行粒子平流。

Description

实现总能量守恒的晶格玻尔兹曼求解器

优先权要求

本申请根据35 USC§119(e)要求于2019年1月10日提交的， 名称为“LATTICEBOLTZMANN SOLVER ENFORMCING TOTAL ENERGY CONSERVE SERVICE(实现总能量守恒的晶格 玻尔兹曼求解器)”的美国临时专利申请第62/790,528号的优先权， 其全部内容通过引用合并于此。

背景技术

本描述涉及诸如物理流体流动的物理过程的计算机模拟。

所谓的“晶格玻尔兹曼方法”(Lattice Boltzmann Method， LBM)是一种在计算流体动力学中使用的有利技术。晶格玻尔兹曼 系统的底层动力学在于动力学理论的基本物理学，其涉及许多粒子根 据玻尔兹曼方程的运动。基础玻尔兹曼动力学系统中有两个基本的动 力学过程——碰撞和平流。碰撞过程涉及遵守守恒定律并弛豫至平衡 的粒子之间的相互作用。平流过程涉及根据粒子的微观速度对粒子从 一个位置到另一个位置的移动进行建模。

对于可压缩流，为了满足质量、动量和总能量守恒，通常使用单 个分布函数。如果使用单个分布函数来求解所有三个量，则用于 LBM的晶格应满足直至八阶的矩(moment)。在这种情况下，随 着矩需求的增加，需要大量的晶格速度，并且所产生的模板尺寸也随 之增加。使用单个粒子分布函数涉及的其他考虑因素包括分布函数稳 定性、动量和能量两者的边界条件、有关体积力(body force)和热 源的考虑因素等等。

发明内容

根据一个方面，一种用于在计算机上模拟流体流动(其中该模拟 实现总能量守恒)的方法包括模拟流体在网格上的活动，该流体的活 动被模拟以便对粒子在网格上的移动进行建模，在计算机可访问的存 储器中为网格中的每个网格位置存储状态向量集，每个状态向量包括 多个条目，这些条目与对应网格位置处的可能动量状态中的特别(particular)动量状态相对应，通过计算机为网格中的网格位置计 算能量值集，通过计算机对于时间间隔执行粒子到后续网格位置的平 流，以及通过计算机通过将特定总能量值添加到平流的粒子的状态并 从在该时间间隔未平流的粒子的状态中减去所述特定总能量值，修正 粒子的状态向量。

各方面还包括计算机程序产品、一个或多个机器可读硬件存储设 备、装置和计算系统。

在本文所述的特征之中，上述方面中的一个或多个可以包括以下 特征中的一个或多个。

在粒子根据修正的状态集平流到后续网格位置之后，该方面包括 在计算矩之前将局部压力项加回到粒子状态以提供适当的压力梯度

局部压力项包括由计算机计算的θ-1个项。模拟流体流动 的活动包括部分地基于第一离散晶格速率集来模拟流体流动。所述方 法还包括部分地基于第二离散晶格速率集来模拟标量的时间演化。第 二离散晶格速率集与第一离散晶格速率集相同。第二离散晶格速率集 的晶格速率的数量与第一离散晶格速率集不同。模拟总能量的时间演 化包括收集对于碰撞和能量算子的来自相邻网格位置的传入分布，对 传入分布进行加权，确定作为碰撞和能量算子乘积的传出分布，并传 播所确定的传出分布。

模拟标量的时间演化包括部分基于碰撞算子模拟标量的时间演化， 该碰撞算子有效地将由碰撞和能量算子的乘积得到的动量通量表达为 一阶项。为了对于任意普朗特(Prandtl)数恢复精确剪切应力，能 量碰撞项由下式给出：

其中，Π＝Σ_ic_ic_i(f_i-f_i ^eq是非平衡分量的经过滤的二阶矩。

该方面还包括应用零净表面通量边界条件，使得传入分布等于所 确定的传出分布。确定传出分布包括确定传出分布以提供零表面标量 通量。标量包括从由温度、浓度和密度组成的组中选择的标量。为网 格中的网格位置计算能量值集包括计算由E_t＝E+v²/2给出的总能 量E_t，其中E是能量，ν是速度。

模拟活动包括由计算机对于特定总能量E_t应用平衡分布和第二分 布函数，其中第二分布被定义为随着流分布f_i平流的特定标量。第二 分布函数考虑了f_i对能量方程的非平衡贡献，以获得边界附近和跨不 同网格分辨率的正确流动行为，其中，用于分布函数的碰撞算子由下 式给出：

其中项Ω_f,Ω_q表示相应的碰撞算子；上述方程中使用的平衡分布 为：

总能量在平流之前添加该项，使得压力被对流，并且为了补偿增 加的能量，将能量从停止状态中去除。去除所增加的能量可以使总能 量守恒并提供正确的压力速度项。通过修正停止状态

以使得

可以提供所增加的能量的去除。

在本文公开的技术中，求解总能量使用被定义为随着流分布平流 的特定标量的第二分布。本文描述的用于模拟流体流动的技术使用晶 格玻尔兹曼(LB)方法和标量求解器传输方程。提供了一种用于可 压缩流的具有总能量守恒和更高稳定性范围的LBM求解器，该求解 器可用于现实世界应用。在碰撞过程中，求解器中采用的平衡分布为 正，在平流过程中进行温度耦合。

根据包括附图和权利要求的以下描述，其它特征和优点将变得清 楚。

附图说明

图1描绘了用于流体流动模拟的系统。

图2-3描绘了流程图，其示出了用于复杂流体流动模拟以及总能 量计算的操作。

图4-6是说明性的基于计算的预期、预测结果的比较。

图7和8图示了两个LBM模型的速度分量(现有技术)。

图9是物理过程模拟系统所遵循的过程的流程图。

图10是微观块(microblock)的透视图(现有技术)。

图11A和11B是图1的系统使用的晶格结构(现有技术)的图 示。

图12和13图示了可变分辨率技术(现有技术)。

图14图示了有助于理解粒子与面元的相互作用的平行六面体。

图15图示了粒子从体元到表面的移动(现有技术)。

图16图示了粒子从表面到表面的移动(现有技术)。

图17是用于执行表面动力学的过程(现有技术)的流程图。

图18是具有与粒子碰撞步骤无关的热力学步骤的用于产生分布 函数的过程(现有技术的流程图)。

具体实施方式

在下面图9中讨论的过程中，描述了使用实现总能量守恒的晶格 玻尔兹曼求解器的流动模拟过程。

在图7，8和10-18中，这些附图均被标记为“现有技术”。这 些附图被标记为现有技术是因为这些附图通常出现在上文所引用的专 利中。然而，这些附图如它们出现在上文引用的专利申请中那样。 但是，这些附图没有考虑对于如下所述的使用实现总能量守恒的晶格 玻尔兹曼求解器的流动模拟的任何修改，这是因为在上文所引用的专 利申请中没有描述LB求解器。

A.求解标量的方法

当完成复杂流体流动模拟时，结合求解流体流动而同时求解标量 (例如温度分布、浓度分布和/或密度)可能是有益的。在本文所述 的系统和方法中，基于基于LBM的物理过程模拟系统，(与矢量相 对的)标量的建模与流体流的建模耦合。可以被建模的示例性标量包 括温度、浓度和密度。

在Pradeep Gopalakrishnan等人的标题为：“Temperature Coupling Algorithmfor Hybrid Thermal Lattice Boltzmann Method (用于混合热晶格玻尔兹曼方法的温度耦合算法)”的共同未决的美 国专利公开US-2016-0188768-A1中详细介绍了将标量建模与模拟流 体流动耦合的技术的具体细节，该美国专利公开的全部内容通过引用 合并于此。使用基于具有相对高的稳定范围的伽利略不变碰撞算子的 碰撞算子。碰撞算子的具体细节在美国专利9,576,087中有详细描述， 该美国专利的全部内容通过引用结合在此。

例如，该系统可被用于确定系统内的对流温度分布。例如，如果 系统(由多个体元所表示的体积形成)包括热源，并且系统内有气流， 则根据气流和与热源的接近程度，系统的某些区域将比其他区域更热。 为了对这种情况进行建模，可以将系统内的温度分布表示为标量，其 中每个体元均具有关联的温度。

在另一个示例中，该系统可被用于确定系统内的化学分布。例如， 如果系统(由多个体元所表示的体积形成)包括污染物源，例如脏弹 或化学物质或其他悬浮在空气或液体中的微粒，并且系统内存在空气 或液体流动，基于流动和与源的接近程度，系统的某些区域将具有比 其他区域更高的浓度。为了对这种情况进行建模，可以将系统内的化 学分布表示为标量，其中每个体元均具有关联的浓度。

在某些应用中，可以同时模拟多个不同的标量。例如，系统可以 模拟系统中的温度分布和浓度分布两者。

标量可被以不同的方式建模。例如，用于求解标量传输方程的晶 格玻尔兹曼(LB)方法可用于间接地求解标量传输。例如，本文描 述的方法可以提供以下二阶宏观标量传输方程的间接解：

除了用于流体流动的晶格玻尔兹曼函数外，还引入了第二组分布 函数用于传输标量。该方法为一个体积中的每个体元分配向量以表示 流体流，并且给该体积中的每个体元分配标量以表示所需的标量变量 (例如，温度、密度、浓度等)。该方法完全恢复了满足精确守恒定 律的宏观标量传输方程。该方法被认为与其他非LBM方法相比，可 以提高所确定的标量的准确性。另外，这种方法提供了增强的解决复 杂边界形状的能力。

可以将这种用于对标量进行建模的方法与基于晶格玻尔兹曼方法 (LBM)的时间显式CFD/CAA解决方案方法(例如可从马萨诸塞 州伯灵顿的Exa公司获得的PowerFLOW系统)结合使用。与基于 离散化宏观连续方程的方法不同，LBM从“介观(mesoscopic)” 玻尔兹曼动力学方程开始，以预测宏观流体动力学。所得的可压缩且 不稳定的求解方法可用于预测各种复杂的流动物理学，例如航空声学 和纯声学问题。下面提供了对于基于LBM的模拟系统的总体讨论， 然后是可以与流体流动模拟结合使用以支持这种建模方法的标量求解 方法的讨论。

对模拟空间建模

在基于LBM的物理过程模拟系统中，流体流动可以由以一组离 散速度c_i进行评估的分布函数值f_i表示。分布函数的动力学是由方程 (I.1)支配的，其中f_i(0)已知为平衡分布函数，被定义为：

其中，

这个方程是描述分布函数f_i的时间演化的众所周知的晶格玻尔兹 曼方程。左手侧表示由于所谓的“流化过程”造成的分布的变化。流 化过程是一块流体在网格位置出发，然后沿其中一个速度向量移动到 下一个网格位置。在那个点，计算“碰撞算子”，即，附近流体块对 出发的流体块的影响。流体只能移动到另一个网格位置，因此速度向 量的适当选择是必要的，使得所有速度的所有分量都是共同速率 (common speed)的倍数。

第一个方程的右手侧是上面提到的“碰撞算子”，其表示由于流 体块之间的碰撞造成的分布函数的变化。这里使用的碰撞算子的特定 形式是Bhatnagar、Gross和Krook(BGK)。它迫使分布函数去向 由第二个方程给出的规定值，其中第二个方程是“平衡”形式。

BGK算子是根据以下物理论点来构造的，即不管碰撞的细节如 何，分布函数都经由碰撞接近由

{f^eq(x,v,t)}给出的良好 定义的局部平衡：

其中参数τ表示经由碰撞达到平衡的特征弛豫时间。处理粒子(例如 原子或分子)时，通常将弛豫时间视为常数。

根据这个模拟，常规的流体变量，诸如质量ρ和流体速度u，作 为方程(I3)中的简单求和而被获得：

其中，ρ，μ和T分别是流体密度、速度和温度，并且D是离散 速度空间的维度(不一定等于物理空间维度)。

由于对称性考虑，速度值集以这样一种方式来选择：该方式使得 它们在配置空间中跨越时形成某种晶格结构。这种离散系统的动力学 遵循具有如下形式的LBE:

f_i(x+c_i,t+1)-f_i(x,t)＝C_i(x,t)

其中碰撞算子通常采取如上所述的BGK形式。通过平衡分布形式的 适当选择，可以从理论上示出，晶格玻尔兹曼方程产生正确的流体动 力学和热-流体动力学。即，从f_i(x,t)导出的流体动力矩在宏观限制 中服从Navier-Stokes方程。这些矩由上述方程(I3)定义。

c_i和w_i的集合值定义LBM模型。LBM模型可以高效地在可扩 展计算机平台上实现并且对于时间不稳定流和复杂的边界条件以最大 的健壮性运行。

从玻尔兹曼方程获得用于流体系统的运动的宏观方程的标准技术 是Chapman-Enskog方法，在该方法中采取全玻尔兹曼方程的逐次 逼近。在流体系统中，密度的小扰动以声速行进。在气体系统中，声 速一般由温度确定。流体中可压缩性的效果的重要性是由特征速度与 声速的比来测量的，该比被称为马赫数。对于常规的基于LBM的物 理过程模拟系统的进一步说明，参考通过引用申请US-2016-0188768- A1而并入的上文。

参考图1，示出了用于模拟(例如，物理对象的表示附近的)流 体流动的系统10。该实现方式中的系统10基于客户端-服务器架构， 并且包括被实现为大规模并行计算系统12的服务器系统12和客户端 系统14。服务器系统12包括存储器18、总线系统11、接口20(例如，用户接口/网络接口/显示器或监视器接口等)和处理设备24。在 存储器18中的是网格准备引擎32和模拟引擎34。系统10访问存储 2D和/或3D网格、坐标系和库的数据存储库38。

虽然图1示出了存储器18中的网格准备引擎32，但网格准备引 擎可以是在与服务器12不同的系统上执行的第三方应用。无论网格 准备引擎32是在存储器18中执行还是在不同于服务器12的系统上 执行，网格准备引擎32都接收用户提供的网格定义30，并且网格准 备引擎32准备网格并将所准备的网格发送到模拟引擎34。模拟引擎 34包括粒子碰撞交互模块、粒子边界建模模块和执行如下所述的平 流操作的平流模块。模拟引擎34还包括总能量求解器，如下面在图 3中讨论的。

现在参考图2，示出了用于模拟物理对象的表示附近的流体流动 的过程40。在文中将描述的示例中，物理对象是翼型体(airfoil)。 但是，翼型体的使用仅仅是说明性的，这是因为物理对象可以是任何 形状，并且特别地可以具有(一个或多个)平坦和/或弯曲的表面。 该过程例如从客户端系统14或通过从数据储存库42检索而接收用于 正模拟的物理对象的网格。在其它实施例中，外部系统或服务器12 基于用户输入生成用于正模拟的物理对象的网格。该过程根据检索到 的网格预先计算44几何量，并使用与检索到的网格对应的预先计算 的几何量执行46动态晶格玻尔兹曼模型模拟。晶格玻尔兹曼模型模 拟包括LBM网格中粒子分布的演化和粒子根据修正的状态向量集到 下一单元

的平流的模拟，该状态向量集是通过将特定总能量值添加 到平流的粒子的状态中并且从在时间间隔内未平流的粒子的状态减去 该特定总能量值而被修正的。

总能量守恒可压缩流求解器

当完成复杂的流体流动模拟时，结合求解流体流动而同时求解标 量(例如温度分布、浓度分布和/或密度)可能是有益的。区分特定 标量与标量分布，标量分布为q(方程10)和h(方程2)之间的差 异，并使用特定标量q为E_t，而h为密度*E_t。讨论的方法是为在 被对流的流的顶部的标量安置，如上面提到的美国专利公开US- 2016-0188768-A1所述。

总能量守恒中压力项的处理

现在参考图3，现在描述具有总能量守恒过程50的可压缩LBM 求解器，该可压缩LBM求解器具有比现有方法更好的稳定性范围， 并且适用于现实世界的应用。需要求解的总能量守恒方程式为：

其中总能量E_t由E_t＝E+v²/2给出，ρ是密度，v是速度，p 是压力，k是导热系数，T是温度，τ是剪切应力。这是用于高速流 的所有计算流体动力学工具都需要求解的一个通用方程。

总能量求解器涉及的复杂性之一是在方程1中添加了压力对流项 “p”。与用于高速流的所有计算流体动力学工具一样，晶格玻尔兹 曼方法(LBM)必须针对总能量守恒来求解该方程。在所有情况下， 尤其是对于压力对流项p，LBM方法在求解该方程式时都有一些缺 点。所公开的方法能够对于复杂问题解决这一点并具有其他优点。

LBM方法存在的压力对流问题可以解释如下：正常的标量传输 具有以下方程

上述方程的左手侧包括对流项，而该方程的右手侧具有扩散项。 与总能量方程不同，上述方程中的所有项仅包含一个标量变量S，这 使得平流更为简单。过程50根据平衡分布的选定形式来模拟52粒子 分布的演化。

除了平衡分布(方程1的选定形式)52之外，还使用54用于计 算特定总能量E_t的第二分布函数。过程50包括LBM中粒子到下一 单元q的平流56。通过这两种分布的乘积实现了满足能量守恒。这 些分布函数的乘积被定义为随着流分布f_i平流的特定标量。与独立使用第二分布函数来求解能量的现有方法不同，该第二分布函数考虑了 f_i对能量方程的非平衡贡献，以获得边界附近和跨不同网格分辨率的 正确流行为。

分布函数的碰撞算子由下式给出：

f_i′(x，t)＝f_i(x，t)+Ω_f[f_i(x，t)，f_i ^eq(x，t)] 方程2

方程2和方程3中的项Ω_f,Ω_q被用于代表各自的碰撞算子。美国 专利第9,576,087号中广泛地讨论了方程5中的碰撞算子Ω_f，该美国 专利的全部内容通过引用合并于此。上面方程中使用的平衡分布为

在用于总能量求解器的上述分布函数方程7中显见的观察结果是 平衡分布不涉及方程1中的压力项(p)和θ-1项，即p＝ρRT，这对 于热流至关重要。第二分布是特定总能量，其是内部能量E＝C_νT和 动能ν²/2的总和。

因此，该分布函数在碰撞操作期间将始终提供正值，因此该分布 函数提供了具有相应的高稳定性的求解器。在方程(2)中碰撞后的 状态按原样平流，方程4的分布函数将仅导致动量中的等温压力梯度， 而没有总能量方程中的压力对流项。压力梯度和压力对流项是由于平 流而产生的，而碰撞不会影响这些项，这与如下所示的使用平衡分布 并且在平流和碰撞两者期间都包括压力项的其他方法形成对比。

术语θ＝RT/T₀对应于压力并使上述平衡为负，因此现有技术高 度不稳定。还应注意，第二分布函数

被独立地用于总能量，这会 导致跨不同网格分辨率的噪声。

处理50通过将特定总能量添加到被平流的状态，同时从停止状 态(非运动状态)中去除特定总能量，来修正58平流状态。

更具体地，在该修正58中，为了恢复正确的动量和能量方程， 仅在平流步骤期间进行压力引入。在平流过程之前修正网格位置处的 运动状态

值，如下面的方程8所示：

但是，将此项(方程5中的第二个方程)添加到q_i也为运动状 态添加了额外的能量。为了补偿此项表示的增加的能量，从停止状态 (粒子不对流的状态)中除去了相同量的能量。通过这样做，该过程 节省了能量并提供了正确的压力速度项。具体地，通过修正停止状态

以使得

来提供总能量守恒。

平流后，在计算矩之前，将局部压力项加回到粒子状态，从而产 生适当的压力梯度

平流之前(方程8)和平流之后(方程6)的修正恢复了用于动 量的正确状态方程(或正确的压力梯度

)，并且还恢复了用于总 能量的正确的压力对流项

由于去除了压力项(方程9)， 因此在碰撞过程中状态变为正，并提高了总能量求解器的稳定性范围。

应用对于平流的总能量修正来对平流进行修正，在该总能量修正 中计算机通过将特定总能量值添加到平流的粒子的状态并从在该时间 间隔内未平流的粒子的状态中减去特定总能量值来修正粒子的状态向 量，如下文的图9中所讨论的。

能量的正则(Regularizd)碰撞算子

在上述美国专利9,576,087中讨论了用于流动状态的碰撞算子Ω_f。 下面讨论用于能量的碰撞算子Ω_q。具有一阶矩的正则碰撞由下式给 出：

方程7是经滤波的对总能量的非平衡贡献。通过这两个状态的相 乘来满足能量守恒，该相乘使得碰撞项、方程(7)的右手侧的第二 项

为：

上项的零矩为零并节省能量。在Chapman-Enskog扩展过程中， 碰撞项的一阶矩涉及热扩散的推导和由切应力完成的功。作为一阶近 似，碰撞项可以表示为：

在上述方程中，动量通量项ρυυ是两个状态相乘

的结果。此项 ρυυ导致针对高速流动的数值不稳定，并且还会导致依赖于马赫数的 扩散。为了减轻这些影响，通过将碰撞算子定义为下式来去除项ρυυ：

上面的表达式省略了动量通量项ρυυ，并提高了对于高马赫流的 数值稳定性。

对于碰撞算子，可以使用许多技术。其中一种技术称为“正则 化”。在Hudong和Raoyang的基于晶格玻尔兹曼的标量的专利 (未决)中，“正则化”与附加的伽利略不变性特征一起使用。但是， 此正则化碰撞算子对于高速流并非理想的。在推导过程中，算子最终具有诸如ρυυ(方程12)之类的高阶项。对于LBM，速度值通常始 终小于“1”(一个)，因此，对于较低的速度流，这些误差项将低 得多，并且可以忽略不计。但是，对于更高速率流，速度值将高达1 阶，因此误差项ρυυ相对较高，并且通过方程10被去除。

任意普朗特数的粘性效应

能量求解器中使用的弛豫时间为τ_q，这将导致粘度等于热扩散率 (τ_q-0.5ρRT)的剪切应力项完成不正确的功。为了对于任意普朗特 数恢复精确剪切应力，以下附加项被添加到能量碰撞项中。

其中Π＝Σ_ic_ic_i(f_i-f_i ^eq)是经滤波的非平衡分量的二阶矩。

总能量求解器能力

总能量求解器可被包含在诸如

的现有流求解器中， 并且该总能量求解器可被用于广泛的工业应用，该应用包括涉及高 亚音速和跨音速流的应用。

本节中讨论了几个基准，以显示所公开的对于可压缩流的总能量 求解器的某些潜在优势，特别是在能量守恒方面。将模拟结果与分析 解决方案和典型的基于PDE的有限差分混合求解器解决方案进行比 较。

图4示出了在马赫速率(Ma)下对于Ma＝0.7，x＝距前 沿1.5m，在平板上方的流的总温度的曲线图。如图4所示，用于在 板上方的流的总能量求解器(实线62)精确地恢复了总温度，而混 合求解器(LBM流求解器和基于有限差分的熵求解器)(虚线64) 不能保持总温度。

图5示出了在马赫速率(Ma)下对于Ma＝0.7，x＝距前沿 1.5m，平板的静态温度的曲线图。该图相对于作为u/u_∞(y轴)的u 来绘制作为T/T_∞(x轴)的温度T。将用实线66表示的本公开的总 能量求解器与混合求解器(LBM流求解器和基于有限差分的熵求解 器)(虚线68)进行比较。该图说明了对于Ma＝0.7，x＝距前 沿1.5m在马赫速率(Ma)下平板的静态温度与沃尔兹方程(星号线 70)相匹配。

参考图6示出针对背压P_b＝0:75P_t的会聚-发散验证(CDV)喷 嘴的温度模拟结果的曲线图。该曲线图相对于x(y轴)绘制温度T 为T/T_∞(x轴)。将以实线72表示的本公开的总能量求解器与混合 求解器(LBM流求解器和基于有限差分的熵求解器)(虚线74)进 行比较。本公开的总能量求解器与分析确定点划线76具有良好的一 致性。静态温度的比较表明，基于有限差分方法的可压缩求解器在冲 击后无法恢复温度，而本公开的总能量求解器实现了预期值。

具有公开的总能量求解器的系统可以用于确定系统内的对流温度 分布。例如，如果系统(由多个体元表示的体积形成)包括热源，并 且系统内有气流，则根据气流和与热源的接近程度，系统的某些区域 将比其他区域更热。为了对这样的系统建模，可以将系统内的温度分 布表示为标量，其中每个体元均都具有关联的温度。

在另一个示例中，该系统可被用于确定系统内的化学分布。例如， 如果系统(由多个体元所表示的体积形成)包括污染物源，例如脏弹 或化学物质或其他悬浮在空气或液体中的微粒，并且系统内存在空气 或液体流动，基于流动和与源的接近程度，系统的某些区域将具有比 其他区域更高的浓度。为了对这种情况进行建模，可以将系统内的化 学分布表示为标量，其中每个体元均具有关联的浓度。

在某些应用中，可以同时模拟多个不同的标量。例如，系统可以 模拟系统中的温度分布和浓度分布两者。

标量可被以不同的方式建模。例如，用于求解标量传输方程的晶 格玻尔兹曼(LB)方法可用于间接地求解标量传输。例如，本文描 述的方法可以提供以下二阶宏观标量传输方程的间接解：

在这样的布置模拟中，除了用于流体流动的晶格玻尔兹曼函数外， 还引入了第二组分布函数用于传输标量。该方法为一个体积中的每个 体元分配向量以表示流体流，并且给该体积中的每个体元分配标量以 表示所需的标量变量(例如，温度、密度、浓度等)。该方法完全恢 复了满足精确守恒定律的宏观标量传输方程。该方法被认为与其他非 LBM方法相比，可以提高所确定的标量的准确性。另外，这种方法 提供了增强的解决复杂边界形状的能力。

可以将这种用于对标量进行建模的方法与基于晶格玻尔兹曼方法 (LBM)的时间显式CFD/CAA解决方案方法(例如可从马萨诸塞 州伯灵顿的Exa公司获得的PowerFLOW系统)结合使用。与基于 离散化宏观连续方程的方法不同，LBM从“介观(mesoscopic)” 玻尔兹曼动力学方程开始，以预测宏观流体动力学。所得的可压缩且 不稳定的求解方法可用于预测各种复杂的流动物理学，例如航空声学 和纯声学问题。下面提供了对于基于LBM的模拟系统的总体讨论， 然后是可以与流体流动模拟结合使用以支持这种建模方法的标量求解 方法的讨论。

参考图7，第一模型(2D-1)100是包括21个速度的二维模型。 在这21个速度中，一个速度(105)表示不移动的粒子；三组四个速 度表示在沿晶格的x或y轴的正方向或负方向以归一化速率(r) (110-113)、归一化速率的两倍(2r)(120-123)或归一化速率的 三倍(3r)(130-133)移动的粒子；以及两组四个速度表示相对于x 和y晶格轴二者以归一化速率(r)(140-143)或归一化速率的两倍 (2r)(150-153)移动的粒子。

也如图8所示，示出了第二模型(3D-1)200——包括39个速 度的三维模型，其中每个速度由图8的箭头之一表示。在这39个速 度中，一个速度表示不移动的粒子；三组六个速度表示在沿晶格的x、 y或z轴的正方向或负方向以归一化速率(r)、归一化速率的两倍(2r)或归一化速率的三倍(3r)移动的粒子；八个速度表示相对于 x、y、z晶格轴当中的全部三个晶格轴以归一化速率(r)移动的粒 子；以及十二个速度表示相对于x、y、z晶格轴当中的两个晶格轴以 归一化速率的两倍(2r)移动的粒子。

也可以使用更复杂的模型，诸如包括101个速度的3D-2模型以 及包括37个速度的2D-2模型。

对于三维模型3D-2，在101个速度中，一个速度表示不移动的 粒子(组1)；三组六个速度表示在沿晶格的x、y或z轴的正方向 或负方向以归一化速率(r)、归一化速率的两倍(2r)或归一化速 率的三倍(3r)移动的粒子(组2、4和7)；三组八个速度表示相 对于x、y、z晶格轴当中的全部三个晶格轴以归一化速率(r)、归 一化速率的两倍(2r)或归一化速率的三倍(3r)移动的粒子(组3、 8和10)；十二个速度表示相对于x、y、z晶格轴当中的两个晶格轴 以归一化速率的两倍(2r)移动的粒子(组6)；二十四个速度表示 相对于x、y、z晶格轴的当中两个晶格轴以归一化速率(r)和归一 化速率的两倍(2r)移动并且不相对于剩下的轴移动的粒子(组5)； 以及二十四个速度表示相对于x、y、z晶格轴当中的两个晶格轴以归 一化速率(r)移动并且相对于剩下的轴以归一化速率的三倍(3r) 移动的粒子(组9)。

对于二维模型2D-2，在37个速度中，一个速度表示不移动的 粒子(组1)；三组四个速度表示在沿晶格的x或y轴的正或负方向 以归一化速率(r)、归一化速率的两倍(2r)或归一化速率的三倍 (3r)移动的粒子(组2、4和7)；两组四个速度表示相对于x晶 格轴和y晶格轴二者以归一化速率(r)或归一化速率的两倍(2r) 移动的粒子；八个速度表示相对于x晶格轴和y晶格轴之一以归一化 速率(r)移动并且相对于另一个轴以归一化速率的两倍(2r)移动 的粒子；以及八个速度表示相对于x晶格轴和y晶格轴之一以归一化 速率(r)移动并且相对于另一个轴以归一化速率的三倍(3r)移动 的粒子。

上述LBM模型在二维和三维中为流的数值模拟提供了特定类的 高效且健壮的离散速度动力学模型。这种类型的模型包括离散速度的 特别集合和与那些速度关联的权重。这些速度与速度空间中的笛卡尔 坐标的网格点相符，这有助于离散速度模型(尤其是被称为晶格玻尔 兹曼模型的种类)的准确和高效实现。利用这种模型，流可以以高保 真度被模拟。

参考图9，物理过程模拟系统根据过程300操作以模拟诸如流体 流动的物理过程。在模拟之前，将模拟空间建模为体元(voxel)的 集合(步骤302)。典型地，使用计算机辅助设计(CAD)程序生成 模拟空间。例如，CAD程序可以用于绘制位于风洞中的微器件。此 后，处理由CAD程序产生的数据，以添加具有适当分辨率的晶格结 构，并考虑模拟空间内的对象和表面。

晶格的分辨率可以基于正在被模拟的系统的雷诺数(Reynolds number)来选择。雷诺数涉及流的粘度(v)、流中的对象的特征长 度(L)和流的特征速度(u)：

Re＝uL/v 方程(I4)

对象的特征长度表示对象的大尺度特征。例如，如果微型设备周 围的流正在被模拟，则该微型设备的高度可以被认为是特征长度。当 对象的小区域(例如，汽车的侧视镜)周围的流是所关心的时，模拟 的分辨率可以增加，或者增加分辨率的区域可以在所关心的区域周围 被采用。体元的维度随着晶格的分辨率增加而减小。

状态空间被表示为f_i(x,t)，其中f_i表示在时间t在由三维向量x 表示的晶格位点处在状态i下的每单位体积的元素或粒子的数量(即， 状态i下的粒子的密度)。对于已知的时间增量，粒子的数量被简称 为f_i(x)。晶格位点的所有状态的组合被表示为f(x)。

状态的数量由每个能级内可能的速度向量的数量来确定。速度向 量由具有三个维度x、y和z的空间中的整数线性速率组成。对于多 种属模拟，状态的数量增加。

每个状态i表示处于特定能级(即，能级零、一或二)的不同速 度向量。每个状态的速度c_i利用其在三个维度当中的每一个维度中的 “速率”指示如下：

c_i＝(c_ix,c_iy,c_iz,). 方程(I5)

能级零状态表示在任何维度都不移动的停止的粒子，即，c_stopped＝(0,0,0)。能级一状态表示在三个维度之一中具有±1速率并且在其 它两个维度中具有零速率的粒子。能级二状态表示在所有三个维度中 都具有±1速率、或者在三个维度之一中具有±2速率并且在其它两 个维度中具有零速率的粒子。

生成三个能级的所有可能的排列给出总共39个可能的状态(一 个能量零状态、6个能量一状态、8个能量三状态、6个能量四状态、 12个能量八状态和6个能量九状态)。

每个体元(即，每个晶格位点)由状态向量f(x)表示。该状态向 量完全定义体元的状态并且包括39个条目。这39个条目对应于一个 能量零状态、6个能量一状态、8个能量三状态、6个能量四状态、 12个能量八状态和6个能量九状态。通过使用这个速度集，系统可以对实现的平衡状态向量产生麦克斯韦-玻尔兹曼统计。对于总能量， 对于第二分配q_i(x,t)使用相同数量的项(39个状态)来表示特定总能 量。连同状态向量f(x)一起，这被用于求解总能量守恒。

为了处理效率，体元被分组在被称为微观块的2x2x2的体积中。 微观块被组织成允许体元的并行处理并且最小化与数据结构关联的开 销。微观块中用于体元的速记符号被定义为N_i(n)，其中n表示微观 块中晶格位点的相对位置并且n{0,1,2,...,7}。微观块在图10中 示出。

参考图11A和11B，表面S(图11A)在模拟空间(图11B)中 被表示为面元F_α的集合：

S＝{F_α} 方程(I6)

其中α是列举特定面元的索引。面元不限于体元的边界，但是通常具 有量级为与该面元相邻的体元的尺寸的尺寸，或者具有稍小于与该面 元相邻的体元的尺寸的尺寸，以使得面元影响相对少量的体元。为了 实现表面动力学，向面元分配属性。具体地，每个面元F_α具有单位 法线(n_α)、表面积(A_α)、中心位置(x_α)和描述面元的表面动 态属性的面元分布函数(f_i(α))。总能量分布函数q_i(α)被以用于面 元的流分布和体元交互相同的方式被处理。

参考图12，不同的分辨率水平可以在模拟空间的不同区域中使 用，以提高处理效率。通常，对象655周围的区域650是最关心的并 且因此利用最高分辨率进行模拟。因为粘度的影响随着离对象的距离 而减小，所以采用降低的分辨率水平(即，扩大的体元体积)来模拟 在离对象655按增加的距离隔开的区域660，665。类似地，如图13 中所示，较低的分辨率水平可以被用来模拟对象775的较不显著特征 周围的区域770，而最高的分辨率水平被用来模拟对象775的最显著 特征(例如，前沿和后缘表面)周围的区域780。边远区域785利用 最低的分辨率水平和最大的体元来模拟。

C.识别受面元(facet)影响的体元

再次参考图9，一旦模拟空间已经被建模(步骤302)，受一个 或多个面元影响的体元就被识别(步骤304)。体元可以以多种方式 受面元影响。首先，被一个或多个面元相交的体元受影响在于：该体 元相对于非相交的体元具有减小的体积。这会发生是因为面元以及在 由该面元表示的表面下面的材料占据了体元的一部分。分数因子P_f(x) 指示体元的不受面元影响的部分(即，可以被流体或为对其模拟流的 其它材料占据的部分)。对于非相交体元，P_f(x)等于1。

通过将粒子传送到面元或者从面元接收粒子而与一个或多个面元 相交的体元也被识别为受面元影响的体元。被面元相交的所有体元都 将包括从面元接收粒子的至少一个状态以及向面元传送粒子的至少一 个状态。在大多数情况下，附加的体元也将包括这种状态。

参考图14，对于具有非零速度向量c_i的每个状态i，面元F_α从 由平行六面体G_iα定义的区域接收粒子，或向其传送粒子，其中平行 六面体G_iα具有由速度向量c_i和面元(|c_in_i|)的单位法线n_α的向量点 积的量值定义的高度以及由面元的表面积A_α定义的基部，使得平行 六面体G_iα的体积V_iα等于：

V_ia＝|c_in_a|A_a 方程(I7)

当状态的速度向量指向面元时(|c_i n_i|<0)，面元F_α从体积V_iα 接收粒子，并且当状态的速度向量指向远离面元时(|c_i n_i|>0)，向 该区域传送粒子。如下面将要讨论的，当另一个面元占据平行六面体 G_iα的一部分时，即，在诸如内角的非凸特征部附近会发生的一种状 况，这个表达式被修改。

面元F_α的平行六面体G_iα可以重叠多个体元的部分或全部。体元 的数量或其中的部分依赖于相对于体元尺寸的面元尺寸、状态的能量 以及面元相对于晶格结构的朝向。受影响的体元的数量随着面元的尺 寸而增大。因此，如上面所指出的，面元的尺寸通常被选择为位于该 面元附近的体元的尺寸的量级或小于位于该面元附近的体元的尺寸。

被平行六面体G_iα重叠的体元N(x)的部分被定义为V_iα(x)。利用 这个项，在体元V_iα(x)和面元F_α之间移动的状态i粒子的通量Γ_iα(x) 等于该体元中状态i粒子的密度(N_i(x))乘以与该体元重叠的区域的 体积(V_iα(x))：

Γ_ia(x)＝N_i(x)+V_ia(x) 方程(I8)

当平行六面体G_iα被一个或多个面元相交时，以下条件为真：

V_ia＝∑V_a(x)+∑V_ia(β) 方程(I9)

其中第一个求和考虑了被G_iα重叠的所有体元并且第二项考虑了 与G_iα相交的所有面元。当平行六面体G_iα不被另一面元相交时，这 个表达式简化为：

V_ia＝∑V_ia(x) 方程(I10)

D.执行模拟

一旦受一个或多个面元影响的体元被识别(步骤304)，定时器 就被初始化，以开始模拟(步骤306)。在模拟的每个时间增量期间， 粒子从体元到体元的移动由说明粒子与表面面元的相互作用的平流阶 段模拟(步骤308-316)。接下来，碰撞阶段(步骤318)模拟每个 体元内粒子的相互作用。其后，定时器递增(步骤320)。如果递增 后的定时器不指示模拟完成(步骤322)，则平流阶段和碰撞阶段 (步骤308-320)重复。如果递增后的计时器指示模拟已完成(步骤 322)，则模拟的结果被存储和/或显示(步骤324)。

1.用于表面的边界条件

为了正确模拟与表面的相互作用，每个面元满足四个边界条件。 首先，由面元接收的粒子的组合质量等于由面元传送的粒子的组合质 量(即，到面元的净质量通量等于零)。第二，由面元接收的粒子的 组合能量等于由面元传送的粒子的组合能量(即，到面元的净能量通 量等于零)。这两个条件可以通过要求在每个能级(即，能级一和能 级二)的净质量通量等于零来满足。

其它两个边界条件与和面元相互作用的粒子的净动量相关。对于 没有皮肤摩擦的表面(在本文被称为滑动表面)，净切向动量通量等 于零并且净法线动量通量等于在面元处的局部压力。因此，与面元的 法线n_α垂直的组合接收和传送的动量的分量(即，切向分量)相等， 而与面元的法线n_α平行的组合接收和传送的动量的分量(即，法线 分量)之间的差等于在该面元的局部压力。对于非滑动表面，表面的 摩擦相对于由面元接收的粒子的组合切向动量将由面元传送的粒子的 组合切向动量减小了与摩擦量相关的因子。

2.从体元到面元的聚集

作为模拟粒子和表面之间的相互作用的第一步，粒子从体元被聚 集并被提供给面元(步骤308)。如以上所指出的，体元N(x)和面元 F_α之间的状态i粒子的通量为：

Γ_iα(x)＝N_i(x)V_iα(x) 方程(I11)

由此看来，对于指向面元F_α的每个状态i(c_in_α<0)，由体元提 供给面元F_α的粒子的数量为：

Γ_iαV→F＝Σ_XΓ_iα(x)＝Σ_X N_i(x)V_iα(x) 方程(I12)

只有其V_iα(x)具有非零值的体元需要被求和。如以上所指出的， 面元的尺寸被选择为使得V_iα(x)仅对于少数体元具有非零值。因为 V_iα(x)和P_f(x)可以具有非整数值，所以Γ_α(x)作为实数被存储和处理。

3.从面元到面元的移动

接下来，粒子在面元之间被移动(步骤310)。如果用于面元F_α 的传入状态(c_in_α<0)的平行六面体G_iα被另一面元F_β相交，则由F_α 接收的状态i粒子的一部分将来自面元F_β。特别地，面元F_α将接收 在前一时间递增期间由面元F_β产生的状态i粒子的部分。这种关系 在图17中示出，其中被面元F_β相交的平行六面体G_iα的部分1000等 于被面元F_α相交的平行六面体G_iβ的部分1005。如以上所指出的， 相交的部分被表示为V_iα(β)。利用这个项，面元F_β与面元F_α之间的 状态i粒子的通量可以被描述为：

Γ_iα(β，t-1)＝Γ_i(β)V_iα(β)/V_iα 方程(I.13)

其中Γ_i(β,t-1)是在前一时间递增期间由面元F_β产生的状态i粒子的测 量。由此看来，对于指向面元F_α的每个状态i(c_in_α<0)，由其它面 元提供给面元F_α的粒子的数量为：

Γ_iαF→F＝Σ_βΓ_iα(β)＝Σ_βΓ_i(β，t-1)V_iα(β)/V_iα 方程(I.14)

并且到面元中的状态i粒子的总通量是：

Γ_iIN(α）＝Γ_iαF→F+Γ_iαF→F＝Σ_xN_i(x)V_iα+Σ_βΓ_i(β，t-1)V_iα(β)/V_iα 方程(I.15)

用于面元的状态向量N(α)(也被称为面元分布函数)具有对应 于体元状态向量的M个条目的M个条目。M是离散晶格速率的数量。 面元分布函数N(α)的输入状态被设置为等于进入那些状态中的粒子 的通量除以体积V_iα：

对于c_i n_α<0，N_i(α)＝Γ_iIN(α)/V_iα 方程(I.16)

面元分布函数是用于从面元生成输出通量的模拟工具，并且不一 定表示实际的粒子。为了生成准确的输出通量，值被分配给分布函数 的其它状态。向外的状态是利用上述用于填充向内的状态的技术来填 充的：

对于c_i n_α≥0，N_i(α)＝Γ_iOTHER(α)/V 方程(I.17)

其中Γ_iOTHER(α)是利用上述用于生成Γ_iIN(α)的技术来确定的，但是将 该技术应用到除传入状态(c_i n_α<0)之外的状态(c_i n_α≥0)。在替代 的方法中，Γ_iOTHER(α)可以利用来自前一时间步长的Γ_iOTHER(α)的值生 成，使得：

Γ_iOTHER(α,t)＝Γ_iOUT(α,t-1) 方程(I.18)

对于平行状态(c_i n_α＝0)，V_iα和V_iα(x)都是零。在用于N_i(α)的表 达式中，V_iα(x)出现在分子中(根据用于Γ_iOTHER(α)的表达式)且V_iα 出现在分母中(根据用于N_i(α)的表达式)。因此，当V_iα和V_iα(x)接 近零时，用于并行状态的N_i(α)被确定为N_i(α)的极限。具有零速度的 状态(即，静止状态和状态(0,0,0,2)和(0,0,0,-2))的值在模拟的开 始基于针对温度和压力的初始条件被初始化。然后，这些值随时间被 调整。

4.执行面元表面动力学

接下来，对每个面元执行表面动力学，以满足以上讨论的四个边 界条件(步骤312)。用于为面元执行表面动力学的过程在图18中 示出。最初，通过确定粒子在面元处的组合动量P(α)来确定到面元 F_α的组合动量(步骤1105)：

对于所有i，

根据这个方程，法线动量P_n(α)被确定为：

P_n(α)＝n_α·P(α) 方程(I.20)

然后，这个法线动量利用推/拉技术被消除(步骤1110)，以产 生N_n-(α)。根据这种技术，粒子以只影响法线动量的方式在状态之间 被移动。推/拉技术在美国专利No.5,594,671中被描述，其通过引用 被结合于此。

其后，N_n-(α)的粒子被碰撞，以产生玻尔兹曼分布N_n-β(α)(步骤 1115)。如下面关于执行流体动力学所描述的，玻尔兹曼分布可以通 过向N_n-(α)应用一组碰撞规则来实现。

然后，基于传入的通量分布和玻尔兹曼分布，确定用于面元F_α 的传出通量分布(步骤1120)。首先，传入的通量分布Γ_i(α)和玻尔 兹曼分布之间的差被确定为：

ΔΓ_i(α)＝Γ_iIN(α)-N_n-βi(α)V_iα 方程(I.21)

使用该差，传出的通量分布是：

对于n_αc_i>0，Γ_iOUT(α)＝N_n-βi(α)V_iα-.Δ.Γ_i*(α) 方程(I.22)

并且其中i*是与状态i具有相反方向的状态。例如，如果状态i是(1, 1,0,0)，则状态i*是(-1,-1,0,0)。为了说明皮肤摩擦及其它因素，传 出通量分布可以进一步细化为：

对于n_αc_i>0，

其中C_f是皮肤摩擦的函数，t_iα是与n_α垂直的第一切线向量，t_2α是与 n_α和t_1α都垂直的第二切线向量，并且ΔN_j,1和ΔN_j,2是对应于状态i和 所指示的切线向量的能量(j)的分布函数。该分布函数根据下式确定：

其中j对于能级1状态等于1并且对于能级2状态等于2。

用于Γ_iOUT(α)的方程的每一项的功能如下。第一项和第二项实施 法线动量通量边界条件至如下程度：碰撞已有效产生玻尔兹曼分布， 但包括切向动量通量异常。第四项和第五项校正这种异常，这种异常 可能由于离散效应或因碰撞不足造成的非玻尔兹曼结构引起。最后， 第三项添加指定量的皮肤摩擦，以实施在表面上的切线动量通量的期 望改变。摩擦系数C_f的生成在下面被描述。需要注意的是，涉及向 量操纵的所有项都是可以在开始模拟之前被计算的几何因子。

由此，切向速度被确定为：

u_i(α)＝(P(α)-P_n(α)n_α)/ρ， 方程(I.25)

其中ρ是面元分布的密度：

和以前一样，传入的通量分布与玻尔兹曼分布之差被确定为：

ΔΓ_i(α)＝Γ_iIN(α)-N_n-βi(α)V_iα 方程(I.27)

于是，传出通量分布变为：

Γ_iOUT(α)＝N_n-βi(α)V_iα-ΔΓ_i*(α)+C_f(n_αc_i)[N_n-βi*(α)-N_n-βi(α)]V_iα 方程(I.28)

其对应于由之前的技术确定的传出通量分布的前两行，但不需要对异 常切向通量的校正。

使用任一种方法，结果产生的通量分布都满足所有的动量通量条 件，即：

其中，p_α是在面元F_α处的平衡压力并且基于向该面元提供粒子的体 元的平均密度和温度值，并且u_α是在该面元处的平均速度。

为了确保质量和能量边界条件被满足，对每个能级j测量输入能 量与输出能量之差：

其中索引j表示状态i的能量。然后，这个能量差被用来生成差异项：

对于c_jin_α>0，

这个差异项被用来修正传出通量，使得该通量变为：

对于c_jin_α>0，Γ_αjiOUTf＝Γ_αjiOUT+δΓ_αji 方程(I.32)

这种操作校正质量和能量通量，同时保留切向动量通量未被更改。如 果流在面元的邻域是大致均匀的并且接近平衡，则这种调整较小。在 调整之后，基于邻域平均属性加上由于邻域的非均匀性或非平衡属性 造成的校正，结果产生的法线动量通量稍微被更改为是平衡压力的值。

5.从体元到体元的移动

再次参考图9，粒子沿三维直线晶格在体元之间移动(步骤 314)。这种体元到体元的移动是对不与面元相互作用的体元(即， 不位于表面附近的体元)执行的唯一移动操作。在典型的模拟中，不 位于足够接近表面以便与表面相互作用的体元构成体元的大多数。

每个分开的状态表示在三个维度x、y和z当中的每一个中以整 数速率沿晶格移动的粒子。整数速率包括：0，±1和±2。速率的符 号指示粒子正在沿对应的轴移动的方向。

对于不与表面相互作用的体元，移动操作在计算上是非常简单的。 状态的整个布居(population)在每个时间增量期间从其目前的体元 移动到目的地体元。同时，目的地体元的粒子从那个体元移动到其自 己的目的地体元。例如，在+1x和+1y方向(1,0,0)上移动的能级1粒 子从其当前的体元移动到在x方向是+1并且对其它方向是0的体元。 粒子以其在移动之前具有的状态(1,0,0)终止在其目的地体元。基于 与其它粒子和表面的局部相互作用，体元内的相互作用将有可能改变 用于那个状态的粒子计数。如果不是，则粒子将继续以相同的速率和 方向沿晶格移动。

对于与一个或多个表面相互作用的体元，移动操作变得稍微更复 杂。这会导致一个或多个分数粒子被传送到面元。这种分数粒子到面 元的传送导致分数粒子留在体元中。这些分数粒子被传送到被面元占 据的体元中。

还示出了对平流的总能量修正315，其中计算机通过将特定的总 能量值添加到平流的粒子的状态并从未在时间间隔内平流的粒子的状 态中减去特定的总能量值来修正粒子的状态向量。

参考图15，当体元905的状态i粒子的一部分900被移动到面元 910时(步骤308)，剩余的部分915被移动到面元910在其中定位 并且状态i粒子从其指向面元910的体元920。因此，如果状态布居 等于25并且V_iα(x)等于0.25(即，体元的四分之一与平行六面体G_iα 相交)，则6.25个粒子将被移动到面元F_α并且18.75个粒子将被移 动到被面元F_α占据的体元。因为多个面元会与单个体元相交，所以 传送到被一个或多个面元占据的体元N(f)的状态i粒子的数量为：

其中N(x)是源体元。

6.从面元到体元的分散

接下来，来自每个面元的传出粒子被分散到体元(步骤316)。 本质上，这一步是将粒子从体元移动到面元的聚集步骤的逆。从F_α 面元移动到体元N(x)的状态i粒子的数量是：

其中P_f(x)说明部分体元的体积减小。由此看来，对于每个状态i，从 面元指向体元N_(x)的粒子的总数为：

在从面元到体元分散粒子之后，将粒子与已经从周围的体元平流 输送的粒子组合，并且将结果整数化，有可能某些体元中的某些方向 可能下溢(变成负)或者上溢(在八位实现方式中超过255)。这将 导致在质量、动量和能量被截断以适应允许的值范围之后这些量的增 益或损失。

为了防止这种情况的发生，出界的质量、动量和能量在违规状态 的截断之前被累加。对于状态所属的能量，等于增益(由于下溢)或 损失(由于上溢)的值的质量的量被加回到随机(或顺序)选择的具 有相同能量并且本身不受制于上溢或下溢的状态。从这种质量和能量 的添加而造成的额外动量被累加并被加到来自截断的动量。通过仅将 质量添加到相同的能量状态，当质量计数器达到零时，质量和能量都 被校正。最后，动量利用推/拉技术被校正，直到动量累加器返回到 零。

7.执行流体动力学

执行流体动力学(步骤318)图9。这个步骤可以被称为微动力 学或体元内操作。类似地，平流过程可以称为体元间操作。以下描述 的微动力学操作也可以被用来在面元碰撞粒子，以产生玻尔兹曼分布。

流体动力学在晶格玻尔兹曼方程模型中通过已知为BGK碰撞模 型的特定碰撞算子来确保。这种碰撞模型模仿真实流体系统中的分布 的动力学。碰撞过程可以通过方程1和方程2的右手侧很好地得到描 述。在平流步骤之后，利用方程3从分布函数获得流体系统的守恒量， 具体而言是密度、动量和能量。从这些量，由方程(2)中的f^eq指出 的平衡分布函数完全由方程(4)规定。速度向量集c_i、权重的选择 都在表1中列出，连同方程2一起确保宏观行为遵循正确的流体动力 学方程。

E.可变分辨率

参考图12，可变分辨率(如以上所讨论的)采用不同大小的体 元，在下文中称为粗糙体元12000和精细体元1205。(以下讨论涉 及具有两种不同大小的体元；应该认识到的是，可以将所描述的技术 应用于三种或更多种不同大小的体元以提供更高的分辨率水平)。粗 糙体元和精细体元的区域之间的接口被称为可变分辨率(VR)接口 1210。

当在表面或表面附近使用可变分辨率时，面元可以与VR接口两 侧的体元相互作用。这些面元被分类为VR接口面元1215(F_αIC)或 VR细面元1220(F_αIF)。VR接口面元1215是位于VR接口的粗糙侧 上并且具有延伸到精细体元中的粗糙平行六面体1225的面元。(粗 糙平行六面体是根据粗糙体元的维度来确定c_i的维度，而精细平行六 面体是根据精细体元的维度来确定c_i的维度)。VR精细面元1220 是定位在VR接口的精细侧上并且具有延伸到粗糙体元的精细平行六 面体1230中的面元。与接口面元相关的处理还可以涉及与粗糙面元1235(F_αC)和精细面元1240(F_αF)的相互作用。

对于这两种类型的VR面元，都以精细的比例执行表面动力学， 并且如上所述进行操作。但是，就粒子向VR面元和从VR面平流的 方式而言，VR面元与其它面元不同。

使用图13中所示的可变分辨率过程1300来处理与VR面元的交 互。此过程的大多数步骤是使用上面针对与非VR面元的交互而讨论 的同等步骤执行的。在粗糙时间步长(即，对应于粗糙体元的时间段) 期间执行过程1300，该粗糙时间步长包括两个阶段，每个阶段对应 于精细时间步长。在每个精细时间步长期间执行面元表面动力学。因 此，VR接口面元F_αIC被认为是两个尺寸相同且朝向相同的精细面元， 分别称为黑色面元F_αICb和红色面元F_αIC。黑色面元F_αICb与粗糙时间 步长内的第一精细时间步长相关联，而红色面元F_αICr与粗糙时间步 长内的第二精细时间步长相关联。

最初，通过第一表面到表面平流阶段，粒子在面元之间移动(平 流)(步骤1302)。粒子以权重因子V_-αβ从黑色面元F_αICb移动到粗 糙面元F_βC，该权重因子对应于从面元F_α延伸出并位于面元F_β后面 的粗糙平行六面体(图12、1225)的未阻塞部分的体积减去从面元 F_α延伸出并位于面元F_β后面的精细平行六面体的未阻塞部分(图12、 1245)。精细体元的c_i的量值是粗糙体元的c_i的量值的一半。如上所 述，针对面元F_α的平行六面体的体积定义为：

V_iα＝|c_in_α|A_α 方程(I.36)

因此，由于体元的表面积A_α在粗糙平行六面体和精细平行六面 体之间不变化，并且由于单位法线n_α的量值始终为1，因此对应于面 元的精细六面体的体积为该面元的对应粗糙平行六面体的体积的一半。

粒子以权重因子V_αβ从粗糙面元F_αC移至黑色面元F_βIC，该权重 因子V_αβ对应于从面元F_α延伸出并位于面元F_β后面的精细平行六面 体的未阻塞部分的体积。

粒子以权重因子V_αβ从红色面元F_αICr移动到粗糙面元F_βC，并且 以权重因子V_-αβ从粗糙面元F_αC移动到红色面元F_βICr。

粒子以权重因子V_αβ从红色面元F_αICr移动到黑色面元F_βICb。在 此阶段，不会发生黑到红平流。此外，由于黑色和红色面元代表连续 的时间步长，所以绝不会发生黑到黑平流(或红到红平流)。出于类 似的原因，该阶段中的粒子以权重因子V_αβ从红色面元F_αICr移动到精 细面元F_βIF或F_βF，并且以相同的权重因子从精细面元F_αIF或F_αF移 动到黑色面元F_αICb。

最后，粒子以相同的权重因子从精细面元F_αIF或F_αF移动到其他 细面元F_βIF或F_βF，然后以权重因子V_Cαβ的粗糙面元F_αC移到其他粗 糙面元F_C，其权重因子V_Cαβ对应于从面元F_α延伸出并位于面元F_β 后面的粗糙平行六面体的未阻塞部分的体积。

在粒子在表面之间平流之后，在第一聚集阶段从体元聚集粒子 (步骤1304-1310)。对于精细面元F_αF使用精细平行六面体从精细 体元聚集粒子(步骤1304)，并对于粗糙面元F_αC使用粗糙平行六面 体从粗糙体元聚集粒子(步骤1306)。然后对于黑色面元F_αIRb以及 对于VR精细面元F_αIF使用精细平行六面体从粗糙和精细体元两者聚 集粒子(步骤1308)。最终，利用粗糙平行六面体和精细平行六面 体之间的差异，对于红色面元F_αIRr从粗糙体元聚集粒子(步骤 1310)。

接下来，将与精细体元或VR面元交互的粗糙体元分解为精细体 元的集合(步骤1312)。将在单个粗糙时间步长内将粒子传输到细 细体元的粗糙体元的状态被分解。例如，未被面元相交的粗糙体元的 适当状态被分解为如图4的微块那样被定向的八个精细体元。被一个 或多个面元相交的粗糙体元的适当状态被分解为与未被任何面元相交 的粗糙体元部分相对应的完全和/或部分精细体元的集合。粗糙体元 和由其分解得到的精细体元的粒子密度N_i(x)相等，但是精细体元的 分数因子P_f可以不同于粗糙体元的分数因子，也可以不同于其他体 元的分数因子。

此后，对精细面元F_αIF和F_αF执行表面动力学(步骤1314)，并 对黑色面元F_αICb执行表面动力学(步骤1316)。使用图11所示的 和上面讨论的过程进行动力学。

接下来，使粒子在包括实际精细体元和由粗糙体元分解得到的精 细体元的精细体元之间移动(步骤1318)。一旦粒子已经移动，粒 子就从精细面元F_αIF和F_αF分散到精细体元(步骤1320)。

颗粒也从黑色面元F_αICb分散到精细体元(包括由粗糙体元分解 得到的精细体元)(步骤1322)。如果体元在当时没有表面存在的 情况下接收到粒子，则粒子会分散到精细体元。特别地，当体元是实 际的精细体元(与由粗糙体元分解得到的精细体元相反)时，当超出 体元N(x)1个速度单位的体元N(x+c_i)为实际的精细体元时，或者当超 出体元N(x)1个速度单位的体元N(x+c_i)为由粗糙体元分解得到的精细 体元时，粒子会分散到体元N(x)。

最后，通过对于精细体元执行流体动力学来完成第一精细时间步 长(步骤1324)。对其执行流体动力学的体元不包括由粗糙体元分 解而产生的精细体元(步骤1312)。

过程1300在第二精细时间步长期间实施类似的步骤。最初，在 第二表面到表面平流阶段中，使颗粒在表面之间移动(步骤1326)。 粒子从黑色面元平流到红色面元，从黑色面元平流到精细面元，从精 细面元平流到红色面元，以及从精细面元平流到精细面元。

在粒子在表面之间平流之后，在第二聚集阶段从体元聚集粒子 (步骤1328-1330)。对于红色面元F_αIRr使用精细平行六面体从精细 体元聚集粒子(步骤1328)。还对于精细面元F_αF和F_αIF使用精细平 行六面体从精细体元聚集粒子(步骤1330)。

此后，如上所述，对于精细面元F_αIF和F_αF(步骤1332)，对于 粗糙面元F_αC(步骤1134)，以及对于红色面元F_αICr(步骤1336) 执行表面动力学。

接下来，使用精细分辨率在体元之间移动粒子(步骤1338)， 使得粒子在精细体元和代表粗糙体元的精细体元之间来回移动。然后， 使用粗糙分辨率在体元之间移动粒子(步骤1340)，使得粒子在粗 糙体元之间来回移动。

接下来，在组合步骤中，将粒子从面元分散到体元，而将代表粗 糙体元的精细体元(即，由粗糙体元分解得到的精细体元)凝聚为粗 糙体元(步骤1342)。在此组合步骤中，粒子通过使用粗糙平行六 面体从粗糙面元分散到粗糙体元，使用精细平行六面体从精细面元分 散到精细体元，使用精细平行六面体从红色面元分散到精细或粗糙体 元，并使用粗糙平行六面体和精细平行六面体之间的差异从黑色面元 分散到粗糙体元。最后，对于精细体元和粗糙体元执行流体动力学 (步骤1344)。

F.标量传输求解器

如上所述，可以将各种类型的LBM应用于求解流体流动，流体 流动作为标量传输的背景载体。在模拟过程中，流体流和标量传输两 者都被模拟。例如，流体流是使用晶格玻尔兹曼(LB)方法模拟的 流，并且使用在本文中被称为标量传输方程的一组分布函数模拟标量 传输。

虽然此处提供了用于模拟流体流动的LBM方法的详细描述，但 以下是可与标量模拟结合使用的一种模拟流体流动的方法的示例：

其中f_i(x,t)(i＝1,...,19)是粒子分布函数，τ是单个弛豫时间，

是具有三阶流体速度展开的平衡分布函数，

其中T₀＝1/3。离散晶格速度c_i为：

对于静止粒子，w₀＝1/3，对于笛卡尔方向的状态，w_i＝1/18，对 于双对角线方向的状态，w_i＝1/36。g_i(x,t)是外部体积力项。流体动 力量和u是粒子分布函数的矩：

如上所述，流体求解器与生成标量传输信息的标量传输求解器结 合使用。因此，除了流体求解器外，还引入了一组单独的分布函数 T_i用于标量传输。因此，对于系统中的每个体元，系统模拟流体流和 标量传输两者，以生成代表流体流的状态矢量和代表标量变量的标量 数量。这些模拟结果作为条目存储在计算机可访问的介质中。

标量传输函数集提供了求解二阶宏观标量传输方程的间接方式。 T_i提供了对标量数量的动态演化进行建模的方程：

T_i是标量分布函数，T是要求解的标量。τ_T是对应于标量扩散 率的弛豫时间。弛豫时间提供了系统弛豫到平衡所需时间的度量。分 别在方程(38)，(39)和(41)中定义了项f_i,ρ,T₀以及u。

用c_i表示的晶格速度集可以是用于标量模拟的离散晶格速率集。 通常，用于标量分布的晶格速率集不必与用于流体分布的晶格速率集 相同，这是因为标量求解器是附加到基本流体求解器的附加系统。例 如，在标量演化的模拟中可以使用比在流体流动的模拟中更少的晶格 速率。只要标量晶格速率集是流体晶格速率集的子集，就可以对于标 量应用不同的标量晶格速率集。例如，当将19速率LB模型用于流 体模拟时，可以将6速率LB模型用于标量模拟。由于19速率LB 模型比6速率LB模型具有更高阶的晶格对称性，因此在下面提供的 示例中使用了相同的用于标量的19速率模型。

标准的众所周知的BGK(例如，如上所述)包括所有阶的非平 衡矩。认为包括所有非平衡矩不必然是各向同性的、流体动力学的、 或物理上有意义的。因此，使用BGK正则化/过滤碰撞算子形式。碰 撞算子Φ_i代表未来的碰撞因子。该碰撞算子仅在相关的受支持阶次 (例如，仅第一阶)提取非平衡标量属性。在与不支持/不希望的高 阶模式相关联的不平衡属性被去除的同时，算子还保留和释放感兴趣 的模式。此投影足以恢复标量传输物理学(例如，对流和扩散)。相 信使用这种未来碰撞算子可以显著降低噪声，提供更好的平流性能 (例如，可以更加伽利略不变)，并且与公知BGK算子的其他解相 比，被认为比其他解更稳定。这样的形式(例如，如方程43所示) 确保了仅一阶不平衡矩对于流体动力范围内的标量扩散有贡献。通过 该碰撞过程，可以滤除所有更高阶的不平衡矩。相信如上所述使用碰撞算子Φ_i提供了益处，包括消除了在BGK中表现出的数字噪声和改 善鲁棒性。标量T充当其自身的平衡，并且不需要标量平衡分布函 数的复杂表达式。碰撞算子Φ_i的总体计算相当高效。相信过滤更高阶 的不平衡矩可以额外地提供减少在更高阶平衡解中可能存在的混叠的 优点。

可以示出(42)中的碰撞遵守了标量守恒定律。在方程(42)的 两侧乘以f_i′(x,t)＝f_i(x+c_i,t+1)并注意：

结果是

其中T_i′(x,t)表示方程(42)的右侧。因此，标量碰撞算子保留 局部ρT，这意味着在标量被视为温度的情况下实现了局部能量守恒。 由于T_i与f_i一起传播，所以在平流过程中能量分布E＝f_iT_i被完全保持。 因此实现了ρT的全局守恒。此外，最值得注意的是，该方案保持了 标量T均匀性的精确不变性。这样直观地看出，如果在各处

则无论背景流场如何，在以后的所有时间 在各处，φ(x,t)＝0且

以前的任何晶格玻尔兹曼标 量模型都没有展示这种基本属性。

使用Chapmann-Enskog展开，可以示出等式(42)恢复了以下 二阶宏观标量传输方程：

其中κ＝(τ_T-1/2)T₀。均匀性不变条件确保ρ在

之外。

边界条件

LBM的一个显著优势是能够处理复杂的几何形状。在用于在任 意几何形状上实现无摩擦/有摩擦边界条件(BC)的广义体积LB表 面算法中，质量守恒，并且边界上的切向和法向动量通量被精确实现。 局部细致平衡已完全满足。作为这种方法的直接扩展，可得出在用于 标量的任意几何形状上的绝热(零标量通量)BC。一旦实现了绝热 BC，就可以实现规定的有限通量BC。

与其他逐点LB不同，边界条件是在离散化的表面元素集上实行 的。这些分段平面元件一起表示弯曲的几何形状。在粒子平流期间， 如图14所示，每个表面元件收集来自其相邻的流体单元的传入颗粒 (步骤1410)。通过来自下层表面元素的平行六面体在粒子移动方 向上与单元的体积重叠来对传入分布

进行加权(步骤1412)。 在接收到输入量之后，应用以下的表面标量算法(步骤1414)。未 来确定表面的传出分布，

其中，

P_i(≡|n·c_i|A)

c_i·nn<0,c_in＝-c_i·P_i＝P_i·n

这里n是指向流体域的表面法线，并且c_i·nn<0,c_in＝-c_i。 P_i(≡|n·c_i|A)是与表面法线n和给定表面元素的面积A相关联的粒 子方向c_i上的平行六面体的体积，显然P_i＝P_i*。根据相同的表面平 流过程，传出分布从表面元素传播回流体单元(步骤1416)。不难 证明，上述表面标量碰撞达到了精确的零表面标量通量。在传出方向 上求和，传出标量通量为：

注意P_i＝P_i*以及方程(49)中Tⁱⁿ的定义，右手侧的第二个求和 项为零。另外，由于质量通量守恒，

总的传出标量通量与总的传入标量通量相同：

因此，在任意几何形状上完全满足零净表面通量(绝热)BC。 如果在表面上指定了外部标量源Q(t)，则可以将源项直接添加到 公式(48)

如果边界条件具有规定的标量数量T_w(例如表面温度)，则可 以相应地计算表面热通量：

Q(t)＝ρC_pκ(T_ω-Tⁱⁿ). 方程(I.52)

数值验证

图15-18显示了四组模拟结果，这些结果证明了LB标量求解器 的在其数值精确性、稳定性、伽利略不变性、网格朝向独立性等方面 的能力。作为比较，还呈现使用两种不同的二阶FD方案，van Leer 型通量限制器方案的和直接混合方案(中央和一阶迎风方案的混合) 的结果。

A.剪切波衰减

第一个测试案例是恒定的均匀流体流携带的温度剪切波衰减。初 始温度分布是均匀的，加上晶格波长L＝16，幅值

的空间正弦变化。T_A是常数。背景平均流的速度 为0.2，热扩散系数κ为0.002，如此低的分辨率和κ可以很好地验证 数值稳定性和精确性。对于没有背景流的温度衰减，LB标量求解器 和有限差分法显示出与理论极佳的一致性。在非零背景平均流的情况 下，LB标量求解器仍能够精确地与理论相比。然而，FD结果显示 出明显的数值误差；在图15中绘制了在晶格时间步长81的温度曲线。 对于通量限制器FD方案，可以清楚地看到数值扩散，而混合FD方 案既不能保持正确的温度曲线，也不能保持其位置。

B.带有热源的倾斜通道

第二个测试案例是模拟具有不同晶格取向的通道流中的温度分布。 通道壁是自由滑动的，因此流体流保持均匀，作为结果U₀＝0.2。通 道宽度为50(晶格间距)，流量Re为2000。热扩散系数κ为0.005。 墙壁上的温度在T_W＝1/3时固定。在体流体域中应用了恒定体积热 源q＝5×10^-6。流在流向上具有周期性边界条件，这在晶格对齐情 况下容易实现。当通道(浅色)倾斜如图16所示地倾斜时，输入和 输出通道边界在坐标方向上完美匹配，从而在流向上再次实现了周期 性边界条件。为了证明晶格独立性，我们选择倾斜角度为26.56度。 像第一个测试案例一样，当通道被晶格对齐时，使用LB标量求解器 和两种FD方案的温度分布与解析解非常匹配。但是，当通道倾斜时， FD方案的结果与理论明显背离。跨倾斜通道的温度分布的模拟结果 在图16中示出.LB结果清楚地示出为与晶格方向无关。FD方法的误差也源于其本质上难以在倾斜边界方向上进行梯度计算，因此引入 了额外的数值伪像。由于LB标量粒子平流因此处介绍的BC是精确 的，因此它能够实现独立于晶格取向的标量演化。

C.倾斜通道中的温度传播

由于在笛卡尔网格上的非晶格对齐的近墙壁区域中缺少邻居信息， 因此很难准确估计局部梯度，这对于基于FD的方法至关重要。此外， 由于强烈依赖于迎风信息，标量平流的计算对于FD方法可能会进一 步被妥协。相反，如上所述，LB标量求解器的边界处理实现了精确 的局部标量通量守恒。这样的近壁区域中的标量平流可以被精确地计 算。作为示例，在倾斜30度的通道中进行了高温对流。自由滑动和 绝热BCs在实体壁上被实施，并且流体速度沿着通道是恒定的U₀＝0.0909。热扩散系数κ为0.002。最初，除入口处的T＝4/9以外，其 他所有地方的温度均为1/3。然后，该温度前沿应由于均匀背景流体 流而在晚些时候在无变形的情况下对流到下游位置。在图17中示出 了计算出的在晶格时间步长2000的跨通道的温度前沿分布。LB标 量求解器的温度前沿保持近似笔直的曲线。另一方面，这两个FD方 案的温度前沿在近壁区域显示出明显的变形。还值得一提的是，LB 标量结果显示出最薄的温度前沿，这意味着LB标量求解器的数值扩 散较小。

D.瑞利伯纳德(Rayleigh-Bernard)自然对流

瑞利伯纳德(RB)自然对流是数值求解器精确性验证的经典基 准。它具有简单的案例设置，但具有复杂的物理现象。当瑞利数Ra 超过某个临界值时，系统会经历从无流动到对流的转变。如图18所 示，当前的研究是在布西内克斯(Boussinesq)近似下进行的，其中 作用在流体上的浮力定义为αρg(T-T_m)，其中α为热膨胀率，g为 重力，T_m为顶部边界和底部边界的平均温度值。由于当R_a超过临界 值R_a＝1707.762时，最不稳定波数为κ_c＝3.117，因此在研究中使用 的分辨率为101×50。此处使用的Pr为0.71。RB对流建立后，两块 板之间的热传递大大增强。Nu＝1+<u_yT>H/κΔT描述了热传递 的增强。

使用上述总能量守恒(而不是常规方法)的流体模拟将提供流体 模拟以及任何定制的对应计算数据输出的类似描述。但是，当对象 (例如，实际物理对象)被建模时，使用上述曲线方法的这种流体模 拟可以比其它方法更快地和/或通过使用比其它方法少的计算资源来 进行流体模拟。

本文讨论的平衡分布不涉及压力项，因此在碰撞操作期间其值始 终为正，从而提供了相对较强的稳定性。碰撞状态按原样对流后， 导致动量中的等温压力梯度，总能量方程中没有压力对流项。总P 能量的守恒是通过修正停止状态来完成的，例如在平流之后，在计算 提供适当压力梯度的矩之前，将局部压力项加回到粒子状态。

已经描述了许多实现方式。但是，将理解的是，在不脱离权利要 求的精神和范围的情况下，可以进行各种修改。因此，其它实现方式 在所附权利要求的范围内。

Claims (20)

1.一种用于在计算机上模拟流体流动的方法，其中该模拟实现总能量守恒，所述方法包括：

模拟流体在网格上的活动，流体的活动被模拟以便对粒子在网格上的移动进行建模；

在计算机可访问的存储器中为网格中的每个网格位置存储状态向量集，每个状态向量包括与对应网格位置处的可能动量状态中的特别动量状态相对应的多个条目；

通过计算机为网格中的网格位置计算能量值集；

通过计算机对于时间间隔执行粒子到后续网格位置的平流；以及

通过计算机通过将特定总能量值添加到平流的粒子的状态并从在该时间间隔未平流的粒子的状态中减去所述特定总能量值，修正粒子的状态向量集。

2.根据权利要求1所述的方法，其中，在粒子根据修正的状态集平流到后续网格位置之后，所述方法还包括：

在计算矩之前将局部压力项加回到粒子状态以提供适当的压力梯度

3.根据权利要求2所述的方法，其中，局部压力项包括由计算机计算的θ-1个项。

4.根据权利要求2所述的方法，其中，模拟流体流动的活动包括部分地基于第一离散晶格速率集来模拟流体流动，并且所述方法还包括：

部分地基于第二离散晶格速率集来模拟标量的时间演化，第二离散晶格速率集与第一离散晶格速率集相同，或者第二离散晶格速率集的晶格速率的数量与第一离散晶格速率集不同。

5.根据权利要求1所述的方法，其中，模拟总能量的时间演化包括：

收集对于碰撞和能量算子的来自相邻网格位置的传入分布；

对传入分布进行加权；

确定作为碰撞和能量算子乘积的传出分布；以及

传播所确定的传出分布。

6.根据权利要求1所述的方法，其中，为了对于任意普朗特(Prandtl)数恢复精确剪切应力，能量碰撞项由下式给出：

其中，Π＝∑_ic_ic_i(f_i-f_i ^eq是非平衡分量的经过滤的二阶矩。

7.根据权利要求5所述的方法，还包括：

应用零净表面通量边界条件，使得传入分布等于所确定的传出分布。

8.根据权利要求1所述的方法，其中，为网格中的网格位置计算能量值集还包括计算由E_t＝E+v²/2给出的总能量E_t，其中E是能量，v是速度。

9.根据权利要求1所述的方法，其中，模拟活动包括

由计算机对于特定总能量E_t应用平衡分布和第二分布函数，其中第二分布被定义为随着流分布f_i平流的特定标量。

10.根据权利要求9所述的方法，其中，第二分布函数考虑了f_i对能量方程的非平衡贡献，以获得边界附近和跨不同网格分辨率的正确流动行为，其中，用于分布函数的碰撞算子由下式给出：

其中项Ω_f,Ω_q表示相应的碰撞算子；上述方程中使用的平衡分布为：

11.一种有形地体现在计算机可读介质中的计算机程序产品，所述计算机程序产品包括指令，所述指令在执行时模拟物理过程流体流动，并且使得计算系统：

模拟流体在网格上的活动，流体的活动被模拟以便对粒子在网格上的移动进行建模；

在计算机存储器中为网格中的每个网格位置存储状态向量集，每个状态向量包括与对应网格位置处的可能动量状态中的特别动量状态相对应的多个条目；

为网格中的网格位置计算能量值集；

对于时间间隔执行粒子到后续网格位置的平流；以及

通过将特定总能量值添加到平流的粒子的状态并从在该时间间隔未平流的粒子的状态中减去所述特定总能量值，修正粒子的状态向量集。

12.根据权利要求11所述的计算机程序产品，其中，总能量在平流之前添加该项，使得压力被对流，并且为了补偿增加的能量，将能量从停止状态中去除。

13.根据权利要求12所述的计算机程序产品，其中，去除所增加的能量可以使总能量守恒并提供正确的压力速度项。

14.根据权利要求13所述的计算机程序产品，其中，通过修正停止状态

以使得

提供所增加的能量的去除。

15.一种用于模拟物理过程流体流动的计算机系统，所述系统包括：

一个或多个处理器设备；

计算机存储器，耦合到所述一个或多个处理器设备；以及

计算机可读介质，存储有指令，所述指令在执行时模拟物理过程流体流动，并且使得计算系统：

模拟流体在网格上的活动，流体的活动被模拟以便对粒子在网格上的移动进行建模；

在计算机存储器中为网格中的每个网格位置存储状态向量集，每个状态向量包括与对应网格位置处的可能动量状态中的特别动量状态相对应的多个条目；

为网格中的网格位置计算能量值集；

对于时间间隔执行粒子到后续网格位置的平流；以及

通过将特定总能量值添加到平流的粒子的状态并从在该时间间隔未平流的粒子的状态中减去所述特定总能量值，修正粒子的状态向量。

16.根据权利要求15所述的系统，其中，所述系统进一步包括指令以进行如下操作：

在粒子根据修正的状态集平流到后续网格位置之后，在计算矩之前将局部压力项加回到粒子状态以提供适当的压力梯度

17.根据权利要求16所述的系统，其中，局部压力项包括由系统计算的θ-1个项。

18.根据权利要求16所述的系统，其中，模拟流体流动的活动的指令包括以下指令：

部分地基于第一离散晶格速率集来模拟流体流动，并且

部分地基于第二离散晶格速率集来模拟标量的时间演化，第二离散晶格速率集与第一离散晶格速率集相同，或者第二离散晶格速率集的晶格速率的数量与第一离散晶格速率集不同。

19.根据权利要求15所述的系统，其中，为网格中的网格位置计算能量值集的指令包括计算由E_t＝E+v²/2给出的总能量E_t的指令，其中E是能量，v是速度。

20.根据权利要求15所述的系统，还包括使得系统进行如下操作的指令：

对于特定总能量E_t应用平衡分布和第二分布函数，其中第二分布被定义为随着流分布f_i平流的特定标量。

Images