《武汉工程大学学报》 2011年11期
95-98
出版日期：2011-11-30
ISSN:1674-2869
CN:42-1779/TQ
激光等离子体羽辉膨胀的流体
动力学模型及数值模拟
0引言激光烧蚀固体靶形成的等离子体羽辉的膨胀特性在脉冲激光沉积薄膜制备［1］、激光等离子体波导［2］、极紫外光刻机光源等离子体碎屑抑制［3］和等离子体光谱［4］等方面有着重要的应用.一般情况下，激光烧蚀形成的等离子体羽辉的空间分布是随时间变化的，对羽辉的时空演化规律主要采用粒子模拟［5］和流体动力学模拟［6］的方法，数值求解过程都较为复杂.本文从流体动力学方程组出发，假设烧蚀形成的等离子体羽辉的空间分布不随时间变化，将偏微分方程组简化为常微分方程组，最后数值求解得到了稳态烧蚀的激光等离子体羽辉的温度、速度和密度的空间分布规律，数值计算得出的规律与实验观测结果吻合的较好.1理论模型当高强度激光辐照固体靶材表面时，表面形成的高温高密的等离子体羽辉会迅速向外膨胀.为了描述激光烧蚀形成的等离子体，考虑到等离子体的德拜长度λD=6.9T/Ne cm，其中T是以K为单位的等离子体电子温度，Ne是以cm-3为单位的等离子体电子数密度［7］.对于CO2激光诱导产生的等离子体，如果电子温度T=105 K，临界电子密度Ne=1019 cm-3，计算出德拜长度为nm量级.我们所研究的激光等离子体系统的特征尺寸一般在mm量级，远大于等离子体的德拜长度，因此等离子体整体显示出电中性，并且可以忽略等离子体的微观运动行为而只需考虑等离子体的集体行为，也就是说，可以采用流体力学描述替代粒子描述.在欧拉描述下，可压缩流体动力学的方程为［8］ρ(v·)v+P+ρvt=0(1)
·(12v2+ε)v-·(βT)+t(12v2+ε)=0(2)·(ρ v)+ρt=0(3)
其中ρ是等离子体质量密度，v是等离子体流体速度，P是等离子体压强，ε是等离子体热力学能，β是等离子体的热导率.在激光等离子体中，由于离子的热导率远小于电子的热导率，离子的温度Ti远小于电子的温度T，因此可以忽略离子的影响.考虑到等离子体的平均电离度Z后，等离子体的行为将主要由电子决定.等离子体的热力学能可写为ε=3/2NekBT，将热力能的表达式代入(2)式，联立(1)式，能量守恒方程可以改写为
(v·)T+Tt+23T(·v)-βNe·(T2.5T)=0(4)考虑烧蚀的激光等离子体羽辉处于稳态，即：各物理量的空间分布不随时间变化，并假设等离子体的平均电离度Z也不随时间变化，羽辉可以看作理想气体，将ρ=Nemi/Z，P=NekBT代入动量守恒方程(1)式得到
miv·12v2+ZkB［T·v-T(·v)］=0(5)
其中kB为玻尔兹曼常数，mi为等离子体中离子的质量.假设激光能量主要在临界面附近被等离子体吸收，并且临界面的烧蚀速度由当地声速决定，即：v0=cs=ZkBT0/mi，其中T0为临界面的等离子体温度.为了方便起见，令u=vv0，τ=TT0，则(5)式可归一化为12u·(u2)+τ·u-τ(·u)=0(6)令n=Ne/N0，N0为激光等离子体临界面电子密度，并考虑在稳态情况下，则(4)式可归一化为
n u·τ+23nτ(·u)-T02.5N0cs·(βτ2.5τ)=0(7)
同理，在稳态假设下，式(3)可归一化为·(n u)=0(8)假设等离子体膨胀过程具有对称性，则三维的等离子体膨胀动力学方程组(6),(7)和(8)式可以写为12urkd(u2)dr-τd(rku)dr+urkdτdr=0(9)
nudτdr+23rknτd(rku)dr-βT02.5N0csddrτ2.5dτdr=0(10)d(rknu)dr=0(11)
其中k=0，1和2分别对应于平面对称、柱对称和球对称的情况.求解(9)式容易得到rknu=rk0n0u0(12)
这里我们假设r0、n0和u0分别为临界面处的位置、密度和速度，则n0=u0=1.不失一般性，我们定义归一化的位置坐标ζ=lnrr0，则ddr=e-ζr0ddζ，此时临界面的位置坐标就移到归一化位置坐标的原点.采用归一化坐标描述，并且联立式(12)，则(9)、(10)式可以化为ududζ-τudvdζ+dτdζ-kτ=0(13)
dτdζ+23τk+1ududζ-Cddζτ2.5e(k-1)ζdτdζ=0(14)
其中C=βT2.50N0r0ν0，非线性常微分方程组(13)、(14)的边界条件为u｜ζ=0=τ｜ζ=0=1，τ｜ζ=∞=0.第11期吴涛，等：激光等离子体羽辉膨胀的流体动力学模型及数值模拟
武汉工程大学学报第33卷
2数值求解与讨论下面数值求解非线性常微分方程组(13)、(14)，方便起见，取常数C=1，并对该方程组采用中心差分离散化后得到
12h(ui+1+ui)-τi+1ui+1+τiui(ui+1-ui)-k2(τi+1+τi)+τi+1-τih=0(15)
1h(τi+1-τi)+13(τi+1+τi)［k+12h1ui+1+1ui(ui+1-ui)］-13(e(k-1)ζi+1+e(k-1)ζi)×［k-12h(τ2.5i+1+τ2.5i)(τi+1-τi)］-13e(k-1)ζi+1+e(k-1)ζi［12h2(τ2.5i+1+τ2.5i)(τi+2-τi+1-τi+τi-1)］-13e(k-1)ζi+1+e(k-1)ζi［54h2(τ1.5i+1+τ1.5i)(τi+1-τi)2］=0(16)
其中h为空间步长，ζi=i×h，ui=u(ζi)和τi=τ(ζi)，i=0，1，2，3…取整数.如上所述，边界条件为u0=τ0=1，代数方程组(15)和(16)式的根即为微分方程的数值解，我们采用牛顿迭代法来求解该非线性差分方程组，这里需要给出u和τ的初始值，而带状的雅克比矩阵的阶数由网格点的个数决定，可以采用高斯消去法求解，迭代的次数与初始值的设定有关，只要设定较合适的初始值，一般迭代五次就能满足精度要求.考虑到激光等离子体的尺度特征，取r0=20 μm，r=3 mm，则0<ζ<5，步长h=0.01，即500个网格点.利用计算机语言编写程序数值求解，图1，图2和图3分别给出了激光等离子体羽辉温度和膨胀速度的数值求解结果，并与实验数据和其它的理论模型计算结果作了对比［911］.图1归一化等离子体温度τ随归一化坐标ζ的变化
关系曲线
Fig.1Normalized plasma temperature τ as a function of
normalized position ζ图2归一化等离子体膨胀速度u随归一化坐标ζ的
变化关系曲线
Fig.2Normalized plasma plume expansion velocity u as
a function of normalized position ζ值得指出的是，对于激光烧蚀平板靶材或者液滴靶材形成的等离子体羽辉的膨胀过程是介于柱对称和球对称之间的，对于平面对称k=0的情况，方程给出的是平凡解.另外，考虑到(12)式及归一化坐标ζ的定义，容易得到n(ζ)=e-kζv(ζ)，图3给出了归一化等离子体电子数密度随归一化坐标的变化规律.图3归一化等离子体密度n随归一化坐标ζ的变化
关系曲线
Fig.3Normalized plasma denstiy n as a function of
normalized position ζ从图1可以看出，归一化的等离子体温度在临界面处最高，随着空间距离的增大而减小直至零，且随k值的增大减小的更快.从图2可以看出，归一化的等离子体膨胀速度随距离的增大而增加，最后趋于一个稳定的渐进值.从图3可以看出，归一化的等离子体电子数密度随距离的增大而减小直至趋于零.实验结果和其它理论模型给出的计算结果也画在图中，结果表明：数值模拟得到的稳态烧蚀膨胀的等离子体羽辉的温度、速度和数密度的变化规律与文献［911］中实验测量结果及Anisimov的理论模型计算结果吻合的较好.3结语根据流体动力学的质量连续性、动量守恒和能量守恒这三个基本方程，建立了激光烧蚀靶材等离子体羽辉稳态膨胀的动力学模型，并根据烧蚀形成的等离子体的特征参量对方程组进行归一化处理，得到了一组非线性的常微分方程.采用有限差分和牛顿迭代法，利用计算机语言编程对方程组进行数值求解，得到了激光等离子体羽辉的温度、速度和密度空间分布的一般规律，得到的数值结果与实验观测结论是一致的.通过物理量的归一化处理，该理论模型将激光等离子体羽辉膨胀的过程标准化，方便和实验数据的对比分析，有助于理解激光等离子体羽辉膨胀的动力学过程.参考文献：