trans.m

% SIMULATION OF TRANSIENT STATE BEHAVIOUR OF A CONTINUOUS CULTURE CHEMOSTAT
%
% Chris A. du Plessis
% Department of Microbiology & Biochemistry
% University of the Orange Free State
% Bloemfontein
% South Africa
% E-mail: [e-mail]
%
% DESCRIPTION
%
% M-file to simulate transient state conditions of a chemostat after
% conditions at steady state have been changed. This will simulate conditions
% from one steady state (or from start-up) until a (new) steady state is
% achieved. 
%
% HOW TO RUN THE FILE
%
% Open this M-file from Matlab. In Matlab type the file name of this M-file.
% You will be prompted for certain inputs.  Give the inputs and press enter.
% The assumption is made that the % reactor is already at steady state 
% (or at start-up).  The inputs you are prompted for are, thus, the
% new reactor conditions and/or parameters.  
%
% OUTPUTS
%
% The outputs are shown in graphical form to show the change in conditions,
% since change of reactor conditions/parameters took place.  Three graphs
% are generated.
% Figure 1: Biomass concentration vs Time
% Figure 2: Residual reactor substrate concentration vs Time
% Figure 3: Specific growth rate vs Time  
% The figures show the transient conditions caused by the change in
% paramters or operating conditions and allow visualisation of the new
% steady state conditions.  If new steady state conditions are not
% evident from figures, re-run the M-file and increase the "Total Time"
% and/or the number of "Time Steps".
%
% EXAMPLE VALUES
% Maximum specific growth rate value Umax in 1/hour: 0.4
% Half-saturation constant Ks in mg/L: 200
% Yield coefficient Y: 0.3
% Pump flow rate in L/hour: 0.4
% Volume of the reactor in L: 1
% Concentration of the feed nutrient in mg/L: 20000
% How many time steps do you want?, >1000 recommended: 2000
% Total time period you would like to simulate in hours: 30
% Initial cell conc in the reactor in mg per L: 10000
% Initial residual substrate conc in the reactor in mg per L: 100
% Initial specific growth rate, i.e. dilution rate at equilibrium in 1/h: 0.35
% 
% CODE:
clear all
clear X
clear S
clear u
umax=input ('Type in the maximum specific growth rate value Umax in 1/hour:  ');
Ks=input ('Type in half-saturation constant Ks in mg/L:  ');
Y=input ('Type in the yield coefficient Y:  ');
F=input ('Type in the pump flow rate in L/hour:  ');
V=input ('Type in the volume of the reactor in L:  ');
C=input ('Type in the concentration of the feed nutrient in mg/L:  ');
nm=input ('How many time steps do you want?, >1000 recommended:  ');
Ttot=input ('Type in the total time period you would like to simulate in hours:  ');
dt=Ttot/nm;
X(1)=input ('Initial cell conc in the reactor in mg per L:  ');
S(1)=input ('Initial residual substrate conc in the reactor in mg per L:  ');
u(1)=input ('Initial specific growth rate, i.e. dilution rate at equilibrium in 1/h:  ');
disp ('  ')
disp ('  ')
disp ('  ')
disp ('  ')
disp ('Please wait for calculations and graphs, this might take a minute!')
T(1)=0;
for i = 1:nm-1
D=F/V;
u(i+1)=umax*S(i)/(Ks+S(i));
X(i+1)=X(i)+(dt*X(i)*(u(i)-D));
S(i+1)=S(i)+(dt*((D*(C-S(i)))-(u(i)*X(i)/Y)));
T(i+1)=T(i)+dt;
end
n=0: nm-1;
out=[T;X;S;u];
plot(T,X,'b')
xlabel('Time, hours');
ylabel('Biomass concentration, mg/L');
hold
figure
plot(T,S,'g')
xlabel('Time, hours');
ylabel('Substrate concentration, mg/L');
hold
figure
plot(T,u,'r')
xlabel('Time, hours');
ylabel('Specific growth rate, 1/h');
hold