周志宏，易 賢，桂業偉，李鳳蔚

（1．中國空氣動力研究與發展中心 空氣動力學國家重點實驗室，四川 綿陽 621000；2．西北工業大學航空學院，陜西 西安 710072）

水滴撞擊特性的高效計算方法

周志宏1，易 賢1，桂業偉1，李鳳蔚2

（1．中國空氣動力研究與發展中心 空氣動力學國家重點實驗室，四川 綿陽 621000；2．西北工業大學航空學院，陜西 西安 710072）

針對拉格朗日方法計算水滴撞擊特性效率低、通用性差等問題，發展了一種水滴撞擊特性的高效計算方法。在求解繞流流場的基礎上，結合逐級結構化管理的邊界信息存儲方式，采用目標擴散追蹤方法對水滴所在網格單元進行快速計算，并插值得到該點處的流場信息，逐個求解水滴運動方程得到各水滴的運動軌跡，從而確定水滴撞擊極限、收集系數等撞擊特性參數。通過對NACA0012翼型、GA-W（1）兩段翼型和某三段翼型的計算得到不同狀態下的水滴撞擊特性，計算結果表明，該方法與傳統方法相比具有計算效率高、結果可靠、通用性好等優點。

結冰；水滴撞擊特性；目標擴散追蹤方法；拉格朗日法

0 引 言

結冰會改變飛機的繞流流場，破壞空氣動力學性能，影響飛機的操縱性和穩定性，危害飛行安全，嚴重時能導致機毀人亡的嚴重事故［1］。我國幅員遼闊，氣象條件非常復雜，飛機結冰現象比較常見，隨著民用航空運輸業的發展、飛機飛行密度的提高以及對飛機全天候飛行的要求，飛行中遭遇到結冰氣象條件的幾率更是大幅提高，結冰已成為飛行安全事故的主要隱患之一［2］。

水滴撞擊特性研究是飛機結冰研究的主要內容之一，是結冰預測以及防／除冰系統設計的基礎［3］，水滴撞擊特性的計算方法主要有拉格朗日和歐拉兩種方法［4］。拉格朗日法是在模擬繞流流場的基礎上，采用差分法逐個求解水滴運動方程，得到各水滴的運動軌跡，進而求得水滴撞擊區域和形體表面上的局部水收集系數等參數［5］，確定水滴的撞擊特性。該方法形式簡單，方法成熟，尤其是模擬SLD（過冷大水滴）成冰過程時［6-7］，它能比歐拉方法更直觀、方便的體現大水滴在運動、撞擊過程中的變形、破碎、飛濺等動力學行為的局部運動細節，因此被廣泛應用于積冰研究中。

拉格朗日方法求解水滴運動方程時，每個時間步都必須計算水滴所屬網格單元，現在的處理方法通常是采用全場網格遍歷的方法搜索水滴所處網格單元［8］，這種方法計算量大，水滴是否撞擊到物面的判斷過程繁雜，對于三維外形和不連通的多段結冰面（如多段翼型）的撞擊判斷尤其難度大、方法通用性差［9-10］。由于對每個水滴在各個時間步上都進行一次全場網格遍歷效率太低，本文提出的目標擴散追蹤方法將水滴當成追蹤目標，沿著水滴飛行的軌跡進行擴散追蹤，采用射線求交法判斷水滴與網格單元的內外關系，大大提高了計算效率；并通過簡潔、高效的結構化邊界信息管理方式，成功解決了傳統方法計算不連通的多段結冰面通用性差、判斷水滴是否與物面相撞困難等問題。本文方法為飛機防除冰系統的設計中飛機部件水滴撞擊特性提供了一種快速、準確的計算方法。

1 計算方法

水滴撞擊特性的拉格朗日計算過程分兩步：首先，用CFD方法求解繞流空氣流場；然后，在流場解的基礎上求解水滴軌跡，進而獲得水滴撞擊極限、水滴收集率等水滴撞擊特性參數。

1.1 空氣流場計算

空氣流場積分形式控制方程為［11］：

采用Van Leer迎風格式進行離散，應用LUSGS隱式格式進行時間推進，結合SA一方程湍流模型獲得繞流流場定常解［12］。

1.2 水滴撞擊特性計算

1.2.1 水滴軌跡運動方程的求解

建立水滴運動方程時，作如下假設：①水滴體積很小，不會影響流場的性質；②水滴在運動過程中不碰撞、不合并、不分裂，水滴尺寸保持不變；③水滴在運動過程中和周圍空氣不發生質量、熱交換，溫度、粘性、密度等介質參數保持不變［13-14］。

考慮作用在水滴上的重力、浮力和阻力，根據牛頓第二定律，水滴軌跡運動方程可以寫成［15］：

每個水滴軌跡運動方程均可當成一階常微分方程的初值問題來求解，從遠處開始計算水滴的運動軌跡，采用四步龍格-庫塔方法，由i時刻的水滴及其繞流信息得到i＋Δi時刻的水滴位置、速度，水滴繞流信息由所處網格單元各節點信息插值得到，當水滴流出邊界或撞擊到部件時開始計算下一水滴軌跡，直到完成所有水滴的軌跡計算。

1.2.2 水滴所處單元的搜索

水滴軌跡計算過程中，每個時間步上都要進行水滴網格單元的搜索，搜索方法的快慢直接影響到結冰模擬算法的效率。特別在過冷大水滴結冰時，為精確捕捉大水滴在運動、撞擊過程中的動力學細節，需計算的水滴數量運比常規小水滴結冰模擬時要多，提高水滴所處單元的計算效率尤為重要。

傳統全場遍歷方法的效率太低，本文提出一種目標擴散追蹤方法來計算水滴所在網格單元，采用射線求交法判斷水滴與網格單元的內外關系，以前一時間步水滴坐標為起點，以本時間步水滴坐標為目標位置進行擴散追蹤以計算水滴所處的網格單元。計算發現，由于步長很小，水滴在一個時間步的移動距離基本都在一至兩個網格單元之間，除了在物面附近極密的粘性網格外，其他絕大部分位置只需要進行一、兩次簡單的位置判定就能獲得水滴所屬網格單元。該方法將傳統方法中的全流場的遍歷過程簡化為幾次簡單的線段相交計算，大大節省了計算時間，極大地提高了搜索效率。

以二維情況下對一個水滴進行擴展目標搜索為例，其過程如下：

1）水滴初始位置坐標為（xA，yA），全場搜索該點所處網格單元（I0，J0）。

2）由網格單元（I0，J0）各節點處的值插值得到水滴繞流信息，求解水滴軌跡方程，得到Δi時間后水滴位置坐標（xB，yB）。

3）計算AB與網格單元（I0，J0）四條邊的交點。若全不相交，則B點仍處于網格單元（I0，J0）中，水滴新位置A的坐標值為：xA＝xB，yA＝yB，返回2）；若與某條邊相交，該邊不屬于塊邊界則進行第4步，該邊屬于塊邊界則進行第5步。

4）水滴流入與該邊相鄰的網格單元中，從相交點往B方向延伸一個很小的量，該點坐標為（xC，yC），水滴新位置A的坐標值為：xA＝xC，yA＝yC，返回3）。

5）判斷該邊的邊界條件：若為對接邊界，則新的網格單元為相鄰網格塊上對應的單元，其序號采用計算空氣流場對接邊界條件時的方式得到，從相交點往B方向延伸一個很小的量，該點坐標為（xC，yC），水滴新位置A的坐標值為：xA＝xC，yA＝yC，返回3）；若為遠場邊界，則水滴流出計算域，結束該點追蹤；若為物面邊界，此時水滴撞擊到形體上，記錄撞擊位置，結束該點追蹤。

1.2.3 水滴撞擊物面的判定

水滴撞擊點的傳統計算方法需要在每個時間步判定水滴與形體的相對關系，并計算水滴軌跡線與所有物面網格線是否有交點，過程十分繁瑣、計算量龐大，且很難用一個通用程序來處理不同復雜幾何構型。采用目標擴散追蹤方法計算水滴所處網格單元時，在上述方法第5步的計算過程中就直接得到了水滴與物面的撞擊點，不需進行額外的特殊處理，而這需要依賴于本文所采用的由上而下逐級結構化管理的邊界信息儲存方式。

將網格塊中各面的邊界信息以“區”為單位存儲，按照“塊”、“面”、“區”順序對各網格塊的邊界信息進行結構化管理，每個存儲單位中的信息包含該“區”所屬網格塊、所處邊界面、該“區”序號、邊界類型、網格點起始點與結束點等信息，將搭接邊界分為左右兩部分，按照相應的塊、起始結束點序號以及搭接方向一一對應。這種逐級結構化的邊界信息管理方式能保證程序不受計算網格拓撲結構的限制、可靈活應用于不同復雜幾何外形。

通過結構化的邊界信息儲存方式，水滴軌跡計算過程中可快速得到與水滴軌跡線相交的網格單元邊的信息，實現高效、準確的水滴碰撞計算。

1.3 水滴撞擊特性的確定

確定水滴運動軌跡之后，可進一步得到防冰部件的水滴撞擊極限、總收集系數、局部水收集系數等水滴撞擊特性參數［16］，用于防冰系統的設計。

2 算例與分析

通過不同外形水滴軌跡的數值模擬及其與相關實驗數據的對比對方法進行了檢驗和驗證。

2.1 方法效率測試

對規模為7880、14820、23760的NACA0012翼型網格，規模為17600的GA-W（1）兩段翼型網格（網格分為三塊），規模為14720的某三段翼型網格（網格分為四塊）分別采用直接搜索方法和本文方法進行300個水滴點的跟蹤計算，直接搜索方法耗時分別為14．4s、26．1s、44．7s、41．3s、35．4s，本文方法耗時分別為1．3s、1．5s、1．7s、1．4s、1．3s。計算結果表明擴展目標追蹤方法能大幅提高計算效率，且計算效率受網格規模、拓撲結構的影響不大，網格越密、拓撲結構越復雜，其優勢就越明顯。

2.2 方法準確性測試

圖1 NACA0012翼型撞擊極限及水滴軌跡Fig.1 Trajectories and impingement limit of water droplets for NACA0012 airfoil

對長度為0．5334的NACA0012翼型進行了水滴軌跡計算，氣象條件如下：M＝0.32，p＝89867Pa，LWC＝0．55g／m3。case1中MVD＝15μm；case2中MVD＝40μm。圖1為采用本文方法計算得到的case1、case2下的水滴軌跡以及撞擊極限，圖2為計算局部收集系數（Beta）與實驗值的對比，s為翼型上的水滴撞擊點到駐點的距離，上表面為正，下表面為負，s＝0處為局部收集系數達到最大值的位置。計算結果與參考文獻［17］中給出的目前公認結冰計算能力最強的軟件Lewice的結果基本上一致，表明本方法的計算結果是準確可靠的。

圖2 水滴收集系數與試驗結果的對比Fig.2 Water collection coefficient compared with the experiment data

2.3 方法通用性測試

在M＝0．27，α＝2．0°，Re＝0．1E＋7，MVD＝20．0μm，p＝1．013×105Pa，T＝262．9K條件下，分別對帶30%子翼的GA-W（1）兩段翼型和某三段翼型進行水滴軌跡計算，水滴運動軌跡及撞擊區域如圖3、圖4所示；各段的水滴收集系數分布規律如圖5、圖6所示，計算結果合理，這表明本方法可對各種不連通幾何外形進行水滴撞擊特性的計算。

圖3 GA-W（1）兩段翼型水滴軌跡及撞擊極限Fig.3 Trajectories and impingement limit of water droplets for GA-W（1）two-element airfoil

圖4 某三段翼型水滴軌跡及撞擊極限Fig.4 Trajectories and impingement limit of water droplets for three-element airfoil

圖5 GA-W（1）兩段翼型水滴收集系數Fig.5 Water collection coefficient for GA-W（1）two-element airfoil

圖6 某三段翼型水滴收集系數Fig.6 Water collection coefficient for three-element airfoil

3 結 論

本文提出的水滴撞擊特性的拉格朗日計算方法基于目標擴散追蹤方法以及邊界信息結構化管理，克服了常用方法計算大量水滴的撞擊特性時計算量大、水滴撞擊判斷過程繁雜、應用于復雜外形的通用性差等缺點，具有計算效率高、計算結果可靠、對復雜外形通用性好等優點，為飛機防除冰系統的設計中飛機部件水滴撞擊特性提供了一種快速、準確的計算方法。

［1］ISAAC G A，et al．Recent Canadian research on aircraft inflight icing［J］．Canadian Aeronautics and Space Journal，2001，47（3）：213-222．

［2］HONSEK R，HABASHI W G，AUBE M S．Eulrian modeling of in-flight icing due to supercooled large droplets［J］．Journal of Aircraft，2008，45（4）：1290-1296．

［3］SUN Z G，ZHU C X，ZHU C L．Development of s of tware for aircraft icing simulation［J］．Computer Simulation，2012，29（4）：104-107．（in Chinese）

孫志國，朱程香，朱春玲．飛機結冰數值仿真軟件開發［J］．計算機仿真，2012，29（4）：104-107．

［4］ZHANG D L，YANG X，ANG H S．Numerical simulation of supercooled water droplets impingement on icing surfaces［J］．Journal of Aerospace Power，2003，18（1）：197-201．（in Chinese）

張大林，楊曦，昂海松．過冷水滴撞擊結冰表面的數值模擬［J］．航空動力學報，2003，18（1）：197-201．

［5］COSSALI G E，COGHE A，MARENGO M．The impact of a single drop on a wetted solid surface［J］．Experiments in Fluids，1997，22（6）：463-472．

［6］ZHANG C，KONG W L，LIU H．An investigation on the breakup model for icing simulation of supercooled large droplets［J］．ACTA Aerodynamica Sinica，2013，31（2）：144-150．（in Chinese）

張辰，孔維梁，劉洪．大粒徑過冷水滴結冰模擬破碎模型研究［J］．空氣動力學學報，2013，31（2）：144-150．

［7］BOURGAULT Y，HABASHI W G，DOMPIERRE J，et al．An Eulerian approach to supercooled droplets impingement calculations［R］．AIAA-97-0176，

［8］YANG Q，CHANG S N，YUAN X G．Study on numerical method for determining the droplet trajectories［J］．ACTA Aeronautica et Astronautica Sinica，2002，23（2）：173-176．（in Chinese）

楊倩，常士楠，袁修干．水滴撞擊特性的數值計算方法研究［J］．航空學報，2002，23（2）：173-176．

［9］YI X，WANG K C，GUI Y W．Study on Eulerian method for icing collection efficiency computation and its application［J］．ACTA Aerodynamica Sinica，2010，28（5）：596-602．（in Chinese）

易賢，王開春，桂業偉．結冰面水滴收集率歐拉計算方法研究及應用［J］．空氣動力學學報，2010，28（5）：596-602．

［10］TRUJILLO M F，MATHEWS W S，LEE C F，et al．Modeling and experiments of impingement and atomization of a liquid spray on a wall［J］．International Journal of Engine Research，2000，1（1）：87-105．

［11］SHIN J，BERKOWITZ B，GHEN H H，et al．Prediction of ice shapes and their effect on airfoil drag［J］．Journal of Aircraft，1994，31（2）：301-316．

［12］ZHOU Z H．Navier-Stokes analysis in complex configurations and icing simulation on airplanes［D］．Xi′an：Northwestern Polytechnical University，2011．

［13］CHANG S N，YANG Q M，LI Y．Quasi-steady numerical simulation of ice accretion on airfoil［J］．ACTA Aerodynamica Sinica，2011，06：302-308．（in Chinese）

常士楠，楊秋明，李延．翼型表面結冰準定常數值模擬［J］．空氣動力學學報，2011，06：302-308．

［14］PAPADAKIS M，RACHMAN A，WONG S C，et al．Waterimpingement experiments on a NACA 23012 airfoil with simulated glaze ice shapes［R］．AIAA 2004-0565，2004．

［15］YI X，ZHU G L，WANG K C，et al．Numerically simulating of ice accretion on airfoil［J］．ACTA Aerodynamica Sinica，2002，20（4）：428-433．（in Chinese）

易賢，朱國林，王開春，等．翼型積冰的數值模擬［J］．空氣動力學學報，2002，20（4）：428-433．

［16］YI X，ZHU G L．Computation of glaze ice accretion on airfoil［J］．ACTA Aerodynamica Sinica，2004，22（4）：490-493．（in Chinese）

易賢，朱國林．考慮傳質傳熱效應的翼型積冰計算［J］．空氣動力學學報，2004，22（4）：490-493．

［17］YVES B．A finite element method study of eulerian dropletes impingement models［J］．International Journal for Numerical Methods in Fluids，1999，29（4）：429-449．

An efficient method to simulate water droplet trajectory and impingement

ZHOU Zhihong1，YI Xian1，GUI Yewei1，LI Fengwei2
（1.State Key Laboratory of Aerodynamics，China Aerodynamics Research and Development Center，Mianyang 621000，China；2.Aeronautics School，Northwestern Polytechnical University，Xi′an 710072，China）

Ice formation on aircraft is of great safety concern because icing may lead the aircraft to a dangerous situation quickly．Numerical method for determining the trajectories of the water droplets is the base of the icing research and the design of the anti-icing system．The traditional Lagrangian method for simulating water droplet trajectory and impingement has some shortages such as lower computation efficiency and poor applicability currency in peculiar geometry．In order to overcome these shortages，an improved Lagrangian method was developed．It′s based on the calculation of the flow field around icing surface，using a method of searching extended target to calculate the droplet location and determining the insert value of flow field information at the droplet location in a grid cell，the boundary of the droplet trajectories are judged based on a manage system of boundary information with framework structure．The method can overcome theses shortages meanwhile maintaines the robustness of original method．Local droplet collection and impingement efficiency at NACA0012 and multi-element airfoils are calculated with this method in order to verify its correctness．The results show that the improved method is very efficient，reliable，and robust．Furthermore，it can be used in different geometries directly．

ice；the trajectories of the water droplets；method for searching extended target；Lagrangian method

V211.3

Adoi:10.7638／kqdlxxb-2012.0179

0258-1825（2014）05-0712-05

2012-10-31；

2013-01-01

國家自然科學基金（11172314）；中國博士后基金（2012M512065）；四川省國際合作計劃（12064HH0042）

周志宏（1981-），男，湖南省漣源人，博士，研究方向：計算流體力學、飛機結冰．E-mail：zhouzhihong029＠163．com

周志宏，易 賢，桂業偉，等．水滴撞擊特性的高效計算方法［J］．空氣動力學學報，2014，32（5）：712-716．

10．7638／kqdlxxb-2012．0179． ZHOU Z H，YI X，GUI Y W，et al．An efficient method to simulate water droplet trajectory and impingement［J］．ACTA Aerodynamica Sinica，2014，32（5）：712-716．