View source for Assignment 6 - 2010

{{OtherSidebar
| due_dates = 01 December 2010
| dates_alt_text = Due date(s)
| questions_PDF = Assignment-6-2010.pdf
| questions_text_alt = Assignment questions
| solutions_PDF = 
| solutions_text_alt = Solutions, prepared by Ali, Elliot and Kevin
| other_instructions = Hand-in at class.
}}

<rst>
<rst-options: 'toc' = False/>
<rst-options: 'reset-figures' = False/>
.. |m3| replace:: m\ :sup:`3`

.. rubric:: To integrate ODEs - both by hand and with computer software.
 	
.. Potential questions:

	Exothermic reaction in a PFR; does temperature exceed 300C over the distance of the pipe?  (See reactor design project)

	Epidemic example, p 557 in Friendly book
	
	Stiff system?
	
	
Question 1 [2]
===============

In assignment 1 you derived the dynamic balance for a waste water treatment system (aerobic heterotrophic growth using Monod kinetics); and in assignment 3 you used Newton's method to find the steady-state solution to the nutrient (:math:`N`) and biomass (:math:`B`) balances:

.. math::
		\frac{dN}{dt} &= N_\text{in}(t)\frac{F(t)}{V} - N(t)\frac{F(t)}{V} - \left(\frac{1}{Y_B}\right)\left(\mu_{\text{max}}\frac{N}{K+N(t)}\right)B(t) \\
		\frac{dB}{dt} &= - B(t)\frac{F(t)}{V} + \left(\mu_{\text{max}}\frac{N(t)}{K+N(t)}\right)B(t)

*	The tank volume is assumed constant at :math:`V = 1600` |m3|
*	The conversion efficiency is constant at :math:`Y_B=0.80`
*	The two Monod kinetics parameters are: :math:`\mu_\text{max} = 5 \,\,\text{day}^{-1}` and :math:`K = 20` g/|m3|
*	The inlet flow rate is :math:`F(t) = 5000` |m3|/day
*	The nutrient inflow is :math:`N_\text{in}(t) = 100` g/|m3|

These dynamic balances describe how nutrients in the inlet (i.e. the :math:`N_\text{in} = 100` g/|m3| of waste in the wastewater) are converted to biomass, B. The stream leaving the well-mixed CSTR should have lower nutrient levels, :math:`N(t)`.

#.	Code these equations in a MATLAB or Python function and use the built-in integrators in either software package to find the steady state solutions when using starting values of :math:`N(t=0) = 100` g/|m3| and :math:`B(t) = 20` g/|m3| and show that you obtain the same steady state values as found in assignment 3.

#.	Most real processes are never at steady state.  In this reactor a (slightly) more realistic input, to account for the diurnal nature of the flow, is: 

	.. math::
	
		F(t) = 5000+900\sin\left(\displaystyle \frac{2\pi}{24}t -\frac{\pi}{5}\right)

	*	Plot 100 hours of process operation at these conditions; use the steady state values from part 1 as your initial conditions.  
	*	Environmental limits prohibit outlet nutrient concentration from exceeding 60 g/|m3|.  Does this system meet those standards?

#.	Using the simulation from part 2 (with the sinusoidal input) for the biomass and nutrient profiles, see what happens if your plant accidentally dumps nutrient into the inlet stream for 10 hours when an upstream unit fails.  Use an increased nutrient influx of :math:`N_\text{in}(t) = 150` g/|m3| between :math:`55 \leq t \leq 65` hours.

	*	Does the plant still manage to treat this accident without violating the 60 g/|m3| outlet limit on :math:`N(t)`?
	*	Would it handle this problem had it occurred between :math:`75 \leq t \leq 85` hours?

.. Solution:  see code/wastewater.py file

Question 2 [3]
===============

.. Reactor system, with recycle loop; disturbance in middle reactor at t = 10; see the disturbance propagate through the system and the recycle loop

.. Series of CSTR's with reentry into middle CSTR

Consider three CSTR's in series, operating at constant volume of :math:`V = 100` L.  A first order reaction, :math:`{\sf A \rightarrow B + C}` takes place in the system, with reaction rate = :math:`r = 0.09 C_{\sf A}`.  The inlet conditions to the first reactor are constant with time: inlet concentration = :math:`C_{{\sf A},0} = 2.7` mol/L at a flow rate of :math:`F_\text{in}^{0}` = 15.8 L/min.  The outlet from the first reactor, :math:`C_{{\sf A},1}`, is the inlet to the second reactor, and the outlet from the second reactor, :math:`C_{{\sf A},2}`, is passed to third reactor.  The outlet from the third reactor, :math:`C_{{\sf A},3}`, can be used to calculate the overall conversion = :math:`\displaystyle \frac{C_{{\sf A},0} - C_{{\sf A},3}}{C_{{\sf A},0}}`.

The dynamic balances in each reactor are:

.. math::
	
	\frac{d C_{{\sf A},1}(t)}{dt} &= \frac{F}{V_1}\bigg(C_{{\sf A},0} - C_{{\sf A},1} \bigg) - kC_{{\sf A},1} \\
	\frac{d C_{{\sf A},2}(t)}{dt} &= \frac{F}{V_2}\bigg(C_{{\sf A},1} - C_{{\sf A},2} \bigg) - kC_{{\sf A},2} \\
	\frac{d C_{{\sf A},3}(t)}{dt} &= \frac{F}{V_3}\bigg(C_{{\sf A},2} - C_{{\sf A},3} \bigg) - kC_{{\sf A},3} \\
	
1.	Describe how you would ordinarily solve the system of equations :math:`f({\bf x}) = {\bf 0}` where :math:`{\bf x} = \big[C_{{\sf A},1}, C_{{\sf A},2}, C_{{\sf A},3} \big]` to find the steady-state values of :math:`\big[C_{{\sf A},1}, C_{{\sf A},2}, C_{{\sf A},3} \big]`.   Now use a computer-based ODE solver to integrate the equations to steady state and report the steady-state values of :math:`{\bf x} = \big[C_{{\sf A},1}, C_{{\sf A},2}, C_{{\sf A},3} \big]` and the overall conversion.

	.. [1.720, 1.095, 0.698], conversion = (2.7 - 0.698)/2.7 = 74.1%

You are tasked with increasing the conversion of the overall system.  You are investigating two possibilities: adding a recycle stream, and operating at larger tank volumes.

2.	What are the new steady state values of :math:`{\bf x} = \big[C_{{\sf A},1}, C_{{\sf A},2}, C_{{\sf A},3} \big]` and the overall conversion when :math:`V_1 = V_2 = V_3 = 150` L?  Why does operating at a larger tank volumes increase conversion?  You can answer this question even if you haven't taken the reactor design course yet - *hint*: interpret the above modelling equations.

	.. [1.456, 0.785, 0.423], conversion = (2.7 - 0.423)/2.7 = 84.3%

3.	Let's take a look at what happens when a disturbance enters the system.  Simulate the case when a step change of +1.2 mol/L in :math:`C_{{\sf A},0}` is made at time :math:`t = 50` (i.e. the inlet concentration is 2.7 mol/L for 50 minutes, then increased to 3.9 mol/L for the rest of the time).  Choose suitable initial conditions so that your system is roughly at steady state before :math:`t=50`, continue the simulation on until the new steady state is reached.

	What do you notice about the shape of the concentration *vs* time trajectories from each reactor (superimpose the 3 trajectories on the same plot).


Adding a recycle stream from the final reactor to the first or second reactor is an exercise that you can complete on your own.  It leads to a counter-intuitive phenomenon: you'd expect an increase in conversion, but ....

Question 3 [2]
================

.. note:: This question is adapted from the textbook by Bradie, listed under the `Suggested readings <http://modelling3e4.connectmv.com/wiki/Suggested_readings>`_ sections of the course website.

A genetic switch is a biological mechanism to activate whether or not a particular protein product from a cell is synthesized.  For example, under certain conditions a pigment protein can be synthesized, which indicates the presence/absence of another chemical (the switch).

This model has been proposed for a genetic switch:

.. math::

	\frac{dg}{dt} = 0.01s - 1.51g + 3.03 \frac{g^2}{1+g^2}\qquad g(0) = 0.0

In this notation, the pigment protein concentration is :math:`g(t)` and :math:`s(t)` is the concentration of the chemical (the switch) that activates the gene that produces the pigment.

The :math:`-1.51g` term indicates the protein's natural decay rate is proportional to its concentration, while the last term describes the positive feedback  the protein exerts on its own formation.

A switch has two important characteristics:

	*	Below a threshold value of :math:`s` the gene concentration must be close to zero (i.e. little/no pigment produced means the switch is off), while above a threshold, the gene concentration, :math:`g`, must be high (i.e. the switch has been turned on).
	
	*	Once turned on (i.e. high values of :math:`g`), the switch must remain on, even if :math:`s` is returned below the threshold value.  Similarly, if the switch was never on to begin with, it must still remain "off" if :math:`s` is below the threshold.  This phenomenon is known as the hysteresis effect.
	

#.	First we wish to find the threshold value of :math:`s` and interpret what :math:`s` means.  Plot and label four trajectories of :math:`g(t)` using :math:`s=20.0, 20.33, 20.67, 21.0` over a time range of 100 minutes, and describe what the plots are showing.  

	*	Be specific: how would you use these plots to design your genetic switch?  
	*	What does :math:`s` mean in your design? (Hint: it has two interpretations)

#.	Now we wish to verify if the hysteresis effect is present, because it is important the switch cannot be turned off.  For each of the four values of :math:`s` in part 1, run a simulation where :math:`s(t>150)=0` (i.e. at time :math:`t=150` the concentration of :math:`s` is set from its current value to zero).  Show your results over a time horizon of 300 minutes and **comment on the results**.

.. Solution:: code/genetic_switch.py

Question 4 [2]
==============

.. note:: This question continues on from question 5, assignment 5.  It was adapted from a question in the final exam, 2009, worth 3+1 out of 50 marks.

Consider a well-mixed tank containing species :math:`{\sf A}` at concentration :math:`C_{\sf A}=1` mol/L. There is no inlet or outlet from the tank, i.e. it is a batch system. A fluid of constant density is considered throughout this question.  The reactor is filled to 5.2 m\ :sup:`3`.

A catalyst is introduced in the tank at time :math:`t=0` minutes to initiate the reaction :math:`{\sf A} \rightarrow {\sf B} + {\sf C}`. Due to quick catalyst deactivation, the dynamic equation describing the evolution of concentration :math:`C_{\sf A}` in the tank is given by:

.. math::

	\frac{{\rm d}C_{\sf A}(t)}{{\rm d}t} = - k\,C_{\sf A}(t)\,e^{-\frac{t}{2}}
	
with :math:`k=1\ \rm min^{-1}`.

#.	Determine the concentration of A in the reactor at time :math:`t=20` minutes using the classical Runge-Kutta method with step size :math:`h=10` minutes.

#.	Compare and discuss your answer with respect to the analytical solution: 

	.. math::
	
		C_{\sf A}(t) = C_{\sf A}(t=0)\ \exp\left[2\,k\left(e^{-\frac{t}{2}}-1\right)\right]
	

.. Solution: see code/catalyst.py

.. raw:: latex

	\newpage

Question 5 [3]
==============

.. From Tutorial 9, 2008

Consider the CSTR shown below, where an irreversible reaction of the form :math:`{\sf A}\rightarrow {\sf B}` takes place.

.. figure:: ../che3e4/Assignments/Assignment-6/images/CSTR_1input.png
	:scale: 68
	:align: center
	
**Model parameters**

* :math:`C_{\sf A}^{\rm in} = 0.79 \rm kmol\ m^{-3}`
* :math:`T^{\rm in} = 352.6 \rm K`
* :math:`F = 0.2 \rm m^3\ min^{-1}` 
* :math:`V = 0.1 \rm m^3` 
* :math:`k_{0} = 7.2\times 10^{10} \rm min^{-1}` 
* :math:`E = 8.314\times 10^{4} \rm kJ\ kmol{-1}` 
* :math:`\Delta H = -4.78\times 10^4 \rm kJ\ kmol^{-1}` 
* :math:`C_{p} = 0.239 \rm kJ\ kg^{-1}\ K^{-1}` 
* :math:`\rho = 1000 \rm kg\ m^{-3}` 
* :math:`R = 8.314 \rm kJ\ kmol^{-1}\ K^{-1}` 

#.	Derive, by hand, the differential equations describing the evolution of the concentration of species :math:`{\sf A}`, :math:`C_{\sf A}(t)`, as well as of the temperature of the outlet stream, :math:`T(t)`, based on the following assumptions:

 	*	The concentrations of :math:`{\sf A}` in the inlet stream is :math:`C_{\sf A}^{\rm in}`, the inlet and outlet volumetric flow rates are constant and both equal to :math:`F`, and the temperature of the inlet stream is :math:`T^{\rm in}`.
	*	The reaction has first-order kinetics, with the dependence of the reaction rate :math:`r` (moles of :math:`{\sf A}` consumed per unit volume and per unit time) on the temperature given by :math:`r=k_0\, e^{-\frac{E}{RT}}\, C_{\sf A}`.
			
		where :math:`E` is the activation energy; :math:`R`, the universal gas constant; :math:`T`, the temperature at which the reaction is taking place; and :math:`k_0`, a constant.
		
	*	The heat of reaction is :math:`\Delta H` and no heat is lost or added to the reactor, other than through the liquid streams.
	
	*	The density and specific heat capacity of the fluid are :math:`\rho` and :math:`C_p` respectively.

#.	For the values of the process parameters listed in the table and the following initial conditions :math:`C_{\sf A}(0) = 0.7\ \rm kmol/m^3` and :math:`T(0)=355\ \rm K`, compute, by hand, the value of :math:`\big[C_{\sf A}(t),\; T(t)\big]` at :math:`t=0.2`, using:

	*	Euler's method with :math:`h=0.1`, and
	*	Heun's predictor-corrector method with :math:`h = 0.1`, and
	*	the fourth-order Runge-Kutta method with :math:`h=0.2`
	
#.	Implement these equations in computer software and compare the integration by hand with the computer version at :math:`t=0.2`.

.. Solution: see code/cstr_integration.py

.. raw:: latex

	\vspace{0.25cm}
	\hrule
	\begin{center}END\end{center}

</rst>