楚晨晖,陈少林

(1.无锡环境科学与工程研究中心,无锡 214100;2.无锡城市职业技术学院 建筑与环境工程学院,无锡 214100;3.南京航空航天大学 土木与机场工程系,南京 211100)

圆柱绕流是流体力学及钝体空气动力学研究的经典问题。早期ACHENBACH[1]、NISHIMURA[2]、ROSHKO[3]、GROVE[4]等通过试验手段分别对不同雷诺数条件下的圆柱绕流进行了深入研究。随着CFD技术的发展和完善,数值模拟技术也被广泛的应用在圆柱绕流的理论研究和工程实践中[5-7]。由于椭圆柱截面具有偏流线型的特点,在流体中具有更好的对抗流体作用,且椭圆柱作为介于圆柱和平板之间的一种更具代表性的钝体,深入研究椭圆柱绕流问题,有助于发现钝体绕流具有的共性规律[8-10]。前人对于圆柱、椭圆柱等类圆柱钝体绕流都做了系统的研究工作。梭形柱长轴为迎风面时也属于钝体绕流范畴,但其具有自身的特点,如表面平缓、边缘尖锐,与圆柱体、椭圆柱体等一般钝体的几何特性明显不同。梭形柱的边缘角度以及表面曲率的不同会改变其表面附近的空气流线进而影响其气动特性。随着建筑形态不断多样化,诸如梭形等类圆柱体的建构筑物不断涌现。这些结构形态给人们在带来优美的审美享受的同时,也考验工程师对于结构抗风的设计经验。然而,对梭形柱的绕流问题相关研究却鲜有报到。因此,有必要对梭形柱体的表面气动力及流场特性进行深入研究,以为理论研究和工程实践提供参考和指导。

基于此,本文应用Fluent软件对超临界态下(Re=7.2×105)的梭形柱体进行CFD数值模拟研究。采用剪切应力运输(SST)模式进行梭形柱绕流数值模拟,并将所得计算参数与已有文献试验数据对比。结果表明采用SST模式的k-ω湍流模型能够很好地模拟超临界态下的钝体绕流。然后,保持湍流模型、网格尺寸以及计算参数不变,通过改变梭形柱长短轴比,分析研究不同尺寸下的梭形柱在超临界态的气动力及流场特性。研究结果可以为类圆形柱绕流的基础研究和工程实践提供参考依据。

1 湍流模型及气动参数计算公式

1.1 湍流模型选取

采用平均雷诺应力方法求解Navier-Stokes方程,需要引入湍流模型方程使其封闭并有解。剪切应力运输(SST)湍流模型是由F R MENTER[11]提出的一种基于k-ω两方程的湍流模型,被认为适合于对钝体结构分离流动的模拟,近年来得到了较多的应用[12-14]。因此选用SST湍流模型模拟超临界条件下的梭形柱绕流试验。

1.2 气动参数计算公式

为描述风压在柱体表面的分布特性,引入无量纲平均风压系数及脉动风压系数,其定义分别为

(1)

(2)

结构阻力系数、升力系数计算公式分别为

(3)

(4)

式中:Cd,Cl分别为结构的无量纲阻力系数、升力系数;Fd,Fl分别为计算获得的结构整体阻力、升力,N;A为结构迎风垂直投影面积,m2。

斯脱罗哈数是描述雷诺数与脱涡频率之间关系的无量纲系数,其计算公式为

(5)

式中:St为无量纲斯脱罗哈数;fs为脱涡频率,Hz;D为结构特征长度,m;U为来流平均速度,m/s。

2 CFD参数设置及数值模拟合理性验证

2.1 圆形截面及梭形截面几何尺寸

采用直径为30 cm、高度为50 cm的圆柱体作为验证模型。然后通过改变短轴尺寸,使圆柱体逐渐转变为梭形柱体。为更好地表示梭形柱体几何尺寸,参考椭圆扁率,定义梭形柱体扁率C:

(6)

式中:C为梭形柱类无量纲参数——扁率;a为梭形柱长轴尺寸,m;b为梭形柱短轴尺寸,m。

梭形柱两端边缘尖锐,而圆形柱、椭圆形柱表面轮廓光滑平缓,明显不同于前者。因此,为更深入地研究尖锐边缘对柱体气动特性的影响,在圆形柱、梭形柱的基础上,单独建立一个扁率C=0.1的椭圆柱作为数值试验对照组。圆形柱、梭形柱、椭圆柱几何形状对比如图1所示。定义来流方向为0°,所有测点沿逆时针均匀布置于柱体中部位置(Z=0.25 m),圆柱及梭形柱尺寸如图2所示。

图1 圆柱、梭形柱、椭圆柱示意(单位:mm)

图2 圆柱及梭形柱尺寸及测点布置

2.2 网格模型及计算参数

采用O型结构化网格对圆柱体及梭形柱体进行建模。SST湍流模型对近壁面的边界层网格要求较高,其无量纲壁面距离y+接近于1。通过试算,并进行网格无关性检验,确定近壁面第1层网格高度为2×10-2mm,增长率为1.2,网格参数如表1所示。梭形柱网格参数参照圆柱,仅改变短轴长度,网格模型如图3所示。为使壁面不影响流体流场和尾流发展,柱体中心距入口为10D,距出口为20D,其中D为结构特征长度。入口速度为35 m/s的均匀流,出口为远场压力条件,四周壁面均采用对称边界条件以防止壁面阻流对内流场造成的干扰。计算域如图4所示。

表1 圆柱几何尺寸及网格数

图3 柱体壁面网格划分

图4 计算域边界条件及尺寸

控制方程采用SIMPLIC,压力项采用standard格式,动量项采用二阶迎风格式。超临界条件圆柱绕流St介于0.2～0.3,估算周期为0.028～0.043 s,为更好地表现脱涡周期内特征,一个周期内至少计算20步,通过步长无关性检验,确定时间步距取0.001 s。

2.3 数值模拟合理性验证

先通过定长计算获得合理初始流场条件,再进行SST瞬态计算,待流场稳定达到统计状态,进行数据后处理。图5给计算后的圆柱体壁面平均y+分布。由图5可知,壁面y+处于0.09～0.92,满足近壁面网格尺寸的要求。

图5 圆柱表面平均y+分布

通过数值模拟得到的气动参数与文献[15]风洞试验数据对比列于图6和表2,可以看出SST计算结果和试验数据吻合较好,验证了网格划分的合理性和数值计算的准确性。

表2 数值模拟与文献试验数据对比

3 超临界条件梭形柱绕流气动特性分析

3.1 阻力、升力系数特性

图7、图8给出了超临界条件下不同扁率梭形柱阻力、升力系数计算结果。

由图7可以明显发现,随着梭形柱扁率增加,平均阻力系数也不断增大,阻力系数均方根也不断增大,说明扁率越大,梭形柱整体阻力变化幅度也不断变大。整体上平均阻力系数关于扁率近似成线性关系。因此,为方便估算梭形柱阻力系数,采用最小二乘法进行一次函数拟合,得到Re=7.2×105条件下的阻力系数估算公式。该拟合函数与拟合测点的相关系数为R=0.9875。

(7)

由图8可知,平均升力系数基本处于0位置,说明梭形柱两端升力随时间对称平衡变化,不存在类似圆柱在临界条件下的“单边泡”效应[15]。而升力波动方面,呈现出小扁率范围内(C≤0.3)升力均方根随扁率增加而增加,大扁率范围内(C>0.3)升力均方根随扁率增加而下降规律。当C=0.3左右时梭形柱升力波动最为剧烈。

仔细观察阻力、升力变化规律可以发现,C=0.1梭形柱长短轴相对差异度只有10%,外形上和圆柱体非常相近,但仅仅压缩较少短轴尺寸,柱体阻力系数就增加近1倍,升力脉动强度也急剧变化。分析原因,很可能是梭形柱两端存在尖锐边缘阻碍了流体在近壁面的附着,使得分离角提前到达,柱体阻力系数产生突变。为了验证上述猜想,通过另外建立扁率C=0.1的椭圆柱体模型,在同条件下计算三者气动参数。计算后的数据如表3所示。从计算结果可以发现,外形相近的3种柱体中,圆形柱和椭圆柱的气动性质较为接近,而扁率相同的椭圆柱和梭形柱,二者气动特性差异显著。椭圆柱两端平缓,阻力增加完全来自于扁率外形的变化。而同扁率条件下,梭形柱阻力系数也明显大于椭圆柱,说明阻力增加不但受类柱体扁率变化影响,同时还受到边缘尖锐程度的影响。证实了前述猜想,即边缘尖锐会明显增加风阻效应。其风阻增加机理将在第4节流场特性分析中进一步探讨。

表3 外形相似的3种柱体平均阻力系数对比

3.2 脱涡频率特征及斯脱罗哈数

由式(5)可知,St是f的无量纲化表示,不同扁率下梭形柱的斯脱罗哈数及脱涡频率对比如图9所示。从图中可以观察出,梭形柱体的脱涡频率在12～26 Hz,整体上随着扁率增加而降低,但在扁率较小(C≤0.3)范围内,脱涡频率起伏下降。进一步地,可以看出C=0～0.5范围内梭形柱体的St变化相对较小,基本处于0.18～0.22区间浮动,近似圆柱脱涡频率。从表3中可以看出,梭形柱与圆柱、椭圆柱的脱涡频率差异也较明显,说明了尖锐边缘也是影响脱涡频率的主要因素。当C≥0.6后,脱涡频率快速下降。大扁率与小扁率的脱涡频率差异明显。

图9 不同扁率梭形柱对应的斯脱罗哈数和脱涡频率

对以上现象,试做如下解释。实际上梭形是由两个圆形相交而得,除圆形半径为0.15 m外,其余梭形柱的相交圆半径分别为:R0.1=0.151 m,R0.2=0.154 m,R0.4=0.170 m,R0.6=0.218 m,R0.8=0.390 m。由半径变化规律可以发现,包括圆柱在内,前三者迎风等效半径相近,均为0.15 m左右,与圆柱半径接近;而后三者迎风等效半径较大。空气脱涡方向是沿着表面切线方向的,等效半径越大意味着迎风表面曲率越小,风速分离壁面后远离迎风中心,再回流形成旋涡的行程变长,频率变小。

3.3 风压系数分布特征

将不同扁率梭形柱的风压系数随角度变化规律绘于图10中。通过观察图10可以得到以下信息:①圆柱分离角大约在110°左右,而梭形柱的分离角基本都处于90°和270°处,即尖锐边缘处。②扁率越大,梭形柱表面随角度正压降到负压的降压速度越慢。③边缘负压系数和背向各测点负压系数均随扁率增加而增大,且背向各测点增加速率快,直至达到甚至超过边缘负压系数。

观察图11,根据脉动风压随扁率的变化特点,可以将梭形柱分为3个区域:0°～90°和270°～360°范围内迎风面风压波动随扁率变化并不剧烈;90°和270°边缘区,脉动风压急剧增大,边缘附近测点脉动风压强度对扁率变化敏感性最高;90～270°范围内背风区内,从梭形边缘延伸至中心线背风点,风压脉动效应对扁率变化的敏感度逐渐降低。

C=0.8对应梭形柱的平均风压和脉动风压分布明显区别于其他梭形柱,其特点主要有:①正风压在边缘处急剧突变到负压区,几乎没有过渡变化,迎风面和背风面的平均风压随角度变化均较小。②边缘风压系数与背风面各点风压系数差异已经不明显。③迎风面风压波动与小扁率梭形柱趋于一致;边缘测点波动强度远高于小扁率的;背风面风压波动规律与其他扁率梭形柱类似。从上述风压分布特征可以看出,随着扁率的不断增大,梭形柱绕流已经逐渐接近于理想无厚度平板,呈现出与类圆柱体明显不同的气动特性。

4 超临界条件梭形柱流场分布特性分析

4.1 速度场分布特征

图12给出了典型扁率下的平均风速分布情况。从图中可以观察到以下风速场特征:①各扁率梭形柱迎风向的风速变化基本相同。②在梭形柱两侧边缘存在明显的成耳朵状的风速增速区,且随着扁率增加“风耳”逐渐变细变薄。图中右上角显示了梭形柱边缘处的风速变化特征。可以看出,随着梭形柱扁率提高,“风耳”根部不断前移,边缘附近风速急剧下降。“风耳”前移、变薄意味着边缘范围内风速变化加快,根据伯努利原理就解释了第3节中边缘区脉动风压系数在边缘处随扁率突变明显的机理。③以来流平均风速35 m/s为参考,可以观察到,尾流区分为风速回升区和风速降速区。随着扁率增加,尾流的风速降速区长度在不断增加。而在大扁率梭形柱(C≥0.6)尾流区出现了逐渐回升至参考风速的风速回升区。风速回升区共有3个,分别对称分布在梭形柱边缘背风区域两侧,以及迎风轴线正后方,且随着扁率增大回升区域也不断扩大。

图12 不同扁率梭形柱的平均速度场分布

4.2 压力场分布特征

图13给出了梭形柱附近的压力场分布情况。可以很明显地从压力场中发现梭形柱压力分布的形成机理。①由速度场及压力场分布图可以分析出,尾流区的风速回升区不断增大导致尾流区逐渐形成了负压区,且随着扁率不断增大,负压区也逐渐变大直至靠近近壁面。可以观察到,正是由于负压区不断地接近背风面,使得背风近壁的负压系数不断增大。②负压区近似于圆形,随着区域不断扩大,其与背风面的接触点是由背风中心点向两边缘扩散,从而使背风中点负压沿壁面向边缘负压逐渐减小。③当负压区增大到完全包络住背风面后,边缘测点的压力突变几乎消失,与其他背风测点同在巨大的负压区内,负压系数趋于一致。

图13 不同扁率梭形柱的平均压力场分布

图13结合图10,可以解释阻力系数随扁率增加的机理。高雷诺数下柱体气动阻力主要来自于正负压差。从图中13中可以定性看出,随着扁率增加正压区变化并不显著,但负压区变化明显,由此带来的压力差增大,从而导致阻力增大。而从图10定量给出的风压分布也可以看出,随着扁率增加,梭形柱体正负压差增加明显。结合第3节猜想,进一步验证结论:梭形柱扁率越大、边缘越尖锐,其风阻效应越强。

4.3 涡量分布特征

考虑RANS方法对流场的时均化处理,故涡量瞬态流场图仅可表现出梭形柱涡量随扁率的变化的定性趋势。图14是超临界条件的梭形柱体三维Z向瞬态涡量。比较圆柱和梭形柱可以发现,圆柱的旋涡基本成垂直柱体状,说明圆柱绕流Z向主要是平面涡,空间相关性较高。而梭形柱Z向涡量呈现两端窄,中间宽形态,空间变化更加丰富,说明了梭形柱的旋涡具有空间分布特性。此外,可以看出随着扁率增加,梭形柱尾流涡的宽度不断发展变大,生成了体积更大形态更饱满的旋涡。形成大涡需要更长的时间,因此从旋涡形态上解释了St随扁率增大而变小的机理。

图14 不同扁率梭形柱的Z向瞬态涡量

为描述扁率对尾流区平均涡量的大小和位置的影响,定义梭形柱中心至平均涡心的距离为涡心距l,定义涡心至涡区最远边缘距离为涡半径r,如图15所示。图中给出了上述定义变量与梭形柱长边之比随扁率的变化规律。从图中可以看出,与圆柱体相比,梭形柱涡心显著逼近柱体,而涡心位置随梭形柱扁率增加变化并不明显。涡半径方面,在圆柱过渡到梭形柱时,涡半径先缩小,之后随扁率增加而显著增加。

5 结论

为研究超临界条件下(Re=7.2×105)梭形柱的绕流特性,基于SST模式的k-ω湍流模型采用Fluent软件对不同扁率梭形柱进行数值计算。通过对计算结果的分析,得到以下主要结论:

1) 梭形柱阻力系数平均值和均方差均随扁率增加而增加,且阻力系数平均值与扁率近似成线性关系。平均升力系数基本维持在0左右,升力脉动效应随扁率增加先增强,后减弱。

2) 在小扁率(C≤0.3)条件下,梭形柱脱涡频率与圆柱体相似,在大扁率(C>0.3)条件下,梭形柱脱涡频率随扁率增加而急剧下降。整体上扁率越大,脱涡频率越小。

3) 随着扁率增加,梭形柱尾流区的负压区逐渐扩大,直至与背风面接触,甚至包络住边缘,导致边缘风压系数与背风面各点风压差异随扁率增加而逐渐缩小。

4) 尾流区涡心到梭形柱中心距离对扁率变化不敏感;旋涡半径随扁率增加显著变大。