View source for Tutorial 4 - 2010 - Solution

{{OtherSidebar
| due_dates = 7 October 2010
| dates_alt_text = Due date(s)
| questions_PDF = Tutorial-4-2010.pdf
| questions_text_alt = Tutorial questions
| solutions_PDF = Tutorial-4-2010-Solution.pdf
| solutions_text_alt = Solutions by Ali Sahlodin
| other_instructions = Hand-in at class.
}}
<rst>
<rst-options: 'toc' = False/>
<rst-options: 'reset-figures' = False/>


.. |m3| replace:: m\ :sup:`3`
.. highlight:: python

.. rubric:: Tutorial objectives

* Derive linear models for actual systems.
* Solve these models by hand, and with computer software.

.. rubric::  Solutions prepared by Ali Sahlodin.

Question 1 [1]
===============

In the second tutorial you derived the set of linear equations for the reaction:

.. math::
	{\sf P_2I_4} + n\,{\sf P_4} + p\,{\sf H_2O} \longrightarrow 4\,{\sf PH_4I} + q\,{\sf H_3PO_4} 

where :math:`n`, :math:`p` and :math:`q` denote the stoichiometric coefficients for the :math:`\sf P_4`, :math:`\sf H_2O` and :math:`\sf H_3PO_4` species respectively.

#.	Solve the linear equations using Gauss elimination, by hand.
#.	Verify your solution using the ``mldivide`` (MATLAB) or ``solve`` (Python) function.  How long does it take to find the solution?  Contrast this with the approach used in tutorial 2.

Solution
-----------

Recall from the second tutorial that we have the following balance equations;

	balance for :math:`{\sf O}:  p- 4q=0`,

	balance for :math:`{\sf P}:  4n - q=2`,

	balance for :math:`{\sf H}:  2p-3q = 16`.


These equations can be written in the form of :math:`Ax=b` where :math:`A` is the matrix of coefficients, and :math:`b` is a column vector containing the right-hand side constants. Also :math:`x=[n\quad p\quad q]^T` is the vector of unknowns. So, we have:

	.. math::
		\left[ \begin{array}{ccc}			0 & 1 & -4 \\			4 & 0 & -1 \\			0 & 2 & -3 		\end{array} \right] \left[ \begin{array}{c} 			x_{1} \\ x_{2} \\ x_{3} 		\end{array} \right] = \left[ \begin{array}{c} 			0 \\ 2 \\ 16		\end{array} \right]

Now we solve this system using the Gauss elimination (:math:`R_i` denotes the :math:`i-th` row of the matrix). Note that your solution procedure (not the final results) may be different if you have arranged the equations and unknowns differently so you have a different :math:`A` and :math:`b` than the above.

**Forward Elimination**

As a first step, :math:`A_{2,1}` must be eliminated using :math:`A_{1,1}` in the first row. However, :math:`A_{1,1}=0` by which we cannot divide. So, we do a partial pivoting by switching the first and second rows in the above matrix. That is:

	.. math::
		\left[ \begin{array}{ccc}			4 & 0 & -1 \\			0 & 1 & -4 \\			0 & 2 & -3 		\end{array} \right] \left[ \begin{array}{c} 			x_{1} \\ x_{2} \\ x_{3} 		\end{array} \right] = \left[ \begin{array}{c} 			2 \\ 0 \\ 16		\end{array} \right]

Interestingly, :math:`A_{2,1}` is already zero with the new rearrangement. This is a lucky case since nothing needs to be done to the second row of the matrix. Now we eliminate :math:`A_{3,2}` using :math:`A_{2,2}`. To do so, multiply :math:`R_2` by :math:`-2`, and add the result to :math:`R_3`:

	:math:`R_3 \leftarrow R_3-2R_2`

which results in

	.. math::
		\left[ \begin{array}{ccc}			4 & 0 & -1 \\			0 & 1 & -4 \\			0 & 0 & 5 		\end{array} \right] \left[ \begin{array}{c} 			x_{1} \\ x_{2} \\ x_{3} 		\end{array} \right] = \left[ \begin{array}{c} 			2 \\ 0 \\ 16		\end{array} \right]

The forward elimination is complete.

**Backward substitution**

From :math:`R_3`: :math:`5x_3=16\rightarrow x_3=q=3.2`. Substituting :math:`x_3` into :math:`R_2`: :math:`x_2-4\times 3.2=0\rightarrow x_2=p=12.8`.
Finally, substituting :math:`x_2` and :math:`x_3` into :math:`R_1`: :math:`4x_1-3.2=2\rightarrow x_1=n=1.3`.

*Solution by software*

	.. twocolumncode::
	    :code1: ../che3e4/Tutorials/Tutorial-4/code/question1.m
	    :language1: matlab
	    :header1: MATLAB code
	    :code2: ../che3e4/Tutorials/Tutorial-4/code/question1.py
	    :language2: python
	    :header2: Python code

The time required to solve the above system using built-in :math:`Ax=b` solvers is much less than that required by the brute-force approach in Tutorial 2.

Question 2 [1.5]
================

When it comes down to the basics, any recipe can be reduced to 4 components: fat, carbohydrates, protein and moisture (usually water).  What differentiates each recipe though is the list of ingredients which are combined, and the ratio of fat:carbohydrate:protein:water.  Of course how the ingredients are mixed and the reactions that take place during cooking are also important, but that affects mainly taste and texture, not nutritional value.

Let's consider a biscuit recipe for this question.  One can use butter for the fat component, use brown sugar and chocolate for the carbohydrate component (chocolate also contributes to the fat component), flour adds to the carbohydrate component also and somewhat to the protein component, while eggs, milk or water make up the rest of the recipe for protein and moisture.

We would like to blend raw materials where the ratio of fat:carbohydrate:protein:water of the uncooked ingredients is 20:50:10:20.  We have the following ingredients available, broken down on a mass percentage basis.

============= ====== ==== ===== ===== ========== === =========
Component     Butter Lard Flour Sugar Whole milk Egg Chocolate
============= ====== ==== ===== ===== ========== === =========
Fat           85     81   0     0     4          10  14
Carbohydrate  0      1    86    98    6          1   75
Protein       1      4    10    0     4          35  6 
Moisture      14     14   4     2     86         54  5 
============= ====== ==== ===== ===== ========== === =========

#.	Write down the 4 mass balance equations, one for each component, in the form :math:`Ax = b`, where :math:`x` is vector that represents the mass of each ingredient to be used, a vector with 7 rows and one column.
#.	What is the technical (mathematical) reason why we cannot solve this system of equations using Gauss elimination?  Translate your mathematical answer to English, explaining it in terms of the recipe.
#.	If we are constrained to using only butter, flour, whole milk and eggs: solve, using LU decomposition on the computer, for the amounts of each to be used, in grams, to obtain 100 grams of uncooked ingredients, in the required ratio.

If you prefer, you can interpret this question in the context of rubber manufacturing, where various oils, polypropylene and existing rubber materials are blended to create a new rubber with desired physical properties.

Solution
----------

#.	For a 100 units of the recipe, we should have 20 units of fat, 50 units of carbohydrate, 10 units of protein, and 20 units of moisture. Given the percentage of these components in each of the ingredients, the four mass balance equations can be written as: 

	==============   ===================================================================================
	Fat              :math:`85x_1+     81x_2+   0x_3 +    0x_4+     4x_5+  10x_6+  14x_7=20`
	Carbohydrate     :math:`0x_1+      1x_2+    86x_3+    98x_4+    6x_5+  1x_6+   75x_7=50`
	Protein          :math:`1x_1+      4x_2+    10x_3+    0x_4+     4x_5+  35x_6+  6x_7=10`
	Moisture         :math:`14x_1+     14x_2+   4x_3+     2x_4+     86x_5+ 54x_6+  5x_7=20`
	==============   ===================================================================================

	The above equations can be written in form :math:`Ax = b` as

	.. math::
			\left[ \begin{array}{ccccccc}				85&81&0&0&4&10&14 \\				0&1&86&98&6&1&75 \\				1&4&10&0&4&35&6\\				14 & 14 & 4 & 2 & 86 & 54 & 5			\end{array} \right] \left[ \begin{array}{c} 				x_{1} \\ x_{2} \\ x_{3} \\ x_{4}\\x_5\\x_6\\x_7			\end{array} \right] = \left[ \begin{array}{c} 				20 \\ 50 \\ 10\\ 20			\end{array} \right]

	where :math:`x_i` represent the mass of the 7 ingredients from butter to chocolate, respectively.

#.	The above system cannot be solved for :math:`x` because the number of unknowns is **not** equal to the number of equations. In other words, the degrees of freedom is not zero, and we must have either more equations or less unknowns to be able to solve this system. In terms of the recipe, there are infinite combinations of the ingredients that yield the desired component ratio. So, there is no unique recipe given only the above information.

#.	If we limit ourselves to the use of only butter, flour, whole milk and eggs, we imply that the mass of the other three ingredients is zero in our recipe. So, we have eliminated three unknowns by fixing them to zero. Now, we have four unknowns in four equations. This system is now solvable. The new system will only have the coefficients and unknowns corresponding to butter, flour, whole milk and eggs as:

	.. math::
		\left[ \begin{array}{ccccccc}			85 & 0  & 4  & 10 \\			0  & 86 & 6  & 1 \\			1  & 10 & 4  & 35\\			14 & 4  & 86 & 54		\end{array} \right] \left[ \begin{array}{c} 			x_{1} \\ x_{3} \\x_5\\x_6		\end{array} \right] = \left[ \begin{array}{c} 			20 \\ 50 \\ 10\\ 20		\end{array} \right]


	Recall from class that :math:`A = LU`, but with a permutation or pivot matrix, it becomes :math:`PA = LU`, or alternatively: :math:`A = P^{-1}LU`

	:math:`Ax = b\rightarrow P^{-1}LUx = b`.  Let us :math:`y=Ux`.

	Then solve :math:`(P^{-1}L)y = b` for :math:`y`. Finally, solve :math:`Ux = y` for :math:`x`. It can be implemented as follows:

	.. twocolumncode::
	    :code1: ../che3e4/Tutorials/Tutorial-4/code/question2.m
	    :language1: matlab
	    :header1: MATLAB code
	    :code2: ../che3e4/Tutorials/Tutorial-4/code/question2.py
	    :language2: python
	    :header2: Python code

.. A = np.array([[85,81,0,0,4,10,14],[0,1,81,98,6,1,75],[1,4,10,0,4,35,6],[14,14,9,1,86,54,5]])
.. A = A[:,[0,2,4,5]]
.. b = np.array([20, 50, 10, 20])
.. linalg.solve(A[:,[0,2,4,5]],b)*100


Question 3 [1.5]
================

In class we derived the molar balance for a species in a similar flowsheet.  Use the slightly modified flowsheet below

.. figure:: images/reactor-network.png
	:scale: 30
	:align: center

and let :math:`\alpha` be the proportion of :math:`m_3` that is sent to reactor RB, and the rest, :math:`(1-\alpha)m_3`, is sent to reactor RC.  Similarly, let  :math:`\beta` represent the proportion of :math:`m_8` that leaves the system, and let :math:`(1-\beta)m_8` be the amount sent to the recycle stream.  Also let :math:`X_A` be the conversion of the species in reactor RA, and similarly for :math:`X_B` and :math:`X_C`.

#.	Write a MATLAB or Python function that will return a matrix :math:`A` given the five inputs, :math:`\alpha, \beta, X_A, X_B`, and :math:`X_C`.  This matrix :math:`A` represents the coefficients in the molar balances on the flowsheet, and these balances can be solved using :math:`Ax = b`.
#.	Given that :math:`m_1` = 10 mol/second, and that the molar conversion of the inlet stream to reactor A is 5%, i.e. :math:`X_A = 0.05`, :math:`X_B = 0.60` and :math:`X_C = 0.85`.  What are the molar flow rates, :math:`m_6, m_7, m_9`, and :math:`m_{10}` when:

	*	:math:`\alpha=0.2`, :math:`\beta=0.7` 
	*	:math:`\alpha=0.2`, :math:`\beta=0.4`
	
	Are your answers reasonable?
		
	You should use either MATLAB or Python's built-in functions  to solve the system of equations.
	
Solution
----------

#.	The feed stream, :math:`m_1`, is the molar flow rate of the species of interest. A molar balance can be written around each reactor, as well as split and each junction point:


        :math:`m_1 = \text{some given value}` mol/s

        :math:`m_2 = m_1 + m_9`

        :math:`m_3 = (1-X_A)m_2`

        :math:`m_4 = \alpha m_3`

        :math:`m_5 = (1-\alpha) m_3`

        :math:`m_6 = (1-X_B) m_4`

        :math:`m_7 = (1-X_C) m_5`

        :math:`m_8 = m_6 + m_7`

        :math:`m_9 = (1-\beta) m_8`

        :math:`m_{10} = \beta m_8`

Rearranging the above equations to have all unknowns on one side, and all known values on the other side gives:

        :math:`m_1 = \text{some given value}` 

        :math:`-m_1 + m_2 - m_9=0`

        :math:`-(1-X_A)m_2 + m_3=0`

        :math:`- \alpha m_3+ m_4=0`

        :math:`- (1-\alpha) m_3 + m_5 =0`

        :math:`- (1-X_B) m_4+ m_6 =0`

        :math:`- (1-X_C) m_5 + m_7 =0`

        :math:`-m_6 - m_7+ m_8 = 0`

        :math:`- (1-\beta) m_8+ m_9 =0`

        :math:`-\beta m_8 + m_{10} = 0`

If we have the values of the parameters :math:`\alpha, \beta, X_A, X_B`, and :math:`X_C`, then we can write a function that accepts these input arguments and whose output is the matrix :math:`A`. This function, and the code that calls it later on can be implemented as shown below.  With the given parameter values, the molar flow rates can be solved easily. First we have to evaluate the matrix :math:`A` by calling the function ``matrix_A``, and providing the values of the parameters. Then, the system is solved for :math:`m_i`. 

.. twocolumncode::
    :code1: ../che3e4/Tutorials/Tutorial-4/code/matrix_A.m
    :language1: matlab
    :header1: MATLAB code
    :code2: ../che3e4/Tutorials/Tutorial-4/code/matrix_A.py
    :language2: python
    :header2: Python code

The solution for the two scenarios above differ only in terms of the :math:`\beta` value.

* :math:`\beta = 0.7` is :math:`x = [ 10.0, 10.6, 10.1, 2.01, 8.06, 0.806, 1.21, 2.01, 0.604, 1.41]`
* :math:`\beta = 0.4` is :math:`x = [10.0,  11.3, 10.7, 2.14, 8.58, 0.858, 1.29, 2.14, 1.29, 0.859]`

It shows that as :math:`\beta` decreases (i.e. we recycle more of our stream :math:`m_8`) that the amount of the species leaving the overall reactor decreases from 1.41 to 0.859.  This is expected, because we know that as we increase the recycle rate, that a greater amount of unconverted material will be reacted.  So even this relatively simplistic model makes engineering sense.

Use the above code and see what happens when you let :math:`\beta = 0.0` or :math:`\beta = 1.0`: do the results still make sense?

.. raw:: latex

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

</rst>