Question 5.PS.1: Non-Isothermal CSTR In Chapter 4, we presented two case stud...

Note that this function takes as inputs not only the vector of unknowns x, but also the parameters β, γ , and Da. For convenience, we unpack the input vector into the concentration, x(1), and the temperature, x(2).

Since we put the unknown vector in the logical order [x_1, x_2], the Jacobian is

J=\begin{bmatrix} −1 −Dae^{−β/x_2} & −Dax_1e^{−β/x_2}βx^{−2}_2 \\Daγ e^{−β/x_2} & −1 + Daγ x_1e^{−β/x_2}βx^{−2}_2\end{bmatrix}                                                             (5.7.22)
which is computed from a similar MATLAB function:

As was the case in getR, we are passing the problem parameters to getJ, and we have unpacked the unknown vector for convenience.

For a given value of the residence time τ , and thus a given Dämkohler number, Da, we need to use a slight modification of Program 3.5 to implement Newton–Raphson:

The obvious difference with the original program is that we are sending the values of the parameters β, γ , and Da. The more subtle difference is in line 10, where we cut off the maximum number of iterations at kmax = 20, whereas our original program allowed this number to go to 100. We will explain the reason for this choice in a moment. This problem has multiple steady states including two turning points that make the continuation method fail. The strategy to find the steady states is as follows.

(1) Fill in the lower part up to the first turning point starting from a low temperature.
(2) Fill in the upper part down to the second turning point starting from a high temperature.
(3) Find an intermediate point (the trickiest part).
(4) Fill up to the second turning point from the intermediate point.
(5) Fill down to the first turning point from the intermediate point.
For each “fill” step, we wrote a program to perform the continuation method until it fails.

The inputs to this function are the initial guesses for the concentration and temperature, contained in the vector x, at the first value of the residence time τ . We also send the parameters γ and β, as well as the step size Δτ that we are going to use for the continuation method. Remember from Section 3.9 that we need to take a small step size if we use a zero-order continuation method, which is the case here.

The main part of continuation_tau is a while loop that executes so long as a new solution was found in the last iteration. Before we try to find this new solution, we first convert τ to Da with a simple program:

You may be wondering why we bothered to write a separate program with such a simple calculation when we could have just put this on line 5 of continuation_tau. The reason is that we decided to work in terms of τ , rather than Da, for the continuation. If at some point we decide to change the problem parameters, in particular the conversion factor in Eq. (5.7.21), we would need to remember to make the change in continuation_tau. You may be very good at remembering the details of your programs, but we are getting old! By separating this calculation into its own subfunction, we are less likely to forget where we put Eq. (5.7.21).

Now let’s go back to continuation_tau. Line 7 calls Newton–Raphson to compute the concentration and temperature at the current value of τ . Recall from our modified version of nonlinear_newton_raphson.m that we will stop the algorithm if we go past 20 iterations or the solution diverges. In either case, the value of flag will be zero. If this happens, we know that Newton–Raphson failed and line 9 stops the continuation method. If Newton–Raphson converged, then we record the value of the temperature and the residence time in the output vectors in lines 11 and 12. To get ready for the next iteration, we update the value of the residence time in line 13 and update our counter for the output vectors in line 14 so that we can record the new values in the right spot. Our algorithm is a zero-order continuation method since, when the while loop executes again, we use the value of x from the last iteration as the initial guess for the next iteration.

Let’s see how well our program works. For the lower branch, we start with a very small residence time of τ = 0.1 and step in increments of Δτ = 0.005. If the residence time is small, then there should be hardly any reaction at all. We thus selected a dimensionless concentration of x^{(0)}_1 = 0.9 and temperature of x^{(0)}_2 = 1, corresponding to a small decrease in concentration with no accompanying temperature rise. As we see in Fig. 5.18, the continuation method works well until τ = 1.870. When we try to solve the problem for τ = 1.875 using the converged solution at τ = 1.870, we cannot get convergence. Thus, it is reasonable to assume that this is the turning point on the lower branch and we need a better initial guess to find a solution.

For the upper branch, we now need to consider a residence time that exceeds the end of the lower branch. We chose to use an initial residence time τ = 2 and the same increment Δτ = −0.005, except that for the upper branch we are trying to march down in residence time. If you think about this case for a moment, you would expect that a large residence time should lead to substantial reaction. Thus, we would expect a high temperature and low outlet concentration. We decided to use x^{(0)}_1 = 0.04 , i.e. 96% conversion, and x^{(0)}_2 = 1.14 as the initial guess. You might not think that a 14% increase in the temperature relative to the feed is very hot, but remember that this is in absolute temperature. So if we feed the reactor near room temperature, say 300 K, then the outlet of the reactor (342 K = 156ºF) exceeds the hottest weather temperature ever recorded on Earth (134ºF in Death Valley, CA on July 10, 2013). Figure 5.18 shows that this was indeed a good initial guess. The continuation method marches down to τ = 0.8550 before failing at τ = 0.8500. We thus assume this is the region for the turning point for the upper branch.

There is a rather peculiar result in our continuation method for the upper branch. As we mentioned above, we set the maximum number of iterations in Newton–Raphson to kmax = 20, and this indeed is a value that works. If we increase the number of iterations to kmax = 100 and let Newton–Raphson wander around the phase space this long near the turning point, we do eventually find a solution. The problem is that the solution is not physical! This highlights an important point about nonlinear equations governing physical processes; they may have many real solutions, so we need to use our knowledge of chemistry and physics to decide whether or not the solution is correct.

We are now ready for the hardest step of this turning point problem, namely finding a middle point to fill in the branch between the two turning points. Fortunately, our solutions for the upper and lower branch give us some very useful insights. First, we know that the two turning points are near τ = 0.85 and τ = 1.87. So we should pick an intermediate value, say τ = 1, that should be on the branch connecting the two turning points. Second, we know the values of the concentration and temperatures at these turning points. We would expect that the other steady state at τ = 1 should have a temperature and concentration intermediate between the upper branch and lower branch. We thus guessed x^{(0)}_1 = 0.4  and  x^{(0)}_2 = 1.08 and fortunately Newton–Raphson converged to the solution x_1 = 0.401 627  and  x_2 = 1.079 781.

We are now in great shape to complete the set of steady states. We can just march up and down in τ from the converged solution, using the continuation method until we reach the turning points we found in Fig. 5.18. Figure 5.19 shows the final result for the steady states of the CSTR. Below the turning point near τ = 0.85 there is a single steady state at a low temperature. Above the turning point near τ = 1.87 there is also a single steady state, but this time we have a hot temperature. The interesting behavior occurs between the two turning points. Here, for a given value of the residence time τ , there appear to be three steady states. Thus, it seems that the reactor could operate at three different temperatures for exactly the same residence time!

We should also try to figure out whether these steady states are stable or not; if the steady state is stable, then when we start up the reactor we know that it will not eventually reach that temperature and stay there. To make this calculation, we made a small modification to our continuation method program to also determine the stability:

The new parts of the program begin on line 11. We first compute the Jacobian from the steady state values of x returned from line 7, and we then call another function to determine whether this Jacobian corresponds to a stable steady state:

For convenience, this function starts by unpacking the entries of the Jacobian into the form

J=\begin{bmatrix} a& b \\ c& d \end{bmatrix}                                         (5.7.23)

We then check to see if the largest real part of the eigenvalues is positive. This is broken up into two steps. First, we see if the eigenvalues are imaginary. If so, then we just have to look at the trace of the matrix to find the real part. If the eigenvalues are real, then the term in the square root of Eq. (5.3.22) is also important. Notice that we only have to consider the positive term, since if the larger real eigenvalue is negative, then the smaller real eigenvalue also must be negative. The program cont_tau_stab returns a vector called stable that has a value of zero if the steady state is unstable and a value of 1 if the steady state is stable. Note that we are not considering the case of neutral stability. From a computational standpoint, round-off error makes it impossible to actually get a value of identically zero from the eigenvalue calculation in check_stability. (It is worth pointing out that MATLAB may return a value of zero because it “thinks” you wanted to work with integers instead of floating point.)

Figure 5.20 reports the results of this stability calculation. As we can see, only the upper and lower branches of the system are stable; the branch connecting the two turning point consists of a number of unstable steady states. Thus, for residence times between the two turning points, the reactor is either hot (the upper branch) or cold (the lower branch). Comparing Fig. 5.20 with Fig. 5.12, we see that this problem exhibits a hysteresis bifurcation.

Phase Plane
If there are two stable steady states between the turning points, how do we know which one is the observed state? After all, a well mixed reactor can only be operating at one temperature at a given point in time. The answer is that the steady state depends on the initial conditions. Some initial conditions will lead to the hot steady state, others will lead to the cold steady state. Since the problem is highly nonlinear, due to the exponential in the reaction rate, it is not obvious that a reactor that starts hot will stay hot and vice versa. Moreover, since there is an unstable steady state at a temperature between the two stable steady states, any trajectory that passes near the unstable state could be deflected.

To look at this more concretely, let’s consider the case where τ = 1. We already computed the unstable steady state at x_1 = 0.401 627  and  x_2 = 1.079 781 as the seed point to create the branch connecting the turning points. For this residence time, the cold steady state corresponds to x_1 = 0.952 803  and  x_2 = 1.006 293 and the hot steady state is x_1 = 0.127 103  and  x_2 = 1.116 383. These three steady states are seen in Fig. 5.20.

The phase plane for these parameters is in Fig. 5.21. We can clearly see the so-called “basin of attraction” for each of the stable steady states, namely the region of the phase space initial concentration and temperature that leads to a particular steady state. While it is true that lower initial temperatures tend to lead to the cold steady state, the dividing line between the two basins of attraction is not at a fixed temperature. Rather, the system selects the hot or cold steady state depending on both the initial temperature of the reactor and the amount of reactant inside the reactor when it starts. Recall that we are looking at an exothermic reaction. Thus, even if the reactor is somewhat cold when it starts, if there is a large amount of reactant inside at the start then the reaction rate will be somewhat fast (due to the concentration), which will increase the temperature. The increased temperature further increases the reaction rate, reducing the concentration. These two effects compete – if the temperature does not increase fast enough relative to the rate at which the reactants are depleted, then the reaction falls into the low temperature steady state. However, if the temperature can go up fast enough, then it can be maintained at a high temperature even as the reactant concentration approaches zero. The reaction must eventually slow down, since the reactant flow into the reactor is very slow compared to the rate of reaction at these temperatures, and eventually the system settles into the hot steady state.

The phase plane in Fig. 5.21 raises some important chemical engineering issues that are discussed at length in reaction engineering textbooks. First, we should note that the high temperature state corresponds to a high conversion. (Recall that the conversion is given by 1 − x_1, since we scaled concentration with the feed concentration.) Thus, if we want the reaction to be efficient, we need to use care to make sure that we actually end up operating at the higher temperature. This is a delicate control problem, since the way we set up the initial conditions here is artificial – it would be extremely difficult to impose a fixed concentration and temperature at t = 0 and then instantaneously start the reaction at t > 0 unless the reaction vessel is empty at the start of the process. Second, some of the initial conditions at high initial concentrations inside the reactor lead to very high maximum temperatures, much larger than the steady state temperature. The dynamics towards the steady state are important if, for example, the reaction temperature or pressure exceed the operating limits of the equipment. Moreover, very high temperatures could activate alternate reaction pathways with higher Gibbs free energies that we could safely ignore at a lower temperature. In short, understanding the steady state is obviously important, but ignoring the transients can literally kill you.

Limit Cycles
Depending on the various physical parameters, the non-isothermal CSTR can exhibit limit cycles of the type seen in Fig. 5.22. Qualitatively, you may be able to rationalize the existence of a limit cycle from the coupling between the reaction rate and the temperature. As the reactor gets hot, the reactants are consumed more quickly. While an exothermic reaction will continue to release heat, eventually there are very few reactants inside the CSTR to generate the heat. The reactants thus need to be replenished by the feed stream. If the feed stream is cold, the fluid flowing into the reactor will lower the temperature. (The temperature can also be lowered by having heat transfer at the boundaries, which we have not included in our model.) Eventually, there may be enough reactants to again send the temperature surging upward, and the cycle repeats again. While this logic sounds like it might produce a limit cycle, it is not necessarily a limit cycle or even cyclic behavior in the first place. For example, the trajectories in Fig. 5.21 seem to show part of the logic (a quick rise in temperature) but there is no cyclic behavior – the trajectories simply cool down to a steady (ignited) state. We could just as easily imagine that the scenario we just described would lead to chaotic behavior where the ups and downs of concentration and temperature are similar on each iteration but do not exactly overlap.

It turns out that Uppal, Ray, and Poore, in a classic paper in 1974, discovered limit cycles by numerical analysis of the non-isothermal CSTR when they also included heat transfer at the surface of the reactor. This additional heat transfer mechanism requires modifying Eq. (5.7.3) to include an addditional term

\frac{d}{dt}(ρC_pVT) = (ρC_pT_F)q − (ρC_pT)q − k_0e^{−E/RT}C_AV ΔH_{rxn} − hA(T − T_c)                                          (5.7.24)

where h is the heat transfer coefficient over the area A of the reactor and T_c is the temperature of the cooling fluid (where the notation “cooling” assumes that T_c < T). When we convert this to dimensionless form, we have

\frac{dx_2}{dθ} = 1 − x_2 + Da γ x_1e^{−β/x_2} − U(x_2 − x_c)                              (5.7.25)

where the dimensionless heat transfer coefficient is

U =\frac{hA}{ρC_pq}                               (5.7.26)

and the cooling fluid temperature is measured relative to the inlet fluid

x_c = T_c/T_f                                        (5.7.27)

Figure 5.22 shows one particular limit cycle for the non-isothermal CSTR corresponding to the case where the cooling fluid is at the same temperature of the feed temperature. We will give you the chance to explore limit cycles in the non-isothermal CSTR in more detail in Problem 5.20. The analysis by Uppal, Ray, and Poore, which includes the conditions we used to produce Fig. 5.22, was a major advance in reaction engineering.