MATLAB Answers

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

39 views (last 30 days)
Christopher Van Horn on 14 Jan 2021
Commented: Alan Stevens on 30 Jan 2021
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 CommentsShowHide 1 older comment
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.

### Answers (1)

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')
##### 2 CommentsShowHide 1 older comment
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')

Sign in to comment.

R2020a

### Community Treasure Hunt

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

Start Hunting!