间断污染点源的四维计算机模拟

下面来介绍一下利用MATLAB对某污染点源进行四维模拟，具体如下：

1、首先在MATLAB的主界面编辑器中分别写入下列两个代码：

主程序代码：

%这是主程序,命名为Main1206.m
clc;clear all;close all;
%本程序研究的对象是瞬时点源污染扩散问题
subplot(1,2,1);
daspect([1 1 1]);
xmin=-500;dx=10;xmax=500;ymin=-100;dy=5;ymax=100;zmin=-10;dz=1;zmax=10;
[cx,cy,ct]=meshgrid(xmin:dx:xmax,ymin:dy:ymax,10:10:50);
[grad_cx,grad_cy,grad_cz,grad_ct]=grad(cx,cy,0,ct);
hcone=coneplot(-500:10:500,-100:5:100,10:10:50,grad_cx,grad_cy,grad_ct,cx,cy,ct,0.9);
xlim([-60 70]);ylim([-50 25]);zlim([0 150]);
grid on;camlight right;lighting phong;camproj perspective;
xlabel('横向/(mg/L*s)','fontsize',18);ylabel('纵向/(mg/L*s)','fontsize',18);
zlabel('时间/s','fontsize',18);view(3);
title('10秒至50秒XOY切面扩散梯度流锥图','fontsize',18);
subplot(1,2,2);
[x,y,z]=meshgrid(xmin:dx:xmax,ymin:dy:ymax,zmin:dz:zmax);
[grad_x,grad_y,grad_z,grad_t]=grad(x,y,z,20);
quiver3(x,y,z,grad_x,grad_y,grad_z,1.1,'r-');
xlabel('横向/(mg/L*s)','fontsize',18);ylabel('纵向/(mg/L*s)','fontsize',18);
zlabel('河流竖直方向(mg/L*s)','fontsize',18);
axis tight;
title('时间T=20秒时3D扩散梯度向量图','fontsize',18);

子程序代码：

% 这是子程序，命名为grad.m
function [grad_x,grad_y,grad_z,grad_t]=grad(x,y,z,t)
%以下各梯度grad_x,grad_y,grad_z,grad_t的表达式可以由此命令用计算机求出：
%syms t M Dx Dy Dz x ux y uy z uz K
%grad_maple=maple('grad(M/(8*(pi*t)^(3*2)*sqrt(Dx*Dy*Dz))*exp(-(x-ux*t)^2/(
%4*Dx*t)-(y-uy*t)^2/(4*Dy*t)-(z-uz*t)^2/(4*Dz*t))*exp(-K*t),vector([x,y,z,t]));')
%先用simplify和collect命令对各变量梯度表达式进行化简和合并
M=200400;K=4.2/(24*60*6);
Dx=50;Dy=5;Dz=2;ux=1.5;uy=0.2;uz=0.1;
grad_x=-1./16.*M./pi.^6./t.^7./(Dx.*Dy.*Dz).^(1./2).*(x-ux.*t)./Dx.*exp(-1./4.*(x-ux.*t).^2./Dx./t-1./4.*(y-uy.*t).^2./Dy./t-1./4.*(z-uz.*t).^2./Dz./t).*exp(-K.*t);
grad_y=-1./16.*M./pi.^6./t.^7./(Dx.*Dy.*Dz).^(1./2).*(y-uy.*t)./Dy.*exp(-1./4.*(x-ux.*t).^2./Dx./t-1./4.*(y-uy.*t).^2./Dy./t-1./4.*(z-uz.*t).^2./Dz./t).*exp(-K.*t);
grad_z=-1./16.*M./pi.^6./t.^7./(Dx.*Dy.*Dz).^(1./2).*(z-uz.*t)./Dz.*exp(-1./4.*(x-ux.*t).^2./Dx./t-1./4.*(y-uy.*t).^2./Dy./t-1./4.*(z-uz.*t).^2./Dz./t).*exp(-K.*t);
grad_t=-3./4.*M./pi.^6./t.^7./(Dx.*Dy.*Dz).^(1./2).*exp(-1./4.*(x-ux.*t).^2./Dx./t-1./4.*(y-uy.*t).^2./Dy./t-1./4.*(z-uz.*t).^2./Dz./t).*exp(-K.*t)+1./8.*M./pi.^6./t.^6./(Dx.*Dy.*Dz).^(1./2).*(1./2.*(x-ux.*t)./Dx./t.*ux+1./4.*(x-ux.*t).^2./Dx./t.^2+1./2.*(y-uy.*t)./Dy./t.*uy+1./4.*(y-uy.*t).^2./Dy./t.^2+1./2.*(z-uz.*t)./Dz./t.*uz+1./4.*(z-uz.*t).^2./Dz./t.^2).*exp(-1./4.*(x-ux.*t).^2./Dx./t-1./4.*(y-uy.*t).^2./Dy./t-1./4.*(z-uz.*t).^2./Dz./t).*exp(-K.*t)-1./8.*M./pi.^6./t.^6./(Dx.*Dy.*Dz).^(1./2).*exp(-1./4.*(x-ux.*t).^2./Dx./t-1./4.*(y-uy.*t).^2./Dy./t-1./4.*(z-uz.*t).^2./Dz./t).*K.*exp(-K.*t);

2、命名保存好至自定义路径下，点击运行。结果如下：

从它的切面扩散梯度流锥图可以看出，随着时间的推移，它的扩散梯度的程度就越发猛烈，但是扩散速度却越来越趋于平稳，而且在达到一定时间后，基本保持稳定状态了，纵方向的相比于横向的要更为明显一些。其次，在3D扩散梯度向量图中可以看出，河流竖直方向的流量与其扩散梯度近似成反比关系，在一定的时间内，污染物的扩散梯度形态是：在扩散原点附近，污染物扩散不活跃，这是因为于污染源附近的污染物浓度过大，污染物浓度已经达到饱和或过饱和状态，各个污染物分子之间相互抑制，不容易发生溶解和水解，所以沿着x、y、z三个方向的扩散速度是先增大后减少，请大家继续关注！！！