99
Mathematics and Computers in Simulation XXIII (1981) 99-106 North-Holland Publishing Company
A PHYSICO-CHEMICAL MODEL AND SIMULATION STIRRED TANK REACTORS *
OF pH-PROCESS IN CONTINUOUS
P. JUTILA, P.J. ORAVA and B. SALMELIN Control Engineering
Laboratory,
Helsinki University of Technology,
1. Introduction The acidity of aqueous solutions is quite an important quantity in many processes in nature as well as in industry. The acidity is customarily expressed by the pH-variable which gives the concentration of the hydrogen ions H+ in a logarithmic scale. The functioning of certain processes requires narrow ranges of pH-values. Examples of this kind of cases can be found in water and waste-water treatment processes. In addition, there exist reliable and simple instruments for the direct measurement of pH. Partly due to these facts, the observation and control of pH is quite an important area in concentration process engineering. The mixing processes in stirred tank reactions form the most common class of (concentration) pH-process systems. The control of pH-processes is often difficult originating in nonlinearities of their models as well as in the changes of these nonlinearities caused by changing composition of chemical species. For example, the control of waste-waters of some industrial processes, consisting of fast varying amounts of concentrated wastes, often is difficult to be realized by conventional methods. On the other hand the pH-processes in clean water treatment form a case of control problems in which the incoming material concentrations do not change very much and fast and the control task is relatively easy. Quite a comprehensive survey of the problems applications of the pH-control is presented by Trevathan [ 11. The dissociation processes, involving the pH-variable as a measure for the acidity, is a classical area in chemistry. However, the quantitative physico-chemical
* Presented at the SIMS-meeting May 19-25,
0378-4754/81/0000-0000/$02.50
1980
@North-Holland
SF-02150
Espoo 15, Finland
approach has only recently been applied to studies of control [2,3]. More lately, a practical linear feedback control algorithm based on the physico-chemical model was presented by Niemi and Jutila [4] and applied to the control of a continuous flow laboratory process. In the present paper, the above-mentioned model is extended by including hypothetical chemiLa1 species, strong and weak ones, both one-valued and twovalued ones. The dissociation constants for the weak species will be evaluated by a static fitting to the titration curves of real liquid samples in each application. Via a certain transformation, it is shown that the unknown concentrations of the hypothetical species can be estimated from the pH-measurements by the aid of a linear Kalman-filter algorithm. The abovementioned control algorithm can then be used also in connection with the extended process model in question. The study of the system and it’s parts by simulations is useful before practical implementations of control.
2. pH of Aqueous Solutions
Dissociation of the water (H, 0) molecule to the hydrogen (Hf) and hydroxyl (OH-) ions is a very fast and reversible process. The dependence of the ion concentrations on each other is unambiguously expressed by the value of the ion product K,. l The concentrations can be used instead of the activities with a good accuracy especially for dilute solutions. Only isothermal systems shall be studied, and the
’ A list of symbols can be found in the Appendix.
P. Jutila et al. / Physico-chemical
100
model
value of the ion product at 25°C is given by K,
= aH+ . aOH- = C+C- = lo-l4
.
(moles/litre)*
(1) Dissolved acids produce H+ ions and bases OH- ions by dissociation. The dissociation reactions of liquid phase species are very fast, too. Strong species dissociate fully and the weak ones partly. Since the liquid is electrically neutral, an equation, the electro-neutrality condition, which consists of the concentrations of all ions present must be valid at any point. For the special case of liquid containing strong acids and bases, one weak acid and base of single step of dissociation (i.e. one-valued) and one weak acid and base with two dissociation steps (i.e. two-valued) we obtain:
-2c,-,--c,-,=o.
(2)
The 1st and 2nd term of (2) are the above concentrations Cf and C-. A strong, i.e. fully dissociated base produces positive cations in a concentration which is equal to the total concentration of the base at the point under inspection (3rd term), and a strong acid negative anions in a concentration equal to the total concentration of the acid (4th term). A two-valued base produces cations of charge one and two (5th and 6th term). Of course the number of ions with charge two must be multiplied by a factor of 2. An one-valued weak base produces cations of charge one only (7th term). A two-valued weak acid produces anions of charge one and two (8th and 9th term). A 1. one-valued weak acid produces ions of charge one (10th term). Combining (2) with the dissociation equilibria of the weak compounds, Eq. (3), one by some effort obtains the final form of the electroneutrality condition, Eq. (4). ,
K,‘,
= 3
Kj2Ccl2
-c+
_
+K,‘,
-0.
(4)
The dissociation constants Ki,, Kfl, Kz2, KL1, Ki,, and concentrations C,, , C,, , C,;, C;r-, CQ2, G2 a29 q2, ($1, cbl? Cll, Cl:, Cp2, Cb2 and Cg2 are presented in (moles/litre). Extension to the case of chemicals with more dissociation steps is straightforward. It is worth mentioning already now, that (4) is linear in Cn - C,, Car, c
G2, C’,.rrand Cp2. The acidity of a solution is easily means of a pH meter. This produces is linear dependent on the pH value but the dependence of the latter on highly nonlinear. pH = -log,,
[a,+/(moles/litre)]
= -log,,
[C+/(moles/litre)]
3. Continuous
pa-reactor
measured by a voltage which of the solution, concentration is
.
(5)
and the control algorithm
For the complete model of pH-process in stirred tank we also need the equations for the tank dynamics. The process scheme is presented in Fig. 1. If
, a2
K&
_
“-5;
,
Ki2
Cbl c,,
=c,r+c,-,
+ca<-,
Cc& = Ca2 + ca< t
c,
-‘-‘ii
,
cb2
(3)
Process
hquid
C,, = C,,, + cb+r + cb’ii , Fig. 1. Stirred
=cb2 -+c& ;
Cc& cp1, Cpz,.
tank reactor.
Symbol
CM refers to vector
CC,,,
P. Jutila et al. / Physico-chemical
the mixing in the reactor is thought to be ideal and only strong control reagents are used, we obtain: - Cl3) =
Q,(cA, -
v&p
= QpCcYPp -
62,
-
cBp)
+ QC>(CA
(Qp
+ QC)
+ QC(cAC -
-
CBC)
CB> ,
C,lp .
(6)
Symbol Cd refers to vector (Colt, Cez, Cpl, Cpz) so that the latter equation includes four scalar differential equations. Eqs. (4) (5) and (6) form now a complete model of the pH-reactor in question. A better model for real mixing process is an ideal tanks in series model. A series of three tanks is shown in Fig. 2. The hydrogen ion concentrations in tanks and the tank volumes are symbolized as shown in Fig. 2. The control principle is now presented in Fig. 3. All the chemical species (control and process chemicals) are flowing through the mixing dynamics and they then determine the hydrogen ion concentration
Flow
%3P ->
101
according to the relation of (4). Let us now suppose that both the pH and all the weak concentrations are measured. The weak concentrations can in principle be measured by ion-selective electrodes. When we now know the weak concentrations Cap and the hydrogen ion concentration C+ (calculated from (5) we can compute the concentration (CA - Cn) from (4). Other possibility to the determination of (CA Cu) is the direct measurement of the concentrations of strong species. This alternative is included in Fig. 3. On the other hand when we know the pH-reference value pHref and the weak concentrations C,, we can calculate from (4) and (5) also the instantaneous value (CA - Cu)ref. (CA - Cu)ref is the difference of the amounts of strong components that are needed to be present in the process output to bring the output PH to pHr,r. When now the difference d = (CA - CB)& - (CA - C’n) is fed to a linear PID-controller which controls the concentration of the strong reagent chemical, the feedback control loop is accomplished. In practice it is difficult to control the reagent concentration and that’s why we control the reagent flow Qc. Because the control chemical concentration is much larger than that of the process chemical, we have Qc << Qn. If for instance the reagent chemical is a base we can approximately write Q&t) Ccs c
Fig. 2. Three ideal tanks in series.
I@*
model
v and
&Q
Fig.
> 2
, *
QP
C+
Ionic
mixing
equilib r,um Eq.
-
CL) -
-. LineOr
controi-
'AP - CBP -
d -*
Meoiu rement
CBC, oc
ler
pH ref
e and
Flow Fig.
2
G Meosu -
~on~ineor
eQb
rements
CO”“t?sions
Fig. 3. The control
loop scheme
with the proposed
control
algorithm.
102
P. Jutila et al. / Physico-chemical
QcSCB(t), where CcS and QCS refer to steady state values. In other words the control flow can be considered to be constant and the concentration to vary. The simulation of this system is going on for the time being. It has been found that the right root C+ of (4) can be found by Newton-method. If the number processing capability of the computer is not good enough, the concentration qualities must be scaled. As a simple case we study now the neutralization of one weak acid, namely acetic acid, with a strong base, potassium hydroxide. Now in (3) KJ,= 1.75 . lo-’ (moles/litre) and all other dissociation constants are zero. The electroneutrality condition can be given, after some algebra, in the polynomial form [4] : C+3 +(CB
+K:,)C+2+(K;,CB-K,C,, -K,)C+
-K;,K,=O.
(7)
Along to the above described algorithm we should now measure pH and the acetate-ion concentration C,, . In practice it is easier to measure the potassiumion concentration CB by ion-selective potassium electrode. When pH (and thus C’) and CB are known the acetate ion concentration C,, can be calculated from the algebraic equation (7), the electroneutrality condition. Thus the control loop remains equivalent. When we now use the dynamics model presented in Fig. 2, we can simulate the system behavior [5]. The presented algorithm performance must then be compared with that of the conventional feedback pH-control. In Fig. 4, simulation results are presented in a case, where the input pH of the process liquid (acetic acid) is changed stepwise from 4.0 to 4.2 and back to 4.0. A computer controlled continuous flow system
model
Fig. 5. Laboratory test results: (a) pH-PI control; (b) the presented algorithm.
has been constructed in the Control Engineering Laboratory, in order to test the pHcontro1 -algorithms. In Fig. 5 two runs with this test system are presented, corresponding to the simulation runs in Fig. 4. The test liquids were made from pure chemicals and ionchanged water. In both cases PI-controllers were used. From the responses it can be seen, that the conventional pH-control reacts slowly when pH deviation is large and tends to oscillate near the neutral point. The new method is better in this respect. The differences between simulation and test results exist because of nonideality of the process liquids (always some amount of carbon dioxide present) and of nonideility of the ion-selective electrode. It has also been stated by simulation, that the presented algorithm is very insensitive to additive constant measurement errors of the ion-selective electrode [6]. Further it can be stated both by simulation and analytically that the algorithm converges to the setpoint, when the controller is tuned so that system is stable. The algorithm suits also to cases, in which the reaction result, the salt, precipitates [7].
4. Model involving hypothetized
Fig. 4. Simulation results: (a) pH-PI control;(b) the presented algorithm.
species
In the previous section we thought, that the concentrations of all weak species of the process liquid could be measured. This is very seldom, if ever, possible. First, we do not usually know all the species present (e.g. the case of waste waters).! Second, the technology of continuous measurement of ions is inadequate. We do not namely have ion-selective electrodes for all ions and besides the existing electrodes are either too unstable (plastic membrane electrodes) or too slow and non-selective to other ions, especially to
P. Jutila et al. / Physico-chemical
hydrogen ions (glass electrodes). Because of these measuring troubles the full accuracy of the model is not of much importance. We can consequently measure only the pH of liquids and the concentrations of the weak components must be estimated. The model must be as simple as possible to guarantee reliable estimation and it must be able to describe the main features of pH-process. In other words the weak species of the model does not need to have accurate counterpart in process liquids and the model need not have accurate counterpart in process liquids and the model still preserves the physico-chemical character of the problem. This kind of thinking thus results in models involving hypothetical species. The dissociation constants of the weak species can be fitted by instantaneous, i.e. static, considerations. The fitting is based on the electroneutrality condition, (4), which determines the static relation for every time instant, the so called titration curve. The fitting of dissociation constants of the hypothetical weak species is now made to typical titration curves received from practical samples. The concentrations of the hypothetical species in on-line control are then estimated simultaneously with the control action for instance with an on-line filtering technique. The control algorithm will be of the same structure as presented in the preceeding section. In Fig. 6 two titration curves presented in literature are shown, [8]. The abscissa scale is changed, being now the concentration scale for 1-valued strong base (e.g. potassium hydroxide, KOH) instead of the flow scale of concentrated 2-valued almost strong
model
103
12 PH 10
8
6
Fig. 7. A model with parameter
values: as given in Table 1.
base (calcium hydroxide, Ca(OH),) in [8]. The titration curves in Fig. 7 are constructed to approximate those in Fig. 6 with a model which consists of one strong acid, two weak one-values bases and two weak one-valued acids. In other words dissociation constants Kil and Kzl of (3) and (4) are zero. The curves of Fig. 7 tit to those in Fig. 6 fairly good between pH-values 4 and 10, but in the neighborhood of the neutrality point they are somewhat too steep. This could be avoided by nonzero values of Kz, or/and Kzi. Of course such a more complex model has better fitting capabilities than a simpler one. In order to yield a linear estimation problem for the concentrations the dissociation constants must be a priori fitted. Simultaneous estimation of dissociation constants would namely lead to a complex nonlinear problem. Table 1 Parameter
values of the model presented
in Fig. 7.
Concentration
PH 10
Dissociation contants (moles/litre)2
8 6
0 12 3 x 2,7~10~3mole/l~tre
L
5 ‘B
Fig. 6. Two practical titration curves to which the fitting of dissociation parameters of the model presented in Fig. 7 is made.
Case b (moles/litre)
KW
lo-l4
Gl
6.3
G
10-4
CCQ
10-3
8
CP1
10-S
10-q
92
10-3
10-S
cA
1.56.
& 2
Case a (moles/like)
e2
G1 G G2
0 10-S 5.10-6 0 5.10-S
. 1O-3
1O-3
1.3.
7.7.
10-Z .10-4
10-3
P. Jutila et al. 1 Physico-chemical
104
5. The estimation of the hypothetized concentrations and the combination of estimation and control The use of the model involving hypothetical species for the control forms a complicated problem of simultaneous estimation and control. In this paper we suppose that the estimation of the concentration variables can be accomplished under the normal control operation by assuming the control action as a known variable. In other words, we suppose that a kind of separation of the estimation and control is valid in this case. As control algorithm we use the method presented in Section 3. As the weak component concentrations are now used the estimated values instead of direct measurement results, which was the case in Section 3. The value of (C, - C,) is now determined in a different way as was discussed in Section 3. Let us now suppose that we use the model of (4) and that the concentrations of the weak species in the last tank (of the three tanks in series model) are Co,p = (Cal, Ccu2, C,, , Cp,). Let us further suppose that the changes of these concentrations are slow compared with the reactor dynamics so that they can also be used as estimates for the inputs of the first two tanks. Now the concentration of the strong species in the input (C,, - C,,) can be computed from (4) and (5), when the input pH is measured. This concentration is then summed to the actuating control concentration C,,Q,/(Q, f Qp) and the sum is fed to the dynamics of the model. Thus we obtain the strong component output (CA - Cn). It has been supposed above the control flow to be so small, that it won’t change the troughput flow considerably. The value of pH, both at the output and the input point, depends very nonlinearly on the concentrations of the species, according to the electroneutrality con-
model
dition (4) and the definition of the pH-variable. However, the values of C+ as well as C+ - K,/C+ in (4) are uniquely determined by the pH-values. By inspection the structure of (4) it is convenient to use the quantity C’ - Kw/C+ in the construction of the final output variable because (4) is linear with respect to the other concentration variables. Naturally, the pHmeters produce measuring errors, so that the calculation of C+ or Ci - Kw/C’ is not accurate. In this paper we suppose that the errors can be interpreted as an additive term u in the quantity C+ - K,/C+. Further, we suppose that there exists an additive 4dimensional noise vector w in the weak concentration components of the vectoral concentration input CM, cfl2k) of the last stirred tank (v) = (&lk, CorZk, cp,k, for modelling the disturbances in the inner dynamics of the process. By these suppositions the structure of the output section of the process can be depicted by the Fig. 8 below. The part enveloped by broken lines in the Fig. 8 represents the memoryless output mapping of the system with the primary output $I, the value given by the pH-meter. The additional output y is n_ow determined from $I by the relation y = (lO-PH - K, IOp”H,(molesllitre) thus giving a noisy estimate for the quantity C+ - K,/C+. The final output 7 is calculated as the quantity y - (CA - C,). As mentioned before the strong component output (C, - C,) is determined by the system dynamics. By electro-neutrality condition (4) the additive error u in 7 is then determined by the relation: 7 = HCap + v, where HC,
=
K%+GI C+2 + K&C+ + KiIK$
-K$ v (CA-CE)k -
-2
v
Fig. 8. The output system of the process, where cab,, Cc&k, +lk,
cflk,
and CC@= (c,,,
c,,
$1,
‘+.$.
= (Calk,
(8)
t K;,K,C+
•tK;,K;,Ci2
K;,K;,P2Cp, K;
t K;,K,C+
+ K&K;,C+2
(9) f&C%32
- Kt2 •tK,
The error u can mathematically be interpreted as a noise although it appears as an output of the system by Fig. 8 (the discrepancy is caused by the ‘backward transformation from j5H to v). The quantity HC, in (8) and (9) is linear in the concentration variables to
.
R Jutila et al. / Physico-chemical
Flow
model
105
and
Fig. 9. The simulation scheme of combined estimation and control.
be estimated,
but the row vector H depends
on the
concentration C+ . An estimate for C+ can be calculated from the measured p”Hvalue. Now the quantity 7 in (9) is thought to be the final output quantity, which will be used for estimation. When the dynamic equations for the concentrations of weak components Co,p are combined with the output (9), a linear differential system is obtained, although the coefficient in (8) depends on the output itself. With a suitably defined quadratic cost function, a proper filtering problem is obtained. It can be shown theoretically, [ 121, that the solution of this problem yields an algorithm, which has the same structure as that of the linear Kalman-fdter. The location of process noise only in input of the last tank reduces the dimension of the Riccati equation considerably. The addition of an exponential decaying weighting term into the cost function makes the algorithm more adaptive. Thus the weak components C$ have only a first order dynamics driven by noise w and only the strong component (C, - Cn) has the full dynamics of the tank series (see Fig. 2). The combined control and estimation simulation scheme is shown in Fig. 9. It can be noted that we can not calculate the exact value j7 but only an estimate of it 7 because we must use an estimate eA- en instead of CA - CB . This originates in the fact that in
practice we must use a noisy (measured input pHvalue in the calculation of the quantity C, - CB.
6. Conclusions A Kalman-filter approach to estimate hypothetical concentrations of a physico-chemical model in connection with a continuous pH control algorithm is presented. This linear control algorithm causes the controlled pH to converge to setpoint value, if the system is tuned to be stable. The method can be extended to cases where the control reagent is weak and manyvalued. On the other hand processes with dissolving kinetics (e.g. the case of solid reagent) cannot accurately be included in this problem formulation. If the throughput flow is time variant, the problem can be maintained time invariant by choosing the sampling time inverse proportional to the flow. The dimension of the model is restricted by the computing capabilities. In addition, the use of high dimensional models may cause uncertainties and singularities in the estimation of parameters and states, as well as in the tuning of algorithm. The presented algorithm can also be combined with any kind of feedforward! control algorithm. One possibility is to use the estimated weak concentrations in the algebraic calculation
106
P. Jutila et al. / Physico-chemical
of feedforward control and combine it with the output of the feedback linear controller. The research of the method is going on in the Control Engineering Laboratory. Parallels will be made to other advanced adaptive feedback algorithms reported in literature (Shinskey [8] ; Boberg, [9] ; Bucholt and Kiimmel [lo]; Gustavsson [ll]).
Appendix ion concentration ion concentration concentration of strong acid(s) concentration of strong base(s) concentrations of weak acids concentration of weak bases concentration vector (C,, , Cor2, Cpl , $2)
f, k
concentrations of undissociated acids concentrations of undissociated bases concentrations of weak anions concentrations of weak cations dissociation constants of weak acids dissociation constants of weak bases ion product of water volumetric flow input noise measurement error signals measured pH values output noise output matrix output variables as subindices of variables refer to process feed and control, respectively as subindices refer to first and second tank of the dynamic model, respectively - the third tank is referred
model
with symbols without subindex as superindex refers to estimated variables
References [ 1] V.L. Trevathan, Advanced Control of pH. ISA preprint, 69-82. [2] T. MacAvoy, Time optimal and Ziegler-Nichols control, Process Des. Develop. 11 (1972) 71-78. [3] I.D. Richter, C.D. Fournier, R.H. Ash and S. Marcikic Waste neutralization control, Instrum. Technol. 21 (1974) 35-40. [4] A. Niemi and P. Jutila, Algorithms for pH control, 2nd Polish-English Seminar on Real-Time Process Control 22-23 May, (1978), Rogow, Poland; also Rev. Tee. Ing., University of Zulia, Maracaibo, Venezuela (1979), Vol. 2, No. l-2. P. Jutila and B. Salmelin, A control method for neutralization of weak acid in stirred tank reactor, poster, (1979) 27th IUPAC Congress. [6 B. Salmelin, Testing of a pH control method and apparatus (in Finnish), M.Sc. Thesis, Helsinki University of Technology. and testing of some spe]7 P.K. Jutila, The implementation cial algorithms of pH control in a computer controlled stirred tank reactor (in Finnish), Lit. Tech. Thesis, Helsinki University of Technology. I F.G. Shinskey, Adaptive pH controller monitors nonlinear process, Control Engrg. 2 (1974) 57-59. [9] L. Boberg, An attempt with adaptive pH-control (in Swedish) Report E 77-13, Chalmers University of Technology. [lo] F. Bucholt and M. Kiimmel, Selftuning control of pHneutralization process, Automatica 15 (1979) 665-671.
is
[ 111 T.K. Gustavsson,
Adaptive pHcontro1 with least squares parameter estimation, Report 79-7, Abe Akademi, Department of Chemical Engineering (1979). [ 121 M.T. Nihtila, Optimal state-linear filtering trough implicit output equation, accepted for publication in preprints of IFA/1981-World Congress, Kyoto, Japan (1981).