clc;close all
%Project to simulate the performance of the reservoir as a sedimentation
%tank
%1/31/2021
%Given
simtime = 15050; %seconds
n = 500; %number of time steps
dt = simtime/(n); %time step
t = linspace (0,15050,n); %time
d = 3; %orifice diameter in ft
alpha = 0.7 ; %inflow concentration coefficient
beta = 0.64 ; %inflow concentration coefficient
r = 9.5*10^-4 ; %rate of sediment deposition
density = 102 ; %lb/ft^3 sediment density
g = 32.2; %gravity acceleration
cd = .65 ; %discharge coefficient
a = pi*(d.^2)/4; %cross sectional AREA of pipe
AREA = 600*600 ; %floor area of storage
%% set up arrays
%Initial conditions
simSize = size(t);
qi = zeros(simSize);
ci = zeros(simSize);
qo = zeros(simSize);
c = zeros(simSize);
s = zeros(simSize);
cs = zeros(simSize);
h = zeros(simSize);
m = zeros(simSize);
%Time loop
for k=2: length(t)
newTime = t(k); %Get the new time
qiNew = dischargeIn(newTime); %Get the new fluid inflow
if newTime > 9000
qi(k) = 0;
else
qi(k) = qiNew; %Insert new fluid inflow into inflow
%array
end
ciNew= flowConcentration(alpha,beta,qiNew); %Concentration in flow
ci(k) = ciNew; %Insert new concentration inflow array
hOld = h(k-1); %get old height
qoNew = dischargeOut(cd,g,a,hOld); %use old height to calc outflow
qo(k) = qoNew; %insert into array
qiOld = qi(k-1);
ds = storageStep(dt,qiOld,qiNew,qoNew); %calc change in storage
sOld = s(k-1); %Get old storage volume
sNew = sOld + ds; % Get new storage volume
s(k) = sNew; % Store new volume
hNew = sNew/AREA; %Get new height
h(k) = hNew; %store new height
csOld = cs(k-1);
ciOld = ci(k-1);
cOld = c(k-1);
qoOld = qo(k-1);
csNew = CSTR(cOld,ciOld,csOld,qiOld,qoOld,dt,r); %calcuate new cs
cs(k) = csNew; %store
cNew = concentration(csNew,sNew);
c(k) = cNew;
dm = ciNew * ds;
m(k) = m(k-1) + dm;
end
%check
mi = trapz(t,qi);
MassIn = sum(mi);
mo = trapz(t,qo);
MassOut = sum(mo);
rcs = r*cs;
rm = trapz(t,rcs);
ResidedMass = sum(rm);
TotalMass = MassIn - MassOut - ResidedMass;
figure(1)
plot(t,ci,t,c)
figure(2)
plot(t,m)
%% functions
function ds = storageStep(dt,qiOld,qiNew,qo)
avg = (qiOld + qiNew)/2;
ds = (avg - qo) * dt;
end
function qo = dischargeOut(cd,g,a,h)
qo = cd*a*sqrt(2*g)*h^.5;
end
function qi = dischargeIn(t)
qi = 750 / pi * (1 - cos ( pi * t /4500 )) ;
end
function ci = flowConcentration (alpha,beta,qi)
ci = alpha * qi^beta;
end
function c = concentration(cs,s)
c = cs/s;
end
function csnew =CSTR(c,ci,cs,qi,qo,dt,r)
csnew = cs + ((ci*qi)- (c*qo)-(r*cs))*dt;
end
