- Email: [email protected]
Dynamic Modelling and Simulation of an Industrial Amidation Reactor
Copyright © IFAC Dvnamics and Control of Chemical Reactors (DYCORD+ '89). Maastricht. The Netherlands. 1989
DY~A~IICS
OF
CO~T1NUOUS
REACTORS
DYNAMIC MODELLING AND SIMULATION OF AN INDUSTRIAL AMIDATION REACTOR Y. Toure*, J. Biston*, C. Darmet*, G. Gilles* and J-F. Thierry** *Laboratoire d'Automatique et de Genie des Proctides, Unit'ersite Lyon I , Villeurbanlle, France **RhOne-Poulenc Recherches, Saint-FOIlS, France
Abstract. In order to get a better control of an amidation industrial process, a dynamic modelling study of the plug-flow reactor has been carried out. An infinitesimal phenomenological balance leads to a set of coupled partial differential equations. The chemical s t a ti c part of the model i s got from an analogy with respect to a batch reactor kinetics. Identifiability results of this non linear model are presented. From the non linear distributed parameter dynamic knowledge model, a static model is drawn out and leads to the writing of a static simulation package in order to find the static profiles and to choose the operating points. Keywords. Textile industry modelling ; identification
chemical industry simulation.
INTRODUCTION
distributed parameter systems
length in which circulates a thermal fluid. These sectlons a re thermally independent from each other. Each one is itself a combination of parallel flow and counterflow heat exchangers, displayed in the reactional mixture and in the double shell of the reactor. A typical section is represented on Fig. 1.
Getting a better quality and a better regu larity in the production of synthetic threads implies a very good control of the amidation plants on the upper part of the manufacturing process. These plants being constituted of a set of rather complex apparatus, a very good knowledge of the static and dynamic behaviours of each of these unitary processes is necessary for the design of a mov"e performant control system. Among these apparatus, the amidation reactor takes a prominant part because it contributes the most to the ex tent of the reaction. This evaporator reactor is a tubular heat exchanger in which circu l ates a thermal fluid and which is the center of coupled thermal and chemical phenomena. An infinitesimal phenomenological balance allows to get a dynamic model under the form of a set of coupl ed partial differential equations. One can derive a static model giving the characteristic variables profiles along the reactor.
A constant flow-rate of a homogeneous reaction mixture composed of water, diamine and diacid is heated in the reactor. The water of the mixture is evaporated while the polyamidation reaction progresses. The kinetic model of this chemical reaction depends upon the water concentration of the medium. KINETIC MODELLING Amidation reaction The polyamide 66 (Nylon 66) can be prepared by linear condensati on of a diamine, the hexamethylene diamine (HMDA) and a diacid, the adipic acid. The global reaction can be represented as follows
After a description of the industrial amidation reactor, the reactional part of the static model is analysed by considering the analogy existing with the reaction kinetics inside a batch reactor. This method presents th e advantage of allowing easier realizations of the experiments on a laboratory batch reactor (Rhone-Poulenc Company), and so of allowing "a validation of the static chemical model without any experiments on the plant. Identifiability results are presented.
n H2N - (CH2)G - NH
+
n HOOC - (CH6)'- COOH
HUm - (CH2)6 - NHCO - (CH2h - COJ
~
OH + (2n-1)HlO
n A simplified approach is to consider only functional groups in the reaction medium : end groups - COOH and - NH" intramolecular groups - CONH and molecules of water, leading to the following representation : - COOH + - NH,
The following part is devoted to the elaboration of the global dynamic model ; then, the static model of the industrial reactor is used in order to crea te a simulation tool which will help a process operator to easily tune the operating points. Then it is shown how to use the dynamic model in order to build a control model which is necessary for the design of a new process control system.
-+ +-
- CONH -
l
+ H2 O
or, in a symbolic form A (acid)
+
B
(base)
1 -+ +-
Y
+
W
l
(amide)
(water)
An equilibrium can be reached, when the polycondensat ion reaction (1) is balanced by the hydrolysis reaction (2). Kinetics
PROCESS DESCRIPTION
From the two following assumptions - all reacting species A, Band Y have the same reactlvlty (Flo r y's equireactivity" assumption) which does not depend on the degree of polymeri-
It is an isobaric continuous reactor. On a thermal point of view, it" can be considered as aheat exchanger composed of several heating sections of equal
III
112
Y. Tour<' fI al.
zation of the molecule - the reaction obeys a first order kinetics with respect to each reacting species (overall second order kinetics), a kinetic equation can then be derived : kl A B - k z Y W Allowance must be made for the highly non ideal character of the medium, especially due to the dependence of activity coefficients on the water content (Steppan, 1987). The general expression of rate "constants" kl and kz is given by : r
~
fla , b , W) exp( - Ei/RT) (i ~ 1,2) i i i The six parameters (ai, bi, Ei) of the model have been estimated from experimental data obtained on a batch laboratory scale reactor. From these expek
~
riments, temperature, pressure and concentrations
were chosen in a range consistent with industrial corresponding conditions. The Marquardt method (Marquardt D.W., 1963 ; Lawson C. L., 1977) was used for parameter identification from the dynamic model of the batch reactor. The sensitivity analysis developed in the appendix allowed to reduce the number of the model parameters inside the considered working section. Fig. 2 shows the comparison between this model and the measurements (+) on the batch reactor.
In this part, we consider all the physical aspects which surround the chemical reaction. One of the results of the chemical kinetics modelling is the quasi non thermality of the reaction . This will allow some further simplifications. Phenomenological balance The main assumption is considering the reactor as an isobaric plug flow reactor. The infinitesimal mass and energy balance are then set inside each section Ri(i ~ 1 to n), considering only one spatial variable along the flow axis of the reaction mixtu r e. Any diffusion along this axis is neglected (J. Villermaux, 1982 ; L . Borel, 1984 ; H.S . Fogler, 1981) .
One of the main phenomena is vaporisation . By means of the energy contribution, the water is mainly evaporated inside the medium, which leads to taking into account the water vaporisation enthalpy which depends upon temperature and pressure : Lv (T p ' Pr) (J. Ficini, 1977 ; M. Bertin, 1982 ; L. Borel, 1984) . So, the mixture flow-rate conti nuously varies along the section Ri and along the whole reactor. This phenomenon is modelled by the medium water flow-rate balance : dF dz
dx
1 [r ___
e
dZ
xe - 1
(1)
x (z,t) is the mixture water mass fraction and Fe(z,t) denotes the mass flow-rate of the reaction m~xture . By means of the reaction factor r , we then get a coupling of the chemical reaction with the thermal aspect because the vaporisation energy can locally be written as dF E
v
~ - L (T v
p'
P) r
........E. dz
dx
aT ........E.
e
dZ
3z
(3)
The equations (1), (2) and (3) describe phenomena that are supposed to be instantaneous (mass conser' vation, equi libr ium) . The main dynamic aspects especially occur on the level of thermal exchanges and of the mass balance for the main reaction. We then get a set of coupled partial differential equations. The heat exchanger We showed (Fig. 1) that the heating system of each section Ri can be decomposed into simple heat ex changers involving only heat exchanges with the rea< tion medium but not between themselves. The heat exchange with the medium is mainly processed by convection. For the reaction mixture, we then ge t the following energy balance equation : aT k sdp(T - T ) + k s (T+ - T Sp Cp Ppp d dp P ep ep e F
af
GLOBAL DYNAMIC MODEL
........E.~
This state function, due to Ogata (Ogata, 1960) h, been reevaluated among Rhone-Poulenc research grot in order to take into account the reaction extent. It relies the medium water concentration with the temperature and the pressure, which allows to writ t the variation of xe(z, t) under the following ex pression :
(2)
By mass conservation, the flow-rate variation of the liquid reaction mixture is equal, in terms of absolute value, to the vapor flow-rate variation . This is due to the fact that the mixture flow-rate is constant at the entrance of the reactor. We consider that this flow-rate variation is processed in order to satisfy the liquid-vapor equilibrium law of the medium water . Thus we introduce the state function of the water-polymer mixture.
+ k
3T s (T - T ) - F C __ P _ E ep ep e p p pp 3z v
(4) The coupling with the chemical reaction is here explicit by means of the Ev term «1) and (2» . For the heating elements, we get the classical equations of the parallel flow or counterf l ow heat exchangers with constant flow rates (Sacadura, 1980 J .C. Friedly, 19 72). Thus, for the R'1 section, we get
-
3T e
at aT
3T v
e
+
e
e
az- -
h (T e e
- Tp ) (5)
3T+
aze -
h (T+ - Tp) e e
at
- v
3Td
d YV d -3z - h d (T d - Tp )
at
e
aT
In the last equation, the Y factor (0 < Y <1) is a distribution coefficient of the double shell energ y between its lower part filled with the liqu id reaction mixture, its upper part (vapor) and the heat losses. Mass balance It is a balance according to the Fluid Mechanics concepts (Vi ll ermaux, 1982) . The vaporisation pheno menon occurs in the balance and we consider that there is a small hexamethylene diamine loss (B) (a fraction of the total vapor flux, with a « 1) 3A
at
3B
3B
at
- r - PpSp Fp
3W
3W
at
a; -
r - PpSp Fp
B PpSp W
az - PpSp
3F
2
3z 3F
........E. 3z
Dvnamic Modelling and Simulation of an Industrial Amidation Reactor (6)
r -
THE STATIC SINULATION TOOL The static model, derived from equations (1) to (6) is written by considering the steady state behaviour of the process (the partial derivative with respect to time is then equal to zero in the equations (4) to (6)). The aim of this simulation tool is double rify the model with respect to the steady behaviour, then and especially to predict tem output variables in order to tune the points. This is schematized on Fig. 3.
: to vestate the sysoperating
A previous work of heat exchange parameter estimation on each section Ri is carried out, taking into account their dependency upon the operating point, and also of estimation of physical parameters which depend upon the thermal fluid and which have been set as a constant in each section. Solving the set of non linear differential equations obviously implies the knowledge of the boundary conditions. These conditions depend upon each Ri section. On a typical structure shown on Fig. 1, the boundary conditions are : T- (1) e
T ei
T~(O)
T+ (0) e
T+ (1) e
Td (1)
T (i) p (0)
T (i-1) p (1)
(input of the heating fluid)
(7)
113
into account the residence of the mixture in the whole reactor, and to tune the operating points. Fig. 4 shows the simulation result for the section of Fig. 1. The measurement points are located at the section boundaries. A sample at the output of the section allowed the chemical model validation (Fig. 4-b). The temperatures (Fig. 4-a) are normalized with respect to the nominal input temperature of the thermal fluid Tei : Tj becomes Tj/Tei. The mass concentrations A , B , Y , Ware expressed in terms of equivalent by kilogram. CONCLUSION AND PROSPECTS A dynamic model of an industrial amidation reactor has been presented as a set of coupled partial differential equations. From this model, a static model has been derived, allowing to build a simulation tool which elaborates the profiles along the reactor and which can give an off-line help and eventually an on-line help to the choice of the reactor operating points. The dynamic model, presently analyzed, can be reduced by means of a transformation of some partial differential equations into ordinary differential equations, in order to obtain a more realistic con-
trol model. By means of intermediate measurements along the reactor, specific distributed parameter control methods (Pohjolainen, 1982) can be used. From an approximated, entirely differential, model we can also apply more classical methods of optimal control under constraints with observer. A comparative analysis must point out the efficiency of the new control systems with regard to the existing system. REFERENCES
This implies that the simulation initial conditions - for example T~(O), T~(O), Td(O) must in fact be boundary condltions because of the cascade structure of the heating elements of the Ri section in the real case. The simulation software must iteratively solve the so defined problem with limits, rapidly converge and be robust with respect to initialization. The numerical solving algorithm which has been chosen is the Newton-Gauss one (J.E. Dennis, 1983).
Steppan, D.D., M.F. Doherty and M.F. Malone (1987). A kinetic and equilibrium model for Nylon 66 Polymerisation. Journal of appl. polymer science vol. 33, 2333-2344. Marquardt, D.W. (1963). An algorithm for leastsquares estimation of nonlinear parameters. SIAM J., 11, 431-441. Lawson, C.L., R.J. Manson (1977). Solving least squares problems. Prentice-Hall, Inc. Villermaux, J. (1982). Conception and functioning of reactors, Lavoisier, Paris.
This leads to solving a system tions : find the x vector such is written under the form : f(x) = RT(x) R(x) that we have example, i~ the-case of system xT
of non linear equaas f(~) = 0., f(~) to minimize. For (7) :
IT - ( 1), T- (0), T+ (0), Td ( 1), T . (0 )] Lee e pl Tpl.(0) - TP'. 1(1) T- (1) - T . e el
R(x) [
T-(O) e
J
T+(O) e
«1) - T (l) d The initialization problem concerns T .(0), T:(O) and Td(O). The minimization algorithmP~s : ~+ = x_ -
D(~+)
T
J(~+)J
-1
J(~+)
T
R(~+)
for the Ri section (i = 1 to n). ~+ and J(~+) (Jacobian of f(x) for x = x+) are numerically evaluated by means-of solving-differential equations (J.P. Nougier, 1987 ; K. Arbenz, 1980 ; J. Friedly, 1972). For the realization, "A.C.S.L." simulation package has been used on I.B.M. - P.C./A.T. computer. This method gives good results, converging after a maximum of 4 iterations. This allows to propose a tool for help to decision if one takes
Borel, L. (1984). Thermodynamic and energetics. Presses Poly techniques Romandes. Suisse Fogler, H.S. (1981). Chemical Reactors. ACS Symposium series 168. Ficini, J., N. Lumbroso, J.C. Depsay(1977). Thermodynamic and chemical equations. Herman. Paris Bertin, M. J.P. Faroux, J. Renault (1982). Thermodynamics. Bordas. Paris Sacadura,J.F. (1980). Initiation to heat transfert. Tech & Doc. Paris. Friedly, J.C. (1972). Dynamic Behaviors of processes. Prentice-Hall, New Jersey. Ogata, N. (1960). Studies on polycondensation reactions of Nylon salt. Makromolekulare Chem. 42, 52. Dennis, J.E., R.B. Schnabel (1983). Numerical methods for unconstrained optimization and nonlinear equations. Prentice Hall, New Jersey. Nougier, J. (1987). Numerical computation methods. Masson. Paris. Arbenz, K., A. Wohlhauser. (1980). Numerical Analysis. Presses Poly techniques Romandes. Suisse. Pohjolainen; S.A. (1982). Robust multivariable PIController for infinite dimentional systems. I.E.E.E. trans. Autom. control., 27, 17.
11-1
Y. Toure NOMENCLATURE
r
reaction rate
Pr
pressure temperature of the j-th element
F.
mass flow-rate
J v. J
F
axial velocity (v
If we assign a confidence degree a for the parameters set, we point out a confidence volume in the parameters space which is an ellipsoid whose main axis are inversely proportional to the singular values of the sensitivity matrix. 68 T (SbT Sb) 68 = -J F(&) p/(N-p) = Q(a-)
g
-)
S. p .
j
J J
h.
k. s. JP JP
h.
inverse of time constant
k. JP
heat transfert coefficient between the j-th heating element and the mixture
s.
exchange area per length unit between the j-th heating element and the mixture
J
JP
J
S.C . p
J PJ g
S.
cross section for the fluid in thej-thelement
C .
heat capacity
p.
density
J
PJ J
(1/.
J ~ the criterion value at the identification end.
T.
J
rl
T with Sb = U . L . v where L is the singular values diagonal matrix if T T 68 (V LT L V ) 68 = Q(&) T Letting y = V 68 ,
p ~
yi/(Q«(i)/ Oi)
i=l The confidence volume is the set of 8 contained in this ellipsoid centererl in 6*. Results The identification is processed on different models using a set of 70 measurements.
e : heating wire ("+" or 11_" with respect to
the circulating direction of the reaction mixture) d
(double) shell
p
reaction mixture
g
thermal fluid
Criterion : .398 Singular values : 6000, 19000, 500, 2.1, .2, 1 . 9 10-' Ellipsoid main axis
molecular weight of W molecular weight of B APPENDIX
Initial model
KINETIC MODEL ANALYSIS AND A POSTERIORI IDENTIFIABILITY
The aim is to calculate a confidence interval for the p estimated parameters with regard to the set of measurementS. Notations S is the N-vector of the measurements realized on the system. M* = M( 8*) is the unknown N- vector which would be the output of the model tuned with the optimal values 8* of the parameters. M= M(e5 is the output N-vector of the model tuned with the estimated values of the parameters.
e
0
7. 10-'
4 10-'
- .01
-.01
- .01
0
0
0
0
0
0
0
0
-2
50
-.02
.6
0
0
.03
.4
0
0
0
0
0
0
3
50
7 10- '
-.01
- .03
0
0
0
The singular values range is enormous (2. 10 7 ), which points out a very important correlation bet ween the parameters a2 and b2 which are badly estimated. The parameter precision is : 30 %,
6b!ibl
18 %,
2000%,
6b2/b2
2000%,
.8 %
72%
Simplified models Assumptions - The model is valid. Without any noise, the vectors H * and S would be equal. - The noise vector is a normal random vector. Each componant is independant, centered, presenting a variance 0 2
•
- Though the model is nonlinear, around the optimum values of the parameters and for small deviations, it can be linearized with respect to the parameters. The sensitivity matrix is denoted Sb : M(k) = M*(k) + aM(k)/a 8i.( 8- 8t) or M = M* + Sb68 Fisher test
After analysis of this matrix, we can look after 2 reductions. The first one consists in grouping the 2 E parameters, which have close values, into one. The second one, moreover, sets to zero the a2 cof ficient which has the worst estimation. The criteria roughly remain the same (.403), the singular values range decreases, the precisions get better and the identification algorithm convergence becares faster. The parameters precision for the last model resumes the improveness : 30 Z,
6bl / bl = 18 Z,
6E/ E = 1 Z,
43 Z.
This linear development generates the TIp affine plane in the RN measurement space shown on Fig. 5. The orthogonal vectors el = M - M* (p components) and e2 = S- M(N- p components) allow to define a Fisher variable with p and N-p freedom degrees F(a) = < el, el>/p < e2, e2>/(N-p)
68T(S~
Sb) M J
N- p
IIM - M* II
P
11 S - M 11
Other simplifications have been done (for example setting b2) but the criterion grew while the precision remained the same.
115
DVllamic Mode lling and Simulation of an Industrial Amidatioll Reactor
R. 1-1
Ri
T
~ e
I===========T=:e:===.'~"
Ri +l
i npu t convey i ng fl u id
~
E
ouput con vey in g
T
d
r OI---------------:--~~ z/L
fluid
Fi g . 1. Section desc ription
1.00 1
1.00 T
. BO •
.80 •
.60 •
. 60 •
.40 •
. 40 •
. 20 •
. 20 TIME (1IVl) ._- -+- -
. 0
50
- - .. . - - - ----t
100 150 200 250 300 350 INITIAL CONCENTRATION ( ~ ) = 30 S T = 210 C
•
0
400
50
TIME (IrI) - - --t- _ "_ _
..... _ .
-
t -
---
---t ___ _ _ t ______ __ .... _ _ - - - f
100 150 200 250 300 350 INITIAL CONCENTRATION ( W ) · ]8 ST · 210 C
400
1.00
,
1.00,
. 80 ;
.80
+
.60 •
.60
.40 ,
.40
.20
.20 •
. 0
I-_-.-J
I - _........L - _ _ + - - _._-+-- ___ t - - - - -+
62 .5
TIME (IrI)
_ _ t-- _ _ _ -+ ____ - - - - i
125 187.5 250 312.5 375 437 .5 INITIAL CONCENTRATION ( W) = 3 0 ST' 200 C
500
TIME (IrI) .0
50
100 15O 200 250 300 J50 INITIAL CONCENTRATION (W ) • 7 ST· 230 C
Fig. 2. Kinetic model responses and experimental data (the vertical axis is the normalized concentration Y/Y max )'
400
lIt)
Y. Toure et ai,
heating set points
composition and temperature mixture input (chemical composition initial telnpcrature)
Reactor
energy comsumption
Fig. 3. Static simulation tool scheme.
o o o
o
I"-
I"-
o o
(Xl
(T)
Ol
Ol
o
o
o
Lf)
Lf)
..,..
I"Ol
I"Ol
I"Ol
I
.
WO
WO
f-
f-
N
N
Ol
Ol
(0
o
o o
(Xl
Ol
+
o o o
I"(Xl
00 0....0 f-
f-
o(Xl
(0
o
o
0
0
Lf)
(Xl ..,..
..,..
Ol
Ol
Ol
I"-
0
0
0
0
Lf)
0.00
0.10
0.20
0.30
0.40
0.50
Z/L Fig. 4a 0 0
0 0
~
~
0
0
0N
0
0.60
0.70
0.80
0.90
1.0
Dynamic Modelling and Simulation of an Industrial Amidation Reactor
Fig. 5. Id entifiabi lity : definition of a Fisher variable by means of the e rr o r vec tor projection.
117
DY~A~IICS
OF
CO~T1NUOUS
REACTORS
DYNAMIC MODELLING AND SIMULATION OF AN INDUSTRIAL AMIDATION REACTOR Y. Toure*, J. Biston*, C. Darmet*, G. Gilles* and J-F. Thierry** *Laboratoire d'Automatique et de Genie des Proctides, Unit'ersite Lyon I , Villeurbanlle, France **RhOne-Poulenc Recherches, Saint-FOIlS, France
Abstract. In order to get a better control of an amidation industrial process, a dynamic modelling study of the plug-flow reactor has been carried out. An infinitesimal phenomenological balance leads to a set of coupled partial differential equations. The chemical s t a ti c part of the model i s got from an analogy with respect to a batch reactor kinetics. Identifiability results of this non linear model are presented. From the non linear distributed parameter dynamic knowledge model, a static model is drawn out and leads to the writing of a static simulation package in order to find the static profiles and to choose the operating points. Keywords. Textile industry modelling ; identification
chemical industry simulation.
INTRODUCTION
distributed parameter systems
length in which circulates a thermal fluid. These sectlons a re thermally independent from each other. Each one is itself a combination of parallel flow and counterflow heat exchangers, displayed in the reactional mixture and in the double shell of the reactor. A typical section is represented on Fig. 1.
Getting a better quality and a better regu larity in the production of synthetic threads implies a very good control of the amidation plants on the upper part of the manufacturing process. These plants being constituted of a set of rather complex apparatus, a very good knowledge of the static and dynamic behaviours of each of these unitary processes is necessary for the design of a mov"e performant control system. Among these apparatus, the amidation reactor takes a prominant part because it contributes the most to the ex tent of the reaction. This evaporator reactor is a tubular heat exchanger in which circu l ates a thermal fluid and which is the center of coupled thermal and chemical phenomena. An infinitesimal phenomenological balance allows to get a dynamic model under the form of a set of coupl ed partial differential equations. One can derive a static model giving the characteristic variables profiles along the reactor.
A constant flow-rate of a homogeneous reaction mixture composed of water, diamine and diacid is heated in the reactor. The water of the mixture is evaporated while the polyamidation reaction progresses. The kinetic model of this chemical reaction depends upon the water concentration of the medium. KINETIC MODELLING Amidation reaction The polyamide 66 (Nylon 66) can be prepared by linear condensati on of a diamine, the hexamethylene diamine (HMDA) and a diacid, the adipic acid. The global reaction can be represented as follows
After a description of the industrial amidation reactor, the reactional part of the static model is analysed by considering the analogy existing with the reaction kinetics inside a batch reactor. This method presents th e advantage of allowing easier realizations of the experiments on a laboratory batch reactor (Rhone-Poulenc Company), and so of allowing "a validation of the static chemical model without any experiments on the plant. Identifiability results are presented.
n H2N - (CH2)G - NH
+
n HOOC - (CH6)'- COOH
HUm - (CH2)6 - NHCO - (CH2h - COJ
~
OH + (2n-1)HlO
n A simplified approach is to consider only functional groups in the reaction medium : end groups - COOH and - NH" intramolecular groups - CONH and molecules of water, leading to the following representation : - COOH + - NH,
The following part is devoted to the elaboration of the global dynamic model ; then, the static model of the industrial reactor is used in order to crea te a simulation tool which will help a process operator to easily tune the operating points. Then it is shown how to use the dynamic model in order to build a control model which is necessary for the design of a new process control system.
-+ +-
- CONH -
l
+ H2 O
or, in a symbolic form A (acid)
+
B
(base)
1 -+ +-
Y
+
W
l
(amide)
(water)
An equilibrium can be reached, when the polycondensat ion reaction (1) is balanced by the hydrolysis reaction (2). Kinetics
PROCESS DESCRIPTION
From the two following assumptions - all reacting species A, Band Y have the same reactlvlty (Flo r y's equireactivity" assumption) which does not depend on the degree of polymeri-
It is an isobaric continuous reactor. On a thermal point of view, it" can be considered as aheat exchanger composed of several heating sections of equal
III
112
Y. Tour<' fI al.
zation of the molecule - the reaction obeys a first order kinetics with respect to each reacting species (overall second order kinetics), a kinetic equation can then be derived : kl A B - k z Y W Allowance must be made for the highly non ideal character of the medium, especially due to the dependence of activity coefficients on the water content (Steppan, 1987). The general expression of rate "constants" kl and kz is given by : r
~
fla , b , W) exp( - Ei/RT) (i ~ 1,2) i i i The six parameters (ai, bi, Ei) of the model have been estimated from experimental data obtained on a batch laboratory scale reactor. From these expek
~
riments, temperature, pressure and concentrations
were chosen in a range consistent with industrial corresponding conditions. The Marquardt method (Marquardt D.W., 1963 ; Lawson C. L., 1977) was used for parameter identification from the dynamic model of the batch reactor. The sensitivity analysis developed in the appendix allowed to reduce the number of the model parameters inside the considered working section. Fig. 2 shows the comparison between this model and the measurements (+) on the batch reactor.
In this part, we consider all the physical aspects which surround the chemical reaction. One of the results of the chemical kinetics modelling is the quasi non thermality of the reaction . This will allow some further simplifications. Phenomenological balance The main assumption is considering the reactor as an isobaric plug flow reactor. The infinitesimal mass and energy balance are then set inside each section Ri(i ~ 1 to n), considering only one spatial variable along the flow axis of the reaction mixtu r e. Any diffusion along this axis is neglected (J. Villermaux, 1982 ; L . Borel, 1984 ; H.S . Fogler, 1981) .
One of the main phenomena is vaporisation . By means of the energy contribution, the water is mainly evaporated inside the medium, which leads to taking into account the water vaporisation enthalpy which depends upon temperature and pressure : Lv (T p ' Pr) (J. Ficini, 1977 ; M. Bertin, 1982 ; L. Borel, 1984) . So, the mixture flow-rate conti nuously varies along the section Ri and along the whole reactor. This phenomenon is modelled by the medium water flow-rate balance : dF dz
dx
1 [r ___
e
dZ
xe - 1
(1)
x (z,t) is the mixture water mass fraction and Fe(z,t) denotes the mass flow-rate of the reaction m~xture . By means of the reaction factor r , we then get a coupling of the chemical reaction with the thermal aspect because the vaporisation energy can locally be written as dF E
v
~ - L (T v
p'
P) r
........E. dz
dx
aT ........E.
e
dZ
3z
(3)
The equations (1), (2) and (3) describe phenomena that are supposed to be instantaneous (mass conser' vation, equi libr ium) . The main dynamic aspects especially occur on the level of thermal exchanges and of the mass balance for the main reaction. We then get a set of coupled partial differential equations. The heat exchanger We showed (Fig. 1) that the heating system of each section Ri can be decomposed into simple heat ex changers involving only heat exchanges with the rea< tion medium but not between themselves. The heat exchange with the medium is mainly processed by convection. For the reaction mixture, we then ge t the following energy balance equation : aT k sdp(T - T ) + k s (T+ - T Sp Cp Ppp d dp P ep ep e F
af
GLOBAL DYNAMIC MODEL
........E.~
This state function, due to Ogata (Ogata, 1960) h, been reevaluated among Rhone-Poulenc research grot in order to take into account the reaction extent. It relies the medium water concentration with the temperature and the pressure, which allows to writ t the variation of xe(z, t) under the following ex pression :
(2)
By mass conservation, the flow-rate variation of the liquid reaction mixture is equal, in terms of absolute value, to the vapor flow-rate variation . This is due to the fact that the mixture flow-rate is constant at the entrance of the reactor. We consider that this flow-rate variation is processed in order to satisfy the liquid-vapor equilibrium law of the medium water . Thus we introduce the state function of the water-polymer mixture.
+ k
3T s (T - T ) - F C __ P _ E ep ep e p p pp 3z v
(4) The coupling with the chemical reaction is here explicit by means of the Ev term «1) and (2» . For the heating elements, we get the classical equations of the parallel flow or counterf l ow heat exchangers with constant flow rates (Sacadura, 1980 J .C. Friedly, 19 72). Thus, for the R'1 section, we get
-
3T e
at aT
3T v
e
+
e
e
az- -
h (T e e
- Tp ) (5)
3T+
aze -
h (T+ - Tp) e e
at
- v
3Td
d YV d -3z - h d (T d - Tp )
at
e
aT
In the last equation, the Y factor (0 < Y <1) is a distribution coefficient of the double shell energ y between its lower part filled with the liqu id reaction mixture, its upper part (vapor) and the heat losses. Mass balance It is a balance according to the Fluid Mechanics concepts (Vi ll ermaux, 1982) . The vaporisation pheno menon occurs in the balance and we consider that there is a small hexamethylene diamine loss (B) (a fraction of the total vapor flux, with a « 1) 3A
at
3B
3B
at
- r - PpSp Fp
3W
3W
at
a; -
r - PpSp Fp
B PpSp W
az - PpSp
3F
2
3z 3F
........E. 3z
Dvnamic Modelling and Simulation of an Industrial Amidation Reactor (6)
r -
THE STATIC SINULATION TOOL The static model, derived from equations (1) to (6) is written by considering the steady state behaviour of the process (the partial derivative with respect to time is then equal to zero in the equations (4) to (6)). The aim of this simulation tool is double rify the model with respect to the steady behaviour, then and especially to predict tem output variables in order to tune the points. This is schematized on Fig. 3.
: to vestate the sysoperating
A previous work of heat exchange parameter estimation on each section Ri is carried out, taking into account their dependency upon the operating point, and also of estimation of physical parameters which depend upon the thermal fluid and which have been set as a constant in each section. Solving the set of non linear differential equations obviously implies the knowledge of the boundary conditions. These conditions depend upon each Ri section. On a typical structure shown on Fig. 1, the boundary conditions are : T- (1) e
T ei
T~(O)
T+ (0) e
T+ (1) e
Td (1)
T (i) p (0)
T (i-1) p (1)
(input of the heating fluid)
(7)
113
into account the residence of the mixture in the whole reactor, and to tune the operating points. Fig. 4 shows the simulation result for the section of Fig. 1. The measurement points are located at the section boundaries. A sample at the output of the section allowed the chemical model validation (Fig. 4-b). The temperatures (Fig. 4-a) are normalized with respect to the nominal input temperature of the thermal fluid Tei : Tj becomes Tj/Tei. The mass concentrations A , B , Y , Ware expressed in terms of equivalent by kilogram. CONCLUSION AND PROSPECTS A dynamic model of an industrial amidation reactor has been presented as a set of coupled partial differential equations. From this model, a static model has been derived, allowing to build a simulation tool which elaborates the profiles along the reactor and which can give an off-line help and eventually an on-line help to the choice of the reactor operating points. The dynamic model, presently analyzed, can be reduced by means of a transformation of some partial differential equations into ordinary differential equations, in order to obtain a more realistic con-
trol model. By means of intermediate measurements along the reactor, specific distributed parameter control methods (Pohjolainen, 1982) can be used. From an approximated, entirely differential, model we can also apply more classical methods of optimal control under constraints with observer. A comparative analysis must point out the efficiency of the new control systems with regard to the existing system. REFERENCES
This implies that the simulation initial conditions - for example T~(O), T~(O), Td(O) must in fact be boundary condltions because of the cascade structure of the heating elements of the Ri section in the real case. The simulation software must iteratively solve the so defined problem with limits, rapidly converge and be robust with respect to initialization. The numerical solving algorithm which has been chosen is the Newton-Gauss one (J.E. Dennis, 1983).
Steppan, D.D., M.F. Doherty and M.F. Malone (1987). A kinetic and equilibrium model for Nylon 66 Polymerisation. Journal of appl. polymer science vol. 33, 2333-2344. Marquardt, D.W. (1963). An algorithm for leastsquares estimation of nonlinear parameters. SIAM J., 11, 431-441. Lawson, C.L., R.J. Manson (1977). Solving least squares problems. Prentice-Hall, Inc. Villermaux, J. (1982). Conception and functioning of reactors, Lavoisier, Paris.
This leads to solving a system tions : find the x vector such is written under the form : f(x) = RT(x) R(x) that we have example, i~ the-case of system xT
of non linear equaas f(~) = 0., f(~) to minimize. For (7) :
IT - ( 1), T- (0), T+ (0), Td ( 1), T . (0 )] Lee e pl Tpl.(0) - TP'. 1(1) T- (1) - T . e el
R(x) [
T-(O) e
J
T+(O) e
«1) - T (l) d The initialization problem concerns T .(0), T:(O) and Td(O). The minimization algorithmP~s : ~+ = x_ -
D(~+)
T
J(~+)J
-1
J(~+)
T
R(~+)
for the Ri section (i = 1 to n). ~+ and J(~+) (Jacobian of f(x) for x = x+) are numerically evaluated by means-of solving-differential equations (J.P. Nougier, 1987 ; K. Arbenz, 1980 ; J. Friedly, 1972). For the realization, "A.C.S.L." simulation package has been used on I.B.M. - P.C./A.T. computer. This method gives good results, converging after a maximum of 4 iterations. This allows to propose a tool for help to decision if one takes
Borel, L. (1984). Thermodynamic and energetics. Presses Poly techniques Romandes. Suisse Fogler, H.S. (1981). Chemical Reactors. ACS Symposium series 168. Ficini, J., N. Lumbroso, J.C. Depsay(1977). Thermodynamic and chemical equations. Herman. Paris Bertin, M. J.P. Faroux, J. Renault (1982). Thermodynamics. Bordas. Paris Sacadura,J.F. (1980). Initiation to heat transfert. Tech & Doc. Paris. Friedly, J.C. (1972). Dynamic Behaviors of processes. Prentice-Hall, New Jersey. Ogata, N. (1960). Studies on polycondensation reactions of Nylon salt. Makromolekulare Chem. 42, 52. Dennis, J.E., R.B. Schnabel (1983). Numerical methods for unconstrained optimization and nonlinear equations. Prentice Hall, New Jersey. Nougier, J. (1987). Numerical computation methods. Masson. Paris. Arbenz, K., A. Wohlhauser. (1980). Numerical Analysis. Presses Poly techniques Romandes. Suisse. Pohjolainen; S.A. (1982). Robust multivariable PIController for infinite dimentional systems. I.E.E.E. trans. Autom. control., 27, 17.
11-1
Y. Toure NOMENCLATURE
r
reaction rate
Pr
pressure temperature of the j-th element
F.
mass flow-rate
J v. J
F
axial velocity (v
If we assign a confidence degree a for the parameters set, we point out a confidence volume in the parameters space which is an ellipsoid whose main axis are inversely proportional to the singular values of the sensitivity matrix. 68 T (SbT Sb) 68 = -J F(&) p/(N-p) = Q(a-)
g
-)
S. p .
j
J J
h.
k. s. JP JP
h.
inverse of time constant
k. JP
heat transfert coefficient between the j-th heating element and the mixture
s.
exchange area per length unit between the j-th heating element and the mixture
J
JP
J
S.C . p
J PJ g
S.
cross section for the fluid in thej-thelement
C .
heat capacity
p.
density
J
PJ J
(1/.
J ~ the criterion value at the identification end.
T.
J
rl
T with Sb = U . L . v where L is the singular values diagonal matrix if T T 68 (V LT L V ) 68 = Q(&) T Letting y = V 68 ,
p ~
yi/(Q«(i)/ Oi)
i=l The confidence volume is the set of 8 contained in this ellipsoid centererl in 6*. Results The identification is processed on different models using a set of 70 measurements.
e : heating wire ("+" or 11_" with respect to
the circulating direction of the reaction mixture) d
(double) shell
p
reaction mixture
g
thermal fluid
Criterion : .398 Singular values : 6000, 19000, 500, 2.1, .2, 1 . 9 10-' Ellipsoid main axis
molecular weight of W molecular weight of B APPENDIX
Initial model
KINETIC MODEL ANALYSIS AND A POSTERIORI IDENTIFIABILITY
The aim is to calculate a confidence interval for the p estimated parameters with regard to the set of measurementS. Notations S is the N-vector of the measurements realized on the system. M* = M( 8*) is the unknown N- vector which would be the output of the model tuned with the optimal values 8* of the parameters. M= M(e5 is the output N-vector of the model tuned with the estimated values of the parameters.
e
0
7. 10-'
4 10-'
- .01
-.01
- .01
0
0
0
0
0
0
0
0
-2
50
-.02
.6
0
0
.03
.4
0
0
0
0
0
0
3
50
7 10- '
-.01
- .03
0
0
0
The singular values range is enormous (2. 10 7 ), which points out a very important correlation bet ween the parameters a2 and b2 which are badly estimated. The parameter precision is : 30 %,
6b!ibl
18 %,
2000%,
6b2/b2
2000%,
.8 %
72%
Simplified models Assumptions - The model is valid. Without any noise, the vectors H * and S would be equal. - The noise vector is a normal random vector. Each componant is independant, centered, presenting a variance 0 2
•
- Though the model is nonlinear, around the optimum values of the parameters and for small deviations, it can be linearized with respect to the parameters. The sensitivity matrix is denoted Sb : M(k) = M*(k) + aM(k)/a 8i.( 8- 8t) or M = M* + Sb68 Fisher test
After analysis of this matrix, we can look after 2 reductions. The first one consists in grouping the 2 E parameters, which have close values, into one. The second one, moreover, sets to zero the a2 cof ficient which has the worst estimation. The criteria roughly remain the same (.403), the singular values range decreases, the precisions get better and the identification algorithm convergence becares faster. The parameters precision for the last model resumes the improveness : 30 Z,
6bl / bl = 18 Z,
6E/ E = 1 Z,
43 Z.
This linear development generates the TIp affine plane in the RN measurement space shown on Fig. 5. The orthogonal vectors el = M - M* (p components) and e2 = S- M(N- p components) allow to define a Fisher variable with p and N-p freedom degrees F(a) = < el, el>/p < e2, e2>/(N-p)
68T(S~
Sb) M J
N- p
IIM - M* II
P
11 S - M 11
Other simplifications have been done (for example setting b2) but the criterion grew while the precision remained the same.
115
DVllamic Mode lling and Simulation of an Industrial Amidatioll Reactor
R. 1-1
Ri
T
~ e
I===========T=:e:===.'~"
Ri +l
i npu t convey i ng fl u id
~
E
ouput con vey in g
T
d
r OI---------------:--~~ z/L
fluid
Fi g . 1. Section desc ription
1.00 1
1.00 T
. BO •
.80 •
.60 •
. 60 •
.40 •
. 40 •
. 20 •
. 20 TIME (1IVl) ._- -+- -
. 0
50
- - .. . - - - ----t
100 150 200 250 300 350 INITIAL CONCENTRATION ( ~ ) = 30 S T = 210 C
•
0
400
50
TIME (IrI) - - --t- _ "_ _
..... _ .
-
t -
---
---t ___ _ _ t ______ __ .... _ _ - - - f
100 150 200 250 300 350 INITIAL CONCENTRATION ( W ) · ]8 ST · 210 C
400
1.00
,
1.00,
. 80 ;
.80
+
.60 •
.60
.40 ,
.40
.20
.20 •
. 0
I-_-.-J
I - _........L - _ _ + - - _._-+-- ___ t - - - - -+
62 .5
TIME (IrI)
_ _ t-- _ _ _ -+ ____ - - - - i
125 187.5 250 312.5 375 437 .5 INITIAL CONCENTRATION ( W) = 3 0 ST' 200 C
500
TIME (IrI) .0
50
100 15O 200 250 300 J50 INITIAL CONCENTRATION (W ) • 7 ST· 230 C
Fig. 2. Kinetic model responses and experimental data (the vertical axis is the normalized concentration Y/Y max )'
400
lIt)
Y. Toure et ai,
heating set points
composition and temperature mixture input (chemical composition initial telnpcrature)
Reactor
energy comsumption
Fig. 3. Static simulation tool scheme.
o o o
o
I"-
I"-
o o
(Xl
(T)
Ol
Ol
o
o
o
Lf)
Lf)
..,..
I"Ol
I"Ol
I"Ol
I
.
WO
WO
f-
f-
N
N
Ol
Ol
(0
o
o o
(Xl
Ol
+
o o o
I"(Xl
00 0....0 f-
f-
o(Xl
(0
o
o
0
0
Lf)
(Xl ..,..
..,..
Ol
Ol
Ol
I"-
0
0
0
0
Lf)
0.00
0.10
0.20
0.30
0.40
0.50
Z/L Fig. 4a 0 0
0 0
~
~
0
0
0N
0
0.60
0.70
0.80
0.90
1.0
Dynamic Modelling and Simulation of an Industrial Amidation Reactor
Fig. 5. Id entifiabi lity : definition of a Fisher variable by means of the e rr o r vec tor projection.
117