基于Messinger擴展模型蒸發對結冰溢流特性的影響研究

摘 要：傳統Messinger擴展模型將蒸發簡化為壁面溫度的線性函數，以便推導冰層和水膜溫度微分方程的原函數。但結冰過程中蒸發與壁面溫度并非線性關系，若引入蒸發非線性計算式，冰層和水膜溫度微分方程將無法推導出原函數。而且有研究者在應用Messinger擴展模型時忽略了蒸發對結冰的影響，在某些狀態下可能造成較大誤差。針對以上問題，本文對Messinger擴展模型做出改進。將冰層和水膜溫度假設為帶未知系數的線性方程，建立未知系數的隱示表達式，采用牛頓迭代對未知系數進行求解，以數值迭代代替方程推導。基于NACA0012翼型，對忽略蒸發以及采用不同蒸發計算式對結冰外形的影響進行了分析。通過計算不同環境溫度下的結冰外形，并與FENSAP-ICE軟件和試驗冰形進行對比，驗證了該改進方法的有效性。

關鍵詞：Messinger模型； Messinger擴展模型； 飛機結冰； 結冰計算； 結冰軟件

中圖分類號：V211.3 文獻標識碼：A DOI：10.19452/j.issn1007-5453.2024.02.007

基金項目： 航空科學基金（2022Z006026005）

數值模擬是研究飛機結冰的重要手段。1953年，B.L. Messinger[1]提出飛機結冰條件下的結冰模型，該模型將飛機結冰過程假設為定常問題，采用代數方程組定義結冰控制方程，容易程序實現，可快速得到計算結果，至今仍然被普遍采用。雖然很多研究者對Messinger模型做出改進[2-3]，但是理論框架基本不變。

在Messinger模型基礎上，Y. Bourgault等[4-5]提出Shallow Water模型，假設結冰表面未凍結的水會形成一層薄水膜，水膜受外部氣流剪切力的作用而流動，建立水膜高度和溫度的偏微分方程組，主要應用于FENSAP-ICE軟件。在2000年前后，T. G. Myers等[6-7]聯立冰層和水膜熱傳導方程，基于線性假設將結冰速率方程簡化為一階常微分方程，并稱之為Messinger擴展模型。考慮結冰過程中的水膜流動特性，Myers將淺水膜流動方程[8-9]耦合進了Messinger擴展模型中[10]，并進一步發展至可預測任意三維表面淺水膜的流動與凍結[11-12]。后續很多研究者基于Myers模型開展了結冰數值模擬研究[13-14]。

國內結冰數值模擬研究早期以Messinger模型為主[15-16]，近些年對Myers模型開展了非常多的研究工作。侯碩等[17]將Myers模型由二維笛卡兒坐標系[12]推廣至正交曲線坐標系，并與Messinger模型計算結果進行了對比。曹廣州等[18-19]將水膜模型擴展至三維復雜表面，分析了水膜流動對結冰模擬的改善作用。雷夢龍等[20]對Myers模型結冰類型的判斷標準進行了修改，并將單步法、多步法計算結果與試驗冰形進行了對比。李浩然等[21]將Messinger模型預估的壁面溫度作為Myers模型的邊界條件，從而解決了Myers模型壁面溫度設置的經驗性問題。

無論是Messinger模型，還是Messinger擴展模型，抑或是引入淺水膜流動方程的Myers模型，都涉及對流換熱、蒸發等熱流項的計算，其中蒸發一般為結冰表面溫度的非線性函數。在使用Messinger模型時，一般采用LEWICE使用手冊提出的蒸發計算式[22]。而為推導冰層和水膜溫度的顯示表達式，Myers在提出Messinger擴展模型時，將非線性蒸發計算式進行線性化處理[10]。國外研究者基本沿用了相同的計算式[13-14]，國內研究者雖然在文獻中未列出蒸發的具體計算式[17，20-21]，但結冰速率方程與Messinger擴展模型是一致的[6-7]，可以推斷蒸發仍然采用線性計算式。

根據現有研究，未見有研究者在使用Messinger擴展模型時采用蒸發非線性計算式。而且參考文獻[18]和[19]忽略了蒸發，其出現液態水過度溢流的現象有可能是忽略蒸發導致的。因此，本文將蒸發非線性計算式引入Messinger擴展模型，以牛頓迭代代替微分方程原函數繁瑣的推導，以提高該模型對熱流任意計算式的適用性。另外，通過對比是否考慮蒸發計算得到的冰形，以驗證蒸發對結冰表面溢流水結冰特性的影響。計算NACA0012翼型在不同環境溫度下的結冰外形，與FENSAP-ICE軟件計算結果和試驗冰形進行對比，以驗證本文研究方法的有效性。

1 Messinger擴展模型

1.1 控制方程

上述方程組的初始條件和邊界條件如下所述，其中空氣/水界面和空氣/冰界面的邊界條件來自Messinger模型中的能量守恒方程。

（1） 在冰/壁界面

當過冷水滴撞擊到溫度在凍結溫度以下的壁面時，凍結過程經歷兩個階段[7]。基于恒定壁溫邊界條件，水滴撞擊到壁面時立即凍結。在該階段只有冰存在，即霜冰。隨著冰層逐漸增厚，當厚度足夠大時，冰層熱傳導速率降低，過冷水滴凍結潛熱將產生足夠的熱量使冰層表面溫度達到凍結溫度，并出現液態水。在該階段冰、水同時存在，即光冰。

1.2 霜冰

2 數值求解

2.1 求解步驟

由以上內容可知，Myers將結冰速率方程（4）轉化為一階常微分方程（20），該方程可通過Runge-Kutta多步迭代方法進行數值求解。下面詳細介紹具體的求解步驟，如圖1所示。該求解步驟只包含單步法，即整個結冰過程只使用一套網格、只計算一次流場和水滴撞擊特性。具體步驟為：（1）設置飛行和結冰條件參數，計算開始，讀取流場和水滴撞擊特性信息；（2）從駐點開始循環結冰表面控制體，每個控制體求解一次結冰模型；（3）根據式（21）和式（22）計算臨界結冰厚度bg和臨界結冰時間tg，根據分類方法判斷結冰生長階段；（4）當結冰處于霜冰階段時，根據式（15）計算霜冰生長厚度b，根據式（14）迭代求解系數a3，根據式（11）得到冰層溫度線性函數，代入冰層厚度b得到冰層表面溫度Ti；（5）當結冰處于光冰階段時，假設當前b的初值為極小值，代入式（17）得到系數a3，根據式（23）計算水膜厚度h，代入式（19），迭代求解系數a5，將a3和a5代入式（20），采用Runge-Kutta多步迭代法求解冰層厚度b；（6）根據式（5）計算流入水m？in；（7）遍歷所有結冰表面控制體；（8）計算結束，輸出計算結果。

2.2 方程系數迭代

根據Messinger擴展模型求解步驟，在求解結冰速率方程（20）之前，需要先求解系數a3和a5（見式（17）和式（19））。由式（19）可知，系數a5的計算式并不是顯示的，且其具體值取決于各熱流項的計算式，（見式（24）～式（32））。Myers將蒸發簡化為線性計算式（31），從而推導出系數a5的顯示計算式，以便于結冰速率方程（20）的求解。在各熱流項當中，蒸發在不同文獻中多出現不同的計算式，若針對每種蒸發計算式推導一種系數a5的顯示計算式，將使函數推導變得煩瑣。若蒸發采用式（29），系數a5的顯示計算式將無法導出。

3 計算結果與分析

空氣流場基于航空工業氣動院自主開發的航空飛行器氣動力仿真軟件飛廉（FeiLian）V1.0[25]進行計算，采用MENTER_SST湍流模型、Roe對流通量離散格式以及Venkatakrishnan限制器。水滴撞擊特性基于航空工業氣動院自主開發的航空飛行器結冰仿真軟件穿云（ChuanYun）V1.0進行計算，采用中心對流通量離散格式。

3.1 臨界結冰厚度分析

臨界結冰厚度bg是判斷霜冰生長轉為光冰生長的重要參數，對結冰生長速率有重要影響。這里忽略蒸發以及基于蒸發計算式（29）、式（31）分析蒸發對臨界結冰厚度的影響，計算狀態見表1，其中LWC為液態水含量。

3.2 蒸發影響分析

基于NACA0012翼型，弦長0.5334m，以表1狀態為基準，環境溫度分別為-4.46°C、-6.13°C、-10°C，忽略蒸發以及使用不同蒸發計算式得到的冰形與試驗冰形[26]的對比如圖7～圖9所示。圖中x/c、 y/c代表量綱一坐標，c代表弦長。

對比不同蒸發計算式得到的冰形，兩者整體上差異很小，原因在于蒸發量占凍結量的比重很小；兩者差異主要體現在結冰后半部分，與前面的分析結論一致，即當撞擊水質量較低、對流換熱系數較高時，蒸發對結冰的影響較大；在環境溫度為-4.46°C時，如圖7所示，與蒸發非線性計算式相比，蒸發線性計算式得到的冰形出現更靠后的溢流結冰。因此，使用簡化的蒸發線性計算式整體上對冰形的影響很小，但在高溫結冰條件下，會產生更多的溢流水結冰，對結冰極限位置的預測造成一定影響。

對比是否考慮蒸發得到的冰形，若忽略蒸發，將產生非常明顯的溢流結冰。在環境溫度為-6.13°C時，如圖8所示，考慮蒸發時上表面結冰極限位置與試驗冰形基本一致，而忽略蒸發時結冰極限位置明顯超出試驗冰形。造成這一現象的主要原因是，雖然蒸發量占凍結量的比重很小，但蒸發熱流占所有冷卻熱流的比重很大，忽略蒸發將降低結冰過程中的冷卻效應，導致凍結系數減小，未凍結的液態水量增大，液態水向下游溢流，超出水滴撞擊極限并凍結。

3.3 翼型結冰

基于NACA0012翼型，弦長0.5334m，以表1狀態為基準，環境溫度分別為-4.46°C、-6.13°C、-7.79°C、-10°C、-13.35°C、-23.35°C，分別采用Messinger模型和Messinger擴展模型進行計算，蒸發采用非線性計算式（29），并與FENSAP-ICE和試驗冰形進行對比，如圖10～圖15所示。

由于FENSAP-ICE采用的Shallow Water模型仍然基于平衡溫度假設，這一點與Messinger模型是相同的，因此穿云V1.0軟件Messinger模型的計算結果與FENSAP-ICE軟件的計算結果呈現出較好的一致性。在圖13和圖14中，FENSAP-ICE在翼型上表面水滴撞擊極限位置的溢流結冰更加明顯，這可能是由Shallow Water模型中采用剪切力驅動水膜流動導致的。

Messinger模型會明顯低估凍結系數[7，17]，導致更多的溢流水產生并在下游凍結，尤其在環境溫度為-4.46°C（見圖10）和-6.13°C（見圖11）條件下，溢流水結冰極限遠超水滴撞擊極限。而且在環境溫度為-23.35°C（見圖15）條件下，Messinger模型在駐點處凍結系數的低估導致駐點處撞擊水未完全凍結，向兩側溢出后立即凍結，導致非正常冰角出現。

與試驗冰形相比，在結冰趨勢上穿云V1.0采用Messinger擴展模型的計算結果要明顯優于Messinger模型和FENSAP-ICE軟件的計算結果。對于翼型下表面結冰極限，Messinger擴展模型計算結果與試驗冰形基本一致；對于翼型上表面結冰極限，除環境溫度為-4.46°C（見圖10）、-6.13°C（見圖11）和-7.79°C（見圖12）條件外，與試驗冰形基本一致。在平均結冰厚度上，Messinger擴展模型計算結果與試驗冰形有較好的一致性，在最大結冰厚度上仍然有一定差距。

由于Messinger擴展模型引入冰層熱傳導，避免Messinger模型對凍結系數的低估，對未凍結液態水的溢流產生抑制作用，因此Messinger擴展模型對結冰極限和平均結冰厚度的預測更加接近試驗冰形。而且在環境溫度為-23.35°C（見圖15）條件下，Messinger擴展模型的計算結果并未出現Messinger模型計算的非正常冰角，這也是由于冰層熱傳導的作用。對于高溫結冰條件，Messinger擴展模型在翼型上表面仍然產生過度的溢流水結冰現象。原因可能是此處結冰表面較為粗糙，現有粗糙度模型低估了此處的對流換熱系數，即使Messinger擴展模型引入冰層熱傳導，也無法抑制溢流水結冰。

4 結論

本文在原有Messinger擴展模型的基礎上，引入蒸發非線性計算式，以牛頓迭代冰層和水膜溫度線性方程未知系數的方法代替微分方程原函數的推導。通過計算與分析，得到如下結論：

（1）以牛頓迭代代替方程推導，使Messinger擴展模型能夠適用于任意形式的蒸發計算式（或其他熱流計算式），牛頓迭代可快速收斂，僅耗費較小計算資源。

（2）使用簡化的蒸發線性計算式，整體上對冰形的影響很小，但在高溫結冰條件下，會產生更多的溢流水結冰，對結冰極限位置的預測造成一定影響。若忽略蒸發，冷卻熱能降低，液態水過度溢流，將導致明顯的溢流水結冰現象。

（3）Messinger模型計算結果與FENSAP-ICE軟件計算結果基本一致，Messinger擴展模型對水滴凍結量的預測更加合理，能夠明顯抑制溢流水結冰，計算結果更加接近試驗冰形。對于高溫結冰條件，在未來工作中通過對下游結冰粗糙區域的對流換熱系數的精準模擬方法進行研究，將能夠進一度提高結冰極限預測精度。

參考文獻

[1]Messinger B L. Equilibrium temperature of an unheated icing surface as a function of airspeed[J]. Journal of the Aeronautical Sciences， 1953， 20（1）： 29-42.

[2] Montreuil E， Chazottes A， Guffond D. Enhancement of prediction capability in icing accretion and related performance penalties[R]. AIAA-2009-3969， 2009.

[3] 戰培國. 美國NASA結冰試驗設備體系綜述[J]. 航空科學技術， 2021，32（5）： 1-6. Zhan Peiguo. Review on the system of icing facilities in NASA[J]. Aeronautical Science Technology，2021，32（5）：1-6.（in Chinese）

[4]Bourgault Y， Habashi W G， Beaugendre， H. Development of a shallow water icing modeling FENSAP-ICE[R]. AIAA-1999-0246， 1999.

[5]Bourgault Y， Boutanios Z， Habashi W G. Three-dimensional Eulerian approach to droplet impingement simulation using FENSAP-ICE， Part 1： Model， algorithm and validation[J]. Journal of Aircraft， 2000， 37（1）： 95-103.

[6]Myers T G， Hammond D W. Ice and water film growth from incoming supercooled droplets[J]. International Journal of Heat and Mass Transfer， 1999， 42： 2233-2242.

[7]Myers T G. Extension to the Messinger model for aircraft icing[J]. AIAA Journal， 2001， 39（2）： 211-218.

[8]Myers T G. Thin films with high surface tension[J]. Society for Industrial and Applied Mathematics， 1998， 40（3）： 441-462.

[9]Myers T G， Thompson C P. Modelling the flow of water on aircraft in icing conditions[J]. AIAA Journal， 1998， 36（6）： 1010-1013.

[10]Myers T G， Charpin J P F， Thompson C P. Slowly accreting glaze ice due to supercooled water impacting on a cold surface[J]. Physics of Fluids， 2002， 14（1）： 240-256.

[11]Myers T G， Charpin J P F， Chapman S J. The flow and solidifi‐cation of a thin fluid film on an arbitrary three dimensional sur‐face[J]. Physics of Fluids， 2002， 14（8）： 2788-2803.

[12]Myers T G， Charpin J P F. A mathematical model for atmospheric ice accretion and water flow on a cold surface[J]. International Journal of Heat and Mass Transfer， 2004， 47： 5483-5500.

[13]？zgen S， Canibek M. Ice accretion simulation on multielement airfoils using extended Messinger model[J]. Heat and Mass Transfer， 2009， 45（3）： 305-322.

[14]？zgen S， Canibek M. In-flight ice formation simulation on finite wings and air intakes[J]. The Aeronautical Journal， 2012， 116： 337-362.

[15]朱程香，孫志國，付斌，等.水滴多尺寸分布對水滴撞擊特性和結冰增長的影響[J]. 南京航空航天大學學報，2010，42（5）： 620-624. Zhu Chengxiang， Sun Zhiguo， Fu Bin， et al. Effects of multidispersed droplet distribution on droplet impingement and ice accretion[J]. Journal of Nanjing University of Aeronautics Astronautics， 2010， 42（5）： 620-624. （in Chinese）

[16]孫志國，朱程香，付斌，等. 二維翼型結冰數值計算[J]. 航空動力學報，2010，25（7）： 1485-1490. Sun Zhiguo， Zhu Chengxiang， Fu Bin， et al. Numerical re‐search of ice accretion on two-dimensional airfoils[J]. Journal of Aerospace Power， 2010， 25（7）： 1485-1490. （in Chinese）

[17]侯碩，曹義華. 基于潤滑理論的二維積冰數值模擬[J]. 北京航空航天大學學報， 2014， 40（10）： 1442-1450. Hou Shuo， Cao Yihua. Numerical simulation of two dimensional ice accretion based on lubrication theory[J]. Journal of Beijing University of Aeronautics and Astronautics， 2014， 40（10）： 1442-1450. （in Chinese）

[18]曹廣州，吉洪湖，胡婭萍，等. 模擬飛機迎風面三維積冰的數學模型[J]. 航空動力學報，2011，26（9）： 1953-1963. Cao Guangzhou， Ji Honghu， Hu Yaping， et al. An icing model for simulating three dimensional ice accretion on the upwind surfaces of a plane[J]. Journal of Aerospace Power， 2011， 26（9）： 1953-1963. （in Chinese）

[19]曹廣州，吉洪湖，斯仁. 迎風面三維結冰過程中水膜流動的計算方法[J]. 航空動力學報，2015，30（3）： 677-685. Cao Guangzhou， Ji Honghu， Si Ren. Computational methodolo‐gy of water film flow in three-dimensional ice accretion on up‐wind surface[J]. Journal of Aerospace Power， 2015， 30（3）： 677-685. （in Chinese）

[20]雷夢龍，常士楠，楊波. 基于Myers模型的三維結冰數值仿真[J]. 航空學報，2018，39（9）： 121952. Lei Menglong， Chang Shinan， Yang Bo. Three-dimensional nu‐merical simulation of icing using Myers model[J]. Acta Aeronau‐tica et Astronautica Sinica， 2018， 39（9）： 121952. （in Chinese）

[21]李浩然，段玉宇，張宇飛，等. 結冰模擬軟件AERO-ICE中的關鍵數值方法[J]. 航空學報，2021，42（S1）： 726371. Li Haoran， Duan Yuyu， Zhang Yufei， et al. Numerical method of ice-accretion software AERO-ICE[J]. Acta Aeronautica et Astronautica Sinica， 2021， 42（S1）： 726371. （in Chinese）

[22]Wright W B. User’s manual for LEWICE version 3.2[R]. NASA-CR-2008-214255， 2008.

[23]Cansdale J T. Helicopter rotor ice accretion and protection research[J]. Vertica， 1981， 5（4）： 357-368.

[24]Verdin P， Charpin J P F， Thompson C P. Multistep results in ICECREMO2[J]. Journal of Aircraft， 2009， 46（5）： 1607.

[25]喬龍，李艷亮，楊思源，等. 基于面向對象的非結構航空CFD軟件體系結構設計[J]. 航空科學技術，2022，33（7）： 66-72. Qiao Long， Li Yanliang， Yang Siyuan， et al. Software architecture of an unstructured aviation CFD solver using object-oriented techniques[J]. Aeronautical Science Technology， 2022， 33（7）： 66-72. （in Chinese）

[26]Shin J， Bond T. Results of an icing test on a NACA 0012 airfoil in the NASA Lewis Icing Research Tunnel[R]. NASA-TM-105374， 1992.

Research on the Influence of Evaporation on Runback Icing Characteristics Based on Extended Messinger Model

Ning Yijun1， Xie Wenlong2， Zhu Dongyu1， Zhang Fukun1

1. Liaoning Provincial Key Laboratory of Aircraft Ice Protection， AVIC Aerodynamics Research Institute， Shenyang 110034， China

2. Nanjing University of Aeronautics and Astronautics， Nanjing 210016， China

Abstract： The traditional extended Messinger model simplifies evaporation as a linear function of the wall temperature in order to derive the original functions of the differential equations of the ice layer and water film temperature. However， the relationship between evaporation and wall temperature is not linear in the process of icing. If the nonlinear equation of evaporation is introduced， the original functions of the differential equations of ice layer and water film temperature cannot be derived. In addition， some researchers ignored the effect of evaporation on icing when applying extended Messinger model， which may cause large errors in some icing conditions. To solve the problems above， the extended Messinger model is improved. The ice layer and water film temperature are assumed to be a linear equation with unknown coefficients， the implicit expressions of unknown coefficients are established， the unknown coefficients are solved by Newton iteration， and the equation derivation is replaced by numerical iteration. Based on NACA0012 airfoil， the effects of ignoring evaporation and using different evaporation formulas on the icing shape were analyzed. The icing shapes at different ambient temperatures were calculated and compared with FENSAP-ICE software and experimental ice shapes， and the effectiveness of the improved method was verified.

Key Words： Messinger model； extended Messinger model； aircraft icing； icing calculation； icing software