Tutorial 2 - 2010 - Solution

Due date(s):	22 September 2010
(PDF)	Tutorial questions (updated)
(PDF)	Solutions, prepared 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

- Some questions to help you feel comfortable deriving model equations for actual chemical engineering systems. - Brute force solving of equation systems. The rest of the course will focus on better ways to solve these equations. - Interpreting source code written on paper.

.. rubric:: Recap of tutorial rules

• Tutorials can be done in groups of two - please take advantage of this to learn with each other.
• Tutorials must be handed in at the start of class on Wednesday. No electronic submissions - thanks!

Question 1 [2]

===

.. note::

Parts 1 and 2 of this question were from the 2006 final exam (slightly modified). It was worth 10% of the 3 hour exam

Consider a mixing tank with fluid volume :math:`V` where fluids A and B with densities :math:`\rho_A` and :math:`\rho_b`, specific heat capacities :math:`C_{p,A}`, :math:`C_{p,B}` and temperatures :math:`T_A` and :math:`T_B` are the inlet streams at volumetric flow rates :math:`F_A` and :math:`F_B`. The outlet stream is at a flow rate :math:`F_A+F_B`. The specific heat capacity and density of the outlet streams is given by :math:`C_p=(C_{p,A}F_A+C_{p,B}F_B)/F` and :math:`\rho=(\rho_AF_A+\rho_BF_B)/F`. The fluid also looses heat from the tank at a rate :math:`q=k_T(T-T_\text{wall})` where :math:`T_\text{wall}` is the constant tank wall temperature, :math:`k_T` is a constant and :math:`T` denotes the current fluid temperature.

1. . Using 3-step modelling approach shown in class, derive a dynamical balance describing the time-dependent exit stream temperature.
2. . Can the steady state exit stream temperature be higher than both :math:`T_A` and :math:`T_B`? Explain.
3. . Calculate, *by-hand*, the steady-state exit temperature, using that

* :math:`V` = 10 |m3| * :math:`\rho_A` = 1200 kg/\ |m3| and :math:`\rho_b` = 950 kg/\ |m3| * :math:`C_{p,A}` = 2440 J/(kg.K) and :math:`C_{p,B}` = 3950 J/(kg.K) * :math:`T_A` = 320 K and :math:`T_B` = 350 K and :math:`T_\text{wall}` = 300K * :math:`F_A = F_B` = 0.05 |m3|/s * :math:`k_T` = 200 W/(m\ :sup:`2`.K) :math:`\times` 24 m\ :sup:`2` = 4800 W/K

Solution

---

• • Step 1: Definition of the problem**

• What are the inputs and outputs?*

In terms of material inputs/output, the process has two inputs as the flow of streams A and B. It also has one output which is a mixture of those two streams.

In terms of energy input/output, we have energy coming in and going out from the material streams as well as energy transfer through the shaft and the wall.

• Lumped vs. distributed, steady-state vs. dynamic?*

Since we have an impeller inside the tank, we can assume that we have a well-mixed environment which is considered a lumped system. That is, there is NO spatial distribution of quantities through the system. As for change of quantities with time, the flow rate is at steady state as the out-flow equals the sum of the in-flows. For the temperature, let us go with the dynamic formulation for the moment as requested by the problem (we will make the steady-state assumption when we want to solve for the steady-state temperature.)

• • Step 2: Controlling mechanisms:**

This is a mixing tank problem, where there is no reaction nor mass diffusion. However, there is heat transfer from the wall. Note that the work input the by the impeller is safely ignored. Also, radiation is negligible at normal temperatures.

Also, we neglect the contributions from the kinetic and potential energies in the energy balance.

• • Step 3: Development of a set of model equations**

Let us denote by the subscript :math:`m` the conditions of the mixed fluid that goes out of the tank. The dynamic heat balance is written as follows.

.. math::

\frac{d(\rho_m V C_{p,m} T(t))}{dt} = \rho_A F_A C_{p,A}(T_A - T(t)) + \rho_B F_B C_{p,B}(T_B - T(t)) - k_T(T(t) - T_{wall})

where :math:`F=F_A+F_B`.

• • Remarks and solutions to other questions**

• I am confused as to the sign of* :math:`k_T(T - T_{wall})` *in the heat balance equation. Should I put a negative before it or not?*

OK. First you do not need to know whether the heat is lost from or gained into the system to determine the sign of :math:`k_T(T-T_{wall})`. Remember that the heat always goes from the higher temperature to the lower one. So, you cannot even tell whether heat gets out of or comes into the system without knowing the value of :math:`T`, which is indeed your unknown. But you have to put a negative sign in this case even if you had not been told about the heat loss. Here is why:

If :math:`T>T_{wall} \rightarrow T-T_{wall}>0`. In this case there is heat loss, which is correctly modelled as :math:`-k_T(T-T_{wall})<0` (negative heat contribution).

If :math:`T<T_{wall} \rightarrow T-T_{wall}<0`. In this case there is heat gain, which is again correctly modelled using :math:`-k_T(T-T_{wall})>0` (positive heat contribution)

• Simplifications:*

The liquid volume remains constant as the input and output flows are equal. Also, this is a liquid system where the density can be considered constant. Also, the specific heat capacity can be considered constant if there is no severe variations in T (for gaseous system we would have to consider a variable specific heat capacity). So, the volume :math:`V`, density :math:`\rho_m`, and specific heat capacity :math:`C_{p,m}` can be taken out of the derivative:

.. math::

\rho_m V C_{p,m} \frac{dT(t)}{dt} = \rho_A F_A C_{p,A}(T_A - T(t)) + \rho_B F_B C_{p,B}(T_B - T(t)) - k_T(T(t) - T_{wall})

• Can the steady state exit stream temperature be higher than both* :math:`T_A` and :math:`T_B`? *Explain.*

In general, it depends. There is no heat generation or consumption (e.g. via reaction) inside the tank. If there was no heat conduction through the wall (adiabatic tank), the steady-state temperature would be between :math:`T_A` and :math:`T_B`, (:math:`T_B>T_A`). However, with the heat conduction, the steady-state temperature depends on the value of :math:`T_{wall}`, and other parameters in the process (e.g. :math:`k_T`). In this case where :math:`T_{wall}<T_A<T_B`, we can say that the steady-state temperature will definitely be less than :math:`T_B`. But, we cannot decide on whether the steady-state temperature is greater than :math:`T_A` or not before solving the model for :math:`T`.

• Calculate, by-hand, the steady-state exit temperature*

Recall the heat balance equation:

.. math::

\rho_m V C_{p,m} \frac{dT(t)}{dt} = \rho_A F_A C_{p,A}(T_A - T(t)) + \rho_B F_B C_{p,B}(T_B - T(t)) - k_T(T(t) - T_{wall})

The term describing the time derivative becomes zero in steady state, so we are left with:

:math:`\rho_A F_A C_{p,A}(T_A - \overline{T}) + \rho_B F_B C_{p,B}(T_B - \overline{T}) - k_T(\overline{T} - T_{wall}) = 0`

By substituting the values given in the problem, and solving for :math:`\overline{T}`, the steady-state temperature is obtained as: :math:`\overline{T}\bf = 336.3` K.

.. note:: A MATLAB code for solving the above equation - not required for full grade. Please read about the ``fsolve`` for details. We will cover this in the 3rd part of the course.

.. code-block:: matlab

% Define the parameters and variables rho_A = 1200; rho_b = 950; C_pA = 2440; C_pB = 3950; T_A = 320; T_B = 350; T_w = 300; F_A = 0.05; F_B = 0.05; k_T = 4800;

% Define the function to be solved, in the form: f(x)=0 LHS = @(T) rho_A*C_pA*F_A*(T_A - T) + rho_b*C_pB*F_B*(T_B - T) - k_T*(T - T_w);

% Solve the function for T; Initial guess= 300 K [T, LHS_val] = fsolve(@(T) LHS(T),300) %LHS_val should be zero at %the solution % Solution % -------- % T = 336.3292, though when you report your answer, you should % use at most 3 significant figures. % LHS_val = -9.1386e-09

Question 2 [2]

=====

.. note:: You don't need to write any code for this question.

Consider 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 :math:`\sf P_4`, :math:`\sf H_2O` and :math:`\sf H_3PO_4` respectively.

1. . Derive the equations necessary to solve for :math:`n, p`, and :math:`q` by equating atoms of :math:`{\sf P}`, :math:`{\sf H}`, and :math:`{\sf O}` on the reactant and product sides.

1. . In the next section of the course we will use Gauss Elimination to solve these equations. For now though, let's describe a brute force approach. First, complete the two lines of this MATLAB function:

.. code-block:: matlab

function total_error = equation_error( n, p, q )

% Given the values of n, p, and q, calculate the error of each balance equation. % Returns the sum of squares of the errors.

error_1 = _______________  % from the P-balance error_2 = 2*p - 3*q - 16;  % from the H-balance error_3 = _______________  % from the O-balance

total_error = (error_1)^2 + (error_2)^2 + (error_3)^2;

end  % end of function

or complete this Python function:

.. code-block:: python

def equation_error(n, p, q ): """ Given the values of n, p, and q, calculate the error of each of the 3 equations. Returns the sum of squares of the errors. """ error_1 = _______________ error_2 = 2*p - 3*q - 16 error_3 = _______________

return error_1**2 + error_2**2 + error_3**2

1. . Since we known that :math:`n, p`, and :math:`q` must be positive, we can construct a set of 3 nested for-loops, as shown below in MATLAB and Python. Describe in plain English what the code does.

**In MATLAB:**

.. code-block:: matlab

smallest = 0.0; largest = 14.9; step_size = 0.1; vector = smallest : step_size : largest;

% How many elements in each vector? num = length(vector); errors = zeros(num, num, num);

index_n = 0; index_p = 0; index_q = 0; for n = smallest : step_size : largest index_n = index_n + 1; for p = smallest : step_size : largest index_p = index_p + 1; for q = smallest : step_size : largest index_q = index_q + 1;

 % Calculate the error at this value of n, p and q: errors(index_n, index_p, index_q) = equation_error(n, p, q);

end index_q = 0; end index_p = 0; end [min_error, min_index] = min(errors(:)) [index_n, index_p, index_q] = ind2sub([num, num, num], min_index); disp(['Solution at ', num2str([vector(index_n), vector(index_p), vector(index_q)])])

**In Python:**

.. code-block:: python

import numpy as np

smallest = 0.0 largest = 15.0 step_size = 0.1 vector = np.arange(smallest, largest, step_size)

# How many elements in each vector? num = len(vector)

errors = np.zeros( (num, num, num) )

for index_n, n in enumerate(vector): for index_p, p in enumerate(vector): for index_q, q in enumerate(vector):

# Calculate the error at this value of n, p and q: errors[index_n, index_p, index_q] = equation_error(n, p, q)

# Which combination had the smallest error? min_index = np.argmin(errors) index_n, index_p, index_q = np.unravel_index(min_index, (num, num, num)) n, p, q = vector[index_n], vector[index_p], vector[index_q] print(n, p, q)

1. . How many times will the function ``equation_error`` be called?

1. . What will this function output be if :math:`(n, p, q) = (1.0, 9.2, 2.5)`?

Solution

---

1. . In the above reaction, component I is already balanced on both sides. Therefore, it remains to write the balance equations for components P, H, and O.

Balance for :math:`{\sf P}: 2 + 4n = q`, balance for :math:`{\sf H}: 2p = 16 + 3q`, balance for :math:`{\sf O}: p = 4q`.

This gives rise to a system of three equations with three unknowns which can be solved easily using tools from linear algebra. It was not required to solve for the unknowns.

1. . The equations can be rearranged and set to set to zero. Each of the balance equations ``error_1``, ``error_2``, and ``error_3`` must be zero in order to obtain the correct :math:`n`, :math:`p`, and :math:`q`.

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

1. . This code will find the minimum of the sum of the squared errors. The errors are the amount by which the balance equations are incorrect given the values of :math:`n`, :math:`p`, and :math:`q`. If an exact solution is found, then ``total_error`` would be zero.

*A more detailed explanation*: The code tries different values of :math:`n`, :math:`p`, and :math:`q` over a specified range (using the for loops), and evaluates the ``total_error`` for each set of values through the function ``equation_error``. The resulting values of the ``total_error`` are stored in a 3D array called ``errors``. The ``index_n``, ``index_p``, and ``index_q`` determine the location of an entry in the 3D array, and correspond to the values of :math:`n`, :math:`p`, and :math:`q` at which the ``total_error`` has been evaluated.

At the end of the for loops, we find the minimum value in the resulting array ``errors``, which means that we find the minimum ``total_error``. Then, using matrix tools we pull out the indices where this minimum ``total_error`` has been stored, and subsequently, we find the corresponding :math:`n`, :math:`p`, and :math:`q` that have minimized our ``total_error``. As mentioned earlier, the ``total_error`` should ideally be zero. Note that the individual errors error_1, error_2, and error_3 are squared because negative and positive values are both deviation and should not cancel out each other.

1. . The equation_error function is called as many times as the most inner loop (``q``-loop) is triggered. The ``q``-loop by itself is run from ``q=0`` to ``q=14.9`` with step size of 0.1. Therefore, it is run 1+(14.9-0.0)/0.1=150 times. However, the ``q``-loop itself is inside an outer loop, the ``p``-loop, which is also run 150 times. These two loops themselves are inside the ``n``-loop that again is run 150 times. Therefore, the total number of calls to the equation_error function will be ``150*150*150`` = **3375000** times.

You can have your code return the total number of calls quite simply as seen in this solution code given in the :ref:`bonus question <bonus-question>`.

1. . To obtain the total_error at :math:`(n, p, q) = (1.0, 9.2, 2.5)`, we can simply call the function ``equation_error`` with input arguments being the values of :math:`n`, :math:`p`, and :math:`q`. In both MATLAB and Python, the function output is **26.9**.

.. _bonus-question:

Bonus question [0.5]

=========

Using the code given in question 2, report what the ``min_index`` variable is and what are the values of :math:`n, p`, and :math:`q` which give minimum error to the set of equations. How long did it take to find the solution of this simple linear equation system?

Solution

---

• • MATLAB**

The ``min_index`` is the index of the entry in the array ``errors`` at which the minimum ``total_error`` has been stored. Note that ``errors`` is converted to a column vector through ``errors(:)``. Therefore, this ``min_index`` returns the index in the resulting vector. Its value is ``min_index = 739214``. In order to locate the corresponding values of :math:`n`, :math:`p`, and :math:`q`, we need to convert this vector-wise index to an array-wise subscript. Recall that we have a 3D array with each dimension corresponding to one of the unknowns. So, if we convert the ``min_index`` to a subscript form of ``(index_n, index_p, index_q)``, we can then find the values of :math:`n`, :math:`p`, and :math:`q` that correspond to the ``min_error``. This is done through the ``ind2sub`` command. The answer will be: :math:`\textbf{n=1.3, p=12.8, q=3.2}`.

The time taken to do the 3 nested loops can be easily found from the ``tic`` and ``toc`` commands. The solution code appended below shows how it is used.

Note the discrepancy in timing between MATLAB and Python. Recent MATLAB versions have technology built-in to accelerate code by vectorization. Transparent to the user, they will rewrite your MATLAB code to speed it up. If you turn this acceleration off, by typing ``feature accel off`` before you run the script, the time taken is roughly 94 seconds instead of 3.8 seconds.

• • Python**

The Python code works exactly the same way, except using a different function, ``np.unravel_index`` to locate the smallest entry in the 3D array. The time taken on my machine was around 54 seconds. It can be reduced by using vectorized calculations.

.. twocolumncode::

   :code1: ../che3e4/Tutorials/Tutorial-2/code/question_2.m
   :language1: matlab
   :header1: MATLAB code
   :code2: ../che3e4/Tutorials/Tutorial-2/code/question-2.py
   :language2: python
   :header2: Python code

.. raw:: latex

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

</rst>