© IETek 2002



FERHAN KAYIHAN              IETek – Integrated Engineering Technologies

                                         5533 Beverly Ave NE

Tacoma, WA 98422-1402

tel 253-925 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 product-end 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 non-trivial 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 co-current 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 co-current 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 counter-current 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 blow-line (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.


Continuous digester modeling has gone through significant evolutions starting with the H-factor 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 plug-flow 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. Liquor-solid 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. 



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 well-known problem on transition front dispersion related to the number of CSTRs used in plug-flow 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, steady-state 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 multi-structured 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 single-phase tubular reactor as seen in Figure 6. For the steady-state solution the time and space coordinates are interchangeable. Therefore, solving the batch reactor problem for the residence time is equivalent to solving the plug-flow problem through space-time. 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 steady-state 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 non-uniform 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, inter-particle solid and liquor density gradients are calculated in a manner similar to Gustafson [2] as compared to the lumped-parameter 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 kAi = kAoi exp (-EAi/RT) and kBi = kBoi exp (-EBi/RT).

Liquor component rates are related through stoichiometric relationships as

Where    RLG  = RS1 + RS2 ,    RC  = RS3 + RS4 + RS5    and   RS  = RLG + RC

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.









Stoichiometric coefficient for mass of effective alkali (EA) consumed/mass of reacting carbohydrate




Stoichiometric coefficient for mass of effective alkali consumed/mass of reacting lignin




Stoichiometric coefficient for mass of hydrosulfide (SH-) consumed/mass of reacting lignin




Diffusivity [m2/min]




Activation energies, i = 1,…,5  [kJ/gmol]

[29.3 115 34.7 25.1 73.3]

[29.3 115 34.7 25.1 73.3]


Activation energies, i = 1,…,5  [kJ/gmol]

[31.4 37.7 41.9 37.7 167]

[31.4 37.7 41.9 37.7 167]


Volume fraction of liquor, compaction




Kappa number = 654*lignin mass/total solid mass.




Pre-exponential factors [m3 liq/min/kgEA]

[0.2809 6.035x1010 6.4509 1.5607 10197]

[0.3954 1.457x1011 28.09 7.075 5826.7]


Pre-exponential factors [m3 liq/min/(kgEA kgEA)1/2]

[9.26 0.489 28.09 10.41 5.7226x1016]

[12.49 1.873 124.9 47.86 3.225 x1016]


Reaction rate effectiveness factor




Diffusion rate effectiveness factor




Liquor component densities, j = 1,…,4 [kg/m3 of liquor], indices are respectively: effective alkali (EA as NaOH), hydrosulfide (HS), dissolved solids and dissolved lignin




Solid component densities, i = 1,…,5 [kg/m3 of solid], indices are respectively: high reactivity lignin, low reactivity lignin, carbohydrate, galactoglucomman and araboxylan




Solid component densities in wood chips [kg/m3 of solid]

[72 280 755 33 360]

[150 225 675 75 375]


Unreactive portion of solid component densities in wood chips [kg/m3 of solid]

[0 0 435 18 0]

[0 0 535 18 0]


Chip (column) pressure [kPa]




Pressure drop [kPa/m]




Gas constant [kJ/gmol oK] = 0.0083144




Reaction rate of liquor component j [(kg/min)/(m3 of liquor)]




Reaction rate of solid component i [(kg/min)/(m3 of solid)]




Temperature [oK]




Relative fluid velocity [m/sec]






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 multi-volume compartments for “larger” chips using spherical geometry with equal total chip volume.





Figure 3. Alternative modeling approaches to moving-bed 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, 228-231,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:87-96, 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:3175-3192, 1997.

[6] F.A. Michelsen. A dynamic mechanistic model and model-based analysis of continuous Kamyr digester. Dr Ing Thesis, 1995 Report no. 95-4-W, 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 two-phase flow in a continuous digester. TAPPI Journal, 122-126, Dec 1987.






For additional simulation results see Advanced Digester Toolbox section.


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.


 Go to TOP of Page


For additional information or questions please contact


5533 Beverly Ave NE, Tacoma WA 98422-1402, USA

Tel: (253) 925-2179,  Fax: (253) 925-5023


© IETek 1996-2002, all rights reserved.

Home ] [ Advanced Digester Model ] Advanced Digester Toolbox ] Batch Digester ] Paper Machine Monitoring & Control Workshop ] ACC 2000 Monitoring Workshop ] DOE Digester Control Project ] Digester Benchmark Model v1 ] Digester Benchmark Toolbox v1 ] Custom Short Courses ] UDel PCMC ] Publications ] Projects ]