Format And Type Fonts CCHHEEMMIICCAALL EENNGGIINNEEEERRIINNGG TTRRAANNSSAACCTTIIOONNSS VOL. 29, 2012 A publication of The Italian Association of Chemical Engineering Online at: www.aidic.it/cet Guest Editors: Petar Sabev Varbanov, Hon Loong Lam, Jiří Jaromír Klemeš Copyright © 2012, AIDIC Servizi S.r.l., ISBN 978-88-95608-20-4; ISSN 1974-9791 DOI: 10.3303/CET1229223 Please cite this article as: Musabekova L. M., Dausheeva N. N. and Jamankarayeva M. A., (2012), Methodology of calculating reaction-diffusion processes with moving boundaries of kinetic zones, Chemical Engineering Transactions, 29, 1333-1338 1333 Methodology of Calculating Reaction-diffusion Processes with Moving Boundaries of Kinetic Zones Leila M. Musabekovaa, Nurjamal N. Dausheevaa, Madina A. Jamankarayevab a State Univ. of South Kazakhstan, Dep. Inform. Technology, Tauke Khan av. 5, 160012, Shymkent, Kazakhstan; tel. 007 725 2210933, e-mail: mleyla@bk.ru b International Humanitarian and Technical University, Department of Computer Science and Mathematics, Baitursynov street 80, 160000, Shymkent, Kazakhstan; e-mail: d_madina08@mail.ru The main factors which promote arising of kinetic zones with moving boundaries in perfect systems have been considered in our work (Musabekova and Brener, 2004). In this work we submit the methodology for calculating the reaction - diffusion processes in non-perfect mixtures. Two cases, namely, the case of instant chemical reactions and the case of reaction with finite velocity have been considered. 1. Methodology of calculating the sorption process accompanied by the instant chemical reaction The mathematical model of the sorption process in the gas-liquid system accompanied by the instant irreversible reaction with moving reaction front describes two time periods. The first period is the time between zero point and the moment when the reaction front became moving. The second period describes the process in the system with moving reaction front. Unlike the previous work (Musabekova and Brener, 2004), the new developed model (Musabekova, 2011) takes into account the intermolecular interaction in reaction- diffusion systems. As the main tools for the numerical solution of the appropriate mathematical model the Krank-Nikolson's and Newton- Rafson's modified methods with sweep separate and bisection methods have been used. For accounting the intermolecular interaction in non-ideal systems, it becomes necessary to use the dependence of diffusion coefficients on concentrations of reaction products, what can be made by the help of the special parameter w (Musabekova L.M., 2011):  XwDD AXi 21 ~  ,   XXAAAXAXw   2 , (1) Where D - coefficient of diffusion of a component in real system; iD - coefficient of diffusion for ideal system; AX XX - energy of interaction between molecules of reagents A and X, A and A, X and X accordingly; - the parameter depending on model of a liquid state. As a result of rearrangements (Musabekova and Brener, 2004) we obtained the equations of convective diffusion of components B and E, valid at Х >0, t>0: . ~~ , ~~ 2222 2222 tCXCDXCD tCXCDXCD EEEEBEB BEBEBBB   (2) The expressions for calculating concentrations of the active component of absorbent B and the product of reaction Е in a liquid layer at time foregoing to the moment of arising of the reaction moving front are: mailto:mleyla@bk.ru mailto:d_madina08@mail.ru 1334                  tSХХerfctSХtSTTRRCC AB 24exp22 ~ 11 2 1212  (3)                        BA CtSХХerfctSХtSRRRTRTC 24exp22 22 2 2122112  ,                   tSХХerfctSХtSRRTTCC AE 24exp22 ~ 11 2 111212  (4)                        tSХХerfctSХtSRRRTRTCA 24exp22 22 2 22122112  . The expressions for concentrations *),( tХCB and *),( tХCE are used as initial data for the calculation of concentrations profiles of components A, B and E with allowing for the influence of inter-molecular interaction between reaction products at the period *tt  . The reaction front becomes moving at the moment of time *tt  . However, under the condition *tt  mathematical model (2), has no analytical decisions. Therefore the adaptive modified method with using the scheme of Krank-Nikolson has been developed and the special numerical experiment has been carried out. We used also both the sweep separate method for definition of profiles of components concentration A, B, E and the modified iterative method for controlling the concentration profiles of reaction products. Let’s consider application of the Krank-Nikolson’s method in this scheme. We had a grid of points at the following levels of time: lt , lll ttt 1 . We introduce the following denotations for components concentration B and E: ),( ji j i txBB  , ),( ji j i txEE  respectively. Then using the finite differences method according to the implicit scheme, we obtain the calculation scheme for a component B:        j i j i j i j i j i j iBB bBBbaaBbBBbaaBD 11 1 1 11 1 )()( ~ (5)    abba t BB bEEbaaEbEEbaaED j j i j ij i j i j i j i j i j iBE                   1 11 1 1 11 1 )()( ~ . Simplifying the expressions we obtain the numerical mathematical model for calculating the absorbent concentrations profiles at time *tt  :                 j i j i j iBE j i j i BB j iBB bEEbaaEDbBB D baaBD 11 1 1 11 1 )( ~ )~( ~  (6)  j i j i j iBE j i j i BB j iBB aEEbabEDaBB D babBD 1111 )( ~ )~( ~             . where: ),( ji j i txAA  , ),( ji j i txEE  are concentrations of components A and E. With using the equation of diffusion (2) and the finite differences method according to the implicit scheme, the scheme for the components A and E at the period looks as follows. For solving this problem both at the first area where the concentration profiles of components A and E are calculated, and at the second area where the profiles for components B and E are calculated, we used sweep separate method. According the approach the finite difference equations (5) - (6) are divided on two groups. The first group is used for calculating the component concentration B. The second group is used for calculating the component concentration E. Equations, which choose in the first group, have been solved by the sweep separate method. The concentration of component E we left to be constant in calculating process, and the appropriate value has been taken from previous layer. When calculating for the first group is finished, the obtained value of the component concentration B is passed to the second group for calculating the component E concentration with the help of the sweep method. The last value is used in external iteration at j+1-th layer for the first group. 1335 At realization of the convergence condition the iteration was finished. The calculated values of components concentration В and Е, А and Е are used on j+2-th layer and etc. For calculating of the sweep coefficients the boundary conditions for the first and second areas were used separately. Calculations of concentration profiles of the component E at the surface of the reaction front show, that the calculated concentrations of E (I) and E (II) didn't be equal nearby the reaction front on the left and on the right. Therefore for increasing the accuracy of calculations we controlled the concentration of reaction product Е so that the following conditions were satisfied:  )()( BAEBA JJJJJ ,  ESEÅ MM 21 , (7) where 2)21( MM EEES  , and AJ , BJ , EJ are the flows of components А, В and Е; and 1Å , 2Å are the concentrations of components; M is the mesh point of the concentration E at surface of the reaction front. While calculating the concentration profiles of the component A we controlled concentration of the absorbed component aC , varying the time step t , so that the following conditions were satisfied: 3 22 10),(    jjsA ttÕC , mol/m 3 (8) Step by time t for every level of time j we calculated by bisection method. The algorithm of sweep separate method consists of follows. The previous concentrations A and B for second area were unknown, therefore we took the calculated values of concentrations from the previous part, i.e. without taking account of influence of reaction product, as initial approximations. At achievement of the convergence condition the iteration cycle has been finished. а) w=0 б) w=0.03 0 0 . 0 0 0 5 0 . 0 0 1 0 . 0 0 1 5 0 . 0 0 2 0 . 0 0 4 7 9 8 6 3 . 1 8 9 6 1 6 3 6 . 3 8 4 0 3 1 1 9 . 5 7 8 4 4 6 Т0 T1 T2 T3 T4 T5 Т6 Длина реак т ора L , м К о н ц е н т р а ц и и С а , С b Figure 1- Profiles of concentrations of components A and B along the length L of the reactor at different time T at а) w=0, б) w=0.03. Ca, Cb - concentrations of components A и B in mol/m 3 0 0.0005 0.001 0.0015 0.002 0.0028849  3.1965956 6.3960761 9.5955566 Т 0 T1 T2 T3 T4 T5 Т 6 Длина реактора L, м Са , С b Length of reactor L, m C o n c e n tr a ti o n С а , С в , m o l/ m 3 Length of reactor L, m C o n c e n tr a ti o n С а , С в , m o l/ m 3 1336 Figure 1 depicts the obtained received results of numerical experiments for values of parameters w = 0 and w = 0.03. The components concentration A and B are shown in different time T0-T6 starting from time of intrusion of reaction front t* with time interval t, where abscissa is a length of the reactor L. We can see that the concentration of the reaction product at the reaction front increases but the time of formation of the reaction front t* decreases with increasing parameter w. So, the time at which the reaction front begins to move into the liquid layer became less. It was established that the conversion and output of the product reaction increase with increasing concentration of component B and the parameter w. We can see also from the results of calculations that the non-ideality of the system essentially influences on the characteristics of the process; therefore the non-ideality must be taking into account at calculating the absorption accompanying by the chemical reactions. The results of carrying out researches can be used in the engineering method of calculating the intensity of chemical reactions and under designing the non-isothermal through-reactors with minimal reaction length. The numerical data we have processed by the methods that are adopted for the analysis and processing of data sets of natural experiments with the assessment of accuracy. The coefficients of the acceleration of mass-transfer at chemisorption on the basis of the well-known Sherwood film model we defined as follows: 0=1+r0S0, 1=1+r1S1, (9) where r0= BA DD ~~ , S0=B/AS, (10) r1= BE DD ~~ , S1=B/Еfr, (11)  XwDD AXi 21 ~  ,   XXAAAXAXw   2 , (12) Characteristic modified number of the Sherwood is: ShМ =  pA hD  ~ (13) Mass transfer coefficient of a captured component in the liquid phase has been determined by the characteristic depth and by the penetration time for the two stages of process: for ttp:   21~~ 4 AAL DerftDa  , (15) where -parameter of the film model. For the first stage of the process the period of the formation of moving reaction front, taking into account of the intermolecular interaction, reads:    2 1 12 2 2112222 12 2* 22 4              S TT S RTRT CRRCt AB  , (16) where: BEBB BEBB DD DD R ~~ ~~ 2 1 1      , EEEB EEEB DD DD R ~~ ~~ 2 1 2      , BEBB DD T ~~ 1 2 1   , EEEB DD T ~~ 1 2 2   . (17) The obtained expressions for the second stage of the process are: 1.The characteristic penetration time is:        52,3 2121018,0132,1   BABBBAABAABp CwDCwDPHCt . (18) 2. The characteristic penetration depth is:         772,017,0 5 21 21 21 21 11069,7                        BBBB AAAA BBBB EEEE p CwD CwD CwD CwD h . (19) 3. The concentration of the absorbed component at the interface is follows:     402,0798,0 369,2 PBBA ttPHCCC  . (20) 4. The average mass-transfer coefficient in the liquid phase at the initial stage of moving reaction front: 1337             HPHCCP BB 811,0 358,2445,2 . (21) In the formulas the following notations are used: BA CC , are the concentrations of the absorbed component and active component in absorbent; EBA DDD , , are the diffusion coefficients of the absorbed component, the active component and the reaction product, respectively;  - the mass- transfer coefficient; H - Henry’s constant; P - pressure in the gas phase. Check of the adequacy of the formulas shown, that they provide a calculation error which does not exceed 20 percent under the entire range of parameters. The third stage of the process in accordance with the film model was calculated. The obtained above expressions were used as the characteristic scales in the assessment of the control parameters of the film model. On the base of the proposed method we can give the reasonable assessments of the applicability of the film model to the fast processes. The process can be calculated by applying the film model only if the characteristic time of the process in an apparatus is much larger than the characteristic time of penetration. Otherwise, we must take into account the stage of the growth of the velocity of the concentration front. 2. Calculation of the heat and mass transfer under autowave process Description of heat and mass transfer in chemical apparatuses under auto wave regimes of phases interaction associates with very great mathematical difficulties. In this part of our paper we investigate the possibility of the appearance of dissipative structures in non- isothermal through-reactors. So, we submit the mathematical model which includes the diffusion-kinetic equations for the two reagents at the availability of a reversible first-order reaction and the heat transfer equation taking into account the heat of reaction in the non-isothermal reaction-diffusion systems in tubular through-reactor (Brener, 2005): BAAAAA CkCkzCSjzCDtC 21 22~  , (22) BABBBB CkCkzCSjzCDtC 21 22~  , (23) pcHzTSjzTtT   22 (24) where AC , BC are the concentration of components A and B, respectively, 0AC is the concentration of reagent A at the reactor input; ÀD ~ , ÂD ~ are diffusion coefficients of reagents; zt , are time and space coordinates, respectively, j is the total consumption of reagents through the reactor, T is the temperature;  is the average coefficient of thermal conductivity;  is average density of reagent mixture ; pc is a average mixture heat; H is the total heat of reaction, S is cross-sectional area of the reactor. The model is designed for calculating heat and mass transfer in non-isothermal tubular chemical reactor, which allows taking into account the existence of regimes of reactor work, accompanied by the appearance the moving wave fronts of the volume of apparatus. We have considered both the case of reactors with heat transfer between apparatus and outside media, and the case of the adiabatic reactor. Stationary regimes and conditions of their stability are investigated by the methods of hydrodynamic stability, and with using of a numerical experiment too. As a result, the control parameters have been determined, including thermal, diffusion and kinetic characteristics for both reagents, taking into account the two stages of reactions (Musabekova, 2008). Here the sequence of analyzing the model has been shown: 1. Input of required data for calculation. 2. Calculating and plotting of the concentrations of components X and Y (A and B in the graphs) and the temperature T along the length of the reactor with taking into account: 1) results of the analysis of stability of stationary regimes in non-isothermal tubular chemical apparatus and conditions for creating the moving concentration and temperature fronts with known rates of forward and 1338 backward stages of chemical conversions in the reactor; 2) wave characteristics of periodic dissipative structures and the concentration and temperature fronts in tubular reactors. 3. Processing of the results of calculation, namely: plotting amplitude of oscillations for components A and B in the time. 4. Carrying out a series of numerical experiments for different values of control parameters. 5. Saving data of calculation. Let’s submit the main formulas for engineering calculation methods (Brener, 2006): 1. The expressions for calculating the wave characteristics of the dissipative structures in the form of moving circulating cells:           L mtA sinexp1 ,           L mt sinexp3 ,           L mtB sinexp2 (25) 2. The minimum length of the reactor, for which can be implemented such cells is defined as follows:         YX DKDK hK L ~~ 21 3 min   . (26) 3. Conclusions The algorithms and the methods for calculating the chemical reactors with fast chemical reactions, which can lead to propagation of moving concentration fronts in apparatuses volume, have been proposed. It was shown that the velocity of moving structures strongly influences on the intensity of transfer processes. Two cases are described, namely: instantaneous chemical reactions and reactions occurring at a finite velocity. The developed method of calculation for the case with instantaneous chemical reaction can be recommended for use in calculations of chemisorptions processes, characterized by the availability of phase fast reaction with the formation of moving reaction front (in particular, the absorption of iodine compounds). Using the proposed method it is possible to get the reasonable assessments of the applicability of the film model to the fast processes. The developed method of calculation for the case of the reaction occurring with a finite rate confirmed the influence of moving auto-wave structures on the intensity of heat and mass transfer in the system. References Brener A.M., Musabekova L.M., 2006, Autowave regimes of heat and mass transfer in the non- isothermal through-reactors, Advanced Computational Methods in Heat transfer IX, Wessex Institute of Technology, Published by WIT Press Ashurst Lodge, Ashurst Southamption SO40 7AA, UK, 181-191. Musabekova L.M., 2011, Method for calculating the reaction-diffusion processes with account of non- ideal systems, International Scientific Conference NERPO, MGOU, Moscow, 294-299 (in Russian). Musabekova L.M., Brener A.M., 2004, The methodology for calculating the process of chemisorption in systems with a moving front of the instantaneous irreversible reaction, Works of the 5th Minsk International Forum on Heat and Mass Transfer MMF, Minsk, T.1, 319-323 (in Russian). Brener A.M., Serimbetov M.A., Musabekova L.M., 2005, Non-local equations for concentration waves in reacting diffusion systems, Twelfth International conference on Computational Methods and Experimental Measurements XII, WIT Press Southampton, Boston, 93-103. Musabekova L.M., Yunusova A.A., Dausheeva N.N., Eskendirov Sh.Z, 2008, Modeling of through- reactors with allowance for the characteristics of phases distribution, 18th International Congress of Chemical and Process Engineering.-Prague, Czech Republic, August 24-28th, 131-132.