OPEN CHANNEL FLOW profiles [Matlab]

Prev: How to use callAllOptimPlotFcns
Next: rayleighchan and mlseeq help....???

From: PiYaL on 15 Nov 2009 00:19

I need a matlab scripts developed for >>>improved euler method,
standard step method and fourth order Runge kutta method and
simultaneous solution method.
USING EACH METHOD >M1, M2, and S3 profiles has to be plotted . Plots
should contain bed elevation , normal depth, critical depth lines.in
case of M1 profile , set the downstream depth as 2*Yn and in case of
M2 profile the downstream depth 1.05*Yc and in case of S3 profile
Upstream depth as 0.05m.

The below sample script done for modified euler method>>
Thanks in advance.

% Calculating the Normal and Critical Depth for Different Open Channel
% Cross Sections.

% Parameters for Trapezoidal, Rectangular, and Triangular Cross
Section
clear all;
Q=100; %Discharge (m^3/s)
s=0.00025; %Bottom slope
n=0.040; %Manning Roughness Coefficient
yinitial=0.1; %Initial Guess to start the iterations
%% Normal and Critical Depths Calculation

%Loop for Normal Depth

for j=1:3
b=[10 50 0]; %Bottom Width (m) for Trap, Rect, Tria. x-
sec. respectively.
z=[1.0 0.0 1.0]; %Side Slope for Trapezoidal,
Rectangular, and Triangular X-sec respectively.
y(1)=yinitial;
i=1;
dy(1)=1e-1;
while (abs(dy(i))>1e-3)
A(i)=b(j)*y(i)+z(j)*(y(i))^2;
P(i)=b(j)+2*(z(j)^2+1)^(0.50)*y(i);
R(i)=A(i)/(P(i));
f(i)=sqrt(s)/n*(A(i)*R(i)^(2/3))-Q; % function
ff(i)=sqrt(s)/n*((-4/3)*R(i)^(5/3)*sqrt(1+z(j)^2)+5/3*R(i)^
(2/3)*(b(j)+2*z(j)*y(i))); % f' Derivative
y(i+1)=y(i)-f(i)/ff(i);
dy(i+1)=-f(i)/ff(i);
i=i+1;
end
normaldepth(j)=y(i);
conv(j)=i;

end

for j=1:3 %Side Slope for Trapezoidal, Rectangular, and
Triangular X-sec respectively.
%Loop for Critical Depth
alpha=1;
theta=0;
yc(1)=yinitial;
ic=1;
dyc(1)=1e-1;
while (abs(dyc(ic))>1e-3)

Ac(ic)=b(j)*yc(ic)+z(j)*(yc(ic))^2;
Bc(ic)=b(j)+2*z(j)*yc(ic);
Pc(ic)=b(j)+2*(z(j)^2+1)^(0.50)*yc(ic);
Rc(ic)=Ac(ic)/(Pc(ic));
Dc(ic)=Ac(ic)/(b(j)+2*z(j)*yc(ic));
fc(ic)=Ac(ic)^(3/2)*Bc(ic)^(-1/2)-Q/sqrt(9.81/alpha); %
function
ffc(ic)=Ac(ic)^(3/2)*(-1/2)*Bc(ic)^(-3/2)*2*z(j)+Bc(ic)^(-1/2)*
(3/2)*Ac(ic)^(1/2)*Bc(ic); % f' Derivative
yc(ic+1)=yc(ic)-fc(ic)/ffc(ic);
dyc(ic+1)=-fc(ic)/ffc(ic);
ic=ic+1;
end
criticaldepth(j)=yc(ic);
convc(j)=ic;
end
%% Water Surface Calculations

dx=2.0; %Interval in channel direction
Y =2.0*normaldepth(1); %Initial Water Depth at the control location
%Y =1.05*criticaldepth(1);
%Y=0.05;
Dir ='U'; %direction of solution : Upstream(U) or
Downstream(D)?', 's');
if Dir == 'U';
dx=-dx;
end
L=40000; %total length of the channel
Nx=round(abs(L/dx));
BedElev(1)=0; %Begining of Bed Elevation
WaterElev(1)=Y(1); %Begining of Water Elevation
b=b(1); % channel bottom width for Trap.
Section
z=z(1); % Side Slope for Trap. Section
for i=1:Nx;
A(i)=b*Y(i)+z*Y(i)^2; % Area m^2
B(i)=b+2*z*Y(i); % Top Width
P(i)=b+2*Y(i)*sqrt(1+z^2); % Wetted Perimeter
R(i)=A(i)/P(i); % Hydr. radius
V(i)=Q/A(i); % Velocity
Sf(i)=n^2*V(i)^2/(R(i)^(4/3)); % Energy slope
f(i)=(s-Sf(i))/(1-Q^2*B(i)/9.81/A(i)^3); % dy/dx=f=(s-sf/(1-fr^2))
yy(i)=Y(i)+0.5*f(i)*dx; % first modified depth or
y(i+1/2)=yold+.5*dx*f'
A(i)=b*yy(i)+z*yy(i)^2; % Calculate new area but
using yy which is ymod1
B(i)=b+2*z*yy(i); % Calculate new top width
but using yy which is ymod1
P(i)=b+2*yy(i)*sqrt(1+z^2); % Calculate new wetted
perimeter but using yy which is ymod1
R(i)=A(i)/P(i);
V(i)=Q/A(i);
Sf(i)=n^2*V(i)^2/(R(i)^(4/3));
ff(i)=(s-Sf(i))/(1-Q^2*B(i)/9.81/A(i)^3); %Calculate new dy/dx but
using yy which is ymod1
Y(i+1)=Y(i)+ff(i)*dx; % Calcualte final y
which is yold+f'(y(i+1/2))*dx
BedElev(i+1)=-dx*(i)*s; % Bed elevation
for plotting purpose
WaterElev(i+1)=BedElev(i+1)+Y(i+1); % Water
surface elevation for plotting purpose
X(i+1)=dx*i; % X: distance along
the channel in (m)
end
CriticalElev=BedElev.*1.0+criticaldepth(1); % get the critical depth
line for plotting purpose
NormalElev=BedElev.*1.0+normaldepth(1); % get the normal depth
line for plotting purpose
Yfinalm=Y(i+1);
figure(1);
plot(X,BedElev,'k-','LineWidth',2); hold;
plot (X,WaterElev,'LineWidth',2);
plot(X,CriticalElev,'-.r')
plot(X,NormalElev,'--m' )
h = legend('Bed','Water Surface','Critical Depth Line','Normal Depth
Line',1);
xlabel('Length(m)');
ylabel('Elevation(m)');
title('Water Surface Profile');
XY=[X' BedElev' WaterElev'];
d = {'X', 'Bed Elev','WaterElev'};
dir=pwd;
xlswrite('water depth.xlsx',d, 'WaterSurface' ,'A1')
xlswrite('water depth.xlsx',XY,'WaterSurface','A2');

From: TideMan on 15 Nov 2009 05:14

On Nov 15, 6:19 pm, PiYaL <yasha...(a)yahoo.com> wrote:
> I need a matlab scripts developed for >>>improved euler method,
> standard step method and fourth order Runge kutta method and
> simultaneous solution method.
> USING EACH METHOD >M1, M2, and S3 profiles has to be plotted . Plots
> should contain bed elevation , normal depth, critical depth lines.in
> case of M1 profile , set the downstream depth as 2*Yn and in case of
> M2 profile the downstream depth 1.05*Yc and in case of S3 profile
> Upstream depth as 0.05m.
>
> The below sample script done for modified euler method>>
> Thanks in advance.
>
> % Calculating the Normal and Critical Depth for Different Open Channel
> % Cross Sections.
>
> % Parameters for Trapezoidal, Rectangular, and Triangular Cross
> Section
> clear all;
> Q=100; %Discharge (m^3/s)
> s=0.00025; %Bottom slope
> n=0.040; %Manning Roughness Coefficient
> yinitial=0.1; %Initial Guess to start the iterations
> %% Normal and Critical Depths Calculation
>
> %Loop for Normal Depth
>
> for j=1:3
> b=[10 50 0]; %Bottom Width (m) for Trap, Rect, Tria. x-
> sec. respectively.
> z=[1.0 0.0 1.0]; %Side Slope for Trapezoidal,
> Rectangular, and Triangular X-sec respectively.
> y(1)=yinitial;
> i=1;
> dy(1)=1e-1;
> while (abs(dy(i))>1e-3)
> A(i)=b(j)*y(i)+z(j)*(y(i))^2;
> P(i)=b(j)+2*(z(j)^2+1)^(0.50)*y(i);
> R(i)=A(i)/(P(i));
> f(i)=sqrt(s)/n*(A(i)*R(i)^(2/3))-Q; % function
> ff(i)=sqrt(s)/n*((-4/3)*R(i)^(5/3)*sqrt(1+z(j)^2)+5/3*R(i)^
> (2/3)*(b(j)+2*z(j)*y(i))); % f' Derivative
> y(i+1)=y(i)-f(i)/ff(i);
> dy(i+1)=-f(i)/ff(i);
> i=i+1;
> end
> normaldepth(j)=y(i);
> conv(j)=i;
>
> end
>
> for j=1:3 %Side Slope for Trapezoidal, Rectangular, and
> Triangular X-sec respectively.
> %Loop for Critical Depth
> alpha=1;
> theta=0;
> yc(1)=yinitial;
> ic=1;
> dyc(1)=1e-1;
> while (abs(dyc(ic))>1e-3)
>
> Ac(ic)=b(j)*yc(ic)+z(j)*(yc(ic))^2;
> Bc(ic)=b(j)+2*z(j)*yc(ic);
> Pc(ic)=b(j)+2*(z(j)^2+1)^(0.50)*yc(ic);
> Rc(ic)=Ac(ic)/(Pc(ic));
> Dc(ic)=Ac(ic)/(b(j)+2*z(j)*yc(ic));
> fc(ic)=Ac(ic)^(3/2)*Bc(ic)^(-1/2)-Q/sqrt(9.81/alpha); %
> function
> ffc(ic)=Ac(ic)^(3/2)*(-1/2)*Bc(ic)^(-3/2)*2*z(j)+Bc(ic)^(-1/2)*
> (3/2)*Ac(ic)^(1/2)*Bc(ic); % f' Derivative
> yc(ic+1)=yc(ic)-fc(ic)/ffc(ic);
> dyc(ic+1)=-fc(ic)/ffc(ic);
> ic=ic+1;
> end
> criticaldepth(j)=yc(ic);
> convc(j)=ic;
> end
> %% Water Surface Calculations
>
> dx=2.0; %Interval in channel direction
> Y =2.0*normaldepth(1); %Initial Water Depth at the control location
> %Y =1.05*criticaldepth(1);
> %Y=0.05;
> Dir ='U'; %direction of solution : Upstream(U) or
> Downstream(D)?', 's');
> if Dir == 'U';
> dx=-dx;
> end
> L=40000; %total length of the channel
> Nx=round(abs(L/dx));
> BedElev(1)=0; %Begining of Bed Elevation
> WaterElev(1)=Y(1); %Begining of Water Elevation
> b=b(1); % channel bottom width for Trap.
> Section
> z=z(1); % Side Slope for Trap. Section
> for i=1:Nx;
> A(i)=b*Y(i)+z*Y(i)^2; % Area m^2
> B(i)=b+2*z*Y(i); % Top Width
> P(i)=b+2*Y(i)*sqrt(1+z^2); % Wetted Perimeter
> R(i)=A(i)/P(i); % Hydr. radius
> V(i)=Q/A(i); % Velocity
> Sf(i)=n^2*V(i)^2/(R(i)^(4/3)); % Energy slope
> f(i)=(s-Sf(i))/(1-Q^2*B(i)/9.81/A(i)^3); % dy/dx=f=(s-sf/(1-fr^2))
> yy(i)=Y(i)+0.5*f(i)*dx; % first modified depth or
> y(i+1/2)=yold+.5*dx*f'
> A(i)=b*yy(i)+z*yy(i)^2; % Calculate new area but
> using yy which is ymod1
> B(i)=b+2*z*yy(i); % Calculate new top width
> but using yy which is ymod1
> P(i)=b+2*yy(i)*sqrt(1+z^2); % Calculate new wetted
> perimeter but using yy which is ymod1
> R(i)=A(i)/P(i);
> V(i)=Q/A(i);
> Sf(i)=n^2*V(i)^2/(R(i)^(4/3));
> ff(i)=(s-Sf(i))/(1-Q^2*B(i)/9.81/A(i)^3); %Calculate new dy/dx but
> using yy which is ymod1
> Y(i+1)=Y(i)+ff(i)*dx; % Calcualte final y
> which is yold+f'(y(i+1/2))*dx
> BedElev(i+1)=-dx*(i)*s; % Bed elevation
> for plotting purpose
> WaterElev(i+1)=BedElev(i+1)+Y(i+1); % Water
> surface elevation for plotting purpose
> X(i+1)=dx*i; % X: distance along
> the channel in (m)
> end
> CriticalElev=BedElev.*1.0+criticaldepth(1); % get the critical depth
> line for plotting purpose
> NormalElev=BedElev.*1.0+normaldepth(1); % get the normal depth
> line for plotting purpose
> Yfinalm=Y(i+1);
> figure(1);
> plot(X,BedElev,'k-','LineWidth',2); hold;
> plot (X,WaterElev,'LineWidth',2);
> plot(X,CriticalElev,'-.r')
> plot(X,NormalElev,'--m' )
> h = legend('Bed','Water Surface','Critical Depth Line','Normal Depth
> Line',1);
> xlabel('Length(m)');
> ylabel('Elevation(m)');
> title('Water Surface Profile');
> XY=[X' BedElev' WaterElev'];
> d = {'X', 'Bed Elev','WaterElev'};
> dir=pwd;
> xlswrite('water depth.xlsx',d, 'WaterSurface' ,'A1')
> xlswrite('water depth.xlsx',XY,'WaterSurface','A2');

I recommend that you get a hold of Henderson's textbook "Open Channel
Flow". It's very old (1966), but it has all you need to do this
homework assignment.
At least, I found it did, back in 1968 when I had the same assignment
as you have. But we didn't have Matlab back then. Only Fortran,
using Hollerith cards on an IBM mainframe.

|
Pages: 1
Prev: How to use callAllOptimPlotFcns
Next: rayleighchan and mlseeq help....???