高Re数圆柱绕流二维RANS模拟适用性分析

祝志文

(湖南大学 风工程试验研究中心，长沙 410082)

圆柱绕流即为经典的流体力学问题，亦为流场配置最简单、最基本的钝体绕流。由于圆柱绕流随流动Re数的变化呈现极复杂的变化特征，因此吸引诸多研究者以期从流动试验或数值模拟上解释其流动形态与机理［1］。以往因受计算资源限制，圆柱体绕流数值模拟通常在较低的Re数下进行［2］。虽CFD模拟可用于桥梁抗风［3-4］，但斜拉桥拉索在高风速下Re数一般在105以上，属于高Re数范围，相关数值模拟报导较少。众所周知，圆柱绕流随流动Re数、来流条件及圆柱表面粗糙度的变化呈现极复杂的变化特征，尤其阻力系数呈现较显著的Re数效应，如图1所示［5］，其中Re＜2×105称为亚临界区，阻力系数的Re数效应并不明显。Re＞5×105称为超临界区，圆柱阻力系数随Re数的增大而缓慢增长。2×105＜Re＜5×105通常被定义为圆柱的临界Re数区。在斜拉索绝大部分工作状态风速下，其流动Re数可能均落在圆柱体的临界Re数区，其阻力系数随Re数的变化而快速变化，合理确定该阻力系数会有较大困难。

图1 圆柱体阻力系数随Re数变化曲线Fig.1 Drag coefficient of circular cylinder Vs.Re number

本文以直径D=0.12 m圆柱为研究对象，通过数值模拟获得圆柱气动力参数与漩涡脱落St数。基于CFD模拟的质量控制要求［6］，开展网格无关与计算时间步长无关检查。通过获得圆柱多个工作风速下的阻力系数与St数，与相关文献研究结果对比，评价二维RANS数值模型在获得高Re数圆柱绕流气动特性上的可行性，澄清CFD研究认识误区，为高Re数圆柱气动特性数值研究提供建议。

1 数值方法

1.1 控制方程

绕圆柱断面非定常二维不可压流动的雷诺时均N-S方程为:

其中:ρ，μ分别为空气密度与分子粘性;p，t分别为静压与时间;ui，u'i和uj，u'j(i=1，2;j=1，2)分别为气流沿坐标轴xi(i=1，2)的平均与脉动速度;为雷诺应力。

上述雷诺应力的引入使控制方程不封闭，需引入湍流模型以求解。若基于涡粘假设，可将雷诺应力表示为:

其中:μt=ρCμk2/ε为湍流粘性;Cμ为经验常数;k，ε分别为湍动能及耗散率，通过求解湍流模型方程确定。

标准k-ε模型会过高估计流动滞点区域的湍动能，不适用于风工程问题的数值模拟［6］。本文采用综合标准k-ε模型与k-ω湍流模型的SSTk-ε湍流模型。该模型结合上述两湍流模型优点，抛弃其缺点。即利用标准k-ε模型适合剪切层模拟而k-ω模型适合近壁区模拟优点，从而通过设定一混合函数，使k-ω模型能在边界层内靠近壁面使用，而边界层外使用kε模型求解。相关研究认为，对分离点附近边界层内的非平衡流动，k-ωSST的模拟结果［7］优于标准k-ε模型。

1.2 计算域与计算参数描述

采用如图2所示形状的计算域，以便能采用分区网格划分方案并形成绕圆柱的高质量正交四边形单元。采用大的计算域，图2中圆柱左侧Z1区为半圆形计算域，该外边界为计算域入口，入口至圆柱中心距离为26D。圆柱右侧计算域由Z2，Z3两个区组成，其中Z2采用与Z1类似的半圆形计算域。Z3区右侧边界为计算域出口，至圆柱中心距离为40D。因此上下边界至圆柱中心距离同样为26D，对应的堵塞度小于2%。上述边界至圆柱中心距离由试算获得，可避免在外边界施加的边界条件对计算域内数值计算影响。

Z1、Z2区均采用结构化贴体四边形网格，见图3。Z3区为非结构四边形网格。圆柱周向等分为140个网格。为进行网格无关性检查，最靠近圆柱表面第一个网格点至物面距离h0分别为0.000 02 m，0.000 01 m，0.0000 05 m，分别称为粗网格 G1、细网格 G2、最细网格G3。沿物面法向，采用1.06的网格高度生长率，以保证在物面等流动变量变化梯度大的位置获得高分辨率网格。Z2、Z3区公共边两侧网格尺度保持接近。对三套网格系统，Z3区网格划分完全相同，区别在于Z1、Z2区网格数量及网格分辨率。因此在圆柱周围及尾迹区的大范围内能获得较高质量的正交四边形网格。三套网格系统各自总单元数N见表1。

图2 计算域分区Fig.2 Partition in computational domain

图3 圆柱周围网格Fig.3 Grids arrangement around cylinder

表1 网格划分参数Tab.1 Grid key parameters

对计算域外边界，在入口处采用速度入口条件，水平向速度等于来流速度，垂直水平向速度等于零，来流湍流度取0.25%。在计算域出口采用流动出口条件，即沿出口边界的法线方向，速度梯度等于零。计算域的上、下边界均采用对称边界条件，即垂直于该边界方向，流动变量梯度为零。在静止圆柱表面，采用无滑移边界条件，速度与湍流度均为零。

图4 G1网格圆柱表面的Y+值分布Fig.4Y+value on circular cylinder surface of G1

对圆柱绕流的非定常计算，控制方程时间离散采用二阶隐式格式，空间离散采用二阶迎风格式。压力方程与动量方程解耦采用SIMPLEC算法及欠松弛迭代。在网格无关及时间步长无关检查中，来流风速设为14.6 m/s，对应的圆柱流动 Re数为1.2 ×105。计算0.000 05 s～0.000 5 s范围内多个时间步大小。对每一时间步数值计算，每一子步进行50次迭代，当动量方程与湍流方程的残差小于1×10-5时，认为此时间步迭代计算已收敛。气动力系数及其它参数的获取，可用足够多时间步进行数值计算，充分剔除初始计算影响，气动力系数时程及趋势基本稳定后开始记录。三套网格对应的圆柱表面最大Yplus(Y+max)(表1)，因此均满足SSTk-ω湍流模型对物面Yplus值要求。其中粗网格系统G1在圆柱表面的Yplus分布见图4。本文数值模拟均在CFD专用软件Fluent 3.2.26中实现。

1.3 网格与时间步长无关检查

风工程数值计算，尤其需重视计算模型的网格无关与计算时间步长无关检查，目的为获得与网格尺度、单元数量、及非定常计算时间步大小无关解，即确定数值模拟结果不再随网格尺度与时间步长大小而变化的计算参数。此为风工程数值模拟质量保证的重要步骤。通常要求至少用三套网格进行绕流模拟，并比较不同网格与时间步长下数值模拟结果［6］。不同网格尺度主要集中在物面网格分辨率与离开物面的网格生长率(一般要求不大于1.2)。

定义二维圆柱阻力系数CD、升力系数CL、St数分别为:

其中:FD，FL分别对应作用在圆柱截面上的阻力和升力，阻力顺来流流向为正，升力垂直来流方向向上为正;U为计算域入口风速;f为二维圆柱漩涡脱落频率。

图5、图6分别为G2网格系统在时间步长dt=0.000 1 s下获得的圆柱截面阻力与升力系数时程，均表现为非纯谐波，说明除主频率成分外，亦有高阶谐波成分。但升力系数谐波特性较阻力系数好。升力系数主频率为阻力系数的二分之一，也对应圆柱的漩涡脱落频率，由此可确定二维圆柱与漩涡脱落有关的St数。因阻力系数时程幅值波动较大，因此不同计算工况下的阻力系数平均值及根方差值可能会受所取时程位置及长度影响，反映涡脱频率特征的St数更具参考价值。图7、图8为三套不同计算网格分别在不同时间步长下获得的圆柱阻力系数平均值、根方差RMS及涡脱St数。

图5 阻力系数时程Fig.5 Drag coefficient time history

图6 升力系数时程Fig.6 Lift coefficient time history

图7 不同工况阻力系数平均值和RMSFig.7 Mean and RMS values of drag coefficients in different simulation cases

图8 不同工况St数Fig.8 Strouhal number in different simulation cases

由图7、图8知，在所有计算工况下，阻力系数RMS值均较小，而阻力系数平均值则呈现差别，随时间步的增大，气动参数差别亦增大。为有效分辨钝体绕流非定常特征，非定常计算时间步长应不小于涡脱周期的1/200～1/500。因此本文取时间步长不大于0.000 1 s的不同网格获得的气动参数结果。可见，三套网格计算所得阻力系数平均值及St数的差别很小，减小时间步长该值无明显变化。因此，当时间步长取不大于0.000 1 s时，三套网格的模拟结果无明显变化，可认为此时三套网格均能获得与网格与时间步长无关的解。基于阻力系数平均值与St数在时间步长大于0.000 2 s后的变化及所用机时，本文确定细网格G2为后续圆柱不同风速下的计算网格。

2 二维圆柱气动特性模拟结果

为获得二维圆柱在不同Re数下的气动特性，据以上数值模型网格无关与时间步无关检查结果，分别对该二维圆柱计算来流风速为10 m/s、15 m/s、20 m/s、25 m/s、30 m/s、35m/s、40m/s 和 45m/s 下的绕流场。考虑时间步长影响，前四个风速的时间步长为0.000 1 s，后四个风速的时间步长为0.000 05 s，所有计算工况下圆柱表面的Y+值均满足SSTk-ω湍流模型要求。图9为计算所得圆柱截面在多个来流风速下的阻力系数时程，可见经初始计算波动及平稳后阻力系数均逐渐呈规律振动。图10为30 m/s来流风速下圆柱周围流场压力等高线，从尾迹压力等高线可观察到圆柱尾迹的漩涡脱落情况。

图9 不同风速下的阻力系数时程Fig.9 Drag coefficients time histories under various wind speeds

图10 圆柱周围的压力等高线云图Fig.10 Pressure contours around cylinder

图11为圆柱绕流在T/4、T/2、3T/4及T四个时刻的尾迹漩涡脱落图画，其中T/4，3T/4对应升力系数时程分别达到负最大值与正最大值时刻，而T/2，T分别对应升力系数时程回到零点及离开零点时刻。可见数值模拟能显示圆柱尾迹漩涡产生、拉长、脱落及漂移过程。

图11 一个周期内圆柱尾迹漩涡脱落Fig.11 Vortex shedding in wake of circular cylinder during one cycle

数值模拟所得二维圆柱阻力系数平均值与涡脱St数随Re数的变化分别见图12、图13。在计算的Re数范围内，圆柱阻力系数平均值及涡脱St数与相关文献［9-10］结果变化趋势相同，即本文数值模型能在一定程度上反映气动参数的Re数效应。数量上，无论阻力系数平均值或涡脱St数，均与相关文献结果存在较大差异。Re＜3.55×105时，本文数值模拟所得阻力系数平均值小于风洞试验结果［10］，而 Re＞3.55 ×105时，本文数值模拟所得可能会大于风洞试验阻力系数平均值，表明本文CFD数值模型在预测圆柱定常气动特性上的不足。另外，从图13可见，本文CFD数值模型明显高估圆柱涡脱St数［9］，即圆柱实际涡脱周期被明显低估了，反映了基于RANS的CFD数值模型在预测圆柱非定常气动特性上的不足，可能与RANS对非定常流动的时均处理，使部分非定常特性被平均化，导致难以准确捕捉圆柱绕流的非定常特性。本文曾采用雷诺应力模型(RSM)代替SSTk-ω湍流模型进行圆柱绕流模拟，研究结果基本相同。篇幅有限，不再赘述。

图12 圆柱阻力系数平均值随Re数的变化Fig.12 Cable drag coefficients vs.Re number

图13 涡脱St数随Re数的变化Fig.13 Strouhal number vs.Re number

3 结论

基于CFD专用软件Fluent3.2.26与严格数值模拟方法及步骤，对二维RANS数值模型在预测圆柱定常与非定常气动特性上展开研究，结论如下:

(1)本文给出的圆柱计算域分区网格划分方法能在圆柱周围及尾迹区大范围内形成高质量正交四边形网格。由此开展的网格无关与时间步长无关检查是本文数值方法重要基础。

(2)二维RANS数值模型能一定程度给出圆柱气动参数的Re数效应，也能给出圆柱尾迹漩涡产生、生长、脱落及漂移的周期性过程。

(3)二维RANS数值模型给出的圆柱阻力系数平均值与涡脱St数，均与相关文献结果存在较大差异。尤其CFD模拟明显高估了圆柱的涡脱St数，反映出二维RANS模型在预测圆柱定常与非定常气动特性上的不足。

采用二维RANS模型难以合理预测圆柱的气动力及涡脱特性。若研究圆柱定常与非定常特性，需尝试采用三维模拟或其它相关方法。

［1］Nishimuraa H，Taniike Y.Aerodynamic characteristics of uctuating forces on a circular cylinder［J］.Journal of Wind Engineering and Industrial Aerodynamics，2001，89(7-8):713-723.

［2］Franke J，Frank W.Large eddy simulation of the flow past a circular cylinder at ReD=3 900［J］.Journal of Wind Engineering and Industrial Aerodynamics，2002，90(10):1191-1206.

［3］祝志文，陈政清，顾 明.桥梁气弹响应分析的时域气动模型［J］.中国公路学报，2007，20(6):43-48.

ZHU Zhi-wen，CHEN Zheng-qing，GU Ming.Time-domain aerodynamic model in aeroelastic response analysis of bridges，China Journal of Highway and Transport，2007，20(6):43-48.

［4］祝志文，陈政清.数值模拟桥梁断面的颤振导数和颤振临界风速［J］.中国公路学报，2004，15(4):41-50.

ZHU Zhi-wen，CHEN Zheng-qing.Numerical simulations for aerodynamic derivatives and critical flutter velocity of bridge deck［J］.China Journal of Highway and Transport，2004，15(4):41-50.

［5］Simiu E，Scanlan R H.Wind effects on structures-an introduction to wind engineering［M］.New York:Wiley，Third Edition，1996.

［6］Casey M，Wintergerste T.Best practice guidelines:ercoftac special interest group on“quality and trust in industrial CFD”［M］.Fluid Dynamics Laboratory，Sulzer Innotec，2000.

［7］ Franke J，Hirsch C，Jensen A G.，et al.Recommendations on the use of CFD in wind engineering［C］.Proceeding of International Conference on Urban Wind Engineering and Building Aerodynamics.Belgium，2004.

［8］Poulin S，Larsen A.Drag loading of circular cylinders inclined in the along-wind direction［J］.Journal of Wind Engineering and Industrial Aerodynamics，2007，95(9-11):1350-1363.

［9］Norberg C.Fluctuating lift on a circular cylinder:review and new measurements［J］.Journal of Fluids and Structures，2003，17(1):57-96.

［10］林志兴，杨立波，李文勃.斜拉索顺桥向风阻系数的试验研究［J］.郑州大学学报(工学版)，2005，26(1):38-41.

LIN Zhi-xing，YANG Li-bo，LI Wen-bo.Experimental study on drag coefficientsofstay-cablescorresponding towind direction along the bridge centrallin［J］.Journalof Zhengzhou University(Engineering Science)，2005，26(1):38-41.