天津大学科学计算选讲第二次作业激光被动调Q速率方程组数值仿真学

天津大学科学计算选讲第二次作业激光被动调Q速率方程组数值仿真学生:徐世林学号:1011202008学院:精密仪器与光电子工程学院专业:光电子技术(博士)日期:2012年5月22日1激光被动调Q速率方程组数值仿真1.物理背景Q值是评定激光器中光学谐振腔质量好坏的指标,称为品质因数。Q值定义为在激光谐振腔内，储存的总能量与腔内单位时间损耗的能量之比,数学表达式为02WQvdWdt。式中W---腔内储存的总能量，dW/dt为光子能量的损耗速率，即单位时间内损耗的能量，0v为激光的中心频率。Q值是评定激光器中光学谐振腔质量好坏的指标，在激光谐振腔内储能不变情况下，损耗越小Q值越大。一般固体脉冲激光器的直接输出时，由于激光器工作阈值附近，输出脉冲是系列尖峰振荡，输出脉宽比较宽，输出功率低。激光调Q是一种广泛用于产生巨脉冲功率激光的运转方式，通过改变Q值(改变阈值)，控制激光产生时间，输出巨脉冲，压缩脉冲宽度，提高峰值功率。激光调节Q值的途径一般是采取改变腔内损耗的办法来调节腔内的Q值。调Q激光器输出的脉冲从时间上分为三阶段：脉冲建立时间、脉冲前沿时间和脉冲后沿时间。从Q值最小变到最大Q值即损耗从最大变到最小需要的时间叫开关时间，开关时间对激光脉冲的影响很大，按开关时间的大小分为快、慢两种类型。被动调Q是利用材料一种可饱和吸收特性做为Q开关，这种可饱和吸收染料是一种非线性吸收物质，把它放在谐振腔内，利用它对光的可饱和吸收特性来改变谐振腔内的吸收损耗，起到Q开关的作用。可饱和吸收特性指材料的透射特性随着能量密度的增大而变透明，当能量密度达到某值时，材料的吸收达到饱和，从而具有很高的透射率。吸收体对入射光强的改变应是这两个过程的综合过程。高能态粒子越多，则吸收光强越少。当高能态粒子数n2和基态粒子数n1相等时，这时从低能态到高能态的受激跃迁的粒子数和从高能态到基态放出光子的粒子数相等，即吸收的和放出的光子相等，达到饱和。常用的可饱和吸收材料有机染料，掺吸收性离子或含有色心（由点缺陷引起）的晶体。激光器速率方程组是表征激光器腔内光子数和工作物质各有关能级上的粒子数随时间变化的微分方程组。对于一般的脉冲激光器，脉冲形成时间长，泵浦、自发辐射、受激跃迁过程都是存在的。采用Q开关技术后，各参量之间关系发生了很大变化，需要根据Q开关的过程分析速率方程变化。由于被动调Q脉冲通常2只有纳秒量级，在脉冲期间可以不考虑抽运源的影响，研究连续脉冲的情况，需要考虑抽运速率、增益介质的上能级寿命、可饱和吸收体的恢复时间等因素，得到下面考虑了可饱和吸收体激发态吸收的连续抽运被动调Q速率方程组为1222(0)(ln)(1)0gsgssessgssrpTagssgsgsgsgsdnlnlnnldttRdnnnRcndtNdnnncndt（1）是激光腔中光子数密度；n是增益介质反转粒子数密度；ngs是可饱和吸收体基态粒子数密度，nes是可饱和吸收体激发态粒子数密度；n0s是可饱和I吸收体总粒子数密度;σ和l分别是增益介质的受激发射截面和长度；gs和σes分别是可饱和吸收体基态和激发态的吸收截面；ls是可饱和吸收体沿光腔轴线的长度；R是输出镜的反射率；是激光器腔体的耗散性损挺；c是光速；r为光在腔中往返一用的时间；是反转因子，对于四能级系统为1三能级系统为2；RP(t)是抽运速率；NT是增益介质的总粒子数密度；a是增益介质的上能级寿命:gs是可饱和吸收体的恢复时间；为了接近实际所添加的因子（1-n/NT）。2.MATLAB中的ode求解器一阶常微分方程的初值，其一般形式为00(,)()dyftydtyty，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB具有丰富的函数，我们将其统称为求解器（solver），其一般格式为：[T,Y]=solver(odefx,t,y0)。该函数表示在区间t=[t0,tf]上，用初始条件y0求解显式常微分方程。solver为命令ode45，ode23，ode113，ode15s，ode23s，ode23t，ode23tb之一，这些命令各有特点。我们这里主要考虑的是常用的ode45求解器，其特点：1.ode45是最常用求解微分方程数值解的命令，对于刚性方程组不宜采用。2.ode45求解器属于变步长的一种，采用Runge-Kutta算法；和他采用相同算法的变步长求解器还有ode23。33.ode45表示采用四阶，五阶Runge-Kutta单步算法,截断误差为(Δx)3，解决的是非刚性常微分方程。4.ode45是解决数值解问题的首选方法，若长时间没结果，应该就是刚性的，可换用ode23试试。ode45函数的调用格式为[T,Y]=de45(odefun,tspan,y0,options)，odefun是函数句柄,可以是函数文件名tspan是区间[t0tf]或者一系列散点[t0,t1,...,tf]，options是求解参数设置,可以用odeset在计算前设定误差,输出参数等。3.被动调Q速率方程组数值求解matlab中编程求解被动调Q的速率方程组，并作出光子数密度、反转粒子数密度以及基态粒子数密度随时间的变化的曲线，各物理量系数的典型取值下表。物理量典型值单位物理量典型值单位5.4×10-23m2gs8.7×10-23m2es2.2×10-23m2NT1.68×1026m3n11.82n21.800.02l1mma750µsls1mmgs3µs1速率方程组（M文件“rate_eq.m”）程序代码:functiondy=rate_eq(t,y,flag,Rp,T0,R)sigma=5.4e-23;%增益介质的受激发射截面sigma_gs=8.7e-23;%可饱和吸收体基态的吸收截面sigma_es=2.2e-23;%可饱和吸收体激发态的吸收截面N_T=1.68e26;%增益介质的总粒子数密度tao_a=750e-6;%增益介质的上能级寿命tao_gs=3e-6;%可饱和吸收体的恢复时间n1=1.82;%Nd3+:YAG的折射率n2=1.80;%Cr4+:YAG的折射率4delta=0.02;%激光器腔体的耗散性损耗l=0.001;%:f:曾1益介质的长度ls=0.001;%可饱和吸收体沿光腔轴线的长度gamma=1;%反转因子，对于四能级系统为1，三能级系统为2c=2.997963e8;%真空中的光速lc=n1*l+n2*ls;%谐振腔等效光程长度tr=lc/c;%光在腔中往返一周的时间n0s=-log(T0)/(sigma_gs*ls);%求可饱和吸收体粒子数密度y(1)=max(y(1),1);%光子数密度的最小值%被动调Q稠合方程组:dy=[y(1)*(2*sigma*y(2)*l-2*sigma_gs*y(3)*ls-2*sigma_es*...(n0s-y(3))*ls-(log(1/R)+delta))/tr;Rp*(1-y(2)/N_T)-gamma*sigma*c*y(1)*y(2)-y(2)/tao_a;(n0s-y(3))/tao_gs-sigma_gs*c*y(1)*y(3)];然后，在MATLAB中调用求解器ode45求解速率方程组（M文件“Qswitch.m”）程序代码:%被动调Q速率方程数值求解clc%Clearcommandwindowclear%ClearvariablesandfunctìonsfrommemorycloseallT0=0.7;%可饱和吸收体初始透射率R=0.8;%输出镜反射率Rp=2e28;%抽运速率y0=[1;0;0];%设定初值tspan=[00.05];%设定计算时间范围tic[t,y]=ode45('rate_eq',tspan,y0,[],Rp,T0,R);%解速率方程组(ode45常微分方程的数值求解)tocy(:,1)=max(y(:,1),1);figuresubplot(3,1,1);plot(t,y(:,3));xlabel(‘时间(s)');ylabel('基态粒子数密度(m^{-3})');figure%将光子数密度和反转粒子数密度随时间变化画于同一图中[AX,H1,H2]=plotyy(t,y(:,1),t,y(:,2));set(H2,'LineStyle','--')xlabel(‘时间(s)');set(get(AX(1),‘Ylabel’),‘String‘,‘光子数密度(m^{-3})')set(get(AX(2)，'Ylabel')，'String'，'反转粒子数密度(m^{-3})')程序运行可得被动调Q速率方程组的数值求解结果500.010.020.030.040.050510x1023时间(s)光子数密度(m-3)00.010.020.030.040.05012x1025时间(s)反转粒子数密度(m-3)00.010.020.030.040.0505x1024时间(s)基态粒子数密度(m-3)图一、光子数密度、反转粒子数密度以及基态粒子数密度随时间变化00.010.020.030.040.050510x1023时间(s)光子数密度(m-3)00.010.020.030.040.05012x1025反转粒子数密度(m-3)图二、光子数密度和增益介质反转粒子数密度随时间变化MATLAB中调用plotyy函数，将光子数密度和增益介质反转粒子数密度随时间变化在同一张图中作出，可以看出在被动调Q过程中，每产生一个光脉冲，增益介质的反转粒子数密度就会急剧减少。