TMA681 Tillämpad matematik K, 2001

STUDIO 4. THE TANK REACTOR: STABILITY.

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.

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

1. Introduction

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

(1)  \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 (1), 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 (1), α 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 (1).

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,

(2)  \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}

Date: April 12, 2001, Niklas Ericsson, Stig Larsson, Mathematics, Chalmers University of Technology.

Exercise 1. Show that the solution to (2) is given by,

(3)  \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 (3). On the other hand, if we instead consider Ū as fixed, we know that X̄ is one solution to (2), 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 2. Modify the file data.m from Studio 2, so that, given X̄ (as above), Ū is computed from (3). 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 3. 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 (1).)

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 set U = Ū in tank2.m. Exercise 4. 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 (1) 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 (2). This is the non-uniqueness mentioned in Section 2.

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 (1). 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 (1) that is close to X̄. With

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

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

(5)  \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 (5),

(6)  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

(7)  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 1. Verify (6) and (7).

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, (5) simplifies to

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

with solution (we assume that A is diagonalizable)

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

where λ_{i}, g_{i} are eigenvalues and eigenvectors 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 eigenvalues of A. Exercise 5. Compute the eigenvalues 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 eigenvalues 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 eigenvalues, one positive and one negative. Because of the positive eigenvalue one of the terms in (9) will grow exponentially with time, and this explains the instability of X̄. Exercise 6. Solve (8) with the same initial perturbations as in Exercise 4. Hint: First write the function file lineartank.m that computes the right-hand side of (8):

  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 (8) 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 (8) 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 4, to X̄ + x(s), with x(s) computed as in Exercise 6: Exercise 7. First solve (1), as in Exercise 4, with S = 1 and X0 = Xbar + [0; 0.05]. Then solve (8), as in Exercise 6, 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 eigenvalues 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.