Here's the task for this cstr model: given the feed temperature(s) and a heat exchange amount, Qset, find the reactor outlet temperature.
This is a simple root finding problem with a scalar function and a scalar independent variable, T. The function, f, that we need to make zero is f = Q - Qset where Q is the heat exchange requirement at the current value of T. The cstr_T routine discussed in the previous post provides Q for the current T guess.
Root finding method
This problem is fairly easy because the operating curve that leads to a value of Q is smooth. There are several good methods for finding a scalar root: I chose the bisection method.
The first task in this method is to find two temperatures that bracket a solution. The first T value chosen produces a value of f that is either positive or negative. The direction to take for the next temperature is the negative of the sign of f. Additional steps are taken until the sign of f changes, meaning the root has been bracketed by the last two points.
After bracketing, the temperature difference is halved repeatedly until the tolerance is achieved. On each iteration, the two closer points that bracket the root are retained and one of the older points is discarded.
Guaranteed stable solutions
This routine will never obtain an unstable solution. This is due partially to the fact that the step direction is determined by the sign of f determined on the first guess of T. The figure below may help show how temperatures will move depending upon the initial T guess.
If a temperature is chosen such that it falls on the red curve in the white zone, the temperature will be increased until the sign of f changes. If the temperature falls on the curve in the yellow zone, the temperature will be decreased until the sign of f changes. These directions always lead to a stable root. Once the root is bracketed, there is no step direction, only step halving.
What happens if initial guess is at the unstable solution?
It is highly unlikely that the initial temperature guess would be equal to the exact root value, but it could be close enough to result in a value of f that would be within the desired tolerance, thereby returning an unstable solution. To prevent that, the routine requires that the bracketing search doesn't begin until an f value is found that does not satisfy the tolerance. Once that criterion is met, the step direction for the rest of the search is computed as described above. Thus, initial guesses at an unstable solution will result in "movement" to a stable solution. Initial guesses at a stable solution will first move away from the solution and then return to it.
How did I achieve solutions at unstable points earlier?
The single phase algorithm that I used prior to creating cstr_Q used the nonlinear equation solver in Mathcad. This method would usually find the closest solution to the starting guess, including an unstable solution. The search algorithm that Mathcad uses had no special information or directions that would avoid unstable solutions.
Will your process simulator give an unstable solution?
Maybe yes, maybe no. Check the documentation: even better, run some tests. Conduct the outlet temperature survey and select a Qset that has an unstable solution. Run the cstr model with that Qset and the T value for that solution as the starting guess. Does the program find another stable solution?
The bisection routine that I created in Mathcad was based on a FORTRAN version in Numerical Recipes, W.H. Press et al, Cambridge U. Press (1986). This book is a major source of programs for the most successful routines needed in a computing toolbox. The programs are given in FORTRAN and Pascal in my (First) edition, but there is a newer edition with C++ programs.
Another source is the web, Wikipedia in particular. Many of the root finding routines are described there in brief mathematical form, but you may need to look further if you want a coded example.