编程实现:
NX=100;
NY=40;
NT=73;
ox=50/100;
ot=1;
P0=75000;
Pa=101325;
A1=(1-0.873/800)*220*0.15*0.1386;
A2=(1-0.873/800)*220*0.55*0.0231;
A3=(1-0.873/800)*220*0.3*0.017328;
%初始条件
for j=1:NX+1
for l=1:NY+1
P(j,l,1)=P0;
P(j,l,2)=P0;
end
end
%顶部边界条件
for m=2:NT
for j=2:NX-1
P(j,1,m)=Pa;
end
end
%左边边界条件
for m=2:NT
for l=1:NY
P(1,l,m)=-2000;
end
end
%底部边界条件
for m=2:NT
for j=2:NX-1
P(j,NY+1,m)=P(j,NY-1,m);
end
end
%右边边界条件
for m=2:NT
for l=2:NY-1
P(NX+1,l,m)=P(NX-1,l,m);
end
end
%内部网格
for m=2:NT
for i=1:10
for j=2:NX-1
for l=2:NY-1
t(l,m-1)=(m-1)*ot+l*ox/2;
B(l,m-1)=0.873*1100*9.8*l*ox*0.365*10^(-6)*exp(-10^(-3)*(m-1)*ot*365);
f(l,m-1)=1.145*(A1*exp(-0.1386*t(l,m-1))+A2*exp(-0.0231*t(l,m-1))+A3*exp(-0.017328*t(l,m-1))+B(l,m-1));
P(j,l,m)=4.3*365*24*3600*((3*P(j+1,l,m)+3*P(j-1,l,m)+P(j,l+1,m)+P(j,l-1,m)-8*P(j,l,m))/(ox*ox))*ot+f(l,m-1)*ot+P(j,l,m-1)
end
end
for j=2:NX-1
t(NY,m-1)=(m-1)*ot+NY*ox/2;
B(NY,m-1)=0.873*1100*9.8*NY*ox*0.365*10^(-6)*exp(-10^(-3)*(m-1)*ot*365);
f(NY,m-1)=1.145*(A1*exp(-0.1386*t(NY,m-1))+A2*exp(-0.0231*t(NY,m-1))+A3*exp(-0.017328*t(NY,m-1))+B(NY,m-1));
P(j,NY,m)=4.3*365*24*3600*((3*P(j+1,NY,m)+3*P(j-1,NY,m)+P(j,NY+1,m)+P(j,NY-1,m)-8*P(j,l,m))/(ox*ox))*ot+f(NY,m-1)*ot+P(j,NY,m-1);
end
for l=2:NY-1
t(l,m-1)=(m-1)*ot+l*ox/2;
B(l,m-1)=0.873*1100*9.8*l*ox*0.365*10^(-6)*exp(-10^(-3)*(m-1)*ot*365);
f(l,m-1)=1.145*(A1*exp(-0.1386*t(l,m-1))+A2*exp(-0.0231*t(l,m-1))+A3*exp(-0.017328*t(l,m-1))+B(NY,m-1));
P(NX,l,m)=4.3*365*24*3600*((3*P(NX+1,l,m)+3*P(NX-1,l,m)+P(NX,l+1,m)+P(NX,l-1,m)-8*P(j,l,m))/(ox*ox))*ot+f(l,m-1)*ot+P(j,l,m-1);
end
end
end