856
Ind. Eng. Chem. Res. 1996, 35, 856-863
Kinetic Model for Kraft Pulping of Hardwood Reinaldo Giudici* and Song Won Park LSCP-Process Control and Simulation Laboratory, Department of Chemical Engineering, University of Sa˜ o Paulo, Polytechnic School, P.O. Box 61548, 05424-970 Sa˜ o Paulo, Brazil
A comprehensive model for the kraft pulping kinetics of Eucalyptus saligna hardwood is presented. Kinetic parameters were estimated by fitting the model to available experimental data taken from the literature over a range of process variables. The model takes into account the effect of hydroxide and sulfide concentration in the liquor as well as the temperature-time history of the cooking. Model predictions were successfully compared with an independent set of bench-scale plant data for lignin and carbohydrate dissolution. The model is able to predict quite well the trends of the process variables. Introduction The kraft wood pulping is a process of lignin removal by alkali and sulfide at high temperature. Still today it is the dominant pulping process worldwide. Therefore, there is a strong economic incentive for better knowledge of the process. The accurate design and control of the pulping process require a reliable model that accounts for both the pulping kinetics and the transport limitations in the wood chips, as well as the flow characteristics in the digester. Further difficulties in the analysis of the pulping processes arise from the complex nature of the wood. The present work focuses on the kinetics of pulping. Reviews on the subject were presented by Wilder and Dalesky (1964) and Kleppe (1970). Chemical aspects of pulping were reviewed by Gierer (1980). More recently, Park and Pinto (1990) analyzed the evolution from empirical to mechanistic models for pulping kinetics. An early well-known attempt to model the kraft pulping kinetics was made by Vroom (1957). He combined the effect of cooking temperature-time history with an Arrhenius type expression for the reaction rate temperature dependence into one variable, named the “H-factor”. The H-factor is still used in many control schemes for the pulping process. However, the correlation based on the H-factor does not take into account the effects of wood composition and concentration of pulping agents. Subsequently, several different semiempirical kinetic models have been proposed in the literature. However, as pointed out by Saltin (1992), all of them have the same general form:
-dW/dt ) (k1[OH-]a + k2[OH-]b[HS-]c)(W - W∞) (1) where W stands for the concentration of a wood component (e.g., lignin, carbohydrate) and W∞ is the wood component concentration after infinite time. In a number of works, the kinetics of the kraft pulping is divided into three periods, namely, initial, bulk, and residual. Equation 1 is then applied for the delignification with different parameters for each period. The transition points between the periods are expressed as a function of lignin content and are identified empiri* Author to whom correspondence should be addressed. E-mail:
[email protected]/
[email protected]. Fax: (5511)8132380.
0888-5885/96/2635-0856$12.00/0
cally using plots of carbohydrate content versus lignin content. The approximation of this plot by three intersecting straight lines considers the rate of carbohydrate removal to be proportional to the delignification rate in each period. It should be noted, however, that the transition points depend on the wood type and, even more significantly, they may vary with the pulping conditions for the same wood. This is evidence of a delignification rate dependence on the morphologic change of the lignin. Although this “stepwise” approach is largely used, it has some disadvantages. From the point of view of modeling, the change of the transition points with the wood type is not a concern, because the kinetic parameters are also wood-dependent. However, the change of transition points with the pulping conditions may represent a serious problem for the prediction capability of such a stepwise approach. An alternative, better approach tries to incorporate mechanistic information into the model. An example is the model presented by Saltin (1992), that considers the lignin as divided into high reactivity or fast lignin (with β ethers) and low reactivity or slow lignin (with γ ethers). Carbohydrate removal can be expressed in terms of the major components (cellulose, glucomannan, xylan), each one following equations similar to eq 1. An even more mechanistic approach was presented by Burazin and McDonough (1988) for softwood. They consider that both carbohydrate components and lignin undergo the transformations by a detailed mechanism involving a number of reaction steps in series and/or in parallel. The final reaction network for lignin involves bulk delignification, condensation, and residual delignification. The best model for carbohydrates involves peeling, stopping, and cleavage reactions, with different parameters for cellulose, glucomannan, and xylan. Three pulping agents were considered, namely, NaOH, NaSH, and anthraquinone. Compared to the literature on softwood pulping kinetics, there is a lack of information about an accurate hardwood pulping kinetics. The objective of the present work is to present a comprehensive kinetic model for kraft pulping of eucalyptus and to obtain its kinetic parameters by fitting the model to experimental data. Kinetic Model The mechanistic model proposed for the kraft pulping kinetics by Burazin and McDonough (1988) was used as a starting point in the development of our model for eucalyptus, with the necessary adaptations. © 1996 American Chemical Society
Ind. Eng. Chem. Res., Vol. 35, No. 3, 1996 857
According to Burazin and McDonough (1988), lignin is a highly branched three-dimensional heterogeneous polymer consisting of three major precursors joined together by several types of bonds. Lignin undergoes degradation reactions, which produce dissolved lignin, while condensation (cross-linking) reactions produce residual lignin, which is less reactive than native lignin but otherwise very difficult to distinguish from native lignin. Carbohydrates are present as cellulose, a high molecular weight linear polymer of glucose, and hemicellulose, low molecular weight polymers of several sugars. The peeling reaction depolymerizes a carbohydrate chain one sugar unit at a time from the reducing (hemiacetal) end of the chain. The stopping reaction oxidizes the reducing end to a metasaccharinic end which does not peel. The cleavage reaction breaks a carbohydrate chain in two, creating a new reducing end which subsequently may peel. The mechanism can be represented by: kb(NaOH,NaSH)
kc
dLR/dt ) kcLD - krLR
(4)
dLD/dt ) kbLN - kcLD
(5)
dLDI/dt ) krLR
(6)
kb ) [NaOH][NaSH]1/2T1/2 exp(k1 - E1Tm)
(7)
kc ) T1/2 exp(k2 - E2Tm)
(8)
kr ) [NaOH]T1/2 exp(k3 - E3Tm)
(9)
and for the carbohydrate reactions
C ) CN + CX
(10)
dCN/dt ) kclCX - (kpl + kst)CN
(11)
dCX/dt ) kstCN - kclCX
(12)
dCD/dt ) kplCN
(13)
kcl ) [NaOH]T1/2 exp(k4 - E4Tm)
(14)
kr(NaOH)
LN 98 LD 98 LR 98 LDI kcl(NaOH)
kpl(NaOH,NaSH)
CX {\ } CN 98 CD k (NaOH)
kpl ) [NaOH][NaSH]1/2T1/2 exp(k5 - E5Tm)
st
In the kinetic model, it is assumed that total lignin (L) is present in two forms, native (LN) and residual lignin (LR). Bulk delignification reactions convert native lignin to dissolved lignin (LD). Condensation reactions convert LD to residual lignin, and residual delignification converts LR to inert dissolved lignin (LDI). For carbohydrates, the model considers that total carbohydrates (C) consist of two fractions, native carbohydrates (CN) and oxidized carbohydrates (CX). Stopping reactions convert CN to CX, and cleavage converts CX back to CN. The only distinction between CN and CX is that CX does not peel. Peeling converts CN to dissolved carbohydrates CD. The main differences between our model and that of Burazin-McDonough’s paper are as follows: (a) we did not consider the effect of anthraquinone, since this pulping agent was not present in the experiments used to fit the model; (b) we model both the bulk delignification and the carbohydrate peeling rates as first-order in NaOH and half-order in NaSH, while they assumed two parallel reactions for each of these steps. We did not attempt to regress the reaction orders as additional parameters because the initial concentrations of the pulping agents were not varied in the available experimental data. Then a number of different values and combinations for these reaction orders were tried, and a satisfactory result was found for the values aforementioned. Regarding this point, our model is analogous to those of Teder and Olm (1981) and Kubo et al. (1983), in contrast to the expressions used by Lemon and Teder (1973) and Saltin (1992). According to these considerations, the rate and mass balance equations for the species involved in the delignification are
L ) LN + LR
(2)
dLN/dt ) -kbLN
(3)
kst ) [NaOH]T1/2 exp(k6 - E6Tm)
(15) (16)
where
Tm )
1 (T1 - 417 )
(17)
In eq 17, T is the absolute temperature (K) and 417 K is an average experimental temperature. This equation represents a typical variable transformation useful to reduce interaction between the kinetic parameters ki and Ei during the parameter estimation procedure (Himmelblau, 1970; Kittrell, 1970; Bates and Watts, 1988). Note that eqs 7-9 and 14-16 were written in terms of logarithms of preexponential parameters rather than preexponential rate constants. This parameter transformation is recommended in the literature (Watts, 1994) because it behaves more linearly regarding the parameter estimation. The temperature variation during a pulping run is given by a ramp followed by a constant temperature period:
T(t) )
{
T0 + at Tp
t e (Tp - T0)/a t g (Tp - T0)/a
(18)
where a is the heating rate (°C/min) and Tp is the “plateau” temperature. The starting temperature was T0 ) 30 °C. This temperature variation is accounted for in the estimation procedure. Experimental Data For the parameter estimation of our eucalyptus pulping model, the experimental data taken from Bugajer’s thesis (Bugajer, 1984) were used. Bugajer’s experiments were carried out with Eucalyptus saligna in a bench-scale 20-L digester with automatic temperature
858
Ind. Eng. Chem. Res., Vol. 35, No. 3, 1996
Table 1. Wood Characterization and Experimental Cooking Conditions Used by Bugajer (1984) Wood Characteristics soluble in cold water soluble in hot water soluble in NaOH (1%) solution soluble in ethanol-benzene soluble in ether cellulose Cross-Bevan pentosanas lignin Klason lignin soluble in acid mineral residue at 600 °C
2.42% 3.54% 14.73% 1.12% 0.11% 58.59% 17.88% 25.32% 3.45% 0.36%
Granulometric Analysis (Circular Hole Screening) sieve diameter (cm) 2.86 2.22 1.59 0.95 0.48 fines active alkali sulfidity effective alkali liquor/wood ratio a
retention (% wt) 12.5 17.8a 26.5a 32.7a 8.4 2.1 Reaction Conditions 105.3 g of Na2O/L 25% 14% (Na2O/dry wood) 4:1
Fractions used in the pulping runs sum up to 77%.
control. The runs were done using different cooking times (25, 50, 75, and 100 min), temperatures (80, 140, 150, 155, 160, 165, and 170 °C), and heating rates (1.85-2.28 °C/min) but keeping constant the initial active alkali (14% Na2O) and sulfidity (25%) in the liquor, as well as the liquor/wood ratio (4:1). Table 1 shows the wood characterization and the pulping conditions used in the experiments. The chip thickness, that is the critical dimension for internal diffusion of the components in the wood chip, was about 3 mm. This value is within the limit in which the chip thickness does not affect the results (Wilder and Daleski, 1965; Pekkala, 1983). No mechanical pretreatment of the chips was reported in Bugajer’s thesis. However, as pointed out by Gustafson et al. (1983), the liquor penetration can be considered complete for this thickness about the time the temperature reaches values high enough so that the pulping reactions occur. In addition, the mass-transfer resistance from the bulk phase to the chips is negligible for laboratory digesters (Gustafson et al., 1983). Therefore, the system can be considered not mass transfer limited. In the present work, the experimental data from Bugajer (1984) were split into two sets. The first one, with 12 runs and 58 data points, was used in the parameter estimation procedure. The second group, with 2 runs and 13 points, was used to check the model predictions with an independent data set. Table 2 reproduces the experimental data of Bugajer’s experiments. Estimation of Alkali and Sulfide Consumption Parameters Bugajer (1984) did not use a high liquor/wood ratio as is usually used in laboratory pulping kinetic studies. Instead, she used the ratio 4:1, which is closer to the typical ratio used in the pulp and paper industry. Due to the low excess of alkali and sulfide, the NaOH and NaSH concentrations may vary significantly during the pulping run. Therefore, it is necessary to account for the consumption of alkali and sulfide in the model.
Table 2. Experimental Data Taken from Bugajer (1984)
run
Tp (°C)
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 1-v1 2-v1 3-v1 4-v1 5-v1 6-v1 1-v3 2-v3 3-v3 4-v3 5-v3 6-v3 7-v3
80.3 80.0 79.8 78.8 81.7 80.9 81.5 81.9 80.5 141.2 142.7 142.9 140.5 141.0 142.2 140.2 141.3 141.0 142.2 149.7 149.0 149.8 149.1 155.7 154.5 154.6 155.0 153.2 160.1 159.5 159.3 160.1 160.2 160.8 161.3 159.1 160.0 159.0 165.7 165.0 165.0 164.8 164.9 165.0 164.8 164.8 164.5 164.3 169.0 169.8 168.6 170.9 169.1 170.0 169.0 170.0 170.3 168.0 80.0 140.3 159.1 158.1 159.3 159.4 80.0 140.0 168.9 169.1 168.5 168.1 168.8
a t (°C/min) (min) 1.80 1.88 1.90 1.94 2.06 2.04 2.03 2.07 2.00 2.28 2.20 2.05 1.84 2.22 2.23 2.20 2.23 2.22 2.04 2.18 2.16 2.17 2.17 2.28 2.08 2.18 2.19 2.16 2.20 2.02 2.16 2.17 2.17 2.11 2.19 2.15 2.17 2.15 2.07 2.08 2.04 2.07 2.08 2.08 2.04 2.06 2.07 1.92 1.85 1.86 1.85 1.88 1.85 1.87 1.85 1.86 1.87 1.84 1.92 2.00 1.99 1.94 1.99 2.05 2.18 2.20 2.11 2.13 2.13 2.12 2.14
0 53 77 101 125 25 50 75 100 125 50 80 105 135 150 50 75 100 125 155 55 80 110 130 55 85 107 132 157 59 89 110 139 162 62 85 110 135 160 63 90 116 140 165 68 90 115 140 170 75 100 125 160 165 75 100 125 150 175 25 55 65 91 115 163 23 50 65 90 115 135 165
L (%) 28.77 24.1 23.2 22.6 23.4 22.9 24.8 23.5 22.4 22.4 15.8 8.97 7.92 4.52 3.62 15.5 11.1 6.79 4.52 3.90 13.9 6.51 2.78 2.33 13.5 4.15 2.48 1.48 1.02 9.50 2.18 1.74 1.50 1.85 11.0 2.41 1.24 1.04 0.82 5.43 1.59 1.11 0.98 1.39 4.60 1.64 0.85 0.97 1.02 2.12 1.42 0.996 0.879 0.964 2.36 1.38 1.08 0.833 0.66 24.0 20.1 12.3 4.86 2.42 1.66 26.5 19.5 9.06 1.92 1.28 1.33 1.26
[NaOH] [NaSH] (mol of (mol of C (%) Na2O/L) Na2S/L) 71.23 64.7 62.1 62.4 61.0 63.3 65.6 64.6 61.8 61.6 55.2 56.3 54.0 52.1 51.6 58.3 54.3 54.3 51.5 54.1 56.7 53.8 51.2 51.2 55.8 50.8 49.3 49.8 50.5 51.2 49.1 49.2 49.3 47.3 56.3 51.1 50.2 48.6 49.4 48.0 48.7 47.4 47.6 47.9 52.0 50.3 49.7 48.7 48.4 48.8 47.4 47.8 46.9 46.8 50.2 47.8 48.1 47.3 48.2 61.1 60.4 54.7 53.5 50.4 51.5 59.7 59.8 53.0 53.2 49.0 48.6 47.3
0.594 0.596 0.570 0.455 0.512
0.1053 0.1001 0.1042 0.0765 0.0877
0.318 0.207 0.162 0.161 0.099 0.276 0.143 0.101 0.075 0.236 0.158 0.096 0.049 0.055
0.0464 0.0279 0 0.0363 0.0092 0.0746 0.0188 0.0297 0.0267 0.0823 0.0381 0.0285 0.0153 0
0.181 0.061 0.054 0 0.016
0.0403 0.0291 0.0195 0.0073 0
0.086 0.061 0.051 0 0
0.0523 0.0174 0.088 0 0
0.077 0.043 0.026 0.017 0
0.0123 0.0217 0.0060 0 0
Due to the complex macromolecular nature of the lignin and carbohydrate structures as well as their
Ind. Eng. Chem. Res., Vol. 35, No. 3, 1996 859
of liquor, [NaSH] is in moles of Na2S per liter of liquor, and C and L are in percent of dry wood. These relationships are illustrated in Figures 1 and 2. The parameter 0.0140 in eq 19 means that on average 0.0140 mol/L of NaOH is consumed when the content of total carbohydrate decreases 1 unit. A similar meaning can be given to the parameter 0.0128, using lignin in place of carbohydrate. The parameter 0.164 in eq 20 represents some selectivity (or, better, reactivity) ratio between sulfide and hydroxide for the pulping reactions. In order to incorporate these expressions into the model, eqs 19 and 20 should be written in their differential form, with their respective initial conditions: Figure 1. Comparison between the experimental data and the predictions by eq 19 for NaOH concentration.
d[NaOH] dC dL ) 0.0140 + 0.0128 dt dt dt [NaOH]|at t)0 ) [NaOH]0 (21) d[NaSH] d[NaOH] ) 0.164 dt dt [NaSH]|at t)0 ) [NaSH]0 (22)
Figure 2. Comparison between the experimental data and the predictions by eq 20 for NaSH concentration.
fragments formed during the pulping, the determination of the stoichiometry of the reactions between NaOH/ NaSH and the wood components can be very difficult. Difficulties in considering the alkali and sulfide consumption in the model may be overcome by using “effective stoichiometric” coefficients that are determined from the experiments. These effective stoichiometric factors relate the amount of NaOH or NaSH consumed for a given amount of lignin and/or carbohydrate reacted. Obviously, they are not true stoichiometric coefficients in the sense of simple reactions. They are more like a yield factor, although yield is usually defined in terms of one product and one reactant. Nonetheless, this approach is practical and useful for taking into account the change of NaOH and NaSH concentrations. In order to avoid increasing the (already large) number of parameters in the estimation procedure, the stoichiometric coefficients for NaOH and NaSH consumption were estimated independently of the kinetic parameters. From the global variation of lignin, carbohydrate, NaOH, and NaSH during each run, one obtains the following relationships, whose coefficients were found by linear regression:
[NaOH] ) -0.67 + 0.0140C + 0.0128L
(19)
[NaSH] ) 0.00705 + 0.164[NaOH]
(20)
where [NaOH] is expressed in moles of Na2O per liter
Note that these equations couple the lignin model equations to the carbohydrate model equations, which otherwise (that is, in the case of large excess of NaOH and NaSH) would be decoupled. Expressions mathematically similar to eq 21 are found in the works of Gustafson et al. (1983) and Saltin (1992). The approach used by Gustafson et al. starts from fitting the carbohydrate content and the alkali concentration as linear functions of the lignin content for both the initial and bulk periods during a kraft cooking. This gives two linear equations which can be solved for the two unknowns, the coefficients in eq 21. This approach is sensitive to the slopes of the initial and bulk periods. In addition, carbohydrate and lignin consumptions are correlated, and both these facts together would explain the tenuous physical significance of the parameters arising from this procedure (for instance, Gustafson et al. obtained a negative coefficient for the lignin term). The approach taken here (and also used by Saltin) seems to be better since it does not require the evaluation of slopes, an operation that usually amplifies experimental errors. Some inaccuracy cannot be avoided, however, since the correlation between lignin and carbohydrate consumption is always present in conventional cooking experiments. Estimation of Kinetic Parameters The model is represented by ordinary differential equations (ODE):
dx(t)/dt ) f(x(t),t,β)
x(t0) ) x0
(23)
and measured variables y are related to the state variables x in the form:
y(t) ) h(t,x(t))
(24)
The maximization of the likelihood function provides a general formulation for the choice of an objective function in a parameter estimation problem. Depending on the error structure, different objective functions can be obtained, as discussed in the paper by Biegler et al.
860
Ind. Eng. Chem. Res., Vol. 35, No. 3, 1996
Table 3. Parameter Summary number of data points n ) 58 × (2 dependent variables) ) 112 number of parameters p ) 12 minimum residue Smin ) 217.5 variance s2 ) Smin/(n - p) ) 2.09 parameter
estimate
standard error
t-value
parameter
estimate
standard error
t-value
k1 k2 k3 k4 k5 k6
-2.676 -5.541 -4.490 -7.341 -3.514 -3.944
0.220 0.157 0.136 0.936 0.567 0.746
12.16 35.29 33.01 7.84 6.20 5.29
E1 E2 E3 E4 E5 E6
7686 3079 20293 26447 8217 6657
735 767 2343 8565 1258 2389
10.46 4.01 8.66 3.09 6.53 2.79
Parameter Correlation k1 E1 k2 E2 k3 E3 k4 E4 k5 E5 k6 E6
k1
E1
k2
E2
k3
E3
k4
E4
k5
E5
k6
E6
1 0.51 0.08 -0.16 -0.65 -0.52 -0.03 0.08 0.21 0.18 0.13 0.06
1 0.08 0.55 -0.09 -0.16 0.04 0.03 0.11 0.15 0.12 0.11
1 0.50 -0.18 -0.08 -0.07 0.22 -0.22 -0.21 -0.18 -0.10
1 0.27 0.18 0.025 0.08 -0.14 -0.09 -0.06 0.006
1 0.21 0.004 -0.10 -0.08 -0.06 -0.07 -0.04
1 0.09 -0.02 0.004 0.05 0.08 0.11
1 -0.85 0.49 0.68 0.84 0.81
1 -0.77 -0.78 -0.81 -0.55
1 0.89 0.79 0.37
1 0.94 0.73
1 0.85
1
(1986). In the present work we used the least squares criterion: n
min S(β) ) β
[y(ti,β) - y(ti)]2 ∑ i)1
(25)
Statistically more rigorous criteria would have been to use either the multiresponse estimation procedure of minimizing the determinant of the sum of squares and cross-products of the observations (Hunter, 1967) or a multifunctional error-in-variables method (Reilly and Patino-Leal, 1981). However, for the present case the simple least-squares criterion was adequate for our purposes. For the model equations 2-18, one has:
vector of state variables x ) [LN,LR,LD,CN,CX]T
(26)
vector of the initial conditions x0 ) [L0,0,0,C0,0]T
(27)
vector of measured variables y ) [L,C]T
(28)
β ) [k1,E1,k2,E2,k3,E3,k4,E4,k5,E5,k6,E6]T
(29)
vector of parameters
In general, methods for solving the estimation parameters in ODE’s may include either the direct integral method or the so-called indirect integral method. The direct method fits the equivalent integral equations, solving the integrals of the experimental points by a numerical quadrature scheme (Himmelblau et al., 1967). Although Vadja et al. (1986) have shown that this method is robust and numerically efficient, it has some disadvantages from the point of view of the statistics, leading to biased estimates of the parameters.
In addition, this technique requires a dense data set in order to evaluate the quadrature; therefore, it is not applicable to the present case. The indirect method fits the solution of the ODE’s, obtained numerically by an adequate step-by-step technique, leading to unbiased estimates, and is applicable to nondense data sets. This approach was used for the parameter estimation in the present work. The system of ODE’s was solved by a standard variable-step fourthorder Runge-Kutta-Gill method. The minimization of criterion (25) was done using Marquardt’s method (Marquardt, 1963) with numerical evaluation of the Jacobian matrix. Results and Discussion Table 3 shows the estimated kinetic parameters, as well as the correlation matrix of the parameters. The kinetic parameters for carbohydrate reactions present higher cross-correlation than those for delignification. Although some parameters are better estimated than others, all estimated parameters are significant at the 95% confidence level. The parity plot in Figure 3 compares experimental and predicted values for total lignin content L and for total carbohydrate content C. Good agreement was achieved. Most deviations between experimental and calculated values are distributed within 1.5% for lignin and within 2.0% for carbohydrate. Also, the deviations present no trends with the dependent variables. Figure 4 shows the distribution of the errors as a function of time (the independent variable). Higher deviations occur for a few points taken at shorter times. Two factors might be used to explain this behavior. First, these points correspond to experimental samples taken when the concentration varies faster, mostly during the temperature ramping. Such points are presumably much more affected by the errors due to the delay in stopping the reactions. Second, the generic term carbohydrate embodies the cellulose as well as the hemicellulose, a component also present in the wood, but hemicellulose undergoes faster dissolution than cellulose at the beginning of the cooking. Nevertheless,
Ind. Eng. Chem. Res., Vol. 35, No. 3, 1996 861
Figure 6. Comparison between model predictions (s) and experimental values of lignin (0) and carbohydrate (4) for run v1. Figure 3. Lignin (0) and carbohydrate (4) contents as calculated by the model versus experimental measurements.
Figure 4. Deviation between calculated and experimental values for lignin (0) and for carbohydrate (4) versus time.
Figure 5. Carbohydrate contents versus lignin contents for runs 1-58: experimental values (0) and model prediction (+).
such deviations are not a concern because from the engineering’s point of view the predictions of industrial interest are those at the end of the cooking. Figure 5 shows the classical relation between carbohydrate and lignin during the pulping. The points follow the well-known trend of three periods: initial, bulk, and residual. At first sight, the fluctuations in the experimental points might be ascribed to experimental errors. However, this is not the case for the prediction points, since they were calculated by the fitted model. This result clearly shows that the model is sensitive enough to respond to the small variations
in heating rate, temperature, and cooking time that occurred from one run to another. In addition, for a given run, the mechanistic model is able to reproduce this behavior smoothly, without imposing “a priori” the three periods, as can be seen in Figures 7 and 9. Delignification and carbohydrate reactions are parallel reactions regarding the consumption of NaOH and NaSH. Incidentally, eq 21 shows that the coefficients for lignin and for carbohydrate are of the same magnitude. It is interesting to note that the lignin model is able to predict an increase of lignin content with time under certain conditions, namely, at low alkali concentration and low temperature. This is due to the fact that the condensation reaction does not depend on NaOH concentration and has lower activation energy than the other lignin reactions. This behavior is experimentally observed in Bugajer’s data. This is also corroborated by Hartler (1978), according to whom the condensation reactions take place when the cooking liquor lignin concentration is high and the alkali concentration low. A low value for the activation energy of the lignin condensation reaction has been found. This is an indication that this reaction should be more diffusioncontrolled than the others. For a different wood (pine), Burazin and McDonough (1988) also observed a low activation energy for the condensation, compared to the other delignification reactions. Those authors decided to fix the activation energy for the condensation reaction to a value expected for diffusion-controlled reactions. From an engineering point of view, a simple fit of kinetic expressions valid over only one of the three periods is enough for several applications. However, this piecewise approach overlooks the fundamentals of the pulping phenomena. Only a mechanistic kinetic model that considers the changes in the reactive species during the pulping may be valid over the whole process and reproduce the transitions between periods smoothly. Finally, the detailed mechanistic kinetic model developed by Burazin and McDonough for pine was able to represent eucalyptus pulping with only minor modifications. These findings suggest that this comprehensive mechanism may represent an improvement in the understanding of the kraft pulping kinetics. Checking the Model Predictions In order to provide a fair check of the quality of the model predictions, we separated a set of two runs of Bugajer (1984), namely, runs v1 and v3, not used in the parameter estimation step. Figures 6-9 illustrate that the simulation results follow the trends of the experi-
862
Ind. Eng. Chem. Res., Vol. 35, No. 3, 1996
Work is now underway to combine the present kinetics with a model for the digester, including the model for diffusion-reaction phenomena inside the wood chips as previously developed by one of the authors (Park, 1988). In addition, the model predictions will be compared with industrial digester data. These applications will be reported in a forthcoming paper of this sequence. Acknowledgment
Figure 7. Carbohydrate contents versus lignin contents for run v1: (s) model; (0) experimental values.
Figure 8. Comparison between model predictions (s) and experimental values of lignin (0) and carbohydrate (4) for run v3.
Figure 9. Carbohydrate contents versus lignin contents for run v3: (s) model; (0) experimental values.
mental data well. This provides an adequate validation test for the kinetic model. Summary The methodology of kinetic parameter estimation in an ordinary differential equation model has been applied to the kraft pulping of eucalyptus. A comprehensive mechanistic kinetic model has been proposed based on the previous one presented by Burazin and McDonough (1988), with some minor modifications. Lignin reactions included bulk delignification, lignin condensation, and residual delignification. The carbohydrate model considers peeling, stopping, and cleavage reactions. Effects of alkali and sulfide, as well as of temperature, have been taken into account. Kinetic parameters have been estimated from fitting the model to available literature data. All estimates have been found significant at 95% confidence. The model has shown to be able to follow quite well the trends of the main process variables.
The authors thank Prof. Frank Quina for revising the manuscript. Nomenclature a ) heating rate (K/s) a, b, c ) parameters in eq 1 C ) total carbohydrate concentration (% wood) CN ) native lignin concentration (% wood) CD ) dissolved carbohydrate concentration (% wood) CX ) oxidized carbohydrate concentration (% wood) Ei ) activation energy divided by gas constant (K) kb ) rate constant for bulk delignification reaction kc ) rate constant for condensation reaction kr ) rate constant for residual delignification kcl ) rate constant for cleavage reaction kpl ) rate constant for peeling reaction kst ) rate constant for stopping reaction ki ) kinetic parameter L ) total lignin concentration (% wood) LD ) dissolved lignin concentration (% wood) LDI ) inert dissolved lignin concentration (% wood) LN ) native lignin concentration (% wood) LR ) residual lignin concentration (% wood) [NaOH] ) hydroxide concentration (mol of Na2O/L) [NaSH] ) hydrosulfide concentration (mol of Na2S/L) n ) number of experimental points S ) least-squares criterion t ) time (s) T ) temperature (K) Tm ) variable defined in eq 17 T0 ) initial cooking temperature (K) Tp ) final cooking temperature (K) x ) vector of state variables y ) vector of measured variables W ) wood component concentration W∞ ) wood component concentration after infinite time β ) vector of parameters
Literature Cited Bates, D. M.; Watts, D. G. Nonlinear regression analysis and its applications; Wiley: New York, 1988. Biegler, L. T.; Damiano, J. J.; Blau, G. E. Nonlinear parameter estimation: a case study comparison, AIChE J. 1986, 32 (1), 29-45. Bugajer, S. Kinetics of delignification reactions in the eucalyptus pulping process (in Portuguese). Ph.D. Thesis, Polytechnic School, University of Sa˜o Paulo, Sa˜o Paulo, Brazil, 1984. Burazin, M. A.; McDonough, J. Building a mechanistic model of kraft-anthraquinone pulping kinetics. Tappi J. 1988, March, 165-169. Froment, G. F.; Hosten, L. H. Catalytic kinetics: modelling. In Catalysis Science and Technology; Anderson, J. R., Boudart, M., Eds.; Springer-Verlag: Berlin, 1981; Vol. 2, Chapter 3, pp 97170. Froment, G. F.; Bischoff, K. B. Chemical Reactor Analysis and Design, 2nd ed., Wiley: New York, 1990. Gierer, J. Chemical aspects of kraft pulping. Wood Sci. Technol. 1980, 14 (4), 241-266. Gustafson, R. R.; Sleicher, C. A.; McKean, W. T.; Finlayson, B. A. Theoretical model of the kraft pulping process. Ind. Eng. Chem. Process Des. Dev. 1983, 22, 87-96.
Ind. Eng. Chem. Res., Vol. 35, No. 3, 1996 863 Hartler, N. Sorption cooking. Sven. Papperstidn. 1978, 81 (14), 457-463. Himmelblau, D. M. Process Analysis by Statistical Methods; Wiley: New York, 1970. Himmelblau, D. M.; Jones, C. R.; Bischoff, K. B. Determination of rate constants for complex kinetic models. Ind. Eng. Chem. Fundam. 1967, 6, 539. Hunter, W. G. Estimation of unknown constants from multiresponse data. Ind. Eng. Chem. Fundam. 1967, 6 (3), 461-463. Kittrell, J. R. Mathematical modeling of chemical reactions. Adv. Chem. Eng. 1970, 8, 97-183. Kleppe, P. J. Kraft pulpingsfeature review. Tappi J. 1970, 53 (1), 35-47. Kubo, M.; Yoshida, H.; Tamao, M.; Ueno, T. A kinetic model of delignification in kraft pulping. Proceedings, 1983 International Symposium on Wood and Pulping Chemistry, Tsukuba, Japan; Japan Tappi: Tsukuba, Japan, 1983, Vol. 2, pp 130-135. Lemon, S.; Teder, A. Kinetics of delignification in kraft pulping. I. Bulk delignification of pine. Sven. Papperstidn. 1973, 76 (11), 407-414. Marquardt, D. W. An algorithm for least squares estimation of nonlinear parameters. J. Soc. Ind. Appl. Math. 1963, 11 (2), 431-441. Park, S. W. An heterogeneous model for kraft pulping of eucalyptus (in Portuguese). M. Sc. Thesis, Polytechnic School, University of Sa˜o Paulo, Sa˜o Paulo, Brazil, 1988. Park, S. W.; Pinto, J. M. Kinetics of kraft pulping: the evolution from empirical to mechanistic models (in Portuguese). Proceedings of the 23rd Pulp and Paper Brazilian Annual Meeting, Sa˜o Paulo, Nov 5-9, 1990; pp 69-91.
Pekkala, O. Some features of residual delignification during kraft pulping of Scots pine. Pap. Puu 1983, 65 (4), 251-263. Reilly, P. M.; Patino-Leal, H. A Bayesian study of the error-invariables model. Technometrics 1981, 23 (3), 221-231. Saltin. J. F. A predictive dynamic model for continuous digesters. Tappi Proc. Pulp. Conf. 1992, 261-268. Teder, A; Olm, L. Extended delignification by combination of modified kraft pulping and oxygen bleaching. Pap. Puu 1981, 63 (4a), 315-326. Vadja, S.; Valko´, P.; Yermakova, A. A direct-indirect procedure for estimation of kinetic parameters. Comput. Chem. Eng. 1986, 10 (1), 49-58 Vroom, K. E. The H-factor: a means of expressing cooking times and temperatures as a single variable. Pulp Pap. Can. 1957, 58 (Convention Issues), 228-231. Watts, D. G. Estimating parameters in nonlinear rate equations. Can. J. Chem. Eng. 1994, 72, 701-710. Wilder, H. D.; Daleski, E. J. Kraft pulping kinetics. Tappi J. 1964, 47 (5), 270-275. Wilder, H. D.; Daleski, E. J. Delignification rate studies. Tappi J. 1965, 48 (5), 293-297.
Received for review June 8, 1995 Revised manuscript received November 27, 1995 Accepted December 6, 1995X IE950341Z X Abstract published in Advance ACS Abstracts, February 1, 1996.
