| 查看: 896 | 回复: 6 | |||
| 当前主题已经存档。 | |||
| 【有奖交流】积极回复本帖子，参与交流，就有机会分得作者 chenfeng79 的 3 个金币 | |||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | |||
[交流]
【求助】一个脉冲热效应的matlab程序不能运行
|
|||
|
本人编写了一个脉冲热效应的matlab 程序，显示结果为边界值，也就是主程序不能运行，即 if=1后面出现了问题，恳求有关高手给予指点。程序如下 lear clc % ************************************* a=1.5;b=5;%晶体参数 h=0.1;ht=5e-4;%计算步长 tol=0.001;maxl=3000;%***误差 kr=0.014;%热导率 absorb=0.91;%晶体吸收系数 w0=0.32;%泵光半径 Q=15;%泵浦功率 dndt=7.3e-6;%dn/dt n0=1.82;%折射率 o=7.5e-6;%热膨胀系数 L=30;%腔长 R=50;%输出镜半径 TOUT=0.1;%输入镜透过率 ls=946e-6;lp=808e-6;lf=1062e-6;%激光波长，泵浦波长，平均荧光波长 oabs=7.9e-24;oem=4e-24;%受激吸收截面，受激发射截面 t=230e-6;%上能级寿命 tl=0.025;%脉冲持续时间 rzxs=4e-4; % when adopt inter must follow inter<=f, inter=f maybe the best c=3e11;%光速 hc=6.62e-34;%普朗克常数 n=fix(a/h)+1;m=fix(b/h)+1;nt=fix(tl/ht)+1; % importent nt=fix(tl/ht)+1; otherwise the number will decrease inter%计算取用格子数 q=zeros(2*n-1,m); % ******************************** h0=0.006;%对流换热系数 %********************************* % 边界条件 Es=hc*c/ls;Ep=hc*c/lp;%光子能量 K=0.006;%光光转化效率 fa=0.074;%玻耳兹曼因子 lr=0.015;%腔内损耗 ql=4.22; %晶体密度 cl=0.59;%晶体比热 ws=((ls/pi)^2*L*(R-L))^(0.25);%激光半径公式 z0=n0*w0^2*pi/lp; for j=1:m wp(j)=w0*(1+j*h/z0)^(0.5);%泵浦光斑半径公式 pp(j)=Q*exp(-absorb*j*h);%泵浦光传播功率 ps=Q*K/TOUT*exp(0.5);%激光传播功率 end %程序主体 U=ones(2*n-1,m,nt); U0=16; U(1,1:m,1:nt)=16;U(2*n-1,1:m,1:nt)=16;U(1:2*n-1,1,1:nt)=16;U(1:2*n-1,m,1:nt)=16;U(1:2*n-1,1:m,1)=16; ****************************************************** for T=2:nt for i=2:2*n-2 for j=1:m q(i,j)=2*Q*(1-(lp/lf+(lp/ls-lp/lf)*... (ps*exp(-2*(wp(j)/2)^2/ws^2)*oem/Es)/(1+(t*ps*2/(pi*ws^2)*exp(-2*(wp(j)/2)^2/ws^2)*oem/Es))*... (t-fa/(pp(j)*2/(pi*wp(j)^2)*exp(-2*(wp(j)/2)^2/ wp(j)^2)*oabs/Ep))))*absorb/(3.14*w0^2*(1-exp(-absorb*b)))*exp(-2*((i-n)*h)^2/w0^2)*exp(-absorb*j*h); if j==1 U(i,j,T)=U(i,j,T-1)+kr*ht/ql/cl*((U(i+1,j,T-1)-2*U(i,j,T-1)+U(i-1,j,T-1))/h^2+(U(i,j,T-1)-U(i-1,j,T-1))/(i-n+1)/h^2+(-2*U(i,j,T-1)-2*h*rzxs*(U(i,j,T-1)-U0)+2*U(i,j+1,T-1)))+q(i,j)*ht/ql/cl elseif j==m U(i,j,T)=U(i,j,T-1)+kr*ht/ql/cl*((U(i+1,j,T-1)-2*U(i,j,T-1)+U(i-1,j,T-1))/h^2+(-U(i-1,j,T-1))/(i-n+1)/h^2+(U(i,j+1,T-1)-2*U(i,j,T-1)+2*h*rzxs*(U(i,j,T-1)-U0)+2*U(i,j-1,T-1)))+q(i,j)*ht/ql/cl else U(i,j,T)=U(i,j,T-1)+kr*ht/ql/cl*((U(i+1,j,T-1)-2*U(i,j,T-1)+U(i-1,j,T-1))/h^2+(U(i,j,T-1)-U(i-1,j,T-1))/(i-n+1)/h^2+(U(i,j+1,T-1)-2*U(i,j,T-1)+U(i,j-1,T-1)))+q(i,j)*ht/ql/cl end end end end % ******************************************** output=U(n+1,5,  ;output=output( ;plot(output) [ Last edited by nono2009 on 2009-12-24 at 22:02 ] |
回复此楼» 猜你喜欢
材料与化工371求调剂
已经有9人回复
-
085600调剂
已经有9人回复
-
346分的生物与医药08600求调剂
已经有6人回复
-
304求调剂（085602，过四级，一志愿985）
已经有13人回复
-
328分调剂
已经有6人回复
-
求调剂
已经有9人回复
-
358求调剂
已经有7人回复
-
生物与医药调剂
已经有3人回复
-
367求调剂
已经有3人回复
-
工科08专硕机械275求调剂
已经有3人回复
3楼2009-12-25 16:26:00
孤鸿飘逸
金虫 (正式写手)
- 应助: 1 (幼儿园)
- 金币: 5861.4
- 散金: 499
- 红花: 17
- 帖子: 617
- 在线: 432.5小时
- 虫号: 864807
- 注册: 2009-10-07
- 性别: GG
- 专业: 机器人学及机器人技术
2楼2009-12-24 23:12:35
孤鸿飘逸
金虫 (正式写手)
- 应助: 1 (幼儿园)
- 金币: 5861.4
- 散金: 499
- 红花: 17
- 帖子: 617
- 在线: 432.5小时
- 虫号: 864807
- 注册: 2009-10-07
- 性别: GG
- 专业: 机器人学及机器人技术
4楼2009-12-26 19:27:16
5楼2009-12-26 20:12:55
