CN112069742A - 使显式数值方案稳定化 - Google Patents

Abstract

本公开涉及使显式数值方案稳定化。公开了用于模拟固体表面周围的流体流动的计算机实现的技术，包括：接收模拟空间的模型，该模拟空间包括被表示为体元的集合的晶格结构和物理对象的表示，其中体元具有用于考虑物理对象的表面的适当的分辨率，模拟流体体积中的粒子的移动，其中粒子的移动引起粒子之间的碰撞，由计算系统识别两个体元之间的面，其中至少一个面违反稳定性条件，使用两个体元附近的空间平均温度梯度计算修正的热通量，其中至少一个面违反稳定性条件，以及由计算系统对粒子执行平流操作以到往后续的体元。

Description

使显式数值方案稳定化

背景技术

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

计算流体动力学是流体力学的分支，其涉及计算机实现的数值分析技术以分析和模拟物理环境中的流体流动。

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

在包括“晶格玻尔兹曼方法”(LBM)的实际计算流体动力学模拟问题中发现的一个共同方面是涉及复杂几何形状的问题，例如用于表面和体积离散化的不规则栅格(grid)。

发明内容

扩散主导(diffusion-dominated)的物理现象的数值模拟通常用于涉及传导性热转移、质量扩散、电传导等的应用。这些现象的控制等式被公式化为一组偏微分等式(PDE)，其包括不稳定(不稳定的)扩散和体积源(volumetric source)项。数值解涉及对感兴趣的空间域进行离散化，并应用时间积分技术以在时间上推进解。空间离散化过程通常使用高度自动化的栅格生成工具来实现，而时间离散化(用于应用时间积分技术以在时间上推进解，即时间步子(time-step))被选择以在可接受的计算成本下确保数值解的稳定性和准确性。

特别地，通常被称为时间推进(time marching)方案的Courant-Friedrichs-Lewy(CFL)约束的稳定性特性确定了可以在不向解引入显著不稳定性的情况下使用的最大时间步子。通常采用两种类型的时间推进方案——隐式和显式。

隐式方法通过构造来满足CFL约束，并且因此可以使用大的时间步子而不使解不稳定(然而，太大的时间步子通常导致不准确的结果)。隐式方法需要矩阵系数的大系统的解，因此使得它们的实际实现既不平凡又在计算上是昂贵的。

另一方面，显式方法实现起来非常简单、(每次迭代)计算成本低并且高度可并行化。然而，显式方法需要满足CFL约束。显式扩散方案的该CFL约束规定了由(κΔ_t)/(Δ_x^2)给出的CFL数小于特定的限制(其是O(1)的)，其中κ是扩散率，Δ_x是最小空间栅格的尺寸，而Δ_t是时间步子。换句话说，如果空间栅格尺寸Δ_x在域中的任何位置以因子F减小，则时间步子Δ_t将必须以因子F²减小以便维持数值稳定性。

因此，对于具有小尺寸元素的空间栅格，显式方法可能需要极小的时间步子，从而严重影响模拟性能。即使这种小尺寸元素的数量在模拟域中非常有限，也是这种情况，因为整个域中的最小元素确定CFL条件并因此确定时间步子。

对于涉及复杂几何形状的实际问题，使用不规则栅格对于表面和体积离散化是不可避免的。在这些栅格上，Δ_x可能会显著变化，并且由于CFL约束所需的极小时间步子，显式方案的使用可能变得非常低效。

因此，显式方案实施者花费大量时间和努力来试图改进空间栅格的质量，试图减轻该问题。即使这样，也几乎不可能从实际空间域的任何离散化中去除所有小尺寸的元素，并且因此，小的时间步子(至少局部地)是使显式解稳定的唯一方式。

根据一个方面，用于模拟固体表面周围的流体流动的系统、计算机实现的方法和计算机程序产品包括接收模拟空间的模型，该模拟空间包括被表示为体元(voxel)的集合的晶格结构和物理对象的表示，其中体元具有适当的分辨率以考虑物理对象的表面，模拟流体体积中的粒子的运动，其中粒子的运动引起粒子之间的碰撞，识别两个体元之间的面(面元(facet))，其中所述面中的至少一个面违反稳定性条件，使用两个体元附近的空间平均温度梯度(temperature gradient)来计算修正的热通量(heat flux)，其中所述面中的该至少一个面违反稳定性条件，对粒子执行平流(advection)操作以到达后续的体元。

当至少一个元素违反约束(例如，Counton-Friedrichs-Lewy(CFL)约束)时，所公开的技术引入对两个相邻元素之间的热通量计算的修正。对热通量计算的这些修正取决于两个元素的材料和几何性质以及元素的紧邻处的感兴趣的量的现有状态，并且帮助使数值解稳定而不管两个元素的尺寸，并且确保空间和时间精度(spatio-temporal accuracy)。当两个相邻元素较大(并且因此满足CFL约束)时，新提出的通量计算为显式方案实现方式，隐含了这种新方法与显式方法是一致的，并且还克服显式方法中的上述缺陷。

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

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

附图说明

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

图2描绘了示出基于晶格玻尔兹曼模型的针对热通量计算的操作模拟的流程图。

图3-图7示出了涉及当元素中的至少一个元素违反约束时对两个相邻元素之间的热通量计算进行修正的方面的流程图。

图8和9示出了在欧几里德空间中表示的两个LBM模型的速度分量(现有技术)。

图10是物理过程模拟系统所遵循的过程的流程图，该物理过程模拟系统使用两个相邻元素之间的修正的热通量计算。

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

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

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

图15示出了粒子的运动(现有技术)。

图16示出了受表面的面元影响的区域(现有技术)。

图17示出了表面动力学(现有技术)。

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

具体实施方式

对模拟空间建模

在基于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的时间演化的众所周知的晶格玻尔兹曼等式。左手侧表示由于所谓的“流化过程”造成的分布的变化。流化过程是一块流体在网格(mesh)位置出发，然后沿其中一个速度向量移动到下一个网格位置的时候。在那个点，计算“碰撞算子”，即，附近流体块对出发的流体块的影响。流体只能移动到另一个网格位置，因此速度向量的适当选择是必要的，使得所有速度的所有分量都是共同速率(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，它们与网格和模拟引擎34一起提供模拟过程30。

在存储器18中的是网格准备引擎32和模拟引擎34。虽然图1示出了存储器18中的网格准备引擎32，但是该网格准备引擎可以是在与服务器12不同的系统(例如，系统14或另一系统)上执行的第三方应用。无论网格准备引擎32是在存储器18中执行还是在与服务器12不同的系统上执行，网格准备引擎32都接收用户供应的网格定义28，并且网格准备引擎32准备网格并将所准备的网格发送到模拟引擎34。

模拟引擎34包括粒子碰撞相互作用模块34a、粒子边界模型模块34b和执行平流操作的平流模块34c。系统10访问存储2D和/或3D网格和库的数据储存库38。平流模块34c包括子模块36，其根据对两个相邻元素之间的通量计算的修正来执行平流操作，如下面讨论的。

在模拟引擎中执行模拟之前，将模拟空间建模为体元的集合。通常，使用计算机辅助设计(CAD)程序来生成模拟空间。例如，CAD程序可以用于绘制位于风洞中的微型器件。此后，处理CAD程序产生的数据，以增加具有适当分辨率的晶格结构，并考虑模拟空间内的对象和表面。

现在参考图2，示出了用于模拟在物理对象的表示附近的流体流动的过程。在本文将讨论的示例中，物理对象是翼型(airfoil)。然而，翼型的使用仅仅是说明性的，因为物理对象可以是任何形状，并且特别地可以具有(一个或多个)平面和/或弯曲表面。该过程例如从客户端系统14或通过从数据仓库检索来接收35a用于正被模拟的物理对象的网格。在其他实施例中，外部系统或服务器12基于用户输入而生成用于被模拟的物理对象的网格。该过程根据检索到的网格预先计算35b几何量，并且使用与检索到的网格对应的预先计算的几何量来执行35c动态晶格玻尔兹曼模型模拟。晶格玻尔兹曼模型模拟包括根据从引擎34b(图1)产生的边界确定(图2中未示出)对LBM网格中的粒子分布34a的演化35d和粒子到下一单元格

的平流35e的模拟。平流35c过程测试37a CFL约束违反，并且如果发现违反，则使用引擎36(图1)应用对通量计算的修正37b。

针对不规则空间栅格上的扩散问题的稳定化显式数值方案

为了本说明书的目的，假定了显式欧拉(Euler)方案和有限体积公式。在以下描述中，感兴趣的量是温度，而支配等式是热传导等式。该数值方案需要计算元素的所有面上的热通量。随后，这些通量被累加并且被用于更新所考虑的元素的温度。考虑共享面的两个相邻元素α和β。根据热传导的傅立叶定律，热通量为：

其中

是在公共面处的热导率，是

正交于该公共面的温度梯度，并且“m”被用于指定在时间步子“m”处求取所述量。通常使用的傅立叶定律的负号被去掉，因为考虑的是热量进入α(而不是热量离开α)。计算所使用的温度梯度以确保粒子平流的平滑度，尤其是在存在不同尺寸的元素的情况下。如果两个相邻元素α和β满足CFL约束，则通过将热通量乘以公共面的面积A^αβ和时间步长尺寸Δ_t，来获得在时间步子m到m+1期间跨越公共面的能量转移，即

在传统方法中，在时间步子结束时元素α的最终温度是根据到达α的净能量转移(来自所有面的能量转移的总和)而计算的：

注意，上面的等式(3)表明，温度变化与净热通量成比例，与元素的尺寸

成反比，即，对于小的元素，相同的净能量转移导致较大的温度变化。

参考图3，通量修正计算过程36包括时间步子的确定。在通量计算过程36中，时间步子被测试36a。当时间步子没有大到足以违反两个元素中的至少一个元素的CFL约束时36b，上述形式(等式3)将不可能导致解的数值不稳定性，因此可以使用上述形式的通量计算。然而，当时间步子足够大以致违反两个元素中的至少一个元素的CFL约束时36c，上述形式(等式3)可能导致解的数值不稳定性。在这种情况下，应用修正的通量计算过程36d。

在不失一般性的情况下，可以如下解释该修正的方法：

假设元素α小于相邻元素β，并且因此至少对于元素α违反了CFL约束。这种数值不稳定性的出现是因为对于元素α(假设其尺寸较小)，假设被用于计算

的温度梯度在整个时间步子Δ_t期间保持在恒定值是不正确的。如上所述，对于相同的净能量转移，小元素的温度变化较大，因此标准的显式时间积分要求减小时间步子以确保恒定温度梯度假设是有效的。对于给定的Δ_t，只要问题中存在该不稳定性，这个问题就会持续，并且仅在每个元素上的所有传入和传出的通量彼此精确平衡时的稳定状态下才会消失。

参考图4，作为修正的方法的一部分，通量修正的通量计算将如等式(1)中定义的项

细分49a为两部分：(1a)施加的通量(applied flux)

49b，其将(在上面的求和中)被用于α的温度演变(temperature evolution)，以及(1b)余额通量(balanceflux，或者称为余额通量)通量

49c，其将被传输到界面αβ的另一侧而不改变元素的温度，并且它们被表达如下：

其中，项

和ΔG由下式给出：

以及

在上述等式中，这些几何特征由在计算温度梯度时使用的距离d^αβ、元素体积

和

和(由相邻元素共享)的公共面的面积A^αβ表示。材料性质由

和

(这里ρ表示密度并且C_p表示比热)和热导率

来说明。求和中的通量项

和

提供了对分别可能存在于元素α和β不同面上的通量的估计。

两个通量项

和

的物理解释如下。所施加的通量

表示(由等式(1)给出的)总通量

中的可以用于元素α的温度演变而不引入数值不稳定性的一部分。这种形式可以从由元素α和β组成的隔离系统的第一原理导出，从而包括对该系统附近的连续演变的温度场的影响的估计。为此，该项示出了对元素的几何/热属性

以及这些元素与周围环境的相互作用(ΔG)二者的依赖性。

应当注意，在等式(4)中，

仍然依赖于在前一时间步子处的在其他面处观察到的通量，以提供对于在当前时间步子期间在那些面处的正在进行的相互作用的估计。对于强烈的瞬态问题，这导致了

和

间的失配。

称为余额通量的第二项

解释了这种失配，该通量被传输穿过α到达界面αβ另一侧上的其相邻元素。该余额通量仅在它被沉积在大到足以满足CFL约束的元素中时才应用于温度演变，在此之前该通量沿通量方向相继传输。

上述方案严格地确保在任何两个元素α和β之间的界面处结合正确量的总通量

同时精确地控制可用于小元素α的温度演变的通量的量。总之，该方案能够保持数值稳定性以及良好的空间和时间精度。最后，应当注意，在稳定状态下，所施加的通量变为等于全部通量，

并且因此余额通量变为等于零，

修正的通量计算过程的特征

如以上描述中所提到的，并且如图5-7中所示，修正的通量计算过程36需要若干算法过程。修正的通量计算过程包括识别要应用热通量的修正定义的面52a。在不违反CFL条件的所有面处，使用热通量的标准定义(等式3)52b。在违反CFL条件的所有面处，例如，在两个元素(其中至少一个违反CFL条件准则)之间的任何面处，使用修正的热通量计算52c。使用所考虑的元素附近的空间平均温度梯度来计算总热通量，以确保解的平滑性。(相反，标准热通量利用基于传统差分形式计算的温度梯度)。修正的热通量计算过程52包括将修正的热通量划分52d成两个部分——施加的通量项49b和余额通量项49c。然后根据该划分执行通量计算52e。

现在参考图6，使用施加的通量计算52f的通量计算(等式3)总是被用于所考虑的元素的温度演变等式中。

现在参考图7，余额通量可以或可以不被用于温度演变，这取决于所涉及的元素的尺寸。如果元素足够小以致于违反CFL约束52g，则计算余额通量52h并且将其在通量方向上传输52i到相邻元素，并且如果相邻元素足够大以致于不违反CFL约束，则余额通量被用于52j温度演变。余额通量沿着通量的方向相继地传输，直到余额通量最终被转移到足够大的元素(足够大以至于不违反CFL约束)，在该元素处余额通量被应用于温度演变。

算法唯一性和优点

已知几种方法来克服在具有变化的元素尺寸的不规则栅格上的扩散问题中的数值不稳定性的问题。最常用的方法是对栅格生成工具施加附加约束以减少这种场景。即使那样，由于该问题不能完全避免，所以通常的做法是使用足够小以确保稳定性的全局时间步子，或者当遇到小的栅格元素时使用局部子循环。第一种方法(小的全局时间步子)实质上增加了计算成本(即使在空间栅格上任何地方存在小尺寸的元素的单次出现也是如此)，而第二种方法(局部子循环)增加了算法及其实现的复杂性。一种替代方法是使用隐式方案而不是显式方案，或者至少将显式方案的使用限制于靠近小元素的(一个或多个)局部区域。这种隐式方案方法受到实施复杂度以及解的非局部性质的影响，解的非局部性质为系统提供了对于大规模并行计算机实施方案而言不方便的等式。

相反，修正的通量计算方法提供了几个明显的优点。修正的通量计算方法允许使用基于时间精度考虑(而不是栅格中最小元素的尺寸)而选择的单个时间步长尺寸。对于每种可想到的情况，这在计算成本和实现的简易性方面是巨大的益处。修正的通量计算方法具有对两个相邻元素的几何性质的依赖性，因此确保修正的通量计算方法将正常工作而不管元素的尺寸如何。因此，可以在很大程度上放宽对网格生成过程的通常约束(栅格质量、尺寸等)。计算各种项的计算成本是合理的，因为该项的数学形式是简单的并且不涉及任何迭代。这与可能具有较高计算成本的现有方法(减少时间步子尺寸或使用混合隐式-显式方案)形成鲜明对比。此外，由于公式的体积性质，该方案保持了精确的守恒性，这在许多应用中是需要的。修正的通量计算本质上仍然是显式的，并且需要来自距所考虑的元素较小距离内的元素的信息。这种显式性意味着需要对现有计算系统进行从原始系统实现方案的最小改变。保留了原始显式方法的并行化特性，因此修正的通量计算方法可以利用大规模并行计算机实现方案。

因此，与现有方法相比，修正的通量计算方法具有优势。虽然以上对通量计算修正方法的描述是基于用于热传导的有限体积公式，但是该方法实际上适用于许多扩散主导的问题。

对模拟空间建模

在基于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的模拟系统的一般讨论，其可以与CAD过程结合使用以进行流体流动模拟。

参考图8，第一模型(2D-1)200是包括21个速度的二维模型。在这21个速度中，一个速度(205)表示不移动的粒子；三组四个速度表示在沿晶格的x或y轴的正方向或负方向以归一化速率(r)(210-213)、归一化速率的两倍(2r)(220-223)或归一化速率的三倍(3r)(230-233)移动的粒子；以及两组四个速度表示相对于x和y晶格轴二者以归一化速率(r)(240-243)或归一化速率的两倍(2r)(250-253)移动的粒子。

参考图9，示出了第二模型(3D-1)260——包括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模型在二维和三维中为流的数值模拟提供了特定类的高效且健壮的离散速度动力学模型。这种类型的模型包括离散速度的特别集合和与那些速度关联的权重。这些速度与速度空间中的笛卡尔坐标的栅格点相符，这有助于离散速度模型(尤其是被称为晶格玻尔兹曼模型的种类)的准确和高效实现。利用这种模型，流可以以高保真度被模拟。

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

物理过程模拟系统根据使用上述修正的通量计算过程的过程270进行操作。晶格的分辨率可以基于正在被模拟的系统的雷诺数(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_ii＝(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个能量九状态。通过使用这个速度集，系统可以对实现的平衡状态向量产生麦克斯韦-玻尔兹曼统计。

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

微观块在图11中示出。

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

S＝{F_α} 等式(I6)

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

参考图13，不同的分辨率水平可以在模拟空间的不同区域中使用，以提高处理效率。通常，对象322周围的区域320是最关心的并且因此利用最高分辨率进行模拟。因为粘度的影响随着离对象的距离而减小，所以采用降低的分辨率水平(即，扩大的体元体积)来模拟在离对象322按增加的距离隔开的区域324，326。

类似地，如图14中所示，较低的分辨率水平可以被用来模拟对象342的较不显著特征周围的区域340，而最高的分辨率水平被用来模拟对象342的最显著特征(例如，前沿和后缘表面)周围的区域344。边远区域346利用最低的分辨率水平和最大的体元来模拟。

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

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

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

参考图15，对于具有非零速度向量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.执行模拟

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

1.用于表面的边界条件

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

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

2.从体元到面元的聚集

作为模拟粒子和表面之间的相互作用的第一步，粒子从体元被聚集并被提供给面元(步骤278)。如以上所指出的，体元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)＝∑_XN_i(x)V_iα(x) 等式(I12)

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

3.从面元到面元的移动

接下来，粒子在面元之间被移动(步骤280)。如果用于面元F_α的传入状态(c_in_α<0)的平行六面体G_iα被另一面元F_β相交，则由F_α接收的状态i粒子的一部分将来自面元F_β。特别地，面元F_α将接收在前一时间递增期间由面元F_β产生的状态i粒子的部分。

现在参考图17，其中示出了在先前时间递增期间由面元F_β产生的状态i粒子的关系。在图17中，示出了被面元F_β相交的平行六面体G_iα的部分380等于被面元F_α相交的平行六面体G_iβ的部分382。如以上所指出的，相交的部分被表示为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.执行面元表面动力学

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

对于所有i，

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

P_n(α)＝n_α·P(α) 等式(I.20)

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

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

然后，基于传入的通量分布、考虑CFL约束违反的修正的通量计算和玻尔兹曼分布，确定用于面元F_α的传出通量分布(步骤398)。首先，传入的通量分布Γ_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)

这种操作校正质量和能量通量，同时保留切向动量通量未被更改。如果流在面元的邻域是大致均匀的并且接近平衡，则这种调整较小。在调整之后，基于邻域平均属性加上由于邻域的非均匀性或非平衡属性造成的校正，结果产生的法线动量通量稍微被更改为是平衡压力的值。如果违反CFL约束，则该过程将修正的通量计算方法应用285于图10的过程中涉及的任何通量计算。

5.从体元到体元的移动

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

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

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

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

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

其中N(x)是源体元。

6.从面元到体元的分散

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

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

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

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

7.执行流体动力学

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

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

可变分辨率

可变分辨率(如US2013/0151221A1中所讨论的)也可以被采用并且将使用不同尺寸的体元，例如粗体元和精细体元。

通过利用基于晶格玻尔兹曼的独特瞬态物理原理，系统可以执行精确预测真实世界条件的模拟。例如，当变化的影响对于设计和预算非常显著时，工程师在构建任何原型之前，在设计过程的早期评估产品性能。该系统可以使用CAD几何结构精确有效地执行空气动力学、气动声学和热管理模拟。该系统可以执行模拟以解决以下应用：空气动力学(空气动力学效率；车辆处理；污染和水管理；面板变形；行驶动力学)，气动声学(温室风噪声；车身底部风噪声；间隙/密封噪声；镜子、汽笛和音调噪音；天窗和车窗抖振；旁路/社区噪音；冷却风扇噪音)，热管理(冷却气流；热保护；制动冷却；行驶循环模拟；关机和浸泡；电子和电池冷却；ROA/进气口)，气候控制(机舱舒适性；HVAC单元和配电系统性能；HVAC系统和风扇噪音；除霜和除雾)，动力总成：(动力传动系统冷却；排气系统；冷却套；发动机缸体)，污水和水管理(柱溢出、污垢和灰尘堆积、轮胎喷雾)。

本说明书中描述的主题和功能操作的实施例可以在数字电子电路、有形体现的计算机软件或固件、计算机硬件(包括本说明书中公开的结构及其结构等同物)或它们中的一个或多个的组合中实现。本说明书中描述的主题的实施例可以实现为一个或多个计算机程序(即，在有形非暂时性程序载体上编码的一个或多个计算机程序指令模块，用于由数据处理装置执行或控制数据处理装置的操作)。计算机存储介质可以是机器可读存储设备、机器可读存储基板、随机或串行存取存储器设备、或它们中的一个或多个的组合。

术语“数据处理装置”指的是数据处理硬件，并且包括用于处理数据的所有类型的装置、设备和机器，包括例如可编程处理器、计算机或多个处理器或计算机。该装置还可以是或进一步包括专用逻辑电路(例如，FPGA(现场可编程门阵列)或ASIC(专用集成电路))。除了硬件之外，该装置还可以包括为计算机程序创建执行环境的代码(例如，构成处理器固件、协议栈、数据库管理系统、操作系统或它们中的一个或多个的组合的代码)。

计算机程序也可以被称为或描述为程序、软件、软件应用、模块、软件模块、脚本或代码，可以用任何形式的编程语言编写，包括编译或解释语言或者声明性或过程性语言，它可以以任何形式部署，包括作为独立程序或作为模块、组件、子例程或适用于计算环境的其他单元。计算机程序可以但不必对应于文件系统中的文件。程序可以存储在保存其他程序或数据的文件的一部分中(例如，存储在标记语言文档中的一个或多个脚本、专用于所讨论的程序的单个文件中、或者在多个协调文件(例如，存储一个或多个模块，子程序或代码部分的文件)中)。可以部署计算机程序，使得程序在一个计算机上或在位于一个站点上或分布在多个站点上并通过数据通信网络互连的多个计算机上执行。

适合于执行计算机程序的计算机可以基于通用或专用微处理器或两者，或任何其他类型的中央处理单元。适用于存储计算机程序指令和数据的计算机可读介质包括介质和存储器设备上的所有形式的非易失性存储器，包括例如半导体存储器设备(例如，EPROM，EEPROM和闪存设备)、磁盘(例如，内部硬盘或可移动盘)、磁光盘、CD ROM和DVD-ROM盘。处理器和存储器可以由专用逻辑电路补充或并入专用逻辑电路中。

本说明书中描述的主题的实施例可以在包括后端组件(例如，作为数据服务器)或包括中间件组件(例如，应用服务器)或包括前端组件(例如，具有图形用户界面的客户端计算机或用户可以通过其与本说明书中描述的主题的实现交互的web浏览器)或者一个或多个这样的后端、中间件或前端组件的计算系统中实现。系统的组件可以通过任何形式或介质的数字数据通信(例如，通信网络)互连。通信网络的示例包括局域网(LAN)和广域网(WAN)(例如，互联网)。

计算系统可以包括客户端和服务器。客户端和服务器通常彼此远离，并且通常通过通信网络进行交互。客户端和服务器的关系借助于在各个计算机上运行并且彼此具有客户端-服务器关系的计算机程序而产生。在一些实施例中，服务器将数据(例如，HTML页面)发送到用户设备(例如，用于向与用户设备交互的用户显示数据和接收用户输入的目的)，用户设备充当客户端。可以在服务器处从用户设备接收在用户设备处生成的数据(例如，用户交互的结果)。

已经描述了主题的特定实施例。其他实施例在以下权利要求的范围内。例如，权利要求中记载的动作可以以不同的顺序执行并且仍然实现期望的结果。作为一个示例，附图中描绘的过程不一定需要所示的特定顺序或连续顺序来实现期望的结果。在某些情况下，多任务处理和并行处理可能是有利的。

Claims (14)

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

由一个或多个计算系统接收模拟空间的模型，所述模拟空间包括被表示为体元的集合的晶格结构和物理对象的表示，其中所述体元具有考虑所述物理对象的表面的适当的分辨率；

由所述一个或多个计算机系统模拟流体体积中的粒子的运动，其中所述粒子的运动引起所述粒子之间的碰撞；

由所述计算系统识别两个体元之间的面，其中所述面中的至少一个面违反稳定性条件；

由所述计算系统使用所述两个体元附近的空间平均温度梯度来计算修正的热通量，其中所述面中的所述至少一个违反所述稳定性条件；

由所述计算系统对所述粒子执行平流操作以到达后续的体元。

2.如权利要求1所述的计算机方法，还包括：

在没有一个面违反所述稳定性条件的相邻面处，计算基于温度差形式而计算的温度梯度。

3.如权利要求2所述的计算机方法，其中，所述温度差形式确定到达所述元素中的每一个元素的净能量转移以计算在时间步子结束时所述元素的最终温度，所述净能量转移对应于来自所述元素的所有面的能量转移的总和。

4.如权利要求1所述的计算机方法，其中计算所述修正的热通量还包括：

由所述计算系统计算施加的通量；以及

由所述计算系统计算余额通量。

5.如权利要求4所述的计算机方法，其中，计算出的施加的通量被用于所考虑的元素的温度演变。

6.如权利要求4所述的计算机方法，其中，所述余额通量取决于所述元素的尺寸而被用于温度演变。

7.如权利要求6所述的计算机方法，其中，当所述元素的尺寸足够小以至于违反约束时，所述余额通量在所述温度演变中被使用。

8.如权利要求4所述的计算机方法，还包括：

由所述计算机系统将所述余额通量在所述通量的方向上传输到一个或多个相邻元素。

9.如权利要求4所述的计算机方法，其中，所述余额通量沿着所述通量的方向相继地传输，直到所述余额通量最终被转移到足够大的元素，在该元素处所述余额通量被应用于该足够大的元素处的温度演变。

10.如权利要求9所述的计算机方法，其中，当所述元素足够大时，所述余额通量被用于温度演变。

11.如权利要求1所述的计算机方法，其中，所述条件包括稳定性特征。

12.如权利要求1所述的计算机方法，其中，所述条件是对时间推进方案的Counton-Friedrichs-Lewy(CFL)约束，所述Counton-Friedrichs-Lewy(CFL)约束确定能够用于维持稳定分布的最大时间步长尺寸。

13.一种用于模拟固体表面周围的流体流动的装置，所述装置包括：

存储器；

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

执行根据权利要求1-10中的任何一项或多项所述的功能。

14.一种有形存储在一个或多个机器可读硬件存储设备上的计算机程序产品，包括可执行指令，所述可执行指令能够由一个或多个处理设备执行以使计算机执行根据权利要求1-10中的任何一项或多项所述的功能。

Images