馬乙楗，柴得林，王強(qiáng)，易賢，*，孔滿昭

1．中國空氣動(dòng)力研究與發(fā)展中心 結(jié)冰與防除冰重點(diǎn)實(shí)驗(yàn)室，綿陽 621000

2．中國空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，綿陽 621000

3．航空工業(yè)第一飛機(jī)設(shè)計(jì)研究院，西安 710089

航空發(fā)動(dòng)機(jī)攝入冰晶易發(fā)生結(jié)冰，進(jìn)而導(dǎo)致發(fā)動(dòng)機(jī)動(dòng)力損失或損壞［1-2］。相比于傳統(tǒng)過冷水結(jié)冰，冰晶結(jié)冰存在冰晶運(yùn)動(dòng)、相變、黏附等物理現(xiàn)象，結(jié)冰機(jī)理復(fù)雜［3-5］。由于研究冰晶結(jié)冰的實(shí)驗(yàn)設(shè)備復(fù)雜、成本昂貴［6］，有必要開展冰晶結(jié)冰數(shù)值計(jì)算方法研究［7］。在冰晶結(jié)冰計(jì)算中，冰晶的運(yùn)動(dòng)相變決定了撞擊壁面的冰晶的形態(tài)與撞擊位置，黏附特性決定了參與結(jié)冰相變的冰晶的質(zhì)量通量，二者直接關(guān)系到壁面冰形的增長與發(fā)展，是準(zhǔn)確模擬冰晶結(jié)冰現(xiàn)象的關(guān)鍵環(huán)節(jié)［8］。因此，為深刻認(rèn)識冰晶結(jié)冰現(xiàn)象，需要開展冰晶運(yùn)動(dòng)相變與黏附計(jì)算方法研究，系統(tǒng)分析冰晶運(yùn)動(dòng)相變與黏附特性。

圍繞冰晶運(yùn)動(dòng)相變以及冰晶黏附過程，國內(nèi)外研究者開展了一系列研究，研究表明冰晶形狀對冰晶運(yùn)動(dòng)以及相變過程有顯著影響［9-11］，Norde等［12］模擬了冰晶運(yùn)動(dòng)過程中的傳熱過程，指出粒子運(yùn)動(dòng)和傳熱強(qiáng)烈依賴于冰晶形狀，在冰晶大量融化之前冰晶形狀通常為非球狀；Nilamdeen和Habashz［13］假設(shè)冰晶為平板狀，開展了冰晶的動(dòng)力學(xué)特性研究，指出冰晶在流場的流動(dòng)特性接近于圓盤狀顆粒；Iuliano等［14］假設(shè)冰晶為圓柱或圓盤狀，開展了冰晶阻力模型和傳熱過程的數(shù)值模擬；Trontin等［10］提出了冰晶的傳熱傳質(zhì)模型，并對自由流中的非球狀冰晶顆粒傳熱傳質(zhì)進(jìn)行了研究，通過實(shí)驗(yàn)所得冰形完善了冰晶黏附模型。郭向東等［11］假設(shè)冰晶為六角板狀和六棱柱狀，采用Villedieu等［15］提出的熱力學(xué)模型研究了冰晶形狀對其運(yùn)動(dòng)過程中換熱和軌跡的影響。針對冰晶黏附效率，Zhang等［16］通過比較冰晶粒徑和結(jié)構(gòu)表面水膜厚度來判斷冰晶是否黏附或破碎，該過程在Fluent軟件上進(jìn)行，計(jì)算得到黏附冰晶質(zhì)量通量，進(jìn)而計(jì)算出結(jié)冰冰形。Villedieu等［15］提出了基于結(jié)構(gòu)表面水膜厚度、粒子表面水膜厚度、粒子等效粒徑和法向動(dòng)能恢復(fù)系數(shù)的冰晶黏附效率計(jì)算方法，該方法計(jì)算冰形與Trontin等［10］實(shí)驗(yàn)結(jié)果的吻合效果并不理想。Trontin等［10］提出了一種基于冰晶粒子融化率、總液態(tài)水含量和總水含量的冰晶黏附效率的計(jì)算方法，并根據(jù)NASA-NRC混合氣象條件結(jié)冰的實(shí)驗(yàn)冰形數(shù)據(jù)校核了該黏附效率計(jì)算模型。綜上，國內(nèi)外學(xué)者已經(jīng)圍繞冰晶運(yùn)動(dòng)相變過程以及冰晶黏附過程建立了一系列數(shù)值計(jì)算模型，但在相關(guān)研究中，僅根據(jù)少數(shù)結(jié)冰實(shí)驗(yàn)、冰形數(shù)據(jù)確定了模型系數(shù)并做了簡單模型驗(yàn)證，未開展冰晶形狀對冰晶運(yùn)動(dòng)相變特性影響的深入研究，也未開展來流參數(shù)和冰晶融化效應(yīng)對冰晶黏附特性影響的系統(tǒng)性分析與研究。

因此，本文建立了拉格朗日框架下冰晶運(yùn)動(dòng)-傳熱傳質(zhì)耦合的數(shù)值計(jì)算方法，分析基于Monte Carlo方法的冰晶撞擊與黏附收集系數(shù)計(jì)算方法；依托NNWICE平臺［17］開發(fā)相應(yīng)的計(jì)算程序，實(shí)現(xiàn)冰晶運(yùn)動(dòng)-傳熱傳質(zhì)耦合與冰晶黏附的數(shù)值計(jì)算。以NACA0012翼型為對象，完成翼面冰晶/混合相結(jié)冰過程中，冰晶運(yùn)動(dòng)相變與黏附過程的模擬，研究冰晶形狀對冰晶運(yùn)動(dòng)相變特性的影響，系統(tǒng)分析不同來流溫度和液態(tài)水含量條件下冰晶的黏附特性，研究冰晶運(yùn)動(dòng)相變與黏附特性的影響規(guī)律。

1 計(jì)算方法

由于撞擊到翼型表面的冰晶不一定全部黏附，黏附冰晶量需在撞擊收集系數(shù)基礎(chǔ)上進(jìn)一步計(jì)算冰晶黏附特性得到。基于上述分析，計(jì)算冰晶黏附特性的流程如圖1所示。首先，求解空氣流場，為冰晶運(yùn)動(dòng)相變、水滴運(yùn)動(dòng)計(jì)算提供背景流場；再通過冰晶運(yùn)動(dòng)-傳熱傳質(zhì)計(jì)算得到冰晶撞擊特性與撞擊表面冰晶物態(tài)；最后，根據(jù)計(jì)算出的冰晶、水滴撞擊特性和冰晶物態(tài)完成冰晶黏附特性計(jì)算。采用NNW-FlowStar軟件完成空氣流場計(jì)算，空間離散采用一階迎風(fēng)格式，通量格式采用Roe格式，湍流模型采用SST（Shear Stress Transport）k-ω模型。

圖1 冰晶黏附特性計(jì)算流程Fig. 1 Calculation process of ice crystal adhesion characteristics

1.1 冰晶運(yùn)動(dòng)-相變計(jì)算

1） 冰晶運(yùn)動(dòng)控制方程

根據(jù)Stokes定理、牛頓第二定律，在拉格朗日框架下的冰晶運(yùn)動(dòng)控制方程為

式中：up為冰晶運(yùn)動(dòng)速度；ua為空氣速度；CD為阻力系數(shù)；Rep為冰晶相對雷諾數(shù)；μa為空氣動(dòng)力黏度；t為時(shí)間；ρa(bǔ)和ρp分別為空氣和冰晶粒子的密度；g為重力加速度；dp為冰晶等效體積直徑；系數(shù)

2） 冰晶形狀參數(shù)與阻力系數(shù)

采用以下形狀參數(shù)表征不同形狀冰晶［16］。

冰晶等效體積直徑（dp）表示具有和冰晶粒子相同體積的球狀冰晶的直徑，即

式中：Vp為冰晶粒子的體積。

冰晶球度（Φ）表示等效球狀冰晶表面積（Ap）與冰晶實(shí)際表面積（A）之比，即

冰晶橫向球度（Φ⊥）表示等效球狀冰晶的橫截面積（Ap，⊥）與冰晶流向投影的橫截面積（A⊥）之比，即

冰晶顆粒縱橫比E（Aspect Ratio）表示冰晶軸向長度b（沿旋轉(zhuǎn)軸）與徑向長度a（垂直于旋轉(zhuǎn)軸）之比，即

通常冰晶形狀以板狀和柱狀為主［11］。為模擬對流云內(nèi)的冰晶顆粒，表1給出了球狀和3種非球狀冰晶粒子的相關(guān)形狀參數(shù)計(jì)算公式［15］。

表1 冰晶形狀參數(shù)計(jì)算公式［15］Table 1 Calculation formulas for ice crystal shape parameters［15］

H?lzer和Sommerfeld［18］提出了一種考慮粒子方向并適用于所有形狀粒子的阻力系數(shù)表達(dá)式，本文采用式（7）計(jì)算式（1）中的冰晶阻力系數(shù)。

3） 冰晶相變模型

根據(jù)Trontin等［10］提出的相變模型，冰晶運(yùn)動(dòng)中的相變過程分為3個(gè)階段：① 過冷冰晶階段，冰晶溫度低于融化溫度，吸熱后溫度逐漸升高，僅發(fā)生升華；② 冰晶融化階段，冰晶溫度保持為融化溫度，吸熱后發(fā)生融化與蒸發(fā)，直至完全融化；③ 水滴蒸發(fā)階段，冰晶完全融化為水滴，溫度高于融化溫度，吸熱后溫度不斷升高，僅發(fā)生蒸發(fā)。

圖2展示了各階段冰晶狀態(tài)，其中Tp表示冰晶粒子溫度，Tm表示冰晶融化溫度。各階段傳熱傳質(zhì)計(jì)算公式見文獻(xiàn)［10］。

4） 冰晶運(yùn)動(dòng)相變計(jì)算方法

冰晶運(yùn)動(dòng)方程與傳熱傳質(zhì)方程采用一階歐拉向前時(shí)間離散格式進(jìn)行求解，對于冰晶運(yùn)動(dòng)控制方程式（1），其離散形式為

圖2 各階段冰晶狀態(tài)示意圖Fig. 2 Schematic diagram of state of ice crystals at each stage

式中：n代表第n個(gè)時(shí)間步；Δt代表時(shí)間步長。而冰晶傳熱傳質(zhì)方程根據(jù)冰晶所處狀態(tài)不同，具體公式見文獻(xiàn)［10］。

通常冰晶運(yùn)動(dòng)伴隨傳熱傳質(zhì)發(fā)生，本文按以下計(jì)算方法求解冰晶運(yùn)動(dòng)-相變耦合模型：① 初始化冰晶參數(shù)（冰晶初始位置、粒徑、形狀、初始溫度）；② 采用徑向基函數(shù)插值方法由冰晶空氣流場變量插值得到冰晶所在位置處物理量；③ 求解冰晶相變方程更新冰晶質(zhì)量、形態(tài)、球度、溫度等參數(shù)與變量；④ 將與求解運(yùn)動(dòng)方程相關(guān)的參數(shù)變量代入運(yùn)動(dòng)方程，計(jì)算更新冰晶位置、速度等變量；⑤ 重復(fù)步驟2~步驟4直至冰晶撞擊壁面/飛出指定區(qū)域或完全蒸發(fā)，完成軌跡求解。

1.2 冰晶黏附特性與收集系數(shù)計(jì)算

1） 冰晶撞擊收集系數(shù)計(jì)算

采用基于統(tǒng)計(jì)原理的Monte Carlo方法［19-20］計(jì)算冰晶撞擊收集系數(shù)，該方法既能適應(yīng)復(fù)雜三維外形，又能克服軌跡交叉、冰晶運(yùn)動(dòng)過程中質(zhì)量損失等問題。

在Monte Carlo方法下，冰晶撞擊收集系數(shù)（βs，MC）的計(jì)算式為

式 中：Ms為 網(wǎng) 格 單 元s從τ0到τtotal時(shí) 刻 的 質(zhì) 量 流率積分；nτ為τ時(shí)刻撞擊到s上的冰晶數(shù)量；mˉp為單個(gè)粒子的質(zhì)量；T為總時(shí)間步數(shù)；As為表面單元s的面積；u∞表示來流速度；IWC表示來流冰水含量。

2） 冰晶融化率與黏附效應(yīng)模型

定義mp，i為冰核質(zhì)量即冰晶中冰的質(zhì)量，則融化率ηm可定義為

式中：m表示冰晶總質(zhì)量。

文獻(xiàn)［21］提出了將黏附效率系數(shù)εs作為冰晶黏附效應(yīng)的數(shù)學(xué)描述，其定義為黏附的冰晶質(zhì)量與撞擊冰晶質(zhì)量之比，并給出了計(jì)算冰晶黏附效率的一種計(jì)算方式。首先定義當(dāng)空氣中只有冰晶以及融化冰晶時(shí)的黏附效率系數(shù)εs，c，其計(jì)算式為

式中：Kc為比例常數(shù)，通常取Kc=2.5。

定義空氣中同時(shí)含有冰晶和水滴時(shí)的黏附效率系數(shù)εs，d，計(jì)算式為

式中：Kd為比例常數(shù)，其計(jì)算式為

式中：Kd*=0.6；Tw為壁面處溫度；Tf為冰晶融化溫度，根據(jù)水的三相圖可知在來流入口壓力為98 000 Pa時(shí)，Tf=273.152 75 K；T*為冰晶完全無法黏附的閾值溫度，取文獻(xiàn)［21］中的建議值273.15 K。

φc表示撞擊冰晶的質(zhì)量分?jǐn)?shù)，φd表示撞擊水滴的質(zhì)量分?jǐn)?shù)，其計(jì)算式為

式中：βimp，c、βd分別表示冰晶、水滴撞擊局部收集系數(shù)； LWC表示來流液態(tài)水含量。

對式（12）~式（15）進(jìn)行推導(dǎo)，可得

根據(jù)εs，c和εs，d，給出冰晶黏附效率εs和冰晶黏附收集系數(shù)βcatch的計(jì)算公式為

此外，式（13）表明壁面溫度也會對冰晶黏附率產(chǎn)生影響。為探究融化率與來流液態(tài)水含量對冰晶黏附收集系數(shù)的影響，本文假設(shè)翼面溫度為恒定值Tw=273.5 K，收集過程均在Tf

2 冰晶運(yùn)動(dòng)相變特性

2.1 冰晶形狀對運(yùn)動(dòng)軌跡的影響

通過球度對冰晶形狀特征進(jìn)行數(shù)學(xué)描述，本節(jié)將討論不同來流條件下，冰晶球度對冰晶運(yùn)動(dòng)軌跡的影響。

以二維NACA0012翼型作為結(jié)冰對象，翼型弦長為1 m，參照文獻(xiàn)［13］，按照表2設(shè)置來流參數(shù)并開展計(jì)算與分析，表中MVD表示冰晶粒子的平均容積直徑。相同最大尺度（軸向長度與徑向長度的較大者）的不同形狀冰晶在運(yùn)動(dòng)過程中的沿程球度與運(yùn)動(dòng)軌跡結(jié)果如圖3~圖6所示。

當(dāng)來流總溫為?5 ℃時(shí)，由圖3可看出，冰晶球度在運(yùn)動(dòng)過程中保持恒定。這是由于冰晶未發(fā)生融化。圖4展示了同一位置釋放的4種形狀冰晶的運(yùn)動(dòng)軌跡，3種非球狀冰晶運(yùn)動(dòng)軌跡較為接近且與球狀差距較大。這是由于球狀冰晶的質(zhì)量與球度都明顯大于非球狀冰晶。圖5給出了來流總溫35 ℃冰晶沿程球度的變化。由圖5可知，非球狀冰晶融化前球度不變，開始融化后球度沿程增大，長橢球狀冰晶球度將逐漸增大至1（融化為水滴）；六角平板狀以及扁橢球狀冰晶球度接近且變化趨勢基本一致。圖6展示了來流總溫35 ℃下同一位置釋放4種冰晶的運(yùn)動(dòng)軌跡，與圖4相似，非球狀冰晶的運(yùn)動(dòng)軌跡接近且與球狀冰晶的差距明顯。算例1、算例2的數(shù)值實(shí)驗(yàn)結(jié)果證明，冰晶形狀將對冰晶運(yùn)動(dòng)軌跡產(chǎn)生顯著影響。

表2 工況參數(shù)Table 2 Running parameters

圖3 來流總溫?5 ℃下冰晶沿程球度Fig. 3 Sphericity of ice crystals along trajectory at total temperature of ?5 ℃

圖4 來流總溫?5 ℃下冰晶運(yùn)動(dòng)軌跡Fig. 4 Trajectories of ice crystals at total temperature of ? 5 ℃

圖5 來流總溫35 ℃下冰晶沿程球度Fig. 5 Sphericity of ice crystals along trajectory at total temperature of 35 ℃

圖6 來流總溫35 ℃下冰晶運(yùn)動(dòng)軌跡Fig. 6 Trajectories of ice crystals at total temperature of 35 ℃

2.2 冰晶形狀對冰晶融化率的影響

采用表2所示的工況參數(shù)，僅改變來流總溫為15 ℃、35 ℃設(shè)置算例3和算例4。計(jì)算得到相同最大尺度的不同形狀冰晶的沿程溫度及融化率變化曲線如圖7、圖8所示。

圖7 來流總溫15 ℃下冰晶粒子溫度及融化率變化曲線Fig. 7 Variation of particle temperature and melting ratio at total temperature of 15 ℃

圖8 來流總溫35 ℃下冰晶粒子溫度及融化率變化曲線Fig. 8 Variation of particle temperature and melting ratio at total temperature of 35 ℃

圖7中不同形狀的冰晶運(yùn)動(dòng)過程中溫度與融化率變化曲線差異明顯，長橢球狀冰晶溫升與融化最快，其次是扁橢球狀和六角平板狀冰晶，球狀最慢；所有冰晶撞擊到壁面時(shí)均未完全融化，而球狀冰晶剛開始融化。這是由于球狀冰晶的質(zhì)量明顯大于其余冰晶，升溫所需外界能量更多。由于扁橢球狀與六角平板狀冰晶形狀和球度類似，故其溫升與融化速率差異較小。

圖8給出了來流總溫35 ℃時(shí)4種冰晶粒子溫度及冰晶融化率變化曲線，冰晶溫度與融化率變化趨勢與算例3基本相同。由于來流總溫更高，冰晶粒子換熱更劇烈，升溫與融化速率顯著增加；相同冰晶撞擊到壁面時(shí)，算例4中融化率較算例3增大；長橢球狀冰晶已完全融化。

算例3和算例4的數(shù)值實(shí)驗(yàn)結(jié)果證明，冰晶形狀對運(yùn)動(dòng)過程中冰晶溫度與融化率具有顯著影響，這將影響冰晶黏附特性及后續(xù)冰晶結(jié)冰相變過程中撞擊到結(jié)構(gòu)表面的冰晶粒子的顯熱、物態(tài)進(jìn)而影響冰形計(jì)算；此外，在最大尺度相同的條件下，相較于扁橢球狀和六角平板狀冰晶，長橢球狀冰晶在運(yùn)動(dòng)過程中換熱、融化現(xiàn)象更顯著。

3 冰晶黏附特性

云霧場通常是指包含由液態(tài)水滴和固態(tài)冰晶形成的云霧的空間區(qū)域，云霧場性能主要由以下云霧參數(shù)衡量：粒子的平均容積直徑（MVD）、液態(tài)水含量（LWC）、冰水含量（IWC）、總水含量（TWC）等，上述參數(shù)是影響翼面結(jié)冰的重要因素。根據(jù)Currie等［22］的實(shí)驗(yàn)結(jié)果，冰晶黏附效率取決于撞擊粒子的液態(tài)水含量與總水含量的比值。通常LWC來自于2部分：冰晶融化與混合相來流中的過冷水滴。為分析冰晶融化以及過冷水滴含量對冰晶收集系數(shù)的影響，將通過改變來流溫度以及云霧參數(shù)開展數(shù)值計(jì)算研究。

3.1 來流溫度對黏附特性的影響

采用150 μm六角平板狀冰晶，參照表2設(shè)置其余工況參數(shù)，并改變來流總溫分別為?5、5、15、35 ℃，對算例5~算例8進(jìn)行數(shù)值實(shí)驗(yàn)。

圖9給出了不同溫度條件下的冰晶撞擊收集系數(shù)和冰晶黏附效率（εs）。圖中S表示翼型上下翼面的點(diǎn)到翼型前緣駐點(diǎn)的距離，L表示特征長度，此處為翼型的弦長。從圖9（a）可看出來流總溫為?5 ℃和5 ℃下的冰晶撞擊收集系數(shù)（βimp，c）曲線接近，隨著來流總溫升高冰晶的撞擊極限和撞擊收集系數(shù)峰值均增大。?5 ℃與5 ℃的冰晶均未進(jìn)入融化階段，質(zhì)量損失均很小，撞擊收集系數(shù)曲線相似。隨著來流總溫升高至冰晶融化，冰晶在靠近翼型表面繞流時(shí)隨流性變?nèi)?，冰晶撞擊收集系?shù)峰值、撞擊極限增大。

圖9 冰晶撞擊收集系數(shù)與黏附效率曲線Fig. 9 Ice crystal impingement coefficient and sticking efficiency

由圖9（b）可知算例5與算例6的黏附效率曲線接近且呈鐘形。算例7與算例8冰晶黏附效率的范圍以及峰值較算例5與算例6顯著增大且維持在某一定值附近。對于算例5與算例6，由于冰晶并未融化，黏附效率完全取決于水滴收集量，曲線呈鐘形。隨著來流總溫的升高冰晶開始融化，這時(shí)黏附效率大小主要取決于融化率（ηm）和水滴收集系數(shù)（βd），其值取水滴對應(yīng)黏附效率（εs，d）與 冰 晶 融 化 對 應(yīng) 黏 附 效 率（εs，c）中 較 大 值。故算例7與算例8中冰晶黏附效率范圍與冰晶撞擊收集系數(shù)范圍一致。此外，算例7中靠近前緣附近水滴導(dǎo)致的黏附效應(yīng)相較于冰晶融化導(dǎo)致的黏附效應(yīng)更明顯，該區(qū)域的黏附效率明顯增高。

圖10給出了?5 ℃、15 ℃和35 ℃下冰晶、水滴收集系數(shù)和冰晶黏附收集系數(shù)（βcatch）曲線，圖10（d）給出了3種來流總溫條件下冰晶黏附收集系數(shù)曲線。選取收集系數(shù)曲線主要特征量——收集系數(shù)峰值和撞擊上極限弧面坐標(biāo)值，并進(jìn)行對比，表3給出了相應(yīng)的數(shù)據(jù)，表中Tt∞為來流總溫，峰值降幅指βcatch，峰值相對于βimp，c峰值的降幅，上極限降幅指黏附上極限弧面坐標(biāo)相對于撞擊上極限弧面坐標(biāo)的降幅。

圖10 混合相條件下冰晶和水滴收集系數(shù)曲線Fig. 10 Ice crystal and droplet collection coefficient under mixed phase conditions

表3 不同來流總溫收集系數(shù)曲線主要特征對比Table 3 Comparison of main characteristics of collection coefficient at different total temperatures

如圖10（a）所示，在來流總溫?5 ℃下冰晶不融化，冰晶僅在液滴撞擊區(qū)域內(nèi)發(fā)生黏附，黏附收集系數(shù)極限與水滴的撞擊極限相同。從圖10（b）可看出來流總溫15 ℃冰晶融化，冰晶撞擊收集系數(shù)峰值與黏附效率均增大，故其黏附收集系數(shù)的峰值增大；冰晶黏附收集系數(shù)極限與冰晶撞擊極限一致。如圖10（c）所示，來流總溫為35 ℃時(shí)，冰晶進(jìn)一步融化，此時(shí)融化導(dǎo)致的黏附效應(yīng)占主導(dǎo)，黏附效率值趨于1，故冰晶撞擊收集系數(shù)曲線與冰晶黏附收集系數(shù)曲線接近。由圖10（d）可知，冰晶黏附收集系數(shù)峰值與極限隨著來流總溫的升高呈增大趨勢。分析表明冰晶融化對冰晶黏附有顯著影響。

3.2 來流IWC和LWC對黏附特性的影響

針對混合相條件下，來流云霧參數(shù)對冰晶黏附效率和黏附收集系數(shù)產(chǎn)生的影響展開研究與討論。

由于不同來流總溫條件下，冰晶融化程度存在差異，這對冰晶黏附收集系數(shù)結(jié)果產(chǎn)生影響。為去除該影響因素，參照表2中算例1設(shè)置工況參數(shù)，來流總溫取?5 ℃保持不變，但改變來流LWC和IWC值 設(shè) 置 算 例9：LWC = 0.3 g/m3，IWC = 0.3 g/m3；算 例10：LWC = 0.3 g/m3，IWC = 0.7 g/m3；算 例11：LWC = 0.7 g/m3，IWC = 0.3 g/m3；算 例12：LWC = 0.7 g/m3，IWC = 0.7 g/m3；并進(jìn)行數(shù)值實(shí)驗(yàn)。

4種工況下計(jì)算結(jié)果中冰晶撞擊上極限弧面坐標(biāo)值均為0.052 68。冰晶黏附上極限弧面坐標(biāo)值均為0.037 2，這是由于冰晶均未融化，冰晶黏附極限均與水滴撞擊極限一致。表4給出撞擊與黏附收集系數(shù)峰值的對比結(jié)果。

表4 不同來流云霧條件下收集系數(shù)曲線主要特征對比（算例9~算例12）Table 4 Comparison of main characteristics of collection coefficient under different conditions of cloud field （Case9-Case12）

?5 ℃來流總溫下冰晶黏附效率和冰晶黏附收集系數(shù)計(jì)算結(jié)果如圖11所示。圖11顯示，隨著LWC/TWC比值增大，冰晶黏附效率和黏附收集系數(shù)峰值逐漸增大。當(dāng)LWC/TWC值一定時(shí)（算例9與算例12），在來流總溫一致的情況下，冰晶將有著同樣的黏附效率和黏附收集系數(shù)。這是因?yàn)楸Р蝗诨瘯r(shí)，黏附效率取決于撞擊壁面處水滴質(zhì)量分?jǐn)?shù)大小。

為進(jìn)一步研究冰晶融化時(shí)LWC和IWC對冰晶黏附效率與冰晶黏附收集系數(shù)的影響，來流其余參數(shù)設(shè)置參照算例9~算例12，改變來流總溫為15 ℃，分別對算例13~算例16進(jìn)行數(shù)值實(shí)驗(yàn)。

算例13~算例16計(jì)算結(jié)果中冰晶撞擊上極限與冰晶黏附上極限的弧面坐標(biāo)值均為0.073 93，這是由于4種工況下冰晶均融化，冰晶黏附極限均與冰晶撞擊極限一致。表5給出了收集系數(shù)曲線特征對比。

圖12給出了15 ℃條件下不同IWC與LWC下冰晶黏附效率和黏附收集系數(shù)的計(jì)算結(jié)果。如圖12（a）所示，除算例15外，其余算例的黏附效率曲線差距很小，整體保持在0.4左右，算例15中靠近翼型前緣位置冰晶黏附效率相較于其他算例明顯增高，由圖12（b）可得相似結(jié)果；此外，圖中各算例的冰晶黏附收集系數(shù)曲線極限均與冰晶撞擊收集系數(shù)極限一致。

圖11 來流總溫?5 ℃不同云霧參數(shù)下黏附效率與黏附收集系數(shù)Fig. 11 Sticking efficiency and catching coefficient with different cloud field parameters of total temperature of incoming flow ?5 ℃

表5 不同來流云霧條件下收集系數(shù)曲線主要特征對比（算例13~算例16）Table 5 Comparison of main characteristics of collection coefficient under different conditions of cloud field （Case13-Case16）

圖12 來流總溫15 ℃不同云霧參數(shù)下黏附效率與黏附收集系數(shù)Fig. 12 Sticking efficiency and catching coefficient with different cloud field parameters of total temperature of incoming flow 15 ℃

由于來流總溫為15 ℃時(shí)冰晶融化，冰晶黏附收集系數(shù)極限范圍與冰晶撞擊收集系數(shù)一致。對于算例15外的算例，撞擊水滴的質(zhì)量分?jǐn)?shù)較小，冰晶融化導(dǎo)致的黏附效應(yīng)占主導(dǎo)，故具有相同的黏附效率曲線和黏附收集系數(shù)曲線；而算例15的來流LWC/TWC值較大，在翼型前緣位置附近收集到的液態(tài)水滴質(zhì)量分?jǐn)?shù)較高，水滴對應(yīng) 黏 附 效 率（εs，d）大 于 冰 晶 融 化 對 應(yīng) 黏 附 效 率（εs，c），故在該區(qū)域冰晶黏附效率和黏附收集系數(shù)相較于其余算例有明顯增高。

4 結(jié) 論

基于NNWICE平臺，對冰晶/混合相結(jié)冰過程中冰晶運(yùn)動(dòng)相變與黏附的計(jì)算方法與影響因素開展了研究，得到以下結(jié)論：

1） 建立了拉格朗日框架下冰晶運(yùn)動(dòng)-傳熱傳質(zhì)耦合的數(shù)值計(jì)算方法，分析了冰晶幾何形狀對冰晶運(yùn)動(dòng)相變特性的影響，發(fā)現(xiàn)同粒徑條件下類似柱狀的長橢球狀冰晶相較于扁橢球狀冰晶和六角平板狀冰晶在運(yùn)動(dòng)過程中換熱、融化現(xiàn)象更顯著。

2） 分析了基于Monte Carlo方法的冰晶/混合相結(jié)冰黏附收集系數(shù)計(jì)算方法。該方法可以在考慮冰晶運(yùn)動(dòng)過程中質(zhì)量損失的情況下，有效模擬冰晶收集系數(shù)。

3） 系統(tǒng)分析了來流溫度等云霧參數(shù)對冰晶黏附特性的影響。來流總溫升高，冰晶融化率增加，撞擊表面微元液態(tài)水質(zhì)量分?jǐn)?shù)增大，冰晶更易黏附；來流LWC/TWC增大，撞擊表面微元液態(tài)水質(zhì)量分?jǐn)?shù)增大，冰晶更易黏附。