[Home] [Journal papers]
peer reviewed logo

Johan Jansson, Freddie Grobler and Erik Dahlquist

Continuous pulp digesters are large reactors with significant residence time. They play a critical role in shaping the pulp qualities achieved at the product-end. For this reason, there is a big drive towards developing fundamental dynamic models to be used for better process understanding and control (especially in the development of model based controllers).

In this paper we discuss the use of physical models for diagnostics and process optimization of a continuous digester. 

Accurate diagnostics are not possible without the use of a pressure-flow net. The pressure-flow net also makes it possible to validate pressure and flow sensors in the circulation around the digester. For control purposes however, the pressure-flow net is not so important.

Two dynamic physical models were built for a continuous digester, one with a pressure-flow net and one without. The model without the pressure-flow net has been tuned and tested with good results at Ngodwana mill. It was not possible to get the model which included the pressure-flow net running with as many digester section blocks as we wanted due to the system becoming stiff. However, for a lower number of sections results were obtained. With an adjustment of the boundary conditions, higher number of sections should be possible to test as well.

Continuous digesters are used in the pulp and paper industries worldwide to remove lignin from wood chips. The most important final pulp quality use for control in a continuous kraft pulp digester is the kappa number. The kappa number is a measure of the residual lignin left in the final product after the cooking process (delignification) is finished.

The delignification process is dependent on several factors. The following are the more important ones:

  • the liquor concentration use for delignification (alkali to wood ratio) and the amounts of wood and liquor fed into the digester (liquor to wood ratio)
  • the temperature and residence time in the digester (indicated by the H-factor)
  • the wood specie and chip quality (dimensions and moisture).

The liquor used for the delignification process (strong white liquor – SWL) consists of two active components namely sodium hydroxide (NaOH) and sodium sulphide (Na2S). The concentrations of these two components are expressed as effective alkalinity and can be defined as:
digester_1 (1)

The ratio between NaOH and Na2S also impact on the extent of delignification. The ratio between NaOH and Na2S is expressed as sulphidity and can be expressed as:
digester_2 (2)

In most modeling exercises the digester is regarded as a continuous stirred tank reactor (CSTR) [Christensen, Smith, Albright and Williams, 1983] [Wisnemski, Doyle and Kayihan, 1997]. To increase the accuracy of the model the number of CSTR's in series can be increased, but this also impacts on the complexity of the model. If a simplified model is assumed, (making use of one CSTR), the following balances would be applicable:

  • total mass balance (if the liquor to wood ratio is constant, as well as the packing factor in the digester, the total filled volume in the digester would not change – making a total mass balance unnecessary)
  • lignin balance (wood balance)
  • EA balance (liquor balance)
  • energy balance

Just focusing on the lignin balance, it is well known that for the delignification of thin wood chips at a constant sulphidity and alkali concentration the following equation can be used [Kerr, 1970]:
digester_3 (3)

where:L = amount of lignin in the pulp mass at time t
Ci = initial EA concentration [Kg/m3]
k = reaction rate constant [Kg/s]

The reaction rate constant is temperature dependent according to the Arrhenius relationship:
digester_4 (4)

where:k0 = constant [Kg/s]
E = activation energy [J/ Kg]
R = universal gas constant [J/ Kg-K]
T = temperature [K]

For the kraft pulping process, this relationship has been determined as:
digester_5 (5)

where :A, B = constant (16113 and 43,2 respectively)

The solution of this time dependent lignin balance can be done provided that all the conservation equations (total mass balance, EA balance and energy balance) are solved simultaneously. 

Most digesters are controlled to keep the temperature and flow as constant as possible, depending on the wood quality and the recipe. It is however difficult for automatic control devices to compensate for changes in wood quality.
Changes in wood quality include differences in wood type, chip size, chip hardness and moisture content. Variations in wood quality result in a pulp production of variable quality and yield. It can also lead to disruptions in the process, such as hangs in the digester or channeling. These disturbances can result in a non-uniform distribution of cooking liquor between the chips and increased residence time.

This article describes the work and tuning result of two mathematical models of a continuous digester that incorporates a pressure-flow net and one without the pressure-flow net, [Aveling A et al 2006].

The model was built on the same basic principle as the Purdue model [Bhartiya et al 2001]. The main difference between this model and the Purdue model is that this model contains a pressure-flow network. The digester is discredited into a two-dimensional model with a number of volume elements both from top to bottom and from the centre and outwards. Each volume element has a specific size and the flow resistance between two volume elements is determined by the size of the open volume between the chips.

Modelica, with the Dymola graphical interface was used to construct the 2-D model. The model contains two types of streams. One stream is of free liquor and the other is of chips and trapped liquor.

The model in Modelica was constructed as two model blocks configured essentially as indicated in figure 1 and figure 2.

Figure 1: Overview of how the digester can be discredited in 2-D model.


Figure 2: Principle model around an extraction screen.

Figure 3 shows the rectangular volume elements within the digester and the valves between them symbolizing the pressure drops between them. The pressures at the boundary limits (BL) are determined at the outset.


Figure 3: Schematic overview of the graphical model

The equations for the chemical reactions in the digester section model are the same as those used and described by Wisnewski et al (1997). We have added hydraulic calculations so that the model takes pressure flow aspects into account [Eborn J, 2001].

In the pressure-flow network model, there is a pressure drop between each volume element. This is expressed as an admittance factor (A), which is equal to the contact area between two volume elements multiplied by the distance from the centre to the interface. The sum of the pressure drop is calculated as the sum of the pressure drop over the half volume element for the two volume elements in contact. So each volume element is sending this "inverse resistance" or admittance factor to four different valves. The admittance factor is independent of the geometry of the volume element and includes an adjustment factor for the chip size distribution. If all the chips are large the volume between them will be large and the pressure drop will be low, whereas if there are a lot of pins or flakes, the pressure drop will be larger.

The model also takes into account the dissolution of lignin. As lignin dissolves and the static pressure increases (as the wood chips move down the digester) the chips maintain their volume but become softer [Harkonen, 1987]. This leads to a higher compression and an increased pressure drop [Michelsen, 1995]. This is modeled as a reduction of the valve opening between volume elements. In the simplest version of the model, it is assumed that the reduction of lignin in the chip is proportional to the reduction of the valve opening (V).

The flow (Q) through the valves is calculated as a function of the pressure difference and pressure drop for the flow from BL to volume element

For element 1 the equations can be written as follows:
digester_6 (6)

digester_7 (7)

PBL1,1 is then applied to volume element 1 and used to calculate the pressure:
digester_8 (8)

P1 can be formulated as a function of the parameters from the connections in all four directions. For the downward connection it becomes:
digester_9 (9)

For the right side connection it becomes:
digester_10 (10)

For the left side connection it becomes:
digester_11 (11)

The variables in these equations are:
Px = pressure [Pa]
ρx = density [kg/m3]
hx = height [m]
g = gravitational constant [m/s2]

This information is combined with the corresponding information for the other valves and volume elements using one of the solvers attached to the Modelica / Dymola environment. The solver that has been used for this model is the "Dassl" solver for differential algebraic problems [Dymola, User Manual, 2004]. The direction of flow between two volume elements must also be determined from the pressure difference between them.

The pressures and the flows in the liquid are calculated from the dynamic and static pressures:
digester_12 (12)

The dynamic pressure calculation must take in to consideration how much the chip column is compacted.

The compaction (  ) is calculated both from how soft the chips are (Eq. 13) and the size of the chips (Eq. 14):
digester_13 (13)

digester_14 (14)

In the above equations:
P = pressure in a volume element [Pa]
Pin = pressure in an adjacent volume element [Pa]
ρl = density of the liquid [kg/m3]
A = area through which liquid is able to flow from the first volume element to the second [m2]
h = difference in height between the centers of gravity of the two volume elements [m]
g = gravitational constant [m/s2]
Fin = flow from the first volume element to the second [m3/s]
η lignin = compaction factor due to the dissolution of lignin from the chips
Zflis = original amount of lignin in the wood chips entering the Digester [kg/m3]
Z lignin = amount of lignin in the chips leaving the volume Element [kg/m3]
η flake = geometric impact of small flakes in between the larger chips causing an increased pressure drop
Z flake = number of flakes in the actual volume element [kg/m3]

Figure 4 schematic representation of the compaction of wood chips, which is dependent on the softness. The softness of the wood chips is dependent on the lignin concentration inside the chips

Figure 4: Schematic representation of the compaction of wood chips

The greater the amount of lignin that is dissolved from the wood chips, the softer the chips become and thus the greater the degree of compaction. This in turn increases the flow resistance in the free space between the chips. Large chips result in larger spaces between the chips than smaller chips due to geometric constraints.


The reaction rates Ri for dissolution of fast lignin, slow lignin and hemicelluloses are calculated by an Arrhenius equation:
digester_15 (15)

where: Ri = reaction rate for i (i represents either fast lignin, slow lignin or hemicelluloses) [Kg/s]
Ai = Arrhenius constant for the component [m/s2]
Ei = specific energy [J/ Kg]
R = gas constant [J/ Kg-K]
T = temperature [K]
Z = represents the concentration of the component i or hydroxide (OH) or hydrogen sulfide (HS) in the chip moving out of a volume element compared to the concentration in the white liquor [kg/m3]
  ε = void volume between the chips
ρ = density of the liquid [kg/m3]
V = volume element [m3]
c = chips

The Arrhenius constants are the constants that are adjusted for different wood species.

The following simplified equation is used to calculate concentrations, mass transfer of liquid into chips and extraction of lignin from the chips into the liquid in the volume elements.
digester_16 (16)

In the above equation dL/dt is the dissolution rate of lignin per time unit [Lindström M, 1977]. This rate is dependent on the concentration of hydroxide (OH), hydrogen sulphide (HS) and temperature (T) given by the Arrhenius expression [Lindgren C T, 1997]. The constants a, b, c, d and f are specific for each wood type and are related to the geometric dimensions, reactivity of the lignin and the diffusion rate into and out of the liquid, which is dependent on the density of the wood chip. Values for these constants for particular wood types may be obtained from the literature or calculated from process data, [Jansson J et al, (2004b)].

Equations 15 and 16 are relatively easy to tune against real data because the power function to  is a linear expression. These two equations are commonly used in the literature, [Grace T.M et al, 1989]

The chemical consumption of hydroxide and hydrogen sulphide is calculated from the reaction rates (R) and the stoichiometric coefficients ( ) for the liquid with respect to wood chips for fast and slow reacting lignin and hemicelluloses:
digester_17 (17)

The temperature (Tc) controls not only the reaction rate but also the diffusion rate RD where  m is the diffusion constant. This is expressed in the following equation:
digester_18 (18)

The diffusion rate is used to calculate the exchange of matter between the chips, the liquor trapped in the chips (Eq. 19), and the free liquor (Eq. 20).
digester_19 (19)

digester_20 (20)

Fc =flow of chips in or out of the volume element [m3/s]
F1 =corresponding liquid flow rate in or out of the volume element [m3/s]

When calculating the temperature we assume that it is the same for flows in and out. The following equation is used for it:

digester_21 (21)

ρ =density [kg/m3]
Cp =heat capacity [J/Kg-K]
T=temperature [K]

for either the liquid (l) or the chip (c), for matter entering (in) or exiting (out) a volume element.

The model with the pressure-flow net was tested by observing whether the chips reacted correctly. The primary focus of these tests was to investigate the changes in flow resistance as a result of compaction as lignin dissolved from the chips, and the effect of reducing free space between the chips by the introduction of additional flakes.

Figure 5 shows the calculation of the downward free liquor flow rate by the model.

A solution for the flow rate was found after 1000 s. This time could be reduced if more appropriate initial values were used. The average size of incoming chips at the top of the digester and the number of flakes was increased after 2000 s. The temperature inside the digester started to increase after 4000 s.

Figure 5 shows that the maximum free liquor flow rate occurred between level 1 and level 2 and that the flow rate decreased downwards through the digester. This decrease in flow rate was due to compaction of chips in the digester, which reduces the free space in the column thereby increasing the flow resistance. The compaction depends on the number of flakes in the chips and the degree of delignification that has occurred.

The initial delignification rate is very high because of the high concentration of lignin in the chips and the high concentration of solvent in the white liquor. The concentration of flakes was reduced after 2000 s, increasing the space between the chips and therefore lowering the flow resistance and increasing the flow rate.

After 4000 s the temperature inside the digester increased, also increasing the rate of delignification. This increased the compaction in the column, thereby reducing the flow. [Kayihan F et al, (2004)]..

So far the result obtained from the model that was developed at Ngodwana rendered reasonably good results of the kinds of in data we used, to get better and more useful results we need to connect the model on-line. The on-line tests will be the next step in the work.

The physical dimensions of the Ngodwana digester was entered into the model and was tuned with average values of data collected over a one year period at a 16 rpm production rate. The Arrhenius constants were adjusted to fit the output data of kappa in the blow line and the residual alkali in the upper and lower extraction. After a proper fit the model was verified at different production rates. The results obtained are given in table 1.

Production Rate (rpm)

C5 Residual
EA (g Na2O/l)

C15,C16 Residual EA (g Na2O/l)


Retention time (hrs)





Data (C15)




































Table 1 :. The tuning and verification result (average data from the PI system)

The result of 13 rpm is not very good, as a limited amount of data is available at this rate. The retention times are very close to that experienced in practice. Also a big variance between the model and the actual data for the C15 / C16 residual can be explained. The way the model is currently setup, one value is given for the lower extraction. This one value can be assumed to be an average value between the C15 and C16 actual values. Unfortunately only the residual on the C15 extraction is measured. From experience it can be assumed that the C16 extraction residual will be a lot lower due to the wash water that is added at the bottom of the digester and this will reduce the average value between C15 and C16 extractions. What is important to note here is the fact that the average data on C15 and the values predicted by the model (average of C15 and C16) are showing the same trend.


Figure 6: Step change in the chips feed with the response in the upper and lower extraction flow and also the kappa in the blow line.

Figure 6 shows the time it takes from the time you make a step change in the chip feed till the time you reach different levels in the digester

In the figure 6, the top graph shows a step change from 15 to 17 rpm in the chip feed. In the second graph it can be seen that after about 20 minutes the upper extraction flow starts to change and it takes about 1.1 hours before stabilizing again. The third graph from top shows the response of the lower extraction flow. The change is relatively small and takes longer to stabilize. This is due to the low solids process employed and the fact that most of the reaction takes place in the top of the digester. The last graph is an indication of the kappa in the blow line.

The response of the model was verified by making a production change from 14 rpm to 16 rpm by making use of the data obtained from PI. Two scenarios with different time constants were created to test the response. Table 2 shows the time constants for the two scenarios used.


Case 1
Time Constants

Case 2
 Time Constants

Chip feed



Liquor feed



C5 extraction



C6 flows



C6 temp



Wash liquor



Table 2: Time constants for changes.

The results of the two cases of production change can be seen in graphs 7 and 8.

Graph 7: The response in the residue alkali in the upper and lower extraction and also the kappa in the blow line for case 1.


Graph 8: The response in the residue alkali in the upper and lower extraction and also the kappa in the blow line for case 2.

From the above it is clear that the model reacts as expected and if tuned properly can be used to optimize production changes.

From the two models generated, only the model without the pressure net flow was tested on the #1 digester at Ngodwana. The reason for this is that the pressure flow net model can only be run with a limited amount of blocks. The more blocks used, the more difficult it becomes to solve the model. For our application at Ngodwana it is not necessary to use a pressure flow net model. The pressure flow net is more important if you use the model as a diagnostic tool, [Jansson J et al, 2004b].

The model without the pressure flow net showed great possibility. Although average data has been used the model shows great potential. Currently this model is only set up for a single specie but can also be used to incorporate swings between different species. For this the physical and chemical differences between the two types of wood need to be taken into account. Hard wood has shorter fibres, less lignin and higher density than soft wood.

Once the reliability of the model has been proven, it can be used for diagnostics, training and control. Problems within the digester can be predicted by comparing measurements to predictions by the model of what can be expected under normal operations. Mathematical models of reactors can also be used for advanced control systems such as MPC, [Jansson J et al, 2004a].

Avelin. A , Jansson. J and Dahlquist. E "Use of Mathematical Models and Simulators for ON-line Applications in Pulp and Paper Industry" Proceedings for the 5th Mathmod Conference in Vienna, February, 2006

Bhartiya, Dufour, Doyle, "Thermal- Hydraulic Modelling of a Continuous Pulp Digester", Proceedings from Conference on Digester modeling in Annapolis, June 2001.

Christensen, T., Smith, C. C., Albright, L. F., Williams, T. J., "Dynamic Modelling of the Kamyr Digester: Normal Operation Including Hardwood-Softwood Swings", Tappi Journal, November 1983, vol. 66, no. 11, p. 65-68.

Dynasim AB, "Dymola Multi-Engineering Modelling and Simulation",2004.

Eborn, J., "On Model Libraries for Thermo-Hydraulic Applications", Doctorial thesis LTH, Lund, dept Automation and control, 2001.

Grace T.M. , "Pulp and Paper Manufacture, Volume 5, Alkaline pulping", Appita, 1989, 3rd edition, p. 49-52.

Härkönen, E. J., "A Mathematical Model for Two-Phase Flow in a Continuous Digester", TAPPI Journal, 1987.

Jansson, J., Lindberg, T., Dahlquist, E., "Process Optimization and Model Based Control in Pulp and Paper Industry", AFCON Conference, Cape Town, 2003.

Jansson, J., Lindberg, T., Dahlquist, E., Persson, U., "Process Optimization and Model Based Control in Pulp Industry", TAPPSA Journal, November 2004.

Jansson, J., Dahlquist, E., "Model Based Control and Optimization in Pulp Industry", Conference Proceedings SIMS 2004 in Copenhagen, Sept 22-24, 2004.

Kayihan F. , "A stochastic continuous digester model to capture transition, compaction and chip size distribution effects" IETek 2002

Kerr, A. J., "The Kinetics of Kraft Pulping – Progress in the Development of a Mathematical Model", Appita, November 1970, vol. 24, no. 3, p. 180-188.

Lindgren, C. T., "Kraft Pulping Kinetics and Modelling; The Influence of HS, OH and Ionic Strength", Doctorial Thesis at KTH, Department of Pulp and Paper Chemistry, 1997.

Lindström, M., "Some Factors Affecting the Amount of Residual Phase Lignin During Kraft Pulping", Doctorial Thesis at KTH, Department of Pulp and Paper Chemistry, 1977.

Michelsen, F. A., "A Dynamic Mechanistic Model and Model-Based Analysis of a Continuous Kamyr Digester", Ph.D. Thesis, University of Trondheim, 1995.

Wisnewski, P. A., Doyle, F. J., Kayihan, F.,. "Fundamental Continuous Pulp-Digester Model for Simulation and Control",. AIChE Journal, December 1997, Vol 43, no 12, p. 3175-3192.

Johan Jansson (1)(3), Freddie Grobler (2) and Erik Dahlquist(3)
1: Sappi, 48 Ameshoff Street, Braamfontein, Johannesburg, South Africa
2: Sappi, N4 Highway, Ngodwana, South Africa
3:Malardalen University, Box 883, SE 721 23 Vasteras, Sweden.

[Home] [Journal papers]