TMA681 Tillämpad matematik K, 2001

STUDIO 1. THE TANK REACTOR: MASS 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. Linear differential equations of first order

(MA2 9.2) Recall that \frac{d}{dt} exp(t) = exp(t) and exp(0) = 1. This means that the exponential function u(t) = exp(t) satisfies the initial-value problem

u′(t) = u(t), t > 0; u(0) = 1. (1)

Let a and u_{0} be real numbers and consider the initial-value problem

u′(t) = au(t), t > 0; u(0) = u_{0}. (2)

The solution is u(t) = u_{0} exp(at). Check this!

Exercise 1. Solve (2) with Matlab. Write the following function file funk1.m, which defines the right-hand side of the differential equation.

  function y=funk1(t,u)
  global a
  y=a*u;

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

  global a
  a=1; T=1; u0=1;

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

  [t,u]=ode45('funk1', [0 T], u0);
  plot(t,u);

Note how the value of a enters into the function via the command global a which is written both in the function file funk1.m and in the main program data.m. Start the computation by typing the following on the Matlab command line:

>> data
>> solve1
>> a=-1
>> T=2
>> solve1

Compute and observe the solutions for various values of a, positive, negative, and zero. If you type » hold on then Matlab will plot several curves in the same figure.

Now let b be another real number and consider the initial-value problem

u′(t) = au(t) + b, t > 0; u(0) = u_{0}. (3)

Remember (MA2 10.5) that the solution is given by u(t) = u_{h}(t) + u_{p}(t), where u_{h}(t) = A exp(at) is a solution of the homogeneous equation u′ − au = 0 and u_{p}(t) = B is a particular solution of the inhomogeneous equation u′ − au = b.

Date: March 28, 2001, Niklas Ericsson, Stig Larsson, Mathematics, Chalmers University of Technology.

Homework 1. Determine the constants A and B and show that the solution is

(4) u(t) = u_{0} exp(at) + \frac{b}{a}(exp(at) − 1) (if a ≠ 0).

What is the solution when a = 0? For which values of a does the solution approach an equilibrium u(t) → ū as t → ∞? Determine ū. Exercise 2. Solve (3) with Matlab. Define b in the file data.m and change the global command to global a b in both places. Observe what the solutions look like for a positive, negative, and zero. For which values of a does the solution approach an equilibrium u(t) → ū as t → ∞?

Now we consider the situation when a = a(t), b = b(t) are not constant:

u′(t) = a(t)u(t) + b(t), t > 0; u(0) = u_{0}. (5)

We use the method of integrating factor (MA2 9.2). Write the equation as u′ − a(t)u = b(t) and multiply it by the integrating factor e^{−A(t)}, where A(t) = ∫_{0}^{t} a(s) ds, so that A′(t) = a(t) and A(0)=0. We get

\frac{d}{dt}(e^{−A(t)}u(t)) = e^{−A(t)}u′(t) − a(t)e^{−A(t)}u(t) = e^{−A(t)}b(t).

We integrate from 0 to T:

[e^{−A(t)}u(t)]_{0}^{T} = ∫_{0}^{T} e^{−A(t)}b(t) dt

and hence

e^{−A(T)}u(T) − e^{−A(0)}u(0) = ∫_{0}^{T} e^{−A(t)}b(t) dt.

Using A(0) = 0, u(0) = u_{0}, multiplying by e^{A(T)}, and replacing t by s, T by t, we finally get

u(t) = u_{0}e^{A(t)} + ∫_{0}^{t} e^{A(t)−A(s)}b(s) ds, where A(t) = ∫_{0}^{t} a(s) ds. (6)

Exercise 3. Solve (5) with Matlab. Use, for example, a = ±1, b(t) = 1 + 0.1 sin(7t). Homework 2. Use constant functions a(t) = a, b(t) = b in (6) and compute the integrals to obtain (4).

2. The tank reactor

In a sequence of studio sessions we will study the ideal mixed tank reactor, see Figure 1. The goal is to design the reactor so that it will operate in a stable way at (c_{f} − c)/c_{f} = 0.5 (“50 % omsättningsgrad”).

We will use the following data:

V = 1.2 m^{3} volume of the tank

c_{p} = 4.19 kJ/(kg K) heat capacity of the fluid

ρ = 1000 kg/m^{3} density of the fluid

q_{ref} = 3.0 · 10^{−4} m^{3}/s reference flux

c_{f} = 5 kmol/m^{3} concentration at inflow of the interesting substance T_{f} = 70^{∘}C temperature at inflow δH = −83.7 kJ/mol heat of reaction

A_{K} = 1.0 m^{2} area of the cooler

κ = 0.58 kJ/(m^{2} s K) heat transfer coefficient of the cooler PSfrag replacements

c
c_{f}
q
T
T_{f}
T_{K}
q_{K}
T_{Kf}
A_{K}
V
V_{K}

Figure 1. The tank reactor.

We introduce the following variables:

c [mol/m^{3}] concentration in the reactor of the interesting substance T [K] temperature in the reactor

q [m^{3}/s] flux through the reactor

T_{K}, T_{Kf} [K] temperatures in the cooler and at the cooler inflow

q_{K} [m^{3}/s] flux through the cooler

V_{K} [m^{3}] volume of the cooler

The reaction is exothermal and of first order with rate of reaction kc [mol/(m^{3}s)]. The rate coefficient depends on the temperature according to the Arrhenius law:

k = k_{0} exp(−E/(RT)) [s^{−1}] (7)

where R [8.31 J/(mol K)] is the gas constant, E [J/mol] is the activation energy and k_{0} [s^{−1}] is the rate constant of the reaction. The following rates have been measured:

\begin{matrix}T [K] & 343 & 353 & 363 & 373 & 383 & 393 & 403 \\ k [s^{−1}] & 2.8 · 10^{−5} & 5.6 · 10^{−5} & 11.2 · 10^{−5} & 22.4 · 10^{−5} & 44.8 · 10^{−5} & 89.6 · 10^{−5} & 179.2 · 10^{−5} \end{matrix}

(Based on “Kemisk reaktionsteknik. Övningsuppgifter”, Kemisk reaktionsteknik, CTH 1993, uppgift 7.1.)

2.1. Mass balance. When we build a mathematical model for the tank reactor it is important not to consider all aspects at once. We therefore begin by making a big simplification: we assume that the rate coefficient k is constant, for example, k =?? corresponding to T =?? (choose a value from the table). Recall that k depends strongly on T, see (7), so this is not very realistic, but it will be a good starting point for our investigation.

The mass balance equation is

V \frac{dc}{dt} = q(c_{f} − c) − ckV. [mol/s]

This equation says that the rate of change of the total amount of the interesting substance is equal to the influx minus the outflux minus the reaction rate. In order to make the equation dimensionless we divide by q_{ref}c_{f} [mol/s]. We get

(8) \frac{V}{q_{ref}} \frac{d}{dt}(\frac{c}{c_{f}} ) = \frac{q}{q_{ref}} (1 − \frac{c}{c_{f}} ) − \frac{c}{c_{f}} k \frac{V}{q_{ref}} .

We introduce τ = V/q_{ref} [s] (“uppehållstid”) and the dimensionless variables

s = t/τ (dimensionless time),
X = c/c_{f} (dimensionless concentration),
U = q/q_{ref} (dimensionless flux).

Then, by the chain rule,

\frac{dX}{ds} = \frac{dX}{dt} \frac{dt}{ds} = τ \frac{dX}{dt} = \frac{V}{q_{ref}} \frac{d}{dt}(\frac{c}{c_{f}} ),

and (8) becomes

\frac{dX}{ds} = U(1 − X) − Xkτ,

that is,

\frac{dX}{ds} = −(kτ + U)X + U; X(0) = X_{0}. (9)

Note that this is of the form (5) with a = −(kτ + U) and b = U. Exercise 4. Change your Matlab programs from part 1 so that they solve (9). Let first U be equal to a constant value Ū. Does the solution X(s) approach an equilibrium X̄ as s → ∞? Hint: Begin the file data.m by

  global ktau Ubar
  V=1.2;   qref=3.0e-4;
  k=?? ;
  tau=V/qref;   ktau=k*tau;
  Ubar=1;  % equilibrium value of U

Change also the function file funk1.m to a file tank.m beginning with

  function y=tank(s,X)
  global ktau Ubar
  U=Ubar;

Exercise 5. Recall that we want the reactor to operate at (c_{f} − c)/c_{f} = 0.5, i.e., at X̄ = 0.5. Determine Ū so that this is achieved. Hint: the equation for X̄, Ū is obtained by setting \frac{dX}{ds} = 0 in (9), i.e.,

−(kτ + Ū)X̄ + Ū = 0. (10)

Find a formula for Ū in terms of kτ and X̄. Explain what happens in the extreme cases X̄ = 0, X̄ = 1. Insert the command Xbar=.5; in the file data.m and the formula for Ubar in the file

solve1.m.

Exercise 6. We now investigate if the operating point X̄ = 0.5 is stable with respect to changes of the initial value X_{0}. Set U = Ū and solve (9) with Matlab with several initial values X_{0}. Is it stable? Exercise 7. We next investigate if the operating point X̄ = 0.5 is stable with respect to changes of the flux U. Set X_{0} = X̄ = 0.5 and solve (9) with Matlab with, for example, U(s) = Ū +0.1 sin(7s). Is it stable? Exercise 8. Finally, compute the required flux q̄ (in dimensional units [m^{3}/s]).

Next week we will include the temperature dependence in k and we will see that the operating point X̄ = 0.5 is not always stable then.