Free surface flows 89
dyc(1)=0.1;
while (abs(dyc(k))>0.0001)
Ac(k)=b(j)*yc(k)+m(j)*(yc(k))^2;
Tc(k)=b(j)+2*m(j)*yc(k);
fc(k)=Ac(k)^(3/2)*Tc(k)^(-1/2)-Q/sqrt(g);
dfc(k)=-m(j)*(Ac(k)/Tc(k))^(3/2)+(3/2)*sqrt(Ac(k)*T
c(k));
%
dfc(k)=-m(j)*Ac(k)^(3/2)*Tc(k)^(-3/2)+
(3/2)*Tc(k)^(-1/2)*Ac(k)^(1/2)*Tc(k);
yc(k+1)=yc(k)-fc(k)/dfc(k);
dyc(k+1)=-fc(k)/dfc(k);
k=k+1;
end
ycr(j)=yc(k);
iter_yc(j)=k;
end
% Water surface profile calculations for the trapezoidal
channel;
Dx=1;
L=4000;
nx=fix(abs(L/Dx));
% Subcritical flow boundary condition;
% Water depth at the control section yc Y =1.5*yn(1);
% Water depth at the control section yc % Y =(yn(1)+ycr(1))/2;
% Supercritical flow boundary condition;
% Water depth at the control section yn % Y=(yn(1)+ycr(1))/2;
% Water depth at the control section y % Y=0.1*yn(1);
if yn-ycr>0
Dx=-Dx;
end
z(1)=0;
h(1)=Y(1);
b=b(1);
m=m(1);
% Main program;
for i=1:nx;
A(i)=b*Y(i)+m*Y(i)^2;
T(i)=b+2*m*Y(i);
Pw(i)=b+2*Y(i)*sqrt(1+m^2);
Rh(i)=A(i)/Pw(i);
u(i)=Q/A(i);
Sf(i)=n^2*u(i)^2/(Rh(i)^(4/3));
% ODE f=dy/dx=(S-Sf)/(1-Fr^2);
f(i)=(S-Sf(i))/(1-Q^2*T(i)/g/A(i)^3);
yy(i)=Y(i)+0.5*f(i)*Dx;
A(i)=b*yy(i)+m*yy(i)^2;