ndiffusion_um.m

用于计算脉冲激光照射下物体的温度分布

%Ndiffusion_um.m计算扩散浓度分布
%空间距离以微米为单位，以cm、m为单位时离散方程系数值十分小。
% 每次运行程序前需
%   1、导入温度分布数据
%   2、检测边界设置rd,zd，尝试求解
%   3、检查tnn,决定温度采样点与时间步进点的倍数
clear;
clc;
format long e;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%参数赋值%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load T060502.mat;
load E:\MATLAB6.5\work\lasertemprature\T05_24_1.mat
load E:\MATLAB6.5\work\lasertemprature\dd05_24_1.mat ;
% load tt06_05_2.mat;


global rd zd R Z r z rn zn td t tn L1 M1 TT;
%T=600*ones(1,10);
%dd=linspace(0,1e0,10);
Nc=1.0e20*1e-12;                    %扩散源初始浓度%单位:um^(-3)
b=1.0e0;                           %扩散源厚度%单位:um
Sc=0;
TT=abs(T);                          %温度%单位：K
tt=dd-dd(1);                        %时间％单位：s
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%空间/时间离散化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%采内点法，圆柱坐系，z方向上有zn个控制容积，r方向上有rn个控制容积，加上
%边界结点，z方向上有zn＋2个结点，r方上有rn＋2个结点。
%时间层划分，须与温度时间分布配合
%global rd zd R Z r z rn zn L1 M1 A N;

%空间划分
rd=60;                             %r方向边界%单位um
zd=1e-2;                             %z方向边界%单位um
rdd=50.;                           %热斑半径，认为近似与激光焦斑半径相等%单位：um
rn=40;                             %r方向上控制容积个数
zn=60;                             %r方向上控制容积个数
zn1=40;

%时间划分
ttt=length(T);                      %由温度的时间分布得到间隔最短的时层划分     

tnn=10;                             %ttt与实际步进的次数的倍数
tn=fix(ttt/tnn);                   %实际步进次数
td=tt(end);                        %时间终点，单位 s
t=tt([1:tnn:ttt]);                 %计算所用到的温度采样点序列

R=linspace(0,rd,rn+1);                                       %R(j)代表r方向上第j个控制容积界面，共有rn+1个
Z1=[linspace(0,zd/4,zn1+1)];
Z=[Z1([1:end-1]),linspace(zd/5,zd,zn-zn1+1)];                  %Z(j)代表z方向上第i个控制容积界面，共有zn+1个
r=[R(1),0.5*[R([1:end-1])+R([2:end])],R(end)];               %r(j)代表r方向上第j个结点的坐标，共有rn+2个
z=[Z(1),0.5*[Z([1:end-1])+Z([2:end])],Z(end)];               %z(i)代表z方向上第i个结点的坐标，共有rn+2个
L1=zn+2;
M1=rn+2;
%dz=inline('z(i)-z(i-1)','z','i');%结点距离
%Dz=inline('Z(i)-Z(i-1)','Z','i');%界面距离
%dr=inline('r(j)-r(j-1)','r','j');%结点距离
%Dr=inline('R(j)-R(j-1)','R','j');%界面距离
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%给定初始浓度分布%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N=zeros(L1*M1,tn);
Ntemp=zeros(L1*M1,1);             %用来存放上一个时层的数据

%%%%%%%%%%%%%%%%%%%%%%%%时层步进%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for l=2:tn,
%l=tn;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%方程离散化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%二维空间，有L1*M1个结点,系数矩阵为L1*M1阶方阵

A=zeros(L1*M1,7);
%对内结点离散方程计算系数
for j=2:M1-1,
    for i=2:L1-1,
        m=L1*(j-1)+i;
        dv=0.5*(r(j)+r(j-1))*(R(j)-R(j-1))*(Z(i)-Z(i-1));
        A(m,2)=r(j)*CoeffD(m,1,l*tnn)*(R(j)-R(j-1))/(z(i+1)-z(i));   %aE
        A(m,3)=r(j)*CoeffD(m,2,l*tnn)*(R(j)-R(j-1))/(z(i)-z(i-1));   %aW
        A(m,4)=r(j+1)*CoeffD(m,3,l*tnn)*(Z(i)-Z(i-1))/(r(j+1)-r(j)); %aN
        A(m,5)=r(j-1)*CoeffD(m,4,l*tnn)*(Z(i)-Z(i-1))/(r(j)-r(j-1)); %aS
        A(m,6)=dv/(tt((l)*tnn)-tt((l-1)*tnn));                       %aP0
        A(m,1)=A(m,2)+A(m,3)+A(m,4)+A(m,5)+A(m,6);                   %aP
        A(m,7)=Sc*dv+A(m,6)*Ntemp(m);                                %b常数项
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%处理边界条件%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% r=0，第二类边界条件，补充边界结点方程
 j=1;
 for i=2:L1,
     m=L1*(j-1)+i;
        A(m,1)=1;%aP
        A(m,2)=0;%aE
        A(m,3)=0;%aW
        A(m,4)=1;%aN
        A(m,5)=0;%aS
        A(m,6)=0;%aP0
        A(m,7)=0;%b常数项
 end
 
 %r=rd ，N=Na,
j=M1;
for i=2:L1,
     m=L1*(j-1)+i;
        A(m,1)=1;%aP
        A(m,2)=0;%aE
        A(m,3)=0;%aW
        A(m,4)=0;%aN
        A(m,5)=0;%aS
        A(m,6)=0;%aP0
        A(m,7)=0;%b常数项
 end
 
 %z=zd,N=Na
 i=L1;
 for j=2:M1,
     m=L1*(j-1)+i;
        A(m,1)=1;%aP
        A(m,2)=0;%aE
        A(m,3)=0;%aW
        A(m,4)=0;%aN
        A(m,5)=0;%aS
        A(m,6)=0;%aP0
        A(m,7)=0;%b常数项
 end
     
     
 %z=0,N=N(r,0,t)
 N1=0;
 for j=2:M1-1,
    for i=2:L1-1,
        m=L1*(j-1)+i;
        dv=0.5*(r(j)+r(j-1))*(R(j)-R(j-1))*(Z(j)-Z(j-1));
        N1=N1+Ntemp(m)*dv;
    end
end
N1=N1/((0.5*rdd^2)*b);
N0=Nc-N1;
i=1;
for j=1:M1,
     m=L1*(j-1)+i;
     %Ntemp(m)=N0;
        A(m,1)=1;%aP
        A(m,2)=0;%aE
        A(m,3)=0;%aW
        A(m,4)=0;%aN
        A(m,5)=0;%aS
        A(m,6)=0;%aP0
        A(m,7)=N0;%b常数项
 end
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%求解离散方程%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Ntemp=ADISolve(A,Ntemp);
 N(:,l)=Ntemp;
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
 
 %l
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%时层步进结束%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% save N06_02_3.mat N;
% save tt06_02_3.mat tt;
Nplot
hold off