周 凯，王 震，陈维山，龙晓军

（1.山东农业大学 机械与电子工程学院，山东 泰安 271018；2.哈尔滨工业大学 机器人技术与系统国家重点实验室，哈尔滨 150001）

0 引 言

格子Boltzmann方法是一种不同于传统数值方法的流体建模和计算方法[1-3]。传统的数值方法的基本思路是将宏观控制方程进行离散，然后利用数值求解的方法去求解该方程，例如有限差分、有限体积、有限元法等。格子Boltzmann方法的是基于分子运动学和统计力学，在微观粒子尺度上建立离散速度模型。在满足质量、动量和能量守恒的基本条件下，通过粒子迁移和碰撞两个微观行为来反应宏观特性，这使得算法简化为对粒子迁移和碰撞的模拟。格子Boltzmann方法采用均匀的正交网格，使得网格算法的工作量大大简化，特别是在处理复杂曲边界时，结合更为准确的边界处理算法和多块网格耦合算法，这种优势更加明显。

多个圆柱的绕流问题在工程中是很常见的，例如海洋平台支撑柱、桥墩和水底管路与流体之间的作用以及它们之间的相互作用。圆柱之间相互作用的存在使得此类问题比单圆柱绕流的情况要复杂得多。这类问题的早期研究多侧重于对实验现象的观察和圆柱力学参数的测量。近代以来数值模拟方法得到了巨大发展，此类问题可以用数值模拟的方法进行研究，并取得了一定研究成果。

Zdravkovich[4-5]等对不同排列方式的双圆柱绕流问题进行了深入的实验研究。针对串列圆柱，研究发现了临界间距的存在。当两圆柱间距小于该临界值时，上游圆柱没有明显的脱涡现象。这一临界值约为3.5倍直径。上下游圆柱尾流场速度分布在该临界间距附近突然发生改变，上游圆柱出现涡脱落现象，两柱之间速度和尾流速度都变大。随着流体力学理论的完善和计算机技术的发展，计算流体力学成为解决流动问题的重要工具。对串列圆柱绕流的数值模拟发现了与实验现象一致的结果，同时获得了更多的流场细节，模拟结果也更接近实验结果[6]。

文章将基于格子Boltzmann方法，对流体的流动和流固耦合作用进行数值模拟。对边界的处理，采用经典的粒子反弹格式，即流体粒子在碰到固体边界后速度反向。针对格子Boltzmann方法的均匀规则网格在处理曲边界问题时的缺陷，采用更为精确的边界处理算法，让流体粒子与固体的碰撞反弹过程在固体格点处发生。为提高计算效率，满足固体周围相对较高的网格密度要求，开发分块网格耦合算法。以单个静止圆柱绕流的已有研究结果为依据，对本文的数值算法进行验证。最后，对不同间距下静止串列双圆柱绕流问题进行数值模拟，研究不同情况下圆柱之间的相互作用以及尾流特征的变化。

1 数值方法

1.1 格子Boltzmann方法

在格子Boltzmann方法中，粒子的运动由波尔兹曼方程描述

其中：u 为流速速度，c=δx/δt，ωα是加权系数，

最终，可以得到粒子的动力学演化方程：

（1）碰撞过程：

（2）迁移过程：

粒子的运动分解为碰撞和迁移两个过程。在迁移过程中，粒子沿eα方向运动，与相邻网格点进行数据交换；在碰撞过程中，粒子与来自不同方向的粒子进行碰撞，在满足质量守恒、动量守恒和能量守恒的条件下，重新构造粒子分布函数。通过（5）式可以看出，粒子的碰撞方程是显式可解的，迁移过程的计算量也较小，算法的编制相对简单。

在粒子的分布函数确定后，流场的宏观量可以表示为

本文采用经典二维九速度模型（D2Q9），将粒子速度离散为9个方向，即α的值有9个，如图1所示。不同方向的速度矢量表达式为

图1 二维9速度（D2Q9）模型Fig.1 The 2D,nine-velocity lattice(D2Q9)model

1.2 多块网格耦合算法

不同于传统的贴体网格处理方法，格子Boltzmann通常采用均匀结构网格来划分流场，网格算法相对更容易实现。然而，在处理流固耦合问题时，固体壁面周围对网格的密度要求较高，而远离壁面处这种要求相对较低，均匀网格显然会造成计算量的浪费。因此，研究多块网格的耦合算法是必要的。本文在保证质量连续和应力连续前提下，采用了多块网格耦合的方法来处理流固耦合中的网格问题。

为简单起见，文章以两块网格的耦合为例进行分析。粗细网格的尺寸比m，而且在网格块的计算中满足 δx=δy=δt，粗细网格交界处的网格结构如图2所示。

前面已经提到，运动粘度和无量纲松弛时间应满足（8）式。根据（8）式，为保证在交界处运动粘度的连续，不同网格块的无量纲松弛时间应该满足（9）式。

图2 不同网格块交界处网格结构Fig.2 Interface stucture between two blocks of different spacing

其中：下标f表示细网格，c表示粗网格。可以将速度分布函数分解为平衡部分和非平衡部分因为在网格交界面上必须满足速度和密度的连续，因此

图3 网格耦合算法计算流程图Fig.3 Flow chart of the computational procedure in the multi-block method

文章采用的多块网格尺寸比m=2，相邻网格会有重合。粗网格和细网格有部分节点是重合的，重合节点处可以利用（18）式和（19）式直接进行数据的交换，其他未知点采用多点插值的方法得出，不同的插值方式对计算结果将产生不同的影响。计算流程如图3所示，n表示当前的时间步，箭头后表示的是该步骤后得到的计算参数。

1.3 曲边界处理方法和作用力的提取

在格子Boltzmann方法中，边界条件的处理方法将直接关系到计算的精度和稳定度。对于静止的固体边界，最常见的处理方法是切向速度为零的无滑移边界条件。具体到格子Boltzmann方法中，入射粒子在碰到边界后，速度方向变为相反方向，这种处理方法也被称为反弹格式。反弹格式在处理静止边界问题上可以得到较好的效果。但是，由于格子Boltzmann方法多采用均匀结构网格，边界处会出现锯齿，只能采用加密网格的方法来提高对边界的拟合精度，因此无法准确拟合曲边界。郭照立等[8-9]结合平直边界的非平衡态外推格式，得到了一种曲边界处理格式，这种格式通过在固体格点位置处执行碰撞过程来确定碰撞后的分布函数，可以获得较好的效果。

曲边界处网格结构如图4所示。以s点代表固体格点，f代表流体格点，w代表固体边界点。本文所采用的方法假定碰撞在s处执行，分析以一个方向的处理为例，其他方向的处理方法类似。在固体点s处，将碰撞前的分布函数分解为平衡态feq和非平衡态fneq两部分

图4 曲边界示意图Fig.4 Schematic description of curved wall boundary

其中，平衡态部分采用粒子的平衡分布函数近似

对于（21）式中的非平衡部分，可表示为

至此，碰撞前固体点的分布函数就确定了。然后，根据LBGK碰撞模型进行标准碰撞，可以得到碰撞后的分布函数

在处理流固耦合问题中，固体的受力参数往往是问题的目标参数，因此流体与固体相互作用力的提取是必要的。本文采用动量交换法来计算流体与圆柱之间的相互作用力，即通过计算流体粒子与壁面发生碰撞前后的动量变化量，根据动量定理来计算粒子所受的流体作用力。在格子Boltzmann方法中，在与壁面碰撞前粒子的动量可以表示为

与壁面碰撞后粒子速度反向，其动量可以表示为

其中的负号表示与α方向的相反方向。

由此，可以计算出流体粒子所受的作用力。根据经典牛顿力学中相互作用力的性质，二者大小相

同而方向相反。因此，通过作用力沿壁面的积分就可以得到固体所受的作用力，具体形式为

进而可以得出圆柱的阻力系数Cd和升力系数Cl：

其中：FD和FL分别是圆柱所受的流体作用力在顺流方向和横流方向的分量，U为来流速度。

2 数值方法验证

为验证算法的可靠性，本文在耦合网格下对Re=200静止单圆柱绕流问题进行了数值模拟，其中雷诺数Re=DU/v，D为圆柱直径，U为来流速度，v是流体的运动粘度。计算区域的划分和网格划分情况如图5所示。在网格密度要求较高的圆柱周围采用尺寸较小的网格块，其他区域采用尺寸相对较大的网格块。 细网格区域取即一倍直径D分为50个格子，粗网格区域网格尺寸是细网格区域的2倍，即网格尺寸比m=2。网格块尺寸间距和时间步大小的确定兼顾了计算的精度和效率。静止单个圆柱绕流问题已经有较多的理论研究、实验测量和数值模拟研究成果，通过与已有结果对比可以验证本文算法的可靠性。

本文数值模拟结果如图6所示。圆柱表面发生了周期性的脱涡现象，正负涡从圆柱上下表面交替脱落。由圆柱升阻力系数曲线可以看出，其升阻力系数曲线呈现周期性震荡，Cd的频率是Cl频率的两倍，这是因为每个周期脱落了一个正涡和一个负涡。

图5 计算区域的划分Fig.5 The division of computational domain

图6 单网格下圆柱绕流计算结果Fig.6 The results of single cylinder with single grid block

表1是本文结果与文献研究成果的比较，可以看出本文结果和文献结果符合得较好。其中的Strouhal数是表征圆柱脱涡情况的重要参数，定义为

其中：f为圆柱的脱涡频率，U为自由来流速度。

表1 结果对比Tab.1 Comparison of results

3 串列双圆柱相互影响研究

对于静止单个圆柱的绕流问题，流场特征和圆柱受力情况主要依赖于雷诺数Re。对于串列双圆柱绕流问题，除了雷诺数之外，两圆柱中心间距比L/D的影响也是关键的（L为前后圆柱中心间距，D为圆柱直径）。

3.1 计算区域和网格划分

本文在耦合网格下对Re=200串列双圆柱绕流问题进行了数值模拟，计算区域的划分如图7所示。计算区域分为两部分，两个圆柱周围采用尺寸较小的细网格，离圆柱较远的区域采用尺寸较大的粗网格。粗网格区域取粗细网格尺寸比 m=2。

图7 计算区域的划分Fig.7 The division of computational domain

3.2 计算结果分析

本文对圆柱中心间距比L/D=1.5、2.0、3.0和4.0的情况进行了数值模拟。不同间距比下的上下游圆柱升力系数曲线、阻力系数曲线、流场涡量等值线和流场流线如图8～11所示。

图8 L=1.5D情况下双圆柱计算结果Fig.8 The results of two cylinders with L=1.5D

图11 L=4.0D情况下双圆柱计算结果Fig.11 The results of two cylinders with L=4.0D

本文将不同间距下上下游圆柱升力系数、阻力系数和Strouhal数进行了汇总，并和文献结果进行了对比，如表2所示。这三个参数是表征圆柱受力情况的重要参数，可以作为判断数值模拟结果可靠性的依据。可以看出，本文计算结果与文献计算值基本相符。

表2 结果对比Tab.2 Comparison of results

综合数值模拟结果可以看出临界间距的存在（3.0＜L/D＜4.0）。当间距小于该临界值时，上游圆柱的分离剪切层附着在下游圆柱上，只有下游圆柱有脱涡现象。此时，下游圆柱的阻力系数较小且小于零。当两圆柱间距超过该临界值时，上下游圆柱都产生明显的涡脱落，且升力系数最大幅值都变大，平均阻力系数也突然变大，特别是对下游圆柱而言，由一个负值变到一个较大的正值。

具体来看，当L/D=1.5、2.0（小于临界值）时，下游圆柱的阻力系数为负值，即下游圆柱并没有受到与来流方向相同的阻力，而是受到了与来流方向相反的推力。L/D=1.5、2.0时，两柱的升力系数幅值都比较小，尤其上游圆柱尾流受到限制，其升力系数几乎为零。当间距比增大到4.0（大于临界值）时，两圆柱阻力系数和升力系数都突然增大。上游柱的升力系数振幅突然增大至0.8左右，下游圆柱的升力系数振幅增至更大的值。下游圆柱的阻力系数由负值变为较大的正值，但其依然比上游柱的阻力系数小，每个圆柱的阻力系数也都小于单柱绕流时的阻力系数。由升阻力系数变化曲线还可以看出，升力振动正负幅值相等，即圆柱所受升力合力为零。另外，上、下游圆柱的涡脱落频率保持相等，但都较单圆柱绕流小。从Strouhal数的变化趋势可以看出，当间距比跨越临界间距值时Strouhal数会突然增大，而随着两圆柱间距的继续增大，Strouhal数接近于单一圆柱绕流情况。

4 结 论

本文基于格子Boltzmann方法，结合多块网格耦合算法和曲边界处理算法，对雷诺数Re=200条件下的静止串列双圆柱绕流进行了数值模拟。数值模拟得到了和已有研究成果一致的结果，验证了临界间距的存在。当两圆柱间距小于该临界值时，上游圆柱没有明显脱涡现象，下游圆柱出现周期性脱涡；当间距大于该临界值时，上游圆柱开始出现周期性脱涡，而下游圆柱所受的升、阻力系数明显提高，并且下游圆柱阻力方向发生了突变。不同间距下两圆柱受力参数和脱涡情况与已有研究成果符合得较好。

[1]Chen S,Doolen G D.Lattice Boltzmann method for fluid flows[J].Annual Review of Fluid Mechanics,1983,30:329-364.

[2]韩文骥,颜 开,褚学森,张 珂.LBM方法数值模拟飞溅流动的影响因素研究[J].船舶力学,2015,19(4):356-362.Han Wenji,Yan Kai,Chu Xueseng,Zhang Ke.Influence factors to numerical simulation of splash flow based on LBM[J].Journal of Ship Mechanics,2015,19(4):356-362.

[3]黄桥高,潘 光.疏水表面流体流动特性的格子Boltzmann方法模拟[J].船舶力学,2016,20(10):1211-1218.Huang Qiaogao,Pan Guang.Lattice Boltzmann simulation of liquid flow characteristics of hydrophobic surfaces[J].Journal of Ship Mechanics,2016,20(10):1211-1218.

[4]Zdravkovich M M.Review of flow interference between two circular cylinders in various arrangements[J].ASME Journal of Fluids Engineering,1977,99:618-633.

[5]Zdrvakoeieh M M.The effets of interference between circular cylinders in a cross flow[J].Journal of Fluids and Structures,1987,l:239-261.

[6]梁亮文.低雷诺数下圆柱横向受迫振荡和涡激运动的数值分析[D].上海:上海交通大学,2009.Liang Liangwen.Numerical analysis of forced oscillation and vortex-induced motion of circular cylinder in cross flow with low Reynolds number[D].Shanghai:Shanghai Jiao Tong University,2009.

[7]Yu D,Mei R,Shyy W.A multi-block lattice Boltzmann method for viscous fluid flows[J].International Journal for Numerical Methods in Fluids,2002,39:99-120.

[8]Guo Z L,Zheng C G,Shi B C.An extrapolation method for boundary conditions in lattice Boltzmann method[J].Physics of Fluids,2002,14:2007-2010.

[9]龚 帅,郭照立.流向振荡圆柱绕流的格子Boltzmann方法模拟[J].力学学报,2011,43(1):11-17.Gong Shuai,Guo Zhaoli.Lattice Boltzmann simulation of the flow around a circular cylinder oscillating streamwisely[J].Chinese Journal of Theoretical and Applied Mechanics,2011,43(1):11-17.

[10]Chan C T,Anastasiou K.Solution of incompressible flows with or without a free surface using the finite volume method on unstructured triangular meshes[J].International Journal for Numerical Methods in Fluids,1995,29:35-57.

[11]Lecointe Y,Piquet J.On the use of several compact methods for the study of unsteady incompressible viscous flow round a circular cylinder[J].Computers&Fluids,1984,12(4):255-280.

[12]Rosenfeld M,Kwak D,Vinokur M.A solution method for the unsteady and in compressible Navier-Stokes equations in generalized coordinate systems[R].AIAA Paper 88-0718,1988.

[13]刘 松,符 松.串列双圆柱绕流问题的数值模拟[J].计算力学学报,2000,17(3):260-266.Liu Song,Fu Song.Numerical simulation of flow past two cylinders in tandem arrangement[J].Chinese Journal of Computational Mechanics,2000,17(3):260-266.