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 c0 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 N2 inlets
in order to generate an inert atmosphere and avoid alkali
reaction with CO2 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 Ai: 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, CA : alkali
concentration of solution in side A of the cell (concentrated
solution), and CB : alkali concentration of solution
in side B of the cell (dilute solution).
By assuming also that CA ( = C 0A : side A solution concentration)
remains constant, and that the initial side B
solution concentration
C0B = 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],
C0A = 1
[M], 2 z = 0.15 [cm], V = 1 [lt.], Ai = 8.41 [cm2 ]. 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 C0A , resulted in negligible changes in the values
of D, which may therefore be considered as a constant
for small changes in C0A .
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, CD (t) : Concentration for diffusion alone, CDR (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: s2 = g2 - 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 c0, 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 CN ( t,0) in advance) calculated in this
way can be observed in Fig. 7 for different values of c0.
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 CD, 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 CD profiles up to t = 2.425 min. is shown.
Light (CD) 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 c0 = 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.
Ct( 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 4th 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 c0, 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.