用Lattice Boltzmann方法模拟二维水力过渡过程

《用Lattice Boltzmann方法模拟二维水力过渡过程》由会员分享，可在线阅读，更多相关《用Lattice Boltzmann方法模拟二维水力过渡过程（10页珍藏版）》请在装配图网上搜索。

1、用Lattice Boltzmann方法模拟二维水力过渡过程程永光，张 慧(武汉大学水利电力学院)摘 要：本文旨在探索用Lattice Boltzmann方法模拟二维水力过渡过程的可行性。应用Chapman Enskog展开技术，导出了一套多尺度方程，将LB方程和其所对应的宏观方程联系起来；根据水力过渡过程基本方程的特点，建立了计算二维水力过渡过程的一个LB模型；应用该模型对一维正弦波的传播、一维断波的传播与反射、二维拐角衍射波等典型算例，以及近似为二维问题的混凝土蜗壳在机组甩负荷和起动工况下的水力过渡过程进行了模拟分析。研究表明用LB方法模拟多维水力过渡过程是可行的，值得进一步探讨。关键词：

2、水力过渡过程；Lattice Boltzmann方法；模拟；二维流场作者简介：程永光(1968-)，男，山西武乡人，博士，副教授，主要从事水电站水力学研究。在水利水电工程中，当流道横向尺度与纵向尺度相当时，水流的二维或三维特性往往比较明显，用一维方法进行水力过渡过程计算一般难以反映实际，多维计算分析势在必行。二十世纪七十年代以来，国内外先后出现过格子法(Latticework)2，5、双特征线法(Bicharacteristics)4、近特征线法(Nearcharacteristics)3、类特征线法(Characteristics like)1，5等多种计算多维流体过渡过程的数值方法，它们计

3、算过一些二维问题，取得了较好结果。这些方法有容易失稳3、格式复杂4和计算网格固定1，5等不足，计算二维过渡过程一般均可，但如果计算三维问题，难度会大大增加。Lattice Boltzmann(LB)方法是近十年来发展起来的一种由运动论原理出发的流动计算技术68，具有算法简单、压力直接算出、并行度高、几何边界易处理等优点，目前在非线性偏微分方程求解、多相流、多孔介质流、特别是大规模三维流场的模拟上显示了较大潜力。能不能发挥LB方法的优点并将之引入水电工程流场计算，值得逐步探讨。文献7建立了Lattice Boltzmann一维水击计算模型，并通过典型算例的计算分析与特征线法进行了比较，结果表明L

4、B方法计算一维水力过渡过程有效可行。本文建立LB二维水力过渡过程模型，计算一些理想流动并尝试模拟混凝土蜗壳内的水力过渡过程，为将来把LB方法应用于实际二维和三维水力过渡过程的模拟打下基础和提供借鉴。1 LB方程及多尺度方程1.1 Lattice Boltzmann方程LB方法的核心是Lattice Boltzmann方程，这是一组可由统计力学中Boltzmann方程导出的，关于固定网格子粒子分布函数的演化式。它的最简单形式是带BGK单松弛碰撞项的LB方程6f(r+e,t+1)-f(r,t)=-1/f(r,t)-f(0)(r,t),=0，1，2，b(1)其中f(r,t)是单粒子分布函数；f(0)

5、(r,t)为局部平衡分布函数，是局部达到平衡态时的粒子分布函数；为无量纲松弛时间，反映f(r,t)趋近于f(0)(r,t)的速度；r 是计算网格节点；e是网格节点连线形成的向量，也就是粒子的运动速度向量；表示给定的粒子运动方向；粒子共有b+1个运动方向(包括e0=0).方程的右边反映运动粒子相互碰撞的影响，左边反映粒子的迁移。LB方法模拟流场的过程实际上就是中观方程(1)在计算网格上的简单迭代，流场的宏观量如密度、流速等按下式由中观量f(r,t)计算得到(2)1.2 Chapman Enskog展开及多尺度方程LB方程的演化对应一定的宏观流动现象。利用Chapman Enskog展开和多尺度分

6、析技术可以给出LB方程和它所对应宏观方程的联系8.设流场离平衡态不远，分布函数f(r,t)对局部平衡分布函数f(0)(r,t)只有微小的偏离，且满足(3)于是可作如下Chapman Enskog展开f=f(0)+f(1)+2f(2)+O(3)(4)其中1是个很小的正数，称为Knudsen数，是粒子运动的平均自由程与宏观特征长度的比值，-1与网格划分的规模相当。这样就有(5)再引进两种时间尺度t1,t2和空间尺度r1，并令t=t1+2t2,r=r1(6)对f(r+e,t+1)在(r,t)作Taylor展开，利用(2)(6)式，经一系列数学变换8，得到如下多尺度方程(7)(8)O(2)：(9)(1

7、0)式中：j,k,m是正交坐标方向，以上采用了张量标记法。2 二维水力过渡过程LB模型的建立2.1 二维水力过渡过程基本方程二维水力过渡过程的基本方程可由典型的Euler方程得到(11)(12)dp=c2d(13)式(11)是连续方程，式(12)是动量方程，式(13)是等熵状态方程，其中，为密度，p是压力，ui为速度分量(i,j=1,2),c是流体波速。对于uc，即Mach数较小的流动，式(12)中的对流项常被略去，这样一来式(11)(13)就变成线性化的水力过渡过程基本方程。2.2 平衡分布函数及对应的宏观方程LB模型的建立，实际上就是确定式(1)中的平衡分布函f(0)。f(0)的确定原则是

8、使它对应的宏观方程与要模拟流场的控制方程相容。对于二维问题，常用的有九点正方形网格模型，它的九个粒子运动向量e=(ex，ey)(=0,1,8)(见图1)，表示为矩阵(14)根据LB方法平衡分布函数的形式，考虑到二维水力过渡过程的基本方程(11)(13)，适当选取局部平衡分布函数f(0)的表达式f(0)=t1+3(eu)+9/2(eu)2-3/2|u|2，t0=4/9;t=1/9,=1,2,3,4;t=1/36,=5,6,7,8 (15)图1 正方形九点模型示意这样的平衡分布函数对应于什么样的宏观方程?将式(14)、(15)代入多尺度方程，进行尺度“粘合”，即(7)+(9)2，(8)+(10)2

9、，并注意到尺度展开式(4)、(6)，则可得到如下宏观方程(16)(17)dp=C2sd(18)其中波速Cs=1/;=2-1/6为运动粘性系数；Sij=1/2(iuj+jui)；M为Mach数。方程(16)(18)是平衡分布函数取式(15)时LB模型对应的宏观方程，可以看出，它们与方程(11)(13)十分一致。仅有的差别是：式(16)和式(17)多了O(2)和O(M3)，因为与计算网格尺度相当，故Mach数较小时，以上LB模型具有二阶空间精度；j(2Sij)对应于NS方程中的粘性项，是实际流体中存在的，计算中可适当调整粘性系数以使计算稳定而又不至于过分耗散；式(18)中的模型波速Cs=1/，对应

10、的式(13)中的c是实际流体的波速，它们一般不相等，因此用LB方法模拟方程(11)(13)需要根据相似原理进行尺度变换8。若要模拟线性化的水力过渡过程基本方程，只须将平衡分布函数表达式(15)中最后两个平方项去掉即可。2.3 模拟方法及边界处理用LB方法模拟流场时，首先要根据模型、流场形状及计算要求将流场划分为计算网络。以往LB只能在规则均匀的网格上运行，近来发展出一种基于插值的LB非均匀网格算法。LB方法的具体计算过程十分简单。在每一时步内的基本计算步骤是：t时刻，将平衡分布函数(15)代入LB方程(1)，将方程(1)在所有的计算网格r上求解，得到f(r+e,t+1)；如果是规则均匀的网格，

11、则r+e也是网格点，若是非均匀网格，则r+e不在网格点上，我们需要由f(r+e,t+1)插值得到网格点上的分布函数f(r,t+1);最后由式(2)计算出各个网格节点上的密度和速度等宏观量。如此反复，所得结果就是宏观流动，即方程(16)(18)的解。计算中，每个时刻所有节点上各个方向的f均需求得。在边界上，总有一些方向的f不能直接按LB方程计算，需根据边界情况确定。方程(2)可以提供3个边界条件方程，计算证明，只要这3个方程得到满足，就能保证边界处理有二阶精度。对于速度边界，压力是未知的，对于压力边界，某速度分量是未知的，这些未知的压力或速度分量与未知的分布函数一道，构成的末知量个数一般大于3.

12、为了将边界方程封闭，需适当增补一些条件。如图1所示的水平线边界，有压力p或切向流速u，以及分布函数f2，f5，f6共四个未知量，常令f2=f4为增补的条件以使方程封闭。3 算例分析3.1 一维波动传播与反射这里用二维模型对二维的极限情况 一维波动的传播和反射进行模拟。第一个算例1，5的计算网格布置和边界条件如图2，上游为速度边界，下游为自由边界，两侧是循环边界。初始流场静止，计算开始后给定上游流速u=0.1sin(t/20)。此例用模拟线性化水力过渡过程的LB模型计算。图2展示了 t=200时刻(时步)流速沿x轴线的分布曲线。可见，除波前略有耗散外，计算值与解析值非常吻合，说明本文LB模型基本

13、上没有频散、耗散也很小，波速也是准确的。第二个算例的网格布置同算例一，但上游为压力边界，下游为速度边界，两侧仍为循环边界。初始时刻整个流场p=1.0，u=0.1.在t=0时刻下游边界流速突然降为零，之后继续为零，而上游边界压力一直保持p=1.0(即=3)不变，这样就得出图3所示的压力波(断波)的传播和反射的典型水击现象。图中实线为压力沿程分布线，虚线为流速沿程分布线，点划线表示下游边界的压力随时间的变化过程。压力波的传播和反射过程完全符合理论分析。图上压力波传播的速度接近理论值Cs=1/0.5773.流速突然变为零所引起的下游边界压力变化值约为0.173，也与理论值p=-uCs=0.13/=0

14、.1732很吻合。传播过程中波锋面虽有坦化趋势，但仍保持了较陡峻的形状，宽度在十个网格左右。该例说明本文的LB模型能够正确反映压力波的传播和反射特性，并较好地模拟断波。图2 压力分布曲线图3 断波的传播与反射3.2 二维波的拐角衍射文献15，9都对一断波迎面撞在一90拐角障碍物上而形成的衍射波作为模拟分析。这里，给定初始静止流场的压力p1=1.0，在距上游边界25的位置布置一个矩形障碍物。t=0时刻后，强加上游边界压力p2=1.1，该压力向下游传播。计算到t=63时步时，得到图4所示的压力分布图。它与文献9的解析结果图案十分相近。图5是拐角两条边上压力分布曲线的比较情况，可见LB方法的结果与其

15、他方法计算的结果和解析解吻合较好。图4 拐角衍射压力分布图5 拐角两条边上的压力分布比较3.3 当作二维问题的混凝土蜗壳内水力过渡过程严格说，用二维方法模拟混凝土蜗壳内的水力过渡过程是不合适的，因为二维模型难以反映蜗壳高度变化的影响。这里把它作二维问题处理，目的是验证LB方法的模拟能力，而非真正计算混凝土蜗壳的水力过渡过程，结果不能作实际应用。该计算取图6所示尺寸的虚拟混凝土蜗壳，包角250。蜗壳进口平均流速Vc=5m/s，引水道及蜗壳内波速400m/s，额定水头Hr=40m.设引水道进口断面处是水位恒定的水库，导叶为孔口出流条件，绝对速度方向角45，其它边界为非滑移固壁边界，导叶直线启闭，启

16、闭时间为Ts=1s.图7是正常运行工况的流速分布图，由于这里的LB模型没有加入湍流模式，计算中边壁取非滑移条件，故显示层流的流速分布规律，同时粘性偏大导致了水力坡降偏大，虽然与实际有差别，但对以压力计算为主的水力过渡过程不会有太明显影响。图8是甩负荷工况下蜗壳内几个特征点(见图6)的压力变化过程，变化规律呈现一维水击计算中的末相水击特性，而且越向蜗壳尾部压力越高，与实际相符。图9展示的是甩负荷工况导叶关闭终了时刻t=Ts=1压力在蜗壳内的分布情况。图10为机组启动工况下特征点的压力变化过程，越向蜗壳尾部压力下降越大。由于稍大，压力波动的衰减稍微快一些。但总的来说能较好反映蜗壳内水力过渡过程的二

17、维特性。图6 混凝土蜗壳的布置和尺寸图7 正常运行工况下的流速向量图8 甩负荷工况不同点压力变化过程图9 甩负荷 工况导叶关闭终了时刻压力分布4 结 语为了尝试把新近发展起来的LB方法应用于二维水力过渡过程计算，本文由LB方法的主要原理和二维水力过渡过程基本方程出发，推导建立了LB的二维水力过渡过程计算模型。通过对一维正弦波的传播、一维断波的传播与反射、二维拐角绕射波等典型算例，以及近似为二维问题的混凝土蜗壳在机组甩负荷和起动工况下的水力过渡过程的初步模拟分，证明LB方法计算二维水力过渡过程是可行的。当然，要将LB真正应用于实际多维水力过渡过程计算，很多问题还有待研究，如LB模型与湍流模式的结

18、合、LB模型计算时的粘性系数的合理取值、滑移边界的恰当处理，能反映第三维影响(如混凝土蜗壳高度变化)的二维模型，以及三维水力过渡过程模型等。图10 机组启动工况不同点压力变化过程参 考 文 献：1 Wylie E B.Linearized two dimensional fluid transientsJ。Journal of Fluids Engineering,ASME,1984,106:227-2322 Wylie E B,Streeter V L.Multidimensional fluid transients by latticeworkJ。Journal of Fluids En

19、gineering,ASME,1980,102:203-2103 Shin Y W,Kot C A K.Two dimensional fluid transient analysis by the method of nearcharacteristicsJ。Journal of Computational Physics,1978,28:211-2314 Shin Y W,Valentin R A.Numerical analysis of fluid hammer waves by the method of characteristicsJ。Journal of Computation

20、al Physics,1976,20:220 2375 Wylie F B,Streeter V L,Suo L S.Fliud transients in systemsM。Prentice Hall,Inc.New Jersey.U.S.A.,19936 Qian Y,dHumieres D,Lallemand P.Lattice BGK models for Navier Stokes equationJ。Europhysics Letters.1992,17:479-4847 程永光，张师华，陈鉴治。用Lattice Boltzmann方法模拟水击J。水利学报，1998(6)：25-308 程永光。Lattice Boltzmann方法及其在流场分析中应用的研究：D。武汉：武汉水利电力大学，1998.69 Keller J B.Diffraction of a shock or an electromagnetic pulse by right angled wedgeJ。Journal of Applied Physics,1952,23:1267 1268