关于降雨入渗模型的计算（实际就是二维扩散方程）
本人正在做本科毕业论文，需要编制一个计算二维降雨入渗模型的计算程序，数学模型已经建立完毕，可是编程序没有头绪，哪位大侠有这方面的资料啊？最好能提供相似的源程序可以参考，matlab或者c++的都可以，拜托啦！ %环境设置，转换后需要去除clear;
clc;
%开始各种变量设置
%------------------------------------------------------------------
%基坑宽度和深度Lx，Ly
Lx=300;
Lz=300;
%计算内的时间段T
T=200;
%定义时间深度和宽度差值dt dz dx
dt=1;
dz=5;
dx=5;
%初始含水率v0
v0=0.28;
%饱和含水率vs
vs=0.45;
%--------------------------------------------------------------------
%定义差分网格节点数
%竖直节点数Nz
Nz=Lz/dz+1;
%水平节点数Nx
Nx=Lx/dx+1;
%时间节点数M
M=T/dt+1;
%方程运算初始化和边界条件设置
%---------------------------------------------------------------------
%含水状态矩阵初始化
%v(k,i,j)表示第k时刻坐标为(i,j)的网格点的含水率
%所有网格点的初始含水率值为v0
for k=1:M
for i=1:Nx
for j=1:Nz
v(k,i,j)=v0;
end
end
end
%定义系数矩阵的初始值:
%方程一：a1(i,j),b1(i,j),c1(i,j),h1(i,j)
%方程二：a2(i,j),b2(i,j),c2(i,j),h2(i,j)
a1=zeros(Nz,Nz);
b1=zeros(Nz,Nz);
c1=zeros(Nz,Nz);
h1=zeros(Nz,Nz);
A1=zeros(Nz,Nz);
H1=zeros(Nz,1);
a2=zeros(Nx,Nx);
b2=zeros(Nx,Nx);
c2=zeros(Nx,Nx);
h2=zeros(Nx,Nx);
A2=zeros(Nx,Nx);
H2=zeros(Nx,1);
%输出文件定义
fileout=fopen('fileout.dat','wt');
%初始条件,t=0
%t=0时所有点的含水率为v0
for i=1:Nx
for j=1:Nz
v(1,i,j)=v0;
end
end
%上边界条件,z=0
%z=0,t>ta,即地表的含水率为vs
%？？？？？？？？这个设置与论文不符 ？？？？？？？？？？
%下边界条件，z=Nz
%所有深度为Nz的网格点的含水率为v0
for k=1:M
for i=1:Nx
v(k,i,Nz)=v0;
end
end
%---------------------------------------------------------------------
%开始差分方程组计算
r1=dt/(dx*dx);
r2=dt/(dz*dz);
r3=dt/(dz*2);
%计算v(k,i,j)
for k=1:M
if mod(k,2)==1;
for j=2:Nz-1
for i=2:Nx-1
%D代表z方向扩散系数，K代表x方向扩散系数
%交替使用两种差分格式以简化运算和消除误差
%时段除以2看是否为1 判断是否为奇数
D1 = 5.2*((v(k+1,i,j)/vs)^3.6+(v(k+1,i-1,j)/vs)^3.6);
D2 = 5.2*((v(k+1,i+1,j)/vs)^3.6+(v(k+1,i,j)/vs)^3.6);
D3 = 5.2*((v(k,i,j)/vs)^3.6+(v(k,i,j-1)/vs)^3.6);
D4 = 5.2*((v(k,i,j+1)/vs)^3.6+(v(k,i,j)/vs)^3.6);
K = 0.02*((v(k,i,j+1)/vs)^2.8-(v(k,i,j-1)/vs)^2.8);
a1(i,j) = r1*D1;
b1(i,j) = -(1+r1*(D1+D2));
c1(i,j) = r1*D2;
h1(i,j) = -r2*D3*v(k,i,j-1)+(r2*(D4+D3)-1)*v(k,i,j)-r2*D4*v(k,i,j+1)+r3*K;
end
%处理特殊的系数值
b1(2,j) = r1*5.2*((v(k+1,2,2)/vs)^3.6+(v(k+1,1,2)/vs)^3.6)-(1+r1*( 5.2*((v(k+1,2,2)/vs)^3.6+(v(k+1,1,2)/vs)^3.6)...
+ 5.2*((v(k+1,2,2)/vs)^3.6+(v(k+1,3,2)/vs)^3.6)));
b1(Nx,j) =r1*5.2*((v(k+1,Nx,2)/vs)^3.6+(v(k+1,Nx-1,2)/vs)^3.6)-(1+r1*( 5.2*((v(k+1,Nx,2)/vs)^3.6+(v(k+1,Nx-1,2)/vs)^3.6)...
+ 5.2*((v(k+1,Nx-1,2)/vs)^3.6+(v(k+1,Nx-2,2)/vs)^3.6)));
%针对固定的k,j进行方程组运算
%定义并求解方程组A1x=H
A1= diag(a1(3:Nx,j),-1)+diag(b1(2:Nx,j))+diag(c1(2:Nx-1,j),1);
H1=h1(2:Nx,j);
v(k+1,2:Nx,j)=crout(A1,H1);
%输出到数据文件
fprintf(fileout,'\n%\n','k='+k);
for i=1:Nx
for j=1:Nz
fprintf(fileout,'%8.7f',v(k,i,j));
end
fprintf(fileout,'\n');
end
fprintf(fileout,'\n%\n','end of k');
end
else
for i=2:Nx-1
for j=2:Nz-1
D1 = 5.2*((v(k,i,j)/vs)^3.6+(v(k,i-1,j)/vs)^3.6);
D2 = 5.2*((v(k,i+1,j)/vs)^3.6+(v(k,i,j)/vs)^3.6);
D3 = 5.2*((v(k+1,i,j)/vs)^3.6+(v(k+1,i,j-1)/vs)^3.6);
D4 = 5.2*((v(k+1,i,j+1)/vs)^3.6+(v(k+1,i,j)/vs)^3.6);
K = 0.02*((v(k+1,i,j+1)/vs)^2.8-(v(k+1,i,j-1)/vs)^2.8);
a2(i,j) = r2*D3;
b2(i,j) = -(1+r2*(D3+D4));
c2(i,j) = r2*D4;
h2(i,j) = -r1*D1*v(k,i-1,j)+(r1*(D1+D2)-1)*v(k,i,j)-r1*D2*v(k,i+1,j)+r3*K;
end
%处理特殊的系数值
h2(i,2) = r1*5.2*((v(k,2,2)/vs)^3.6+(v(k,1,2)/vs)^3.6)*v(k,1,2)+(r1*( 5.2*((v(k,2,2)/vs)^3.6+(v(k,1,2)/vs)^3.6+...
((v(k,2,2)/vs)^3.6+(v(k,3,2)/vs)^3.6))-1)*v(k,2,2) - r1*5.2*((v(k,2,2)/vs)^3.6+(v(k,3,2)/vs)^3.6)*v(k,3,2)+...
r3*0.02*((v(k+1,2,3)/vs)^2.8-(v(k+1,2,1)/vs)^2.8)-r2*5.2*((v(k+1,2,2)/vs)^3.6+(v(k+1,2,1)/vs)^3.6))*v(k+1,2,1);
h2(i,Nz)=r1*5.2*((v(k,2,Nz-1)/vs)^3.6+(v(k,1,Nz-1)/vs)^3.6)*v(k,1,Nz-1)+(r1*( 5.2*((v(k,2,Nz-1)/vs)^3.6+(v(k,1,Nz-1)/vs)^3.6+...
((v(k,2,Nz-1)/vs)^3.6+(v(k,3,Nz-1)/vs)^3.6))-1)*v(k,2,Nz-1) - r1*5.2*((v(k,2,Nz-1)/vs)^3.6+(v(k,3,Nz-1)/vs)^3.6)*v(k,3,Nz-1)+...
r3*0.02*((v(k+1,2,Nz)/vs)^2.8-(v(k+1,2,Nz-2)/vs)^2.8)-r2*5.2*((v(k+1,2,Nz-1)/vs)^3.6+(v(k+1,2,Nz-2)/vs)^3.6))*v(k+1,2,Nz);
%针对固定的k,j进行方程组运算
%定义并求解方程组A1x=H
A2= diag(a2(i,3:Nz),-1)+diag(b2(i,2:Nz))+diag(c2(i,2:Nz-1),1);
H2=h2(i,2:Nz)';
v(k+1,i,2:Nz)=crout(A2,H2)';
%输出到数据文件
fprintf(fileout,'\n%\n','k='+k);
for i=1:Nx
for j=1:Nz
fprintf(fileout,'%8.7f',v(k,i,j));
end
fprintf(fileout,'\n');
end
fprintf(fileout,'\n%\n','end of k');
end
end
for j=1:Nz-1
if v(k+1,1,j)>=vs
v(k+1,1,j)=vs;
end
end
end
%计算结束
%---------------------------------------------------------------------
%开始输出结果图像
z=0:dz:Lz;
x=0:dx:Lx;
%time=input('请输入想要查看的时刻k！\n','d')
%if(time >M || time <0)
%output('输入错误！')
%end
=mesh(x,z);
surf(X,Z,v(time,:,:))
grid on
%保存文件
fclose(fileout);
%清理运行环境
clear;
clc
怎么没人回复呢，这是编的一个程序，麻烦高手给看看啦~ 都本科毕业了，还这水平？连matlab都不会，临阵抱佛脚，不好。
另外，你的方程，都没提供。
你给出的代码，也忘了去掉勾选表情符号。
你先想到的是有没有参考的源代码，而不愿意花半天，去了解一下matlab 或者 C 怎么编程。 我专业又不是数学或者应用信计的，俺是土木工程滴，把一个土木工程的实际问题搞到这种程度还可以啦，不要打击人呢 我专业又不是数学或者应用信计的，俺是土木工程滴，把一个土木工程的实际问题搞到这种程度还可以啦，不要打击人呢 楼主好会找理由，学土木的怎么啦，我也是学土木的，matlab已经成为一个通用软件，理工科的大学生都应该掌握，不需要理由！ 6# liuyanchang
我没有说过我不会啊~
那您能给我看看上面的那个程序吗？另外给你传一个完善后的程序，麻烦您把求出来的数据绘出图来好吗？二维降雨入渗模型，结果是不同时刻的体积含水率，v(k,i,j)，固定k值，画一个，固定i值画一个图，固定j值画一个图，谢谢啦，小弟实在愚笨，本科的时候最怕的就是编程
页:
[1]