资源描述
单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,单击此处编辑母版标题样式,#,环境流体力学,Matlab,编程,题目分析,由于有害物质在池底均匀分布,池底及池壁对有害物质不完全吸收,则可按,瞬时平面源一维扩散,处理。,已知条件:,池底面积,200m200m,水深,50m,有害物质质量,4000kg,扩散系数,1.0cm,2,/s,1.,解析法,1.1,计算公式,因底部不完全吸收,由底部沿垂向的扩散浓度应为,式中,x,为自池底量起的距离,,kg/m,扩散系数,D=1.0cm,2,/s=8.64m,2,/,天。,2.,有限差分法,2.1,网格剖分,令,x,=,xi,=,ih,,,h,=1m,,,0,i,50,t,=,tk,=,k,,,=1,天,,0,k,365,一维扩散的扩散方程为:,2.2,古典隐格式,稳定性,和,收敛性,差分格式为(与解析法区分,差分法用,u,表示,C,):,令步长比,表示成矩阵形式,A,B,3.,Matlab,程序,clc;clear;close all;%,清除屏幕上的内容,清除工作区里的内容,关闭画的图,x=-8.64*ones(1,48);,A=18.28*eye(49)+diag(x,-1)+diag(x,1);,InvA=inv(A);%,求矩阵,A,的逆矩阵,M=0.1;%,参数赋值,D=8.64;,r=8.64;,l=50;,B=zeros(49,365);,for t=1:365,B(1,t)=r*2*M/sqrt(4*pi*D*t);,B(49,t)=r*2*M/sqrt(4*pi*D*t)*exp(-l2/4/D/t);,end,U=zeros(49,365);%,用来记录,U,的数值,横坐标为,1365,天,纵坐标为,149m,。初始化全为,0,。,for t=1:365,if t=1,U(:,t)=InvA*B(:,t);,else U(:,t)=InvA*(U(:,t-1)+B(:,t);,end,end,figure%,水池污染物浓度作图,hold on;,grid;,plot(U(5,:),k),plot(U(10,:),b),plot(U(15,:),y),plot(U(20,:),r),plot(U(25,:),g),plot(U(30,:),m),plot(U(35,:),c),plot(U(40,:),r:),plot(U(45,:),b:),axis(0 370 0 0.014);%,横坐标取,0-370,,纵坐标取,0-0.014,title(,水池污染物浓度,);%,加图形标题,xlabel(,时间,/d);%,加,x,轴说明,ylabel(,浓度,/kg.m-3);%,加,y,轴说明,legend(i=5,i=10,i=15,i=20,i=25,i=30,i=35,i=40,i=45);%,加图例,C=zeros(50,365);%,解析法计算各点的值,for k=1:50,for t=1:365,C(k,t)=2*M/sqrt(4*pi*D*t)*exp(-k2/4/D/t);,end,end,figure%,不同点处解析法计算与有限差分法计算比较,subplot(3,3,1),hold on;,grid;,plot(U(5,:),r);,plot(C(5,:),b);,title(i=5);,xlabel(,时间,/d);,ylabel(,浓度,/kg.m-3);,图,1,有限差分法污染物浓度随时间变化的曲线,图,2,有限差分法和解析法的计算结果比较,谢谢!,
展开阅读全文