黃 輝,洪延姬,李 倩,曹正蕊,馮海兵

（1.裝備指揮技術(shù)學院基礎部,北京 101416;

2.裝備指揮技術(shù)學院研究生院,北京 101416）

1 引 言

將高能脈沖激光聚焦在空氣中,當輻射強度超過一定的閾值時,會發(fā)生空氣擊穿現(xiàn)象,產(chǎn)生等離子體。自從1963年J.F.Ready首次報道激光導致空氣擊穿的現(xiàn)象后[1],由于該現(xiàn)象在高功率激光應用中的重要性,自此引起了廣泛關(guān)注[2-7]。

在模擬激光點火時,能量瞬時沉積模型[2-5]中不考慮空氣擊穿和激光能量的沉積過程,直接把激光脈沖能量作為流體守恒方程的初始條件。龔平等[2]采用前期點爆炸自模擬解和后期高分辨率PPM格式相結(jié)合的方法;H.Yan等[3]采用1維球?qū)ΨQ模型,初始時刻光斑區(qū)域的溫度滿足高斯分布;類似地,I.G.Dors等[4]考慮了激光能量沉積的非對稱性和電離、離解效應對流場的影響,將等離子計算模型引入連續(xù)輸運方程,數(shù)值計算得到的流場在100 μs時出現(xiàn)了渦環(huán)。R.Rozman等[6]模擬了激光誘導等離子形成過程,研究了逆韌制吸收、光致電離等各種吸收機制對等離子流場特性的影響。倪曉武等[7]結(jié)合激光等離子體的特性對強激光致空氣擊穿的過程進行數(shù)值模擬,但是其中描述光學擊穿的電子數(shù)密度模型中僅考慮了電離和電子的損失,與實際物理過程差別比較大。

本文中將采用光線追蹤法,建立簡化的物理模型來研究空氣擊穿現(xiàn)象?？諝獗粨舸┖?溫度瞬時升高到一定的值,入射激光能量通過逆韌致吸收機制沉積下來,進而完成隨后的物理過程。對激光等離子體流場進行數(shù)值模擬,給出空氣擊穿后爆轟波形成和演化過程的物理圖像。

2 物理模型和數(shù)值模擬方法

2.1 光路和光線傳輸模型

入射激光束通過透鏡聚焦到焦點區(qū)域,光路如圖1所示。由于光的衍射,光線的實際匯聚區(qū)域不是一個點,而是具有一定大小的光斑,叫做艾里（Airy）斑。艾里斑半徑ra的表達式為

式中:λ為激光波長,f為聚焦透鏡的焦距,D為最小孔徑大小;β為衍射極限倍數(shù),需要用實驗來確定,一般來說,β的取為3～8,甚至更高。

圖1 透鏡聚焦系統(tǒng)Fig.1 Lens focusing system

在模擬激光束時使用幾何光學,假設激光束由獨立的光線組成,按照一定的方法將激光束離散成若干條光線,每條光線除了有確定的直線方程,還攜帶一定的功率。光線傳播過程中,輻射強度的變化規(guī)律遵循輻射輸運方程

式中:c為光速;I為輻射強度,對一條光線來說,就是它攜帶的功率;l是沿激光傳輸方向的傳輸距離;μv為氣體對激光的吸收系數(shù),由高溫氣體狀態(tài)方程得到。

如果方程（2）中μv=0,即不考慮傳播過程中能量的損失,艾里斑上的平均功率密度為

式中:ra為艾里斑的半徑。數(shù)值求解方程（2）時,需要知道光線在網(wǎng)格中的傳播路徑。光線在網(wǎng)格中傳播時,從1個網(wǎng)格單元進入下1個網(wǎng)格單元,有2種情形,如圖2所示。圖2（a）中,第k條光線穿過1條邊進入（i,j）網(wǎng)格,圖2（b）中,光線通過1個網(wǎng)格點進入網(wǎng)格。

圖2 光線從1個網(wǎng)格進入下1個網(wǎng)格Fig.2 The ray goes from one cell into another cell

2.2 空氣的擊穿閾值

導致空氣被激光擊穿產(chǎn)生等離子體的機制,已被人們普遍接受的有2種,分別是多光子電離過程和級聯(lián)吸收過程。根據(jù)N.Kroll等[8]的研究,空氣的擊穿閾值計算公式為

式中:擊穿閾值Ibd的單位為W/cm2,波長λ單位為 μm,脈寬tp單位為s,氣壓p單位為 Pa。

純凈空氣的擊穿閾值對于CO2激光約為108～1010W/cm2,具體數(shù)值隨脈沖寬度、焦斑直徑等因素而異,當空氣中含有氣溶膠粒子時,擊穿閾值將降低2～3個數(shù)量級。

當焦平面上的功率密度大于擊穿閾值,即滿足

時,焦平面會出現(xiàn)擊穿現(xiàn)象,對應的時刻記為ta,式中Ia的意義參見式（3）。ta總是處于激光脈沖波形的上升前沿范圍內(nèi)。由式（5）僅能夠判斷焦平面上是否被擊穿,不能對整個流場進行研究。

2.3 基于光線追蹤法的擊穿模型

要研究流場的擊穿情況,先求出所有網(wǎng)格的功率密度。在圖2（a）的情形下,光線k對（i,j）網(wǎng)格功率密度的貢獻為

式中:Pk為光線進入（i,j）網(wǎng)格時攜帶的功率,SAB為（i,j）網(wǎng)格AijBij邊的面積。

在圖2（b）的情形下,光線k對（i,j）網(wǎng)格功率密度的貢獻為

與式（6）一樣,式（7）中的SBC為邊BijCij的面積。（i,j）網(wǎng)格的功率密度是每條光線的貢獻之和,即

為得到流場的功率密度分布,在擊穿模型中需要2次根據(jù)式（8）計算網(wǎng)格的功率密度。后1次計算與前1次不一樣,在沿著光線路徑計算的過程中,如果碰到第1次所得功率密度大于擊穿閾值的網(wǎng)格,就假設光線不再繼續(xù)傳播,認為光線被這個網(wǎng)格單元屏蔽掉。這樣假設和物理過程相符,因為被擊穿的區(qū)域,空氣高度電離,光量子要被強烈地吸收,這些區(qū)域?qū)饩€來說,是不透明的。

如果流場中存在網(wǎng)格（i,j）,其功率密度Iij＞Ibd,則流場被激光擊穿,對應的時刻記為tRtm,tRtm是對整個流場進行分析得到的結(jié)果。擊穿模型的研究對象是網(wǎng)格,可以很方便地和輻射流體動力學計算耦合起來。此外,模型可以考慮非均勻介質(zhì)的擊穿,當介質(zhì)非均勻時,各處的擊穿閾值不相同,光線在非均勻介質(zhì)中傳播時,方向會發(fā)生改變,介質(zhì)可能會吸收一部分的激光能量。

2.4 數(shù)值模擬條件和方法

計算區(qū)域為0.1 m×0.05 m的矩形,采用2維軸對稱模型,除矩形的1條長邊為對稱軸邊界條件外,其他3個邊界為超聲速出口邊界條件,在構(gòu)造網(wǎng)格時采用了矩形網(wǎng)格。激光束從左側(cè)入射,光路是旋轉(zhuǎn)軸對稱的,對稱軸為透鏡的光軸,圖1為光路沿軸線的剖面圖,激光束半徑R=0.05 m,透鏡焦距f=0.25 m,焦點在流場中的坐標為（0.05 m,0 m）。由于計算模型是2維軸對稱的,式（6）和式（7）中的面積S是柱坐標下,相應網(wǎng)格邊繞軸線旋轉(zhuǎn)1周得到的圓臺的側(cè)面積。激光脈沖波形對空氣擊穿過程有很大的影響,實驗測得的波形拖尾段有很多振蕩,為方便理論分析和數(shù)值計算的進行,將激光脈沖波形擬合為計算公式

式中:功率P（t）的單位為 W,時間t的單位是 s,常數(shù)因子a=0.366 666,b=-7.041 07×10-7,c=3.666 66×10-13;對應的單脈沖能量E0=32.5 J,脈沖寬度tp=10 μs,峰值功率時刻th=1 μs。標準大氣壓下,應用公式（4）,得到空氣的擊穿閾值約為1.43×109W/cm2。

當激光強度變化不很劇烈,或者傳輸?shù)奶卣骶嚯x不是很大時,激光穿過特征距離所需要的時間很短,方程（2）左邊第1項可以被忽略,數(shù)值求解輻射疏運方程變得非常簡單。得到的沉積能量作為源項耦合進入流體動力學方程組。流體計算使用有限體積法,選用Roe格式,并且通過使用具有最小模值通量限制器的MUSCL方法達到空間2階精度,用預估-校正法進行時間積分。

3 結(jié)果與分析

3.1 擊穿模型的驗證

基于光線追蹤法建立的擊穿模型,計算功率密度時用到了網(wǎng)格邊的面積,模型可能依賴于網(wǎng)格。為驗證擊穿模型的合理性,使用表1中的3種計算網(wǎng)格,對擊穿模型進行驗證。

表1 驗證擊穿模型使用的3種網(wǎng)格Table 1 Three computational grids for qualifying the breakdown model

表2中給出了針對不同的衍射極限倍數(shù)β,使用上述3種網(wǎng)格,擊穿模型的計算結(jié)果,表中tRtm為擊穿模型得到的擊穿時刻,ta為求解式（5）得到的擊穿時刻,n是艾里斑穿過的網(wǎng)格數(shù),以網(wǎng)格3作為標準。表中只列出了擊穿時刻,擊穿位置總是焦平面處最靠近對稱軸的那個網(wǎng)格單元。

比較表2中不同網(wǎng)格得到的擊穿時刻可知,隨著衍射極限倍數(shù)β的增加,網(wǎng)格對擊穿時刻tRtm的影響越來越小。此外,當焦點處的網(wǎng)格尺寸比艾里斑半徑小得多時,利用光線追蹤法對整個流場進行分析得到的擊穿時刻與直接判斷焦平面擊穿得到的結(jié)果相差比較小。文獻[7]中模擬空氣擊穿,焦點處的網(wǎng)格也比較密,光斑半徑跨越了10個網(wǎng)格點。

表2 不同網(wǎng)格和衍射極限的計算結(jié)果Table 2 Results for different grids and diffraction limits

當流場中有網(wǎng)格被擊穿時,擊穿網(wǎng)格的溫度被賦予1個給定的值,該溫度不能夠太低,否則,入射激光能量不能夠沉積下來。表2給出了常溫常壓下,按照Y.B.Zeldovich等[9]的理論,各種溫度對應的輻射自由程lf,僅考慮逆韌致吸收。

表3 密度為1.29 kg/m3的情況下不同溫度對應的輻射自由程Table 3 Radiation free path for different temperatures andρ=1.29 kg/m3

表1中計算網(wǎng)格的典型尺寸為0.1～1 mm,由表3可知,當擊穿區(qū)域的溫度高于15000 K時,對應的輻射自由程等于或者小于網(wǎng)格典型尺寸,激光能量以約60%的比例沉積到焦點區(qū)域,轉(zhuǎn)化成氣體的能量。而當擊穿區(qū)域溫度為14 000 K或更低時,能量沉積效率為0。文獻中已經(jīng)說明吸氣式激光推進中等離子體的溫度為10 000～100 000 K[4,7,9],從量級上可以看出,本文所得結(jié)果與實際情況較吻合。

3.2 空氣擊穿后的流場演化

采用表1中的網(wǎng)格3,衍射極限倍數(shù)β=12,擊穿區(qū)域的溫度為15000 K。由表2,t=0.44 μs時刻,空氣被擊穿,圖3為擊穿形成的等離子,在t=0.94、1.44和2.44 μs時刻流場溫度等高線圖,這些等高線較細致地描述了等離子體流場的演化過程。

圖3 流場溫度等高線圖Fig.3 Temperature contours of the field

圖4是使用高功率T EA脈沖CO2激光器,在大氣條件下得到的實驗陰影照片,入射激光束參數(shù)與數(shù)值模擬條件相同。對比數(shù)值計算結(jié)果和實驗結(jié)果,發(fā)現(xiàn)兩者所描述的過程基本一致;數(shù)值方法得到的等高線圖中,迎著激光入射方向傳播的激波不如實驗照片中的尖銳,相對要圓滑一些。圖片的顯示方式是一個原因,計算程序是基于2維軸對稱情形,繪制出的等高線是剖面圖,而實驗所拍攝的照片是3維立體火球的平面投影。

圖4 激光擊穿空氣后的流場序列陰影照片（以空氣擊穿的時間為0時刻）Fig.4 Shadowing photographs of field after breakdown by taking breakdown time as zero time

4 結(jié) 論

采用光線追蹤法,對空氣擊穿模型做了較為詳細的分析和討論。透鏡聚焦情形下的計算結(jié)果表明,當焦平面處的艾里斑半徑比該處網(wǎng)格尺寸大得多時,模型對計算網(wǎng)格的依賴性比較小,能夠較好地描述流場的擊穿現(xiàn)象。利用本文建立的擊穿模型,得到了等離子流場爆轟波的形成和演化過程,計算結(jié)果與實驗結(jié)果基本一致。

降低激光在擊穿空氣產(chǎn)生等離子體環(huán)節(jié)所消耗的能量,提高激光能量的利用效率有著積極的實際意義,本文中建立的擊穿模型可為能量沉積的研究奠定基礎,也可為輻射流體計算點火模型的確定提供參考及研究思路。此外,數(shù)值計算中考慮了高溫真實氣體效應,使擊穿模型更能反映出等離子體的產(chǎn)生和膨脹等物理過程。

需要指出的是,本文模型中計算功率密度的方法存在一些局限性,當光線以大的入射角進入網(wǎng)格時,比如在2次反射的環(huán)聚焦情形下,這種計算方法誤差比較大,下一步應進行適當改進,以減小模型的誤差。

[1] Ready J F.Effect of High-power Laser Radiation[M].New York:Academic Press Inc,1971:12-57.

[2] 龔平,唐志平.大氣呼吸模式激光推進的機理分析及數(shù)值模擬[J].爆炸與沖擊,2003,23（6）:501-508.GONG Ping,TANG Zhi-ping.Mechanism analysis and numerical simulation of air-breathing model for laser propulsion[J].Explosion and Shock Waves,2003,23（6）:501-508.

[3] Yan H,Adelgren R,Boguszko M,et al.Laser energy deposition in quiescent air[J].AIAA Journal,2003,41（10）:1988-1995.

[4] Dors I G,Parigger C G,Lewis J W L.Fluid dynamics effects following laser-induced optical breakdown[C]//Proceedings of the 38th Aerospace Sciences Meeting and Exhibit.AIAA 2000-0717,2000.

[5] Phuoc T X.An experimental and numerical study of laser-induced spark in air[J].Optics and Lasers in Egnineering,2005,43（2）:112-129.

[6] Rozman R,Grabec I,Govekar E.Influence of absorption mechanisms on laser-induced plasma plume[J].Applied Surface Science,2008,254（11）:3295-3305.

[7] 倪曉武,王文中,陸建,等.強激光致空氣擊穿過程的數(shù)值模擬[J].兵工學報,1998,19（2）:134-138.

NI Xiao-wu,WANG Wen-zhong,LU Jian,et al.Numerical simulation of laser-induced air plasma shock wave[J].Acta Armamentarii,1998,19（2）:134-138.

[8] Kroll N,Watson K M.Theoretical study of ionization of air by intense laser pulses[J].Physical Review A,1972,5（4）:1883-1889.

[9] 澤爾道維奇,萊依捷爾.激波和高溫流體動力學現(xiàn)象物理學:上冊[M].張樹材,譯.北京:科技出版社,1980:260.