空气光学击穿的光线追踪模型*

黄 辉,洪延姬,李 倩,曹正蕊,冯海兵

（1.装备指挥技术学院基础部,北京 101416;

2.装备指挥技术学院研究生院,北京 101416）

1 引 言

将高能脉冲激光聚焦在空气中,当辐射强度超过一定的阈值时,会发生空气击穿现象,产生等离子体。自从1963年J.F.Ready首次报道激光导致空气击穿的现象后[1],由于该现象在高功率激光应用中的重要性,自此引起了广泛关注[2-7]。

在模拟激光点火时,能量瞬时沉积模型[2-5]中不考虑空气击穿和激光能量的沉积过程,直接把激光脉冲能量作为流体守恒方程的初始条件。龚平等[2]采用前期点爆炸自模拟解和后期高分辨率PPM格式相结合的方法;H.Yan等[3]采用1维球对称模型,初始时刻光斑区域的温度满足高斯分布;类似地,I.G.Dors等[4]考虑了激光能量沉积的非对称性和电离、离解效应对流场的影响,将等离子计算模型引入连续输运方程,数值计算得到的流场在100 μs时出现了涡环。R.Rozman等[6]模拟了激光诱导等离子形成过程,研究了逆韧制吸收、光致电离等各种吸收机制对等离子流场特性的影响。倪晓武等[7]结合激光等离子体的特性对强激光致空气击穿的过程进行数值模拟,但是其中描述光学击穿的电子数密度模型中仅考虑了电离和电子的损失,与实际物理过程差别比较大。

本文中将采用光线追踪法,建立简化的物理模型来研究空气击穿现象。空气被击穿后,温度瞬时升高到一定的值,入射激光能量通过逆韧致吸收机制沉积下来,进而完成随后的物理过程。对激光等离子体流场进行数值模拟,给出空气击穿后爆轰波形成和演化过程的物理图像。

2 物理模型和数值模拟方法

2.1 光路和光线传输模型

入射激光束通过透镜聚焦到焦点区域,光路如图1所示。由于光的衍射,光线的实际汇聚区域不是一个点,而是具有一定大小的光斑,叫做艾里（Airy）斑。艾里斑半径ra的表达式为

式中:λ为激光波长,f为聚焦透镜的焦距,D为最小孔径大小;β为衍射极限倍数,需要用实验来确定,一般来说,β的取为3～8,甚至更高。

图1 透镜聚焦系统Fig.1 Lens focusing system

在模拟激光束时使用几何光学,假设激光束由独立的光线组成,按照一定的方法将激光束离散成若干条光线,每条光线除了有确定的直线方程,还携带一定的功率。光线传播过程中,辐射强度的变化规律遵循辐射输运方程

式中:c为光速;I为辐射强度,对一条光线来说,就是它携带的功率;l是沿激光传输方向的传输距离;μv为气体对激光的吸收系数,由高温气体状态方程得到。

如果方程（2）中μv=0,即不考虑传播过程中能量的损失,艾里斑上的平均功率密度为

式中:ra为艾里斑的半径。数值求解方程（2）时,需要知道光线在网格中的传播路径。光线在网格中传播时,从1个网格单元进入下1个网格单元,有2种情形,如图2所示。图2（a）中,第k条光线穿过1条边进入（i,j）网格,图2（b）中,光线通过1个网格点进入网格。

图2 光线从1个网格进入下1个网格Fig.2 The ray goes from one cell into another cell

2.2 空气的击穿阈值

导致空气被激光击穿产生等离子体的机制,已被人们普遍接受的有2种,分别是多光子电离过程和级联吸收过程。根据N.Kroll等[8]的研究,空气的击穿阈值计算公式为

式中:击穿阈值Ibd的单位为W/cm2,波长λ单位为 μm,脉宽tp单位为s,气压p单位为 Pa。

纯净空气的击穿阈值对于CO2激光约为108～1010W/cm2,具体数值随脉冲宽度、焦斑直径等因素而异,当空气中含有气溶胶粒子时,击穿阈值将降低2～3个数量级。

当焦平面上的功率密度大于击穿阈值,即满足

时,焦平面会出现击穿现象,对应的时刻记为ta,式中Ia的意义参见式（3）。ta总是处于激光脉冲波形的上升前沿范围内。由式（5）仅能够判断焦平面上是否被击穿,不能对整个流场进行研究。

2.3 基于光线追踪法的击穿模型

要研究流场的击穿情况,先求出所有网格的功率密度。在图2（a）的情形下,光线k对（i,j）网格功率密度的贡献为

式中:Pk为光线进入（i,j）网格时携带的功率,SAB为（i,j）网格AijBij边的面积。

在图2（b）的情形下,光线k对（i,j）网格功率密度的贡献为

与式（6）一样,式（7）中的SBC为边BijCij的面积。（i,j）网格的功率密度是每条光线的贡献之和,即

为得到流场的功率密度分布,在击穿模型中需要2次根据式（8）计算网格的功率密度。后1次计算与前1次不一样,在沿着光线路径计算的过程中,如果碰到第1次所得功率密度大于击穿阈值的网格,就假设光线不再继续传播,认为光线被这个网格单元屏蔽掉。这样假设和物理过程相符,因为被击穿的区域,空气高度电离,光量子要被强烈地吸收,这些区域对光线来说,是不透明的。

如果流场中存在网格（i,j）,其功率密度Iij＞Ibd,则流场被激光击穿,对应的时刻记为tRtm,tRtm是对整个流场进行分析得到的结果。击穿模型的研究对象是网格,可以很方便地和辐射流体动力学计算耦合起来。此外,模型可以考虑非均匀介质的击穿,当介质非均匀时,各处的击穿阈值不相同,光线在非均匀介质中传播时,方向会发生改变,介质可能会吸收一部分的激光能量。

2.4 数值模拟条件和方法

计算区域为0.1 m×0.05 m的矩形,采用2维轴对称模型,除矩形的1条长边为对称轴边界条件外,其他3个边界为超声速出口边界条件,在构造网格时采用了矩形网格。激光束从左侧入射,光路是旋转轴对称的,对称轴为透镜的光轴,图1为光路沿轴线的剖面图,激光束半径R=0.05 m,透镜焦距f=0.25 m,焦点在流场中的坐标为（0.05 m,0 m）。由于计算模型是2维轴对称的,式（6）和式（7）中的面积S是柱坐标下,相应网格边绕轴线旋转1周得到的圆台的侧面积。激光脉冲波形对空气击穿过程有很大的影响,实验测得的波形拖尾段有很多振荡,为方便理论分析和数值计算的进行,将激光脉冲波形拟合为计算公式

式中:功率P（t）的单位为 W,时间t的单位是 s,常数因子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。

当激光强度变化不很剧烈,或者传输的特征距离不是很大时,激光穿过特征距离所需要的时间很短,方程（2）左边第1项可以被忽略,数值求解辐射疏运方程变得非常简单。得到的沉积能量作为源项耦合进入流体动力学方程组。流体计算使用有限体积法,选用Roe格式,并且通过使用具有最小模值通量限制器的MUSCL方法达到空间2阶精度,用预估-校正法进行时间积分。

3 结果与分析

3.1 击穿模型的验证

基于光线追踪法建立的击穿模型,计算功率密度时用到了网格边的面积,模型可能依赖于网格。为验证击穿模型的合理性,使用表1中的3种计算网格,对击穿模型进行验证。

表1 验证击穿模型使用的3种网格Table 1 Three computational grids for qualifying the breakdown model

表2中给出了针对不同的衍射极限倍数β,使用上述3种网格,击穿模型的计算结果,表中tRtm为击穿模型得到的击穿时刻,ta为求解式（5）得到的击穿时刻,n是艾里斑穿过的网格数,以网格3作为标准。表中只列出了击穿时刻,击穿位置总是焦平面处最靠近对称轴的那个网格单元。

比较表2中不同网格得到的击穿时刻可知,随着衍射极限倍数β的增加,网格对击穿时刻tRtm的影响越来越小。此外,当焦点处的网格尺寸比艾里斑半径小得多时,利用光线追踪法对整个流场进行分析得到的击穿时刻与直接判断焦平面击穿得到的结果相差比较小。文献[7]中模拟空气击穿,焦点处的网格也比较密,光斑半径跨越了10个网格点。

表2 不同网格和衍射极限的计算结果Table 2 Results for different grids and diffraction limits

当流场中有网格被击穿时,击穿网格的温度被赋予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中计算网格的典型尺寸为0.1～1 mm,由表3可知,当击穿区域的温度高于15000 K时,对应的辐射自由程等于或者小于网格典型尺寸,激光能量以约60%的比例沉积到焦点区域,转化成气体的能量。而当击穿区域温度为14 000 K或更低时,能量沉积效率为0。文献中已经说明吸气式激光推进中等离子体的温度为10 000～100 000 K[4,7,9],从量级上可以看出,本文所得结果与实际情况较吻合。

3.2 空气击穿后的流场演化

采用表1中的网格3,衍射极限倍数β=12,击穿区域的温度为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激光器,在大气条件下得到的实验阴影照片,入射激光束参数与数值模拟条件相同。对比数值计算结果和实验结果,发现两者所描述的过程基本一致;数值方法得到的等高线图中,迎着激光入射方向传播的激波不如实验照片中的尖锐,相对要圆滑一些。图片的显示方式是一个原因,计算程序是基于2维轴对称情形,绘制出的等高线是剖面图,而实验所拍摄的照片是3维立体火球的平面投影。

图4 激光击穿空气后的流场序列阴影照片（以空气击穿的时间为0时刻）Fig.4 Shadowing photographs of field after breakdown by taking breakdown time as zero time

4 结 论

采用光线追踪法,对空气击穿模型做了较为详细的分析和讨论。透镜聚焦情形下的计算结果表明,当焦平面处的艾里斑半径比该处网格尺寸大得多时,模型对计算网格的依赖性比较小,能够较好地描述流场的击穿现象。利用本文建立的击穿模型,得到了等离子流场爆轰波的形成和演化过程,计算结果与实验结果基本一致。

降低激光在击穿空气产生等离子体环节所消耗的能量,提高激光能量的利用效率有着积极的实际意义,本文中建立的击穿模型可为能量沉积的研究奠定基础,也可为辐射流体计算点火模型的确定提供参考及研究思路。此外,数值计算中考虑了高温真实气体效应,使击穿模型更能反映出等离子体的产生和膨胀等物理过程。

需要指出的是,本文模型中计算功率密度的方法存在一些局限性,当光线以大的入射角进入网格时,比如在2次反射的环聚焦情形下,这种计算方法误差比较大,下一步应进行适当改进,以减小模型的误差。

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

[2] 龚平,唐志平.大气呼吸模式激光推进的机理分析及数值模拟[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] 倪晓武,王文中,陆建,等.强激光致空气击穿过程的数值模拟[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] 泽尔道维奇,莱依捷尔.激波和高温流体动力学现象物理学:上册[M].张树材,译.北京:科技出版社,1980:260.