濕空氣凝結(jié)對(duì)ONERA M6機(jī)翼氣動(dòng)特性影響的數(shù)值研究

孫秀玲 ,李 亮 ,李國君

(西安交通大學(xué)能源與動(dòng)力工程學(xué)院,陜西 西安 710049)

0 引 言

目前有許多求解Navier-Stokes方程對(duì)翼型繞流進(jìn)行數(shù)值模擬的程序,它們對(duì)航空工業(yè)的發(fā)展起了巨大的推動(dòng)作用。這些計(jì)算程序基本上都假設(shè)流動(dòng)工質(zhì)為完全氣體,流動(dòng)過程中不發(fā)生相變。然而在實(shí)際的大氣飛行中,各種氣象條件下空氣中總是包含一定數(shù)量的水蒸汽,機(jī)翼周圍壓力、溫度的變化在一定條件下會(huì)引起水蒸汽發(fā)生非平衡凝結(jié)現(xiàn)象。非平衡凝結(jié)是指水蒸汽膨脹越過飽和線后并不立即凝結(jié),而是在熱力學(xué)狀態(tài)參數(shù)達(dá)到一定非平衡程度時(shí)才突然發(fā)生的凝結(jié)現(xiàn)象。圖1顯示了Tornado戰(zhàn)斗機(jī)機(jī)翼表面附近凝結(jié)形成的水霧[1]。

大氣飛行條件下1kg空氣中平均的水蒸汽含量約為0.01kg,凝結(jié)放熱量Q與CpT∞(Cp為來流空氣比熱,T∞為來流空氣溫度)的比值約為10%,這個(gè)量值的熱量加入到跨聲速氣流中,不僅會(huì)顯著改變機(jī)翼表面附近的壓力、溫度和超音速區(qū)大小等流動(dòng)參數(shù),甚至還會(huì)引起穩(wěn)定或振蕩的激波,從而導(dǎo)致機(jī)翼升力/阻力特性的顯著改變[2-5]。

另一方面,凝結(jié)放熱導(dǎo)致的激波與邊界層相互作用,引起邊界層位移厚度與動(dòng)量厚度的顯著增加,是誘發(fā)流動(dòng)分離的原因之一[6],對(duì)機(jī)翼表面附近渦的產(chǎn)生和發(fā)展產(chǎn)生影響。此外,在飛機(jī)機(jī)動(dòng)過程中,攻角、偏航角等參數(shù)的變化影響凝結(jié)過程,凝結(jié)過程反過來又會(huì)對(duì)機(jī)翼氣動(dòng)性能產(chǎn)生影響。這些因素綜合作用,使得濕空氣凝結(jié)現(xiàn)象對(duì)飛機(jī)在各種飛行姿態(tài)下的機(jī)動(dòng)性能產(chǎn)生不可忽視的影響。此外,大氣飛行中機(jī)翼等部位的結(jié)冰問題也與濕空氣凝結(jié)現(xiàn)象緊密相關(guān),濕空氣凝結(jié)是形成冰晶的重要機(jī)制之一[7]。

圖1 Tornado飛機(jī)機(jī)翼表面的凝結(jié)現(xiàn)象[1]Fig.1 Condensation around Tornado[1]

從上個(gè)世紀(jì)60年代起,不斷有研究機(jī)構(gòu)和學(xué)者對(duì)濕空氣凝結(jié)流動(dòng)問題開展研究。NASA在其8英尺風(fēng)洞中研究了濕空氣凝結(jié)對(duì)機(jī)翼氣動(dòng)性能的影響[8];Schnerr對(duì)NACA0012翼型和某圓弧翼型進(jìn)行了一系列研究[3,5,6];Zierep[9],Rusak[10]等學(xué)者從凝結(jié)起始馬赫數(shù)和凝結(jié)放熱量等方面對(duì)凝結(jié)流動(dòng)之間的相似性進(jìn)行了理論分析,取得了一些進(jìn)展。這些研究表明:跨聲速流動(dòng)與凝結(jié)過程相耦合,凝結(jié)過程與來流馬赫數(shù)、攻角、空氣濕度等參數(shù)之間存在強(qiáng)烈依賴關(guān)系,濕空氣凝結(jié)對(duì)翼型升力/阻力特性有顯著影響。盡管如此,在濕空氣凝結(jié)流動(dòng)研究方面,仍有一些重要的問題沒有解決,目前尚沒有一般性的結(jié)論指出在什么樣的條件下(大氣溫度、壓力、相對(duì)濕度、飛行速度、攻角、飛機(jī)幾何特征參數(shù)等)發(fā)生凝結(jié);另外,凝結(jié)影響機(jī)翼氣動(dòng)性能的一般規(guī)律也并不清楚。解決這些問題,對(duì)提高飛機(jī)在復(fù)雜大氣條件和各種飛行姿態(tài)下的氣動(dòng)性能具有重要意義。

本文發(fā)展了濕空氣非平衡凝結(jié)流動(dòng)數(shù)學(xué)物理模型和計(jì)算方法,并對(duì)ONERA M6機(jī)翼周圍的濕空氣凝結(jié)流動(dòng)進(jìn)行了分析。

1 濕空氣凝結(jié)流動(dòng)數(shù)學(xué)物理模型

將濕空氣凝結(jié)流視為一個(gè)兩相系統(tǒng):氣相由干空氣和水蒸汽即濕空氣構(gòu)成,凝結(jié)相為彌散在濕空氣中所有凝結(jié)物的集合。在大氣飛行條件下,由于水蒸汽相變發(fā)生的溫度低于三相點(diǎn)溫度,水蒸汽在跨聲速條件下非平衡相變的產(chǎn)物可能是液態(tài)微小水滴或固態(tài)冰晶。相變過程影響機(jī)翼氣動(dòng)特性的本質(zhì)是相變潛熱對(duì)氣流加熱產(chǎn)生的影響。文獻(xiàn)[11]指出,為了計(jì)算這種影響,可以假設(shè)相變產(chǎn)物為液態(tài)水,并在成核模型中對(duì)平面水的表面張力進(jìn)行適當(dāng)修正來計(jì)算成核過程中微小水滴的表面張力。因此本文假設(shè)濕空氣相變產(chǎn)物為微小水滴。

由于跨聲速濕空氣流動(dòng)中非平衡凝結(jié)形成的水滴直徑在1μm以下,水滴和空氣流之間的速度滑移可以忽略不計(jì),因而可以建立如下歐拉/歐拉型的數(shù)值模型。描述氣相即濕空氣流動(dòng)的質(zhì)量、動(dòng)量及能量守恒方程為

其中,ρa(bǔ)v為濕空氣密度,ρ為濕空氣與水滴混合物的密度,U為速度矢量,П為應(yīng)力張量,pav為濕空氣靜壓,Eav為濕空氣總能,ht為濕空氣總焓,hfg為凝結(jié)潛熱。引入源項(xiàng)-ρ﹒m、-ρ﹒mU、-ρ﹒m(ht-hfg)考慮流動(dòng)過程中氣液兩相之間的相互作用對(duì)濕空氣流動(dòng)的影響,這里﹒m為質(zhì)量凝結(jié)速率,其表達(dá)式為:

其中N為水滴數(shù)目,r為水滴半徑,Y為凝結(jié)出的液態(tài)水在濕空氣與水滴混合物中的質(zhì)量百分?jǐn)?shù),J為成核率,dr/dt為水滴生長(zhǎng)率,ρl為液態(tài)水密度,rc為臨界半徑。表達(dá)式(4)中兩項(xiàng)分別考慮了成核和水滴生長(zhǎng)對(duì)質(zhì)量凝結(jié)速率的貢獻(xiàn)。在求得液態(tài)水質(zhì)量百分?jǐn)?shù)和水滴數(shù)后,水滴半徑由下式確定

為了使方程組封閉,補(bǔ)充濕空氣的狀態(tài)方程

如前所述,由于假設(shè)凝結(jié)出的水滴與氣相之間無速度滑移,因此對(duì)凝結(jié)出的液相只需描述水滴數(shù)目N和液態(tài)水質(zhì)量份額Y。凝結(jié)形成的水滴實(shí)際上是具有不同半徑的水滴族,如果對(duì)這些水滴族一一計(jì)算,會(huì)導(dǎo)致計(jì)算過程極其復(fù)雜。依據(jù)兩相流理論,凝結(jié)形成的水滴遵循統(tǒng)計(jì)規(guī)律[13],因而可以采用體積平均的水滴數(shù)N進(jìn)行計(jì)算。

另外在濕空氣流中,隨著水蒸汽不斷凝結(jié)為水滴,濕空氣中的含濕量也在不斷減少,水蒸汽的分壓隨之減小。為了計(jì)算水蒸汽的當(dāng)?shù)胤謮毫v,還需要知道濕空氣中含濕量d的分布。令

由含濕量定義可得ρv=ρa(bǔ)vD,其中ρv為水蒸汽密度。由守恒定律可以得到描述N,Y及D分布的控制方程為

成核率J和水滴生長(zhǎng)率dr/dt由經(jīng)典成核理論計(jì)算,本文采用的表達(dá)式為[12]

其中,ω和ν為兩個(gè)修正系數(shù),qc為凝結(jié)系數(shù),σ為水滴表面張力,mm為水分子質(zhì)量,k為熱傳導(dǎo)系數(shù),T為溫度,Kn為玻爾茲曼常數(shù),Prg為普朗特?cái)?shù),過冷度ΔT由水蒸汽的當(dāng)?shù)胤謮毫v和溫度T確定。

其中Ts為水蒸汽當(dāng)?shù)貕毫?pv對(duì)應(yīng)的飽和溫度。這樣,由方程(1-3),方程(5、6)以及方程(8-10)就構(gòu)成了求解跨聲速濕空氣非平衡凝結(jié)流動(dòng)的全部方程。

需要說明的是,Schnerr[11]對(duì)濕空氣凝結(jié)流動(dòng)建立的模型中采用了面積平均的水滴半徑,相應(yīng)需要求解四個(gè)偏微分方程來描述水滴的大小和數(shù)目的分布。本文采用體積平均的水滴數(shù)N水滴半徑r,由此由方程(8、9)即可確定水滴的分布,從而減少了計(jì)算量。由兩相流理論[13],水滴的分布滿足一定的統(tǒng)計(jì)規(guī)律,因此采用面積平均或體積平均的水滴半徑本質(zhì)上是相同的,計(jì)算結(jié)果一致,這一點(diǎn)由本文第3節(jié)的對(duì)比也可看到。

2 濕空氣凝結(jié)流動(dòng)數(shù)值計(jì)算方法

在濕空氣兩相非平衡凝結(jié)流動(dòng)計(jì)算中,氣液兩相間的耦合方式及解法也是一個(gè)重要問題。二維氣液兩相流動(dòng)的控制方程組聯(lián)立求解時(shí)由于方程組的雅可比系數(shù)矩陣剛性較大,不易收斂,對(duì)三維流動(dòng),這個(gè)問題更加突出。本文處理方法是液相和氣相各自作為相對(duì)獨(dú)立的開口系,依靠在控制方程中添加考慮兩相間相互作用的源項(xiàng)實(shí)現(xiàn)耦合計(jì)算。兩相耦合關(guān)系如圖2所示。氣相即濕空氣流動(dòng)求解模塊向液相求解模塊傳遞濕空氣密度ρa(bǔ)v、壓力pav、溫度T及速度矢量U,液相模塊由凝結(jié)模型(11、12)計(jì)算水滴成核率和生長(zhǎng)率,并求解方程(8-10),由此得到液相凝結(jié)及水滴的分布情況;同時(shí),液相模塊向氣相流動(dòng)求解模塊傳遞質(zhì)量凝結(jié)速率﹒m和凝結(jié)潛熱hfg,從而可以計(jì)算方程(1-3)中的源項(xiàng)。采用這種方式,可以在大量已有的單相流動(dòng)計(jì)算程序基礎(chǔ)上發(fā)展?jié)窨諝饽Y(jié)流動(dòng)的計(jì)算程序。

圖2 兩相耦合方式Fig.2 Interphase coupling

在機(jī)翼的濕空氣非平衡凝結(jié)繞流計(jì)算中,一方面機(jī)翼表面附近可能存在氣動(dòng)激波,另一方面,由于非平衡凝結(jié)發(fā)生的突然性,凝結(jié)流場(chǎng)中液相參數(shù)變化劇烈,而且凝結(jié)過程對(duì)當(dāng)?shù)亓鲃?dòng)參數(shù)的變化很敏感。因此,在數(shù)值求解中需要采用高分辨率的差分格式。在本文計(jì)算中,對(duì)氣液兩相方程組的求解均采用了Harten提出的TVD格式[14]。

這樣氣相方程(1-3)和液相方程(8、9、14)均可在直角坐標(biāo)系下分別寫為如下矢量形式

對(duì)于NS方程計(jì)算,紊流模型采用SA模型,可以與氣相方程(1-3)統(tǒng)一寫成式(15)的形式。直角坐標(biāo)系下三個(gè)方向的離散方法相同。以方向?yàn)槔?令A(yù)=?F/?Q=RΛ R-1,其中 R 為A 的特征矩陣,R-1為R的逆陣,Λ為A的特征值構(gòu)成的對(duì)角陣。方程(15)離散為,

求解中采用顯式時(shí)間推進(jìn),紊流模型采用SA模型。氣相方程組求解時(shí),遠(yuǎn)場(chǎng)邊界給定濕空氣的壓力、溫度、相對(duì)濕度或含濕量,翼面給定粘性固體邊界。液相方程求解中,遠(yuǎn)場(chǎng)及翼面邊界條件均由內(nèi)點(diǎn)外推。求解中為了保證計(jì)算穩(wěn)定,濕空氣流動(dòng)計(jì)算的CFL數(shù)取為干空氣流動(dòng)計(jì)算時(shí)的0.1倍。

3 濕空氣凝結(jié)流動(dòng)模型驗(yàn)證

首先對(duì)NACA0012翼型周圍的濕空氣凝結(jié)流動(dòng)進(jìn)行了無粘計(jì)算。由于缺乏濕空氣凝結(jié)流動(dòng)中翼型表面壓力系數(shù)的試驗(yàn)數(shù)據(jù),因此計(jì)算結(jié)果與文獻(xiàn)[11]中的結(jié)果進(jìn)行了比較,以驗(yàn)證本文所建立的模型。計(jì)算中的流動(dòng)條件為:弦長(zhǎng)c=0.1m,壓力p∞=65600Pa,溫度 T∞=259K,馬赫數(shù) M∞=0.8,相對(duì)濕度 φ∞=50%,攻角 α=0°。遠(yuǎn)場(chǎng)邊界距翼型約30倍弦長(zhǎng)。

計(jì)算得到的翼型表面壓力系數(shù)分布如圖3所示?？梢钥吹?本文結(jié)果與文獻(xiàn)[11]中的結(jié)果吻合良好。在相對(duì)濕度φ∞=50%的流動(dòng)中,由于水蒸汽凝結(jié)放熱對(duì)氣流加熱,激波前馬赫數(shù)明顯減小,激波位置也向上游移動(dòng),導(dǎo)致翼型表面壓力系數(shù)發(fā)生顯著變化。

圖3 NACA0012翼型表面壓力系數(shù)分布Fig.3 Pressure coefficient of NACA0012 airfoil

4 機(jī)翼濕空氣凝結(jié)流動(dòng)分析

利用本文發(fā)展的數(shù)值方法對(duì)ONERA M6機(jī)翼[15]周圍的濕空氣凝結(jié)流動(dòng)進(jìn)行了分析。計(jì)算中流動(dòng)條件為壓力 p∞=129000Pa,溫度 T∞=293.15K,馬赫數(shù)Ma∞ =0.84,攻角 α=3.06°,相對(duì)濕度 φ∞ =50%,計(jì)算網(wǎng)格總數(shù)約為30萬。

沿ONERA M6機(jī)翼展向三個(gè)截面的壓力系數(shù)如圖4所示。干空氣(φ∞=0%)流動(dòng)計(jì)算清晰地捕獲到了機(jī)翼表面的激波;而在濕空氣(φ∞=50%)流動(dòng)中,由于凝結(jié)潛熱對(duì)跨聲速氣流的加熱作用,導(dǎo)致激波消失。另外,與干空氣(φ∞=0%)流動(dòng)相比,濕空氣(φ∞=50%)凝結(jié)流動(dòng)中機(jī)翼下表面壓力系數(shù)變化不大,而上表面壓力系數(shù)則有顯著變化。這是由于下表面氣流膨脹速率低,凝結(jié)過程不明顯,凝結(jié)放熱對(duì)流動(dòng)的影響較小;而上表面的跨聲速流動(dòng)中,由于發(fā)生明顯的濕空氣凝結(jié)放熱,對(duì)跨聲速流動(dòng)造成了顯著影響,從而使機(jī)翼上表面的壓力系數(shù)發(fā)生明顯變化。這一點(diǎn)也可以從下文圖7機(jī)翼表面附近凝結(jié)出的液態(tài)水質(zhì)量份額看到。

為了分析凝結(jié)放熱對(duì)氣流參數(shù)的影響,圖5給出干空氣流動(dòng)和濕空氣凝結(jié)流動(dòng)條件下機(jī)翼表面附近的馬赫數(shù)分布?？梢钥吹?凝結(jié)流動(dòng)中機(jī)翼上表面激波前的馬赫數(shù)減小,并且沿著翼展方向減小的幅度越來越大,這與圖4中兩種流動(dòng)條件下壓力系數(shù)的變化相對(duì)應(yīng)。

圖4 ONERA M6機(jī)翼表面壓力系數(shù)分布Fig.4 Pressure coefficient of ONERA M6 wing

圖6給出ONERA M6機(jī)翼上表面壓力等值線分布。可以看到干空氣流動(dòng)中上表面附近形成了清晰的λ形激波;而在濕空氣凝結(jié)流動(dòng)中,由于凝結(jié)放熱使馬赫數(shù)減小,λ形激波結(jié)構(gòu)消失,激波位置向下游移動(dòng),并且激波強(qiáng)度也明顯減弱。

機(jī)翼表面壓力系數(shù)的變化是由于空氣中水蒸汽凝結(jié),相變潛熱對(duì)跨聲速氣流加熱,導(dǎo)致機(jī)翼表面附近流動(dòng)參數(shù)發(fā)生改變的結(jié)果。因此機(jī)翼氣動(dòng)特性的改變與機(jī)翼表面附近的凝結(jié)過程直接相關(guān)。圖7給出了65%翼展位置機(jī)翼表面附近凝結(jié)出的液態(tài)水的質(zhì)量份額Y,沿翼展其它位置處Y的分布與此類似。在機(jī)翼上表面,由于凝結(jié)發(fā)生在前緣,整個(gè)上表面附近的液態(tài)水質(zhì)量份額約為1%;下表面凝結(jié)發(fā)生的位置向后移動(dòng),在50%弦長(zhǎng)以后才可以凝結(jié),機(jī)翼下表面附近的液態(tài)水質(zhì)量份額約為0.5%。

圖5 ONERA M6機(jī)翼表面附近馬赫數(shù)等值線分布(等值線間隔0.05)Fig.5 Contours of Mach number of ONERAM6 wing(ΔM=0.05)

圖6 ONERA M6機(jī)翼上表面壓力等值線分布Fig.6 Upper surface pressure of ONERA M6 wing

圖7 65%展向位置處ONERA M6表面附近液態(tài)水質(zhì)量份額Y(Δ Y=0.002)Fig.7 Mass fraction of waterY around ONERA M6 wing at 65%span(Δ Y=0.002)

5 結(jié)論與展望

濕空氣中的水蒸汽凝結(jié)對(duì)機(jī)翼氣動(dòng)特性會(huì)造成顯著影響。本文建立了濕空氣非平衡凝結(jié)流動(dòng)的數(shù)值模型,并對(duì)ONERA M6機(jī)翼進(jìn)行了濕空氣非平衡凝結(jié)流動(dòng)的初步研究。在攻角為3.06°和空氣相對(duì)濕度50%時(shí)的ONERA M6機(jī)翼濕空氣凝結(jié)流動(dòng)與干空氣流動(dòng)相比,機(jī)翼表面壓力系數(shù)有顯著變化。造成機(jī)翼氣動(dòng)特性顯著變化的原因在于:濕空氣中水蒸汽凝結(jié)放熱對(duì)跨聲速氣流加熱,導(dǎo)致機(jī)翼表面附近的流速、壓力與流場(chǎng)結(jié)構(gòu)發(fā)生了顯著變化。

濕空氣凝結(jié)流動(dòng)研究涉及復(fù)雜的相變過程,如果考慮大氣中塵埃、離子等各種微小懸浮顆粒對(duì)凝結(jié)過程的影響,模型將更為復(fù)雜。但濕空氣凝結(jié)對(duì)機(jī)翼氣動(dòng)特性的影響顯著,有必要在試驗(yàn)和計(jì)算方面進(jìn)一步開展深入的研究工作。

[1]http://www.steehouwer.com[DB/OL].

[2]CAM PBELL J F,CHAMBERS J R,AND RUMSEY C L.Oberservation of airplane flow fields by natural condensation effects[J].J.Aircraft,1989,26(7):593-604.

[3]SCHNERR G H,MUNDINGER G.Similarity,drag and lift in transonic flow with given internal heat addition[J].Eur.J.Mech.,B/Fluids,1993,12(5):597-612.

[4]LI L,SUN X,FENG Z and LI G.Transonic flow of moist air around a NACA 0012 airfoil with non-equilibrium condensa-tion[J].Progress in Natural Science,2005,15(9):838-842.

[5]SCHNERR G H.Transonic aerodynamics including strong effects from heat addition[J].ComputersFluids,1993,22(2-3):103-116.

[6]SCHNERR G H,LI P.Shock wave/boundary layer interaction with heat addition by non-equilibrium phase transition[J].Int.J.Multiphase Flow,1993,19(5):737-749.

[7]WILLIAM L W,GLENN G,THOMAS J H,et al.AircRAFT-PRODUCED ICE PArticles(APIPs):additional results and further insights[J].J.Applied Meteorology,2003,42(5):640-652.

[8]JORDAN F L.Investigation at near-sonic speed of some effects ofhumidity on the longitudinal aerodynamic characteris-tics of a NASA supercritical wing research airplane model[R].NASA-TM-X-2618,1972.

[9]ZIEREP J,LIN S.Bestimmung des kondensation-beginns bei der entspannung ǜfeuchter luft in berschalldǜsen[J].For-schung im Ingenieurwesen,1967,33(6):169-172.

[10]RUSAK Z,LEE J C.Transonic flow of moist air around a thin airfoil with non-equilibrium and homogeneous condensa-tion[J].J.Fluid Mechanics,2000,43(1):173-199.

[11]SCHNERR G H,DOHRMANN U.Transonic flow around airfoils with relaxation and energy supply by homogeneous condensation[J].AIAA Journal,1990,28(7):1187-1193.

[12]GUHA A.,YOUNG J.B.Time-marching prediction of unsteady condensation phenomena due to supercritical heat addi-tion[J].IMechE,1991,C423/057.

[13]MOORE M J,SIEVERDING C H.Two-phase steam flow in turbines and separators[M].Hemisphere Publishing Corpora-tion,1976.

[14]HARTEN A.High resolution schemes for hypersonic conservation laws[J].Journalof Computational Physics,1983,49:357-393.

[15]http://www.grc.nasa.gov/WWW/wind/valid/archive.html[DB/OL].