基于表面粗糙度修正的冷卻風扇數值模擬方法研究

傅佳宏，張 宇，肖寶蘭，張旭方，左 強

(浙江大學城市學院 工程學院， 杭州 310015)

軸流式冷卻風扇在工業領域具有廣泛的應用，對其性能進行精準的數值模擬研究，在冷卻系統匹配、節能減排、氣動噪聲抑制等方面具有重要的工程應用意義[1-2]。風扇旋轉時，通過具有不同壓力曲線的風扇葉片與空氣、護風罩之間復雜的雙向流固耦合[3-4]及周期性作用，形成旋轉流場與風扇前后的靜壓升。

目前，工程領域的主流風扇性能預測模型是多重參考坐標系(MRF)計算模型，該模型將風扇流動區域定義在旋轉坐標系下，其它區域定義在靜止參考系下，通過坐標系的不斷重建來模擬風扇的旋轉。由于模型中包含詳細的風扇葉片幾何信息，可以直接模擬出壓力階躍以及旋轉流場等風扇特性，在預測精度及計算資源的消耗方面具有相當的優勢[5-6]。由于采用了定常的近似求解，如何對MRF方法進行修正，進一步提高其模擬精度仍然面臨著挑戰。Wang等[7]研究了風扇旋轉區的選擇對其性能的影響，認為風扇旋轉區應該為護風罩所包含區域，而不僅僅是風扇直徑與厚度形成的圓柱區。Shankar等[8]采用MRF模型建立了風扇的數字風洞實驗室用以進行實際風扇的修正。Gullberg等[9-11]認為對于風扇旋轉流體區域的選擇要包括與風扇有劇烈相互作用的區域和風扇附近非對稱旋轉部件，并提出將風扇轉速提高14%所得的數值模擬結果可很好地與實驗值進行匹配。Sengupta等[12]指出，通過MRF模型旋轉流體區域、紊流模型、風扇網格質量進行一定的適應性修正后，完全可以達到工程實際的需求。耿麗珍等[13]采用MRF進行了風扇噪聲分析方法的數值仿真研究。肖紅林等[14]研究指出可以采用減小風扇和導風罩之間的間隙的方法來提高風扇模擬的精度。倪計民等[15]采用MRF方法研究了風扇與導風罩之間相對位置關系對于冷卻風扇性能的影響，得到不同風扇具有不同最佳安裝參數的結論。石海民等[16]進行了多風扇之間的MRF建模方法研究。

在上述的MRF計算模型修正中，都將風扇表面作為水力光滑面處理，而實際風扇表面粗糙度對于風扇表面渦的流動阻力以及壁面函數都有一定的影響。文中結合壁面函數，通過第一層網格離壁面的無因次距離y+，在湍流計算中加入粗糙度函數對壁面律做出修正，使得風扇表面更貼近實際情況，進一步提高風扇的數值模擬精度。

1 風扇表面粗糙度修正

通常粗糙度由表征粗糙顆粒類型的粗糙度常數Cs與粗糙度高度Ks決定，共同組成了粗糙度影響因子ΔB，對于風扇葉片表面，影響公式為[17]

(1)

(2)

其中,ν為流體運動粘度。粗糙度影響因子ΔB的定義式為

(3)

引入無量綱參數u+與y+，分別為第一層網格離壁面的無因次速度與距離：

(4)

(5)

其中：u是流體的時均速度；Δy是到壁面的距離；μ是流體的動力粘度。

對于黏性底層有u+=y+，從而得到風扇葉片表面摩擦速度uτ為

(6)

聯立式(5)與式(6)可得到粗糙度高度Ks與風扇葉片第一層網格高度Δy的迭代方程組：

(7)

通過對風扇表面粗糙度的修正，可以讓風扇的數值模擬更接近實際情況，結合網格、湍流模型等因素進一步分析修正后風扇數值模擬結果，并進行實驗驗證。

2 風扇數值模型的構建

2.1 MRF模型

通過風扇三維實體模型建立風扇數值仿真模型，采用MRF模型進行風扇性能仿真，其主要方法是在風扇流體區域建立多重參考坐標系，坐標系能夠隨著流體的旋轉而不斷重建。固定坐標系與移動坐標系中的速度矢量按照如下關系轉換[18]：

Vr=v-(Vt+w×r),

(8)

式中，Vr為移動坐標系參考速度;V為絕對坐標系速度;Vt為坐標系移動速度;W為旋轉角速度;r為風扇半徑。

根據以上變換，在風扇旋轉流體區域中，控制方程組為

(9)

2.2 邊界條件的設置

以某半徑為762 mm的車用冷卻風扇為例，采用風扇周期性邊界條件，分析紊流模型、葉片表面粗糙度、近壁面網格處理、網格無關性等對于風扇數值模擬結果的影響，從而提高風扇數值仿真精度。模擬風扇共有10片葉片，風扇直徑為762 mm，轉鼓直徑為245 mm，給定風扇轉速為2 100 r/min，如圖1所示。

圖1 風扇幾何模型示意圖

設置初始迭代條件為y+，得到Δy=0.000 892 m，Ks=0.000 046 4 m，因此，在設置邊界層網格時，取邊界層網格厚度為0.9 mm，葉片表面摩擦度高度為0.046 4 mm。

風扇采用周期性旋轉邊界建立風扇數值模擬風洞，風扇為十葉片均布，風扇的旋轉偏移角為36°，設置風扇入口段半徑與風扇半徑一致，出口段為了使風扇的流動能夠充分發展，取其半徑與長度皆為入口段的2倍，如圖2所示。

圖2 風扇數值模型示意圖

風扇流場采用MRF模型進行穩態數值計算，由于沒有熱交換器參與熱交換，忽略能量方程，設置流動介質空氣的物性參數為常物性；入口與出口分別設為壓力入口與壓力出口，表壓設置為0 Pa；采用基于壓力修正的SIMPLEC算法進行流場計算，采用二階迎風格式進行網格的離散，當殘差小于等于10-4次方時認為計算收斂。數值模擬在20核，32 G內存，2 T硬盤高性能工作站上進行，每個case根據網格數量的多少大約在1.5～3.0 CPU時收斂。

3 實驗驗證

按照如圖3所示的工業通風機標準化風道上進行數值模型的試驗驗證。試驗時通過改變調速電機的轉速來控制風扇轉速，通過改變節流加載板的孔隙率來改變風道的進氣阻力，從而得到風扇在不同轉速n下的風扇流量qv與靜壓H曲線。具體可參照文獻[19]。

注：1.錐形集流器；2.節流加載板；3、5.整流器；4.試驗風筒；6.空氣流量計；7.風扇；8.扭矩轉速傳感器；9.調速電機

風扇全壓ptp為冷卻空氣通過風扇后總壓的升高量，包括冷卻空氣靜壓與動壓之和，可表示為

(10)

風扇靜壓H表示為

(11)

其中：pout，pin分別為風扇出口和進口處靜壓，通常pout與大氣環境相連，可認為pout=0；cin、cout分別為風扇進出口軸向速度，由于風扇進出口截面積大致相等，可認為cin=cout；ρ為空氣密度；n為風扇轉速，r/min；qv為風扇流量，m3/s。

4 數值仿真結果分析

4.1 網格無關性分析

不同網格密度下得到的風扇網格無關性分析，如圖4所示。圖4(a)為當風扇轉速為2 100 r/min，流量為0.5 m3/s時，風扇的靜壓值隨著網格密度變化的規律。圖4(b)為風扇轉速為2 100 r/min時，不同網格密度下得到的風扇性能曲線。

圖4 風扇數值模型網格無關性分析

可以看出，隨著風格的增加，靜壓值不斷增加，當風扇旋轉流體區域網格增加到57萬時，隨著網格的增加，風扇靜壓值幾乎不變。在各個網格密度下，風扇性能曲線趨勢基本一致，主要在曲線的兩端距試驗值有一定的偏差，當網格數量小于30萬時，數值仿真結果偏差較大，隨著網格數量的增加，偏差逐漸減小；然而當網格密度達到87萬時，數值仿真結果并沒有顯著的改善，反而在流量超過5 m3/s時，靜壓值上出現了明顯的流量波動。因此，最后選取風扇旋轉流體區域網格為57萬。

4.2 湍流模型對比分析

在計算風扇旋轉流場等有強旋流和有彎曲壁面的流動時，可能出現時均應變率特別大的情形，采用標準k-ε模型時，針對此情形有可能產生負的正應力，從而導致流動失真，因此，在對具有旋轉流動的流場進行數值模擬時，通常采用考慮旋轉流場的RNG k-ε模型或者能夠反映主流時均應變率的Realizable k-ε模型。

相同風扇旋轉流體區域下，分別采用RNG k-ε模型和Realizable k-ε模型進行風扇數值模擬，不同的紊流模型對于風扇數值模擬精度的影響如圖5所示。為了減小網格因素對于仿真精度的影響，采用加密網格進行數值計算，加密后，風扇網格數量為59萬，計算工況為風扇轉速2 100 r/min。

圖5 不同紊流模型下風扇數值模擬結果對比

從圖5中可以看出，相比Realizable k-ε模型，RNG k-ε模型具有更高的仿真精度，在風扇工作范圍內(2～5 m3/s)，仿真值與試驗值吻合較好，誤差小于5%，但在低流量區域(0～2 m3/s)誤差較大。因為在低流量、高靜壓工況下，湍流模型自身的局限性及風扇葉片壓力面和吸力面之間的高壓力梯度，RNG k-ε紊流模型在極端情況下不能很好地模擬出風扇的流動特性，但在風扇常用工況點，該模型具有足夠的精度，具有一定的工程實用價值。

Realizable k-ε模型在整體趨勢上能夠模擬出風扇的氣動特性，但是存在較為明顯的過度模擬現象，主要表現為在中低流量工況(0～4 m3/s)下，仿真靜壓值要明顯小于試驗值，但是在大流量工況(>4 m3/s)下，仿真靜壓值又要明顯大于試驗值，因為數值仿真將光滑風扇葉片進行了人為的離散，無法完全模擬出光滑表面，在局部可能出現高曲率的情況，導致應變率出現極值點，影響了整個風扇的數值仿真精度。

4.3 粗糙度修正結果分析

采用近壁面粗糙度修正與僅采用簡單壁面函數的風扇數值仿真結果與試驗值的對比如圖6所示。可以明顯看出，在增加了近壁面處理后，數值仿真結果更接近試驗值，因為其流動更符合實際情況，雖然在低流量區域靜壓值與試驗值還是存在100 Pa左右的偏差，但經過近壁面處理后風扇性能曲線更為平順，改善了未經過近壁面處理時產生的一些畸點，且在風扇常用工況范圍內(2～5 m3/s)，試驗值與仿真值的偏差在3%以內，完全可以滿足工程應用的需要。

圖6 壁面處理前后計算結果對比

圖7所示為比較得到的相對較好的57萬網格數量，RNG k-ε紊流模型，表面粗糙度修正之后，風扇在各個轉速下的數值模擬結果與試驗值的對比，從中可以看出，除了高壓低速區，試驗值與模擬值都能較好的吻合，可見文中提出的風扇數值模擬方法可以提高全工況范圍內風扇性能的數值模擬精度。

圖7 風扇性能試驗值與理論計算值對比

5 結 論

1)分析比較了Realizable k-ε紊流模型與RNG k-ε紊流模型對于風扇數值仿真精度的影響， 發現后者在常用工況點具有足夠的精度與計算效率。

2)常用風扇MRF計算模型基礎上增加了風扇葉片表面摩擦度修正，提高了風扇數值仿真精度。

3)通過全轉速下的風扇模擬與試驗驗證，證明了該方法的普適性，可以采用風扇相似定理得到不同轉速下的風扇性能。