Ethanol Bioreactor
Bioreactors create a biologically active environment for the production of chemicals. A bioreactor is a vessel where organisms are grown and preferential products are produced by controlling the feed and temperature. Bioreactors can either be aerobic or anaerobic.
One particular batch bioreactor is fed at varying rates throughout the batch process to produce ethanol. Algebraic equations define the heat generation, rate constants, volume, and other factors used in the model. Symbols and units are defined in the source code.
Growth Rate
$$\mu=\mu_{\max} \frac{{S}}{{K}_{{SX}}+{S}} \frac{{O}_{{liq}}}{{K}_{{OX}}+{O}_{{liq}}}\left(1-\frac{{P}}{{P}_{\max}}\right) \frac{1}{1+\exp (-(100-{S}))}$$
$$\mu_{\max}=\left[\left({a}_1\left({~T}-{k}_1\right)\right)\left(1-\exp \left({b}_1\left({~T}-{k}_2\right)\right)\right)\right]^2$$
$${P}_{\max}={P}_{\operatorname{maxb}}+\frac{{P}_{\operatorname{maxT}}}{1-\exp \left(-{b}_2\left({~T}-{k}_3\right)\right)}$$
$${q}_{{P}}={a}_{{P}} \mu+{b}_{{P}}$$
Non-Growth Ethanol Production
$${b}_{{P}}={c}_1 \exp \left(-\frac{{A}_{{P} 1}}{{~T}}\right)-{c}_2 \exp \left(-\frac{{A}_{{P} 2}}{{~T}}\right)$$
Ethanol Consumption Rate
$${q}_{{S}}=\frac{\mu}{{Y}_{{XS}}}+\frac{{q}_{{P}}}{{Y}_{{PS}}}$$
Oxygen Consumption Rate
$${q}_{{O}}=\frac{{q}_{{O}, \max}}{{Y}_{{Xo}}} \frac{{O}_{{liq}}}{{K}_{{OX}}+{O}_{{liq}}}$$
Biological Cell Mass Deactivation Rate
$${K}_{{d}}={K}_{{db}}+\frac{{K}_{{dT}}}{1+\exp \left(-{b}_3\left({~T}-{k}_4\right)\right)}$$
Oxygen Saturation Concentration
$${O}^*=\frac{{zO}_{{gas}} {RT}}{{K}_{{H}}}$$
Oxygen Mass Transfer Coefficient
$${k}_{l} {a}=\left({k}_{l} {a}\right)_0(1.2)^{{T}-20}$$
Liquid and Vapor Volumes
$${V}_{{l}}+{V}_{{g}}={V}$$
Differential equations are from material, species, and energy balances. The differential equations relate the input flow and substrate concentrations to ethanol production.
Liquid Volume
$$\frac{d V_1}{d t}=Q_{i n}-Q_e$$
Total Cell Mass
$$\frac{{dX}}{{dt}}=\frac{{Q}_{{in}}}{{V}_1}\left({X}_{{t}, {in}}-{X}_{{t}}\right)+\mu {X}_{{v}}$$
Total Biologically Active Cell Mass
$$\frac{d X_v}{d t}=\frac{Q_{i n}}{V_1}\left(X_{v, i n}-X_v\right)+\left(\mu-K_d\right) X_v$$
Glucose (Substrate) Concentration
$$\frac{{dS}}{{dt}}=\frac{{Q}_{{in}}}{{V}_{{l}}}\left({S}_{{in}}-{S}\right)-{q}_{{S}} {X}_{{v}}$$
Ethanol (Product) Concentration
$$\frac{{dP}}{{dt}}=\frac{{Q}_{{in}}}{{V}_1}\left({P}_{{in}}-{P}\right)+{q}_{{P}} X_{{v}}$$
Liquid Oxygen Concentration
$$\frac{{dO}_{liq}}{{dt}}=\frac{{Q}_{{in}}}{{V}_{{l}}}\left({O}^*-{O}_{{liq}}\right)+\left({k}_{{l}} {a}\right)\left({O}^*-{O}_{{liq}}\right)-{qo}_{{o}} {X}_{{v}}$$
Gas Oxygen Concentration
$$\frac{d O_{{gas}}}{d t}=\frac{F_{{air}}}{V_g}\left(O_{{gas}, { in}}-O_{{gas}}\right)-\frac{V_l\left(k_l a\right)}{V_g}\left(O^*-O_{{liq}}\right)+O_{{gas}} \frac{Q_{{in}}-Q_e}{V_g}$$
Bioreactor Temperature
$$\frac{{dT}}{{dt}}=\frac{{Q}_{{in}}}{{V}_{{l}}}\left({T}_{{in}}-{T}\right)-\frac{{T}_{{ref}}}{{V}_{{l}}}\left({Q}_{{in}}-{Q}_{{e}}\right)+\frac{{q}_{{o}} {X}_{{v}} \Delta {H}}{{MW}_{{O}} \rho {C}_{{p}, {br}}}-\frac{{K}_{{T}} {A}_{{T}}\left({T}-{T}_{{c}}\right)}{{V}_{{l}} \rho {C}_{{p}, {br}}}$$
Cooling Agent Temperature
$$\frac{{dT}_{{c}}}{{dt}}=\frac{{F}_{{c}}}{{V}_{{cj}}}\left({T}_{{c}, {in}}-{T}_{{c}}\right)+\frac{{K}_{{T}} {A}_{{T}}\left({T}-{T}_{{c}}\right)}{{V}_{{cj}} \rho_{{c}} {C}_{{p}, {c}}}$$
The model simulates biomass, ethanol, and by-product production. The following chart show total biomass and viable (living) biomass in the reactor.
The model can be used to investigate advanced process control to maximize performance subject to a feeding strategy and measured disturbances.
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
m = GEKKO(remote=False)
# Create time vector: t=[0, 0.1, 0.2,...,36.9,37], [hours]
tm = np.linspace(0,37,371)
# Insert smaller time steps at the beginning
tm = np.insert(tm,1,[0.001,0.005,0.01,0.05])
m.time = tm
nt = len(tm)
# Define constants and parameters
#################################
# Kinetic Parameters
a1 = 0.05 # Ratkowsky parameter [oC-1 h-0.5]
aP = 4.50 # Growth-associated parameter for EtOh production [-]
AP1 = 6.0 # Activation energy parameter for EtOh production [oC]
AP2 = 20.3 # Activation energy parameter for EtOh production [oC]
b1 = 0.035 # Parameter in the exponential expression of the maximum specific growth rate expression [oC-1]
b2 = 0.15 # Parameter in the exponential expression of the maximum specific growth rate expression [oC-1]
b3 = 0.40 # Parameter in the exponential expression of the specific death rate expression [oC-1]
c1 = 0.38 # Constant decoupling factor for EtOh [gP gX-1 h-1]
c2 = 0.29 # Constant decoupling factor for EtOh [gP gX-1 h-1]
k1 = 3 # Parameter in the maximum specific growth rate expression [oC]
k2 = 55 # Parameter in the maximum specific growth rate expression [oC]
k3 = 60 # Parameter in the growth-inhibitory EtOH concentration expression [oC]
k4 = 50 # Temperature at the inflection point of the specific death rate sigmoid curve [oC]
Pmaxb = 90 # Temperature-independent product inhibition constant [g L-1]
PmaxT = 90 # Maximum value of product inhibition constant due to temperature [g L-1]
Kdb = 0.025 # Basal specific cellular biomass death rate [h-1]
KdT = 30 # Maximum value of specific cellular biomass death rate due to temperature [h-1]
KSX = 5 # Glucose saturation constant for the specific growth rate [g L-1]
KOX = 0.0005 # Oxygen saturation constant for the specific growth rate [g L-1]
qOmax = 0.05 # Maximum specific oxygen consumption rate [h-1]
# Metabolic Parameters
YPS = 0.51 # Theoretical yield of EtOH on glucose [gP gS-1]
YXO = 0.97 # Theoretical yield of biomass on oxygen [gX gO-1]
YXS = 0.53 # Theoretical yield of biomass on glucose [gX gS-1]
# Physicochemical and thermodynamic parameters
Chbr = 4.18 # Heat capacity of the mass of reaction [J g-1 oC-1]
Chc = 4.18 # Heat capacity of cooling agent [J g-1 oC-1]
deltaH = 518000 # Heat of reaction of fermentation [J mol-1 O2]
Tref = 20 # Reference temperature [oC]
KH = 200 # Henry's constant for oxygen in the fermentation broth [atm L mol-1]
z = 0.792 # Oxygen compressibility factor [-]
R = 0.082 # Ideal gas constant [L atm mol-1 oC-1]
kla0 = 100 # Temperature-independent volumetric oxygen transfer coefficient [-h]
KT = 360000 # Heat transfer coefficient [J h-1 m-2 oC-1]
rho = 1080 # Density of the fermentation broth [g L-1]
rhoc = 1000 # Density of the cooling agent [g L-1]
MO = 32.0 # Molecular weight of oxygen [g mol-1]
# Bioreactor design data
AT = 1 # Bioreactor heat transfer area [m2]
V = 1800 # Bioreactor working volume [L]
Vcj = 50 # Cooling jacket volume [L]
Ogasin = 0.305 # Oxygen concentration in airflow inlet [g L-1]
# Define variables
##################
mi = m.Var(name='mi',lb=0)
# Qin is a step function:
# Qin = Qin0 + 15H(t-5) + 5H(t-10) - 6H(t-20) - 14H(t-35)
# where H(t-t0) heaviside function
Qin_step = np.zeros(nt)
Qin_step[np.where(tm>=5)] += 15
Qin_step[np.where(tm>=10)] += 5
Qin_step[np.where(tm>=20)] -= 6
Qin_step[np.where(tm>=35)] -= 14
Qin = m.Param(value=Qin_step, name='Qin')
# Fixed variables, they are constant throughout the time horizon
Xtin = m.FV(value=0, name='Xtin')
Xvin = m.FV(value=0, name='Xvin')
Qe = m.FV(value=0, name='Qe')
Sin = m.FV(value=400, lb=0, ub=1500)
Pin = m.FV(value=0, name='Pin')
Fc = m.FV(value=40, name='Fc')
Fair = m.FV(value=60000, name='Fair')
Tin = m.FV(value=30, name='Tin')
Tcin = m.FV(value=15, name='Tcin')
Vl = m.Var(value=1000, name='Vl') # lb=-0.0, ub=0.75*V
Xt = m.Var(value=0.1, name='Xt') # lb=-0.0, ub=10
Xv = m.Var(value=0.1, name='Xv') # lb=-0.0, ub=10
S = m.Var(value=50, name='S') # lb=+0.0, ub=10000
P = m.Var(value=0, name='P')
Ol = m.Var(value=0.0065, name= 'Ol')
Og = m.Var(value=0.305, name='Og')
T = m.Var(value=30, name='T') # lb=20, ub=40
Tc = m.Var(value=20, name='Tc') # lb=0, ub=30
Sf_cum = m.Var(value=0, name='Sf_cum')
# Define algebraic equations
############################
# Specific growth rate of cell mass
mimax = m.Intermediate(((a1*(T-k1))*(1-m.exp(b1*(T-k2))))** 2)
Pmax = m.Intermediate(Pmaxb + PmaxT/(1-m.exp(-b2*(T-k3))))
m.Equation(mi == mimax * (S / (KSX+S)) * (Ol/(KOX + Ol)) \
* (1 - P/Pmax) * (1 / (1+m.exp(-(100-S)))))
# Specific production rate of EtOH
bP = m.Intermediate(c1*m.exp(-AP1/T) - c2*m.exp(-AP2/T))
qP = m.Intermediate(aP*mi + bP)
# Specific consumption rate of glucose
qS = m.Intermediate(mi/YXS + qP/YPS)
# Specific consumption rate of oxygen
qO = m.Intermediate(qOmax*Ol/YXO/(KOX+Ol))
# Specific biological deactivation rate of cell mass
Kd = m.Intermediate(Kdb + KdT/(1+m.exp(-b3*(T-k4))))
# Saturation concentration of oxygen in culture media
Ostar = m.Intermediate(z*Og*R*T/KH)
# Oxygen mass transfer coefficient
kla = m.Intermediate(kla0*1.2**(T-20))
# Bioreactor phases equation
Vg = m.Intermediate(V - Vl)
# Define differential equations
###############################
m.Equation(Vl.dt() == Qin - Qe)
m.Equation(Vl*Xt.dt() == Qin*(Xtin-Xt) + mi*Vl*Xv)
m.Equation(Vl*Xv.dt() == Qin*(Xvin-Xv) + Xv*Vl*(mi-Kd))
m.Equation(Vl*S.dt() == Qin*(Sin-S) - qS*Vl*Xv)
m.Equation(Vl*P.dt() == Qin*(Pin - P) + qP*Vl*Xv)
m.Equation(Vl*Ol.dt() == Qin*(Ostar-Ol) + Vl*kla*(Ostar-Ol) - qO*Vl*Xv)
m.Equation(Vg*Og.dt() == Fair*(Ogasin-Og) - Vl*kla*(Ostar-Ol) + Og*(Qin-Qe))
m.Equation(Vl*T.dt() == Qin*(Tin-T) - Tref*(Qin-Qe) \
+ Vl*qO*Xv*deltaH/MO/rho/Chbr - KT*AT*(T-Tc)/rho/Chbr)
m.Equation(Vcj*Tc.dt() == Fc*(Tcin - Tc) + KT*AT*(T-Tc)/rhoc/Chc)
m.Equation(Sf_cum.dt() == Qin*Sin)
# solve ODE
m.options.SOLVER= 1
m.options.IMODE = 7
m.options.NODES = 3
# m.open_folder()
m.solve(disp=False)
# Plot results
plt.figure(1)
plt.title('Total & Viable Cellular Biomass')
plt.plot(m.time, Xv.value, label='Xv')
plt.plot(m.time, Xt.value, label='Xt')
plt.legend()
plt.ylabel('Biomass concentration [g/L]')
plt.xlabel('Time [h]')
plt.grid()
plt.minorticks_on()
plt.ylim(0)
plt.xlim(m.time[0],m.time[-1])
plt.tight_layout()
plt.figure(2)
plt.title('Substrate (S) & Product (P) concentration')
plt.plot(m.time, S.value, label='S')
plt.plot(m.time, P.value, label='P')
plt.legend(); plt.grid()
plt.ylabel('Conc [g/L]')
plt.xlabel('Time [h]')
plt.minorticks_on()
plt.ylim(0)
plt.xlim(m.time[0],m.time[-1])
plt.tight_layout()
plt.figure(3)
plt.title('Bioreactor & Cooling jacket temperature')
plt.plot(m.time, T.value, label='T')
plt.plot(m.time, Tc.value, label='Tc')
plt.legend()
plt.ylabel('Temperature [oC]')
plt.xlabel('Time [h]')
plt.grid()
plt.minorticks_on()
plt.ylim(0)
plt.xlim(m.time[0],m.time[-1])
plt.tight_layout()
fig4, ax = plt.subplots()
ax.title.set_text('Dissolved & Gaseous Oxygen concentration')
lns1 = ax.plot(m.time, Ol.value, label='[Oliq]', color='c')
ax.set_xlabel('Time [h]')
ax.set_ylabel('Oliq [g/L]', color='c')
ax.minorticks_on()
ax2 = ax.twinx()
lns2 = ax2.plot(m.time, Og.value, label='[Ogas]', color='y')
ax2.set_ylabel('Ogas [g/L]', color='y')
ax2.minorticks_on()
lns = lns1 + lns2
labs = [l.get_label() for l in lns]
ax.legend(lns, labs, loc='best')
ax.grid()
fig4.tight_layout()
plt.figure(4)
plt.figure(5)
plt.title('Feeding Policy')
plt.plot(m.time, Qin.value, label='Qin')
plt.legend()
plt.ylabel('Qin [L/h]')
plt.xlabel('Time [h]')
plt.grid()
plt.minorticks_on()
plt.ylim(0)
plt.xlim(m.time[0],m.time[-1])
plt.tight_layout()
plt.figure(6)
plt.title('Check >=0 Constraints')
plt.plot(tm,bP.value,label='bP')
plt.plot(tm,mi.value,label='mi')
plt.legend(); plt.grid()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
m = GEKKO(remote=False)
# Create time vector: t=[0, 0.1, 0.2,...,36.9,37], [hours]
tm = np.linspace(0,37,371)
# Insert smaller time steps at the beginning
tm = np.insert(tm,1,[0.001,0.005,0.01,0.05])
m.time = tm
nt = len(tm)
# Define constants and parameters
#################################
# Kinetic Parameters
a1 = 0.05 # Ratkowsky parameter [oC-1 h-0.5]
aP = 4.50 # Growth-associated parameter for EtOh production [-]
AP1 = 6.0 # Activation energy parameter for EtOh production [oC]
AP2 = 20.3 # Activation energy parameter for EtOh production [oC]
b1 = 0.035 # Parameter in the exponential expression of the maximum specific growth rate expression [oC-1]
b2 = 0.15 # Parameter in the exponential expression of the maximum specific growth rate expression [oC-1]
b3 = 0.40 # Parameter in the exponential expression of the specific death rate expression [oC-1]
c1 = 0.38 # Constant decoupling factor for EtOh [gP gX-1 h-1]
c2 = 0.29 # Constant decoupling factor for EtOh [gP gX-1 h-1]
k1 = 3 # Parameter in the maximum specific growth rate expression [oC]
k2 = 55 # Parameter in the maximum specific growth rate expression [oC]
k3 = 60 # Parameter in the growth-inhibitory EtOH concentration expression [oC]
k4 = 50 # Temperature at the inflection point of the specific death rate sigmoid curve [oC]
Pmaxb = 90 # Temperature-independent product inhibition constant [g L-1]
PmaxT = 90 # Maximum value of product inhibition constant due to temperature [g L-1]
Kdb = 0.025 # Basal specific cellular biomass death rate [h-1]
KdT = 30 # Maximum value of specific cellular biomass death rate due to temperature [h-1]
KSX = 5 # Glucose saturation constant for the specific growth rate [g L-1]
KOX = 0.0005 # Oxygen saturation constant for the specific growth rate [g L-1]
qOmax = 0.05 # Maximum specific oxygen consumption rate [h-1]
# Metabolic Parameters
YPS = 0.51 # Theoretical yield of EtOH on glucose [gP gS-1]
YXO = 0.97 # Theoretical yield of biomass on oxygen [gX gO-1]
YXS = 0.53 # Theoretical yield of biomass on glucose [gX gS-1]
# Physicochemical and thermodynamic parameters
Chbr = 4.18 # Heat capacity of the mass of reaction [J g-1 oC-1]
Chc = 4.18 # Heat capacity of cooling agent [J g-1 oC-1]
deltaH = 518000 # Heat of reaction of fermentation [J mol-1 O2]
Tref = 20 # Reference temperature [oC]
KH = 200 # Henry's constant for oxygen in the fermentation broth [atm L mol-1]
z = 0.792 # Oxygen compressibility factor [-]
R = 0.082 # Ideal gas constant [L atm mol-1 oC-1]
kla0 = 100 # Temperature-independent volumetric oxygen transfer coefficient [-h]
KT = 360000 # Heat transfer coefficient [J h-1 m-2 oC-1]
rho = 1080 # Density of the fermentation broth [g L-1]
rhoc = 1000 # Density of the cooling agent [g L-1]
MO = 32.0 # Molecular weight of oxygen [g mol-1]
# Bioreactor design data
AT = 1 # Bioreactor heat transfer area [m2]
V = 1800 # Bioreactor working volume [L]
Vcj = 50 # Cooling jacket volume [L]
Ogasin = 0.305 # Oxygen concentration in airflow inlet [g L-1]
# Define variables
##################
mi = m.Var(name='mi',lb=0)
# Qin is a step function:
# Qin = Qin0 + 15H(t-5) + 5H(t-10) - 6H(t-20) - 14H(t-35)
# where H(t-t0) heaviside function
Qin_step = np.zeros(nt)
Qin_step[np.where(tm>=5)] += 15
Qin_step[np.where(tm>=10)] += 5
Qin_step[np.where(tm>=20)] -= 6
Qin_step[np.where(tm>=35)] -= 14
Qin = m.Param(value=Qin_step, name='Qin')
# Fixed variables, they are constant throughout the time horizon
Xtin = m.FV(value=0, name='Xtin')
Xvin = m.FV(value=0, name='Xvin')
Qe = m.FV(value=0, name='Qe')
Sin = m.FV(value=400, lb=0, ub=1500)
Pin = m.FV(value=0, name='Pin')
Fc = m.FV(value=40, name='Fc')
Fair = m.FV(value=60000, name='Fair')
Tin = m.FV(value=30, name='Tin')
Tcin = m.FV(value=15, name='Tcin')
Vl = m.Var(value=1000, name='Vl') # lb=-0.0, ub=0.75*V
Xt = m.Var(value=0.1, name='Xt') # lb=-0.0, ub=10
Xv = m.Var(value=0.1, name='Xv') # lb=-0.0, ub=10
S = m.Var(value=50, name='S') # lb=+0.0, ub=10000
P = m.Var(value=0, name='P')
Ol = m.Var(value=0.0065, name= 'Ol')
Og = m.Var(value=0.305, name='Og')
T = m.Var(value=30, name='T') # lb=20, ub=40
Tc = m.Var(value=20, name='Tc') # lb=0, ub=30
Sf_cum = m.Var(value=0, name='Sf_cum')
# Define algebraic equations
############################
# Specific growth rate of cell mass
mimax = m.Intermediate(((a1*(T-k1))*(1-m.exp(b1*(T-k2))))** 2)
Pmax = m.Intermediate(Pmaxb + PmaxT/(1-m.exp(-b2*(T-k3))))
m.Equation(mi == mimax * (S / (KSX+S)) * (Ol/(KOX + Ol)) \
* (1 - P/Pmax) * (1 / (1+m.exp(-(100-S)))))
# Specific production rate of EtOH
bP = m.Intermediate(c1*m.exp(-AP1/T) - c2*m.exp(-AP2/T))
qP = m.Intermediate(aP*mi + bP)
# Specific consumption rate of glucose
qS = m.Intermediate(mi/YXS + qP/YPS)
# Specific consumption rate of oxygen
qO = m.Intermediate(qOmax*Ol/YXO/(KOX+Ol))
# Specific biological deactivation rate of cell mass
Kd = m.Intermediate(Kdb + KdT/(1+m.exp(-b3*(T-k4))))
# Saturation concentration of oxygen in culture media
Ostar = m.Intermediate(z*Og*R*T/KH)
# Oxygen mass transfer coefficient
kla = m.Intermediate(kla0*1.2**(T-20))
# Bioreactor phases equation
Vg = m.Intermediate(V - Vl)
# Define differential equations
###############################
m.Equation(Vl.dt() == Qin - Qe)
m.Equation(Vl*Xt.dt() == Qin*(Xtin-Xt) + mi*Vl*Xv)
m.Equation(Vl*Xv.dt() == Qin*(Xvin-Xv) + Xv*Vl*(mi-Kd))
m.Equation(Vl*S.dt() == Qin*(Sin-S) - qS*Vl*Xv)
m.Equation(Vl*P.dt() == Qin*(Pin - P) + qP*Vl*Xv)
m.Equation(Vl*Ol.dt() == Qin*(Ostar-Ol) + Vl*kla*(Ostar-Ol) - qO*Vl*Xv)
m.Equation(Vg*Og.dt() == Fair*(Ogasin-Og) - Vl*kla*(Ostar-Ol) + Og*(Qin-Qe))
m.Equation(Vl*T.dt() == Qin*(Tin-T) - Tref*(Qin-Qe) \
+ Vl*qO*Xv*deltaH/MO/rho/Chbr - KT*AT*(T-Tc)/rho/Chbr)
m.Equation(Vcj*Tc.dt() == Fc*(Tcin - Tc) + KT*AT*(T-Tc)/rhoc/Chc)
m.Equation(Sf_cum.dt() == Qin*Sin)
# solve ODE
m.options.SOLVER= 1
m.options.IMODE = 7
m.options.NODES = 3
# m.open_folder()
m.solve(disp=False)
# Plot results
plt.figure(1)
plt.title('Total & Viable Cellular Biomass')
plt.plot(m.time, Xv.value, label='Xv')
plt.plot(m.time, Xt.value, label='Xt')
plt.legend()
plt.ylabel('Biomass concentration [g/L]')
plt.xlabel('Time [h]')
plt.grid()
plt.minorticks_on()
plt.ylim(0)
plt.xlim(m.time[0],m.time[-1])
plt.tight_layout()
plt.figure(2)
plt.title('Substrate (S) & Product (P) concentration')
plt.plot(m.time, S.value, label='S')
plt.plot(m.time, P.value, label='P')
plt.legend(); plt.grid()
plt.ylabel('Conc [g/L]')
plt.xlabel('Time [h]')
plt.minorticks_on()
plt.ylim(0)
plt.xlim(m.time[0],m.time[-1])
plt.tight_layout()
plt.figure(3)
plt.title('Bioreactor & Cooling jacket temperature')
plt.plot(m.time, T.value, label='T')
plt.plot(m.time, Tc.value, label='Tc')
plt.legend()
plt.ylabel('Temperature [oC]')
plt.xlabel('Time [h]')
plt.grid()
plt.minorticks_on()
plt.ylim(0)
plt.xlim(m.time[0],m.time[-1])
plt.tight_layout()
fig4, ax = plt.subplots()
ax.title.set_text('Dissolved & Gaseous Oxygen concentration')
lns1 = ax.plot(m.time, Ol.value, label='[Oliq]', color='c')
ax.set_xlabel('Time [h]')
ax.set_ylabel('Oliq [g/L]', color='c')
ax.minorticks_on()
ax2 = ax.twinx()
lns2 = ax2.plot(m.time, Og.value, label='[Ogas]', color='y')
ax2.set_ylabel('Ogas [g/L]', color='y')
ax2.minorticks_on()
lns = lns1 + lns2
labs = [l.get_label() for l in lns]
ax.legend(lns, labs, loc='best')
ax.grid()
fig4.tight_layout()
plt.figure(4)
plt.figure(5)
plt.title('Feeding Policy')
plt.plot(m.time, Qin.value, label='Qin')
plt.legend()
plt.ylabel('Qin [L/h]')
plt.xlabel('Time [h]')
plt.grid()
plt.minorticks_on()
plt.ylim(0)
plt.xlim(m.time[0],m.time[-1])
plt.tight_layout()
plt.figure(6)
plt.title('Check >=0 Constraints')
plt.plot(tm,bP.value,label='bP')
plt.plot(tm,mi.value,label='mi')
plt.legend(); plt.grid()
plt.show()
import numpy as np
from scipy.integrate import odeint
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
#Simulate fed-batch operation
# Specify the simulation time (hrs)
tspan = np.linspace(0,37,371); t=tspan
# Specify values of control variables [Qin0 Qin Xtin Xvin Qe Sin Fc Fair Tin Tcin]
u0 = [0.0,0.0,0.0,0.0,0.0,400,40,60000,30,15]
# Specify initial conditions [Xt Xv S P Oliq Ogas T Tc Vl Gloss MP]
x0 = [0.1,0.1,50,0,0.0065,0.305,30,20,1000,0.0,0.0]
ux0 = tuple(u0 + x0)
#Qin = 15*heaviside(t-5) + 5*heaviside(t-10)
# - 6*heaviside(t-20) - 14*heaviside(t-35)
Qin = np.zeros_like(tspan)
Qin[np.where(tspan>=5)] += 15
Qin[np.where(tspan>=10)] += 5
Qin[np.where(tspan>=20)] -= 6
Qin[np.where(tspan>=35)] -= 14
QinInterp = interp1d(tspan,Qin,bounds_error=False)
def ethanol(x,t,Qin0,Qin,Xtin,Xvin,Qe,Sin,Fc,Fair,
Tin,Tcin,Xt0,Xv0,S0,P0,Oliq0,Ogas0,T0,
Tc0,Vl0,Sf_cum0,Time0):
## #Initial Conditions
## Xt0 = u[10] # Initial total cellular biomass, [g L-1]
## Xv0 = u[11] # Initial viable cellular biomass, [g L-1]
## S0 = u[12] # Initial substrate/Glucose concentration, [g L-1]
## P0 = u[13] # Initial product/Ethanol concentration, [g L-1]
## Oliq0 = u[14] # Initial Dissolved oxygen concentration, [g L-1]
## Ogas0 = u[15] # Initial Gas phase oxygen (bubbles) in the fermentation broth, [g L-1]
## T0 = u[16] # Initial Temperature in the bioreactor, [oC]
## Tc0 = u[17] # Initial Temperature of the cooling agent in the jacket, [oC]
## Vl0 = u[18] # Initial Culture volume in the bioreactor, [L]
## Sf_cum0 = u[19] # Initial Cumulative substrate/glucose fed to the bioreactor, [g]
## Time0 = u[20] # Initial batch time, [h]
##
## #Control variables
## Qin0 = u[0] # Volumetric inflow rate, [l/h-1]
## Qin = u[1] # Volumetric inflow rate, [l/h-1]
## Xtin = u[2] # Total biomass concentration in the bioreactor feed, [g L-1]
## Xvin = u[3] # Viable biomass concentration in the bioreactor feed, [g L-1]
## Qe = u[4] # Volumetric outflow rate, [l/h-1]
## Sin = u[5] # Substrate/Glucose concentration in bioreactor feed, [g L-1]
## Fc = u[6] # Cooling agent inlet volumetric flowrate, [L h-1]
## Fair = u[7] # Airflow inlet volumetric flowrate, [L h-1]
## Tin = u[8] # Temperature of bioreactor feed, [oC]
## Tcin = u[9] # Temperature of cooling agent inlet, [oC]
# 1D Interpolation for Qin
Qin = QinInterp(t)
#Definition of model parameters
#Kinetic parameters
a1 = 0.05 # Ratkowsky parameter [oC-1 h-0.5]
aP = 4.50 # Growth-associated parameter for ethanol production, [-]
AP1 = 6.0 # Activation energy parameter for ethanol production, [oC]
AP2 = 20.3 # Activation energy parameter for ethanol production, [oC]
b1 = 0.035 # Parameter in the exponential expression of the maximum specific growth rate np.expression, [oC-1]
b2 = 0.15 # Parameter in the exponential expression of the growth inhibitory ethanol concentration np.expression, [oC-1]
b3 = 0.40 # Parameter in the exponential np.expression of the specific death rate expression,[oC-1]
c1 = 0.38 # Constant decoupling factor for ethanol production, [gP gX-1 h-1]
c2 = 0.29 # Constant decoupling factor for ethanol production, [gP gX-1 h-1]
k1 = 3.00 # Parameter in the maximum specific growth rate expression, [oC]
k2 = 55.0 # Parameter in the maximum specific growth rate expression, [oC]
k3 = 60.0 # Parameter in the growth-inhibitory ethanol concentration expression, [oC]
k4 = 50.0 # Temperature at the inflection point of the specific death rate sigmoid curve, [oC]
Pmaxb = 90 # Temperature-independent product inhibition constant, [g L-1]
PmaxT = 90 # Maximum value of product inhibition constant due to temperature, [g L-1]
Kdb = 0.025 # Basal specific cellular biomass death rate, [h-1]
KdT = 30.00 # Maximum value of specific cellular biomass death rate due to temperature, [h-1]
KSX = 5 # Glucose saturation constant for the specific growth rate, [g L-1]
KOX = 0.0005 # Oxygen saturation constant for the specific growth rate, [g L-1]
qOmax = 0.05 # Maximum specific oxygen consumption rate, [h-1]
#Metabolic parameters
YPS = 0.51 # Theoretical yield of ethanol on glucose, [gP gS-1]
YXO = 0.97 # Theoretical yield of biomass on oxygen, [gX gO-1]
YXS = 0.53 # Theoretical yield of biomass on glucose, [gX gS-1]
#Physicochemical and thermodynamic parameters
Chbr = 4.18 # Heat capacity of the mass of reaction, [J g-1 oC-1]
Chc = 4.18 # Heat capacity of the cooling agent, [J g-1 oC-1]
DeltaH = 518000 # Heat of reaction of fermentation, [J mol-1 O2]
Tref = 20 # Reference temperature, [oC]
KH = 200 # Henry's constant for oxygen in the fermentation broth, [atm L mol-1]
z = 0.792 # Oxygen compressibility factor, [-]
R = 0.082 # Ideas gas constant, [L atm mol-1 oC-1]
kla0 = 100 # Temperature-independent volumetric oxygen transfer coefficient, [h-1]
KT = 360000 # Heat transfer coefficient, [J h-1 m-2 ??C-1]
rho = 1080 # Density of the fermentation broth, [g L-1]
rhoc = 1000 # Density of the cooling agent, [g L-1]
MO = 32.0 # Molecular weight of oxygen (O2), [g mol-1]
#Bioreactor design data
AT = 1.0 # Bioreactor heat transfer area, [m2]
V = 1800 # Bioreactor working volume, [L]
Vcj = 50 # Cooling jacket volume, [L]
Ogasin = 0.305 # Oxygen concentration in airflow inlet, [g L-1]
#Definition of model variables
#State variables
Xt = x[0] # Total cellular biomass, [g L-1]
Xv = x[1] # Viable cellular biomass, [g L-1]
S = x[2] # Substrate/Glucose concentration, [g L-1]
P = x[3] # Product/Ethanol concentration, [g L-1]
Oliq = x[4] # Dissolved oxygen concentration, [g L-1]
Ogas = x[5] # Gas phase oxygen (bubbles) in the fermentation broth, [g L-1]
T = x[6] # Temperature in the bioreactor, [oC]
Tc = x[7] # Temperature of the cooling agent in the jacket, [oC]
Vl = x[8] # Culture volume in the bioreactor, [L]
Sf_cum = x[9] # Cumulative amount of substrate/glucose fed to the bioreactor, [g]
Time = x[10] # Batch time, [h]
# Definition of model equations
# Kinetic rates
# -----------------------------
# Specific growth rate, [h-1]
mmax = ((a1*(T-k1))*(1-np.exp(b1*(T-k2))))**2
Pmax = Pmaxb + PmaxT/(1-np.exp(-b2*(T-k3)))
m1 = mmax * S/(KSX + S) * Oliq/(KOX + Oliq) * (1 - P/Pmax) * 1/(1+np.exp(-(100-S)/1)) # Specific growth rate, [h-1]
if m1 >= 0:
m = m1
else:
m=0.0
# Non-growth-associated ethanol specific production rate, [h-1]
if S > 0:
bP = c1 * np.exp(-AP1/T) - c2 * np.exp(-AP2/T) # Non-growth-associated ethanol specific production rate, [h-1]
else:
bP = 0.0
qP = aP*m + bP
# Ethanol consumption specific rate
qS = m/YXS + qP/YPS
# Oxygen consumption specific rate
qO = qOmax*Oliq/YXO/(KOX + Oliq)
# Specific biological deactivation rate of cell mass
Kd = Kdb + KdT/(1+np.exp(-b3*(T-k4)))
# Saturation concentration of oxygen in culture media
Osat = z*Ogas*R*T/KH
# Oxygen mass transfer coefficient
kla = kla0*1.2**(T-20)
# Volume of the gas phase in the bioreactor
Vg = V - Vl
#Material balances
#-----------------
# Volume of liquid culture
dVl = Qin - Qe
# Total cell mass
dXt = m*Xv + Qin/Vl*(Xtin-Xt)
# Total mass of biologically active cells
dXv = (m-Kd)*Xv + Qin/Vl*(Xvin-Xv)
# Glucose concentration
dS = Qin/Vl*(Sin-S) - qS*Xv
# Ethanol concentration
dP = Qin/Vl*(-P) + qP*Xv
# Disolved oxygen
dOliq = Qin/Vl*(Osat - Oliq) + kla*(Osat-Oliq) - qO*Xv
# Oxygen gas phase
dOgas = Fair/Vg*(Ogasin-Ogas) - Vl*kla/Vg*(Osat - Oliq) + Ogas*(Qin-Qe)/Vg
# Energy balances
#---------------
# Bioreactor temprature
dT = Qin/Vl*(Tin-T) - Tref/Vl*(Qin-Qe) + qO*Xv*DeltaH/MO/rho/Chbr - KT*AT*(T-Tc)/Vl/rho/Chbr
# Cooling agent temperature
dTc = Fc/Vcj*(Tcin-Tc) + KT*AT*(T-Tc)/Vcj/rhoc/Chc
# Yields & Productivity
#---------------------
# Cumulative amount of glucose fed to the bioreactor
dSf_cum = Sin*Qin
dTime = 1
# Definition of state derivatives vector
# State derivatives
dxdt = [dXt,dXv,dS,dP,dOliq,dOgas,dT,dTc,dVl,dSf_cum,dTime]
# [dxdt,mmax,Pmax,bP,m,Kd,Qin]
return dxdt
# test function
print(ethanol(x0,0.0,*ux0))
# Simulate the bioreactor operation until the selected time tf
x = odeint(ethanol,x0,tspan,args=ux0)
#plots Results
#Total and Viable Cellular Biomass
plt.figure(figsize=(7,4))
plt.plot(tspan,x[:,0])
plt.plot(tspan,x[:,1])
plt.title('Total & Viable Cellular Biomass')
plt.ylabel('Biomass concentration [g/L]')
plt.xlabel('t [h]')
plt.grid()
plt.legend(['Xt','Xv'])
plt.tight_layout()
plt.figure(figsize=(7,4))
plt.title('Substrate (S) & Product (P) concentration')
plt.plot(tspan,x[:,2], label='S')
plt.plot(tspan,x[:,3], label='P')
plt.legend(); plt.grid()
plt.ylabel('Conc [g/L]')
plt.xlabel('Time [h]')
plt.minorticks_on()
plt.ylim(0)
plt.xlim(t[0],t[-1])
plt.tight_layout()
plt.figure()
plt.title('Bioreactor & Cooling jacket temperature')
plt.plot(tspan,x[:,6], label='T')
plt.plot(tspan,x[:,7], label='Tc')
plt.legend()
plt.ylabel('Temperature [oC]')
plt.xlabel('Time [h]')
plt.grid()
plt.minorticks_on()
plt.ylim(0)
plt.xlim(t[0],t[-1])
plt.tight_layout()
fig4, ax = plt.subplots()
ax.title.set_text('Dissolved & Gaseous Oxygen concentration')
lns1 = ax.plot(t,x[:,4], label='[Oliq]', color='c')
ax.set_xlabel('Time [h]')
ax.set_ylabel('Oliq [g/L]', color='c')
ax.minorticks_on()
ax2 = ax.twinx()
lns2 = ax2.plot(t,x[:,5], label='[Ogas]', color='y')
ax2.set_ylabel('Ogas [g/L]', color='y')
ax2.minorticks_on()
lns = lns1 + lns2
labs = [l.get_label() for l in lns]
ax.legend(lns, labs, loc='best')
ax.grid()
fig4.tight_layout()
plt.figure()
plt.title('Feeding Policy')
plt.plot(tspan, Qin, label='Qin')
plt.legend()
plt.ylabel('Qin [L/h]')
plt.xlabel('Time [h]')
plt.grid()
plt.minorticks_on()
plt.ylim(0)
plt.xlim(tspan[0],tspan[-1])
plt.tight_layout()
plt.show()
from scipy.integrate import odeint
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
#Simulate fed-batch operation
# Specify the simulation time (hrs)
tspan = np.linspace(0,37,371); t=tspan
# Specify values of control variables [Qin0 Qin Xtin Xvin Qe Sin Fc Fair Tin Tcin]
u0 = [0.0,0.0,0.0,0.0,0.0,400,40,60000,30,15]
# Specify initial conditions [Xt Xv S P Oliq Ogas T Tc Vl Gloss MP]
x0 = [0.1,0.1,50,0,0.0065,0.305,30,20,1000,0.0,0.0]
ux0 = tuple(u0 + x0)
#Qin = 15*heaviside(t-5) + 5*heaviside(t-10)
# - 6*heaviside(t-20) - 14*heaviside(t-35)
Qin = np.zeros_like(tspan)
Qin[np.where(tspan>=5)] += 15
Qin[np.where(tspan>=10)] += 5
Qin[np.where(tspan>=20)] -= 6
Qin[np.where(tspan>=35)] -= 14
QinInterp = interp1d(tspan,Qin,bounds_error=False)
def ethanol(x,t,Qin0,Qin,Xtin,Xvin,Qe,Sin,Fc,Fair,
Tin,Tcin,Xt0,Xv0,S0,P0,Oliq0,Ogas0,T0,
Tc0,Vl0,Sf_cum0,Time0):
## #Initial Conditions
## Xt0 = u[10] # Initial total cellular biomass, [g L-1]
## Xv0 = u[11] # Initial viable cellular biomass, [g L-1]
## S0 = u[12] # Initial substrate/Glucose concentration, [g L-1]
## P0 = u[13] # Initial product/Ethanol concentration, [g L-1]
## Oliq0 = u[14] # Initial Dissolved oxygen concentration, [g L-1]
## Ogas0 = u[15] # Initial Gas phase oxygen (bubbles) in the fermentation broth, [g L-1]
## T0 = u[16] # Initial Temperature in the bioreactor, [oC]
## Tc0 = u[17] # Initial Temperature of the cooling agent in the jacket, [oC]
## Vl0 = u[18] # Initial Culture volume in the bioreactor, [L]
## Sf_cum0 = u[19] # Initial Cumulative substrate/glucose fed to the bioreactor, [g]
## Time0 = u[20] # Initial batch time, [h]
##
## #Control variables
## Qin0 = u[0] # Volumetric inflow rate, [l/h-1]
## Qin = u[1] # Volumetric inflow rate, [l/h-1]
## Xtin = u[2] # Total biomass concentration in the bioreactor feed, [g L-1]
## Xvin = u[3] # Viable biomass concentration in the bioreactor feed, [g L-1]
## Qe = u[4] # Volumetric outflow rate, [l/h-1]
## Sin = u[5] # Substrate/Glucose concentration in bioreactor feed, [g L-1]
## Fc = u[6] # Cooling agent inlet volumetric flowrate, [L h-1]
## Fair = u[7] # Airflow inlet volumetric flowrate, [L h-1]
## Tin = u[8] # Temperature of bioreactor feed, [oC]
## Tcin = u[9] # Temperature of cooling agent inlet, [oC]
# 1D Interpolation for Qin
Qin = QinInterp(t)
#Definition of model parameters
#Kinetic parameters
a1 = 0.05 # Ratkowsky parameter [oC-1 h-0.5]
aP = 4.50 # Growth-associated parameter for ethanol production, [-]
AP1 = 6.0 # Activation energy parameter for ethanol production, [oC]
AP2 = 20.3 # Activation energy parameter for ethanol production, [oC]
b1 = 0.035 # Parameter in the exponential expression of the maximum specific growth rate np.expression, [oC-1]
b2 = 0.15 # Parameter in the exponential expression of the growth inhibitory ethanol concentration np.expression, [oC-1]
b3 = 0.40 # Parameter in the exponential np.expression of the specific death rate expression,[oC-1]
c1 = 0.38 # Constant decoupling factor for ethanol production, [gP gX-1 h-1]
c2 = 0.29 # Constant decoupling factor for ethanol production, [gP gX-1 h-1]
k1 = 3.00 # Parameter in the maximum specific growth rate expression, [oC]
k2 = 55.0 # Parameter in the maximum specific growth rate expression, [oC]
k3 = 60.0 # Parameter in the growth-inhibitory ethanol concentration expression, [oC]
k4 = 50.0 # Temperature at the inflection point of the specific death rate sigmoid curve, [oC]
Pmaxb = 90 # Temperature-independent product inhibition constant, [g L-1]
PmaxT = 90 # Maximum value of product inhibition constant due to temperature, [g L-1]
Kdb = 0.025 # Basal specific cellular biomass death rate, [h-1]
KdT = 30.00 # Maximum value of specific cellular biomass death rate due to temperature, [h-1]
KSX = 5 # Glucose saturation constant for the specific growth rate, [g L-1]
KOX = 0.0005 # Oxygen saturation constant for the specific growth rate, [g L-1]
qOmax = 0.05 # Maximum specific oxygen consumption rate, [h-1]
#Metabolic parameters
YPS = 0.51 # Theoretical yield of ethanol on glucose, [gP gS-1]
YXO = 0.97 # Theoretical yield of biomass on oxygen, [gX gO-1]
YXS = 0.53 # Theoretical yield of biomass on glucose, [gX gS-1]
#Physicochemical and thermodynamic parameters
Chbr = 4.18 # Heat capacity of the mass of reaction, [J g-1 oC-1]
Chc = 4.18 # Heat capacity of the cooling agent, [J g-1 oC-1]
DeltaH = 518000 # Heat of reaction of fermentation, [J mol-1 O2]
Tref = 20 # Reference temperature, [oC]
KH = 200 # Henry's constant for oxygen in the fermentation broth, [atm L mol-1]
z = 0.792 # Oxygen compressibility factor, [-]
R = 0.082 # Ideas gas constant, [L atm mol-1 oC-1]
kla0 = 100 # Temperature-independent volumetric oxygen transfer coefficient, [h-1]
KT = 360000 # Heat transfer coefficient, [J h-1 m-2 ??C-1]
rho = 1080 # Density of the fermentation broth, [g L-1]
rhoc = 1000 # Density of the cooling agent, [g L-1]
MO = 32.0 # Molecular weight of oxygen (O2), [g mol-1]
#Bioreactor design data
AT = 1.0 # Bioreactor heat transfer area, [m2]
V = 1800 # Bioreactor working volume, [L]
Vcj = 50 # Cooling jacket volume, [L]
Ogasin = 0.305 # Oxygen concentration in airflow inlet, [g L-1]
#Definition of model variables
#State variables
Xt = x[0] # Total cellular biomass, [g L-1]
Xv = x[1] # Viable cellular biomass, [g L-1]
S = x[2] # Substrate/Glucose concentration, [g L-1]
P = x[3] # Product/Ethanol concentration, [g L-1]
Oliq = x[4] # Dissolved oxygen concentration, [g L-1]
Ogas = x[5] # Gas phase oxygen (bubbles) in the fermentation broth, [g L-1]
T = x[6] # Temperature in the bioreactor, [oC]
Tc = x[7] # Temperature of the cooling agent in the jacket, [oC]
Vl = x[8] # Culture volume in the bioreactor, [L]
Sf_cum = x[9] # Cumulative amount of substrate/glucose fed to the bioreactor, [g]
Time = x[10] # Batch time, [h]
# Definition of model equations
# Kinetic rates
# -----------------------------
# Specific growth rate, [h-1]
mmax = ((a1*(T-k1))*(1-np.exp(b1*(T-k2))))**2
Pmax = Pmaxb + PmaxT/(1-np.exp(-b2*(T-k3)))
m1 = mmax * S/(KSX + S) * Oliq/(KOX + Oliq) * (1 - P/Pmax) * 1/(1+np.exp(-(100-S)/1)) # Specific growth rate, [h-1]
if m1 >= 0:
m = m1
else:
m=0.0
# Non-growth-associated ethanol specific production rate, [h-1]
if S > 0:
bP = c1 * np.exp(-AP1/T) - c2 * np.exp(-AP2/T) # Non-growth-associated ethanol specific production rate, [h-1]
else:
bP = 0.0
qP = aP*m + bP
# Ethanol consumption specific rate
qS = m/YXS + qP/YPS
# Oxygen consumption specific rate
qO = qOmax*Oliq/YXO/(KOX + Oliq)
# Specific biological deactivation rate of cell mass
Kd = Kdb + KdT/(1+np.exp(-b3*(T-k4)))
# Saturation concentration of oxygen in culture media
Osat = z*Ogas*R*T/KH
# Oxygen mass transfer coefficient
kla = kla0*1.2**(T-20)
# Volume of the gas phase in the bioreactor
Vg = V - Vl
#Material balances
#-----------------
# Volume of liquid culture
dVl = Qin - Qe
# Total cell mass
dXt = m*Xv + Qin/Vl*(Xtin-Xt)
# Total mass of biologically active cells
dXv = (m-Kd)*Xv + Qin/Vl*(Xvin-Xv)
# Glucose concentration
dS = Qin/Vl*(Sin-S) - qS*Xv
# Ethanol concentration
dP = Qin/Vl*(-P) + qP*Xv
# Disolved oxygen
dOliq = Qin/Vl*(Osat - Oliq) + kla*(Osat-Oliq) - qO*Xv
# Oxygen gas phase
dOgas = Fair/Vg*(Ogasin-Ogas) - Vl*kla/Vg*(Osat - Oliq) + Ogas*(Qin-Qe)/Vg
# Energy balances
#---------------
# Bioreactor temprature
dT = Qin/Vl*(Tin-T) - Tref/Vl*(Qin-Qe) + qO*Xv*DeltaH/MO/rho/Chbr - KT*AT*(T-Tc)/Vl/rho/Chbr
# Cooling agent temperature
dTc = Fc/Vcj*(Tcin-Tc) + KT*AT*(T-Tc)/Vcj/rhoc/Chc
# Yields & Productivity
#---------------------
# Cumulative amount of glucose fed to the bioreactor
dSf_cum = Sin*Qin
dTime = 1
# Definition of state derivatives vector
# State derivatives
dxdt = [dXt,dXv,dS,dP,dOliq,dOgas,dT,dTc,dVl,dSf_cum,dTime]
# [dxdt,mmax,Pmax,bP,m,Kd,Qin]
return dxdt
# test function
print(ethanol(x0,0.0,*ux0))
# Simulate the bioreactor operation until the selected time tf
x = odeint(ethanol,x0,tspan,args=ux0)
#plots Results
#Total and Viable Cellular Biomass
plt.figure(figsize=(7,4))
plt.plot(tspan,x[:,0])
plt.plot(tspan,x[:,1])
plt.title('Total & Viable Cellular Biomass')
plt.ylabel('Biomass concentration [g/L]')
plt.xlabel('t [h]')
plt.grid()
plt.legend(['Xt','Xv'])
plt.tight_layout()
plt.figure(figsize=(7,4))
plt.title('Substrate (S) & Product (P) concentration')
plt.plot(tspan,x[:,2], label='S')
plt.plot(tspan,x[:,3], label='P')
plt.legend(); plt.grid()
plt.ylabel('Conc [g/L]')
plt.xlabel('Time [h]')
plt.minorticks_on()
plt.ylim(0)
plt.xlim(t[0],t[-1])
plt.tight_layout()
plt.figure()
plt.title('Bioreactor & Cooling jacket temperature')
plt.plot(tspan,x[:,6], label='T')
plt.plot(tspan,x[:,7], label='Tc')
plt.legend()
plt.ylabel('Temperature [oC]')
plt.xlabel('Time [h]')
plt.grid()
plt.minorticks_on()
plt.ylim(0)
plt.xlim(t[0],t[-1])
plt.tight_layout()
fig4, ax = plt.subplots()
ax.title.set_text('Dissolved & Gaseous Oxygen concentration')
lns1 = ax.plot(t,x[:,4], label='[Oliq]', color='c')
ax.set_xlabel('Time [h]')
ax.set_ylabel('Oliq [g/L]', color='c')
ax.minorticks_on()
ax2 = ax.twinx()
lns2 = ax2.plot(t,x[:,5], label='[Ogas]', color='y')
ax2.set_ylabel('Ogas [g/L]', color='y')
ax2.minorticks_on()
lns = lns1 + lns2
labs = [l.get_label() for l in lns]
ax.legend(lns, labs, loc='best')
ax.grid()
fig4.tight_layout()
plt.figure()
plt.title('Feeding Policy')
plt.plot(tspan, Qin, label='Qin')
plt.legend()
plt.ylabel('Qin [L/h]')
plt.xlabel('Time [h]')
plt.grid()
plt.minorticks_on()
plt.ylim(0)
plt.xlim(tspan[0],tspan[-1])
plt.tight_layout()
plt.show()
Save as reactor.m
clear; clc; close all;
%Simulate fed-batch operation for xx hours
% Specify the simulation time
tspan = [0 37];
% Specify values of control variables [Qin0; Xtin; Xvin; Qe; Sin; Fc; Fair; Tin; Tcin]
u0 = [0.0; 0.0; 0.0; 0.0; 400; 40; 60000; 30; 15];
% Specify initial conditions [Xt; Xv; S; P; Oliq; Ogas; T; Tc; Vl; Gloss; MP]
x0 = [0.1; 0.1; 50; 0; 0.0065; 0.305; 30; 20; 1000; 0.0; 0.0];
ux0 = [u0;x0];
% Simulate the bioreactor operation until the selected time tf
[t,x] = ode15s(@EthanolModel,tspan,x0,[],ux0);
N=length(t);
for i=1:size(t,1)
[~,mmax(i,:),Pmax(i,:),bP(i,:),m(i,:),Kd(i,:),Qin(i,:)] = EthanolModel(t(i,:),x(i,:),ux0);
end
% Loss of substrate at the end of the bioreactor operation-Cumulative ethanol mass produced
Gloss = x(N,3)*x(N,9);
%plots Results
%Total and Viable Cellular Biomass
figure(1)
f=plot(t,x(:,1), t,x(:,2));
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid
title('Total & Viable Cellular Biomass')
ylabel('Biomass concentration [g/L]','FontSize',12,'FontWeight', 'bold')
xlabel('t [h]','FontSize', 12, 'FontWeight','bold')
legend('Xt','Xv')
%Substrate(S) and Production(P) concentration
figure(2)
f=plot(t,x(:,3),'y', t,x(:,4),'g');
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid
title('Substrate (S) & Product (P) concetration')
ylabel('Concentration [g/L]','FontSize',12,'FontWeight', 'bold')
xlabel('t [h]','FontSize', 12, 'FontWeight','bold')
legend('S','P')
% Dissolved and gaseous oxygen concentration
figure(3)
title('Dissolved and Gaseous oxygen concentration')
hold on
yyaxis left
f=plot(t,x(:,5),'color',[0.3010 0.7450 0.9330]);
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid on
ylabel('Oliq [g/L]','FontSize',12,'FontWeight', 'bold', 'color',[0.3010 0.7450 0.9330])
yyaxis right
f=plot(t,x(:,6),'color',[250/255 128/255 114/255]);
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid on
ylabel('Ogas [g/L]','FontSize',12,'FontWeight', 'bold','color',[250/255 128/255 114/255])
xlabel('t [h]','FontSize', 12, 'FontWeight','bold')
legend('Oliq','Ogas')
hold off
% Bioreactor & Cooling jacket temperature
figure(4)
f=plot(t,x(:,7), t,x(:,8));
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid
title('Bioreactor & Cooling jacket temperature')
ylabel('Temperature [oC]','FontSize',12,'FontWeight', 'bold')
xlabel('t [h]','FontSize', 12, 'FontWeight','bold')
legend('T','Tc')
% Specific growth and death rate
figure(5)
f=plot(t,mmax, t,m, t,Kd);
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid
title('Specific growth and death rate')
ylabel('Specific growth and death rate [1/h]','FontSize',12,'FontWeight', 'bold')
xlabel('t [h]','FontSize', 12, 'FontWeight','bold')
legend('mimax','mi','Kd')
% Production rate
figure(6)
title('Production rate')
hold on
yyaxis left
f=plot(t,bP);
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid on
ylabel('bP [1/h]','FontSize',12,'FontWeight', 'bold')
yyaxis right
f=plot(t,Pmax);
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid on
ylabel('Pmax','FontSize',12,'FontWeight', 'bold')
xlabel('t [h]','FontSize', 12, 'FontWeight','bold')
legend('bP','Pmax')
hold off
% Culture Volume & Total substrate fed to the BR
figure(7)
title('Culture Volume & Total substrate fed to the BR')
hold on
yyaxis left
f=plot(t,x(:,9));
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid on
ylabel('Vl [L]','FontSize',12,'FontWeight', 'bold')
yyaxis right
f=plot(t,x(:,10));
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid on
ylabel('Sf_cum [g]','FontSize',12,'FontWeight', 'bold')
xlabel('t [h]','FontSize', 12, 'FontWeight','bold')
legend('Vl','Sf_cum')
hold off
% Feeding Policy
figure(8)
title('Feeding policy')
f=plot(t,Qin);
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid on
ylabel('Qin [L/h]','FontSize',12,'FontWeight', 'bold')
xlabel('t [h]','FontSize', 12, 'FontWeight','bold')
legend('Qin')
%Simulate fed-batch operation for xx hours
% Specify the simulation time
tspan = [0 37];
% Specify values of control variables [Qin0; Xtin; Xvin; Qe; Sin; Fc; Fair; Tin; Tcin]
u0 = [0.0; 0.0; 0.0; 0.0; 400; 40; 60000; 30; 15];
% Specify initial conditions [Xt; Xv; S; P; Oliq; Ogas; T; Tc; Vl; Gloss; MP]
x0 = [0.1; 0.1; 50; 0; 0.0065; 0.305; 30; 20; 1000; 0.0; 0.0];
ux0 = [u0;x0];
% Simulate the bioreactor operation until the selected time tf
[t,x] = ode15s(@EthanolModel,tspan,x0,[],ux0);
N=length(t);
for i=1:size(t,1)
[~,mmax(i,:),Pmax(i,:),bP(i,:),m(i,:),Kd(i,:),Qin(i,:)] = EthanolModel(t(i,:),x(i,:),ux0);
end
% Loss of substrate at the end of the bioreactor operation-Cumulative ethanol mass produced
Gloss = x(N,3)*x(N,9);
%plots Results
%Total and Viable Cellular Biomass
figure(1)
f=plot(t,x(:,1), t,x(:,2));
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid
title('Total & Viable Cellular Biomass')
ylabel('Biomass concentration [g/L]','FontSize',12,'FontWeight', 'bold')
xlabel('t [h]','FontSize', 12, 'FontWeight','bold')
legend('Xt','Xv')
%Substrate(S) and Production(P) concentration
figure(2)
f=plot(t,x(:,3),'y', t,x(:,4),'g');
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid
title('Substrate (S) & Product (P) concetration')
ylabel('Concentration [g/L]','FontSize',12,'FontWeight', 'bold')
xlabel('t [h]','FontSize', 12, 'FontWeight','bold')
legend('S','P')
% Dissolved and gaseous oxygen concentration
figure(3)
title('Dissolved and Gaseous oxygen concentration')
hold on
yyaxis left
f=plot(t,x(:,5),'color',[0.3010 0.7450 0.9330]);
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid on
ylabel('Oliq [g/L]','FontSize',12,'FontWeight', 'bold', 'color',[0.3010 0.7450 0.9330])
yyaxis right
f=plot(t,x(:,6),'color',[250/255 128/255 114/255]);
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid on
ylabel('Ogas [g/L]','FontSize',12,'FontWeight', 'bold','color',[250/255 128/255 114/255])
xlabel('t [h]','FontSize', 12, 'FontWeight','bold')
legend('Oliq','Ogas')
hold off
% Bioreactor & Cooling jacket temperature
figure(4)
f=plot(t,x(:,7), t,x(:,8));
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid
title('Bioreactor & Cooling jacket temperature')
ylabel('Temperature [oC]','FontSize',12,'FontWeight', 'bold')
xlabel('t [h]','FontSize', 12, 'FontWeight','bold')
legend('T','Tc')
% Specific growth and death rate
figure(5)
f=plot(t,mmax, t,m, t,Kd);
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid
title('Specific growth and death rate')
ylabel('Specific growth and death rate [1/h]','FontSize',12,'FontWeight', 'bold')
xlabel('t [h]','FontSize', 12, 'FontWeight','bold')
legend('mimax','mi','Kd')
% Production rate
figure(6)
title('Production rate')
hold on
yyaxis left
f=plot(t,bP);
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid on
ylabel('bP [1/h]','FontSize',12,'FontWeight', 'bold')
yyaxis right
f=plot(t,Pmax);
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid on
ylabel('Pmax','FontSize',12,'FontWeight', 'bold')
xlabel('t [h]','FontSize', 12, 'FontWeight','bold')
legend('bP','Pmax')
hold off
% Culture Volume & Total substrate fed to the BR
figure(7)
title('Culture Volume & Total substrate fed to the BR')
hold on
yyaxis left
f=plot(t,x(:,9));
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid on
ylabel('Vl [L]','FontSize',12,'FontWeight', 'bold')
yyaxis right
f=plot(t,x(:,10));
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid on
ylabel('Sf_cum [g]','FontSize',12,'FontWeight', 'bold')
xlabel('t [h]','FontSize', 12, 'FontWeight','bold')
legend('Vl','Sf_cum')
hold off
% Feeding Policy
figure(8)
title('Feeding policy')
f=plot(t,Qin);
set(gca,'linewidth',1,'fontsize',12)
set(f,'linewidth',3)
grid on
ylabel('Qin [L/h]','FontSize',12,'FontWeight', 'bold')
xlabel('t [h]','FontSize', 12, 'FontWeight','bold')
legend('Qin')
Save as EthanolModel.m
function [dxdt,mmax,Pmax,bP,m,Kd,Qin] = EthanolModel(t,x,u)
%Definition of model parameters
%Kinetic parameters
a1 = 0.05; % Ratkowsky parameter [oC-1 h-0.5]
aP = 4.50; % Growth-associated parameter for ethanol production, [-]
AP1 = 6.0; % Activation energy parameter for ethanol production, [oC]
AP2 = 20.3; % Activation energy parameter for ethanol production, [oC]
b1 = 0.035; % Parameter in the exponential expression of the maximum specific growth rate expression, [oC-1]
b2 = 0.15; % Parameter in the exponential expression of the growth inhibitory ethanol concentration expression, [oC-1]
b3 = 0.40; % Parameter in the exponential expression of the specific death rate expression,[oC-1]
c1 = 0.38; % Constant decoupling factor for ethanol production, [gP gX-1 h-1]
c2 = 0.29; % Constant decoupling factor for ethanol production, [gP gX-1 h-1]
k1 = 3.00; % Parameter in the maximum specific growth rate expression, [oC]
k2 = 55.0; % Parameter in the maximum specific growth rate expression, [oC]
k3 = 60.0; % Parameter in the growth-inhibitory ethanol concentration expression, [oC]
k4 = 50.0; % Temperature at the inflection point of the specific death rate sigmoid curve, [oC]
Pmaxb = 90; % Temperature-independent product inhibition constant, [g L-1]
PmaxT = 90; % Maximum value of product inhibition constant due to temperature, [g L-1]
Kdb = 0.025; % Basal specific cellular biomass death rate, [h-1]
KdT = 30.00; % Maximum value of specific cellular biomass death rate due to temperature, [h-1]
KSX = 5; % Glucose saturation constant for the specific growth rate, [g L-1]
KOX = 0.0005; % Oxygen saturation constant for the specific growth rate, [g L-1]
qOmax = 0.05; % Maximum specific oxygen consumption rate, [h-1]
%Metabolic parameters
YPS = 0.51; % Theoretical yield of ethanol on glucose, [gP gS-1]
YXO = 0.97; % Theoretical yield of biomass on oxygen, [gX gO-1]
YXS = 0.53; % Theoretical yield of biomass on glucose, [gX gS-1]
%Physicochemical and thermodynamic parameters
Chbr = 4.18; % Heat capacity of the mass of reaction, [J g-1 oC-1]
Chc = 4.18; % Heat capacity of the cooling agent, [J g-1 oC-1]
DeltaH = 518000; % Heat of reaction of fermentation, [J mol-1 O2]
Tref = 20; % Reference temperature, [oC]
KH = 200; % Henry's constant for oxygen in the fermentation broth, [atm L mol-1]
z = 0.792; % Oxygen compressibility factor, [-]
R = 0.082; % Ideas gas constant, [L atm mol-1 oC-1]
kla0 = 100 ; % Temperature-independent volumetric oxygen transfer coefficient, [h-1]
KT = 360000; % Heat transfer coefficient, [J h-1 m-2 ??C-1]
rho = 1080 ; % Density of the fermentation broth, [g L-1]
rhoc = 1000 ; % Density of the cooling agent, [g L-1]
MO = 32.0; % Molecular weight of oxygen (O2), [g mol-1]
%Bioreactor design data
AT = 1.0; % Bioreactor heat transfer area, [m2]
V = 1800; % Bioreactor working volume, [L]
Vcj = 50; % Cooling jacket volume, [L]
Ogasin = 0.305; % Oxygen concentration in airflow inlet, [g L-1]
%Definition of model variables
%State variables
Xt = x(1); % Total cellular biomass, [g L-1]
Xv = x(2); % Viable cellular biomass, [g L-1]
S = x(3); % Substrate/Glucose concentration, [g L-1]
P = x(4); % Product/Ethanol concentration, [g L-1]
Oliq = x(5); % Dissolved oxygen concentration, [g L-1]
Ogas = x(6); % Gas phase oxygen (bubbles) in the fermentation broth, [g L-1]
T = x(7); % Temperature in the bioreactor, [oC]
Tc = x(8); % Temperature of the cooling agent in the jacket, [oC]
Vl = x(9); % Culture volume in the bioreactor, [L]
Sf_cum = x(10); % Cumulative amount of substrate/glucose fed to the bioreactor, [g]
Time = x(11); % Batch time, [h]
%Initial Conditions
Xt0 = u(10); % Initial total cellular biomass, [g L-1]
Xv0 = u(11); % Initial viable cellular biomass, [g L-1]
S0 = u(12); % Initial substrate/Glucose concentration, [g L-1]
P0 = u(13); % Initial product/Ethanol concentration, [g L-1]
Oliq0 = u(14); % Initial Dissolved oxygen concentration, [g L-1]
Ogas0 = u(15); % Initial Gas phase oxygen (bubbles) in the fermentation broth, [g L-1]
T0 = u(16); % Initial Temperature in the bioreactor, [oC]
Tc0 = u(17); % Initial Temperature of the cooling agent in the jacket, [oC]
Vl0 = u(18); % Initial Culture volume in the bioreactor, [L]
Sf_cum0 = u(19); % Initial Cumulative substrate/glucose fed to the bioreactor, [g]
Time0 = u(20); % Initial batch time, [h]
%Control variables
Qin0 = u(1); % Volumetric inflow rate, [l/h-1]
Qin= Qin0 + 15*heaviside(t-5) + 5*heaviside(t-10) - 6*heaviside(t-20) - 14*heaviside(t-35);
Xtin = u(2); % Total biomass concentration in the bioreactor feed, [g L-1]
Xvin = u(3); % Viable biomass concentration in the bioreactor feed, [g L-1]
Qe = u(4); % Volumetric outflow rate, [l/h-1]
Sin = u(5); % Substrate/Glucose concentration in bioreactor feed, [g L-1]
Fc = u(6); % Cooling agent inlet volumetric flowrate, [L h-1]
Fair = u(7); % Airflow inlet volumetric flowrate, [L h-1]
Tin = u(8); % Temperature of bioreactor feed, [oC]
Tcin = u(9); % Temperature of cooling agent inlet, [oC]
% Definition of model equations
% Kinetic rates
% -----------------------------
% Specific growth rate, [h-1]
mmax = ((a1*(T-k1))*(1-exp(b1*(T-k2))))^2;
Pmax = Pmaxb + PmaxT/(1-exp(-b2*(T-k3)));
m1 = mmax * S/(KSX + S) * Oliq/(KOX + Oliq) * (1 - P/Pmax) * 1/(1+exp(-(100-S)/1)); % Specific growth rate, [h-1]
if m1 >= 0
m = m1;
else
m=0.0;
end
% Non-growth-associated ethanol specific production rate, [h-1]
if S > 0
bP = c1 * exp(-AP1/T) - c2 * exp(-AP2/T) ; % Non-growth-associated ethanol specific production rate, [h-1]
else
bP = 0.0;
end
qP = aP*m + bP;
% Ethanol consumption specific rate
qS = m/YXS + qP/YPS;
% Oxygen consumption specific rate
qO = qOmax*Oliq/YXO/(KOX + Oliq);
% Specific biological deactivation rate of cell mass
Kd = Kdb + KdT/(1+exp(-b3*(T-k4)));
% Saturation concentration of oxygen in culture media
Osat = z*Ogas*R*T/KH;
% Oxygen mass transfer coefficient
kla = kla0*1.2^(T-20);
% Volume of the gas phase in the bioreactor
Vg = V - Vl;
%Material balances
%-----------------
% Volume of liquid culture
dVl = Qin - Qe;
% Total cell mass
dXt = m*Xv + Qin/Vl*(Xtin-Xt);
% Total mass of biologically active cells
dXv = (m-Kd)*Xv + Qin/Vl*(Xvin-Xv);
% Glucose concentration
dS = Qin/Vl*(Sin-S) - qS*Xv;
% Ethanol concentration
dP = Qin/Vl*(-P) + qP*Xv;
% Disolved oxygen
dOliq = Qin/Vl*(Osat - Oliq) + kla*(Osat-Oliq) - qO*Xv;
% Oxygen gas phase
dOgas = Fair/Vg*(Ogasin-Ogas) - Vl*kla/Vg*(Osat - Oliq) + Ogas*(Qin-Qe)/Vg;
% Energy balances
%---------------
% Bioreactor temprature
dT = Qin/Vl*(Tin-T) - Tref/Vl*(Qin-Qe) + qO*Xv*DeltaH/MO/rho/Chbr - KT*AT*(T-Tc)/Vl/rho/Chbr;
% Cooling agent temperature
dTc = Fc/Vcj*(Tcin-Tc) + KT*AT*(T-Tc)/Vcj/rhoc/Chc;
% Yields & Productivity
%---------------------
% Cumulative amount of glucose fed to the bioreactor
dSf_cum = Sin*Qin;
dTime = 1;
% Definition of state derivatives vector
% State derivatives
dxdt = [dXt; dXv; dS; dP; dOliq; dOgas; dT; dTc; dVl; dSf_cum; dTime];
return
%Definition of model parameters
%Kinetic parameters
a1 = 0.05; % Ratkowsky parameter [oC-1 h-0.5]
aP = 4.50; % Growth-associated parameter for ethanol production, [-]
AP1 = 6.0; % Activation energy parameter for ethanol production, [oC]
AP2 = 20.3; % Activation energy parameter for ethanol production, [oC]
b1 = 0.035; % Parameter in the exponential expression of the maximum specific growth rate expression, [oC-1]
b2 = 0.15; % Parameter in the exponential expression of the growth inhibitory ethanol concentration expression, [oC-1]
b3 = 0.40; % Parameter in the exponential expression of the specific death rate expression,[oC-1]
c1 = 0.38; % Constant decoupling factor for ethanol production, [gP gX-1 h-1]
c2 = 0.29; % Constant decoupling factor for ethanol production, [gP gX-1 h-1]
k1 = 3.00; % Parameter in the maximum specific growth rate expression, [oC]
k2 = 55.0; % Parameter in the maximum specific growth rate expression, [oC]
k3 = 60.0; % Parameter in the growth-inhibitory ethanol concentration expression, [oC]
k4 = 50.0; % Temperature at the inflection point of the specific death rate sigmoid curve, [oC]
Pmaxb = 90; % Temperature-independent product inhibition constant, [g L-1]
PmaxT = 90; % Maximum value of product inhibition constant due to temperature, [g L-1]
Kdb = 0.025; % Basal specific cellular biomass death rate, [h-1]
KdT = 30.00; % Maximum value of specific cellular biomass death rate due to temperature, [h-1]
KSX = 5; % Glucose saturation constant for the specific growth rate, [g L-1]
KOX = 0.0005; % Oxygen saturation constant for the specific growth rate, [g L-1]
qOmax = 0.05; % Maximum specific oxygen consumption rate, [h-1]
%Metabolic parameters
YPS = 0.51; % Theoretical yield of ethanol on glucose, [gP gS-1]
YXO = 0.97; % Theoretical yield of biomass on oxygen, [gX gO-1]
YXS = 0.53; % Theoretical yield of biomass on glucose, [gX gS-1]
%Physicochemical and thermodynamic parameters
Chbr = 4.18; % Heat capacity of the mass of reaction, [J g-1 oC-1]
Chc = 4.18; % Heat capacity of the cooling agent, [J g-1 oC-1]
DeltaH = 518000; % Heat of reaction of fermentation, [J mol-1 O2]
Tref = 20; % Reference temperature, [oC]
KH = 200; % Henry's constant for oxygen in the fermentation broth, [atm L mol-1]
z = 0.792; % Oxygen compressibility factor, [-]
R = 0.082; % Ideas gas constant, [L atm mol-1 oC-1]
kla0 = 100 ; % Temperature-independent volumetric oxygen transfer coefficient, [h-1]
KT = 360000; % Heat transfer coefficient, [J h-1 m-2 ??C-1]
rho = 1080 ; % Density of the fermentation broth, [g L-1]
rhoc = 1000 ; % Density of the cooling agent, [g L-1]
MO = 32.0; % Molecular weight of oxygen (O2), [g mol-1]
%Bioreactor design data
AT = 1.0; % Bioreactor heat transfer area, [m2]
V = 1800; % Bioreactor working volume, [L]
Vcj = 50; % Cooling jacket volume, [L]
Ogasin = 0.305; % Oxygen concentration in airflow inlet, [g L-1]
%Definition of model variables
%State variables
Xt = x(1); % Total cellular biomass, [g L-1]
Xv = x(2); % Viable cellular biomass, [g L-1]
S = x(3); % Substrate/Glucose concentration, [g L-1]
P = x(4); % Product/Ethanol concentration, [g L-1]
Oliq = x(5); % Dissolved oxygen concentration, [g L-1]
Ogas = x(6); % Gas phase oxygen (bubbles) in the fermentation broth, [g L-1]
T = x(7); % Temperature in the bioreactor, [oC]
Tc = x(8); % Temperature of the cooling agent in the jacket, [oC]
Vl = x(9); % Culture volume in the bioreactor, [L]
Sf_cum = x(10); % Cumulative amount of substrate/glucose fed to the bioreactor, [g]
Time = x(11); % Batch time, [h]
%Initial Conditions
Xt0 = u(10); % Initial total cellular biomass, [g L-1]
Xv0 = u(11); % Initial viable cellular biomass, [g L-1]
S0 = u(12); % Initial substrate/Glucose concentration, [g L-1]
P0 = u(13); % Initial product/Ethanol concentration, [g L-1]
Oliq0 = u(14); % Initial Dissolved oxygen concentration, [g L-1]
Ogas0 = u(15); % Initial Gas phase oxygen (bubbles) in the fermentation broth, [g L-1]
T0 = u(16); % Initial Temperature in the bioreactor, [oC]
Tc0 = u(17); % Initial Temperature of the cooling agent in the jacket, [oC]
Vl0 = u(18); % Initial Culture volume in the bioreactor, [L]
Sf_cum0 = u(19); % Initial Cumulative substrate/glucose fed to the bioreactor, [g]
Time0 = u(20); % Initial batch time, [h]
%Control variables
Qin0 = u(1); % Volumetric inflow rate, [l/h-1]
Qin= Qin0 + 15*heaviside(t-5) + 5*heaviside(t-10) - 6*heaviside(t-20) - 14*heaviside(t-35);
Xtin = u(2); % Total biomass concentration in the bioreactor feed, [g L-1]
Xvin = u(3); % Viable biomass concentration in the bioreactor feed, [g L-1]
Qe = u(4); % Volumetric outflow rate, [l/h-1]
Sin = u(5); % Substrate/Glucose concentration in bioreactor feed, [g L-1]
Fc = u(6); % Cooling agent inlet volumetric flowrate, [L h-1]
Fair = u(7); % Airflow inlet volumetric flowrate, [L h-1]
Tin = u(8); % Temperature of bioreactor feed, [oC]
Tcin = u(9); % Temperature of cooling agent inlet, [oC]
% Definition of model equations
% Kinetic rates
% -----------------------------
% Specific growth rate, [h-1]
mmax = ((a1*(T-k1))*(1-exp(b1*(T-k2))))^2;
Pmax = Pmaxb + PmaxT/(1-exp(-b2*(T-k3)));
m1 = mmax * S/(KSX + S) * Oliq/(KOX + Oliq) * (1 - P/Pmax) * 1/(1+exp(-(100-S)/1)); % Specific growth rate, [h-1]
if m1 >= 0
m = m1;
else
m=0.0;
end
% Non-growth-associated ethanol specific production rate, [h-1]
if S > 0
bP = c1 * exp(-AP1/T) - c2 * exp(-AP2/T) ; % Non-growth-associated ethanol specific production rate, [h-1]
else
bP = 0.0;
end
qP = aP*m + bP;
% Ethanol consumption specific rate
qS = m/YXS + qP/YPS;
% Oxygen consumption specific rate
qO = qOmax*Oliq/YXO/(KOX + Oliq);
% Specific biological deactivation rate of cell mass
Kd = Kdb + KdT/(1+exp(-b3*(T-k4)));
% Saturation concentration of oxygen in culture media
Osat = z*Ogas*R*T/KH;
% Oxygen mass transfer coefficient
kla = kla0*1.2^(T-20);
% Volume of the gas phase in the bioreactor
Vg = V - Vl;
%Material balances
%-----------------
% Volume of liquid culture
dVl = Qin - Qe;
% Total cell mass
dXt = m*Xv + Qin/Vl*(Xtin-Xt);
% Total mass of biologically active cells
dXv = (m-Kd)*Xv + Qin/Vl*(Xvin-Xv);
% Glucose concentration
dS = Qin/Vl*(Sin-S) - qS*Xv;
% Ethanol concentration
dP = Qin/Vl*(-P) + qP*Xv;
% Disolved oxygen
dOliq = Qin/Vl*(Osat - Oliq) + kla*(Osat-Oliq) - qO*Xv;
% Oxygen gas phase
dOgas = Fair/Vg*(Ogasin-Ogas) - Vl*kla/Vg*(Osat - Oliq) + Ogas*(Qin-Qe)/Vg;
% Energy balances
%---------------
% Bioreactor temprature
dT = Qin/Vl*(Tin-T) - Tref/Vl*(Qin-Qe) + qO*Xv*DeltaH/MO/rho/Chbr - KT*AT*(T-Tc)/Vl/rho/Chbr;
% Cooling agent temperature
dTc = Fc/Vcj*(Tcin-Tc) + KT*AT*(T-Tc)/Vcj/rhoc/Chc;
% Yields & Productivity
%---------------------
% Cumulative amount of glucose fed to the bioreactor
dSf_cum = Sin*Qin;
dTime = 1;
% Definition of state derivatives vector
% State derivatives
dxdt = [dXt; dXv; dS; dP; dOliq; dOgas; dT; dTc; dVl; dSf_cum; dTime];
return
Optimize Ethanol Production
The simulation model of an ethanol bioreactor can be used to optimize the production of ethanol by providing a virtual platform to test different scenarios and strategies without the need for costly and time-consuming experimental trials. Follow these steps to optimize the ethanol production:
- Define the parameters: Once the model is built, the next step is to define the parameters that affect the production of ethanol. These may include variables such as temperature, substrate concentration, batch time, and other relevant variables.
- Run simulations: Using the defined parameters, run simulations to determine the best conditions for maximizing ethanol production. The simulations can be run repeatedly, each time changing one or more variables to determine the optimal combination.
- Analyze the results: After running the simulations, analyze the results to identify the optimal conditions for maximizing ethanol production. The analysis can be done using statistical techniques or data visualization tools to help identify trends and patterns.
- Optimize the bioreactor: Based on the results of the analysis, optimize the bioreactor operation by adjusting the parameters to the identified optimal conditions. This may involve adjusting the reactor temperature, substrate concentration, or other variables.
- Validate the model: test the simulation model by comparing the results obtained from the simulation with experimental data obtained from the actual bioreactor. This ensures that the results of the optimization are accurate.
The optimization strategy is used to determine the optimal Tcin (temperature to the cooling jacket) and Fair (air flow).
Maximize Ethanol Production with Cooling Temperature
Maximize Ethanol Production with Air Flow Rate
import numpy as np
from scipy.integrate import odeint
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
#Simulate fed-batch operation
# Specify the simulation time (hrs)
tspan = np.linspace(0,37,371); t=tspan
# Specify values of control variables [Qin0 Qin Xtin Xvin Qe Sin Fc Fair Tin Tcin]
u0 = [0.0,0.0,0.0,0.0,0.0,400,40,60000,30,15]
# Specify initial conditions [Xt Xv S P Oliq Ogas T Tc Vl Gloss MP]
x0 = [0.1,0.1,50,0,0.0065,0.305,30,20,1000,0.0,0.0]
ux0 = (u0 + x0)
#Qin = 15*heaviside(t-5) + 5*heaviside(t-10)
# - 6*heaviside(t-20) - 14*heaviside(t-35)
Qin = np.zeros_like(tspan)
Qin[np.where(tspan>=5)] += 15
Qin[np.where(tspan>=10)] += 5
Qin[np.where(tspan>=20)] -= 6
Qin[np.where(tspan>=35)] -= 14
QinInterp = interp1d(tspan,Qin,bounds_error=False)
def ethanol(x,t,Qin0,Qin,Xtin,Xvin,Qe,Sin,Fc,Fair,
Tin,Tcin,Xt0,Xv0,S0,P0,Oliq0,Ogas0,T0,
Tc0,Vl0,Sf_cum0,Time0):
## #Initial Conditions
## Xt0 = u[10] # Initial total cellular biomass, [g L-1]
## Xv0 = u[11] # Initial viable cellular biomass, [g L-1]
## S0 = u[12] # Initial substrate/Glucose concentration, [g L-1]
## P0 = u[13] # Initial product/Ethanol concentration, [g L-1]
## Oliq0 = u[14] # Initial Dissolved oxygen concentration, [g L-1]
## Ogas0 = u[15] # Initial Gas phase oxygen (bubbles) in the fermentation broth, [g L-1]
## T0 = u[16] # Initial Temperature in the bioreactor, [oC]
## Tc0 = u[17] # Initial Temperature of the cooling agent in the jacket, [oC]
## Vl0 = u[18] # Initial Culture volume in the bioreactor, [L]
## Sf_cum0 = u[19] # Initial Cumulative substrate/glucose fed to the bioreactor, [g]
## Time0 = u[20] # Initial batch time, [h]
##
## #Control variables
## Qin0 = u[0] # Volumetric inflow rate, [l/h-1]
## Qin = u[1] # Volumetric inflow rate, [l/h-1]
## Xtin = u[2] # Total biomass concentration in the bioreactor feed, [g L-1]
## Xvin = u[3] # Viable biomass concentration in the bioreactor feed, [g L-1]
## Qe = u[4] # Volumetric outflow rate, [l/h-1]
## Sin = u[5] # Substrate/Glucose concentration in bioreactor feed, [g L-1]
## Fc = u[6] # Cooling agent inlet volumetric flowrate, [L h-1]
## Fair = u[7] # Airflow inlet volumetric flowrate, [L h-1]
## Tin = u[8] # Temperature of bioreactor feed, [oC]
## Tcin = u[9] # Temperature of cooling agent inlet, [oC]
# 1D Interpolation for Qin
Qin = QinInterp(t)
#Definition of model parameters
#Kinetic parameters
a1 = 0.05 # Ratkowsky parameter [oC-1 h-0.5]
aP = 4.50 # Growth-associated parameter for ethanol production, [-]
AP1 = 6.0 # Activation energy parameter for ethanol production, [oC]
AP2 = 20.3 # Activation energy parameter for ethanol production, [oC]
b1 = 0.035 # Parameter in the exponential expression of the maximum specific growth rate np.expression, [oC-1]
b2 = 0.15 # Parameter in the exponential expression of the growth inhibitory ethanol concentration np.expression, [oC-1]
b3 = 0.40 # Parameter in the exponential np.expression of the specific death rate expression,[oC-1]
c1 = 0.38 # Constant decoupling factor for ethanol production, [gP gX-1 h-1]
c2 = 0.29 # Constant decoupling factor for ethanol production, [gP gX-1 h-1]
k1 = 3.00 # Parameter in the maximum specific growth rate expression, [oC]
k2 = 55.0 # Parameter in the maximum specific growth rate expression, [oC]
k3 = 60.0 # Parameter in the growth-inhibitory ethanol concentration expression, [oC]
k4 = 50.0 # Temperature at the inflection point of the specific death rate sigmoid curve, [oC]
Pmaxb = 90 # Temperature-independent product inhibition constant, [g L-1]
PmaxT = 90 # Maximum value of product inhibition constant due to temperature, [g L-1]
Kdb = 0.025 # Basal specific cellular biomass death rate, [h-1]
KdT = 30.00 # Maximum value of specific cellular biomass death rate due to temperature, [h-1]
KSX = 5 # Glucose saturation constant for the specific growth rate, [g L-1]
KOX = 0.0005 # Oxygen saturation constant for the specific growth rate, [g L-1]
qOmax = 0.05 # Maximum specific oxygen consumption rate, [h-1]
#Metabolic parameters
YPS = 0.51 # Theoretical yield of ethanol on glucose, [gP gS-1]
YXO = 0.97 # Theoretical yield of biomass on oxygen, [gX gO-1]
YXS = 0.53 # Theoretical yield of biomass on glucose, [gX gS-1]
#Physicochemical and thermodynamic parameters
Chbr = 4.18 # Heat capacity of the mass of reaction, [J g-1 oC-1]
Chc = 4.18 # Heat capacity of the cooling agent, [J g-1 oC-1]
DeltaH = 518000 # Heat of reaction of fermentation, [J mol-1 O2]
Tref = 20 # Reference temperature, [oC]
KH = 200 # Henry's constant for oxygen in the fermentation broth, [atm L mol-1]
z = 0.792 # Oxygen compressibility factor, [-]
R = 0.082 # Ideas gas constant, [L atm mol-1 oC-1]
kla0 = 100 # Temperature-independent volumetric oxygen transfer coefficient, [h-1]
KT = 360000 # Heat transfer coefficient, [J h-1 m-2 ??C-1]
rho = 1080 # Density of the fermentation broth, [g L-1]
rhoc = 1000 # Density of the cooling agent, [g L-1]
MO = 32.0 # Molecular weight of oxygen (O2), [g mol-1]
#Bioreactor design data
AT = 1.0 # Bioreactor heat transfer area, [m2]
V = 1800 # Bioreactor working volume, [L]
Vcj = 50 # Cooling jacket volume, [L]
Ogasin = 0.305 # Oxygen concentration in airflow inlet, [g L-1]
#Definition of model variables
#State variables
Xt = x[0] # Total cellular biomass, [g L-1]
Xv = x[1] # Viable cellular biomass, [g L-1]
S = x[2] # Substrate/Glucose concentration, [g L-1]
P = x[3] # Product/Ethanol concentration, [g L-1]
Oliq = x[4] # Dissolved oxygen concentration, [g L-1]
Ogas = x[5] # Gas phase oxygen (bubbles) in the fermentation broth, [g L-1]
T = x[6] # Temperature in the bioreactor, [oC]
Tc = x[7] # Temperature of the cooling agent in the jacket, [oC]
Vl = x[8] # Culture volume in the bioreactor, [L]
Sf_cum = x[9] # Cumulative amount of substrate/glucose fed to the bioreactor, [g]
Time = x[10] # Batch time, [h]
# Definition of model equations
# Kinetic rates
# -----------------------------
# Specific growth rate, [h-1]
mmax = ((a1*(T-k1))*(1-np.exp(b1*(T-k2))))**2
Pmax = Pmaxb + PmaxT/(1-np.exp(-b2*(T-k3)))
m1 = mmax * S/(KSX + S) * Oliq/(KOX + Oliq) * (1 - P/Pmax) * 1/(1+np.exp(-(100-S)/1)) # Specific growth rate, [h-1]
if m1 >= 0:
m = m1
else:
m=0.0
# Non-growth-associated ethanol specific production rate, [h-1]
if S > 0:
bP = c1 * np.exp(-AP1/T) - c2 * np.exp(-AP2/T) # Non-growth-associated ethanol specific production rate, [h-1]
else:
bP = 0.0
qP = aP*m + bP
# Ethanol consumption specific rate
qS = m/YXS + qP/YPS
# Oxygen consumption specific rate
qO = qOmax*Oliq/YXO/(KOX + Oliq)
# Specific biological deactivation rate of cell mass
Kd = Kdb + KdT/(1+np.exp(-b3*(T-k4)))
# Saturation concentration of oxygen in culture media
Osat = z*Ogas*R*T/KH
# Oxygen mass transfer coefficient
kla = kla0*1.2**(T-20)
# Volume of the gas phase in the bioreactor
Vg = V - Vl
#Material balances
#-----------------
# Volume of liquid culture
dVl = Qin - Qe
# Total cell mass
dXt = m*Xv + Qin/Vl*(Xtin-Xt)
# Total mass of biologically active cells
dXv = (m-Kd)*Xv + Qin/Vl*(Xvin-Xv)
# Glucose concentration
dS = Qin/Vl*(Sin-S) - qS*Xv
# Ethanol concentration
dP = Qin/Vl*(-P) + qP*Xv
# Disolved oxygen
dOliq = Qin/Vl*(Osat - Oliq) + kla*(Osat-Oliq) - qO*Xv
# Oxygen gas phase
dOgas = Fair/Vg*(Ogasin-Ogas) - Vl*kla/Vg*(Osat - Oliq) + Ogas*(Qin-Qe)/Vg
# Energy balances
#---------------
# Bioreactor temprature
dT = Qin/Vl*(Tin-T) - Tref/Vl*(Qin-Qe) + qO*Xv*DeltaH/MO/rho/Chbr - KT*AT*(T-Tc)/Vl/rho/Chbr
# Cooling agent temperature
dTc = Fc/Vcj*(Tcin-Tc) + KT*AT*(T-Tc)/Vcj/rhoc/Chc
# Yields & Productivity
#---------------------
# Cumulative amount of glucose fed to the bioreactor
dSf_cum = Sin*Qin
dTime = 1
# Definition of state derivatives vector
# State derivatives
dxdt = [dXt,dXv,dS,dP,dOliq,dOgas,dT,dTc,dVl,dSf_cum,dTime]
return dxdt
T_val = np.linspace(7,40,21)
P_val = np.empty_like(T_val)
ux1 = ux0
for i,T in enumerate(T_val):
ux1[9] = T
x = odeint(ethanol,x0,tspan,args=tuple(ux1))
P_val[i] = x[-3,3]
plt.figure(figsize=(6,3))
plt.plot(T_val, P_val, label='Maximize Product')
plt.legend(); plt.grid()
plt.ylabel('Conc [g/L]')
plt.xlabel('Tcin')
plt.tight_layout()
plt.savefig('bioreactor_Tcin.png',dpi=300)
F_val = np.linspace(10,1000,41)
P_val = np.empty_like(F_val)
ux2 = ux0
for i,F in enumerate(F_val):
ux2[7] = F
x = odeint(ethanol,x0,tspan,args=tuple(ux2))
P_val[i] = x[-3,3]
plt.figure(figsize=(6,3))
plt.plot(F_val, P_val, label='Maximize Product')
plt.legend(); plt.grid()
plt.ylabel('Conc [g/L]')
plt.xlabel('Airflow [L/hr]')
plt.tight_layout()
plt.savefig('bioreactor_Fair.png',dpi=300)
plt.show()
from scipy.integrate import odeint
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
#Simulate fed-batch operation
# Specify the simulation time (hrs)
tspan = np.linspace(0,37,371); t=tspan
# Specify values of control variables [Qin0 Qin Xtin Xvin Qe Sin Fc Fair Tin Tcin]
u0 = [0.0,0.0,0.0,0.0,0.0,400,40,60000,30,15]
# Specify initial conditions [Xt Xv S P Oliq Ogas T Tc Vl Gloss MP]
x0 = [0.1,0.1,50,0,0.0065,0.305,30,20,1000,0.0,0.0]
ux0 = (u0 + x0)
#Qin = 15*heaviside(t-5) + 5*heaviside(t-10)
# - 6*heaviside(t-20) - 14*heaviside(t-35)
Qin = np.zeros_like(tspan)
Qin[np.where(tspan>=5)] += 15
Qin[np.where(tspan>=10)] += 5
Qin[np.where(tspan>=20)] -= 6
Qin[np.where(tspan>=35)] -= 14
QinInterp = interp1d(tspan,Qin,bounds_error=False)
def ethanol(x,t,Qin0,Qin,Xtin,Xvin,Qe,Sin,Fc,Fair,
Tin,Tcin,Xt0,Xv0,S0,P0,Oliq0,Ogas0,T0,
Tc0,Vl0,Sf_cum0,Time0):
## #Initial Conditions
## Xt0 = u[10] # Initial total cellular biomass, [g L-1]
## Xv0 = u[11] # Initial viable cellular biomass, [g L-1]
## S0 = u[12] # Initial substrate/Glucose concentration, [g L-1]
## P0 = u[13] # Initial product/Ethanol concentration, [g L-1]
## Oliq0 = u[14] # Initial Dissolved oxygen concentration, [g L-1]
## Ogas0 = u[15] # Initial Gas phase oxygen (bubbles) in the fermentation broth, [g L-1]
## T0 = u[16] # Initial Temperature in the bioreactor, [oC]
## Tc0 = u[17] # Initial Temperature of the cooling agent in the jacket, [oC]
## Vl0 = u[18] # Initial Culture volume in the bioreactor, [L]
## Sf_cum0 = u[19] # Initial Cumulative substrate/glucose fed to the bioreactor, [g]
## Time0 = u[20] # Initial batch time, [h]
##
## #Control variables
## Qin0 = u[0] # Volumetric inflow rate, [l/h-1]
## Qin = u[1] # Volumetric inflow rate, [l/h-1]
## Xtin = u[2] # Total biomass concentration in the bioreactor feed, [g L-1]
## Xvin = u[3] # Viable biomass concentration in the bioreactor feed, [g L-1]
## Qe = u[4] # Volumetric outflow rate, [l/h-1]
## Sin = u[5] # Substrate/Glucose concentration in bioreactor feed, [g L-1]
## Fc = u[6] # Cooling agent inlet volumetric flowrate, [L h-1]
## Fair = u[7] # Airflow inlet volumetric flowrate, [L h-1]
## Tin = u[8] # Temperature of bioreactor feed, [oC]
## Tcin = u[9] # Temperature of cooling agent inlet, [oC]
# 1D Interpolation for Qin
Qin = QinInterp(t)
#Definition of model parameters
#Kinetic parameters
a1 = 0.05 # Ratkowsky parameter [oC-1 h-0.5]
aP = 4.50 # Growth-associated parameter for ethanol production, [-]
AP1 = 6.0 # Activation energy parameter for ethanol production, [oC]
AP2 = 20.3 # Activation energy parameter for ethanol production, [oC]
b1 = 0.035 # Parameter in the exponential expression of the maximum specific growth rate np.expression, [oC-1]
b2 = 0.15 # Parameter in the exponential expression of the growth inhibitory ethanol concentration np.expression, [oC-1]
b3 = 0.40 # Parameter in the exponential np.expression of the specific death rate expression,[oC-1]
c1 = 0.38 # Constant decoupling factor for ethanol production, [gP gX-1 h-1]
c2 = 0.29 # Constant decoupling factor for ethanol production, [gP gX-1 h-1]
k1 = 3.00 # Parameter in the maximum specific growth rate expression, [oC]
k2 = 55.0 # Parameter in the maximum specific growth rate expression, [oC]
k3 = 60.0 # Parameter in the growth-inhibitory ethanol concentration expression, [oC]
k4 = 50.0 # Temperature at the inflection point of the specific death rate sigmoid curve, [oC]
Pmaxb = 90 # Temperature-independent product inhibition constant, [g L-1]
PmaxT = 90 # Maximum value of product inhibition constant due to temperature, [g L-1]
Kdb = 0.025 # Basal specific cellular biomass death rate, [h-1]
KdT = 30.00 # Maximum value of specific cellular biomass death rate due to temperature, [h-1]
KSX = 5 # Glucose saturation constant for the specific growth rate, [g L-1]
KOX = 0.0005 # Oxygen saturation constant for the specific growth rate, [g L-1]
qOmax = 0.05 # Maximum specific oxygen consumption rate, [h-1]
#Metabolic parameters
YPS = 0.51 # Theoretical yield of ethanol on glucose, [gP gS-1]
YXO = 0.97 # Theoretical yield of biomass on oxygen, [gX gO-1]
YXS = 0.53 # Theoretical yield of biomass on glucose, [gX gS-1]
#Physicochemical and thermodynamic parameters
Chbr = 4.18 # Heat capacity of the mass of reaction, [J g-1 oC-1]
Chc = 4.18 # Heat capacity of the cooling agent, [J g-1 oC-1]
DeltaH = 518000 # Heat of reaction of fermentation, [J mol-1 O2]
Tref = 20 # Reference temperature, [oC]
KH = 200 # Henry's constant for oxygen in the fermentation broth, [atm L mol-1]
z = 0.792 # Oxygen compressibility factor, [-]
R = 0.082 # Ideas gas constant, [L atm mol-1 oC-1]
kla0 = 100 # Temperature-independent volumetric oxygen transfer coefficient, [h-1]
KT = 360000 # Heat transfer coefficient, [J h-1 m-2 ??C-1]
rho = 1080 # Density of the fermentation broth, [g L-1]
rhoc = 1000 # Density of the cooling agent, [g L-1]
MO = 32.0 # Molecular weight of oxygen (O2), [g mol-1]
#Bioreactor design data
AT = 1.0 # Bioreactor heat transfer area, [m2]
V = 1800 # Bioreactor working volume, [L]
Vcj = 50 # Cooling jacket volume, [L]
Ogasin = 0.305 # Oxygen concentration in airflow inlet, [g L-1]
#Definition of model variables
#State variables
Xt = x[0] # Total cellular biomass, [g L-1]
Xv = x[1] # Viable cellular biomass, [g L-1]
S = x[2] # Substrate/Glucose concentration, [g L-1]
P = x[3] # Product/Ethanol concentration, [g L-1]
Oliq = x[4] # Dissolved oxygen concentration, [g L-1]
Ogas = x[5] # Gas phase oxygen (bubbles) in the fermentation broth, [g L-1]
T = x[6] # Temperature in the bioreactor, [oC]
Tc = x[7] # Temperature of the cooling agent in the jacket, [oC]
Vl = x[8] # Culture volume in the bioreactor, [L]
Sf_cum = x[9] # Cumulative amount of substrate/glucose fed to the bioreactor, [g]
Time = x[10] # Batch time, [h]
# Definition of model equations
# Kinetic rates
# -----------------------------
# Specific growth rate, [h-1]
mmax = ((a1*(T-k1))*(1-np.exp(b1*(T-k2))))**2
Pmax = Pmaxb + PmaxT/(1-np.exp(-b2*(T-k3)))
m1 = mmax * S/(KSX + S) * Oliq/(KOX + Oliq) * (1 - P/Pmax) * 1/(1+np.exp(-(100-S)/1)) # Specific growth rate, [h-1]
if m1 >= 0:
m = m1
else:
m=0.0
# Non-growth-associated ethanol specific production rate, [h-1]
if S > 0:
bP = c1 * np.exp(-AP1/T) - c2 * np.exp(-AP2/T) # Non-growth-associated ethanol specific production rate, [h-1]
else:
bP = 0.0
qP = aP*m + bP
# Ethanol consumption specific rate
qS = m/YXS + qP/YPS
# Oxygen consumption specific rate
qO = qOmax*Oliq/YXO/(KOX + Oliq)
# Specific biological deactivation rate of cell mass
Kd = Kdb + KdT/(1+np.exp(-b3*(T-k4)))
# Saturation concentration of oxygen in culture media
Osat = z*Ogas*R*T/KH
# Oxygen mass transfer coefficient
kla = kla0*1.2**(T-20)
# Volume of the gas phase in the bioreactor
Vg = V - Vl
#Material balances
#-----------------
# Volume of liquid culture
dVl = Qin - Qe
# Total cell mass
dXt = m*Xv + Qin/Vl*(Xtin-Xt)
# Total mass of biologically active cells
dXv = (m-Kd)*Xv + Qin/Vl*(Xvin-Xv)
# Glucose concentration
dS = Qin/Vl*(Sin-S) - qS*Xv
# Ethanol concentration
dP = Qin/Vl*(-P) + qP*Xv
# Disolved oxygen
dOliq = Qin/Vl*(Osat - Oliq) + kla*(Osat-Oliq) - qO*Xv
# Oxygen gas phase
dOgas = Fair/Vg*(Ogasin-Ogas) - Vl*kla/Vg*(Osat - Oliq) + Ogas*(Qin-Qe)/Vg
# Energy balances
#---------------
# Bioreactor temprature
dT = Qin/Vl*(Tin-T) - Tref/Vl*(Qin-Qe) + qO*Xv*DeltaH/MO/rho/Chbr - KT*AT*(T-Tc)/Vl/rho/Chbr
# Cooling agent temperature
dTc = Fc/Vcj*(Tcin-Tc) + KT*AT*(T-Tc)/Vcj/rhoc/Chc
# Yields & Productivity
#---------------------
# Cumulative amount of glucose fed to the bioreactor
dSf_cum = Sin*Qin
dTime = 1
# Definition of state derivatives vector
# State derivatives
dxdt = [dXt,dXv,dS,dP,dOliq,dOgas,dT,dTc,dVl,dSf_cum,dTime]
return dxdt
T_val = np.linspace(7,40,21)
P_val = np.empty_like(T_val)
ux1 = ux0
for i,T in enumerate(T_val):
ux1[9] = T
x = odeint(ethanol,x0,tspan,args=tuple(ux1))
P_val[i] = x[-3,3]
plt.figure(figsize=(6,3))
plt.plot(T_val, P_val, label='Maximize Product')
plt.legend(); plt.grid()
plt.ylabel('Conc [g/L]')
plt.xlabel('Tcin')
plt.tight_layout()
plt.savefig('bioreactor_Tcin.png',dpi=300)
F_val = np.linspace(10,1000,41)
P_val = np.empty_like(F_val)
ux2 = ux0
for i,F in enumerate(F_val):
ux2[7] = F
x = odeint(ethanol,x0,tspan,args=tuple(ux2))
P_val[i] = x[-3,3]
plt.figure(figsize=(6,3))
plt.plot(F_val, P_val, label='Maximize Product')
plt.legend(); plt.grid()
plt.ylabel('Conc [g/L]')
plt.xlabel('Airflow [L/hr]')
plt.tight_layout()
plt.savefig('bioreactor_Fair.png',dpi=300)
plt.show()
References
- Bioreactor Simulation for Ethanol Production using GEKKO, Stack Overflow, Nov 9, 2022.
