相变传热与流体流动数值分析第114讲课件

1、相变传热与流体流动数值分析(第11-14讲)格子Boltzmann方法Lattice Boltzmann Method主要内容 14.1 启发式边界条件 14.2 动力学边界条件 14.3 外推边界条件 14.4 曲面边界条件概述 14.5 LBM模拟的基本步骤本章重点；边界条件处理 上章回忆；理论推导边界条件处理如何进行LBM模拟单相模型(D2Q9模型)多相和多组分模型（SC模型）14. LBM边界处理每个时步之后，内部流场节点上的分布函数均已经获得，但是边界节点上的部分分布函数是未知的。需要根据已知的宏观边界条件确定出边界节点上相应的分布函数的取值。根据边界处理格式的特性，对边界处理的方式

2、进行分类，主要分为启发式格式、动力学格式、外推格式以及其他复杂边界处理格式。原理： 根据边界上诸如周期性、对称性、充分发展等宏观物理学特性，通过微观粒子的的运动规则直接确定边界节点上的未知分布函数。优点： 不需要较复杂的数学推导和公式求解类型： 周期性边界、对称边界、充分发展以及用于固壁边界的反弹格式、镜面反弹格式、反弹与镜面反弹混合格式等14.1 启发式格式如果流场在空间呈现周期性变化或在某个方向无穷大，常常将周期性单元取出作为模拟区域，并在相应边界上采用周期性边界。其边界处理格式是指，当流体粒子从一侧边界离开流场时，在下一个时步就会从流场的另侧边界进入流场。周期性边界能严格保证整个系统的质

3、量和动量守恒。14.1.1 周期性边界处理格式周期格式可以表示为：14.1.1 周期性边界处理格式14.1.2 对称边界处理格式由于对称性问题，为了节省计算资源，可以取物理模型的一半作为模拟的区域并在对称轴上采取对称边界处理当流体在通道内流动达到充分发展后，密度与速度等物理量在主流方向上不再发生变化，即他们的空间导数为0其处理格式可以表示为：这是格子Boltzmann方法中处理充分发展边界最简单和最常用的方法14.1.3 充分发展边界处理格式处理静止无滑移壁面的一种常用格式，包括标准反弹格式、半步长反弹格式及修正反弹格式标准反弹格式：对于静止固体边界，常用的处理方法就是对边界上的粒子进行弹回处

4、理，即假设粒子与壁面碰撞后沿粒子原方向的反方向逆转14.1.4 反弹格式标准反弹格式： 操作简单，能够严格保证系统的质量和动量守恒； 一阶精度，降低了格子玻尔兹曼方法的精度； LBM格子波尔兹曼方法在内节点上具有二阶精度 -降低方法精度 修正反弹格式或者半步长反弹格式-二阶精度与标准反弹格式的最大区别：在碰撞过程中，边界节点与流体区域的节点一样，所有的离散速度方向的粒子均参与碰撞。得到的分布函数在进入下一次的碰撞和迁移的过程，这一格式称为修正反弹格式。修正反弹格式： 14.1.4 反弹格式 半步长反弹格式与标准反弹格式的反弹原理相同。不同的是：固体边界位于两个格点中间。 半步长反弹格式： 物理

5、意义： 在 格子节点处碰撞，速度指向壁面的粒子经过 时间到达固体壁面，然后与壁面碰撞后反弹，速度逆转，再经过 时间沿原路返回到 格子节点处。反弹边界的优点在于操作简单，对于处理复杂的不规则边界（如多孔介质）具有很大的优势;14.1.4 反弹格式对于光滑壁面，如果壁面对于流体没有摩擦作用，即壁面与流体之间没有动量交换，则通常采用镜面反射格式来实现自由滑移边界。流体节点(i-1,2)入射粒子 迁移到壁面节点(i,1)后，如果沿原路返回即 方向，即为反弹处理；如果以壁面法线方向 反射，则为镜面反射处理。 14.1.5 镜面反射格式壁面法向速度为零；切向速度分量不为零!适合既不能用简单的反弹格式处理，

6、也不能用镜面反射格式的情况，如微通道中的气体流动；将反弹和镜面反射相综合起来，准确地实现真实的气体与固体之间的相互作用。因此需要定义一个弹回系数 （ ），用来表述粒子在和壁面作用时沿原路弹回所占的比例。14.1.6 反弹与镜面反射混合格式 当 时，表示纯反弹处理当 时，表示纯镜面反射处理当 时，表示理想的漫反射 越小，壁面滑移速度越大；对任意 ，壁面法向速度为零；14.1.6 反弹与镜面反射混合格式 动力学格式主要利用边界上宏观物理量的定义，直接求解边界点上未知分布函数的方程组以获得边界节点上待定的分布函数，动力学格式主要包括。 Nobel格式 非平衡反弹格式 反滑移格式 质量修正格式14.2

7、 动力学格式 14.2.1 Nobel格式通常情况下，速度或压力其中一个已知，以速度边界条件为例：未知量：3个方程1995年，Nobel等人使用内能作为补充条件方程封闭，联立求解1995年，Inamuro等提出；边界节点上的未知分布函数可以用一个新的平衡态分布函数来替代为了克服壁面滑移，需要在这个平衡态分布函数中加入一个反滑移速度进行修正，该反滑移速度可以抵消掉采用反弹格式处理无滑移边界时所产生的滑移速度，从而使壁面上的流体速度与避免速度相等。设未知分布函数形式如下：其中 ， ， ，都是未知参数， 表示壁面的速度。14.2.2 反滑移格式通用性差1997年，Zou与He给出了另外一种补充条件：

8、假设分布函数的非平衡部分在垂直于边界的方向上仍然满足反弹格式速度边界/压力边界14.2.3 非平衡态反弹格式流场内部对于速度边界条件：已知量：未知量：四个未知量，需要四个方程The macroscopic density formula is one equation:14.2.3 非平衡态反弹格式The macroscopic velocity formula gives two equations:x-direction:y-direction:Components of ea are all unit vectorsAssuming ux = 014.2.3 非平衡态反弹格式假设分布函数

9、的非平衡部分在垂直于边界的方向上仍然满足反弹格式14.2.3 非平衡态反弹格式第四个方程，方程封闭，联立求解Two equations have the directional density unknowns f4, f7 and f8 in common, so rewrite them with those variables on the left hand side:14.2.3 非平衡态反弹格式Equating them gives: 14.2.3 非平衡态反弹格式Solving for r:Solving the bounceback equation for f4: Solvi

10、ng for f7:14.2.3 非平衡态反弹格式Solving for f8: 14.2.3 非平衡态反弹格式/ Zou and He velocity BCs on north side. for( i=0; ini; i+) fi = ftempnj-1i; rho0 = ( fi0 + fi1 + fi3 + 2.*( fi2 + fi5 + fi6) / ( 1. + uy0); ru = rho0*uy0; fi4 = fi2 - (2./3.)*ru; fi7 = fi5 - (1./6.)*ru + (1./2.)*( fi1-fi3); fi8 = fi6 - (1./6.)

11、*ru + (1./2.)*( fi3-fi1); 14.2.3 非平衡态反弹格式压力（密度）边界Dirichlet boundary conditions constrain the pressure/density at the boundariesSolution is closely related to that for velocity boundaries A density r0 is specified and velocity is computedSpecifying density is equivalent to specifying pressure since t

12、here is an equation of state (EOS) relating them directlyFor single component D2Q9 model, the relationship is simply P = RTr with RT = 1/3. More complex EOS applies to single component multiphase modelsWe assume that velocity tangent to the boundary is zero and solve for the component of velocity no

13、rmal to the boundary. 14.2.3 非平衡态反弹格式Assume that velocity tangent to the boundary is zero and solve for the component of velocity normal to the boundary Need to solve for v, f4, f7 and f8 压力（密度）边界14.2.3 非平衡态反弹格式压力（密度）边界14.2.3 非平衡态反弹格式压力（密度）边界14.2.3 非平衡态反弹格式/ Zou and He pressure boundary on north sid

14、e. for( i=0; ini; i+) fi = ftempnj-1i; uy0 = -1. + ( fi0 + fi1 + fi3 + 2.*( fi2 + fi5 + fi6) / rho0; ru = rho0*uy0; fi4 = fi2 - (2./3.)*ru; fi7 = fi5 - (1./6.)*ru + (1./2.)*( fi1-fi3); fi8 = fi6 - (1./6.)*ru + (1./2.)*( fi3-fi1); 压力（密度）边界14.2.3 非平衡态反弹格式对于充分发展流动而言，如果边界条件不能使流体在整个流动区域内满足质量守恒要求，则数值结果的精度

15、就会降低，而且数值计算的收敛速度以及稳定性都会受到影响；如果在处理充分发展边界条件时，采取一定的方法，使得流体在整个流动区域内的总体质量守恒要求得到满足，那么数值结果的精度以及计算过程的收敛状况就会有较大程度的提高；何雅玲等基于此思路，提出了格子Bolzmann方法处理充分发展流动的质量修正格式。所谓质量修正，就是通过计算进口和出口的质量流率，确定一个质量修正系数，以此系数修正出口速度，并用修正后的速度更新边界节点上的未知粒子分布函数。14.2.4 质量修正格式具体做法如下：假设出口截面上各点速度x方向分量的相对变化率为常数，即然后，考虑整体质量守恒，也即进出口的质量流率相等，得质量修正系数为

16、：利用质量修正系数即可得到出口速度分布，利用该速度计算出口边界的粒子分布函数。14.2.4 质量修正格式启发式格式和动力学格式的局限性依赖使用的格子Boltzmann模型， 通用性差；需要对边界处的一些物理性质（如速度、密度、压力等）作假设，因而很难推广到更一般的边界处理上。对一些非常规边界（包含梯度信息的边界）存在误差借鉴传统计算流体力学方法中的边界处理方法来构造格子Bolzman方法边界处理格式；外推格式非平衡外推格式14.3 外推格式1996年Chen等人提出；假定在物理边界外有一层虚拟边界，物理边界和虚拟边界格点上的分布函数值都按照标准的演化步骤进行演化，但虚拟边界格点上指向流场的分布

17、函数值使用外推方法确定。14.3.1 外推格式流体层物理边界虚拟边界在碰撞过程中在物理边界上使用给定的边界条件（速度或压力）来计算平衡态分布函数。二阶精度优点普适性好，可以根据这种方法容易地设计出一般边界条件的格式；计算量也比较小，容易实现。在确定未知分布函数时，不需要额外的假设条件，这是外推格式的技巧所在；具有二阶精度；缺点 数值稳定性较差（在外推格式中，因为有负的系数出现，所以数值误差有可能随时间增加，从而导致数值不稳定的发生）14.3.1 外推格式受外推方法和非平衡态反弹方法的启发，吸收两种方法的优点，2002年Guo等人提出了一种新的边界处理方法，即非平衡态外推方法。其基本思想是，将边

18、界节点上的分布函数分解为平衡态和非平衡态两部分，并根据具体的边界条件定义新的平衡态分布来近似平衡态部分，而非平衡态部分则使用非平衡态外推（一阶精度）确定，所得的分布函数的整体逼近精度为二阶精度。14.3.2 非平衡外推格式流场内部节点： G 边界节点： O14.3.2 非平衡外推格式流体内部物理边界流场外部对于非平衡态部分，用邻近的流体格子节点G出的非平衡态部分近似得到，即14.3.2 非平衡外推格式对于速度边界：u已知，密度未知空间精度是二阶的，时间精度也是二阶的，这与格子Boltzmann方法的整体精度是一致的。由于低阶外推格式的数值稳定性一般要优于高阶外推格式的数值稳定性，所以非平衡态外

19、推格式的数值稳定性要好于根据线性外推得到的分布函数外推格式的数值稳定性。同时，与外推格式一样，非平衡态外推格式也具有适用范围广、计算简单，容易实现的优点。14.3.2 非平衡外推格式当所要处理的物理区域并不具有规则的几何形状时适体网格非结构化网格在结构化的直角正交网格，并在适当的位置采用阶梯逼近或者差值处理，以保证满足物理边界上的条件。 阶梯逼近可与诸如反弹格式等前面提到的多种边界处理格式配合使用，实施简单，但是计算精度较低，常用的处理复杂边界处理格式为Filippova与Hanel格式、Bouzidi格式、Lallemand与Luo格式、Guo格式等14.4 复杂边界处理格式1998年， F

20、ilippova与Hanel提出构建虚拟平衡态分布函数，并用其进行线性插值的方法来计算由壁面弹回流体的分布函数值。14.4.1 Filippova与Hanel格式定义变量q，用来表示流体节点中离壁面最近的点到实际物理边界的相对距离，即式中， 为实际物理边界所在的位置。考虑到算法的稳定性， Filippova与Hanel提出用以下关系式确定 和 当 时，当 时，14.4.1 Filippova与Hanel格式2001年，Bouzidi等提出了一种结合反弹格式和空间插值的边界处理方法，用以处理任意曲率的复杂几何边界。该格式的基本思想是利用在弹回方向上用空间插值来获得由边界进入流体的分布函数值。考虑

21、如下图一般的二维情况 ，则线性插值格式给定如下：当 时当 时14.4.2 Bouzidi格式14.4.2 Bouzidi格式在Bouzidi格式中，差值公式包含了在不同时间层次上的分布函数。2003年, Lallemand与Luo针对运动边界对Bouzidi格式进行了改进，在引用差值公式时做了变化，把所有分布函数统一到同一个时间层次上，同时考虑了碰撞和迁移过程，提出了Lallemand与Luo格式。2002年，Guo等提出了一种结合非平衡外推格式和空间差值曲线边界处理方法，其基本思想是，将边界节点的分布函数分为平衡态和非平衡态两部分，平衡态部分由一个虚拟的平衡态分布函数代替，非平衡态部分则由相

22、邻流体节点的非平衡部分插值获得，提出了Guo格式。14.4.3 Lallemand与Luo格式、Guo格式14.5 LBM模拟的基本步骤Define velocity vectors012345678D2Q9 e1e2e3e4e5e6e7e8%define esex(0)= 0; ey(0)= 0ex(1)= 1; ey(1)= 0ex(2)= 0; ey(2)= 1ex(3)=-1; ey(3)= 0ex(4)= 0; ey(4)=-1ex(5)= 1; ey(5)= 1ex(6)=-1; ey(6)= 1ex(7)=-1; ey(7)=-1ex(8)= 1; ey(8)=-114.5 LB

23、M模拟的基本步骤Problem definitionLY=10LX=20tau = 1( relaxation time)g=0.00001(gravity acceleration)%set solid nodesis_solid_node=zeros(LY,LX)for i=1:LX is_solid_node(1,i)=1 is_solid_node(LY,i)=1endLXLY% if is_interior_solid_node(j,i)14.5 LBM模拟的基本步骤Initialize density and fs (assuming zero velocity)%define i

24、nitial density and fsrho=ones(LY,LX);f(:,:,1) = (4./9. ).*rho; f(:,:,2) = (1./9. ).*rho;f(:,:,3) = (1./9. ).*rho;f(:,:,4) = (1./9. ).*rho;f(:,:,5) = (1./9. ).*rho;f(:,:,6) = (1./36.).*rho;f(:,:,7) = (1./36.).*rho;f(:,:,8) = (1./36.).*rho;f(:,:,9) = (1./36.).*rho;14.5 LBM模拟的基本步骤/ Computing macroscopi

25、c density, rho, and velocity, u=(ux,uy). for( j=0; jLY; j+) for( i=0; iLX; i+) u_xji = 0.0; u_yji = 0.0; rhoji = 0.0; if( !is_solid_nodeji) for( a=0; a9; a+) rhoji += fjia; u_xji += exa*fjia; u_yji += eya*fjia; u_xji /= rhoji; u_yji /= rhoji; 14.5 LBM模拟的基本步骤14.5 LBM模拟的基本步骤Periodic BoundariesOn bound

26、ary, neighboring point is on opposite boundary ip = ( i0 )?( i-1):( LX-1); jp = ( j0 )?( j-1):( LY-1);/ Compute the equilibrium distribution function, feq. f1=3.; f2=9./2.; f3=3./2.; for( j=0; jLY; j+) for( i=0; iLX; i+) if( !is_solid_nodeji) rt0 = (4./9. )*rhoij; rt1 = (1./9. )*rhoij; rt2 = (1./36.

27、)*rhoij; ueqxij = uxij; /Can add forcing here; see MATLAB code ueqyij = uyij; uxsq = ueqxij * ueqxij; uysq = ueqyij * ueqyij; uxuy5 = ueqxij + ueqyij; uxuy6 = -ueqxij + ueqyij; uxuy7 = -ueqxij + -ueqyij; uxuy8 = ueqxij + -ueqyij; usq = uxsq + uysq; 14.5 LBM模拟的基本步骤 feqij0 = rt0*( 1. - f3*usq); feqij1

28、 = rt1*( 1. + f1*ueqxij + f2*uxsq - f3*usq); feqij2 = rt1*( 1. + f1*ueqyij + f2*uysq - f3*usq); feqij3 = rt1*( 1. - f1*ueqxij + f2*uxsq - f3*usq); feqij4 = rt1*( 1. - f1*ueqyij + f2*uysq - f3*usq); feqij5 = rt2*( 1. + f1*uxuy5 + f2*uxuy5*uxuy5 - f3*usq); feqij6 = rt2*( 1. + f1*uxuy6 + f2*uxuy6*uxuy6 - f3*usq); feqij7 = rt2*( 1. + f1*uxuy7 + f2*uxuy7*uxuy7 - f3*usq); feqij8 = rt2*( 1. + f1*uxuy8 + f2*uxuy8*uxuy8 - f3*usq); 14.5 LBM