翅片陣列對翼面的沖擊換熱影響數值模擬研究

張天倫， 王克辰， 張 栩， 周文武

(上海交通大學 機械與動力工程學院， 上海 200240)

當飛機穿過含有過冷水滴的低溫云層時，機翼表面尤其是前緣會形成積冰，破壞翼型的氣動外形，引起飛行性能下降甚至造成機毀人亡.Petty 等[1]總結得出：1975—1988年有超過800起由于飛機機翼結冰造成的安全事故.因此，高效的機翼防除冰技術對飛行安全至關重要.在眾多的防除冰方法中，熱氣膜防除冰是一種常見有效的防冰手段，從發動機引入的高溫空氣先經調壓閥進入笛形管(熱氣供應腔)，而后通過噴射孔沖擊到機翼內表面從而達到防除冰目的.整個過程中，氣膜防除冰效果與機翼沖擊換熱的性能緊密相關，因此提升曲面沖擊換熱性能一直是相關領域的研究重點.

國內外學者針對平板沖擊換熱開展了大量的研究.前期研究表明：射流孔的幾何形狀、射流孔距靶面的距離、靶面的幾何結構對射流的流動和換熱都有顯著的影響.Jambunathan等[2]總結了不同雷諾數、不同噴嘴形狀的平板局部和區域平均傳熱的實驗結果和經驗關聯式.Lytle等[3]利用激光-多普勒測速和壁面壓力測量研究了平板射流沖擊的流動結構.Narayanan等[4]通過實驗研究了3種流態沖擊平板時的流場、表面壓力和換熱率，結果表明換熱的二次峰位置與近壁的流動脈動速度變化相對應.Loureiro等[5]考慮了11個不同的實驗數據集，總結了經驗表達式，其中包括一個分段努塞爾數表達式.Siddique等[6]通過數值模擬確定了不同熱流輸入邊界條件下的局部努塞爾數.除了平板流動和換熱機理的研究，部分學者也探索了在平板表面添加粗糙單元(如三角形的翅片[7-9]、帶矩形槽道[8-9]、針肋[10-12]、微型肋加糙[13]、球形突起[14])對射流沖擊換熱的影響.

除了平板模型，半圓形曲面沖擊也是重要的研究對象.在改變圓形曲面的幾何結構方面，Singh[15]利用ANSYS Fluent 17.2對射流沖擊復合凹坑和凸凹表面進行了數值模擬，研究了凹面上凹陷和凸起的交替位置對沖擊傳熱的強化作用.不僅如此，Singh等[16]還建立了一個基于人工神經網絡的替代模型，以確定半圓形凸出物的最佳位置，最大限度地提高凹形表面的傳熱.Gau等[17]通過實驗研究了表面曲率對凹、凸表面沖擊冷卻流動和換熱過程的影響.Lyu等[18]測試了4個凹面(包括一個半圓柱凹面和3個具有不同寬深比的拋物線凹面)的射流沖擊換熱實驗結果，探究了曲率對凹表面單排人字形噴嘴射流沖擊換熱的影響.此外，劉錫晨等[19]進行了圓形和冠齒形射流口的射流實驗研究，實驗結果表明冠齒形噴嘴可以增大半圓形凹面的努塞爾數.

針對曲面沖擊在翼型防除冰的應用研究，Hannat等[20]通過ANSYS CFX數值模擬了三維防冰腔內部的流動與換熱，并與測試案例進行了結果對比.張靖周等[21]采用數值模擬的方法研究了射流孔直徑、射流孔間距等參數對凹腔表面溫度分布的影響.郭之強等[22]研究了光滑表面與具有長方體突起結構對內凹面的流動和換熱特性的影響，模擬結果表明凸起結構可以更加均勻地發散壁面射流.歸曉燁[23]計算曲面槽道內部強制對流，運用 Gnielinski 公式對機翼防冰系統槽道進行了分析，得到了槽道截面長度和寬度變化對換熱系數的影響.關于沖擊彎曲表面的數值模擬湍流模型的使用，Seyedein等[24-25]通過研究認為二方程的湍流模型只能模擬出部分的射流流動特性，改進k-ε和k-ωSST模型可以較好地擬合實際的射流沖擊曲面的換熱和流動特性.

雖然，針對射流沖擊機理以及結防冰應用開展了較多研究，然而之前研究大多以平板或半圓形曲面為主，這與真實的翼型曲面相差較大.因此，本文基于真實翼型結構，在幾乎不改動翼型設計的基礎上，數值探討如何在內壁面添加微小翅片陣列實現曲面沖擊換熱性能提高.以平板和NACA0012對稱翼型為研究模型，首先探索在平板上添加導流翅片陣列的流動換熱性能，而后推廣到曲面翼型.相關研究成果能為進一步提高機翼防除冰性能提供方法借鑒和設計參考.

1 射流沖擊流動結構

平板沖擊射流是一種被廣泛研究的流動現象，沖擊到平板后從駐點向外擴散的流動通常可以分為3個區域(見圖1)：第1個區域為自由射流區，在自由射流區段，來自周圍的流體被卷入射流，從而射流速度u0不斷降低，其垂直速度方向的射流面積逐漸增大；第2個區域為沖擊區，在該區域，垂直平板方向速度快速下降并偏轉，靜壓則快速上升，沖擊區在平板上的尺寸大約為兩個噴嘴直徑大小；第3個區域為壁面射流區，射流完全轉向，速度平行壁面，隨著流體的發展，逐漸出現速度邊界層，壁面射流表現出比平行流動更高的傳熱水平，這是由于壁面射流和周圍空氣之間的剪切所產生的湍流，在換熱表面被輸送到邊界層引起的.

圖1 射流流動分區Fig.1 Division of jet flow

在沖擊區和壁面射流區的過渡段處增加翅片陣列，起到翅片前擾流、翅片段整流的作用，一方面增大駐點附近的湍流水平，提高了駐點區域的換熱；另一方面增大遠離駐點區域的換熱面積和流動速度，提高了遠離駐點區域的換熱.最終實現整體區域換熱水平的提升.

2 平板射流沖擊

2.1 平板模型

圖2為平板計算域模型及網格劃分情況.平板沖擊模型和帶有的翅片陣列結構如圖2(a)和2(b)所示，考慮機翼防冰的笛形腔內的噴嘴長度與噴嘴直徑相當，因此噴嘴長度和噴嘴直徑d均為4 mm，噴嘴出口距靶面的高度為5d，平板厚度為0.3d.為保證出口流動對沖擊區域無影響，整體的計算域尺寸需足夠大，設置為60 mm×60 mm.靶面上的幾何結構除了基準的無翅片以外，還包括8片、12片直翅片和12片彎翅片3種，其中直翅片的長度、高度和厚度分別為d、0.5d和 0.25d，離駐點的距離為1.25d.彎翅形狀為1/4圓弧，圓弧半徑為噴嘴直徑大小.

圖2 平板計算域模型與網格Fig.2 Flat computing domain model and grid

2.2 數值計算

2.2.1計算域與邊界條件設置 平板沖擊的計算域同模型，采用ANSYS ICEM進行非結構化網格劃分，包括流體域和固體域，流體域與固體域的交界面節點一一對應.從圖2(c)和2(d)中可以看出，在射流區域、翅片附近進行了網格加密，在靶面的流體一側增加了壁面棱柱層，棱柱層的第1層厚度為0.005 mm，對應的無量綱壁面距離y+<1，網格質量在0.3以上.生成網格后在Fluent中進行邊界條件的設置和湍流模型的選擇，為了保證入口雷諾數與真實防冰系統射流雷諾數(Re≈10 000)匹配，入口速度按照不同的雷諾數(Re=5 500，8 600，10 000，15 000)分別設置為u0=23.70，36.518，42.61，63.91 m/s，溫度為335 K.流體域的四周為壓力出口，平板固體區域的下底面為恒熱流密度面，熱流密度為 -2 000 W/m2，四周面為絕熱壁面，設置參考Dobbertean[9]的模擬過程.空氣的物性參數求解使用的是理想氣體狀態方程.流體域計算所使用的湍流模型為Realizablek-ε，Enchanted Wall Function模型，該模型作為改進的k-ε算法，能夠模擬射流撞擊，特別是圓形射流，并且為求解靶面近壁處的流動換熱，采用增強壁面函數型的近壁面處理.固體域的溫度分布求解為邊界條件已知的控制體導熱微分方程的求解，材料設置為鋁.采用二階迎風格式的Coupled算法，對整個計算域進行了穩態的數值模擬.

2.2.2CFD驗證 共設計3套不同數目的網格進行網格獨立性的驗證，網格數分別為120萬、225萬和322萬.網格獨立性驗證的結果表明，當網格數大于225萬時，壁面上的溫度隨著網格加密的變化在0.3 K以內，綜合考慮計算精度和模擬效率，最終選擇網格數目225萬的網格.

為了確保網格和湍流模型的準確性，本文將數值模擬得到的駐點位置處的的局部努塞爾數大小與Garimella等[26]通過實驗歸納的實驗關聯式計算值進行了比較.局部努塞爾數的定義式如下：

(1)

式中：h0為局部的對流換熱系數；λ為空氣的導熱系數；qint為交界面處的熱流密度；Tint和Tj分別為交界面固體的溫度和射流的入口溫度.

Garimella等[26]歸納的實驗關聯式如下：

Nu0=0.462Re0.585Pr0.4(H/d)0.024

(2)

式中：Pr為空氣的普朗特數；H為垂直方向上離靶面的高度.式(2)的適用范圍為Re=4 000～23 000，H/d=1～5.

本文數值模擬得到4種雷諾數下平板駐點位置處的Nu0與實驗關聯式所得的結果對比如表1所示.表中：Er為本文模擬所得Nu0與實驗關聯式所得Nu0之間的相對誤差.由表可見，Er在10%以內，其中低雷諾數(Re=5 500)下數值模擬所得Nu0大于實驗關聯式計算的Nu0，高雷諾數下模擬所得Nu0小于實驗關聯式的Nu0.這是由于隨著雷諾數的增大， 流動再循環和混合更強，下游的熱空氣回流到駐點附近，導致Nu0的數值隨著表面熱流密度的變化有輕微的改變[27]，而本文在流固交界面處的熱流密度與文獻[26]中實驗的熱流密度不同.

表1 不同雷諾數下駐點位置的努塞爾數Tab.1 Nusselt numbers of stagnation point at different Reynolds numbers

2.2.3結果分析 為了更加綜合地評估增加翅片之后的換熱效能，流固交界面處的局部努塞爾數是不適用的，參考Dobbertean[9]對有厚度平板的努塞爾數的定義，定義綜合的努塞爾數：

(3)

式中：h為綜合定義的對流換熱系數；qcon為固體外壁面的熱流密度，該值在平板各處為一個恒定值；To為固體外壁面的溫度.該式從研究對象來看將整個平板的換熱進行了考慮，包括外表面的熱量流失、交界面的對流換熱、交界面到外表面的熱傳導，因此是一個綜合的換熱效能.

針對新定義的Nu，首先對比了4種結構在Re=5 500 時的結果，如圖3所示.圖中：x、y分別為以射流孔圓心為原點，X和Y方向下的距離.由圖可見，加8片、12片直翅片相比于不加翅片均可以提升換熱能力，并且8片的換熱效果最好.這是因為在12片直翅片的結構中，射流在向四周擴散時發生了堵塞，12片相較于8片周向面積減小了17.1%.加12片彎翅片相對于無翅片的換熱能力反而減弱，這是由于彎翅片不僅存在直翅片的堵塞問題，而且翅片是彎曲的，氣流在每個翅片的吸力面后緣均會出現明顯的分離，使得氣流更難向整個平板擴展出去，翅片換熱面積增大所增加的換熱小于因為射流無法擴散而導致的熱量堆積，所以換熱能力反而減弱了.

針對圖3分析得到的結果，后續對不同雷諾數換熱性能的分析，僅使用無翅片和增強換熱效果最優的8片直翅片的結果進行比較，如圖4所示.從圖4(a)、4(b)和4(c)、4(d)中最小和最大兩個雷諾數工況可以看出，安裝直翅片相較于無翅片的平板，Nu都是有提高的，這也意味著恒熱流外壁面的溫度有提升.

圖3 4種結構的努塞爾數分布云圖Fig.3 Nusselt number distribution of 4 kinds of structures

圖4 不同雷諾數下的努塞爾數分布云圖Fig.4 Nusselt number distribution at different Reynolds numbers

從(0，0)駐點位置處沿與X軸正方向22.5°方向取一條線上的Nu，得到一系列點線圖，即不同雷諾數下的Nu分布圖，如圖5所示.圖中：r為沿22.5°線上的點到靶面中心的距離；η為相對于無翅片時，有翅片Nu的提升比例.由圖可以直觀地看出換熱效能的提升效果，隨著雷諾數的增大，Nu隨之增大.相同雷諾數下，在徑向距離5d以內，有翅片的Nu均高于無翅片平板的Nu.在低雷諾數(Re=5 500)下，η最大，平均為5.5%左右，其余3個雷諾數下的η在1.3%左右.這是由于在低雷諾數下換熱效果由翅片面積帶來的提升與翅片增強局部對流換熱系數帶來的提升都比較明顯，而隨著雷諾數的增大，后者帶來的提升逐漸占據主導地位，所以整體提升效果沒有低雷諾數時那么明顯.此外，由于翅片的存在導致向外擴散的周向面積減小，從而向外擴散的射流流速增快，換熱增強，徑向距離在翅片之外的區域換熱也有提升.

圖5 綜合努塞爾數對比Fig.5 Comprehensive Nusselt number comparison

3 翼型曲面射流沖擊

3.1 翼型計算模型與網格

熱氣膜防除冰沖擊的壁面是機翼前緣的內表面，在平板上確認增加翅片對于換熱有提升，將該結構應用到翼型曲面上進行驗證.在翼型曲面模型中，整體計算域的寬度、射流孔距駐點的高度以及翼型曲面板的厚度保持不變，在曲面模型中研究了翅片長度對換熱的影響，在原先1倍噴嘴直徑長度的翅片基礎上，增加了2倍噴嘴直徑長度的翅片模型.整體網格如圖6所示，圖中：S為沿著外壁面上與葉片型線相切的方向.同樣在近壁面和翅片附近進行了加密，邊界條件和湍流模型的設置與平板的模型相同.

圖6 曲面計算域模型與網格Fig.6 Curved surface computing domain model and grid

3.2 翼型曲面沖擊換熱結果

圖7為無翅片、8片1倍噴嘴直徑長度翅片和2倍噴嘴直徑長度翅片的曲面板外壁面在Re=5 500 工況下的綜合努塞爾數分布，這里將翼型曲面按照圖6所示的型線S方向進行展平.圖中：s為翼型外表面上的點沿型線S方向距翼型前緣的距離；z為以射流孔圓心為原點，Z方向的距離.增加翅片后紅色的高溫區域明顯增多，并且2倍的翅片的效果要優于1倍的.因此，后續不同雷諾數下的流動換熱的分析在無翅片和2倍噴嘴直徑長度翅片的模型之間展開.

圖7 3種模型的綜合努塞爾數Fig.7 Comprehensive Nusselt number distribution of three models

翼型曲面板模型同樣研究了4種雷諾數下的換熱性能.圖8對比了不同雷諾數下的無翅片和有翅片的Nu分布，中間兩個雷諾數類似，因此選擇最大、最小兩個雷諾數的結果進行展示.在4種雷諾數的工況下，加翅片陣列模型的Nu在各個位置均要高于無翅片模型的Nu.在高雷諾數Re=15 000的時候，換熱效果的提升最為明顯，其現象與平板的結果相反.

圖8 不同雷諾數下的努塞爾數分布云圖(曲面)Fig.8 Nusselt number distribution at different Reynolds numbers (curved surface)

圖9顯示了從(s,z)=(0, 0)駐點位置處延展的兩條曲線上的綜合努塞爾數及其有2倍噴嘴直徑長度翅片時相對無翅片的提升比例η.從圖中可以看出，隨著雷諾數的增大，綜合努塞爾數相應提升，并且帶翅片的換熱始終要優于不帶翅片的.整體的提升比例在4%～10%左右，比平板的提升更為明顯.

結合圖8和圖9觀察到的現象，在翼型曲面上，隨著雷諾數的增大，綜合努塞爾數變化規律與平板不同，平板在小雷諾數下提升效果更明顯，而翼型曲面在高雷諾數下效果更好；此外，整體提升效果也不同，平板提升比例明顯小于翼型曲面.造成這種現象的原因有：① 曲面的曲率改變了射流的流場結構，從而對對流換熱的效果產生了影響；② 在翼型模型中采用的是2倍噴嘴直徑長度的翅片，而在平板上為1倍噴嘴直徑長度的翅片，翅片的長度從增大換熱面積上有不同的提升效果.

圖9 曲面綜合努塞爾數對比Fig.9 Curved surface comprehensive Nusselt number comparison

3.3 翼型曲面沖擊流場

翅片陣列對流場的影響對于平板和曲面的模型是類似的，均可以提高沖擊壁面附近的湍流度.湍流水平的提升可以導致傳熱率的增大，為了進一步說明換熱得到改善的內在機理，對Re=8 600 的工況進行了瞬態的數值仿真計算.圖10給出了距離曲面板內壁面0.5 mm處的近壁曲面上的湍動能強度k以及X、Y、Z方向上的速度u、v、w的均方根誤差(RMSE)，RMSE越大，代表速度的波動越大，即在數值模擬計算中速度的脈動量越大.翅片的存在使得k增大了1倍多；并且主要改變了Y、Z兩個方向上的速度脈動量，這兩個方向也正是平行于曲面的兩個方向.說明翅片的結構在駐點附近對沖擊射流的流場起到了擾流和增強湍流度的作用.

圖11為翅片方向上的時均流線圖.圖中：|V| 為流場合速度的大小.由放大的流線圖可見，射流在到達翅片時與翅片發生了撞擊，速度方向發生了改變，翅片前緣同時也是四周擴散流動的駐點位置，存在高水平的換熱.不過整體流場的結構沒有太大的變化，結合圖10的信息，可知翅片影響的是每一個時刻的流場，特別是平行于翼型曲面方向上的速度脈動量，而翅片對整體時均流場的影響比較小.

圖10 貼近內壁面的速度均方根誤差及湍動能(Re=8 600)Fig.10 RMSE of velocity and turbulent kinetic energy close to inner wall (Re=8 600)

圖11 截面流線圖Fig.11 Cross-section streamline diagram

4 結論

為了提高防除冰目標曲面的射流沖擊換熱，研究了翅片結構對平板和翼型曲面板上射流沖擊換熱的影響，對比不同雷諾數下的Nu，研究結果表明：

(1) 直翅片陣列的結構在平板可以提高射流沖擊的換熱強度，8片的提升效果最好，12片次之，8片1倍噴嘴直徑長度直翅片相較于無翅片可以提升1%～6%.

(2) 翅片的結構在翼型板上同樣可以提高射流沖擊的換熱強度，8片2倍噴嘴直徑長度長直翅片在曲面板上的提升為4%～10%，換熱水平的提升可以提高防除冰效果.雖然增加翅片在平板和翼型曲面上均可以提升換熱強度，不過，由于曲面曲率的影響，射流在曲面上的流動傳熱和平板是有區別的，后續的研究將針對真實的翼型進行探索.

(3) 翅片結構不僅可以增大駐點以及周邊區域的湍動能強度從而提升換熱性能(湍動能的提升達到了1倍左右)，也可以從增加換熱面積的角度進一步提高換熱.

(4) 低雷諾數下，翅片通過改變流動結構增加的換熱量與增大換熱面積增加的換熱量都是比較明顯的；高雷諾數下，則以改變流動結構的提升為主.