杜晓庆，顾李敏，吴葛菲，王玉梁，孙雅慧

(1.上海大学 土木工程系,上海 200444；2.上海大学 风工程和气动控制研究中心,上海 200444)

在实际工程中，圆柱群结构得到广泛应用，如桥梁缆索、输电导线、冷却塔群等[1-2].错列双圆柱是最常见的多圆柱结构布置形式，受到上游圆柱的尾流作用，下游圆柱易发生尾流激振现象，且尾流激振常常发生在小风向角条件下[3-5].尾流激振现象与圆柱的脉动升力联系紧密，但以往针对脉动升力特性的研究很少，其流场机理尚未澄清[6].

受试验条件的限制，以往研究大多集中于平均气动力，对脉动升力的研究相对较少[6].Arie等[7]在Re=1.51×105下，通过风洞试验研究了不同间距下串列双圆柱的脉动气动力，但并未分析其流场机理；Alam等[8-10]通过风洞试验在Re=6.5×104下，研究了错列双圆柱的脉动升力系数随间距布置的变化规律，但对小风向角下的双圆柱的脉动升力的研究仍不够系统，仅研究了0°、10°、25°三个风向角下的脉动升力变化规律.

同时，以往研究者基于流迹显示试验对双圆柱绕流场进行流态分类[6].Zdravkovich[11]将错列双圆柱划分为尾流干扰区、邻近干扰区、邻近和尾流干扰区及无干扰区.Alam等[10]、Gu等[12]、Sumner[13]对不同间距布置下的双圆柱的流场结构做了进一步的细分.但上述学者均是基于低雷诺数下的流迹显示试验对双圆柱的流态进行分析.如Gu等[12]的流迹显示试验是在Re=5.6×103下进行的，Sumner[13]和Alam等[10]分别在Re=850～1 900及Re=350下进行的.但双圆柱绕流有强烈的雷诺数效应[15]，在高雷诺数下，双圆柱的流场特性及干扰机理与低雷诺数下存在明显差异[15].同时在强/台风下，实际工程中的双圆柱结构常工作于Re>105的环境中，因而有必要在高雷诺数下研究双圆柱绕流问题.

本文使用商业软件Fluent里的大涡模拟方法，在高亚临界雷诺数下(Re=1.4×105)，研究了圆心间距P/D=1.5～4、风向角β=0°～90°的错列双圆柱的脉动升力系数，重点研究了β=0°～30°范围内双圆柱表面脉动风压分布特性，结合功率谱、瞬时涡量以及脉动风压场，探讨了上、下游圆柱之间脉动气动力和流场之间的内在联系，讨论了上、下游圆柱之间相互干扰的流场机理.

1 计算模型与参数

1.1 控制方程与亚格子模型

在湍流运动中，来流的运动可以看成由多种尺度的涡叠加形成.为模拟湍流运动，许多学者采用不同的数值模拟方法来进行研究.大涡模拟(LES)方法放弃对全部尺度的旋涡的求解，而是直接通过Navier-Stokes方程求解湍流中的大尺度旋涡，对小尺度的旋涡的求解则采用亚格子尺度模型(SGS)求解.该方法可节约计算资源，提高计算能力，同时也可收集到丰富的流场信息.因此，本文采用大涡模拟方法对错列双圆柱周围的流场进行数值模拟.

通过滤波函数的过滤，得到关于大尺度涡的不可压缩的Navier-Stokes运动方程：

(1)

(2)

本文的亚格子应力τij采用Fluent中的Smagorinsky-Lilly的亚格子尺度模型：

(3)

(4)

式中μt代表亚格子尺度的湍动黏度，且

(5)

(6)

(7)

式中：Δ表示过滤的网格尺度，Cs为Smagorinsky常数，本文取0.1.

1.2 计算模型和计算参数

图1为本文计算所采用的双圆柱的布置及工况示意图.两个圆柱的直径均为D，间距为P，间距比P/D分别为1.5、2、3和4，风向角β范围为0°～90°.由圆柱直径D以及来流风速U0计算得到的雷诺数为1.4×105.

本文采用O型计算域，计算域的直径为46D， 阻塞率在2.2%～4.3%范围内，模型展向长度为2D.本次计算采用速度入口边界条件及自由出口边界条件.圆柱展向两端采用周期性边界条件，圆柱壁面为无滑移边界条件.计算采用的求解器为压力分离求解器.计算模型采用结构化网格，网格总数约为270万～320万.圆柱近壁面最小网格尺度为0.000 1D，经换算得到壁面的y+≈1.无量纲时间步Δt*=0.005(Δt*=ΔtU0/D，其中Δt为实际计算时间步，U0为来流风速).

圆柱表面的风压系数定义为

(8)

式中：p为当地风压，p0为远前方风压，ρ代表来流空气密度.在下文中，C′p均代表脉动风压系数.

阻力系数CD和升力系数CL的定义为：

(9)

(10)

式中：FD和FL分别代表单位长圆柱表面的升力与阻力.具体方向见图1(a)，图中小标“1”和“2”分别表示上游圆柱和下游圆柱.

图1 计算模型示意

1.3 结果验证

对于单圆柱的流场绕流，作者对周向网格数、展向长度、无量纲时间步长等相关因素进行了研究，并与文献结果进行了详细的比较，具体的网格独立性检验和结果验证见文献[17-18].图2给出了本文单圆柱的平均和脉动风压系数与文献中风洞试验结果的比较，结果显示本文单圆柱的结果与文献结果吻合良好.图3、4分别为本文串列双圆柱的平均阻力系数和脉动升力系数与文献结果的比较.从总体上看，本文结果与文献值的变化趋势相同，数值接近，吻合较好.值得指出的是，从图4可看出，不同文献结果在P/D=3～4之间的离散性较大，这主要是因为在此间距比内，串列双圆柱的绕流场会从剪切层再附流态转变为双涡脱流态[14]，从而导致脉动升力系数发生突变.而引起脉动升力系数突变的间距比对雷诺数、来流湍流度和圆柱表面粗糙度等参数非常敏感.因而不同文献的脉动升力系数的差异很可能是因为不同研究者的研究条件的差异造成的.

图2 单圆柱风压系数分布[18]

图3 串列双圆柱平均阻力系数随间距比的变化

图4 串列双圆柱脉动升力系数随间距比的变化

2 结果分析

2.1 脉动升力系数

图5给出了上、下游圆柱在不同间距比下，脉动升力系数随风向角的变化曲线.由图5(a)可知，当P/D=1.5～2时，在全风向角下，上游圆柱的脉动升力系数均远小于单圆柱；在β=5°附近时，上游圆柱脉动升力达到极小值，并随着风向角的增大而增大，且在β=90°时达到极大值.结合下文的流场分析可得，此时下游圆柱的存在抑制了上游圆柱的旋涡脱落，导致上游圆柱的脉动升力系数始终小于单圆柱.当P/D=3～4时，上游圆柱的脉动升力随风向角的变化较为剧烈；在P/D=3时，上游圆柱脉动升力系数仅在β=10°～20°及90°时大于单圆柱，同样在β=5°时取得极小值；在P/D=4，β<45°时，上游圆柱的脉动升力系数均大于单圆柱.

由图5(b)可知，下游圆柱脉动升力系数随风向角的变化比上游圆柱更为剧烈.当P/D=1.5～2时，在β=5°～20°范围内，下游圆柱的脉动升力系数仍远小于单圆柱；当P/D=1.5时，下游圆柱的脉动升力系数随风向角增大存在先增大后减小的趋势，而当P/D=2时，脉动升力系数则始终呈增大的趋势，这可能与P/D=1.5时，双圆柱间存在的偏向流有关.当P/D=3时，下游圆柱的脉动升力仅在β=0°～5°时小于单圆柱，但当β=10°时，下游圆柱脉动升力瞬间增大，说明在该风向角下存在流态转变；当P/D=4时，下游圆柱的脉动升力系数在β=0°时取得极大值，且仅在β=90°时略小于单圆柱，表明随着间距比和风向角的增大，下游圆柱受上游圆柱尾流干扰逐渐减小.

图5 脉动升力系数随风向角的变化

为了更直观地展示上、下游圆柱脉动升力系数随间距比和风向角的变化趋势，图6给出了上、下游圆柱脉动升力系数的等值线图，图例中的黑点代表本文单圆柱的脉动升力系数.由图6可知，在小间距小风向角下，上、下游圆柱的脉动升力系数均远小于单圆柱；而在较大的间距比下，在β=0°～30°范围内，脉动升力系数出现极大值.该图与文献[6]结果吻合较好.针对小间距和较大间距，以下分别以P/D=1.5和3为例分析其在β=0°～30°范围内的气动力特性.

图6 脉动升力系数等值线图

2.2 脉动风压系数

图7、8分别给出了间距比P/D=1.5及3两种间距典型风向角下上、下游圆柱表面的脉动风压系数分布，还给出了单圆柱的脉动风压系数分布作为对比.

由图7可知，当P/D=1.5时，上游圆柱的脉动风压系数远小于单圆柱；与下游圆柱相比，上游圆柱上、下侧表面的脉动风压呈现较好的对称性.下游圆柱的脉动风压随风向角的变化有明显差异.当β=0°时，下游圆柱的脉动风压系数基本呈对称分布且其值远小于单圆柱的值.在β=10°的情况下，下游圆柱的脉动风压系数呈明显的不对称分布，在其迎风面下侧(θ=270°～360°)出现明显的峰值且远大于单圆柱的值.随着风向角的继续增大，下游圆柱的脉动风压系数逐渐接近单圆柱的分布.

当P/D=3时，由图8可看出，在β=5°的情况下，上游圆柱的脉动风压远小于单圆柱.而在其他风向角下，上游圆柱的脉动风压分布与单圆柱基本一致，上游圆柱受下游圆柱的干扰较小.下游圆柱的脉动风压系数分布随风向角发生显著的变化.当β=5°时，下游圆柱的脉动风压系数分布与P/D=1.5，β=10°相近，在下游圆柱迎风面下侧(θ=270°～360°)均出现了脉动风压极值；与P/D=1.5，β=10°不同的是，此时下游圆柱上表面(θ=0°～180°)的脉动风压系数明显增加.在β=10°的情况下，下游圆柱的脉动风压系数整体大于单圆柱的值，且呈现出明显的不对称性，并在其迎风面下侧(θ=270°～360°)处出现峰值.当风向角β增大到20°时，在下游圆柱的下表面(θ=180°～360°)出现两个峰值，下文结合流场进一步分析.

2.3 升力系数功率谱

为了进一步解释上、下游圆柱在β=0°～30°下不同间距的脉动气动力系数变化原因及相互干扰的流场机理，图9、10分别给出了β=0°～30°下的上、下游圆柱的升力系数功率谱.由图9可知，在P/D=1.5的情况下，上游圆柱升力系数功率谱峰值始终较小，此时下游圆柱的存在抑制了上游圆柱的旋涡脱落.对于下游圆柱而言，在风向角β=0°～20°时，其功率谱的峰值较弱，这说明此时下游圆柱旋涡脱落同样受到抑制；当β=20°～30°时，下游圆柱的升力功率谱中St数峰值增加，下游圆柱旋涡脱落强度有所增强，但其涡脱频率仍比单圆柱的要小，此时下游圆柱仍受到较大干扰.

图9 P/D=1.5时升力系数功率谱随风向角的变化

图10 P/D=3时升力系数功率谱随风向角的变化

当P/D=3时，由图10可得，上、下游圆柱的升力系数功率谱与P/D=1.5时有较大差别.在β=0°及5°的情况下，上游圆柱的升力系数功率谱没有出现明显的峰值.随着风向角的进一步增大，上游圆柱St数峰值明显增加，这表明随着风向角的增大，上游圆柱受下游圆柱的干扰减弱，其流态逐渐向单圆柱靠近.反观下游圆柱，其升力系数的功率谱波动更为剧烈.在P/D=3的情况下，除β=5°以外，下游圆柱均存在明显的St数，这表明下游圆柱自身一直存在较强的旋涡脱落.值得指出的是，当P/D=3，β=20°时，下游圆柱的功率谱中出现两个St数，结合下文流场特性进一步分析可知：首先，对于下游圆柱而言，其自身会有一个旋涡脱落频率，对应的St数为0.2；其次，从上游圆柱脱落的旋涡会与下游圆柱产生相互作用，并导致下游圆柱产生另一个涡脱频率，对应的St数为0.17.

2.4 干扰流态分类

图11、12分别给出了β=0°～30°内，上、下游圆柱在各个间距下的典型时刻涡量图(单位S-1)以及脉动风压场，图11还给出了对应涡量时刻的表面风压.依据上、下游圆柱瞬时流场特性以及上文的气动力特性，并结合参考文献[17-18]对平均流场的分析，将上、下游圆柱间的干扰流态主要分为5种干扰流态：尾流干扰、剪切层干扰、邻近干扰、旋涡撞击以及涡的相互作用.

1)尾流干扰流态.当P/D=1.5～3且β=0°时，上、下游圆柱的间距较近，上游圆柱的剪切层会再附到下游圆柱表面，使得下游圆柱完全沉浸在上游圆柱的尾流中，这导致下游圆柱涡脱受到抑制，脉动升力及脉动风压远小于单圆柱.对于上游圆柱而言，下游圆柱的存在也完全抑制了上游圆柱的涡脱，使得上游圆柱的脉动气动力远小于单圆柱.

图11 错列双圆柱的瞬态涡量和表面风压分布

图12 错列双圆柱脉动风压场

2)剪切层干扰流态.当P/D=1.5～2且β=5°～30°、P/D=3且β=5°时，上游圆柱上侧剪切层会紧贴下游圆柱的迎风面并向下运动，随后与下游圆柱下侧剪切层相互作用，使得下游圆柱迎风面下侧出现脉动风压极值.同时，该相互作用导致下游圆柱下侧剪切层没有发生卷起，进而减弱了下游圆柱上、下侧剪切层之间的相互影响，使得下游圆柱脉动升力较小.此外，由于下游圆柱的影响，上游圆柱的上、下侧剪切层并未在圆柱背风面交替卷起形成旋涡，导致其脉动风压以及脉动升力远小于单圆柱.

3)邻近干扰流态.当P/D=3～4且β=20°～30°时，上、下游圆柱背风面均有旋涡脱落，此时上、下游圆柱之间的干扰较弱，两个圆柱的流态都与单圆柱接近，脉动升力系数接近单圆柱.

4)旋涡撞击流态.当P/D=4且β=0°～15°、P/D=3且β=10°～15°时，上游圆柱上下侧剪切层交替卷起形成旋涡，旋涡会撞击在下游圆柱迎风面，造成撞击处的风压产生剧烈的波动，这种旋涡的不断撞击是下游圆柱表面出现脉动风压的极值以及脉动升力极值的主要原因，下游圆柱极易出现尾流激振现象.同时，在这种流态下，上游圆柱基本不受下游圆柱的干扰，其流态与单圆柱的相接近.

5)涡的相互作用流态.当P/D=3～4且β=20°时，上游圆柱上侧剪切层卷起形成的旋涡不再撞击到下游圆柱的表面上，而是与下游圆柱下侧(间隙侧)的旋涡耦合产生相互作用，造成下游圆柱下表面的脉动风压产生极值；同时下游圆柱背风侧的旋涡脱落会造成另一个脉动风压极值的出现.这种上下游圆柱之间的涡的相互作用造成下游圆柱出现较大脉动升力系数，使得下游圆柱极易发生尾流激振.反观上游圆柱，此时的上游圆柱基本不受下游圆柱的干扰，流态与单圆柱的相接近.

值得指出的是，不同的干扰流态下双圆柱的气动性能有很大差异，因此在多圆柱结构设计风荷载取值时需要注意这一问题.此外，当双圆柱发生“涡的相互作用流态”和“旋涡撞击流态”时，下游圆柱的脉动升力较单圆柱有明显增大，从而导致下游圆柱易发生尾流激振，这也是多圆柱结构在设计时需要关注的另一个重要问题.

3 结 论

在高雷诺数(Re=1.4×105)下，本文对错列双圆柱绕流问题进行大涡模拟研究，分析了双圆柱的脉动升力随间距和风向角的变化规律，着重探讨了β=0°～30°范围内错列双圆柱的脉动气动力与流场之间关系，并讨论了相关的流场干扰机理.主要结论如下：

1)上、下游圆柱的脉动升力系数随双圆柱的间距布置变化剧烈.当P/D=1.5～2时，上游圆柱的脉动升力系数在全风向角范围内始终远小于单圆柱，而下游圆柱的脉动升力系数则分别在β=0°～30°及β=70°～90°范围内出现极小值；当P/D=3～4时，上、下游圆柱的脉动升力系数明显增大，随着风向角的增大，上、下游圆柱脉动升力系数均呈先增大后减小的趋势，但下游圆柱脉动升力系数大于上游圆柱.

2)在β=0°～30°范围内，随着风向角及间距比的变化，根据绕流场的脉动风压及典型时刻的涡量图，将双圆柱间相互干扰的流态分为5种：尾流干扰、剪切层干扰、邻近干扰、旋涡撞击及涡的相互作用流态.

3)在尾流干扰流态下，下游圆柱的存在抑制了上游圆柱的旋涡脱落，而下游圆柱处于上游圆柱尾流中，涡脱强度较弱，故上、下游圆柱的脉动升力系数均小于单圆柱；在剪切层干扰流态下，上游圆柱上侧剪切层与下游圆柱下侧剪切层相互耦合，削弱了下游圆柱的旋涡脱落，导致下游圆柱脉动升力系数较小.在邻近干扰流态下，上、下游圆柱之间的干扰较弱，脉动升力系数接近单圆柱.

4)在旋涡撞击流态下，上游圆柱尾流中的涡撞击到下游圆柱迎风面，导致下游圆柱的脉动升力系数远大于单圆柱.而在涡的相互作用流态中，从上游圆柱脱落的旋涡会与下游圆柱的下侧旋涡相互耦合，使得下游圆柱的脉动升力明显增大.值得指出的是，在这两种干扰流态中，受到上游圆柱的干扰，下游圆柱的脉动升力较单圆柱有明显增大，易使得下游圆柱发生尾流激振，这也是多圆柱结构在设计时需要关注的另一个重要问题.