IFAC
Copyright c:> IFAC Dynamics and Control of Process Systems, Jejudo Island, Korea, 200 I
D
0
[>
Publications www.elsevier.coml1ocatelifac
NONLINEAR NON-STIFF MODELS OF REACTIVE DISTILLATION COLUMNS WITH TWO-TIME-SCALE DYNAMICS Nishith Vora *,1 Prodromos Daoutidis *,2,3
* Department of Chemical Engineering and Materials Science
421 Washington Ave. SE, University of Minnesota Minneapolis, MN 55455
Abstract: In this paper, we address the dynamic modeling of reactive distillation processes, where the presence of fast reaction or momentum transport induces two time scale dynamics. We present a systematic modeling framework that allows identifying the source of time S<'.a.le multiplicity and deriving nonlinear non-stiff models of the slow dynamics. Copyright © 20011FAC Keywords: NOnlinE'M Model Reduction, Reactive Distillation, Singular Perturbation, High Index DAEs
1. INTRODt'CTION
Reactive distillation, the process of carrying chemical reaction and multi-stage distillation simultaneously, offers several advantages over conventional reactor-separator configurations (Doherty and Buzad (1992». Motivated by these, the last decade has witnessed a flurry of research on design (e.g. Ung and Doherty (1995», steady state multiplicities (e.g. Ciric and Miao (1994); Giittinger and Morari (1997); Kumar and Daoutidis (l999b», dynamic modeling and simulation (e.g. Alejski and Duprat (1996); Kumar and Daoutidis (1999b); Ruiz et al. (1995» and control (e.g. Balasubramhanya and Doyle III (1998); Kumar and Daoutidis (l999b); Sneesbyet al. (1997); Vora and Daoutidis (2001». In the present work, we address the dynamic modeling of reactive distillation processes. The simultaneous occurrence of various phenomena such as reaction, flow and mass transfer in a reactive distillation column can potentially lead to stiff 1 Current affiliation with General Electric, Corporate Research & Development, Schenectady, NY 12301 2 To whom all correspondence be addressed 3 The authors acknowledge the partial financial support from the Graduate School, University of Minnesota
(ill-conditioned) dynamic models. It is well documented that such stiff dynamic models present difficulties in dynamic simulation, and also result in ill-conditioned model-based controllers (Kumar and Daoutidis (1999a); Vora and Daoutidis (2001». Motivated by these considerations, we focus on the development of a dynamic modeling framework for reactive distillation columns that allows the systematic identification of various sources of stiffness, and the subsequent development of nonlinear non-stiff models suitable for simulation and control.
2. DYNAMIC MODELING A detailed tray-by-tray dynamic model for a multicomponent system, under the assumptions of well-mixed liquid and vapor phases in thermodynamic phase equilibrium and negligible vapor holdup, comprises of differential equations accounting for material and energy balances on each tray, and algebraic equations arising from thermodynamic relations and tray hydraulics correlations (see e.g. Alejski and Duprat (1996); Ruiz et al. (1995»; specifically, in the case of a single reaction the model in a stage i has the form:
119
M/i
= Fi +Vi+1-Vi+Li - 1-L i+(L Vj)ri VIi
(1)
j
Xj,i
= (Fi(xfj,i -
Xi,i)
+ Vi+l(Yj,i+1
- Xi,i)
+Li-1(xi,i-l - Xi,i) - Vi(Yi,i - Xi,i) +(vi - L VjXl,i)ri Vli)/M1i i 1'i = (Fi(Tli - T;) + Vi+1 (Ti+1 - Ti )
+ L i - 1 (Ti - 1 - T i ) t!.Hv + ( Vi+l - Vi)-c
(2)
(8)
t!.Hro - riVli--' )/M/i(3)
p
er,
0= Yj,i
(5)
t!.Pi
(6)
Vi i -
Pi R:I'i
JT
We initially examine the fundamental structure of the algebraic equations that determine the reaction rate (see Eq.5), the vapor flowrate (see Eq.6) and the liquid flowrate (see Eq.7). Specifically, consider the reaction rate expression of Eq.5; it can be expressed as:
(4)
(7)
where Mli is the molar liquid holdup, Xi,i is the mole fraction of component j in the liquid phase, and Ti denotes the temperature in K. Fi denotes the molar feed rate on feed stages, x fj,i is the mole fraction of component j in the feed and Tli denotes the feed temperature in K. t!.Hri is the heat of reaction in J Imol and Vj is the stoichiometric coefficient of species j. The heat capacity cp (in J /mol-K) and latent heat of vaporization t!.Hv (in J /mol) are assumed to be constant and equal for all components. Yj,i denotes the mole fraction of component j in the vapor phase (see Eq.4). Pi and PJ,i denote the total pressure and the saturation pressure of component j, respectively in N/m 2 , 1i,i denotes the activity coefficient accounting for the nonideality in the liquid phase and
where kin".., is the reaction rate constant at the nominal steady state conditions, and Cli = k;/klnoTft • In the same vein, consider the pressuredrop correlation of Eq.6; defining Vi = t!.P;/{, it can be written as:
J
(9)
where P~om is the vapor phase molar density at the nominal steady state conditions, and
(10)
where Pnom is the average molar density at the nominal steady state conditions, and ~i = p;/ Pnom. Substituting the expressions of Eq.8,9,10 to Eqs.l-7, and expressing these in a matrix form, we obtain: Xi = h(x)
+ P~om bvi(x) kvi(x) + kl
+Pnom b/i(x) kli(X)
n •m
bri(X) kri(x)
0= Yi,i rPj,j Pi - 1i,i Xi,i PJ,i ri = klnom Cli Ti Vi = P~om
(11)
where kri(X) = ri, kVi(X) = [ViIVi+11 T and k1i(X) = [Li-1ILif, and the corresponding bri(X), bvi (x) and b/i (x) are appropriately defined. In the model of Eq.ll, we have isolated the parameters (k lnom , P~om and Pnom), which can be relatively large in case of fast reaction or fast liquid and vapor flowrate, and potentially cause stiffness (time scale multiplicity) in the model.
p,.
where
pi = .m,j
denotes the vapor phase molar
density, and ~ is a characteristic parameter of the correlation. The liquid flOv.Tate L i , is a function of the volume of liquid over the weir (~Vli = VliV~'i). The Francis weir formula is used to calculate the liquid flowrate (see Eq.7), where Pi denotes the average molar density in mol/m 3 , VIi and Vwi denote the volume of liquid and the weir volume on stage i in m3 , respectively, and CL denotes the characteristic hydraulic coefficient. The above equations are appropriately modified for the total condenser and the reboiler stage.
3. DERIVATION OF NON-STIFF MODELS In this section, we perform an analysis of the dynamic behavior of reactive distillation columns based on the dynamic model of Eq.ll, when fast reaction or fast liquid and vapor flowrate is present. For illustration purposes, we consider a reactive distillation column with 7 equilibrium stages, comprising of 5 trays (stages 2 to 6), one total condenser (stage 1) and one reboiler (stage 7). We consider a quaternary system of components A, B, C and D, where the following
120
liquid phase reversible reaction between A and B occurs to produce C and D:
effect of a chemical reaction near equilibrium, small pressure drop and small volume of liquid over the weir on the dynamic behavior of the column. In all the case studies, the representative plots demonstrating the transient response of a variable are for stages 2, 4 and 6, which correspond to the top, middle and bottom trays in the column, respectively. Also, the transient response is shown by the solid line, and the steady state is indicated by the dashed line. The reboiler holdup and the condenser holdup act as pure integrators and must be stabilized for such an analYSiS, and this was done in all the case studies by using a perfect controller that manipulated the bottom and distillate fiowrates. In addition, the pressure at the top is also controlled to ensure stability of the column in presence of open-loop disturbances, and moreover to minimize the effect of changes in phase equilibria on column dynamics. This was achieved by using a perfect controller that manipulated the condenser heat duty.
A+B <-> C+D The phase equilibrium is assumed to be ideal, and the relative volatilities of the components in decreasing order are D > A > B > C. The feed streams FA,6 and FB,2 of reactants A and B are fixed at stage 6 and 2 respectively. The products D and C are withdrawn from the top as the distillate stream, and from the bottom stream respectively. The reaction rate r i is given as:
where kJi and "'i represent the reaction rate constant and the reaction equilibrium constant, respectively, and are given by: k Ji = k'} exp( -
:~i)
_ k'} EJ -Eb "'i - kg exp( - R Ti )
3.1 Case I: Fast reaction
where k'}, kg are the pre-exponential factors, and E f' Eb denote the activation energies for the forward and backward reactions respectively. The column specifications and physical properties of the components are listed in Table 1.
In this case study, we analyze the dynamics of a reactive distillation column with a fast reaction. The activation energies for the forward and back, ward reaction are 75 K J /mol and 90 K J /mol, respectively. The characteristic coefficient in the pressure drop correlation is ~ = 70000 (see Eq.6). We initially consider the open-loop dynamic behavior of the above reactive distillation column, and consider a small perturbation (±1%, with a random determination of the sign) in the liquid phase compositions from their steady state values. Fig.1 shows the responses of the composition of A, the liquid holdup and the temperature. An initial
Table 1. Specifications of the column Qh Qc D R B
Pi A Cp
kO kb
6
FA,6 FB,2 TA TB PA PB PC PD
Vw We!! b.H"
Description reboiler heat duty (MW) condenser heat duty (MW) distillate flowrate (mol / s) reflux fiowrate (mol/s) bottom fiowrate (mol/s) Pressure at the top (bar) cross-sectional area (m 2 ) molar heat capacity (J/mole K) pre-exponential factor pre-exponential factor feed fiowrate (mol/s) feed fiowrate (mol/s) feed temperature of A (K) feed temperature of B (K) molar density of A (kmol/m 3 ) molar density of B (kmol/m 3 ) molar density of C (kmol/m 3 ) molar density of D (kmol/m 3 ) weir volume on each tray (m 3 ) effective weir width (m) latent heat (kJ /mol)
Value 1.40 1.32 32.59 87.93 67.41 4.23 2.0 80 1.0 x lO10 1.53 x 10 12 50 50 300 300 31 29 28 32 0.3 0.798 10
The approach pursued is to initially study the open-loop dynamic behavior in the presence of small perturbations in the initial conditions, with the objective to observe and isolate various sources of stiffness, and the associated two time scale behavior. Then we systematically derive the non-stiff dynamic model in the slow time scale using singular perturbation arguments. We consider three case studies to separately analyze the
lQ I
I
I" c
~Il
~0
.
~I
~I
i-
• • • • • •,I, Time, • •••
Time,
n~
I
8
..!!I
~II
e... 1
~I
0
i
I
8
I
..!!.
8
I
I
•••••
Time,
8
11 ~ I
0
~.
~Il
l
, • • • • • •, • • • • • Time, Time, I
11
I
I
. i
I
l
, • • • • • •, Time, ••••• Time, ,I
.
, Time, ••••• 8
I
.•
10
IS
I
I
8
~
I
I
•••••
Time,
8
Fig. 1. Composition of A, liquid holdup and temperature in stages 2, 4 and 6, in the slow time scale for ±1 % perturbation in liquid compositions in case I
121
fast transient in the order of magnitude of seconds and a slower approach to steady state in the order of magnitude of several hours can be observed for all state variables. Clearly, the system exhibits transients in two different time scales. Note however that the fast and slow dynamics cannot be attributed to specific state variables such as compositions, liquid holdups and temperatures. Now consider the transient responses of the reaction rate, vapor flowrate and liquid flowrate shown in Fig.2, for the same perturbation in the initial condition, in the fast time scale. As can be seen, the reaction rates approach their steady state values almost instantaneously, whereas the internal flowrates exhibit much slower dynamics. This indicates that the reaction rates are a possi-
Fig. 3. Reaction rate constant and condition of chemical reaction equilibrium in stages 2, 4 and 6, for ±1 % perturbation in liquid compositions in case I
·.10 ."fd"B ..... r .~ "
,f
/ ... .
~
~ •............_..•..........
'-l
• ........•..• _...•••..... _..
0= Yj ,i Pi - Xj,i PJ,i 1
;:.
I-
I
I
I
I
I
'I
I
I
"
I
I
r
I
I
I
"
I
I
ri
t
•
t[O}~}ra I
I
I
I
I
I
I
I
I
I
I
I
I
lIt
I
I
"
I
I
I
I
I
I I
'
r
I
I
I
I
I
"
I
I
I
I
,
I
I
Vi = P~om (7i iii Li = Pnom ~i Li
8
I~
Time,
I~
Time ,
I
,i
= 2, ... , 6
(13)
Note that the l/t:r term appears in all the differential equations of the above model; this is indicative of the fact that the state variables cannot be separated into fast and slow (see Fig.1 and the related discussion). Equivalently, the model of Eq.13 is in the non-standard singularly perturbed form (Kokotovic et al. (1986)). A non-stiff model of the slow dynamics of the system of Eq.13 can be obtained by setting t:r -+ O. In this limit, the algebraic constraints kri(X) = ri = o must be satisfied in the slow time scale. For the model of Eq.13, these constraints correspond to (XA ,i XB,i - «l-XA,i - XB,i - XD,i)XD,i!Ki)) = 0, which represents the condition of reaction equilibrium. Let Zri = limE~-+o ri denote the reaction
1
t[cl }r;~~·-;1!1~1 Time,
= -t:r O!i ri
I
.i
Fig. 2. Reaction rate, liquid flowrate and vapor flowrate in stages 2, 4 and 6, in the fast time scale for ±1% perturbation in liquid compositions in case I ble source of time scale multiplicity in the column. Fig.3 shows the profiles of kli and ri , i.e. the two factors involved in the calculation of r i , in the fast time scale (see Eq.8). One can observe (i) the fast approach of ri to steady state, and (H) the significant difference in the orders of magnitude of these terms, with kli being large (0(10- 1 )), and ri being small (0(10- 3 )) . The above documents that the reaction rates are indeed a source of time scale multiplicity in the column. This time scale multiplicity, and the associated stiffness caused by the large parameter kino", in the process model can be systematically addressed on the basis of the dynamic model in Eq.ll. Specifically, defining t:r = l/klno ", in Eq.ll, as a small parameter, the dynamic model of Eq.ll takes the form:
t:r
rate in the slow time scale. Then the dynamic model of the slow dynamics can be expressed as:
x = lex) + P~om bvi(x) kvi(x) +Pnom
bli(X) kli(X)
+ bri(X) Zri(X)
0= Yj,i Pi - Xj ,i PJ,i
O=ri Vi = P~om (7i iii L i = Pnom ~i Li
,i = 2, . .. , 6
(14)
l'iote that the index of the above DAE system is greater than one. A direct differentiation in time of the reaction equilibrium constraint in this case results in an equation that is solvable in Zri , hence the index of the DAE system describing the slow dynamics of the process is two.
122
3.2 Case 11: Fast vapor flow In the present case study, we consider a reactive distillation column with fast vapor flow. The characteristic coefficient in the pressure drop correlation is { = 700 (see Eq.6). The activation energies for the forward and backward reaction are 120 K J / mol and 135 K J / mol, respectively. Note that the larger values of the activation energy in comparison with case I imply a smaller reaction rate constant and a slower reaction in the column. We again study the open-loop dynamic behavior of the above reactive distillation column, and consider a small perturbation (±1%) in the liquid phase compositions from their steady state values. Similar to the earlier case study, the state variables such as compositions, liquid holdups and temperatures exhibit an initial fast transient followed by a slower approach to steady state. Fig.4 shows the plots of reaction rate, liquid flowrate and vapor flowrate profiles in stages 2, 4 and 6, respectively, in the fast time scale. As can be seen, the vapor flowrates approach
Fig. 5. Molar density in vapor phase and condition of constant pressure in stages 2, 4 and 6, for ±1% perturbation in liquid compositions in case IT ting tv -4 o. In this limit, the algebraic constraints kv.(x) = 0 are obtained that must be satisfied in the slow time scale. These constraints correspond to tlPi = Pi - Pi-1 = 0, i.e. the condition of constant pressure throughout the column. Notice that in this limit, the pressure-drop correlation Vi (see Eq.6) for calculating the vapor flowrate becomes indeterminate and the vapor flowrate Vi, is not explicitly specified by an algebraic equation but is implicitly fixed by the_ constant pressure
f8 f fg ·_············ J~E ~
...~ _..•...• _.................. I
I
I
I
I
a
~
J
I
I
,
I
J
I
i --
I
I
I
I
I
J
J
constraint. Let Zv ,·
J
~lC=l ~~ ~EJTime .. I
1
I
J
I
J
I
.
I
I
I
J
I
I
I!,
I
J
I
~~-................ Se. \ ~ I
I
,
I
Time,
J I
J
I
I'
I
I
Time,
J I
I
~
-;: J
+Pnom
1
I
Time,
I
I
fv
denote the vapor
bu(x) ku(x)
+ kl
nom
bri(X) kri(X)
0= Yj.i Pi - Xj.' PJ.i
r··" --
I'
40....!
x = f(x) + bvi(X) ZVi(X) I
;;bJTime. ~E1 -.-~::~.~.~..... ;;EJTime," 8
Vi
v
flowrate in the slow time scale. Then the model of the slow dynamics can be expressed as:
~r=:J}~ ~ r-I
= lim t
r. = kl
ai T. O=Pi -Pi - 1 nom
I
Li =
I
Pnom 6i
Li
,i = 2, ... ,6
(15)
Note that again the index of the above DAE system is greater than one. A d.4-ect differentiation in time of the constant pressure (equivalently, no pressure drop) constraint in this case results in an equation that is solvable in Zv; hence the index of the DAE system describing the slow dynamics of the process is two.
Fig. 4. Reaction rate, liquid flowrate and vapor flowrate in stages 2, 4 and 6, in the fast time scale for ±1% perturbation in liquid compositions in case IT their steady state values very quickly, whereas the liquid flowrates and reaction rates exhibit much slower dynamics, indicating that the vapor flowrates are the source of time scale multiplicity in the column. Fig.5 shows the profiles of pf and Vi, i.e. the two factors involved in the calculation of Vi, in the fast time scale. Note (i) the fast approach of "Ci to their steady-state values, and (ii) the significant difference in the orders of magnitude of these terms, with pf being large (0(10 3 )), and Vi being small (0(10 1 )). A non-stiff model of the slow dynamics of the system of Eq.ll in this case can be obtained by defining tv = 1/ P~om as the small parameter, and set-
3.3 Case Ill: Fast liquid flow
In this case study, we analyze the dynamics of a reactive distillation column with fast liquid flow. The characteristic coefficient in the pressure drop correlation is ~ = 70000. The activation energies for the forward and backward reaction are 120KJ/mol and 135KJ/mol, respectively. Note that this results in slower reaction and moderate pressure drop. Similar to the earlier case studies, again we can observe two time scale dynamics in
123
the column. Fig.6 shows the plots of reaction rate, liquid flowrate and vapor flowrate, respectively, for a small perturbation (± 1%) in the liquid phase compositions and the liquid holdups from their steady state values in the fast time scale. As can be seen, the liquid flowrates approach
Again note that the index of the above DAE system is two as a direct differentiation in time of the constant volumetric holdup constraipt in this case results in an equation that is solvable in Li • 4. CONCLUDING REMARKS
-~-U" -fC]
"'" ...-............................. ......... Q Q e e
t
.........
Q
In a reactive distillation column, the simultaneous occurrence and interaction of reactions, flow and mass transfer results in complex dynamics, which often occur in multiple time scales. This study presented a modeling framework that structurally isolates the potential sources of time scale multiplicity in a reactive distillation column, and allows derivation of non-stiff DAE models of index two under the quasi-steady state assumptions of reaction equilibrium, constant pressure throughout the column and constant volumetric holdup in a tray. State-space realizations of these models suitable for simulation and control can be derived using standard index reduction methods for DAEs (Kumar and Daoutidis (1999a)) .
~
S
111'1111
Time,
.. _............ _................
e
Time,
8
IIII
Time,
I
11
-r===J -rr==I-~
~12J}~jo I
I
Time,
I
•
I
I
I
Time,
•
•
I
8
,
Time,
•
•
I
~D · · · · · . · · · · · ·~lO~§ e e
(j
e
'0
.i
~
I
I
Time,
I
•
.......................... _......
e
.j
I
I
I
Time.
I
•
8
I
I
Time,
I
•
I
Fig. 6. Reaction rate, liquid flowrate and vapor flowrate in stages 2, 4 and 6, in the fast time scale for ±1 % perturbation in liquid compositions and holdups in case lIT
References K. Alejski and F. Duprat.
Dynamic simulation of the multicomponent reactive distillation. Chem. Eng. Sd., 51:4237-4252, 1996. L. S. Balasubramhanya and F. J. Doyle Ill. Nonlinear model-based control of a batch reactive distillation column. In Proc. of DYCOPS, 123-128, Corfu, Greece, 1998. A. R. Ciric and P. Miao. Steady state multiplicities in an ethylene glycol reactive distillation column. Ind. Eng. Chem., Res., 33:2738-2748, 1994. M. F . Doherty and G. Buzad. Reactive distillation by design. 'lrans. IChemE, 70:448-458, 1992. T. E. Giittinger and M. Morari. Predicting multiple steady states in distillation: Singularity analysis and reactive systems. Comput. Chem. Engng., 21 S:S995-S1OOO, 1997. P. V. Kokotovic, H. K. Khalil, and J. O'Reilly. Singular Perturbations in Control: Anal!lsis and Design. Academic Press, London, 1986. A. Kumar and P. Daoutidis. Control of Nonlinear Differential Algebraic Eqtiation S!lstems. Research Notes in Mathematics. Chapman & Hall/CRC, 1999a. A. Kumar and P . Daoutidis. Modeling, analysis and control of ethylene glycol reactive distilation column. AIChE J, 45(1):51~8, 1999b. C. A. Ruiz, M. S. Basualdo, and N. J. Scenna. Reactive distillation dynamic simulation. 'lrans. IChemE, 73: 363-378, 1995. M. G. Sneesby, M. O. Tade, R. Dutta, and T. N. Smith. ETBE synthesis via reactive distillation. 2. dynamic simulation and control aspects. Ind. Eng. Chem., Res., 36:1870-1881, 1997. S. Ung and M. F . Doherty. Synthesis of reactive distillation systems with multiple equilibrium chemical reactions. Ind. Eng. Chem., Res., 34:2555-2565, 1995. N. Vora and P. Daoutidis. Dynamics and control of an ethyl acetate reactive distillation column. Ind. Engng. Chem. Res., 40(3) :833- 849, 2001.
their steady state values very quickly, whereas the vapor flowrates and reaction rate exhibit slower dynamics, indicating that the liquid flowrates are the source of time scale multiplicity in the column. This time scale multiplicity can be addressed similar to case I and IT on the basis of Eq.ll, by defining €/ = 1/ Pnom and considering the limit as €/ -t 0, the algebraic constraints k/i(x) = o are obtained, that must be satisfied in the slow time scale. These constraints correspond to A Vii = Vii - Vwi = 0, which corresponds to the condition of constant volumetric holdup in a tray. Notice that in this limit, the Francis weir formula (see Eq.7) for calculating the liquid flowrate Li becomes indeterminate and the liquid flowrate, is not explicitly specified by an algebraic equation but is implicitly fixed by the constant volumetric . Let Z/, = I'zmE1-+O -Li denote hoIdup constramt. €l
the liquid flowrate in the slow time scale. Then the dynamic model of the slow dynamics can be expressed as: :i; =
!(x) + P~om bvi(x) kvi(x) +bli(X) Zli(X) + k fnom bri(x) Zri(X)
0= Yj,i Pi - Xj,i PJ.i ri = k fn om 0i fi
Vi = P~om (Ti ~ o= Vii - Vwi ,i = 2, ... ,6
(16)
124