TMA681 Tillämpad matematik K, 2001

STUDIO 2. THE TANK REACTOR: HEAT BALANCE.

Important: write a readable report of your work in each studio session. You will need this at the end of the course when you prepare for the written exam; several exam questions will be based directly on the studio work which is only documented in your own notes and the instructions.

Some exercises are called homework in order to save time in the studio classroom.

1. Systems of differential equations

Recall that \frac{d}{dt} cos(t) = − sin(t), \frac{d}{dt} sin(t) = cos(t), cos(0) = 1, and sin(0) = 0. This means that the trigonometric functions u_{1}(t) = cos(t) and u_{2}(t) = sin(t) satisfy the initial-value problem

u′_{1}(t) = −u_{2}(t), t > 0,
(1)  u′_{2}(t) = u_{1}(t), t > 0,
u_{1}(0) = 1, u_{2}(0) = 0.

This is a system of (linear) differential equations of first order. Exercise 1. Solve (1) with Matlab. Hint: Write the following function file trig.m, which defines the right-hand side of the system of differential equations.

  function y=trig(t,u)
  y = zeros(2,1); % trig(t,u) must return a column vector
  y(1) = -u(2);
  y(2) = u(1);

Then write the following script file trigdata.m, which gives default values to the variables that are used.

  T = 2*pi;
  u0 = [1; 0];

Finally write the following script file solve2.m, which solves (1) and plots the solution.

  [t,u] = ode45('trig', [0 T], u0);
  subplot(2,1,1) % breaks the figure into a 2-by-1 matrix; selects top half
  plot(t, u(:,1)); % plots first component of u versus t
  hold on
  plot(t, u(:,2), '--'); % dashed line
  hold off
  title('Solid: u_1   Dashed: u_2'), xlabel('t')
  subplot(2,1,2) % selects bottom half
  plot(u(:,1), u(:,2)) % plots second versus first component of u
  title('Phase portrait'), xlabel('u_1'), ylabel('u_2')

Start the computation by typing the following on the Matlab command line:

>> trigdata
>> solve2
>> T=pi
>> solve2

Compute and observe the solutions for various values of T. Consider in particular the phase plot of u_{2} versus u_{1}, where t can be viewed as a parameter along the curve. What point on the curve corresponds to t = 0? t = \frac{π}{2} ? t = 2π?

The hyperbolic functions, defined by cosh(t) = \frac{exp(t)+exp(−t)}{2} and sinh(t) = \frac{exp(t)−exp(−t)}{2}, have the properties \frac{d}{dt} cosh(t) = sinh(t), \frac{d}{dt} sinh(t) = cosh(t), cosh(0) = 1, and sinh(0) = 0. Homework 1. Verify these properties.

Date: May 7, 2001, Niklas Ericsson, Stig Larsson, Mathematics, Chalmers University of Technology.

Exercise 2. Set u_{1}(t) = cosh(t) and u_{2}(t) = sinh(t). Derive an initial-value problem for u_{1}(t) and u_{2}(t) similar to (1). Compute the hyperbolic functions by solving this problem with Matlab. Hint: You need only slightly modify the trigonometric case above.

2. The tank reactor

Last week we considered the simplified case where we assumed that the rate coefficient k is constant, i.e., we neglected the temperature dependence of the reaction rate. Today we will include this dependence (the Arrhenius law) to improve our mathematical model of the tank reactor. Therefore, we need to consider also the temperature, T, in the reactor, and start by deriving an equation for it by considering heat balance in the tank reactor. In doing so we will need to know the temperature, T_{K}, in the cooler tank. Here we will once again make a simplifying assumption, namely, that T_{K} is constant. Physically, this can be expected to be the case when the cooler tank is large compared to the reactor tank.

2.1. Heat balance. The heat balance equation is

(2)  ρc_{p}V \frac{dT}{dt} = ρc_{p}q(T_{f} − T) + (−δH)V ck_{0} exp(−E/(RT)) − κA_{K}(T − T_{K}). [J/s]

Recall that ρ and c_{p} are the density and the heat capacity of the fluid in the reactor, respectively. Since V is the volume of the tank reactor, the left-hand side expresses the rate of change of the internal energy of the fluid in the reactor. Three different mechanisms contribute to this:

• The first term on the right-hand side expresses the net rate of change of the energy due to the temperature difference, T_{f} − T, between the inflow and the outflow (q is the flux through the tank).

• The second term on the right-hand side expresses the increase rate of the energy due to heat released in the (exothermal) reaction. Recall that δH is the heat of reaction, c is the concentration in the reactor of the reacting substance, and k = k_{0} exp(−E/(RT)) is the rate coefficient of the reaction. Note: since the reaction rate, ck, depends both on the concentration and on the temperature, there is a (non-linear) coupling between the equation for the concentration, that we studied last week, and the equation for the temperature.

• The third term on the right-hand side expresses the decrease rate of the energy due to heat transferred from the reactor to the cooler. This rate is proportional to the temperature difference, T − T_{K}, between the reactor and the cooler. Recall that κ and A_{K} are the heat transfer coefficient and the area of the cooler, respectively.

Note that (2) is a non-linear differential equation for the temperature, T, in the tank reactor. In order to make the equation dimensionless we divide by ρc_{p}q_{ref}T_{f} [J/s]. We get

(3) \frac{V}{q_{ref}} \frac{d}{dt}(\frac{T}{T_{f}} ) = \frac{q}{q_{ref}} (1 − \frac{T}{T_{f}} ) + \frac{(−δH)c_{f}}{ρc_{p}T_{f}} \frac{c}{c_{f}} \frac{V}{q_{ref}} k_{0} exp (−\frac{E}{RT_{f}} \frac{T_{f}}{T} ) − \frac{κA_{K}}{ρc_{p}q_{ref}} (\frac{T}{T_{f}} − \frac{T_{K}}{T_{f}} ).

Recall that last week we introduced the dimensionless time s = t/τ, where τ = V/q_{ref} [s]. We now also introduce

X_{1} = c/c_{f} (dimensionless concentration),
X_{2} = T/T_{f} (dimensionless reactor temperature),
U_{1} = q/q_{ref} (dimensionless flux),
U_{2} = T_{K}/T_{f} (dimensionless cooler temperature).

Note that X_{1} is the quantity last week referred to as X, and U_{1} is the quantity last week referred to as U. To express the left-hand side in (3) in the non-dimensional variables we note that, by the chain rule,

\frac{dX_{2}}{ds} = \frac{dX_{2}}{dt} \frac{dt}{ds} = τ \frac{dX_{2}}{dt} = \frac{V}{q_{ref}} \frac{d}{dt}(\frac{T}{T_{f}} ). (4)

We also express the terms on the right-hand side in (3) in the non-dimensional variables:

• The first term on the right-hand side simply becomes

U_{1}(1 − X_{2}). (5)

• For the second term on the right-hand side we introduce the three dimensionless numbers

α = \frac{(−δH)c_{f}}{ρc_{p}T_{f}} , γ = \frac{E}{RT_{f}} , δ = k_{0}τ exp(−γ),

where γ is referred to as the Arrhenius number. This term then becomes

(6)  αX_{1}τk_{0} exp (−\frac{γ}{X_{2}} ) = αX_{1}k_{0}τ exp(−γ) exp(γ−\frac{γ}{X_{2}} ) = αX_{1}f(X_{2}),

where

f(X_{2}) = δ exp(γ−\frac{γ}{X_{2}} ).

Note that we have re-written the Arrhenius law as kτ = f(X_{2}) = δ exp(γ−\frac{γ}{X_{2}} ). The dimensionless numbers γ and δ can be thought of as non-dimensional counterparts to the activation energy E and the rate constant k_{0}. Next week we will estimate γ and δ by fitting the function kτ = δ exp(γ−\frac{γ}{X_{2}} ) to measured reaction rates, but as for today you may use the values γ = 30 and δ = 0.1. The dimensionless number α, on the other hand, can be directly computed from known data.

• Introducing the dimensionless number

β = \frac{κA_{K}}{ρc_{p}q_{ref}}

the third term on the right-hand side becomes

−β(X_{2} − U_{2}). (7)

From (4)–(7), (3) becomes

(8)  \frac{dX_{2}}{ds} = U_{1}(1 − X_{2}) + αX_{1}f(X_{2}) − β(X_{2} − U_{2}),

and replacing X by X_{1}, U by U_{1}, and kτ by f(X_{2}), in the equation for the concentration that we derived last week, we arrive at the system of differential equations,

\frac{dX_{1}}{ds} = U_{1}(1 − X_{1}) − X_{1}f(X_{2}), s > 0,
(9)  \frac{dX_{2}}{ds} = U_{1}(1 − X_{2}) + αX_{1}f(X_{2}) − β(X_{2} − U_{2}), s > 0,
X_{1}(0) = X_{1,0}, X_{2}(0) = X_{2,0}.

Exercise 3. Extend your Matlab programs from last week so that they solve (9). Let first U_{1}

and U_{2} be equal to constant values Ū_{1} and Ū_{2}. Does the solution X(s) = [\begin{matrix}X_{1}(s) \\ X_{2}(s) \end{matrix}] approach an

equilibrium X̄ = [X̄_{1} ; X̄_{2}] as s → ∞? In fact, you should be able to find two equilibrium points by

choosing different initial values X_{0}, say, X_{0} = [\begin{matrix}0.5 \\ 1 \end{matrix}] and X_{0} = [\begin{matrix}0.5 \\ 1.1 \end{matrix}].

Hint: The file data.m could for example look like

  global alpha beta gamma delta Ubar
  dH = -83.7e3; % heat of reaction
  cf = 5e3; % concentration at inflow
  % here you initiate the rest of the variables needed to compute alpha
  alpha = -dH*cf/(rho*cp*Tf);
  % here you initiate the variables needed to compute beta
  % you must choose a value for AK, for instance AK = 1
  beta = ...; % insert the expression for beta
  gamma = 30;

delta = 0.1;
Ubar = [1; 0.97];
S = 10; % final time
X0 = [0.5; 1]; % initial value

Change the function file tank.m to a file tank2.m beginning with

function y=tank2(s,X)
global alpha beta gamma delta Ubar
U(1) = Ubar(1);
U(2) = Ubar(2);
y = zeros(2,1);
y(1) = U(1)*(1-X(1)) - ... ; % complete the expression for dX1/ds
y(2) = ... ; % insert the expression for dX2/ds

(If you want, you can write a separate function file f.m that implements the “Arrhenius function” f(X_{2}) = δ exp(γ−\frac{γ}{X_{2}} ). Then you can call f from tank2. Just don’t forget to declare the variables gamma and delta as global inside f.)

You should finally be able to use the file solve2.m from Exercise 1 with obvious modifications. In a later exercise we will look for an equilibrium at X̄_{1} = 0.5 and analyze the stability of this

desired operating point.

Finally we give a some reasons for writing the mathematical model in dimensionless form:

• the dimensionless equations (9) contain fewer constants and variables than the original equations and they are therefore more convenient to work with;

• it gives a possibility of scaling the variables so that they are not very big or very small;

• it is useful for scaling an apparatus from laboratory size to factory size; this should be done so that the dimensionless constants are the same.