- 积分
- 0
- 注册时间
- 2009-4-18
- 仿真币
-
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2009-4-21 16:10:03
|
显示全部楼层
来自 北京海淀
%环境设置,转换后需要去除
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
[X,Z]=mesh(x,z);
surf(X,Z,v(time,:,:))
grid on
%保存文件
fclose(fileout);
%清理运行环境
clear;
clc
怎么没人回复呢,这是编的一个程序,麻烦高手给看看啦~ |
|