CN111382497A - 任意坐标系中网格上物理流体的计算机模拟 - Google Patents

Abstract

公开了任意坐标系中网格上物理流体的计算机模拟。具体公开了一种计算机实现的用于模拟实体表面附近的流体流动的技术，包括：接收用于表示与实体表面相符的曲线网格的坐标系；利用晶格速度集模拟粒子在流体的体积中的传输，该传输引起粒子之间的碰撞；执行用于粒子的传输的分布函数，该分布函数包括粒子碰撞确定以及与曲线网格相关联的粒子分布的变化；由计算系统在应用于粒子动量值的约束下在坐标系中执行平流操作；以及由计算机系统通过将在坐标系中确定的粒子动量值和空间坐标转换为曲线空间中的动量和空间值而将从模拟得到的值映射到曲线网格上。

Description

任意坐标系中网格上物理流体的计算机模拟

背景技术

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

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

在标准LBM模型中，粒子速度具有一组离散的常数值，并且该常数值在与三维(3D)均匀立方体笛卡尔网格对应的简单布拉维(Bravais)晶格上形成从一个晶格位点到其相邻位点的精确链路。

已经进行了各种尝试以将LBM扩展到基于任意坐标系的网格(任意(一个或多个)网格)。主要的代表性方法之一是放松一对晶格位点之间的一对一平流映射。在这样的任意网格上，从其原始网格位点平流之后的粒子通常不会降落在单个相邻位点上。在那些解决方案中，粒子的位置由那些网格位点通过插值来表示。

发明内容

虽然在许多应用中很重要，但均匀立方体笛卡尔网格对其它应用却构成了基本限制。例如，在现实的流体模拟中，模拟的通常是实体几何弯曲表面。笛卡尔网格不呈现与实体几何弯曲表面的平滑相符性。此外，现实的物理流通常具有集中在某些空间区域和方向中的小结构。例如，在所谓的湍流边界层中，垂直于壁的流动比例远小于切线方向或流体区域的大部分。因此，在边界层内部的壁的法线方向上，对空间分辨率的要求明显更高。立方体笛卡尔网格无法提供在不同方向上具有不同空间分辨率的灵活性。

根据一方面，一种计算机实现的用于模拟实体表面附近的流体流动的方法，包括：由计算系统接收用于表示与实体表面相符的曲线网格的坐标系；利用晶格速度集，模拟粒子在流体的体积中的传输，该传输引起粒子之间的碰撞；执行用于粒子的传输的分布函数，该分布函数包括粒子碰撞确定以及与曲线网格相关联的粒子分布的变化；由计算系统在应用于粒子动量值的约束下在坐标系中执行平流操作；以及由计算机系统通过将在坐标系中确定的粒子动量值和空间坐标转换为曲线空间中的动量和空间值，来将从模拟得到的值映射到曲线网格上。

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

本文公开的方法将在笛卡尔网格上预测的当前基于LBM的模拟扩展到曲线空间中的非笛卡尔网格框架。本文公开的方法基于体积方程，使得满足质量和动量守恒。所得的质量连续性方程将在曲线坐标中具有正确的形式，因此该方法无需引入任何人工的质量源项来校正所得的流体动力学中的伪影。此外，如在流形(manifold)上的连续动力学理论中，在所公开的扩展的LBM中，唯一的外部源项是由于曲线空间而引起的惯性力的存在。该惯性力项不对质量产生贡献，并且是在不依赖于连续动力学理论中的分析形式的情况下被构造的。

离散空间和时间中的该惯性力项在流体动力学极限中渐近地恢复了连续动力学理论中的力。另外，在离散时空LBM模型中，惯性力为底层欧几里得空间强制执行了精确的动量守恒。该力项被构造成使得该力项以适当的离散时间矩向系统增加动量，以便产生粘性阶次的正确结果的Navier-Stokes流体动力学。

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

附图说明

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

图2描绘了示出用于在曲线坐标中制定晶格波尔兹曼模型的操作的流程图。

图3描绘了示出使用与所模拟的物理对象相符的存储的曲线网格进行操作的流程图。

图4描绘了示出使用以曲线坐标表达的晶格玻尔兹曼模型的模拟操作的流程图。

图5描绘了示出应用于分布函数的约束操作的流程图。

图6A-6B描绘了具有在曲线空间和笛卡尔空间中表达的形状的物理对象的表示。

图7A-7B描绘了曲线空间和笛卡尔空间中粒子的移动的表示。

图8A-8B描绘了曲线空间和笛卡尔空间中粒子的移动的另一种表示。

图9和图10图示了在非欧几里得空间中表示的两个LBM模型的速度分量。

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

图12是微观块(microblock)的透视图。

图13A-13B是图1的系统使用的非欧几里得空间中的晶格结构的图示。

图14和15图示了可变分辨率技术。

图16图示了粒子的移动。

图17图示了受表面的面元影响的区域。

图18图示了表面动力学的流程图。

图19图示了粒子从体元到表面的移动。

图20是用于执行表面动力学的过程的流程图。

图21是表示流体模拟的屏幕截图。

具体实施方式

对模拟空间建模

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

f_i(x+c_i，t+1)＝f_i(x，t)+C_i(x，t) 方程(I.1)

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

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

BGK算子是根据(不管碰撞的细节)分布函数都经由碰撞接近由{f^eq(x,v,t)}给出的良好定义的局部平衡的物理论点来构造的：

其中参数τ表示经由碰撞达到平衡的特征弛豫时间。

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

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

f_i(x+c_i，t+1)＝f_i(x，t)+C_i(x，t)

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

ρ(x，t)＝∑_if_i(x，t)；ρ(x，t)u(x，t)＝∑_ic_if_i(x，t) 方程(I.3)

其中，ρ和u分别是流体密度和流体速度。

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。

虽然图1示出了存储器18中的网格准备引擎32，但网格准备引擎可以是在与服务器12不同的系统上执行的第三方应用。无论网格准备引擎32是在存储器18中执行还是在不同于服务器12的系统上执行，网格准备引擎32都接收用户提供的网格定义30，并且网格准备引擎32准备网格并(存储或)将所准备的网格发送到模拟引擎34。模拟引擎34包括粒子碰撞交互模块、粒子边界建模模块和平流操作。系统10访问存储2D和/或3D网格(笛卡尔和/或曲线)、坐标系和库的数据储存库38。

现在参考图2，示出了用于模拟物理对象的表示附近的流体流动的过程40。在这个示例中，物理对象是翼型体(airfoil)。翼型体的使用仅仅是说明性的，因为物理对象可以是任何形状，并且特别地具有(一个或多个)弯曲的表面。该过程例如从客户端系统14或通过从数据储存库38检索而接收42网格，例如与所模拟的物理对象相符的曲线网格。在其它实施例中，基于用户输入的外部系统或服务器12生成与正被模拟的物理对象相符的曲线网格。该过程根据检索到的曲线网格预先计算44几何量，并使用与检索到的曲线网格对应的预先计算的几何量执行46动态晶格波尔兹曼模型模拟。

现在参考图3，下面解释曲线晶格的构造50的方面以及对预先计算的几何量的确定。选择曲线坐标系(与要被模拟的物理对象相符)，使得通过q唯一地定义x＝x(q)，其中q≡(q¹,q²,q³)是沿着参数化基础曲线的三个非共线全等形(non-co-linear congruencies)的坐标值。

系统10以3D曲线网格开始52，其中曲线网格具有在给定空间中的布局，使得其在空间中的所有顶点位置是已知的并且被指定。令x为3维(3D)欧几里得空间中的任何空间点。仅在曲线网格的位点(即，顶点)上定义54空间点x。一旦给定曲线网格(即，“原始”曲线网格{x})，就知道网格{x}上的所有空间位置x_i。

系统10选择56坐标系{q}，该坐标系{q}是空间点x与曲线网格上的位点(q)之间的一对一映射。对于曲线网格{x}上的任何位点x，存在与x相关联的唯一值q，即，x＝x(q)，使得x是唯一定义的。系统10如下在网格上构造q的坐标值：

对于任何位点x(q)，在曲线网格上，该位点沿着第α个(α＝1,2,3)坐标曲线在正方向或负方向上最近的相邻位点都是空间点x_±i。由于这种唯一的映射x_±α＝x(q_±α)，其中q_±α是相邻位点的唯一坐标值，因此完全有可能选择q的空间变化，使得q_±α和q的第α个坐标分量值仅相差恒定距离(d₀)，

并且

其中(α,β＝1,2,3)。在这个示例中，常数d₀被选择为单位元素(1)，而不失一般性。

在这种构造下，坐标值{q}提供了简单的统一的3D立方体笛卡尔晶格结构，其中晶格间距为单位元素(值d₀)。这种变体的“笛卡尔”晶格{q}是由欧几里得空间中原始曲线网格{x}的变形(弯曲、扭曲和拉伸/压缩)产生的。因此，“笛卡尔”晶格{q}的拓扑结构与原始曲线网格{x}相同，但是所得的笛卡尔晶格在非欧几里得空间上。

当提供曲线网格时，指定网格上所有顶点{x}的空间位置，并且也完全确定了曲线网格上从任何一个顶点到另一个顶点的距离。从位点x(q)到该位点的邻居之一x(q_±α)(α＝1；2；3)的距离向量D_±α(q)被定义为：

D_±α(q)≡x(q_±α)-x(q)；α＝1,2,3 (方程1)

例如，D_±1(q)≡x(q¹±1,q²,q³)-x(q¹,q²,q³)。

由于一般曲线网格的空间不均匀性，因此从一个网格位点到其最近邻居位点的空间距离通常随位置而变化。换句话说，D_±α(q)是q的函数。此外，沿着第α个坐标曲线的正方向上的距离值通常不等于负方向上的距离值。

明确地说，就距离向量而言，D_α(q)≠-D_-α(q)。例如，根据由(方程1)给出的定义，

D_-1(q)＝x(q¹-1,q²,q³)-x(q¹,q²,q³)＝-D₁(q¹-1,q²,q³)≠-D₁(q)＝-(x(q¹+1,q²,q³)-x(q¹,q²,q³))

方程2

该不等式仅在当曲线网格是均匀的笛卡尔晶格(使得|D_α|＝常数，与空间坐标值q无关)时才变得处处相等。

系统10可以为所选择的曲线网格存储58所构造的坐标系。

依据基础微分几何学的概念，系统10沿着每个坐标方向构造60正切基础向量

g_α(q)≡[D_α(q)-D_-α(q)]/2；α＝1,2,3 (方程3)

通过这种构造，实现了奇偶对称性，使得g_α(q)＝-g_-α(q)。与欧几里得空间中的笛卡尔坐标系不同，曲线坐标系的基础切向量g_α(q)(α＝1,2,3)通常不是正交的。即，

因此，系统10构造62对应的度量张量(metric tensor)，该度量张量基于上述基础正切向量被定义为：

g_αβ(q)≡g_α(q)·g_β(q),α,β＝1,2,3 (方程4)

以及以x(q)为中心的单元的体积J被定义为：

J(q)≡(g₁(q)xg₂(q))·g₃(q) (方程5)

并且选择适当的“旋向性(handedness)”，使得J(q)＞0，其中J(q)是均匀布拉维晶格的空间中的常数。(如本文所使用的“旋向性”是指如通常在向量分析中提到的所得向量的方向约定(通常为右手规则))。它可以通过以下来验证：

g(q)＝det[g_αβ(q)]＝J²(q) (方程6)

其中det[g_αβ(q)]是度量[g_αβ(q)]张量的行列式。余切基础向量g^α(q)(α＝1,2,3)被构造64为：

g¹(q)≡g₂(q)xg₃(q)/J(q)

g²(q)≡g₃(q)xg₁(q)/J(q)

g³(q)≡g₁(q)xg₂(q)/J(q) (方程7)

基础正切向量和余切向量彼此正交，其中

是Kronecker德尔塔(delta)函数。

类似地，逆度量张量被定义为：

g^αβ(q)≡g^α(q)·g^β(q),α,β＝1,2,3 (方程8)

并且逆度量张量是度量张量的逆，[g^αβ(q)]＝[g_αβ(q)]^-1或

和det[g^αβ(q)]＝1/det[g_αβ(q)]。

在上面定义的基本几何量的情况下，将晶格玻尔兹曼速度向量引入66在一般曲线网格上，类似于标准笛卡尔晶格上的速度向量。

常数

是正整数或负整数或零，并且它是三维坐标数组

的第α个分量值。例如，在所谓的D3Q19中，希腊语索引的范围是从0到18，

c_i∈{(0，0，0)，(±1，0，0)，(0±1，0)，(0，0，±1)，(±1，±1，0)，(±1，0，±1)，(0，±1，±1)}

满足一组矩各向同性和归一化条件，以便恢复正确的完整Navier-Stokes流体动力学。这些是，当存在适当的一组恒定权重{ω_ii＝1,...b}时，该组晶格分量向量允许上至6阶的矩各向同性，即

其中T₀是恒定温度并且取决于一组晶格向量的选择，并且

是Kronecker德尔塔函数。注意的是，三维数组为c_i，而e_i(q)为方程9的晶格向量。

定义了用于与LBM模型一起使用的一组特定几何量，如下：

当网格是均匀笛卡尔网格时，项

消失。

因此，一旦指定曲线网格，就可以完全确定上述所有几何量，并且可以存储68与要模拟的物理对象相符的曲线网格，并且因为几何量0被完全确定，因此可以在开始动态LBM模拟之前对其进行预先计算(参见图2)。

参考图4，模拟过程46(图2)根据例如适于曲线空间的修改后的晶格玻尔兹曼方程(LBE)来模拟80粒子分布的演变。该过程在LBM非欧几里得空间中执行粒子到下一个单元q的平流82、通过将非欧几里得空间中的坐标系的动量和空间坐标转换到欧几里得空间中的曲线系统来将从模拟得到的结果值映射84到原始曲线网格上，并且将欧几里得空间中的映射结果值渲染86到显示器设备上，等等。

参考图5，示出了模拟过程46的更多细节。模拟过程46(图2)根据例如适于曲线空间的修改后的晶格玻尔兹曼方程(LBE)来模拟80(图4)粒子分布的演变，参见方程12。该过程在LBM非欧几里得空间中执行粒子到下一个单元q的平流82(图4)。该过程对粒子平流施加约束83。约束过程83约束83a由于与一般曲线网格的曲率和不均匀性相关联的有效外力而导致的粒子分布的变化、对动量通量施加83b约束、使用定义的83c平衡分布函数，该平衡分布函数恢复流体动力学极限中曲线坐标中正确的欧拉方程和Navier-Stokes方程。因此，该过程通过将非欧几里得空间中的坐标系统中的动量和空间坐标转换成欧几里得空间中的曲线系统来将从模拟得到的结果值映射84(图4)到原始曲线网格上，并将欧几里得空间中的映射结果值渲染86(图4)到显示器设备上，等等。

现在参考图6A，示出了在欧几里得空间(现实世界空间)中的一般3D曲线网格93中的翼型体92，其中仅图示了二个维度(为了清楚起见)。如图所示，网格93大致与翼型体92的形状相符，具体而言是与翼型体92的外表面相符。相符性的程度可以比所示的程度更大或者比图6A中所示的程度更小。

现在参考图6B，示出了非欧几里得空间中的笛卡尔网格95中的(图6A的)翼型体92的表示92’(仅图示了二维)。在这个示例中，翼型体92的表示92’被图示为矩形实心体(仅为了说明目的)。如图所示，笛卡尔网格95通常与翼型体的表示92’的形状相符，具体而言是与其外表面相符。相符性的程度可以比所示的程度更大或更小。

现在参考图7A、图7B，在欧几里得空间(现实世界空间，但为了清楚起见仅图示了二维)中的一般3D曲线网格93中示出了表示流体的粒子的动量的向量96。如图所示，向量96表示沿着穿过曲线空间96的直线而移动的流体的粒子。但是，在变体的笛卡尔非欧几里得空间中，表示该粒子的向量96’被视为具有弯曲的运动。

现在参考图8A、图8B，与图7A，图7B相反，在欧几里得空间(现实世界空间，但为了清楚起见仅图示了二维)中的一般3D曲线网格93中示出了表示流体的粒子的动量的向量97。如图所示，向量97表示以弯曲运动移动通过曲线空间96的流体的粒子。但是，在变体的笛卡尔非欧几里得空间中，表示该粒子的向量97’被视为具有直线运动。现在将描述这些过程的具体细节。

曲线晶格上的体积晶格波尔兹曼建模

所描述的是体积晶格玻尔兹曼建模方法，虽然该方法一般而言适用于各种公式化，但作为说明性示例，将讨论所谓的“等温LBM”的公式化。

修改后的分布类似于粒子分布的演变的标准分布。曲线晶格上的体积晶格波尔兹曼建模被提供为：

N_i(q+c_i，t+1)＝N_i(q，t)+Ω_i(q，t)+δN_i(q，t) 方程12

其中N_i(q,t)是在时间t处对于t的单位元素增量在单元q中属于离散方向c_i的粒子的数量，而不失一般性。(方程12)中的项Ω_α(q,t)是满足局部质量和动量守恒方程的碰撞项：

方程12中的额外项δN_i(q,t)表示由于与一般曲线网格的曲率和不均匀性相关联的有效外力而引起的粒子分布的变化。这个额外项在笛卡尔晶格上的标准LBM中消失。

粒子密度分布函数通过以下方程与N_i(q,t)相关：

J(q)f_i(q，t)＝N_i(q，t) 方程14

如前面所定义的，其中J(q)是以q为中心的单元的体积。基本流体量由标准流体动力矩给出，

其中ρ(q,t)和u(q,t)是位置q处和时间t的流体密度和流体速度。

使用(方程9)中的关系，上面的速度矩被重写为：

并且曲线坐标系中的速度分量值由以下方程给出：

为了简化表示，三维流体速度阵列被定义为U(q,t)≡(U¹(q,t),U²(q,t)U³(q,t))，并且方程17等效地被表达为：

方程18具有与基于标准笛卡尔晶格的LBM相同的流体速度形式。类似地，方程13中碰撞项的动量守恒也可以被写为：

在LBM中，碰撞项通常采用线性化形式，其中M_ij和f_i(q,t)分别表示碰撞矩阵和平衡分布函数。特别地，所谓的Bhatnagar-Gross-Krook(BGK)形式对应于:

其中τ是碰撞弛豫时间。为了恢复正确的Navier-Stokes流体动力学，除了方程13和方程19之外，碰撞矩阵需要满足附加条件。

BGK形式满足这样的附加属性。如上所述，方程12中的额外项δN_i(q,t)表示由于与一般曲线网格的曲率和不均匀性相关联的有效外力引起的粒子分布的变化。这个额外项在笛卡尔网格上的标准LBM中消失。平流过程是从曲线网格中的一个位点到标准LBM中的另一个位点的精确一对一跃点(hop)，即

N_i(q+c_i,t+1)＝N′_i(q,t) (方程22)

其中N′_i(q,t)是(q,t)处的碰撞后分布。由于是曲线网格，所以虽然从单元q平流的粒子的数量与到达单元q+c_i的粒子的数量完全相等(参见方程22)，但动量在平流期间发生变化。

一般而言，

(e_i(q+c_i)N_i(q+c_i,t+1))≠e_i(q)N′_i(q,t)

在上述方程中，不等式的左侧是平流过程结束时的动量值，而右侧是该过程开始时的动量值。存在不等式是因为粒子的路径由于曲线网格而是弯曲(以及拉伸或压缩)的，使得其在平流结束时的速度从其原始值改变。这与在均匀笛卡尔晶格上的情况根本不同，因为粒子在整个平流过程中具有恒定的速度。

因此，总动量值存在以下不等式，

其中(23)中不等号的右侧表示从所有相邻单元中平流出的动量总量，而左侧是沿着弯曲路径平流之后到达单元q的总动量。因此从方程22和方程23中看出，通过平流从所有相邻单元到单元q的净动量变化由以下方程给出：

同样，通过平流出单元q而到单元q的所有相邻单元的净动量变化由以下方程给出：

随后，如果如下对δN_i(q,t)施加约束：

则对于底层欧几里得空间，质量守恒被保留并且精确的动量守恒被恢复。这里

χ(q，t)＝[χ^I(q，t)+χ^o(q，t)]/2

x(q,t)＝[x^I(q,t)+x^o(q,t)]/2

更具体而言，方程26中的第一约束没有通过δN_i(q,t)引入质量源。方程26中的第二约束引入了“惯性力”，其恰好等于在任何晶格位点和时间t在底层欧几里得空间中满足动量守恒的量。该机制类似于弯曲空间中的连续动力学理论描述。以坐标分量形式写为F^α(q,t)＝χ(q,t)·g^α(q)，其中通过直接放置前面小节中的符号，提供：

其中几何函数

在(方程11)中定义。从(27)和(11)中可以看出，如果曲线网格是规则的均匀笛卡尔网格，则F^α(q,t)消失，如常规LBM中那样。(26)中的第二约束也可以用坐标分量形式表达为：

为了恢复完整的粘性Navier-Stokes方程，还需在下面对动量通量施加附加约束，

其中

的具体形式被选择为：

该方程满足方程26、28和29的力矩约束。定义平衡分布函数的形式是为了在流体动力极限中在曲线坐标中恢复正确的欧拉方程以及Navier-Stokes方程。特别地，以下流体动力矩的基本条件是：

其中，在上面，

这些基本条件通过以下平衡分布形式被满足：

上面的平衡分布形式类似于麦克斯韦-波尔兹曼分布的低马赫数展开形式，但以曲线坐标表达。实际上，如果曲线网格是具有g^αβ＝δ^αβ的均匀笛卡尔网格，则上面的平衡分布形式会缩减到标准LBM平衡分布。通过上面定义的所有量和动力学，在(非欧几里得)均匀笛卡尔晶格上{q}上模拟的晶格玻尔兹曼方程12在曲线坐标的情况下服从Navier-Stokes流体动力学。因此，可以通过下面的简单转换将结果值映射到原始曲线网格上，

ρ(x,t)＝ρ(q,t),

u(x,t)＝U^α(q,t)g_α(q) 方程34

非欧几里得速度空间中的笛卡尔坐标

由欧几里得空间中原始曲线网格{x}变形(弯曲、扭曲和拉伸/压缩)产生的变体“笛卡尔”晶格{q}可以用于以与LBM的常规笛卡尔晶格(x)相同的方式模拟物理主体附近的流体流动，这里假设在“笛卡尔”晶格{q}的非欧几里得空间的拓扑结构中平流时，约束被应用于粒子动量，如以上所讨论的，以返回欧几里得空间中的一般曲线网格。

参考图9，第一模型(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)移动的粒子。

也如图10所示，示出了第二模型(3D-1)200——包括39个速度的三维模型，其中每个速度由图9的箭头之一表示。在这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模型在二维和三维中为流的数值模拟提供了特定类的高效且健壮的离散速度动力学模型。这种类型的模型包括离散速度的特定集合和与那些速度关联的权重。这些速度与速度空间中的(非欧几里得空间中的)笛卡尔坐标的网格点相符，这便于离散速度模型(尤其是被称为晶格玻尔兹曼模型的种类)的准确和高效实现。利用这种模型，流可以以高保真度被模拟。

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

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

Re＝uL/v 方程(I.5)

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

状态空间被表示为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,). 方程(I.6)

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

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

每个体元(即，每个晶格位点)由状态向量f(x)表示。该状态向量完全定义体元的状态并且包括39个条目。这39个条目对应于一个能量零状态、6个能量一状态、8个能量三状态、6个能量四状态、12个能量八状态和6个能量九状态。通过使用这个速度集，系统可以对实现的平衡状态向量产生麦克斯韦-玻尔兹曼统计。

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

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

S＝{F_α} 方程(I.7)

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

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

C.识别受面元影响的体元

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

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

参考图16，对于具有非零速度向量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 方程(I.8)

当状态的速度向量指向面元时(|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) 方程(I.9)

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

V_ia＝∑V_a(x)+∑V_ia(β) 方程(I.10)

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

V_ia＝∑V_ia(x) 方程(I.11)

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) 方程(I.12)

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

Γ_iαV→F＝∑_XΓ_iα(x)＝∑_XN_i(x)V_iα(x) 方程(I.13)

只有其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.14)

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

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

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

Γ_iIN(α)＝Γ_iαF→F+Γ_iαF→F＝∑_xN_i(x)V_iα+∑_βΓ_i(β，t-1)V_iα(β)/V_iα

方程(I.16)

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

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

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

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

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

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

对于平行状态(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.21)

然后，这个法线动量利用推/拉技术被消除(步骤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.22)

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

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

并且其中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.26)

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

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

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

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

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

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

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

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

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

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

对于c_jin_α＞0，

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

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

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

5.从体元到体元的移动

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

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

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

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

参考图19，当体元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)图11。这个步骤可以被称为微动力学或体元内操作。类似地，平流规程可以称为体元间操作。以下描述的微动力学操作也可以被用来在面元碰撞粒子，以产生玻尔兹曼分布。

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

E.可变分辨率

参考图20，可变分辨率(如以上所讨论的)采用不同大小的体元，在下文中称为粗糙体元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面元与其它面元不同。

图21示出了流体模拟的屏幕截图。使用上述曲线方法(而不是常规的笛卡尔方法或用于曲线网格的可能的其它方法)的流体模拟将提供流体模拟以及任何定制的对应计算数据输出的类似描述。但是，当具有曲线表面的对象(例如，实际物理对象)被建模时，使用上述曲线方法的这种流体模拟可以更快地和/或通过使用比其它方法少的计算资源来进行流体模拟。

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

Claims (20)

1.一种计算机实现的用于模拟实体的表面附近的流体流动的方法，所述方法包括：

由计算系统接收用于表示曲线网格的坐标系，所述曲线网格与所述实体的表面相符；

利用晶格速度集，模拟粒子在流体的体积中的传输，其中所述传输引起所述粒子之间的碰撞；

执行用于所述粒子的传输的分布函数，其中所述分布函数包括粒子碰撞确定以及与所述曲线网格相关联的粒子分布的变化；

在应用于粒子动量值的约束下，由计算系统在所述坐标系中执行平流操作；以及

由计算机系统通过将在所述坐标系中确定的粒子动量值和空间坐标转换为曲线空间中的动量和空间值，将从模拟得到的值映射到所述曲线网格上。

2.如权利要求1所述的计算机方法，其中，实体表面是弯曲实体表面。

3.如权利要求1所述的计算机方法，其中，所述网格是欧几里得空间中的曲线网格。

4.如权利要求1所述的计算机方法，其中，所述坐标系是在黎曼几何中定义的非欧几里得空间。

5.如权利要求1所述的计算机方法，其中，一般曲线坐标系具有到所选择的坐标系的一对一映射，所选择的所述坐标系是非欧几里得空间中的笛卡尔坐标系。

6.如权利要求3所述的计算机方法，其中，粒子平流被执行，使得在每个时间步长中，粒子分布函数从一个体元(晶格位点)唯一地平流到另一个体元。

7.如权利要求1所述的计算机方法，其中，粒子分布的所述变化使所述计算系统模拟在所述曲线网格中的每个时间增量施加在每个体元中的力的施加，以在所述坐标系中强制实施精确的动量守恒。

8.如权利要求3所述的计算机方法，其中，体积表示被使用，使得一般曲线网格中的每个体元具有定义的体积，并且质量守恒被强制实施，并且质量和动量密度被定义。

9.如权利要求1所述的计算机方法，其中，映射还包括：

确定由于非欧几里得空间中的笛卡尔坐标系中的网格和实体表面的表示而受到约束的动量状态。

10.一种用于模拟实体的表面附近的流体流动的装置，所述装置包括：

存储器；

一个或多个处理器设备，被配置为：

接收用于表示曲线网格的坐标系，所述曲线网格与所述实体的表面相符；

利用晶格速度集，模拟粒子在流体的体积中的传输，其中所述传输引起所述粒子之间的碰撞；

执行用于所述粒子的传输的分布函数，其中所述分布函数包括粒子碰撞确定以及与所述曲线网格相关联的粒子分布的变化；

在应用于粒子动量值的约束下，在所述坐标系中执行平流操作；以及

通过将在所述坐标系中确定的粒子动量值和空间坐标转换为曲线空间中的动量和空间值，将从模拟得到的值映射到所述曲线网格上。

11.如权利要求10所述的装置，其中，实体表面是弯曲实体表面。

12.如权利要求10所述的装置，其中，所述网格是欧几里得空间中的曲线网格。

13.如权利要求10所述的装置，其中，所述坐标系是在黎曼几何中定义的非欧几里得空间。

14.如权利要求10所述的装置，其中，一般曲线坐标系具有到所选择的坐标系的一对一映射，所选择的所述坐标系是非欧几里得空间中的笛卡尔坐标系。

15.如权利要求12所述的装置，其中，粒子平流被执行，使得在每个时间步长中，粒子分布函数从一个体元唯一地平流到另一个体元。

16.如权利要求10所述的装置，其中，粒子分布的所述变化使所述计算系统模拟在所述曲线网格中的每个时间增量施加在每个体元上的力的施加，以在所述坐标系中强制实施精确的动量守恒。

17.如权利要求12所述的装置，其中，体积表示被使用，使得一般曲线网格中的每个体元具有定义的体积，并且质量守恒被强制实施，并且质量和动量密度被定义。

18.如权利要求10所述的装置，其中，映射还包括：

确定由于非欧几里得空间中的笛卡尔坐标系中的网格和实体表面的表示而受到约束的动量状态。

19.一个或多个机器可读的硬件存储设备，其存储能够由一个或多个处理设备执行以使计算机执行以下操作的指令：

接收用于表示曲线网格的坐标系，所述曲线网格与所述实体的表面相符；

利用晶格速度集，模拟粒子在流体的体积中的传输，其中所述传输引起所述粒子之间的碰撞；

执行用于所述粒子的传输的分布函数，其中所述分布函数包括粒子碰撞确定以及与所述曲线网格相关联的粒子分布的变化；

在应用于粒子动量值的约束下，在所述坐标系中执行平流操作；以及

通过将在所述坐标系中确定的粒子动量值和空间坐标转换为曲线空间中的动量和空间值，将从模拟得到的值映射到所述曲线网格上。

20.如权利要求19所述的一个或多个机器可读的硬件存储设备，其中，所述网格是欧几里得空间中的曲线网格。

Images