- Email: [email protected]
Simultaneous estimation of kinetic and heat transfer parameters of a wall-cooled fixed-bed reactor
Pergamon PII:
Chemical Enoineerin O Science, Vol. 51, No. 21, pp. 4791 4800, 1996 Copyright © 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved S0009-2509(96)00330-2 0009-2509/96 $15.00 + 0.00
S I M U L T A N E O U S ESTIMATION OF KINETIC AND HEAT TRANSFER PARAMETERS OF A WALL-COOLED FIXED-BED REACTOR Z H E N - M I N CHENG and WEI-KANG YUAN* UNILAB Research Center of Chemical Reaction Engineering, East China University of Science and Technology, Shanghai 200237, Peoples' Republic of China
(Received 10 July 1995) Abstract--Heat transfer and kinetic parameters of wall-cooledfixed-bedreactors are traditionally obtained through separate heat transfer and kinetic experiments. In spite of its simplicity,the uncertainties resulting from the low parametric sensitivity and high correlationship among the parameters inevitably make the estimation unsuitable for safe extension to actual reacting conditions. In this paper, an effort was made to estimate simultaneously both heat transfer and kinetic parameters under reacting conditions in a single tube wall-cooled fixed-bed reactor, and a two-stage parameter estimation procedure was developed. The correlation and confidence region analysis performed in this paper verifies both theoretical and experimental feasibility of simultaneous estimation of heat transfer and kinetic parameters under real reacting conditions. Copyright © 1996 Elsevier Science Ltd
Keywords: Fixed-bed reactor, parameter estimation, correlation analysis, confidence region. INTRODUCTION The difficulties encountered in fixed-bed reactor design can be ascribed to two reasons. Firstly, models describing reactor performances are usually oversimplified, particularly when the experimental ratio of dr~dr is relatively low (less than 10), leading to enormous influence of particle size on gas velocity distribution, local packing density of solids and gas-solid external transport coefficients. Besides, the randomness of catalyst packing, which causes heat transfer coefficients to differ by as many as three times under the same operating conditions for different experimental runs (Wijngaarden and Westertcrp, 1992), makes decisive modeling impractical. Secondly, the fixed-bed reactor is a very sensitive system due to the exponential dependence of reaction rate on temperature and the low thermal capacity of the system, and hencc any slight change in operating variables will result in large variations in temperature and concentration, even causing reactor runaway. Different opinions have been expressed in designing a fixed-bed reactor. Pcterson and Carberry (1983) urged that as many parameters as possible be independently measured since simultaneous estimation of too many parameters would make thc results meaningless because of their strong correlation. According to this point of view, the heat loss parameters and the reaction rate parameters should be obtained separately. This approach has the advantage of being relatively simple, since only a small number of parameters is estimated at a time. However, the reactor model must be very accurate. Clement and Jorgensen * Corresponding author.
(1983), on the other hand, insisted on the simultaneous estimation of all the parameters as a reliable attempt especially when a simplified model was used, for example, the plug-flow model. It should also be pointed out that good results cannot be achieved from separate estimations as expected. It is well known that the correlation between radial thermal conductivity and wall heat transfer coefficient is very high, as can be observed from the well-known relationship 1
1
dt
U - hw + 8Ker Strong correlation is also frequently encountered between the pre-exponential factor and the activation energy in the Arrhenius expression, e.g. as high as 0.9999, as was noted by Pritchard and Bacon (1978). Further investigation has also shown that heat transfer parameters under reactions could differ considerably from those under reaction-free conditions (Hofmann et al., 1984). It should also be noted that the low parametric sensitivity of heat transfer parameters under non-reacting conditions cannot meet the requirement in the accuracy for those under reacting conditions. In consideration of the overall inherent problems existing in separate parameter estimations, simultaneous estimation of heat transfer and kinetic parameters may be a better recourse for reactor design. EXPERIMENTAL
Determination of heat transfer parameters of a fixed-bed reactor under non-reacting conditions has been conductcd by measuring radial temperaturc
4791
4792
Z H E N - M I N CHENG a n d W E I - K A N G YUAN
profiles at different bed depths by progressively adding catalyst pellets up to different levels, as proposed by Coberly and Marshall (1951). Under reacting conditions, such a method was also used by Valstar et al. (1975). It can be imagined, however, that the method of increasing reactor length by adding catalyst would correspond to a reactor with varying length. According to Danckwerts'-type boundary conditions, the combination of all the exit concentrations and temperatures of a length-varied bed cannot represent the real situation of the whole reactor. To overcome such a drawback, a method was proposed, which is an improvement of Hofmann's (1979) contribution. The reactor, made of 4 cm o.d. stainless-steel tube, consists of six separate sections, is shown in Fig. 1. The lengths of individual sections, counted from the bottom to the top, are 10, 15, 15, 20, 20 and 30 cm, totalling 110 cm. Compared to the catalyst diameter of 3-5 mm, the reactor is long enough to make end effect negligible on axial dispersion. The reactor tube is provided with a jacket (details shown in Fig. 3) for coolant circulation in a 4 mm channel• Conduct oil (SD-320), which can resist a temperature as high as 320°C, is recycled by a high-temperature centrifugal oil pump at a flow rate of 1.7 m3/h between the jackets of the six reactor sections (in parallel, shown in Fig. 1) and an oil tank equipped with a controlled 3 kW electric heating coil. Catalyst pellets are packed only in the four middle sections (5, 15, 20 and 20 cm, totalling 60 cm) while the top and bottom sections are packed with inert ceramic packings. Between every two neighboring sections is sandwiched a specially designed composite probe (Fig. 2), each equipped with 5 gas sampling tubes and 9 thermocouples, located at different radial positions, for radial concentration and temperature measurement. The 9 thermocouples are held in place by a ceramic monolith, which also functions to suppress gas mixing. Temperatures are acquired through a computer which is also used for controlling the operating conditions. A parallel-sequential exothermic reaction, i.e. catalytic oxidation of ethanol by air over CuO catalyst supported on 7-A1203 pellets (3-5 mm in diameter), is used as the working system: (A)
P = E(AkAk T) F
~k2,
cov(ka, k2)
•.-
"-I
cov(kl, kp)[ /
= I COY(k2, kl)
322
! [_cov (kp, k0
cov(kp, k2)
•"
__I
(2) The parameter correlation matrix C, whose element C~j is defined as coy (kl, k j)
Co
6k,6~
can be easily obtained from matrix P. This knowledge, concerned only with parameter variance and parameter correlation, is, however, not enough to describe the actual parameter confidence region. It is known that for a p-dimensional parameter vector, the confidence region is a hyperellipsoid centered at the estimated value k in the p-dimensional parameter space (Hosten, 1974), and the size of the confidence region is equal to the square root of the determinant of matrix P. Complete information concerning the parameter estimation should include the following aspects: (1) (2) (3) (4)
k, the parameter estimation; 6kl, the standard deviation of ki; CO, the correlation between k~ and ks; x//~l, v~22 . . . . . x/~p, the axial lengths of the pdimensional hyperellipsoid; (5) IPI l/2, the size of the confidence region which is equal to 2x/~x2~ ..- 2p. Given the above-mentioned values, the linearized individual confidence interval for a parameter under
2
, C2H40 + (-AH1) +02
The task for parameter estimation is to find the optimal value of a parameter and its confidence region. Assume the estimated error of a parameter vector k is Ak, and its error distribution is represented by a covariance matrix P. Then the relationship between P and Ak can be established (Seinfeld and Lapidus, 1974):
(c)
(B) 1
C2HsOH
PARAMETER ESTIMATION
3
+ Oz
) CO
2 +
H20 + (-AHa).
(1)
T
-~-O2 The mechanism of this reaction is the redox type. According to the Langmuir-Hinshelwood mechanism (Smith, 1981), the rate expression under high oxygen concentration (21%) and low reactant concentration (1%) can be reduced to a simple Arrhenius first-order form. Experiments were carried out under the operating conditions shown in Table 1.
95% level according to the normal distribution is readily found to be: kl ± 1.966ki. The heat and mass balances for the two-dimensional plug-flow model are described in dimensionless forms as follows:
o¢ =~'\-Y7
p ep/+'~'~"
(3)
4793
Simultaneous estimation of kinetic and heat transfer parameters to M.S.
oil out
sampling tubes
thermocouples ,°
. . . . . . . . . . . . . ,
° . . . . . . . . . .
oil pipe
\
o,j
16 way ~ valve
converter
recorded 1 and displayed by computer
adiabatic section
recorded by mass spectrometer
inert packing I
premixed gas
Fig. 1. Overall schematic structure of the reactor and the analysis arrangements.
i I I I
thermocouple I
I
I
I
cerami, monoli
III
sampling tube
:111II I
_1
t
catalyst pellets,I reacting gas I
[1II
reacting~ ~e probe Fig. 3. Interior detailed structure of the reactor.
Fig. 2. The composite temperature and concentration measurement station.
The reaction rates in eqs (3)-(5) are expressed as
~X C
:~2X C
~0
(~2o
l ~Xc~
(4)
~ n = y,4oYo~[kl(l - Xn - X c ) - k z X n ] ~ c = YaoYo2 [(k2Xn + k3(1 - X s - X c ) ]
1 ~0\ (5)
(9) (10)
and the reaction rate constants k~, k2 and k3 are assumed to be of the Arrhenius type:
with the boundary conditions
~=0,
xB=Xc=O, eXR
p=0,
p=
~p
~Xc
(6)
o=1 ~0
- . . . . ~p ~p
(7)
0
-~-~
,
i = 1,2,3.
(8)
8p
(I1)
The dimensionless parameters involved in eqs (3)-(8) are defined as 4Ldp
~0 - - = -- B i ( O - - 1 ) .
1,
ki=kioexp
~m
="
2
d, Pe~,
4Ldp '
~h
~-
d~Peh
Table 1. Operating conditions and the maximum observed temperature rises No.
YAo(%)
Yao (%)
Yco (%)
T~, (K)
To (K)
G (kg/m 2 s)
AT=,, (K)
1 2 3
1.044 0.94 0.86
0.0 0.0 0.0
0.0 0.0 0.0
498 503 507
498 503 507
0.295 0.369 0.369
25 33 57
(12)
4794 /~1 =
ZHEN-MIN CHENG
L UCAo '
~
=
L(-AHt) GCp
-
,
/h
=
hwdt Bi = - 2Ker
L(-AH3) GCp (13) (14)
Since the pre-exponential factor ko and the activation energy E appear as a product in the Arrhenius kinetics, this will result in a correlation of the two parameters higher than 99.9% and thus leading to a serious illness in the objective function. According to the theory of regression analysis, the reaction rate constant k should be reparameterized as k = ~exp[~(1 - 1/0)]
(15)
where ~t = koexp(-E/RTr=f)
and
y = E/RTror. (16)
The reference temperature T, ee in eq. (16) is set equal to the inlet temperature To. Generation of initial values of parameters The initial value of a parameter can initially be generated either by the random search method or by an empirical relation. However, because of the high non-linearity of the process, poor initial value may lead to failure in optimization. Understandably, a good initial guess is preferably obtained from experimental data. As can be found from eqs (3)-(5), it is the reaction rates ~n and Ytc which make the three partial differential equations and thereby all the model parameters coupled together through their functional relations with 0, Xn and Xc. Therefore, ~B and 9tc have to be known a priori so as to make all the three equations independent. It seems that this requirement can only be fulfilled by employing the differential method (Hosten, 1979; Dimitrov and Kamenski, 1991; Kamenski and Dimitrov, 1993). It is observed here that if the axial and radial derivatives are substituted by their numerical values, the values of ~B and ~ c can be obtained directly from eqs (3) and (4), and the three equations involving the two heat transfer parameters and six kinetic parameters can thus be decoupled. An efficient initial parameter value generating procedure is outlined according to Cheng and Yuan (to appear):
(1) Interpolating the values of 0, XB and Xc at dimensionless radial positions p = 0.577 and 1 through profile fitting. (2) Evaluating the radial derivatives by orthogonal collocation with one collocation point. (3) Using the Pad6 rational function R(() =
ao + aa~ + a2~2 + a3~ 3 1 + bl( + b2( 2 + b 3 ( 3
to approximate the measurement of 0, Xn and Xc at p = 0.577 in the axial direction from ( = 0 to ( = 1. The axial derivative is then calculated through analytical differentiation of the rational function.
and WEI-KANGYUAN (4) Evaluating the reaction rates ~n and ~ c from the two mass balance equations (3) and (4), and then substituting them into the heat balance equation (5) to evaluate the heat transfer parameter cth. (5) Obtaining Bi from the boundary condition (8). (6) Kinetic parameters of the three reactions can finally be obtained from reaction rates ~ e and ~Rc at various values of 0, Xn and Xc through non-linear least-squares regression. Precise determination of the model parameters The Gauss-Newton-Levenberg-Marquardt method (Seinfeld and Lapidus, 1974) is employed to produce the final precise estimation. Orthogonal collocation with four interior points is used to convert the three partial differential equations into a set of ordinary differential equations. The parameter vector to be estimated is composed of two heat transfer and six kinetic parameters, and in dimensionless form, is written as
k
= (¢th, B i , ~ t , ~ l , ~ 2 , ~ 2 , ~ 3 , ' ~ 3 ) T .
(17)
There are therefore 12 state equations and 96 parametric sensitivity equations involved, and these 108 differential equations have to be integrated from ~ = 0 to ( = 1 simultaneously, for which it is found that the Runge-Kutta-Merson algorithm is suitable. RESULTSAND DISCUSSION Initial estimates With the experimental temperature and concentration data listed in Table 2, the differential method is first used to produce initial estimates of the eight parameters, with the results shown in Table 3. It is observed from Table 3 that the preliminary estimates for different conditions are in good agreement. However, the objective functions SSRs differ by orders of magnitude: under operating condition no. 3 the SSR is about 1000 times of that under condition no. 1. However, the small SSR value for condition no. 1 does not definitely mean more accurate parameter estimation, since it may probably be caused by low parametric sensitivity. Precise estimates Simultaneous parameter estimation of heat transfer and kinetic parameters has been touched in the literature (Beck and Hassel, 1971; Peuider et al., 1971; Emig and Hosten, 1974; Clement and Jorgensen, 1983; Hofmann et al., 1984). In Clement and Jorgensen's work, the set of partial differential equations of the system were transformed into algebraic equations by double collocation, resulting in 2 and 15 collocation points in the radial and axial directions, respectively. In order to improve convergence, a parameter transformation method suggested by Michelsen (1972) was employed. It is obvious that double collocation has serious limitations in its application, because for example, if 4 radial points and 15 axial points are adopted, 180 algebraic equations will result, along with 180 x 180
Simultaneous estimation of kinetic and heat transfer parameters
4795
Table 2. Experimental data of temperature and concentration at different axial and radial positions (1) Operation no. 1 Temperature (K)
Concentration
zi\ri
rl
r2
r3
r4
r5
r6
r7
r8
r9
zl z2 z3 z,,
506 516 522 510
506 516 522 510
506 514 520 509
505 513 520 509
505 512 519 507
504 511 516 506
504 5ll 515 505
503 509 512 504
503 506 509 504
Zi\ri
YA
Zl
YB YC
z2
YA YB Yc YA
Z3
YB YC
Z4
Ya yn YC
(mol%)
rt
r2
F3
r4
F5
0.9672 0.0616 0.0197
0.9608 0.0623 0.0199
0.9644 0.0618 0.0198
0.9677 0.9693 0.0594 0.0574 0.0189 0.0182
0.6434 0.2736 0.1261
0.6463 0.2744 0.1260
0.6525 0.2678 0.1208
0.6818 0.7008 0.2513 0.2420 0.1088 0.1024
0.3326 0.4031 0.3047
0.3377 0.4039 0.3034
0.3425 0.4019 0.2933
0.3795 0.3924 0.3960 0.3925 0.2722 0.2615
0.2132 0.4159 0.4160
0.2101 0.4161 0.4156
0.2146 0.4169 0.4108
0.2218 0.2308 0.4182 0.4188 0.4005 0.3952
(2) Operation no. 2 Temperature (K)
Concentration (mol%)
z~kr~
rl
r2
r3
r4
rs
r6
r7
rs
r9
zl z2 z3 z,
512 531 537 524
512 528 537 523
512 527 534 521
511 525 532 520
511 523 530 519
510 521 528 516
510 519 523 513
508 517 520 512
507 514 516 510
zi\ri
Z1
Ya YB YC
Z2
YA yn
Z3
YA Ya
YC
YC
YA Z4
YS YC
rl
r2
r3
r4
r5
0.8648 0.0541 0.0175
0.8680 0.0547 0.0177
0.8700 0.0544 0.0176
0.8706 0.8748 0.0525 0.0509 0.0169 0.0163
0.5816 0.2410 0.1128
0.5865 0.2417 0.1127
0.6002 0.2359 0.1079
0.6290 0.6375 0.2214 0.2131 0.0971 0.0912
0.2819 0.3582 0.2927
0.2900 0.3594 0.2907
0.3013 0.3579 0.2787
0.3382 0.3418 0.3528 0.3497 0.2544 0.2422
0.1753 0.3622 0.4043
0.1783 0.3626 0.4036
0.1766 0.3640 0.3976
0.1848 0.1963 0.3665 0.3677 0.3850 0.3787
(3) Operation no. 3 Temperature (K)
Concentration (mol%)
zi'kri
rl
r2
r3
ra
rs
r6
rv
rs
r9
zl z2 z3 z,
517 540 565 526
517 538 565 524
517 535 562 524
516 532 555 523
515 530 553 520
513 527 550 518
513 524 543 515
511 521 534 514
509 516 523 511
zi\ri
Z,
YA YB YC
z2
YA Ys Yc
z3
Ya ya Yc
z4
ya YB Yc
rl
r2
r3
r4
rs
0.7745 0.0548 0.0175
0,7699 0,0525 0,0178
0.7803 0.0491 0.0166
0.7824 0.7800 0.0469 0.0450 0.0158 0.0151
0.5057 0.5104 0.2276 0.2236 0.1301 0.1293
0.5381 0.2021 0.1117
0.5679 0.5811 0.1956 0.1862 0.0995 0.0969
0.1311 0.2960 0.4324
0.1368 0.2972 0.4309
0.1910 0.3028 0.3632
0.2444 0.2613 0.3019 0.3085 0.3175 0.2971
0.0690 0.2407 0.5513
0.0694 0.2527 0.5499
0.0838 0.2534 0.5221
0.0991 0.0104 0.2626 0.2669 0.4963 0.4907
(4) Positions of measurement probes Axial positions of sampling tubes and thermocouples: zj = 5, 20, 40, 60 cm from the inlet. Radial positions of the 5 sampling tubes: r~ = 0, 3, 7, 13 and 18 mm, where 0 mm is the center of the bed. Radial positions of the 9 thermocouples are listed in the underlying table in "ram", where 0 mm is the centre of the bed.
ZHEN-MIN CHENG and WEI-KANGYUAN
4796
Table 2. (contd.)
zi\r~
rl
r2
r3
r,
rs
r6
r7
rs
r9
z~ Z2 Z3 Z4
2.5 3.1 2.4 3.3
5.1 6.0 3.8 4.9
7.2 7.4 7.2 8.2
9.1 10.6 8.8 9.5
11.0 11.7 10.2 11.3
12.9 13.4 12.0 13.5
14.0 14.2 14.2 14.8
15.1 15.0 15.4 15.9
16.1 16.7 16.9 16.5
Table 3. Preliminary results of parameter estimation Operating condition
cth
Bi
0tI
Yl
~2
~)2
0t3
Y3
SSR
l 2 3
3.96 3.45 3.20
3.53 3.19 3.88
- 1.82 -1.54 - 1.26
27.47 26.52 21.88
- 2.27 -1.94 - 1.71
23.28 21.40 20.37
- 2.97 -2.82 -2.49
27.20 31.73 27.60
0.00036 0.0058 0.3588
Table 4. Precise estimation according to Gauss-Newton-Levenberg-Marquardt Operating condition
~h
Bi
~q
71
~2
Yz
a3
73
SSR
1 2 3
3.24 3.14 2.67
5.18 3.99 4.99
- 1.81 - 1.43 -1.40
26.56 24.08 22.80
-2.51 - 2.09 -1.63
30.04 25.05 20.78
-2.91 - 2.60 -2.61
23.84 26.13 28.15
0.00004 0.00003 0.0005
working matrices. In comparison, the largest matrix involved in this work is only 12 x 8. Actual iteration shows that SSR can be reduced to a relatively small value after only about 4 iterations, and the final result can be achieved within 7 to 9 iterations (see Table 4).
Statistical properties of estimated parameters Individual confidence limits. The deviation of the estimated values from the optimal estimates is the most desired information concerning the accuracy of parameter estimation. Standard deviation of parameter ki, defined as 6ki, can be readily obtained from the square root of the ith diagonal element of P. The individual confidence region at the selected probability level of ( 1 - ~) is given by the well-known relation (Seinfeld and Lapidus, 1974) ;,, - l~rk, <. k* <~ i, + [~6~,
(18)
where/~ depends on the confidence level. At the 95% level/~ = 1.96. The individual confidence region for each dimensionless and physical parameter is shown in Table 5.
Joint confidence region. By assuming that the model can be adequately linearized in the vicinity of the parameter estimate k, and if 62 is not known, the joint confidence region can be defined by F-distribution: (~ _ k , ) T p - 1 ( ~ _ k*) = 62pF~ _,,(p, m - - p ) (19) where ~ 2 = S S R m i , / ( m - p). One notices that there will be p variables on the left-hand side subject to the
same eq. (19). In this work, the hyperellipsoid confidence region cannot be represented by the solution of eq. (19) since the parameter number is eight. Instead, the confidence region can be approximately described by the length of each axis as well as the size of the confidence region (see Table 6). The relative ratio of the longest to the shortest axis ranges from 1118 to 14,425, indicating an extremely elongated hyperellipsoid. It shows further that the shape is more elongated under operating condition no. 1 than under conditions nos 2 and 3. Besides, the size of the confidence region under the no. 1 condition is larger by 10, 000 times than that under no. 3 condition. These findings are promising, since a precise estimate can thereby be obtained under sensitive conditions like condition no. 3. Further description of parameter estimation should be referred to the correlation matrix, as can be found in Table 7. The correlation matrices under the first two conditions are similar. The correlation between the two heat transfer parameters ~h and Bi is about 0.8, a substantial reduction compared with 0.95 when no reaction occurs. The correlation between the pre-exponential factor ~i and activation energy ?i of the same reaction is about 0.95, larger than that between different reactions. The most exciting result is observed from the correlation coefficient between one of the six kinetic parameters, ~1, 71, (z2, ?2, 0c3 and Y3, and one of the two heat transfer parameters, ~h and Bi, with an average value of about 0.3, the largest values being 0.488 for condition no. 1 and 0.517 for condition no. 2. F o r condition no. 3, the correlation between the two heat transfer parameters, cth and Bi, has decreased
Simultaneous estimation of kinetic and heat transfer parameters
O +1 -H +1 ",~" eq '~1"
-H -H +1 o~tt~
eq eq ¢-q
¢.q t-,q ¢-.1
X
X
X
-H -H -H O
¢-q ¢-q ¢-q
I
+1 +1 +1
4797
to 0.2, but positively correlated, and additionally the correlation between the kinetic parameters has also decreased considerably. Another important finding is that the average value for the correlation between one of the six kinetic parameters ~a, ~1,0~2,~2,5 3 and ~'3, and one of the two heat transfer parameters, ~h and Bi, still remains at 0.3. Although the reactor behaviors under the three conditions are different in many respects, one thing however is similar: the correlation between a kinetic parameter and a heat transfer parameter remains almost unchanged at an average value of 0.3. This value implies that the correlation between the kinetic and the heat transfer parameters is weak, and therefore simultaneous estimation of transport and kinetic parameters is theoretically feasible.
O~ tt% 1"--
+1 +1 -H
-I-I +1 -H ee~
~,'5 ¢N exl
o~,o
.~- t-q ¢'q
X
X
X
+1 +1 +1 ¢'q ¢'-,i ~
+1 +1 +1
I I I
e,q
e.
e-q
+l +l +l
+1 + l + t
8
i'-I o ¢'q ¢-q eq
I
w; o
~
¢.i i~.1 ¢-i
t~ eO
X
O
+1 +1 +1
X
X
+l +l +1 I
[--
7771
• r..: q-~ t"q ~ O +1 +1 +1
e~ +1 + l +1
~'~1 .-.-~ tt~ ¢'.q + i +1 +1
+1 +1 + l p... ~ tt% "2. "n: 'Q I
e6 ,A e.i e-,
o
..
Comparison with reactor performance Gas velocity and bed voidage fluctuate across the tube diameter; however, the measured radial concentration and temperature profiles are smooth, as can be found in Valstar et al. (1975), Hofmann (1979) and Daszkowski and Eigenberger (1992), as well as in our experiment• Perhaps local fluctuations in temperature and concentration have been damped by mass and heat diffusions• Therefore, when designing wallcooled fixed-bed reactors, the two-dimensional plugflow model permits a reasonable reactor performance description• A comparison is made in Figs 4 and 5 between plug-flow model predictions and experimental findings. It is found that fitting of concentration is very good (Fig. 5), but for temperature (Fig. 4), a substantial disagreement is observed within a distance of about one particle diameter to the wall. It is believed that such a result is mainly because of the different influence of the wall on the heat and mass transport process• The existence of the wall is advantageous in damping concentration differences since it confines the mass transport within the tube. Therefore, in most cases, the concentration profile is usually relatively flat or S-shaped, as shown in Fig. 5. On the contrary, the influence of wall on heat transport is different, causing a decreased thermal conductivity, a higher bed voidage and a shorter gas-solid contacting time in the 'wall region'• All these effects combine to make the radial temperature profile near the wall steeper than the plug-flow prediction, probably because the plug-flow simplification assumes a constant radial effective thermal conductivity and a uniform voidage even up to the wall (Cheng and Yuan, 1993).
o CONCLUSIONS
O
(1) A two-step parameter estimation method is presented in this work. In the first step, an initial value is generated from experimental data without making any assumptions, and in the second step, the Gauss-Newton procedure is employed for further precise estimation. Although 108
ZHEN-MIN CHENG and WEI-KANG YUAN
4798
Table 6. The confidence region under different operating conditions Operating condition 1 2 3
~
~
~33
~
x/~5
x~6
~
10.0 5.87 1.54
3.50 1.97 1.03
2.12 1.36 0.84
1.14 0.52 0.13
0.085 0.058 0.038
0.031 0.019 0.011
0.0041 0.0076 0.0034
x~s
IPI
0.0023 0.0053 0.00058
2.06x 10 -6 3.53 x 10-7 1.43 x 10- to
Table 7. Correlation matrix for all the eight parameters under different operating conditions (1) Operating condition no. 1 a~, 1.00000
ah Bi
--0.81915 -0.40933 0.48811 -0.31009 0.36463 0.05445 -0.08988
at 71 a2 72 a3 ?3
Bi
at
Yl
a2
?2
a3
Y3
1.00000 0.35096 -0.21491 0.31197 -0.24795 -0.14301 0.23354
1.00000 -0.91410 0.60761 -0.65237 -0.70207 0.68397
1.00000 -0.54895 0.68487 0.61478 -0.58653
1.00000 -0.96408 -0.71644 0.70475
1.00000 0.73865 -0.72951
1.00000 -0.97794
1.000~
(2) Operating condition no. 2 ah Bi
al 71 a2 ?2 a3 73
ah
Bi
at
?t
a2
?2
a3
1.00000 --0.75367 -0.40692 0.51703 -0.36376 0.44763 0.13693 --0.16799
?3
1.00000 0.32876 -0.18089 0.42442 -0.33805 -0.25147 0.37214
1.00000 -0.90719 0.64322 -0.70664 -0.76505 0.72490
1.00000 -0.55108 0.72031 0.64709 -0.58556
1.00000 -0.95468 -0.76995 0.76226
1.00000 0.78753 -0.76587
1.00000 -0.97317
Bi
at
Yz
a2
?2
~+
1.00000 0.22007 0.20241 0.60427 -0.47872 -0.28875 0.35252
1.00000 -0.69684 0.47067 --0.59403 -0.75131 0.52766
1.00000 0.05695 0.32318 0.36685 -0.20584
1.00000 -0.89146 -0.67678 0.34740
1.00000 0.73943 -0.37199
1.00000 -0+83327
1.00000
(3) Operating condition no. 3 ah 1.00000
ah Bi
0.20189 -0.11682 0.64973 0.14422 0.23299 -0.19265 0.45371
al ?t a2 Y2 a3 ?3
?3
1.00000
,540 "Z=SCm ÷ Z = 2 0 ~ 530
""z = 5 m -='z ÷Z =40~ +Z== 20cm WelII
-*z = 4Oan + z = 60cm
i
; 5 3 0 ............
.... 51 (~ ......~
"
". .
I
I
3
6
"*Z=SCO1 ÷z=2Ocnl
5 7 0 f. ........
* z = 40cm + z = 60crn
s6o
...
....
.
"'..
{
550
s,o
iii'. l
3o..
[
i\" .....
....
, ....
, 510]
0
)
9
12
15
18
r, m m
A: operating condition No. 1
5000
. . . .
~
~
;
12
15
18
r, mm
B: operating condition No.2
0
3
6
9
12
15
r, mm
C: operating condition No.3
Fig. 4. Radial temperature distributions: difference between measurement and plug-flow prediction.
18
Simultaneous estimation of kinetic and heat transfer parameters • 1=5¢m
4799
nzm2Ocnl
,z=Scm
• zm4OCnl • z-8Ocm O.B
0.41-
" Z=SCm
0.1~
~ Z = 2(] ¢1~
0.8
o at
* Z = 40cra " Z = 60Cm
xz=2Oom
o z : 4 O c r l t v z:l~O(:Im
_:. ___---~
06 ~ - ' ~ - ' - - ~ "
--..____~_~
o.21
o.4
~
0.1
0.2
.
•
0.4:
t
0.2 --
~t I
% a
,t
.t I
6
I
9
,¢--
I
I
12
15
11m
O.
0
.
.
3
.
6
r, m m
.
.
9
~
12
15
1
0
~t
3
r, m m
[ .
[
6
9
2&
1
5
1
r, m m
Fig. 5. Radial concentration distributions (operating condition no. 3).
differential equations have to be integrated simultaneously when four collocation points are used, it is still an easy task using a personal computer, since CPU storage and processing time can be greatly saved as compared to other methods. (2) Since the parameter space is highly elongated, a variety of gradientless methods, such as Simplex or Powell, are less effective in minimizing the objective function. (3) A positive answer can be given to the availability of simultaneous estimation of heat transfer and kinetic parameters. Mathematically, the average correlation between heat transfer and kinetics parameters is only about 0.3. On the other hand, they are physically different: the heat transfer parameters are related to physical properties of the fixed bed as well as the catalyst and the fluid, whereas the kinetic parameters are only related to the chemical reaction, (4) Parameter estimation performed under reacting conditions shows that the size of confidence region under conditions of different parameter sensitivity can differ up to 10,000 times, which implies that it is necessary to do parameter estimation under relatively sensitive conditions.
1 ct point of the F-distribution mass velocity, kg/m2s wall heat transfer coefficient, J/m 2 s K pre-exponential factor parameter vector axial dispersion for heat effective thermal conductivity, J / m s K catalyst bed height, m number of experimental data points numbers of parameters covariance matrix of parameter Peclet number for radial heat transfer Peh Peclet number for radial mass transfer Peru radial position, m or mm r gas constant, cal/mol K R reaction rate, kmol/kgcat hr .:~ sum of squares of residues SSR t l -(~/2) t value at 1 - ~ confidence level temperature, K T linear velocity of reactant, m/s u heat transfer coefficient, J/m 2 s K U measurement covariance matrix V yield of a product X mole fraction y axial position, m z
F1-~ G hw ko k Ke~ Ke, L M P P
Greek letters Acknowledgement---This work is supported by the National Natural Science Foundation of China, the State Education Commission of China, and the China Petrochemical Corporation.
C(
6 AH 7
NOTATION
ai, bi
coefficients in rational function
Bi
Biot number \ 2 - ~ , ]
coy c~ C d
covariance specific heat, kJ/mol K correlation matrix diameter, mm axial diffusion for mass activation energy, kcal/mol
Oea E
2 0 P
Pn 62
dimensionless pre-exponential factor coefficient of reaction rate difference from a standard value reaction heat, kcal/mol dimensionless activation energy eigenvalue of a matrix dimensionless temperature dimensionless radial position packing density of catalyst bed dimensionless axial position square deviation
Supercripts T /x
transpose of a matrix estimated value
;8
ZHEN-MIN CHENG and WEI-KANG YUAN
4800
Subscripts 0 values at z = 0 1, 2, 3 sequential number A C2H5OH B C2H40 C CO2 and H 2 0 h heat transfer rn mass transfer p particle t tube w the wall
REFERENCES
Beck, J. and Hassell, H. L., 1971, calculation of differential rates from experiments with packed tubular reactors. Proceedings 4th European Symposium Chemical Reaction Engineering, Brussels, 1968, p. 229. Pergamon, NY. Cheng, Z. M., 1993a, Measurements on the heat transfer and kinetics parameters in a fixed-bed reactor. Ph.D. thesis, East China University of Science and Technology. Cheng, Z. M. and Yuan, W. K., 1993b, Prediction of velocity distribution in a packed column. J. Chem. Ind. Engng (China) 44, 624-628. Cheng, Z. M. and Yuan, W. K., 1995, Experimental investigation on the diffusion-reaction problems in a fixed-bed reactor. Chem. React. Engng Tech. (China) 11, 173. Cheng, Z. M. and Yuan, W. K., to appear, Initial estimation of heat transfer and kinetic parameters of a wall-cooled fixed-bed reactor. Comput. Chem. Engng. Cheng, Z. M. and Yuan, W. K., to appear, Lumped kinetics of fixed-bed reactors. Prog. Natural Sci. (China). Clement, K. and J~rgensen, S. B., 1983, Experimental investigation of axial and radial thermal dispersion in a packed bed. Chem. Engng Sci. 38, 835. Coberly, C. A. and Marshall, W. R., 1951, Temperature gradients in gas streams flowing through fixed granular beds. Chem. Engng Prog. 47, 141-151. Daszkowski, T. and Eigenberger, G., 1992, A reevaluation of fluid flow, heat transfer and chemical reaction in catalyst filled tubes. Chem. Engng Sci. 47, 2245. Dimitrov, S. D. and Kamenski, D. I., 1991, A parameter estimation method for rational functions. Comput. Chem. Engng 15, 657.
Emig, G. and Hosten, L. H., 1974, On the reliability of parameter estimates in a set of simultaneous nonlinear differential equations. Chem. Engng Sci., 29, 475. Hofmann, H., 1979, Progress in modeling of catalytic fixedbed reactors. Ger. Chem. Engng 2, 258. Hofmann, H., Emig, G. and R6der, W., 1984, The use of an integral reactor with sidestream-analysis for the kinetic investigation of complex reactions. Proceedings of the 8th International Symposium on Chemical Reaction Engineering, 419. Hosten, L. H., 1974, A sequential experimental design procedure for precise parameter estimation based on the shape of the joint confidence region. Chem. Enong Sci. 29, 2247. Hosten, L. H., 1979, A comparative study of short-cut procedures for parameter estimation in differential equations. Comput. Chem. Engng 3, 117. Kamenski, D. I. and Dimitrov, S. D., 1993, Parameter estimation in differential equations by application of rational functions. Comput. Chem. Engng 17, 643. Michelsen, M. L., 1972, A least-squares methods for residence time distribution analysis. Chem. Engng J. 4, 171. Paterson, W. R. and Carberry, J. J., 1983, Fixed-bed reactor modelling: the heat transfer problem. Chem. Engng Sci. 38, 175. Peuider, V., Cerny, J. and Pasck, J., 1971, A contribution to kinetic studies performed on a nonisothermal pilot-plant reactor. Proceedings 4th European Symposium Chemical Reaction Engineering, Brussels, 1968, p. 239. Pergamon, NY. Pritchard, D. J. and Bacon, D. W., 1978, Prospects for reducing correlations among parameter estimates in kinetic models. Chem. Engng Sci. 33, 1539. Seinfeld, J. H. and Lapidus, L., 1974, Mathematical Methods in ChemicalEngineering,Vol. 3: Process Modeling, Estimation, and Identification. Prentice-Hall, Englewood Cliffs, NJ. Smith, J. M., 1981, Chemical Engineering Kinetics, 3rd Edn. McGraw-Hill, Inc., NY. Valstar, J. M., Van Den Berg, P. J. and Oyserman, J., 1975, Comparison between two dimensional fixed-bed reactor: calculations and measurements. Chem. Engng Sci. 30, 723. Wijngaarden, R. J. and Westerterp, K. R., 1992, The statistical character of packed-bed heat transport properties. Chem. Engng Sci. 47, 3125.
Chemical Enoineerin O Science, Vol. 51, No. 21, pp. 4791 4800, 1996 Copyright © 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved S0009-2509(96)00330-2 0009-2509/96 $15.00 + 0.00
S I M U L T A N E O U S ESTIMATION OF KINETIC AND HEAT TRANSFER PARAMETERS OF A WALL-COOLED FIXED-BED REACTOR Z H E N - M I N CHENG and WEI-KANG YUAN* UNILAB Research Center of Chemical Reaction Engineering, East China University of Science and Technology, Shanghai 200237, Peoples' Republic of China
(Received 10 July 1995) Abstract--Heat transfer and kinetic parameters of wall-cooledfixed-bedreactors are traditionally obtained through separate heat transfer and kinetic experiments. In spite of its simplicity,the uncertainties resulting from the low parametric sensitivity and high correlationship among the parameters inevitably make the estimation unsuitable for safe extension to actual reacting conditions. In this paper, an effort was made to estimate simultaneously both heat transfer and kinetic parameters under reacting conditions in a single tube wall-cooled fixed-bed reactor, and a two-stage parameter estimation procedure was developed. The correlation and confidence region analysis performed in this paper verifies both theoretical and experimental feasibility of simultaneous estimation of heat transfer and kinetic parameters under real reacting conditions. Copyright © 1996 Elsevier Science Ltd
Keywords: Fixed-bed reactor, parameter estimation, correlation analysis, confidence region. INTRODUCTION The difficulties encountered in fixed-bed reactor design can be ascribed to two reasons. Firstly, models describing reactor performances are usually oversimplified, particularly when the experimental ratio of dr~dr is relatively low (less than 10), leading to enormous influence of particle size on gas velocity distribution, local packing density of solids and gas-solid external transport coefficients. Besides, the randomness of catalyst packing, which causes heat transfer coefficients to differ by as many as three times under the same operating conditions for different experimental runs (Wijngaarden and Westertcrp, 1992), makes decisive modeling impractical. Secondly, the fixed-bed reactor is a very sensitive system due to the exponential dependence of reaction rate on temperature and the low thermal capacity of the system, and hencc any slight change in operating variables will result in large variations in temperature and concentration, even causing reactor runaway. Different opinions have been expressed in designing a fixed-bed reactor. Pcterson and Carberry (1983) urged that as many parameters as possible be independently measured since simultaneous estimation of too many parameters would make thc results meaningless because of their strong correlation. According to this point of view, the heat loss parameters and the reaction rate parameters should be obtained separately. This approach has the advantage of being relatively simple, since only a small number of parameters is estimated at a time. However, the reactor model must be very accurate. Clement and Jorgensen * Corresponding author.
(1983), on the other hand, insisted on the simultaneous estimation of all the parameters as a reliable attempt especially when a simplified model was used, for example, the plug-flow model. It should also be pointed out that good results cannot be achieved from separate estimations as expected. It is well known that the correlation between radial thermal conductivity and wall heat transfer coefficient is very high, as can be observed from the well-known relationship 1
1
dt
U - hw + 8Ker Strong correlation is also frequently encountered between the pre-exponential factor and the activation energy in the Arrhenius expression, e.g. as high as 0.9999, as was noted by Pritchard and Bacon (1978). Further investigation has also shown that heat transfer parameters under reactions could differ considerably from those under reaction-free conditions (Hofmann et al., 1984). It should also be noted that the low parametric sensitivity of heat transfer parameters under non-reacting conditions cannot meet the requirement in the accuracy for those under reacting conditions. In consideration of the overall inherent problems existing in separate parameter estimations, simultaneous estimation of heat transfer and kinetic parameters may be a better recourse for reactor design. EXPERIMENTAL
Determination of heat transfer parameters of a fixed-bed reactor under non-reacting conditions has been conductcd by measuring radial temperaturc
4791
4792
Z H E N - M I N CHENG a n d W E I - K A N G YUAN
profiles at different bed depths by progressively adding catalyst pellets up to different levels, as proposed by Coberly and Marshall (1951). Under reacting conditions, such a method was also used by Valstar et al. (1975). It can be imagined, however, that the method of increasing reactor length by adding catalyst would correspond to a reactor with varying length. According to Danckwerts'-type boundary conditions, the combination of all the exit concentrations and temperatures of a length-varied bed cannot represent the real situation of the whole reactor. To overcome such a drawback, a method was proposed, which is an improvement of Hofmann's (1979) contribution. The reactor, made of 4 cm o.d. stainless-steel tube, consists of six separate sections, is shown in Fig. 1. The lengths of individual sections, counted from the bottom to the top, are 10, 15, 15, 20, 20 and 30 cm, totalling 110 cm. Compared to the catalyst diameter of 3-5 mm, the reactor is long enough to make end effect negligible on axial dispersion. The reactor tube is provided with a jacket (details shown in Fig. 3) for coolant circulation in a 4 mm channel• Conduct oil (SD-320), which can resist a temperature as high as 320°C, is recycled by a high-temperature centrifugal oil pump at a flow rate of 1.7 m3/h between the jackets of the six reactor sections (in parallel, shown in Fig. 1) and an oil tank equipped with a controlled 3 kW electric heating coil. Catalyst pellets are packed only in the four middle sections (5, 15, 20 and 20 cm, totalling 60 cm) while the top and bottom sections are packed with inert ceramic packings. Between every two neighboring sections is sandwiched a specially designed composite probe (Fig. 2), each equipped with 5 gas sampling tubes and 9 thermocouples, located at different radial positions, for radial concentration and temperature measurement. The 9 thermocouples are held in place by a ceramic monolith, which also functions to suppress gas mixing. Temperatures are acquired through a computer which is also used for controlling the operating conditions. A parallel-sequential exothermic reaction, i.e. catalytic oxidation of ethanol by air over CuO catalyst supported on 7-A1203 pellets (3-5 mm in diameter), is used as the working system: (A)
P = E(AkAk T) F
~k2,
cov(ka, k2)
•.-
"-I
cov(kl, kp)[ /
= I COY(k2, kl)
322
! [_cov (kp, k0
cov(kp, k2)
•"
__I
(2) The parameter correlation matrix C, whose element C~j is defined as coy (kl, k j)
Co
6k,6~
can be easily obtained from matrix P. This knowledge, concerned only with parameter variance and parameter correlation, is, however, not enough to describe the actual parameter confidence region. It is known that for a p-dimensional parameter vector, the confidence region is a hyperellipsoid centered at the estimated value k in the p-dimensional parameter space (Hosten, 1974), and the size of the confidence region is equal to the square root of the determinant of matrix P. Complete information concerning the parameter estimation should include the following aspects: (1) (2) (3) (4)
k, the parameter estimation; 6kl, the standard deviation of ki; CO, the correlation between k~ and ks; x//~l, v~22 . . . . . x/~p, the axial lengths of the pdimensional hyperellipsoid; (5) IPI l/2, the size of the confidence region which is equal to 2x/~x2~ ..- 2p. Given the above-mentioned values, the linearized individual confidence interval for a parameter under
2
, C2H40 + (-AH1) +02
The task for parameter estimation is to find the optimal value of a parameter and its confidence region. Assume the estimated error of a parameter vector k is Ak, and its error distribution is represented by a covariance matrix P. Then the relationship between P and Ak can be established (Seinfeld and Lapidus, 1974):
(c)
(B) 1
C2HsOH
PARAMETER ESTIMATION
3
+ Oz
) CO
2 +
H20 + (-AHa).
(1)
T
-~-O2 The mechanism of this reaction is the redox type. According to the Langmuir-Hinshelwood mechanism (Smith, 1981), the rate expression under high oxygen concentration (21%) and low reactant concentration (1%) can be reduced to a simple Arrhenius first-order form. Experiments were carried out under the operating conditions shown in Table 1.
95% level according to the normal distribution is readily found to be: kl ± 1.966ki. The heat and mass balances for the two-dimensional plug-flow model are described in dimensionless forms as follows:
o¢ =~'\-Y7
p ep/+'~'~"
(3)
4793
Simultaneous estimation of kinetic and heat transfer parameters to M.S.
oil out
sampling tubes
thermocouples ,°
. . . . . . . . . . . . . ,
° . . . . . . . . . .
oil pipe
\
o,j
16 way ~ valve
converter
recorded 1 and displayed by computer
adiabatic section
recorded by mass spectrometer
inert packing I
premixed gas
Fig. 1. Overall schematic structure of the reactor and the analysis arrangements.
i I I I
thermocouple I
I
I
I
cerami, monoli
III
sampling tube
:111II I
_1
t
catalyst pellets,I reacting gas I
[1II
reacting~ ~e probe Fig. 3. Interior detailed structure of the reactor.
Fig. 2. The composite temperature and concentration measurement station.
The reaction rates in eqs (3)-(5) are expressed as
~X C
:~2X C
~0
(~2o
l ~Xc~
(4)
~ n = y,4oYo~[kl(l - Xn - X c ) - k z X n ] ~ c = YaoYo2 [(k2Xn + k3(1 - X s - X c ) ]
1 ~0\ (5)
(9) (10)
and the reaction rate constants k~, k2 and k3 are assumed to be of the Arrhenius type:
with the boundary conditions
~=0,
xB=Xc=O, eXR
p=0,
p=
~p
~Xc
(6)
o=1 ~0
- . . . . ~p ~p
(7)
0
-~-~
,
i = 1,2,3.
(8)
8p
(I1)
The dimensionless parameters involved in eqs (3)-(8) are defined as 4Ldp
~0 - - = -- B i ( O - - 1 ) .
1,
ki=kioexp
~m
="
2
d, Pe~,
4Ldp '
~h
~-
d~Peh
Table 1. Operating conditions and the maximum observed temperature rises No.
YAo(%)
Yao (%)
Yco (%)
T~, (K)
To (K)
G (kg/m 2 s)
AT=,, (K)
1 2 3
1.044 0.94 0.86
0.0 0.0 0.0
0.0 0.0 0.0
498 503 507
498 503 507
0.295 0.369 0.369
25 33 57
(12)
4794 /~1 =
ZHEN-MIN CHENG
L UCAo '
~
=
L(-AHt) GCp
-
,
/h
=
hwdt Bi = - 2Ker
L(-AH3) GCp (13) (14)
Since the pre-exponential factor ko and the activation energy E appear as a product in the Arrhenius kinetics, this will result in a correlation of the two parameters higher than 99.9% and thus leading to a serious illness in the objective function. According to the theory of regression analysis, the reaction rate constant k should be reparameterized as k = ~exp[~(1 - 1/0)]
(15)
where ~t = koexp(-E/RTr=f)
and
y = E/RTror. (16)
The reference temperature T, ee in eq. (16) is set equal to the inlet temperature To. Generation of initial values of parameters The initial value of a parameter can initially be generated either by the random search method or by an empirical relation. However, because of the high non-linearity of the process, poor initial value may lead to failure in optimization. Understandably, a good initial guess is preferably obtained from experimental data. As can be found from eqs (3)-(5), it is the reaction rates ~n and Ytc which make the three partial differential equations and thereby all the model parameters coupled together through their functional relations with 0, Xn and Xc. Therefore, ~B and 9tc have to be known a priori so as to make all the three equations independent. It seems that this requirement can only be fulfilled by employing the differential method (Hosten, 1979; Dimitrov and Kamenski, 1991; Kamenski and Dimitrov, 1993). It is observed here that if the axial and radial derivatives are substituted by their numerical values, the values of ~B and ~ c can be obtained directly from eqs (3) and (4), and the three equations involving the two heat transfer parameters and six kinetic parameters can thus be decoupled. An efficient initial parameter value generating procedure is outlined according to Cheng and Yuan (to appear):
(1) Interpolating the values of 0, XB and Xc at dimensionless radial positions p = 0.577 and 1 through profile fitting. (2) Evaluating the radial derivatives by orthogonal collocation with one collocation point. (3) Using the Pad6 rational function R(() =
ao + aa~ + a2~2 + a3~ 3 1 + bl( + b2( 2 + b 3 ( 3
to approximate the measurement of 0, Xn and Xc at p = 0.577 in the axial direction from ( = 0 to ( = 1. The axial derivative is then calculated through analytical differentiation of the rational function.
and WEI-KANGYUAN (4) Evaluating the reaction rates ~n and ~ c from the two mass balance equations (3) and (4), and then substituting them into the heat balance equation (5) to evaluate the heat transfer parameter cth. (5) Obtaining Bi from the boundary condition (8). (6) Kinetic parameters of the three reactions can finally be obtained from reaction rates ~ e and ~Rc at various values of 0, Xn and Xc through non-linear least-squares regression. Precise determination of the model parameters The Gauss-Newton-Levenberg-Marquardt method (Seinfeld and Lapidus, 1974) is employed to produce the final precise estimation. Orthogonal collocation with four interior points is used to convert the three partial differential equations into a set of ordinary differential equations. The parameter vector to be estimated is composed of two heat transfer and six kinetic parameters, and in dimensionless form, is written as
k
= (¢th, B i , ~ t , ~ l , ~ 2 , ~ 2 , ~ 3 , ' ~ 3 ) T .
(17)
There are therefore 12 state equations and 96 parametric sensitivity equations involved, and these 108 differential equations have to be integrated from ~ = 0 to ( = 1 simultaneously, for which it is found that the Runge-Kutta-Merson algorithm is suitable. RESULTSAND DISCUSSION Initial estimates With the experimental temperature and concentration data listed in Table 2, the differential method is first used to produce initial estimates of the eight parameters, with the results shown in Table 3. It is observed from Table 3 that the preliminary estimates for different conditions are in good agreement. However, the objective functions SSRs differ by orders of magnitude: under operating condition no. 3 the SSR is about 1000 times of that under condition no. 1. However, the small SSR value for condition no. 1 does not definitely mean more accurate parameter estimation, since it may probably be caused by low parametric sensitivity. Precise estimates Simultaneous parameter estimation of heat transfer and kinetic parameters has been touched in the literature (Beck and Hassel, 1971; Peuider et al., 1971; Emig and Hosten, 1974; Clement and Jorgensen, 1983; Hofmann et al., 1984). In Clement and Jorgensen's work, the set of partial differential equations of the system were transformed into algebraic equations by double collocation, resulting in 2 and 15 collocation points in the radial and axial directions, respectively. In order to improve convergence, a parameter transformation method suggested by Michelsen (1972) was employed. It is obvious that double collocation has serious limitations in its application, because for example, if 4 radial points and 15 axial points are adopted, 180 algebraic equations will result, along with 180 x 180
Simultaneous estimation of kinetic and heat transfer parameters
4795
Table 2. Experimental data of temperature and concentration at different axial and radial positions (1) Operation no. 1 Temperature (K)
Concentration
zi\ri
rl
r2
r3
r4
r5
r6
r7
r8
r9
zl z2 z3 z,,
506 516 522 510
506 516 522 510
506 514 520 509
505 513 520 509
505 512 519 507
504 511 516 506
504 5ll 515 505
503 509 512 504
503 506 509 504
Zi\ri
YA
Zl
YB YC
z2
YA YB Yc YA
Z3
YB YC
Z4
Ya yn YC
(mol%)
rt
r2
F3
r4
F5
0.9672 0.0616 0.0197
0.9608 0.0623 0.0199
0.9644 0.0618 0.0198
0.9677 0.9693 0.0594 0.0574 0.0189 0.0182
0.6434 0.2736 0.1261
0.6463 0.2744 0.1260
0.6525 0.2678 0.1208
0.6818 0.7008 0.2513 0.2420 0.1088 0.1024
0.3326 0.4031 0.3047
0.3377 0.4039 0.3034
0.3425 0.4019 0.2933
0.3795 0.3924 0.3960 0.3925 0.2722 0.2615
0.2132 0.4159 0.4160
0.2101 0.4161 0.4156
0.2146 0.4169 0.4108
0.2218 0.2308 0.4182 0.4188 0.4005 0.3952
(2) Operation no. 2 Temperature (K)
Concentration (mol%)
z~kr~
rl
r2
r3
r4
rs
r6
r7
rs
r9
zl z2 z3 z,
512 531 537 524
512 528 537 523
512 527 534 521
511 525 532 520
511 523 530 519
510 521 528 516
510 519 523 513
508 517 520 512
507 514 516 510
zi\ri
Z1
Ya YB YC
Z2
YA yn
Z3
YA Ya
YC
YC
YA Z4
YS YC
rl
r2
r3
r4
r5
0.8648 0.0541 0.0175
0.8680 0.0547 0.0177
0.8700 0.0544 0.0176
0.8706 0.8748 0.0525 0.0509 0.0169 0.0163
0.5816 0.2410 0.1128
0.5865 0.2417 0.1127
0.6002 0.2359 0.1079
0.6290 0.6375 0.2214 0.2131 0.0971 0.0912
0.2819 0.3582 0.2927
0.2900 0.3594 0.2907
0.3013 0.3579 0.2787
0.3382 0.3418 0.3528 0.3497 0.2544 0.2422
0.1753 0.3622 0.4043
0.1783 0.3626 0.4036
0.1766 0.3640 0.3976
0.1848 0.1963 0.3665 0.3677 0.3850 0.3787
(3) Operation no. 3 Temperature (K)
Concentration (mol%)
zi'kri
rl
r2
r3
ra
rs
r6
rv
rs
r9
zl z2 z3 z,
517 540 565 526
517 538 565 524
517 535 562 524
516 532 555 523
515 530 553 520
513 527 550 518
513 524 543 515
511 521 534 514
509 516 523 511
zi\ri
Z,
YA YB YC
z2
YA Ys Yc
z3
Ya ya Yc
z4
ya YB Yc
rl
r2
r3
r4
rs
0.7745 0.0548 0.0175
0,7699 0,0525 0,0178
0.7803 0.0491 0.0166
0.7824 0.7800 0.0469 0.0450 0.0158 0.0151
0.5057 0.5104 0.2276 0.2236 0.1301 0.1293
0.5381 0.2021 0.1117
0.5679 0.5811 0.1956 0.1862 0.0995 0.0969
0.1311 0.2960 0.4324
0.1368 0.2972 0.4309
0.1910 0.3028 0.3632
0.2444 0.2613 0.3019 0.3085 0.3175 0.2971
0.0690 0.2407 0.5513
0.0694 0.2527 0.5499
0.0838 0.2534 0.5221
0.0991 0.0104 0.2626 0.2669 0.4963 0.4907
(4) Positions of measurement probes Axial positions of sampling tubes and thermocouples: zj = 5, 20, 40, 60 cm from the inlet. Radial positions of the 5 sampling tubes: r~ = 0, 3, 7, 13 and 18 mm, where 0 mm is the center of the bed. Radial positions of the 9 thermocouples are listed in the underlying table in "ram", where 0 mm is the centre of the bed.
ZHEN-MIN CHENG and WEI-KANGYUAN
4796
Table 2. (contd.)
zi\r~
rl
r2
r3
r,
rs
r6
r7
rs
r9
z~ Z2 Z3 Z4
2.5 3.1 2.4 3.3
5.1 6.0 3.8 4.9
7.2 7.4 7.2 8.2
9.1 10.6 8.8 9.5
11.0 11.7 10.2 11.3
12.9 13.4 12.0 13.5
14.0 14.2 14.2 14.8
15.1 15.0 15.4 15.9
16.1 16.7 16.9 16.5
Table 3. Preliminary results of parameter estimation Operating condition
cth
Bi
0tI
Yl
~2
~)2
0t3
Y3
SSR
l 2 3
3.96 3.45 3.20
3.53 3.19 3.88
- 1.82 -1.54 - 1.26
27.47 26.52 21.88
- 2.27 -1.94 - 1.71
23.28 21.40 20.37
- 2.97 -2.82 -2.49
27.20 31.73 27.60
0.00036 0.0058 0.3588
Table 4. Precise estimation according to Gauss-Newton-Levenberg-Marquardt Operating condition
~h
Bi
~q
71
~2
Yz
a3
73
SSR
1 2 3
3.24 3.14 2.67
5.18 3.99 4.99
- 1.81 - 1.43 -1.40
26.56 24.08 22.80
-2.51 - 2.09 -1.63
30.04 25.05 20.78
-2.91 - 2.60 -2.61
23.84 26.13 28.15
0.00004 0.00003 0.0005
working matrices. In comparison, the largest matrix involved in this work is only 12 x 8. Actual iteration shows that SSR can be reduced to a relatively small value after only about 4 iterations, and the final result can be achieved within 7 to 9 iterations (see Table 4).
Statistical properties of estimated parameters Individual confidence limits. The deviation of the estimated values from the optimal estimates is the most desired information concerning the accuracy of parameter estimation. Standard deviation of parameter ki, defined as 6ki, can be readily obtained from the square root of the ith diagonal element of P. The individual confidence region at the selected probability level of ( 1 - ~) is given by the well-known relation (Seinfeld and Lapidus, 1974) ;,, - l~rk, <. k* <~ i, + [~6~,
(18)
where/~ depends on the confidence level. At the 95% level/~ = 1.96. The individual confidence region for each dimensionless and physical parameter is shown in Table 5.
Joint confidence region. By assuming that the model can be adequately linearized in the vicinity of the parameter estimate k, and if 62 is not known, the joint confidence region can be defined by F-distribution: (~ _ k , ) T p - 1 ( ~ _ k*) = 62pF~ _,,(p, m - - p ) (19) where ~ 2 = S S R m i , / ( m - p). One notices that there will be p variables on the left-hand side subject to the
same eq. (19). In this work, the hyperellipsoid confidence region cannot be represented by the solution of eq. (19) since the parameter number is eight. Instead, the confidence region can be approximately described by the length of each axis as well as the size of the confidence region (see Table 6). The relative ratio of the longest to the shortest axis ranges from 1118 to 14,425, indicating an extremely elongated hyperellipsoid. It shows further that the shape is more elongated under operating condition no. 1 than under conditions nos 2 and 3. Besides, the size of the confidence region under the no. 1 condition is larger by 10, 000 times than that under no. 3 condition. These findings are promising, since a precise estimate can thereby be obtained under sensitive conditions like condition no. 3. Further description of parameter estimation should be referred to the correlation matrix, as can be found in Table 7. The correlation matrices under the first two conditions are similar. The correlation between the two heat transfer parameters ~h and Bi is about 0.8, a substantial reduction compared with 0.95 when no reaction occurs. The correlation between the pre-exponential factor ~i and activation energy ?i of the same reaction is about 0.95, larger than that between different reactions. The most exciting result is observed from the correlation coefficient between one of the six kinetic parameters, ~1, 71, (z2, ?2, 0c3 and Y3, and one of the two heat transfer parameters, ~h and Bi, with an average value of about 0.3, the largest values being 0.488 for condition no. 1 and 0.517 for condition no. 2. F o r condition no. 3, the correlation between the two heat transfer parameters, cth and Bi, has decreased
Simultaneous estimation of kinetic and heat transfer parameters
O +1 -H +1 ",~" eq '~1"
-H -H +1 o~tt~
eq eq ¢-q
¢.q t-,q ¢-.1
X
X
X
-H -H -H O
¢-q ¢-q ¢-q
I
+1 +1 +1
4797
to 0.2, but positively correlated, and additionally the correlation between the kinetic parameters has also decreased considerably. Another important finding is that the average value for the correlation between one of the six kinetic parameters ~a, ~1,0~2,~2,5 3 and ~'3, and one of the two heat transfer parameters, ~h and Bi, still remains at 0.3. Although the reactor behaviors under the three conditions are different in many respects, one thing however is similar: the correlation between a kinetic parameter and a heat transfer parameter remains almost unchanged at an average value of 0.3. This value implies that the correlation between the kinetic and the heat transfer parameters is weak, and therefore simultaneous estimation of transport and kinetic parameters is theoretically feasible.
O~ tt% 1"--
+1 +1 -H
-I-I +1 -H ee~
~,'5 ¢N exl
o~,o
.~- t-q ¢'q
X
X
X
+1 +1 +1 ¢'q ¢'-,i ~
+1 +1 +1
I I I
e,q
e.
e-q
+l +l +l
+1 + l + t
8
i'-I o ¢'q ¢-q eq
I
w; o
~
¢.i i~.1 ¢-i
t~ eO
X
O
+1 +1 +1
X
X
+l +l +1 I
[--
7771
• r..: q-~ t"q ~ O +1 +1 +1
e~ +1 + l +1
~'~1 .-.-~ tt~ ¢'.q + i +1 +1
+1 +1 + l p... ~ tt% "2. "n: 'Q I
e6 ,A e.i e-,
o
..
Comparison with reactor performance Gas velocity and bed voidage fluctuate across the tube diameter; however, the measured radial concentration and temperature profiles are smooth, as can be found in Valstar et al. (1975), Hofmann (1979) and Daszkowski and Eigenberger (1992), as well as in our experiment• Perhaps local fluctuations in temperature and concentration have been damped by mass and heat diffusions• Therefore, when designing wallcooled fixed-bed reactors, the two-dimensional plugflow model permits a reasonable reactor performance description• A comparison is made in Figs 4 and 5 between plug-flow model predictions and experimental findings. It is found that fitting of concentration is very good (Fig. 5), but for temperature (Fig. 4), a substantial disagreement is observed within a distance of about one particle diameter to the wall. It is believed that such a result is mainly because of the different influence of the wall on the heat and mass transport process• The existence of the wall is advantageous in damping concentration differences since it confines the mass transport within the tube. Therefore, in most cases, the concentration profile is usually relatively flat or S-shaped, as shown in Fig. 5. On the contrary, the influence of wall on heat transport is different, causing a decreased thermal conductivity, a higher bed voidage and a shorter gas-solid contacting time in the 'wall region'• All these effects combine to make the radial temperature profile near the wall steeper than the plug-flow prediction, probably because the plug-flow simplification assumes a constant radial effective thermal conductivity and a uniform voidage even up to the wall (Cheng and Yuan, 1993).
o CONCLUSIONS
O
(1) A two-step parameter estimation method is presented in this work. In the first step, an initial value is generated from experimental data without making any assumptions, and in the second step, the Gauss-Newton procedure is employed for further precise estimation. Although 108
ZHEN-MIN CHENG and WEI-KANG YUAN
4798
Table 6. The confidence region under different operating conditions Operating condition 1 2 3
~
~
~33
~
x/~5
x~6
~
10.0 5.87 1.54
3.50 1.97 1.03
2.12 1.36 0.84
1.14 0.52 0.13
0.085 0.058 0.038
0.031 0.019 0.011
0.0041 0.0076 0.0034
x~s
IPI
0.0023 0.0053 0.00058
2.06x 10 -6 3.53 x 10-7 1.43 x 10- to
Table 7. Correlation matrix for all the eight parameters under different operating conditions (1) Operating condition no. 1 a~, 1.00000
ah Bi
--0.81915 -0.40933 0.48811 -0.31009 0.36463 0.05445 -0.08988
at 71 a2 72 a3 ?3
Bi
at
Yl
a2
?2
a3
Y3
1.00000 0.35096 -0.21491 0.31197 -0.24795 -0.14301 0.23354
1.00000 -0.91410 0.60761 -0.65237 -0.70207 0.68397
1.00000 -0.54895 0.68487 0.61478 -0.58653
1.00000 -0.96408 -0.71644 0.70475
1.00000 0.73865 -0.72951
1.00000 -0.97794
1.000~
(2) Operating condition no. 2 ah Bi
al 71 a2 ?2 a3 73
ah
Bi
at
?t
a2
?2
a3
1.00000 --0.75367 -0.40692 0.51703 -0.36376 0.44763 0.13693 --0.16799
?3
1.00000 0.32876 -0.18089 0.42442 -0.33805 -0.25147 0.37214
1.00000 -0.90719 0.64322 -0.70664 -0.76505 0.72490
1.00000 -0.55108 0.72031 0.64709 -0.58556
1.00000 -0.95468 -0.76995 0.76226
1.00000 0.78753 -0.76587
1.00000 -0.97317
Bi
at
Yz
a2
?2
~+
1.00000 0.22007 0.20241 0.60427 -0.47872 -0.28875 0.35252
1.00000 -0.69684 0.47067 --0.59403 -0.75131 0.52766
1.00000 0.05695 0.32318 0.36685 -0.20584
1.00000 -0.89146 -0.67678 0.34740
1.00000 0.73943 -0.37199
1.00000 -0+83327
1.00000
(3) Operating condition no. 3 ah 1.00000
ah Bi
0.20189 -0.11682 0.64973 0.14422 0.23299 -0.19265 0.45371
al ?t a2 Y2 a3 ?3
?3
1.00000
,540 "Z=SCm ÷ Z = 2 0 ~ 530
""z = 5 m -='z ÷Z =40~ +Z== 20cm WelII
-*z = 4Oan + z = 60cm
i
; 5 3 0 ............
.... 51 (~ ......~
"
". .
I
I
3
6
"*Z=SCO1 ÷z=2Ocnl
5 7 0 f. ........
* z = 40cm + z = 60crn
s6o
...
....
.
"'..
{
550
s,o
iii'. l
3o..
[
i\" .....
....
, ....
, 510]
0
)
9
12
15
18
r, m m
A: operating condition No. 1
5000
. . . .
~
~
;
12
15
18
r, mm
B: operating condition No.2
0
3
6
9
12
15
r, mm
C: operating condition No.3
Fig. 4. Radial temperature distributions: difference between measurement and plug-flow prediction.
18
Simultaneous estimation of kinetic and heat transfer parameters • 1=5¢m
4799
nzm2Ocnl
,z=Scm
• zm4OCnl • z-8Ocm O.B
0.41-
" Z=SCm
0.1~
~ Z = 2(] ¢1~
0.8
o at
* Z = 40cra " Z = 60Cm
xz=2Oom
o z : 4 O c r l t v z:l~O(:Im
_:. ___---~
06 ~ - ' ~ - ' - - ~ "
--..____~_~
o.21
o.4
~
0.1
0.2
.
•
0.4:
t
0.2 --
~t I
% a
,t
.t I
6
I
9
,¢--
I
I
12
15
11m
O.
0
.
.
3
.
6
r, m m
.
.
9
~
12
15
1
0
~t
3
r, m m
[ .
[
6
9
2&
1
5
1
r, m m
Fig. 5. Radial concentration distributions (operating condition no. 3).
differential equations have to be integrated simultaneously when four collocation points are used, it is still an easy task using a personal computer, since CPU storage and processing time can be greatly saved as compared to other methods. (2) Since the parameter space is highly elongated, a variety of gradientless methods, such as Simplex or Powell, are less effective in minimizing the objective function. (3) A positive answer can be given to the availability of simultaneous estimation of heat transfer and kinetic parameters. Mathematically, the average correlation between heat transfer and kinetics parameters is only about 0.3. On the other hand, they are physically different: the heat transfer parameters are related to physical properties of the fixed bed as well as the catalyst and the fluid, whereas the kinetic parameters are only related to the chemical reaction, (4) Parameter estimation performed under reacting conditions shows that the size of confidence region under conditions of different parameter sensitivity can differ up to 10,000 times, which implies that it is necessary to do parameter estimation under relatively sensitive conditions.
1 ct point of the F-distribution mass velocity, kg/m2s wall heat transfer coefficient, J/m 2 s K pre-exponential factor parameter vector axial dispersion for heat effective thermal conductivity, J / m s K catalyst bed height, m number of experimental data points numbers of parameters covariance matrix of parameter Peclet number for radial heat transfer Peh Peclet number for radial mass transfer Peru radial position, m or mm r gas constant, cal/mol K R reaction rate, kmol/kgcat hr .:~ sum of squares of residues SSR t l -(~/2) t value at 1 - ~ confidence level temperature, K T linear velocity of reactant, m/s u heat transfer coefficient, J/m 2 s K U measurement covariance matrix V yield of a product X mole fraction y axial position, m z
F1-~ G hw ko k Ke~ Ke, L M P P
Greek letters Acknowledgement---This work is supported by the National Natural Science Foundation of China, the State Education Commission of China, and the China Petrochemical Corporation.
C(
6 AH 7
NOTATION
ai, bi
coefficients in rational function
Bi
Biot number \ 2 - ~ , ]
coy c~ C d
covariance specific heat, kJ/mol K correlation matrix diameter, mm axial diffusion for mass activation energy, kcal/mol
Oea E
2 0 P
Pn 62
dimensionless pre-exponential factor coefficient of reaction rate difference from a standard value reaction heat, kcal/mol dimensionless activation energy eigenvalue of a matrix dimensionless temperature dimensionless radial position packing density of catalyst bed dimensionless axial position square deviation
Supercripts T /x
transpose of a matrix estimated value
;8
ZHEN-MIN CHENG and WEI-KANG YUAN
4800
Subscripts 0 values at z = 0 1, 2, 3 sequential number A C2H5OH B C2H40 C CO2 and H 2 0 h heat transfer rn mass transfer p particle t tube w the wall
REFERENCES
Beck, J. and Hassell, H. L., 1971, calculation of differential rates from experiments with packed tubular reactors. Proceedings 4th European Symposium Chemical Reaction Engineering, Brussels, 1968, p. 229. Pergamon, NY. Cheng, Z. M., 1993a, Measurements on the heat transfer and kinetics parameters in a fixed-bed reactor. Ph.D. thesis, East China University of Science and Technology. Cheng, Z. M. and Yuan, W. K., 1993b, Prediction of velocity distribution in a packed column. J. Chem. Ind. Engng (China) 44, 624-628. Cheng, Z. M. and Yuan, W. K., 1995, Experimental investigation on the diffusion-reaction problems in a fixed-bed reactor. Chem. React. Engng Tech. (China) 11, 173. Cheng, Z. M. and Yuan, W. K., to appear, Initial estimation of heat transfer and kinetic parameters of a wall-cooled fixed-bed reactor. Comput. Chem. Engng. Cheng, Z. M. and Yuan, W. K., to appear, Lumped kinetics of fixed-bed reactors. Prog. Natural Sci. (China). Clement, K. and J~rgensen, S. B., 1983, Experimental investigation of axial and radial thermal dispersion in a packed bed. Chem. Engng Sci. 38, 835. Coberly, C. A. and Marshall, W. R., 1951, Temperature gradients in gas streams flowing through fixed granular beds. Chem. Engng Prog. 47, 141-151. Daszkowski, T. and Eigenberger, G., 1992, A reevaluation of fluid flow, heat transfer and chemical reaction in catalyst filled tubes. Chem. Engng Sci. 47, 2245. Dimitrov, S. D. and Kamenski, D. I., 1991, A parameter estimation method for rational functions. Comput. Chem. Engng 15, 657.
Emig, G. and Hosten, L. H., 1974, On the reliability of parameter estimates in a set of simultaneous nonlinear differential equations. Chem. Engng Sci., 29, 475. Hofmann, H., 1979, Progress in modeling of catalytic fixedbed reactors. Ger. Chem. Engng 2, 258. Hofmann, H., Emig, G. and R6der, W., 1984, The use of an integral reactor with sidestream-analysis for the kinetic investigation of complex reactions. Proceedings of the 8th International Symposium on Chemical Reaction Engineering, 419. Hosten, L. H., 1974, A sequential experimental design procedure for precise parameter estimation based on the shape of the joint confidence region. Chem. Enong Sci. 29, 2247. Hosten, L. H., 1979, A comparative study of short-cut procedures for parameter estimation in differential equations. Comput. Chem. Engng 3, 117. Kamenski, D. I. and Dimitrov, S. D., 1993, Parameter estimation in differential equations by application of rational functions. Comput. Chem. Engng 17, 643. Michelsen, M. L., 1972, A least-squares methods for residence time distribution analysis. Chem. Engng J. 4, 171. Paterson, W. R. and Carberry, J. J., 1983, Fixed-bed reactor modelling: the heat transfer problem. Chem. Engng Sci. 38, 175. Peuider, V., Cerny, J. and Pasck, J., 1971, A contribution to kinetic studies performed on a nonisothermal pilot-plant reactor. Proceedings 4th European Symposium Chemical Reaction Engineering, Brussels, 1968, p. 239. Pergamon, NY. Pritchard, D. J. and Bacon, D. W., 1978, Prospects for reducing correlations among parameter estimates in kinetic models. Chem. Engng Sci. 33, 1539. Seinfeld, J. H. and Lapidus, L., 1974, Mathematical Methods in ChemicalEngineering,Vol. 3: Process Modeling, Estimation, and Identification. Prentice-Hall, Englewood Cliffs, NJ. Smith, J. M., 1981, Chemical Engineering Kinetics, 3rd Edn. McGraw-Hill, Inc., NY. Valstar, J. M., Van Den Berg, P. J. and Oyserman, J., 1975, Comparison between two dimensional fixed-bed reactor: calculations and measurements. Chem. Engng Sci. 30, 723. Wijngaarden, R. J. and Westerterp, K. R., 1992, The statistical character of packed-bed heat transport properties. Chem. Engng Sci. 47, 3125.