宽高比为0.62的矩形柱绕流特性研究

李 鳌，程文明，王书标，杜 润

(1. 西南交通大学机械工程学院，四川 成都 610031；2..轨道交通运维技术与装备四川省重点实验室，四川 成都 610031)

1 引言

钝体绕流作为风工程领域基础问题，一直深受国内外学者广泛关注。在风洞试验和数值仿真方面有大量的研究资料。Lyn，Einav等[1]应用LDV(激光多普勒风速测量技术)对方柱绕流进行风洞试验。Norberg[2]在低速闭环风洞对不同宽高比矩形柱绕流进行了系统研究。秦浩，肖姚等[3]结合PIV(粒子图像测速法)和POD(正交分解法)对方柱绕流尾流场进行了解析和研究，结果与Lyn，Einav等[1]的吻合较好，同时应用γ-Reθ和SST k-ω数值模拟也得到了不错的结果。数值研究方面，国内外学者就方柱绕流[4-8]运用LES(大涡模拟)仿真，均获得不错效果。Sohankar[9]采用LES对雷诺数为Re=105的宽高比为0.4～4的矩形进行仿真分析，得到的结果和Norberg[2]实验结果相一致。Tian，Ong等[10]采用SST k-ω对宽高比0.05-1的矩形进行分析，但对宽高比小于0.6时结果不理想。Islam，Zhou等[11]应用ILBM(不可压缩亚格子玻尔兹曼模型)对低雷诺数(Re≤250)矩形进行分析，得到与高雷诺数不同的结论。Mannini，Soda等[12]对宽高比5.0的矩形讨论了风向角和雷诺数的影响，并与风洞结果吻合较好。

实际生活中，矩形柱和类矩形柱的结构很多，例如起重机箱型梁结构，针对其展开研究意义明显，而试验和仿真也均表明：矩形柱绕流在阻力系数上存在一个临界值，即宽高比0.62附近时阻力系数最大。而对于其机理产生，也作了各种解释和分析[13， 14]。但是单独对其进行分析和研究的讨论略显不够。

因此，本文主要针对该临界情况进行深入分析。选择Norberg[2]文中雷诺数Re=19000作为研究条件，采用大涡模拟方法作为主要手段，通过矩形柱气动特征，流场及尾涡特性等等进行分析和讨论，研究阻力和结构关系和机理，并进行详细地分析与讨论。

2 湍流模型与计算方法

2.1 大涡模拟模型

对不可压缩流体控制方程N-S方程而言，经过滤波得到LES的连续性和动量方程可以表示为

(1)

(2)

(3)

其中，νSGS为亚格子模型涡粘系数，表达式为

(4)

(5)

2.2 计算流场及算法

针对R=0.62矩形柱的数值计算仿真的流场域及相应的边界条件如图1所示(针对二维模拟展向L为0)。

计算域在空间各方向无量纲长度(即长度与特征高度的比值)Lu，LD，h，H，L分别为12，30，1，20，4。本文展向长度方向选取的是L/b=6.45(L/h=4)。针对展向长度L/b的选取，国内外学者对此也进行了很多相关性实验和仿真研究。Vickery[15]，McLean & Gartshore[16]等通过实验数据验证了方柱展向长度相关性长度大约在L/b=5.0，Tamura，Miyaji等[17]通过仿真计算也得到展向长度在L/b=4.0～10.0时结果差异不大的结论，周强，廖海黎等[7]也针对方柱展向长度进行分析和讨论。

计算边界条件的选取如下：入口处采用inlet速度入口边界，均匀来流速度条件：u=U， v=w=0；出口采用零压力出口：拟压为0；柱体壁面采用无滑移壁面：No-slip；上下端面采用对称边界条件：Symmetry；展向边界采用周期边界条件：Periodic。

图1 计算流场域和边界条件

在空间离散上，采用有界的二阶中心差分格式，时间离散上，采用有界的二阶隐式法，对压力速度耦合采用PISO求解算法；收敛精度设置为10-5，无量纲时间步长为ΔtU / D=2.78×10-3，以确保库朗数小于1，且迭代在20次内完成。

流场域选用六面体网格进行空间离散。三维流域在XY平面采用混合网格，即柱体表面 5h 范围内采用O型网格，外围采用正交型网格，并沿着展向z方向进行延伸。贴体O型网格采用较低网格增长率(1.05)以考虑边界层效应；外围则采用较大网格增长率(1.1)在保证计算合理基础上减少计算成本。网格在XY方向示意图如图2所示。二维流场及网格设置。

图2 计算流场域内XY平面内网格划分示意图

3 网格无关性验证

在进行绕流分析前一般需要进行网格无关性验证，这里选用3套疏密不同的网格，网格差异主要在于表面首层网格高度的不同，Case1、Case2 和Case3沿着展向方向解析度为Δz/h=0.1，首层网格高度设置有所差异，具体参数如表1所示。将三组不同网格计算得到的平均积分分量结果与其余相关实验和数值结果进行了对比，结果如表1所示。对比三组网格Case1、Case2 和Case3，结果间差异不大，且与参考值间吻合较好。综合考虑计算精度和成本，选用Case2结果作为参照。此外，表1中还列出了采用二维LES计算得到的结果，即Case4数据作为对比。

表1 计算结果对比

4 结果与讨论

以高雷诺数Re=19000时，临界宽高比R=0.62下的矩形柱为研究对象，结合通用的LES湍流模型对矩形柱结构下的气动特征，绕流特性及其产生机理等进行研究和分析，同时对二维和三维LES进行对比研究。

4.1 气动特性研究

表1对比了计算所得三分力系数及频率参数等与前人试验和仿真所得结果以及方柱绕流对比情况，所得三分力系数可直接通过计算所得时程曲线(如图3所示LES计算时程曲线)取平均值得到，频率参数St 数需要对升力系数结果进行FFT变换得到涡脱频率 f，再根据计算公式St=fh/U 得到(h 为特征高度，U为来流速度)。

图3 三维LES升力系数 Cd 和阻力系数 Cl 时程曲线

由表1结果可知，三维LES计算得到Case2的平均阻力系数 Cd与参考仿真值[9， 10， 18]间的差异在1%，与实验值[2， 13]吻合也很好； St 数与实验值[2， 13]均比较吻合，相对误差也在2%左右；而升力系数均方根 Clrms与LES仿真结果[9]吻合，远小于二维RANS仿真[10， 18]结果。而Case4为二维LES计算所得，计算所得Cd和St 数与三维LES结果相似，均和参考值间吻合很好；而升力系数均方根 Clrms高于Sohankar[9]采用三维LES所得结果(相对偏差28.4%)，同时又小于RANS计算[10， 18]所得结果(相对偏差分别为51.5%和39.1%)。

表1对比结果说明气动力参数的预测在阻力系数和频率上，采用各模型计算所得结果都比较好；而在升力系数均方根上，采用三维LES结果最好，其次是二维LES，采用二维RANS计算预测效果最差。

图4为柱体截面表面平均压力系数 Cp的分布情况，同时与方柱表面压力分布的实验[19]及仿真结果[7， 8]进行了对比。由图可知，柱体前端面AB面和上端面BC面，矩形柱和方柱压力分布差异不大，差异主要发生在后端面CD面：矩形柱后端面压力分布明显小于方柱，且在D处差值最大。而后端面压力分布反映的是绕流阻力的大小，也表明了该矩形柱阻力系数更大，与表1中Cd结果相呼应。

同时，BC段上端面压力系数分布可知，三维LES结果大于二维LES，即二维LES上端面负压更大，而上下端面压力分布反映的是绕流升力系数相关参数，表明二维仿真二维LES升力参数估值过大，与表1中Clrms结果相呼应。也说明了Clrms是一个与三维相关的参数，与周强等[7]方柱绕流结果相吻合。

图4 矩形柱中截面表面压力系数Cp分布

4.2 流场特性研究

图5给出并对比了中心平面(y=0平面)上平均流向速度的分布情况(横坐标X-position=x/h，下同)。有图可以看出，矩形柱和方柱一样，流向速度在柱体后壁面速度为零，然后逐渐减小到最小速度Umin，再单调增加到与来流速度 U 相同为止。相对于方柱，矩形柱最小速度比方柱更小，最小值到结构中心截面(x=0平面)距离更短(定义为绕流回转长度)，说明矩形柱后的尾涡距离更短，即矩形柱背压更大(如图4)，因此阻力系数更大(如表1)。同时，尾流区内，矩形柱流向速度在 x/h > 6时约为0.4U小于方柱的0.4U，说明矩形柱尾流区内尾涡区域更长。

对比二维LES和三维LES结果发现，二维LES计算结果所得的回转长度略小，说明二维LES计算阻力系数更大，与表1计算所得气动参数相对应；最小流速 umin/U 相较于三维LES偏大，说明展向长度对umin/U有较大影响，对回转长度影响较小，所得结果存在一定偏差。

图5 矩形柱中心线平均流向速度分布

图6和图7分别为尾流区内距离柱体后壁面分别为x/h=1，2和3.5时流向速度，竖向速度的平均值及均方根，以及对应的雷诺应力。通过对比不同位置处速度及其脉动量的差异可以反映尾流区内尾涡变化情况，同时对比方柱和矩形柱在相同位置处速度及其脉动值差异可以反映不同结构对尾流区的影响，也可以对比二维和三维仿真在尾流区的差异情况。

图6(a)为尾流区平均流向速度不同位置速度分布，速度分布以y=0处为中心呈现“U型”的对称分布，且随着远离尾流区逐渐变缓；图6(b)为平均竖向速度分布，分布以y=0处为中心呈现“N型”的反对称分布。与方柱试验及仿真结果对比发现，矩形柱在y=0附近速度相较于方柱更小，整体上“U型”更窄，“N型”更缓，这些说明了尾涡区内矩形柱的尾涡宽度更窄，间距更短，也即印证了矩形柱阻力系数更大，和表1结果相呼应。同时，对比二维LES和三维LES结果二者平均速度差异不大，结果基本一致，说明展向方向的长度对平均速度流场影响不大。

图6 尾流区不同位置平均速度分布(—，Present-Rectangular-3DLES；---，Present-Rectangular-3DLES；，Exp-Square(Lyn et，al，Re21400)；--×--，LES-Square(Zhou et，al，Re22000))

速度平均值反映的是尾涡中结构的时均特性，而脉动值和雷诺应力则可以反映结构的湍流特性，如图7为尾流区内的速度均方根及雷诺应力分布。图7(a)为来流速度均方根分布，由于卡门涡街的存在，存在两个峰值，呈“M型”对称分布于中心点两侧，同时随着远离后壁面而变缓并消失，说明旋涡在下游逐渐减弱；图7(b)为竖向速度均方根分布，在中心线处呈现“∩型”对称分布；图7(c)为雷诺应力分布，呈“N型”反对称分布。

对比图7中方柱和矩形柱脉动参数可以发现，矩形柱的流向速度均方根分布在x/h=3.5时“M型”双峰依旧明显，说明旋涡在流向速度方向相较于方柱衰减更慢；竖向速度方向“∩型”相较于方柱在靠近后壁面时峰值稍大，而远离柱体后趋于一致；雷诺应力上矩形柱和方柱在各个位置上表现差异不明显。以上说明方柱和矩形柱在湍流特性上的差异主要体现在流向方向上，矩形柱相较于方柱尾涡更大，且更接近于后壁面，同时尾流区更长。

图7 尾流区不同位置速度均方根分布

同时，图7中二维LES和三维LES结果对比发现，二维LES的仿真结果变化幅值更大，趋势更明显，说明二维LES仿真在尾涡模拟上会进行一定程度上的放大，也说明展向长度对速度脉动值模拟，即尾涡湍流特征的描述上有一定影响。这也解释了表1中二维LES得到升力系数均方根偏大的原因。

通过对宽高比0.62的矩形柱的三分力系数，表面压力系数分布，流向速度分布，尾流区不同位置时均速度，湍流特性等等结果的与方柱间的对比分析和讨论，解释说明了该临界宽高比下矩形柱的特殊性和相应的产生机理。而二维和三维计算结果的对比也说明进行三维计算的必要性。

5 结论

针对高雷诺数Re=19000下，临界宽高比为0.62的矩形柱，采用大涡模拟(LES)方法，就二维和三维模型进行绕流数值模拟。在验证网格及数值结果准确性的基础上，详细分析了矩形柱的气动力和流场特性，并与方柱结果进行分析对比，得到以下结论：

1)采用三维大涡模拟算法，对宽高比为0.62的矩形柱进行绕流分析，得到气动参数等和相关实验和仿真数据吻合很好。同时对比矩形柱表面压力系数，中心线来流速度分布，以及不同位置处平均速度以及湍流特性的分布，可以很好地解释临界宽高比0.62气动参数的特殊性。

2)对比二维和三维大涡模拟计算结果，可以发现：二维大涡模拟计结果在升力系数均方根，柱体上下壁面压力系数分布，尾流区最小速度，速度脉动值分布等的预测上偏差较大，说明这些参数对绕流结构的展向长度较敏感，三维效应明显。

3)考虑到计算效率和计算成本，如果在只关心结构阻力系数以及涡脱频率时，选用二维计算更合理。

至于采用三维计算时，展向长度对矩形柱绕流特性的影响程度和规律，还需要后续深入研究探讨。