Microsoft Word - TOC_R.doc HUNGARIAN JOURNAL OF INDUSTRIAL CHEMISTRY VESZPRÉM Vol. 36(1-2) pp. 149-153 (2008) DETERMINATION OF SAFETY OPERATING REGIMES BASED ON THE ANALYSIS OF CHARACTERISTIC EQUATION OF STATE-SPACE MODEL T. VARGA , G. HORVÁTH, J. ABONYI University of Pannonia, Department of Process Engineering, H-8200 Veszprém Egyetem street 10., HUNGARY E-mail: vargat@fmt.uni-pannon.hu Determining safe regions of operation is important at the design, control and optimization of process systems. Process alarms and safety systems should be designed based on such knowledge. In classical control systems these regions are represented as independent constraints on process variables. In case of a highly exothermic reaction takes place in any kind of reactor there is a high risk to develop reactor runaway if the produced heat can’t be removed. Reactor runaway phenomena is a serious problem in many chemical industrial processes, hence the most of reactions are exothermic. Reactor runaway means a sudden and considerable change in process variables such as reactor temperature. To detect and forecast the development of runaway a model based criteria is introduced in this manuscript. The criteria is based on indirect Ljapunov’s stability analysis of a process model. The aim of this work is to compare methods for calculation of eigenvalues. Eigenvalues of the Jacobian matrix were calculated with a numerical method and there were also computed analytically. Stability of two-phases fed-batch reactor with highly exothermic reactions are analyzed with the investigated stability analysis methods. Analytical solution is more sensitive and of course it is more accurate than the numerical one still the algorithm based on the analytical solution is a lot more faster. Keywords: operating regime, process model, simulation fed-batch reactor, stability analysis. Introduction Safe and at the same time optimal operation is always a key issue in the wide spectrum of chemical engineering problems. E.g. the design of optimal feeding trajectory of a fed-batch reactor where a highly exothermic reaction takes place. So to apply any kind of optimization methods some constraints on operating variables must be determined, i.e. the safe operating regimes must be characterized. In this paper the first step of this problem solution is investigated. It is introuced, how the safe operating regimes can be extracted from a process model. The most important phenomena from the safety of operation was investigated closely. It is called reactor runaway. It has contributed to industrial chemical accidents, most notably the 1984 explosion of a Union Carbide plant in Bhopal, India that produced methyl isocyanate. Thermal runaway is also a concern in hydrocracking, an oil refinery processes. Reactor runaway means a sudden and considerable change in process variables. It is a serious problem in many chemical industrial technologies, like oxidation processes and polymerization technologies [1-3]. E.g. in case of a highly exothermic reaction thermal runaway occurs when reaction rate increases due to an increase in temperature. It causes a further increase in temperature and it further increases the reaction rate. The temperature increasing will be stoped only if all reactants are depleted and reaction rate is getting zero. Detection of runaway has two main important aspects. On one hand runaway forecast has a safety aspect, since it is important for avoiding the damage of constructional material of reactor or in the worst case the reactor explosion. On the other hand it has a technology aspect, since the forecast of the runaway can be used for avoiding the development of hot spots in catalytic bed, which speeds up the ageing of catalyst [4; 5]. The avoid of runaway can be important in decreasing of the produced amount of byproducts, e.g. in synthesis of 2-octanone from 2-octanol [1; 6; 7]. A control system which is able to modify operating conditions of reactor in appropriate time decreases the costs and increases the safety of operation. The first step to develop such control system is the development of a reliable runaway criterion. Most of runaway criteria found in literature can be classified as data- or model-based criteria [8-10]. To apply data-based criterion it is necessary to have some measured data. It means the use of a data-based criterion has some restrictions on the forecasting. Other problem with data-based methods is found in measurement conditions, e.g. measurement noise can result in false forecast. Model-based criteria require parameter sensitivity and/or stability analysis, so for the application of these kinds of criteria it is necessary to have exact process model with correct model parameters. In the design process of a reactor in which a reversible reaction takes place the first step must be the calculation of equilibrium temperature. Reactor temperature can’t cross this line. After everything is 150 known about the equilibrium the calculation of optimal temperature profile or trajectory can be performed. It can be generated easily with differentiating of the rate of reaction with respect to reactor temperature. Hence the most of reactions in chemical industry are exothermic the generated heat must be removed to keep the operation in stable and controllable operating regime. In our earlier studies Ljapunov’s stability analysis method proved to be appropriate to detect reactor runaway [11; 12]. The applied stability analysis method based on the analysis of eigenvalues [13]. The aim of this article is the investigation of the relability of numerically calculated eigenvalues. To check the applicability of the choosen calculation method the results are compared with the analytically calculated eigenvalues. Next to relability the calculation time is the second important property of the proposed tool since it can be a part of an Operator Support System (OSS) [14-15]. After a short introduction about the analyzed chemical engineering problem and the possible detection technics the investigated case study is presented, than results obtained by applying the developed programs are compared and finally some conclusions are stated and some possible future steps are described. Stability analysis In advanced model based techniques reactor runaway is detected through the stability analysis of the process model, e.g. by Ljapunov’s indirect method [16]. The state space model of the process: ( ) ,xf d xd = τ where: x – state variable; τ – independent variable (e.g. time). Generation of the Jacobian-matrix of first principles model is the first step in application of Ljapunov’s indirect method to investigate stability: ( ) , x xf J 0x ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ = where J – Jacobian matrix; x0 – investigated work point. It is followed by the examination of eigenvalues of the Jacobian-matrix. ,0IJ =λ− where λ – eigenvalues; I – unit matrix. In case all of eigenvalues are negative than the model is stable but if one of eigenvalues is over zero than model is unstable at the investigated operating point of reactor. To calculate the eigenvalues of the Jacobian the determinant method can be applied. In case there are m correlations and n state variables the size of Jacobian matrix is nxm. To apply determinant method it is necessary that the number of correlations is equal with the number of state variables so the Jacobian matrix must be a square matrix. Eigenvalues can be calculated from the following characteristic polinomial: ( ) , JJ JJ IJdet mmm1m m1111 λ− λ− =λ− K MOM L Ljapunov’s stability analysis is suitable in detection of the development of runaway and the algorithm based on this stability analysis can separate the cases whether the runaway occurs in the reactor or not [11; 12]. However, it is quite difficult to implement this criterion in an industrial environment. Furthermore, it is applicable to detect when the runaway has been occurred, which is interesting from the analysis of historical process data point of view, but it is not usable for on-line process monitoring and control, where the operator is interested in the region of the safe operation and the last time instant when the runaway can be avoided by the control system. For the predictive analysis of the process not only a detailed (accurate) process model is needed, but also a process simulator that is able to estimate the trajectories of the process variables in case of normal and abnormal operations. This simulator should be also able to model the dynamical behavior control system, including its safety elements. This knowledge is extremely important since these elements of the control system define and sometimes ”widen” the region of safe operation. The most common method for calculating eigenvalues of a small Jacobian matrix is the calculation of roots of characteristic polynomial. It can be performed with applying numerical root-finding algorithms or in case the order of characteristic polynomial is not higher than four the roots can be calculated analytically. Next to root-finding the second common method is the power method or vector iteration method. In this paper the accuracy, reliabilty and the speed of the developed algorithms based on numerical and analytical root-finding are analyzed and compared in stability analysis of a complex system. Case study To illustrate the applicability of numerically calculated Ljapunov’s indirect stability analysis a well-stirred fed- batch reactor with jacket cooling was investigated, where 2-octanone is produced from 2-octanol in a highly exothermic oxidation reaction. 151 Reactor model The main product, 2-octanone is produced in a two- phases reaction system in a fed-batch reactor from 2- octanol. The 2-octanol is continuously fed into the organic phase, which does not exist before the start of the feeding. Hence, the organic phase is the dispersed and the aqueous nitric acid phase where all the reaction steps take place is the continuous phase. Two phases are connected by the mass and heat transfer phenomena. Based on Castellan et al.’s work [17] which shows that the oxidation processes in room temperature mostly take place according to ionic mechanism, Woezik and Westerterp [6] determined a reaction scheme describing the way of oxidation of 2-octanol to 2-octanone in nitric acid, which is proved to be feasible. In the first step of this reaction mechanism the nitrosonium ion forms. This reaction has a very long induction time, which can be shortened by adding a small amount of an initiator like NaNO2 to generate the necessary nitrous acid faster. The nitrosonium ion forms only in the aqueous phase. It is followed by the oxidation of 2-octanol forming 2- octanone and the further oxidation producing different carboxylic acids. This process can be considered as an advanced benchmark problem in the field of process engineering. In [1] the control of this process is analyzed. In [6] the safety issues related to the operation is discussed. A report how the simulator of this process can be used in process engineering education can be read in [7]. From the viewpoint of our paper it is the most relevant work where the Hazard and Operability Analysis (HAZOP) of this process is developed and also based on the application of process simulator. Mathematical model of the system was worked out by Woezik and Westerterp [6]. The introduced mechanistic model is based on a simplified form of the above introduced mechanism to predict the dynamic behavior of the reactor. The simplified reaction mechanism: XBP B2PBA ⇒+ +⇒+ where A is 2-octanol, B is nitrosonium ion, P is 2- octanone and X represents all the byproducts. In the next part just a brief introduction will be given about this model; a well detailed description can be found in [9-10]. The structure of the dynamic model can be found in Fig. 1. In the interconnection of organic and aqueous phases are considered only the mass transfer process and the temperature of the mixed phases are the same. The reactor insight and the jacket are coupled by only the heat transfer through the reactor wall. The main and the side reaction takes place in aqueous phase so at first 2- octanol have to get through there from the feeded organic phase. To intensificate the process a well-stirred reactor can be applied with as huge heat transfer area as it can be built in because of the highly exothermic reactions. Reactor Reactor insight Jacket Aqueous phase A P B X A P X Mass transport Reaction(s) Heat transport Liquid phase W Organic phase N Figure 1: Structure of Woezik and Westerterp’s reactor model Based on the assumed reaction mechanism the overal reaction rates are calculated with the following correlations: ( ) ( ) ,cck1r cck1r aq B org P2,eff org 2 aq B org A1,eff org 1 ⋅⋅⋅ε−= ⋅⋅⋅ε−= where εorg – void fraction of the organic phase in the reactor; keff,- – effective reaction rate constant, since in both equations the concentration of reagents are considered in different phases; c- org – concentration of component indexed in subscript in organic phase; c- aq – concentration of component indexed in subscript in aqueuos phase. The effective reaction rate constants are depend on the temperature of the reactor insight and the acidity of the aqueous phase which is based on the concentration of solved nitric acid: ,ekk 0,0HR ,eff,A Hm TR E 0 ,eff,eff ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ⋅− ⋅ − −− − − ⋅= where k0eff,- – effective reaction rate constant; EA,eff,- – activation energy; R – ideal gas constant; TR – temperature of reactor insight; H0 – Hammett-acidity function. To describe the component mass trajectories during the operation the following component mass balances which are ordinary differential equations must be solved. The A component is continously feeding into the reactor during the first part of the operation so volume of reaction medium is increasing too: ( ) ,rVcB dt cVd 1 Rin A in,R org A R ⋅−⋅= ⋅ where VR – volume of reaction medium; BR,in – flow rate of reagent feed; cA in – concentration of A component in feed. 152 Because only the A component is being fed into the reactor during the operation the term which is considered this is missing from the other component mass balances: ( ) ( )21R aq B R rrV dt cVd −⋅= ⋅ ( ) ( )21R org P R rrV dt cVd −⋅= ⋅ ( ) 2 R org X R rV dt cVd ⋅= ⋅ ( ) ( ) ,rrV dt cVd 21 R aq N R +⋅−= ⋅ where cN aq – concentration of N component (nitric acid) in aqueous phase. To simulate the effect of producing heat from reactions and the jacket cooling on the temperature of the reactor a heat balance must be solved: ( ) ,QQQQ HC 1 dt dT RARCin,R rR R −−+⋅= where HCR – heat capacity of reaction medium; Qr – heat produced by reactions; QR,in – heat produced by the feeding; QRC – heat flux from the reactor insight to the jacket; QRA – heat flux from the reactor insight to ambient. Finally to calculate the modifications in jacket temperature: ( ) ,QQ HC 1 dt dT RCin,C C C +⋅= where TC – temperature of jacket; HCC – heat capacity of cooling medium; QC,in – heat produced by the feeding of jacket. 0 2 4 6 8 10 12 14 16 18 0 1 2 3 t [h] n A [kmol] n P [kmol] n X [kmol] n B [kmol] 0 2 4 6 8 10 12 14 16 18 16 17 18 19 t [h] n N [kmol] 0 2 4 6 8 10 12 14 16 18 260 270 280 t [h] TR [K] TC [K] 0 2 4 6 8 10 12 14 16 18 0 1 2 3 t [h] n A [kmol] n P [kmol] n X [kmol] n B [kmol] 0 2 4 6 8 10 12 14 16 18 16 18 20 t [h] n N [kmol] 0 2 4 6 8 10 12 14 16 18 260 270 280 t [h] TR [K] TC [K] a) Normal reactor operation – analytically b) Normal reactor operation – numerically 0 2 4 6 8 10 12 14 16 18 0 1 2 3 t [h] n A [kmol] n P [kmol] n X [kmol] n B [kmol] 0 2 4 6 8 10 12 14 16 18 14 16 18 t [h] n N [kmol] 0 2 4 6 8 10 12 14 16 18 300 350 400 t [h] TR [K] TC [K] 0 2 4 6 8 10 12 14 16 18 0 1 2 3 t [h] n A [kmol] n P [kmol] n X [kmol] n B [kmol] 0 2 4 6 8 10 12 14 16 18 14 16 18 20 t [h] n N [kmol] 0 2 4 6 8 10 12 14 16 18 300 350 400 t [h] TR [K] TC [K] c) Reactor runaway occurs – analytically d) Reactor runaway occurs – numerically Figure 2: Calculated trajectories of state variables 153 Results and discussions The model was solved in MATLAB with a numerical solver based on Runge-Kutta method. To illustrate the complex dynamical behavior of process the following two experiments are performed. Trajectories of state- variables in Fig. 2a-b show the normal operation of the reactor when 2-octanone is the main product and the total quantity of byproducts is low. A small, only 3 K change in the inlet temperature of the jacket results the development of reactor runaway (see in Fig. 2c-d). In this case byproducts are mainly generated during the operation while the conversion of 2-octanol is significantly decreased at the end of the operation. In Fig. 2 grey areas represent the instability regimes of the process. These areas are determined by using Ljapunov’s indirect stability analysis. It can be seen that in a small time interval the reactor becomes instable still in normal reactor operation (Fig. 2a) applying analytic calculations. This is not thermal instability, which is a synonym of reactor runaway, but it can be originated from the autocatalytic reaction scheme. The main difference between analytically and numerically calculated eigenvalues is that while in the analytically calculation all the eigenvalues are arranged from the Jacobian-matrix of the developed process model, in the numerically calculation only the Jacobian- matrix is generated and eigenvalues are determined by using the solver of MATLAB. In Fig. 2 on the left hand side can be found the result of analytical analysis while the numerical is shown on the right hand side. Comparing results seen in Fig. 2 it can be seen that the sensitivity of the analytical calculation is more higher than the numerical technic. It means that runaway detection based on stability analysis and calculating the eigenvalues of Jacobian matrix detects the development of runaway sooner than using a numerical solver to calculate eigenvalues. Other very important thing that the algorithm based on analytical solution is three times faster than the other algorithm. This difference in calculation times makes the analytical technic more applicable in on-line monitoring system or in an OSS. Conclusions and future work Process operators usually have only temperature measurements to follow the operation and to check and keep the safe way of operation. Therefore a tool based on process models and stability analysis can be very useful in checking the appropriate operation and helping their work. The paper investigates that numerical analysis of characteristic equation of a dynamic process model how reliable in detection of reactor runaway. Analytical solution is not just more accurate but much more faster than the other numerical one. Therefore the message of this work is that with a bit more use of mathematic a lot more reliable and applicable tool can be developed to forecast or simply to detect the development of reactor runaway. Further research will focus on how the extracted knowledge can be transformed into a set of constrains on the process variables and how these constrains can be applied in batch-to-batch and in feeding trajectory optimization. ACKNOWLEDGEMENTS The authors would like to acknowledge the financial support of the Cooperative Research Centre (VIKKK, project 2004-III/2), the Hungarian Research Found (OTKA 49534), the Bolyai Janos fellowship of the Hungarian Academy Science, and the Öveges fellowship. REFERENCES 1. WOEZIK B. A. A., WESTERTERP K. R.: Chemical Engineering and Processing 41 (2001) 59-77 2. ALBERT J., LUFT G.: Chemical Engineering and Processing 37 (1998) 55-59 3. KAO C. S., HU K. H.: Journal of Loss Prevention in the Process Industries 15 (2002) 213-222 4. HENDA R., MACHAC A., NILSSON B.: Chemical Engineering Journal 143 (2008) 195-200 5. VÉCHOT L., BIGOT J. P., TESTA D., KAZMIERCZAK M., VICOT P.: Journal of Loss Prevention of Process Industries 21 (2008) 359-366 6. WOEZIK B. A. A., WESTERTERP K. R.: Chemical Engineering and Processing 39 (2001) 521-537 7. EIZENBERG S., SHACHAM M., BRAUNER N.: Journal of Loss Prevention in the Process Industries 19 (2006) 754-761 8. BARKELEW C. H.: Chemical Engineering Progress Symposium Series 25 (1959) 37-46 9. ADLER J., ENIG J. W.: Combustion and Flame 8 (1964) 97-103 10. LACEY A. A.: International Journal of Engineering Science, 21 (1983) 501-515 11. VARGA T., ABONYI J., SZEIFERT F.: Acta Agraria Kaposvariensis 10 (2006) 121-133 12. VARGA T., SZEIFERT F., ABONYI J.: Acta Agraria Kaposvariensis 11 (2007) 175-186 13. SASTRY S.: Nonlinear systems: Analysis stability and control, Springer (1999) 14. VARDE P. V., SANKAR S., VERMA A. K.: Reliability Engineering and System Safety 60 (1998) 53-69 15. GÓES A. G. A, ALVARENGA M. A. B., MELO P. F. F.: Reliability Engineering and System Safety 87 (2005) 149-161 16. LYAPUNOV A. M.: International Journal of Control 55 (1992) 531-773 17. CASTELLAN A., BART J. C. J., CAVALLARO S.: Catalysis Today 9 (1991) 255-283