邓 雄,刘 强,王秋旺,罗振兵,程 盼

(1. 西安交通大学 能源与动力工程学院, 陕西 西安 710049; 2. 国防科技大学 空天科学学院, 湖南 长沙 410073)

射流冲击冷却是流体在一定压差下通过喷孔高速喷出冲击换热表面进行冷却,其具有较高的局部换热效率,在内燃机冷却、金属热处理、航空发动机叶片冷却以及电子器件冷却等领域具有广泛的应用[1-2]。合成射流无须流体供应装置和管路输送系统,通过流场环境本身中的流体来“合成”周期性的非定常射流,是一种零质量射流技术[3]。合成射流的本质是一系列非定常涡结构的形成、传播和相互作用。涡结构强烈的卷吸作用和湍流度使得合成射流冲击热表面时,加速了热边界层不断破坏和重建,能有效强化换热[4-5],与传统连续射流冲击冷却相比,相同雷诺数条件下,其散热量是连续射流的3倍多[6]。因此,合成射流在强化换热领域[7-8],尤其是受限空间电子器件散热方面,具有较好的应用前景。美国通用电气公司2012年在拉斯维加斯国际照明展中展示了采用合成射流与翅片组合冷却技术的Energy Smart®LED灯泡,被认为是解决LED器件散热问题的重大机遇[9]。

合成双射流技术是一种性能更优的合成射流技术,其核心作动器结构是两个射流出口、两个腔体和一个振动膜,振动膜在往复振动过程中不断地压缩和膨胀两侧腔体内的流体,在两个出口的剪切作用下交替形成两股相位差为180°的合成射流,两股射流逐渐向下游传播最终融合成一股合成双射流[10-12]。合成双射流除了具有常规合成射流的优越性能外,能量利用效率和射流频率都提高了1倍,控制能力和环境适应性更强,而且合成双射流是两股合成射流融合而成的,其包含的涡结构更丰富,对于强化换热更有利。另外,合成双射流冲击冷却具有独特的射流矢量控制特性,能实现对大面积电子器件的高效、均匀散热和动态“热点”去除,应用于电子器件散热更具优势[13-14]。

合成射流冲击平板属于对流传热范畴,其散热冷却性能取决于冲击射流的流场特性。因此,了解合成双射流冲击平板的流场结构,对于认识其传热机理和强化传热性能至关重要。McGuinn等[15]采用高速粒子图像测速和热线风速仪研究了非对称合成射流冲击光滑平板的流场特性,并根据流场拓扑结构和尾迹流强度,将冲击合成射流随斯托克斯数的增加分为四个不同区。Xu和Wang[16-17]通过激光诱导荧光和粒子图像测速方法研究了合成射流涡环撞击平板的流场结构演化过程,揭示了斯托克数和雷诺数对合成射流涡环冲击平板的影响。Greco等[18]研究了斯特罗哈数和冲击距离对合成射流冲击流场的影响,得出斯特努哈数反映了涡环强度和尾迹流强度的相对大小,冲击距离主导涡环发展和固壁引起的逆压梯度。Silva-Llanca等[19-20]分析了合成射流涡动力学特性,揭示了合成射流冲击冷却强化换热机制。Zhao等[21]采用锁相粒子图像测速技术研究了矩形出口合成射流体流场结构的演化过程,获得了三维涡结构和表征流体与壁面相互作用的扫掠涡结构。Lü等[22]研究了活塞式合成射流冲击平板的传热性能,分析了花瓣型、弧形射流出口的数量和长宽比对冲击换热的影响。

本文采用大涡模拟方法对合成双射流冲击平板进行数值模拟,基于有限时间Lyapunov指数方法对其流场的拉格朗日涡系结构进行识别,并结合本征正交分解方法,分析流场的模态特征,为合成双射流冲击冷却布局设计和机理研究奠定基础。

1 计算模型与处理方法

1.1 计算模型

本文采用大涡模拟方法进行数值计算分析,亚格子模型为Smagorinsky-Lily模型,离散格式采用PISO(pressure implicit with splitting of operators),时间项采用二阶Crank-Nicolson格式进行推进求解,压力项和对流项分别采用二阶空间离散格式和具有三阶精度的MUSCL格式。图1为合成双射流冲击平板物理计算模型二维示意图。顶部为冲击面,由于本文着重于对流场结构进行研究,不考虑冲击散热特性,因而设置为绝热无滑移壁面。下部为合成双射流激励器,其压电陶瓷膜直径为0.046 m,两出口宽为0.002 m、深0.004 m,出口间隔0.005 m。坐标原点设为膜片的底部,以膜片振动方向为x方向,垂直方向为y方向。

图1 计算模型Fig.1 Computational model

采用准二维假设,即在三维展向方向上对网格进行拉伸计算,计算域的长宽高分别取109 mm×5 mm×30 mm。对于上部计算域,开展了三组不同数量的网格进行验证,分别是粗网格(501×11×251)、中等网格(601×11×301)和细网格(601×11×351),并对壁面和出口附近进行局部加密。图2给出了三组网格下时均流场中左侧出口中心线上的速度分布。从图中可以看出,对于粗网格,其在冲击面附近的速度值偏低,而中等网格和细网格的速度曲线能够很好地吻合。考虑计算代价,采用中等网格开展研究,总网格数约200万。所用计算方法已在参考文献[23]得到了充分验证。式(1)给定激励器左右膜片速度入口条件:

图2 网格无关性验证结果Fig.2 Validation of the grid independence

(1)

即激励器驱动频率为100 Hz。计算域顶部为平板,四个侧面给定压力出口边界条件,下边界采用无滑移固壁条件,数值计算时间步长取2×10-5s。

1.2 有限时间Lyapunov指数方法

涡系结构的识别与分析手段较多,总的来说,分为基于欧拉体系和拉格朗日体系两种辨识方法。欧拉体系下涡结构的提取分析是基于当地速度梯度张量,要求具有较高空间分辨精度的三维速度场,并且具有一定的当地性,在物理客观性方面存在一定的问题。拉格朗日体系的涡结构辨识是以流体质点在流场中对流特性的时间积分为研究对象,而不是使用速度梯度张量,有限时间Lyapunov指数(finite-time Lyapunov exponents,FTLE)方法[24-26]即是其中一种。

(2)

其最大线性指数增长率即为流体粒子最大的轨迹拉伸率。

1.3 本征正交分解方法

本征正文分解(proper orthogonal decomposition, POD)的原理是在优化目标场变量均方差的基础上提取模态,是一种有效的识别显性特征和时间的方法。近些年,POD方法已在合成射流流场分析中得到了广泛的应用,如合成射流控制圆柱绕流后驻点的涡旋分析[27]、受合成射流激励的后台阶分离流分析[28]、合成双射流流场特征结构[29]等。

考虑到物理量Ω∈2存在边界∂Ω,假定有以下随时空变化的目标场u(x,t)∈L2(Ω),POD的基本思想则是在有限集中寻找一组最优正交基ai(t)和φi(x)来表示目标场[30],即

(3)

式中,φi(x)是空间模态,ai(t)是与每个空间模态对应的时间系数。寻找最优正交基的问题可以进一步转化为求解目标场矩阵的N个特征值问题,也就是Meyer等发展的snapshots POD方法[31]。该方法将N个时间所对应的目标场(速度或者压力)信息构建为矩阵,求解目标场的自协方差矩阵的特征值,则可求出POD的最优分解。

2 结果与讨论

2.1 合成双射流冲击平板流场特性

为了清晰地认知合成双射流冲击平板的流场特性,图3给出了其时间平均流场。由图可见,整个流场大致关于激励器两出口中心轴线呈反对称分布,两侧均存在一对正负涡。激励器两出口(x为-0.004 5～-0.002 5 m和0.002 5～0.004 5 m)附近的流动比较明显,由于两出口同时产生的是相位差为180°的两股射流,即一出口“吹程”对应另一出口的“吸程”,则射流会受到卷吸作用而向两出口中心轴线偏移。

图3 合成双射流冲击平板流动时间平均流场结果Fig.3 Time-average flow field of dual synthetic jets impinging plate

图4给出了合成双射流冲击平板流动一个周期内8个不同相位下的流线与展向涡量ωz云图,图中以激励器左出口的涡结构刚从出口处脱落的时刻作为0°相位。在相位Φ=0°时,激励器振动膜向左收缩,左侧腔体内的流体受到压缩形成射流,射流在出口剪切作用下形成涡对并脱落;随着时间的推移,射流涡逐渐向下游迁移(Φ=45°),同时受激励器右侧出口的“吸程”作用,整体涡结构略向右侧倾斜;在Φ=90°时激励器左侧腔体处于最大压缩状态,对腔内流体的压缩量也达到最大;此后振动膜逐渐向右运动,左侧腔体压缩量逐渐减小但仍处于压缩状态,射流涡持续向下游迁移(Φ=135°);在Φ=180°相位点,振动膜回到中间位置并开始向右收缩,激励器右侧腔体开始压缩,射流涡即将在右侧出口脱落;射流涡将经历与左侧射流涡同样的过程(Φ=225°,Φ=270°,Φ=315°),而左侧射流涡在自诱导速度下进一步向下游发展并在平板附近分解为许多小结构涡,此时合成双射流冲击平板完成一个周期。

(a) Φ=0°

由图4还可见,在远离合成双射流冲击平板核心区左右两侧分别存在一个稳定的涡,而核心区的涡结构极其复杂,不同结构大小的射流涡相互掺混,并随着时间不断迁移;且在合成双射流激励器工作过程中,核心区平板附近不断有涡的融合与破碎,加强了气流的掺混,对于强化换热十分有利。

2.2 拉格朗日流场结构

对计算得到的合成双射流冲击平板等时间间隔速度场计算向后向积分的FTLE,积分时长为T=0.01 s,时间间隔为t=0.02 ms,得到FTLE场随时间的变化。

将激励器一个工作周期分为8个相位点,如图5所示是FTLE场在一个周期内随时间的变化情况。同样取激励器左出口(x为-0.004 5～-0.002 5 m)射流涡即将脱落时刻为0°相位,右出口(x为0.002 5～0.004 5 m)射流涡即将脱落时刻为180°相位。在Φ=0°时,激励器振动膜开始向左侧运动,其左侧腔体开始处于压缩状态,射流涡从左出口脱落;随着振动膜继续向左运动,激励器左侧腔体压缩量持续增大,射流涡继续向下游迁移(Φ=45°);在Φ=90°时刻,激励器左侧腔体处于最大压缩状态,同样右侧腔体处于最大扩张状态,受激励器右出口的“吸程”影响,射流向右侧偏移;此后,振动膜开始右侧运动,直到位于激励器腔体中心(Φ=180°),在此过程中射流涡继续向下游迁移并略向右侧偏移;在Φ=180°时刻,激励器右侧出口开始有射流涡脱落,射流涡经历与左侧射流涡相同的经历(Φ=225°,Φ=270°,Φ=315°),不断向下游迁移的过程中略向左侧偏移;而此前激励器左侧出口产生的射流涡在此过程中,继续向平板迁移,不断发生融合与破碎。

(a) Φ=0°

在激励器一个工作周期内,前半周期和后半周期,在远离其射流核心区的左右区域,皆存在一个较稳定的涡系结构。射流核心区的涡系结构较为复杂,尤其在平板附近,涡的迁移与融合破碎较为频繁。在激励器的工作过程中,激励器出口产生的射流涡不断迁移到平板附近,射流核心区平板附近的涡则向两侧迁移,破碎融入两侧的稳定涡系中。由图5可见,合成双射流所产生的涡系结构丰富,变化迅速,流场中的气流存在大量的掺混与快速运动。

上述结果表明,图5利用拉格朗日体系下的FTLE方法对流场涡系结构的识别能与图4中的涡量云图很好地对应,证实了FTLE识别方法的客观性与识别精度。

2.3 POD结果分析

为了更深入地了解合成双射流冲击平板流场的结构特征,采用POD方法对流场进行模态分解,将激励器的一个工作周期等分为500个相位点,取三个工作周期作为原始流场数据,共1 500个流场“快照”数据。

图6是前30阶模态的相对能量和累积能量结果。模态1占据了总能量的35%,前6阶模态占据大约80%的能量,剩余较少的能量分散在其余的模态中。这是由于POD是基于特征值大小(也就是能量高低)对分解得到的模态进行降阶排序,大尺度的旋涡结构能量较高,对应的特征值大,也就占据了前几阶模态,小尺度的、破碎的涡结构则表征在后面的低阶模态中。

(a) 相对能量(a) Relative mode energy

由于前6阶模态的能量贡献率达到了80%左右,图7给出了流场经过POD分解后的前6阶模态。结果表明,模态1关于激励器两出口中心轴线大致呈对称分布,主要展现了由两股射流相互卷吸干扰在流场中形成的大尺度涡结构;模态2和模态3则主要展现了激励器两出口附近的涡结构,且左右出口分别交叉存在一对正负涡,关于出口中心轴线呈反对称分布,这是由于激励器左右出口工作相位差为180°,“吹程”和“吸程”交叉进行;模态4和模态6则主要体现了流场中结构相对较小的涡结构,且与模态1大致相似,关于两出口中心轴线大致呈对称分布;模态5展现了射流经出口喷射出去后继续向下游迁移发展的现象。

(a) 模态1(a) Mode 1

POD分解的本质就是对不同特征流动依据其含能大小进行降阶排序,而涡结构的能量往往与涡结构的空间尺度呈正相关。根据图7模态分解的结果来看,模态1和模态2中的涡结构尺度较大且相对集中;而随着射流向下游迁移发展并与壁面发生冲击作用,大尺度的涡结构由于黏性耗散和冲击作用逐渐破碎分解为能量较低的小尺度涡结构,且分布更加分散,其随着能量水平的减小逐步呈现在后几阶模态中。

图8给出了前6阶模态的时间系数功率谱。前2阶模态时间系数的功率谱密度值比后4阶要高一个数量级,前2阶体现了流场的主要特征,后4阶是对流场的细节补充。模态1的时间系数的频率为50 Hz,是激励器工作频率的一半,这是因为模态1主要体现流场中的大尺度涡结构,冲击平板距离激励器出口较近,受冲击平板的影响,流场中的大尺度涡结构的变化没有随激励器的工作周期瞬时响应,而是变化相对较慢,前一个工作周期射流的影响持续时间较长,而后4阶模态则出现了二阶谐频。

图8 前6阶模态时间系数功率谱Fig.8 Power spectrum of the first 6 order time coefficients

根据POD分解结果,合成双射流冲击平板流场的大尺度含能涡系结构主要是在两股射流卷吸与相互作用较强的射流核心区以及激励器出口附近,流场中涡结构大致关于激励器出口中心轴线呈对称或者反对称分布,表明合成双射流的流场具有高度的对称性。

3 结论

采用大涡模拟方法对合成双射流冲击平板流动进行了数值模拟,基于FTLE方法对流场的拉格朗日涡结构进行了识别,并对流场开展了本征正交分析,得出以下结论:

1)在合成双射流冲击平板过程中,激励器左右出口进行周期性“吸吹”工作过程,射流核心区尤其平板附近涡结构较为复杂,左右两侧远离核心区域存在稳定的尺度较大的涡结构,核心区复杂的涡结构与已发表研究结果[32]一致;

2)利用FTLE方法识别的拉格朗日涡结构能较清晰地反映流场的涡系结构特性与演化特性,且能与涡量较好对应;

3)POD分析结果表明,模态1占据总能量的35%,前6阶模态占据了大约总能量的80%,模态阶数越高,所表征的流动结构尺度越小;

4)POD分解结果的前6阶模态所反映的流场特性,表明合成双射流冲击平板的流场具有高度的对称性,该结果与已发表研究[33]一致。

研究结果进一步加深了对合成双射流冲击平板流动特性的了解,为下一步开展散热特性以及合成喷雾流场特性的研究奠定了基础。