CN110175342A - 用于高速流的基于格点玻尔兹曼的求解器 - Google Patents

用于高速流的基于格点玻尔兹曼的求解器 Download PDF

Info

Publication number
CN110175342A
CN110175342A CN201910126500.2A CN201910126500A CN110175342A CN 110175342 A CN110175342 A CN 110175342A CN 201910126500 A CN201910126500 A CN 201910126500A CN 110175342 A CN110175342 A CN 110175342A
Authority
CN
China
Prior art keywords
entropy
lattice point
state
particle
fluid
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN201910126500.2A
Other languages
English (en)
Inventor
P·乔帕拉克里斯纳恩
张绕阳
陈沪东
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Dassault Systemes Americas Corp
Original Assignee
Dassault Systemes Simulia Corp
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Dassault Systemes Simulia Corp filed Critical Dassault Systemes Simulia Corp
Publication of CN110175342A publication Critical patent/CN110175342A/zh
Pending legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/25—Design optimisation, verification or simulation using particle-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/30—Circuit design
    • G06F30/36—Circuit design at the analogue level
    • G06F30/367—Design verification, e.g. using simulation, simulation program with integrated circuit emphasis [SPICE], direct methods or relaxation methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/08—Probabilistic or stochastic CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Algebra (AREA)
  • Software Systems (AREA)
  • Operations Research (AREA)
  • Databases & Information Systems (AREA)
  • Microelectronics & Electronic Packaging (AREA)
  • Computing Systems (AREA)
  • Fluid Mechanics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Computer Graphics (AREA)

Abstract

描述了用于在计算机上模拟流体流的技术，其涉及稳定的熵求解器。所述技术包括：跨网格模拟流体的活动，流体的活动被模拟以跨网格对粒子的移动进行建模；在计算机能够访问的存储器中存储网格中每个网格位置的状态矢量的集合，状态矢量中的每一个状态矢量包括多个条目，该多个条目与对应网格位置处的可能的动量状态中的特定动量状态对应；通过以下操作来模拟流的熵的时间演变：针对碰撞操作，从相邻的网格位置收集传入的分布集合；由计算机计算每个位置中的标量值；将输出分布确定为碰撞操作和热源添加的乘积；以及通过计算机针对一时间间隔执行粒子到后续网格位置的平流来修改流。

Description

用于高速流的基于格点玻尔兹曼的求解器
背景技术
Lattice Boltzmann方法(LBM，格点玻尔兹曼方法)是用于流体模拟的一类计算流体动力学(CFD)方法。不是求解Navier-Stokes方程，而是用诸如Bhatnagar-Gross-Krook(BGK)之类的碰撞模型求解离散Boltzmann方程以模拟牛顿流体(Newtonian fluid)的流。通过模拟有限数量的粒子上的流动和碰撞过程，内在的粒子交互表明了适用于更大质量的粘性流行为的微观世界。
发明内容
根据一个方面，一种用于在计算机上模拟流体流的方法包括跨网格模拟流体的活动，流体的活动被模拟以跨网格对粒子的移动建模；在计算机可访问的存储器中存储网格中每个网格位置的状态矢量集合，这些状态矢量中的每个状态矢量包括多个条目，这些条目与对应网格位置处的可能的动量状态中的特定动量状态对应；通过针对碰撞操作从相邻网格位置收集传入的分布集合、由计算机计算每个位置的标量值、将输出分布确定为碰撞操作和热源添加的乘积，来模拟流的熵的时间演变；以及通过计算机针对一时间间隔执行粒子到后续网格位置的平流来修改该流。
以上方面中可以包括下列特征中的一个或多个。
模拟流体流的活动包括部分地基于离散格点速率的第一集合来模拟流体流，并且该方法还包括部分地基于离散格点速率的第二集合来模拟熵量的时间演变。离散格点速率的第二集合与离散格点速率的第一集合具有相同的格点速率。该方法还包括由计算机计算来自传入的格点集(lattice set)矢量的高阶误差项、确定高阶误差项的平均值、以及从碰撞算子中减去高阶误差项的平均值。计算附加热源项并将其加到第二状态，并且该方法还包括由计算机计算通过流体粘度加热和通过流体传导加热的效应以及由计算机计算熵扩散并该熵扩散从附加热源项移除。该方法还包括由计算机计算网格中网格位置的物理量的集合。该方法提供了Lattice Boltzmann熵求解器，其避免了二阶速度项。碰撞算子涉及非平衡(non-equilibrium)计算而没有速度的任何二阶项。该方法还包括通过提供格点矢量qi的附加集合来求解熵，以表示具体的熵s，并且这些状态的时间演变由下式给出：
其中qi是格点矢量，x是方向，ci是状态的速度，Δt是时间t的改变，是平衡时的格点矢量，τq是弛豫时间，V是是非平衡贡献，p是压力，并且Qs是附加热源。该方法中碰撞算子与下式相关：
其它方面包括计算机程序产品和装置，诸如计算机和数据处理系统。
上述各方面中的一个或多个可以提供以下优点中的一个或多个。熵求解器避免项a(ρvv)，在其他方式下该项引入误差(例如，特别是对于高速流的不稳定性)，因为该项与Ma2(速率平方)成比例。因此，该熵求解器对于高速应用是稳定的。在Lattice Boltzmann方法的背景下，基于LB的熵求解器对于解决需要高度准确的瞬态结果连带可压缩性影响的高速应用是有效的。具有高温度比的应用也受益于LB熵求解器，这是由于其独特的碰撞模型提供的增强的稳定性。
本申请依据35U.S.C.§119要求于2018年2月20日提交的标题为“LatticeBoltzmann based solver for high speed flows”的美国临时专利申请序列No.62/632,584的优先权，其全部内容通过引用结合于此。
附图说明
图1描绘了用于模拟流体流的系统，该系统包括用于可压缩流的熵求解器。
图2描绘了示出制定具有熵求解器的Lattice Boltzmann模型模拟的操作的流程图。
图3描绘了使用具有熵求解器的Lattice Boltzmann模型的模拟操作的流程图。
图4A、4B针对流体动力学和基于Lattice Boltzmann的技术分别描绘了所模拟的流。
图5是声压级(SPL)与频率的关系图。
图6A-6B是模拟的方向角与频率和声压级的关系图。
图7和8图示了在欧几里德空间中表示的两个LBM模型的速度分量(现有技术)。
图9是由物理过程模拟系统遵循的过程的流程图。
图10是微块的透视图(现有技术)。
图11A-11B是图1的系统使用的格点结构的图示(现有技术)。
图12和13图示了可变分辨率技术(现有技术)。
图14图示了粒子的移动(现有技术)。
图15图示了受表面的面元影响的区域(现有技术)。
图16图示了表面动力学的流程图(现有技术)。
图17是用于执行表面动力学的过程的流程图(现有技术)。
图18图示了粒子从体素(voxel)到表面的移动。
图19是代表流体模拟的屏幕截图(现有技术)。
在下面的附图和描述中阐述了本发明的一个或多个实施例的细节。根据描述和附图以及权利要求，本发明的其它特征、目的和优点是清楚的。
具体实施方式
用于模拟流体流的一种方法是所谓的Lattice Boltzmann模型(LBM)。在基于LBM的物理过程模拟系统中，流体流由分布函数值表示，使用众所周知的、描述分布函数的时间演变的Lattice Boltzmann方程(见下面的方程1)在一组离散速度下进行求值。该分布函数涉及两个过程：流动过程和碰撞过程。
其中ξ和u分别是无量纲格点矢量和速度，并且其中fi(x,t)是在时间t处在位置x处沿着格点方向ci的粒子分布函数，其中是平衡分布，并且τ是特征弛豫时间。
流动过程是当流体袋从网格位置开始然后沿着速度矢量之一移动到下一个网格位置。在那个时候，计算“碰撞因子”，即，附近的流体袋对起始流体袋的影响。流体只能移动到另一个网格位置，因此速度矢量的恰当选择是必要的，以便所有速度的所有分量都是公共速率的倍数。
第一个方程的右边是前面提到的“碰撞算子”，它表示由于流体袋之间的碰撞引起的分布函数的改变。碰撞算子的特定形式是Bhatnagar、Gross和Krook(BGK)算子。碰撞算子迫使分布函数到达由第二个方程给出的规定值，它是“平衡”形式。
根据构造BGK算子，使得分布函数接近明确定义的局部平衡，而不管碰撞的细节如何：
其中参数τ表示经由碰撞到达平衡的特征弛豫时间。处理粒子(例如，原子或分子)时，弛豫时间通常取为恒定值。
根据这个模拟，常规流体变量(诸如质量ρ和流体速度u)作为方程2的简单求和被获得。
其中ρ和u分别是上面提到的流体密度、速度。
由于对称性的考虑，以这样一种方式选择速度值的集合，使得它们在配置空间中跨越时形成某些格点结构。这种离散系统的动力学服从具有以下形式的LatticeBoltzmann方程(LBE)
fi(x+ciΔt，t+Δt)-fi(x，t)＝Ci(x，t)
其中碰撞算子通常采用BGK形式，如上所述。
通过恰当选择平衡分布形式，理论上可以证明Lattice Boltzmann方程产生了正确的流体动力学。即，从fi(x,t)导出的流体动力矩服从宏观极限中的Navier-Stokes方程。这些矩由上面的方程(3)定义。
ci和wi的集体值定义LBM模型。该LBM模型可以在可缩放的计算机平台上高效地实现，并且对于时间不稳定的流和复杂边界条件以强稳健性运行。
从Boltzmann方程获得流体系统的运动的宏观方程的标准技术是Chapman-Enskog方法，在Chapman-Enskog方法中采用完全Boltzmann方程的连续近似。在流体系统中，密度的小扰动以声速行进。在气体系统中，声速一般由温度确定。通过特征速度和声速之比来测量流中的可压缩性的影响的重要性，该比被称为Mach数。对于常规的基于LBM的物理过程模拟系统的进一步解释，参考上面通过引用并入的申请US-2016-0188768-A1。
热Lattice Boltzmann模型
由方程1和方程2描述的LBM仅适用于求解具有均匀温度的低速流。对于热流，通常引入单独的求解器来求解能量方程，并且将计算出的温度耦合回到Lattice Boltzmann方法中。温度方程可以通过有限差分方法或者通过引入附加的格点矢量集合来求解。
采用独特的基于LBM的温度求解器。在这种基于LBM的温度求解器中，使用附加的格点矢量集合来表示具体标量，即，伴随流矢量并从一个网格点移动到另一个网格点的温度。它具有若干优点，并且在公开的专利申请US20130151221A1中可获得更详细的描述，该申请的全部内容通过引用并入本文。
本文使用类似的运输机制用于熵求解器，如下面所讨论的。
一旦知道温度，就将温度与流体LBM求解器耦合。传统方法在LBM(方程1)的时间演变中对热流采用以下平衡分布(方程4)。
其中ξ、u和θ分别是无量纲格点矢量、无量纲速度和无量纲压力。然而，对于高温度比的应用，项θ-1将变为负，这使得求解器变得不稳定。
在Pradeep Gopalakrishnan等人的标题为“Temperature Coupling Algorithmfor Hybrid Thermal Lattice Boltzmann Method”的同样未决的美国专利公开US-2016-0188768-A1中详细描述了针对相对宽范围的流条件碰撞期间使平衡分布为正的耦合技术的具体细节，该专利的全部内容通过引用并入本文。
在方程1中阐述的BGK碰撞算子具有相对低的稳定范围。虽然方程1的算子对于低Mach(Mach<0.6)和低温度比应用是有用的，但是方程1的该算子在较高Mach(Mach>0.6)和较高温度比应用中趋于变得不稳定。适用于航空航天领域的高亚音速和跨音速应用的温度求解器需要基于伽利略不变(Galilean-invariant)算子的更稳定的碰撞算子，其具有相对高的稳定性范围。这种碰撞算子的具体细节在美国专利No.9,576,087中详细描述，该专利的全部内容通过引用并入本文。
高速可压缩流的能量守恒
对于高速流，特别是Mach数高于0.3的流，由于可压缩性和加热改变流体的粘度造成的对流体的能量的影响是不可忽略的。由方程5给出的可压缩流的能量方程解释了这种可压缩性影响。
其中e是内能，V是速度，k是导热系数，并且Φ是耗散函数，表示机械能由于粘度而变成热能的比率。对于不可压缩流，压力功项、方程5左侧的第三项和粘性耗散项退出，并且方程5变成相对容易求解的典型的对流-扩散标量方程。
然而，由于存在发散速度项，使用有限差分方法或其它方法整体上求解方程5导致不稳定性。为了避免这个项，通过更可行的选项(即，通过求解下面给出的熵方程(方程6))来强制执行能量守恒：
其中s是系统的熵。
可压缩流的基于Lattice Boltzmann的熵求解器
熵方程(方程6)可以使用像有限体积或有限差分方法之类的传统方法来求解。对于复杂问题，网格具有多个分辨率级别并且具有复杂的边界条件。传统方法在计算这些情况的梯度时存在困难，其引入了强的网格依赖性和数值伪像。
本文描述的是基于Lattice Boltzmann(LB)的熵求解器34a(图1)，其可以处理这些复杂问题。求解熵使用附加的格点矢量qi的第二集合，它们被引入以表示具体的熵“s”。这些状态的时间演变由下式给出：
其中p是压力，并且Qs是附加的源(例如，附加的热源)。方程7与方程1一起将熵项从一个网格点携带到另一个网格点，并且获得作为两个方程的乘积的熵守恒方程(方程8)。
其中S是总熵。这个方法将总熵作为平流的流上的标量项来运输，类似于如美国专利公开US-2016-0188768-A1中所描述的为基于LBM的温度求解器所采用的机制。这种用于运输熵的方法有几个优点，包括在不同网格级别附近和复杂边界附近的相对平滑的过渡。
方程7右侧的第二项是所谓的正则化碰撞算子。在BGK碰撞算子中(由方程1给出)，计算实际状态与平衡状态之间的差异，称为“非平衡部分”。这个非平衡部分被放松(relax)，使得实际状态将放松到平衡状态。
然而，该算子涉及所有阶数的非平衡效应，这对于守恒是不必要的。此外，涉及所有阶数的非平衡效应可能导致不稳定。另一方面，方程7中使用的算子仅考虑非平衡贡献的一阶效应(通过乘以格点矢量计算出的一阶矩)，如下：
然后将这个经滤波的非平衡部分乘以伽利略不变算子ci-V。
现在参考图1，描述包括用于高速和可压缩流的基于LB的熵求解器34b的系统10。这个实现方案中的系统10基于客户端-服务器或基于云的体系架构，并且包括实现为大规模并行计算系统12(独立或基于云)的服务器系统12和客户端系统14。服务器系统12包括存储器18、总线系统11、接口20(例如，用户接口/网络接口/显示器或监视器接口等)和处理设备24。在存储器18中的是网格准备引擎32和模拟引擎34。
虽然图1示出了网格准备引擎32在存储器18中，但网格准备引擎可以是在与服务器12不同的系统上执行的第三方应用。不管网格准备引擎32是在存储器18中执行还是在与服务器12不同的系统上执行，网格准备引擎32都接收用户提供的网格定义30，并且网格准备引擎32准备网格并将准备好的网格发送(和/或存储)到模拟引擎34。模拟引擎34包括碰撞交互模块34a、用于高速流的基于LB的求解器模块34b、边界模块34c以及平流粒子碰撞交互模块34d。系统10访问存储2D和/或3D网格(笛卡尔和/或曲线)、坐标系统和库的数据储存库38。
现在参考图2，示出了用于模拟围绕物理对象的表示的流体流的过程40。在本文将讨论的示例中，物理对象是翼型(airfoil)。然而，翼型的使用仅仅是说明性的，因为物理对象可以是任何形状，并且特别地可以具有(一个或多个)平面和/或弯曲表面。过程40例如从客户端系统14接收42或从数据储存库38检索用于被模拟的物理对象的网格(或栅格)。在其它实施例中，外部系统或者服务器12基于用户输入生成用于被模拟的物理对象的网格。该过程从检索到的网格预先计算44几何量，并使用与检索到的网格对应的预先计算的几何量执行动态的Lattice Boltzmann模型模拟46。Lattice Boltzmann模型模拟包括对粒子分布演变的模拟46a、熵求解器的求值46b，以及粒子和熵到LBM网格中的接下来的单元格的平流46c。
参考图3，模拟过程46根据(例如，适用于用于高速流的基于LB的求解器模块34b(图1)的)经修改的Lattice Boltzmann方程(LBE)来模拟粒子分布的演变。过程46(参见图2)执行碰撞操作46a、熵的时间演变的模拟46b、随后是粒子到LBM空间中的接下来的单元的平流46c。
模拟熵的时间演变46b包括从碰撞操作收集52来自相邻网格位置的传入的分布集合，并且收集54表示具体标量(例如，熵)的格点矢量的附加集合，其在流矢量上并且从一个栅格点移动到另一个栅格点。模拟熵的时间演变还包括由计算机计算56网格中每个位置的一组标量值。计算机随后将传出分布确定58为碰撞操作方程7和附加热源Qs的乘积。
通过计算机计算60来自传入的格点集合矢量的高阶误差项和从法向碰撞算子中减去62的对应平均值来稳定碰撞算子。附加热源项被计算64并且被加到66第二状态，并且是计算机计算附加热源通过流体粘度加热流体和通过流体传导加热流体的效应的结果。通过计算机计算熵扩散，并且从附加热源项中移除68计算出的熵扩散的结果。
由于扩散引起的热源的处理
在没有任何附加热源的情况下，基于LBM的标量求解器将产生典型的标量对流-扩散方程，如下：
其中值ks取决于τq中所使用的弛豫时间。附加热源项(考虑下面给出的方程6与方程10之间的差异)将恢复熵方程的正确行为。
用于高速流的增强的熵碰撞算子
如上面提到的美国专利9,576,087中所提到的，描述了用于流状态的特殊碰撞算子。它具有宽范围的稳定性，并且可以被应用于高Mach应用。
然而，在高Mach数流下，用于熵的碰撞算子变得不稳定。
Chapman-Enskog展开理论被应用于流和熵状态(方程1和方程7)的演变，以恢复Navier-Stokes方程和熵标量方程。扩散项由于碰撞操作而产生，并且其中的误差可能使熵求解器不稳定。
考虑非平衡分量的一阶计算(方程9)，在这种情况下，通过使用泰勒级数展开并使用LBM的基本格点要求，一阶非平衡分量可以被近似为：
在上面的表达式中，乘以压力(pδ)的熵梯度是唯一需要的分量，而附加项(ρvv)引入误差。项(ρvv)是导致特别是高速流不稳定的主要原因，因为该项与Ma2(速率平方)成比例。
对上述非平衡计算的修改如下进行，以移除这个项。
上述非平衡分量移除了高阶误差项，并提供了对高速应用而言稳定的熵求解器。将方程(11)和方程(13)代入方程(7)，熵演变的完整形式由下式给出：
基于LB的熵求解器能力
上面讨论的LB熵求解器对于求解需要高度准确的瞬态结果连带可压缩性影响的高速应用是有效的。具有高温度比的应用也受益于该LB熵求解器，这是由于其独特的碰撞模型提供的增强的稳定性。
本节介绍了一项基准研究，以示出基于LB的熵求解器对可压缩流(尤其是与高速、可变网格分辨率、复杂几何形状和高温度比流相关)的潜在优势。将模拟结果与实验结果和典型的基于的有限差分FD熵求解器进行比较。
现在参考图4A和4B，阐述了在喷射喷嘴的温度比大气温度高2.702倍时对于处于0.901Mach的喷嘴流实施的模拟的结果。为了预测由喷嘴产生的准确噪声，有必要获得更准确的瞬态结果。
图4A针对FD方法的示出了压力振荡50a，并且图4B针对基于LB的方法示出了压力振荡50b。
FD方法在喷嘴52a出口附近产生高水平的噪声54a，其涉及由于可变分辨率而影响喷嘴52a的流下游56a的复杂边界条件。基于LB的方法对于喷嘴52b出口具有更干净(即，更低)的噪声54b水平，并且在喷嘴下游56b也保持精细的结构。
现在参考图5，描绘了以分贝为单位的声压级(SPL)与以Hz为单位的频率的关系图。该图是针对FD和(具有熵求解器的)LB方法的模拟声压级与实验结果的比较。由于数值伪像，FD方法预测噪声水平远高于实验结果(尤其是在1KHz以上的频率)，而基于LB的方法与实验结果相比性能相对好。
现在参考图6A和6B，示出了与实验相比基于LB的求解器在预测声压级图中的准确度。
参考图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个速度的三维模型，其中每个速度由图6中的箭头之一表示。在这39个速度中，一个速度表示不移动的粒子；三组六个速度表示沿着格点的x、y或z的正方向或者负方向上以归一化速率(r)、归一化速率的两倍(2r)或归一化速率的三倍(3r)移动的粒子；八个速度表示相对于所有三个x、y、z格点轴以归一化速率(r)移动的粒子；以及十二个表示相对于x、y、z格点轴中的两个以归一化速率的两倍(2r)移动的粒子。
也可以使用更复杂的模型，诸如包括101个速度的3D-2模型和包括37个速度的2D-2模型。
对于101个速度的三维模型3D-2，一个速度表示不移动的粒子(组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模型提供了一类具体的高效且稳健的离散速度动力学模型，以用于二维和三维流的数值模拟。这种模型包括离散速度和与这些速度相关联的权重的特定集合。速度与速度空间中笛卡尔坐标(欧几里德空间中)的网格点一致，这便于准确和高效地实现离散速度模型，特别是称为Lattice Boltzmann模型的类型。使用此类模型，可以高保真地对流进行模拟。
参考图9，物理过程模拟系统根据过程300操作，以模拟诸如流体流的物理过程。在模拟之前，模拟空间被建模为体素的集合(302)。通常，使用计算机辅助设计(CAD)程序生成模拟空间。例如，CAD程序可以用于绘制位于风洞中的微设备。此后，处理由CAD程序产生的数据，以添加具有适当分辨率的格点结构并考虑模拟空间内的对象和表面。
可以基于被模拟的系统的Reynolds数来选择格点的分辨率。Reynolds数与流的粘度(ν)、流中的对象的特征长度(L)以及流的特征速度(u)相关：
Re＝uL/v. 方程(I.5)
对象的特征长度表示对象的大的尺度特征。例如，如果正在模拟微设备周围的流，那么微设备的高度可以被认为是特征长度。当关注对象的小区域(例如，汽车的侧镜)周围的流时，可以增加模拟的分辨率，或者可以在感兴趣的区域周围采用增加分辨率的区域。随着格点的分辨率增加，体素的尺寸减小。
状态空间被表示为fi(x，t)，其中fi表示在时间t在由三维矢量x表示的格点位点处状态i下每单位体积的元素或粒子的数量(即，状态i下的粒子密度)。对于已知的时间增量，粒子的该数量简称为fi(x)。格点位点的所有状态的组合表示为f(x)。
状态的数量由每个能级内可能的速度矢量的数量确定。速度矢量由具有三个维度x、y和z的空间中的整数线性速率组成。对于多品种模拟，状态的数量增加。
每个状态i表示在特定能级(即，能级零、一或二)处的不同速度矢量。每个状态的速度ci以其在三个维度中的每个维度中的“速率”表示如下：
ci＝(cix，ciy，ciz) 方程(I.6)
能级零状态表示不在任何维度上移动的停止的粒子，即，c停止＝(0,0,0)。能级一状态表示在三个维度中的一个维度中具有a±1速率并且在另外两个维度中具有零速率的粒子。能级二状态表示在所有三个维度上具有a±1速率、或者在三个维度中的一个维度中具有a±2速率并且在其它两个维度中具有零速率的粒子。
生成三个能级的所有可能的排列给出总共39种可能的状态(一个能量零状态，6个能量一状态，8个能量三状态，6个能量四状态，12个能量八状态，以及6个能量九状态)。
每个体素(即，每个格点位点)由状态矢量f(x)表示。状态矢量完全定义体素的状态并且包括39个条目。这39个条目与一个能量零状态、6个能量一状态、8个能量三状态、6个能量四状态、12个能量八状态和6个能量九状态对应。通过使用这个速度集，系统可以产生用于所实现的平衡状态矢量的Maxwell-Boltzmann统计。
为了处理效率，体素被分组为称为微块的2x2x2体积。微块被组织成允许体素的并行处理并且使得与数据结构相关联的开销最小化。微块中体素的简写符号定义为Ni(n)，其中n表示格点位点在微块内的相对位置，并且n∈{0,1,2，...，7}。
微块在图10中示出。
还参考图11A和11B，表面S(图11A)在模拟空间(图11B)中表示为面元(facet)Fα的集合：
S＝{Fα} 方程(1.7)
其中α是枚举特定面元的索引。面元不限于体素边界，而是通常尺寸为与面元相邻的体素的尺寸的数量级或略小于该体素的尺寸，使得面元影响相对少量的体素。为了实现表面动力学，将特性指派给面元。具体而言，每个面元Fα具有单位法向(nα)、表面积(Aα)、中心位置(xα)以及描述面元的表面动态特性的面元分布函数(fi(α))。
参考图12，可以在模拟空间的不同区域中使用不同级别的分辨率以提高处理效率。通常，对象655周围的区域650是最感兴趣的，因此以最高分辨率进行模拟。因为粘度的效应随着距对象的距离而减小，所以采用逐步降低的分辨率级别(即，扩展的体素体积)来模拟与对象655间隔距离增加的区域660、665。
类似地，如图13中所示，较低级别的分辨率可以用于模拟对象775的较不重要的特征周围的区域770，而最高分辨率级别用于模拟对象775的最重要的特征周围的区域780(例如，前沿和后沿表面)。使用最低分辨率级别和最大体素来模拟外围区域785。
识别受面元影响的体素
再次参考图9，一旦模拟空间已被建模(302)，就识别受一个或多个面元影响的体素(304)。体素可以以多种方式受面元的影响。首先，与一个或多个面元相交的体素受到影响，因为该体素相对于非相交的体素具有减小的体积。发生这种情况是因为面元以及由面元表示的表面下面的材料占据了体素的一部分。分数因子Pf(x)指示体素的不受面元影响的部分(即，可以被为其模拟流的流体或其它材料占据的部分)。对于非相交的体素，Pf(x)等于1。
通过将粒子转移到面元或从面元接收粒子而与一个或多个面元交互的体素也被识别为受面元影响的体素。与面元相交的所有体素将包括从面元接收粒子的至少一个状态以及将粒子转移到面元的至少一个状态。在大多数情况下，附加的体素也将包括这些状态。
参考图14，对于具有非零速度矢量ci的每个状态i，面元Fα从由平行六面体Giα定义的区域接收粒子或将粒子转移到该区域，其中平行六面体Giα具有高度，该高度由速度矢量ci和面元的单位法向nα的矢量点积量(|cini|)的幅度以及由面元的表面积Aα定义的基数定义，使得平行六面体Giα的体积Viα等于：
Via＝|cina|Aa 方程(I.8)
当状态的速度矢量指向面元(|ci ni|<0)时，面元Fα从体积Viα接收粒子，并且当状态的速度矢量指向远离面元(|ci ni|>0)时，面元Fα将粒子转移到该区域。如下面将讨论的，当另一个面元占据平行六面体Giα的一部分时(可能发生在诸如内角之类的非凸特征附近的条件)，必须修改这个表达式。
面元Fα的平行六面体Giα可以与多个体素中的部分体素或全部体素重叠。体素或其部分的数量取决于面元相对于体素的尺寸的尺寸、状态的能量、以及面元相对于格点结构的朝向。受影响的体素的数量随着面元的尺寸而增加。因而，如上所述，面元的尺寸通常被选择为位于面元附近的体素的尺寸或小于该体素的尺寸。
由平行六面体Giα重叠的体素N(x)的部分被定义为Viα(x)。使用这个项，在体素N(x)和面元Fα之间移动的状态i粒子的通量Γiα(x)等于该体素中的状态i的密度(Ni(x))乘以与该体素重叠的区域的体积(Viα(x))：
Γiα(x)＝Ni(x)Viα(x) 方程(I.9)
当平行六面体Giα与一个或多个面元相交时，满足以下条件：
Viα＝∑Vα(x)+∑Viα(β) 方程(I.10)
其中第一个总和考虑与Giα重叠的所有体素，第二项考虑与Giα相交的所有面元。当平行六面体Giα未与另一个面元相交时，这个表达式将精简为：
Viα＝∑Viα(x) 方程(I.11)
执行模拟
再次参考图9，一旦识别出受一个或多个面元影响的体素(304)，就初始化定时器以开始模拟(306)。在模拟的每个时间增量期间，通过平流阶段(308-316)来模拟粒子从体素到体素的移动，该平流阶段考虑粒子与表面面元的交互。接下来，碰撞阶段(318)部分地基于上面讨论的熵求解器的使用来模拟每个体素内的粒子的交互。此后，计时器递增(320)。如果递增的计时器未指示模拟完成(322)，那么重复平流和碰撞阶段(308-320)。如果递增的定时器指示模拟完成(322)，那么存储和/或显示模拟的结果(324)。
如图9中所示，“平行路径”用于熵状态的平流。这个并行路径一般包括并行动作308'-318'，然后是用熵求解器34b对流求值。熵状态沿着三维直线格点在体素之间移动314'，类似于上面讨论的体素到体素移动的314。只要递增的定时器不指示模拟完成，该过程就重复(322)。
1.用于表面的边界条件
为了正确模拟与表面的交互，每个面元必须满足四个边界条件。首先，由面元接收的粒子的组合质量必须等于由面元转移的粒子的组合质量(即，到面元的净质量通量必须等于零)。其次，由面元接收的粒子的组合能量必须等于由面元转移的粒子的组合能量(即，到面元的净能量通量必须等于零)。通过要求每个能级(即，能级一和二)的净质量通量等于零，可以满足这两个条件。
另外两个边界条件与粒子与面元交互的净动量相关。对于没有表面摩擦的表面(在本文称为滑动表面)，净切向动量通量必须等于零，并且净法向动量通量必须等于面元处的局部压力。因此，垂直于面元的法向nα的组合的接收和转移动量的分量(即切向分量)必须相等，而平行于面元的法向nα的组合的接收和转移动量的分量(即，法向分量)之间的差异必须等于面元处的局部压力。对于非滑动表面，表面的摩擦按照与摩擦量相关的因子、相对于由面元接收的粒子的组合切向动量减小了由面元转移的粒子的组合切向动量。
2.从体素到面元的搜集
作为模拟粒子与表面之间的交互的第一步，从体素搜集粒子并将粒子提供给面元(308)。如上所述，体素N(x)和面元Fα之间的状态i粒子的通量为：
Γiα(x)＝Ni(x)Viα(x) 方程(I.12)
由此，对于指向面元Fα的每个状态(cinα<0)，由体素提供给面元Fα的粒子的数量是：
ΓiαV→F＝∑X Γiα(x)＝∑X Ni(x)Viα(x) 方程(I.13)
只有Viα(x)具有非零值的体素必须求和。如上所述，面元的尺寸被选择成使得Viα(x)仅对少量体素具有非零值。因为Viα(x)和Pf(x)可以具有非整数值，所以Γα(x)作为实数被存储和处理。
3.从面元移动到面元
接下来，粒子在面元之间移动(310)。如果用于面元Fα的传入状态(cinα<0)的平行六面体Giα与另一个面元Fβ相交，那么由面元Fα接收的状态i粒子的一部分将来自面元Fβ。具体而言，面元Fα将接收在前一时间增量期间由面元Fβ产生的状态i粒子的一部分。
这个关系在图15中示出，其中与面元Fβ相交的平行六面体Giα的一部分1000等于与面元Fα相交的平行六面体Giβ的一部分1005。如上所述，相交的部分表示为Viα(β)。使用这个项，面元Fβ和面元Fα之间的状态i粒子的通量可以被描述为：
Γiα(β,t-1)＝Γi(β)Viα(β)/Viα 方程(I.14)
其中Γi(β,t-1)是在前一时间增量期间由面元Fβ产生的状态i粒子的测度。由此，对于指向面元Fα的每个状态(cinα<0)，由其它面元提供给面元Fα的粒子的数量为：
ΓiαF→F＝∑βΓiα(β)＝∑β Γi(β，t-1)Viα(β)/Viα 方程(1.15)
并且状态i粒子进入面元的总通量为：
用于面元的状态矢量N(α)，也称为面元分布函数，具有与体素状态矢量的M个条目对应的M个条目。M是离散格点速率的数量。面元分布函数N(α)的输入状态被设置为，对于cinα＜0，等于粒子进入那些状态的通量除以体积Viα：
Ni(α)＝ΓiIN(α)/Viα 方程(I.17)
面元分布函数是用于从面元生成输出通量的模拟工具，并且不一定代表实际粒子。为了生成准确的输出通量，将值指派给分布函数的其它状态。使用上述的用于填充向内状态的技术来填充向外状态：
Ni(α)＝ΓiOTHER(α)/Viα 方程(I.18)
对于cinα≥0，其中使用上述用于生成ΓiIN(α)的技术确定ΓiOTHER(α)，但是将该技术应用于除传入状态(ci nα＜0)之外的状态(ci nα≥0)。在另一种方法中，可以使用前一次的ΓiOUT(α)值生成ΓiOTHER(α)，使得：
ΓiOTHER(α，t)＝ΓiOUT(α，t-1). 方程(I.19)
对于平行状态(cinα＝0)，Viα和Viα(x)均为零。在用于Ni(α)的表达式中，Viα(x)出现在分子中(根据用于ΓiOTHER(α)的表达式)，而Viα出现在分母中(根据用于Ni(α)的表达式)。因而，用于平行状态的Ni(α)表示被确定为当Viα和Viα(x)逼近零时Ni(α)的极限。具有零速度的状态(即，静止状态和状态(0，0，0，2)和(0，0，0，-2))的值基于温度和压力的初始条件在模拟开始时被初始化。然后随时间调整这些值。
4.执行面元表面动力学
接下来，对每个面元执行表面动力学，以满足上面图9(312)讨论的四个边界条件。
在图16中图示了用于执行面元的表面动力学的过程。首先，对于所有i，通过确定面元处粒子的组合动量P(α)来确定与面元Fα垂直的组合动量(1105)：
由此，法向动量Pn(α)被确定为：
Pn(α)＝nα·P(α). 方程(I.21)
然后使用推送/拉取技术消除这个法向动量(1110)以产生Nn-(α)。根据这个技术，粒子以仅影响法向动量的方式在状态之间移动。推送/拉取技术在美国专利No.5,594,671中描述，该专利通过引用并入。
此后，Nn-(α)的粒子碰撞以产生Boltzmann分布Nn-β(α)(1115)。如下面关于执行流体动力学所描述的，可以通过将碰撞规则集应用于Nn-(α)来实现Boltzmann分布。
基于传入的通量分布和Boltzmann分布确定(1120)面元Fα的传出通量分布。首先，传入的通量分布Γi(α)与Boltzmann分布之间的差异被确定为：
ΔΓi(α)＝ΓiIN(α)-Nn-βi(α)Viα 方程(I.22)
利用这个差异，对于nαci>0，传出的通量分布为：
ΓiOUT(α)＝Nn-βi(α)Viα-.Δ.Γi*(α), 方程(I.23)
其中i＊是具有与状态i相反的方向的状态。例如，如果状态i是(1，1，0，0)，那么状态i＊是(-1，-1，0，0)。为了考虑皮肤摩擦和其它因素，对于nαci＞0，传出的通量分布可以进一步细化为：
ΓiOUT(α)＝Nn-Bi(α)Viα-ΔΓi*(α)+Cf(nα·ci)-[Nn-βi*(α)-Nn-βi(α)]Viα+(nα·ci)(t1α·ci)ΔNj，1Viα+(nα·ci)(t2α·ci)ΔNj，2Viα
方程(I.24)其中Cf是皮肤摩擦的函数，tiα是垂直于nα的第一切向矢量，t2α是与nα和t1α都垂直的第二切向矢量，并且ΔNj，1和ΔNj，2是与状态i的能量(j)和指示的切向矢量对应的分布函数。这些分布函数根据下式确定：
其中j对于能级1状态等于1，并且对于能级2状态等于2。
ΓiOUT(α)的方程的每个项的功能如下。第一项和第二项强制执行法向动量通量边界条件，使得碰撞在产生Boltzmann分布时是有效的，但包括切向动量通量异常。第四项和第五项纠正这种异常，其可能是由于由于碰撞不充分导致的离散效应或非Boltzmann结构而产生的。最后，第三项添加指定量的皮肤分数以强制执行表面上的切向动量通量的期望改变。下面描述摩擦系数Cf的生成。要注意的是，涉及矢量操纵的所有项都是可以在开始模拟之前计算的几何因子。
由此，切向速度被确定为：
ui(α)＝(P(α)-Pn(α)nα)/ρ， 方程(I.26)
其中ρ是面元分布的密度：
ρ＝∑iNi(α) 方程(I.27)
如上所述，传入的通量分布与Boltzmann分布之间的差异被确定为：
ΔΓi(α)＝ΓilN(α)-Nn-βi(α)Viα 方程(I.28)
于是，传出的通量分布变为：
ΓiOUT(α)＝Nn-βi(α)Viα-ΔΓi*(α)+Cf(nαci)[Nn-βi*(α)-Nn-βi(α)]Viα 方程(I.29)
这与由前一种技术确定的传出通量分布的前两行对应，但不需要校正异常的切向通量。
使用任一种方法，结果得到的通量分布都满足所有动量通量条件，即：
其中pα是面元Fα处的平衡压力，并且基于向该面元提供粒子的体素的平均密度和温度值，并且uα是面元处的平均速度。
为了确保满足质量和能量边界条件，测量每个能级j的输入能量与输出能量之间的差异：
其中索引j表示状态i的能量。然后对于cjinα>0，使用这个能量差来生成差异项：
这个差异项用于修改传出的通量，使得对于cjinα>0，该通量变为：
ΓαjiOUTf＝ΓαjiOUT+δΓαji 方程(I.33)
这个操作校正质量和能量通量，同时保持切向动量通量不变。如果流在面元附近大致均匀并接近平衡，那么该调整是小的。调整之后得到的法向动量通量被略微更改为基于邻域平均特性的平衡压力加上由于邻域的非均匀性或非平衡特性引起的校正的值。
5.从体素到体素移动
再次参考图9，粒子沿着三维直线格点在体素之间移动(314)。这种体素到体素移动是在不与面元交互的体素(即，不位于表面附近的体素)上执行的唯一移动操作。在典型的模拟中，不位于表面足够近而不与表面交互的体素构成体素的大部分。
每个独立的状态表示粒子在三个维度(x、y和z)中的每一个中沿着格点以整数速率移动。整数速率包括：0、±1和±2。速率的符号指示粒子正沿着对应轴移动的方向。
对于不与表面交互的体素，移动操作在计算上非常简单。在每个时间增量期间，状态的整个群体从其当前体素移动到其目的地体素。同时，目的地体素的粒子从该体素移动到它们自己的目的地体素。例如，在+1x和+1y方向上(1,0,0)中移动的能级1粒子从其当前体素移动到在x方向上为+1并且在其它方向上为0的体素。该粒子最终在其目的地体素处于与移动前相同的状态(1,0,0)。体素内的交互将有可能基于与其它粒子和表面的局部交互来改变那个状态的粒子数。如果不是这样，那么粒子将继续沿着格点以相同的速率和方向移动。
对于与一个或多个表面交互的体素，移动操作变得稍微复杂一些。这可以导致一个或多个分数粒子被转移到面元。这种分数粒子到面元的转移导致分数粒子保留在体素中。这些分数粒子被转移到由面元占据的体素中。
参考图17，当用于体素905的状态i粒子的一部分900移动到面元910时(308)，剩余部分915移动到(面元910所在并且状态i的粒子从其指向面元910的)体素920中。因此，如果状态群体等于25并且Viα(x)等于0.25(即，体素的四分之一与平行六面体Giα相交)，那么6.25个粒子将移动到面元Fα并且18.75个粒子将移动到由面元Fα占据的体素。因为多个面元可以与单个体素相交，所以转移到由一个或多个面元占据的体素N(f)的状态i粒子的数量是：
其中N(x)是源体素。
6.从面元到体素的散射
接下来，来自每个面元的传出粒子被散射到体素(316)。从本质上讲，这与粒子从体素移动到面元的搜集相反。从面元Fα移动到体素N(x)的状态i粒子的数量是：
其中Pf(x)考虑部分体素的体积减小。由此，对于每个状态i，从面元指向体素N(x)的粒子的总数量是：
在将粒子从面元散射到体素之后，将它们与从周围体素平流进来的粒子组合，并使结果整数化，某些体素中的某些方向可能下溢(变为负)或者上溢(在八位实现方案中超过255)。在截断质量、动量和能量以适应允许的值范围之后，这将导致这些量的增益或者损失。为了防止这种情况发生，超出边界的质量、动量和能量在违规状态的截断之前被累加。对于该状态所属的能量，等于(由于下溢)增益或(由于溢出)丢失的值的质量的量被添加回具有相同能量并且本身不经受上溢或下溢的随机(或顺序)选定状态。从这种质量和能量的添加所产生的附加动量被累加并添加到截断所得的动量上。通过仅向相同的能量状态添加质量，当质量计数器达到零时，质量和能量二者都被校正。最后，使用推送/拉取技术校正动量，直到动量累加器返回到零。
7.执行流体动力学
执行流体动力学(318)图9。这个步骤可以被称为微动力学或体素内操作。类似地，平流过程可以被称为体素间操作。下面描述的微动力学操作也可以用于碰撞面元处的粒子以产生Boltzmann分布。
通过称为BGK碰撞模型的特定碰撞算子在Lattice Boltzmann方程模型中确保流体动力学。这个碰撞模型模仿真实流体系统中的分布的动力学。碰撞过程可以通过方程I.1和方程I.2的右侧很好地描述。在平流步骤之后，使用方程I.3从分布函数获得流体系统的守恒量，具体而言是密度、动量和能量。从这些量中，由方程(I.2)中的feq表示的平衡分布函数由方程(I.4)完全指定。速度矢量集ci、权重的选择(都在表1中列出)与方程I.2一起确保宏观行为服从正确的流体动力学方程。
参考图18，系统生成1200粒子运输的分布，其中熵的确定与碰撞操作分开。在这个示例中，分布函数包括平流，并且通过利用熵确定增强平流(而不是利用熵确定增强粒子碰撞)来将熵确定包括在分布部分中。一般而言，平流包括粒子的运输(例如，在水平方向上从一个区域到另一个区域)。
在操作中，系统在格点速度集中模拟(1202)一定体积的流体中的粒子的运输，其中运输造成粒子之间的碰撞。该系统还基于用于粒子的运输的分布函数生成(1204)粒子分布，其中分布函数包括热力学步骤和粒子碰撞步骤，并且其中热力学步骤基本上独立于粒子碰撞步骤并且与粒子碰撞步骤分离。
在图12的示例中，系统还确定(1204)确定在特定时间t处流体体积中的特定位置x处的碰撞的碰撞后分布函数fi′(x,t)，其中fi′(x,t)＝fi(x,t)+Ci(x,t)，其中Ci是碰撞算子，并且fi是碰撞前粒子的分布函数。
该系统在1206根据方程7通过使用格点矢量qi的附加集合来求解熵和运输，以表示特定的熵。方程7和方程1将熵从一个网格点带到另一个网格点，从而有效地提供熵作为与流体流的平流一起平流的标量，如美国专利公开US-2016-0188768-A1中对于诸如温度之类的标量的平流所讨论的那样。这个技术在不同网格级别附近和复杂边界附近提供平滑过渡。
可变分辨率
也可以采用可变分辨率(如US 2013/0151221 A1中所讨论的)并且将使用不同尺寸的体素，例如粗略体素和精细体素。
图19示出了流体模拟的屏幕截图。使用了上述熵求解器方法的流体模拟将提供流体模拟的类似描述以及任何常规的对应计算的数据输出。然而，在具有曲线表面的对象(例如，实际物理对象)被建模时，使用了上述熵求解器方法的这种流体模拟可以更快地进行流体模拟，和/或使用比其它方法更少的计算资源，并且还可以用于对高Mach数流体流中的对象进行建模同时为需要高度准确的瞬态结果连带可压缩性影响的应用维持相对稳定的熵求解器。
硬件和软件实现
系统10一般可以被实现为各种不同的电气或电子计算或处理设备中的任何一个，并且可以执行上面讨论的各种步骤的任意组合以控制所公开的热管理系统的各种部件。
系统10一般可以并且可选地包括处理器(或多个处理器)、存储器、存储设备和输入/输出设备中的任何一个或多个。可以使用系统总线互连这些部件中的部分或全部，如图1中所示。处理器能够处理指令以供执行。在一些实施例中，处理器是单线程处理器。在某些实施例中，处理器是多线程处理器。通常，处理器能够处理存储在存储器中或存储设备上的指令，以在输入/输出设备上显示用于用户界面的图形信息，并执行上面讨论的各种监视和控制功能。适用于本文公开的系统的处理器包括通用和专用微处理器，以及任何类型的计算机或计算设备的唯一处理器或多个处理器之一。
存储器存储系统内的信息，并且可以是计算机可读介质，诸如易失性或非易失性存储器。存储设备能够为系统10提供大容量存储。一般而言，存储设备可以包括被配置为存储计算机可读指令的任何非瞬态有形介质。例如，存储设备可以包括计算机可读介质和相关联的部件，包括：磁盘，诸如内部硬盘和可移除盘；磁光盘；以及光盘。适于有形地体现计算机程序指令和数据的存储设备包括所有形式的非易失性存储器，包括例如半导体存储器设备，诸如EPROM、EEPROM和闪存设备；磁盘，诸如内部硬盘和可移除盘；磁光盘；以及CD-ROM和DVD-ROM盘。本文公开的系统的处理器和存储器单元可以由ASIC(专用集成电路)补充或并入其中。
输入/输出设备为系统10提供输入/输出操作，并且可以包括键盘和/或定点设备。在一些实施例中，输入/输出设备包括用于显示图形用户界面和系统相关信息的显示单元。
本文描述的特征，包括用于执行各种测量、监视、控制和通信功能的部件，可以在数字电子电路系统中实现，或者在计算机硬件、固件中或它们的组合中实现。方法步骤可以在有形地体现在信息载体中(例如，在机器可读存储设备中)的计算机程序产品中实现，用于由(例如，系统10的)可编程处理器执行，并且特征可以由执行这种指令的程序以执行上述任何步骤和功能的可编程处理器执行。适于由一个或多个系统处理器执行的计算机程序包括可以直接或间接地被使用的指令集，以使执行指令的处理器或其它计算设备执行某些活动，包括上面讨论的各种步骤。
适于与本文公开的系统和方法一起使用的计算机程序可以用任何形式的编程语言编写，包括编译或解释的语言，并且可以以任何形式部署，包括作为独立程序或作为模块、部件、子例程，或适于在计算环境中使用的其它单元。
已经描述了许多实现方案。然而，应该理解的是，在不脱离权利要求的精神和范围的情况下，可以进行各种修改。因此，其它实现方案在以下权利要求的范围内。

Claims (21)

1.一种用于在计算机上模拟流体流的方法，所述方法包括：
跨网格模拟流体的活动，流体的所述活动被模拟以跨所述网格对粒子的移动进行建模；
在计算机能够访问的存储器中存储所述网格中每个网格位置的状态矢量的集合，所述状态矢量中的每一个状态矢量包括多个条目，所述多个条目与对应网格位置处的可能的动量状态中的特定动量状态对应；
通过以下操作来模拟所述流的熵的时间演变：
针对碰撞操作，从相邻的网格位置收集传入的分布集合；
由计算机计算每个位置中的标量值；
将输出分布确定为所述碰撞操作和热源添加的乘积；以及
通过计算机针对一时间间隔执行所述粒子到后续网格位置的平流来修改所述流。
2.如权利要求1所述的方法，其中模拟流体流的活动包括部分地基于离散格点速率的第一集合来模拟所述流体流，并且所述方法还包括：
部分地基于离散格点速率的第二集合来模拟熵量的时间演变。
3.如权利要求1所述的方法，其中离散格点速率的所述第二集合与离散格点速率的所述第一集合是相同的格点速率。
4.如权利要求1所述的方法，还包括：
由计算机根据传入的格点集合矢量计算高阶误差项；
确定所述高阶误差项的平均值；以及
从碰撞算子中减去高阶误差项的所述平均值。
5.如权利要求1所述的方法，其中附加热源项被计算并被加到第二状态，并且所述方法还包括：
由计算机计算通过流体粘度加热和通过流体传导加热的效应，
由计算机计算熵扩散并将所述熵扩散从所述附加热源项移除。
6.如权利要求1所述的方法，还包括：
由计算机计算所述网格中的网格位置的物理量的集合。
7.如权利要求1所述的方法，其中Lattice Boltzmann熵求解器避免二阶速度项。
8.如权利要求1所述的方法，其中碰撞算子涉及非平衡计算，而没有速度的任何二阶项。
9.如权利要求1所述的方法，还包括：
通过提供格点矢量qi的附加集合来求解熵，以表示具体的熵s，并且这些状态的时间演变由下式给出
其中qi是格点矢量，x是方向，ci是状态的速度，Δt是时间t的改变，是平衡时的格点矢量，τq是弛豫时间，V是是非平衡贡献，p是压力，并且Qs是附加热源。
10.如权利要求9所述的方法，其中所述碰撞算子与下式相关：
11.一种被配置为模拟流体流的计算系统，所述计算系统包括：
一个或多个处理器设备；
存储器，可操作地耦合到一个或多个处理器设备，所述存储器存储计算指令以使所述一个或多个处理器设备：
跨网格模拟流体的活动，流体的所述活动被模拟以跨所述网格对粒子的移动进行建模；
在所述存储器中存储所述网格中每个网格位置的状态矢量的集合，所述状态矢量中的每一个状态矢量包括多个条目，所述多个条目与对应网格位置处的可能的动量状态中的特定动量状态对应；
通过用于以下操作的指令模拟所述流的熵的时间演变：
针对碰撞操作，从相邻的网格位置收集传入的分布集合；
由计算机计算每个位置中的标量值；
将输出分布确定为所述碰撞操作和热源添加的乘积；以及
通过计算机针对一时间间隔执行所述粒子到后续网格位置的平流来修改所述流。
12.如权利要求11所述的系统，其中用于模拟流体流的活动的指令包括用于部分地基于离散格点速率的第一集合来模拟流体流的指令；
部分地基于离散格点速率的第二集合来模拟熵量的时间演变。
13.如权利要求11所述的系统，其中离散格点速率的所述第二集合与离散格点速率的所述第一集合是相同的格点速率。
14.如权利要求11所述的系统，还包括用于以下操作的指令：
根据传入的格点集合矢量计算高阶误差项；
确定所述高阶误差项的平均值；以及
从碰撞算子中减去高阶误差项的所述平均值。
15.如权利要求11所述的系统，其中附加热源项被计算并被加到第二状态，并且所述方法还包括：
计算通过流体粘度加热和通过流体传导加热的效果，
由计算机计算熵扩散并将所述熵扩散从所述附加热源项移除。
16.如权利要求11所述的系统，还包括用于以下操作的指令：
计算所述网格中的网格位置的物理量的集合。
17.如权利要求11所述的系统，其中Lattice Boltzmann熵求解器避免二阶速度项。
18.如权利要求11所述的系统，其中所述碰撞算子涉及非平衡计算而没有速度的任何二阶项。
19.如权利要求11所述的系统，还包括用于以下操作的指令：
通过提供格点矢量qi的附加集合来求解熵，以表示具体的熵s，并且这些状态的时间演变由下式给出
其中qi是格点矢量，x是方向，ci是状态的速度，Δt是时间t的改变，是平衡时的格点矢量，τq是弛豫时间，V是是非平衡贡献，p是压力，并且Qs是附加热源。
20.如权利要求19所述的系统，其中所述碰撞算子与下式相关：
21.一种有形地存储在非瞬态硬件存储设备上的计算机程序产品，所述计算机程序产品包括用于将计算系统配置为执行以下操作的可执行指令：
跨网格模拟流体的活动，流体的所述活动被模拟以跨所述网格对粒子的移动进行建模；
在存储器中存储所述网格中每个网格位置的状态矢量的集合，所述状态矢量中的每一个状态矢量包括多个条目，所述多个条目与对应网格位置处的可能的动量状态中的特定动量状态对应；
通过用于以下操作的指令模拟所述流的熵的时间演变：
针对碰撞操作，从相邻的网格位置收集传入的分布集合；
由计算机计算每个位置中的标量值；
将输出分布确定为所述碰撞操作和热源添加的乘积；以及
通过计算机针对一时间间隔执行所述粒子到后续网格位置的平流来修改所述流。
CN201910126500.2A 2018-02-20 2019-02-20 用于高速流的基于格点玻尔兹曼的求解器 Pending CN110175342A (zh)

Applications Claiming Priority (4)

Application Number Priority Date Filing Date Title
US201862632584P 2018-02-20 2018-02-20
US62/632,584 2018-02-20
US16/274,403 2019-02-13
US16/274,403 US20190258764A1 (en) 2018-02-20 2019-02-13 Lattice Boltzmann Based Solver for High Speed Flows

Publications (1)

Publication Number Publication Date
CN110175342A true CN110175342A (zh) 2019-08-27

Family

ID=65529372

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910126500.2A Pending CN110175342A (zh) 2018-02-20 2019-02-20 用于高速流的基于格点玻尔兹曼的求解器

Country Status (4)

Country Link
US (1) US20190258764A1 (zh)
EP (1) EP3528146A1 (zh)
JP (1) JP7296216B2 (zh)
CN (1) CN110175342A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112687002A (zh) * 2021-03-15 2021-04-20 四川省公路规划勘察设计研究院有限公司 一种三维地质模型网格优化方法

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR102123254B1 (ko) * 2018-08-28 2020-06-16 한국과학기술연구원 짤리스 엔트로피 및 레이니 엔트로피를 이용한 유체유동 시뮬레이션 방법
US20210133374A1 (en) * 2019-10-30 2021-05-06 Dassault Systemes Simulia Corp. Computer System for Simulating Physical Process Using Lattice Boltzmann based Scalar Transport Enforcing Galilean Invariance for Scalar Transport
CN111489447B (zh) * 2020-04-14 2022-04-29 西北工业大学 一种适用于格子Boltzmann方法的直角网格自适应建模方法
US11847391B2 (en) * 2020-06-29 2023-12-19 Dassault Systemes Simulia Corp. Computer system for simulating physical processes using surface algorithm
CN112100099B (zh) * 2020-09-28 2021-06-08 湖南长城银河科技有限公司 一种面向多核向量处理器的格子玻尔兹曼优化方法
CN112765841B (zh) * 2020-12-31 2024-05-07 西北工业大学深圳研究院 一种用以解决流体变物性计算的预处理格子玻尔兹曼方法
CN113792432B (zh) * 2021-09-15 2024-06-18 沈阳飞机设计研究所扬州协同创新研究院有限公司 基于改进型fvm-lbfs方法的流场计算方法
CN116306279B (zh) * 2023-03-15 2024-06-07 重庆交通大学 一种水动力自由面lb模拟方法、系统及存储介质

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1482434A2 (en) * 2001-08-10 2004-12-01 Xencor Protein design automation for protein libraries
US20130116997A1 (en) * 2011-11-09 2013-05-09 Chenghai Sun Computer simulation of physical processes
CN103425833A (zh) * 2013-08-07 2013-12-04 湖南大学 一种基于熵格子波尔兹曼模型的并行流体计算实现方法
US20150356217A1 (en) * 2013-07-24 2015-12-10 Exa Corportion Lattice Boltzmann Collision Operators Enforcing Isotropy and Galilean Invariance
CN106650133A (zh) * 2016-12-28 2017-05-10 中南大学 加热炉内导热流体仿真方法
WO2018026327A1 (en) * 2016-08-04 2018-02-08 Nanyang Technological University An apparatus for transferring heat from a heat source to a heat sink

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5377129A (en) 1990-07-12 1994-12-27 Massachusetts Institute Of Technology Particle interaction processing system
US10360324B2 (en) 2011-12-09 2019-07-23 Dassault Systemes Simulia Corp. Computer simulation of physical processes
US10515159B2 (en) * 2013-03-06 2019-12-24 Dassault Systemes Simulia Corp. Flow-induced noise source identification
JP6434512B2 (ja) 2013-07-31 2018-12-05 エクサ コーポレイション ハイブリッド熱格子ボルツマン法のための温度結合アルゴリズム
CN106529063A (zh) 2016-11-14 2017-03-22 宜兴八达流体技术有限公司 基于cfd技术的流体系统及其设计方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1482434A2 (en) * 2001-08-10 2004-12-01 Xencor Protein design automation for protein libraries
US20130116997A1 (en) * 2011-11-09 2013-05-09 Chenghai Sun Computer simulation of physical processes
US20150356217A1 (en) * 2013-07-24 2015-12-10 Exa Corportion Lattice Boltzmann Collision Operators Enforcing Isotropy and Galilean Invariance
CN103425833A (zh) * 2013-08-07 2013-12-04 湖南大学 一种基于熵格子波尔兹曼模型的并行流体计算实现方法
WO2018026327A1 (en) * 2016-08-04 2018-02-08 Nanyang Technological University An apparatus for transferring heat from a heat source to a heat sink
CN106650133A (zh) * 2016-12-28 2017-05-10 中南大学 加热炉内导热流体仿真方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
寇伟元;白志刚;: "基于格子玻尔兹曼方法的多型绕流数值模拟", 天津理工大学学报, no. 02 *
邵卫东;李军;: "格子玻尔兹曼方法中的特征无反射边界条件研究", 工程热物理学报, no. 02 *
黄兵方;闻炳海;邱文;赵琬玲;陈燕雁;: "基于晶格Boltzmann方法的接触角实时测量研究", 广西师范大学学报(自然科学版), no. 01 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112687002A (zh) * 2021-03-15 2021-04-20 四川省公路规划勘察设计研究院有限公司 一种三维地质模型网格优化方法

Also Published As

Publication number Publication date
JP2019215848A (ja) 2019-12-19
JP7296216B2 (ja) 2023-06-22
EP3528146A1 (en) 2019-08-21
US20190258764A1 (en) 2019-08-22

Similar Documents

Publication Publication Date Title
CN110175342A (zh) 用于高速流的基于格点玻尔兹曼的求解器
JP6434512B2 (ja) ハイブリッド熱格子ボルツマン法のための温度結合アルゴリズム
JP6561233B2 (ja) 等方性及びガリレイ不変性を強化した格子ボルツマン衝突演算子
JP6562307B2 (ja) 層流乱流遷移のモデル化を含む物理的プロセスのコンピュータシミュレーション
US11847391B2 (en) Computer system for simulating physical processes using surface algorithm
JP2010500654A (ja) 物理プロセスのコンピュータシミュレーション
JP7334125B2 (ja) 任意座標系のメッシュにおける物理的流体のコンピュータシミュレーション
US20220188485A1 (en) Lattice Boltzmann Solver Enforcing Total Energy Conservation
JP2021072123A (ja) スカラ輸送についてのガリレイ不変を強制する格子ボルツマンベースのスカラ輸送を使用して物理的プロセスをシミュレートするためのコンピュータシステム
JP7496049B2 (ja) 全エネルギー保存を実施する格子ボルツマンソルバ
JP7402125B2 (ja) 陽的数値拡散問題を安定化させる不規則空間グリッドによる物理流体のコンピュータシミュレーション
Hertel et al. Using a moving mesh PDE for cell centres to adapt a finite volume grid
Yang et al. A fast iterated orthogonal projection framework for smoke simulation
Hong et al. Intuitive Control of Deformable Object Simulation Using Geometric Constraints.

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
TA01 Transfer of patent application right
TA01 Transfer of patent application right

Effective date of registration: 20240319

Address after: Massachusetts USA

Applicant after: Dassault Systemes Americas Corp.

Country or region after: U.S.A.

Address before: Rhode Island USA

Applicant before: Dassault Systemes Simulia Corp.

Country or region before: U.S.A.