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

摘要

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

说明书

技术领域

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

背景技术

目前，LBM方法作为一种严格物理真实的计算流体动力学仿真方法，在多个流体仿真 领域都有着重要的应用。等温环境或无热环境下LBM方法已经有了大量成熟且稳定的研究 和应用，而热流体动力学(Thermohydrodynamics)作为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)的LatticeBhatnagar-gross-krook (LBGK)控制方程和表示热流体运动的热能方程，采用n维离散空间的m速度(DnQm) 多维离散网格模型，进行流场建模。建立流体仿真物理模型，包括流体求解区域、边界条件 和初始条件，流体边界条件有动态开放边界和封闭边界两种类型，初始条件包括流体密度、 初始速度、温度、流体类型等。连续的Boltzmann方程表示为：

$\frac{f}{t}+\xi ·▿f=\Omega \left(f\right)---\left(1\right)$

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

$J_{BGK}\left(f\right)=\frac{1}{{\tau }_c}[f^{eq}-f]---\left(2\right)$

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

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

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

$f_i^{\left(eq\right)}={\omega }_i\rho [1+\frac{c_i·u}{{\mathrm{RT}}_0}+\frac{{(c_i·u)}^2}{2{\left({\mathrm{RT}}_0\right)}^2}-\frac{u^2}{2{\mathrm{RT}}_0}]---\left(3\right)$

$h_i^{\left(eq\right)}={\omega }_i{p}_0[\frac{c_i·u}{{\mathrm{RT}}_0}+{\left(\frac{c_i·u}{{\mathrm{RT}}_0}\right)}^2-\frac{u^2}{{\mathrm{RT}}_0}-\frac{1}{2}(\frac{c_i^2}{{\mathrm{RT}}_0}-D)]+{\mathrm{Ef}}_i^{\left(eq\right)}---\left(4\right)$

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

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

$f_i(x+c_i{\Delta }_t,t+{\Delta }_t)=f_i(x,t)-{\omega }_f[f_i(x,t)-f_i^{eq}(x,t)]+{\Delta }_t(1-\frac{{\omega }_f}{2})F_i(x,t)---\left(5\right)$

$h_i(x+c_i{\Delta }_t,t+{\Delta }_t)-h_i(x,t)=-{\omega }_h[h_i(x,t)-h_i^{eq}(x,t)]---\left(6\right)$

其中，ω_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₀，可得：

$R=\frac{pV}{nT}=\frac{p_0}{{\mathrm{\rho T}}_0}---\left(8\right)$

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

$p=\frac{p_0}{{\mathrm{\rho T}}_0}\frac{nT}{V}---\left(9\right)$

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

$p=\alpha ·\frac{p_0}{T_0}T---\left(10\right)$

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

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

$F_i(x,t)={\omega }_i\rho [\frac{c_i·a}{{\mathrm{RT}}_0}+\frac{(c_i·a)(c_i·u)}{{\left({\mathrm{RT}}_0\right)}^2}-\frac{a·u}{2{\mathrm{RT}}_0}]---\left(11\right)$

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

$\left|a\right|=ds·p=ds·\alpha ·\frac{p_0}{T_0}T---\left(12\right)$

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

$ds·\alpha ={\omega }_i·4\pi ·{\left(\frac{\Delta x}{2}\right)}^2---\left(13\right)$

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

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

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

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

ρ＝∑f_i

ρu＝∑c_if_i(14)

$\rho E=\underset{i}{\Sigma }{h}_i$

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

E＝c_vT+u²/2(15)

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

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