二维方柱绕流阻塞效应的大涡模拟

郜 阳， 全 涌， 顾 明

(同济大学 土木工程防灾国家重点实验室， 上海 200092)

在对建筑结构进行风洞试验研究时，阻塞效应可能导致风洞内的流场特性偏离真实流场，从而影响试验结果的准确性.深入研究阻塞比(Blockage Ratio，以下简称R)对流场以及柱体气动力的影响对保证风洞试验结果的正确性是至关重要的.

研究人员[1-4]已采用风洞试验的方法，就阻塞比对模型风压分布、气动力等气动特性的影响展开了大量的研究.大多数的学者，如黄剑等[3-4]，均是通过对相同截面形式、不同缩尺比的模型进行风洞试验，来研究阻塞效应.而少数学者[1]考虑到模型尺寸对风洞试验结果的影响，通过改变在风洞内安装的隔板的位置实现不同阻塞比，保证了模型的尺寸不变，从而对阻塞效应展开研究.尽管已有大量试验研究，但对阻塞效应的定量认识仍然比较模糊.

利用计算流体力学(computational fluid dynamics，以下简称CFD)开展阻塞效应的研究最早开始于1980年代，Stafford[5]通过数值方法验证了风洞洞壁形状的改变对降低阻塞效应的可行性.之后，Davis[6]对两种阻塞比(16.7%和25%)的矩形柱进行数值模拟，研究表明平均阻力系数以及Strouhal数随着阻塞比的增大而增加.Okajima[7]、Sharma[8]和Patil[9]等通过数值模拟方法考察了方柱阻塞效应，研究了平均阻力系数和脉动升力系数随阻塞比的变化规律.值得注意的是，这些关于阻塞效应的数值模拟均是在二维平面内建立计算域，并且局限于低雷诺数范围内.相比于二维平面，三维空间的数值模拟更为复杂，但也更接近于实际.遗憾的是，到目前为止，仅有少数学者针对三维空间建模，系统地研究阻塞比的影响.基于Nakagawa[2]的方形截面柱风洞试验，Kim[10]在三维空间建立流场，研究在约束条件下的流场特性，并与试验结果吻合较好，但该文献仅给出了阻塞比为20%一种阻塞比工况.

方柱作为研究钝体绕流的常用模型，在工程中也是广泛采用的截面形式.针对二维方柱绕流，无论是风洞试验还是数值模拟都有很多已发表的文献，尤其是在雷诺数Re=2.2×104条件下[11-16].本文选取同样的截面形式及相近的雷诺数，采用三维大涡模拟(LES)，研究阻塞比(柱体迎风面面积与计算域截面面积之比)对方柱绕流的影响.在下文中，首先介绍了数值计算采用的湍流模型、边界条件、求解器设置以及网格划分等参数；然后通过与以往试验及数值模拟的结果进行对比，验证本文数值模拟计算结果的合理性；最后详细分析了阻塞比对方柱绕流气动力及流场的影响.

1 计算参数设置

1.1 湍流模型

本文采用大涡模拟进行数值模拟计算.通过对流体连续性方程及动量方程空间上的滤波，得到过滤后的Navier-Stocks控制方程，表达式如下：

(1)

(2)

(3)

式中：ui(i=1,2,3)是指3个方向的风速分量；p表示的压力；v为运动粘性系数；带有上划线的量表示空间滤波后的流场变量；τij代表亚格子应力，在动态Smagorinsky亚格子应力模型中[17]，利用线性涡粘模型将亚格子应力的各项异性部分与最小求解尺度上的应变率联系起来，如式(4)所示.

(4)

1.2 计算域及边界条件

为了研究不同的阻塞比对柱体气动力以及风场特性的影响，设置了5种阻塞比工况，不同工况之间只是计算域的宽度发生变化，计算域的宽度从24D减小到4D(D为模型迎风面宽度)，对应的阻塞比由4.2%增大25%，各工况的阻塞比以及对应的计算域尺寸设置如表 1所示.图 1所示为计算域尺寸及边界条件，上游入口到柱体迎风面距离为6D，柱体背风面到下游出口的距离为16D，计算域顺流向长度总计为23D.以往的数值模拟，多数采用了上游4.5D～8D以及下游15D～16.5D的计算域长度[13-14,18-19]，本文的长度介于上述范围.为了形成三维湍流结构，展向长度取4D[13,20].计算域宽度B的取值决定了阻塞比的大小，见表1.采用图 1所示的网格区域划分方法：在计算域横向上划分了多个不同的区域，其中I区为柱体附近区域，其宽度为2D，在此区域内网格比较密集，用于提高柱体附近区域涡旋捕捉的精度；II区为靠近约束边界的区域，宽度为0.5D，在计算域的两侧设置了无滑移边界条件，网格局部加密；通过改变III区的宽度改变计算域的宽度，得到不同的阻塞比.

表1 各阻塞比工况计算域尺寸

a 侧视图

b 俯视图

对于5种阻塞比工况， I区和II区两个区域内的网格划分完全一致，以提高计算结果的可对比性.I区是柱体所在区域，在柱体表面周向均匀布置 200个网格.对于LES数值模拟，需要保证障碍物附近区域有足够精度的网格分辨率，在柱体近壁面第1层网格高度约为D/2000，要小于0.1D/Re0.5=6.7×10-4D，向外延伸的增长率为1.13，计算得到的所有工况柱体近壁面Y+<2(Y+=uτy/v，uτ为壁面摩擦速度；y为近壁面第1层网格高度).柱体展向网格长度为0.1D.网格划分全部采用结构化，图 2为工况1的整体网格划分以及柱体附近区域网格局部放大示意图.由于精确控制每一个分区网格的尺寸，在分区分界线两侧的网格尺寸基本一致，避免相邻网格尺寸突变过大而导致的计算误差.

图2 工况1计算域网格划分示意图

图1也给出了计算域的各边界条件.入口边界条件设置为速度入口(velocity-inlet)，u=(U0,0,0)，湍流强度为0，雷诺数Re=2.2×104(根据入口平均风速以及柱体宽度尺寸计算得到)；柱体表面以及计算域的两侧边界均采用无滑移边界条件(Wall)；展向的边界设置为周期边界条件(Periodic)；出口采用自由出流边界条件(Outflow).

数值计算采用有限体积法，采用SIMPLEC格式求解压力速度耦合方程组，对流项采用数值耗散低的二阶中心差分格式(bounded central differencing)，时间离散项都采用二阶全隐式方法(bounded second-order implicit).为了保证计算最大柯朗数(Courant number)小于1，量纲为一计算时间步长Δt*为0.003 2(Δt*=Δt·U0/D，Δt为实际计算时间步，U0为来流风速).为加快计算的收敛速度，在进行LES计算之前，先采用雷诺平均定常计算(RANS)方法对流场进行初始化.在LES计算稳定之后开始采样，采样时长均大于40个旋涡脱落周期，所有计算工况均采用以上求解设置.

2 结果分析

2.1 计算结果验证

在进行不同阻塞比工况计算之前，针对工况1进行了两组不同数量网格的计算，如表 2所示.同时在表 2中给出了各气动力参数的计算结果.可以看出，对网格进行加密之后，各气动力参数并没有明显的变化.因此本文以下工况采用网格2的计算结果，并且其他阻塞比工况在I区和III区域工况1保持一致的网格.

表2 工况1网格收敛性验证

为了验证本文数值模拟的准确性，将本文工况1工况计算结果与以往的试验[11,21-22]结果以及数值模拟[13-16,23-24]结果进行对比.与本文不同的是，这些数值模拟中，两个侧面的边界均不是无滑移边界条件.表 3列举了风洞试验、LES模拟以及直接模拟(DNS)的结果，对比项包括：平均阻力系数(CD,mean)、脉动阻力系数(CD,rms)、脉动升力系数(CL,rms)、斯托罗哈数(St)、回流区长度(Lf)等.由表 3可见，本文的平均阻力系数、脉动阻力系数相比前人的研究结果略微偏大，而斯托罗哈数和回流区长度与各数值模拟以及试验结果都比较接近.表 3中也给出了文献中相应的模型长宽比(Ar)以及来流的湍流强度(Iu).

图3给出了工况1沿柱体中心截面环向的风压系数分布图，并将DNS[15]和LES[16,24]模拟结果以及风洞试验结果[21-22]一并示于图中.图 3a为平均风压系数分布，可以看出在迎风面，本文结果与风洞试验结果以及前人数值模拟结果吻合较好；在侧面与背风面，本文结果与风洞试验结果略有不同，但与Trias[15]和Cao[16]的数值模拟结果吻合较好.

图3b为脉动风压系数的分布，在迎风面以及背风面，本文计算结果与试验结果、数值模拟结果吻合较好；在侧面，本文计算结果要略大于试验结果，而与其他数值模拟结果吻合较好.

由表3及图3可见，不同的风洞试验及数值模拟结果存在明显的差异，这主要是由不同风洞实验室的试验条件、模型尺寸、来流湍流度等以及数值模拟的湍流模型、边界条件、网格划分等因素引起.而采样时间不足也会引起结果的偏差[15].但整体上，本文的数值模拟结果与前人试验结果及数值模拟结果均较为接近，是可靠有效的.

表3 整体气动力特性参数对比

a 平均风压系数

b 脉动风压系数

2.2 阻塞效应分析

2.2.1阻塞比对柱体气动力系数的影响

当阻塞比较大时，在方柱前缘边角分离的气流以及尾流会受到边界的约束作用，从而对柱体的气动力产生影响.柱体气动力系数的统计值随阻塞比的变化规律如图 4所示.其中，图 4a给出了平均阻力系数(CD,mean)、迎风面风力系数均值(Cf,mean)以及背风面风力系数均值(Cr,mean)随阻塞比的变化规律.从中可以看出，随着阻塞比R由4.2%增大到25%，平均阻力系数由2.35增加到3.0，增加了26.7%，同时迎风面平均风力系数略有减小，背风面风力系数(主要是吸力)均值绝对值则单调增加.这主要是由于阻塞比增大，气流流通面积被压缩，柱体尾流的流速增加，使得柱体背风面的吸力增加.图 4a还给出了Kim[10]在阻塞比为20%时方柱绕流的LES计算结果，其平均阻力系数与本文工况4(R=20%)的结果完全一致.

图 4b中给出了脉动阻力系数(CD,rms)、背风面风力系数脉动值(Cr,rms)和脉动升力系数(CL,rms)随阻塞比的变化规律.可以看出，脉动阻力系数值与背风面风力系数脉动值基本重合，这表明，阻力的脉动特性主要是来自背风面吸力的脉动特性的贡献.这种脉动特性会随着阻塞比的增大略微有所增加.脉动升力系数在阻塞比由4.2%增大到25%的过程中，由1.56逐渐增加到2.24，增加了43.6%.在图4b给出了Kim[10]在阻塞比为20%条件下方柱绕流的LES数值模拟结果，其脉动升力系数也与本文工况4(R=20%)的结果基本一致.

a 气动力均值

b 气动力脉动值

2.2.2阻塞比对Strouhal数的影响

柱体Strouhal数随阻塞比的变化规律如图 5所示，当阻塞比从4.2%增大到25%，Strouhal数从0.133增加到了0.175，增大了31%.这是由于洞壁对流场的约束作用加强，钝体绕流的流速增加所导致的.对于阻塞比较小的两组工况(R=4.2%和R=10%)，Strouhal数值相差很小，相差3%；同样平均阻力系数以及脉动升力系数值也很接近，相差在3%以内，如图 4所示.

图5 阻塞比对Strouhal数的影响

图6为各柱体阻力系数功率谱曲线.可以看出，对于阻塞比较小的前两个工况1、工况2(R=4.2%和10%)，阻力系数的功率谱曲线的值偏小，没有出现峰值，这时的阻力系数的脉动比较低.当阻塞比增大到一定值(R≥16.7%)，阻力系数的频谱曲线会在大约2倍Strouhal数附近出现明显的峰值，尤其工况5(R=25%)，其阻力系数频谱曲线的峰值非常突出.黄剑等[25]对不同阻塞比的三维矩形柱进行风洞试验时，也发现当阻塞比过大时，阻力系数频谱曲线会在大约2倍Strouhal数处出现峰值.

图6 阻塞比对阻力系数功率谱的影响

2.2.3阻塞比对风压系数的影响

图 7为各工况的柱体中心横截面(z=0)的平均风压系数以及脉动风压系数分布.如图 7所示，在柱体迎风面，随着阻塞比的增大，除靠近边缘区域的风压系数均值略微减小而脉动略微增大外，内部区域的风压系数几乎没有变化.在两侧面，随着阻塞比的增大，负风压系数的均值绝对值以及脉动值均显著增加(见图 7a、7b).同时可以看出，柱体两侧面，负风压系数均值绝对值和脉动风压系数都在靠近后缘附近出现极大值.极大值的位置对应柱体侧面附近回流区涡心的位置.随着阻塞比的增大，该极大值愈加明显，回流区涡心附近区域的流速在增加，且涡心位置逐渐向上游移动.

a 平均风压系数

b 脉动风压系数

在紧邻柱体背风面下游形成一对反向旋转的回流涡，由此在背风面形成稳定的负压.随着阻塞比的增大，负风压系数均值绝对值和脉动值，均有较小的增加.而阻塞比R≤10%的两个工况，洞壁对流场的约束不强烈，其平均风压系数及脉动风压系数的分布均比较接近.

图 8为各工况沿计算域中轴线分布的时间平均风压系数.随着阻塞比的增大，尾流区的负风压系数均值会显著减小.在x/D=1附近风压系数极小值位置处，当阻塞比R由4.2%增大到25%时，极小值由-2.23减小到-2.97，出口处的风压系数均值由-0.33减小到-0.89.同时可以看出，对于阻塞比较小(R≤10%)的两个工况，其沿轴线的平均风压系数分布比较接近，阻塞比的影响比较小.

图8 流场中轴线平均风压系数分布随阻塞比的变化规律

2.2.4阻塞比对柱体附近流场的影响

图9a、9b所示为顺流向平均速度分布，随着横向y的增加，顺风向速度由零先减小，在y/D=0.6附近取得极小值之后迅速增加，在y/D=0.8附近增加到极大值，之后再逐渐减小，直到在边界附近迅速减小到零(图中未画出).随着阻塞比的增大，极大值的位置会逐渐向远离柱体方向移动，且数值有明显的增加.相比于上游(x/D=-0.125)，下游位置(x/D=0.125)的顺流向平均速度的极大值出现的位置要更远离柱体.通常定义边界层厚度为壁面附近速度为远场速度的99%的位置[26]，由此可见随着阻塞比的增大，柱体壁面边界层的厚度略有减小.

顺流向雷诺正应力沿横向分布如图 9c、9d所示.对于上游位置处，随着远离柱体，雷诺正应力先增加后减小.而对于下游位置，雷诺正应力沿横向分布的规律要更为复杂一些，要经过两次先增加后减小的过程，最后逐渐减小到接近于零.随着阻塞比的增大，雷诺正应力也在逐渐变大.对于上游位置处，在y/D<0.7范围内，阻塞效应较为明显，在y/D>0.7范围内，不同阻塞比工况的雷诺正应力分布比较接近；对于下游位置处，这个分界点为y/D=0.8.

雷诺剪切应力沿横向分布如图 9e、9f所示.在上游的y/D<0.8以及下游的y/D<0.9区域内，随着阻塞比的增大，雷诺剪切应力逐渐减小；在上游的y/D>0.8以及下游的y/D>0.9区域内，随着阻塞比的增大，雷诺剪切应力有逐渐增加，随着与柱体距离的增加，阻塞效应会逐渐减弱.另外当阻塞比R≤10%时，流场特性比较接近，阻塞比对流场特性的影响较小.

a x/D=-0.125b x/D=0.125

c x/D=-0.125d x/D=0.125

e x/D=-0.125f x/D=0.125

2.2.5阻塞比对流场瞬时涡量的影响

图10给出了4种不同阻塞比(10%、16.7%、20%、25%)工况的瞬时涡量云图.图中所示平面为流场中轴面，瞬时涡量采用量纲为一的形式(ω*=ωD/U0).以逆时针旋转为正，实线代表的是正向涡旋，虚线代表的是反向涡旋.图中所有量纲为一涡量等值线最大值、最小值取值分别为±10，间隔为1.

对比4组不同阻塞比工况，可以看出随着阻塞比的增大，尾流流场涡旋的复杂性在增强、脱落的旋涡间距减小.如图 10a所示，当R≤10%时，洞壁对柱体绕流的约束作用较小，在柱体后缘脱落的旋涡向尾流耗散；由于洞壁给定的是无滑移边界条件，因此在上、下边界表面会形成一定厚度的边界层.当阻塞比增大到R=16.7%时，洞壁对流场的约束效应已经比较明显，柱体背风面附近的涡旋也更加复杂，尾流涡旋的耗散进一步减弱，在x/D=12处依旧可以看见比较明显的大尺寸旋涡；值得注意的是，此时，上、下洞壁表面不仅有边界层的出现，同时在此边界层内部也会出现一些小的涡旋，并向下游脱落.当阻塞比进一步增大到20%和25%时，尾流涡旋的耗散进一步减弱，同时在洞壁边界层内形成的涡旋的尺寸也在随之增大，并与柱体脱落的旋涡相互掺混.尤其是R=25%时，在上侧洞壁表面x/D=5以及x/D=13位置处有明显的大尺寸涡旋.Kim[10]通过被动粒子追踪模拟，也观察到沿上、下表面释放的粒子会形成二次旋涡，并向下游移动，这种现象会周期性地在上、下壁面交替出现.Okajima[22]在低雷诺数条件下，进行不同阻塞比的二维方柱的数值模拟也发现了相同的现象.

a 工况2(R=10%)

b 工况3(R=16.7%)

c 工况4(R=20%)

d 工况5(R=25%)

3 结论

对4.2%、10%、16.7%、20%和25%这5种不同阻塞比条件下的方柱绕流进行了CFD大涡模拟，将计算域横向边界设置为“wall”边界条件，用于模拟洞壁对流场的约束作用.通过与以往风洞试验及数值模拟结果的对比，验证了本文计算结果的合理性.在此基础上，分析了阻塞比对柱体气动力、表面风压分布、柱体附近流场特性以及尾流瞬时涡量的影响，得到了以下结论：

(1) 当阻塞比4.2%逐渐增大到25%时，方柱的平均阻力系数增大了27.6%，脉动升力系数增大了43.6%，Strouhal数增大了31%，且阻力系数频谱曲线会在2倍Strouhal数处逐渐出现明显的峰值.

(2) 随着阻塞比的增大，方柱迎风面的风压系数几乎不变，而背风面和侧面的负风压系数均值绝对值以及脉动值均明显增加；尾流中轴线的负风压系数均值绝对值随着阻塞比的增大明显增大.

(3) 随着阻塞比的增大，在分离剪切层内部、紧邻柱体表面区域的气流回流速度以及分离剪切层外部的顺流向流速均会增加；阻塞比对雷诺应力的影响主要表现在分离剪切层内部区域，并且下游受到的影响明显要强于上游位置处.

(4) 当阻塞比增大到一定程度(本文为R=16.7%)后，洞壁对流场约束作用的增强，在尾流区洞壁上、下表面附近会出现二次旋涡；随着阻塞比的增大，该涡旋的尺寸增加，与方柱上脱落的涡旋掺混现象也会越加明显.

本文重点考察当阻塞比在较大范围内变化时，柱体气动力及流场特性的变化.为了明确给出允许范围内的有价值的阻塞度建议值以及具有普适性的修正方法，需要更多的计算模拟，同时必须进行大量的模型试验，进行更深入地研究.