Ind. Eng. Chem. Process Des. Dev. 1081, 20, 1-1 1
1
REVIEW Chemical Reaction Invariants and Variants and Their Use in Reactor Modeling, Simulation, and Control Kurt V. Waller' and Perttl M. Maklla Process Control Laboratory, Department of Chemical Engineering, A bo Akademi, 20500 A BO,Finland
In numerical analysis of complicated systems transformation of variables is frequently a useful tool to simplify the treatment. This is indeed the case for chemical reaction systems, where the primary state variables naturally transform to reaction invariants and reaction variants, as previously shown by Asbjornsen and co-workers. Nonlinearities and possible instability caused by the reactions are completely described in the space expressed by the variants. Reaction invariants have their origin in such basic laws as those expressing conservation of atoms, charge, and energy, and they are fundamental, whereas there is great freedom in the choice of reaction variants. In contrast to previous work in the field, the present paper advocates that the variants generally should be chosen nonorthogonal to the invariants. By choosing nonorthogonal variants, e.g., the number of state variables that have to be measured or estimated in multivariable control systems can be kept at a minimum. The utilization of the concepts is illustrated through applications. In a simulation of an industrial cement kiln the computational burden is significantly reduced through the transformation. The dynamics of acid-base systems can be treated in a linear reaction invariant space, all nonlinearities being included in an output relation and possibly the flow pattern description. In control of chemical reactors the dimensionality of the problem c a n be considerably reduced by working in the space of the reaction variants only. On the other hand, the space of reaction invariants is the natural starting point for systematic design of pH-control systems.
Introduction For obvious reasons, variables in physical systems are not often defined so as to facilitate numerical treatment of the system. Transformation of variables is then frequently an efficient means to make the system analysis easier, and it may even be the key to making a quantitative treatment possible. Chemical reaction systems lend themselves extremely well to variable transformations. Here the methods of linear algebra are especially useful, as pointed out previously by Aris and Mah (1963). Beek's review (1962) on the design of packed bed reactors is among the first papers containing material along these lines. A more general transformation approach can be based on the observation that the basic state variables, such as species concentrations, transform naturally into reaction invariants and reaction variants. This has been pointed out by Asbjornsen and co-workers in a series of papers and texts on reactor dynamics and control (Asbjornsen and Fjeld, 1970; Fjeld, 1971; Asbjornsen, 1972a,b, 1974; Fjeld et al., 1974). In these texts systematic methods were developed using only linear algebra for the transformations. The results are directed toward handling large chemical systems on computers. Although the methods mentioned are used in some simulation programs (Asbjornsen, 1975) for investigation of large industrial reactor systems (Villadsen et al., 1976; Asbjornsen, 1978),not much Use of the concepts is reported on in the literature. The papers are not even mentioned in a recent extensive review of chemical reactor theory (Lapidus and Amundson, 1977). A partial explanation for this might be that the treatment in the mentioned papers 0196-4305/81/ 1120-0001$0 1.OO/0
mainly has had its starting point in the differential equations describing the reactor and that the concepts have been illustrated largely by analysis of local stability. In addition to giving a review of the various ways suggested in the literature to decompose the state vector into invariants and variants, the present paper stresses the fact that the invariants derive from such fundamental principles as the conservation of atoms, charge, energy, etc., and have been used by chemists and chemical engineers for a long time in one way or another in simple systems. In this respect the paper is tutorial and complementary to previous papers. Further, it is advocated that the variants should be chosen on other grounds than orthogonality. Some applications illustrate the advantage of using variants which are nonorthogonal to the invariants, as well as of the concepts in general, with special emphasis on control. Reaction Invariants-Chemical Background Throughout the paper, the term chemical reaction invariants will be used for (state) variables that are not affected by the chemical reactions taking place in the system. For example, in the system A B, the variable CA + CB is unaffected by the reaction, or equivalently, the reaction is unobservable in the variable CA + cB. This (state) variable, CA + cB, is then a reaction invariant in the system mentioned. It should be emphasized that the invariants so defined need not be and usually are not invariant with respect to time or space. One way of expressing the stoichiometry of the system is actually through the chemical reaction invariants. Knowledge of the invariants means that the stoichiometry can be directly calculated as well as knowledge of the
-
0 1980 American Chemical Society
2
Ind. Eng. Chem. Process Des. Dev., Vol. 20, No. 1, 1981
stoichiometry means that the invariants are readily computed. The fact that the invariants express chemically trivial properties does not mean that the introduction of the invariants in the quantitative reactor treatment would not introduce new and valuable tools in the analysis. What makes the invariants a most valuable tool is that they are defined in such a way that they utilize the constraints expressed by the stoichiometry in a very efficient way for numerical computations. However, behind both the invariants and the stoichiometry there are well-known conservation laws governing all chemical reactions and restricting the number of independent chemical reactions possible in a system. Examples used in the sequel are the laws of conservation of the number of atoms of each chemical element in the system and the electric charge in the system. To summarize, the chemical reaction invariants, the stoichiometry, and the conservation laws are three various ways of expressing the same phenomena, in this case constraints on the system. The relations between the three concepts are quantitatively treated in the sequel. Consider a system with a set (mjp= 1, ...,p ) of distinguishable chemical species. Let (a$ = 1, ..., q ) be the set of all (atomic) elements in the compositions of ml,..., mp. Then define a matrix A with the following entries: A, = number of atoms of the element ai in the species mi;Ap+lj = electric charge of the species mi. Since the matrix A so defined contains information about the conservation of atoms and charge in the system, it is suitably called an invariance matrix. The chemical reactions in the system can be compactly written as
NTm e o (1) where 1\7r is a stoichiometric matrix and m = [ml,..., mpIT is the vector of chemical symbols. Since the generation of chemical elements and of electric charge in any chemical reaction is zero, it follows that AN=O (2) See, e.g., Aris and Mah (1963) for a more extensive discussion. Thus, the rows of the invariance matrix are orthogonal to the columns of N or the rows of the stoichiometric matrix NT. One consequence of this is that rank(N) 5 dim(m) - rank(A) (3) i.e., at most [dim(m) - rank(A)] linearly independent reactions can be observed in the system. The invariance matrix A contains the information necessary for writing down the possible chemical reactions in the system. However, to make the reaction formulas acceptable to the chemist, linear transformations may have to be made, as illustrated by Aris (1970). The invariance matrix A as defined above is made up of information about invariants in the system and may be directly used for calculation of chemical reaction invariants, as previously pointed out by Asbjornsen (1976a,b). The method is chemically most appealing but has the following shortcoming. The difference between the number of chemical species in the system and the rank of A gives the maximum number of (linearly) independent chemical reactions possible in the system, as expressed by eq 3. Which reactions actually take place in a specific system is, however, an empirical fact not revealed by the basic conservation laws. Each system can have specific constraints which restrict the number of independent reactions actually
observed in the system. The resulting invariants may be due to reaction passivation expressed by inert molecular groups or caused by the reaction mechanism itself without the existence of inert groups (Bjornbom, 1975). Reaction invariants can also be caused by kinetic constraints imposing relations between reaction rates. Thus, it is known that enantiomorphs of an asymmetric molecule (optical isomers) are produced with equal rates from symmetric molecules in a symmetric environment, i.e., by ordinary synthetic reactions. To conclude, some of the chemical reaction invariants of a system are directly given by the invariance matrix A expressing the conservation of atoms and charge. These invariants may or may not be all the reaction invariants in the system. The actual number of invariants is an empirical quantity and the only general way to obtain all the invariants is to start from the empirical facts of the system, expressed, e.g., as the stoichiometry. To put it another way: the chemical reaction formulas (1) express the conservation of atoms, charge, inert molecular groups and other constraints on the system. Assume now that the identification of the invariants in the system by atoms, charge, inert groups etc. has been successful. The invariance matrix A can then be extended so as to include the extra invariants mentioned above. If NT now denotes the matrix expressing the observed stoichiometry, then also in this case AN = 0 holds. A can now be called an extended invariance matrix. It is easy to understand that rank(A) = dim(m) - rank(N), i.e., that the total number of possible reaction invariants is given by the difference between the number of chemical species in the system and the number of linearly independent chemical reactions. Two simple examples illustrate the general concepts discussed above. In the first example there are six chemical species, and the vector of chemical symbols is given by m = [CH3CH=CHz, CHJHO, CH20,03, HzO, H2O2lT The invariance matrix (with al = C, a2 = H, a3 = 0)is then given by eq 3a. In this case the conservation of charge A=
3 6 0 0
[
2 4 1 0
1 2 1 0
1
0 0 0 carbon 0 2 2 hydrogen 3 1 2 oxygen 0 0 0 charge
(3a)
is trivial as the six species are all uncharged. The rank of A is three. Then, according to eq 3, rank(N) can be at most three; i.e., at most three independent reactions can take place between the six species. However, the actual rank of N observed to be sufficient in describing the system depends on the reaction conditions. Thus, e.g., the theoretically possible reaction CH3CH=CH2 + O3 ~ t 3CH20 . can be passivated, i.e., not observed under the actual reaction conditions. The conditions can be chosen so that the following single overall reaction between the six chemical species is the only significant reaction observed CH3CH=CHz + O3 + HzO CH3CHO + CHBO + HzOz In this case the molecular groups CH3 and CH2 remain inert in the reaction. If the invariance matrix is extended with these inert groups, i.e., the following two rows are added to A
Ind. Eng. Chem. Process Des. Dev., Vol. 20, No. 1, I981 3
w expresses the reaction invariants if and only if DN=O
(5)
then the rank of the so-extended invariance matrix becomes five, expressing the fact that only one overall reaction is observed between the six species. The second example contains three chemical species, namely m = [C12,C1+, C1-IT
i.e., the rows of D are orthogonal to the rows of the stoichiometric matrix NT. The vectors w and v are now given by w = Dx (6)
Now there is only one chemical element present. Looking only a t the species, this would allow two linearly independent reactions between the three species. However, because the electric charge must also be conserved, only one reaction is possible. This is seen from the invariance matrix
with
which has rank two. In the latter example the two variables w1 = 2cCl2+ cc1+ + ccl- (concentration of element chlorine)
w2 = ccl+
- eel-
(charge density)
are unaffected by any chemical reaction in the system, and no information of reaction rates is gained through measuring the variables w1and w2. As stated in the first paragraph of this chapter, such variables are called chemical reaction invariants. Variables which are affected by a chemical reaction are then naturally called reaction uariants. In the latter example ccI,, ccl+, and ccl- are all reaction variants. The reaction invariants are state variables which have an appealingly simple chemical behavior. It is natural to utilize this fact by employing, instead of the variables cc12,ccI+,and ccI-, the set wl,w2,and, e.g. (as the only reaction variant needed), cCl2as state variables in the analysis of the system. The original state variables can then be reconstructed from the new ones with a linear transformation (and vice versa). Matrix methods make possible a systematic treatment of these transformations. These computational methods are considered in the following section. Choice a n d Computation of Reaction Invariants a n d Variants The transformation of the variables describing a chemical reactor into a reaction invariant part and a reaction variant part can be accomplished by a linear transformation of variables. Let x denote the state vector or, more precisely, the vector of thermodynamic state variables for the reactor. It should be noted that the treatment to follow also is true in the case when x contains temperature and/or some other variable(s) related to the energy balance(s) for the reactor if the stoichiometric matrix is extended correspondingly. Thus, the whole state vector can be manipulated with essentially the same matrix methods. This will be illustrated later on under the heading Generalizations. The aimed linear transformation x = Pw
[:I
+ T v E [P,T]
-
(4)
has the property that [P,T] is a nonsingular square matrix, i.e., [P,T]-' exists. Furthermore, if we denote D [P,T]-l =
[
1
v = Lx
(7)
dim(w) = dim(x) - rank(N) dim(v) = rank(N)
(8)
The vector [wT,vTITthus contains the transformed state variables, w being the reaction invariant part and v the variant part. For the numerical computation of the invariants and variants a number of methods can be used. The essential point of the transformation (4) is the orthogonality property (5) of the matrices D and N, and any method to compute a matrix of appropriate size orthogonal to N can be used to obtain the invariants. Asbjornsen (1974) has given a matrix partitioning procedure to find the reaction invariants w. (For simplicity, we assume in the sequal full rank of N i.e., all the chemical reactions treated are linearly independent. There is no difficulty to extend the treatment to the case when N does not have full rank-several of the references do-but it is probably equally simple to reduce the equation system first so that all reactions considered are (linearly) independent before the reaction invariants and variants are computed.) The chemical species are arranged so that N can be partitioned as
.=[?I
where N1 is a nonsingular square matrix. Then D in eq 6 is given by
D = [-NZN1-',1]
(10)
which is seen to be orthogonal to N and to have full rank. Fjeld et al. (1974) have used a transformation given by x=[+Nv (11) where t contains the reaction invariants and v, as before, are the variants. This transformation does not, however, fulfill the conditions defining the transformation (4) because the number of transformed state variables (t and v) is larger than the number of original state variables (x). The transformation (11) is obtained from (4) by putting 5 = Pw andT = N Transformation (11) has the following inverse relations [ = (I - NN']x
(12)
v = Ntx (13) where Nf = (NW)-lNT is the pseudo-inverse of N. The matrix (I- NN') does not have full rank and thus contains redundant information, as also pointed out by Fjeld et al. (1974). Since the empirical facts of the reaction system-usually expressed through the stoichiometric matrix N-must be the starting point for computation of the invariants, and the invariants are orthogonal to the columns of N,every general method to compute the invariants starts from N.
4
Ind. Eng. Chem. Process Des. Dev., Vol. 20, No. 1, 1981
Furthermore, any orthogonalization method can be used, not only the two mentioned above. For example, in simple systems a basis set for the invariants can be picked just be inspection of N as illustrated, e.g., by Hammarstrom (1979). The other principal way to obtain the reaction invariants is to use the invariance matrix A. As explained, this method, however, does not always give all the invariants. On the other hand, it is a most simple and illustrative approach. Asbjornsen (1976a,b) has suggested the use of the invariance matrix in the following way. The decomposition of the state vector according to eq 4 is done as
X = ATw + Nv Then D and L in eq 6 and 7 are given by D = (AT)+
species-or another variable desirable for the use as a variant-increases the rank of the matrix [DTiLT]. The procedure goes on until the proper number of variants have been found, Le., rank([DTLT]) = dim(x). A remark on the computation of [P,T] is appropriate. If the variants are directly chosen as elements of x, which may be suitable, e.g., in control applications, then by proper numbering of the elements in x, eq 7 can be written v = [I,O]x (21) If [D/L] is partitioned according to eq 21a where now D2
(14) is a nonsingular square matrix, [P,T] according to eq 20 is given by eq 22. Thus, only a matrix of size (dim(x),(15)
and
L = Nt
(16) In the cases when the invariance matrix does not have full rank, the matrix A in eq 14 represents only a full rank portion of the complete invariance matrix. To the above suggestions from the literature, the following can be added. D in eq 6 is simply chosen as the invariance matrix D=A (17) This choice was used in the last illustration example in the preceding section. Further, a simple choice for the variants is
rank(N)) has to be inverted to compute the matrix [P,T]. An illustration of the various approaches discussed is given by the following example. A Simple Example for Illustration. Consider the following reaction system in liquid phase A + B e R + S rl
r2 (234 where rl and r2 denote the reaction rates. Writing the system in the standard form (1)gives eq 23b. At constant A+RieS+D
L = NT (18) Then the state vector x is reconstructed from the invariants w and the variants v according to
+
x = AT(AAT)-lw N(NTN)-'v
(19) To conclude the discussion of the invariants, it should be noted that the reaction invariants in fact define a unique linear space. Thus all methods give invariants which dif€er only by linear transformations. As stated above, the invariants are always orthogonal to the rows of the stoichiometric matrix, i.e., DN = 0. The variants, on the other hand, can be more freely chosen. In the methods described above to compute the variants they were chosen to be orthogonal to the invariants. However, this choice of variants, in the sequel called orthogonal reaction uariants, is unnecessarily unflexible. There is actually great freedom in the choice of reaction variants, and we suggest that the variants are chosen on other grounds than orthogonality. Combining eq 4,6, and 7, it is seen that, [P,T] is given by
[P,T] =
[
Any choice of reaction variants making the inverse of the right hand side of eq 20 exist is permissible. For example, as possible reaction variants, v can be chosen directly to be the state variables of a subset (with rank(N) elements) of the set (mjb = 1, ...,p ) of chemical species in the system. By a proper choice of the subset the inverse in eq 20 exists. In small systems a proper set of variants is usually found simply by inspection. In large systems, however, inspection becomes unreliable, and a proper set of variants has to be found by computation. This poses no principal difficulty: the computer tests whether the addition of a desired
temperature and pressure the state vector is given, e.g., by
x=
[CAY CB, CR, C S ,
%IT
where ci expresses the concentration of the chemical species miin the solvent. In order not to confuse the fundamentals of invariants and variants with a complicated flow pattern, consider a batch reactor. The dynamics of the reactor is described by the five equations (together with the appropriate initial conditions) in eq 24, where a dot above a symbol denotes C A = -rl - r2 CB
= -rl
CR
= rl - r2 or c = N r
Cs = r1
+ r2
CD = r2 (24) ita time derivative. It is seen that the effect of the matrix N is to spread out the two reaction rates rl and r2 to all five equations. Now we illustrate what happens in the analysis if we transform the state variables into reaction invariants and variants. Out of the various methods to calculate the invarianta we use the invariance matrix. Other approaches are illustrated later in connection with the applications. The number of linearly independent invarianta is given by dim(x) - rank(N), which in this example is three. It is convenient to denote A = XU,B = Z,R = XZ,S = Y, and D = X 2 Z . Then a possible choice for the invariance
Ind. Eng. Chem. Process Des. Dev., Vol. 20, No. 1, 1981 5
matrix is given by eq 25. Equations 6 and 17 then give A B R S D
[ : ::3:
A=
0 1
1 0
the invariants
O Y
w2
= CA
+ CR + 2cD + CS
w3
= CB
+ CR + CD
w1 = CA
(25)
(26)
Consider now the two variants. They can be chosen in a number of ways. Consider first the orthogonal variants, previously dominating in the literature. A simple choice is given by eq 18. Ula
= -CA - CB + CR + CS = -CA - CR + cs + CD
(274 The orthogonal variants suggested by Fjeld et al. (1974), eq 13, are Ulb = 1/15(-3CA - 4cB + 5cR + 3cs - CD) U2a
uzb
= '/15(-3CA + CB - 5cR + 3cs + 4cD)
system. They are consequences of such fundamental laws as those expressing the conservation of atoms and charge in a chemical reaction. Naturally, in this respect also such phenomena as phase changes, adsorption, desorption, etc. can be viewed as reactions in a broad sense, something touched upon also by Deans et al. (1970). It is then obvious that the invariant space is extended when use is made of other analogous conservation laws such as the conservation of energy. One consequence of this is that in an adiabatic continuous flow reactor the energy is reaction invariant. If the reactor is liquid phase, the enthalpy is (with good approximation) a reaction invariant. On the other hand, the approximation of invariant enthalpy may be quite inaccurate in a vapor phase reactor (Makila and Waller, 1980) where the difference between enthalpy and internal energy is larger. To illustrate how the invariant space is extended, the simple example of the last chapter can be used. The extension is made by including an energy balance in the analysis. For the batch reactor with liquid phase reactions the state vector then becomes, e.g.
(27b)
Now, if we consider other choices than the orthogonal variants, then the variants can be chosen simply as species concentrations, e.g.
and the reactor state equations are given by
L '
Actually, the variants can be chosen as any two species concentrations except the combination CA and CQ. This last combination would mean that the matrix on the right-hand side of eq 20 would be singular as also is evident from eq 24. It is also easy to see from the reaction system 23 or from eq 24 that knowledge of only CA and cs does not give information enough for computation of rl and r2. Now, the reactor can be described by the following equations instead of eq 24 W = 0
irla = 4r1 C2a = rl
+ r2
+ 4r2
(284
or W = 0
Lilb
= rl
ljZb =
r2
or W = 0
J
where the term U(Tj- 77 represents the heat transfer between the reactor jacket and the reaction mixture, and vTr expresses the effects of the reactions (exothermal or endothermal) on the temperature of the mixture. The matrix
can be considered an extended stoichiometric matrix. The reaction invariants are also now given by eq 6 w = Dx where DN, = 0 and rank(D) = dim(x) - rank(N,). Also, in this case, the invariants for the batch reactor follow the equation W = 0
i.e., they are constant throughout the reactor operation time, thus also invariant with respect to time. In the case of adiabatic operation, i.e., U = 0, enthalpy is a reaction invariant. If the chemical species are arranged so that eq 9 and 10 can be used, the reaction invariants are given by eq 32. Thus, in addition to the reaction
ulc = -rl - r2 dzc = -rl
(28~) It is seen that the suggestion of Fjeld et al. (1974), eq 28b, means that in batch reactors the variants are the so-called extents of reaction. After the transformation of the reactor equations utilizing invariants and variants, the analysis of the reactor is made with only two equations inatead of the original five in eq 24. If they are of interest, the original state variables are then calculated by means of eq 4. Generalizations I t has been emphasized that the chemical reaction invariants are just one way of expressing constraints on the
invariants we have in the isothermal case, there is now also the invariant (expressing the reaction invariance of the enthalpy) wq = -VTN~-'C~ +T (33) The same matrix methods as used in the isothermal case can be used for the analysis of invariants and variants in the general case when the chemical reactions bring about temperature effects as well as concentration effects in the reactor.
6
Ind. Eng. Chem. Process Des. Dev., Vol. 20, No. 1. 1981
It is also obvious that the inuariant space is reduced by phenomena or actions analogous to chemical reactions. If, for example, the reactor is controlled by manipulation of the concentration of a species, in principle this change in concentration has the same effect as a change due to a chemical reaction. Analogously, cooling or heating of the reactor means that the energy balance is affected. This view was adopted by Asbjornsen (1974),who extended the reaction rate vector to a reaction and control rate vector. He also provides systematic procedures for computation of t h e reaction and control invariants. For the simple example discussed above, the enthalpy is not reaction invariant in the nonadiabatic case. This natural result is seen from the matrix Ne,eq 31, since if UZO rank(N,) = rank(N) + 1 This means that the number of reaction variants in the nonadiabatic case is one more than in the adiabatic case. The generalization of invariants and variants to include also such phenomena as control does not change the general principles and approaches previously discussed, only the details of the treatment. For example, also in the general case the variants can be chosen orthogonal or nonorthogonal to the invariants. Application 1: Simulation of an Industrial Cement Kiln Westerlund (1977) has employed the projection operator method given by Fjeld et al. (1974) and Asbjornsen (1974) to find a set of state variables which would make possible the numerical solution of a minimum number of partial differential equations in a simulation study of multivariable control strategies for an industrial dry cement kiln. The dynamic material balances describing the kiln were expressed through a plug flow model and were of the form (34)
- -
The reaction system given in eq 35a was used to describe CaC03 CaO + COz Si02 + 2Ca0 (CaO)2.SiOz 4Ca0 + Alz03+ Fez03 (Ca0)4.A1203-Fez03 3Ca0 + A1203 (Ca0)3.A1z03 (CaO)2.SiOz+ CaO (Ca0)3.Si02 (35a) the formation of clinker in the kiln. This reaction system is compactly written in the form (1)G m o where the stoichiometric matrix G and the vector of chemical symbols m are given by eq 35b. Since the matrix NT has full -+
-
E l l- 2 0 1 -0 1 10
NT= 0
-4 0 0 - 3 0 0 - 1 0 0
0 0
0
0
o o0 o 0l 0
-1 --1 1 0 0 - 1 0 0 1 0 -10 0 0 0 1
rank, the suggested reactions are (linearly) independent. To transform the material balances into a reaction invariant part and a variant part, Westerlund (1977) used eq 12 and 13. The material balances for the reaction invariants are then given by
with the boundary conditions at z = 0 tf = (I - NNt)cf
(37) If tf can be considered to be constant, then according to (36) I is constant over the whole reactor length. Due to effective control of the concentration of the feed to the cement kiln, the boundary conditions at z = 0 were considered constant in Westerlund (1977). Hence the material balances (34), a system of 10 partial differential equations (PDEs), could be reduced to a set of five PDEs describing the reaction variants
The kiln model is completed by combining the material balances (37) and (38) with the energy balances for the kiln system (Westerlund, 1977). Thus, the computation time for the analysis could be significantly reduced by taking into account the special properties of the reaction invariants in the system in question. In the work of Westerlund (1977), the kinetics for the system (35a) was taken to be pseudolinear and of the form r = Kc =
h s , &F,
[ ~ C C ~ C O ~ ,
b A ,
k c d T (39a)
where the following abbreviations have been used: C' = CaO, S = SiOz,A = Alz03,F = Fe203. The strongly temperature-dependent k values were estimated from steady-state kiln operation data. By use of eq 11, eq 39a becomes r = Kt + KNv (39b) Inserting eq 39 into eq 38 gives the form of eq 40 for the a"
-av_ - - F - + at
az
r
k , ( t , - VI) k z ( b - uz) k3(t7
-
u3)
k,(t,
-
u, -
k,(E2
+
U, -
u4)
]
(40)
2 ~ -, 4 ~ , - 3 ~ , us)
material balances, now very appealing for numerical computation. Transformation back into the species concentrations is then simply made by use of eq 11 c=t+Nv (11') Thus, for the case when the invariants t are constant also with respect to time and space changes in the species concentrations are very easily computed from the changes in the variants. The example is well suited for an illustration of different methods of decomposing the concentration vector into invariants and variants. The invariants and the variants v were given by Westerlund (1977) as eq 40a and 40b. The symmetric matrix PD has rank = dim(m)- rank(N) = 10 - 5 = 5; i.e., only five of the ten componentxi in 5 are linearly independent. As seen, this decomposition results in quite complex expressions for both the invariants and the variants. The decomposition in which an invariance matrix A is constructed using the conservation of atoms, charge, and other inert groups is chemically more appealing. (This decomposition is, however, probably not numerically more
Ind. Eng. Chem. Process Des. Dev., Vol. 20, No. 1, 1981 7
-
= PDc=
4 8 -4 -13 3 -15 -4 -8 4 8 4 8 5 10
48 -4 52 6 -1 8 4 -4 -4 -5
52 4 48 -6 1 -8
-
-6 -13 6 54 29 25 13 -13 -13 16
1 3 -1 29 34 -5 -3 3 3 37
-8 -15 8 25 -5
70 -5 5 25 -20
-4 -8 4 13 -3 -5 68 32 -28 -10
4 4 8 8 -4 -4 -13 -13 3 3 5 -5 32 -28 68 28 28 48 10 10
-.
5 10 -5 16 37 -20 -10 10 10 47 -
av at
-= -F-
-6 -46 13 13 -17
48 6 4 4 5
1 29 -3 -3 -37
-8 25 -5 -25 20
-4 13 -32 28 10
4 -13 32 -28 -10
4 -13 -28 52 -10
1
5 17 -10 c -10 53 (40b)
efficient in this example than the method described above.) Using the invariance matrix, the reaction invariants are given by w = Ac. In this example the matrix A, taken with only linearly independent rows, can be conveniently chosen as shown in eq 40c. A is the atom matrix for Ca,
C, Si, Al, and Fe in the state vector. The reaction invariants are then simply given by eq 41, where the same abW =
-cCaCO,
+
CC'
+
2CC'?S
%',AF
+ 3CC',A
CCaCO, + CCO,
cs + CC',S +
CC',S
CA + CC',AF
+
CC'~A
CF + W,AF
1
3cC',S
(41)
breviations as above have been used. In addition to giving simple expressions for the reaction invariants, the invariance matrix in this case gives as extra information that the maximum number of linearly independent reactions allowed by the conservation of chemical elements is equal to five (= number of chemical species minus number of linearly independent elementary atomic and charge invariants). Since also rank(N) = 5 for the suggested reaction system (35), the stoichiometry of the system is sufficient to describe all theoretically possible (linearly) independent reactions between the ten species. In this method no matrix inversion is needed to construct the reaction invariants. Simple expressions for orthogonal variants would be obtained by use of v = NTc. However, in order to reconstruct c from w = Ac and v = Wc, matrix inversion would be needed as shown by eq 19. However, considering the assumed reaction kinetics for the system in question, eq 39a, a most natural choice for the variants, in this case nonorthogonal to the invariants, would be v=
[cCaCO3,
CS, CF, CA, CC'lT
+
av LNr = -Faz
+
LN
av -F - + az
10-2
4 -13 -8 -8 -10
av
az
C
(40a)
v = Lc =
on the material balances (34) gives eq 43, which is very
(42)
Then, letting the operator L according to v = Lc operate
similar to eq 40 with the basic difference, however, that no invariants appear in eq 43. In this case this advantage is not great, since the invariants were assumed to be constant in time and space. In a more general case the advantage is greater. Application 2: Dynamic Models for Acid-Base Systems The literature shows several examples of difficulties to develop dynamic models for acid-base systems. When acids and bases react in water solution, the system is a chemical reactor. However, it is usually assumed that the system is in equilibrium. This means that the system model can be obtained from a reactor model by letting the reaction rates approach infinity. The difficulties encountered in such a derivation can be circumvented by a direct elimination of the reaction rates. Some attempts in this direction have been made for simple systems (McAvoy et al., 1972; Orava and Niemi, 1974),but no general treatment has been found in the literature. However, already the definition of chemical reaction invariants indicates that this concept will solve the problem in a straightforward way, and as will be shown, this approach does not merely result in a dynamic model, but the systematic treatment reveals the basic features of the process. This gives a natural starting point for a systematic design of control systems. The presentation below is largely based on the work of Gustafsson (1978) and Gustafsson and Waller (1979). Consider an acid-base reaction process involving water solutions of strong acids and bases, 1-, 2-, and 3-basic weak acids as well as corresponding bases and ampholytes. The stoichiometry may be described by the reactions H20
* H++ OH-
a e h + b b s h + d d*h+e
(44)
where a, b, d, e, and h are vectors with chemical symbols as elements. a denotes undissociated weak acids, b, d, and e denote 1-,2-, and 3-acidic bases or ampholytes, respectively. Further, hT = [H+,Hf, ..., H+]. The dimension of all symbol vectors is equal, dim(a) = a. If not all weak acids in reaction system (44) are 3-basic, the symbol vectors e and d include purely hypothetical species. This is no limitation for the description (44), because the dissociation constants of the conjugate acids may be chosen to be zero. Thus the concentrations of such hypothetical species will always be zero. The choice of equal dimensions for all symbol vectors in (44)is made only to reduce the notation in the general treatment of stoichiometry and reaction invariance in the sequel. For a specific case there is no need to construct such hypothetical species since the derivation of reaction invariants in the
8
Ind. Eng. Chem. Process Des. Dev., Vol. 20, No. 1, 1981
needed in modeling the pH process when the relative change in the H 2 0 concentration is small. Then it is convenient to consider the following set of chemical species {OH-,H+,ai,bi, di,ei;i = 1, ...,a]. The reaction H+ + OH+ H20is then written as H+ + OH- F? 0. The zero follows because the product, H20, of the reaction has been excluded from the set of chemical species given above. Even when an incomplete set of chemical species is employed, the reaction system can be written in the general matrix form W m + o and the reaction invariants of the system can be found through Asbjornsen’s matrix partitioning procedure (1974). Generally the method of forming an invariance matrix to find the reaction invariants can be used only on a complete set of chemical species. However, the reaction invariants for an incomplete system can be found from those of the corresponding complete system. In other words, the reaction invariants for an incomplete system form a subset of the reaction invariants for the complete system. The invariants for the incomplete system are obtained simply by choosing such linear combinations of the reaction invariants of the complete system that in their expressions do not contain any of the concentrations of the species belonging to the complete set but not to the incomplete system. Thus, if in the acid-base system (44), H 2 0 is not considered, the conservation of oxygen (or hydrogen) cannot be employed. However, the conservation of charge and the total concentrations of the weak acids apply. Consequently, the reaction invariants of the in; they are complete system are given by w1 and w , + ~i.e., one less than in the complete system. All the reactions of the system (44) essentially involve a proton transfer. Many ionic reactions of this kind have a very fast relaxation to the equilibrium. If it is assumed that all the reactions in (44) have this property, it can be assumed that the reactions are in equilibrium at any point and every moment in the system. The Gibbs phase rule (Denbigh, 1964) then gives the number of degrees of freedom at any point in the reactor for the acid-base system (44) to be a + 3. If, furthermore, it is assumed that temperature and pressure are constant, there remain a + 1 degrees of freedom at any point of the reactor at any time. Thus, the linear space of reaction invariants gives the chemically natural and elegant state variables w1 and w,+~to model the acid-base system for any flow pattern in the chemical reactor. The following examples (Gustafsson and Waller, 1979) of the invariants w1and wd1 of (47) provide an illustration. (1) For strong acids and bases mT = [OH-, H+]
sequel is valid for all combinations of 1-, 2-, and 3-basic weak acids. The vector of chemical symbols m should be chosen so as to give a stoichiometric matrix of appropriate structure. Such a choice is m = [OH-, bT, dT, eT,aT,H+, H201T
(45)
If the reaction system (44) is written in the matrix form (l),NTm i=? 0 , the matrix N is given by eq 46. Now
r1
o~l
o~
o~
dim(m) = 4a + 3, dim(N) = (4a+ 3)-(3a + l),rank(N) = 3a + 1. The matrix N can be partitioned as shown above so that N1 is a nonsingular square matrix. The linear space of reaction invariants has the dimension dim(m) - rank(N) = a + 2 and, e.g., Asbjornsen’s matrix partitioning method (1974) can be used to give a set of a + 2 linearly independent variables w. These are obtained by use of eq 10.
w = [-N2N1-11 I]c Carrying out the matrix operations gives w1 = c , + c b + c d + c, a
wm+1
=
CH
-
COH
-
r=l wa+Z
=
CHzO
a
ccb,i
-
(10’)
a
2ccd,i - 3 c c , i i=l
i=l
+ COH
(47)
As previously stated, a chemically more appealing way to obtain the reaction invariants is to construct an invariance matrix. From the reaction system (44), a natural choice of inert chemical groups is the atomic groups which are the elements of the vector e. These atomic groups are the 3-acidic bases, which are assumed not to yield further protons. The atomic groups of the vectors d, b, and a are obtained from the corresponding bases of group e through addition of one, two, and three protons, respectively. Thus e is a vector of transparent groups which are contained in the corresponding elements of d, b, and a in the ratio 1:l:l. Combining this with the conservation of oxygen and charge in the chemical reactions taking place gives the invariance matrix (eq 48) where k; = number of oxygen atoms in ai, j ; = charge of a,, and iT = [l,1, ..., 11. The rows of A are linearly independent. (A possible row for hydrogen atoms is left out in A since such a row would be linearly dependent on the rows for oxygen atoms and electric charge.) This gives a set of a + 2 independent reaction invariants according to w = Ac. A simple linear transformation gives (47). Thus, the obtained reaction invariants have their chemical interpretation in the conservation of the total concentrations of weak acids, of oxygen atoms, and of electric charge in the reactions 44. It is not always convenient or necessary to use a complete set of chemical species to describe a reaction system. For the reaction system (44) the solvent, water, is not OHA=
w = CH - COH (2) For the system m T = [OH-, Ac-, HAC, H+] = C H A ~+ CAC 102
w1 w2
eT
I kT
I kT
I kT
iT)
(jT
- 2iT)
COH - CAc
=
cH2c03
+
‘Hc03
+
CO?-, HzCO3, H+] ‘COa
CH - COH - CHCO3 - 2cCO3
(4) For the system mT = [OH-, H2P04-,HC03-, HPOd2-, C032-,P043-, H3P04,H2C03, H+l
dT
-
-
(3) For the system mT = [OH-, HCO;,
bT
(jT
= CH
aT
H + H,O group e oxygen
Ind. Eng. Chem. Process Des. Dev., Vol. 20, No. 1, 1981 9 w1
=
CH3P04
w2 = W3
+ CHzP04 + ‘HPO4 + ‘Po4 +
C H ~ C O ~ CHCO,
+ CCO,
= CH - COH - CHzP04 - CHCO, - 2cHPo4 - 2CC03 - 3cP04
Now, to obtain a dynamic model for the acid-base system, the flow pattern has to be considered. For example, the dynamic model of a perfectly mixed CSTR, conveniently written in the standard form (Asbjornsen, 1974) 1 c = -(cf - c) + Nr 7
becomes for an acid-base system with very fast reactions
together with the appropriate initial conditions. Thus,the open-loop model of the acid-base CSTR is linear and decoupled, and the uncontrolled reactor is asymptotically stable. This is a very natural starting point for design of control systems, as illustrated below in Application 3. Orava and Niemi (1974) have derived a linear model for strong acids and bases (example 1 above) and McAvoy et al. (1972) a model for the system acetic acid and sodium hydroxide (example 2 above). Actually, the modeling approach they suggest for pH in CSTR‘s is essentially a special case of the reaction invariant method. The basic properties of these models have their natural interpretation in the properties of reaction invariants. Design of Control Systems for Chemical Reactors In several ways the reaction invariants and variants can be utilized to promote the design and subsequent operation of control systems for chemical reactors. One of the key properties of a control system is its stability. The primary causes for instability are the chemical reactions and the control actions. Therefore, for local stability it is often necessary to consider the reaction-or-control-variantpart of the system only. Consequently, in many systems it is natural to consider the reaction-and-control-invariantsas disturbances or fdters for disturbances and to compute control laws for the reaction-or-control-variant part of the system only. This approach has, e.g., the advantage of reducing the dimension of the problem. The approach has been suggested by Asbjornsen (1974) and illustrated by Hammarstrom (1979). Both have chosen the reaction-or-control-variant part of the system to be orthogonal to the reaction-and-controlinvariant part. However, the choice of orthogonality has the important drawback that it “would require all states to be available” (Asbjornsen, 1974). This fact was illustrated above, where eq 27a and 27b show that orthogonal variants imply that all species concentrations are measured or estimated. The disadvantage was discussed by Hammarstrom (1979), who proposed a method to modify the control laws, calculated for the control-or-reaction-variant part of the system, so that the modified suboptimal control laws could be realized using feedback of only as many states as there are control-or-reaction-variants. The disadvantage is eliminated also if the reaction-orcontrol-variants are chosen so that they are closer connected to the variables employed in the implementation of the control laws. The approach is illustrated above under the heading “a simple example for illustration”, where the relaxation of the orthogonality demand was demonstrated to reduce the measurement vector from all five components for orthogonal variants to two for nonorthogonal variants.
Thus, both the design and implementation of the control systems can be very significantly facilitated if the synthesis and implementation is confiied to the space of the variants only. This space contains the possible instability and also the possible nonlinearities caused by the chemical reactions. Thus, by use of the chemical reaction invariants the state space of main concern in control can be considerably reduced, something that facilitates the design, implementation, and operation of control systems. Another utilization of the concepts of reaction invariants and variants for control is illustrated in the following section. Application 3: Reaction Invariant Control of pH In Application 2 a model for acid-base systems was derived. This has already solved a problem which has caused some trouble: many pH-control systems have been designed without a proper model of the system. The model derived has a nice form: expressed by reaction invariants it consists of a number of decoupled subsystems. If the process is a CSTR, the subsystems are of first order. Let us now consider feedback control of the system. For the case of a CSTR with one feed inlet, one control stream, and one effluent stream, the system model takes the simple form dw Ff Fc Fe - = -wf + -w, - -w V V dt V In pH-control systems pH is usually measured, and thus a relation between pH and the state variables w of the system is needed. Using the equilibrium expressions c ~ c O H = K,; cH’cb,i = Ka,i-ca,i;etc., and the expression pH = -log [cH/(mol/L)]gives pH as an implicit function of the reaction invariants w,+1 = f(PH, w1) (51) The output, pH, is a strongly nonlinear function of the state (WI, W,+I). For economical and environmental reasons the control stream in most cases contains only a strong acid or base. The control vector is then simply
w, = [O, 0,
... 0, %+‘IT 9
and the control affects only the last partial system of eq 50. The dynamics of the control loop may then be presented by the scalar equation
Thus, for control purposes only the last equation of the model need be considered. The a first equations in (50) act only as first-order filters for disturbances. Through introduction of the control action, the reaction invariant w,+‘ becomes a reaction-or-control-variant(of course, it is still a reaction invariant). Thus also in this application the control system is designed in the space of reaction-or-control-variants. In a linear control loop F, is manipulated through a linear feedback of A W , + ~ .When the control stream contains only strong acids or bases, the reaction invariant variable wc,,+1 of this stream is wc,,+l = c,,H - c,,OH. Thus, if c c a>> c,,OH or c,B
