请教关于2D-FD扩散方程的问题

leonidas

这是一个2D-FD扩散方程结合Crank-Nicolson法，总是得不到正确的解，解总是零，条件：表面浓度设为1，表面有掩模板，浓度分布应该是表面没掩模板的地方为1，向着x轴y 轴方向递减，边界条件按实验情况，
%输入数据
xdim=input('Enter xdim ');
ydim=input('Enter ydim ');
delx=input('Enter delx ');
dely=input('Enter dely ');
mask=input('Enter mask (um) ');

%划分网格
mask=round(mask/delx);
d=round((xdim-mask)/2);
C=zeros(ydim,xdim);
%在表面加一层作为扩散源和掩模板
C=[zeros(1,xdim);C];
C=apply_bound(C,mask,d);

%变量
Da=6e-16; %自扩散系数
M=0.72;
t_total=300; %时间
delt=10;  
q=1.6e-19;  
T=523;   %温度
maxit=xdim*ydim;

%坐标
YCoords=0:dely:dely*(ydim-1);
ycoords=YCoords.*1e6;
XCoords=0:delx:delx*(xdim-1);
XCoords=XCoords-max(XCoords)/2;
xcoords=XCoords.*1e6;

%循环
time=0;
while time<t_total,
   
% 第一步
dr=1-(1-M).*C(2:end-1,2:end-1);
diff=(C(2:end-1,3:end)-C(2:end-1,1:end-2))./(4*delx^2);

b=obtain_b(Da,delt,delx,dr);
a=obtain_a(Da,1-M,delt,delx,dr,diff);
c=obtain_c(Da,1-M,delt,delx,dr,diff);
g=obtain_g(C(:,2:end-1),Da,1-M,delt,dely,dr);

for n=2:ydim
    x=zeros(xdim,1);
    x=L(a(n-1,:),b(n-1,:),c(n-1,:),g(n-1,:));
    x=x';
    C(n,:)=x;
end

C=(10*eps).*round(C./(10*eps));
C=apply_bound(C,mask,d);

time=time+delt;

%第二步
C=C';
dr=1-(1-M).*C(2:end-1,2:end-1);
diff=(C(2:end-1,3:end)-C(2:end-1,1:end-2))./(4*dely^2);

b=obtain_b(Da,delt,dely,dr);
a=obtain_a(Da,1-M,delt,dely,dr,diff);
c=obtain_c(Da,1-M,delt,dely,dr,diff);
g=obtain_g(C(:,2:end-1),Da,1-M,delt,delx,dr);
g(:,1)=g(:,1)-C(2:end-1,1).*a(:,1);

for m=2:xdim-1,
    y=zeros(ydim+1,1);
    y=L(a(m-1,:),b(m-1,:),c(m-1,:),g(m-1,:));
    y=y';
    C(m,:)=y;
end

C=C';
C=apply_bound(C,mask,d);
C=(10*eps).*round(C./(10*eps));

time=time+delt/2;
max(max(C));
end %finish while

C=0.5.*(C+abs(C));
C=C(2:end,:);
image(C);

function C=apply_bound(C,mask,d);
C(1,1+d:mask+d)=1;
C(2,:)=C(1,:);
function a=obtain_a(D,alpha,dt,h,dr,diff);
a=((D/h)./dr).*(1/h-(alpha/2).*(diff./dr));
function b=obtain_b(D,dt,h,dr);
b=(-2/dt)-((2*D/h^2)./dr);
function c=obtain_c(D,alpha,dt,h,dr,diff);
c=((D/h)./dr).*(1/h+(alpha/2).*(diff./dr));
function g=obtain_g(C,D,alpha,dt,h,dr);
d=(C(3:end,:)-C(1:end-2,:))./(2*h);
d2=(C(3:end,:)+C(1:end-2,:)-2.*C(2:end-1,:))./(h^2);
g=-2.*C(2:end-1,:)./dt-(D./dr).*((alpha.*d.^2)./dr+d2);

function d=L(a,b,c,g)
a=[a,0];
b=[0,b,0];
c=[0,c];
g=[0,g,0];
A=diag(a,-1);
B=diag(b);
C=diag(c,1);
h=A+B+C;
d=h/g;

[ 本帖最后由 eight 于 2007-9-4 20:10 编辑 ]