激光維持等離子體多物理場耦合模型與仿真*

張東荷雨 劉金寶 付洋洋

(清華大學(xué)電機(jī)工程與應(yīng)用電子技術(shù)系,北京 100084)

1 引言

激光維持等離子體(laser-sustained plasma,LSP)是一種利用激光與電離氣體相互作用產(chǎn)生的等離子體.LSP 可以從激光輻射場中吸收能量并維持在接近熱力學(xué)平衡狀態(tài),通常也被稱為連續(xù)光放電(continuous optical discharge,COD)[1,2].LSP具有高溫、高密度、高輻射強(qiáng)度等特點[3],在激光推進(jìn)[4]、光學(xué)檢測[5]、等離子體化學(xué)和材料加工[6]等領(lǐng)域展現(xiàn)了重要的應(yīng)用價值.

在集成電路制造領(lǐng)域,LSP 已經(jīng)發(fā)展成為半導(dǎo)體晶圓光學(xué)檢測中的主流光源技術(shù)[7].晶圓缺陷檢測貫穿半導(dǎo)體工藝的各個環(huán)節(jié),以保證芯片產(chǎn)品的質(zhì)量和良率.其按技術(shù)原理可分為電子束檢測和光學(xué)檢測兩大類,前者主要用于缺陷的精細(xì)檢測,而后者基于光學(xué)成像原理對表面缺陷進(jìn)行快速定位.鑒于晶圓表面缺陷尺寸小、介質(zhì)種類多,檢測光源須同時具備高亮度、寬光譜范圍等特點.與傳統(tǒng)的電弧光源相比,LSP 光源體積小、能量沉積效率高,同等功率下LSP 光源發(fā)光強(qiáng)度較短弧氙燈可提高幾個數(shù)量級[3],因此更易滿足半導(dǎo)體晶圓檢測應(yīng)用的需求.此外,LSP 穩(wěn)定工作時無需高壓電極長時間參與放電,電極受到的損傷相對較輕、壽命更長,等離子體狀態(tài)及發(fā)光質(zhì)量受到電極材料濺射和蒸發(fā)的影響較小,可長時間穩(wěn)定工作.

圖1 所示為一種用于產(chǎn)生LSP 光源的實驗裝置.預(yù)填充了高氣壓稀有氣體的密閉腔體內(nèi)置電極,通過在上下兩極之間外施高壓脈沖產(chǎn)生初始等離子體.腔體外的激光束經(jīng)過擴(kuò)束、聚焦等處理后,被照射于兩電極之間的等離子體區(qū)域,此處為能量集中交換區(qū)域.如圖2 所示,對于中心波長為10.6 μm 的 CO2激光,其能量主要通過逆韌致輻射方式注入等離子體,同時等離子體內(nèi)也不斷進(jìn)行著自由-自由(free-free,ff)、自由-束縛(free-bound,fb)和束縛-束縛(bound-bound,bb)多種輻射過程,向外輻射能量并在電子躍遷過程中發(fā)出特定波長的光,另有部分能量通過熱擴(kuò)散等過程被損失.在實際應(yīng)用場景下,LSP 輻射出的光會被收集并進(jìn)行均質(zhì)化處理,隨后被照射到待檢測物體表面.圖1 右側(cè)簡要展示了LSP 光源在半導(dǎo)體晶圓缺陷檢測中的應(yīng)用,極高亮度的LSP 輻射光在待檢測晶圓表面不斷移動,通過收集并分析不同位置反射或折射的光信號,可實現(xiàn)對晶圓表面結(jié)構(gòu)缺陷的檢測和定位.

圖1 LSP 實驗裝置結(jié)構(gòu)示意圖Fig.1.Schematic diagram of the LSP experimental setup.

圖2 CO2 激光條件下LSP 中主要能量傳遞機(jī)制,包括: 逆韌致輻射(ff)吸收; 自由-自由(ff)/自由-束縛(fb)輻射; 束縛-束縛(bb)輻射; 熱擴(kuò)散Fig.2.Main energy transfer processes of plasma sustained by CO2 laser including bremsstrahlung absorption,free-free(ff)/free-bound (fb) emission,bound-bound (bb) emission,and thermal diffusion.

基于LSP 原理,日本Energetiq 公司[8]研發(fā)了低功率光源LDLS (laser-driven light source),可以提供輻射波長從深紫外(170 nm)到紅外波段(2.1 μm)的高強(qiáng)度光照.LDLS 光源采用近紅外連續(xù)半導(dǎo)體激光器作為激勵源,可產(chǎn)生溫度高于10000 K、微米尺寸的Xe LSP,在各波段的輻照強(qiáng)度均明顯高于傳統(tǒng)氣體放電光源.經(jīng)測試,50 WLDLS 光源在紫外至可見光波段的光照亮度是75 W 短弧氙燈的5—10 倍[9].美國科磊公司(KLATencor)基于LSP 光源技術(shù)開發(fā)了商用的寬光譜等離子圖案晶圓缺陷檢測系統(tǒng),如KLA-39xx、KLA-29xx 系列等 .此外,針對如何提高LSP 光源亮度與穩(wěn)定性,科磊公司嘗試在腔體方案中引入反向渦流氣體腔[10]、漩渦軸[11]等設(shè)計,以改善高功率條件下封閉腔體內(nèi)可能出現(xiàn)的局部過熱現(xiàn)象.雖然國外已經(jīng)出現(xiàn)了較為成熟的商用LSP 光源設(shè)備,但由于LSP 涉及流場、溫度場和激光能量場等多物理場的耦合,且相關(guān)實驗診斷難度大、機(jī)理復(fù)雜度高,因此仍存在大量物理現(xiàn)象和底層物理機(jī)制尚未充分闡明.

本文首先簡要回顧了LSP 在實驗與模擬方面的研究進(jìn)展,對相關(guān)研究結(jié)果進(jìn)行梳理與總結(jié).在此基礎(chǔ)上,詳細(xì)介紹LSP 二維多物理場耦合流體模型,給出各種氣體物性參數(shù)與輸運(yùn)參數(shù)的計算方法.所建模型可用于研究氣體種類、腔體結(jié)構(gòu)、激光參數(shù)(功率、模式)等因素對LSP 時空特性及性能的影響機(jī)制,滿足不同應(yīng)用場景下LSP 優(yōu)化設(shè)計的需求.基于仿真模型,觀察到了LSP 的穩(wěn)態(tài)結(jié)構(gòu)、演化過程及振蕩現(xiàn)象,研究結(jié)果將為LSP 參數(shù)調(diào)控、系統(tǒng)性能優(yōu)化以及相關(guān)技術(shù)難題的解決提供理論指導(dǎo).

2 激光維持等離子體的實驗與模擬

自連續(xù)高功率 CO2激光器出現(xiàn)以后,蘇聯(lián)科學(xué)院力學(xué)問題研究所Generalov 等[12]首次在聚焦激光束焦點附近產(chǎn)生并維持了高溫等離子體,成功驗證了Raizer[1]關(guān)于LSP 的理論設(shè)想.1972 年,阿夫科-埃弗里特實驗室(AERL) Kantrowitz[4]首次提出了激光推進(jìn)的概念(其中LSP 是激光推進(jìn)中實現(xiàn)能量轉(zhuǎn)化的重要方式),即推進(jìn)劑通過吸收遠(yuǎn)程高功率激光的輻射在噴管中被加熱膨脹并產(chǎn)生推力.相較于常見的化學(xué)推進(jìn)方法,激光推進(jìn)可以突破燃料最高燃燒溫度的限制,能量吸收效率(50%—86%)與熱能轉(zhuǎn)化效率更高[13–15].實驗中發(fā)現(xiàn),激光光束聚焦射入高壓氣體后,穩(wěn)態(tài)的LSP會出現(xiàn)在焦點上游,其核心位置和尺寸相對穩(wěn)定,此時LSP 吸收的入射激光強(qiáng)度恰好與熱輻射、傳導(dǎo)和對流損失相平衡.平衡位置與尺寸大小受到激光入射功率、激光模式[16,17]、聚焦能力[6,18–20]、氣體壓強(qiáng)及對流流速[6,21,22]等參數(shù)的影響.現(xiàn)有的實驗與模擬研究結(jié)果顯示,特定條件下可以觀察到呈現(xiàn)雙峰形態(tài)的LSP 溫度分布[20,23,24].改變實驗條件可以調(diào)控出現(xiàn)高溫核心的個數(shù): Jeng 和Keefer[25]的模擬結(jié)果表明低流速、高激光功率、高氣壓、高F數(shù)的情況下,LSP 更傾向于出現(xiàn)雙核心結(jié)構(gòu).

已有研究表明,LSP 的各項參數(shù)會顯著影響其發(fā)光亮度與輻射特性.2016 年,俄羅斯科學(xué)院力學(xué)問題研究所Zimakov 等[19]與北京航空航天大學(xué)Shi 等[16]的實驗結(jié)果均表明,更高的功率會增加LSP 的長度與核心溫度,進(jìn)而獲得輻射強(qiáng)度更高的光源.同時降低激光F數(shù)會顯著提高激光閾值功率,產(chǎn)生的LSP 長度與長寬比減小,同時LSP輻射光譜中波長小于 400 nm 的部分占總光譜強(qiáng)度的比例提升.華中科技大學(xué)Hu 等[20]通過在實驗中改變LSP 的腔內(nèi)氣壓,確定了在相同激光功率下,高氣壓條件更容易獲得小體積、高溫、高亮度的LSP 光源.

早在1970 年,Raizer[1]提出了亞音速條件下激光火花(放電)傳播的一維計算模型,用來確定連續(xù) CO2激光維持等離子體的可行性及功率閾值.1974 年,美國佛羅里達(dá)大學(xué)Batteh 與Keefer[26]將模型擴(kuò)展至準(zhǔn)二維,增加了對徑向熱傳導(dǎo)損失和激光聚焦特性的描述.在此基礎(chǔ)上,美國伊利諾伊大學(xué)Glumb 與Krier[27]引入收斂的激光束及更準(zhǔn)確的氣體特性對模型進(jìn)一步拓展,但缺少對流體速度變化的考量.1985 年,美國賓夕法尼亞州立大學(xué)Merkle 等[28,29]首次提出了較完備的LSP 二維穩(wěn)態(tài)模型,包含了激光光束的聚焦作用和截面徑向能量分布,并考慮了氣體物理參數(shù)隨溫度的變化,但該模型將吸收系數(shù)假設(shè)為常數(shù),存在較大不足.美國田納西大學(xué)Jeng 等[30]在此基礎(chǔ)上對二維模型進(jìn)行修正與改進(jìn),通過與實驗結(jié)果[22]進(jìn)行對比驗證了穩(wěn)態(tài)模型的準(zhǔn)確性.1996 年,Conrad 等[31]提出了含時的LSP 二維流體模型,但該模型對激光吸收系數(shù)以及LSP 輸運(yùn)系數(shù)、物性參數(shù)的計算較為簡化,盡管背景氣體為空氣,但未考慮原子譜線引起的熱輻射損耗.2009 年,土耳其中東技術(shù)大學(xué)Rafatov 等[32]引入了能夠更精確描述激光輻射傳播過程的準(zhǔn)光學(xué)近似,假設(shè)光束在等離子體中沿拋物線傳播,模擬結(jié)果能夠與早期實驗研究結(jié)果相符合[33,34],并總結(jié)了激光功率、F數(shù)[35]、氣流流速對LSP 位置、溫度的影響.同年,美國賓夕法尼亞州立大學(xué)Akarapu 等[17]引入了激光功率場,以代替光線追跡法對LSP 中的功率沉積進(jìn)行模擬,給出了不同激光模式下LSP 能量吸收效率和形態(tài)區(qū)別,驗證了環(huán)形激光模式下LSP 吸收效率低但體積更小的實驗結(jié)論.該模型雖然可以在二維上較好地模擬LSP 的穩(wěn)態(tài)溫度分布,但所解流體方程中未考慮含時項,不能對LSP 的暫態(tài)發(fā)展過程進(jìn)行研究.

隨著對LSP 特性研究的深入,實驗中觀察到的LSP 較長時間尺度不穩(wěn)定現(xiàn)象引起了越來越多的關(guān)注.Zimakov 等在實驗中利用紋影法[36]、發(fā)射光譜頻譜分析[37]等手段觀察到了LSP 的周期波動規(guī)律,總結(jié)了氣體壓強(qiáng)對波動頻率的影響.在此基礎(chǔ)上,Lavrentyev 等[38]通過類比火焰的波動規(guī)律提出了定量估計LSP 振蕩周期頻率的計算公式,Kotov 等[39]則對該研究提出的半經(jīng)驗公式進(jìn)行了補(bǔ)充,通過自然對流模型與實驗現(xiàn)象的對比分析,驗證了LSP 波動產(chǎn)生只與流場有關(guān)的猜想.

目前,理論計算模型能夠較為合理地描述穩(wěn)態(tài)下LSP 溫度分布、核心位置和長度等狀態(tài)參數(shù),能夠更加準(zhǔn)確地預(yù)測實驗結(jié)果,為掌握LSP 中的復(fù)雜過程、關(guān)鍵影響因素、調(diào)控手段等提供了有力研究工具.然而,現(xiàn)有研究尚存在一些不足: 目前LSP的模擬研究多為穩(wěn)態(tài)模型,僅能得到LSP 的最終穩(wěn)定狀態(tài),缺少對激光輻射發(fā)展過程進(jìn)行描述的含時模型; 實驗中觀察到的LSP 不穩(wěn)定波動現(xiàn)象尚未在模擬研究中體現(xiàn); 已有研究多聚焦于單一氣體種類下的LSP 模型,可推廣性不強(qiáng).為彌補(bǔ)以上不足,本文建立了LSP 二維含時流體模型,較為精確地考慮了激光吸收與輻射損失功率密度項,給出了計算Xe 等離子氣體的物性參數(shù)及輸運(yùn)系數(shù)的方法.

3 計算模型

根據(jù)Raizer[3]的研究,LSP 的中心溫度可以達(dá)到 10000—20000 K,甚至更高,通常滿足局部熱力學(xué)平衡(local thermodynamic equilibrium,LTE).基于LTE 假設(shè),本文在模型中認(rèn)為所有粒子(電子、離子和中性粒子)分布在統(tǒng)計上滿足玻爾茲曼關(guān)系,且各能量自由度上具有相同的溫度T,即T=Te=Ti=Tg.LSP 氣壓在幾bar 到十幾bar范圍(1 bar=105Pa),因此可以采用流體模擬的方法.在流速小于音速,以及較高的溫度梯度條件下,可采用層流可壓(laminar,compressible)方法.層流可壓流體可由質(zhì)量守恒(連續(xù))方程、動量守恒方程和能量守恒方程描述.質(zhì)量守恒(連續(xù))方程為

其中,ρ 是質(zhì)量密度,u是流體速度.Navier-Stokes方程(動量方程)為

等式左側(cè)是流體動量變化速率,右側(cè)是流體受到的力.其中,p是氣體壓強(qiáng),fb是體積力,一般可以忽略,是黏滯力張量,采用牛頓流體計算公式表示為,其中,μv是黏性系數(shù),是單位張量.能量方程(忽略氣體動能)為

其中,左側(cè)是流體焓變速率,T是溫度,Cp是定壓比熱容.右側(cè)是流體的加熱或損失功率密度,單位均為W/m3,包括擴(kuò)散功率密度?·(kT?T)、激光加熱功率密度Pabs,以及輻射功率密度 ?Prad,kT是熱導(dǎo)率.

(1)式—(3)式與理想氣體狀態(tài)方程,即

構(gòu)成關(guān)于ρ,u,p和T的完備方程組.(4)式中Mmol是摩爾質(zhì)量,R是摩爾氣體常數(shù).為求解式(1)式—(4)式,還需確定參數(shù)Cp,kT,μv,Pabs,Prad,下面給出這些參數(shù)的具體計算方法.

3.1 參數(shù)計算

LTE 等離子體的μv,Cp和kT與氣體的組分有關(guān),而氣體的平衡組分隨著T和p的不同而變化,因此首先需要確定等離子體體系的組成[40,41].LTE 等離子體中各物種處于熱力學(xué)平衡狀態(tài),此時系統(tǒng)的吉布斯自由能最小[42].系統(tǒng)吉布斯自由能G可表示為

其中,μi是物種i的化學(xué)勢,Ni是物種i的粒子數(shù),kB是玻爾茲曼常數(shù),qi是物種i的配分函數(shù),Nspe為物種數(shù).本文計算中考慮的組分包括: 電子,Xe,Xe+,Xe2+,Xe3+與Xe4+六種粒子,即物種數(shù)Nspe=6.粒子的配分函數(shù)qi可拆分為關(guān)于平動能的qi,tr與關(guān)于內(nèi)能的qi,int兩部分,即

其中,

其中,h是普朗克常量,V是體積,m是質(zhì)量,gj和Ej是能級j對應(yīng)的簡并度和能量.除以上約束條件外,Ni還要滿足元素守恒,此時求得使G達(dá)到最小的Ni值,即可得到LTE 等離子體各組分的數(shù)密度,求解這一約束優(yōu)化問題可采用拉格朗日乘子法.圖3 給出了不同氣壓條件下穩(wěn)定狀態(tài)的Xe 等離子體中各粒子數(shù)密度隨溫度的變化.

圖3 不同壓強(qiáng)的平衡狀態(tài)下Xe 氣體組分Fig.3.Composition of Xe plasma at different pressures.

3.2 激光與等離子體相互作用

當(dāng)激光頻率大于等離子體頻率時,激光可以進(jìn)入等離子體并在等離子體中傳播.不考慮激光在等離子體中傳播產(chǎn)生的電磁場影響,基于幾何光學(xué)假設(shè),根據(jù)Beer-Lambert 定律,沿激光傳播路徑s上激光強(qiáng)度I的變化可表示為

其中,s表示激光傳播路徑上的單位矢量;α 表示吸收系數(shù),即激光經(jīng)過單位長度被等離子體吸收的強(qiáng)度百分比.當(dāng)采用的激光束中心波長不同時,激光對等離子體的能量注入(加熱)機(jī)制存在一定差異.如用中心波長為10.6 μm 的CO2激光器照射 H2,Ar,Xe 等種類氣體實現(xiàn)LSP,等離子體主要通過逆韌致輻射(inverse bremsstrahlung)過程吸收激光能量,即電子從光子中吸收能量,再經(jīng)過碰撞過程將能量傳遞給重粒子(電離和激發(fā)),進(jìn)而提高等離子體溫度.若采用高功率近紅外(near-infrared,NIR)激光器(常見中心波長包括1064,1070,1080 nm等)作為光源,則還需要考慮束縛-束縛(bb)吸收形式.

3.2.1 ff 吸收

對于逆韌致輻射(ff)吸收,吸收系數(shù)αff的計算式為[17,43]

其中,Z,ni為離子電荷數(shù)和數(shù)密度;e,me,ne,Te分別為電子的電荷量(絕對值)、質(zhì)量、數(shù)密度和溫度;c為光速;ν 和νp=分別為激光頻率和等離子體頻率.lnΛ為庫侖對數(shù),Λ等于vT/(2πνppmin)與vT/(2πνpmin) 中的較小值,vT為 電子的熱速度,pmin為電子-離子碰撞過程中的最小瞄準(zhǔn)距離.

3.2.2 bf 吸收

其中,n是i粒子j能級的主量子數(shù),Zn是非氫原子的有效電荷數(shù),Pn是原子或離子在n層軌道的電子占據(jù)數(shù)(occupation number),Uij是i粒子j能級的電離能.當(dāng)光子能量hν大于j能級原子(或離子)進(jìn)一步電離所需要的能量Uij時,光電離反應(yīng)截面＞0,否則=0.

3.2.3 bb 吸收

對于束縛-束縛(bb)吸收,吸收系數(shù)與激光頻率有關(guān),首先給出特定頻率下的吸收系數(shù)計算公式:

其中,A21為電子從能級2 躍遷到能級1 的自發(fā)輻射躍遷概率,下標(biāo)2 代表高能級、下標(biāo)1 代表低能級,n2為處于高能級的電子數(shù)密度,?(ν) 為等離子體吸收光譜的線型函數(shù).

若計算某一溫度下等離子體的bb 吸收系數(shù),還需要考慮入射激光的光譜形狀,并在相關(guān)頻率波段范圍內(nèi)進(jìn)行積分:

其中f(ν) 為入射激光光譜的線型函數(shù),ν1和ν2分別為等離子體能夠吸收激光的最小頻率與最大頻率.

3.2.4 總吸收系數(shù)

LSP 的總激光吸收系數(shù)等于以上三種吸收系數(shù)之和:

在入射CO2激光中心波長為 10.6 μm 的條件下,不同壓強(qiáng)下Xe 的吸收系數(shù)隨溫度的變化曲線如圖4 所示.

圖4 10.6 μm CO2 激光條件下不同壓強(qiáng)Xe 等離子體的吸收系數(shù)隨溫度變化Fig.4.Absorption coefficient of Xe plasmas with temperature at different pressures under the condition of 10.6 μm CO2 laser.

3.3 等壓比熱容

等壓比熱容Cp描述了焓H與溫度T之間的關(guān)系.首先定義計算系統(tǒng)的質(zhì)量密度ρ,即

其中,ni和mi分別為物種i的數(shù)密度和粒子質(zhì)量.

焓H與配分函數(shù)等參數(shù)的關(guān)系為

其中,Hf,i為物種i的合成焓.

定壓比熱容Cp與焓H、溫度T的關(guān)系為

不同壓強(qiáng)下Xe 等離子體定壓比熱容隨溫度變化曲線如圖5(a)所示.

圖5 不同壓強(qiáng)下Xe 等離子體各參數(shù)隨溫度的變化 (a)恒壓比熱容; (b)黏滯系數(shù); (c)熱導(dǎo)率; (d)各項輻射損失功率密度Fig.5.Parameters of Xe plasmas scaling with temperature at different pressures: (a) Specific heat capacity; (b) viscosity; (c) thermal conductivity; (d) radiative loss power density.

3.4 碰撞積分

根據(jù)Chapman 和Cowling[45]的理論,黏滯系數(shù)μv與熱導(dǎo)率kT等輸運(yùn)參數(shù)可通過構(gòu)建碰撞積分計算得到.計算以上兩項輸運(yùn)系數(shù)需用到的碰撞積分種類在表1 中列出.

表1 輸運(yùn)系數(shù)與用于計算的碰撞積分Table 1.Transport coefficients and the collision integrals used for calculation.

物種i與物種j相互作用的碰撞積分可表示為

其中,b,χij分別為瞄準(zhǔn)距離和碰撞偏角;εij=是由組分碰撞時相對速度g定義的初始動能; 上標(biāo)l代表碰撞積分的類型,s代表碰撞積分的階數(shù);χij是碰撞參數(shù)b、動能ε和粒子間相互作用勢V(r) 的函數(shù);r為相互作用粒子間的距離.

不同粒子對的相互作用需采用合適的勢能函數(shù)進(jìn)行描述,對應(yīng)的碰撞積分計算方法也不同.本文中用到的Xe 等離子體中包含 e-Xe,Xe-Xe,e-,Xe-Xe+與共5 種類型的粒子對,描述這些粒子對之間相互作用的碰撞積分計算方法見表2.下面介紹幾種描述粒子間相互作用的碰撞積分計算方法.

表2 Xe 等離子體中粒子對相互作用對應(yīng)的碰撞積分計算方法Table 2.Calculation of collisional integrals corresponding to particle pair interactions in Xe plasmas.

3.4.1 指數(shù)排斥勢

指數(shù)排斥勢(exponential repulsive potential)常被用于描述中性原子之間或中性原子與離子之間的相互作用勢[50,51],它的表達(dá)式為

其中,A和ρ 兩個參數(shù)描述了指數(shù)排斥勢的形狀.其對應(yīng)的碰撞積分計算公式為

其中β=ln[A/(kBT)],I(l,s)為對應(yīng)不同l和s組合的特定值,一般可通過查表得到I(l,s)與β 的關(guān)系[52].

3.4.2 蘭納-瓊斯勢

另一種用于描述中性原子之間或中性原子與帶電粒子之間相互作用勢的勢能函數(shù)類型為蘭納-瓊斯勢(Lennard-Jones potential),其使用較廣泛的表達(dá)形式為

其中,ε0和σ 兩個參數(shù)描述了該勢能函數(shù)形狀.其對應(yīng)的折合碰撞積分通過查表可得[53],并由折合碰撞積分(reduced collision integrals)可以計算蘭納-瓊斯勢對應(yīng)的碰撞積分.

按照Chapman 的處理方式,折合碰撞積分定義為其他勢能的碰撞積分與硬球碰撞積分的比值,即

3.4.3 電荷交換相互作用

電荷交換相互作用的實驗測定截面Qex,ij可以參數(shù)化表示為

其中,g為相對速度,A、B為擬合參數(shù)(見表2 第三行).將(26)式代入(19)式,可以得到對應(yīng)的碰撞積分計算式為[54]

3.4.4 庫侖勢與屏蔽庫侖勢

用于描述帶電粒子間相互作用的庫侖勢能可以表示為

其中,Zi為帶電粒子i的凈電荷數(shù).

由于等離子體的聚集效應(yīng),須對庫侖勢做修正,即采用屏蔽庫侖勢計算帶電粒子之間的相互作用.修正后的屏蔽庫侖勢公式可表示為

其中,V0=(Z1Z2e2)/(4πε0λD),λD為德拜長度.可通過查表得到對應(yīng)碰撞積分與溫度的關(guān)系[49].

3.5 黏滯系數(shù)

根據(jù)Chapman 和Cowling[45]的理論,多組分系統(tǒng)的黏滯系數(shù)可表示為

其中,列向量x={x1,x2,···}T中xi為各物種的體積分?jǐn)?shù);I為單位列向量;H為二維矩陣,其對角元素Hii和非對角元素Hij分別為

式中,μij為二元黏性系數(shù),其計算式為

不同壓強(qiáng)下Xe 等離子體黏滯系數(shù)隨溫度變化曲線如圖5(b)所示.

3.6 熱導(dǎo)率

熱導(dǎo)率由幾部分共同組成,即

其中,kT,hp為重粒子貢獻(xiàn)的熱導(dǎo)率,kT,e為電子貢獻(xiàn)的熱導(dǎo)率.重粒子貢獻(xiàn)kT,hp表達(dá)式為

其中,xT與(30)式中定義相同,I為單位列向量.矩陣L的對角元素Lii和非對角元素Lij分別為

其中,Nhea為重粒子的種類數(shù),二元熱導(dǎo)率kij可表示為

不同類型碰撞積分的比例Aij和Bij定義為

Devoto[55]給出了電子貢獻(xiàn)kT,e的三階計算公式:

其中,qxy與平均截面有關(guān),計算方法如下[55]:

不同壓強(qiáng)下Xe 等離子體熱導(dǎo)率隨溫度變化曲線如圖5(c)所示.

3.7 輻射功率

如圖2 所示,輻射損失包含了ff,fb 和bb 過程.其中ff 和fb 過程產(chǎn)生連續(xù)譜,bb 過程產(chǎn)生線狀譜.不同輻射機(jī)制引起的輻射損失功率密度求和,得到輻射損失功率密度Prad.

3.7.1 ff 輻射

式中α 為精細(xì)結(jié)構(gòu)常數(shù),a0為Bohr 半徑,Ry為Rydberg 常量,Zeff為有效凈電荷數(shù).α,a0和Ry計算式(均為國際單位制)為

有效凈電荷數(shù)Zeff定義式為[57]

其中,Zi為第i種離子的凈電荷數(shù),Nhp為重粒子種類數(shù).

由空間上的各向同性,立體角積分即為系數(shù) 4π.

3.7.2 fb 輻射

由Kirchhoff 定律,LTE 狀態(tài)下等離子體的輻射與吸收相同,且輻射系數(shù)jν與吸收系數(shù)αν之間存在如下關(guān)系:

式中Bν為黑體輻射系數(shù),

3.7.3 bb 輻射

bb 輻射損失功率密度為

式中h為普朗克常量,Nhp為全部重粒子(原子和離子)種類數(shù),lines 為粒子i的全部能級躍遷譜線,νi,21為重粒子i從高能級(2)躍遷到低能級(1)過程中輻射出的光子頻率,Ai,21為重粒子i從高能級到低能級的自發(fā)輻射的Einstein 系數(shù)(自發(fā)輻射幾率),ni,2為處于高能級的重粒子i數(shù)密度.

3.7.4 總輻射損失功率

綜上,總輻射損失功率為

不同壓強(qiáng)下Xe 等離子體輻射功率隨溫度變化曲線如圖5(d)所示.

4 結(jié)果與分析

4.1 計算域與邊界條件設(shè)置

本文研究基于連續(xù) CO2激光器(中心波長λ0=10.6 μm)與Xe 等離子體的相互作用,在柱坐標(biāo)系下建立模型,以激光入射等離子體方向為z軸,與其垂直方向為r軸.根據(jù)LSP 的形態(tài)特性,θ 方向各物理量近似均勻分布,因此可建立二維圓柱軸對稱模型,并采用COMSOL Multiphysics?軟件[58]進(jìn)行建模計算.如圖6 所示,設(shè)定圓柱形區(qū)域長度為 200 mm(z軸坐標(biāo)范圍為?100—100 mm),直徑為 70 mm(r軸坐標(biāo)范圍為 0—35 mm),激光從左側(cè)入射,在中心處聚焦.根據(jù)幾何光學(xué)假設(shè)對光束半徑的定義,確定了流體區(qū)域與激光區(qū)域的分界.計算區(qū)域劃分為矩形網(wǎng)格,與激光傳播方向平行,并在激光焦點處適當(dāng)加密,網(wǎng)格的最小尺寸約為 5 μm,總網(wǎng)格數(shù)約為 104個.采用變時間步長算法,相對誤差(relative tolerance)設(shè)置為 0.5%,對應(yīng)的時間步長約為 0.1 μs.

圖6 LSP 模擬計算域設(shè)置Fig.6.Setup of the LSP computational domain.

模型采用的控制方程及對應(yīng)邊界條件見表3,具體設(shè)置為: 計算區(qū)域的上邊界設(shè)定為無滑移壁面,即控制u=0,下邊界為對稱邊界; 考慮流場的作用,左側(cè)入口設(shè)置為速度邊界,且僅在軸向上存在一固定速度,徑向速度為 0,使得氣流方向與激光傳播方向一致; 右側(cè)出口處為壓強(qiáng)邊界,設(shè)置此處壓強(qiáng)大小為常數(shù); 傳熱場計算中,在對稱軸外的所有邊界上設(shè)置恒定溫度條件,邊界溫度為 500 K.為模擬起始時刻等離子體的產(chǎn)生,設(shè)置演化的初始溫度分布為中心位于激光焦點的二維高斯分布.

表3 模型控制方程及對應(yīng)邊界條件Table 3.Control equations and boundary conditions of the model.

4.2 穩(wěn)定狀態(tài)下溫度與速度分布

本文采用瞬態(tài)計算模型模擬LSP 的時空發(fā)展過程.當(dāng)觀察到核心位置及溫度分布在幾個 ms 內(nèi)均不再發(fā)生明顯變化時,判定LSP 的發(fā)展達(dá)到穩(wěn)態(tài),一般從激光注入至穩(wěn)態(tài)所需時間為10—15 ms.圖7 為穩(wěn)定狀態(tài)下LSP 的等溫線(圖7(a))與速度流線(圖7(b)).圖7—圖10 結(jié)果的計算區(qū)域與4.1 節(jié)中描述相同,并根據(jù)高溫核心的位置和尺寸進(jìn)行了局部放大.圖7 給出達(dá)到穩(wěn)態(tài)后演化時間為 10 ms 時的結(jié)果,計算條件為:壓強(qiáng)為5 atm(1 atm=1.013×105Pa),激光入射方向平行的氣體流速為 10 m/s,入射高斯激光總功率為 500 W,F數(shù)為f/7.從圖7(a)中顯示的溫度分布可以看到,等離子體核心處最高溫度超過 15000 K,穩(wěn)定核心位于靠近光束焦點偏上游位置.LSP 高溫區(qū)域形狀整體呈現(xiàn)“水滴形”,即靠近激光入口的頭部溫度梯度更高,而尾部溫度梯度較低,存在明顯的拖尾現(xiàn)象.

圖7 仿真得到典型的穩(wěn)態(tài)Xe-LSP 參數(shù)二維分布 (a)等溫線分布; (b)速度流場分布Fig.7.Two-dimensional (2D) profiles of the steady-state Xe-LSP from the simulation results: (a) Temperature isotherm profile; (b) velocity streamline.

圖7(b)則對稱展示了對應(yīng)位置的速度流線圖.與初始的軸向流速相比,LSP 核心的頭部前沿面處速度出現(xiàn)了徑向分量.這是由于等離子體頭部快速吸收激光能量導(dǎo)致該區(qū)域產(chǎn)生熱阻力,形成熱阻塞效應(yīng),迫使部分原本流向等離子體的軸向氣流改變方向,從周圍的冷區(qū)通過,路徑發(fā)生彎曲[59].

4.3 不同激光模式的影響

為探究不同激光模式對LSP 特性的影響,在模型中通過改變左側(cè)邊界處的激光功率分布函數(shù)設(shè)置了不同的入射激光模式,如高斯激光模式(TEM00)與環(huán)形激光模式(TEM01?).高斯模式激光光強(qiáng)I的徑向分布可表示為

其中,P是入射激光功率,w是高斯光束的特征半徑,I(w) 是中心光強(qiáng)的 1/e2(≈13.5%) .對于高斯光束,r越大,I(r) 越小.本文用到的環(huán)形模式激光光強(qiáng)徑向分布可表示為

對于環(huán)形模式,隨著r增大,I先增大后減小,在r=0 處光強(qiáng)為0,在r=處最強(qiáng).將r=2w處定義為激光半徑RL.在r∈(0,RL) 范圍內(nèi)高斯和環(huán)形激光的功率分別為總功率的 99.98% 和99.70%,定義r∈(0,RL) 范圍為激光區(qū).

對比圖8(a),(b)中同等 500 W 功率下兩種激光模式的LSP 溫度分布,可以看到由環(huán)形激光產(chǎn)生的LSP 高溫核心徑向尺寸更大、長寬比更小,整體溫度更高,且高溫核心區(qū)域更靠近光束焦點處,這一差異在Shi 等[16]的實驗中也已經(jīng)得到證實.以上現(xiàn)象可以從能量分布與吸收的過程給出解釋.如圖8(a)所示,與環(huán)形模式相比,高斯光束的激光強(qiáng)度分布更集中在中心對稱軸上,且相同功率下其強(qiáng)度幅值高于環(huán)形激光.因此在注入等離子體過程中,z軸附近的光線經(jīng)過相同距離后剩余功率更高,可以加熱更多的下游氣體,進(jìn)而形成的高溫核心具有更大的長寬比.與此同時,激光強(qiáng)度分布的形狀與穩(wěn)定狀態(tài)下的核心形狀具有相關(guān)性,如環(huán)形激光強(qiáng)度極值出現(xiàn)在r=的圓環(huán)上,導(dǎo)致在此模式下產(chǎn)生的LSP 高溫核心位置發(fā)生了徑向偏移,不在對稱軸上.

4.4 不同時刻的演化過程

為探究初始階段在激光作用下等離子體的演化過程,圖9 展示了達(dá)到穩(wěn)態(tài)前不同時刻激光束焦點附近的溫度分布變化.激光注入起始時,由于焦點 (z=0 mm)處激光強(qiáng)度最大,此處的初始等離子體具有最高溫度,因此會率先吸收較多的激光能量并在t=0.03 ms 時形成高溫核心.t=0.03—1.0 ms時,焦點處形成的高溫核心逐漸沿著激光入射方向(z軸)向上游移動.約t=1 ms 后,核心的位置和形態(tài)均達(dá)到相對穩(wěn)定的狀態(tài),頭部位于焦點上游約 6 mm 處,最高溫度可達(dá)到約 23000 K,說明激光能量得到了等離子體的有效吸收.

圖9 不同時刻LSP 的溫度分布圖.模擬條件為 5 atm 的Xe 環(huán)境下,入射激光功率為 500 W,F 數(shù)為 f/7,氣體流速為10 m/sFig.9.Temperature distributions of LSP at different times during laser propagation.Background gas: Xe,pressure:5 atm,laser power: 500 W,F-number:f/7,velocity: 10 m/s.

在之后的發(fā)展過程中(t=1.0—10.0 ms),由于軸向氣體流速的存在(本例中流速設(shè)置為 10 m/s),LSP 的尾部逐漸向后發(fā)展并不斷延長,使得LSP的整體長度和體積增加,但靠近焦點前側(cè)的頭部高溫核心體積及位置均已達(dá)到穩(wěn)定,說明此時已經(jīng)形成了小體積、穩(wěn)定的高溫度LSP 核心.

4.5 不穩(wěn)定性探究

在采用高功率環(huán)形激光、且流場速度較大的條件下,容易觀察到LSP 尾部的不穩(wěn)定波動現(xiàn)象.圖10 給出了不同條件下形成穩(wěn)態(tài)LSP 在不同時刻的溫度分布,其中圖10(a)—(c)為環(huán)形光束下產(chǎn)生的LSP 在3,4 和 5 ms 的形態(tài).可以觀察到環(huán)形激光條件下的LSP 形狀不規(guī)則,且隨時間存在著形態(tài)上的變化.Kotov 等[39]將LSP 實驗中觀察到的類似現(xiàn)象類比為層流振蕩火焰在燃燒過程中會出現(xiàn)的Kelvin-Helmholtz 不穩(wěn)定現(xiàn)象.此類不穩(wěn)定性通常發(fā)生在密度和速度不同的流體分界面上.由于LSP 內(nèi)部溫度遠(yuǎn)高于環(huán)境氣體,導(dǎo)致其密度遠(yuǎn)低于外部,因此LSP 與外部氣體組成了一個具有密度差的“雙流”系統(tǒng).當(dāng)入射激光功率較高且外部流場速度較大時,兩流體在切向(z軸方向)上存在明顯速度差,產(chǎn)生對流.此時,界面附近法向方向的微小擾動可引起溫度在時間和空間分布上的Kelvin-Helmholtz 波動.

圖10 高功率、高氣體流速下產(chǎn)生的不穩(wěn)定LSP.圖中結(jié)果的氣體壓強(qiáng)均為 10 atm,采用環(huán)形激光,功率為 5000 W,F 數(shù)為f/3.5.圖(a)—(c)展示了入口流速設(shè)置為 10 m/s 產(chǎn)生的LSP 在 3.6,4.0 和 4.4 ms 的二維溫度分布; 圖(d)—(f)展示了入口流速設(shè)置為 20 m/s 產(chǎn)生的LSP 在 3.6,4.0 和 4.4 ms 的二維溫度分布Fig.10.Unstable LSP generated under the condition of high laser power and high gas velocity.The gas pressure for all the results is set to 10 atm.The laser power is 5000 W and F-number equals to f/3.5.Panels (a)–(c) show the two-dimensional temperature distribution of the LSP when inlet velocity is set to 10 m/s at 3.6,4.0,and 4.4 ms .Panels (d)–(f) show the two-dimensional temperature distribution of the LSP when inlet velocity is set to 20 m/s at 3.6,4.0,and 4.4 ms.

5 結(jié)論

本文基于流體計算方法建立了激光維持等離子體多物理場耦合的含時模型,應(yīng)用該模型研究了激光維持等離子體的穩(wěn)態(tài)特性、能量注入機(jī)制、初始演化過程及流體不穩(wěn)定性.模擬結(jié)果表明,入射激光的模式、功率、F數(shù),以及氣體種類、壓強(qiáng)、流速等參數(shù)條件會對LSP 的穩(wěn)態(tài)性質(zhì)產(chǎn)生影響.在相同功率和F數(shù)條件下,高斯模式與環(huán)形模式激光產(chǎn)生的LSP 具有不同的形狀和位置.本文對LSP的形成過程進(jìn)行了探究,給出了激光注入至到達(dá)穩(wěn)態(tài)過程中不同時刻的LSP 溫度分布.由于激光功率密度在焦點處最高,在激光注入起始時很容易在焦點處率先產(chǎn)生LSP 高溫核心,在接下來的演化過程中,LSP 核心會沿著激光入射方向向上游移動,最終穩(wěn)定在能量吸收與耗散達(dá)到平衡的位置.此外,在高功率環(huán)形激光入射、且流速較大的條件下,可以在模擬結(jié)果中觀察到LSP 的不穩(wěn)定現(xiàn)象.這一現(xiàn)象在實驗中也已經(jīng)被報道,但具體的物理機(jī)制仍待進(jìn)一步研究.本文提出的LSP 瞬態(tài)模型可與實驗探究互為補(bǔ)充,彌補(bǔ)實驗因診斷難度大、參數(shù)調(diào)控范圍小對LSP 微觀機(jī)制探究的不足.后續(xù)工作將考慮等離子體對光線的折射作用以及等離子體對自身輻射能量的自吸收過程.同時將基于改進(jìn)模型繼續(xù)深入研究LSP 中產(chǎn)生的不穩(wěn)定波動現(xiàn)象及臨界條件,給出相關(guān)物理機(jī)制的合理解釋.