科学网

 找回密码
  注册
搜索
查看: 6819|回复: 3

求教对流扩散方程的MATLAB求解问题

[复制链接]
发表于 2012-1-5 08:28:15 | 显示全部楼层 |阅读模式
是一个二维的对流扩散问题,对流在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<x<L (表示在表面(y = 0)存在药品(扩散物质)的区域(0<x<L)浓度是饱和的C0)
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<x<2 and  0<y<0.1
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)

回复

使用道具 举报

发表于 2012-2-4 03:07:23 | 显示全部楼层
回复 lygenius 的帖子

对流扩散问题应采用迎风格式,你可检查一下。
回复 支持 反对

使用道具 举报

 楼主| 发表于 2012-2-8 04:30:32 | 显示全部楼层
回复 scientifique 的帖子

这贴子以为要沉了呢,感动啊,这个问题其实现在我已经解决了,除了差分格式,还要注意步长的问题,不然很容易发散。但依然感谢这位“scientifique”仁兄提醒相助,辛苦,辛苦
回复 支持 反对

使用道具 举报

发表于 2012-10-12 21:19:01 | 显示全部楼层
这个话题有意义!学习了,楼主把正确的发上来看看:)
回复 支持 反对

使用道具 举报

Archiver|手机版|科学网 ( 京ICP备14006957 )

GMT+8, 2017-11-22 09:44

Powered by ScienceNet.cn

Copyright © 2007-2017 中国科学报社

快速回复 返回顶部 返回列表