易 賢，王開春，馬洪林，朱國林

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

0 引 言

風(fēng)力機(jī)在寒冷的氣候（溫度低于0℃）環(huán)境下運(yùn)行時(shí)，如果遇到含有過冷水滴的氣流，或者遭遇雨雪天氣，在風(fēng)力機(jī)表面（包括葉片、輪轂、機(jī)艙、塔架、風(fēng)速風(fēng)向儀等風(fēng)力機(jī)所有部件表面）會(huì)發(fā)生結(jié)冰現(xiàn)象。結(jié)冰會(huì)改變?nèi)~片的氣動(dòng)外形，破壞風(fēng)力機(jī)的流場(chǎng)特性，使風(fēng)能利用系數(shù)降低，減小出力；同時(shí)，結(jié)冰還會(huì)引起額外的過載及振動(dòng)，冰載荷的不平衡會(huì)增加機(jī)組部件的疲勞，輕者，使風(fēng)力機(jī)停止運(yùn)行，帶來經(jīng)濟(jì)損失，重者，造成葉片損壞甚至導(dǎo)致整個(gè)風(fēng)機(jī)倒塌，危及人們的生命財(cái)產(chǎn)安全。鑒于風(fēng)力機(jī)結(jié)冰的危害，近年來人們?cè)絹碓街匾曔@方面的研究［1］。

風(fēng)力機(jī)結(jié)冰包括沉積結(jié)冰（Precipitation Icing，主要由凍雨和濕雪引起）和過冷水滴結(jié)冰（In－Cloud Icing，主要由過冷水滴撞擊在風(fēng)力機(jī)表面所引起）兩種形式［2］，其中過冷水滴結(jié)冰是低溫氣象條件下常見的一種結(jié)冰形式。水滴收集率（也稱水滴局部收集系數(shù)）定義為風(fēng)力機(jī)某局部區(qū)域?qū)嶋H所收集的水量與該區(qū)域可能收集的水量最大值之比，它表征了部件表面的水滴撞擊范圍以及撞擊區(qū)域內(nèi)水量的分布。水滴收集率是結(jié)冰研究中最重要的參數(shù)之一，只有獲得水滴收集率，才能進(jìn)行結(jié)冰預(yù)測(cè)、探測(cè)以及防／除冰系統(tǒng)設(shè)計(jì)。對(duì)于過冷水滴結(jié)冰，常用的研究手段有工程估算、結(jié)冰風(fēng)洞試驗(yàn)和數(shù)值計(jì)算。其中試驗(yàn)研究是早期主要采用的方法［3－4］，近年來，開始陸續(xù)出現(xiàn)用數(shù)值計(jì)算的手段開展研究的報(bào)道［1，5－7］，由于旋轉(zhuǎn)條件下風(fēng)力機(jī)流場(chǎng)的復(fù)雜性，這些計(jì)算多采用拉格朗日法計(jì)算水滴收集率，或者采用葉片的截面翼型進(jìn)行二維簡化計(jì)算。本文在前期發(fā)展的飛機(jī)結(jié)冰三維計(jì)算方法基礎(chǔ)上［8－10］，建立了適合于大型風(fēng)力機(jī)結(jié)冰過程水滴收集率計(jì)算的三維數(shù)值方法和計(jì)算程序，并對(duì)某1.5MW級(jí)水平軸風(fēng)力機(jī)的水滴收集率進(jìn)行了計(jì)算，獲得了不同直徑水滴在風(fēng)力機(jī)表面的撞擊特征以及水滴收集率在葉片上的分布規(guī)律，為開展風(fēng)力機(jī)結(jié)冰防護(hù)研究打下了較好基礎(chǔ)。

1 水平軸風(fēng)力機(jī)空氣流場(chǎng)計(jì)算方法

風(fēng)力機(jī)的空氣繞流場(chǎng)是典型的包含旋轉(zhuǎn)機(jī)械的流場(chǎng)，對(duì)于這類流場(chǎng)計(jì)算，有單旋轉(zhuǎn)坐標(biāo)系方法、多參考坐標(biāo)系（Multiple Reference Frame，簡稱 MRF）方法［11］、滑移網(wǎng)格方法、動(dòng)態(tài)重疊網(wǎng)格方法等，其中，當(dāng)流場(chǎng)中同時(shí)存在旋轉(zhuǎn)部件和非旋轉(zhuǎn)部件（如轉(zhuǎn)子／定子）時(shí)，MRF方法由于其魯棒、簡便、高效的特點(diǎn)，在工程中得到了廣泛的使用，一些CFD商業(yè)軟件如FLUENT、STAR－CD、CFX等，都集成了該方法。

MRF方法的基本思路為：將計(jì)算區(qū)域分為旋轉(zhuǎn)區(qū)和非旋轉(zhuǎn)區(qū)，在非旋轉(zhuǎn)區(qū)，流動(dòng)的控制方程為非旋轉(zhuǎn)坐標(biāo)系（如慣性坐標(biāo)系）下的方程，在旋轉(zhuǎn)區(qū)域，流動(dòng)的控制方程為旋轉(zhuǎn)坐標(biāo)系下的控制方程。對(duì)不同區(qū)域的流動(dòng)分別進(jìn)行求解，在區(qū)域的交界面上，要求流場(chǎng)的參數(shù)完全一致。

本文在原有的非旋轉(zhuǎn)流場(chǎng)計(jì)算研究的基礎(chǔ)上，建立了計(jì)算風(fēng)力機(jī)流場(chǎng)的MRF方法，現(xiàn)概述如下。

1.1 空氣流場(chǎng)計(jì)算的控制方程

非旋轉(zhuǎn)區(qū)域空氣流場(chǎng)計(jì)算的控制方程為低速粘流的時(shí)均N－S方程，其通用形式為：

方程（1）中，φ為輸運(yùn)變量，ρa(bǔ)為空氣密度，ua為空氣速度，Γφ為擴(kuò)散系數(shù)，qφ為源項(xiàng)，φ、Γφ和qφ取不同的值，可代表空氣的連續(xù)性方程、動(dòng)量方程和其他標(biāo)量（如湍動(dòng)能等）的輸運(yùn)方程，方程中各項(xiàng)的物理意義及具體表達(dá)式可參見文獻(xiàn)［8］，此處不再詳述。

在旋轉(zhuǎn)區(qū)域內(nèi)，控制方程也可以寫成輸運(yùn)方程的形式：

方程（2）中，uar為空氣在旋轉(zhuǎn)坐標(biāo)系中的速度，qφr為旋轉(zhuǎn)坐標(biāo)系中的源項(xiàng)。當(dāng)φ＝1、qφr＝0時(shí)，方程（2）為連續(xù)性方程，當(dāng)φ＝uar時(shí)，方程（2）為動(dòng)量方程，此時(shí)源項(xiàng)的表達(dá)式為：

其中，ω為旋轉(zhuǎn)角速度，r為流體質(zhì)點(diǎn)在旋轉(zhuǎn)坐標(biāo)系中的矢徑。

在旋轉(zhuǎn)區(qū)域和非旋轉(zhuǎn)區(qū)域的交界面上，要求流場(chǎng)參數(shù)一致，兩個(gè)坐標(biāo)系內(nèi)的速度關(guān)系為：

1.2 控制方程的有限體積離散

采用有限體積法離散求解方程（1）、（2）。以點(diǎn)P為中心點(diǎn)的六面體為例，方程中各項(xiàng)的離散方法為：

（1）對(duì)流項(xiàng)的離散。

式（5）中，nb＝e，w，n，s，t，b，分別代表以 P點(diǎn)為中心的控制體的六個(gè)面，Snb為各面對(duì)應(yīng)的面積，對(duì)于空氣相的控制方程有：

控制體界面上的變量值fnb用迎風(fēng)插值和線性插值相結(jié)合的方法計(jì)算，表達(dá)式為：

其中為迎風(fēng)插值的量，為線性插值的量，ε為混合因子，0≤ε≤1。

（2）源項(xiàng)的離散。輸運(yùn)方程的源項(xiàng)對(duì)不同的方程都有不同的表達(dá)式，為了使其離散化表達(dá)式盡可能逼近源項(xiàng)本身，加強(qiáng)代數(shù)方程的主對(duì)角優(yōu)勢(shì)，提高代數(shù)方程組求解的穩(wěn)定性，對(duì)源項(xiàng)采用線性化方式做進(jìn)行處理，即令

則源項(xiàng)的離散形式為：

式中δV為控制體的體積。

（3）時(shí)間項(xiàng)離散。方程（1）、（2）的時(shí)間項(xiàng)采用如下二階精度的隱式離散：

其中，上標(biāo)n＋1、n和n－1分別代表t＋Δt時(shí)刻、t時(shí)刻和t－Δt時(shí)刻的值。

2 水滴項(xiàng)的控制方程及求解方法

水滴收集率計(jì)算在空氣流場(chǎng)計(jì)算完成之后進(jìn)行，主要有兩種方法［12－13］，一種是拉格朗日法，即用拉格朗日的觀點(diǎn)，計(jì)算跟蹤流場(chǎng)中單個(gè)水滴的運(yùn)動(dòng)，并判斷其是否與物面碰撞，從而得到物面各處的水滴收集量及其分布；另一種是歐拉法，該方法將水滴看作與空氣單向作用的、在空間連續(xù)分布的液態(tài)相，通過流體力學(xué)的觀點(diǎn)求解水滴相的控制方程，得到流場(chǎng)中水滴容積分?jǐn)?shù)分布，進(jìn)而獲得水滴收集率。鑒于拉格朗日法在三維計(jì)算中的不足，本文發(fā)展了三維水滴收集率計(jì)算的歐拉方法。

2.1 水滴相流場(chǎng)計(jì)算的控制方程

水滴相方程的求解統(tǒng)一在慣性坐標(biāo)系中進(jìn)行，引入水滴容積分?jǐn)?shù)α，其定義為空間微團(tuán)中水滴相所占的體積比例，則可以建立水滴相的控制方程，包括連續(xù)方程和動(dòng)量方程，分別為［8］：

方程（12）、（13）中，ud為水滴速度，ρd為水滴密度，g為重力加速度，K為慣性因子，其表達(dá)式為：

式（14）中，μa為空氣動(dòng)力粘性系數(shù)，d為水滴直徑，CD為水滴阻力系數(shù)，Re為相對(duì)雷諾數(shù)，其表達(dá)式為：

水滴運(yùn)動(dòng)過程中，阻力隨相對(duì)雷諾數(shù)變化而變化，本文采用如下公式確定阻力［8］：

與空氣控制方程類似，水滴相控制方程可以統(tǒng)一寫成輸運(yùn)方程的形式：

式中，qφ為源項(xiàng)，φ取1、ud、vd或wd分別代表連續(xù)方程和x、y、z方向的動(dòng)量方程。

2.2 水滴相控制方程的離散方法

對(duì)于方程（17），對(duì)流項(xiàng)和源項(xiàng)的離散方法與方程（1）、（2）一致，時(shí)間項(xiàng)采用一階顯式離散，即：

2.3 壁面邊界條件

與空氣流場(chǎng)計(jì)算在壁面采用無滑移壁面邊界條件不同，對(duì)于水滴相計(jì)算，采用壁面吸入邊界條件，即如果水滴與物面碰撞，則認(rèn)為水滴從碰撞點(diǎn)流出。

3 水滴收集率的計(jì)算

由于水滴容積分?jǐn)?shù)α較?。?0－6量級(jí)），可以認(rèn)為空氣和水滴是單向作用，即只考慮空氣對(duì)水滴的作用，忽略水滴對(duì)空氣的作用。因此，水滴收集率的計(jì)算步驟可概括為：首先，計(jì)算空氣場(chǎng)，本文采用SIMPLE方法計(jì)算空氣場(chǎng)，湍流模型為標(biāo)準(zhǔn)k－ε湍流模型［14］；其次，在得到空氣流場(chǎng)分布的基礎(chǔ)上，求解水滴相控制方程；最后，水滴收集率β可在獲得當(dāng)?shù)厮稳莘e分?jǐn)?shù)α和水滴與物面的相對(duì)速度udr之后，由以下表達(dá)式得到：

其中，α∞為遠(yuǎn)場(chǎng)水滴容積分?jǐn)?shù)，u∞為遠(yuǎn)場(chǎng)來流速度，n為物面碰撞點(diǎn)處的單位法線向量。

由于udr＝ud－ω×r，式（19）可以寫成：

從式（20）可以看出，旋轉(zhuǎn)狀態(tài)下的水滴收集率β，由兩部分組成，一部分是風(fēng)力機(jī)旋轉(zhuǎn)速度為0時(shí)的水滴收集率β0，另一部分是由于旋轉(zhuǎn)而產(chǎn)生的水滴收集率βr。在慣性坐標(biāo)系下，水滴收集率的物理意義定義為：物面某局部區(qū)域所收集的水量與該區(qū)域所收集的水量最大值之比，因此β0值應(yīng)該在［0，1］區(qū)間，而βr的值則與物面某區(qū)域的旋轉(zhuǎn)半徑有關(guān)，半徑越大，βr的值也越大。因此，在旋轉(zhuǎn)條件下，由式（19）計(jì)算得到的β值將會(huì)出現(xiàn)大于1的情況。

4 大型風(fēng)力機(jī)水滴收集率計(jì)算

根據(jù)以上數(shù)值方法，開發(fā)了相應(yīng)的計(jì)算程序，并針對(duì)某1.5MW級(jí)的水平軸風(fēng)力機(jī)進(jìn)行了應(yīng)用計(jì)算。

4.1 計(jì)算外形及網(wǎng)格

計(jì)算外形選用某1.5MW級(jí)的水平軸風(fēng)力機(jī)，半徑R＝41m，忽略塔架的影響。如圖1，坐標(biāo)原點(diǎn)取為風(fēng)力機(jī)輪轂中心，坐標(biāo)軸方向?yàn)椋簒軸與遠(yuǎn)場(chǎng)來流方向一致，y軸垂直向上，z軸按右手系確定。計(jì)算采用多塊對(duì)接網(wǎng)格，網(wǎng)格規(guī)模為1800萬。圖2顯示了風(fēng)力機(jī)表面輪轂附近的網(wǎng)格。

圖1 計(jì)算外形及坐標(biāo)方向Fig.1 Wind turbine configuration and frame

圖2 輪轂附近表面網(wǎng)格 Fig.2 Surface grid near the hub

4.2 風(fēng)力機(jī)空氣流場(chǎng)計(jì)算結(jié)果

本文計(jì)算了兩種來流條件下風(fēng)力機(jī)的氣動(dòng)特性，并與設(shè)計(jì)值進(jìn)行了對(duì)比，計(jì)算條件為：

計(jì)算條件1：來流速度8m／s，風(fēng)輪轉(zhuǎn)速16.7r／min，漿距角0°，葉尖速比為9。

計(jì)算條件2：來流速度11m／s，風(fēng)輪轉(zhuǎn)速17.4r／min，漿距角3.2°，葉尖速比為6.8。

圖3給出了速度為8m／s時(shí)葉片上的壓力計(jì)算結(jié)果。圖3（a）給出了葉片上壓力分布云圖，圖3（b）～圖3（f）顯示的是葉片上R＝4m、13m、22m、31m和40m五個(gè)不同位置截面上的壓力系數(shù)曲線。計(jì)算顯示，在靠近輪轂附近（R＝4m），由于葉片截面還不是翼型形狀，迎風(fēng)面和背風(fēng)面的壓力呈現(xiàn)相反的變化趨勢(shì)，在30%弦長處，吸力面在迎風(fēng)和背風(fēng)面之間發(fā)生轉(zhuǎn)換；葉片中部（R＝22m）一直到葉尖附近（R＝40m），壓力分布呈相似的曲線輪廓。

圖3 葉片表面壓力分布（v＝8m／s）Fig.3 Pressure distribution on the blade（v＝8m／s）

表1給出了兩種速度條件下，本文計(jì)算的軸功率與設(shè)計(jì)值的比較，可見本文計(jì)算的功率與設(shè)計(jì)值接近，計(jì)算值比設(shè)計(jì)值只分別低0.6%和1.3%。

表1 計(jì)算的軸功率與設(shè)計(jì)值的對(duì)比Table 1 Comparison of computed power and designed power

采用固定8m／s的來流速度、改變風(fēng)輪轉(zhuǎn)速的方式，計(jì)算了不同葉尖速比條件下風(fēng)輪的氣動(dòng)特性。圖4（a）～圖4（c）分別給出了λ＝9、4和2三種葉尖速比條件下，吸力面上的極限流線分布情況。圖中顯示，葉尖速比λ＝9時(shí)，吸力面的流動(dòng)為附著流，未出現(xiàn)分離，整個(gè)葉片表面近似為二維弦向流動(dòng)，葉尖速比λ＝4時(shí)，在葉片的中部首先出現(xiàn)分離，當(dāng)葉尖速比降至2時(shí)，整個(gè)吸力面上的流動(dòng)已完全分離。

圖4 吸力面極限流線分布Fig.4 Limiting streamlines distribution on suction surface

4.3 風(fēng)力機(jī)水滴收集率計(jì)算結(jié)果

在風(fēng)力機(jī)空氣流場(chǎng)計(jì)算的基礎(chǔ)上，開展了水滴收集率計(jì)算，對(duì)比了不同直徑水滴的收集率在風(fēng)力機(jī)表面的分布情況。計(jì)算條件為：來流速度8m／s，密度1.225kg／m3，環(huán)境壓力101325Pa，風(fēng)輪轉(zhuǎn)速16.7r／min，漿距角0°，葉尖速比為9，水滴直徑選取了三種，從大到小分別為100μm、40μm和20μm。

圖5～圖8給出了水滴直徑d＝100μm時(shí)的計(jì)算結(jié)果。圖5為風(fēng)力機(jī)表面的水滴收集率分布云圖，由圖可見，水滴主要撞擊在輪轂表面和葉片前緣；圖6為收集率云圖在輪轂附近的放大，圖7為單個(gè)葉片表面的分布云圖，圖8給出的是葉片表面不同展向位置的水滴收集率曲線。計(jì)算顯示：（1）從葉片根部到葉尖，水滴收集率逐漸增大，葉尖附近（R＝40m）的最大收集率達(dá)到7.5以上，而葉片根部附近（R＝4m）的收集率僅為0.2左右；（2）輪轂上的水滴收集率極?。ú坏?.1），小于葉片根部的收集率。

根據(jù)以上水滴收集率計(jì)算結(jié)果，可以對(duì)輪轂表面的結(jié)冰厚度進(jìn)行估算，物體表面某處結(jié)冰的厚度h可以表達(dá)為［8］：

圖5 風(fēng)力機(jī)表面水滴收集率分布云圖（d＝100μm）Fig.5 Contours of collection efficiency（d＝100μm）

圖6 輪轂附近的水滴收集率分布云圖（d＝100μm）Fig.6 Collection efficiency near hub（d＝100μm）

圖7 葉片表面的水滴收集率分布云圖（d＝100μm）Fig.7 Collection efficiency on blade（d＝100μm）

圖8 葉片上不同位置水滴收集率曲線（d＝100μm）Fig.8 Collection efficiency distribution on blade（d＝100μm）

其中f為凍結(jié)系數(shù)（0≤f≤1），β為當(dāng)?shù)氐乃问占剩琇WC為空氣中的液態(tài)水含量［15］，V為來流速度，dt為結(jié)冰時(shí)間，ρi為冰的密度。

根據(jù)式（20），輪轂表面β的最大值取為0.1，取f為1，V取為8m／s，LWC取0.5g／m3，ρi取700kg／m3，則1h的結(jié)冰厚度約為0.2cm，連續(xù)50h的結(jié)冰厚度也僅為10cm。由于結(jié)冰主要影響風(fēng)力機(jī)的氣動(dòng)特性和葉片的載荷強(qiáng)度及疲勞特性，因此在實(shí)際分析時(shí)，可以忽略輪轂的結(jié)冰，主要考慮葉片結(jié)冰的影響。

圖9、圖10給出了水滴直徑d＝40μm時(shí)，風(fēng)力機(jī)表面的水滴收集率分布；圖11、圖12給出的是水滴直徑為20μm的結(jié)果?？梢钥闯觯沃睆阶冃≈?，水滴撞擊特性的變化規(guī)律為：（1）水滴收集率隨水滴直徑變小而減?。唬?）不同直徑水滴在風(fēng)力機(jī)表面的撞擊區(qū)域以及水滴收集率在葉片上的分布特征保持類似。

圖9 輪轂附近的水滴收集率分布云圖（d＝40μm）Fig.9 Collection efficiency near hub（d＝40μm）

圖10 葉片上不同位置水滴收集率曲線（d＝40μm）Fig.10 Collection efficiency distribution on blade（d＝40μm）

圖11 輪轂附近的水滴收集率分布云圖（d＝20μm）Fig.11 Collection efficiency near hub（d＝20μm）

圖12 葉片上不同位置水滴收集率曲線（d＝20μm）Fig.12 Collection efficiency distribution on blade（d＝20μm）

以上結(jié)果說明對(duì)于過冷水滴在風(fēng)力機(jī)葉片表面的結(jié)冰，其基本規(guī)律有兩個(gè)，一是冰主要結(jié)在葉片前緣，二是葉尖結(jié)冰遠(yuǎn)比葉根處嚴(yán)重。因此，在進(jìn)行結(jié)冰分析、設(shè)計(jì)結(jié)冰探測(cè)方法和防除冰方案時(shí)，葉片前緣，尤其是葉尖附近的葉片前緣，需要重點(diǎn)關(guān)注。

5 結(jié) 論

建立了適合于大型風(fēng)力機(jī)結(jié)冰過程水滴收集率計(jì)算的三維數(shù)值方法和計(jì)算程序，并對(duì)某1.5MW級(jí)水平軸風(fēng)力機(jī)的水滴收集率進(jìn)行了計(jì)算，得到如下結(jié)論：

（1）本文的工作基于前期發(fā)展的飛機(jī)結(jié)冰三維計(jì)算方法，旋轉(zhuǎn)和非旋轉(zhuǎn)區(qū)的空氣流場(chǎng)、水滴相的控制方程均采用統(tǒng)一的有限體積離散方法，方便了計(jì)算程序的編制，對(duì)風(fēng)力機(jī)流場(chǎng)和水滴收集率的計(jì)算結(jié)果初步表明本文發(fā)展的算法以及相應(yīng)的計(jì)算程序是有效的。

（2）對(duì)某1.5MW級(jí)水平軸風(fēng)力機(jī)結(jié)冰過程水滴收集率的計(jì)算表明，水滴主要撞擊在輪轂表面和葉片前緣，輪轂上的水滴收集率極小，從葉片根部到葉尖，水滴收集率逐漸增加，且葉尖處收集率遠(yuǎn)比葉根處高。

（3）在進(jìn)行結(jié)冰分析、設(shè)計(jì)結(jié)冰探測(cè)方法和防除冰方案時(shí)，可以忽略輪轂和葉片根部附近結(jié)冰的影響，葉片前緣，尤其是葉尖附近的葉片前緣，需要重點(diǎn)關(guān)注。

（4）風(fēng)力機(jī)結(jié)冰過程涉及流動(dòng)、傳質(zhì)、傳熱和相變等多種物理現(xiàn)象的耦合，對(duì)該問題進(jìn)行三維數(shù)值模擬是一項(xiàng)復(fù)雜的工作，本文發(fā)展的風(fēng)力機(jī)流場(chǎng)和水滴收集率求解方法，為進(jìn)行三維結(jié)冰計(jì)算打下了較好的基礎(chǔ)。

［1］FORTIN G ，PERRON J.Wind turbine icing and deicing［R］.AIAA－2009－0274，2009

［2］FROHBOESE P，ANDERS A.Effects of icing on wind turbine fatigue loads［J］.JournalofPhysics：ConferenceSeries，2007，（75）：012061.

［3］JASINSKI W J，NOE S C，SELIG M S，et al.Wind turbine performance under icing conditions［R］.AIAA－1997－0977，1997.

［4］HOCHART C，F(xiàn)ORTIN G，PERRON J.Icing simulation of wind－turbine blades［R］.AIAA－2007－1373，2007.

［5］MAKKONEN L，LAAKSO T，MARJANIEMI M，et al.Modeling and prevention of ice accretion on wind turbines［J］.WindEngineering，2001，25（1）：3－21.

［6］BATTISTI L，F(xiàn)EDRIZZI R，RIALTI M，et al.A model for the design of hot－air based wind turbine ice prevention system［A］.Proceedings of the WREC［C］.Aberdeen，2005.

［7］朱程香，王瓏，孫志國，等.風(fēng)力機(jī)葉片翼型的結(jié)冰數(shù)值模擬研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2011，29（4）：522－528.

［8］易賢，王開春，桂業(yè)偉，等.結(jié)冰面水滴收集率歐拉計(jì)算方法研究及應(yīng)用［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2010，28（5）：596－601.

［9］易賢，桂業(yè)偉，朱國林.飛機(jī)三維結(jié)冰模型及其數(shù)值求解方法［J］.航空學(xué)報(bào)，2010，31（11）：2152－2158.

［10］易賢，桂業(yè)偉，朱國林，杜雁霞.運(yùn)輸機(jī)翼型結(jié)冰的計(jì)算和實(shí)驗(yàn)［J］.航空動(dòng)力學(xué)報(bào)，2011，26（4）：808－813.

［11］LUO J Y，ISSA R I，GOSMAN A D.Prediction of impeller－induced flows in mixing vessels using multiple frames of reference［J］.InIChemESymposiumSeries，1994，136：549－556.

［12］BRAGG M B，GREGOREK G M，SHAW R J.An analytical approach to airfoil icing［R］.AIAA－81－0403，1981.

［13］BOURGAULT Y，HABASHI W G，DOMPIERRE J，et al.An eulerian approach to supercooled droplets impingement calculations［R］.AIAA－97－0176，1997.

［14］ZHU G L，KRONAST M.The calculation of ground effect on a car flow field using two dimensional Navier－Stokes equations［J］.ACTAAerodynamicSinica，1993，11（1）：35－40.

［15］易賢.飛機(jī)積冰的數(shù)值計(jì)算與積冰試驗(yàn)相似準(zhǔn)則研究［D］.綿陽：中國空氣動(dòng)力研究與發(fā)展中心，2007.