请教关于2D-FD扩散方程的问题
这是一个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=;
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=;
b=;
c=;
g=;
A=diag(a,-1);
B=diag(b);
C=diag(c,1);
h=A+B+C;
d=h/g;
[ 本帖最后由 eight 于 2007-9-4 20:10 编辑 ]
页:
[1]