基于數值模擬的NSDBD等離子體激勵器防冰特性

賈韞澤，桑為民，蔡旸

西北工業大學 航空學院，西安 710072

飛行器在積冰氣象條件下飛行時，由于過冷水滴的撞擊在其表面形成的積冰會使氣動性能下降。在積冰條件下，防冰系統保障在飛行安全方面尤為重要。

對于飛機積冰的研究，通常包括：積冰預測、防/除冰措施研究、積冰對空氣流場影響等。在積冰預測研究中水滴的收集情況是重要的一部分，Durst等[1]使用拉格朗日方法和歐拉方法求解水滴碰撞情況，計算發現兩種方法得到的結果基本一致。基于控制體內的質量與能量平衡，Messinger[2]建立了積冰質量、能量平衡方程。基于Messinger模型,Macarthur[3]考慮了溢流水的影響，對霜冰與光冰條件下的冰型進行了預測，與試驗對比良好。隨著研究的增多，出現了許多商用的積冰預測軟件，如：LEWICE、ONERA 2D/3D、FENSAP-ICE等。在防冰措施方面，熱防冰方法已經得到運用。通過熱源的不同，現有熱防冰主要分為氣熱防冰與電熱防冰。氣熱防冰主要是通過引入噴氣式發動機的熱空氣到防冰區，對防冰區進行加熱。電防冰是直接通過布置好的發熱元件進行加熱。電熱防冰已在波音787客機的機翼上得到了使用。在數值模擬方面，Al-Khalil[4]、Morency[5]、Silva[6-7]等對氣熱與電熱防冰進行了大量研究。隨著近年來對等離子體流動控制的研究不斷增多[8]，阻擋介質放電(DBD)等離子體激勵器對流場溫度的影響已經被注意到[9-11]。這意味著等離子體也可以作為一種熱源加入流場中。目前等離子流動控制方面的研究有：流動分離控制[12-13]、激波控制[14-15]、預混空氣丙烷點火[16]、升阻控制[17]等。對等離子體激勵器進行數值研究的常用方法有：① 粒子單元直接蒙特卡羅方法[18]；② 基于集總電路的方法[19]； ③ 流場方程與泊松電場方程耦合的方法[20-21]； ④ 基于試驗與理論結果的唯象學方法[22-24]。唯象學模型不關注放電細節, 僅抓住激勵對流場的直接效應, 將激勵以動量和能量源項形式耦合到流體方程,節約計算耗時[23]。從目前公開文獻來看，等離子體防冰方面研究開展較少，文獻[25]中研究人員在圓柱表面布置阻擋介質等離子體激勵器，在結冰條件下通過試驗研究了等離子體在防/除冰方面的影響。

在本文中，首先建立基于Messinger方法的積冰模型，對于典型光冰霜冰進行驗證計算；其次，使用唯象學等離子體模型，計算了納秒脈沖阻擋介質放電(NSDBD)作用下的空氣流場，并對NSDBD等離子體激勵器的防冰特性進行了數值研究。

1 數值模擬方法及計算模型

1.1 空氣流場計算

采用非定常雷諾平均Navier-Stokes(URANS)方程[24]對流場進行求解，NSDBD等離子體激勵器對空氣流場的影響作為動量和能量源項與URANS方程耦合，忽略流場對等離子體的影響。

(1)

(2)

(3)

式中:ρ為空氣密度；t為時間;ui為速度；xi為空間坐標；μ為空氣黏度；cp為空氣比熱；p為壓力；T為溫度；τij為切應力張量；h和H分別為焓和總焓；E為總能；Pr為普朗特數；SM為動量源項，不考慮，即SM=0；SE=ERv(x,y)/τR為能量源項，其中τR為熱源作用時間，ERv(x,y)為等離子體能量密度，mJ/cm3，將由唯象學等離子體模型給出；符號“-”“″”“～”分別表示參數的雷諾平均值、波動分量和Favre平均值。空氣流場求解采用二階迎風格式、SST(Shear Stress Transport)k-ω湍流模型進行計算。

1.2 水滴軌跡與水滴收集系數

水滴軌跡的求解采用拉格朗日方法，軌跡方程為

(4)

圖1為碰撞到翼型表面的水滴軌跡示意圖，圖中dy′為水滴釋放位置的間距，ds為翼型上水滴碰撞位置的間距。Su、Sl分別為上下水滴碰撞極限位置，y′為碰撞極限位置對應的水滴釋放位置的間距。在求得水滴軌跡的情況下,局部水滴收集系數β定義為

(5)

圖1 水滴軌跡示意圖Fig.1 Sketch of droplet trajectories

1.3 積冰熱力學模型

本文的積冰熱力學模型基于Messinger模型[3,26]考慮一個控制體內的質量與能量守恒，圖2為控制體內能量與質量守恒示意圖[27]。質量與能量守恒方程分別為

(6)

(7)

在控制體內定義凍結系數F為

(8)

當F=1時表示當前控制體內水全部凍結為冰；當F=0時表示沒有任何凍結發生。

流出控制體的水量可表示為

(9)

(10)

式中：A為換熱面積；hc為對流換熱系數；Ts為表面平衡溫度；Tairsurf為與控制體對流換熱的空氣溫度。

圖2 控制體內質量與能量守恒示意圖[27]Fig.2 Sketch of mass and energy balance of control volume[27]

1.4 唯象學等離子體模型

(11)

文獻[28]中的試驗也表明能量輸入Eenergy(f,U)與激勵頻率f(Hz)有關:

(12)

同時能量輸入Eenergy(f,U,p0)也與大氣壓p0(單位為torr，1 torr≈133.3 Pa)相關,其表達式為

Eenergy(f,U,p0)=

(13)

式(12)～式(13)中：a、b、c、m、n為唯象學模型的常量，取值見表1。

最后根據能量的空間分布，能量密度分布ERv(x,y)可表示為

(14)

式中:Emission(x,y)與Emissiontotal為能量分布函數,詳見文獻[23];η為用于加熱氣體的能量占總能量的百分比。

表1 唯象學模型常量Table 1 Constants of phenomenological model

2 經典算例測試及結果分析

2.1 積冰預測

選取NACA0012翼型，其弦長為0.533 4 m，積冰時間t為360 s，圖3、圖4給出了霜冰預測結果，計算條件見表2，圖5、圖6給出了光冰預測結果，計算條件見表3。表2和表3中，AOA為迎角，MVD為平均水滴直徑，LWC為液態水含量。

霜冰的計算中選取每60 s更新一次網格進行流場的重新計算。圖3顯示了每一積冰時間步冰型的增長情況(圖中c為弦長,X′為翼型的弦線方向，Y′為與弦線垂直的方向，翼型前緣點為坐標原點，圖例clean代表未結冰的干凈翼型)。圖4中將本文計算結果與文獻計算[29]和試驗[30]結果進行了對比，發現基于Messinger方法建立的結冰模型能良好地預測霜冰冰型。

圖3 各時間步霜冰冰型Fig.3 Shape of rime ice (per time step)

圖4 霜冰冰型計算與試驗結果的對比Fig.4 Comparison of shapes of rime ice between calculated and experimental results

表2 霜冰計算條件Table 2 Calculation condition of rime ice

ParameterValueAOA/(°)4Velocity/(m·s-1)67.05Temperature/K253.69MVD/μm20LWC/(g·m-3)1.0

光冰的計算中選取每60 s更新一次網格進行流場的重新計算。圖5顯示了每一積冰時間步冰型的增長情況。圖6中將計算結果與試驗對比發現，預測的光冰冰型與試驗[30]大體趨勢一致，有光冰典型的冰角，光冰前緣冰型厚度在X′方向上基本與試驗一致。但與試驗結果在上冰角有所差距，對比其他計算結果[29-30]發現，與試驗均有一定區別。綜合霜冰與光冰的計算結果，可以認為本文積冰模型是有效的。

圖5 各時間步光冰冰型Fig.5 Shape of glaze ice (per time step)

圖6 光冰冰型計算與試驗結果的對比Fig.6 Comparison of shapes of glaze ice between calculated and experimental results

表3 光冰計算條件Table 3 Calculation condition of glaze ice

ParameterValueAOA/(°)4Velocity/(m·s-1)67.05Temperature/K265.36MVD/μm20LWC/(g·m-3)1.0

2.2 唯象學等離子模型測試

這部分將應用唯象學等離子體模型與流場耦合求解進行計算，驗證算例都為單個等離子體激勵器加載在平板上，計算等離子體激勵后流場的變化。NSDBD等離子體激勵器結構示意圖如圖7所示,激勵器由陰陽電極、阻擋介質組成，等離子體在陰陽電極之間產生。

圖7 NSDBD等離子體激勵器結構Fig.7 Structure of NSDBD plasma actuator

2.2.1 算例1

本算例中計算條件與文獻[21]中的計算條件相同,峰值電壓為14 kV，氣壓為一個標準大氣壓，空氣溫度為300 K。選取從激勵器開始工作時,與文獻[21]相同時間點t=4,8,16,25 μs的壓力云圖進行對比分析。

圖8 等離子體激勵器作用下不同時刻的壓力云圖Fig.8 Contours of pressure at different time with plasma actuator

圖8給出了等離子體激勵器作用下不同時刻的壓力云圖。圖中Δp是相對大氣壓的壓力值，從圖中可以發現在等離子體能量的注入下，流場中形成了圓弧形狀的沖擊波。其中低壓區的發展緊隨著高壓區。對比圖9文獻[21]的計算結果，發現本文計算所得波的強度和位置與之基本一致。

圖9 文獻[21]中不同時刻的壓力云圖Fig.9 Contours of pressure at different time in Ref. [21]

2.2.2 算例2

為了進一步驗證唯象學等離子體激勵器模型，同樣布置等離子體激勵器在平板上，并參照文獻[31]中的紋影試驗。計算條件與試驗相同,峰值電壓為50 kV，氣壓為一個標準大氣壓，空氣溫度為300 K。圖10為t=16 μs時刻流場的情況，左側為紋影試驗[31],右側為計算得到的無量綱密度(ρ/ρ0)云圖。

圖10 紋影試驗[31]與計算結果對比(t=16 μs) Fig.10 Comparison of shadow image between experimental[31] and calculated results (t=16 μs)

圖11 波的位置計算與試驗[31]結果對比Fig.11 Comparison of wave positions between calculated and experimental[31] results

圖12 不同時刻在X=-0.001 5 m處沿Y方向的 無量綱密度分布 Fig.12 Distribution of dimensionless density for different time at X=-0.001 5 m in Y direction

從圖10中可以發現在t=16 μs時刻，計算的圓弧形沖擊波的形狀與紋影試驗結果[31]基本一致。圖11取X=-0.001 5 m處波在Y方向的位置，與試驗結果[31]相對比。發現波的傳播速度與試驗相一致，這意味著波的強度得到了比較好的預測。圖12給出了0.5 μs≤t≤16 μs時間范圍內X=-0.001 5 m處波在Y方向上無量綱密度的分布。觀察圖像可以發現有1道壓縮波從壁面處生成，其后跟著1道比較弱的膨脹波。在0.5 μs≤t≤2 μs時間范圍內壓縮波的寬度較窄，且強度較大。隨著波的傳播，膨脹波逐漸趕上壓縮波，壓縮波強度逐漸降低。

綜上兩個等離子體的算例，可以說明唯象學模型能較好地反應出等離子體能量的注入對流場的影響。

3 NSDBD等離子體激勵器防冰特性

在積冰模型和唯象學等離子體模型的基礎上，本節將結合這兩種模型，對NSDBD等離子體激勵器防冰特性進行分析。

3.1 NSDBD等離子體激勵器分布

首先將等離子激勵器安置在NACA0012翼型前緣，內置電極、外漏電極按照流向分布，布置位置見圖13。根據唯象學模型僅抓住激勵對流場的直接效應, 將激勵以能量源項形式耦合到流體方程，按前述式(14)進行計算。由唯象學模型考慮將等離子體對空氣流場的影響耦合進能量源項，得到的無量綱能量密度分布見圖14，記為狀態A。根據圖14可以看能量密度主要分布在等離子體激勵器外漏電極與內埋電極之間，這與等離子體激勵器放電過程中產生等離子體的位置一致。根據文獻[30]中的試驗研究，其環境溫度在244.75～268.7 K之間。為驗證等離子體防冰的有效性，選擇較為極端的低溫環境(244.75 K)進行數值模擬，其他計算條件詳見表4。

圖13 NACA0012翼型上NSDBD等離子體激勵器的 分布(狀態A)Fig.13 Distribution of NSDBD plasma actuator on NACA0012 airfoil (State A)

圖14 NACA0012翼型上無量綱能量密度分布Fig.14 Distribution of dimensionless energy density on NACA0012 airfoil

表4 防冰計算條件Table 4 Calculation condition of anti-icing

ParameterValueAOA/(°)4Velocity/(m·s-1)67.05Temperature/K244.75MVD/μm20LWC/(g·m-3)1.0

在本節防冰數值計算時等離子體激勵器參數為：峰值電壓30 kV, 激勵器頻率4 kHz, 環境壓力101 300 Pa。文獻[25]對阻擋介質等離子體激勵器防冰進行了試驗研究，試驗環境溫度為263.75 K，選取的激勵器電壓為15 kV。在本文中計算環境溫度為244.75 K，因此依照此環境溫度，最后選擇更高的峰值電壓30 kV。

3.2 NSDBD等離子體開啟時的流場結果

為了分析NSDBD等離子體激勵器對流場的影響，首先計算了不加載等離子體時的流場，壓力云圖和溫度云圖見圖15(a)與圖16(a)。采用URANS方法對等離子體作用下的流場進行計算。選取的激勵頻率為4 kHz，周期則為250 μs。圖15與圖16給出了不同時刻的壓力與溫度云圖。圖15中Δp是相對于環境壓力(101 300 Pa)的壓力值。

圖15 不同時刻的壓力云圖與流線分布Fig.15 Contours of pressure and distribution of streamlines at different time

圖16 不同時刻的溫度云圖Fig.16 Contours of temperature at different time

根據圖15和圖16，可以發現等離子對翼型周圍流場的影響主要分為兩方面：① 形成沖擊波；② 加熱翼型表面空氣。對比圖15(a)～圖15(d)可以發現。等離子體引起翼型表面的壓力變化形成沖擊波，改變了周圍流動，流線有種被“沖散”的現象。隨著時間推移，圖15(c)～圖15(d)所對應的時間段內，可以發現這種沖擊效應逐漸減弱，變得不明顯，由翼型表面向遠場擴散。對比圖16(a)～圖16(d)的溫度云圖可以發現，在NSDBD等離子體激勵器作用下機翼表面被高溫氣體所覆蓋。在一個激勵周期內(250 μs)翼型表面的氣體一直保持在一個較高的溫度。圖17給出了翼型表面氣體在一個激勵周期內時間平均溫度(Tmean)與不加等離子體時的對比。圖17中橫軸為翼型表面距離駐點的位置(S)除以弦長得到的無量綱長度。在NSDBD等離子體激勵器開啟時，翼型前緣表面的時間平均空氣溫度要遠高于不開啟時，且氣流溫度在273.15 K之上。圖18給出了駐點處空氣溫度在一個激勵周期內的變化，駐點處的空氣被快速加熱到一個較高的溫度，在等離子體放電結束之后空氣與周圍發生對流后溫度逐漸降低。可見這些在翼型周圍保持較長時間的高溫氣體將會起到防冰作用。

圖17 翼型前緣時間平均空氣溫度對比Fig.17 Comparison of time-averaged air temperature around airfoil at leading edge

圖18 駐點附近空氣溫度在一個激勵周期內的變化Fig.18 Variation of air temperature around stagnation point in a pulse period

3.3 NSDBD等離子體激勵器防冰結果

圖19 凍結系數對比Fig.19 Comparison of freezing fraction

圖20 t=60 s時的冰型對比Fig.20 Comparison of ice shapes at t=60 s

在3.2節計算結果的基礎上，考慮到NSDBD等離子體激勵器對防冰在溫度方面的影響，將激勵器一周期的時間平均翼型前緣空氣溫度耦合進積冰模型的式(10)中進行計算。由于目前尚未考慮積冰表面對等離子體激勵器的影響。因此積冰模型中，僅由干凈翼型向前推進一個時間步，計算經過Δt=60 s后的情況。圖19與圖20分別給出了表4積冰條件下的凍結系數與推進一個時間步的冰型。圖19中可以看出，當不開啟等離子體激勵器時，凍結系數F基本上為1，這意味著水滴碰撞到機翼表面立即發生凍結。開啟等離子體激勵器時，凍結系數基本為0，這表示沒有發生凍結形成積冰。圖20更直觀地顯示出在開啟等離子體激勵器時，前緣防冰區沒有積冰出現。可見等離子體影響下在翼型前緣形成的高溫氣體起到了防冰作用。

3.4 NSDBD等離子體激勵器防冰特性

按照前述的方法，對等離子體激勵器參數對其防冰特性的影響進一步開展了研究。流場環境條件與3.3節中的相同，具體見表4。

3.4.1 峰值電壓對防冰效果的影響

積冰時間為60 s，等離子體激勵器的布置與3.1節中相同，激勵器頻率為4 kHz,選取不同峰值電壓(U=15，20，27，30 kV)進行計算，結果見圖21。

圖21 不同峰值電壓下的防冰特性Fig.21 Anti-icing property at different peak voltages

圖21(b)為不同峰值電壓下冰型的對比，在30 kV時，沒有積冰產生。在15～27 kV時有積冰出現，且隨著峰值電壓降低，積冰越嚴重。由圖21(a)可知，在峰值電壓為15 kV時，凍結系數有點類似于光冰，即碰撞到翼型表面的水滴沒有立即發生凍結，而是有部分向后溢流。15 kV時的冰型有點類似于光冰，有冰角出現。根據圖21(c)，可以看到15 kV時的前緣空氣平均溫度在273.15 K附近，而其他峰值電壓下，平均溫度基本都高于273.15 K。這種溫度的差異導致在15 kV時積冰較嚴重，其他所計算的峰值電壓下，雖還有積冰產生，但并不多。這種峰值電壓引起的溫度上的差異，通過分析唯象學模型式(11)可以看到，峰值電壓越高，等離子體能量輸入就越高，導致了高電壓產生較高的溫度。并且由于等離子體激勵器分布位置都與圖13中的一致，這幾個峰值電壓下的溫度圖像有著相同的趨勢。總體上電壓的大小對等離子體防冰效果有著較大影響。從能量的角度考慮，存在最優的防冰電壓取值。可以看到27 kV時尚未完全防止結冰產生，30 kV時已看不到冰層出現。在27～30 kV之間，選取更多的電壓(U=28,29 kV)進行計算，在計算過程中得到的冰型見圖22，最后可以看出最優防冰電壓區間為29～29.25 kV。

圖22 不同峰值電壓對應的冰型Fig.22 Ice shapes at different peak voltages

3.4.2 激勵器頻率對防冰效果的影響

積冰時間為60 s，等離子體激勵器的布置能量分布與3.1節相同，峰值電壓選取為30 kV,選取不同激勵頻率(f=1，2，4 kHz)進行計算，結果見圖23。

圖23 不同激勵頻率下的防冰特性Fig.23 Anti-icing property at different plus frequencies

圖23(b)給出了峰值電壓為30 kV時不同激勵頻率下，等離子體激勵器開啟時翼型前緣積冰情況。在4 kHz時前緣沒有積冰，在1 kHz與2 kHz 時仍有積冰產生。圖23(c)中更高的激勵頻率對應著較高的溫度。雖然激勵器分布位置相同，但在不同的激勵頻率下溫度曲線上沒有顯現出相同的趨勢。分析唯象學模型式(12)，激勵頻率越高，等離子體能量越大，使得翼型表面空氣溫度越高。由于圖23(c)中溫度為一個激勵周期的時間平均溫度，激勵頻率越低，則激勵周期越長。也就是說激勵頻率較低時，翼型表面的高溫氣體與周圍環境中的冷空氣有較長的熱交換時間。導致激勵頻率低時溫度曲線表現得更加平滑，如圖23(c)中f=1 kHz的曲線。在圖23(c)中也能發現高溫氣體隨著熱交換時間的增加，在翼型表面覆蓋的范圍更廣，有延翼型表面向下游發展的趨勢。綜上，在較低的激勵頻率下翼型表面的氣溫較低，有大范圍低于273.15 K的情況，使得低頻率下NSDBD等離子體激勵器防冰效果較差。根據上述計算，在其他給定條件不變的情況下，從能量的角度，大致可以看出最優激勵頻率選擇范圍應在2～4 kHz。

3.4.3 激勵器的布置方式對防冰效果的影響

若將NSDBD等離子體激勵器外漏電極與內埋電極位置交換，那么其能量密度的分布也將反向。按照圖13的等離子體激勵器位置，將其電極位置交換，得到新的電極分布如圖24所示，得到的反向能量密度分布見圖25。記這個新的能量密度分布為狀態B。積冰條件見表4，積冰時間為60 s。等離子體激勵器參數為：峰值電壓30 kV，激勵頻率4 kHz。計算結果與按圖13中的能量分布方式所得結果的對比見圖26。

在圖26中，State A曲線代表按照圖13中的能量密度分布計算所得的結果；State B曲線代表按照圖25反向分布的能量密度計算所得結果。觀察圖26(b)可以發現B狀態下有積冰產生，防冰效果沒有原狀態好。由于峰值電壓與激勵頻率以及環境氣壓都相同，根據唯象學模型，這兩種等離子體激勵器布置方式注入的總能量是相同的。觀察圖25(b)可以發現在駐點附近能量密度處于較低的水平。從圖26(c)中狀態B的溫度分布也可以看到駐點附近的溫度處于較低狀態(低于273.15 K)，加上駐點附近的熱交換較為劇烈，駐點附近發生水滴凍結的現象，從圖26(a)中也可以看到駐點附近凍結系數較大。可以發現等離子激勵器分布對防冰特性的影響需要具體分析，即使注入總能量相同的條件下，防冰特性也會不同。

圖24 NACA0012翼型上NSDBD激勵器的分布 (狀態B)Fig.24 Distribution of NSDBD plasma actuator on NACA0012 airfoil (State B)

圖25 NACA0012翼型無量綱能量密度反向分布Fig.25 Distribution of dimensionless energy density in opposite direction on NACA0012 airfoil

圖26 不同能量密度分布下的防冰特性Fig.26 Anti-icing property at different distribution of energy density

3.4.4 NSDBD等離子激勵器防冰能耗分析

根據唯象學模型，可以得到等離子體激勵器的能量密度分布，進而可以求得等離子體的功率密度。激勵頻率取4 kHz時不同峰值電壓下對應的能量密度與等離子體的功率密度見圖27。需要注意的是由于使用唯象學模型進行計算，因此這里提到的等離子體功率密度，僅為等離子體影響范圍內等離子體提供的功率，并不是等離子體激勵裝置總體的能量消耗。

圖27 不同峰值電壓下的等離子能量密度與功率密度Fig.27 Energy density and power density of plasma at different peak voltages

據3.4.1節的計算分析可知等離子體激勵器最優電壓取值區間在29～29.25 kV，對應的功率密度大小為3.24～3.31 kW/m2。根據文獻[6-7]中的電熱防冰計算，其來流速度為89.4 m/s，環境溫度為251.55 K，液態水含量為0.55 g/m3，水滴直徑為20 μm，防冰電加熱能量密度選擇在18.6～43.3 kW/m2的范圍內。雖然文獻中的條件與本文等離子體防冰計算時的有差別，但從功率密度的角度可以看到等離子功率密度要小于電加熱功率密度，從整體能耗方面考慮，這對降低防冰系統總能耗是有幫助的。未來的工作中將進一步完善計算模型并分析等離子體防冰裝置的總能耗，將其與電熱裝置總能耗的差異進行研究。

4 結 論

1) NSDBD等離子體激勵器對空氣流場的影響通過唯象學模型計算，發現其能快速加熱周圍氣體，使翼面附近溫度升高，峰值電壓對其影響較大。在本文環境條件下存在的最優防冰效果峰值電壓的取值區間為29～29.25 kV。

2) NSDBD等離子體激勵器激勵頻率對防冰效果的影響有兩方面：① 激勵頻率越大，等離子體注入能量越大，翼面氣流溫度越高，防冰效果越好；② 激勵頻率較低時，激勵周期較長，一個周期內，翼面高溫氣體與環境冷空氣熱交換時間較長，此時氣流溫度進一步降低，防冰效果不佳。

3) NSDBD等離子體激勵器布置方式對防冰效果的影響需要具體分析。即使等離子體輸入的總能量相同，在不同的布置方式下也可能會有不同的防冰效果。

4) 根據唯象學模型計算出等離子體的功率密度，將其與典型電熱防冰功率密度比較，等離子體功率密度較低，對降低防冰裝置整體能耗方面是有利的。

[1] DURST F, MILOIEVIC D, SCH?NUNG B. Eulerian and lagrangian predictions of particulate two-phase flows: A numerical study[J]. Applied Mathematical Modelling, 1984, 8(2): 101-115.

[2] 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.

[3] MACARTHUR C D. Numerical simulation of airfoil ice accretion: AIAA-1983-0112[R]. Reston, VA: AIAA, 1983.

[4] AL-KHALIL K M, HORVATH C, MILLER D R, et al. Validation of NASA thermal ice protection computer codes. Ⅲ—The validation of ANTICE: AIAA-1997-0051[R]. Reston, VA: AIAA, 1997.

[5] MORENCY F, BRAHIMI M T, TEZOK F, et al. Hot air anti-icing system modelization in the ice predict ion code CANICE: AIAA-1998-0192[R]. Reston, VA: AIAA, 1998.

[6] SILVA G A L, SILVARES O M, ZERBINI E J G J. Numerical simulation of airfoil thermal anti-ice operation: Part Ⅰ—Mathematical modeling[J]. Journal of Aircraft, 2007, 44(2): 627-634.

[7] SILVA G A L, SILVARES O M, ZERBINI E J G J. Numerical simulation of airfoil thermal anti-ice operation: Part Ⅱ—Implementation and results[J]. Journal of Aircraft, 2007, 44(2): 635-641.

[8] ROTH J R, SHERMAN D M, WILKINSON S P.Boundary layer flow control with a one atmosphere uniform glow discharge surface plasma: AIAA-1998-0328[R]. Reston, VA: AIAA, 1998.

[9] WINKEL R, CORREALE G, KOTSONIS M. Effect of dielectric material on thermal effect produced by ns-DBD plasma actuator: AIAA-2014-2119[R]. Reston, VA: AIAA, 2014.

[10] ERFANI R, HALE C, KONTIS K. The influence of electrode configuration and dielectric temperature on plasma actuator performance: AIAA-2011-0955[R]. Reston, VA: AIAA, 2011.

[11] NUDNOVA M, KINDUSHEVA S, ALEKSAHDROV N, et al. Rate of plasma thermalization of pulsed nanosecond surface dielectric barrier discharge: AIAA-2010-0465[R]. Reston, VA: AIAA, 2010.

[12] ROUPASSOV D V, NIKIPELOV A A, NUDNOVA M M, et al. Flow separation control by plasma actuator with nanosecond pulsed-periodic discharge[J]. AIAA Journal, 2009, 47(1):168-185.

[13] LITTLE J, TAKASHIMA K, NISHIHARA M, et al. Separation control with nanosecond-pulse-driven dielectric barrier discharge plasma actuators[J]. AIAA Journal, 2012, 50(2): 350-365.

[14] SHANG J S, MENART J, KIMMEL R, et al. Hypersonic inlet with plasma induced compression: AIAA-2006-0764[R]. Reston, VA: AIAA, 2006.

[15] NISHIHARA M, TAKASHIMA K, RICH J W, et al. Mach 5 bow shock control by a nanosecond pulse surface DBD: AIAA-2011-1144[R]. Reston, VA: AIAA, 2011.

[16] GALLEY D, PILLA G, LACOSTE D, et al. Plasma-enhanced combustion of a lean premixed air-propane turbulent flame using a nanosecond repetitively pulsed plasma: AIAA-2005-1193[R]. Reston, VA: AIAA, 2005.

[17] MIZOKAMI T, NOGUCHI D, FUKAGATA K. Lift and drag control using dielectric barrier discharge plasma actuators installed on the wingtips: AIAA-2013-2456[R]. Reston, VA: AIAA, 2013.

[18] FONT G I, JUNG S, ENLOE C L, et al. Simulation of the effects of force and heat produced by a plasma actuator on neutral flow evolution: AIAA-2006-0167[R]. Reston, VA: AIAA, 2006.

[19] ORLOV D M , CORKE T C, PATEL M P. Electric circuit model for aerodynamic plasma actuator: AIAA-2006-1206[R]. Reston, VA: AIAA, 2006.

[20] BOUEF J P, LAGMICH Y, CALLEGARI T, et al. Electro hydro dynamic force and acceleration in surfaces discharges: AIAA-2006-3574[R]. Reston, VA: AIAA, 2006.

[21] UNFER T, BOEUF J P. Modelling of a nanosecond surface discharge actuator[J]. Journal of Physics D: Applied Physics, 2009, 42(19): 194017-194028.

[22] SUZEN Y B, HUANG P G, JACOB J D, et al. Numerical simulations of plasma based flow control applications: AIAA-2005-4633[R]. Reston, VA: AIAA, 2005.

[23] 趙光銀， 李應紅， 梁華， 等. 納秒脈沖表面介質阻擋等離子體激勵唯象學仿真[J].物理學報, 2015, 64(1): 015101-015111.

ZHAO G Y, LI Y H, LIANG H, et al. Phenomenological modeling of nanosecond pulsed surface dielectric barrier discharge plasma actuation for flow control[J]. Acta Physica Sinica, 2015, 64(1): 015101-015111 (in Chinese).

[24] CHEN Z L, HAO L Z, ZHANG B Q. A model for nanosecond pulsed dielectric barrier discharge actuator and its investigation on the mechanisms of separation control over an airfoil[J]. Science China Technological Sciences, 2013, 56(5): 1055-1065.

[25] CAI J S, TIAN Y Q, MENG X S, et al. An experimental study of icing control using DBD plasma actuator[J]. Experiments in Fluids, 2017, 58(8): 102.

[26] TRAN P, BRAHIMI M T, PARASCHIVOIU I, et al. Ice accretion on aircraft wings with thermodynamics effects: AIAA-1994-0605[R]. Reston, VA: AIAA, 1994.

[27] 易賢, 桂業偉, 朱國林. 飛機三維結冰模型及其數值求解方法[J]. 航空學報, 2010, 31(11): 2152-2158.

YI X, GUI Y W, ZHU G L. Numerical method of a three-dimensional ice accretion model of aircraft[J]. Acta Aeronautica et Astronautica Sinica, 2010, 31(11): 2152-2158 (in Chinese).

[28] TAKASHIMA K, ZUZEEK Y, LEMPERT W R, et al. Characterization of surface dielectric barrier discharge plasma sustained by repetitive nanosecond pulses: AIAA-2010-4764[R]. Reston, VA: AIAA, 2010.

[29] FORTIN G, ILINCA A, LAFORTE J, et al. Prediction of 2D airfoil ice accretion by bisection method and by rivulets and beads modeling: AIAA-2003-1076[R]. Reston, VA: AIAA, 2003.

[30] SHIN J, BOND T H. Results of an icing test on a NACA 0012 airfoil in the NASA lewis icing research tunnel: AIAA-1992-0647[R]. Reston, VA: AIAA, 1992.

[31] STARIKOVSKII A Y, ROUPASSOV D V, NIKIPELOV A A, et al. Acoustic noise and flow separation control by plasma actuator: AIAA-2009-0695[R]. Reston, VA: AIAA, 2009.