賈韞澤，桑為民，蔡旸

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

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

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

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

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

1.1 空氣流場(chǎng)計(jì)算

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

(1)

(2)

(3)

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

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

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

(4)

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

(5)

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

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

本文的積冰熱力學(xué)模型基于Messinger模型[3,26]考慮一個(gè)控制體內(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)水全部?jī)鼋Y(jié)為冰；當(dāng)F=0時(shí)表示沒(méi)有任何凍結(jié)發(fā)生。

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

(9)

(10)

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

圖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)與激勵(lì)頻率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)典算例測(cè)試及結(jié)果分析

2.1 積冰預(yù)測(cè)

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

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

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

圖4 霜冰冰型計(jì)算與試驗(yàn)結(jié)果的對(duì)比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)行流場(chǎng)的重新計(jì)算。圖5顯示了每一積冰時(shí)間步冰型的增長(zhǎng)情況。圖6中將計(jì)算結(jié)果與試驗(yàn)對(duì)比發(fā)現(xiàn)，預(yù)測(cè)的光冰冰型與試驗(yàn)[30]大體趨勢(shì)一致，有光冰典型的冰角，光冰前緣冰型厚度在X′方向上基本與試驗(yàn)一致。但與試驗(yàn)結(jié)果在上冰角有所差距，對(duì)比其他計(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é)果的對(duì)比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é)等離子模型測(cè)試

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

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

2.2.1 算例1

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

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

圖8給出了等離子體激勵(lì)器作用下不同時(shí)刻的壓力云圖。圖中Δp是相對(duì)大氣壓的壓力值，從圖中可以發(fā)現(xiàn)在等離子體能量的注入下，流場(chǎng)中形成了圓弧形狀的沖擊波。其中低壓區(qū)的發(fā)展緊隨著高壓區(qū)。對(duì)比圖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é)等離子體激勵(lì)器模型，同樣布置等離子體激勵(lì)器在平板上，并參照文獻(xiàn)[31]中的紋影試驗(yàn)。計(jì)算條件與試驗(yàn)相同,峰值電壓為50 kV，氣壓為一個(gè)標(biāo)準(zhǔn)大氣壓，空氣溫度為300 K。圖10為t=16 μs時(shí)刻流場(chǎng)的情況，左側(cè)為紋影試驗(yàn)[31],右側(cè)為計(jì)算得到的無(wú)量綱密度(ρ/ρ0)云圖。

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

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

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

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

3 NSDBD等離子體激勵(lì)器防冰特性

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

3.1 NSDBD等離子體激勵(lì)器分布

首先將等離子激勵(lì)器安置在NACA0012翼型前緣，內(nèi)置電極、外漏電極按照流向分布，布置位置見圖13。根據(jù)唯象學(xué)模型僅抓住激勵(lì)對(duì)流場(chǎng)的直接效應(yīng), 將激勵(lì)以能量源項(xiàng)形式耦合到流體方程，按前述式(14)進(jìn)行計(jì)算。由唯象學(xué)模型考慮將等離子體對(duì)空氣流場(chǎng)的影響耦合進(jìn)能量源項(xiàng)，得到的無(wú)量綱能量密度分布見圖14，記為狀態(tài)A。根據(jù)圖14可以看能量密度主要分布在等離子體激勵(lì)器外漏電極與內(nèi)埋電極之間，這與等離子體激勵(lì)器放電過(guò)程中產(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等離子體激勵(lì)器的 分布(狀態(tài)A)Fig.13 Distribution of NSDBD plasma actuator on NACA0012 airfoil (State A)

圖14 NACA0012翼型上無(wú)量綱能量密度分布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í)等離子體激勵(lì)器參數(shù)為：峰值電壓30 kV, 激勵(lì)器頻率4 kHz, 環(huán)境壓力101 300 Pa。文獻(xiàn)[25]對(duì)阻擋介質(zhì)等離子體激勵(lì)器防冰進(jìn)行了試驗(yàn)研究，試驗(yàn)環(huán)境溫度為263.75 K，選取的激勵(lì)器電壓為15 kV。在本文中計(jì)算環(huán)境溫度為244.75 K，因此依照此環(huán)境溫度，最后選擇更高的峰值電壓30 kV。

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

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

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

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

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

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

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

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

3.4 NSDBD等離子體激勵(lì)器防冰特性

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

3.4.1 峰值電壓對(duì)防冰效果的影響

積冰時(shí)間為60 s，等離子體激勵(lì)器的布置與3.1節(jié)中相同，激勵(lì)器頻率為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)為不同峰值電壓下冰型的對(duì)比，在30 kV時(shí)，沒(méi)有積冰產(chǎn)生。在15～27 kV時(shí)有積冰出現(xiàn)，且隨著峰值電壓降低，積冰越嚴(yán)重。由圖21(a)可知，在峰值電壓為15 kV時(shí)，凍結(jié)系數(shù)有點(diǎn)類似于光冰，即碰撞到翼型表面的水滴沒(méi)有立即發(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)生，但并不多。這種峰值電壓引起的溫度上的差異，通過(guò)分析唯象學(xué)模型式(11)可以看到，峰值電壓越高，等離子體能量輸入就越高，導(dǎo)致了高電壓產(chǎn)生較高的溫度。并且由于等離子體激勵(lì)器分布位置都與圖13中的一致，這幾個(gè)峰值電壓下的溫度圖像有著相同的趨勢(shì)?？傮w上電壓的大小對(duì)等離子體防冰效果有著較大影響。從能量的角度考慮，存在最優(yōu)的防冰電壓取值。可以看到27 kV時(shí)尚未完全防止結(jié)冰產(chǎn)生，30 kV時(shí)已看不到冰層出現(xiàn)。在27～30 kV之間，選取更多的電壓(U=28,29 kV)進(jìn)行計(jì)算，在計(jì)算過(guò)程中得到的冰型見圖22，最后可以看出最優(yōu)防冰電壓區(qū)間為29～29.25 kV。

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

3.4.2 激勵(lì)器頻率對(duì)防冰效果的影響

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

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

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

3.4.3 激勵(lì)器的布置方式對(duì)防冰效果的影響

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

在圖26中，State A曲線代表按照?qǐng)D13中的能量密度分布計(jì)算所得的結(jié)果；State B曲線代表按照?qǐng)D25反向分布的能量密度計(jì)算所得結(jié)果。觀察圖26(b)可以發(fā)現(xiàn)B狀態(tài)下有積冰產(chǎn)生，防冰效果沒(méi)有原狀態(tài)好。由于峰值電壓與激勵(lì)頻率以及環(huán)境氣壓都相同，根據(jù)唯象學(xué)模型，這兩種等離子體激勵(lì)器布置方式注入的總能量是相同的。觀察圖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ù)較大。可以發(fā)現(xiàn)等離子激勵(lì)器分布對(duì)防冰特性的影響需要具體分析，即使注入總能量相同的條件下，防冰特性也會(huì)不同。

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

圖25 NACA0012翼型無(wú)量綱能量密度反向分布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等離子激勵(lì)器防冰能耗分析

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

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

據(jù)3.4.1節(jié)的計(jì)算分析可知等離子體激勵(lì)器最優(yōu)電壓取值區(qū)間在29～29.25 kV，對(duì)應(yīng)的功率密度大小為3.24～3.31 kW/m2。根據(jù)文獻(xiàn)[6-7]中的電熱防冰計(jì)算，其來(lái)流速度為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í)的有差別，但從功率密度的角度可以看到等離子功率密度要小于電加熱功率密度，從整體能耗方面考慮，這對(duì)降低防冰系統(tǒng)總能耗是有幫助的。未來(lái)的工作中將進(jìn)一步完善計(jì)算模型并分析等離子體防冰裝置的總能耗，將其與電熱裝置總能耗的差異進(jìn)行研究。

4 結(jié) 論

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

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

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

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

[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ì)阻擋等離子體激勵(lì)唯象學(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è)偉, 朱國(guó)林. 飛機(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.