©
IETek 2002 FERHAN
KAYIHAN
IETek – Integrated Engineering Technologies
5533 Beverly Ave NE Tacoma,
WA 984221402 tel
253925 2179
Continuous pulp digesters are
large moving bed reactors with significant residence time. They play a critical
role in shaping the pulp and paper qualities achieved at the productend of a
mill and significantly influence energy balance and environmental burden of the
operations. There is renewed interest and much progress towards developing
fundamental dynamic models to be used for better process understanding and
especially in the development of model based controllers. One of the challenging areas in
fundamental model development is the incorporation of grade transition with
dynamic compaction (in packed bed of chips) and the transport/conversion effects
due to chip size distribution. The challenge arises not from the theoretical
understanding of how to include these factors into the fundamental model, but
from the considerable size of the resulting equations that can easily exceed
70,000 ODEs representing the dynamic behavior. This work is focused on the numerical aspects of reducing the model order to a tractable size so that process dynamics can be captured with a fast computational approach without any loss of information. Overall utility of this high fidelity stochastic dynamic model is in its ability to capture concealed process variability that could not be quantified with a distributed parameter approach constructed with hidden lumped parameters. Continuous digesters are very
complex vertical tubular reactors, used in the pulp and paper industry to remove
lignin from wood chips. Aqueous solution of sodium hydroxide and hydrosulfide,
called white liquor, is used to react with porous and wet wood chips. Usually,
continuous digesters are separated into multiple reaction and extraction zones
to carry out a specific process sequence. Depending on the production needs of a
pulp mill and on the state of the art of digester design at the time of
installation, there may be numerous differences between digesters. However,
common to all is the general sequence of transport and reaction processes that
govern the overall operation. Due to the complexities of these physical and
chemical phenomena and the fact that wood chips are nonuniform and constantly
changing, regulating product quality in a digester is a nontrivial task. The particular digester design
chosen for this article is the dual vessel EMCC (extended modified continuous
cooking) arrangement as shown in Figure 1. A brief description of the process is
provided to familiarize the reader with the basics. Detailed analysis and
descriptions are available in textbooks and PhD theses on the subject. Wet chips are steamed to remove
air in the pores and fed into the impregnation vessel (IV) together with white
liquor. In the impregnation vessel, white liquor penetrates into the chips and
equilibrates with initial moisture for about 30 minutes depending on the
production rate. In the IV, both chips and liquor move in the cocurrent
downward direction. From the IV, the chips are carried
into the top section of the digester with hot liquor that brings the mixture to
the desired reaction temperature. The top section of the digester, referred to
as the cook zone, is a cocurrent section where the main reactions take place.
Chips react from inside out owing to the significant internal pore volume and
associated surface area. Therefore, overall reaction rates depend on the
concentration levels of entrapped liquor and the diffusion rates from free
liquor that replenishes the active ingredient holdup in the pores. Spent liquor
saturated with dissolved solids at the end of the cook zone is extracted for
chemical recovery elsewhere in the mill. Chips follow into the MCC (modified
continuous cooking) and the EMCC zones, now countercurrent to fresh dilute
white liquor that simultaneously continue mild delignification reactions and
extract valuable inorganic solids from the pores of chips. As packed reactors, digesters are very unique in that the packing (main ingredient of the process) is continuously in motion, nonuniform in size and undergo through variable compaction both with respect to conversion and differential head pressure. Extent of reaction, defined through the blowline (exit) Kappa number, is the major performance measurement. Other important factors are the yield of the process and the fiber properties of the final product. Although various operating conditions may yield the same Kappa number, important fiber properties like strength are reaction path dependent. FEATURES
OF THE NEW GENERATION MODEL Continuous digester modeling has
gone through significant evolutions starting with the Hfactor approximation of
Vroom [1] through the fundamental developments from University of Washington by
Gustafson [2], Purdue University by Butler [3], Christensen [4] and Wisnewski
[5], and University of Trondheim by Michelsen [6]. Although the transport
phenomena related fundamental issues (e.g. diffusion and reaction) are modeled
with sufficient detail in published models we do not have all of the mechanistic
features formulated with sufficient rigor to handle transition dynamics
accurately. Recent efforts by Puig [7] at University of Delaware introduced a
mixing parameter approach to approximate transition behavior. Owing to current
keen industrial interest in managing transition dynamics to handle multiple
species and frequent rate changes in existing continuous digesters, we have
undertaken the task of developing a new fundamental model that provides the
necessary dynamic characteristics needed to track both species and transition
behavior without severe simplifying assumptions. In addition to this central
objective there is also certain interest to capture the strengths of all the
previous models into the new one to assure rigorous and high fidelity dynamic
predictions within the limits of our chemical and hydrodynamic understanding of
the process. The new model that will be summarized here contributes to the field
by (a) removing some of the limiting assumptions on solid phase behavior and (b)
introducing a novel and efficient numerical procedure. Main features of the new model are
the following: 1.
Solid phase moves in distinct “plugs” downward in the
digester satisfying the experimental evidence of plugflow behavior. Plug size
can be as small as desired for numerical accuracy. For example, about 200 solid
plugs in a 6hr residence time digester capture approximately 360/200 = 1.3 min
flow intervals. 2.
Chip size distribution can be represented by discrete size
mass fraction allocations. Usually 7 size cut categories are sufficient to
represent practical differences between chip thickness variations. Each plug has
its own representative mass fractions for predefined size cut dimensions. 3.
Solid, entrapped liquor and free liquor are the three
distinct phases in transport and reaction calculations. Liquorsolid reactions
occur between solid segments of each chip size and the entrapped liquor
contained within. Diffusion between entrapped liquor segments and free liquor
determine liquor density (composition) gradients and replenishment via free
liquor. Solid volume segments are represented in equivalent spherical shapes and
multiple discrete “shells” are used to capture diffusion and reaction
interactions within large solid particles (see Figure 2). 4.
Each plug entering the digester system carries its own
distinct set of physical, chemical and transport properties as model parameters.
Thus, there may be as many simultaneous wood species in the digester at any time
as there are distinct plugs. 5.
Solid compaction is a computed dynamic state for each plug
that depends on chemical and hydrodynamic properties of the digester.
SOLID
PLUGS AND THE CINEMATIC APPROACH Perhaps the most distinguishing
feature of the new model is the departure from the traditional CSTR
approximation of the tubular reactor behavior. Figure 3 shows the difference
between the stationary coordinate reference that is inherent in the CSTR
approach and the moving coordinates of the moving plugs with the new approach.
Although there are obvious simplifying advantages of the CSTR
compartmentalization of the plug flow reactor, for the transition problem we are
focusing on, there are certain undesirable dynamic side effects. Figure 4
quantifies the wellknown problem on transition front dispersion related to the
number of CSTRs used in plugflow approximation. Even with a large number CSTR
units like 2500 there remains approximately 30 minutes of artificial
“filtering” of the signal by the time it propagates through the digester. We
are very much interested in simulating transition fronts without numerically
generated forgiveness in process dynamics or having to employ impractically
large number of model states. In the simplest classical form,
steadystate behavior of a plug flow reactor with a uniform single phase can be
treated equivalently as either a stationary coordinate or a moving coordinate
problem with space or time as the independent variable respectively. However,
for dynamic considerations both space and time need to be considered
simultaneously and the equivalent partial differential equations need to be
solved. For single section reactors with explicit boundary conditions, and most
importantly with axial dispersion assumption, we can solve the dynamic problem
efficiently by employing orthogonal collocation methods. However, for the
problem at hand we have three severe limitations against using orthogonal
collocation approximations for numerical convenience: (1) the multistructured
nature of the reactor compartments and the mix of co/counter current flow
schemes as seen in Figure 5 make boundary conditions implicit and therefore
introduce severe (if not prohibitive) difficulty in implementing established
solution procedures; (2) requirement of reasonable axial dispersion in solid
phase to avoid numerical instabilities, which in tern would reintroduce the
artificial dynamic forgiveness we are trying to avoid; and (3) by expanding the
model to include (a) multiple chip sizes and thickness as additional independent
space variables and (b) dynamic compaction of the solid phase, we essentially
force the problem to stretch beyond the established scope of orthogonal
collocation procedures as applied to tubular reactors. The procedure used to accommodate
all of the model features into a practical numerical algorithm is the cinematic
approach. To explain it better let us concentrate on the simple problem of
singlephase tubular reactor as seen in Figure 6. For the steadystate solution
the time and space coordinates are interchangeable. Therefore, solving the batch
reactor problem for the residence time is equivalent to solving the plugflow
problem through spacetime. Furthermore, batch problem can be solved in multiple
small incremental time sequences without loosing any accuracy. Each one of these
snapshots of time can be considered as a sequence of instantaneous move in space
followed by a holding time equivalent to travel time for that space during which
reactions and other transport phenomena occur. Now consider the dynamic
conditions where both time and space are independent. The cinematic approach
will be applicable if one follows multiple discrete batch “plugs” instead of
just one as in the steadystate case. Here, physical transport through space is
an algebraic calculation while reaction and diffusion operations are time
dependent only. Obviously, space related approximations and the resulting
accuracy are controlled by the discretization of the plug size. Applying the cinematic approach to
the more complex problem of digester dynamics where multiple phases are involved
in a variety of flow direction sequences through a nonuniform reactor diameter
is a bit more complicated concern, but not impossible if the details of the
material balance conditions are carefully implemented. A significant advantage
of transforming physical transport (or convective flow) related model states
into algebraic computations is the complete elimination of ODE stiffness.
Otherwise, had we stayed with the CSTR approach and fully implemented material
balances to track compaction and flowrate dynamics, we would have a major issue
with the solution of stiff and sparse set of large number of ODEs. Basic steps of the numerical
algorithm can be stated as: 1.
Initially fill digester and IV with solid plugs and liquor. 2.
For a period of “sampling time = Δt” carry out
transport and reaction calculations (between phases in contact with each other
within the height boundaries of each plug) like a batch reactor. 3.
Compute solid compactions and rearrange solid effective
volumes and locations. As total digester volume is invariable allow for free
liquor movement to accommodate solid plug reshaping. All time dependent model
states are computed and available at this point providing full profile values. 4.
To initiate the next sampling time move necessary number of
plugs into, through and out of the digester (and IV) depending on specified
flowrates. Entrapped liquor by definition moves with the corresponding plugs.
Each new plug carries with it a unique set of physical and chemical properties
that stays with the plug through its residence in the reactor system. Total
number of plugs in the digester (and IV) may change depending on compaction and
level control. 5.
Move free liquor into, through and out of the digester (and
IV) depending on specified flowrates. Time dependent liquor properties are
introduced with new feed. 6. Go to step 2. Governing material and energy
balance equations of the new generation model are similar to those reported by
Wisnewski [5]. Details will be kept to a minimum here due to space limitations.
It is worth pointing out that, to track chip size and diffusion rate
limitations, interparticle solid and liquor density gradients are calculated in
a manner similar to Gustafson [2] as compared to the lumpedparameter approach
of [5]. Furthermore, dynamic compaction calculations are done with similar
assumptions and correlations as in Michelsen [6], which are originally proposed
by Harkonen [8]. Reaction rates for solid densities
are specified as
with reaction rate constants k_{Ai}
= k_{Aoi} exp (E_{Ai}/RT) and k_{Bi}
= k_{Boi} exp (E_{Bi}/RT). Liquor component rates are related
through stoichiometric relationships as
Where R_{LG } = R_{S1} + R_{S2 }, R_{C } = R_{S3} + R_{S4} + R_{S5} and R_{S } = R_{LG} + R_{C } Entrapped
liquor density diffusion coefficient is
Chip compaction under pressure is
given by the Harkonen [8] correlations as
Nomenclature
and parametric values are listed separately as a table. A
new generation continuous digester model is developed to provide high fidelity
simulation capabilities specifically suitable for specie and rate transition
dynamics. Undesirable “filtering” effects of CSTR approximation and
numerical difficulties associated stiff ODEs are avoided by using moving solid
plugs with a cinematic algorithm for numerical solutions. Dynamic compaction and
chip size distribution effects are easily incorporated in the model. Simulation
examples and additional model features will be demonstrated during the workshop.
Figure 1. Dual vessel EMCC
continuous digester.
Figure 2. Chip size distribution
and the approximation of diffusion effects on delignification rates by discrete
treatment of multivolume compartments for “larger” chips using spherical
geometry with equal total chip volume.
Figure 3. Alternative modeling
approaches to movingbed tubular reactors as applied to the digester problem.
Figure 4. CSTR approximation to
plug flow with 6hr residence time and the accuracy of transition front.
Figure 5. Traditional tubular
reactor structures and comparison to the digester model.
Figure 6. Application of cinematic
discetization of time for a simple tubular reactor.
[1] K.E.
Vroom, The H factor: a means of expressing cooking times and temperatures as a
single variable, Pulp and Paper Magazine of Canada, Convention issue,
228231,1957. [2] R.R. Gustafson, C.A. Sleicher,
W.T. McKean and B.A. Finlayson, Theoretical model of a Kraft pulping process, Ind.
Eng. Chem. Process Des. Dev., 22:8796, 1983. [3] A.C. Butler and T.J. Williams. A Description and User's
Guide for the Purdue Kamyr Digester Model. Technical Report 152, Purdue
University, PLAIC, Purdue Engineering, West Lafayette, IN 47907, December 1988 [4] T. Christensen, L.F. Albright, and T.J. Williams. A
Mathematical Model of the Kraft Pulping Process. Technical Report 129, Purdue
University, PLAIC, Purdue University, West Lafayette, IN 47907, May 1982. [5] P.A. Wisnewski, F.J. Doyle III, and F. Kayihan.
Fundamental continuous pulp digester model for simulation and control. AIChE
J., 43:31753192, 1997. [6] F.A. Michelsen. A dynamic mechanistic model and
modelbased analysis of continuous Kamyr digester. Dr Ing Thesis, 1995 Report
no. 954W, University of Trondheim [7] L.J. Puig, F.J.
Doyle III, and F. Kayihan. Reaction profile control of grade transitions in the
continuous pulp digesters. Control Systems 2000, Victoria BC, Canada, May 2000. [8] E.J. Harkonen, A
mathematical model for twophase flow in a continuous digester. TAPPI Journal,
122126, Dec 1987.
A copy of GUI showing nominal operating conditions for Softwood.
A copy of GUI showing dynamic response to 3C change in Cook, MCC and EMCC temperatures. Notice the changes taking place in chip pressure and compaction. Due to additional cooking, the chip column is getting lighter.
A copy of GUI showing dynamic response to 3C change in Cook, MCC and EMCC temperatures, about 6 hours after temperature changes.
A copy of GUI showing
dynamic response to 3C change in Cook, MCC and EMCC temperatures, about 16 hours
after temperature changes. Detail of the Kappa history shows that the blow Kappa
transient takes more than 12 hours (2xresidence time of plug flow reactor) to
settle. This prolonged transient is due to the countercurrent flow of liquor and
its nonlinear effects on process behavior.
A copy of GUI showing
dynamic response to 2% change (drop) in chip moisture. Response is captured
about 20 hours after moisture change.
A copy of GUI showing
dynamic response to 6 kg/m3 change (drop) in EA strength (measured as equivalent
NaOH). Response is captured about 20 hours after moisture change.
For additional information or questions please contact IETek 5533 Beverly Ave NE, Tacoma WA 984221402, USA Tel: (253) 9252179, Fax: (253) 9255023 © IETek 19962002, all rights reserved. 
