郭向東，王梓旭，李明，劉蓓

1.中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點實驗室，綿陽 621000 2. 中國空氣動力研究與發(fā)展中心 飛行器結(jié)冰與防除冰重點實驗室，綿陽 621000

飛機結(jié)冰廣泛存在于飛行實踐中，并嚴(yán)重威脅飛行安全[1-2]。研究飛機結(jié)冰的主要途徑包括數(shù)值模擬、結(jié)冰風(fēng)洞試驗和飛行試驗3種，其中結(jié)冰風(fēng)洞試驗顧名思義是利用結(jié)冰風(fēng)洞在地面模擬飛機結(jié)冰過程，驗證飛機防除冰系統(tǒng)，開展飛機結(jié)冰問題研究[3]。相對于數(shù)值模擬和飛行試驗，結(jié)冰風(fēng)洞試驗具有結(jié)冰條件易控、試驗成本相對較低、結(jié)果可靠等優(yōu)點，是目前研究飛機結(jié)冰的主要手段。為模擬高空低溫低壓的云霧環(huán)境，結(jié)冰風(fēng)洞主要包括制冷系統(tǒng)、高度模擬系統(tǒng)和噴霧系統(tǒng)三大部分，其中噴霧系統(tǒng)利用噴嘴在風(fēng)洞穩(wěn)定段內(nèi)生成大量懸浮液滴，隨后液滴由低溫氣流攜帶進入試驗段，形成試驗所需的過冷水滴云霧[4]。結(jié)冰試驗中，氣流靜溫低于冰點，為防止噴霧系統(tǒng)中液態(tài)水結(jié)冰，須采用加熱方式提高液態(tài)水溫度，這樣便導(dǎo)致穩(wěn)定段內(nèi)生成的液滴溫度高于氣流靜溫，通常情況下，液滴會與氣流發(fā)生充分的熱交換，進而在試驗段內(nèi)達到過冷狀態(tài)(液滴溫度與試驗段靜溫一致)，形成試驗所需的過冷水滴[5]。但是，由于水滴與氣流間的熱交換過程十分復(fù)雜，涉及對流傳熱、相變傳質(zhì)和三維收縮效應(yīng)等多物理過程的耦合作用，試驗參數(shù)的變化可能導(dǎo)致試驗段內(nèi)液滴偏離過冷狀態(tài)，進而改變積冰生長過程和冰形特征，影響結(jié)冰風(fēng)洞試驗結(jié)果，因此探索結(jié)冰風(fēng)洞中液滴傳熱過程，揭示試驗參數(shù)的影響規(guī)律，評估液滴過冷狀態(tài)，對結(jié)冰風(fēng)洞試驗準(zhǔn)確模擬結(jié)冰過程具有重要意義。

針對結(jié)冰風(fēng)洞液滴過冷問題，Willbanks和Schulzt[6]率先發(fā)展了基于拉格朗日法的運動液滴傳質(zhì)傳熱計算方法，并對發(fā)動機高空試驗臺內(nèi)的液滴傳熱過程開展了研究，隨后Miller等[7]采用這一方法，研究了NASA Glenn IRT結(jié)冰風(fēng)洞內(nèi)的液滴過冷問題，發(fā)現(xiàn)初始液滴溫度對液滴過冷狀態(tài)無顯著影響，此外Bellucci等[8]采用試驗方法對CIRA結(jié)冰風(fēng)洞高速試驗段和低速試驗段構(gòu)型的液滴過冷狀態(tài)進行了評估。但是這些研究主要關(guān)注工程層面上液滴過冷狀態(tài)的評估，缺乏對結(jié)冰風(fēng)洞中液滴傳熱傳質(zhì)過程的規(guī)律性認(rèn)識。

近年來，隨著國內(nèi)大型結(jié)冰風(fēng)洞——3 m×2 m結(jié)冰風(fēng)洞的建成，液滴過冷問題逐漸引起了一些國內(nèi)學(xué)者的關(guān)注[9-10]。其中郭向東等[10]發(fā)展了基于Euler法的氣液兩相耦合流動計算方法，研究了3 m×2 m結(jié)冰風(fēng)洞主試驗段構(gòu)型液滴傳熱過程，并評估了液滴過冷狀態(tài)，研究結(jié)果表明結(jié)冰風(fēng)洞中液滴傳熱過程可以分為準(zhǔn)一維傳熱和三維收縮傳熱兩個階段，其中三維收縮效應(yīng)對液滴過冷狀態(tài)起決定性作用，但是他們的研究中未考慮液滴相變傳質(zhì)對液滴傳熱過程的影響，導(dǎo)致對結(jié)冰風(fēng)洞液滴過冷狀態(tài)的評估結(jié)果較為保守。

因此，本文在郭向東等的計算方法基礎(chǔ)上，進一步引入Hill矩方法[11]，發(fā)展了基于Euler法的氣液兩相傳質(zhì)傳熱耦合流動計算方法，實現(xiàn)了液滴相變傳質(zhì)過程的模擬。進而利用該方法，針對3 m×2 m結(jié)冰風(fēng)洞主試驗段構(gòu)型，考察了液滴相變效應(yīng)對液滴傳熱過程的影響，揭示了相變效應(yīng)的影響規(guī)律，開展了參數(shù)影響研究，分析了相對濕度、試驗段氣流速度和液滴尺寸對液滴傳質(zhì)傳熱過程的影響，評估了液滴過冷狀態(tài)，為國內(nèi)飛機結(jié)冰地面模擬試驗數(shù)據(jù)應(yīng)用評估提供支撐。

1 計算方法

1.1 控制方程和物理模型

本文借鑒了文獻[10]中基于Euler法的氣液兩相耦合流動計算方法，進一步耦合Hill矩函數(shù)方程組[12-14]，實現(xiàn)了相變傳質(zhì)過程的模擬。為獲得簡化物理模型，對氣液兩相進行以下假設(shè)[15]：①氣相為理想氣體，遵循理想氣體法則；②液滴為球形，不發(fā)生變形和破碎，液滴間不發(fā)生碰撞；③液滴內(nèi)溫度均勻分布；④忽略重力效應(yīng)；⑤ 黏性和熱傳導(dǎo)效應(yīng)僅發(fā)生在氣液兩相間；⑥不考慮液滴凍結(jié)過程以及水蒸氣凝結(jié)成核過程。

基于以上假設(shè)，簡化后的控制方程組為

1) 氣相方程組

(1)

(2)

(3)

(4)

2) 液相方程組

(5)

(6)

(7)

3) 矩函數(shù)方程組

(8)

(9)

(10)

式中：下標(biāo)g和d分別代表氣相和液相;ρ、p、V和E分別為密度、壓力、速度和單位質(zhì)量總能量；氣相能量Eg和液相能量Ed的表達式分別為

(11)

式中：cv為氣相定容比熱；cl為液滴比熱；T為溫度。

σ和α分別為液滴有效密度和水蒸氣質(zhì)量分?jǐn)?shù)，其中σ對應(yīng)液態(tài)水含量(LWC)；n階矩函數(shù)定義為

(12)

其中：r為液滴半徑；f為粒子尺寸分布函數(shù)；Q0、Q1和Q2為零階、一階和二階矩函數(shù)，根據(jù)定義可見Q0表示單位體積內(nèi)液滴總個數(shù)。

S0為質(zhì)量源項，表示為

(13)

式中：ρl為液滴物理密度；dr/dt為液滴半徑變化率。

SD和SM分別為兩相動量交換源項和兩相間傳質(zhì)引起的動量傳遞源項，表示為

SD=Q0FD,SM=S0Vd

(14)

式中：FD為單個液滴受到的氣相作用力。

SEg和SEd分別為氣相和液相能量源項，表示為

(15)

式中：qheat與qrad分別為單個液滴與氣相間的對流傳熱率和輻射傳熱率；qgas,mass和qliquid,mass分別為單位質(zhì)量傳質(zhì)引起的氣相和液相能量變化率，表示為[16-17]

(16)

其中：cp為水蒸氣定壓比熱；L為汽化潛熱。

S1和S2分別為一階和二階矩函數(shù)方程源項，表示為

(17)

單個液滴受到的氣相作用力FD和單個液滴與氣相間的對流傳熱率qheat表示為[18]

(18)

qheat=Nukgπd(Tg-Td)

(19)

式中：d為液滴直徑；μg為氣相黏性系數(shù)；kg為氣相熱傳導(dǎo)系數(shù)；Red為液滴相對雷諾數(shù)，表示為

(20)

CD為液滴阻力系數(shù)，針對球型液滴，采用Schiller-Naumann圓球阻力模型，表示為

(21)

針對球型液滴，Nu采用Ranz-Marshall半經(jīng)驗關(guān)系式，表示為

(22)

式中：Pr為Prandtl數(shù)。

單個液滴與氣流間的輻射傳熱率qrad表示為[19]

(23)

式中：ε為熱輻射發(fā)射率，根據(jù)文獻[19]，ε取0.95;σSB為Stefan-Boltzmann常數(shù)。

(24)

式中：DAB為二元質(zhì)量擴散系數(shù)；Rv為水蒸氣氣體常數(shù)；pv,s為液滴表面飽和水氣壓；Hu為相對濕度，表示為

(25)

式中：pv,∞為氣相環(huán)境中水蒸氣分壓。

Sh為Sherwood數(shù)，針對球型液滴，Sh采用Ranz-Marshall半經(jīng)驗關(guān)系式，表示為

(26)

式中：Sc為Schmidt數(shù)。根據(jù)文獻[16，21]，Ranz-Marshall半經(jīng)驗關(guān)系式適用于Red<137、溫度低于0 ℃時球型液滴的傳質(zhì)傳熱過程，匹配本文的研究范圍(尤其對于高風(fēng)速工況下依然成立)，因此本文計算模型中Nu和Sh采用該關(guān)系式。

1.2 數(shù)值方法

本文采用與文獻[10]相同的數(shù)值方法，控制方程組采用有限體積方法離散，非定常項采用二階Euler格式離散，空間輸運項采用二階迎風(fēng)型的Roe格式離散，源項采用隱式格式[10]。

氣液兩相邊界條件采用不同形式處理：氣相入口條件采用壓力入口，出口條件采用壓力出口，壁面采用無黏滑移壁面條件(1.1節(jié)假設(shè)⑤)；液相入口條件采用速度入口，出口采用出流條件，壁面處認(rèn)為液滴可以直接穿過，不考慮壁面對液相的作用。

1.3 驗證算例

選取文獻[21]中的試驗結(jié)果驗證本文計算方法。文獻[21]采用懸掛液滴試驗方法研究了液滴凍結(jié)過程中的液滴溫度變化特性，文中將液滴凍結(jié)過程分為過冷(Supercooling Stage)、復(fù)輝(Recalescence Stage)、凍結(jié)(Freezing Stage)和冷卻(Cooling Stage)4個階段，其中過冷階段對應(yīng)本文研究的內(nèi)容，因此選取文獻中過冷階段的試驗結(jié)果驗證本文計算方法。驗證結(jié)果如圖1所示，圖中：Va和Ta分別表示氣流速率和氣流靜溫。從圖中可以看出：本文計算結(jié)果與試驗結(jié)果吻合得較好，不同工況下液滴溫度誤差均小于1%，因此本文計算方法可以合理模擬液滴過冷階段的傳熱過程。針對以上的驗證過程，應(yīng)該指出的是考慮到懸掛液滴法對液滴尺寸和試驗風(fēng)速的限制，該方法僅獲得了大液滴(毫米量級)低風(fēng)速(小于1 m/s)下的試驗結(jié)果，而缺乏小液滴(微米量級)高風(fēng)速的試驗結(jié)果。進一步，文獻[21]指出，盡管缺乏真實噴霧的試驗數(shù)據(jù)，但是試驗數(shù)據(jù)中液滴相對雷諾數(shù)范圍(36

2 結(jié)果與分析

本節(jié)針對3 m×2 m結(jié)冰風(fēng)洞主試驗段收縮構(gòu)型，采用本文發(fā)展的計算方法，考察了基準(zhǔn)工況下液滴相變效應(yīng)對液滴傳熱過程的影響，揭示了相變效應(yīng)的影響規(guī)律，開展了參數(shù)影響研究，分析了相對濕度、試驗段氣流速度和液滴尺寸對液滴傳質(zhì)傳熱過程的影響，評估了試驗段中心處液滴過冷狀態(tài)。

3 m×2 m結(jié)冰風(fēng)洞主試驗段構(gòu)型如圖2所示，該構(gòu)型包括穩(wěn)定段、收縮段和試驗段3部分，其中構(gòu)型入口截面為高8 m×寬11 m，位于穩(wěn)定段噴霧耙處，出口截面為高2 m×寬3 m，位于試驗段中心處，構(gòu)型全長16 m，收縮比為14.67。

根據(jù)文獻[10]的研究結(jié)果可知，試驗段氣流速度和液滴尺寸是影響液滴過冷狀態(tài)的主要參數(shù)，因此本文選取相對濕度、試驗段氣流速度和液滴尺寸3個參數(shù)開展參數(shù)研究，計算工況矩陣如表1所示，表中LWCi、di、Tdi、T0i、Udi、UTS和Hui分別代表初始液態(tài)水含量、初始液滴直徑、初始液滴溫度、初始?xì)饬骺倻?、初始液滴速率、試驗段氣流速度、初始相對濕度。Case 1為基準(zhǔn)工況，考察典型工況下結(jié)冰風(fēng)洞中相變對液滴傳熱過程的影響；Case 2選取了70%、80%、90%和100%這4個典型的相對濕度狀態(tài)，考察相對濕度的影響；Case 3選取了55 、77、130、164 m/s這4個典型試驗段風(fēng)速，考察氣流速度的影響，其中各風(fēng)速對應(yīng)的入口總壓分別為2 000、4 000、12 000、20 000 Pa，參考壓力為105Pa；Case 4選取了40、60、80、100 μm這4個典型液滴直徑，考察液滴尺寸的影響，值得指出的是由于假設(shè)液滴為球型，不考慮液滴變形和破碎等動力學(xué)特性，因此僅選取了小尺寸液滴開展研究。

表1 計算工況矩陣Table 1 Calculation condition matrix

2.1 基準(zhǔn)工況

基準(zhǔn)工況選取FAR Part25 附錄C中典型結(jié)冰條件[22]，為考察相變對液滴傳熱過程的影響，下文通過對比構(gòu)型內(nèi)液滴和顆粒的傳熱過程，進而揭示液滴相變效應(yīng)的影響規(guī)律，其中針對顆粒傳熱過程的模擬，本文采用文獻[10]中的計算方法來實現(xiàn)。

圖3給出了基準(zhǔn)工況下構(gòu)型中心線處氣流靜溫和液滴溫度變化曲線，從圖中可以看出，類似于文獻[10]中顆粒傳熱過程，整個液滴傳熱過程近似以構(gòu)型6 m處截面為界分成準(zhǔn)一維傳熱和三維收縮傳熱兩個階段。在準(zhǔn)一維傳熱階段內(nèi)(入口至6 m處截面，包括穩(wěn)定段和收縮段入口區(qū)域)，風(fēng)洞型面變化較緩，氣流以準(zhǔn)一維方式低速流動，氣流速度和靜溫變化較小，接近入口流速和總溫，液滴以準(zhǔn)一維穩(wěn)態(tài)方式傳熱；在三維收縮傳熱階段內(nèi)(6 m處截面至出口，包括收縮段后部區(qū)域和試驗段)，風(fēng)洞型面變化較大，三維收縮效應(yīng)較強，氣流場在此階段內(nèi)發(fā)生了顯著變化，氣流速度加速增大，氣流靜溫則等熵下降，出口氣流速度越大，出口氣流靜溫越低，三維收縮效應(yīng)(氣流可壓縮性)顯著影響液滴傳熱過程，液滴以三維動態(tài)方式傳熱。對比液滴與顆粒的溫度變化曲線可見：在準(zhǔn)一維傳熱階段，液滴溫度下降趨勢快于顆粒溫度，當(dāng)液滴溫度和顆粒溫度趨于穩(wěn)定時(約在0.6 m處)，液滴溫度(約-5.4 ℃)低于顆粒溫度(-5 ℃)，此時液滴溫度稱為液滴濕球溫度[16]；在三維收縮傳熱階段，液滴溫度下降趨勢慢于顆粒，最終在構(gòu)型出口處(試驗段中心)，液滴溫度(-7.41 ℃)高于顆粒溫度(-7.77 ℃)。

圖4給出了構(gòu)型中心線處液滴直徑和相對濕度變化曲線，從圖中可以看出：結(jié)冰風(fēng)洞中液滴經(jīng)歷了先蒸發(fā)后凝結(jié)兩個階段，在準(zhǔn)一維傳熱階段，氣流中水蒸氣未達到飽和狀態(tài)，液滴處于蒸發(fā)狀態(tài)，液滴直徑不斷減小，同時蒸發(fā)產(chǎn)生的水蒸氣則使氣流的相對濕度不斷增大；在三維收縮階段，氣流靜溫迅速下降，氣流飽和水汽壓降低，而水蒸氣分壓變化較小(適用于低速流動)，導(dǎo)致氣流的相對濕度迅速增大(對應(yīng)式(25))，氣流中水蒸氣從未飽和狀態(tài)轉(zhuǎn)換為過飽和狀態(tài)，液滴則從蒸發(fā)狀態(tài)轉(zhuǎn)換為凝結(jié)狀態(tài)。

進一步理論分析相變效應(yīng)對液滴傳熱的影響。首先液相能量方程式(7)可以化簡為

(27)

式中:由于液滴過冷階段輻射傳熱遠(yuǎn)小于對流傳熱和相變傳熱，因此忽略輻射傳熱項，同時聯(lián)立式(19)和式(24)，則式(27)可化簡為

(28)

ξheat+ξmass

(29)

式中：ξheat為對流傳熱引起的液滴溫度空間變化率；ξmass為液滴相變引起的液滴溫度空間變化率。需要指出的是，式(29)是在一維定常流動條件下得到的，但是對于三維收縮階段，構(gòu)型中心線處的氣流沿流向加速流動，仍可以近似認(rèn)為氣流運動為準(zhǔn)一維，因此式(29)適用于構(gòu)型中心線處整個液滴傳熱過程的分析。

圖5給出了構(gòu)型中心線處ξmass和ξheat的變化曲線，從圖中可以看出，在準(zhǔn)一維傳熱階段，當(dāng)溫度較高的液滴進入未飽和的低溫氣流(Hui<100%)中時，液滴與氣流間發(fā)生對流傳熱，液滴向氣流傳遞熱能(ξheat<0)，同時蒸發(fā)效應(yīng)從液滴吸收所需的汽化潛熱(ξmass<0)，而顆粒與氣流間僅發(fā)生對流傳熱過程，因此液滴溫度下降速度快于顆粒的速度(0～0.6 m)；當(dāng)液滴溫度下降至氣流靜溫時[23]，由于蒸發(fā)效應(yīng)仍持續(xù)從液滴吸收熱量，因此液滴溫度會持續(xù)下降，但對于顆粒而言，其溫度此時將會與氣流靜溫保持一致，顆粒與氣流間的對流傳熱過程停止，顆粒達到穩(wěn)定狀態(tài)；進一步，當(dāng)液滴溫度下降至液滴濕球溫度時，蒸發(fā)效應(yīng)所需的汽化潛熱(ξmass<0)與氣流向液滴傳遞的熱能(ξheat>0)達到平衡(|ξmass|=|ξheat|)，液滴溫度趨于穩(wěn)定，此時ξmass與ξheat沿流向(x軸)對稱分布(0.6～6 m)，液滴進入穩(wěn)定蒸發(fā)階段。在三維收縮階段，三維收縮效應(yīng)導(dǎo)致氣流靜溫迅速下降，液滴穩(wěn)定蒸發(fā)過程被破壞，ξmass與ξheat不再對稱分布，液滴溫度下降(ξmass+ξheat<0)；此階段內(nèi)對流傳熱的能量傳遞方向發(fā)生了轉(zhuǎn)換，即ξheat>0轉(zhuǎn)換為ξheat<0，促進了液滴溫度下降，同時氣流中水蒸氣從未飽和狀態(tài)轉(zhuǎn)換為過飽和狀態(tài)，液滴則從蒸發(fā)狀態(tài)轉(zhuǎn)換為凝結(jié)狀態(tài)，即ξmass<0轉(zhuǎn)換為ξmass>0，則抑制了液滴溫度下降，而顆粒與氣流間僅存在對流換熱過程，缺少了凝結(jié)過程的抑制效應(yīng)，因此可見凝結(jié)效應(yīng)是三維收縮階段內(nèi)液滴溫度下降趨勢慢于顆粒的主要原因，最終導(dǎo)致構(gòu)型出口處(試驗段中心)，液滴溫度高于顆粒溫度。

2.2 相對濕度影響

首先考察相對濕度的影響，在結(jié)冰風(fēng)洞試驗中，為避免氣流濕度過低導(dǎo)致液滴大量蒸發(fā)，會采用加濕器提高風(fēng)洞回路內(nèi)氣流濕度，因此此處選取較高的相對濕度參數(shù)開展研究。

圖6(a)和圖6(b)給出了不同初始相對濕度下構(gòu)型中心線處液滴直徑和相對濕度變化曲線，從圖中可以看出：在準(zhǔn)一維傳熱階段，隨著初始相對濕度的增加，液滴尺寸的減小趨勢減慢，表明液滴蒸發(fā)效應(yīng)減弱，尤其當(dāng)濕度100%時，液滴尺寸在溫度穩(wěn)定階段不再發(fā)生變化，蒸發(fā)效應(yīng)消失；與此相反，在三維收縮階段，隨著初始相對濕度的增加，水蒸氣過飽和程度增強，液滴尺寸的增大趨勢加快，表明凝結(jié)效應(yīng)增強。

圖6(c)給出了不同初始相對濕度下構(gòu)型中心線處液滴溫度變化曲線，從圖中可以看出：在準(zhǔn)一維傳熱階段，隨著初始相對濕度增加，液滴蒸發(fā)效應(yīng)減弱，液滴溫度下降趨勢減慢(如放大圖所示)，當(dāng)液滴趨于穩(wěn)定蒸發(fā)狀態(tài)時(約在0.6 m處)，初始相對濕度越大則液滴濕球溫度越高，尤其當(dāng)相對濕度為100%時，蒸發(fā)效應(yīng)消失，液滴溫度與顆粒溫度趨于一致；在三維收縮傳熱階段，隨著初始相對濕度增加，液滴凝結(jié)效應(yīng)增強，導(dǎo)致構(gòu)型出口處(試驗段中心)液滴溫度不斷增大，尤其當(dāng)相對濕度100%時，液滴凝結(jié)效應(yīng)最強，出口液滴溫度最高，比顆粒溫度高近1 ℃。

圖7給出了不同初始相對濕度下構(gòu)型出口中心處液氣溫度差，從圖中可以看出：隨著初始相對濕度增加，液氣溫差不斷增大，其中初始相對濕度為70%時，液滴溫差與顆粒的一致，約為0.1 ℃，而當(dāng)初始相對濕度增加至100%時，液滴溫差增至約0.87 ℃，比顆粒溫度差高近0.8 ℃。由此可見，相對濕度會影響試驗段內(nèi)液滴過冷狀態(tài)，但影響程度較弱，本文計算工況下最大液氣溫差(相對濕度為100%)小于1 ℃，液滴仍認(rèn)為處于過冷狀態(tài)。

2.3 試驗段氣流速度影響

考察試驗段氣流速度的影響，根據(jù)3 m×2 m結(jié)冰風(fēng)洞典型試驗工況，選取55、77、130、164 m/s這4個試驗段氣流速度開展研究，其中各試驗段氣流速度(構(gòu)型出口氣流速度)對應(yīng)的穩(wěn)定段氣流速度(構(gòu)型入口氣流速度)分別為4.11、5.70、9.15及11.00 m/s。

圖8給出了不同試驗段氣流速度下構(gòu)型中心線處液滴直徑和相對濕度變化曲線，從圖中可以看出：在準(zhǔn)一維傳熱階段，隨著試驗段氣流速度的增加，穩(wěn)定段氣流速度增大，液滴相變時間減少，液滴尺寸的減小趨勢減慢，表明蒸發(fā)效應(yīng)減弱；在三維收縮階段，隨著試驗段氣流速度的增加，盡管液滴相變時間不斷減少，但是在三維收縮效應(yīng)的影響下，氣流靜溫顯著降低，水蒸氣過飽和度顯著增大，導(dǎo)致液滴尺寸的增大趨勢加快，進而表明凝結(jié)效應(yīng)增強。

圖9給出了不同試驗段氣流速度下構(gòu)型中心線處氣流靜溫和液滴溫度變化曲線，從圖中可以看出，在準(zhǔn)一維傳熱階段，隨著試驗段氣流速度的增加，穩(wěn)定段氣流速度增大，但是由于穩(wěn)定段速度增大幅度較小(從約4 m/s增至11 m/s)，氣流可壓縮性較弱，因此各工況下的氣流靜溫和液滴溫度變化趨勢無顯著差異，氣流靜溫接近氣流總溫，并且液滴溫度穩(wěn)定后，液滴濕球溫度也基本一致。在三維收縮階段，隨著試驗段氣流速度的增加，氣流可壓縮性顯著增強，試驗段氣流靜溫顯著下降，具體的當(dāng)試驗段氣流速度從55 m/s增至164 m/s時，試驗段氣流靜溫約從-7 ℃降至-18 ℃，因此在結(jié)冰風(fēng)洞試驗中，氣流可壓縮性導(dǎo)致的氣流靜溫變化是不可忽略的；同時，在這一階段隨著試驗段氣流速度的增加，液滴溫度的下降趨勢加快，但是相對于顆粒，由于凝結(jié)效應(yīng)的不斷增強，抑制了液滴溫度的下降趨勢，導(dǎo)致出口處(試驗段中心)液滴與顆粒溫差不斷增大。

圖10給出了不同試驗段氣流速度下構(gòu)型出口中心處液氣溫度差，從圖中可以看出：隨著試驗段氣流速度的增加，液氣溫差不斷增大，并且在相變效應(yīng)的影響下，其增大趨勢強于顆粒，其中試驗段氣流速度為55 m/s時，液滴溫差約為0.3 ℃，而顆粒的溫差則趨近于0 ℃，當(dāng)試驗段氣流速度增加至164 m/s時，液滴溫差增至約3.3 ℃，比顆粒的溫差高近2 ℃。由此可見，氣流速度會顯著影響試驗段內(nèi)液滴過冷狀態(tài)，同時相變效應(yīng)則會顯著增強影響程度，其中小尺寸液滴(40 μm≤di≤100 μm)在高風(fēng)速時(UTS≥164 m/s)將偏離過冷狀態(tài)(液氣溫差超過3 ℃)。

2.4 液滴尺寸影響

最后，考察液滴尺寸的影響，由于采用球型液滴阻力模型，不考慮液滴變形和破碎效應(yīng)，因此在We數(shù)的約束下[24]，僅選擇尺寸較小的液滴，以較好地滿足球型液滴假設(shè)。

圖12給出了不同初始液滴尺寸下構(gòu)型中心線處液滴溫度變化曲線，從圖中可以看出：在準(zhǔn)一維傳熱階段，隨著初始液滴尺寸增加，液滴熱容量增大，液滴溫度下降趨勢減慢，液滴溫度穩(wěn)定后，液滴濕球溫度無顯著變化，液滴與顆粒溫差基本一致；在三維收縮階段，隨著初始液滴尺寸增加，液滴溫度下降趨勢減慢(與文獻[10]的結(jié)論一致)，但是相對于顆粒，凝結(jié)效應(yīng)不斷減弱，促進了液滴溫度的下降趨勢，導(dǎo)致出口處(試驗段中心)液滴與顆粒溫差不斷減小。

圖13給出了不同初始液滴尺寸下構(gòu)型出口中心處液氣溫度差，從圖中可以看出：隨著初始液滴尺寸增加，液氣溫差不斷增大，但是在相變效應(yīng)的影響下，其增大趨勢弱于顆粒的溫度差，其中液滴直徑為40 μm時，液滴溫度差約為0.6 ℃，而顆粒溫度差約為0.1 ℃，當(dāng)液滴直徑增加至100 μm時，液滴溫度差增至約0.85 ℃，比顆粒的溫差低近0.05 ℃。由此可見，液滴尺寸會影響試驗段液滴過冷狀態(tài)，但相變效應(yīng)卻會減弱影響程差，其中在低風(fēng)速(UTS≤77 m/s)試驗工況下，小尺寸液滴(di≤100 μm)仍可以認(rèn)為處于過冷狀態(tài)(液氣溫差小于1 ℃)，這與文獻[10]的評估結(jié)果一致。

3 結(jié) 論

本文發(fā)展了基于Euler法的氣液兩相傳質(zhì)傳熱耦合流動計算方法，模擬了結(jié)冰風(fēng)洞中氣液兩相傳質(zhì)傳熱耦合流動過程，研究了結(jié)冰風(fēng)洞中液滴相變效應(yīng)。

1) 結(jié)冰風(fēng)洞中液滴經(jīng)歷了先蒸發(fā)后凝結(jié)兩個階段。在準(zhǔn)一維傳熱階段，液滴為蒸發(fā)狀態(tài)，蒸發(fā)效應(yīng)加快了液滴溫度下降趨勢，并且促使液滴溫度趨于濕球溫度；在三維收縮階段，液滴則從蒸發(fā)狀態(tài)轉(zhuǎn)換為凝結(jié)狀態(tài)，凝結(jié)效應(yīng)減慢了液滴溫度下降趨勢，導(dǎo)致構(gòu)型出口處(試驗段中心)液氣溫差增大，影響液滴過冷狀態(tài)。

2) 初始相對濕度增加，會減弱蒸發(fā)效應(yīng)而增強凝結(jié)效應(yīng)，進而影響試驗段內(nèi)液滴過冷狀態(tài)。但其影響程度較弱，本文計算工況下最大液氣溫差(相對濕度為100%)小于1 ℃，液滴仍認(rèn)為處于過冷狀態(tài)。

3) 氣流速度會顯著影響試驗段內(nèi)液滴過冷狀態(tài)，隨著氣流速度的增加，蒸發(fā)效應(yīng)減弱而凝結(jié)效應(yīng)增強，進而顯著增強了影響程度，其中小尺寸液滴(40 μm≤di≤100 μm)在高風(fēng)速時(UTS≥164 m/s)將偏離過冷狀態(tài)(液氣溫差超過3 ℃)。

4) 初始液滴尺寸增加，會導(dǎo)致蒸發(fā)效應(yīng)和凝結(jié)效應(yīng)均減弱，進而減弱了液滴尺寸對試驗段內(nèi)液滴過冷狀態(tài)的影響程度，對于小尺寸液滴(di≤100 μm)在低風(fēng)速(UTS≤77 m/s)試驗工況下，仍可以認(rèn)為處于過冷狀態(tài)(液氣溫差小于1 ℃)。

需要指出的是，實際結(jié)冰風(fēng)洞中，氣液兩相流動過程十分復(fù)雜，涉及動量傳遞、質(zhì)量傳遞和能量傳遞等多物理過程的強耦合作用，而本文發(fā)展的基于Euler法的氣液兩相傳質(zhì)傳熱耦合流動計算方法是基于一定假設(shè)下的簡化分析方法，因此開展計算方法驗證試驗則是下一步急需開展的工作。

此外，針對液滴直徑超過100 μm的大尺寸液滴，由于其存在變形和破碎等動力學(xué)效應(yīng)，超出了本文采用的球形液滴數(shù)學(xué)模型適用范圍，因此建立大尺寸液滴傳質(zhì)傳熱模型，則是今后需要開展的研究。

參 考 文 獻

[1] 杜雁霞， 李明， 桂業(yè)偉， 等. 飛機結(jié)冰熱力學(xué)行為研究綜述[J]. 航空學(xué)報， 2017， 38(2)： 520706.

DU Y X，LI M, GUI Y W, et al. Review of thermodynamic behaviors in aircraft icing process[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(2): 520706 (in Chinese).

[2] 桂業(yè)偉， 周志宏， 李穎暉， 等. 關(guān)于飛機結(jié)冰的多重安全邊界問題[J]. 航空學(xué)報， 2017， 38(2): 520723.

GUI Y W, ZHOU Z H, LI Y H, et al. Multiple safety boundaries protection on aircraft icing[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(2): 520723 (in Chinese).

[3] 易賢， 王斌， 李偉斌， 等. 飛機結(jié)冰冰形測量方法研究進展[J]. 航空學(xué)報， 2017, 38(2): 520700.

YI X, WANG B, LI W B, et al. Research progress on ice shape measurement approaches for aircraft icing[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(2): 520700 (in Chinese).

[4] RAGNI A, ESPOSITO B, MARRAZZO M, et al. Calibration of the CIRA IWT in the high speed configuration[C]//43th AIAA Aerospace Science Meeting and Exhibit. Reston, VA: AIAA, 2005.

[5] STEEN L E, IDE R F, van ZANTE J F, et al. NASA Glenn icing research tunnel: 2014 and 2015 cloud calibration procedures and results: NASA/TM-2015-218758 [R]. Washington, D. C. : NASA, 2015.

[6] WILLBANKS C E, SCHULZT R J. Analytical study of icing simulation for turbine engines in altitude test cells[J]. Journal of Aircraft, 1975, 12(12): 960-967.

[7] MILLER D R, ADDY H E, IDE R F. A study of large droplet ice accretions in the nasa glenn irt at near-freezing conditions: NASA TM-1996-107142-REV1 [R]. Washington, D. C. : NASA, 1996.

[8] BELLUCCI M, ESPOSITO B M, MARRAZZO M, et al. Calibration of the CIRA IWT in the low speed configuration[C]∥45th AIAA Aerospace Science Meeting and Exhibit. Reston, VA： AIAA，2007.

[9] 易賢, 馬洪林, 王開春, 等. 結(jié)冰風(fēng)洞液滴運動及傳質(zhì)傳熱特性分析[J]. 四川大學(xué)學(xué)報, 2012, 44(2): 132-135.

YI X, MA H L, WANG K C, et al. Analysis of water droplet movement and heat / mass transfer in an icing wind tunnel[J]. Journal of Sichuan University, 2012, 44(2): 132-135 (in Chinese).

[10] 郭向東，王梓旭，李明，等．結(jié)冰風(fēng)洞中液滴過冷特性數(shù)值研究[J]．航空學(xué)報, 2018, 39(1): 121254.

GUO X D, WANG Z X, LI M, et al. Numerical study of droplet supercooling in an icing wind tunnel [J]. Acta Aeronautica et Astronautica Sinica, 2018, 39(1): 121254 (in Chinese).

[11] HILL P G. Condensation of water vapor during supersonic expansion in nozzles[J]. Journal of Fluid Mechanics, 1966, 25(3): 593-620.

[12] CHENG W, LUO X S, van DONGEN M E H. On condensation-induced waves[J]. Journal of Fluid Mechanics, 2010, 651(1): 145-164.

[13] 郭向東, 黃生洪, 吳穎川, 等. 氫燃料燃燒加熱風(fēng)洞中水蒸氣相變效應(yīng)的數(shù)值研究[J]. 推進技術(shù), 2017, 38(4):932-941.

GUO X D, HUANG S H, WU Y C, et al. Numerical evaluation of water vapor phase transition effects in a hydrogen-fueled combustion-heated wind tunnel[J]. Journal of Propulsion Technology, 2017, 38(4): 932-941 (in Chinese).

[14] CAO Y, LUO X S. Numerical study on Homogeneous condensation in a vortex[J]. Procedia Engineering, 2015, 126: 607-611.

[15] LUO X S, WANG G, OLIVIER H. Parametric investigation of particle acceleration in high enthalpy conical nozzle flows for coating applications[J]. Shock Waves, 2008, 17(5): 351-362.

[16] KERSEY J, LOTH E, LANKFORD D. Effects of Evaporating Droplets on Shock Waves[J]. AIAA Journal, 2010, 48(9): 1975-1986.

[17] RUSSO E, KUERTEN J G M, van der Geld C W M, et al. Water droplet condensation and evaporation in turbulent channel flow[J]. Journal of Fluid Mechanics, 2014, 749: 666-700.

[18] SAITO T. Numerical analysis of dusty-gas flows[J]. Journal of Computational Physics, 2002, 176(1): 129-144.

[19] FUJITA A, KUROSE R, KOMORI S. Experimental study on effect of relative humidity on heat transfer of an evaporating water droplet in air flow[J]. International Journal of Heat and Mass Transfer, 2010, 36(3): 244-247.

[20] MILLER R S, HARSTAD K, BELLAN J. Evaluation of equilibrium and non-equilibrium evaporation models for many-droplet gas-liquid flow simulations[J]. International Journal of Multiphase Flow, 1998, 24(6): 1025-1055.

[21] HINDMARSH J P, RUSSEL A B, CHEN X D. Experimental and numerical analysis of the temperature transition of a suspended freezing water droplet[J]. International Journal of Heat and Mass Transfer, 2003, 46(7): 1199-1213.

[22] COBER S, BERNSTEIN B, JECK R, et al. Data and analysis for the development of an engineering standard for supercooled large drop conditions: DOT/FAA/AR-09/10[R]. Washington, D. C.: FAA, 2009.

[23] KUROSE R, FUJITA A, KOMORI S. Effect of relative humidity on heat transfer across the surface of an evaporating water droplet in air flow[J]. Journal of Fluid Mechanics, 2009, 624: 57-67.

[24] LUXFORD G. Experimental and modelling investigation of the deformation, drag and break-up of drizzle droplets subjected to strong aerodynamic forces in relation to SLD aircraft icing[D]. Bedfordshire: Cranfield University, 2005: 5-7 .