是一个二维的对流扩散问题，对流在X方向， 扩散作用在Y 方向，所以方程为 dC/dt=D*d2C/dy2-V*dC/dx (注：
D= Diffusion coefficient, C=concentration, V=
流体速度，并假设它与表面的垂直方向Y成正比，即：V=a*Y. d2C/dy2表示浓度对Y的二阶导数。 ) 边界条件为：
1 ,C = C0 when y = 0 and 0
(表示在表面(y =
0)存在药品(扩散物质)的区域(0
2, C = 0 when y> 0.1cm for all x
(表示对难容物质来说，由于流动的因素，离表面一定距离后不存在溶质)
3, C = 0 when x = 0 for all y (因为在存在药品区域之前是 没有浓度的)
4, C = 0 when y = 0 and x> L (表示
在不存在药品的表面也没有浓度)
本人用论文中介绍的显示差分算法在MATLAB 编程 [C(i , j , k+1)- C(i , j , k)]/dt =
D/(dy)^2 * [ C(i , j-1 , k) - 2*C(i , j , k) + C(i , j+1 , k)] -
V/dx * [C(i , j , k) - C(i -1, j , k)]
下面是我的代码：
问题是，虽然可以得到仿真(步长设置不好还容易发散)，但是明显没有对流项的作用，画出的图应该是从左到右被“冲”的感觉，实在找不到问题出在哪儿了，恳请大侠赐教，不胜感谢！
%%%%%%%%%%%%%%%%%%
clear
clc
D = [0 ,2 , 0 , 0.1 ] %%%% region of
0
0
Mx =200
My =100
N = 1000
T = 10
dif = 9.86*10^(-6) %%%%diffusion coefficient
alfa = 4.25
visco = 0.89
dx = (D(2) - D(1))/Mx;
dy = (D(4) - D(3))/My;
dt = T/N; t = [0:N]*dt;
%initialization and boundary condition
for j = 1:Mx + 1
for i = 1:My + 1
u(i,j) = 0 ;
end
end
for j = 1:Mx /2
u(1,j) = 1;
u(My+1,j) = 0;
end
ry = dif*dt/(dy*dy);
for k = 1:N
u_1 = u; t = k*dt;
for i= 2:My+1
vx = alfa*i*D(4)/My; %%%%to be
rx = vx*dt/dx;
rx1 = 1 - rx;
for j =2 : Mx
u(i , j )=ry*(u_1(i,j-1)-2*u_1(i , j)+u_1(i , j+1)) +rx1*u_1(i ,
j)+rx*u_1(i-1 , j);
end
end
for i = 1: My + 1
u(i , Mx + 1)=u(i , Mx);
end
end
i = 1:Mx+1;
j= 1: My+1;
[I, J]=meshgrid(i,j);
surf(I , J , u)
figure;
mesh(I, J , u)
