Estimating pure diffusion contributions in alkaline pulping processes

Vicente Costanza ^† and Pedro Costanza ^‡

^† Grupo Tecnología de la Madera, INTEC (CONICET-UNL), Güemes 3450, 3000 Santa Fe, Argentina.
tsinoli@ceride.gov.ar
^‡ Facultad de Ingeniería Química, Universidad Nacional del Litoral, Santiago del Estero 2654, 3000 Santa Fe, Argentina
pedroc@intec.unl.edu.ar

Abstract — A model that predicts isothermal alkali diffusion and reaction with acetyl groups in moist wood chips was derived and approximated. System parameters were estimated from unsteady-state experimental data. Simulation results reinforce the idea that the diffusion effect is not fully exploited in pulping processes. Traditionally, digestion is conducted at high temperature, where delignification reaction kinetics is enhanced and the reaction effect is predominant. This approach is being reviewed by modern industry since energy and environmental savings associated with low temperature operation might compensate for high-yield productivity. The concentration of alkali at the center of the chip is a measure of the completeness of wood deacetylation, which translates into the aptitude of the final product for pulping purposes. This concentration is predicted here from the solution to a pair of coupled ODE's. Since alternatives combining both low and high-temperature processes are being studied, the results in this paper provide basic data for optimization analysis.

Keywords — reaction-diffusion, wood pulping, alkaline pulping, impregnation, transport phenomena.

I. INTRODUCTION

A comprehensive model for diffusion and reaction of alkali in wood chips is important for optimizing usual pulping processes in the paper industry. Alkaline treatment is increasingly replacing sulfite chemimechanical pulping due to the enforcement of environmental regulations in many countries. The alkali generates carboxylic groups in contact with cellulose and hemicellulose (an effect which is not desired for pulping purposes), and deacetylates the wood, both mechanisms resulting in a softening and swelling of the solid material needed for ulterior defibration. According to previous work (Zanuttini and Marzocchi, 1997) poplar wood (Populus deltoides f. angulata and carolinensis) swelling caused by alkaline treatment is related to by alkaline treatment is related to deacetylation rather than to carboxylic group content.
Experiments are time-consuming and expensive, and industrial conditions may be difficult to reproduce in the laboratory. An accurate theoretical model, on the other hand, provides a quick and inexpensive means for studying, understanding, and monitoring the pulping process and mechanisms. Our model's usefulness for optimization purposes will depend on its accuracy and robustness when temperature and bulk concentration are allowed to vary within a wide range, since these are the main operation conditions that can be manipulated.
Traditionally, pure alkaline treatment is conducted at temperatures well above 100 °C. High temperature results in an overall diffusion-controlled process with flat reaction profiles in agreement with the unreacted core model; for short reaction times, this leads to mostly unreacted chip cores. Ideally, reaction should occur as soon as the alkali reaches acetyl groups. It is therefore desired that the penetration of the alkali into the wood be reaction-driven, which happens only at high reaction rates. Perhaps one of the reasons why some researchers (Jimenez et al., 1989a, 1989b; Zanuttini and Marzocchi, 1997; etc.) adopted the flat front model was the familiarity with pressure treatments in wood impregnation, where diffusion is practically ignored (Costanza and Miyara, 1988, 2000). At low temperatures alkali diffusion becomes substantial, the process is no longer diffusion- controlled, and there is a reaction gradient inside the chip, so reaction occurs in all points of the material, leaving no unreacted core. The results obtained in this work suggest that the pulp and paper industries may be misusing the low-temperature diffusion effect. In other words, taking diffusion into account and allowing for a pretreatment of the chips would result in overall energy savings. The high temperatures used in common onestep reaction processes would no longer be required. Moreover, high temperatures cause important material loss together with high alkali consumption. The combination of alkaline pretreatment and sulfite final diges tion has recently been proposed as a compromise (Zanuttini et al., 1998).
Previous experimental work on the pure diffusion mechanism deserves criticism. Lönnberg and Robertsen, 1992, did not distinguish between diffusion alone and diffusion combined with chemical reaction, so they obtained not really pure but combined coefficients. Talton and Cornell, (1987), measured diffusion in the reverse direction, i.e. from saturated wood to dilute solution, without considering eventual distortions due to hystheresis effects. In our work an experimental procedure is set up to obtain reliable diffusion coefficients for alkali solutions in poplar wood chips in the radial direction. Extensions to other types of impregnating solutions, wood species, and cutting angles are immediate.
It should also be noticed that the model described here would not respond accurately for thicker wood samples due to wood inhomogeneity, moisture content, air entrapped in the vessels, etc. (Kelso et al., 1963; Siau, 1984). Even under our simplifying assumptions a differential mass balance had to be solved (Slattery, 1972), and some parameters had to be precisely defined and determined. One of the parameters considered was the pure diffusion coefficient. Also, an alternative expression for the reaction rate was introduced, although this implied that a new value for the kinetic constant had to be determined. The mass balance was solved by establishing an analogy with unsteady heat conduction in a finite solid slab (Bird et al., 1960). A computer program was developed for simulation purposes, but it led to time-consuming calculations. To cope with this inconvenience an approximation reducing the mathematical model to ODE's was obtained and compared against experimental results.

II. METHODS

A. Geometry of the problem and basic modeling

Extensive studies have proven that chip thickness is the most critical dimension in alkaline pulping processes. It has been shown that increasing chip thickness results in a lower pulping rate (higher screening rejects); chip length and width having minimal influence on delignification (Akhtaruzzaman and Virkola, 1979). Pulp chips have a length-to-thickness ratio of about 5:1. In alkaline pulping the diffusivities in the three primary directions are very similar (Rydholm, 1965; Hartler, 1962); therefore, chip thickness is the critical geometrical parameter with respect to diffusion.
Figure 1 shows the basic chip geometry adopted hereafter for modeling purposes. At early stages of the diffusion- reation process the expression of the reaction rate derived by Zanuttini and Marzocchi, (1997): , where A: acetyl concentration in [% acetyl / wood orig.], S: NaOH concentration in [g / l], a » 2, b »1.35, may be approximated by , where c: NaOH concentration in [mole / l].

This expression is proposed because of its simplicity; it still reflects the dependence on both the alkali concentration, regarded as the active agent, as well as on the limited capacity of the acetyl groups for accepting reaction. The form of the reaction term remains nonlinear but it is more suitable for simulation purposes, since it depends on only one variable c. Also, it shows a similar qualitative behavior: when the concentration of alkali in wood c rises, then the term (c0 - c) tends to 0 (thus curbing the reaction rate), and so decreases the acetyl concentration A with the extent of reaction.
This new expression of the reaction rate leads to the mass balance:

(1)

Dimensionless variables are defined

and then the governing equation takes the form

(2)

where

(3)

The term K c₀ in expression (3) may be regarded as a dimensionless number reflecting the relative weights of diffusion and reaction.

B. System parameters: (1) Diffusion coefficient

An experimental set up was designed in order to estimate the diffusion coefficient of alkali solutions in poplar wood chips in the radial direction. It is well known that this is the most frequent direction for chip thickness in the pulping industry, due to the way knives cut the wood during the chipping process.
A diffusion cell was built similar to a prototype described by Lönnberg and Robertsen, 1992. It consists of two chambers connected by a 40 x 40 mm channel, in which wood slices varying in thickness from 1 to 10 mm fit. The only diffusion contact between the two chambers occurs through the mounted piece of wood. The 1-liter chambers are both fitted with stirrers to eliminate concentration gradients within the chambers. The dilute side chamber is equipped with a 1/cm conductivity measurement cell and a temperature compensation probe, connected to a YSI model 34 conductivity meter. Both chambers are also equipped with N₂ inlets in order to generate an inert atmosphere and avoid alkali reaction with CO₂ from the surroundings. Figure 2 shows a scheme of the diffusion cell.

Assuming that there is molecular diffusion through the wood chip, that there it occurs in only one direction (radial), that the chemical reaction is already completed, and neglecting convective flow, then for each value of the time variable t:

(4)

where J : alkali molar diffusive flow velocity, ÑC : concentration gradient in the wood chip, D : radial alkali diffusion coefficient in wood.
Since J is by definition where A_i: interface area, and n : number of alkali moles, and n = C V, i.e. alkali concentration multiplied by cell volume, and assuming that V is constant in time, we can then write

Standard simplifying assumptions lead to

(5)

with z : half the thickness of the wood chip, C_A : alkali concentration of solution in side A of the cell (concentrated solution), and C_B : alkali concentration of solution in side B of the cell (dilute solution).
By assuming also that C_A ( = C ⁰_A : side A solution concentration) remains constant, and that the initial side B solution concentration C⁰_B = 0, then

(6)

Integrating (5) we obtain:

(7)

where and .
Conversely, derivation of ΔC in (6) with respect to t , for t = 0, results in:

(8)

From equation (8) and experimental values of , the diffusion coefficient is finally calculated,

(9)

Experimental procedure

• The wood samples used were radial poplar wood slices. The cuttings were made with a circular saw and afterwards polished with sandpaper, which resulted in probes of nearly 30 x 30 x 1.5 [mm].
• Working Conditions: Temperature = 25 [°C], C⁰_A = 1 [M], 2 z = 0.15 [cm], V = 1 [lt.], A_i = 8.41 [cm² ]. For additional details see Costanza, (1998).
• Three consecutive experiences were made with the same poplar wood sample, yielding the results shown in Figures 3 to 5.

The slope shown in Fig. 3 is lower than those in Figures 4b and 5. This can be explained by assuming that in this first experience there is chemical reaction together with alkali diffusion. The second experience leads to Figures 4a and 4b. Figure 4b expands the initial data of Fig. 4a. In the experience leading to Fig. 4a, diffusion was given enough time to confirm the exponential growth behavior derived analytically.
The slopes in Figures 4b and 5 are equal, which shows that the diffusion-only assumption was correct. By discarding initial noisy data and applying equation (9) the experimental diffusion coefficient is finally obtained

This value for D was found acceptable after comparison against data previously obtained and quoted in the bibliography (Lönnberg and Robertsen, 1992; Siau, 1984; Talton and Cornell, 1987).
Other experiments carried out with the same wood chips and in the same operational conditions, but at different values of C⁰_A , resulted in negligible changes in the values of D, which may therefore be considered as a constant for small changes in C⁰_A .

C. System parameters: (2) Reaction rate kinetic constant

As announced, a modification to the expression for the reaction rate given in Zanuttini and Marzocchi, 1997, is used in what follows. The new expression is more suitable for simulation under present conditions, but in order to make it compatible with the original one the following steps were taken:

• Calculation of the concentration values for conductivity data from Figures 3 (diffusion + chemical reaction) and 4 (diffusion alone).
• Calculation of the approximate slopes

where C'(t) = r: reaction rate, C_D (t) : Concentration for diffusion alone, C_DR (t): Concentration for diffusion +
chemical reaction, Δt= 1 [min.], and

• Estimation of the corresponding values for k (Figure 6) by dividing , which resulted in k » 0,7 for the early stages.

It can be seen from Fig. 6 that for initial times, k is nearly constant and does not depend on the value of C.

D. Mathematical treatment of the reaction-diffusion model

Equation (2) is first attacked in the pure diffusion case. This involves making K · c · C = 0 in (2) and solving Poisson's equation by separating variables, analogously to the solution of Laplace´s equation for the case of heat conduction in a finite slab (see Bird et al., 1960, page 352, for details).
By applying the initial condition C(0,n) = 1 and the orthogonality property of the basis functions the solution takes the form

(10)

The next stage in the mathematical treatment is basically an artifice to extend this approach to the case when there is also chemical reaction. By introducing the term K · c · C in (2) and proposing separate variables as if nonlinearity could be neglected, the result would be .  Assuming K · c to be approximately constant, and calling and , then: s² = g² - K · c, which leads to:

(11)

Recalling the initial condition and the orthogonality property of the basis functions:

(12)

This equation was numerically approximated. Concentration C profiles were obtained as a function of position x for increasing values of t. Problems arise when trying to find a numerical solution to equation (12) because
both sides of the equation are functions of C, so an iterative method was required. Results provided a first insight of the real process and its evolution, which turned out to be very useful for ulterior simulations.

E. Numerical treatment of the full equations

A finite difference scheme combined with iteration of profiles, much in the style of predictor-corrector methods, was used to treat equation (3) as an alternative to previously described approximation attempts. In this approach the quadratic reaction term was kept intact. With adopted values of D, k and c₀, and using small time increments, values of C at h = x = 0 were initially predicted by the solution of equation (13) (see section G below). As noted before, the fact that makes it necessary to iterate until both members of (12) converge to a single value of C( t,0). This was achieved by using an improved Newton-Raphson root-finding method. Values of C( t,0) (called C_N ( t,0) in advance) calculated in this way can be observed in Fig. 7 for different values of c₀.

The boundary condition C( t,1) = 0 t implies and then, from equation (3), . Also, symmetry requires . This information, together with the predicted values of C( t,0) and an iterative method that modifies C profiles until the terms of equation (3) are satisfied within a small tolerance, were used to generate solutions depicted in Fig. 8.

Both pure diffusion C_D, and diffusion-reaction C profiles, are plotted as soon as a fixed tolerance is achieved. In Fig. 8 a resulting plot for accumulative C and C_D profiles up to t = 2.425 min. is shown.
Light (C_D) profiles correspond to the case in which the only operative transport mechanism is molecular diffusion, while dark (C) profiles represent the combined effect of diffusion and chemical reaction.
Figure 8 is obtained by working with n from 0 to 50, 1 by 1, x from 0 to 0.075, with 0.001 increment, t from 0.003 to 2.425 min., with 0.003 increment, and c₀ = 1 [M].

F. Experiments and model validation

Wood slices similar to those used to determine the diffusion coefficient were placed in a water bath and vacuumed up to 30 mmHg. Under these conditions the slices sank in the bath in just a few minutes, and released air bubbles for several hours afterwards. They were given enough time for all the air to be evacuated. Before alkaline treatment, the four thin faces were sealed so as to prevent alkaline attack, taking special care in preventing the wood from drying out. Alkaline treatments were then conducted under the following conditions: T = 298 K (25°C); C = 0.25, 0.5, and 1 M; t = 3, 4.5, and 6 minutes.
In each case the final time was just the necessary to allow for model validation. Once each treatment time was reached, the chips were immersed into liquid nitrogen for at least two hours to quell the reaction. The chips were then stored in a freezer at -20 ºC. For chemical analysis 100-µm-thick slices were cut from the middle of the frozen chips, except in those cases in which the material had been excessively softened, where thickness was increased to 200 µm to preserve slices' integrity. The slices were immediately weighed, then submerged in 20-mL water and quantitatively neutralized to determine their alkaline load with 0.005 M HCl; this process carried on during enough time to permit alkali diffusion into the neutralizing environment, which was indicated by the phenolphthalein coloration. The dry weighing of the slices allowed their liquid content to be determined and expressed as g liquid/g treated wood, and the alkali concentration as g NaOH/L liquid. Alkali concentration experimental results compared with model predictions are shown in Fig. 9. Experimental data show striking agreement with experimental results. Only experiments at t = 3, 4.5 and 6 min. were performed due to the complexity of the technique employed and considering that the model is intended only for small times. Under present restrictions, the model appears fairly accurate.

Non-entire slices had to be discarded. The present method is limited because differences in concentration corresponding to distances lower than 100 or 200 mm (slices' thickness) could not be appreciated.

G. Simplifications and reduction to ODEs

A numerical solution to equation (3) through finite differences, as described in the previous section, has the main disadvantage that, in order to predict a diffusionreaction profile at a given dimensionless time t i.e. C_t( h) = C(t ,h), 0 < h £ 1 the whole scheme has to be run from  t = 0 to the desired  t value, using small Dt increments to guarantee accurate values for partial time derivatives and to avoid growing perturbations due to the discontinuous initial condition. Also, timeconsuming iterations arise in correcting concentration profiles at each time-step.
To minimize these inconveniences an approximation setup using ODE's was developed and contrasted against the results of the finite differences solution. This scheme is based in the particular forms of the original PDE's terms. By calling f ( t) º C ( t, 0), , and substituting these expressions in equation (3) for ( t,0) the following ODE is obtained:

(13)

Since the function h is still unknown, a series expansion for C(t,h)is used to generate some relation between f and h, which is necessary for the initial value problem to be properly defined.
As , then also .
This allows us to eliminate odd powers from the expansion and write:

(14)

Taking into account only the first two terms in the power series, and solving numerically, the results shown in Fig. 10 are obtained.

Qualitative comparison between Figures 7 and 10 indicates that a better approximation for the power series should be attempted. An independent ODE for h is then generated as follows. Notice that by differentiating (3) with respect to h
. Differentiating (3) twice with respect to h gives:

Let us now call , and since , then the additional relation appears

(15)

Taking now into account the first three terms of the second member of expression (14), the following differential equation system is finally obtained

(16)

The results for solving this system with a 4^th order Runge-Kutta method are shown in Fig. 11.

A new comparison, this time between Figures 7 and 11, allows us to conclude that an acceptable approximation of the alkali (NaOH) concentration in the middle of the chip as a function of time may be achieved through the method just described.
In order to further verify the profiles obtained in Figure 8, the following procedure may be followed. Integration of expression (3) with respect to h between 0 and 1 gives

(17)

where we call
For the pure diffusion case

(18)

Now defining , and , and subtracting expressions (17) and (18):

(19)

The form of each C profile depicted in Figure 8 and the boundary conditions of the problem admit a third order polynomial in h as an acceptable approximation

(20)

where C(t, 0) = nondimensional concentration in the middle of the chip.
Working with boundary conditions gives: a = 0, b = -3× g, and
Replacing all this in (20) allows approximating equation (19) by

where . Using now the estimations

the following ODE is obtained for d

(21)

which is solved together with (16). Results of numerical integration are shown in Fig. 12.

They are consistent with numerical integration of the first member of equation (19).
Most important are the qualitative features of the curves in Fig. 12, namely: For usual bulk concentrations, the curves grow from 0 to a maximum, and then decay, in accordance with the natural threshold implicit in the reaction rate expression.
The differences between pure diffusion and diffusion plus reaction grow with c₀, which is logical since the influence of concentration is essential in the case of reaction but marginal with regard to diffusion.
High values of d for all times mean that diffusion is unduly neglected.

III. CONCLUSIONS

Evidence generated in the laboratory compelled to distinguish between pure diffusion coefficients and combined diffusion-reaction coefficients previously published. Clarifying this issue and obtaining reliable values for the main process parameters was fundamental to pursue further modeling and simulation.
The relation between these coefficients indicates that diffusion and reaction mechanisms are both significant in alkaline pulping. Their relative importance depends mainly on temperature, which affects the diffusion coefficient and the kinetic reaction constant to different extents. In this paper simple numerical tools have been devised to predict the dynamics of the combined process accurately enough.
The concentration of alkali at the center of the chip is a measure of the completeness of wood deacetylation, which translates into the suitability of the final product for pulping purposes. This concentration can be predicted by integrating a pair of coupled ODE's with temperature dependent coefficients, resulting from a novel approach to the original diffusion-reaction PDE's.
Since alternatives combining both low- and hightemperature processes are currently being explored, the results in this paper provide useful equations for static and dynamic optimization analysis.

ACKNOWLEDGMENTS

The authors are indebted to Prof. Alberto J. Miyara, Facultad de Ingeniería - Universidad Nacional de Entre Ríos, for the valuable advice he provided during the edition of this work.

REFERENCES

1. Akhtaruzzaman, A. F. and Virkola, N. E., Paperi ja Puu 61, 629 (1979).
2. Bird, R.B., Stewart, W., and Lightfoot, E., Transport Phenomena, John Wiley & Sons, New York (1960).
3. Costanza, P., "Informe de Beca de Iniciación en la Investigación", Univ. Nac. del Litoral, (1998).
4. Costanza, V. and Miyara, A. J., "Modelización Parcial del Proceso de Impregnación de Maderas," Actas VI Congreso Forestal Argentino III, 591-594 (1988).
5. Costanza, V. and Miyara, A.J., "Dynamics of Hardwood Impregnation," Holzforschung 54, 183-188 (2000).
6. Donaldson, L. A., "Distribution of Inorganic Elements in Sub-Fossil Totara Wood," Holzforschung 47, 451-457 (1993).
7. Hartler, N., Paperi ja Puu 44, 365 (1962).
8. Jimenez, G., Gustafson, R., and McKean, W., "Modeling Incomplete Penetration of Kraft Pulping Liquor," J. of Pulp and Paper Sci. 15, 110-115 (1989).
9. Jimenez, G., Gustafson, R. McKean, W., and Chian, D., "The role of penetration and diffusion in nonuniform pulping of softwood chips," Tappi J. 72, 163-167 (1989).
10. Kelso, W. C., Gertjejanson, C., and Hossfield, R., "The Effect of Air Blockage Upon the Permeability of Wood to Liquids," Agr. Exp. Sta., University of Minnesota, Tech. Bull 242, 1-40 (1963).
11. Lönnberg, B and Robertsen, L., "Chemical Diffusion in Wood," Proc. AIChE Forest Products Sym., 49-55 (1992).
12. Rydholm, S., Pulping Processes, Wiley, New York (1965).
13. Siau, J.F., Transport Processes in Wood, Springer- Verlag, New York (1984).
14. Slattery, J.C., Momentum, Energy, and Mass Transfer in Continua, McGraw-Hill, New York (1972).
15. Talton, J. H. Jr. and Cornell, R. H., "Diffusion of sodium hydroxide in wood at high pH as a function of temperature and the extent of pulping," Tappi J. 70, 115-118 (1987).
16. Zanuttini, M. and Marzocchi, V., "Kinetics of alkaline deacetylation of poplar wood," Holzforschung 51, 251-256 (1997).
17. Zanuttini, M., Citroni, M., Martinez, M. J. and Marzocchi, V., "Chemimechanical Pulping of Poplar Wood, Alkaline Pretreatment at Low Temperature," Holzforschung 52, 405-409 (1998).

Received: October 24, 2000.
Accepted for publication: August 13, 2001.
Recommended by Subject Editor R. Piacentini.