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

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

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

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

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

1 風(fēng)扇表面粗糙度修正

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

(1)

(2)

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

(3)

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

(4)

(5)

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

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

(6)

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

(7)

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

2 風(fēng)扇數(shù)值模型的構(gòu)建

2.1 MRF模型

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

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

(8)

式中，Vr為移動坐標系參考速度;V為絕對坐標系速度;Vt為坐標系移動速度;W為旋轉(zhuǎn)角速度;r為風(fēng)扇半徑。

根據(jù)以上變換，在風(fēng)扇旋轉(zhuǎn)流體區(qū)域中，控制方程組為

(9)

2.2 邊界條件的設(shè)置

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

圖1 風(fēng)扇幾何模型示意圖

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

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

圖2 風(fēng)扇數(shù)值模型示意圖

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

3 實驗驗證

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

注：1.錐形集流器；2.節(jié)流加載板；3、5.整流器；4.試驗風(fēng)筒；6.空氣流量計；7.風(fēng)扇；8.扭矩轉(zhuǎn)速傳感器；9.調(diào)速電機

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

(10)

風(fēng)扇靜壓H表示為

(11)

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

4 數(shù)值仿真結(jié)果分析

4.1 網(wǎng)格無關(guān)性分析

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

圖4 風(fēng)扇數(shù)值模型網(wǎng)格無關(guān)性分析

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

4.2 湍流模型對比分析

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

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

圖5 不同紊流模型下風(fēng)扇數(shù)值模擬結(jié)果對比

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

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

4.3 粗糙度修正結(jié)果分析

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

圖6 壁面處理前后計算結(jié)果對比

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

圖7 風(fēng)扇性能試驗值與理論計算值對比

5 結(jié) 論

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

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

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