%
% Exercise 7: Carbon Removal Process
%
% - three constant volume tanks
% - well mixed
% - Monod kinetics
% - constant biomass
%
% Bob Newell, May 1996
%
%-----------------------------------------------------------
% parameters
muhat = 3.733;
Ks = 20;
Koh = 0.2;
Yh = 0.6;
bh = 0.4;
fp = 0.1;
a = 332;
b = 16;
So_sat = 10;
%
vol1 = 4;
vol2 = 4;
vol3 = 4;
air1 = 6.75;
air2 = 1.60;
air3 = 0.93;
flow = 100;
cod = 75;
%
f = input( 'Flowrate to Tank 1 [100] > ' );
if ~isempty(f), flow = f; end;
f = input( 'COD to Tank 1 [75] > ' );
if ~isempty(f), cod = f; end;
%
air = input( 'Air to Tank 1 [6.75] > ' );
if ~isempty(air), air1 = air; end;
air = input( 'Air to Tank 2 [1.60] > ' );
if ~isempty(air), air2 = air; end;
air = input( 'Air to Tank 3 [0.93] > ' );
if ~isempty(air), air3 = air; end;
%
kla1 = a * ( 1 - exp( - air1 / b ) );
kla2 = a * ( 1 - exp( - air2 / b ) );
kla3 = a * ( 1 - exp( - air3 / b ) );
%
global p
p = [ muhat Ks Koh Yh bh fp a b So_sat ...
             vol1 vol2 vol3 kla1 kla2 kla3 ...
             flow cod ];
%-----------------------------------------------------------
% initial states
Ss1 = 12.18;
So1 = 3;
Ss2 = 2.49;
So2 = 3;
Ss3 = 1.44;
So3 = 3;
%
xo = [ Ss1 So1 Ss2 So2 Ss3 So3 ];
%-----------------------------------------------------------
% integrate the state equation
to = clock;
days = 6/24;
[ t, x ] = ode45( 'Carbon_Plant', [0 days], xo );
secs = etime( clock, to )
%-----------------------------------------------------------
% evaluate the outputs
y = Carbon_Plant_Outputs( t, x );
%-----------------------------------------------------------
% plot the results
figure;
subplot( 2, 1, 1 );
plot( 24*t, y(:,5:7) );
axis( [ 0 24*days 0 80 ] );
ylabel( 'COD (mg/l)' )
legend( 'Tank1', 'Tank2', 'Tank3', -1 )
%
subplot( 2, 1, 2 );
plot( 24*t, y(:,8:10) );
axis( [ 0 24*days 0 8 ] );
xlabel( 'Time (hrs)' )
ylabel( 'DO (mg/l)' )
legend( 'Tank1', 'Tank2', 'Tank3', -1 )
%
[ n, m ] = size( y );
fprintf( 1, 'COD = %7.1f%7.1f%7.1f\n', ...
            y(n,5), y(n,6), y(n,7) );
fprintf( 1, 'DO  = %7.1f%7.1f%7.1f\n', ...
            y(n,8), y(n,9), y(n,10) );
%
%-----------------------------------------------------------
% the end!   
