CN105278346B - 一种基于离散格子Boltzmann双分布模型的热流体仿真方法 - Google Patents

Abstract

本发明公开了一种基于离散格子Boltzmann双分布模型的热流体仿真方法，涉及基于物理的流体仿真技术领域，以实现包含热传导的流体仿真为目的，基于物理真实的离散格子Boltzmann方法(LBM)的双分布模型，在LBM的流体仿真基础上，围绕含热传导的流体仿真展开研究。首先利用基于欧拉网格的LBM方法进行流场建模，构建流体仿真的主体部分，然后利用无黏性热耗散和压缩功的LBM双分布函数计算流体的热传导，并通过提出的耦合算法计算LBM双分布函数中热量与动量的相互转化，进而可以对诸如风的产生、热流体流动、流体热交换等基于物理的热流体现象进行计算机仿真。

Description

一种基于离散格子Boltzmann双分布模型的热流体仿真方法

技术领域

本发明涉及热流体仿真的技术领域，具体涉及一种基于离散格子Boltzmann双分布模型的热流体仿真方法。

背景技术

目前，LBM方法作为一种严格物理真实的计算流体动力学仿真方法，在多个流体仿真领域都有着重要的应用。等温环境或无热环境下LBM方法已经有了大量成熟且稳定的研究和应用，而热流体动力学(Thermo hydrodynamics)作为LBM领域的一个基础问题，一直以来都受到众多学者的青睐，许多学者一直在探索数值稳定、精确性好、算法结构简单且能反映热流动基本特性的方法，然而至今未有一个健全的模型能够完美的解决热流动问题。

一般来说，热流动的LBM模型可以分为三类：多速模型，双分布模型和近些年最新提出的与差分方法相结合的混合模型。多速模型是等温模型的直接推广，密度、速度、温度均由速度分布函数f_i求解得到。然而，多速模型往往需要使用比等温模型更多的离散速度来得到温度的宏观演化方程。与等温模型相比，多速模型的平衡态分布函数包含更高阶的速度项。早先多速模型的一个主要缺点是模型的Prandtl数不可调，此后虽有学者从不同角度发展了Prandtl数可调的模型，但热量方程中的黏性耗散项的输运系数可能不正确。多速模型的另一个不足是数值稳定性差，这严重制约了该类模型的普遍应用。

比较而言，后来提出的双分布函数LBM模型和混合方法可有效克服多速模型的两个基本缺点。双分布模型的基本思想是使用两个分布函数，分别用于速度场和温度场的描述。此类模型具有良好的数值稳定性，而且可以模拟比多速模型温度变化范围更大的流场，同时格子结构也比较简单。双分布模型按照是否可以耦合又分为两种：非耦合的双分布函数模型和耦合的双分布函数模型。非耦合的双分布函数模型中速度分布函数影响热量或温度的分布函数，而热量或温度分布函数并不影响速度分布函数，反映在流动中即为流场影响温度场，而温度场不反作用于流场，这是非耦合双分布模型的主要缺点。而耦合的双分布函数模型常用于不可压缩流的模拟，通过Boussinesq假设将速度分布函数与热量分布函数耦合起来。何雅玲、王勇在著作<<格子Boltzmann方法的理论及应用>>中将双分布函数模型和多速模型结合起来，提出了可用于可压缩流动的耦合的双分布函数模型，该模型满足完全气体状态方程，并且可以用于可压缩流体。该模型的缺点是模型相较于双分布模型的建模过程过于复杂，计算量很大，且容易继承多速模型数值稳定性差的缺点。

混合模型作为一种新型的LBM热流体模型，与双分布模型比较类似，对热能方程也是单独进行处理，但是使用的方法是差分方法求解，混合模型的显著优点是可以使压力与温度耦合起来，是一种耦合方法。

发明内容

本发明要解决的技术问题是：克服了现有模拟热流体运动模型的缺点与不足，提出了一种基于离散格子Boltzmann双分布模型的热流体仿真方法，方法在双分布模型框架内实现了热能与动能的耦合，解决了非耦合模型不能实现热能向动能转化的缺点，同时比现有耦合模型具有更小的复杂度和计算量，满足了物理真实的热流体运动的模拟仿真。

本发明采用的技术方案为：一种基于离散格子Boltzmann双分布模型的热流体仿真方法，包括以下四个步骤：

步骤(1)、仿真流体的建模：利用求解纳维-斯托克斯方程(Navier-Stokes，N-S)的Lattice Bhatnagar-gross-krook(LBGK)控制方程和表示热流体运动的热能方程，采用n维离散空间的m速度(DnQm)多维离散网格模型，进行流场建模。建立流体仿真物理模型，包括流体求解区域、边界条件和初始条件，流体边界条件有动态开放边界和封闭边界两种类型，初始条件包括流体密度、初始速度、温度、流体类型等；

步骤(2)、求解LBM速度分布函数：双分布函数模型包括两个部分，分别表征流体流动的速度分布和热量传导的热能分布，求解LBM双分布函数可得到流体流动的速度变化和热能变化。在一个时间步长Δt内，根据步骤(1)中的初始条件首先求解流体的速度分布函数，得到流体的速度分布；

步骤(3)、求解LBM热能分布函数：时间步长Δt不改变，利用LBM双分布模型的热能分布函数，实现步骤(2)中速度分布向热能分布的转化。首先，由步骤(2)得到流体在仿真流场中的速度分布，根据流体的速度分布确定流体的速度矢量和密度值；然后，根据得到的流体速度矢量、密度等求解热能分布函数，得到流体的热能分布；接下来，根据所得到的热能分布确定在下一个时间步长Δt+1时由热能向动能的转化量大小；

步骤(4)、时间步长Δt+1，进入下一个仿真的时间步内，重复步骤(2)、(3)，实现包含热传导的流体的连续仿真。

本发明的原理在于：

本发明提出了一种单纯基于LBM双分布模型的耦合方法，耦合性可理解为：热量(本发明中针对总能形式)分布函数演化方程中存在的速度分布函数项，以及平衡态热量分布函数中存在的密度和速度项体现了流场对温度的影响；而温度场对流场的影响则是通过气体状态方程发生的，本发明推导了由气体状态方程向压强项转化的物理过程，提出了由温度向压力项的转化公式，使得LBM双分布模型中的速度分布函数中含有与温度有关的作用力项。对可压缩Navier-Stokes方程组来说，假定压力给定，那么连续方程和动量方程4个标量方程含4个未知量(密度及3个速度分量)，方程组封闭，密度场、速度场可以确定，而后温度场亦可以确定。于是，可以通过状态方程由密度和温度确定压力场，并开始新一轮的演化。

对比单分布函数可压缩模型与耦合的双分布函数可压缩模型，从恢复同一层次的宏观方程来看，如果确定平衡态分布函数的方法相同，则由于双分布函数对约束条件的降阶，双分布函数模型的平衡态分布函数通常较单分布函数模型中的简单。从调节比热容比和普朗特数来看，双分布函数模型具有较大的优势。在双分布函数模型中，只需要在平衡态能量分布函数中引入一个与比热容比相关的参数便可以调节比热容比，这比单分布函数模型中通过采用势能或者多能级来调节比热容比方便不少。

本发明与现有技术相比的优点在于：

1、本发明提出的基于离散格子Boltzmann双分布模型的热流体仿真方法，应用于计算机动画和计算机仿真领域，相对于现有的流体仿真方法物理更加真实，算法更加严谨、精确。

2、本发明的方法可以在可压缩热流体仿真中实现热能向动能的转化，进而对诸如空气流通(如风的生成)、流体热传导等场景进行模拟。

3、相较于现有计算流体动力学领域的LBM热流体仿真方法的非耦合模型，本发明的方法可以实现热能和动能的双向转化，弥补了现有非耦合方法不能实现热能向动能转化的不足。

4、相较于现有计算流体动力学领域的LBM热流体仿真方法的耦合模型，在满足了物理真实的同时，模型更加简洁，具有更小的算法复杂度和计算量。

附图说明

图1为本发明一种基于离散格子Boltzmann双分布模型的热流体仿真方法流程图；

图2为二维LBM示意图；

图3为本发明的三维热能分布示意图；

图4为本发明的二维热能分布示意图；

图5为热能向动能转化结果示意图；

图6为本发明提出的双分布模型三维仿真结果。

具体实施方式

图1给出了一种离散格子Boltzmann双分布模型的热流体仿真方法的总体处理流程，下面结合其他附图及具体实施方式进一步说明本发明。

本发明提出供一种离散格子Boltzmann双分布模型的热流体仿真方法，具体实施为基于双分布LBM模型热流体仿真，主要步骤介绍如下：

1、热流体仿真先期建模

利用求解纳维-斯托克斯方程(Navier-Stokes，N-S)的Lattice Bhatnagar-gross-krook(LBGK)控制方程和表示热流体运动的热能方程，采用n维离散空间的m速度(DnQm)多维离散网格模型，进行流场建模。建立流体仿真物理模型，包括流体求解区域、边界条件和初始条件，流体边界条件有动态开放边界和封闭边界两种类型，初始条件包括流体密度、初始速度、温度、流体类型等。连续的Boltzmann方程表示为：

ξ是粒子速度，Ω(f)是表示碰撞过程的碰撞因子，▽是哈密顿算子。f是连续方程的平衡态分布函数，平衡态分布函数保证了流体局部的质量与动量守恒。f线性依赖于流体密度ρ和速度u。根据上式，可通过求数值解的方法模拟流体的宏观运动。然而碰撞项Ω(f)并非分布函数的线性方程，只与分子作用力有关，因此求解起来具有难度，这使得Boltzmann方程在实际应用中受到了很大限制。改进了的Bhatnager-Gross-Krook(BGK)模型假设流体运动具有平衡态，分子之间的碰撞过程会促使分布函数ff趋近于平衡态。由于碰撞导致的变化量和ff与平衡态的差值成正比，即：

τ_c为松弛时间，1/τ_c称之为平均碰撞频率。

从宏观运动的角度，粒子速度ξ和分布函数ff也相应的被离散到n维空间，这意味着在t时刻，某个格点x处将有沿各个离散方向的分布函数。格子Boltzmann方程(LBE)是BGK方程的特殊离散格式，从空间、时间和速度的角度将连续矢量离散成沿各个特定的方向的标量值。以二维流场为例，LBM模型示意图如图2所示。

本发明涉及表示热流体运动的LBM双分布函数模型，此模型中，流体运动被看作为两个步骤——碰撞(Collision)和迁移(Stream)。其算法核心包括两组方程，平衡态分布函数(也叫平衡态方程)和碰撞方程，平衡态方程如下：

其中，表示平衡态的速度分布函数，表示平衡态的热能分布函数，eq是equilibrium的缩写，i表示离散速度的序号，ω_i是权系数，c_i表示离散速度，u表示网格格点处的宏观速度，R是气体普适常量，T₀是参考温度，p₀是压强，p₀＝ρRT₀，若仿真流体为气体时，D表示气体分子的自由程，E表示格点处的内能。LBM模型的基本思想认为流体在运动过程中有趋向于平衡态分布的趋势，因此平衡态分布函数是求解LBM流体的核心。

另外，流体运动过程中的“碰撞”过程直接决定了流体的运动状态和运动规律，碰撞方程是LBM双分布模型的第二组核心方程，表示如下：

其中，ω_h＝2Δ_t/(2τ_h+Δt)，ω_f＝2Δ_t/(2τ_f+Δ_t)，τ_f是速度的松弛时间，τ_h是内能的松弛时间。x是当前格子位置，t表示当前时间。流场不受外力的影响，只跟流场内部热量转化有关，在速度分布函数中增加一项系统内力F_i(x，t)。

2、热能向动能转化

求解LBM热能分布函数的过程中，速度分布函数的平衡态对热量分布函数的平衡态求解有直接影响，如方程(3)、(4)所示，与成正相关，而的求解并不依赖于由此根据气体状态方程建立热能向动能转化的关系式。具体过程如下：

理想气体状态方程为：

pV＝nRT (7)

其中，p为气体压强，n为物质的量，T为气体温度，V是气体体积，R是气体普适常量。由权利要求2已知在热量的平衡态分布函数中，压强p₀＝ρRT₀，可得：

由(7)、(8)得到压强p的表达式：

对于基于LBM方法的流体仿真而言，气体体积V表示格点处抽象流体粒子团的体积，而n表示格点处流体粒子的物质的量，因此有ρV·α＝n，α是与仿真流体规模相关的可控参量，压强表达式简化为：

由上述推导过程可知，流场内流体的速度分布可以通过温度的改变量在迁移过程中施加影响。

考虑到流场不受外力的影响，只跟流场内部热量转化有关，因此压强项最终可转化为方程(5)中的作用力项F_i(x，t)。F_i(x，t)的表达式为：

上式中只有a为未知力，a的方向与速度u方向一致，大小为|α|＝ds·p，ds为压力作用的微观面积，则：

而∑_ids·α可以近似看作为仿真空间中以格子步长Δx/2作为半径的球面积，那么：

最后将α带入方程(11)，由此，以上过程建立起完整的由热能向动能转化的物理机制。

3、求解热流体运动方程步骤

求解LBM速度分布函数：双分布函数模型包括两个部分，分别表征流体流动的速度分布和热量传导的热能分布，求解LBM双分布函数可得到流体流动的速度变化和热能变化。求解顺序如下：

在一个时间步长Δt内，根据初始模型求解流体的速度分布态分布函数(3)，得到流体的速度平衡态分布函数；将得到的速度平衡态分布函数带入方程(4)，得到热能平衡态分布函数；然后根据得到的流体速度矢量、密度等求解热能分布函数，得到流体的热能分布；接下来，根据所得的热能分布确定由热能向动能的转化量大小，即求解F_i(x，t)；再然后求解速度分布函数和热能分布函数的迁移，即求解方程(5)、(6)。根据求得的分布函数f_i和h_i，计算宏观量，得到宏观密度、速度和内能，宏观量的累加过程如下：

ρ＝∑f_i

ρu＝∑c_if_i (14)

根据所得内能E与温度T的关系即可求得温度场：

E＝c_vT+u²/2 (15)

其中，c_v为所求流体的定容比热。最后，时间步长Δt+1，进入下一个仿真的时间步内，重复以上步骤，实现包含热传导的流体的连续仿真。

如图3-6所示，为证明本方法在计算机动画和计算机仿真领域的正确性，设计了一个封闭边界的三维流场场景，底端正中心有一个占底部面积1/4的开放边界作为流体的输入，热能或动能都通过此开放边界影响流场。图3表示输入热能和动能，双分布模型热传导的热能分布示意图，图4表示输入热能和动能，在三维温度场中截取一个平面后的热能分布示意图。通过这两个图可以明显看出，本发明的方法能够真实有效的对热流体运动过程中的热传导进行仿真。图5表示静止流场内，在热能驱动下的空气对流产生，即实现了诸如热流驱动产生风的现象模拟。图6是本发明方法涉及方法的综合仿真结果。需要说明的是，图3、图4、图6示意图中颜色深浅代表内能的大小。图5中矢量长短代表速度的大小。

Claims (1)

1.一种基于离散格子Boltzmann双分布模型的热流体仿真方法，其特征在于包括以下四个步骤：

步骤(1)、仿真流体的建模：利用求解纳维-斯托克斯方程(Navier-Stokes，N-S)的Lattice Bhatnagar-gross-krook(LBGK)控制方程和表示热流体运动的热能方程，采用n维离散空间的m速度(DnQm)多维离散网格模型，进行流场建模，建立流体仿真物理模型，包括流体求解区域、边界条件和初始条件，流体边界条件有动态开放边界和封闭边界两种类型，初始条件包括流体密度、初始速度、温度和流体类型；

步骤(2)、求解LBM速度分布函数：双分布函数模型包括两个部分，分别表征流体流动的速度分布和热量传导的热能分布，求解LBM双分布函数可得到流体流动的速度变化和热能变化，在一个时间步长Δt内，根据步骤(1)中的初始条件首先求解流体的速度分布函数，得到流体的速度分布；

步骤(3)、求解LBM热能分布函数：时间步长Δt不改变，利用LBM双分布模型的热能分布函数，实现步骤(2)中速度分布向热能分布的转化，首先，由步骤(2)得到流体在仿真流场中的速度分布，根据流体的速度分布确定流体的速度矢量和密度值；然后，根据得到的流体速度矢量、密度求解热能分布函数，得到流体的热能分布；接下来，根据所得到的热能分布确定在下一个时间步长Δt+1时由热能向动能的转化量大小；

步骤(4)、时间步长Δt+1，进入下一个仿真的时间步内，重复步骤(2)、(3)，实现包含热传导的流体的连续仿真；

步骤(1)中所述的流体建模过程，该方法使用空间场内的均匀网格表示流场，利用DnQm网格模型表示流体在网格格点处的属性，将纳维－斯托克斯(Navier-Stokes，N-S)和流体的热传导方程转化为只与时间相关的离散形式的LBM双分布函数形式求解，基于物理的LBM具体过程为假设流体运动具有平衡态，分子之间的碰撞过程会促使流体速度和热能分布函数都趋近于平衡态，在流体趋于平衡态的过程中，从空间、时间和速度的角度将连续矢量离散成沿各个特定的方向的标量值，将连续的流体运动看作为两个离散步骤——碰撞(Collision)和迁移(Stream)，通过离散时间Δt驱动流体在网格格点处的碰撞或迁移过程以实现流体模拟过程， 碰撞过程用平衡态分布函数求解，以含黏性热耗散和压缩功的双分布函数的总能形式为例，碰撞过程的双分布函数具体表示如下：

其中，f_i ^(eq)表示平衡态的速度分布函数，表示平衡态的热能分布函数，eq是equilibrium的缩写，i表示离散速度的序号，ω_i是权系数，c_i表示离散速度，u表示网格格点处的宏观速度，R是气体普适常量，T₀是参考温度，p₀是压强，p₀＝ρRT₀，若仿真流体为气体时，D表示气体分子的自由程，E表示格点处的内能；

在步骤(3)中所述求解LBM热能分布函数的过程中，速度分布函数的平衡态f_i ^(eq)对热量分布函数的平衡态求解有直接影响，具体体现为与Ef_i ^(eq)成正相关，而f_i ^(eq)的求解并不依赖于由此根据气体状态方程建立热能向动能转化的关系式，具体过程如下：

理想气体状态方程为：

pV＝nRT

其中，p为气体压强，n为物质的量，T为气体温度，V是气体体积，R是气体普适常量，在热量的平衡态分布函数中，压强p₀＝ρRT₀，可得：

进而得到压强p的表达式：

对于基于LBM方法的流体仿真而言，气体体积V表示格点处抽象流体粒子团的体积，而n表示格点处流体粒子的物质的量，因此有ρV·α＝n，α是与仿真流体规模相关的可控参量，压强表达式简化为：

由上述推导过程，可知流场内流体的速度分布可以通过温度的改变在迁移过程中施加影响，迁移过程由显式的LBM方程求解，可用f_i和h_i表示，具体方程如下：

其中，ω_h＝2Δ_t/(2τ_h+Δ_t)，ω_f＝2Δ_t/(2τ_f+Δ_t)，τ_f是速度的松弛时间，τ_h是内能的松弛时间，x是当前格子位置，t表示当前时间，考虑到流场不受外力的影响，只与流场内部热量转化有关，在速度分布函数中增加一项系统内力F_i(x，t)，F_i(x，t)不是流场所受外力，而是内能向动能转化时影响速度分布的流场内力，F_i(x，t)的表达式为：

上式中只有a为未知力，a的方向与速度方向一致，大小为|a|＝ds·p，ds为压力作用的微观面积，则：

而∑_ids·α可以近似看作仿真空间中以格子步长Δx/2作为半径的球面积，那么

由此，以上过程建立起完整的由热能向动能转化的物理机制。