非对称凹腔橫向阵列压电风扇强化冷却特性研究

2022-08-30张冬冬谭晓茗张靖周鹿世化李鑫郡

南京航空航天大学学报 2022年4期

张冬冬，谭晓茗，张靖周，鹿世化，李鑫郡,3

（1.南京师范大学能源与机械工程学院，南京 210042；2.南京航空航天大学能源与动力学院，南京 210016；3.南京航空航天大学航空动力系统重点实验室，南京 210016）

压电风扇作为一种新型的强化换热装置，一般主要由压电材料和柔性叶片组成［1］。由于压电材料本身所固有的逆压电效应，当对其施加交流电压时，电能能够持续高效地转化为机械能。在交流电压作用下，压电材料以与该输入电信号相同的频率交替膨胀与收缩，并驱动附着其上的柔性叶片同频振荡。随着柔性叶片的振荡，将在自由端产生伪射流并导致较高局部对流传热能力［2］。此外，由于具有结构简单、可靠性高、功耗低和噪声小等突出优点，压电风扇作为一种潜在的强化换热器件在电子冷却应用［3-4］、可再生能源［5］以及能量收集［6］等方面均得到了广泛的研究。

早期研究者主要致力于揭示压电风扇激励的流动特性。大量研究表明，压电风扇激励的局部流动具有涡旋特征，且最大射流速度与振动参数（振动频率、振幅和振型）［7-8］、几何参数［9-10］以及运行环境［11-13］相关。特别是Ebrahimi 等［14］成功捕捉到了从压电风扇后缘和侧缘位置形成的瞬态3D 涡旋结构，及其向下游脱落和射流形成过程。基于此，许多研究者则将压电风扇潜在的冲击强化换热性能作为研究重点［15-17］。例如，Kimber 等［15］量化了单一参数（包括频率、振幅和几何形状）的影响及其对局部传热性能的相对影响。结果表明，受热平面上局部传热系数的分布将从扇尖到表面间隙较小的叶状转变为中间间隙处的近圆形，然后在较大间隙处进一步变为椭圆形。此外，存在特定间隙使得传热效率最佳。Liu 等［16］研究了压电风扇排布方式对传热影响，发现垂直和水平布置的换热系数量级相同，但水平热源布置显著优于垂直情况。此外，Zhou 等［17］采用粒子图像测速技术和感温涂料技术研究证实压电风扇在小间隙和高雷诺数下同样具有强化换热的能力。

然而，随着研究的深入，仅依靠单风扇强化平面的散热方式已达上限，难以通过参数优化使其效率进一步提升［18］。因此，为了进一步探索提升使用压电风扇散热的可能性，单个压电风扇结合被动散热策略（各类扩展表面，如散热器、圆柱形受热面和异形表面）［19-20］、耦合多压电风扇［21］和耦合多风扇与被动策略相结合的方式开始引起人们的关注［22-25］。Ma 等［22］以及Li 和Wu［23］分别研究了由双压电风扇冷却的板翅式和针翅式散热器的传热。在他们的实验中，双压电风扇垂直安装在散热器的正上方，阐述了压电风扇的相位差、配置和高度以及散热器尺寸对散热器热性能的影响。此外，Sufian 和Abdullah［24］设计了一种更紧凑的组合结构，通过将多个振动梁嵌入到翅片阵列中来增强翅片阵列传热并减少冷却器体积。Abdullah 等［25］通过使用实验设计（Design of experiment，DOE）方法，对3 个压电风扇的叶尖间隙和定向角的优化进行了研究，以最大限度地提高翅片散热器的散热性能。

从上述讨论看出，尽管目前国内外学者针对压电风扇换热特性开展了大量研究，但应用对象主要集中在具有对称外形的热源，针对非对称热源的研究还非常缺乏，特别是针对非对称凹腔这类经常出现在涡轮叶片、电子设备散热等场景中的冷却结构目前基本没有涉及［26-27］。因此，本文主要针对利用多压电风扇系统冷却非对称凹腔问题，开展非稳态数值研究，重点关注振动相位、相对曲率等因素对腔内流动换热特性的影响。相关研究成果对于指导非对称凹腔型结构冷却系统设计具有一定的参考意义。

1 计算模型与数据处理方法

1.1 计算模型

本文所使用的特定压电风扇如图1（a）所示，主要包括压电陶瓷片和不锈钢膜片两部分。其中压电陶瓷片尺寸为18 mm×8 mm×1 mm，不锈钢膜片尺寸为28 mm×8 mm×1 mm，一阶共振频率为67.3 Hz。图1（b）所示为压电风扇简化计算模型，其长度LPF为46 mm，宽度WPF为8 mm。叶尖振幅AP设为4 mm，叶尖沿振动方向最大位移APP为8 mm。

图1 压电风扇结构示意图Fig.1 Schematic of the piezoelectric fan

为了研究多风扇系统对非对称凹腔的冷却换热性能，本文构建了如图2（a）所示的三维计算模型。该模型由一个类似于涡轮叶片前缘的非对称凹腔和垂直排列在其顶部的多压电风扇系统组成。模型长度与高度分别为120 mm 和48 mm，宽度则随着凹腔曲率而变化。多压电风扇系统沿展向垂直安装于凹腔中心线上方，并通过固定底部边缘（x=0）位置来限制其仅沿y方向振动。非对称凹面与压电风扇叶尖间隙高度G固定为2 mm，多压电风扇间距P为2 mm。图2（b）为非对称凹腔剖面结构。其中右侧半圆半径RR恒为48 mm，左侧半圆半径RL分别取8 mm、16 mm 和24 mm，于是得到3 种非对称凹腔结构。表1 使用分段函数精确描述了凹腔的剖面几何形状。

表1 非对称凹面轮廓Table 1 Profiles of asymmetrical concave surfaces

图2 计算模型Fig.2 Basic geometry of computational domainst

为了定量描述非对称凹腔的几何特征，定义两个无量纲参数：基于压电风扇振动尺度的半圆形表面的曲率K以及两侧半圆形表面的相对曲率Kr，公式如下

式中：R为圆弧的半径；DPF为压电风扇叶尖振动包络区特征长度；AP和WPF为压电风扇叶尖振幅与宽度，如图1 所示。

此外，计算域的边界条件如下：压电风扇简化为无厚度绝热壁面，其振动轨迹由用户自定义函数（User defined function，UDF）依据测振结果定义；非对称凹面为无滑移边界，其恒热流密度为q=1 600 W/m2；其余边界均设为压力边界条件。

1.2 计算方法

压电风扇作为周期性运动器件，其运动轨迹可由时间-位移函数描述。依据压电风扇一阶谐振频率下振动测试结果，文中压电风扇运动轨迹方程如下

式中：x表示叶片到固定端距离；t表示压电风扇运行时间；Y（x）表示压电风扇某处位移峰值。式（5）中的相关系数为：p1=-1.856×10-7；p2=2.347×10-4；p3=-7.25×10-3；p4=6.281×10-2；p5=-5.124×10-2。

本文中计算模型的建立基于以下假设：（1）流动不可压缩，具有温度相关特性；（2）重力影响忽略不计；（3）辐射效应被忽略。因此，选择3D 非定常雷诺平均Navier-Stokes 方程（RANS）作为控制方程。

式中：ρ为流体密度；p为流体压力；τij为黏性应力；gi为i方向的重力加速度。

式中：cp为空气比热容；T为温度；k为热导率。

此外，湍流模型选择SST（Shear stress transport）k-ω两方程模型，并选用SIMPLEC 算法进行压力速度耦合，其中压力、动量、湍动能耗散率和能量采用二阶迎风格式进行离散。

1.3 网格独立性验证

由于采用动态网格技术描述压电风扇振动过程时，需要在每个计算时间步中对压电风扇叶尖处的网格进行局部重构和加密处理，因此消除系统的网格敏感性至关重要。选择非对称凹腔中心线（C=0）处的时均局部对流传热系数作为关键参数。在一个完整的周期内，时均对流传热系数可以使用如下表达式进行计算。

式中：t和dt分别为每个周期的时间和积分时间步长；h为壁面瞬时对流换热系数，其计算公式如下

式中：q、Tw和Ta分别为加热表面热流密度、瞬时壁面温度和瞬时凹腔内平均空气温度。

图3 是以多风扇反相振动为例，通过数值模拟得到的网格数对非对称凹腔中心线（C=0）处时均局部对流传热系数的影响情况。可以看出，当网格数由150 万增加至210 万时，其时均对流换热系数的最大差异仅在1.5%以内，因此最终非对称凹腔模型的计算网格数选取为150 万。

图3 网格独立性验证Fig.3 Grid independence verification

1.4 实验验证

为了评估本文数值计算结果的可靠性，搭建如图4（a）所示的验证性实验平台。多风扇系统固定于坐标架之上，其中叶尖到凹腔的间隙可通过坐标架进行调节。利用信号发生器和功率放大器为多风扇系统提供激励信号，并通过输出信号的调节使得风扇振幅为4 mm。其中，风扇叶尖振幅的测量采用LK-G3000 激光位移传感器。非对称凹腔的轮廓由胶木块制成，并将厚度仅为0.05 mm 的加热箔膜黏附于凹腔内壁。箔膜由直流电源提供加热热流。热流均匀性则依靠箔膜两端的条形铜片保证。如图4（b）所示，共有9 个K 型热电偶用于监控箔膜下方关键位置瞬时温度信息，同时监测腔内瞬时环境温度。K 型热电偶所测得的瞬时温度信息由数据采集卡收集，采集频率为10 Hz。其中，由于加热箔膜非常薄，因此每次响应时间均不超过10 s，平均每组实验总采样时间不超过15 s。

图4 实验装置示意图Fig.4 Schematic diagram of experimental device

图5 所示为实验值和模拟值在凹面中心线处局部对流传热系数分布的比较情况。可以看出，模拟结果与实验结果虽然具有差异，但整体偏差并不显著，最大偏差仅为7.1%。这些小的偏差很可能是由于数值模型中辐射的遗漏和实验中通过胶木块的热传导所引起。此外，实验中的环境空气温度是由单一热电偶测得，而在模拟中则使用域内的平均流体温度。

图5 沿S 方向实验与数值结果对比Fig.5 Comparison of numerical and experimental hav at centerline along S-direction

1.5 数据处理方法

压电风扇激励强化传热的原理本质上与冲击射流相似，是形成的局部扰动涡结构对热边界层的破坏，进而增强对流换热性能。因此可通过识别风扇振动周期内涡结构的演化过程，来分析三维流场的运动情况。本文采用λ2判据［28］进行涡识别，将流场速度梯度张量J分解为应变张量S和旋转张量Ω。

该判据通过计算两者组合张量S2+Ω2的3 个特征值（λ1≥λ2≥λ3），认为压力达到截面最小的充要条件为λ2＜0，其中λ2＜0 的点即属于涡核空间位置。

为了研究非对称凹面沿展向S和弦向C的平均对流传热系数分布规律，本文分别定义了沿展向和弦向的两个横向平均对流换热系数，其表达式如下

式中L为积分的宽度。沿展向和弦向L分别取风扇最大位移1APP或宽度1WPF。

2 结果与分析

图6 和图7 分别为压电风扇同相以及反相振动时非对称凹腔表面时均对流换热系数分布情况。图中黑色虚框表示压电风扇的振动包络区，黑色实线表示压电风扇的平衡位置，红色虚线表示凹腔上相邻曲面边界。整体看来，时均对流传热系数分布的不对称性主要是由于翼展方向相邻风扇的相互作用和弦向表面不对称共同造成的。从图6 和图7可以看出，Kr与φ对于时均对流传热系数的分布特征具有明显的耦合作用。

对于同相振动（φ=0°）风扇，如图6（a）所示，当非对称凹腔两侧相对曲率较大时（Kr=6），由于侧壁对于腔内流动限制作用差异明显，因此表面时均对流换热系数呈现出明显的非对称分布特征。而且相邻风扇间隙区域内的换热能力得到显著强化，这是由于同相振动（φ=0°）增强了间隙区域内的局部扰动所致。而当相对曲率Kr的值从6 降至2 时，如图6（c）所示，可以看出此时时均换热系数分布的不对称性几乎完全消失，且振动包络以及相邻风扇间隙内的传热均得到显著提升。这说明相对曲率对于非对称凹腔内的扰动具有抑制作用。

图6 压电风扇同相振动时的时均换热系数分布(φ=0°)Fig.6 Distribution of time-averaged local convective heat transfer coefficient on asymmetrical concave surfaces for vibrating in-phase (φ=0°)

图7 给出了多风扇反相振动时（φ=180°）相应的局部时均传热系数分布情况。与同相振动（φ=0°）相比，反相振动的多风扇系统（φ=180°）的强化换热区域主要集中在包络区内，而在相邻风扇间隙中的作用并不明显。这是由于反相振动的相邻风扇激励的局部流动方向始终相反，因而局部扰动被抑制所致。与同相振动（φ=0°）相同的是，随着凹腔相对曲率（Kr）的降低，反相振动时（φ=180°）的换热系数分布不对称性逐渐消失，换热能力同样得到一定程度的提高。

图7 压电风扇反相振动时的时均换热系数分布(φ=180°)Fig.7 Distribution of time-averaged local convective heat transfer coefficient on asymmetrical concave surfaces for vibrating out-of-phase (φ=180°)

为了进一步定量分析不同振动相位（φ）和相对曲率（Kr）对于多风扇系统冷却非对称凹腔的影响，图8 和图9 分别给出了同相及反相振动时非对称凹腔表面积分平均对流传热系数沿展向和弦向的分布情况。图9 中黑色实线代表压电风扇侧缘位置，黑色虚线代表压电风扇振动的最大偏移位置。

图8 同相振动时的积分平均对流传热系数分布(φ=0°)Fig.8 Distribution of laterally-averaged convective heat transfer coefficients for vibrating in-phase (φ=0°)

图9 反相振动时的积分平均对流传热系数分布(φ=180°)Fig.9 Distribution of laterally-averaged convective heat transfer coefficients for vibrating out-of-phase (φ=180°)

如图8（a）所示，当多风扇系统同相振动时，展向分布的峰值havC出现在相邻风扇的间隙区域，并且在越过包络区后迅速降低，这是因为凹腔在展向上对流动阻碍非常小。具有最小相对曲率（Kr=2）的表面产生最高的havC值，其中Kr=6 和Kr=2时凹腔的峰值之间的差异约为30%。类似的规律也出现在沿弦向分布中，如图8（b）所示。然而，不同之处在于沿弦向分布时峰值havS则几乎位于包络区中心。当多风扇系统反相振动时，如图9（a）所示，沿展向分布的峰值havC均匀分布于各个风扇包络区内，且峰值大小随着相对曲率的增加而迅速降低。此时Kr=6 和Kr=2 的表面峰值之差大于50%。弦向分布（图9（b））也具有相似的特征。

综上，可以看出振动相位主要影响峰值出现的位置，而相对曲率则主要决定峰值的大小。

为了揭示多风扇系统所激励的非对称凹腔内瞬时流动特征，以Case 2（Kr=3）为例，图10 和图11 分别给出了同相和反相振动时一个周期内4 种典型风扇位置时λ2=-3×104的瞬态等值面涡结构的演化过程以及受热凹腔瞬时温度分布情况。

图10 多风扇同相振动的瞬时温度和λ2=-3×104瞬态等值面(Case 2)Fig.10 Instantaneous temperature and iso-surface with λ2=-3×104 for multi-piezoelectric fans in-phase vibration (Case 2)

如图10（a）和图10（c）所示，当同相振动的多风扇系统运动至其中间位置时，风扇振动速度达到最大值，此时由于压电风扇叶片侧缘及叶尖对周围空气具有强烈的剪切作用，因此在叶尖及两侧侧缘形成明显的涡结构，并且由于临近风扇同相振动促进了风扇间隙的扰动，进而形成了较强的间隙涡；而当多风扇系统运动至最大位置时，如图10（b）和图10（d）所示，由于此时风扇振动速度减小至0，由于惯性作用使得涡结构从叶片脱落，形成局部射流撞击凹腔表面，进而降低了叶片附近的凹腔表面温度。

对比图11（a）和图10（a），可以发现反相振动的多压电风扇系统在运行至平衡位置时同样会在叶尖及两侧缘形成较为明显的涡结构，但与同相振动相比，其尺度有所降低，而且由于临近风扇运动方向相反此时间隙涡完全消失。这一现象可以解释不同相位时（φ=0°和180°）多风扇间隙区域内的时均换热系数分布特点。对比图11（d）和图10（d），可以看出反相振动时形成的脱落涡尺度明显小于同相振动的情况，因而同相振动的整体换热效率高于反相振动情况。

图11 多风扇反相振动的瞬时温度和λ2=-3×104瞬态等值面(Case 2)Fig.11 Instantaneous temperature and iso-surface with λ2=-3×104 for multi-piezoelectric fans out-of-phase vibration (Case 2)

图12 所示为风扇运行至0 位置时凹腔Case 1（Kr=6）与Case 3（Kr=2）中λ2=-3×104的瞬态等值面涡结构以及受热凹腔瞬时温度分布情况。对比图12（a）、图10（a）与图12（c），可以看出随着凹腔Kr值的减小，同相振动压电风扇所激励的间隙涡具有向中间位置风扇振动包络区内逐渐收缩的趋势，这表明间隙涡的收缩效应可能是引起中间位置压电风扇包络区内换热系数上升的主要原因；而对比图12（b）、图11（a）与图12（d），同样可以发现，压电风扇反相振动时激励脱落涡的尺度会随着Kr值的下降而逐渐上升，同时涡系结构向各个风扇的聚集程度也会得到强化，因此包络区内时均换热系数随着Kr值的下降而逐渐增强。

图12 多风扇同相、反相振动至0 相位瞬时温度和λ2=-3×104瞬态等值面(Case 1 与Case 3)Fig.12 Instantaneous temperature and iso-surface with λ2=-3×104 for multi-piezoelectric fans in-phase and out-of-phase vibration under 0 (Case 1 and Case 3)