傅佳宏，張 宇，肖寶蘭，張旭方，左 強(qiáng)

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

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

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

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

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

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

(1)

(2)

其中,ν為流體運(yùn)動(dòng)粘度。粗糙度影響因子ΔB的定義式為

(3)

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

(4)

(5)

其中：u是流體的時(shí)均速度；Δy是到壁面的距離；μ是流體的動(dòng)力粘度。

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

(6)

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

(7)

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

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

2.1 MRF模型

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

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

(8)

式中，Vr為移動(dòng)坐標(biāo)系參考速度;V為絕對(duì)坐標(biāo)系速度;Vt為坐標(biāo)系移動(dòng)速度;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)格無(wú)關(guān)性等對(duì)于風(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)格時(shí)，取邊界層網(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)扇的流動(dòng)能夠充分發(fā)展，取其半徑與長(zhǎng)度皆為入口段的2倍，如圖2所示。

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

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

3 實(shí)驗(yàn)驗(yàn)證

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

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

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

(10)

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

(11)

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

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

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

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

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

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

4.2 湍流模型對(duì)比分析

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

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

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

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

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

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

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

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

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

圖7 風(fēng)扇性能試驗(yàn)值與理論計(jì)算值對(duì)比

5 結(jié) 論

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

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

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