林海燕, 向 阳, 张 斌, 刘 洪
(上海交通大学 航空航天学院, 上海 200240)
涡结构广泛地存在于自然界及工程应用中,由于涡环是3维流场中最基本的涡系结构,众多学者开展了对其特征性能的研究.Shariff等[1]指出研究涡环可以揭示涡动力学的一些物理特性及基本原理.针对不可压缩涡环,研究者从涡环动力学的特征[2-3]、极限生长原理[4-5]、多涡环相互作用[6]等多个方面开展了大量研究.但是,相比于不可压缩涡环,可压缩涡环的研究十分有限且大多局限于研究可压缩涡环的流动结构特征,如针对涡环夹止和尾迹射流[7]及由激波气泡相互作用产生的涡环的生长演化过程[8]的研究,几乎均未揭示可压缩性对涡环物理特征的影响.对于可压缩流动,一般而言可以分为亚声速、跨声速和超声速3个区域,其流场存在不同的压缩性特点.不同的压缩性会使可压缩涡环产生不同的结构特征.
一般采用激波管产生轴对称的可压缩涡环,而通过调节激波管驱动段和被驱动段之间的压力比可以改变入射激波的马赫数(Mashock),即激波到达激波管出口时的马赫数,进而得到具有不同可压缩性的涡环.Elder等[9]在1952年首次通过纹影实验观测到了由激波管产生的可压缩涡环向下游传播的过程.Arakeri等[10]利用粒子图像测速(PIV)技术研究当Mashock=1.1,1.2,1.3 时的可压缩涡环,发现涡环完全形成后的传播速度是激波到达激波管出口时波后速度的0.7倍.Dora等[11]同样通过PIV技术研究Mashock=1.27,1.37时的可压缩涡环,获得了涡环结构参数的定量数据,并验证了由实验测得的涡环传播速度与理论计算公式求得的结果相符.此外,部分学者还关注了可压缩涡环结构中存在的激波.Baird[12]通过差分干涉测量法对可压缩涡环进行研究的结果表明,当Mashock=1.5时,涡环的径向区域存在嵌入后向激波,但该激波在强烈黏性剪切的作用下终止,并不能传播至涡核中心.Brouillette 等[13]通过电花影图和纹影拍摄技术研究不同条件下激波管产生的可压缩涡环,发现能够产生嵌入激波的Mashock的最小值为1.34.进而,Brouillette等[14]根据激波的结构特征对可压缩涡环进行了更细致的分类:利用纹影拍摄技术发现在临界驱动长度条件下,当Mashock<1.43时,将会产生没有嵌入激波的涡环;当Mashock≥1.43时,则产生带有嵌入激波的涡环;当Mashock=1.6时,嵌入激波的逆压梯度将产生二次涡,这种二次涡被称为反向旋转涡环.随后,又有学者针对反向旋转涡环开展了更为细致的研究,如:Kontis 等[15]通过高速纹影拍摄技术观察到了当Mashock=1.63时,主涡环下游存在多重反向旋转涡环;Thangadurai等[16]通过石蜡油烟雾显示技术同样观测到了反向旋转涡环并详细地研究了其形成机理及演变过程.
上述研究结果均表明,对于可压缩涡环,可压缩性的强弱对涡环流动结构的影响有很大的不同.然而,有关涡环自身的可压缩性描述及可压缩性对涡环自身结构参数影响的较为系统的研究尚未见报道.激波结构的存在与否是一个可对涡环结构的稳定性产生重要影响的因素,而涡环自身的结构参数(如涡核结构,涡环半径,传播速度)对于研究由涡环结构引起的夹带、混合、燃烧等现象均十分重要.本文通过有限体积法对激波管产生的可压缩涡环进行数值计算,以研究涡环自身的可压缩性及其在不同可压缩性条件下对涡环半径、涡核半径、涡环传播速度以及涡核内涡量分布的影响.
采用课题组自行研发的软件进行数值计算;采用有限体积法对3维可压缩Navier-Stokes方程进行空间离散;采用5阶WENO(Weighted Essentially Non-Oscillatory)格式重构对流项;时间推进采用3阶Runge-Kutta法;针对激波管产生涡环这个物理现象,为了方便验证,采用与文献[17]相同的物理模型.激波管计算模型示意图及各组成部分的结构尺寸如图1所示.激波管由驱动段A和被驱动段B组成,内径为64 mm;C和D为环境区域,即涡环的形成区域及运动区域;被驱动段和环境区域的初始压力为101.325 kPa;整个计算区域的初始温度为300 K;改变驱动段的压力即可得到不同驱动段与被驱动段的压力比(Rp),进而得到不同的Mashock.定义入射激波到达激波管出口时时间t=0;激波管出口处x=0.
通过改变涡环生长演化最主要区域(0 mm 表1 网格信息 为了验证计算方法的正确性,在Rp=10及驱动段长度L=230 mm的条件下进行计算TS1-1,并将计算结果与文献[18]中相同条件下算例的计算结果(TS2)进行比较,涡环生长过程中不同时刻涡环截面密度梯度云图对比结果如图3所示.由图3可见:TS1-1与TS2的涡核中心位置以及嵌入激波位置均较为一致(见图3(a));TS1-1与TS2在对称轴附近产生的滑移线位置也较为一致(图3(b));TS1同样得到了反向旋转涡环结构(图3(c)).进而在Rp=7,L=165 mm的条件下进行计算(TS1-2),并将获得的计算结果与文献[17]的计算结果(TS3)进行对比.当t=560 μs时,流场中当地流体沿x轴横向穿过涡环中心线上的径向速度(v)的分布规律如图4所示,其中TS3-m为TS3的实测值.由图4可见,TS1-2与TS3的实验及计算结果无论是趋势还是拐点数据均较为吻合,从而验证了所提计算方法的准确性. 图2 不同网格条件下Γ随t的变化规律 图3 TS1与TS2的涡环截面密度梯度云图对比结果 为了研究不同可压缩性条件下产生的涡环,设计5个具有不同Rp值的算例,每个算例对应的Mashock(由激波到达出口时根据正激波关系计算获得)取值如表2所示. 图4 径向速度分布规律对比 表2 5个算例的激波马赫数 因此,在涡环的随体坐标系下,流函数ψ0为 (1) 式中:u0为涡环传播速度,由各时刻涡量最大值点的所在位置确定.为了消除取涡量最大值所在位置时的测量误差,对时间及位置曲线进行高斯拟合[21],得到u0的平均值. ψ0的计算步骤如下: (1) 定义在对称轴(r=0)上,ψ0=0;以ψ0=0的点为起点,根据速度场(u,v)进行积分,从而求得整个计算平面的ψ0值. (2) 获得整个区域ψ0值的曲线分布,可知ψ0=0(不包括对称轴r=0)形成一个封闭曲线,该曲线即为涡环边界.根据算例5计算得到的流函数分布如图5所示.从图中可以清晰地看到涡环边界ψ0=0.鉴于涡环的对称性,在ψ0=0的上半部分及r=0所包围的区域计算环量. 图5 流函数分布示意图 涡核内的涡量分布是一个对理论分析(数学角度)及涡环的特征描述(物理角度)非常重要的参数.目前,已有多个用于描述涡核内涡量分布的模型被提出.Saffman[22]的研究表明,在真实流体中不可压缩涡环的涡核内涡量分布满足高斯分布. 对涡量分布和位置的无量纲计算公式如下: ωn=ω/ω0 (2) yn=(y-y0)/rc (3) 式中:ω为涡量;ω0为峰值涡量;y为ω对应的径向位置;y0为ω0对应的径向位置;rc为涡核半径. 可压缩涡环由于自诱导速度而向下游传播.Moore[23]提出了该传播速度的理论计算公式: (4) 式中:Γ可通过对由流函数方法确定的边界区域的涡量进行积分得到;rr为涡环半径;Mavortex为涡马赫数,可用于表征涡环的可压缩性, (5) c∞为无穷远处的声速. 激波离开激波管出口后向下游运动的同时发生衍射,导致激波管出口压力减小,产生向激波管上游传播的膨胀波,使得激波管内流体速度增大;被加速的流体在激波管出口形成剪切层后卷起形成涡环,涡环不断生长并向下游传播;涡环完成生长后,从激波管流出的流体会形成尾迹射流.尾迹射流的长度及持续时间取决于激波管驱动段的长度和压力. 图6 当Rp=3时 的演化过程云图 图7 当Rp=6时 的演化过程云图 可压缩涡环根据结构特征可分为不含嵌入激波的涡环及含有嵌入激波的涡环.Kontis等[15]提出根据Mashock是否大于1.43来定量判断是否存在嵌入激波,但Mashock只能作为判断涡环是否会产生的一种参考值,对于表征涡环可压缩性的相对强弱并不具有物理意义.Baird[12]指出把含有嵌入激波的涡环称为超声速涡环,并根据由涡环对称轴上最大轴向速度计算得到的局部马赫数(Malocal)定量区分亚声速和超声速涡环.此外,Moore[23]通过理论推导提出Mavortex也可用于表征涡环的可压缩性. 传统的空气动力学把Ma<0.8称为亚声速,Ma=0.8~1.2称为跨声速,Ma>1.2则称为超声速.而在可压缩涡环的研究中,并没有系统地研究过涡环自身定量的可压缩性.在不同Rp条件下,3种表征可压缩性的马赫数(Mashock,Malocal,Mavortex)的取值变化趋势如图9所示.由图9可见:①Mashock与Rp近似成正比,且都大于1,表明都为可压缩涡环.② 当Rp≥4时,Malocal均大于1.根据文献[12]的定义,Rp≥4产生的皆为超声速涡环,但当Rp=4时Malocal只是略大于1,结合图8可知,此时产生的涡环中并不存在嵌入激波,说明产生的不是传统意义上的超声速涡环,而更趋向于跨声速涡环.③ 当Rp=4时,Mavortex略小于1,表明此时的涡环处于跨声速涡环的范围;当Rp=5时,结合图8可知该条件下的涡环存在嵌入激波,Malocal也远大于1,表明此时产生的是超声速涡环,但其Mavortex的值反而小于1.导致这种现象的原因可能是涡核受到了激波的影响,轴向被拉长,径向被压缩,使得涡核的半径产生变化,导致Mavortex的数值受到了影响,但这种影响在Rp=6之后就不再显著. 通过以上分析,为了使涡环可压缩性的定量表征能够与传统的流动可压缩性相对应,将涡环也区分为亚声速,跨声速和超声速3种特征,通过Malocal和Mavortex可以很明显地区分具有弱压缩性的亚声速特征涡环和具有强压缩性的超声速特征涡环.在跨声速区域,由于从是否存在嵌入激波这一现象上无法区分可压缩涡环的特征,只能通过Malocal和Mavortex的值作近似判断.由于Malocal的计算较为简便,故可通过其对涡环压缩性的强弱进行预判;而Mavortex是基于整个涡环的结构参数计算获得的,因此更能够代表涡环自身的可压缩性,同时也更具有实用意义. 图8 不同Rp下主涡环ω及 的分布云图 图9 3种表征涡环压缩性的Ma随Rp的变化 本节分析Mavortex与可压缩涡环的相关结构参数.描述涡环的重要参数包括ω,rc,rr及u0,下面将对这些参数进行定量分析. 图10 不同Mavortex条件下ωn分布及其相应的高斯分布对比 涡核内无量纲涡量沿yn的分布(ωn)规律以及在对应条件下的高斯分布情况如图10所示.其中,yn=0为涡核中心.由图10可见:随着Mavortex取值的增加,涡量的聚集程度也有所增加;可压缩涡环中涡核内涡量的分布不满足高斯分布,并且随着Mavortex的增加,涡核内涡量更加集中,也愈加偏离高斯分布. Norbury[24]根据涡环的无量纲半径对涡环进行分类.可压缩涡环直径的定义如图11所示.涡环直径取涡量绝对值最大的两点之间的距离,涡环半径则取为该距离的1/2;涡核直径由穿越涡核中心线(平行于y轴)上的轴向速度分布中最大值和最小值的所在位置确定(见图12), 涡核半径则取为涡核直径的1/2. 图11 通过涡环截面涡量分布定义的涡环直径 图12 通过涡核中心u的分布定义的涡核直径 涡环半径和涡核半径随Mavortex的变化趋势如图13所示.由图13可见:随着Mavortex的增加,涡环半径呈逐渐增加的趋势;涡核半径在Mavortex较小时呈逐渐增加的趋势,而在Mavortex>1.2时,即出现相对较强的嵌入激波后,由于激波的拉伸作用,使得涡核径向半径略微变小;在Mavortex接近1的区域,Mavortex对两个半径的影响突然加剧,表明可压缩性开始对涡环产生更剧烈的影响,也可以认为涡环进入超声速特征的区域. 图14 u0及随Mavortex的变化情况 本文结合有限体积法求解3维Navier-Stockes方程,对激波管产生的可压缩涡环进行数值模拟,计算结果与文献的相关实验结果及计算结果相符合;分析涡环的演化过程及其相关结构;给出了定量表征涡环可压缩性的参数,并定量分析可压缩性对涡环参数的重要影响;揭示涡环在可压缩性作用下自身物理特征的变化情况.研究主要得到以下结论: (1) 可压缩涡环的可压缩性可以通过Malocal及Mavortex定量表征.由于Malocal的计算方法较为简单,所以可通过其对涡环压缩性的强弱进行预判;由于Mavortex更能代表涡环自身的可压缩性,所以更具有实用意义. (2) 激波管产生的可压缩涡环能够表现出亚声速、跨声速和超声速3种特征,根据Malocal和Mavortex的取值可以很明显地区分具有弱压缩性的亚声速特征涡环和具有强压缩性的超声速特征涡环.在跨声速区域,由于从嵌入激波这一现象上无法区分可压缩涡环的特征,只能通过Malocal和Mavortex的值作近似判断;而Mashock只能表征产生涡环的激波的强度,对于可压缩涡环的特征判定不具备实际意义. (3) 随着可压缩性的增加,涡环半径呈上升趋势,涡核半径由于嵌入激波的存在使得其在超声速特征区域略微有所减小,而涡环的传播速度逐渐增加;在3种可压缩性特征区域计算得到的涡环传播速度和理论公式计算结果相符合,表明该理论计算公式可同时适用于3种特征区域,从而可以更为准确地预测涡环的传播速度;可压缩涡环涡核内涡量的分布不符合高斯分布,随着可压缩性的增加,涡量更加集中. 致谢感谢上海交通大学高性能计算中心超级计算机Л提供的计算帮助.1.2 流函数确定涡边界法
1.3 涡核内涡量分布的无量纲化
1.4 可压缩涡环传播速度计算
2 结果与讨论
2.1 可压缩涡环的形成及其演化过程
2.2 涡环可压缩性的定量表征
3 可压缩性对涡环物理特征的影响
3.1 涡核内涡量分布
3.2 涡环半径和涡核半径
3.3 涡环传播速度
4 结论