Water Tank input that changes based on volume.

Need to plot the volume of water and concentration of pollutant in a tank of water. Initial inflow is 1.25 m^3/s, and outflow is 1.1. When the tank's volume reaches 200 m^3, the inflow is 0, and returns to 1.25 when the volume drops to 150. How can inflow (Qin) be changed to 0 at V=200, and then reset at V=150? if statement?

Here is code for a prior part to the project that assumes a constant rate of fill:

V0 = 50;             %initial volume (m^3)
c0 = 0;              %initial concentration
Qin = 1.25 * 60;   %inflow (m^3/hour)
Qout = 1.1 * 60;     %outflow (m^3/hour)
cin = 100;           %pollutant concentration (mg/L)
t_final = 50;        %hours
dt = 0.5;            %time step
T = [dt:0.5:50];
t = 0;
%
k = 0;
%
%Initialize Loop Variables
cold = c0;
vold = V0;

for n = [1 : (t_final/dt)];
    t = t + dt;
    k = 0;
    
    alpha = (Qin / (vold + (Qin - Qout) * t)) + k;
    beta = (Qin*cin)/(vold+(Qin - Qout) * t);
    
    c(n) = cold - dt * alpha * cold + dt * beta;
    V(n) = vold + (Qin - Qout) * dt;
    
    cold = c(n);
    vold = V(n);
end


yyaxis left
plot(T,V)
ylabel('Volume (m^3)')
xlabel('Time (hours)')

yyaxis right
plot(T,c)
ylabel('Concentration (mg/L)')

disp('      t        V           c')
format compact
disp([T' V' c'])
more(10)

Answer 1

[Subhamoy Saha]

Please check with the following

clear, clc
V0 = 50;        %initial volume (m^3)
c0 = 0;     %initial concentration
Qin = 1.25 * 60;    %inflow (m^3/hour)
Qout = 1.1 * 60;    %outflow (m^3/hour)
cin = 100;          %pollutant concentration (mg/L)
t_final = 50;       %hours
dt = 0.5;           %time step
T = [dt:0.5:50];
t = 0;
%
k = 0;
%
%Initialize Loop Variables
cold = c0;
vold = V0;
inflow=1; % inflow status variable

for n = [1 : (t_final/dt)]
    t = t + dt;
    k = 0;
    
    if vold>=200 || inflow==0
        inflow=0; % say inflow stopped
        if vold<=150
            inflow=1;  % say inflow start
        end
    end
    
    if inflow==0
        Qin=0;
    else
        Qin=1.25 *60;
    end
    
    alpha = (Qin / (vold + (Qin - Qout) * t)) + k;
    beta = (Qin*cin)/(vold+(Qin - Qout) * t);

c(n) = cold - dt * alpha * cold + dt * beta;
V(n) = vold + (Qin - Qout) * dt;
 
cold = c(n);
vold = V(n);

end
yyaxis left
plot(T,V)
ylabel('Volume (m^3)')
xlabel('Time (hours)')
yyaxis right
plot(T,c)
ylabel('Concentration (mg/L)')
disp('      t        V           c')
format compact
disp([T' V' c'])
more(10)

[Karl Nilsson]

Thanks. I figured it out eventually but this is much more efficient.