LBM相变传热与流体流动数值分析13

相变传热与流体流动数值分析(第11-14讲)格子Boltzmann方法LatticeBoltzmannMethod主要内容

13.1

LBGK模型13.2

MRT模型

13.3多相和多组分模型13.4

LBM热模型概述本章重点；LBM根本模型

连续Boltzmann方程LBGK方程BGK近似泰勒级数展开Chapman-Enskog回归上章回忆；理论推导单相模型(D2Q9模型)多相和多组分模型〔SC模型〕根据不同的划分规那么，LBM的模型包括：单松弛(LBGK)模型

等温模型非等温模型多松弛(MRT)模型不同的松弛时间温度是否变化?

单组分单相模型多组分多相模型多相，多组分?格子Boltzmann-BGK〔LBGK〕模型的演化方程：13.1

LBGK模型代表性的模型：Qian等人1992年提出的DnQb模型和Guo等人2006年提出的DnGb模型。〔n为空间维数，b是离散速度数〕LBGK模型单松弛模型DnQb模型DnGb模型伪不可压模型不可压模型012345678e1e2e3e4e5e6e7e8f1f2f3f4f5f6f7f8f0Histogramviewofthedistributionfunction,f.f1f5125aff2f3f434678f6f7f8Unfold13.1

LBGK模型DnQb模型13.1

LBGK模型碰撞过程：迁移过程：DnQb模型平衡态分布函数:这两个参数是决定LBGK

模型的关键参数，其取值依赖于选用的格子类型。DnQb模型13.1

LBGK模型D1Q3:D1Q5:DnQb模型13.1

LBGK模型格子速度D2Q7:D2Q9:DnQb模型13.1

LBGK模型D3Q15:DnQb模型13.1

LBGK模型D3Q19:DnQb模型13.1

LBGK模型D3Q19:D2Q9:DnQb模型13.1

LBGK模型通过和Chapman-Enskog

展开方法类似的多尺度分析方法，对演化方程进行处理，可以得到LBGK

模型对应的宏观方程如果流体的密度变化不大，即r≈r0，上式就是标准的不可压N-S方程。Note：在推导过程中，假设Mach数充分小，因次LBGK模型的上述结果仅限于低Mach数流动〔Ma≤0.3〕。DnQb模型13.1

LBGK模型对于宏观物理量，可由以下方程得到实际上,LBGK方法势求解不可压N-S方程的一种人工压缩方法。DnQb模型13.1

LBGK模型压力入口压力出口

D2Q9模型----Poiseuilleflow13.1

LBGK模型多松弛模型和单松弛模型的主要区别在于它的碰撞过程包含多个松弛时间〔Multiple-Relaxation-Time〕。MRT模型含有更多的可调参数，剪切黏性和体黏性都是可调的；MRT可以克服普通LBGK模型的明显缺乏，如不受因单松弛参数影响而保证Prandtl数固定不变的限制；MRT模型有更好的稳定性，并且在粘性较高时可减小人工压缩性的影响。优点：缺点：

公式复杂，计算量较大。碰撞矩阵13.2

MRT模型〔〕上式描述了离散速度分布函数在速度空间的时空演化过程。通过b个b维基向量，可以建立速度空间和矩空间之间的关系，即其中，是粒子速度的多项式函数，且线性无关。定义b个矩，其中，即，13.2

MRT模型关键问题：确定碰撞矩阵演化过程，仍然分为碰撞和迁移两步：碰撞：迁移：因为碰撞步只涉及局部计算，写成矢量形式上式两端乘以变换矩阵M，得其中，为矩空间的平衡态函数13.2

MRT模型以标准D2Q9模型为例：变换矩阵：对应的矩为矩空间的平衡态其中a和b为自由参数，当a＝1，且b＝－3，时，该平衡态函数和LBGK模型的平衡态分布函数一致。13.2

MRT模型模型的松弛参数为：剪切黏性和体黏性系数分别为：其中：13.2

MRT模型Multi-ComponentMultiphaseMiscibleFluids/Diffusion(NoInteraction)ImmiscibleFluidsSingleComponentMultiphaseSinglePhase(NoInteraction)NumberofComponentsInteractionStrengthNatureofInteractionAttractiveRepulsive13.3多相和多组分模型颜色模型〔Chromodynamicsmodel〕自由能模型〔Freeenergymodel〕伪势模型〔SCmodel，pseudo-potentialmodel)

基于动理学理论的LBM模型13.3多相和多组分模型Gunstensen等人于1991年提出的第一个多组分多相模型。基于Rothman和keller提出的LGA两相流模型。Grunau等又将此模型推广到密度和粘性变化的两相流系统中。

用不同颜色区分不同相态的流体；不同流体之间的相互作用通过引入颜色梯度来实现，并根据它来调整流体粒子的运动趋势，实现流体的别离或混和。特点：颜色模型

以两相流为例：引入两个分布函数，分别表示红色相和蓝色相流体。混和流体的分布函数的演化方程为其中，表示由流体粒子之间碰撞引起变化，可以用BGK模拟，而表示界面张力引起的扰动。每相和混和流体的宏观流动变量为颜色模型ciWl是碰撞算子的与黏性系数相关的特征值。

所以定义为其中：，A是控制表面张力s的参数，而

显然，在单相区内颜色梯度为0，因此，也就是说，表面张力引起的分布函数变化只在界面处起作用，这是符合物理事实的。引入描述两相差别的局部序参数（orderparameter）并定义局部的颜色梯度颜色模型确定了后，混和分布函数在每一个时间步内根据方程〔〕进行演化，红蓝两相的分布函数那么通过重新标色〔Recoloring〕过程来获得。这一过程是求解下面的极大化问题并满足如下的限制条件其中，是重新标色后的分布函数，而是碰撞前的量。这一极大化的实质是使颜色通量颜色梯度G的方向尽可能一致，并包装每个格点处各相质量及每一方向的质量守恒。颜色模型综上，Gunstensen等人的颜色模型可分为三步进行演化：1.单相碰撞：a)计算混和分布函数；

b)计算和平衡态分布函数；

c)按照单相碰撞规则计算碰撞后分布函数。2.两相碰撞：a)将表面张力扰动加入碰撞后分布函，；

b)对重新标色，得到。3.迁移：NOTE上述模型演化过程中的碰撞对象是混和分布函数！颜色模型

其中分别为每一相的碰撞算子，采用BGK模型颜色模型的改进：Grunau等提出了改进的颜色模型，允许不同颜色的分布函数参与碰撞过程，其演化方程为：

其中平衡态分布依赖每相的密度。颜色模型仍然是外表张力引起的扰动，形式为

其中是与表面张力相关的参数，是与模型相关的常数。

重新标色过程如下：首先通过极大化，使红色流体的动量与颜色梯度方向尽可能一致；之后，令蓝色流体的分布函数为,

其中是重新标色前的混和分布函数。颜色模型

局限性：外表张力与界面走向相关〔各向异性〕；在相界面处附近会产生非物理现象；“颜色能量”的极小化过程的计算量较大；不容易考虑热动力学影响；颜色模型1993年和1994年，Shan和Chen提出了一种能直接地刻画粒子之间相互作用的LBE模型。伪势模型单组分多相模型(SCMP)(e.g.H2O,liquidandvapor)

多相多组分模型(MCMP)(e.g.H2Oandoil)SCMP模型

伪势模型TheprincipaldistinguishingcharacteristicofsinglecomponentmultiphaseLBMsistheincorporationofanattractiveforcebetweenfluid‘particles’.理想气体状态方程：Pispressure(ATM)Visvolume(L)nisnumberofmolsRisgasconstant(0.0821Latmmol-1K-1)Tistemperature(K)if1/Vm(molL-1)isdensity,cs2=RT‘a’termduetoattractiveforcesbetweenmolecules[atmL2mol-2]‘b’termduetofinitevolumeofmolecules[Lmol-1]vanderWaalsEOS:伪势模型vanderWaalsgaslaw(EOS):Vmisthevolumeoccupiedbyonemoleofsubstance.SCMP模型

CO2:P-VSpace伪势模型SCMP模型

Liquid-VaporCoexistence:CO2,P-VSpaceLiquidMolarVolumeVaporMolarVolumeVaporPressurea=3.592L6atmmol-2b=0.04267L3mol-1伪势模型SCMP模型

IncludingDirectionalH-Bonds:Water,298K:P-rSpaceSpinodalSCMP模型

伪势模型AnattractiveforceFbetweennearestneighborfluidparticlesisinducedasfollows:伪势模型SCMP模型----D2Q9Fluid-FluidInteraction假设流体粒子之间存在非局部的相互作用，相应的势函数为：根据相互作用势的定义，可以得到两相间的作用力为GistheinteractionstrengthYistheinteractionpotential:

Y0andr0arearbitraryconstantsG<0forattractionbetweenparticles.Forceisstrongerwhenthedensityishigherwiis1/9fori=1,2,3,4andis1/36fori=5,6,7,8假设只考虑相邻节点间的作用力伪势模型SCMP模型----D2Q9Fluid-FluidInteraction伪势模型SCMP模型----D2Q9ForceimplementationNon-idealComponent伪势模型SCMP模型----D2Q9PressureexpressionRealisticEOSforwater:FollowsidealgaslawatlowdensitycompressibilityofwaterathighdensityNorepulsivepotentialinLBmodel伪势模型SCMP模型----D2Q9Pressureexpression伪势模型SCMP模型----D2Q9Surfaceadhesion(adsorption)润湿性Martys等提出Gwistheinteractionstrength.“—”or“＋”sisa‘switch’thattakesonvalue1ifthesiteatx+eiDtisasolidandis0otherwiseWeseemtohaveflexibilityinthechoiceofthepre-sumfactor;thepapersciteduserorY伪势模型SCMP模型----D2Q9Surfaceadhesion(adsorption)伪势模型SCMP模型----D2Q9Surfaceadhesion(adsorption)如果有外力作用，如重力等伪势模型SCMP模型----D2Q9Phaseseparationr伪势模型SCMP模型----D2Q9Young-Laplaceequation(1805)SCMP模型----D2Q9SurfaceTensionmeasurement)伪势模型其中，k

表示不同的组分。伪势模型中假设流体粒子之间存在非局部的相互作用，相应的势函数为其中为Green函数，决定了组分之间相互作用的强度。与组分的密度有关，表示组分k

的有效密度。不同的表达方式，所得的状态方程不同。伪势模型MCMP模型

参数的绝对值决定了组分之间相互作用的强度，符号决定了两者之间是相互吸引还是相互排斥，即负值是吸引力而正值为排斥力。同样地，假设只考虑临近格点的影响，例如，在D2Q9中根据相互作用势，可得到k类流体粒子受到的粒子间的作用力为伪势模型伪势模型MCMP模型

同样的，通过改变平衡态速度来实现的影响，即混合流体的速度混合流体的密度随后，Shan和Doolen对根本模型进行了改进，重新定义了保证没有作用力时碰撞过程动量守恒，即伪势模型伪势模型MCMP模型

在Shan-Doolen模型中，混合流体的速度为宏观流体变量，可由下面方程确定压力黏度伪势模型伪势模型MCMP模型

参数的绝对值决定了组分和壁面之间相互作用的强度。符号决定了作用力是吸附还是排斥，即负值是吸引力，表明壁面对流体时润湿的，而正值为排斥力，表面壁面对流体而言是非润湿壁面。

s为壁面指数，当格子特性为固体时，其值为1，当格子特性为流体时，其值为0。伪势模型伪势模型MCMP模型----壁面润湿性

同样地，在D2Q9中壁面吸附力的影响，也通过改变平衡态速度来实现。当有外力作用在流体上时，eg.重力，同样可以通过改变平衡态速度的方法来实现伪势模型伪势模型MCMP模型----D2Q9forceimplementation

不同时刻的两组分的别离过程参数选择：伪势模型MCMP模型----D2Q9Phaseseparation伪势模型MCMP模型----D2Q9Phaseseparation

液滴在不同时刻的形状演变模拟条件：流体2被流体1环绕；不考虑外力；密度比和黏度比均为1；伪势模型MCMP模型----D2Q9Young-LaplaceLaw

处两组分沿x轴密度分布图处两组分沿x轴压力分布图压力差与曲率的关系界面张力与相间力参数的关系伪势模型MCMP模型----D2Q9Young-LaplaceLaw

相间力参数对最大的伪速度和互溶密度的影响时，稳态时速度分布图相间力参数的取值实际上受两方面的约束:伪速度越小越好，互溶密度越小越好!伪势模型MCMP模型----D2Q9Young-LaplaceLawLaplaceLaw伪势模型MCMP模型----D2Q9Interfacialtension(asopposedtosurfacetensionbetweenaliquidanditsownvapor)

伪势模型MCMP模型----D2Q9Poiseuille解析解流体266×132个格子；半步长反弹格式、周期边界；恒定的体积力；密度均为1；运动黏度均为0.16；

数值解与解析解的比较伪势模型MCMP模型----D2Q9Poiseuille优点：自动实现相别离，界面自动生成跟踪；直接对微观相互作用力进行描述，能够反映多相/多组分流体动力学的物理本质，应用比较很广泛。缺点：当相互作用力中的有效密度函数取指数形式时，该模型才与热力学相关理论一致。无法模拟具有大密度比和大粘度比的流体系统。〔有研究者指出改变有效密度的表达方式，可以模拟密度比为1000的情况，但是大粘度比〔M>10〕的流体系统仍然无法模拟〕。伪势模型颜色模型和伪势模型都是基于界面现象的唯象模型，含有一些人为的假设。1995年，Swift等人直接从多相/多组分流体的自由能理论出发，构造了与热动力学理论一致的多相和多组分流体模型。根本思想：根据自由能函数构造LBE的平衡态分布函数，通过引入一个非理想流体的热力学压力张量，使得系统的总能〔包括动能、内能和外表能〕守恒得以满足。单组分非理想流体模型两组分非理想流体模型自由能模型根据vanderWaals理论，对于包含相界面的单组分非理想流体系统，当处于局部平衡时，其自由能泛函为单组分非理想流体模型其中，是体相区自由能密度函数，是与外表张力相关的由密度梯度诱导的自由能。一个常用的密度梯度自由能形式为其中，是与外表张力相关的一个参数。自由能模型〔〕根据自由能泛函定义一个非局部压力单组分非理想流体模型其中，是状态方程，表示变分运算。由此可得包含界面奉献的热力学压力张量这里的不包含黏性应力部分。NOTE自由能模型〔〕为构造符合上述热力学理论的模型，Swift等建议在LBGK模型中使用一个包含密度梯度的平衡态分布函数单组分非理想流体模型其中满足条件为确定的具体形式，将其表示为u的多项式，例如对于FHP模型自由能模型根据限定条件式〔〕，可以确定平衡态分布函数的参数单组分非理想流体模型流体宏观速度为自由能模型经过Chapmann-Enskog分析可得式〔〕对应的宏观方程为单组分非理想流体模型其中动量方程右端最后一项在单相区为0，此时方程为一般的N-S方程。但是在相界面附近，密度梯度大，该项不能忽略。自由能模型使用两类分布函数（组分分布函数之和及之差）来描述系统的整体演化，即双组分非理想流体模型自由能模型平衡态分布函数和都包含流体的非均匀特性。仍然满足约束条件〔〕，但此时密度和速度都是混合流体的密度和速度。其中，是两组分的化学势之差，是迁移系数。这里热力学压力张力和化学势差可以根据两组分流体的自由能泛函得到。考虑相互排斥的两组分流体，自由能泛函为双组分非理想流体模型平衡态分布函数则满足以下条件自由能模型双组分非理想流体模型其中是表征相互作用强度的参数。当时，两相会发生分离。考虑相互排斥的两组分流体，自由能泛函为自由能模型双组分非理想流体模型其中根据这一自由能泛函，可得化学势差和热力学压力张力的张量表达式自由能模型双组分非理想流体模型根据上述的化学势差和热力学压力张力的张量表达式，由限制条件（13.3.8）可得平衡态分布函数和的表达式。对FHP模型，的展开系数为自由能模型双组分非理想流体模型相应的平衡态分布函数为其中自由能模型基于上述平衡态分布函数，可得FHP模型对应的宏观方程为双组分非理想流体模型其中和由式（13.3.4）给出（用演化方程中的松弛时间代替）。方程中的输运系数是。是演化方程中的松弛时间。自由能模型局限性：当界面附近存在较大密度梯度时，不满足伽利略不变性，会导致一些非物理现象。例如，在一个匀速运动的流场内，一个最初是圆形的液滴或气泡会随着时间的推移变成椭圆形。--------Swift等人，采用两个可调参数来改变平衡态分布函数的二阶距，可以局部克服这一局限性。--------Inamuro等人基于渐进分析方法，建立了一个满足伽利略不变性的D2Q9自由能模型。--------Kalarakis等人在Swift等人提出的7速模型根底上，通过修正平衡态分布函数系数，得到满足伽利略不变性的自由能模型。自由能模型局限性：无法模拟具有大密度比的流体系统--------Inamuro等人提出了改进的方法，可以模拟流体密度为1000的流体系统。主要思想是采用近似C-H(Cahn-Hilliard)方程的扩散方程来追踪相界面。--------Zheng等人基于自由能模型也提出了改进方法，可以模拟大密度比的流体系统。他们采用近似的C-H方程去追踪和定义界面，在界面的定义上没有人工的剪切。这在物理背景比其他方法更接近朗道的平均场理论。自由能模型13.3.4基于动力学的LBM模型基于Enskog方程的Luo模型简单流体的LBE模型可以从Boltzmann方程得到。同样，从多相或多组分流体的介观动理学出发，也可得到相应的LBE模型。多组分单相流体的LBE模型单组分多相流体的LBE模型基于有效作用力的He-Shan-Doolen模型不可压多相流体的He-Chen-Zhang模型基于Sirovich理论的两组分LBE模型基于拟平衡态理论的多组分LBE模型基于Hamel