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.

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.0m^{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.

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,
(11)  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 9. Solve (11) 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 (11) 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 3. Verify these properties.

Exercise 10. 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 (11). 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

(12) ρ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 (12) 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

(13) \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 (13) 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}} ). (14)

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

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

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

• 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

(16)  α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}). (17)

From (14)–(17), (13) becomes

(18)  \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,
(19)  \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 11. Extend your Matlab programs from last week so that they solve (19). 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 9 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.

STUDIO 3 The tank reactor: Arrhenius’ law.

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. The method of least squares

(TLA 4.9) Consider the linear system of equations

(20)  Ax = b,

where A ∈ R^{m×n}, x ∈ R^{n}, b ∈ R^{m}. If m > n (more equations than unknowns), then the system is “overdetermined” and such a system has no solution in general. Geometrically, the reason for this is that (in general) the vector b lies outside the range (“värderummet”) of A,

R(A) = {y ∈ R^{m} : y = Ax for some x ∈ R^{n}}.

Since the dimension of R(A) is ≤ n and m > n, we realize that the space R(A) does not “fill out” the whole space R^{m}. Therefore it is likely that a given vector b ∈ R^{m} will lie outside R(A), see Figure 1, and then Ax cannot be equal to b and (20) has no solution.

PSfrag replacements

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

b − Ax̂

Ax̂ Av R(A)

Figure 1. Orthogonal projection onto the range of A.

In this situation we seek an approximate solution which makes the residual

b − Ax

as small as possible. More precisely, we seek a vector x̂ ∈ R^{n}, which minimizes the square of the norm (length) of the residual:

(21)  f(x̂) = min f(x), f(x) = ∥b − Ax∥^{2}.

Recall the scalar product ⟨x, y⟩ = y^{t}x and the norm ∥x∥ = \sqrt{⟨x, x⟩} = \sqrt{x^{t}x} of column vectors. From (TLA 6.9) we know that there is a unique vector ŷ = Ax̂ ∈ R(A) such that the distance ∥b − Ax∥ is minimal, i.e., ∥b − Ax̂∥ ≤ ∥b − Ax∥ for all x. The vector ŷ = Ax̂ is the orthogonal projection of b onto R(A), see Figure 1. It is characterized by the condition that b − ŷ = b − Ax̂ is orthogonal to all vectors Av ∈ R(A). This means that

0 = ⟨b − Ax̂, Av⟩ = (Av)^{t}(b − Ax̂) = v^{t}A^{t}(b − Ax̂) = v^{t}(A^{t}b − A^{t}Ax̂).

Since this holds for all v ∈ R^{n}, we may take v = A^{t}b − A^{t}Ax̂ to get

(22)  A^{t}Ax̂ − A^{t}b = 0.

Thus, we can compute x̂ by solving the linear system

(23)  A^{t}Ax = A^{t}b.

Note that the coefficient matrix A^{t}A is n × n and symmetric. The system (23) has at least one solution (namely x̂).

In order to see that the minimization problem (21) is equivalent to solving the linear system (23), we write x = x̂ + v, v = x − x̂, and compute

\begin{aligned}f(x) & = f(x̂ + v) = ∥(b − Ax̂) − Av∥^{2} = ⟨(b − Ax̂) − Av, (b − Ax̂) − Av⟩ \\ & = ⟨b − Ax̂, b − Ax̂⟩ − 2⟨b − Ax̂, Av⟩ + ⟨Av, Av⟩ \\ & = ∥b − Ax̂∥^{2} + 2(Av)^{t}(Ax̂ − b) + ∥Av∥^{2} \\ & = f(x̂) + 2v^{t}A^{t}(Ax̂ − b) + ∥Av∥^{2} \\ & = f(x̂) + 2v^{t}(A^{t}Ax̂ − A^{t}b) + ∥Av∥^{2}. \end{aligned} ; (24)

Taking (22) into account we get

f(x) = f(x̂) + ∥Av∥^{2} ≥ f(x̂).

This shows that x minimizes f(x), if and only if ∥Av∥ = 0, in which case Ax = Ax̂ + Av = Ax̂ and A^{t}Ax = A^{t}Ax̂ = A^{t}b. Therefore, x minimizes f(x) if and only if x is a solution of (23).

We can also interpret this in terms of the general minimization problem in (MA3 2.9–10). Recall Taylor’s formula:

(25)  f(x) = f(x̂ + v) = f(x̂) + v^{t}f′(x̂) + \frac{1}{2}v^{t}f′′(x̂)v + R(x).

Noting that ∥Av∥^{2} = ⟨Av, Av⟩ = (Av)^{t}(Av) = v^{t}(A^{t}A)v, we re-write (24) as

(26)  f(x) = f(x̂ + v) = f(x̂) + 2v^{t}(A^{t}Ax̂ − A^{t}b) + v^{t}(A^{t}A)v.

Comparing (25) with (26), we identify the Jacobi matrix (gradient vector) f′(x̂) = 2(A^{t}Ax̂−A^{t}b), the Hesse matrix f′′(x̂) = 2A^{t}A, and the remainder R(x) = 0. Recall that stationary points are given by the system of equations f′(x) = 2(A^{t}Ax − A^{t}b) = 0, which is the same as (23). Note also that the Hesse matrix is constant (with respect to x) and positive semidefinite: v^{t}f′′(x)v = 2v^{t}(A^{t}A)v = 2∥Av∥^{2} ≥ 0.

Exercise 12. Suppose that the variables y and x are related by y = kx + m. In order to determine the coefficients k and m we make measurements of y and x:

\begin{matrix}x & 5 & 6 & 7 & 8 & 9 & 10 \\ y & 19.5888 & 23.4043 & 25.5754 & 29.1231 & 31.9575 & 35.8116 \end{matrix}

This leads to an overdetermined system of the form

kx_{1} + m = y_{1}
⋮
kx_{6} + m = y_{6}

or, in matrix form Av = y,

[\begin{matrix}x_{1} & 1 \\ ⋮ & ⋮ \\ x_{6} & 1 \end{matrix}] [\begin{matrix}k \\ m \end{matrix}] = [y_{1} ; ⋮ ; y_{6}] .

Solve this system by the least squares method in Matlab. Hint: set up the column vectors x, y and the matrix A=[x ones(size(x))], then form the matrices B=A'*A and g=A'*y. Solve the system Bv = g by the command v=B\g.

Plot the data points (x_{i}, y_{i}) and the fitted function y = kx + m in the same figure. The

following commands are useful: plot(x,y,'or'), fplot('ykxm',[x(1) x(6)]). Here ykxm.m is

a function file that implements the function y = kx + m. Don’t forget to declare global k m both inside the function file and in the main program.

Actually, Matlab’s backslash command v=A\y automatically uses the least squares method when the system Av = y is overdetermined. Try this also!

2. The tank reactor

The rate coefficient depends on the temperature according to the Arrhenius law:

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

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}

The task is now to determine the coefficients k_{0} and E by fitting the rate law (27) to these data. Last week we wrote (27) in dimensionless form

(28)  kτ = δe^{γ(1−1/X)}, where γ = \frac{E}{RT_{f}} , δ = k_{0}τe^{−γ}, X = \frac{T}{T_{f}} , τ = \frac{V}{q_{ref}} .

Introducing new variables r = kτ and ξ = 1 − 1/X we get

(29)  r = δe^{γξ}.

The task is now to fit this function to the given data points (ξ_{i}, r_{i}).

Exercise 13. (Linear least squares method.) Form the logarithm of (29) so that you get a linear relation of the form y = kx + m, namely, (see also (TLA 5.2))

(30)  log(r) = γξ + log(δ).

(Note that the natural logarithm is denoted log(x) in English and in Matlab, but ln(x) in Swedish.) Solve for γ and δ by using the least squares method as in Exercise 12. Begin by forming column vectors X, r, xi and so on. Plot the data points (X_{i}, r_{i}) and the fitted function r = δe^{γ(1−1/X)} in the same figure. Finally, determine k_{0} and E.

Homework 4. (A nonlinear least squares method.) Alternatively, we can form the residual r − δe^{γξ} from the nonlinear relation (29) and minimize the square of its norm

(31)  g(δ, γ) = ∑_{i} (r_{i} − δe^{γξ_{i}})^{2}.

Write a Matlab function that implements this function and use Matlab’s program fminsearch to minimize it. Why does this method give a slightly different result? Hint: the Matlab function norm may be useful for computing the right side of (31).

Exercise 14. Insert the new values for δ and γ in your Matlab programs from Studio 2. Repeat all the computations. Let 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}].

12 3. THE TANK REACTOR: ARRHENIUS’ LAW.

Next week we will look for an equilibrium at X̄_{1} = 0.5 and analyze the stability of this desired operating point.

STUDIO 4 The tank reactor: Stability.

Theory: S. Larsson, “Kompletterande föreläsningsanteckningar och övningar” (K).

1. Introduction

(K 1.1) Recall from Studio 2 the system of differential equations,

(32)  \begin{aligned}\frac{dX_{1}}{ds} & = U_{1}(1 − X_{1}) − X_{1}f(X_{2}) = F_{1}(X, U), s > 0, \\ \frac{dX_{2}}{ds} & = U_{1}(1 − X_{2}) + αX_{1}f(X_{2}) − β(X_{2} − U_{2}) = F_{2}(X, U), s > 0, \\ X_{1}(0) & = X_{1,0}, X_{2}(0) = X_{2,0}, \end{aligned}

which is our mathematical model for the dynamics of the tank reactor. Recall, also, the state variables X_{1} = X_{1}(s) (dimensionless concentration) and X_{2} = X_{2}(s) (dimensionless reactor temperature), depending on s (dimensionless time). These two variables, that describe the state of the tank reactor, are the ones that we compute by solving (32), i.e., X_{1} and X_{2} are output data.

There are two kinds of input data. First we have the initial data, X_{0} = [\begin{matrix}X_{1,0} \\ X_{2,0} \end{matrix}], then the control

variables U_{1} = U_{1}(s) (dimensionless flux) and U_{2} = U_{2}(s) (dimensionless cooler temperature).

Recall, finally, that kτ = f(X_{2}) = δe^{γ(1−1/X_{2})} is the dimensionless rate coefficient given by the Arrhenius law. In Studio 3 you determined the dimensionless numbers γ and δ, which can be thought of as non-dimensional counterparts to the activation energy and the rate constant of the reaction, by fitting the rate law to data. As for the other two dimensionless numbers appearing in (32), α is also reaction dependent since it is proportional to the heat of reaction, whereas β is proportional to the area and the heat transfer coefficient of the cooler.

Our final goal is to design the tank reactor in such a way that it runs in a stable manner at

a specified, desired, operating point X̄ = [X̄_{1} ; X̄_{2}]. In order to achieve this goal, we first determine

corresponding values of the control variables Ū = [Ū_{1} ; Ū_{2}], for which X̄ is a stationary point. Then,

we analyze the stability of X̄ with respect to perturbations of input data and, if necessary, change the value of one or more design parameters. This could, for instance, physically mean varying the area of the cooler, i.e., varying the value of β in our mathematical model (32).

2. Stationary points

(K 1.2) Our first task is, given a desired operating (stationary) point X̄ = [X̄_{1} ; X̄_{2}], to determine

corresponding values of the control variables Ū = [Ū_{1} ; Ū_{2}] by solving the system of equations,

(33)  \begin{aligned}0 & = Ū_{1}(1 − X̄_{1}) − X̄_{1}f(X̄_{2}), \\ 0 & = Ū_{1}(1 − X̄_{2}) + αX̄_{1}f(X̄_{2}) − β(X̄_{2} − Ū_{2}). \end{aligned}

Exercise 15. Show that the solution to (33) is given by,

(34)  \begin{aligned}Ū_{1} & = \frac{X̄_{1}}{1 − X̄_{1}} f(X̄_{2}), \\ Ū_{2} & = X̄_{2} − \frac{1}{β} (\frac{X̄_{1}}{1 − X̄_{1}} (1 − X̄_{2})f(X̄_{2}) + αX̄_{1}f(X̄_{2})). \end{aligned}

We here stress a subtle point. Given X̄, clearly, Ū is uniquely defined by (34). On the other hand, if we instead consider Ū as fixed, we know that X̄ is one solution to (33), however, it might not be the only solution, i.e., there may exist more than one stationary point corresponding to Ū. We will return to this later.

To concretize, let us now specify the values of the state variables at the desired operating point: our objective is to design the tank reactor to operate in a stable manner at (c_{f} −c̄)/c_{f} = 0.5 (“50%

omsättningsgrad”) and at reactor temperature T̄ = 99^{∘}C, i.e., at X̄ = [X̄_{1} ; X̄_{2}] = [\begin{matrix}0.5 \\ (99 + 273.15)/T_{f} \end{matrix}].

Exercise 16. Modify the file data.m from Studio 2, so that, given X̄ (as above), Ū is computed from (34). At the same time, check that you have changed the old values of γ (= 30) and δ (= 0.1) that we used in Studio 2 to the new ones that you determined in Studio 3. Also, check that you have set A_{K} = 1 m^{2}. Hint: You need to change the line

  Ubar = [1; 0.97];

into

  Xbar = [0.5; (99 + 273.15)/Tf]; % Tf = 70 + 273.15
  Ubar = zeros(2,1); % initialize (column vector) Ubar
  Ubar(1) = ...; % insert the expression for Ubar(1)
  Ubar(2) = ...; % insert the expression for Ubar(2)

Exercise 17. Check data.m by calling tank2.m that you wrote in Studio 2:

>> data
>> Xprime = tank2(0, Xbar)

What should the result be? (Note that the value of the first argument may be given arbitrarily, since there is no explicit time dependence in the right-hand side of (32).)

3. Instability of the operating point

We now perform a first stability check of the operating point X̄. We do this by introducing small initial perturbations, i.e., small initial deviations from X̄, in X. In this test, we do not consider perturbations in the control variables, i.e., we take U = Ū in tank2.m.

Exercise 18. Assuming that solve2.m is the name of your script file from Studio 2, from which the call to ode45 is made and the solution is plotted, solve (32) by giving the following Matlab commands:

>> data
>> S = 20;
>> X0 = Xbar + [0; 0.05];
>> solve2
>> X0 = Xbar - [0; 0.05];
>> solve2

Is X̄ stable with respect to these perturbations? Also try some other initial perturbations.

As you have just seen, a small deviation from X = X̄ causes the tank reactor to depart from the desired operating point. Since these kinds of perturbations are inevitable in practice, the reactor will not remain in the desired state, which is therefore not stable. Rather, it will (depending on the initial perturbation) reach one of two other equilibrium points, which seem to be stable ones. These two are also stationary points, corresponding to Ū, i.e., they are also solutions to (33). This is the non-uniqueness mentioned in Section 33.

4. Linear stability analysis

(K 1.3) In order to learn how to “adjust” the tank reactor so that it will operate in a stable way at X̄, we need to systematically study the stability of solutions to (32). We will perform a linear stability analysis based on the assumption of small perturbations.

Let X(s) with input data X_{0}, U(s) be a solution to (32) that is close to X̄. With

X(s) = X̄ + ΔX(s), X_{0} = X̄ + ΔX_{0}, U(s) = Ū + ΔU(s), (35)

we may consider ΔX(s) as a perturbation in X(s) caused by the perturbations ΔX_{0} and ΔU(s) in input data.

If ΔX(s) and ΔU(s) are small, we obtain the linear system

(36)  \begin{matrix}x′(s) = Ax(s) + Bu(s), s > 0, \\ x(0) = x_{0}, \end{matrix}

for the approximate perturbation x(s) ≈ ΔX(s) caused by the perturbations in input data x_{0} = ΔX_{0} and u(s) = ΔU(s). In (36),

(37)  A = [\begin{matrix}\frac{∂F_{1}}{∂X_{1}} (X̄, Ū) & \frac{∂F_{1}}{∂X_{2}} (X̄, Ū) \\ \frac{∂F_{2}}{∂X_{1}} (X̄, Ū) & \frac{∂F_{2}}{∂X_{2}} (X̄, Ū) \end{matrix}] = [\begin{matrix}−Ū_{1} − f(X̄_{2}) & −X̄_{1}f′(X̄_{2}) \\ αf(X̄_{2}) & −Ū_{1} + αX̄_{1}f′(X̄_{2}) − β \end{matrix}] ,

where f′(X̄_{2}) = \frac{γ}{X̄_{2}^{2}} f(X̄_{2}), and

(38)  B = [\begin{matrix}\frac{∂F_{1}}{∂U_{1}} (X̄, Ū) & \frac{∂F_{1}}{∂U_{2}} (X̄, Ū) \\ \frac{∂F_{2}}{∂U_{1}} (X̄, Ū) & \frac{∂F_{2}}{∂U_{2}} (X̄, Ū) \end{matrix}] = [\begin{matrix}1 − X̄_{1} & 0 \\ 1 − X̄_{2} & β \end{matrix}] ,

are called Jacobi matrices of F (X, U) = [\begin{matrix}F_{1}(X, U) \\ F_{2}(X, U) \end{matrix}] at X̄, Ū.

Homework 5. Verify (37) and (38).

5. Stability with respect to perturbations of initial data

(K 1.4) In this section we consider the case u(s) = 0, i.e., we only consider perturbations in initial data. In this case, (36) simplifies to

(39)  \begin{matrix}x′(s) = Ax(s), s > 0, \\ x(0) = x_{0}, \end{matrix}

with solution (we assume that A is diagonalisable)

(40)  x(s) = c_{1}e^{λ_{1}s}g_{1} + c_{2}e^{λ_{2}s}g_{2},

where λ_{i}, g_{i} are eigen-values and eigen-vectors of A, and the c_{i} are constants depending on x_{0}. Clearly, the growth of x(s) (and accordingly the stability of X̄) depends on the eigen-values of A.

Exercise 19. Compute the eigen-values of A using Matlab. Hint: First write the function file jacobianA.m that computes A:

  function A = jacobianA(Xbar)
  global alpha beta gamma delta Ubar
  A = zeros(2,2); % initialize (2x2 matrix) A
  A(1,1) = ...; % insert the expression for A(1,1)
  A(1,2) = ...; % insert the expression for A(1,2)
  A(2,1) = ...; % insert the expression for A(2,1)

  A(2,2) = ...; % insert the expression for A(2,2)

Then you can compute the eigen-values of A by typing:

>> data
>> global A % we declare A as global because we will need this later
>> A = jacobianA(Xbar)
>> eig(A)

As you (hopefully!) noticed, A has two real eigen-values, one positive and one negative. Because of the positive eigen-value one of the terms in (40) will grow exponentially with time, and this explains the instability of X̄.

Exercise 20. Solve (39) with the same initial perturbations as in Exercise 18. Hint: First write the function file lineartank.m that computes the right-hand side of (39):

  function y = lineartank(s,x)
  global A % this is the reason we declared A as global
  y = A*x;

Then modify solve2.m into linearsolve.m. (Just replace tank2 by lineartank in the call to ode45, and X by x everywhere.) Now you can solve (39) by typing:

>> figure % opens a new figure
>> data
>> S = 1;
>> x0 = [0; 0.05];
>> linearsolve
>> x0 = [0; -0.05];
>> linearsolve

One clearly sees the perturbation growth. Note that (39) was derived on the basis of the assumption of small perturbations and that it is not valid if x(s) becomes “too” large. So there is no point in computing much further than to S = 1.

It is instructive to compare X(s), computed as in Exercise 18, to X̄+x(s), with x(s) computed as in Exercise 20:

Exercise 21. First solve (32), as in Exercise 18, with S = 1 and X0 = Xbar + [0; 0.05]. Then solve (39), as in Exercise 20, with S = 1 and x0 = [0; 0.05]. Now compare the first solution, X, to Xbar + x, where x is the second solution. You can do the comparison by writing and running the following script:

  clf % clear current figure
  plot(X(:,1), X(:,2)) % plots second versus first component of X
  hold on
  plot(Xbar(1) + x(:,1), Xbar(2) + x(:,2), '--')
  title('Phase portraits:   Solid: X   Dashed: Xbar + x')
  xlabel('X_1, Xbar_1 + x_1')
  ylabel('X_2, Xbar_2 + x_2')
  hold off

Note how the two curves successively diverge, and how the linear approximation fails to find the stable equilibrium point.

Next week we will conclude the exercise on the tank reactor by “adjusting” it in such a way that the operating point X̄ becomes stable. The idea is to try to “move” the eigen-values of A so that their real parts get the right (negative!) sign. We will also briefly consider stability with respect to perturbations in the control variables.