differential equations of the type Q(t,h(t) where h(t) is based on volume

11 views (last 30 days)
Problem 1.
1. Write a differential equation describing mass conservation in a reservoir for a generic
streamflow input Qin and outlet flow determined by a pipe at the bottom of the reservoir and an
emergency spillway (see figure and equations below).
where c1 is orifice coefficient, Ac is orifice cross-sectional area, c2 is weir coefficient; Lw is the
length of the weir crest. Assume that the outflow is equal to the inflow when the reservoir is full.
Also, the relation between the water volume in the reservoir and the water surface elevation is as
follows:
2. Solve the equation in Matlab using as the input hydrograph a rectangular pulse with flow rate
of 50 m3/s and duration of 10 hours, and let hR-Max=6m, circular orifice size=0.6 m, c1=0.82 (assume
a short pipe), c2=1 (assume a broad-crested weir), VR-Max=1,000,000 m3, Hspill=3m, Lw=3m, and
assume that at time zero h(0)=1m. Plot the input and the output hydrographs.
I have been at this for about 4 months. I have read all the mathworks files and several University files on this subject. Yet every time I try to set dv/dt = ht/dt I get tspan errors. I even asked support for help and got notta. The files attached are not all of my attempts just the clossest attempts and the original pdf question. Please if you cite a file location accompony it with actual code that relates to my question. I am willing to bet that I have already read it and do not understand it. I am a nontraditional student and a disabled veteran.
  2 Comments
Christopher Van Horn
Christopher Van Horn on 20 Jan 2021
mass conservation says that in = out. with an initial height of 1 meter Vr0 = Vr_max*(h0/hr_dam)^0.3 == 5.8419e+05
which is the volume in the reservor at t=0. we also know that dv/dt = the give 50 m^3/s in flow minus the q(t,h(t)) of the pipe and the wier. dh/dt is solvable by taking Vr0 and changing t for dh which produces dh = nthroot((VR+Q_in-Qout)/Vr0, 0.3) so now for any volume we have the coresponding h value.
%% Reservoir Mass Conservation Equation (below spillway)
% dV/dt = Vr-max(h/hr_dam)^0.3
Vr0 = Vr_max*(h0/hr_dam)^0.3
dh = h0;
for t=0:.01:10
Vr = Vr_max*(dh/hr_dam)^0.3;
Qpipe =diff(C_1*Ac*sqrt(2*g*dh));
Qwier =diff((C_1*Ac*sqrt(2*g*dh)+C_2*L_w*(dh-Hspill).^3/2));
Q_out = Qpipe + Qwier;
dh = nthroot((Vr+Q_in-Q_out)/Vr0, 0.3);
end
is where I am at right now It does not work. I am trying to understand how to equate all the givens. Any help would of course be apriciated
Christopher Van Horn
Christopher Van Horn on 29 Jan 2021
still trying to figure out how to tie these together any help would be wonderful
clc
clear
reset(symengine)
% syms h real
% syms t real
% syms t clear
% syms h clear
Radius = 0.3;
Height = 0.6;
alpha = 0.0;
Vr_max = 1000000; % m^3
hr_dam = 6; % m
C_1 = 0.82; %
g = 9.81; %m^2/s
C_2 = 1.00;
Hspill = 3; % m
L_w = 3; % m
hr_max = 6; % Max water depth in meters
h0 = 1;
hspan = [1 6];
tspan = [0 10];
Time= 10*60*60; % converts hours to seconds
Ac = sect_area_cylinder(Radius,Height,alpha);
Vr = Vr_max*(h0/hr_dam)^0.3
Qin = 50
syms h t
Qpipe = C_1*Ac*sqrt(2*g*h)
Qp = diff(Qpipe)
Qwier = C_1*Ac*sqrt(2*g*h)+C_2*L_w*(h-Hspill)^1.5
Qw = diff(Qwier)
QW = matlabFunction(diff(Qwier))
QP = matlabFunction(diff(Qpipe))
% Qoutp = QP(1:.01:3)
% Qoutw = QW(3.01:.01:6)
% Qout = piecewise(h<=0, QP, 3<h<=6, QW)
% Ht =matlabFunction (Qout(h)) %,'vars',[h]
Qout = piecewise(h<=0, C_1*Ac*sqrt(2*g*h), 3<h<=6, C_1*Ac*sqrt(2*g*h)+C_2*L_w*(h-Hspill)^1.5)
QoutH=diff(Qout, h)
fplot(Qout)
fplot(QoutH)
matlabFunction(QoutH,'file','C:\Users\cdrlj\Documents\Spring 2021\final\testMatrix.m')

Sign in to comment.

Accepted Answer

Alan Stevens
Alan Stevens on 30 Jan 2021
Does this help:
g = 9.81; % m/s^2
c1 = 0.82; c2 = 1;
Hspill = 3; % m
Lw = 3; % m
Ac = pi*0.6^2/4; % m^2
Vrmax = 10^6; % m^3
hrmax = 6; % m
Qin = 50; % m^3/s
h0 = 1; % m
V0 = Vrmax*(h0/hrmax)^0.3; % m^3
tend = 10*3600; % secs
dt = 1; % secs
t = 0:dt:tend; % secs
h = zeros(1,numel(t));
h(1) = h0;
Q = zeros(1,numel(t));
Q(1) = c1*Ac*sqrt(2*g*h0);
for i = 2:numel(t)
V = min(V0 + i*Qin*dt,Vrmax);
h(i) = hrmax*(V/Vrmax)^(1/0.3);
if h(i)<=Hspill
Q(i) = c1*Ac*sqrt(2*g*h(i)); % m^3/s
elseif h(i)>Hspill && V<Vrmax
Q(i) = c1*Ac*sqrt(2*g*h(i)) + c2*Lw*(h(i) - Hspill)^(3/2);
else
Q(i) = Qin;
end
end
plot(t/3600,Q),grid
xlabel('time [hours]'),ylabel('flow rate [m^3/s]')
figure
plot(t/3600,h,[0 tend/3600],[Hspill Hspill],'--'),grid
xlabel('time [hours]'),ylabel('height [m]')
legend('h','Hspill')
  3 Comments
Alan Stevens
Alan Stevens on 30 Jan 2021
The following explains my conservation of mass reasoning (though I did forget to incorporate the exit flow correctly!). The time integration I used was essentially a simple Euler integration.
My corrected code is as follows:
h0 = 1; % m
V0 = Vrmax*(h0/hrmax)^0.3; % m^3
tend = 10*3600; % secs
dt = 1; % secs
t = 0:dt:tend; % secs
h = zeros(1,numel(t));
h(1) = h0;
Q = zeros(1,numel(t));
Q(1) = c1*Ac*sqrt(2*g*h0);
V = V0;
for i = 2:numel(t)
H = h(i-1);
if H<=Hspill
Q(i) = c1*Ac*sqrt(2*g*H); % m^3/s
elseif H>Hspill && V<Vrmax
Q(i) = c1*Ac*sqrt(2*g*H) + c2*Lw*(H - Hspill)^(3/2);
else
Q(i) = Qin;
end
dV = (Qin-Q(i))*dt; % Volume change in time dt
V = min(V + dV, Vrmax); % New volume
h(i) = hrmax*(V/Vrmax)^(1/0.3); % New height
end
plot(t/3600,Q),grid
xlabel('time [hours]'),ylabel('flow rate [m^3/s]')
figure
plot(t/3600,h,[0 tend/3600],[Hspill Hspill],'--'),grid
xlabel('time [hours]'),ylabel('height [m]')
legend('h','Hspill')
Christopher Van Horn
Christopher Van Horn on 14 May 2021
Sorry it took so long to get back to this. I appriciate your patience and assistance. Steel structural annalysis gave me some fits..lol.. this semester.

Sign in to comment.

More Answers (0)

Categories

Find more on Historical Contests in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!