存在横流时双压电风扇激励传热特性

李鑫郡，张靖周，谭晓茗

南京航空航天大学 能源与动力学院，南京 210016

对于压电风扇激励流动和传热性能的研究大多以单个压电风扇作为研究对象。Kim等[4-5]采用PIV(Particle Image Velocimetry)锁相测试和烟迹显示方法对自由空间压电风扇激励的瞬时流场特征进行了研究，揭示出在每一个振动周期，振动膜片叶尖附近形成一对相互逆转的大尺度涡结构，并且在这两个相互逆转的涡之间形成高速流动；Kimber等[6]对压电风扇诱导的压力场和拍动射流流动速率进行了实验测试，研究表明可达流动速度与膜片叶尖速度近似呈现抛物线分布关系，当膜片经历其上下最大振幅位置时上游泵吸能力加强；谭蕾[7-8]和孔岳[9]等采用数值模拟方法研究了压电风扇诱导的非定常流场特征。Acikalin[10-12]和Kimber[13-14]等开展了一系列的压电风扇传热研究，较为系统地揭示了激励参数对于换热特性的影响；Liu等[15]实验研究了水平和垂直于壁面放置的单个压电风扇传热特性；Lin[16]基于振动测试获得了谐振风扇位移函数，运用动网格技术研究了压电风扇激励的非定常射流冲击壁面的流动和换热瞬时特征，之后采用同样的方法研究了强迫流动环境中圆柱曲面上压电风扇的作用效果[17]，研究发现在强迫对流通道中，压电风扇激励射流与横流外掠圆柱形成的尾迹流相互作用，能够形成改善圆柱表面对流换热的效果，但在某些情形下也可能削弱圆柱表面对流换热，与主流和压电风扇射流特征速度密切关联；Jeng和Liu[18]实验研究了轴向流通道中热沉上游安置压电风扇的传热特性，研究表明，在通道低雷诺数流动条件下，压电风扇的谐振可以增强主流的湍流强度，起到改善热沉传热能力的作用。

本文关注双压电风扇系统在存在横流环境中的传热特性及其影响机制。基于横流环境中压电风扇振动测试获得风扇位移规律，结合动网格技术数值模拟横流环境中压电风扇激励的三维非定常流场。本文研究的双压电风扇采用并列和串列两种方式，横流速度与压电风扇的特征速度相当，以期揭示双风扇系统在横流环境中的流动传热机理。

1 计算物理模型

1.1 压电风扇

本文研究以某个特定结构的压电风扇激励器作为研究对象，该型压电风扇主要由压电陶瓷片(PZT)、不锈钢膜片以及底座组成，如图1所示。压电陶瓷片单侧粘附于不锈钢膜片，其长×宽×厚(Lp×W×tp)为28.0 mm×25.0 mm×0.5 mm。柔性膜片伸出长度×厚度(Lb×tb)为38.0 mm×0.1 mm。记风扇叶尖沿振动方向处于前后极限位置时的最大位移为App，该位移是叶尖振幅Ap的2倍。

图1 压电风扇结构示意图Fig.1 Geometry schematic of piezoelectric fan

1.2 计算域和边界条件

图2 数值计算模型示意图Fig.2 Schematic of numerical computational model

为了研究双压电风扇系统在横流环境中，排布方式和振动相位差φ对于流场和加热壁面换热能力的影响，分别建立如图2(a)和图2(b)所示的双风扇并列、双风扇串列布置的物理模型，对于并列布置模式，两压电风扇之间的无量纲间距设为P/W=1.2 (其中P为两风扇中心线之间的距离)；对于串列布置模式，两压电风扇之间的无量纲间距也设为P/W=1.2。同时为了对比分析，也建立相同尺寸的单风扇物理模型，如图2(c)所示。根据已有的压电风扇传热研究[13-14]，压电风扇具有局部传热强化的特征，压电风扇叶尖距离加热壁面的间距对于流动和传热状况存在重要的影响，过小的间距由于会引发风扇流和壁面的相互作用而导致风扇振幅衰减，因此本文选择的风扇叶尖距离加热壁面的间距G=3 mm。在三维非定常流场数值计算中，将压电风扇简化处理为均匀的振动膜片，其运动按照特定的位移函数设定。加热表面(Heat foil)设置为无滑移速度边界和恒热流密度(q=1 700 W/m2)热边界条件，外掠加热表面横流入口设为均匀速度进口，横流出口设为压力出口，其余边界均设为压力边界条件。根据已有的针对单个压电风扇振动特性的测试研究[19]，该型压电风扇在一阶谐振频率下振动的叶片位移最大速度或特征速度约为2.3 m/s，因此，本文计算模型中的横流入口速度设定为uCF=2.3 m/s。

2 计算方法

2.1 压电风扇位移函数

压电风扇的类悬臂梁振动过程所激发的流动过程是典型的非稳态运动，且射流速度较低，因此采用非定常不可压缩Navier-Stokes方程模拟。将压电风扇柔性膜片的运动轨迹y′(z,t)使用如下振型函数描述：

y′(z,t)=Y(z)sin(2πft)

(1)

式中：Y(z)为风扇膜片位移；f为振动频率；t为模拟运动时间。压电风扇的位移函数采用文献[19]的测试拟合多项式，即

Y(z)=p1z4+p2z3+p3z2+p4z+p5

(2)

式(2)的单位为mm。多项式系数依次为p1=-1.886×10-6,p2=2.447×10-4,p3=-7.31×10-3,p4=6.078×10-2,p5=-5.092×10-2。

2.2 计算网格及独立性验证

进行压电风扇三维非定常流场数值模拟时，将压电风扇位移函数通过用户自定义函数(UDF)输入，作为压电风扇周期性运动条件，每个周期划分为200个时间步长。计算过程中使用动网格技术描述压电风扇的运动，利用局部重构方法实现计算过程中网格再生成[19]。计算区域采用四面体非结构网格，在压电风扇叶尖附近振动位移最剧烈的区域以及换热表面附近进行适当的网格加密及自适应处理设置。图3所示为压电风扇一个周期内最大位移时刻的局部网格示意图。

计算网格数经过网格独立性验证确定。图4是以单风扇模型为例，通过数值模拟得到的计算网格数对加热表面时均对流换热系数的影响情况。计算过程中通过监测每个循环中的壁面平均温度随时间的变化过程，可以发现在经历约20个周期以后，壁面平均温度随时间呈周期性变化。选择一个周期进行分析，时均对流换热系数采用瞬时值积分平均确定:

(3)

式中：Δt为每个周期的时间；dt为时间步长。

可以看出，当计算网格数在100万和160万之间时，时均对流换热系数的差异仅在3%以内，因此最终单风扇模型的计算网格数选取为132万。运用同样的方式，确定双风扇模型的计算网格数为160万。

图3 不同相位时局部网格重构Fig.3 Re-meshing of local grids at different phases

图4 网格独立性验证Fig.4 Verification of grid independency

2.3 计算方法和验证

应用三维非定常雷诺时均控制方程(RANS)进行数值模拟。湍流模型选择SST(Shear Stress Transport)k-ω两方程湍流模型，该湍流模型被众多的研究者应用于压电风扇流场计算[16-17, 19]。

数值模拟采用Fluent软件，在计算过程中需要采用欠松弛加以控制以确保非定常计算过程的稳定性，当3个连续的风扇振动周期内壁面温度的相对误差小于10-4时，认为数值解达到准稳定的周期性状态[16-17]。

图5 无横流单个风扇时均局部对流换热系数分布Fig.5 Distributions of time-averaged local convective heat transfer coefficients for single fan without cross flow

为验证本文的数值计算方法，选用Kimber等[13]的单个压电风扇在无横流影响下的表面对流换热特性实验进行验证。图5为不同风扇叶尖距离加热壁面的间距下数值计算得到的时均局部对流换热系数分布情况，图中虚线框代表风扇振动时所扫过的区域，实线代表风扇停振时的位置。可以看到压电风扇作用下的加热表面对流换热系数在叶尖包络区两侧呈现哑铃状分布，与Kimber等[13]实验揭示的分布特征相符。

图6为在不同风扇叶尖距离加热壁面的间距下数值计算得到的驻点努塞尔数与Kimber等[13]的实验数据对比。此处的压电风扇雷诺数和努塞尔数分别定义为

(4)

式中：uPF为压电风扇的特征速度；DPF为风扇振动叶尖包络区的特征长度；k为导热系数；ν为运动黏度。uPF、DPF的表达式分别为

(5)

比较表明，本文数值计算结果与Kimber等[13]的实验数据相吻合。

图6 驻点努塞尔数随无量纲距离G/Ap的变化Fig.6 Stagnation Nusselt number vs non-dimensional gap G/Ap

3 结果与分析

3.1 瞬态流场特征

压电风扇振动诱导一系列非定常的拍动涡环并不断地聚合，从而形成风力相对集中的“准连续”射流。为了揭示压电风扇射流与横流相干的涡系结构，借鉴Jeong和Hussain[20]提出的λ2涡识别判据进行分析。该判据是将流场的速度梯度张量J分解为对称部分S(应变率张量)和非对称部分Ω(旋转张量)，并通过计算组合张量S2+Ω2的3个特征值(λ1≥λ2≥λ3)，认为压力达到截面最小的充要条件为λ2<0，其中λ2为负值的点即属于涡核空间位置。J、S、Ω三者之间的关系式为

(6)

图7 单风扇不同振动相位瞬时温度和λ2=-3×104 瞬态等值面Fig.7 Instantaneous temperature and iso-surface for λ2=-3×104 at different vibrating phases for single fan

图7为单个压电风扇在一个运动周期内4种典型振动相位下λ2=-3×104的瞬态等值面发展演变情况。当压电风扇运动至其平衡位置时，如图7(b)和图7(d)所示，其叶尖位移速度达到最大值，依靠风扇叶片对气流的剪切机制在叶尖以及叶片两侧缘形成较为强烈的马蹄涡，而当风扇运动至其最大振幅位置时，由于风扇位移速度逐渐减小，涡结构显著减弱，此时的横流效应相应增强，挟带叶尖脱落的涡向下游迁移。从叶尖脱落的涡对壁面形成冲击作用，从而导致压电风扇振动的包络区域附近的局部对流换热强化；由于周期性的风扇振动，叶尖脱落涡在包络区下游形成周期性扫掠，同时流经叶尖包络区的横流也受到压电风扇的脉动激励，因而在压电风扇下游形成温度相对较低的条带状区域。

图8 双风扇并列反相振动的瞬时温度和λ2=-3×104 瞬态等值面Fig.8 Instantaneous temperature and iso-surface for λ2=-3×104 for dual out-of-phase fans in face-to-face mode

图8和图9分别为双风扇并列模式下反相和同相振动时λ2=-3×104瞬态等值面以及加热壁面瞬时温度分布情况。与单个风扇相比，由于双风扇的作用区域增大，因此加热表面的冷却效果得到增强。同时相邻风扇间的振动对于涡系的演变也具有一定的影响，对于并列布置的双压电风扇，振动相位引起的相邻风扇之间涡结构演化特征具有较大的差异。并列双风扇反相振动时，如图8所示，相邻两个风扇相向运动至各自的平衡位置时，叶尖冲击脱落涡之间形成融合和发展，并对其他时刻的涡系向下游的发展也产生影响，从而在压电风扇下游形成温度更低、范围更宽的影响区域。而对于并列双风扇同相振动，虽然存在着相邻风扇间的相互作用，但每个风扇诱导的涡结构与单个风扇较为接近。

图10和图11分别为双风扇串列模式下反相和同相振动时的涡核瞬态等值面以及加热壁面的瞬时温度分布情况。由于两个风扇串列于横流中，上游的压电风扇诱导的涡在横流作用下向下游迁移，与下游压电风扇激励作用形成耦合，因此，相对于单一风扇，串列双压电风扇叶尖脱落的涡对下游表面对流换热影响作用的区域更长。对比双风扇串列模式下反相和同相振动时的涡结构演变，可以看出，同相振动情形下相邻风扇的脱落涡之间的融合和发展更具优势。

图9 双风扇并列同相振动的瞬时温度和λ2=-3×104瞬态等值面Fig.9 Instantaneous temperature and iso-surface of λ2=-3×104 for dual in-phase fans in face-to-face mode

图10 双风扇串列反相振动的瞬时温度和λ2=-3×104瞬态等值面Fig.10 Instantaneous temperature and iso-surface for λ2=-3×104 for dual out-of-phase fans in edge-to-edge mode

图11 双风扇串列同相振动的瞬时温度和λ2=-3×104瞬态等值面Fig.11 Instantaneous temperature and iso-surface for λ2=-3×104 for dual in-phase fans in edge-to-edge mode

3.2 传热时均结果分析

图12为不同压电风扇作用下的一个周期内恒热流壁面时均局部对流换热系数分布情况。对于单个风扇，如图12(a)所示，受横流的影响，风扇叶尖扫掠的包络区的对流换热系数分布呈现水滴状，与单个压电风扇在无横流环境中的哑铃形分布不同(图5(a))，显然这是横流对风扇激励射流冲击的影响所致。对于并列双风扇，如图12(b)和图12(c)所示，风扇叶尖扫掠的包络区附近局部对流换热系数依然呈现水滴状分布，在包络区的时均换热能力与单风扇几乎相同，但相对于单风扇，水滴状在流向上有所延伸，这是相邻风扇脱落涡之间的融合和发展所致，注意到风扇之间的干涉作用主要体现在下游区域，反相振动时风扇间隙下游区域的换热能力明显得到强化。对于串列双风扇排布，如图12(d)和图12(e)所示，反相振动时上游风扇包络区的换热能力与单风扇情况基本一致，下游风扇包络区的换热能力有所弱化，但下游风扇尾迹区域的换热能力明显高于单风扇的情况；同相振动时，上游风扇与下游风扇包络区内的换热能力均高于单风扇情况，两风扇间隙区域的换热能力也得到明显加强，综合而言，在横流中串列双风扇同相振动时，局部换热能力可以得到更大程度的强化。

为了进一步对比分析各种排布方式和不同振动相位下壁面平均换热能力的差异，沿横流方向定义侧向积分平均对流换热系数为

(7)

式中：Ly为侧向积分的参考长度。对于单风扇和串列双风扇，Ly分别取App和2App；对于并列双风扇，Ly取为双风扇间距P。

图13为沿横流方向的侧向积分平均对流换热系数。图中竖实线、竖虚线分别表示压电风扇中心线位置、风扇侧缘位置。对于单个风扇，如图13(a)所示， 在单风扇振动包络区， 按照App和2App平均的对流换热系数存在较大的差异，而在单风扇振动包络区下游1.5倍风扇宽度(1.5W)之后，按照App和2App平均的对流换热系数相等，说明压电风扇的影响趋于消失。对于并列双风扇，如图13(b)所示，同相和反相振动对上游以及包络区的对流换热影响并无明显差异，但是在下游1W之后反相振动的平均对流换热系数明显高于同相振动；对于串列双风扇，如图13(c)和图13(d)所示，同相振动始终优于反相振动情况，在反相振动时，相邻两风扇之间的对流换热显著低于同相振动的方式。

图12 有横流时加热表面时均局部对流换热系数分布Fig.12 Distributions of time-averaged local convective heat transfer coefficients on heat foil in presence of cross flow

图13 侧向积分平均对流换热系数沿横流方向的分布Fig.13 Distribution of laterally-averaged convective heat transfer coefficient along cross flow direction

4 结 论

本文采用动网格技术对横流环境中垂直加热壁面的单个压电风扇、并列双风扇以及串列双风扇三维非定常流动和传热特性进行了数值模拟，结果表明：

1) 存在横流时，单个风扇叶尖扫掠包络区的对流换热系数分布呈现水滴状，与无横流影响下的哑铃形分布不同。

2) 对于并列布置的双压电风扇，反相振动有利于相邻风扇脱落涡之间的融合和发展，双风扇振动包络区下游的换热能力优于同相振动情形。

3) 对于串列布置的双压电风扇，上游的压电风扇诱导的涡动射流在横流作用下向下游迁移，与后一个压电风扇激励作用形成耦合。反相振动时相邻两风扇之间局部区域的对流换热显著低于同相振动的情形。

[1] YOO J H, HONG J I, CAO W. Piezoelectric ceramic bimorph coupled to thin metal plate as cooling fan for electronic devices[J]. Sensors Actuators A, 2000, 79(1): 8-12.

[2] CHUNG H C, KUMMARI K L, CROUCHER S J, et al. Development of piezoelectric fans for flapping wing application[J]. Sensors and Actuators A: Physical, 2009, 149(1): 136-142.

[3] GILSON G M, PICKERING S J, HANN D B, et al. Piezoelectric fan cooling: A novel high reliability electric machine thermal management solution[J]. IEEE Transactions on Industrial Electronics, 2013, 60(11): 4841-4851.

[4] KIM Y H, WERELEY S T, CHUN C H. Phase-resolved flow field produced by a vibrating cantilever plate between two endplates[J]. Physics in Fluids, 2004, 16(1): 145-162.

[5] KIM Y H, CHUN C H, WERELEY S. Flow field around a vibrating cantilever: Coherent structure education by continuous wavelet transform and proper orthogonal decomposition[J]. Journal of Fluid Mechanics, 2011, 669: 584-606.

[6] KIMBER M, SUZUKI K, KITSUNAI N, et al. Pressure and flow rate performance of piezoelectric fans[J]. IEEE Transactions on Components and Packaging Technologies, 2009, 32(4): 766-775.

[7] 谭蕾， 谭晓茗， 张靖周. 压电风扇激励非定常流动和换热特性数值研究[J]. 航空学报， 2013, 34(6): 1277-1284.

TAN L, TAN X M, ZHANG J Z. Numerical investigation on unsteady flow and heat transfer characteristics of piezoelectric fan[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(6): 1277-1284 (in chinese).

[8] TAN L, ZHANG J Z, TAN X M. Numerical investigation of convective heat transfer on a vertical surface due to resonating cantilever beam[J]. International Journal of Thermal Sciences, 2014, 80: 93-107.

[9] 孔岳， 李敏， 吴蒙蒙. 压电风扇非定常流场速度分布的数值研究[J]. 工程力学， 2016, 33(1): 209-216.

KONG Y, LI M, WU M M. Numerical investigation on the velocity of unsteady flow field induced by piezoelectric fan[J]. Engineering Mechanics, 2016, 33(1): 209-216 (in Chinese).

[10] ACIKALIN T, WAIT S M, GARIMELLA S V, et al. Experimental investigation of the thermal performance of piezoelectric fans[J]. Heat Transfer Engineering, 2004, 25(1): 4-14.

[11] ACIKALIN T, GARIMELLA S V, RAMAN A, et al. Characterization and optimization of the thermal performance of miniature piezoelectric fans[J]. International Journal of Heat and Fluid Flow, 2007, 28(4): 806-820.

[12] ACIKALIN T, GARIMELLA S V. Analysis and prediction of the thermal performance of piezoelectrically actuated fans[J]. Heat Transfer Engineer, 2009, 30(6): 487-498.

[13] KIMBER M, GARIMELLA S V, RAMAN A. Local heat transfer coefficients induced by piezoelectrically actuated vibrating cantilevers[J]. ASME Journal of Heat Transfer, 2007, 129(9): 1168-1176.

[14] KIMBER M, GARIMELLA S V. Measurement and prediction of the cooling characteristics of a generalized vibrating piezoelectric fan[J]. International Journal of Heat and Mass Transfer, 2009, 52(19-20): 4470-4478.

[15] LIU S F, HUANG R T, SHEU W J, et al. Heat transfer by a piezoelectric fan on a flat surface subject to the influence of horizontal/vertical arrangement[J]. International Journal of Heat and Mass Transfer, 2009, 52(11-12): 2565-2570.

[16] LIN C N. Analysis of three-dimensional heat and flow induced by piezoelectric fan[J]. International Journal of Heat and Mass Transfer, 2012, 55(11-12): 3043-3053.

[17] LIN C N. Enhanced heat transfer performance of cylindrical surface by piezoelectric fan under forced convection conditions[J]. International Journal of Heat and Mass Transfer, 2013, 60: 296-308.

[18] JENG T M, LIU C H. Moving-orientation and position effects of the piezoelectric fan on thermal characteristics of the heat sink partially filled in a channel with axial flow[J]. International Journal of Heat and Mass Transfer, 2015, 85: 950-964.

[19] 李鑫郡， 张靖周， 谭晓茗. 单个压电风扇传热特性研究[J]. 航空学报， 2017, 38(7): 120982.

LI X J, ZHANG J Z, TAN X M. Characteristics of heat transfer with a single piezoelectric fan[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(7): 120982 (in Chinese).

[20] JEONG J, HUSSAIN F. On the definition of a vortex[J]. Journal of Fluid Mechanics, 1995, 285: 69-94.