毛曉東，卜雪琴，趙國(guó)昌，王建明

(1.沈陽(yáng)航空航天大學(xué) 航空航天工程學(xué)部(院)，沈陽(yáng) 110136； 2.北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京 100191)

?

基于UDS框架的水滴撞擊特性數(shù)值計(jì)算方法

毛曉東1，卜雪琴2，趙國(guó)昌1，王建明1

(1.沈陽(yáng)航空航天大學(xué) 航空航天工程學(xué)部(院)，沈陽(yáng) 110136； 2.北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京 100191)

摘要：為了預(yù)測(cè)二維/三維復(fù)雜表面在可壓/不可壓流動(dòng)狀態(tài)下的水滴撞擊特性，提出了一種基于UDS框架的水滴撞擊特性數(shù)值計(jì)算方法。水滴相和空氣相認(rèn)為單向耦合，空氣相流場(chǎng)由N-S方程預(yù)先求得。水滴相控制方程基于FLUENT中UDS輸運(yùn)方程框架，自定義編程迭代求解。為了排除壁面處水滴，給出一種水滴撞擊壁面邊界條件處理方法。引入具有自適應(yīng)性的人工數(shù)值擴(kuò)散項(xiàng)以增加計(jì)算穩(wěn)定性，避免局部水滴容積分?jǐn)?shù)異常大導(dǎo)致解的奇異。驗(yàn)證計(jì)算了二維圓柱和三維球面，通過(guò)與公開(kāi)文獻(xiàn)中試驗(yàn)數(shù)據(jù)的對(duì)比，證明了本文計(jì)算方法的準(zhǔn)確性和有效性，能很好的解決復(fù)雜表面在各種流動(dòng)狀態(tài)下的水滴撞擊特性預(yù)測(cè)。

關(guān)鍵詞：UDS框架；水滴撞擊；結(jié)冰；二相流；數(shù)值計(jì)算

當(dāng)飛機(jī)在某些特殊的氣象條件下飛行時(shí)，飛機(jī)表面以及航空發(fā)動(dòng)機(jī)相關(guān)部件表面，均有可能發(fā)生結(jié)冰現(xiàn)象。結(jié)冰將改變飛機(jī)繞流流場(chǎng)以及發(fā)動(dòng)機(jī)進(jìn)氣道流場(chǎng)特性，影響飛機(jī)的操縱性和穩(wěn)定性以及發(fā)動(dòng)機(jī)正常工作特性，危害飛行安全。如果冰層脫落并吸入發(fā)動(dòng)機(jī)進(jìn)氣道，可能直接導(dǎo)致發(fā)動(dòng)機(jī)停機(jī)，最終發(fā)生機(jī)毀人亡的嚴(yán)重事故[1]。水滴撞擊特性計(jì)算，是飛機(jī)及發(fā)動(dòng)機(jī)結(jié)冰研究的首要內(nèi)容，是結(jié)冰冰形預(yù)測(cè)及防除冰系統(tǒng)設(shè)計(jì)的前提和基礎(chǔ)。

水滴撞擊特性計(jì)算方法主要分為拉格朗日法和歐拉法。拉格朗日法把水滴看作離散相，在求解空氣相繞流流場(chǎng)的基礎(chǔ)上，采用有限差分法追蹤水滴在流場(chǎng)中的運(yùn)動(dòng)軌跡，從而得到水滴撞擊極限及局部水滴收集系數(shù)等參數(shù)[2-3]。該方法對(duì)二維或幾何形狀簡(jiǎn)單的表面比較簡(jiǎn)便，但對(duì)于三維或復(fù)雜外形，由于粒子釋放位置較難確定，處理起來(lái)比較復(fù)雜。歐拉法把水滴看作連續(xù)相介質(zhì)，通過(guò)求解水滴相的連續(xù)方程和動(dòng)量方程，得到計(jì)算域各網(wǎng)格節(jié)點(diǎn)處的水滴容積分?jǐn)?shù)和速度矢量。相對(duì)于拉格朗日法，歐拉法無(wú)需確定粒子的釋放位置，特別適合如三維機(jī)翼[4-5]、多段翼[6]、發(fā)動(dòng)機(jī)進(jìn)氣道[7]、發(fā)動(dòng)機(jī)整流帽罩[8]和發(fā)動(dòng)機(jī)整流支板[9]等處的三維復(fù)雜表面的水滴撞擊特性計(jì)算。目前FLUENT商用軟件中的歐拉模型求解器被廣泛應(yīng)用于水滴撞擊特性的數(shù)值求解。歐拉模型是非常復(fù)雜的多相流模型，采用完全耦合，計(jì)算量大，收斂速度慢，容易發(fā)散。而最重要的是它不支持模擬可壓縮流動(dòng)所需的壓力遠(yuǎn)場(chǎng)邊界條件，因此直接影響到結(jié)冰模擬中馬赫數(shù)超過(guò)0.3的情況。本文提出了一種基于FLUENT中UDS框架的水滴撞擊特性計(jì)算方法，該方法計(jì)算快速、準(zhǔn)確，并且普遍適用于可壓縮流及不可壓縮流。

1數(shù)學(xué)模型

1.1模型假設(shè)

對(duì)于空氣相和過(guò)冷水滴相所組成的兩項(xiàng)流動(dòng)，空氣相通過(guò)阻力、湍流效應(yīng)等方式影響水滴相。但由于水滴相的微粒濃度很小，故認(rèn)為水滴對(duì)空氣相無(wú)影響，即單向耦合[7]。在建立水滴相控制方程前，先作如下假設(shè)：(1)水滴為大小統(tǒng)一的球形剛體，忽略變形和破裂；(2)水滴之間沒(méi)有碰撞和融合，在撞擊到壁面后不發(fā)生反彈和飛濺；(3)忽略水滴與空氣間的傳熱和傳質(zhì)；(4)忽略水滴相動(dòng)量方程中的粘性項(xiàng)和壓力項(xiàng)；(5)僅考慮水滴所受到的阻力和重力，忽略非穩(wěn)態(tài)力；(6)阻力認(rèn)為是穩(wěn)態(tài)的；(7)不考慮空氣對(duì)水滴的湍流作用。

1.2控制方程

由于單向耦合，故空氣流場(chǎng)仍由Navier-Stokes方程描述，且可以預(yù)先計(jì)算收斂，作為水滴相控制方程求解的已知條件。水滴的連續(xù)性方程和動(dòng)量方程如式(1)和式(2)所示。

(1)

(2)

式中，ρ為水滴密度，α為水滴容積分?jǐn)?shù)，u為水滴速度矢量，ua為空氣速度矢量，F(xiàn)為作用于水滴上的除阻力外的外力，K是空氣和水滴的動(dòng)量交換系數(shù)，由式(3)定義。

(3)

其中，μa為空氣動(dòng)力粘度，dp為水滴直徑，f是阻力函數(shù)，由于假設(shè)水滴為球形剛體，阻力函數(shù)如式(4)。

(4)

其中，CD為阻力系數(shù)，Re為相對(duì)雷諾數(shù)，分別由式(5)和式(6)給出。

(5)

(6)

式中，ρa(bǔ)為空氣密度。

2基于FLUENT的UDS框架

2.1UDS輸運(yùn)方程

UDS(User Defined Scalar)，即用戶自定義標(biāo)量。FLUENT商用軟件提供了求解UDS輸運(yùn)方程的相關(guān)功能，只需待求解方程滿足輸運(yùn)方程一般形式，即可通過(guò)宏定義進(jìn)行迭代求解。FLUENT中UDS輸運(yùn)方程的一般形式如式(7)。

(7)

式中，φk為用戶自定義的一個(gè)標(biāo)量，F(xiàn)i，Гk，Sφk分別為輸運(yùn)方程的對(duì)流項(xiàng)、擴(kuò)散項(xiàng)和源項(xiàng)系數(shù)。以三維問(wèn)題為例，根據(jù)水滴相控制方程，需要定義4個(gè)UDS，分別是水滴容積分?jǐn)?shù)α，水滴速度分量u，v，w。對(duì)照式(1)、式(2)和式(7)，即可得到4個(gè)UDS輸運(yùn)方程所對(duì)應(yīng)的系數(shù)如表1。

表1　輸運(yùn)方程系數(shù)表

根據(jù)以上定義，利用FLUENT宏函數(shù)進(jìn)行編程和調(diào)用，對(duì)流項(xiàng)、擴(kuò)散項(xiàng)和源項(xiàng)的定義宏分別為DEFINE_UDS_FLUX、DEFINE_DIFFUSIVITY、DEFINE_SOURCE。以預(yù)先計(jì)算收斂的空氣相流場(chǎng)作為已知條件，耦合水滴相邊界條件，進(jìn)行迭代求解，即可獲得各網(wǎng)格點(diǎn)的水滴容積分?jǐn)?shù)及水滴速度矢量。

2.2邊界條件處理

飛機(jī)和發(fā)動(dòng)機(jī)結(jié)冰模擬中，一般為繞流流動(dòng)，因此邊界條件有遠(yuǎn)場(chǎng)邊界條件及壁面邊界條件。對(duì)于遠(yuǎn)場(chǎng)邊界，水滴容積分?jǐn)?shù)α∞可由液態(tài)水含量計(jì)算得到，如式(8)。

(8)

遠(yuǎn)場(chǎng)處認(rèn)為水滴根隨空氣一同運(yùn)動(dòng)，無(wú)相對(duì)速度，故水滴速度即為空氣速度。

u∞=ua，∞

(9)

對(duì)于壁面邊界，當(dāng)過(guò)冷水滴撞擊到壁面處時(shí)，將不會(huì)繼續(xù)以水滴形態(tài)留存于流體計(jì)算域中，而是會(huì)發(fā)生結(jié)冰，或者變成水膜、水珠等方式附著在壁面處。由假設(shè)忽略水滴撞擊壁面處的反彈和飛濺，則固壁面對(duì)于水滴而言，相當(dāng)于一個(gè)單向出口，只許流出計(jì)算域，而不能流入計(jì)算域。因此必需將撞擊到壁面處的水滴予以排除，否則會(huì)產(chǎn)生非物理解。

水滴在物體表面的繞流，根據(jù)水滴速度矢量方向與表面法線向量的關(guān)系，可將壁面分為撞擊區(qū)和非撞擊區(qū)，如圖1所示。對(duì)于撞擊區(qū)的面網(wǎng)格，其UDS輸運(yùn)方程的對(duì)流項(xiàng)按表1所示進(jìn)行計(jì)算，即允許水滴流出計(jì)算域。而對(duì)于非撞擊區(qū)的面網(wǎng)格，對(duì)流項(xiàng)設(shè)定為返回0值，即表示此處沒(méi)有水滴流入計(jì)算域。

圖1　物體表面水滴速度示意圖

2.3數(shù)值擴(kuò)散項(xiàng)

計(jì)算中發(fā)現(xiàn)，在某些局部區(qū)域，一般在繞流物體后方上下兩股氣流交匯處，會(huì)出現(xiàn)水滴容積分?jǐn)?shù)異常增大的情況，體現(xiàn)在某個(gè)網(wǎng)格節(jié)點(diǎn)的水滴容積分?jǐn)?shù)值可能比它周?chē)W(wǎng)格值大100倍以上。Tony把這種現(xiàn)象解釋為兩個(gè)沿不同方向運(yùn)動(dòng)的水滴顆粒交匯后，變成以同一個(gè)速度沿同一方向運(yùn)動(dòng)了，如圖2所示[8]。這使得交匯點(diǎn)的密度出現(xiàn)了一個(gè)無(wú)限大的階躍，這個(gè)無(wú)限大的密度脈沖，就可能帶來(lái)解的奇異。

圖2　水滴交匯后的情況

為了及時(shí)清除水滴容積分?jǐn)?shù)異常的區(qū)域，在水滴連續(xù)性方程(1)中加入人工數(shù)值擴(kuò)散項(xiàng)。加入數(shù)值擴(kuò)散項(xiàng)后的連續(xù)性方程變?yōu)槭?10)。

(10)

數(shù)值擴(kuò)散系數(shù)Г的選取具有一定的自適應(yīng)性，在正常區(qū)域，Г應(yīng)趨于0，以保證解的精度；在水滴容積分?jǐn)?shù)異常區(qū)域，Г應(yīng)足夠大，使之可以迅速清除奇異區(qū)。基于以上要求，本文構(gòu)造數(shù)值擴(kuò)散系數(shù)如式(11)、式(12)所示。

Γ=a(eb-1.0)

(11)

(12)

式中，下標(biāo)P代表當(dāng)前網(wǎng)格，Ni代表與之相鄰網(wǎng)格，b表征周?chē)W(wǎng)格和當(dāng)前網(wǎng)格水滴容積分?jǐn)?shù)的偏離程度，a是比例因子，可根據(jù)經(jīng)驗(yàn)進(jìn)行調(diào)整。由式(11)可知，當(dāng)前網(wǎng)格容積分?jǐn)?shù)正常時(shí)，b值趨于0，擴(kuò)散系數(shù)也趨于0。而對(duì)于異常區(qū)域，偏離程度越大，b值越大，指數(shù)型函數(shù)也可以確保擴(kuò)散系數(shù)足夠大，從而達(dá)到快速清除奇異區(qū)的目的。綜上所述，2.1節(jié)中表1所列UDS(α)方程的擴(kuò)散項(xiàng)系數(shù)Гk應(yīng)變?yōu)棣学ぁ?/p>

3算例驗(yàn)證

為了驗(yàn)證建立模型的準(zhǔn)確性，與相關(guān)試驗(yàn)結(jié)果進(jìn)行對(duì)比。由于本文提出的方法同時(shí)適用于二維及三維物體，故分別選取二維圓柱及三維球面作為對(duì)比算例。對(duì)比參數(shù)為局部水滴收集系數(shù)β，其計(jì)算式為式(13)。

(13)

式中，u為壁面處水滴速度矢量，n為壁面法向向量，V∞為遠(yuǎn)場(chǎng)速度，αn為水滴的相對(duì)容積分?jǐn)?shù)，定義為壁面處水滴容積分?jǐn)?shù)與遠(yuǎn)場(chǎng)水滴容積分?jǐn)?shù)之比。由模型假設(shè)可知，它基于單一水滴直徑，實(shí)際上，云層中的水滴不可能是同一種直徑大小。為了提高模擬精度，使用呈某種分布的若干種直徑來(lái)代替實(shí)際情況，最常見(jiàn)的是Langmuir-D分布[9]，由7種水滴直徑按一定比例構(gòu)成。計(jì)算中首先針對(duì)每種水滴直徑(dp)i求解一次，得到相應(yīng)水滴直徑的局部水滴收集系數(shù)βi，最后由式(14)計(jì)算得到綜合收集系數(shù)。

(14)

式中pi為第i種直徑水滴所占液態(tài)水含量比重。

3.1二維圓柱

二維圓柱作為典型的二維情況，在多篇文獻(xiàn)中都作為驗(yàn)證模型使用[10-11]。本文計(jì)算選取與試驗(yàn)相同的條件：圓柱直徑10.16 cm，自由流速80 m/s，入口空氣密度1.097 kg/m3，環(huán)境壓力89 867 Pa。計(jì)算7種水滴直徑下的收集系數(shù)，平均容積直徑(MVD)為16 μm，滿足Langmuir-D分布。圖3顯示了7種水滴直徑下的局部水滴收集系數(shù)以及由式(14)計(jì)算得到的Langmuir-D分布下的綜合收集系數(shù)，可以看出水滴直徑越大，局部水滴收集系數(shù)峰值和撞擊極限均增加。這是由于水滴直徑越大，質(zhì)量越大，慣性力越大，受氣流影響越小，因此更容易撞擊到物體表面。

圖3　二維圓柱不同直徑水滴局部收集系數(shù)

圖4　二維圓柱水滴收集系數(shù)預(yù)測(cè)值與試驗(yàn)結(jié)果對(duì)比

圖4顯示了二維圓柱局部水滴收集系數(shù)的計(jì)算結(jié)果與試驗(yàn)值的比較。圖中紅色圓圈表示使用單一平均容積直徑MVD的計(jì)算結(jié)果，藍(lán)色三角表示采用相應(yīng)Langmuir-D分布的7種水滴直徑的綜合收集系數(shù)，陰影部分為試驗(yàn)重復(fù)性范圍。由圖4可知，2種計(jì)算結(jié)果與試驗(yàn)值均吻合較好，單一MVD計(jì)算得到的局部水滴收集系數(shù)峰值較Langmuir-D分布的綜合收集系數(shù)略大，但基本都處于試驗(yàn)范圍之內(nèi)。而對(duì)于撞擊極限預(yù)測(cè)，使用Langmuir-D分布的預(yù)測(cè)結(jié)果更接近試驗(yàn)值，使用單一MVD水滴直徑的撞擊極限預(yù)測(cè)值偏小。

3.2三維球面

第二個(gè)驗(yàn)證算例選自Bidwell等人[12]對(duì)三維球面所做的水滴撞擊試驗(yàn)。試驗(yàn)對(duì)象為直徑為15.04 cm的球體，自由流速度75 m/s，環(huán)境靜溫度7 ℃，環(huán)境壓力95 840 Pa，水滴平均容積直徑18.6 μm，滿足Langmuir-D分布。圖5顯示了不同直徑下的局部水滴收集系數(shù)，以及據(jù)此計(jì)算出的Langmuir-D分布下綜合收集系數(shù)。圖6給出了預(yù)測(cè)值與試驗(yàn)數(shù)據(jù)的對(duì)比，可以看出本文的預(yù)測(cè)結(jié)果與試驗(yàn)數(shù)據(jù)非常吻合，驗(yàn)證了本文方法在三維情況下的有效性。

圖5　三維球面不同直徑水滴局部收集系數(shù)

圖6　三維球面水滴收集系數(shù)預(yù)測(cè)值與試驗(yàn)結(jié)果對(duì)比

4結(jié)論

本文提出了一種基于UDS框架的水滴撞擊特性計(jì)算方法，通過(guò)對(duì)水滴撞擊邊界進(jìn)行處理從而將撞擊水從計(jì)算域中排除，引入數(shù)值擴(kuò)散項(xiàng)避免了解的奇異性。通過(guò)二維及三維算例的試驗(yàn)對(duì)比，驗(yàn)證了本文方法的準(zhǔn)確性和有效性，能夠滿足工程應(yīng)用需要。

參考文獻(xiàn)(References)：

[1]VENKATARAMANI K S,PLYBON R C,HOLM R G,et al.Aircraft engine icing model[R].AIAA 2008-440.

[2]楊倩,常士楠,袁修干.水滴撞擊特性的數(shù)值計(jì)算方法研究[J].航空學(xué)報(bào),2002,23(2):173-176.

[3]楊倩,常士楠,袁修干.發(fā)動(dòng)機(jī)進(jìn)氣道水滴撞擊特性分析[J].北京航空航天大學(xué)學(xué)報(bào),2002,28(3):362-365.

[4]韓鳳華,張朝民,王躍欣.飛機(jī)機(jī)翼表面水滴撞擊特性計(jì)算[J].北京航空航天大學(xué)學(xué)報(bào),1995，21(3)：16-21.

[5]張強(qiáng),胡利,曹義華.過(guò)冷水滴撞擊三維機(jī)翼的數(shù)值模擬[J].航空動(dòng)力學(xué)報(bào),2009,24(6):1345-1350.

[6]楊勝華,林貴平,申曉斌.三維復(fù)雜表面水滴撞擊特性計(jì)算[J].航空動(dòng)力學(xué)報(bào),2010,25(2):284-290.

[7]申曉斌,林貴平,楊勝華.三維發(fā)動(dòng)機(jī)進(jìn)氣道水滴撞擊特性分析[J].北京航空航天大學(xué)學(xué)報(bào),2011,37(1):1-5.

[8]吳孟龍,常士楠,冷夢(mèng)堯,等.基于歐拉法模擬旋轉(zhuǎn)帽罩水滴撞擊特性[J].北京航空航天大學(xué)學(xué)報(bào),2014,40(9):1263-1267.

[9]胡劍平,劉振俠,張麗芬.發(fā)動(dòng)機(jī)整流支板大尺寸過(guò)冷水滴撞擊特性[J].航空學(xué)報(bào),2011,32(10):1778-1785.

[10]CROWE C T.Review-Numerical models for dilute gas-particle[J].Journal of Fluids Engineering,1982(104):297-303.

[11]TONG X L,LUKE E A.Eulerian simulations of icing collection efficiency using a singularity diffusion model[R].AIAA 2005-1246.

[12]PAPAKAKIS M,WONG S,RACHMAN A A,et al.Large and small droplet impingement data on airfoils and two simulated ice shapes[R].NASA/TM-2007-213959.

[13]BOURGAULT Y,BOUTANIOS Z,HABASHI W G.Three-dimensional eulerian approach to droplet impingement simulation using fensap-ice,part1:model,algorithm,and validation[J].Journal of Aircraft,2000,37(1):95-103.

[14]WIROGO S,SRIRAMBHATLA S.An eulerian method to calculate the collection efficiency on two and three dimensional bodied[R].AIAA 2003-1073.

[15]BIDWELL C S,MOHLER S R.Collection efficiency and ice accretion calculations for a sphere,a swept MS(1)-317 wing,a swept NACA-0012 wing tip,and axisymmetric inlet,and a Boeing 737-300[R].AIAA-95-0755.

(責(zé)任編輯：宋麗萍英文審校：隋華)

Calculation of water droplet impingement characteristics based on UDS frame

MAO Xiao-dong1，BU Xue-qin2，Zhao Guo-chang1，Wang Jian-ming1

(1.Faculty of Aerospace Engineering,Shenyang Aerospace University,Shenyang 110136，China;2.School of Aeronautic Science and Engineering,Beijing University of Aeronautics and Astronautics,Beijing 100191,China)

Abstract：In order to predict the water droplet impingement characteristics of 2-dimensional and 3-dimensional surfaces under compressible and uncompressible flow conditions,a calculation method based on UDS frame was approached.The air-droplet two-phase flow was considered as one-way coupled.The airflow field could be solved previously and independently using N-S equations.The droplet control equations were user-programed and iterated under fluent UDS frame.To remove the droplet from wall boundaries,a treatment of boundary conditions on the impingement surfaces was given.A self-adaptation artificial numerical diffusion model was proposed to enhance the iteration stability and avoid the singularity of droplet volume fraction in localized region.Validation cases were performed for typical 2-dimensional cylinder and 3-dimensional sphere,through the comparison with the published experimental results,the approached calculation method was proved to be accurate and effective in predicting the droplet impingement characteristics on complex surfaces under various flow conditions.

Key words：UDS frame;water droplet impingement;icing;two-phase flow;numerical calculation

doi:10.3969/j.issn.2095-1248.2016.01.002

中圖分類(lèi)號(hào)：V211.3

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：2095-1248(2016)01-0008-05

作者簡(jiǎn)介：毛曉東(1984-)，男，遼寧大連人，副教授，主要研究方向：飛機(jī)與發(fā)動(dòng)機(jī)防除冰，E-mail：mxdbh@163.com。

基金項(xiàng)目：國(guó)家自然科學(xué)基金(項(xiàng)目編號(hào)：51206008)

收稿日期：2015-09-17