賈韞澤，桑為民，蔡旸

西北工業(yè)大學(xué) 航空學(xué)院，西安 710072

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

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

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

1 數(shù)值模擬方法及計(jì)算模型

1.1 空氣流場計(jì)算

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

(1)

(2)

(3)

式中:ρ為空氣密度；t為時(shí)間;ui為速度；xi為空間坐標(biāo)；μ為空氣黏度；cp為空氣比熱；p為壓力；T為溫度；τij為切應(yīng)力張量；h和H分別為焓和總焓；E為總能；Pr為普朗特?cái)?shù)；SM為動量源項(xiàng)，不考慮，即SM=0；SE=ERv(x,y)/τR為能量源項(xiàng)，其中τR為熱源作用時(shí)間，ERv(x,y)為等離子體能量密度，mJ/cm3，將由唯象學(xué)等離子體模型給出；符號“-”“″”“～”分別表示參數(shù)的雷諾平均值、波動分量和Favre平均值?？諝饬鲌銮蠼獠捎枚A迎風(fēng)格式、SST(Shear Stress Transport)k-ω湍流模型進(jìn)行計(jì)算。

1.2 水滴軌跡與水滴收集系數(shù)

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

(4)

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

(5)

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

1.3 積冰熱力學(xué)模型

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

(6)

(7)

在控制體內(nèi)定義凍結(jié)系數(shù)F為

(8)

當(dāng)F=1時(shí)表示當(dāng)前控制體內(nèi)水全部凍結(jié)為冰；當(dāng)F=0時(shí)表示沒有任何凍結(jié)發(fā)生。

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

(9)

(10)

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

圖2 控制體內(nèi)質(zhì)量與能量守恒示意圖[27]Fig.2 Sketch of mass and energy balance of control volume[27]

1.4 唯象學(xué)等離子體模型

(11)

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

(12)

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

Eenergy(f,U,p0)=

(13)

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

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

(14)

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

表1 唯象學(xué)模型常量Table 1 Constants of phenomenological model

2 經(jīng)典算例測試及結(jié)果分析

2.1 積冰預(yù)測

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

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

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

圖4 霜冰冰型計(jì)算與試驗(yàn)結(jié)果的對比Fig.4 Comparison of shapes of rime ice between calculated and experimental results

表2 霜冰計(jì)算條件Table 2 Calculation condition of rime ice

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

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

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

圖6 光冰冰型計(jì)算與試驗(yàn)結(jié)果的對比Fig.6 Comparison of shapes of glaze ice between calculated and experimental results

表3 光冰計(jì)算條件Table 3 Calculation condition of glaze ice

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

2.2 唯象學(xué)等離子模型測試

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

圖7 NSDBD等離子體激勵器結(jié)構(gòu)Fig.7 Structure of NSDBD plasma actuator

2.2.1 算例1

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

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

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

圖9 文獻(xiàn)[21]中不同時(shí)刻的壓力云圖Fig.9 Contours of pressure at different time in Ref. [21]

2.2.2 算例2

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

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

圖11 波的位置計(jì)算與試驗(yàn)[31]結(jié)果對比Fig.11 Comparison of wave positions between calculated and experimental[31] results

圖12 不同時(shí)刻在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中可以發(fā)現(xiàn)在t=16 μs時(shí)刻，計(jì)算的圓弧形沖擊波的形狀與紋影試驗(yàn)結(jié)果[31]基本一致。圖11取X=-0.001 5 m處波在Y方向的位置，與試驗(yàn)結(jié)果[31]相對比。發(fā)現(xiàn)波的傳播速度與試驗(yàn)相一致，這意味著波的強(qiáng)度得到了比較好的預(yù)測。圖12給出了0.5 μs≤t≤16 μs時(shí)間范圍內(nèi)X=-0.001 5 m處波在Y方向上無量綱密度的分布。觀察圖像可以發(fā)現(xiàn)有1道壓縮波從壁面處生成，其后跟著1道比較弱的膨脹波。在0.5 μs≤t≤2 μs時(shí)間范圍內(nèi)壓縮波的寬度較窄，且強(qiáng)度較大。隨著波的傳播，膨脹波逐漸趕上壓縮波，壓縮波強(qiáng)度逐漸降低。

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

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

在積冰模型和唯象學(xué)等離子體模型的基礎(chǔ)上，本節(jié)將結(jié)合這兩種模型，對NSDBD等離子體激勵器防冰特性進(jìn)行分析。

3.1 NSDBD等離子體激勵器分布

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

圖13 NACA0012翼型上NSDBD等離子體激勵器的 分布(狀態(tài)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 防冰計(jì)算條件Table 4 Calculation condition of anti-icing

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

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

3.2 NSDBD等離子體開啟時(shí)的流場結(jié)果

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

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

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

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

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

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

3.3 NSDBD等離子體激勵器防冰結(jié)果

圖19 凍結(jié)系數(shù)對比Fig.19 Comparison of freezing fraction

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖24 NACA0012翼型上NSDBD激勵器的分布 (狀態(tài)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等離子激勵器防冰能耗分析

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

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

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

4 結(jié) 論

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

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

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

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

[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] 趙光銀， 李應(yīng)紅， 梁華， 等. 納秒脈沖表面介質(zhì)阻擋等離子體激勵唯象學(xué)仿真[J].物理學(xué)報(bào), 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] 易賢, 桂業(yè)偉, 朱國林. 飛機(jī)三維結(jié)冰模型及其數(shù)值求解方法[J]. 航空學(xué)報(bào), 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.