赵红晓

（同济大学航空航天与力学学院，上海 200092）

0 引言

圆柱的绕流现象在工程中广泛存在，如管排［1-2］、桥墩［3-4］、电缆和高耸建筑［5］等结构周围的流动问题。1940 年11 月7 日，美国华盛顿州塔克马海峡大桥因风振致毁，事故原因正是空气绕流形成的卡门涡街引起的桥梁共振。作为典型钝体的一种，近代学者对圆柱绕流的研究一直没有停息过。

1911 年，冯·卡门观察到在一定条件下的定常来流绕过某些物体时，物体两侧会周期性地脱落出旋转方向相反、并排列成有规则的双列线涡，并对涡街的诱导速度、稳定性和阻力等做了分析［6］。Coutanceau等［7］揭示了圆柱后尾流形态和雷诺数之间的对应关系。Zdravkovich［8］研究了双圆柱串列在流体流动速度很小时，不同的间距比L/D（圆柱间距与圆柱直径之比）对两柱间隙中的速度分布和下游柱尾流场中的速度分布具有明显影响。Le 等［9］发现双圆柱绕流在间距比L/D=2 时尾流产生的涡街迅速消失，间距比为4.5 时两条尾流的形状随时间随机发生变化，出现涡脱落方向相同、相反等4 种典型的尾流形状。Yokoi等［10］研究了双柱并列（非等直径），间距比为0.25 ～1.5 之间时流体形态发生交替偏斜，间距比大于1.5时，两柱后分别形成各自的涡脱落。张爱社等［11］利用有限元法对等边三角形布置的等直径3 圆柱绕流问题进行了研究，表明小圆柱的位置、大小、个数以及雷诺数变化等对主圆柱有影响。Kuo 等［12］研究了2 个并列小圆柱置入尾流区，可有效地抑制涡街的产生，并使主圆柱的升力系数降低了70% ～80%，压力、阻力减小也很明显。张洪军等［13］研究了附属小圆柱对并列双圆柱尾流流场的影响。Muddada 等［14］研究了放置于主圆柱后的旋转小圆柱可以有效地抑制涡街的生成。

由于流动的复杂性，目前以多圆柱群为对象的研究相对较少，工程情况中圆柱以圆柱群的形式出现的案例却十分常见，本文选取不等直径多圆柱群实验装置进行仿真研究。

1 数值模拟

1.1 控制方程

本文采用Fluent软件进行数值模拟，选用Simplec算法（半隐式方法Semi-Implicit Method for Pressure Linked Equations）进行计算。本文研究中，低雷诺数下水的流动，可近似为不可压缩流体，则二维连续性方程可简化为：

动量方程为动量守恒定律在流体力学中的表示形式：

式中：u、v、w为3 个速度分量；p为压强；ρ为流体密度；μ为流体动力黏度。

升力系数Cl和阻力系数Cd是流体力学中的无因次量，用来表示物体在流体中的升力和阻力大小

式中：Fl、Fd分别为物体受到的升力和阻力；v∞为无限远场的来流速度；S为物体在来流方向上的投影面积。

雷诺数

式中，v、L分别为流场的特征速度和特征长度。雷诺数物理上表示惯性力和黏性力量级的比。

本文均设置左侧边界为速度入口，右侧边界为压力出口，远场壁面为对称边界条件，圆柱壁面采用无滑移壁面。流体介质为水，密度ρ=1000 kg/m3，动力黏度μ=1.0 ×10-3kg/（m·s）。

1.2 计算模型

多圆柱绕流研究实验装置如图1 所示，以还原实验装置中的流动规律，划分模型进行多组仿真计算。

图1 流态演示装置断面示意图

为避免过于庞大的模型造成不必要的计算时间浪费，需对模型进行适当的简化和划分。

如图2 所示，模型保留了实验装置中近壁面绕流和多圆柱群绕流的流域，并在入口和出口处做了适当的加长。与实验室现象进行对比，观察能否正确模拟实际流动，保留了实验装置中的栅格，多圆柱排列方式、尺寸和编号如图3 所示。

图2 实验装置物理模型及尺寸

图3 多圆柱群尺寸及编号示意

1.3 网格划分及计算方法验证

采用结构化网格的方法进行网格划分，建立oblock与圆柱面以及栅格前缘、尾缘的圆角关联，将其余边与壁面和栅格长边关联，并在湍流情况中，运用式（7）估算第1 层网格到柱面的距离Δy进行相应加密。最终得到的网格具有较好的质量和长宽比，网格划分如图4 所示。

图4 实验装置网格划分

网格划分与层流验证模型时相同。为保证圆柱面上第1 层网格满足y+≈1，估算确定第1 层网格到柱面的距离［15］

代入特征长度L=D=0.01 m，Re=3900，得到Δy=0.03 mm，加密后网格如图5 所示。

图5 圆柱周围及局部网格加密图

2 层流模型和紊流模型验证

参考文献［16］，使用相同的物理模型进行层流验证。左侧为速度入口，右侧为压力出口，远场壁面为对称边界条件，圆柱壁面采用无滑移壁面。计算区域20D×25D，采用结构化网格，计算参数见表1，计算结果见表2。

表1 层流验证模型仿真参数

表2 层流模型验证结果与文献数据对比

由表2 可知，计算结果与文献［15］基本一致，可使用层流模型对小雷诺数的情况进行模拟。

参考文献［15］中使用相同的物理模型进行紊流验证。左侧为速度入口，右侧为压力出口，远场壁面为对称边界条件，圆柱壁面采用无滑移壁面。计算区域20D×40D，靠近圆柱加密网格，采用结构化网格，计算参数见表3，计算结果见表4。

表3 湍流验证模型仿真参数

表4 湍流模型验证结果与文献数据对比

计算结果与文献［15］基本一致，可使用分离涡模型进行模拟。

3 仿真计算

分别采用以上验证计算模型和划分网格的方法分别对多圆柱绕流模型（见图3（a））进行数值模拟，表5中雷诺数均为小圆柱雷诺数，参数设置见表5。

表5 实验装置中多圆柱绕流仿真参数

由速度云图图6（a）可见，在Re=200 时，以大圆柱的前缘为中心，水流在驻点受阻后，被迫向两边流动，在靠近圆柱尾部的边界上，主流开始与圆柱体分离，边界层分离后，在分离区的下游形成回流区，即为尾涡区。后流经栅格区域，在流经栅格的各狭缝中，由于截面积的减少，流动皆有一定的加速。流体流经多列小圆柱时，列与列之间的尾涡出现掺混而紊动更强烈，涡的长度由于受到限制变短。第4 列圆柱则直接在后方产交替脱落涡。

由图6（b）速度云图可见，当Re=5000 时，大圆柱尾涡出现明显的两列涡，且脱落形成长度明显变短，涡街的宽度变窄，尾流呈复杂的湍流特征涡。小圆柱的第1 列圆柱限制了后面第2、3 列圆柱分离剪切层的卷起，使其不能产生强烈的旋涡脱落，第3、4 列圆柱的涡在后方混合，形成更为复杂的尾流。各圆柱后方仍有一定范围的回流泡，旋涡脱落现象在下游第4 列圆柱后尾涡脱落成偏斜的尾流模式，随着雷诺数的增大，第4 列后的紊动就越强烈。

当Re=10000 时，从图6（c）的速度云图可观察到其流动与Re=5000 时十分相似，但紊动更强烈，尾流模式和圆柱受力不具有周期性，是由于尾涡区的长度和紊动强度与来流的雷诺数有关，雷诺数越大，紊动越强烈。

当Re=200 时，由图7（a）、（b）可见，大圆柱受到最大阻力系数最大，小圆柱中第1 列（从左向右看）圆柱阻力次之，第3 列圆柱所受阻力略小于第1 列圆柱，位于中心线上的圆柱12、32 仍相比同列受到更大的阻力。第2 和第4 列圆柱所受阻力比第3 列略小，第4列圆柱所受阻力略大于第2 列。第4 列圆柱升阻力系数的振幅比前几列大。

图7 多圆柱绕流阻力系数随时间变化曲线

当Re=5000 时，由图7（c）、（d）和（e）可见，大圆柱的最大阻力系数仍是最大，但比Re=200 时小。圆柱12 阻力系数次之，第4 列圆柱的阻力系数比第3 列大，位于中心线上的圆柱12、32 仍相比同列其他圆柱受到更大的阻力，各圆柱的阻力系数均有一定幅度的振荡，后方的圆柱较前方有更大的振荡幅度。

当Re=10000 时，由图7（c）、（d）和（e）可见，各圆柱的阻力系数均有一定幅度的下降和振荡，后方的圆柱较前方有更大的振荡幅度。同样位于中心线上的圆柱12、32 仍相比同列受到更大的阻力。

由图7 中不同雷诺数下的阻力系数可见，大圆柱的升力系随着雷诺数的增加而降低，小圆柱的绕流阻力系数总体上也随着雷诺数的增加而降低，随着雷诺数的增加，流动由层流变为紊流，湍动更强烈，各阻力系数均随时间产生微小振荡，振幅逐列增加。

当Re=200 时，由图8（a）～（d）可见，第4 列小圆柱升力系数大约在±4.5 这一范围内振荡，其次是圆柱31、33，圆柱11、13。圆柱21、22 升力系数在零上下振荡，对称轴上的大圆柱和小两圆柱12、32 所受升力为零。第1、3 列和第2、4 列同侧圆柱所受升力仍具反向性。

图8 多圆柱绕流升力系数随时间变化曲线

当Re=5000 时，由图8（e）～（h）可见，升力系数的振荡幅度也逐列增加，第4 列圆柱的升力系数振幅较大，最大值为其他圆柱的2 倍以上，第3 列升力系数次之，第1 列最小。不过前2 列圆柱的振荡较小，除去中心线上的圆柱12、32，其他圆柱受力方向也没发生改变，而大圆柱和小圆柱第4 列圆柱随着尾流的复杂化，其升力系数开始较大幅度的振荡，且方向会随时间改变，说明紊动强度增加，能量损失更大。

当Re=10000 时，由图8（i）～（l）可见，规律和雷诺数5000 时相似，升力系数的振荡幅度较之雷诺数5000 时更大，流体紊动强度继续增加。

对比Re=5000 时的速度云图和实验装置的流动图像（见图9），可发现两者较好符合，认为仿真结果与实际流动是很接近的，可用数值模拟的方法替代实验装置来演示紊流的流态。

图9 实验装置流动图像对比

4 结语

本文对实验装置简化模型在3 种雷诺数下进行了数值计算和流态、阻力系数和升力系数分析，可知，多圆柱绕流在雷诺数Re=5000 时的速度云图和实验装置的流动图像符合，认为仿真结果与实际流动是很接近的，可采用模拟的方法去演示流体的紊动机理，为流体力学虚拟仿真开发提供有益的探索。