三角形布置圆柱体群绕流特性与流动机理研究

涂佳黄，吕海宇，胡 刚，2，王贯宇

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

0 引 言

在实际工程中，相较于单、双柱体形式结构，多柱体结构群形式运用更广泛，如海洋平台立柱、输电线路和桥桥墩柱等。柱体结构群的绕流问题会更复杂，既要考虑结构与流体的作用，还需考虑结构群之间复杂的互扰效应。以往的研究大多集中于串、并列双柱体结构形式[1-3]。Sumner[1]和Alam[2]研究串列双圆柱结构绕流时，发现了三种流型：拓展体流、再附流和共同脱涡流，三种流型中，下游柱对上游柱的影响在间距比L/D<5.0 时较显著[3]。在对18≤Re≤100 范围内的并列双圆柱结构研究中，Liu 等[4]将其近尾流区分为稳定与不稳定模式。

部分学者对低雷诺数下多圆柱群绕流问题进行了研究[5-10]。Bao 等[5]采用数值计算方法研究了间距比和来流角度对等边三角布置圆柱群绕流特性的影响；Zheng等[6]模拟三圆柱群的绕流，发现小间距比工况下，上游布置并列圆柱比下游布置并列圆柱时的流场流动模式更加稳定；同时，Kang 和Su[7]也分析了两种布置方式下三圆柱群涡脱落频率和流体力系数的变化，也发现上游布置并列圆柱时的尾流呈周期性变化流场更加稳定；杨枭枭等[8]详细研究间距比L/D（1.0≤L/D≤6.0）、Re=100 下的等边三角形布置圆柱的尾流模式，发现6 种模式，按间距比分单体（1.0≤L/D≤1.4）、偏斜（1.4

另外一些学者用实验的方法去研究高雷诺数工况下不同参数对三圆柱群绕流特性的影响[11-15]。Pouryoussefifi[11]测量了Re=6.08×104和1.26×104工况下等边三角布置三圆柱表面压力分布；Igarashi 和Suzuki[12]实验研究了串列三圆柱的绕流特性，将上游圆柱剪切层对下游圆柱影响方式划为六种；Lam和Cheung[13]在2.1×103≤Re≤3.5×103范围内对三圆柱群绕流开展了实验研究，在1.27

综上所述，柱体绕流的研究使用了不同的数值方法或者实验方法，相关的参数涉及了不同的Re、间距比、布置方式等，探究了其尾流模式、流体力特性、涡脱落特征，取得了诸多成果。本文运用浸没边界的多松弛格子Boltzmann（IB-MRT-LBM）数值方法研究Re=100 下等边三角形布置三圆柱的流场特性，该方法可以满足计算高效率、高准确度和较好适用性的要求，同时数值模拟方法可以弥补实验方法对低雷诺数的模拟能力的缺陷，便于更加精细刻画流场的动态过程；主要分析不同来流角度和间距比两个参数（其中来流角度α选取0°、15°、30°、45°和60°，间距比Kd范围为1.2≤Kd≤6.0）对等边三角形布置三圆柱体群尾流模式和流体力系数的影响，阐明流场的转变过程。

1 数值方法和验证

1.1 数值方法

1.1.1 多松弛格子Boltzmann方程

格子Boltzmann 方程在速度、时间、空间上表现出特殊的离散形式，可将微观粒子的速度ξ离散成有限维的速度空间{e0,e1,...,eN}，N表示粒子编号，分布函数f也可进行离散{f0,f1,...,fN}，其中fα=fα(r,eα,t),α= 0,1,...,N。在粒子碰撞过程中，可在不同方向上的分布函数中采用不同的松弛时间进行求解，即多松弛的格子Boltzmann模型（MRT-LBM）：

式中：m是mk组成的向量，mk=φk f,k=1,2,...,b，b等于DrQs模型中速度模型的向量个数s；φk是粒子速度的线性无关多项式；M是φk确定的变换矩阵；meq=Mfeq,m=Mf；feq为平衡态的密度分布函数；S为松弛矩阵，引入D2Q9 模型后，松弛矩阵转换为S=MΩM-1=dig(0,se,sε,0,sq,0,sq,sυ,sυ)-1，其中Ω是碰撞矩阵，由无量纲松弛时间τ=τ0/δt组成。

MRT-LBM算法中，不同于实际有量纲的物理系统参数，密度ρ、u和p均与分布函数有关，即

式中，cs是格子声速，定义为cs=c/ 3，在D2Q9模型中，有cs=1/ 3。运动粘度的定义也与实际的物理系统定义不同：

式中，δt=1/Nt，Nt为时间步数，sυ为松弛矩阵中关于运动粘度的一项。

实际流场中存在固体物时，流固边界的流体质点受到边界力作用使得流场密度分布函数发生变化，郭照立等[16]将额外作用力施加到流场中，结合不含外力项的MRT-LBM方程，得到新的方程：

式中，δt为时间步长，F̂为边界作用力分布函数，I是单位矩阵。

1.1.2 浸没边界-格子Boltzmann方法对外力项的速度进行修正，可以改善传统IB-LBM 方法在流线穿透方面的不足，从而严格满足无滑移边界的要求，使流场演化过程中分布函数精度得到提高[17]。本文将浸没边界与格子Boltzmann 方法相结合，将节点额外作用力加入格子Boltzmann方程，即公式（4）。按照文献[18]的计算思路，使用二阶拉格朗日插值法对修正速度进行求解，提高作用力f的精度。如图1所示，曲线代表实际的固体边界，固体采用拉格朗日坐标系描述，流体则采用欧拉坐标系描述，P为拉格朗日节点，和欧拉节点A之间的距离为Δ1，和欧拉节点B之间的距离为Δ2，通过A1、B1、P点对A和B速度进行修正。拉格朗日节点P的速度（up，vp）为已知，由三点拉格朗日插值可计算出节点A修正后的速度：

图1 流体和固体边界节点示意图Fig.1 Schematic diagram of fluid and solid boundary nodes

同理，欧拉节点B的速度可以由点A1、B1、P进行修正。竖向欧拉节点C和D同样需要D1、C1、Q三点插值修正，如果同一个节点被修正了两次，如图重合的节点B和节点C，那么最终修正速度取两次修正的平均值：

通过对修正速度的求解，可以求出作用力f：

式中，uˉ为修正速度，u*为中间过渡速度。

本文的IB-MRT-LBM 计算流程如下：

（a）设定流场的初始值ρ和u，利用公式（1）计算出初始的平衡态密度分布函数meq；

（b）求出流场域各格点的宏观量ρ和u*，进行边界处理；

（c）由修正速度法，求出边界周围欧拉点的修正速度uˉ，根据公式（7）计算流体与固体边界发生碰撞后的作用力f；

（d）通过作用力f求出F̂，代入公式（4）求得下一时刻新的分布函数；

（e）重复上述步骤（b）与（c），直到结果收敛。

1.2 参数有效性验证

为了验证算法参数设置的准确性，选取单圆柱模型做绕流模拟，对模型的边界处理格式与网格分辨率进行验证。无量纲参数如阻力系数CD、升力系数CL与斯托罗哈数St分别定义为

式中，Fx与Fy分别是浸入边界所受到的阻力与升力，fs是圆柱漩涡脱落频率。选取Re=100，无量纲的入口速度u0=0.05，边界条件以及流场域大小设置如图2所示。

图2 单圆柱绕流计算域以及边界条件示意图Fig.2 Schematic diagram of computational domain of flow around a single cylinder and the boundary conditions

对单位D内包含的格子数分别为20、30 和40 三种工况的单圆柱绕流问题进行模拟，并将其流体力系数和St的计算结果与文献[19-20]进行对比，结果见表1。计算结果之间的误差主要是由于边界格式的精度、网格数量与计算域大小等方面共同影响产生的。另外，在其它参数设置不变情况下，流场网格的分辨率会对数值结果产生一些影响，本节对不同格子数工况下的单圆柱绕流问题进行了计算分析。格子数20工况的CD-mean、CL-rms和St值相对于格子数40工况的误差均明显大于格子数为30工况的数值，但是格子数为40工况计算所需时间会明显增加，因而单位直径D包含的格子数选择为30，可以较好地解决计算精确性与效率之间的问题。

表1 网格分辨率验证Tab.1 Grid resolution verification

1.3 问题描述

1.3.1 计算模型和边界条件

均匀来流作用下等边三角形布置三圆柱绕流的计算模型以及边界条件如图3 所示，三个圆柱分别记为圆柱1（C1）、圆柱2（C2）、圆柱3（C3），其圆心围成的等边三角形的重心为坐标原点O，三个圆柱的直径均为D，任意两圆柱圆心之间的距离与圆柱直径的比为Kd=L/D，OO2是C2的圆心O2与流场域坐标原点O的连线，来流角度为来流在入口处与流场域上下边界形成的角度，大小等于直线OO2与X轴的夹角α。计算区域选取为[-25D，40D]×[-20D，20D]。模型参数选择如下：Re=100，1.5≤Kd≤6.0，α=0°、15°、30°、45°、60°。为保证计算的精确性，设置无量纲速度ux=u0=0.05，uy=0，选择非平衡外推边界格式为流场区域边界，流场出口的密度ρ=1，对应压力p==1/3。

图3 等边三角形布置三圆柱绕流计算域以及边界条件示意图Fig.3 Schematic diagram of the computational domain and boundary conditions for the flow around three cylinders in equilateral triangle arrangement

1.3.2 算例验证

为验证IB-MRT-LBM 算法在模拟等边三角布置三圆柱绕流问题的正确性，本小节选择Re=100，1.5≤Kd≤6.0，α=0°工况进行对比分析。图4展示了等边三角形布置三圆柱绕流计算的圆柱阻力系数平均值CD-mean与升力系数脉动值CL-rms随间距比的变化规律，与文献结果进行了比较，发现本文的计算结果与文献的结果基本吻合，存在一些误差可能是计算方法不同导致的。通过以上验证，说明本文算法在求解等边三角形布置三圆柱群绕流问题方面具有正确性与适应性。

图4 等边三角布置三圆柱群流体力系数计算结果对比Fig.4 Comparison of calculation results of fluid force coefficients of three cylinders in equilateral triangle arrangement

2 计算结果与分析

2.1 尾流模式

本节对不同来流角度与间距比下的三圆柱体群的尾流模式进行了总结，见表2，主要模式有七种：单漩涡模式（SV）、单漩涡体向双漩涡过渡模式（SV-DV）、不规则的双漩涡模式（IDV）、规则的双漩涡模式（RDV）、双漩涡向三漩涡过渡模式（DV-TV）、不规则的三漩涡模式（ITV）和规则的三漩涡模式（RTV）。

表2 各工况下三圆柱群尾流模式分布Tab.2 Wake pattern distribution of three-cylinder group at different cases

小间距比工况下，间隙流发展较弱，间隙内部不存在剪切层，三圆柱群上下侧的剪切层以稳定频率脱落漩涡，因而各来流角度工况下尾流模态均为单漩涡模式（SV）。

当Kd=1.5 时，间隙流流量增多，如图5（a）所示，除了自由流侧有明显漩涡脱落之外，C2 后方负号剪切层和C3自由流侧的剪切层融合后脱落漩涡，呈现单漩涡向双漩涡过渡模式（SV-DV）。间距比增大到2.0，当α=45°时，C1尾流和C2尾流发生了融合，如图5（b）所示，C1和C2间隙侧没有涡脱落，间隙流被夹在中间，流动方向不断变化，C2 和C3 间隙侧的尾流被压缩，呈粗细不均的形态，间隙侧剪切层有漩涡脱落但周期性不明显，为不规则双漩涡模式（IDV）。然而，其它来流角度下，流场会转变为规则的双漩涡模式（RDV），如图5（c）所示，C1 尾流在间隙内部伸长，影响C2、C3 涡脱落，由于没有充足的动量提供，C2和C3间隙侧尾流产生的漩涡被拉伸成条状，但整体的下游涡街正负交替、规则有序。

图5 不同尾流模式下流场瞬时涡量图Fig.5 Instantaneous vorticity diagram of the flow field for different wake patterns

中间距比工况下，C1 自由流侧的剪切层有漩涡脱落趋势，但间隙侧剪切层还未脱落漩涡，如图5（d）所示，层间互扰作用以及下游漩涡变形融合显著，C3 下游负号漩涡向右下方移动，在对流过程中正负融合后，正号漩涡消失，整个流场下游的漩涡形态极不规则，尾流模式是双漩涡向三漩涡脱落的过渡模式（DV-TV）。另外一种尾流模式如图5（e）所示，C1 尾流受间隙流影响很小，C1 下游可以找到负号的漩涡，表明C1尾流完全具备了脱落漩涡的能力，自由流侧的剪切层在脱涡过程中和C2剪切层融合，但此时C2 下游漩涡变形程度显著，使得圆柱受力不稳定性增加，为不规则三漩涡脱落模式（ITV）。

大间距比工况下，如图5（f）所示，各圆柱剪切层均发生摆动，C1 和C2 剪切层产生的漩涡发生融合，与C3 后方的漩涡形成两列正负交替的独立涡街，彼此不受干扰，流场脱涡模式具有周期性的特点，因此尾流模式为规则的三漩涡脱落模式（RTV）。

2.2 流体力系数

2.2.1 时均特性

来流角度的变化会使得流体流动的方向出现偏差，尾流区漩涡也将受此影响，从而柱体间存在不同的绕流特性，然而不同角度下升阻力时均系数（CD-mean，CL-mean）随间距比的变化规律类似，如图6和图7 所示。当α≤30°时，C1 的CD-mean在Kd=1.5 处取得最大值后逐渐减小后稳定，如图6（a）所示。当α=45°和60°时，C1的CD-mean分别在Kd=2.5与1.5处取得最小值。之后CD-mean均随着Kd增大而趋于平稳。需注意的是，在α≤30°工况下C1的CD-mean均会大于α=45°和60°工况，由于当α=15°时C1和C3在横流向接近并列排布，柱体距离较大，C1周边流速会增大，如图8（a）所示，于是C1后方形成了低压区，流线在后方较远处才开始汇合，如图8（b）所示，但是该特征在ɑ=45°时并不明显。同时，与α=45°工况相比，α=15°时C1的剪切层范围会缩小，表明C1后方流体只在较小范围内做旋转运动，其余动量较大的流体则径直向下游运动，见图8（c），导致α=15°工况C1承受了更大的阻力。

图6 不同来流角度下三圆柱体群阻力系数平均值随间距比的变化Fig.6 Variation of time-averaged value of the drag coefficient for three-cylinder group with the spacing ratio in different incoming flow angles

图7 不同来流角度下三圆柱体群升力系数平均值随间距比的变化Fig.7 Variation of time-averaged value of the lift coefficient for three-cylinder group with the spacing ratio in different incoming flow angles

图8 ɑ=15o，Kd=4.0工况下三圆柱体群时均流场特性图Fig.8 Time-averaged flow field characteristic diagram of three-cylinder group for ɑ=15o,Kd=4.0

在小间距比（Kd=1.2、1.5）下，当来流角度较大（α=45°、60°）时，C2前表面偏上方的区域存在一个压力驻点，CD-mean远大于其余工况下的结果，如图6（b）所示。α=30°时，C1 和C2 串列，C2 所受阻力会减小。当2.0≤Kd≤3.0与α=30°时，C2的CD-mean会逐渐减小。随间距比和来流角度改变，CD-mean会出现突变：当Kd=3.0 时，由0.4（α=30°）增加至0.82（α=45°），另外当α=30°时，由0.4（Kd=3.0）增大至0.82（Kd=3.5），此处的机理是流线的分布改变导致C2 的受力方式也发生变化。进一步对其时均流线和压力图进行分析，如图9（a）所示，C1 和C2 的自由流侧流线重合，部分的流线偏转产生了小范围的回流，向下流动过程中与具有较高动量的间隙流相遇，受流体间剪切力作用回流流体方向再次偏转，呈“Z”状的流线，最后流线汇合，C1 对C2 屏蔽作用最明显，C2 所受的时均阻力最小。然而，ɑ=45°工况下，C2 后方形成“S”形流线，前方自由流侧形成一个明显的压力驻点，如图9（b）所示；另外Kd=3.5 时的间隙流量增大，C2后方流体回流后再附着显著减弱，如图9（c）所示，造成了突变的形成。

图9 不同工况下三圆柱体群时均流线和压力图Fig.9 Time-averaged streamline and pressure diagram of three-cylinder group at different cases

随着来流角度的增加，C3 逐渐从上游圆柱转变为下游圆柱，受力发生变化。当Kd≤2.0 时，随着来流角度的增大，C3 的CD-mean会减小，如图6（c）所示。当Kd≥2.5、α≥30°时CD-mean均有明显的上升幅度。特别是α=45°工况下，当Kd=3.5 和4.0 时，C3 的CD-mean达到最大值，由于C2 和C3 在横流向近似并列排布，圆柱体间距离大，来流在经过C2和C3间隙时获得了充足的动量，因而间隙流速度较大，在C3周围形成低压区，如图10（a）～（b）所示。另外，与ɑ=15°工况相比，C3 的剪切层长度会减小，如图10（c）所示，因此C3受较大的时均阻力作用。Kd≥4.5后变化逐渐平稳，α=0°和α=15°所受的CD-mean较小。

图10 α=45°，Kd=4.0工况下三圆柱体群时均流场特性图Fig.10 Time-averaged flow field characteristic diagram of three-cylinder group for α=45°,Kd=4.0

对于升力系数而言，C1 的升力系数平均值变化趋势受来流角度的影响较小。随着间距比的增大，其CL-mean值会逐渐减小并趋近于0。在Kd=1.2 处不同来流角度工况CL-mean值相差较大，如图7（a）所示。然而，当α=60°时，由于三圆柱体群关于来流中心线对称，此时CL-mean维持在0。

在小间距比工况下（Kd≤1.5），随着来流角度的增大，C2 的CL-mean会增加，如图7（b）所示，自由流侧的涡脱落对C2 的影响更加明显。然而，在小来流角度工况下，C2 会被自由流侧的剪切层包裹，C2 的下游涡街远离其表面使得CL-mean减小。在α=15°和α=30°工况下，由于自由流侧脱落的漩涡强度增大，因而Kd=1.5～2.0 范围内CL-mean会增大。在α=45°工况下，CL-mean出现负值，这可能是由于C2 下游涡街较C3更加杂乱，导致涡强度损失大，从而C2受C3下游涡街影响显著。但Kd=4.5时，由于C3下游涡街也受干扰，致使其对C2的影响减弱，使C2的CL-mean减小。

当Kd=1.2时，间隙流没有充分发展，外侧自由流对C3的影响占主要地位，分离剪切层的作用使得C3 承受向下的吸力，CL-mean的绝对值较大，如图7（c）所示。同时，随着来流角度的减小，CL-mean会增大。随着间距比的增大，C3间隙侧也出现分离的剪切层，CL-mean有所减小。然而，当α=60°时，在2.0≤Kd≤4.0范围内，C3又受到向上的流体力作用，使得CL-mean为正值。

2.2.2 脉动特性

本节分析三圆柱体所受流体力系数脉动特性的变化情况。图11 和图12 展示了不同来流角度工况下三圆柱的流体力系数均方根值（CD-rms，CL-rms）随间距比的变化特性。随来流角度的增大，C1 的CD-rms会减小，如图11（a）所示。另外，随着间距比的增大，C1 的CD-rms逐渐减小趋近于0，当Kd=1.2 时，C1 的CD-rms会达到最大值，当来流角度较小时，自由流侧剪切层受到C2 的干扰较小，较少附着在C2 的表面，导致发展越充分，C1的阻力系数波动范围会增大。

图11 不同来流角度下三圆柱体群阻力系数均方根值随间距比的变化Fig.11 Variation of root mean square value of the drag coefficient for three-cylinder group with the spacing ratio in different incoming flow angles

图12 不同来流角度下三圆柱体群升力系数均方根值随间距比的变化Fig.12 Variation of root mean square value of the lift coefficient for three-cylinder group with the spacing ratio in different incoming flow angles

在小间距比（Kd≤1.5）下，α=0°工况流场关于顺流向中心线对称，下游涡街分布规则，如图13（a）与13（b）所示，因而C2的CD-rms大于其它角度下的数值。然而，当α=45°时，Kd增加，间隙侧剪切层形成，漩涡脱落位置远离C2，如图13（d）所示，对C2 的影响会减弱，此时CD-rms减小。随着Kd的进一步增加，CD-rms逐渐增大，再次达到极大值后小幅回落。当Kd≥3.0 时，由于尾流模式的改变，C2 的CD-rms出现突变。在α=30°与Kd=3.5 工况下，尾流呈ITV 模式。另外，在α=60°与Kd=4.0 工况下，C2 的CD-rms发生陡降。间隙流进一步发展，切断了间隙内动量以及涡量的输送，如图13（e）所示，C1 尾流拉伸较长而未卷起，间隙流由于层间的相互挤压速度加快，打断了剪切层的相互作用，脱落的漩涡变形程度大，强度损失较大，因而单看C2 或者C3 下游涡街正负交替特征并不明显，同时C2 和C3 涡脱落相位差约为180°，呈反相模式[9]。然而，当Kd增大至4.5后，C1尾流区有了足够的空间，如图13（f）所示，能够形成完整的涡结构，C2 和C3 涡脱落同相，呈同相模式[9]，尾流模式从反相转变为同相时，也会导致C2 的CL-rms发生较大变化。

图13 不同工况下的流场瞬时涡量图Fig.13 Instantaneous vorticity diagrams of the flow field at different cases

当Kd=1.2时，各来流角度下尾流模式均为SV模式，涡街距C3较近，且漩涡脱落最规则，造成C3的CD-rms最高，Kd增至1.5 后，间隙侧剪切层初步发展，尾流模式转变，漩涡位置距C3 较远，导致其CD-rms减小，Kd≥2.0后变化幅度均较小。

升力系数方面，C1的CL-rms在Kd=1.2处取得最大值，逐渐减小至0附近后趋于稳定。随着间距比的增加，尾流模式有所转变，导致CL-rms发生突变。在α=30°与Kd=3.0 工况下CL-rms曲线的跃升，这是由尾流模式从DV-TV 转变为ITV所引起，C1、C2剪切层融合产生更大的漩涡，导致CL-rms达到0.289，之后的CL-rms均未出现显著变化情况。

如图12（b）所示，α=30°与Kd=3.5工况下C2的CL-rms达到1.05。图13（g）和（h）从机理上展示了突变的原由，Kd=3.0 增至3.5，C1 剪切层卷起后脱涡，C2 后方漩涡在下游形成了更加明显的涡街，对柱体受力影响显著。分析两个工况C2 的CD与CL相位，如图14（a）所示，Kd=3.0 时，相位关系紊乱，CD、CL整体变化幅度较小，CL在0 附近变化。然而，当Kd=3.5 时，尾流模式从DV-TV 转变为ITV，其相位曲线呈规则的“月牙形”，CD、CL的变化幅度显著增加。类似地，在α=45°与Kd=4.5 工况下的CL-rms达到0.75，比较图13（i）和图15（j），C1 剪切层区别很大，Kd=4.0 工况下C1 剪切层一直被拉伸而不脱落，而Kd=4.5 工况下C1 剪切层能够脱落漩涡，受此影响，C2、C3 的剪切层被扭转和拉伸，脱落的漩涡变成长条状，尾流模式也由DV-TV 转变为ITV。从相位分析，如图14（b）所示，Kd=4.0时，CD、CL相位变化呈长条状，变化幅度较小，CL均小于0，表明C2只承受向下的流体力作用。另外，CD增大时CL呈减小趋势，CD减小时则相反，反映了C2受力机理，即CD增大时C2下游的涡吸力效应占主导，CL增大时C3下游涡吸力效应占主导。当Kd=4.5时，由于受下游涡街的影响显著，相位曲线呈“半月形”，CL近似关于零点对称，变化幅度显著增大。

图14 不同工况下C2的CD、CL相位图Fig.14 Phase diagram of CD and CL for C2 at different cases

图15 不同工况下三圆柱群斯托罗哈数随间距比的变化情况Fig.15 Variation of the Strouhal number for three-cylinder group with the spacing ratio at different cases

如图12（c）所示，不同来流角度工况下，C3 的CL-rms呈先减小再增大至稳定的变化规律，其中α=0°和α=60°下C3的CL-rms分别与C1、C2一致。值得关注的是，在α=45°与Kd=4.5工况下，C3的CL-rms突然减小，这是由于随间距比增加，下游圆柱对上游C1尾流发展的抑制效果会减弱，三圆柱群间尾流互扰效应导致C3的绕流发生变化，涡街的演化受到影响。

综上所述，剪切层间的互扰效应影响了漩涡的脱落，同时在间隙流的干扰下，形成了不同的漩涡脱落模式，对圆柱体的脉动力系数产生了不同的影响。

2.3 涡脱落频率

图15 为三圆柱体斯托罗哈数随间距比和来流角度的变化情况。当间距比较小时（Kd≤1.5），流场漩涡数量较少，尾流模式为SV 和SV-DV 模式，三圆柱体的涡脱落频率均会逐渐减小。当α≤15°时，三圆柱体结构的涡脱落频率在Kd=2.0处达到最小，然而，当α>15°时St最小值出现在Kd=1.5处，随着Kd的增大，三圆柱体群结构的涡脱落得到加强，St值均会有较大幅度的提升。图16展示了五个工况的能量谱密度（PSD），在图16（a）～（c）中，C3主峰对应的St值变化不大（分别为0.19、0.185、0.175），表明C3受C1和C2的影响较小。在α=30°与Kd=3.0工况下，C2出现了两个较大的PSD 主峰，其St值分别为0.127与0.19（和C1 相同），α=45°时C2 的两个PSD 主峰接近，并且其中一个St仍和C1 相同，但C1 峰值PSD接近于0，在α=30°与Kd=3.5工况下C1和C2都仅有一个主峰，但C1仍接近于0，α=30°与Kd=3.5工况和α=45°与Kd=3.0 工况下，C1 和C2 间隙内不会发生漩涡脱落。α=45°、Kd=4.5 时C3 的St出现一次突变，发生该现象的原因是C3和C2尾流之间的互扰效应致漩涡无法正常脱落。当α=60°时，Kd=4.0与4.5时C2 和C3 主峰对应的St几乎一致，如图16（d）～（e）所示，Kd从4.0 变化至4.5 后，C1 剪切层尾部能够卷起，以C2、C3 相同的频率脱落漩涡，同时，三圆柱主峰所对应的能量均提升较大，表明当Kd=4.5 时C1尾流能够干扰圆柱的漩涡脱落方式。随Kd继续增大时，剪切层间与下游涡相互作用逐渐减弱，后方涡街更加规则，呈现RTV 模式，涡脱落频率最终趋于稳定。α=0°工况下涡街演化更规则有序，圆柱间的互扰效应以及流场下游的漩涡融合变形最弱，因此在α=0°与Kd≥3.5下，C1和C3的St值均较大。除α=60°工况，C2 始终位于顺流向的下游，其和上游圆柱的相互作用更强，相较于C1 和C3，C2 的涡脱落更容易受到间距比和来流角度的影响。

3 结 论

本文基于IB-MRT-LBM 研究了Re=100下等边三角布置三圆柱体结构群的绕流特性，分析了不同来流角度和间距比对三圆柱体群尾流模式及其流体力系数的影响，阐述了其流动机理。主要结论如下：

（1）不同来流角度下尾流模式分为SV、SV-DV、RDV、IDV、DV-TV、RTV、ITV 七种模式。α=30°时下游涡街更加规则，且尾流模式转变更快，α=45°和α=60°时尾流模式则转变较慢。各来流角度下尾流模式均存在从不规则向规则转变的过程，Kd≥5.5后尾流模式均为RTV。

（2）流体力系数时均值受间隙流流速影响较大，除了小间距比，各来流角度下随间距比变化趋势基本一致，α=45°时在较大间距比范围内出现了升力系数均值为负的情况。在Kd≥3.5 后，尾流模式的转变对C1 的流体力系数均方根值有明显影响。同时，剪切层的摆动、融合等互扰效应也会改变流场的脉动特性。

（3）不同来流角度工况下，在Kd≥3.5 时C1 和C3 的St受Kd的影响较小，数值变化相对稳定。α=0°工况下，当Kd≥3.5时，C1和C3的St值均较大。除α=60°工况，C2的涡脱落更容易受到间距比和来流角度的影响。