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

張東荷雨 劉金寶 付洋洋

(清華大學電機工程與應用電子技術系,北京 100084)

1 引言

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

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

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

圖1 LSP 實驗裝置結構示意圖Fig.1.Schematic diagram of the LSP experimental setup.

圖2 CO2 激光條件下LSP 中主要能量傳遞機制,包括: 逆韌致輻射(ff)吸收; 自由-自由(ff)/自由-束縛(fb)輻射; 束縛-束縛(bb)輻射; 熱擴散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]研發了低功率光源LDLS (laser-driven light source),可以提供輻射波長從深紫外(170 nm)到紅外波段(2.1 μm)的高強度光照.LDLS 光源采用近紅外連續半導體激光器作為激勵源,可產生溫度高于10000 K、微米尺寸的Xe LSP,在各波段的輻照強度均明顯高于傳統氣體放電光源.經測試,50 WLDLS 光源在紫外至可見光波段的光照亮度是75 W 短弧氙燈的5—10 倍[9].美國科磊公司(KLATencor)基于LSP 光源技術開發了商用的寬光譜等離子圖案晶圓缺陷檢測系統,如KLA-39xx、KLA-29xx 系列等 .此外,針對如何提高LSP 光源亮度與穩定性,科磊公司嘗試在腔體方案中引入反向渦流氣體腔[10]、漩渦軸[11]等設計,以改善高功率條件下封閉腔體內可能出現的局部過熱現象.雖然國外已經出現了較為成熟的商用LSP 光源設備,但由于LSP 涉及流場、溫度場和激光能量場等多物理場的耦合,且相關實驗診斷難度大、機理復雜度高,因此仍存在大量物理現象和底層物理機制尚未充分闡明.

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

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

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

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

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

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

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

3 計算模型

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

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

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

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

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

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

3.1 參數計算

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

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

其中,

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

圖3 不同壓強的平衡狀態下Xe 氣體組分Fig.3.Composition of Xe plasma at different pressures.

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

當激光頻率大于等離子體頻率時,激光可以進入等離子體并在等離子體中傳播.不考慮激光在等離子體中傳播產生的電磁場影響,基于幾何光學假設,根據Beer-Lambert 定律,沿激光傳播路徑s上激光強度I的變化可表示為

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

3.2.1 ff 吸收

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

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

3.2.2 bf 吸收

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

3.2.3 bb 吸收

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

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

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

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

3.2.4 總吸收系數

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

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

圖4 10.6 μm CO2 激光條件下不同壓強Xe 等離子體的吸收系數隨溫度變化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之間的關系.首先定義計算系統的質量密度ρ,即

其中,ni和mi分別為物種i的數密度和粒子質量.

焓H與配分函數等參數的關系為

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

定壓比熱容Cp與焓H、溫度T的關系為

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

圖5 不同壓強下Xe 等離子體各參數隨溫度的變化 (a)恒壓比熱容; (b)黏滯系數; (c)熱導率; (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 碰撞積分

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

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

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

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

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

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

3.4.1 指數排斥勢

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

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

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

3.4.2 蘭納-瓊斯勢

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

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

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

3.4.3 電荷交換相互作用

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

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

3.4.4 庫侖勢與屏蔽庫侖勢

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

其中,Zi為帶電粒子i的凈電荷數.

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

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

3.5 黏滯系數

根據Chapman 和Cowling[45]的理論,多組分系統的黏滯系數可表示為

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

式中,μij為二元黏性系數,其計算式為

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

3.6 熱導率

熱導率由幾部分共同組成,即

其中,kT,hp為重粒子貢獻的熱導率,kT,e為電子貢獻的熱導率.重粒子貢獻kT,hp表達式為

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

其中,Nhea為重粒子的種類數,二元熱導率kij可表示為

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

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

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

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

3.7 輻射功率

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

3.7.1 ff 輻射

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

有效凈電荷數Zeff定義式為[57]

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

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

3.7.2 fb 輻射

由Kirchhoff 定律,LTE 狀態下等離子體的輻射與吸收相同,且輻射系數jν與吸收系數αν之間存在如下關系:

式中Bν為黑體輻射系數,

3.7.3 bb 輻射

bb 輻射損失功率密度為

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

3.7.4 總輻射損失功率

綜上,總輻射損失功率為

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

4 結果與分析

4.1 計算域與邊界條件設置

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

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

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

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

4.2 穩定狀態下溫度與速度分布

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

圖7 仿真得到典型的穩態Xe-LSP 參數二維分布 (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)則對稱展示了對應位置的速度流線圖.與初始的軸向流速相比,LSP 核心的頭部前沿面處速度出現了徑向分量.這是由于等離子體頭部快速吸收激光能量導致該區域產生熱阻力,形成熱阻塞效應,迫使部分原本流向等離子體的軸向氣流改變方向,從周圍的冷區通過,路徑發生彎曲[59].

4.3 不同激光模式的影響

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

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

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

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

4.4 不同時刻的演化過程

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

圖9 不同時刻LSP 的溫度分布圖.模擬條件為 5 atm 的Xe 環境下,入射激光功率為 500 W,F 數為 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.

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

4.5 不穩定性探究

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

圖10 高功率、高氣體流速下產生的不穩定LSP.圖中結果的氣體壓強均為 10 atm,采用環形激光,功率為 5000 W,F 數為f/3.5.圖(a)—(c)展示了入口流速設置為 10 m/s 產生的LSP 在 3.6,4.0 和 4.4 ms 的二維溫度分布; 圖(d)—(f)展示了入口流速設置為 20 m/s 產生的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 結論

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