串列双锥柱绕流流动特性的大涡模拟研究1)

刘 健 邹 琳 陶 凡 左红成 徐汉斌

(武汉理工大学机电工程学院,武汉 430070)

引言

风能作为一种重要的可再生清洁能源,对其开发利用受到国内外学者的广泛关注[1-2].对于风能的获取,如涡激振动发电[3]、风力压电效应发电[4],增振就尤其重要.

圆柱作为绕流的经典研究对象,通常通过改变圆柱几何形状来改变周围流动状态以实现流动控制,Zhang 等[5]通过改变波浪外形圆柱的波长比以及波幅比,发现波幅比为0.2 时波浪圆柱受到的时均阻力系数显著降低.Zhang 等[6]通过对比研究无限长和一端自由的圆柱,发现有限长圆柱的时均阻力系数和脉动升力系数远小于无限长圆柱.赵桂欣[7]通过对比四种不同端面以及侧面的两端自由的有限长圆柱,发现改变外形对柱体阻力系数均有降低作用,在特定雷诺数下,脉动升力系数较单直圆柱有所提升.赵萌等[8]、乔永亮[9]的研究发现有限长圆柱长径比会明显改变圆柱的流场结构,存在一个临界长径比使得有限长圆柱的流动模式由自由端来流主导的无序涡脱模式转向侧面来流主导的卡门涡街涡脱模式.杨耀宗[10]的研究发现,有限长锥柱在斜率k=0.05 时锥柱振动性能有显著提升.由此可见,对于有限长圆柱,通过改变其表面形状对脉动升力系数有较大提升.

双圆柱较单圆柱的流动模式更为复杂.Alam等[11-12]通过实验,研究不同布置形式的双圆柱,对两个圆柱之间的尾流形态进行分类,研究发现在间距比为2.4～ 3.0,交错角为10°时的错列双圆锥脉动阻力出现最大值.Zdravkovich[13]研究两个串列圆柱间距比小于5 的情况,发现间距比小于3 时,上、下游圆柱之间不会产生明显涡脱,在间距比大于3 时,上下游圆柱之间逐渐产生涡脱.杜晓庆等[14-16]对不同布置状态下的双圆柱进行了数值仿真,按照上、下游圆柱相互作用的不同形式将流态进行分类,解释了下游圆柱在上游圆柱尾流的不同状态下的升阻力特性,发现下游圆柱被上游包裹状态下,下游圆柱的阻力系数显著降低,脉动升力系数反而受到一定抑制,但是相对于单圆柱,下游圆柱脉动升力系数均有不同程度增幅.文献[17-19]的研究发现,串列双圆柱随间距比增加,脉动升力系数和时均阻力系数先是缓慢降低,但是在间距比为3.0～ 4.0 附近发生突增跳跃到一个较高的水平,上游圆柱对下游圆柱的影响由剪切层包裹的抑制状态逐渐变为涡脱撞击状态,使得下游圆柱在超过临界间距比后脉动升力相对于单圆柱大幅提升.对于改变圆柱外形的串列双圆柱,唐涛等[20]数值模拟了间距比为3 的不同长度比的梯形截面串列柱,发现两圆柱之间的间隙流会随着梯形长度比的改变而呈现不同的耦合模式,导致升力发生改变.因此上游圆柱的尾流结构对下游圆柱脉动升力的提升有着决定性的作用,但是对于有限长串列双锥柱,其流动模式和升阻力特性缺乏研究.

综上所述,下游柱体脉动升力系数有可能相对于单柱体有大幅提升,同时表面的形状改变对于柱体脉动升力也有很大影响.根据课题组前期研究,发现锥形圆柱在斜率k=0.05 时相对于直圆柱增振最明显[10],因此,本文将以提升脉动升力系数为目的,利用商用软件Fluent 中的大涡模拟,在亚临界雷诺数(Re=3900)下,对斜率k=0.05 的串列双锥柱在间距比L/Dm=2～ 10 下进行数值模拟,通过分析上、下游锥柱之间的流动机理,探索下游锥柱时均阻力系数和脉动升力系数变化规律,为风力俘能结构群的列阵布局提供理论支持.

1 数值方法

1.1 控制方程

大涡模拟(large eddy simulation,LES)[21-22]采用空间滤波技术,通过截止尺度和滤波函数将涡识别为大涡和小涡,大涡直接解析,小涡则通过亚格子尺度应力模型(sub-grid-scale stress,SGS)封闭求解,过滤后的不可压缩黏性流体Navier-Stokes 方程为

式中,xi,yi为笛卡儿坐标,i,j=x,y,z,ui,uj为滤波速度矢量,为流体所受压力,ρ为流体密度,ν为流体运动黏度,τij为亚格力应力.

根据Smagorinsky[23]的基础SGS 模型,设应力表达式

式中,Δi为网格在i方向的尺度,Cs为Smagorinsky尺度,本文取0.2[24].

1.2 控制方程离散

本文计算控制方程离散均采用有限体积法,压力速度耦合方式为PISO,离散化均采用二阶格式[6].

2 计算域及网格

2.1 模型简介

基于多数人研究的直圆柱,本文加入斜率k,将直圆柱转化为锥柱,锥柱如图1 所示,几何表达式如下

图1 锥柱结构示意图及表面网格划分Fig.1 Schematic diagram of conical cylinder structure and meshing

式中,Dz为锥柱在高度为z处的横截面直径,Dm=0.01 m 为锥柱的平均直径,k为锥柱的斜率,圆柱高度比H/Dm=7,Dmax和Dmin分别为锥柱最大直径和最小直径.

本文中出现的相关参数定义如表1.

表1 相关参数定义Table 1 Parameter definition

2.2 计算域及网格划分

本文研究的串列双锥柱计算域大小为(30Dm+L) × 20Dm× 10Dm,X轴正方向为顺流向,Y轴为横流向,Z轴为圆柱展向,笛卡尔坐标系圆心位于上游锥柱中心处.上游锥柱中心距离计算域入口10Dm,下游锥柱中心距离计算域出口20Dm,锥柱的上下端面距离计算域顶部和底部均为1.5Dm,如图2 所示.入口采用速度入口(velocity-inlet),U∞=5.772 m/s,出口采用压力出口(pressure-outlet),背压设置为0 Pa,计算域两个侧面以及顶面采用对称边界(symmetry),计算域底面和圆柱体表面为无滑移壁面(noslip-wall).流体介质为空气,密度ρ为1.225 kg/m3,运动黏性系数ν为1.48 × 10-5m2/s.

图2 串列双锥柱计算域及网格示意图Fig.2 Schematic diagram of calculation domain and grid of tandem two conical cyliners

本文使用ICEM 对计算域进行结构化网格划分,利用O-block 技术对圆柱周围以及顶部网格进行加密,同时避免了片状低质量网格的生成.大涡模拟对于边界层的要求控制y+≤1,每个圆柱周围有3Dm×3Dm指数形式增长的加密区域,同时也对两个圆柱之间以及下游锥柱尾流区的网格进行了加密,如图2 所示.

2.3 计算模型验证

本文首先验证使用单直圆柱计算模型,探究时间步长Δt、网格参数(圆周节点数,第一层边界层高度)对圆柱的时均阻力系数(Cdmean),脉动升力系数(Clrms)和Strouhal 数(St)的影响,计算结果如表2.综合考虑计算成本和计算精度的影响,本文选用Case2 的网格精度和时间步长进行下面仿真.

表2 时间步长、圆周节点数和第一层高度对时均阻力系数,脉动升力系数和Strouhal 数的影响Table 2 The influence of time step and grid parameters on the Cdmean,Clrms and St

由于自由端数目以及长径比不同,本文单圆柱数值结果与仿真数据(num.)[25]以及实验数据(exp.)[26]有一定误差,为了进一步验证本次仿真所用的算法,网格精度以及时间步长能够适用于串列双锥柱,对串列双锥柱的流动结构进行烟线实验,其实验装置如图3 所示.风洞为直流式风洞(CSUWT2450-40),烟线控制器(SW-2) 通过充电放电加热涂有石蜡油的铜丝,石蜡油在瞬间高温下产生白色烟雾,风洞试验段均匀来流经过铜丝时,白烟在风的作用下绕过固定在风洞低端串列布置的锥柱,通过相机(Canon EOS 850 D)捕捉串列双波浪锥柱之间白烟形成的流场.对串列双波浪锥柱四个截面Y=0,S1(z/H=3/4),S2(z=0),S3(z/H=-3/4)与同截面下仿真结果进行对比.

图3 串列双锥柱实验台架Fig.3 Experimental of two conical cylinders in tandem arrangement

实验结果如图4 所示,串列双锥柱实现显示,Y=0 截面可以看到,间距比L/Dm=5 的串列双锥柱两个自由端来流撞击在下游波浪锥柱迎风面,这与仿真的结论一致,观察S1,S2,S3 截面发现,实验显示锥柱后方回流区发展完成,剪切层已经收缩,尾流撞击在下游锥柱表面,仿真结果与实验结论基本一致.因此可以认为本文验证中Case2 的算法,网格精度及时间步长可用于串列双锥柱仿真研究.

图4 间距比L/Dm=5 时串列双锥形圆柱烟线实验图与仿真流线图Fig.4 Smokeline diagrams and numerical streamline diagrams of two conical cylinders in tandem arrangement at spacing ratio L/Dm=5

3 结果分析

3.1 时均阻力系数和脉动升力系数

时均阻力系数和脉动升力系数是风力俘能结构获取风能的两个重要参数.本文与文献[18]上游柱和下游柱及Case2 中单圆柱数据对比结果,如图5所示,其中“UC”表示上游柱,“DC”表示下游柱.

图5 串列双锥柱时均阻力系数和脉动升力系数随间距比变化Fig.5 Cdmean and Clrms of two conical cylinders in tandem arrangement change with the spacing ratios

本文上、下游锥柱与文献中上、下游圆柱的Cdmean和Clrms随间距比变化规律基本一致,但是数值大小存在很大差异,可能是因为本文锥柱自由端的影响,下文将通过流场结构分析其原因.本研究发现,Cdmean和Clrms在L/Dm=4～ 5 时存在突变,在这个区间可能发生流动状态的改变;在临界间距比之前上、下游锥柱Clrms呈下降趋势,而Cdmean则随间距比先增加后减小;在间距比超过临界值之后,上游锥柱Clrms在临界值之后基本保持不变,且与单直圆柱基本一致.下游锥柱Clrms在间距比L/Dm=5～ 6有一个略微降低然后突增的过程,在L/Dm=7 的时候出现峰值,此时较单直圆柱增长约20.7 倍,此阶段可能存在流动结构改变.上游锥柱Cdmean在低于单直圆柱一定值处波动,下游锥柱Cdmean由一个较小值逐渐增大,趋近于上游锥柱,下游锥柱相对于单直圆柱Cdmean最大减小约90.7%,在L/Dm=7 时,减小约19.8%.

3.2 表面压力系数

表面压力系数Cp能够直接反映出锥柱表面受力情况,为了方便观察表面Cp分布情况,如图6(a)所示,将锥柱的表面沿母线剪开,并以角度“θ”和高度“Z/H”作为横、纵坐标铺平,表面时均Cp分布如图6(b)所示.

图6 间距比L/Dm=2～ 9 下串列双锥柱时均压力系数分布图Fig.6 Distribution diagram of time-average Cp of two conical cylinders in tandem arrangement with spacing ratio L/Dm=2～ 9

从整体来看,无论上游锥柱还是下游锥柱在不同间距比下的时均Cp分布基本关于180°对称,说明来流对上、下游锥柱两侧的影响在时间平均上基本相同.根据云图分布的不同,可将上游圆柱的Cp分为两类,第一类为低压区(压力系数小于-0.6 的区域)集中在锥柱顶部(DIS1),如L/Dm=2,3,5,第二类分为高压区(压力系数大于0.6 的区域)主要集中在锥柱底部(DIS2),如L/Dm=4,6,7,8,9.可以看到上游锥柱除L/Dm=2,3,4,低压区强度相对较大以外,其他间距比分布基本相近,这是Cdmean在临界区域前端逐渐降低的主要原因;由于下游锥柱受到上游锥柱尾流的影响,在小间距比下压力分布强度较弱,随着间距比增加,高压区和低压区范围和强度都有所增加,并且分布形式逐渐接近上游圆锥,这也是下游锥柱阻力系数随间距比增加而增大的主要原因(临界区域除外).锥柱的表面压力分布主要是由其流场决定,下面从流场的本质对压力分布进行分析.

3.3 时均流线图

三维流场较为复杂,但串列双锥柱的流场分布具有一定对称性,下面取Y/Dm=0,S1,S2,S3 截面(见图3(b))对串列双锥柱的时均流线进行分析,如图7 所示.

图7 间距比L/Dm=2～ 9 下串列双锥柱在Y=0,S1,S2,S3 截面时均流线图Fig.7 Time-average streamline diagram of two conical cylinders in tandem arrangement with spacing ratio L/Dm=2～ 9 at Y=0,S1,S2,S3 sections

由于自由端的存在,锥柱上、下端部流体分别存在上洗(流经下自由端的空气向上运动)和下洗(流经上自由端的空气向下运动)现象,上洗、下洗作用会抑制侧面来流在锥柱后方形成周期性涡脱落,从而导致脉动升力系数降低,这也是串列双锥柱脉动升力远小于串列无限长圆柱的主要原因[27-28].空气流经上游锥柱在其后方形成一对Y截面的展向回流区(arch-vortex),同时空气流过柱体侧面,与arch-vortex 相互作用,在柱体后方形成Z截面的横向回流区(side-vortex,以下简称回流区),除此之外,在柱体端面处能观察到顶涡(tip-vortex)[29-30],如图7(g).课题组前期研究发现[10],由于锥柱斜率的存在,导致自由端上洗作用强于下洗,下端形成的arch-vortex 挤压上端形成的arch-vortex 使其变为椭圆形,对于串列双锥柱,在L/Dm≥ 7 时,上游锥柱后方能够明显观察到一对arch-vortex,但在间距比L/Dm=2～ 6 时,受到下游锥柱影响,上游锥柱自由端形成的涡结构发生明显变化,根据上洗和下洗在两个圆柱之间形成的arch-vortex 不同,将两锥之间的流场形式分为两类,第一类是由下洗作用占主导,流体经上自由端向下运动一部分锥柱上端形成一个涡,另一部分与下自由端的上洗流相作用在锥柱下端形成一个涡,如L/Dm=2,3,5,第二类则与第一类相反,如L/Dm=4,6,7,8,9.这两种分类方式正好也解释上一节中DIS1 和DIS2 的分布规律,占主导作用的一端,空气流速更快,上游锥柱背风面对应位置压力也就越低.而无限长串列圆柱上游柱背风面不存在自由端影响,背风面形成的回流区主要来自侧面,压力系数展向对称分布[6].串列双锥柱间流动结构随间距比变化可分为三种状态:剪切层包裹状态,过渡状态,尾流撞击状态.剪切层包裹状态,上游锥柱回流区受下游锥柱影响并未完全发展,剪切层完全包裹住下游锥柱,如图7 中L/Dm=2,3;尾流撞击状态,上游锥柱当前位置回流区已经完全发展并且大小随间距比基本不发生改变,且下游锥柱已位于回流区之外,不再被上游剪切层包裹,如图7 中间距比L/Dm=7,8,9;过渡状态,即上游锥柱自由端来流作用到下游锥柱表面,而侧面剪切层部分包裹下游锥柱,如图7 中间距比L/Dm=4,5,6.随着三种状态变化,上游锥柱的上洗和下洗直接作用在下游锥柱表面的范围逐渐减小,同时下游锥柱后方的回流区随着间距比增加而逐渐扩大,这也是下游锥柱随间距比增加,其高压区和低压区范围逐渐扩大并增强的主要原因.

3.4 涡结构图

涡的识别方法可以分为三代[31],其中第二代中的Q方法运用较为广泛,表达式如下

式中,A为对称张量,B为反对称张量,为矩阵的Frobenius 范数.下面对串列双锥柱的涡采用Q方法识别,如图8 所示,每个间距比提供两个视角view1,view2.

图8 间距比L/Dm=2～ 9 下串列双锥柱瞬时涡量图(Q=1 × 104),view1 为整体视图,view2 为俯视图Fig.8 Instantaneous vorticity diagram of two conical cylinders in tandem arrangement with spacing ratio L/Dm=2～ 9 (Q=1 × 104).View1 is the general view;view2 is the top view

从涡量图中可以看到,前面提到的两种分类,DIS1 在这里显示的是大量涡“下泄”集中在两个锥柱之间的中下部,如图8(a)(b)(d),DIS2 在这里反应的是大量涡“上泻”集中在两个锥柱之间的中上部,如图8(c)～图8(h).从涡的结构上来看,间距比较小时,上游锥柱脱落的涡碎且无序,与下游圆柱产生的涡相互作用,在下游圆柱后方形成“肋状涡”(rib vortex),如图8(a)中所示,当间距比逐渐增大时,上游锥柱脱落的涡逐渐变大并且有周期特性,下游锥柱后方“肋状涡”也逐渐变大变长,如图8(e)～图8(h)所示.而对于无限长圆柱,不同间距比下均能观察到卡门涡街现象[19],所以其脉动升力系数能明显高于本文有限长双串列锥柱.上游锥柱脱落的尾涡作用在下游锥柱表面,与来自下游锥柱端部以及侧面来流相互作用,产生较大的升力,相较于低间距比串列双锥柱,当间距比为L/Dm=7～ 9 时,下游锥柱处于上游锥柱回流区临界位置或回流区下游,上游锥柱尾流能够充分发展,形成类似卡门涡街特征,周期性的漩涡作用在下游锥柱表面,使下游圆锥表面更容易产生大的升力.间距比L/Dm=7 时,上游锥柱尾流发展最充分,其尾流对下游锥柱的作用更为强烈,这可能是下游锥柱Clrms相对于间距比L/Dm=8,9 更大的原因.下面通过分析几个截面的涡量,来具体分析上下游锥柱之间的流动情况.

3.5 Z 截面涡量图

为了方便观察串列双锥柱之间的流动形态,下面取S1,S2,S3(见图3(b)) 三个截面的涡量图,如图9 所示.尾流之间的干扰可分为很多不同形式,根据锥柱截面直径以及流场特性,参照文献[12-13]对流场的分析,将Z截面涡量分为两类,V1:上游锥柱尾流包裹下游锥柱,如图9 中双点划红线外框中涡量图,上、下游锥柱共同形成的剪切层极不稳定,在下游锥柱后方形成的回流区相较于单锥柱[8]很小,在后方形成非周期性的涡脱;V2:下游锥柱位于上游锥柱回流区外端或恰好位于其回流区边界处,如图7 中虚线黑线外框中涡量图,上游锥柱尾涡得到充分发展,尾流撞击在下游锥柱表面,结合下游锥柱产生的涡在下游锥柱后方产生大量不规则碎涡,相较于V1 为下游锥柱提供更大的升力.结合图8 中涡量图显示,上游锥柱的上洗和下洗作用强度不同,使得作用在下游锥柱前端的涡集中位置不同,如图8(f),上游锥柱尾涡集中在锥柱中下部,导致图9 中S1 截面两个圆柱相互作用较弱,这也是图6 中间距比L/Dm=4 的下游圆柱顶部时均压力系数高的主要原因.图9 中,当间距比L/Dm=2,3 时,三个截面流动模式均为V1,即4.3 节提到的剪切层包裹状态,间距比L/Dm=4,5,6 中,三个截面中V1,V2 模式共同存在,即过渡状态,而当间距比L/Dm=7,8,9 时,三个截面流动模式全部变为V2,即尾流撞击状态,这是间距比L/Dm=7,8,9 有更高的脉动升力值主要原因.从图9 中可以观察到,图9(f) 相对于9(g) 和9(h)每个截面下游锥柱距离上游锥柱回流区更近,上游锥柱发展的周期性尾涡作用在下游锥柱表面更剧烈,这是下游锥柱在L/Dm=7 取得Clrms幅值的主要原因.

图9 间距比L/Dm=2～ 9 下串列双锥柱在S1,S2,S3 截面的涡量图Fig.9 Z-vorticity of two conical cylinders in tandem arrangement with spacing ratio L/Dm=2～ 9 at S1,S2,S3 sections

4 结论

本文在Re=3900 下利用大涡模拟对间距比为L/Dm=2～ 10 的串列双锥柱进行了固定绕流仿真.首先,验证了Case2 时间步长以及网格达到收敛,通过烟线实验证明本次网格精度,时间步长和算法可用于串列双锥柱仿真;其次,分析了其升阻力系数随间距比变化的特性;再次,分析了其时均流场信息分析其阻力变化的主要原因;最后,通过涡量图分析了脉动升力系数变化的主要原因,得出以下结论.

(1)随间距比增加,上游锥柱Clrms先减小逐渐趋近于单直圆柱,在L/Dm=4 之后基本保持不变;上游锥柱Cdmean先减小,在L/Dm=4 之后保持低于单直圆柱的波动状态.

(2)下游锥柱Clrms随间距比基本保持增加后减小趋势(除临界区间),在L/Dm=7 时,取得峰值,约为单直圆柱的21.7 倍;下游锥柱Cdmean由一个较小值随间距比增加逐渐接近上游锥柱的Cdmean,下游锥柱相对于单直圆柱Cdmean最大减小约90.7%,在L/Dm=7 时,减小约19.8%.

(3)相较于无限长串列双圆柱,串列双锥柱的上游锥柱自由端的上洗和下洗作用能显著影响侧面来流在锥柱后方形成的涡脱的发展,抑制其产生规则卡门涡街现象;上洗和下洗作用强的一端会使上游锥柱背风面形成大且强的低压区,同时使同端下游锥柱迎风面的高压区减小.

(4)随间距比增加,下游锥柱与上游锥柱回流区的关系有剪切层包裹状态、过渡状态、尾流撞击状态,其中尾流撞击状态时,上游锥柱尾流能够充分发展,形成交替脱落的肋状涡,撞击在下游锥柱表面,使下游锥柱产生较大的脉动升力系数,间距继续增加,这种相互作用的强度会随之减弱,脉动升力系数随之减小.本文的研究能为风力俘能结构的列阵布局提供理论支持.