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

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

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

1 引 言

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

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

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

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

2.1 光路和光線傳輸模型

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

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

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

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

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

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

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

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

2.2 空氣的擊穿閾值

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)果與分析

3.1 擊穿模型的驗(yàn)證

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

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

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

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

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

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

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

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

3.2 空氣擊穿后的流場(chǎng)演化

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

圖3 流場(chǎng)溫度等高線圖Fig.3 Temperature contours of the field

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

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

4 結(jié) 論

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

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

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

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

[2] 龔平,唐志平.大氣呼吸模式激光推進(jìn)的機(jī)理分析及數(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] 倪曉武,王文中,陸建,等.強(qiáng)激光致空氣擊穿過(guò)程的數(shù)值模擬[J].兵工學(xué)報(bào),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] 澤爾道維奇,萊依捷爾.激波和高溫流體動(dòng)力學(xué)現(xiàn)象物理學(xué):上冊(cè)[M].張樹材,譯.北京:科技出版社,1980:260.