非對稱凹腔橫向陣列壓電風扇強化冷卻特性研究

張冬冬，譚曉茗，張靖周，鹿世化，李鑫郡,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)

3 結 論

本文主要采用數值方法針對非對稱凹腔內陣列排布的多壓電風扇系統激勵的三維非定常流動與傳熱特性進行研究，重點討論了振動相位和凹腔相對曲率的影響，得出結論如下：

（1）當非對稱凹腔相對曲率較大時（Kr=6），由于兩側側壁對于腔內弦向流動限制作用差異顯著，因此無論同相與反相振動時其表面時均對流換熱系數均呈現出明顯的非對稱分布形態。而當相對曲率Kr降低到2 時，非對稱分布特征幾乎完全消失。

（2）同相振動時換熱最強的區域出現在鄰近風扇的間隙，而反相振動時包絡區內換熱最強。

（3）同相振動多風扇系統激勵的脫落渦結構尺度明顯大于反相振動時的情況，使得同相振動平均換熱能力優于反相振動；且由于同相振動時激勵的間隙渦結構，也進一步增強了鄰近風扇間隙區域的換熱能力。