喻明浩

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

1 引 言

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

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

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

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

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

圖1 ICP風(fēng)洞系統(tǒng)結(jié)構(gòu)布局Fig.1.Schematic diagram of the ICP wind tunnel system.

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

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

2 控制方程與數(shù)值方法

本文所研究的100 kW級(jí)ICP風(fēng)洞幾何形狀與計(jì)算網(wǎng)格如圖2所示.圖2(a)為遠(yuǎn)場(chǎng)電磁場(chǎng)和ICP炬流場(chǎng)的計(jì)算網(wǎng)格,電磁場(chǎng)的計(jì)算域?yàn)楱C120 mm ≤x≤ 360 mm,0 ≤y≤ 206 mm,由101×66個(gè)網(wǎng)格節(jié)點(diǎn)組成,其中ICP炬(ICP torch)的流場(chǎng)網(wǎng)格由101×38個(gè)網(wǎng)格節(jié)點(diǎn)組成,為了準(zhǔn)確計(jì)算線圈附近的電磁場(chǎng)與進(jìn)氣口處氣流速度,計(jì)算網(wǎng)格在感應(yīng)線圈區(qū)和進(jìn)氣口處進(jìn)行了加密處理.圖2(b)為ICP炬的結(jié)構(gòu)示意圖,感應(yīng)線圈圍繞ICP炬管外壁面纏繞三圈,線圈的內(nèi)徑和外徑分別為47 mm和55 mm,起始線圈的中心坐標(biāo)為(x,y)=(51,51) mm,相鄰感應(yīng)線圈的中心線間隔為16.5 mm.

ICP風(fēng)洞正常工作時(shí),感應(yīng)線圈上通有高頻交變電流,電流頻率為1.78 MHz,工作氣體如空氣或氮?dú)鈴谋诿娓浇莫M窄開口注入,經(jīng)電火花高能點(diǎn)火后,工作氣體開始發(fā)生離解和電離反應(yīng)并產(chǎn)生自由電子.此后,帶電粒子在交變電流產(chǎn)生的交變電磁場(chǎng)作用下進(jìn)一步發(fā)生劇烈的碰撞電離、副電離等反應(yīng),并釋放出大量能量而形成溫度高達(dá)10000 K的ICP氣流,最終該等離子體氣流從ICP炬出口流入大空間試驗(yàn)腔內(nèi)進(jìn)行飛行器氣動(dòng)性能測(cè)試、熱防護(hù)材料燒蝕試驗(yàn)等實(shí)驗(yàn)研究.

圖2 ICP炬計(jì)算網(wǎng)格和幾何結(jié)構(gòu) (a) 電磁場(chǎng)與流場(chǎng)計(jì)算網(wǎng)格; (b) 等離子體炬幾何結(jié)構(gòu)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,質(zhì)量流率m=1.8 g/s,工作氣壓p=10 kPa,電流頻率f=1.76 MHz,工作氣體為空氣.對(duì)于高功率ICP風(fēng)洞,由于存在較大的熱能損失(比如,冷卻水引起的熱損失,電路上的電阻抗引起的能量損失以及輻射損失),沉積到等離子體流動(dòng)中的凈輸入功率通常與從供電系統(tǒng)讀取的面板輸入功率有所不同[36],Fujita等[28,37]對(duì)此進(jìn)行了實(shí)驗(yàn)與仿真研究,研究結(jié)果表明: 對(duì)于110 kW ICP風(fēng)洞其熱效率通常為1/3.因此,本研究ICP風(fēng)洞的熱效率h也設(shè)定為1/3.被等離子體實(shí)際吸收的沉積功率等于熱效率h與面板輸入功率的乘積.仿真計(jì)算時(shí),為了確保能量守恒,每一步流場(chǎng)迭代均使所有流體微元的焦耳加熱率體積積分之和等于等離子體的沉積功率

2.1 流場(chǎng)方程組

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

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

高溫氣體的壓強(qiáng)p由下式計(jì)算:

氣體的總內(nèi)能E定義為

氣體粒子的平動(dòng)能(Etr)、轉(zhuǎn)動(dòng)能(Erot)、振動(dòng)能(Evib)和電子內(nèi)能(Ee)分別由下式計(jì)算:

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

2.2 電磁場(chǎng)方程組

ICP的電磁場(chǎng)分布通常可通過求解以下磁矢勢(shì)方程得到:

其中A為磁矢勢(shì)和i分別表示線圈電流的角頻率、電流密度、真空介電常數(shù)和復(fù)數(shù)單位為線圈電流驅(qū)動(dòng)頻率.此外,由于磁矢勢(shì)A與磁場(chǎng)強(qiáng)度H和電場(chǎng)強(qiáng)度E有如下關(guān)系[20]:

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

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

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

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

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

2.3 空氣化學(xué)反應(yīng)模型

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

式中的化學(xué)反應(yīng)系數(shù)Cr,n和θr取自文獻(xiàn)[33,35].

逆向化學(xué)反應(yīng)速率kb,r由平衡常數(shù)Keq進(jìn)行計(jì)算:

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

式中,nsb,r和nsf,r是組分s在進(jìn)行第r個(gè)化學(xué)反應(yīng)前后的化學(xué)計(jì)量系數(shù),Ms代表物質(zhì)s的摩爾質(zhì)量,下標(biāo)s和j代表物質(zhì)組分編號(hào).

表1 空氣化學(xué)反應(yīng)模型Table 1.Chemical reaction model of air.

2.4 內(nèi)部能量交換模型

由于電子、分子和原子之間存在彈性或非彈性碰撞,碰撞時(shí)會(huì)產(chǎn)生平動(dòng)能、轉(zhuǎn)動(dòng)能、振動(dòng)能或電子能量的交換,為了準(zhǔn)確模擬這些內(nèi)能交換過程,根據(jù)相應(yīng)的數(shù)學(xué)模型計(jì)算粒子內(nèi)能交換率,并將計(jì)算得到的轉(zhuǎn)動(dòng)能交換率Sint,rot、振動(dòng)能交換率Sint,vib和電子能量交換率Sint,e作為能量源項(xiàng)添加到(2)式中的轉(zhuǎn)動(dòng)、振動(dòng)和電子能量守恒方程中.轉(zhuǎn)動(dòng)、振動(dòng)和電子能量交換率Sint,rot,Sint,vib,Sint,e的具體計(jì)算方法如下:

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

2.5 湍流運(yùn)輸方程

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

湍流黏度由下式計(jì)算:

上述方程中的模型常數(shù)如下:

模型函數(shù)表示為

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

2.6 邊界條件和數(shù)值求解方法

2.6.1 邊界條件

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

2.6.2 數(shù)值求解方法

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

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

3 結(jié)果與討論

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

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

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

圖3 等離子體炬內(nèi)氣體流線和速度矢量(上半部分)以及電子溫度云圖(下半部分)分布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給出了軸向和徑向洛倫茲力分布.洛倫茲力作為能量守恒方程的慣性力源項(xiàng)參與到整個(gè)流場(chǎng)的迭代求解中.由圖5可見: 軸向洛倫茲力在感應(yīng)線圈的第一與第二圈之間為正值,即其方向向右.然而在感應(yīng)線圈第二圈與第三圈之間變?yōu)樨?fù)值,即軸向洛倫茲力的方向反向.而對(duì)于徑向洛倫茲力(圖5下半部分),其值始終為負(fù)值,即其方向始終由壁面指向中心軸線,這表明因趨膚效應(yīng)在炬壁附近生成的自由電子將始終受到指向等離子體炬中軸線的慣性力,且徑向洛倫茲力的最大值比軸向洛倫茲的最大值高3.95倍,說明等離子體的徑向動(dòng)量傳遞比軸向動(dòng)量傳遞強(qiáng)烈的多.

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

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

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

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

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

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

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

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

圖9 感應(yīng)線圈中心(x=68 mm)空氣粒子徑向摩爾分?jǐn)?shù)分布Fig.9.Mole fraction of air species along the radial direction at the coil center x=68 mm.

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

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

4 結(jié) 論

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

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

2)電場(chǎng)強(qiáng)度虛部EI的最大值比電場(chǎng)強(qiáng)度實(shí)部ER的大2.9倍,負(fù)的電場(chǎng)虛部EI主要是由外加高頻電流產(chǎn)生,而正的電場(chǎng)虛部則是由等離子體內(nèi)部的自由電子產(chǎn)生,在第二圈感應(yīng)線圈下方距離等離子體炬內(nèi)壁面5.5 mm處自由電子的數(shù)密度達(dá)到最大值.

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

感謝國(guó)家超級(jí)計(jì)算廣州中心提供的天河二號(hào)超級(jí)計(jì)算機(jī)計(jì)算服務(wù)支持.感謝日本宇宙航空研究開發(fā)機(jī)構(gòu)Kazuhiko Yamada等的建議與支持.