基于速度修正IB-LBM的圆柱强迫运动特性研究*

2023-10-30马欢欢王贯宇涂佳黄

湘潭大学自然科学学报 2023年5期

马欢欢，胡刚，2，王贯宇，涂佳黄

(1.湘潭大学土木工程学院，湖南湘潭 411105；2.重庆市城市建设发展有限公司，重庆 400010)

0 引言

当一个圆柱体结构群长期暴露在流体中时，流体长期的流动可能诱发其产生流致振动和流致不稳定性问题，降低结构物在长期荷载下的稳定性，最终导致结构的崩溃，影响整个系统的安全.早期圆柱绕流问题多以模型试验为主，随着计算机性能的快速提高和数值模拟准确性的提高，有限元、有限差分法、有限体积法等离散方法可以更高效地应用于计算流体力学领域，由于其周期短、成本低，使得它更容易进行大规模的数值计算和分析.

格子玻尔兹曼(lattice Boltzmann method，LBM)方法随即出现，目前被广泛用来模拟多组分流等流动问题[1-2]。与传统流体力学的计算方法比，LBM方法有着无须网格划分、高并行性、可处理复杂边界等突出优势[3].赫万恒等[4]使用LBM方法对均匀来流中匀速转动圆柱的绕流问题进行模拟，不同参考系之间流线涡量连续且模拟结果与实验结果吻合，结果说明LBM方法对于模拟旋转圆柱有很好的应用性.张慎等[5]运用Boltzmann方法对城市建筑群进行小尺度模拟分析，并对其可靠性与准确性进行了验证.刘闯等[6]采用格子玻尔兹曼-大涡模拟(LBM-LES)方法模拟三维圆柱绕流，发现串联多圆柱绕流现象存在一个临界间距比，超过此临界间距比时才会产生明显的涡脱落现象.杨红兵[7]采用多松弛格子Boltzmann方法对三圆柱结构在二维横向流中的绕流进行了数值模拟，发现尾流形状和涡脱落的形状受间距比的作用机理.

目前关于平面均匀来流作用下的圆柱绕流问题，除了单圆柱以及串列、并列的多圆柱体系，人们也陆续开始对圆柱的布置方式进行了实验.Rahman等[8]使用LBM方法，采用两种不同的排布方式，对雷诺数Re=150时三角形排列方柱绕流结构进行数值模拟，分析结构流动特性的影响因素，发现流场特性受排布方式和间距比的共同影响.为了模拟实际工程中更加复杂的流固耦合现象，浸没边界-格子玻尔兹曼方法(IB-LBM)被广泛运用.Wang等[9]采用IB-LBM研究了不同间距比与不同雷诺数下等边三角形布置的3个方柱结构的绕流情况，根据间隙比将尾流分为单钝体、偏转、翻转、反相、反相与同相共存、完全发展的同相流等6种流型.郭春雨等[10]基于IB-LBM对3种雷诺数下旋转圆柱进行数值模拟，研究了不同工况下圆柱尾部流场涡量场的变化，发现圆柱流场特性与雷诺数的关系，随着雷诺数的增加，涡的扭曲现象会更严重.Li等[11]对平面壁面附近圆柱体周围尾流涡旋进行数值研究，发现当间距比L/D增加时，相同符号的单行相干结构变为偶极涡旋散发.Chen等[12]采用浸没边界法对Re=100时不同工况串列三圆柱的涡激振动数值模拟，发现在间距比较大时，下游两圆柱发生大幅振动，而上游圆柱的振动响应类似于孤立圆柱.在实际工程中，结构会受到多种环境变量的影响，通过对多个圆柱体串联或并联排列的流场进行控制，以提高结构群的受力特性.工程结构变化的影响不仅是某个单一特性，还会导致系统的流固耦合特性更加复杂.

综上所述，本文运用具有扩展性强、计算速度快、可模拟多相流体等优势的基于速度修正IB-LBM方法研究雷诺数Re=100时均匀来流作用下等边三角形布置旋转三圆柱的流场特性，主要分析不同的来流角度(α)、旋转速度(ω)和间距比(Kd)3个参数，其中α分别选取0°、30°、60°，ω范围为0≤ω≤2.0，Kd范围为2.0≤Kd≤4.0，对圆柱结构流体力系数和流场特性的内在影响机理，并讨论其转变过程.

1 控制方程和数值方法

1.1 传统的IB-LBM基本方程

将浸没边界法(IBM)与格子玻尔兹曼方法(LBM)组合形成浸没边界-格子玻尔兹曼方法(IB-LBM)，基本方程为：

(1)

本文主要研究的是二维情况，采用常用的D2Q9模型，其平衡分布函数的描述为

(2)

1.2 基于速度修正的IB-LBM基本方程

本文采用了Guo等[13]提出的非平衡外推格式模型，通过对外力项的考虑，解决反弹格式和非平衡反弹格式的局限性，整个Boltzmann方程以下列公式表达出来为：

(3)

(4)

(5)

F=Fxi+Fyj

(6)

式中：外力项Fα为离散的边界作用力分布函数；修正力项F为x方向的分量Fx与y方向的分量Fy的矢量和形式.

设通过密度分布函数求解的速度为u*(过渡速度)：

ρu*=∑αeαfα

(7)

(8)

作用力F通过修正速度解得：

(9)

A修正后的速度可由拉格朗日插值法计算：

(10)

(11)

2 验证和计算模型

2.1 参数有效性验证

图1 串列或并列布置双圆柱绕流算例验证模型示意图：(a)串列旋转双圆柱绕流；(b)并列旋转双圆柱绕流

流场中上、下游圆柱的升、阻力系数平均值随圆柱旋转速度的变化规律呈正对称.将其流体力系数计算结果与文献[14-15]进行对比，如图2与图3所示，结果表明圆柱旋转速度对计算结果差值的影响很小，且旋转速度越大，结果越接近，可证实该数值方法的准确性与适用性.

图2 串列布置双圆柱绕流计算结果与文献的对比情况：(a)阻力系数平均值；(b)升力系数平均值

图3 并列布置双圆柱绕流计算结果与文献的对比情况：(a)阻力系数平均值；(b)升力系数平均值

2.2 计算模型

均匀来流作用下旋转等边三角形布置三圆柱绕流的计算模型以及边界条件如图4所示，坐标原点为3个圆柱的圆心所围等边三角形的重心，3个圆柱的直径均为D，C2的圆心O2与坐标原点O的连线为OO2，来流角度α为入口处与流场域上下边界形成的角度，即直线OO2与x轴之间形成的夹角，雷诺数Re=100，间距比2.0≤Kd≤4.0，来流角度0°≤α≤60°.由于C2在逆时针旋转和较高的旋转速度情况下仍然可以形成强度较高的涡流，流场特性明显，因此设置旋转速度为0≤ω1=-ω2=ω3≤2.0.计算区域为[15D，40D]×[20D，20D]，入口速度为ux=u0=0.05，选择非平衡外推边界格式为流场区域边界，流场的初始密度为ρ=1，初始压力为P=1/3.

图4 旋转等边三角形布置三圆柱绕流的计算模型以及边界条件示意图

流场网格的分辨率会对数值结果有明显的影响，为验证不同网格分辨率对计算结果的影响，分别选取D=20、30、40进行单圆柱绕流模拟分析，结果如表1所示，D越大计算结果越精确，但同时计算效率也越低，综合计算精确性与计算效率，本文选择网格分辨率D为30.

表1 网格分辨率验证

3 计算结果分析

3.1 流场分析

本节重点对旋转三圆柱流场瞬时涡量图随圆柱旋转速度的变化情形展开剖析.

图5展示了间距比为2.0时，不同来流角度与不同旋转速度工况下流场瞬时涡量图的变化特性.圆柱间的间隙流在此时充分地发展，圆柱的漩涡脱落随ω增大逐渐减弱直到消失.

图5 间距比为2.0时，不同来流角度与不同旋转速度工况下流场瞬时涡量图：(a)α=0°，ω=0.5；(b)α=0°，ω=1.5；(c)α=0°，ω=2.0；(d)α=30°，ω=0.5；(e)α=30°，ω=1.5；(f)α=60°，ω=0.5；(g)α=60°，ω=1.5

旋转速度ω较小时(ω≤1.0)，圆柱的旋转使得C1下侧的剪切层附着于C2上，脱离出强度和体积很大的漩涡，侵占C2下侧剪切层漩涡脱落的空间，拉长C3剪切层的长度，从而削弱了漩涡的能量，使脱落的漩涡呈长条状，如图5(a)所示；随着来流角度α的增大，漩涡强度逐渐减弱至消散，后尾流场呈现两行平行的涡街，如图5(f)所示.当ω增大至1.5后，C3的漩涡脱落几乎消失，与C2的下侧剪切层脱落的漩涡在流场后尾流区形成了强度很弱的漩涡带，后尾流场不再呈现两行平行的涡街，如图5(e)与图5(g)所示.特别地，α=0°，ω达到2.0后，C1上下侧的剪切层全部进入圆柱的间隙中，C3被下侧的剪切层全部包住，3个圆柱均未产生漩涡脱落，在后尾流区形成两条平行的剪切层，如图5(c)所示.

图6展示了间距比为3.0时，不同来流角度与不同旋转速度工况下流场瞬时涡量图的变化特性.随着Kd的增大，此时圆柱间的互扰作用有所减弱.

图6 间距比为3.0时，不同来流角度与不同旋转速度工况下流场瞬时涡量图：(a)α=0°，ω=0.5；(b)α=0°，ω=1.5；(c)α=30°，ω=0.5；(d)α=30°，ω=1.5；(e)α=60°，ω=0.5；(f)α=60°，ω=1.5

旋转速度ω较小时，由3个圆柱上侧产生的漩涡在后尾流区逐渐融合形成一个范围很大的正极漩涡，与负极漩涡一起在后尾流区域形成两排平行的涡街，如图6(a)所示；随着来流角度α的增大C2的漩涡脱落减弱，下游圆柱间隙剪切层对下游圆柱内侧的漩涡脱落影响较大，从下游圆柱内侧剪切层分离出的漩涡呈长条状且强度较弱，如图6(e)所示.ω增大至1.5后，受到圆柱旋转和剪切层共同作用，C2下游剪切层的长度变长，漩涡分离十分微弱，圆柱C3受其他圆柱的影响较小，其漩涡脱落的作用逐渐消失，如图6(b)所示；同时随着α的增大，由C2脱落的负极漩涡从长条状逐渐向圆形转变，旋转速度的增大并未削弱其漩涡脱落作用，C2仍能向下游分离出强度较大的漩涡，在尾流区形成两排平行的涡街，如图6(f)所示.

图7展示了间距比为4.0时，不同来流角度与不同旋转速度工况下流场瞬时涡量图的变化特性.圆柱之间的互扰作用大大减弱，各圆柱的剪切层均能够自由地发展并形成漩涡脱落.旋转速度ω较小时，C1下侧剪切层与C2上侧剪切层在C2的后下方相切，削弱了由C1下侧剪切层脱落的漩涡的强度，C1和C3脱离的漩涡强度均较低，在下游不远处就消散，如图7(a)所示；当α=60°，C1的剪切层全部进入下游圆柱的间隙后向下游发展，通过间隙的剪切层受下游圆柱的反相旋转作用被拉扯得很长，阻断下游圆柱剪切层的相互作用，使下游圆柱反相向下游分离漩涡，如图7(e)所示.ω增大至1.5后通过间隙的长条状尾流拉长了C2下侧和C3上侧剪切层的长度，使得其向下游分离出的漩涡强度较低且呈长条形，后尾流区呈现两行平行涡街，如图7(b)所示；而来流角度越大，C2向下游分离漩涡强度越大，旋转速度对其影响越小，如图7(f)所示.

图7 间距比为4.0时，不同来流角度与不同旋转速度工况下流场瞬时涡量图：(a)α=0°，ω=0.5；(b)α=0°，ω=1.5；(c)α=30°，ω=0.5；(d)α=30°，ω=1.5；(e)α=60°，ω=0.5；(f)α=60°，ω=1.5

3.2 流体力系数分析

本节重点对旋转三圆柱的流体力系数变化情形展开剖析.图8展示了不同工况下圆柱的流体力系数平均值(CD-mean)、均方根值(CD-rms)随转动速度的变化特性.

图8 不同间距比、不同来流角度工况下圆柱流体力系数随旋转速度的变化：(a)阻力系数平均值；(b)阻力系数均方根值；(c)升力系数平均值；(d)升力系数均方根值

图8(a)是不同Kd下的阻力系数平均值，Kd=2.0时，随着转动速度ω的增大，C1的CD-mean先减小再增大；C2的CD-mean随着ω增大而减小，且来流角度越大，C2的CD-mean越大；C3的CD-mean受ω影响较小，又C1的尾流对C3的作用最强，导致此时C3的CD-mean相对较小.Kd=3.0时，根据来流角度不同，C1的CD-mean随ω的变化出现两种趋势，α≤15°时，C1的CD-mean与ω负相关直至ω达到1.5后出现小幅增大；当α>15°时，C1的CD-mean与ω正相关.Kd=4.0时，随着ω的增大，C1的CD-mean约为1.3，受ω、α的影响很小；C2的CD-mean逐渐减小，且与α成正比关系；圆柱C3的CD-mean逐渐减小，且与α成反比.

图8(b)是不同Kd下的阻力系数均方根值，Kd=2.0时，随着转动速度ω增大，α越大，C1、C2的CD-rms增幅越大；C3的CD-rms受ω和α的影响较小，当ω=1.0时，C3的CD-rms随着ω的增大达到最大值0.452后逐渐趋于稳定.Kd=3.0时，C1和C3的CD-rms受ω和α的影响较小，均在0附近波动.需要注意的是，α=0°时，由于C1的尾流撞向C2的表面，促进它的涡脱作用，使得此时取得最大值0.529.Kd=4.0时，3个圆柱的CD-rms均在0附近波动.

对于升力系数而言，C1和C3的升力系数平均值(CL-mean)总体均与ω呈正相关，但C2的CL-mean与ω呈负相关，各工况下C1的CL-mean均小于圆柱C3的CL-mean，当圆柱的转动速度为零时，CL-mean几乎为0，且3个圆柱的CL-mean变化趋势受来流角度影响较小，如图8(c)所示.而3个圆柱的升力系数均方根值(CL-rms)受转动速度ω与来流角度α的影响则较大，如图8(d)所示.ω≥0.5后，ω和α越大，圆柱C2的CL-rms越大；来流角度α>30°后，圆柱C3的涡脱受C1尾流的激励变得明显，其CL-rms不断增大.Kd=3.0时，圆柱C2的CL-rms受α和ω的影响较圆柱C1更大，圆柱C2的CL-rms分别在ω=0.5、1.5处达到最大值后逐渐减小.Kd=4.0时，当α<45°，圆柱C1的CL-rms与ω负相关，而当α≥45°，其CL-rms逐渐接近于0；圆柱C2的CL-rms结果相对更大；圆柱C3的漩涡脱落作用随ω增大而减弱，CL-rms与α成正比关系.

3.3 顺/横流向速度分析

以x/D=0处(即三圆柱的圆心所构成的等边三角形的重心)作为记录点，图9展示了不同来流角度与间距比下重心位置处的顺/横流向瞬时速度分布随圆柱旋转速度的变化特征.

图9 重心位置处顺/横流向瞬时速度分布随圆柱旋转速度的变化：(a)α=0°，Kd=2.0；(b)α=30°，Kd=2.0；(c)α=60°，Kd=2.0；(d)α=0°，Kd=3.0；(e)α=30°，Kd=3.0；(f)α=60°，Kd=3.0；(g)α=0°，Kd=4.0；(h)α=30°，Kd=4.0；(i)α=60°，Kd=4.0

同样间距比下，在α=0°时，ω对顺/横流向速度的影响较小，圆柱静止时，顺/横流向速度分布分别呈“U”型分布与幅度很小的反对称分布，当圆柱开始转动时，流场顺/横流向速度分布立即呈竖直分布，受C1的尾流向圆柱C2方向靠近的影响，在y/D为0～1范围内的顺流向速度随ω的增大而减小，如图9(a)与图9(d)所示.随着α增大，同工况下ω对顺/横流向速度分布的影响逐渐变大，且ω越大，顺/横流向速度分布的极值越大，受C1的剪切层发展方向影响，顺流向速度增大，且C1尾流的发展方向更加靠近C3，在重心位置处的顺流向速度增大，使得顺流向速度分布的“凹槽”逐渐向y/D=-1处靠近；由于剪切层在y方向上的速度分量为负，故此时横流向速度越来越小，逐渐接近-1到1的斜线式分布，如图9(c)与图9(f)所示.当间距比达到4.0，α=0°时，在重心位置处的横流向速度均在1附近，且呈直线分布.而α=30°、α=60°时顺/横流向瞬时速度与间距比较小时分布情况类似，特别地，当ω=2.0时，顺流向速度分布出现的“凹槽”靠近y/D=-1处，顺流向速度远大于其他工况，横流向速度在y/D为0～-1.0范围内小于0，在y/D为0～1.0范围内大于0，如图9(i)所示.