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

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.

[Christopher Van Horn]

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]

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')

Answer 1

[Alan Stevens]

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')

[Christopher Van Horn]

Allen that is a very interesting theory on how to solve a mass balance equation. however if Im not mistaken dv/dt = a flow rate (m^3/s) this necesitates a differential equation for your Q(i) to make it a flow rate out to match the Qin flowrate of 50 (m^3/s). I do like the h(i) = hrmax*(V/Vrmax)^(1/0.3); that solved how to get around using nthroot which whas giving me problems. Here is what I have resolved from your work above. I am also including to sets of theory using sims as explination of how out put flows. However, that your wonderful code above does not solve my height volume problem with sims or ode. I displayed my current code and myode.m and attached the theory in sim for at the end with and without variables. I deply appreciate your imput and encourage you to keep at it.

clc
clear

definitions
% h(t) = depth of water
% Hspill = hight of spillway
% Qin = amount of water into reservoir
% Qout = amout of water exiting reservoir
% g = gravity
% t = time record of flow (Q)
% c1 = orifice coefficient ?
% c2 = weir coefficient
% Ac = orifice cross sectional area
% Lw = length of the wier crest.
% hr_max = height of wier
% Qin = Qout when reservoir is full ((dVr/dt))
% Vr(t) = volume of reservoir
% Vr_max = max volume of reservoir
% r = radius of the pipe
% 
% The following relation between the reservoir inflow, outflow, and a change in water level holds:
% 
%                              dV/dt = Qin(t)-Qb(H)-Qf(H)
% 
% where V is the reservoir water volume (m3), t is time (s), Qin is the reservoir inflow (m3/s), Qb is
% the breach outflow (m3/s) (in the following text, Qbo denotes overtopping discharge and Qbp piping
% discharge), Qf is the outflow through appurtenant works (outlets, spillways, etc.) in (m3/s), and H is
% the elevation of the water level.
% 
% The breach outflow during overtopping QbO was calculated using an overflow formula for a
% rectangular broad crested weir 
% 
%                             QbO = m*b*sqrt(2g)*(H-Z)^1.5
% 
% where m is the weir coeficient, b is the width of the breach opening (m), g is the acceleration due to
% gravity (m/s2), and Z is the breach bottom elevation .
% 
% Using the overtopping discharge obtained from Equation (2), the water depth hf and velocity vf
% along the slope was derived from the Chezy equation :
% 
%                             hf =((Qb0*n/(b*sqrt(sin(beta)))))^.6
% 
%                             vf = (sqrt(sin(beta)/n)*hf^0.67
% 
% where n is the roughness of the downstream slope and beta is the
% downstream slope angle

Pipe
% For the piping model, a circular opening with a constant leakage pipe diameter and constant
% shear stress and hydraulic gradient along the pipe were assumed. The input parameters
% were the length of the pipe L, the initial reservoir water level H0, the elevation of the pipe outflow
% Houtflow, the critical shear stress along the walls of the pipe taoc, the soil erodibility Ce, the initial pipe
% diameter D0, and absolute roughness delta.

% The outflow through the pipe QbP and average water velocity vpip in the pipe were determined
% using pipe hydraulics
% 
%                 QbP = (1/(sprt(alpha + lambda*(L/D))))*S*sqrt(2*g*(H-Houtflow))
% 
%                                       Vpip = Qbp/S
% 
% where alpha is the kinetic energy (Coriolis) coefficient, L is the length of the pipe (m), D is the pipe diameter
% (m), lambda is the friction loss coefficient, S is the pipe cross-sectional area (m2), and Houtflow is the pipe
% outflow elevation
% 
%                             IE = lambda *(1/D)*((Vpip^2)/(2*g)
%
%                                     Tao = rho*g*(D/4)*IE
%
% using the Newton method with a very small-time step of delta t = 1 s.

Creat Arrays
% Time= [0:10*60*60];
% Vr = zeros(size(Time));
% Qout = zeros(size(Time));
% Qin = zeros(size(Time));
% dh = zeros(size(Time));

Variables
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
hspan = [1 6];
tspan = [0 10];

Ac = sect_area_cylinder(Radius,Height,alpha);


Initial conditions

% 
h0 = 1;
% % Qout(1) = double(Qpipe);
Vr = Vr_max*(h0/hr_dam)^0.3
Qin = 50
dVr(1) = Vr+Qin-C_1*Ac*sqrt(2*g*h0)
dh(1) = nthroot(dVr/Vr, 0.3)
Vr0 = Vr_max*(h0/hr_dam)^0.3;

% 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 
z0 = [Vr;h0]

[t,dz] = ode45(@(t,z) myode(t,z),tspan,z0)


figure(1)
plot(t,dz(:,1))
figure(2)
plot(t,dz(:,2))
 
% % 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
% 
% fun = @(t) Qin(t)-Qwier(t)-Qpipe(t);
% dV = integral(fun,0,Time) % dt


function  [dz] = myode(t,z)
%%  definitions

%   vector z is equivalent to (Vout,h)
%   vector g is equivalent to [g1;g2] ~ dz/dt
%   g1 is equivalent dV0/dt
%   g2 is equivalent dh/dt
%   Angular pluse is dv/dt = 50 (m^3/s) for tspan or 10 hours
%% Imput numbers
        g = 9.81; % meters per sec squared
        C_1 = .82;
        Radius = .3; % meters
        Height = .6; % meters
        alpha = 0; 
        C_2 = 1;
        L_w = 3; % meters
        Hspill = 3; % meters
        Ac = sect_area_cylinder(Radius,Height,alpha);
        Qin = 50; % ~ dv/dt on the left hand side derives volumetric flow (m^3/s)
        h0 = 1;
        Vr_max = 1000000; % m^3
        hr_dam = 6; % m
        Vr0 = Vr_max*(h0/hr_dam)^0.3;
        
        
%          dVr(1) = Vr+Qin-C_1*Ac*sqrt(2*g*h0)
%       [t,y,te,ye,ie] = ode45(odefun,tspan,y0,options)  
%         Vr0 = Vr_max*(h0/hr_dam)^0.3;
%         H0  = nthroot(Vr0/Vr0, 0.3);
%%  calculations
    h = z(2);
    Vr = Vr_max*(h/hr_dam)^0.3;
    % The mass conservation statements of the model
   
             if  (0<h)&&(h<=3)
                    g1 =  C_1*Ac*sqrt(2*g*h);
             elseif (3<h)&&(h<=6)
                    g1 = C_1*Ac*sqrt(2*g*h) + C_2*L_w*(h-Hspill).^3/2;
             else
                    g1 = 50;
             end
             if (h<6)
                
                g2 = hr_dam*(Vr/Vr_max).^(1/0.3);
             else
                g2 = 0;
             end
     dz = [g1;g2];
   
end

[Alan Stevens]

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]

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.