当前位置:首页 >> 电力/水利 >>

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


用 Lattice Boltzmann 方法模拟二维水力过渡过程 程永光,张 慧 (武汉大学水利电力学院)
摘 要: 本文旨在探索用 Lattice Boltzmann 方法模拟二维水力过渡过程的可行性。 应用 Chapman Enskog 展开技术,导出了一套多尺度方程,将 LB 方程和其所对应的宏观方程联系起来;根据水 力过渡过程基本方程的特点,建立了计算二维水力过

渡过程的一个 LB 模型;应用该模型对一维 正弦波的传播、一维断波的传播与反射、二维拐角衍射波等典型算例,以及近似为二维问题的混 凝土蜗壳在机组甩负荷和起动工况下的水力过渡过程进行了模拟分析。研究表明用 LB 方法模拟 多维水力过渡过程是可行的,值得进一步探讨。 关键词:水力过渡过程;Lattice Boltzmann 方法;模拟;二维流场 作者简介:程永光(1968-),男,山西武乡人,博士,副教授,主要从事水电站水力学研究。

在水利水电工程中,当流道横向尺度与纵向尺度相当时,水流的二 维或三维特性往往比较明显, 用一维方法进行水力过渡过程计算一般难 以反映实际,多维计算分析势在必行。二十世纪七十年代以来,国内外 先后出现过格子法(Latticework)
[4] [2,5]

、双特征线法
[3]

(Bicharacteristics) 、近特征线法(Nearcharacteristics) 、类特 征线法(Characteristics like)
[1,5]

等多种计算多维流体过渡过程的数

值方法,它们计算过一些二维问题,取得了较好结果。这些方法有容易 失稳 、 格式复杂 和计算网格固定
[3] [4] [1,5]

等不足, 计算二维过渡过程一般

均可,但如果计算三维问题,难度会大大增加。Lattice Boltzmann(LB) 方法是近十年来发展起来的一种由运动论原理出发的流动计算技术
8] [6~

,具有算法简单、压力直接算出、并行度高、几何边界易处理等优点,

目前在非线性偏微分方程求解、多相流、多孔介质流、特别是大规模三 维流场的模拟上显示了较大潜力。能不能发挥 LB 方法的优点并将之引

入水电工程流场计算,值得逐步探讨。文献[7]建立了 Lattice Boltzmann 一维水击计算模型,并通过典型算例的计算分析与特征线法 进行了比较,结果表明 LB 方法计算一维水力过渡过程有效可行。本文 建立 LB 二维水力过渡过程模型,计算一些理想流动并尝试模拟混凝土 蜗壳内的水力过渡过程,为将来把 LB 方法应用于实际二维和三维水力 过渡过程的模拟打下基础和提供借鉴。 1 LB 方程及多尺度方程 1.1 Lattice Boltzmann 方程 LB 方法的核心是 Lattice Boltzmann 方 程,这是一组可由统计力学中 Boltzmann 方程导出的,关于固定网格子 粒子分布函数的演化式。它的最简单形式是带 BGK 单松弛碰撞项的 LB 方程
[6]

fα (r+eα ,t+1)-fα (r,t)=-1/τ [fα (r,t)-fα (0)(r,t)],α =0, 1,2,?b
(0)

(1)

其中 fα (r,t)是单粒子分布函数;fα (r,t)为局部平衡分布函数,是局 部达到平衡态时的粒子分布函数;τ 为无量纲松弛时间,反映 fα (r,t) 趋近于 fα (r,t)的速度;r 是计算网格节点;eα 是网格节点连线形成 的向量,也就是粒子的运动速度向量;α 表示给定的粒子运动方向;粒 子共有 b+1 个运动方向(包括 e0=0).方程的右边反映运动粒子相互碰撞 的影响,左边反映粒子的迁移。
(0)

LB 方法模拟流场的过程实际上就是中观方程(1)在计算网格上的简 单迭代,流场的宏观量如密度、流速等按下式由中观量 fα (r,t)计算得 到
(2)

1.2 Chapman Enskog 展开及多尺度方程 LB 方程的演化对应一定的宏观 流动现象。利用 Chapman Enskog 展开和多尺度分析技术可以给出 LB 方 程和它所对应宏观方程的联系 . 设流场离平衡态不远,分布函数 fα (r,t)对局部平衡分布函数 fα
(0) [8]

(r,t)只有微小的偏离,且满足
(3)

于是可作如下 Chapman Enskog 展开
fα =fα (0)+ε fα (1)+ε 2fα (2)+O(ε 3) (4)

其中ε 《1 是个很小的正数,称为 Knudsen 数,是粒子运动的平均自由 程与宏观特征长度的比值,ε
-1

与网格划分的规模相当。这样就有
(5)

再引进两种时间尺度 t1,t2 和空间尺度 r1,并令 =ε t1+ε
2 t2

t

,

r



r1

(6)

对 fα (r+eα ,t+1)在(r,t)作 Taylor 展开,利用(2)~(6)式,经一系列 数学变换 ,得到如下多尺度方程
(7) (8)
[8]

O(ε ):

2

(9)

(10)

式中:j,k,m 是正交坐标方向,以上采用了张量标记法。 2 二维水力过渡过程 LB 模型的建立 2.1 二维水力过渡过程基本方程 二维水力过渡过程的基本方程可由典 型的 Euler 方程得到
(11) (12) dp=c2dρ (13)

式(11)是连续方程, 式(12)是动量方程, 式(13)是等熵状态方程, 其中, ρ 为密度,p 是压力,ui 为速度分量(i,j=1,2),c 是流体波速。 对于 u《c,即 Mach 数较小的流动,式(12)中的对流项常被略去, 这样一来式(11)~(13)就变成线性化的水力过渡过程基本方程。

2.2 平衡分布函数及对应的宏观方程 LB 模型的建立,实际上就是确定 式(1)中的平衡分布函 fα 。fα 的确定原则是使它对应的宏观方程与 要模拟流场的控制方程相容。 对于二维问题,常用的有九点正方形网格模型,它的九个粒子运动向量 eα =(eα x,eα y)(α =0,1,?8)(见图 1),表示为矩阵
(0) (0)

(14) 根据 LB 方法平衡分布函数的形式, 考虑到二维水力过 渡过程的基本方程(11)~(13),适当选取局部平衡分布函 数 fα (0)的表达式 fα (0)=tα ρ [1+3(eα ·u)+9/2(eα ·u)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)×ε , (8)×ε +(10)×ε , 并注意到尺度展开式(4)、 (6), 则可得到如下宏观方程
(16) (17) dp=C2sdρ 其中波速 Cs=1/ ;ν =2τ -1/6 为运动粘性系数;Sij=1/2(
i

2

2

(18) uj+
j

ui);M 为 Mach

数。 方程(16)~(18)是平衡分布函数取式(15)时 LB 模型对应的宏观方程,可以 看出,它们与方程(11)~(13)十分一致。仅有的差别是:①式(16)和式(17)多了 O(ε 2)和 O(ε M3),因为 ε 与计算网格尺度相当,故 Mach 数较小时,以上 LB 模

型具有二阶空间精度; ②

j

(2ν ρ Sij)对应于 NS 方程中的粘性项, 是实际流体中

存在的,计算中可适当调整粘性系数 ν 以使计算稳定而又不至于过分耗散;③ 式(18)中的模型波速 Cs=1/ , 对应的式(13)中的 c 是实际流体的波速,它们一

般不相等,因此用 LB 方法模拟方程(11)~(13)需要根据相似原理进行尺度变换 [8] 。

若要模拟线性化的水力过渡过程基本方程, 只须将平衡分布函数表 达式(15)中最后两个平方项去掉即可。 2.3 模拟方法及边界处理 用 LB 方法模拟流场时,首先要根据模型、流 场形状及计算要求将流场划分为计算网络。以往 LB 只能在规则均匀的 网格上运行,近来发展出一种基于插值的 LB 非均匀网格算法。 LB 方法的具体计算过程十分简单。在每一时步内的基本计算步骤 是:①t 时刻,将平衡分布函数(15)代入 LB 方程(1),将方程(1)在所有 的计算网格 r 上求解,得到 fα (r+eα ,t+1);②如果是规则均匀的网格, 则 r+eα 也是网格点,若是非均匀网格,则 r+eα 不在网格点上,我们需 要由 fα (r+eα ,t+1)插值得到网格点上的分布函数 fα (r,t+1);③最后由 式(2)计算出各个网格节点上的密度和速度等宏观量。如此反复,所得 结果就是宏观流动,即方程(16)~(18)的解。 计算中,每个时刻所有节点上各个方向的 fα 均需求得。在边界上, 总有一些方向的 fα 不能直接按 LB 方程计算,需根据边界情况确定。 方程(2)可以提供 3 个边界条件方程,计算证明,只要这 3 个方程 得到满足,就能保证边界处理有二阶精度。对于速度边界,压力是未知

的,对于压力边界,某速度分量是未知的,这些未知的压力或速度分量 与未知的分布函数一道, 构成的末知量个数一般大于 3.为了将边界方程 封闭,需适当增补一些条件。如图 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 模型基本上没有频散、耗散也很小,波速也是 准确的。 第二个算例的网格布置同算例一,但上游为压力边界,下游为速度边界,两 侧仍为循环边界。初始时刻整个流场 p=1.0,u=0.1.在 t=0 时刻下游边界流速突 然降为零,之后继续为零,而上游边界压力一直保持 p=1.0(即 ρ =3)不变,这样 就得出图 3 所示的压力波(断波)的传播和反射的典型水击现象。 图中实线为压力 沿程分布线, 虚线为流速沿程分布线,点划线表示下游边界的压力随时间的变化 过程。 压力波的传播和反射过程完全符合理论分析。图上压力波传播的速度接近 理论值 Cs=1/ ≈0.5773.流速突然变为零所引起的下游边界压力变化值约为 =0.1732 很吻合。传播过程中波锋

0.173,也与理论值 Δ p=-Δ uρ Cs=0.1×3/

面虽有坦化趋势,但仍保持了较陡峻的形状,宽度在十个网格左右。该例说明本 文的 LB 模型能够正确反映压力波的传播和反射特性,并较好地模拟断波。

图 2 压力分布曲线

图 3 断波的传播与反射

3.2 二维波的拐角衍射 文献[1~5,9]都对一断波迎面撞在一 90°拐角障碍物 上而形成的衍射波作为模拟分析。这里,给定初始静止流场的压力 p1=1.0,在距 上游边界 25 的位置布置一个矩形障碍物。 t=0 时刻后, 强加上游边界压力 p2=1.1, 该压力向下游传播。计算到 t=63 时步时,得到图 4 所示的压力分布图。它与文 献[9]的解析结果图案十分相近。 图 5 是拐角两条边上压力分布曲线的比较情况, 可见 LB 方法的结果与其他方法计算的结果和解析解吻合较好。

图 4 拐角衍射压力分布

图 5 拐角两条边上的压力分布比较

3.3 当作二维问题的混凝土蜗壳内水力过渡过程 严格说,用二维方法模拟混凝 土蜗壳内的水力过渡过程是不合适的, 因为二维模型难以反映蜗壳高度变化的影 响。这里把它作二维问题处理,目的是验证 LB 方法的模拟能力,而非真正计算 混凝土蜗壳的水力过渡过程, 结果不能作实际应用。该计算取图 6 所示尺寸的虚 拟混凝土蜗壳,包角 250°。蜗壳进口平均流速 Vc=5m/s,引水道及蜗壳内波速 400m/s,额定水头 Hr=40m.设引水道进口断面处是水位恒定的水库,导叶为孔口 出流条件,绝对速度方向角 45°,其它边界为非滑移固壁边界,导叶直线启闭, 启闭时间为 Ts=1s.图 7 是正常运行工况的流速分布图,由于这里的 LB 模型没有 加入湍流模式,计算中边壁取非滑移条件,故显示层流的流速分布规律,同时粘 性偏大导致了水力坡降偏大, 虽然与实际有差别,但对以压力计算为主的水力过 渡过程不会有太明显影响。图 8 是甩负荷工况下蜗壳内几个特征点(见图 6)的压 力变化过程, 变化规律呈现一维水击计算中的末相水击特性,而且越向蜗壳尾部 压力越高,与实际相符。图 9 展示的是甩负荷工况导叶关闭终了时刻 t=Ts=1 压 力在蜗壳内的分布情况。图 10 为机组启动工况下特征点的压力变化过程,越向 蜗壳尾部压力下降越大。由于稍大,压力波动的衰减稍微快一些。但总的来说能 较好反映蜗壳内水力过渡过程的二维特性。

图 6 混凝土蜗壳的布置和尺寸

图 7 正常运行工况下的流速向量

图 8 甩负荷工况不同点压力变化过程

图 9 甩负荷 工况导叶关闭终了时刻压 力分布

4 结 语
为了尝试把新近发展起来的 LB 方法应用于二 维水力过渡过程计算,本文由 LB 方法的主要原理 和二维水力过渡过程基本方程出发,推导建立了 LB 的二维水力过渡过程计算模型。通过对一维正 弦波的传播、一维断波的传播与反射、二维拐角 绕射波等典型算例,以及近似为二维问题的混凝 土蜗壳在机组甩负荷和起动工况下的水力过渡过 程的初步模拟分,证明 LB 方法计算二维水力过渡 过程是可行的。 当然,要将 LB 真正应用于实际多维水力过渡 过程计算,很多问题还有待研究,如 LB 模型与湍图 10 机组启动工况不同点压力变化 流模式的结合、 LB 模型计算时的粘性系数 ν 的合 过程 理取值、滑移边界的恰当处理,能反映第三维影 响(如混凝土蜗壳高度变化)的二维模型,以及三 维水力过渡过程模型等。 参 考 文 献:
[1] Wylie E B.Linearized two dimensional fluid transients[J] 。Journal of Fluids Engineering,ASME,1984,106:227-232 [2] Wylie E B,Streeter V L.Multidimensional fluid transients by latticework[J] 。Journal of Fluids Engineering,ASME,1980,102:203-210 [3] Shin Y W,Kot C A K.Two dimensional fluid transient analysis by the method of nearcharacteristics[J]。Journal of Computational Physics,1978,28:211-231[4] Shin Y W,Valentin R A.Numerical analysis of fluid hammer waves by the method of characteristics[J]。Journal of Computational Physics,1976,20:220 237

[5] Wylie F B,Streeter V L,Suo L S.Fliud transients in systems[M]。Prentice Hall,Inc.New Jersey.U.S.A.,1993 [6] Qian Y,d'Humieres D,Lallemand P.Lattice BGK models for Navier Stokes equation[J] 。 Europhysics Letters.1992,17:479-484 [7] 程永光,张师华,陈鉴治。用 Lattice Boltzmann 方法模拟水击[J]。水利学报,1998(6): 25-30 [8] 程永光。Lattice Boltzmann 方法及其在流场分析中应用的研究:[D]。武汉:武汉水利电 力大学,1998.6 [9] Keller J B.Diffraction of a shock or an electromagnetic pulse by right angled wedge[J]。Journal of Applied Physics,1952,23:1267 1268


相关文章:
s格子Boltzmann方法及其应用 安明论文
计算流程图 用格子 Boltzamnn 方法模拟二维间题的...二维水力过渡过程、溃坝波以及明渠非恒定流的模拟 ...Thorne,Jr. Lattice Boltzmann Modeling[M]. Florida...
更多相关标签:
lattice boltzmann | 水力过渡过程 | 水电站水力过渡过程 | 水力机组过渡过程 | 二维过渡金属硫化物 | 二维过渡金属碳化物 | boltzmann分布 | boltzmann |