Free surface flows 103
else u(i)=u(i)+Dt*g*(S-(h(i)-h(i-1))/
Dx-u(i)*abs(u(i))/C^2*2/(h(i)+h(i-1)));
end
Q(i)=u(i)*b(i)*(h(i)+h(i-1))/2;
end
% Downstream boundary condition;
Q(nx)=b(nx)*h(nx-1)*sqrt(g*h(nx-1));
% Estimation of the water depths:
for i=1:nx-1
h(i)=h(i)-Dt/Dx/(b(i)+b(i+1))*2*(Q(i+1)-Q(i));
end
% Data stored at different time intervals for plotting;
if k == round(nt/25)
for m=1:nx
h1t(m)=h(m);
Q1t(m)=Q(m);
end
end
if k == round(nt/5)
for m=1:nx
h3t(m)=h(m);
Q3t(m)=Q(m);
end
end
if k == round(nt)
for m=1:nx
h5t(m)=h(m);
Q5t(m)=Q(m);
end
end
end
for j=1:ntot
Qinflow(j)=Qin-Qin*j/ntot;
end
plot(1:ntot,Qinflow,'k','Linewidth',1.5)
xlabel('Time*Dt [s]')
ylabel('Inflow hydrograph [m^3/s]')
axis([0,51000,0,2100])
text(10000,200, 'Td=Flood duration: 5*10^4*Dt [s]')
figure, plot(1:nx,Q1t,'b','Linewidth',1.5)
hold on
plot(1:nx,Q3t,'r','Linewidth',1.5)
plot(1:nx,Q5t,'g','Linewidth',1.5)
xlabel('Longitudinal distance x 1000 [m]')
ylabel('Water discharge [m^3/s]')
axis([0,47,0,2000])
legend('t=Td/25','t=Td/5','t=Td')
figure, plot(1:nx,h1t,'b','Linewidth',1.5)
hold on