2516
Ind. Eng. Chem. Res. 2003, 42, 2516-2524
Modeling of Complex Organic Solid-Liquid Reaction Systems in Stirred Tanks Tapio Salmi,*,† Esko Tirronen,‡ Juha Lehtonen,§ Erkki Paatero,| and Daniel Valtakari† Process Chemistry Group, Laboratory of Industrial Chemistry, Åbo Akademi, FIN-20500 Turku/Åbo, Finland, Espoo Research Centre, Kemira, Box 44, FIN-22071 Espoo, Finland, Neste Engineering, Box 310, FIN-06191 Porvoo, Finland, and Laboratory of Industrial Chemistry, Center of Separation Technology, Lappeenranta University of Technology, Box 20, FIN-53851 Lappeenranta, Finland
Mathematical modeling of gas-solid reactions is very common in chemical engineering, but liquid-solid reactions have received much less attention. Modeling of liquid-solid reactions is not studied very much, but it is gaining importance because the principles of chemical reaction engineering have made a breakthrough in industries dealing with fine and specialty chemicals. General principles for the modeling of complex organic liquid-solid reaction systems in backmixed batch and semibatch reactors are presented. The ingredients of the model are complex kinetics, interfacial mass transfer, and models for solid particles and the reactor. Previously developed theories for solid particles, such as shrinking-particle and porous-particle models, can be applied on liquid-solid systems, but the kinetics is complex and the reactions are carried out (semi)batchwise, which implies that numerical solution algorithms are used. The approach is illustrated with two experimental case studies representing shrinking and porous particles, Claisen condensation and carboxyalkylation of cellulose. Following the modeling principles, it was possible to describe the experimental data very well and predict the behavior of solid particles and the reactor performance. Detailed modeling of the liquid-solid reactions will be used extensively in process scale-up in the future. Introduction Gas-solid reactions are studied very much in chemical engineering, particularly because of their importance in combustion processes. The progress of a gas-solid reaction can be followed by measuring the weight change of the solid reactant. The weight change is related to the reactant conversion. The reaction time vs conversion data can be compared with existing models for reactive solid particles. Many of the models, such as shrinking-particle models and product (ash)layer models, are presented in the classical textbook of Levenspiel.1 Later on, more detailed models, such as grain models, have been developed.2 All of these approaches are based on the concept of a two-dimensional reaction surface in the solid particle. It is, however, possiblesbut not very usualsthat the entire porous particle is under reaction conditions simultaneously. This implies that theories developed for simultaneous reaction and diffusion for porous catalyst particles become applicable for the system. Chemical systems where solid particles react with organic liquids have received much less attention than gas-solid systems, but they are gaining more and more importance because the principles of chemical engineering have made a breakthrough in the production of fine and specialty chemicals, such as production of pharma* To whom correspondence should be addressed. Fax: +3582-2154479. E-mail:
[email protected]. † Åbo Akademi. ‡ Kemira. § Neste Engineering. | Lappeenranta University of Technology.
ceuticals and modification of alimentary products. Liquid-solid reaction systems have some characteristic features to be included in mathematical modeling. Because organic compounds are involved, complex reaction mechanisms play a key role which leads to too complex rate expressions. Thus, model simplification is a vital part of the work to reduce the number of adjustable kinetic parameters in the rate expressions. The activity of the solid material might change during the reaction because of surface restructuring or sterical hindrances on the solid surface due to substitution reactions. To optimize the amount of the solid material in the chemical reactor, reliable models for the solids are needed. On this point, all of the theories developed for gas-solid reactions are helpful. The development of microscopic techniques enables in some cases the model verification for the solid material, but different models often have to be fitted to the experimental data to reveal which one gives the best description of experimental data. Liquid-solid reactions are mostly carried out in stirred-tank reactors operating in batch or semibatch mode. This implies that the concentrations in the liquid phase change continuously and a reactor model is needed. The simple test plots of conversion vs time for solid particles are not directly applicable because they are based on constant concentrations in the bulk phase surrounding the particles. Until the present moment, the description of stirred tanks is based on the concept of complete backmixing, which is achievable on a laboratory scale but not always on a large industrial scale. In the future, a lot of effort should be put into the application of computational fluid dynamics (CFD)
10.1021/ie0206296 CCC: $25.00 © 2003 American Chemical Society Published on Web 01/31/2003
Ind. Eng. Chem. Res., Vol. 42, No. 12, 2003 2517
on reactive liquid-solid systems.3 For the description of laboratory-scale data and for first considerations of large-scale reactors, the simple approach based on complete backmixing will persist. The models are dynamic (time-dependent) and usually highly nonlinear, which implies that nice analytical test plots cannot be applied on experimental data but efficient numerical strategies are needed. The purpose of the present paper is to give a simple overview of the treatment of complex organic liquidsolid reactions in completely backmixed batch and semibatch reactors and to illustrate the modeling concepts with a few case studies. Complex Reaction Kinetics
ri )
∑νijrj
A(l) + S(s) ) A1(l) + P1(l/s)
(a)
A1 + B(l) ) A2(l) + P2(l/s)
(b)
A2 + C(l) ) P3(l/s)
(c)
where A, B, and C are liquid-phase reactants, P1-P3 reaction products, and A1 and A2 short-life intermediates. The intermediates have such low concentrations that they cannot be detected by standard analytical techniques such as GC and HPLC. The overall reaction directly visible for the researcher is
A + B + C + S ) P1 + P2 + P3 The concentrations of the intermediates are eliminated by applying quasi-steady-state or quasi-equilibrium hypotheses. Recently, a more general approach, asymptotic analysis, has been proposed by Haario et al.4 The use of a quasi-steady-state hypothesis implies that the generation rates of the intermediates are set to zero.
(1)
Here, for example, rA1 ) r1 - r2 ) 0 and rA2 ) r2 - r3 ) 0. Such a treatment leads to a very general rate expression, provided that each reaction step is linear with respect to the intermediates:
∏a1a2...aSR - ∏a-1a-2...a-SR)/(∑(∏akalam...)
r)(
(2)
where the numerator represents the mass-action term and the denominator consists of a sum contribution of selected forward and backward reaction steps. The details are explained in the appendix. The concept is illustrated for substitution reactions between solid and liquid, namely, Claisen condensation. For nonlinear reaction steps, i.e., when two or more intermediates react in the same step, an analytical solution for the reaction rate becomes mostly very cumbersome or even impossible.
(3)
Mass Balances for Liquid-Phase Components The mass balances for liquid-phase components will be discussed briefly in this section. We assume that the liquid phase is completely backmixed, but some of the components are fed into the reactor during the experiment. This leads de facto to a semibatch reactor model. The mass balance for an arbitrary component is given by
dni/dt ) n′0i + NiA + riV
Organic systems, such as substitution and elimination reactions, typically involve several steps, where shortlived intermediates play a key role. Typical examples of reaction intermediates are radicals and ionic species, e.g., carbenium ions and carbanions. A simple sequence of reaction steps can be sketched, for instance, as follows:
rAj ) 0
The generation rates of the components are obtained from the reaction rates and the stoichiometry in a straightforward manner:
(4)
where the symbols are defined in the Notation section. It should be noticed that the liquid volume changes because of the feed. Provided that the liquid density remains virtually constant, the volume can be updated with the simple formula
V ) V0 +
∫0t V′ dt f
(5)
where tf is the feed time. After the feeding has been completed, the reactor operates batchwise. As indicated by eqs 4 and 5, it is wise to use an amount of the substance (n) or some dimensionless variation of it in calculations to keep the differential equations explicit. The concentrations needed in the rate expressions are then calculated from ci ) ni/V. The expression for the interfacial flux (Ni) depends on the particular case. Here we consider two examples relevant for complex organic systems: shrinking solid particles and porous particles with simultaneous reaction and diffusion. Mass Balances for Solid Components: Nonporous Shrinking Particles For the sake of simplicity, particles of equal sizes are considered in the sequel, but the treatment can be easily generated to particle size distributions. The chemical reactions are assumed to take place in the liquid phase at the outer surface of the particle. The solid material has a limited solubility in the liquid. The mass balance of the dissolved solid material in the backmixed tank reactor is
dnS/dt ) rSV + NSA
(6)
On the other hand, the mass balance of the solid material itself can be written in the form
dnSs/dt ) -NSA
(7)
A simplified treatment becomes possible when the balances are added together:
dnSs/dt + dnS/dt ) rSV
(8)
Often the solid component is very sparingly soluble and a pseudo-steady-state hypothesis is applied on it, which implies that dnSs/dt . dnS/dt. Thus, the amount of solid
2518
Ind. Eng. Chem. Res., Vol. 42, No. 12, 2003
is directly related to the generation rate of the component: s
dnS /dt ) rSV
(11)
The symbols are explained in the Notation section. For numerical treatment, it is better to introduce the transformation
z ) xa
(12)
Differentiation of eq 12 gives dz/dt ) axa-1 dx/dt, and finally we get
dz/dt ) (MS/m0S)rSV
(13)
The dimensionless radius is then obtained from x ) r/R ) z1/a. De facto z is the dimensionless mass of the solid material in the reactor [z ) mS/m0S ) (nS/MS)/m0S], which is revealed by comparison of eqs 9 and 13. Mass Balances for Solid Components: Porous Particles with Constant Size Particles with equal sizes are considered, and the reaction is assumed to take place everywhere inside the wetted pores of the particle. The liquid-phase species are mobile and transported via diffusion, while the solid components are immobile on the surfaces of the pores. For the liquid-phase components, the mass balance can be written for an infinitesimal volume element (∆VLP) in the particle:
(NiAP)in + ri∆VLP ) (NiAP)out + dni/dt
(14)
After introducing the concentration in the volume element (ni ) ci∆VLP) and letting the volume element approach zero, we get
dci/dt ) -d(NiAP)/dVLP + ri
(15)
The further treatment depends on the detailed geometry and the diffusion model applied. After introduction of a characteristic dimension and shape factor analogous to eq 10
a ) (AP/VP)R
(16)
and assuming that the porosity of the particle (P ) ∆VLP/∆VP) is constant, it can be shown that the balance
(17)
As a simple case, let us consider Fick’s law for diffusion, Ni ) -Dei dci/dr. Equation 17 becomes
dci/dt ) -(Dei/P)[d2ci/dr2 + (2/r) dci/dr] + ri (18) The shape factor has the following values for ideal geometries: a ) 1 for slab, a ) 2 for long cylinder, and a ) 3 for sphere. However, noninteger values can be used for nonideal geometries. For the solid particle, i.e., the immobile functional groups in the pores, the diffusion term in eqs 17 and 18 is zero and the simple relation
dcS/dt ) rS
(10)
it is possible to show that the dimensionless particle radius x ) r/R follows the relation
axa-1 dx/dt ) (MS/m0S)rSV
dci/dt ) -P-1r1-a d(Nira-1)/dr + ri
(9)
It should be remembered that the treatment presented above makes sense only when the turbulence (stirring) around the particle is so vigorous that the concentrations at the particle surface and liquid bulk are equal. If this is not the case, both balances (6) and (7) remain, and a mathematical expression is inserted for the flux (NS). The further treatment of the mass balance depends on the detailed particle geometry. However, by introduction of a characteristic initial dimension (R) of the particle along with a shape factor (a)
a ) (A0P/V0P)R
equation takes a very general form:
(19)
is valid. Numerical Approach The mathematical models presented above are typically used for two purposes: determination of kinetic parameters from experimental data and prediction of the reactor performance. Kinetic parameters are estimated by nonlinear regression, i.e., minimization of an objective function of the type
Q)
∑ωi,t[yi,t(exp) - yi,t]2
(20)
where yi,t(exp) and yi,t denote the variables obtained experimentally and predicted by the model, e.g., the molar amounts given by eq 4. Commonly known optimization methods, such as simplex, gradient, and Levenberg-Marquardt methods,5 are used in the minimization of the objective function, and the theories developed in mathematical statistics are used in the evaluation of the goodness of the parameters.6 Two kinds of serious obstacles appear very often: the minimization method is not able to find the global optimum but remains at a local minimum, and the evaluation of the confidence intervals of the parameters with standard methods are unreliable because of the nonlinearity of the model. If the number of parameters to be determined by regression analysis is high compared to the number of measured variables, huge confidence intervals are predicted for at least some of the parameters. Therefore, graphical considerations of the parameters are very preferable. The parameters can be checked pairwise with contour plots, and the sensitivity of the objective function can be investigated with respect to each parameter by calculating its value as a function of one parameter and simultaneously keeping the other parameters at their optimal values. The approach is illustrated in Figure 1. In the determination of kinetic constants from experimental data, it is numerically useful to introduce an orthogonal transformation of Arrhenius’ law, k ) k0 exp[-(Ea/RG)(T-1 - T0-1)], where T0 is the average value of experimental temperatures. The statistical accuracy of the rate constant is improved by this transformation. Furthermore, it is much easier to give an initial guess for the rate constant at the average temperature (k0) than for the frequency factor. Most models consist of a set of ordinary differential equations (ODEs) and initial value problems (IVPs) with
Ind. Eng. Chem. Res., Vol. 42, No. 12, 2003 2519
Figure 1. Examples of contour plots (a) and parameter sensitivity analysis (b).
respect to the reaction time. Because of the stiffness of the ODE system, which is caused by the chemical complexity, i.e., there is a coexistence of slow and rapid reaction steps and spatial discretization. Thus, implicit and semiimplicit algorithms for stiff ODEs are highly preferable. Several codes have been developed: those based on the backward difference (BD) method of Henrici7-9 as well as codes based on semiimplicit Runge-Kutta methods.10-12 A specific aspect appears when diffusion limitations exist inside the solid particles. To solve the IVP, the spatial coordinate has to be removed by discretization, for instance, by approximation functions, such as orthogonal collocation,13 or by finite differences.14 The finite difference method was used here. Because the first- and second-order derivatives originate from diffusion, central difference formulas were applied.15,16 Typically, a five- or seven-point formula gave a satisfactory accuracy. Asymmetric central difference formulas were used at the outer surface of the particle and close to it. Discretization converts the partial differential equations to a set of IVPs, which were solved with stiff algorithms. The Claisen condensation case was solved with the BD method,8 while the carboxyalkylation case was solved with semiimplicit Runge-Kutta methods.11,12 The parameter estimation problem was solved with the Levenberg-Marquardt method5 implemented in parameter estimation software.17,18 In the first attempts, when then parameters were still totally unknown, the estimation was commenced with the more robust simplex method.
Figure 2. Shrinking-particle model. Case: reaction scheme for Claisen condensation.
Experimental Section Claisen Condensation. Twelve Claisen condensations were carried out at 45-85 °C in an isothermal jacketed glass reactor (2 L) equipped with a turbine stirrer (450 rpm) and a thermocouple. Chlorobenzene was used as the solvent. One of the reactants (A; see Figure 2) was fed into the reactor at a constant speed by a peristaltic pump. The temperature and feed data of A were stored on a PC. The experiment was started by purging the reactor with nitrogen, after which the ester (C) dissolved in chlorobenzene was added. Stirring was switched on, and sodium methoxide (B) was added under a nitrogen flow to prevent atmospheric contamination. The feeding time of A was typically 0.5-1 h. The molar excesses of A and B were 1.1-2.5 (A) and 1.4 (B). Samples were withdrawn from the liquid phase and analyzed with HPLC (C and D) and GC (A and E). Carboxyalkylation of Cellulose. Carboxymethylation and -ethylation of sodium cellulose were studied in an isothermal jacketed glass reactor (0.7 L) equipped with a turbine stirrer (1800 rpm) and a thermocouple. Preliminary experiments with different stirring velocities indicated that the stirring speed selected was high enough to remove external mass-transfer limitations. A cooling condenser was placed on top of the reactor to prevent the escape of the solvent. Monochloroacetic acid (carboxymethylation) and monochloropropanoic acid (carboxyethylation) were used as reagents and 2-propanol as a solvent. The temperatures of the carboxymethylation and -ethylation experiments were 30-80 and 60-95 °C. The reactor was filled with pure cellulose (12 g), and the solvent (500 mL) was added under mixing. A 50 wt % NaOH solution was added to swell the
2520
Ind. Eng. Chem. Res., Vol. 42, No. 12, 2003
Figure 3. Claisen condensation. Examples of the model fit to primary data. Feeding time of A: 1 h (in all data sets displayed in the figure). Excess of A: 1.1 (data sets 9 and 10), 1.8 (data set 11), and 2.5 (data set 12). Excess of B: 1.4 (in all data sets).
cellulose and to bring it to the sodium form (mercerization). The temperature was maintained at 20 °C. Mercerization was allowed to continue for 2 h, after which the temperature was raised to the carboxyalkylation temperature and monochloroacetic or monochloropropanoic acid was added instantaneously. Typical molar ratios of cellulose/NaOH/monochlorocarboxylic acid were 1:8:4. Samples were withdrawn from the liquid phase, HCl was added to stop the reaction, and carboxymethylcellulose (CMC) or carboxyethylcellulose (CEC) was precipitated by washing the samples with ethanol. The average degree of substitution (DS) was determined from the precipitated product by titration with sulfuric acid.19 From some samples, the side reactions of the chlorocarboxylic acid were determined by proton NMR. Claisen Condensation: Reactive Shrinking Particle The main reactions taking place in Claisen condensation are displayed in Figure 2. The scheme follows exactly the general pattern presented by reaction formulas a-c in the section Complex Reaction Kinetics. In the present case the reaction steps are formulated as follows:
A + B f A1 + E A1 + C f A2 A2 f A3 + B A3 + B f D + E which give the overall reaction A + B + C f D + 2E; i.e., the final condensation product (D) and methanol (E) are formed. The task is to predict the concentration profiles in the liquid phase along with the consumption of the solid material, methoxide (B). The rate equation was simplified by assuming the overall reaction to be
Table 1. Kinetic Parameters Claisen Condensationa k at 418 K/[(dm3)2/mol/min] Ek/(kJ/mol) K′
1740b 56.1 3.04
Carboxymethylation of Cellulosec T/K
k0c0/(1/min)b
303 313 333
0.0036 0.0049 0.0114 D′ ) 0.0071 (T/313 K) (1/min)
a
Where Ek is the activation energy of k. Parameter K′ turned out to be temperature-independent in practice. b The standard deviations of the parameters were a few percentages only in both cases. c Where c0 is the initial concentration of hydroxyl groups and D′ is the internal diffusion parameter (effective diffusion coefficient to particle radius).
irreversible. The hypothesis was verified by experimental observations: a complete conversion of the limiting reactant (C) was always achieved. The final form of the rate equation used in regression analysis became
R ) kcAcBcC/(cC + K′cE)
(21)
It should be kept in mind that both k and K′ represent merged parameters, which contain products and sums of fundamental rate and equilibrium constants. The results of parameter estimation are briefly summarized in Table 1, where the merged parameters are also explained in detail. An example of the fit of the model to primary data is displayed in Figure 3. As the figure reveals, the model provides a perfect description of the data. The accuracy of the parameters was good too, and the standard deviations of the modified frequency factor and the activation energy were a few percent only. An important aspect is the behavior of the solid material. First, it is necessary to assume a numerical
Ind. Eng. Chem. Res., Vol. 42, No. 12, 2003 2521
Figure 5. Porous-particle model. Case: reaction scheme for carboxyalkylation of cellulose.
Figure 6. Carboxyalkylation of cellulose. Side reactions of chlorocarboxylic acid.
Figure 4. Claisen condensation. The average size of the solid particles during the reaction.
value for the generalized shape factor (a). Images obtained by scanning electron microscopy revealed that the particles are best represented by cylinders with spherical ends (Figure 4). This gave a value of a ) 2.5. Test calculations showed that the results were not very sensitive with respect to the a value. Thus, an average a value is used throughout the calculations. After these considerations, the consumption of solid material was computed as a function of the reaction time. An example is provided in Figure 4. The figure predicts the consumption and shows the excess of B remaining in the reactor. The simulated picture is in agreement with experimental observations. The model developed here has been used in the optimization of industrial production units. Carboxyalkylation of Cellulose: Reactive Porous Particle The principal reaction scheme for carboxyalkylation of cellulose is displayed in Figure 5. The chlorocarboxylic
acid anion has a partial positive charge in the substituted carbon atom, which makes it susceptible to a nucleophilic attack of the ionized hydroxyl groups of cellulose. In principle, all of the hydroxyl groups of the anhydroglucose unit (OH-2, OH-3, and OH-6) are reactive, which leads to a maximal DS ) 3. The multiple substitution gives, in practice, a complex reaction system with mono- (2, 3, 6), di- (23, 26, 36), and trisubstituted (236) units as products (Figure 5). In addition, the chlorocarboxylic acid undergoes side reactions in the liquid phase. Sodium cellulose provocates elimination of the chlorine atom in longer chain (C3 and longer) carboxylic acids. Besides this side reaction, the hydroxyl and alkoxide ions present in the liquid phase originating from dissolved alkali (NaOH) and solvent alcohol (usually 2-propanol) initiate substitution and elimination reactions with the chlorocarboxylic acid. The side reactions are summarized in Figure 6. They are highly undesirable because they destroy the valuable reagent, the chlorocarboxylic acid. The main and side reactions are comprised of bimolecular steps. Therefore, it is reasonable to assume second-order reaction kinetics. The generation rates of substituted (i) and unsubstituted hydroxyl (0i) groups are obtained from
ri ) kic0icHA
(22)
r0i ) -kic0icHA
(23)
where HA denotes the chlorocarboxylic acid (monochloroacetic or monochloropropanoic acid). The rate of
2522
Ind. Eng. Chem. Res., Vol. 42, No. 12, 2003
chlorocarboxylic acid consumption is expressed with
∑kic0i + k′1c0R′ + k′2cOH)cHA
rHA ) -(
(24)
It should be noted that k′1 and k′2 represent lumped constants which include the effects of both substitution and elimination reactions: k′1 ) k′1N + k′1E and k′2 ) k′2N + k′2E. The DS is obtained from the concentrations of OH-2, OH-3, and OH-6 groups.
DS ) ci/c0
ci ) c0i, i ) 2, 3, 6
(25)
In the case of diffusion resistance inside the particle, DS is a local quantity. The average DS is obtained by numerical integration over the entire particle (x):
∫
DS ) a (ci/c0)xa-1 dx
(26)
In the present case, we assumed a constant porosity for the solid material because only a part of the anhydroglucose units were substituted by the monochlorocarboxylic acid. In addition, the substituted product remained partially in the solid material. Thus, there was no real evidence to account for the change of the porosity. The approximative concept introduced by Flory,20 namely, the independence of the reactivity of a functional group of its environment, was adopted in the present case. This implies that just three basic rate constants are needed, for each hydroxyl group in the anhydroglucose unit. The ratio of the rate parameters (k2, k3, and k6) was fixed a priori, based on separate analysis data from the hydrolysis of substituted cellulose.21 Thus, the rate constants obtained the form
ki ) Rik0
Figure 7. Carboxyalkylation of cellulose. Examples of the model fit to primary data: carboxymethylation (a); carboxyethylation (b).
(27)
where Ri denotes the reactivity ratio. The following reactivity ratio was used in this work: R2:R3:R6 ) 1:0.3:1. Typical results of carboxymethylation and carboxyethylation are displayed in parts a and b of Figure 7, respectively. The model parameters were determined with regression analysis according to the principles explained in the modeling sections. As revealed by Figure 7a,b, the model gives a good overall fit to primary data, the average DS. Furthermore, the model predicts the concentration profiles of monochlorocarboxylic acid (HA) as a function of the spatial coordinate and reaction time (Figure 8). Analogous concentration profiles are obtained for the concentrations of substituted (OH-2, OH-3, and OH-6) groups (Figure 9). Figures 8 and 9 show clearly that diffusion resistance plays a role in substitution of cellulose and has to be included in modeling, either with a real reaction-diffusion model presented here or with a semiempirical model described in a previous paper.22 According to the reaction scheme (Figure 5), the basic model gives also the distribution of mono-, di-, and trisubstituted units, which is important for the product quality of substituted cellulose. After estimation of the kinetic and diffusion parameters (Table 1), it is possible to use the reaction scheme (Figure 5) and compute the concentrations of mono-, di-, and trisubstituted units. The stoichiometric details are not listed here because they are trivial and directly based on Figure 5. A simulation example is presented in Figure 10, in which di- and monosubstituted units dominate in the final
Figure 8. Carboxyalkylation of cellulose. Concentration profiles of chlorocarboxylic acid (HA) during the reaction inside the porous particle.
reaction product of carboxymethylation but also trisubstituted products exist in the mixture. In carboxyethylation, monosubstituted anhydroglucose units strongly dominate in the product. Simulation results displayed in Figure 10 are useful in the optimization of the substitution process: depending on the specific application of CMC, mono- or disubstitution is preferred. Conclusions General principles for the modeling of reactive solid materials in organic liquids were presented. Two models for solid particles were considered: shrinking particles and porous reactive solids. The models for complex reaction kinetics and mass transfer were coupled to reactor models, i.e., mass balances for completely backmixed batch and semibatch reactors. The modeling concept was applied to two case studies: Claisen
Ind. Eng. Chem. Res., Vol. 42, No. 12, 2003 2523
Appendix. General Rate Equations for Complex Kinetics
( ( ( (
)
r1 ) k1 cAcB -
cA1cE ) a1 - a-1cA1 ) r K1
r2 ) k2 cA1cC -
cA2 ) a2cA1 - a-2cA2 ) r K2
r3 ) k3 cA2 -
) )
cA3cB ) a3cA2 - a-3cA3 ) r K3
r4 ) k4 cA3cB r)
)
cDcE ) a4cA3 - a-4 ) r K4
a1a2a3a4 - a-1a-2a-3a-4 a-1a-2a-3 + a-1a-2a4 + a-1a3a4 + a2a3a4
Application to Claisen Condensation.
a1 ) k1cAcB a2 ) k2cC a3 ) k3 Figure 9. Carboxyalkylation of cellulose. Concentration profiles of substituted hydroxyl groups (OH-2, OH-3, and OH-6) during the reaction inside the porous particle.
a4 ) k4cB a1a2a3a4 ) k1k2k3k4cAcB2cC cE a-1 ) k1 K1 a-2 )
k2 K2
cB a-3 ) k3 K3 cDcE a-4 ) k4 K4
Figure 10. Carboxyalkylation of cellulose. Predicted distribution of mono- (2,3,6) (I), di- (23, 26, 36) (II), and trisubstituted (236) (III) anhydroglucose units during the reaction.
a-1a-2a-3a-4 ) condensation (shrinking particle) and carboxyalkylation of cellulose. The case studies illustrated how the models can be used for prediction of the reactor performance. Modeling of the liquid-solid systems starting from the basic principles is necessary for a rational process scaleup. Acknowledgment This work is a part of the activities at the Åbo Akademi Process Chemistry Group within the Finnish Centre of Excellence Program (2000-2005) by the Academy of Finland. The authors are grateful to Prof. B. Holmbom, Dr. Antti Vuori, Dr. Elias Suokas, Dr. Rainer Sjo¨holm, Mr. Paavo Mansikkama¨ki, Ms. O. Gro¨nfors, and Ms. K. Kaljula for their participation in fruitful discussions and experimental work.
k1k2k3k4 2 c c c K1K2K3K4 E B C
Notation A ) liquid-solid interfacial area AP ) outer particle surface area a ) shape factor aj, a-j, al ) factors in rate expressions, eq 2, and the appendix c ) concentration De ) effective diffusion coefficient DS ) degree of substitution Ea ) activation energy K ) equilibrium parameter k ) rate constant M ) molar mass m ) mass N ) flux n ) amount of substance n′ ) molar flow rate Q ) objective function in regression
2524
Ind. Eng. Chem. Res., Vol. 42, No. 12, 2003
R ) particle radius RG ) gas constant r ) radial coordinate ri ) generation rate of component j rj ) rate of reaction step j t ) time V ) volume V′ ) volumetric flow rate x ) dimensionless particle radius y ) dependent variable in regression R ) reactivity ratio in carboxyalkylation P ) particle porosity ν ) stoichiometric coefficient ω ) weight factor in regression Subscripts and Superscripts i ) component j ) reaction k, l, m ) arbitrary component, eq 2 L, l ) liquid phase P ) particle S, s ) solid SR ) total number of reactions t ) time 0 ) initial, inlet, or average quantity Abbreviations A, B, P, ... ) analytically detectable chemical components Ai (i ) 1, 2, ...) ) reaction intermediates CEC, CMC ) carboxyethyl- and carboxymethylcellulose GC ) gas chromatography HA ) monochlorocarboxylic acid HPLC ) high-performance liquid chromatography
Literature Cited (1) Levenspiel, O. Chemical reaction engineering; John Wiley & Sons: New York, 1972 (2nd ed.), 1999 (3rd ed.). (2) Sohn, H. Y.; Szekely, J. A structural model for gas-solid reactions with a moving boundary. Chem. Eng. Sci. 1972, 27, 763778. (3) Baldyga, J.; Bourne, J. R. Turbulent mixing and chemical reactions; John Wiley & Sons: New York, 1999. (4) Haario, H.; Kalachev, L.; Salmi, T.; Lehtonen, J. Asymptotic analysis of chemical reactions. Chem. Eng. Sci. 1999, 54, 11311143. (5) Marquardt, D. W. An algorithm for least squares estimation on nonlinear parameters. SIAM J. 1963, 11, 431-441.
(6) Green, J.; Margerison, D. Statistical treatment of experimental data; Elsevier: Amsterdam, The Netherlands, 1978; Chapter 5. (7) Henrici, P. Discrete variable methods in ordinary differential equations; John Wiley & Sons: New York, 1962. (8) Hindmarsh, A. ODEPACKsA systematized collection of ODE-solvers. In Scientific computing; Stepleman, R., et al., Eds.; IMACS/North-Holland: Amsterdam, The Netherlands, 1983; pp 55-64. (9) Buzzi Ferraris, G.; Manca, D. BzzOde: A new C++ class for the solution of stiff and non-stiff ordinary differential equation systems. Comput. Chem. Eng. 1998, 22, 1595-1621. (10) Caillaud, J. B.; Padmanabhan, L. An improved semiimplicit Runge-Kutta method for stiff systems. Chem. Eng. J. 1971, 2, 227-232. (11) Michelsen, M. L. An efficient general purpose method for the integration of stiff ordinary differential equations. AIChE J. 1976, 22, 594-597. (12) Kaps, P.; Wanner, G. A study of Rosenbrock-type methods of high order. Numer. Math. 1981, 38, 279-298. (13) Villadsen, J.; Michelsen, M. L. Solution of differential equation models by polynomial approximation; Prentice Hall: Englewood Cliffs, NJ, 1978. (14) Schiesser, W. E. The numerical method of lines; Academic Press: San Diego, 1991. (15) Romanainen, J. J.; Salmi, T. Numerical strategies in solving gas-liquid reactor models, Part III. Steady-state bubble columns. Comput. Chem. Eng. 1995, 19, 139-154. (16) Wa¨rnå, J.; Salmi, T. Dynamic modelling of catalytic threephase reactors. Comput. Chem. Eng. 1996, 20, 39-47. (17) Vajda, S.; Valko´, P. REPROCHEsRegression program for chemical engineers; Manual, European committee for computers in chemical engineering education, Budapest, 1985. (18) Haario, H. MODEST User’s Guide; Profmath Oy: Helsinki, Finland, 1994. (19) Wilson, K. A modified method for determination of active agent and degree of substitution in carboxymethyl cellulose (CMC). Sven. Papperstidn. 1960, 63, 714-715. (20) Flory, P. J. Kinetics of polyesterification: a study of the effects of molecular weight and viscosity on reaction rate. J. Am. Chem. Soc. 1939, 61, 3334-3340. (21) Niemela¨, K.; Sjo¨stro¨m, E. Characterization of hardwoodderived carboxymethylcellulose by gas-liquid chromatography and mass spectrometry. Polym. Commun. 1989, 31, 254-256. (22) Salmi, T.; Valtakari, D.; Paatero, E.; Holbom, B.; Sjo¨holm, R. Kinetic study of carboxymethylation of cellulose. Ind. Eng. Chem. Res. 1994, 33, 1454-1459.
Received for review August 10, 2002 Revised manuscript received December 3, 2002 Accepted December 5, 2002 IE0206296
