基于格子Boltzmann方法和有限体积法的方柱绕流特性对比分析

分析

史冬岩;李红群;王志凯

【摘 要】基于黏性流体理论,分别采用格子Boltzmann方法(LBM)和有限体积法(FVM)建立了黏性流场中方柱绕流模型.探究LBM在非光滑曲面钝体绕流方面的应用;并结合FVM进行对比分析.在FVM模型中,采用局部加密的方法对钝体边界进行处理,而在LBM模型中,除了传统的Half-way边界处理方法,还结合了拐角边界处理方法.为获得较好的可对比数据,根据已发表文献中的理论及UDF编译码技术,分别对两模型的进出口边界条件进行了讨论和设置.对比分析了两模型下的速度云图以及获得的升、阻力系数,Strouhal数.结果发现方柱上游压力不受涡脱落影响,雷诺数对其影响也较小;两种方法下的速度、无量纲参数吻合较好;但两者最适进出口边界不同,且相同条件下,LBM比FVM数值模拟能更快达到稳定状态. 【期刊名称】《科学技术与工程》 【年(卷),期】2015(015)028 【总页数】7页(P96-102)

【关键词】方柱绕流;格子Boltzmann方法;有限体积法;对比分析 【作 者】史冬岩;李红群;王志凯

【作者单位】哈尔滨工程大学机电工程学院,哈尔滨150001;哈尔滨工程大学机电工程学院,哈尔滨150001;哈尔滨工程大学机电工程学院,哈尔滨150001 【正文语种】中 文

【中图分类】O357.1

方柱由于其结构简单，因此方柱绕流一直以来是钝体绕流中最具典型、最基础的范例之一，工程上如房屋建筑、海洋立管等很多的现象也可以用方柱绕流模型简化研究[1，2]。

目前分析柱体绕流的数值方法主要有基于有限体积法(FVM)、有限单元法(FEM)、有限差分法(FDM)和格子Boltzmann方法(LBM)等，国内外学者利用传统方法和新兴方法在这方面有很多的研究。Sohankar[3]探究了低雷诺数下，以不同入射角绕流的方柱涡脱落形态以及阻塞比、边界条件对流场的影响。Geller[4]用格子Boltzmann方法、有限单元法和有限体积法对层流流动进行了基本的计算，比较了LBM与FEM、FVM所得的结果，讨论了LBM算法的有效性。刘天成等[5]利用LB并行算法对方柱绕流进行了计算，探究了雷诺数对流场参数的影响。Arumuga[6，7]用格子Boltzmann方法计算了二维不可压缩流稳态和非稳态下的钝体绕流，研究了雷诺数、阻塞比和流场长度对方柱和圆柱绕流的影响。Regulski等[8]用格子Boltzamnn方法和谱元法对比研究了钝体绕流流场中涡街特性。陈静涛[9]用Fluent软件计算了二维湍流模型，通过各项流体参数分析流场中圆柱的尾涡特性。虽然学者们分析了阻塞比、立柱位置及流场长度对流场的影响，以及用CFX和LBM计算了钝体绕流。但是基于FVM讨论流场阻塞比，立柱位置和流场长度对绕流的影响的研究较少，也很少具体的对比分析LBM与FVM的结果；且在用Fluent进行计算的时候，入流边界的速度大都采用的是恒定的值。本文用格子Boltzmann方法和有限体积法分别对方柱绕流进行了数值模拟计算。根据非平衡态外推格式对方柱角点进行LB拐角边界处理，与以往Fluent设置的恒定速度不同，计算采用能降低流场影响的抛物线形式的来流速度；并基于有限体积法探究了不同边界条件对方柱绕流结果的影响，通过阻力系数的对比，可以看出传统方法数值模拟要求的流场边界与LBM计算要求的流场边界条件有一定的区别。

最后，将两种方法所得的结果进行了对比分析。 1.1 控制模型

有限体积法是基于N-S方程，计算的是流体的宏观行为，其控制方程[10]为： 而基于分子动理论的格子Boltzmann方法属于介观方法，考虑的是大量微观粒子的平均行为，计算时采用具有单松弛时间格式的二维九速度(D2Q9)模型，通过Boltzmann-BGK方程[11]计算流体质点从而获得整个流场宏观物理量。 fα(ri+eαδt，t+δt)=

式(2)中，ri表示某一时刻流体质点的位置，eα表示9个粒子速度；τ为松弛时间，与黏度υ有关[2]，τ=(6υ+1)/2；为fα对应的平衡态分布函数，表示为： 式(3)中，cs为格子声速，有，c取1；ωi为权系数。从而可以知道流体模型的宏观密度、速度为

格子Boltzmann-BGK方程可以通过Chapman-Enskog展开导出N-S方程[12]。LBM算法采用的是格子单位，而FVM用到的是物理单位，本文中根据何雅玲教授[10]提到的单位转换方法进行转换。其中，无量纲参数阻力系数Cd、升力系数Cl和斯托劳哈数St的计算均是不需要转换，它们分别[7]为： ；；

式(6)中，Umax为来流速度最大值；D为特征值尺寸；T为涡脱落周期，Fx、Fy分别为立柱受到的阻力和升力。 1.2 LB拐角边界处理

格子Boltzmann模型的边界采用具有二阶精度的Half-way反弹格式，并结合Machado[13]提出的边界拐角处理方式，这种处理能准确计算出流场分布函数，提高数值精度。由于进出口边界是开放的，只有非曲面的方柱存在拐角，因此，只需在方柱的4个拐角处进行处理，采用内插值和外插值格式。如图1所示，我们知道此时刻拐角处的流体质点的速度uc和分布函数f1、f2和f5。因此通过邻近

粒子插值得到拐角点处粒子的密度ρ，这与平直边界处粒子密度的计算不同。然后根据fα的对称性以及非平衡态反弹格式处理得到余下几个方向的分布函数。为了得到宏观参数和LB粒子群之间较为准确的结果，设f1=0，并满足 即可算出拐角点处粒子的宏观参数。 2.1 几何模型

绕流中，钝体与流场之间有一层很薄的边界层，一些流体参数在该边界层影响下剧烈的变化。因此，在对流场模型进行网格划分的时候，采用分块网格画法，并对方柱周围进行加密，减少边界层影响，保证结果的可靠性和准确性。方柱几何模型如图2所示，其中区域长度为L，高度为H，方柱边长为D，阻塞率B=D/H，立柱中心距上边界为H/2，立柱位置为S。

水流速度定义为来流速度Ux和垂直于来流速度的Uy。流场入口边界采用入流速度形式，出口边界采用自由流形式，流场的上下边界、钝体的边界及内部均为固壁无滑移边界。在自由流中，选择合适的来流速度形式和阻塞比有助于减少流场的不稳性带来的影响[14]。因此，选择抛物线形式U=4Umax(H-y)y/H2的来流速度对雷诺数Re=100进行FVM数值计算，采用UDF编译码导入流体计算软件Fluent中，并与横定值进行了比较，所得结果如图3，表1所示。

当来流的速度为一恒定值时，方柱的Cd值为1.5，抛物线时为1.3，从表1可看出，后者与已有文献值更吻合，说明了抛物线形式的来流速度比恒定的来流速度更具有正确性。

2.2 进出口边界模拟结果分析

为了研究基于FVM下流场边界条件对钝体绕流的影响，现借助Fluent软件对单方柱进行数值模拟，分别研究模型的横向规模L、阻塞比B、与入流边界距离S对流场的影响，所得结果如图4、图5所示。

控制流场模型横向长度为L=50D，方柱中心与入流边界的距离S=12.5D，探究阻

塞比B=1/8，1/10，1/15，1/20，1/25下的方柱阻力系数和升力系数。图5(a)是方柱平均阻力系数随不同阻塞比的变化曲线。从图中可以看出，同一雷诺数下，阻塞比对柱体绕流的影响较大，且随着阻塞比的减小，阻力系数先减小后增大，这一变化规律与文献[17]中实验研究变化规律一致。当阻塞比B=1/8时，方柱升力系数(图4)震荡得比较厉害，B=1/10的时候相对其他的阻塞比很早就开始震荡，因此较高的阻塞比有利于抑制流场的不稳定性[7]。

图5(b)是控制流场模型横向长度为L=50D，阻塞比B=1/10，模拟方柱中心与入流边界的距离S=5.5D，10D，12.5D，25D不同工况下的方柱绕流。从图5(b)中Cd变化曲线可以看出，S=5.5D情况下的方柱阻力系数最大，S=10D、12.5D、25D时的阻力系数大小十分接近，但能看出S=12.5D时的阻力系数最小，即所研究中所得的最小值，并且能有效减小模型尺寸，节省CPU计算时间。

图5(c)是控制阻塞比B=1/10，方柱中心与入流边界的距离S=12.5D，分别取不同区域长度L=30D，40D，50D，60D，得到的不同工况下方柱绕流的阻力系数，获得了计算域大小对阻力系数的影响。从图中Cd变化曲线可以看出，模型的横向规模L对流场的影响不是很大。由于考虑到模型边界对绕流的影响，在取模型的时候，尽量往大的尺寸取，但是太大又会增加计算的时间，降低了工作效率。 3.1 压力等值线分析

从图5讨论的结果来看，所讨论的几种情况下L=50D、B=1/10、S=12.5D这组边界条件更能有效减少FVM模拟下边界对流场的影响。而已有文献[14]表明B=1/8的LBM更能提高流场的稳定性。根据相关研究文献[6]，本文选择LB模型尺寸40D×8D，抛物线式的来流速度，最大速度ULmax=0.1，计算所得压力如图6和图7所示。

结果显示整个流场中方柱上游部分的压力最大，且压力等值线近似椭圆形分布，如图6所示，取y1=H/2，得某一时刻下、一个周期内和同一时刻不同雷诺数下y1

线上压力的变化情况[图7(a)、(b)]。可以看出y1线上的压力随x增加呈指数增加，在前驻点上达到最大值；在一个周期内压力也几乎相等，变化趋势一致；同时，随着雷诺数的增加，压力稍微有所增加，但变化趋势一致，在x=300～360之间压力几乎相同。因此可以知道，当流场达到稳定状态之后，立柱上游的压力不会受到涡脱落影响，且雷诺数对压力较大值的影响也不是很大。 3.2 速度云图对比分析

为研究流场不同位置下的速度，特选取三个位置的速度变化，比较两种方法所得的结果，如图8所示。图8给出了LBM和FVM两种方法在某一时刻所得的速度云图，该时刻一漩涡正从方柱的上边缘脱落，且两种方法模拟下的涡街形态保持一致。在方柱的中心位置、涡正脱落处以及涡已脱落处，我们分别标定x=x1，x2，x3。在这三个位置处我们得到了速度随y向距离的变化曲线和y1=H/2处速度在x方向上的速度变化，如图9所示。

从图9(a)中我们可以看出在图8中x1处，此时刻下方柱迎流面的上下角点处有两个涡正在分离。因此，在方柱的内部速度最小(velocity=0)，而在方柱的上下边缘速度达到最大值(velocity>0.005 m/s)，且速度变化曲线关于流场中轴线(y1=H/2)对称。该现象符合方柱绕流流场变化规律。在x=x2处，关于y1=H/2上下流场的对称性不如x=x1处对称性好，因此该位置处速度变化较为复杂，在流场下半部分漩涡处的速度最大，而在上半部分涡脱落处再次达到第二最大速度。因此我们可以总结出，在漩涡处的速度最大，而在流场其它位置处速度保持均匀，在x=x3处的速度变化也正好吻合这一现象。从图9(b)可以知道中轴线上速度随x向距离的增加先减少后增大，最后趋于动态稳定。结合图8，发现以前驻点为中心向前扩展的椭圆形内的速度呈抛物线型减小，在中心处减到最小。在椭圆至来流入口处的速度保持一致，为入口最大速度(velocity=0.005 m/s2)。 3.3 无量纲参数的对比分析

图10分别是基于FVM与LBM两种方法模拟所获得的单方柱的升阻力系数以及相应的升力频谱图，从图10(a)中可以看出，同一雷诺数、相同抛物线型的来流速度下的数据变化曲线不是严格的正余弦变化，而FVM的数据变化就较为均匀，平稳。这主要是因为LBM基于的是分子动理论，是通过计算大量分子的平均行为得到的宏观量，而FVM直接计算的是宏观参数，因此更具有连续性、均匀性。从图10(b)中发现两种方法所得的方柱升力系数积分值均接近零，但LBM获得的方柱涡脱落升力系数振荡幅值稍小于FVM所得的结果，且前者能更快达到平衡，并在图10(c)中的频谱图中得到了验证。通过表2的数据可以看出，两种方法计算所获得Cd、St值与已有文献对比，结果是比较吻合的，从而证明了两种方法在模拟钝体绕流是准确的。

基于微观分子动力论与粘性流体理论，本文分别用格子Boltzmann方法和有限体积法模拟了雷诺数Re=100下的方柱绕流。用有限体积法研究了不同来流速度形式对流场的影响并探讨了不同边界条件对结果的影响，用格子Boltzmann方法探讨了方柱上游压力变化情况；并比较了两种方法所得的结果。结果显示： (1)抛物线形式的来流速度下的流场更为精确，有限体积法在阻塞比B=1/10时边界对流场影响较小，而格子Boltzmann方法是在B=1/8。

(2)流场达到稳定后，方柱上游的压力变化呈指数分布，且不受涡脱落的影响，雷诺数对压力影响不大。

(3)方柱附近漩涡上的速度最大，大于入口来流的最大速度，在远离方柱的涡街上的速度越来越小，该变化与低雷诺数下的流场涡脱落变化规律一致。 (4)基于有限体积法的软件模拟所得结果的精度要受模型网格的影响，而格子Boltzmann则不需要划分网格，因此在网格处理上占有优势。且在相同的时间步内，后者的计算效率更高。

【相关文献】

1 Pavlov A N，Sazhin S S，Fedorenko R P，et al.A conservative finite difference method and its application for the analysis of a transient flow around a square prism.International Journal of Numerical Methods for Heat & Fluid Flow，2000； 10(1): 6—47

2 Perumal D A，Dass A K.Multiplicity of steady solutions in two-dimensional lid-driven cavity flows by lattice Boltzmann method.Computers & Mathematics with Applications，2011； 61(12): 3711—3721

3 Sohankar A，Norberg C，Davidson L.Low-Reynolds-number flow around a square cylinder at incidence: study of blockage，onset of vortex shedding and outlet boundary condition.International Journal for Numerical Methods in Fluids，1998； 26(1): 39—56 4 Geller S，Krafczyk M，Tölke J，et al.Benchmark computations based on lattice-Boltzmann，finite element and finite volume methods for laminar flows.Computers & Fluids，2006； 35(8): 888—7

5 刘天成，葛耀君，曹丰产，等.基于 Lattice Boltzmann 方法的方柱绕流模拟.同济大学学报: 自然科学版，2008； 36(8):1028—1033

Liu T C，Ge Y J，Cao F C，et al.Simulation of flow around square cylinder basedon lattice boltzmann method.Journal of Tongji University(Natural Science)，2008； 36(8):1028—1033

6 Perumal D A，Kumar G V S，Dass A K.Numerical simulation of viscous flow over a square cylinder using lattice Boltzmann method.ISRN Mathematical Physics，2012； 2012:1—16

7 Kumar G V S，Perumal D A，Dass A K.Numerical simulation of viscous flow over a circular cylinder using lattice Boltzmann method.Proceedings of the 37th National & 4th International Conference on Fluid Mechanics and Fluid Power，India: IIT Madras，2010:16—18

8 Regulski W，Szumbarski J.Numerical simulation of confined flows past obstacles-the comparative study of lattice Boltzmann and spectral element methods.Archives of Mechanics，2012； (4): 423—456

9 陈静涛.圆柱绕流的二维数值模拟和尾迹分析.计算机辅助工程，2013； (6):1—6 Chen J T.2D numerical simulation and wake analysis on flow around circular cylinder.Computer Aided Engineering，2013； (6):1—6

10 何雅玲，王 勇，李 庆.格子Boltzmann方法的理论及应用.北京: 科学出版社，2009:26—29 11 Hou S，Zou Q，Chen S，et al.Simulation of cavity flow by the lattice Boltzmann method.Journal of Computational Physics，1995； 118(2): 329—347

12 Yu D，Mei R，Luo L S，et al.Viscous flow computations with the method of lattice Boltzmann equation.Progress in Aerospace Sciences，2003； 39(5): 329—367

13 Machado R.On pressure and corner boundary conditions with two lattice Boltzmann construction approaches.Mathematics and Computers in Simulation，2012； 84: 26—41 14 Breuer M，Bernsdorf J，Zeiser T，et al.Accurate computations of the laminar flow past a square cylinder based on two different methods: lattice-Boltzmann and finite-volume.International Journal of Heat and Fluid Flow，2000； 21(2):186—196 15 Sahu A K，Chhabra R P，Eswaran V.Two-dimensional unsteady laminar flow of a power law fluid across a square cylinder.Journal of Non-Newtonian Fluid Mechanics，2009； 160(2): 157—167

16 Turki S，Abbassi H，Nasrallah S B.Effect of the blockage ratio on the flow in a channel with a built-in square cylinder.Computational Mechanics，2003； 33(1):22—29 17 Davis R W，Moore E F，Purtell L P.A numerical-experimental study of confined flow around rectangular cylinders.Physics of Fluids，1984； 27(1):46—59