I've made progress on building a process simulator. (If you ask "why", then you must have missed the last post.) My initial process consists of a boiling CSTR and a reflux condenser. This process has one recycle stream, which I chose to be the liquid reflux from the partial condenser. I've successfully created the CSTR reactor, condenser, and written a loop that converges. In this post, I will focus on the multiphase CSTR model.
Two models required
I built the CSTR model in two parts. First I created a cstr_T routine that solves the multiphase mass balances for a given outlet temperature. Then I wrote a cstr_Q routine that finds the outlet temperature needed to satisfy the energy balance for a given heat exchange, Q. This routine calls the cstr_T routine in a loop until it converges on the specified Q.
I thought the cstr_T model would be easy because it didn't need to find the root of the heat balance curve and the heat exchange line. To my surprise, it proved to be the most difficult program to write. The difficulty was due to the need to iterate between a flash routine and the reaction mass balances. In my earlier work with a single phase, a Mathcad solve block was used to find the outlet composition. I could not incorporate the PPP flash routine into that solve block so I took the iteration route.
The cstr_T convergence method
Using the current composition guesses for the outlet vapor and liquid phases, the reactions are computed and the total outlet composition is determined. I then calculate the distance between the new and old total outlet composition points. If this distance is not below the desired tolerance, then the outlet stream with the new composition is flashed at T, P and the vapor and liquid compositions are determined. A damped substitution method is used to update these compositions for the next iteration.
I am not very satisfied with this current convergence method because the substitution must be extremely damped. By that I mean that the new phase compositions are not allowed to move very far from the previous guess. Most likely, I will return to this problem later and explore other convergence methods.
Some features of the model
I've allowed for an unlimited number of feed streams. The stream numbers for all of the feeds are supplied as a vector. This avoids having a fixed number of input streams. It also might be useful for exchanging mass between a cstr and a distributed phase such as a gas stream rising through the cstr liquid.
Reactions may take place in one or both phases in the reactor. Each reaction in the network uses either a liquid or a vapor composition in its rate expression. The appropriate reactor volume is used for reactions in each phase.
Reaction rate expressions can be of any form, and reversible reactions may be included. The complete flexibility is due to the fact that the rate expressions are written as Mathcad functions prior to their use in the reactor program. This avoids the practice in most simulators that requires the reaction rates to conform to certain expression forms allowed by the input data forms.
In anticipation that a process might require more than one reactor, I have set up the reaction rates and reaction stoichiometry in nested arrays. This allows me to create reaction sets, one for each reactor. The reaction set is then specified in the call statement for a given reactor.
The figure below demonstrates how the nested array for the stoichiometric matrix is created. In this example, the stoichiometric matrices for two networks have been created. These are manually input as the 0 and 1 elements of the overall matrix. A network can be selected by referring to the desired element of the overall matrix. The figure also shows that an overall matrix created with the Mathcad stack command doesn't produce the desired result.
In a similar manner, the reaction rate function, rset, is created from rate functions for each reactor.
The component production rates for a given reaction set, rsn, are computed with the Rxn function as defined below:
I'm using the hydrogenation of trimethyl pentenes to isooctane as a development example. There are two irreversible reactions and one reversible reaction in the network. In order to demonstrate cooling via boilup, I have selected n-pentane as the solvent. See June 17, 2013 post for more details about this reaction system.
The outlet temperature survey
Using cstr_T, I conducted a case study with the outlet temperature as the variable. The operation curve is shown below.
This plot is much like one I presented in an earlier post (Aug/7/2014). The earlier plot was for the same example, but vaporization was not considered. The unstable region in the earlier plot was between 330 and 360 K. With cooling provided by vaporization, the curve has moved down about 10 K so that the stable region is now between 315 and 350 K.
Next: the cstr_Q routine.