非平衡感應耦合等離子體流場與電磁場作用機理的數值模擬*

喻明浩

(西安理工大學機械與精密儀器工程學院,西安 710048)

1 引 言

再入飛行器完成太空任務重返地球大氣層時,由于其飛行速度很快,在飛行器前方會形成很強的弓形激波,激波層內氣體被急劇壓縮和加熱形成溫度高達幾千攝氏度的等離子體氣流,飛行器頭部材料將被等離子體氣流急劇加熱而發(fā)生分解反應,因此,若不在飛行器頭部表面加載熱防護材料,飛行器很可能被等離子體氣流損壞或燒毀.

為了開發(fā)耐高溫且質量輕的飛行器熱防護材料,近年來,很多國家建立了高溫高焓風洞以開展飛行器防熱材料燒蝕試驗,如電弧加熱風洞、激波風洞、感應耦合等離子體(inductively coupled plasma,ICP)風洞.由于ICP風洞能夠產生連續(xù)、高純度、接近真實飛行環(huán)境的高溫、高焓氣流,因此它被廣泛用于開發(fā)耐燒蝕的熱防護材料[1]、研究飛行器表面氮化反應機理[2]、研制航天器新型薄膜減速傘材料[3]等.ICP除了用于航天領域外,也可以用于研究甲醛氣體協(xié)同凈化處理[4]、超微金屬粉末球化處理[5]、微納米顆粒制備[6]等.在進行上述工業(yè)應用研究時,如果想對這些實驗研究進行準確、細致地分析,就必須先準確掌握ICP的流動特性與形成機理.因此,為了準確獲得ICP的流動特性與形成機理,我們需要對ICP的形成過程、高頻放電特性以及電磁場與流場的作用機理等方面進行深入、細致的研究.

ICP風洞的系統(tǒng)結構通常由三部分組成: 進氣部分、高頻放電部分和試驗部分,這三部分的結構布局和組成單元如圖1所示.ICP風洞的核心工作部分由高頻電源、感應線圈和石英管組成,由于等離子體流動率先形成于石英管內,因此石英管通常又被稱為等離子體炬.ICP的形成過程是一個復雜的物理、化學過程,它涉及高頻放電、電磁加熱、趨膚效應、非平衡內能交換、氣體電離/離解等過程.ICP的具體形成過程如下: 常溫氣體通過進氣道(如圖1右上角所示)進入到石英管前端進行充分預混合[7],當氣體流經感應線圈下方時,經電火花點火,氣體介質將在感應線圈產生的高頻交變電磁場作用下發(fā)生感應放電,此時由于趨膚效應和焦耳熱效應的作用,一部分電能將轉化為分子熱能.隨著氣體溫度不斷升高,氣體分子將發(fā)生電離和離解反應并釋放出大量的熱,進而形成由電子、原子、分子和陽離子組成的準電中性等離子體流動.由于等離子體中存在大量的自由電子和活性陽離子,它們在交變電磁場的作用下會產生感應磁場.由于該感應磁場與感應線圈產生的外加磁場存在一定的相位差和耦合關系,因此所生成的等離子體常被稱為ICP.由于ICP的形成過程較為復雜,且流動溫度也較高,使用實驗測量工具對其流動參數進行全面診斷難度很大,因此開展數值模擬研究成為當前研究ICP的一種重要方法.

自第一臺ICP風洞問世以來[8],世界各國學者紛紛開展了相關的實驗與數值模擬研究工作[9?18].1976年,Barnes和Nikdel[19]使用一維電磁場模型和能量平衡方程聯(lián)立求解的方法研究了氮氣ICP的溫度場和速度場的變化規(guī)律.Mostaghimi和Boulos[20]基于麥克斯韋方程組和畢奧-薩伐爾定律提出了一種二維磁矢勢模型,并成功地將該模型應用到了氬氣ICP的數值模擬中.Punjabi等[21,22]利用FLUENT軟件對不同功率ICP炬內等離子體阻抗、吸收功率、耦合效率等進行了仿真分析,研究發(fā)現: 隨著ICP風洞輸入功率的增大,其等離子體火焰體積隨之增大,等離子體最高電子溫度、最大速度和最大電子數密度均隨之增大[17,21]; 對于不同工作介質(如氮氣、氧氣、氬氣和空氣) ICP流動,在相同的工況條件下(功率、頻率、氣壓、體積流量均相同),氧氣ICP的等離子體火焰體積最小,氬氣、氮氣和空氣ICP的等離子體阻抗隨輸入功率增大而增大,但氧氣ICP的阻抗隨功率增大而呈下降趨勢[21]; 此外,無論以上何種氣體介質,ICP流動的最大溫度、速度和電子數密度均隨電源放電頻率的增大而減小[23].但以上研究均是在局域熱化學平衡假設條件下開展的,只有當工作氣壓較高時,該平衡假設才被認為是有效的.因此,當工作氣壓較低時,將熱化學非平衡模型引入數值模擬中是非常有必要的.對于非平衡ICP研究,Lei等[10]利用COMSOL軟件對低壓ICP流動進行了非平衡仿真,Stewart等[24]使用二溫度模型研究了氬氣ICP的熱非平衡放電過程,Munafò等[25]、Zhang等[26]和Lei等[10]也構建了熱非平衡和化學非平衡模型,并對氬氣ICP流動的非平衡特性進行了研究.

基于氬氣熱化學非平衡模型的發(fā)展,Degrez等[27]采用空氣化學反應動力學模型研究了空氣電離、離解反應模型對流動速度、溫度分布的影響,研究發(fā)現,由于氣體粒子的擴散效應,高頻放電區(qū)氧原子和氮原子出現了分層現象.然而,該研究沒有考慮熱力學非平衡對其仿真結果的影響.Sumi等[28]也采用雙溫度模型研究了ICP的熱非平衡特性以及溫度分布規(guī)律.Morsli和Proulx[29]基于Dunn-Kang的32化學反應模型和雙溫度模型構建了空氣ICP熱化學非平衡磁流體力學模型,并對超聲速空氣ICP的流動規(guī)律進行了研究,研究表明空氣化學反應模型對準確模擬空氣粒子的組分濃度起著至關重要的作用.此外,本課題組之前通過對電磁場方程組進行簡化、對電子輸運系數計算方法進行完善、對四溫度模型和粒子內能交換模型進行補充等,也對非平衡ICP的仿真研究做了些許工作[11,30?32].在前期研究中,本課題組提出了一種收斂速度快、相對簡單的焦耳熱源模型,該模型可以替代感應線圈區(qū)電磁場方程組的求解,并且能較準確地模擬ICP風洞試驗腔內等離子體氣流的溫度分布[11]; 在此基礎上,我們還通過耦合求解遠場電磁場和非平衡流體力學方程組研究了10 kW級氮氣與空氣ICP的流場差異及其熱化學非平衡特性[17],并對10 kW ICP風洞內氮氣離子體的超聲速流動進行了瞬態(tài)模擬[30].除此之外,我們還對空氣ICP的熱化學非平衡特性與氣壓的非線性關系進行了研究,研究發(fā)現當氣壓大于19 kPa時,10 kW空氣ICP流動趨于局域熱力學平衡態(tài)[31],并發(fā)現當電導率的計算精度從一階提高到三階時,仿真得到的氣體最高溫度會下降680 K,所以最終整理提出了一種計算空氣和氮氣等離子體三階精度電子輸運系數的方法,該方法在ICP仿真研究中具有很好的應用潛力[32].

圖1 ICP風洞系統(tǒng)結構布局Fig.1.Schematic diagram of the ICP wind tunnel system.

盡管上述研究工作對ICP數值模擬技術都起到了積極的促進作用,然而目前尚缺少對高功率非平衡ICP電磁場與流場作用機理的深入研究,關于洛倫茲力、焦耳加熱率、電子數密度、電子溫度等參數之間的相互影響關系目前尚不清晰.此外,我們前期研究主要是針對10 kW級中小功率ICP風洞展開的,對于100 kW級高功率ICP風洞來說,由于其輸入功率較大,其電磁場強度、氣體溫度和電子數密度必然將更高,因此其電磁場、流場與熱力場之間的耦合也將變得更加緊密、更加復雜.而且,由于100 kW級ICP風洞的結構參數、輸入功率、放電頻率與10 kW級ICP風洞存在諸多差異,因此本文擬對高功率ICP風洞的流場與電磁場相互作用機理進行研究,以揭示高功率ICP的高頻放電規(guī)律和非平衡流動特性.

本文以100 kW級ICP風洞為研究對象,基于流場-電磁場-熱力場-化學場-湍流場多場耦合求解技術對ICP電磁場與流場的相互作用機理進行研究,以揭示電磁場強度對流場參數的影響機理.為了準確描述高溫空氣微觀粒子的熱運動現象,本文采用含有最新電子振動弛豫時間算法的四溫模型,并使用三階精度電導率和電子導熱率進行仿真計算.由于高頻感應放電在線圈區(qū)發(fā)揮著重要作用,因此通過求解麥克斯韋方程組來計算線圈電流產生的外加電場和ICP產生的感應磁場,并使用11組分空氣、32種化學反應來模擬高溫空氣粒子的電離、離解與電量交換等化學反應過程[33,34].在進行空氣ICP仿真時,Park等[35]所提出的化學反應動力學模型常被用來模擬空氣的電離和離解過程,但Park等并未考慮N2,O2和NO分子的副電離反應由于空氣在溫度超過4000 K和9000 K時將很容易發(fā)生離解和電離反應,故為了準確模擬高溫空氣的化學反應過程,本文在數值模擬中將考慮這些副電離反應.最終,基于上述模型與二維黏性可壓縮Navier-Stokes方程組的耦合求解來對ICP的形成過程與感應放電機理進行描述.

2 控制方程與數值方法

本文所研究的100 kW級ICP風洞幾何形狀與計算網格如圖2所示.圖2(a)為遠場電磁場和ICP炬流場的計算網格,電磁場的計算域為–120 mm ≤x≤ 360 mm,0 ≤y≤ 206 mm,由101×66個網格節(jié)點組成,其中ICP炬(ICP torch)的流場網格由101×38個網格節(jié)點組成,為了準確計算線圈附近的電磁場與進氣口處氣流速度,計算網格在感應線圈區(qū)和進氣口處進行了加密處理.圖2(b)為ICP炬的結構示意圖,感應線圈圍繞ICP炬管外壁面纏繞三圈,線圈的內徑和外徑分別為47 mm和55 mm,起始線圈的中心坐標為(x,y)=(51,51) mm,相鄰感應線圈的中心線間隔為16.5 mm.

ICP風洞正常工作時,感應線圈上通有高頻交變電流,電流頻率為1.78 MHz,工作氣體如空氣或氮氣從壁面附近的狹窄開口注入,經電火花高能點火后,工作氣體開始發(fā)生離解和電離反應并產生自由電子.此后,帶電粒子在交變電流產生的交變電磁場作用下進一步發(fā)生劇烈的碰撞電離、副電離等反應,并釋放出大量能量而形成溫度高達10000 K的ICP氣流,最終該等離子體氣流從ICP炬出口流入大空間試驗腔內進行飛行器氣動性能測試、熱防護材料燒蝕試驗等實驗研究.

圖2 ICP炬計算網格和幾何結構 (a) 電磁場與流場計算網格; (b) 等離子體炬幾何結構Fig.2.Computational mesh and geometry of the inductively coupled plasma torch: (a) Computational mesh of electromagneticand flow-field; (b) geometric construction of the ICP torch.

本文所開展的仿真算例工況條件如下: 面板輸入功率P=90 kW,質量流率m=1.8 g/s,工作氣壓p=10 kPa,電流頻率f=1.76 MHz,工作氣體為空氣.對于高功率ICP風洞,由于存在較大的熱能損失(比如,冷卻水引起的熱損失,電路上的電阻抗引起的能量損失以及輻射損失),沉積到等離子體流動中的凈輸入功率通常與從供電系統(tǒng)讀取的面板輸入功率有所不同[36],Fujita等[28,37]對此進行了實驗與仿真研究,研究結果表明: 對于110 kW ICP風洞其熱效率通常為1/3.因此,本研究ICP風洞的熱效率h也設定為1/3.被等離子體實際吸收的沉積功率等于熱效率h與面板輸入功率的乘積.仿真計算時,為了確保能量守恒,每一步流場迭代均使所有流體微元的焦耳加熱率體積積分之和等于等離子體的沉積功率

2.1 流場方程組

為了簡化計算,本研究假設等離子體是電中性和二維軸對稱的; 流體微元處于熱化學非平衡態(tài)狀態(tài),從微觀上氣體溫度可分為平動溫度Ttr、轉動溫度Trot、振動溫度Tvib和電子溫度Te.仿真計算時,考慮洛倫茲力、湍流以及等離子體產生的感應磁場.最終所建立的非平衡磁流體力學方程組系統(tǒng)由總質量、動量、總能量守恒方程,11組分粒子質量守恒方程,分子振動、轉動和電子能量守恒方程,電磁場方程、湍流方程和32種空氣化學反應動力學模型組成.該方程組系統(tǒng)的矢量形式如下:

(1)式中,守恒向量Q、通量矢量F、Fυ和源項矢量W的表達式分別如下:

高溫氣體的壓強p由下式計算:

氣體的總內能E定義為

氣體粒子的平動能(Etr)、轉動能(Erot)、振動能(Evib)和電子內能(Ee)分別由下式計算:

(3)—(10)式中的氣體輸運系數(如氣體黏性μ、平動導熱系數λtr、振動導熱系數λvib和轉動導熱系數λrot)均由Yos公式[38,39]進行計算; 擴散系數Ds由Curtiss和Hirschfelder[40]給出的公式計算; 氣體電導率s和電子導熱系數λe由三階Sonine多項式和時新的電子-重粒子碰撞截面數據進行計算[32,41,42].三階精度電導率s和電子導熱系數λe的計算公式如下:

2.2 電磁場方程組

ICP的電磁場分布通?？赏ㄟ^求解以下磁矢勢方程得到:

其中A為磁矢勢和i分別表示線圈電流的角頻率、電流密度、真空介電常數和復數單位為線圈電流驅動頻率.此外,由于磁矢勢A與磁場強度H和電場強度E有如下關系[20]:

因此可通過求解磁矢勢來得到電場和磁場分布,進而得到等離子體的焦耳加熱率和洛倫茲力.

對于圓柱形ICP炬,線圈電流可以被認為是由若干平行環(huán)狀電流微元組成,因此磁矢勢可以被認為只有切向分量[20],即考慮到線圈電流產生的電磁場與等離子體產生的電磁場之間存在相位差,故切向矢量勢可表示為復數,即最終,(13)式可寫成如下形式:

式中,電流密度Jc可由計算,其中Rc為線圈半徑,I為線圈電流.待求解出磁矢勢方程后,ICP的焦耳加熱率SJoule可由下式計算得到:

軸向和徑向洛倫茲力FLx和FLy可分別由下式進行計算:

式中,s為等離子體的電導率,焦耳加熱率和洛倫茲力將被作為能量守恒和動量守恒方程的源項參與到流場方程組的迭代求解中.ICP風洞的電磁場與流場將通過焦耳加熱率和洛倫茲力進行耦合聯(lián)接.

2.3 空氣化學反應模型

為了準確模擬空氣在高溫條件下發(fā)生的離解、電離、電量交換、分子復合等化學反應,本文將文獻[33,34]給出的由11種組分(N2,O2,NO,N2+,O2+,NO+,N,O,N+,O+,e–)和32種化學反應組成的動力學模型應用到本研究的數值計算中.該32種化學反應如表1所列,化學反應r的正向反應速率kf,r由Arrhenius公式進行計算[44]:

式中的化學反應系數Cr,n和θr取自文獻[33,35].

逆向化學反應速率kb,r由平衡常數Keq進行計算:

平衡常數由曲線擬合公式進行計算[35].最終,任一空氣組分s因發(fā)生化學反應而產生的質量變化率凈質量生成率可由下式計算:

式中,nsb,r和nsf,r是組分s在進行第r個化學反應前后的化學計量系數,Ms代表物質s的摩爾質量,下標s和j代表物質組分編號.

表1 空氣化學反應模型Table 1.Chemical reaction model of air.

2.4 內部能量交換模型

由于電子、分子和原子之間存在彈性或非彈性碰撞,碰撞時會產生平動能、轉動能、振動能或電子能量的交換,為了準確模擬這些內能交換過程,根據相應的數學模型計算粒子內能交換率,并將計算得到的轉動能交換率Sint,rot、振動能交換率Sint,vib和電子能量交換率Sint,e作為能量源項添加到(2)式中的轉動、振動和電子能量守恒方程中.轉動、振動和電子能量交換率Sint,rot,Sint,vib,Sint,e的具體計算方法如下:

式中,平動能-轉動能交換率QT-R由Park[45]給出的方法計算; 轉動能-振動能量交換率QR-V由Millikan和White[46]給出的方法計算; 轉動能-電子能交換率QR-e由Lazdinis和Petrie[47]給出的方法計算; 平動能-振動能量交換率QT-V由Millikan和White[46]和Park[48]給出的公式計算; 電子能-振動能交換率Qe-V由 Landau-Teller 方程計算[49],其中氮分子和電子的振動能-電子能馳豫時間根據Kim等[50]及Bourdon和Vervisch[51]給出的方法進行計算; 平動能-電子能交換率QT-e由Appleton和Bray[52]給出的公式計算; 空氣分子因離解和電離反應產生的能量損失由 Gnoffo等[53]給出的方法分別進行計算.

2.5 湍流運輸方程

本文采用低雷諾數湍流模型—Abe-Kondoh-Naganok-e模型[54]來考慮湍流對ICP流動的影響,包含湍流動能k和耗散率e的湍流輸運方程如下:

湍流黏度由下式計算:

上述方程中的模型常數如下:

模型函數表示為

式中,y?=(νε)1/4ywd/ν,Rt=k2/(νε) ,n表示分子運動黏度,ywd為與內壁面的距離.

2.6 邊界條件和數值求解方法

2.6.1 邊界條件

對于流場方程組: 1)在等離子體炬入口處,工作氣體由距壁面約2.4 mm的開口注入,氣體質量流率和初始溫度T0=300 K作為已知輸入參數.2)在ICP炬出口處,工作氣壓設為定值p=10 kPa,其他參數由相鄰的內部網格點線性插值得到.3)在等離子體炬壁面處,認為壁面上無滑移、無催化效應、無壓力梯度,壁面溫度由導熱方程?Tw/?n=qjmax進行計算,其中qjmax為內壁面處的熱流量密度.4)在中心軸上,采用軸對稱邊界條件,即Qi,j=Qi,j+1,Q代表守恒變量.對于磁矢勢方程組,電磁場計算區(qū)域的外邊界(圖2(a)中的x=–120 mm,x=360 mm和y=206 mm)設置的離感應線圈足夠遠以使磁矢勢AR和AI在這些外邊界上數值為零,由于100 kW級高功率ICP風洞的電磁場強度較大,其影響范圍也更廣,故與10 kW ICP風洞的磁場外邊界(電磁場趨近于0的邊界)相比,本文的遠場電磁外邊界向外側進行了推移.中心軸線上x=0處我們采用軸對稱邊界條件:Ai,0=Ai,1.

2.6.2 數值求解方法

對于流場和湍流方程組,本文采用有限體積法對其進行離散化,通過時新的低耗散迎風差分法計算其對流項,并通過MUSCL (monotonic upstream-centered scheme for conversation laws)格式將其計算精度提高到二階精度; 黏性項采用二階精度中心差分法進行計算.對于非定常物理流動,為獲得其定常流場,除了要在空間上進行推進求解外,也需要在時間方向上進行推進求解,故本文采用點隱式LU-SGS (lower and upper symmetric Gauss-Seidel)時間推進算法對流場方程組進行求解,采用廣義最小殘差法(generalized minimal residualal method,GMRES)對湍流輸運方程組進行求解.

對于電磁場方程組,本文采用有限差分法對其進行離散化,采用二階中心差分法計算其數值通量,采用亞松弛迭代法快速、穩(wěn)定地求解離散后的方程,松弛因子設定為0.2,當磁矢勢第n步和第n+1步數值解的相對誤差小于10–5時,認為計算收斂.對于化學反應方程組,本文采用有限體積法對方程組進行離散化,利用托馬斯追趕法進行迭代求解.

3 結果與討論

下面將通過對100 kW級空氣ICP風洞在典型工況下(P=90 kW,m=1.8 g/s,p=10 kPa,f=1.76 MHz)的速度、電子溫度、壓強、洛倫茲力、焦耳加熱率、電場強度、電子數密度、粒子摩爾分數、平動溫度的分布規(guī)律進行分析與總結,明確ICP的流場與電磁場分布特性,并揭示其流場與電磁場的相互作用機理.

ICP炬中氣流的速度矢量、流線和電子溫度分布云圖如圖3所示,由于感應線圈區(qū)存在很強的外加電磁能,從入口流入的常溫新鮮空氣在此區(qū)域經點火后迅速發(fā)生放電反應并吸收大量的焦耳熱,因此氣體速度和溫度在此區(qū)域都迅速升高.由于空氣中的氮分子與氧分子在此區(qū)域被迅速電離和離解并釋放大量的熱,因此由電子、陽離子和原子組成的ICP溫度在感應線圈區(qū)很高,最大電子溫度約10500 K.此外,從圖3中的流線分布還可以發(fā)現,在進氣口附近形成了較強的氣流回旋現象.該回流的形成與感應線圈區(qū)的負壓強梯度和電磁加熱現象有很大關系.

圖4為等離子體流動的流線與壓力分布云圖.可見,最大氣壓處在第二和第三圈線圈之間,即在第二個線圈上游存在負壓力梯度,這是線圈上游靠近入口處產生渦流的原因之一.此外,在動量方程中,速度和壓力的分布都受到洛倫茲力的影響,而且在能量守恒方程中焦耳加熱率對流動速度也有影響.因此,線圈下方產生的大面積渦流是由感應線圈區(qū)域的負壓力梯度、洛倫茲力和焦耳加熱現象共同作用產生的.

圖3 等離子體炬內氣體流線和速度矢量(上半部分)以及電子溫度云圖(下半部分)分布Fig.3.Distributions of streamlines and velocity vector (upper),and electron temperature (lower) in the torch.

圖4 等離子體流線(上半部分)和壓力云圖(下半部分)分布Fig.4.Distributions of streamlines (upper) and pressure contour (lower).

圖5給出了軸向和徑向洛倫茲力分布.洛倫茲力作為能量守恒方程的慣性力源項參與到整個流場的迭代求解中.由圖5可見: 軸向洛倫茲力在感應線圈的第一與第二圈之間為正值,即其方向向右.然而在感應線圈第二圈與第三圈之間變?yōu)樨撝?即軸向洛倫茲力的方向反向.而對于徑向洛倫茲力(圖5下半部分),其值始終為負值,即其方向始終由壁面指向中心軸線,這表明因趨膚效應在炬壁附近生成的自由電子將始終受到指向等離子體炬中軸線的慣性力,且徑向洛倫茲力的最大值比軸向洛倫茲的最大值高3.95倍,說明等離子體的徑向動量傳遞比軸向動量傳遞強烈的多.

圖5 軸向洛倫茲力(上半部分)和徑向洛倫茲力(下半部分)的分布Fig.5.Distributions of axial Lorentz force (upper) and radial Lorentz force (lower).

圖6給出了徑向洛倫茲力(上半部分)和焦耳加熱率(下半部分)分布.可以看出,焦耳加熱率的分布形狀和位置與徑向洛倫茲力的非常相似.這表明,對于空氣ICP流動,進行焦耳加熱現象的位置主要由徑向洛倫茲力控制.此外,徑向洛倫茲力還對等離子體最大溫度和速度的位置有一定影響.另一方面,由于感應線圈上通有高頻交變電流,線圈區(qū)域會發(fā)生氣體放電和焦耳加熱現象,由于焦耳加熱率由氣體電導率和電場E共同控制,且等離子體炬壁上始終通有冷卻水,壁面溫度被限制在1000 K以下,導致壁面附近的流動溫度較低,即導電率很小,因此,盡管電場強度E在等離子體炬內壁面上很大,但最大焦耳加熱率并不出現在內壁面上,而是出現在第二個感應線圈下方距離內壁面3.5 mm處.

圖6 徑向洛倫茲力(上半部分)和焦耳加熱率(下半部分)的分布Fig.6.Distributions of Joule heating rate(lower) and radial Lorentz force (upper).

圖7給出了等離子體炬內電場強度虛部EI(上半部分)和實部ER(下半部分)的分布云圖.可知,電場強度虛部的最大值比電場強度實部大2.9倍,所以由磁矢勢計算得到的電場強度虛部EI是總電場強度的主要部分.電場強度實部ER在等離子體炬中始終為負值; 而虛部EI在靠近壁面處為負值,在第二個線圈下方靠近軸線位置處變?yōu)檎?負的電場虛部主要是由外加的高頻電流產生,而正電場則是由等離子體內部電子和陽離子產生的感應電場.

圖7 電場強度分布(虛部EI (上半部分)實部ER(下半部分))Fig.7.Distribution of electric-field intensity (imaginary part EI (upper) and real part ER (lower)).

圖8給出了等離子體炬中電場強度EI(上半部分)和電子數密度ne(下半部分)分布.在強電場的作用下,電子在徑向洛倫茲力作用下向中軸線附近移動,在第二個感應線圈下方距離等離子體炬內壁面5.5 mm處電子數密度達到最大值.大量的自由電子聚集在感應線圈區(qū),電子數密度ne=1.0×1020m–3所包絡的區(qū)域與正的電場強度EI所處位置近似相同,這表明正的電場強度EI由這些自由電子產生,它是這些自由電子在等離子體炬中自由流動時形成的感應電流產生的感應電場.

圖9為軸向位置x=68 mm處11組分空氣的摩爾分數分布.如圖9所示,氧分子被完全解離和電離成氧原子和離子,91%的氮分子在軸線附近被離解成N原子,氮原子N和氧原子O在徑向位置y≤ 30 mm之前為最主要的化學組分.因電離反應生成的自由電子和離子在本算例中仍然較少,電子e–和氧離子O+的摩爾分數約為10–3,二者電荷數趨于一致,說明等離子體呈準電中性.沿等離子體炬半徑方向,隨著平動溫度的降低,氮分子N2和NO分子的摩爾分數逐漸增大; 由于等離子體炬壁附近(30 mm

圖8 電場強度Ei(上)和電子數密度ne(下)的分布Fig.8.Distribution of electric field intensity EI (upper) and electron number density ne (lower).

圖9 感應線圈中心(x=68 mm)空氣粒子徑向摩爾分數分布Fig.9.Mole fraction of air species along the radial direction at the coil center x=68 mm.

圖10為等離子體平動溫度和電子溫度分布云圖,最大電子溫度(9921 K)和最大平動溫度(8507 K)均出現在第二個感應線圈下方靠近等離子體炬壁處.在60 ≤x≤ 85 mm,20 ≤y≤ 30 mm區(qū)域,電子溫度(8500 ≤T≤ 9500 K)比平動溫度高約1000 K,即在該區(qū)域等離子體處于熱力學非平衡態(tài).然而,在靠近中軸線附近,電子溫度與平動溫度基本相等,氣體溫度T≈ 6500 K,說明該區(qū)域空氣重粒子(如N,O)與電子之間的能量交換基本達到了熱力學平衡,此時等離子體流動趨于局域熱力學平衡態(tài).

圖10 等離子體炬內平動溫度(上半部分)和電子溫度(下半部分)分布云圖Fig.10.Distributions of translational (upper) and electronic temperatures (lower) in the torch.

4 結 論

本文基于多物理場耦合流動仿真研究了空氣ICP流場與電磁場的分布特性及其相互作用機理,通過數值模擬獲得了100 kW級ICP炬內等離子體流動的電子溫度、平動溫度、粒子摩爾分數、速度、壓強、洛倫茲力、電場強度、焦耳加熱率等參數的分布規(guī)律,研究發(fā)現:

1)在進氣口與第二圈感應線圈之間出現了較強的渦流,該渦流與感應線圈區(qū)的負壓強梯度和電磁加熱現象有很大關系.徑向洛倫茲力方向始終為負說明因趨膚效應在壁面附近生成的自由電子始終受到指向ICP炬中軸線的洛倫茲力作用; 且徑向洛倫茲力的最大值比軸向洛倫茲的高3.95倍,這表明動量傳遞以徑向為主; 空氣粒子的焦耳熱效應也受徑向洛倫茲力的影響.

2)電場強度虛部EI的最大值比電場強度實部ER的大2.9倍,負的電場虛部EI主要是由外加高頻電流產生,而正的電場虛部則是由等離子體內部的自由電子產生,在第二圈感應線圈下方距離等離子體炬內壁面5.5 mm處自由電子的數密度達到最大值.

3)空氣在感應線圈下方發(fā)生劇烈的離解和電離反應,氧分子被完全解離和電離成氧原子或離子;91%的氮分子在中軸線附近被離解成N原子,N和O原子在徑向位置y≤ 30 mm之前為最主要的化學組分.本次模擬得到的空氣ICP最大電子溫度(9921 K)和最大平動溫度(8507 K)均出現在第二圈感應線圈下方靠近等離子體炬壁面處.在60 mm ≤x≤ 85 mm且 20 mm ≤y≤ 30 mm區(qū)域,電子溫度比平動溫度高大約1000 K,流動處于熱力學非平衡狀態(tài).然而在靠近中軸線附近,電子溫度與平動溫度基本相等(約為6500 K),在該區(qū)域空氣重粒子與電子之間的能量交換處于局域熱平衡.

感謝國家超級計算廣州中心提供的天河二號超級計算機計算服務支持.感謝日本宇宙航空研究開發(fā)機構Kazuhiko Yamada等的建議與支持.