NUMERICAL INVESTIGATION OF INDUSTRIAL CONTINUOUS DIGESTERS By Yaoguo Fan B.Sc, Huazhong University of Science and Technology, China, 1994 M.Sc, Huazhong University of Science and Technology, China, 1997 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIRMENT FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Mechanical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA September 2005 © Yaoguo Fan, 2005 Abstract Digesters, either of batch or continuous type, are large chemical cooking vessels used to convert wood chips into fiber solution or pulp, from which paper products are eventually produced. In the continuous digesters, the complex liquor flow, which usually has a number of circulation loops, interacts with the downward moving wood chips through mass, momentum, and energy exchanges and therefore has a major effect on the chip motion and the cooking reactions. A set of comprehensive differential equations governing the two-phase flow, heat transfer, and chemical reaction is developed, with the purpose of computationally investigating the cooking process in continuous digesters. A new cooking reaction model is also derived, based on extensive experimental data. To model the chip phase motion, the chip phase is assumed to be a non-Newtonian fluid with a slip boundary condition on the digester walls. A number of sub-models are examined and used to account for the compressibility of chip column and the inter-phase friction force. Based on the governing equations and mathematical models, a general two-phase three-dimensional CFD code has been developed to calculate the flow fields for both phases, solid fraction, temperature, and the concentrations of chemicals. The developed CFD code also has the capability of solving many other potential multiphase reactive flow problems. Liquid flow in a laboratory digester filled with compacted chips is simulated by the developed code. The computations compare reasonably well with the experimental observations by ERT (Electrical Resistance Tomography). Currently, the code cannot predict the tracer motion. Therefore detailed comparison between CFD simulation and ERT images cannot be made. Simulations are also performed in a full-scale industrial digester. The predictions of yield and kappa number at the digester exit are accurate for various cooking conditions. The distributions of pressure and chip phase volume fraction are discussed, the chip ii phase velocity profiles at different digester levels are plotted, and the kappa number variation in the digester exit plane is also addressed. Simulations also show that an increased production rate results in a higher kappa number and a lower solid pressure at the bottom exit, and therefore the digester becomes more difficult or impossible to operate. The impact of variations of chip type and cold blow flow rate on cooking performance is also investigated computationally. The extensive investigation of the parameter variation impact on the process has shown that the model developed constitutes a reliable tool to predict the complex process in industrial digesters. iii Table of Contents Abstract ii Table of Contents iv List of Tables vii List of Figures viii List of Symbols ; x Acknowledgements xiv Chapter 1 Introduction 1 1.1 Paper making and kraft pulping 1 1.2 Continuous digesters 2 1.2.1 Batch cooking and continuous cooking 2 1.2.2 Modifications of continuous digesters 4 1.3 Objective of this thesis 6 Chapter 2 Mathematical Models of the Pulping Process 8 2.1 Literature review of digester modelling 8 2.2 Modelling assumptions 12 2.3 Governing equations 15 2.3.1 Continuity equation 15 2.3.2 Momentum equations 16 2.3.3 Species transport equation 17 2.3.4 The energy equation 19 2.4 Sub-models of the liquid phase 20 2.4.1 Flow resistance or inter-phase forces model 20 2.4.2 Compressibility equation for the chip mass 23 2.4.3 Concept of "equivalent diameter" of wood chips 25 2.5 Special treatment of the chip phase 27 2.5.1 Non-Newtonian fluid 27 2.5.2 Slippage on the wall boundary 30 2.6 Thermal properties of liquor 32 2.6.1 Density 32 2.6.2 Viscosity 32 2.6.3 Specific heat capacity and thermal conductivity 33 i v Chapter 3 The Kraft Pulping Reaction Model 35 3.1 Representative reaction models 35 3.1.1 Single variable model 35 3.1.2 Three phase model: Gustafson's model 36 3.1.3 "Two lignin" model 38 3.2 Motivation for deriving a new reaction model 38 3.3 Regression method 39 1. Take full advantage of the available experimental data 39 2. Select suitable regression variables 39 3. Avoid using global variables 40 4. Provide sufficient data to ensure the model accuracy 40 5. Abandon the "three stage" regression form 41 3.4 Regression procedures and the pulping model 41 3.5 Model validation 46 3.6 Conclusion 48 Chapter 4 Computational Methods 49 4.1 Solving Navier-Stokes equations by SIMPLE method 49 4.1.1 Navier-Stokes equations 50 4.1.2 Collocated grids 51 4.1.3 Discretization of the equations by control volume method 52 4.1.4 Solution procedures for steady-state problems 57 4.1.5 Curvilinear grids 57 4.2 Additional treatment of two-phase flow equations 60 4.2.1 Evaluate the mass flux 60 4.2.2 Discretization of the inter-phase interaction terms 60 4.3 Discretization of scalar equations 61 4.4 Imposing boundary conditions 63 4.5 The code and its structure 64 4.6 Code validation 66 4.6.1 Flow in a cavity 66 4.6.2 Liquid flow through porous media along a flat plate 71 4.6.3 Power-law non-Newtonian fluid flowing through a tube 73 Chapter 5 Modelling of Liquid Flow in a Laboratory Digester 76 5.1 Introduction 76 5.2 The model digester and electrical resistance tomography 77 v 5.2.1 The laboratory digester 77 5.2.2 Electrical Resistance Tomography 78 5.3 Computational method 80 5.4 Computational and experimental results 82 5.4.1 Approximation of boundary conditions 82 5.4.2 Liquid flow in the digester 84 5.5 Conclusions 89 Chapter 6 Simulation of an Industrial Continuous Digester 90 6.1 Introduction 91 6.1.1 Two-vessel hydraulic system in HSPP 91 6.1.2 Mill records of digester operations 93 6.2 Blended chip species 93 6.3 Computational domain 94 6.4 Boundary conditions 95 6.4.1 Boundary conditions for velocity 95 6.4.2 Boundary values for species equation (Lignin, Carbohydrates and EA).. . 102 6.4.3 Boundary conditions for the energy equation 107 6.5 Simulation results 107 6.5.1 Boundary values and model constants 107 6.5.2 The mesh 109 6.5.3 Computational results 110 6.6 Discussion of simulation results 118 6.6.1 Mesh convergence check 118 6.6.2 Predictions of kappa number and yield 118 6.6.2 Variation of production rate 119 6.6.3 Cooking performance affected by the chip type 120 6.6.4 Effect of the cold blow rate 125 6.7 Three dimensional simulation 126 6.8 Summary 132 Chapter 7 Summary and Conclusion 134 Chapter 8 Recommendations for Future Work 136 Bioliography 137 Appendix A Mill record of digester operation at HSPP* 145 Appendix B Average mill records of digester operation at HSPP 151 v i List of Tables Table 2.1 Flow resistance through packed wood chips 24 Table 2.2 The equivalent diameters of tested wood chips 26 Table 2.3 Coefficients for viscosity equations 33 Table 3.1 Constants fitted to pulping reactions 37 Table 3.2 Cooking conditions of cited experimental data3 44 Table 6.1 Average fiber length of different wood species 93 Table 6.2 Specific gravity of some common woods 94 Table 6.3 Selection of recorded data from HSPP digester (09-MAR-2004) 98 Table 6.4 Calculation of solid/liquor phase flow rate through the top inlet 100 Table 6.5 Chemical composition of some soft wood species 104 Table 6.6 Chemical composition of some soft wood in North America 104 Table 6.7 Composition of white liquor in HSPP (g/L) 105 Table 6.8 Data relevant to the modified cooking circulation (09-MAR-2004) 107 Table 6.9 Computed boundary values for the case run on March 9, 2004 108 Table 6.10 Computed kappa # and yield versus mill records 119 Table 6.11 Typical combinations of incoming wood chips in HSPP 121 vii List of Figures Figure 1.1 two-vessel hydraulic continuous cooking system (conventional) 3 Figure 1.2 Two-vessel hydraulic Kamyr EMCC continuous digester 5 Figure 2.1 Chip and free liquor phase in a REV 13 Figure 2.2 Inter-phase forces acting on different phases 16 Figure 2.3 Pressure drop over wood chips [7] 22 Figure 2.4 Solid fraction versus solid pressure [7] 23 Figure 2.5 Flow curves for a Bingham plastic and Newtonian fluids 28 Figure 2.6 Bingham plastic shear stress approximated by the Cross model [27] 29 Figure 2.7 Bingham plastic shear stress approximated by Pougatch et al's model [6].... 30 Figure 2.8 No-slip and slip conditions at the bottom wall 31 Figure 3.1 Effective alkali and sulfide consumptions vs lignin removal (h=5) 42 Figure 3.2 Fitted effective alkali and sulfide profiles 43 Figure 3.3 Lignin profiles by model predictions and experiments 46 Figure 3.4 Lignin profiles by model predictions and experiments 47 Figure 3.5 Model prediction versus experiments 47 Figure 4.1 Colocated grids in 2-D situation 52 Figure 4.2 Two-dimensional curvilinear grids 57 Figure 4.3 Illustration of symmetrical boundary conditions 63 Figure 4.4 Flow chart of the developed CFD code 65 Figure 4.5 Geometry and boundary conditions for 2D cavity flows 66 Figure 4.6 Flow field for 2D lid driven problem (Re=1000) 67 Figure 4.7 Velocity profiles in 3D lid driven problem (Re=100, mesh size: 25x25) 68 Figure 4.8 Velocity profiles in 3D lid driven problem (Re=400, mesh size: 25x25) 68 Figure 4.9 Stream function and temperature contours from Chen et al. (Re=100) 69 Figure 4.10 Stream function and temperature contours by this work (Re=100) 70 Figure 4.11 Streamlines in 2D lid-driven cavity with both side walls inclined at 45° 71 Figure 4.12 Streamlines in 2D lid-driven cavity by Ferziger and Peric's 71 Figure 4.13 Boundary layer flow passing through porous medium along a flat plate 72 Figure 4.14 Velocity profiles by code prediction and analytical solutions 73 Figure 4.15 Velocity distributions for power law fluids 74 vm Figure 5.1 Configuration of the model digester 78 Figure 5.2 The laboratory digester and the ERT sensor planes 80 Figure 5.3 Flow resistance over chip column, for reference pine chips 81 Figure 5.4 Plane jet: water injected from triple holes, m0=0.33kg/s 82 Figure 5.5 Plane jet: water injected from a single port, ra0=0.33kg/s 83 Figure 5.6 Computational domain, mesh, inlet and outlets 84 Figure 5.7 Flow field given by ERT image and CFD modelling 85 Figure 5.8 Flow field given by ERT image and CFD modelling 86 Figure 6.1 Two-vessel hydraulic digester system in Howe Sound Pulp and Paper 91 Figure 6.2 Illustration of the computational domain of the digester 95 Figure 6.3 Illustration of the liquor flow at the digester bottom 101 Figure 6.4 The extraction/injection ports of the modified cooking circulation 106 Figure 6.5 The structured 2-D mesh shown in reduced size for better viewing 109 Figure 6.6 Computational results of the operating case from HSPP I l l Figure 6.7 The liquor phase velocity and the chip phase pressure 112 Figure 6.8 The liquor phase velocity and the chip pressure at the extraction port 113 Figure 6.9 Velocity distribution of the chip phase 114 Figure 6.10 Chip phase velocity profiles at different 115 Figure 6.11 Kappa number variation at the bottom exit of the digester 116 Figure 6.12 Chip phase residence time versus chip locations 117 Figure 6.13 Computed downward liquor phase velocities using different meshes 118 Figure 6.14 Cooking performance affected by variation in production rate 120 Figure 6.15 Cooking performance affected by the chip type 121 Figure 6.16 Cooking affected by equivalent diameter of wood chips 122 Figure 6.17 Chip phase pressure affected by the equivalent diameter 124 Figure 6.18 Chip phase volume fraction affected by the equivalent diameter 124 Figure 6.19 Cooking performance affected by the cold blow rate 125 Figure 6.20 Liquor extraction shifted between screen sectors .' 126 Figure 6.21 The structured mesh shown in reduced size (original size:30x20xl80) 127 Figure 6.22 Velocity and pressure distributions for the liquor phase 128 Figure 6.23 Liquor phase streamline and the chip phase pressure 129 Figure 6.24 The velocity of the chip phase 130 Figure 6.25 Computed distributions of scalar variables 131 i x List of Symbols ^li ' ' '' coefficients in discrete equations constants for liquor viscosity A , B empirical constant in modified Ergun equation coefficients of discrete equations A x , A y , A z surface area of control volume Ck empirical constant C / " CPJ heat capacity of cooking liquor, kJ/(kg-°C) empirical constant in equation (2.22), MJ/(kg-m3) C mass fraction of carbohydrates cx convective flux in x direction Z) coefficient of diffusion rate; diffusivity, m /s 2 molecular diffusion coefficient, m /s Ddls hydraulic dispersion, m /s D e q equivalent diameter, meter active energy per mol for lignin, kJ/mol d ,exp //c// diffusive flux on e-face approximated by explicit scheme inter-phase friction force, Newton the constant of gravity, 9.8 m/s the hydrosulfide ion concentration, g/L k relative reaction rate of lignin k constant parameter in cross model species related constants in cooking reactions L, Lr mass fraction of lignin, % X characteristic slip length, in meter m, n power index in non-Newtonian models mn mass transfer rate between chip phase and free liquor phase mm mass flow rate of the entrapped liquor phase, kg/s OH hydroxyl ion concentration OH' hydroxyl ion concentration Pc pressure in chip phase, Pa Pf pressure in free liquor phase, Pa Pi prediction pressure, Pa p correction pressure, Pa Q,Q, heat generated or supplied, source term in i direction QA, QR volume flow rate in axi- and radial directions respectively R gas constant, kJ/mol/K Rt,R2 constants in flow resistance equation R, the generation rate of species i REV representative element volume [OH] mole concentration of hydroxide, mol/L [SH] mole concentration of sulfide, mol/L th heating time, s tc cooking time, s T time, s T temperature, °C Tc cooking temperature, °C velocity of chip phase, m/s velocity of free liquor phase, m/s u] prediction velocity, m/s xi Uj correction velocity, m/s Uj calculated velocity without pressure gradient, m/s ( * ] interpolated prediction velocity on the east face of a control volume, m/s Vo superficial velocity, m/s Vmv the volume of a REV, m Vfm tigmr volume of free liquor, m 3 VsoM ma, volume of solid material, m 3 Usiip, wsiip slip velocity at the wall, m/s Wj molecular weight of species i , g/mol Xj denotes x, y, z when 1=1,2,3 X dry solids concentration, % [Xj] mole concentration of species i , mol/L Y yield Yj mass fraction of species i Yj species mass fraction in cooking liquor Greek letters a dispersivity of porous media, meter, equation (2.17) a empirical constants in A-S equation (2.26) a dimensionless parameter in equation (2.44) acarb mass fraction of carbohydrates based on oven-dry wood ctijg mass fraction of lignin based on oven-dry wood a u over-relaxation parameter in Equation (4.18) xii p empirical constants in A-S equation p dimensionless parameter in equation (2.44b) K kappa number V the kinematic viscosity, m Is molecular conductivity of cooking liquor and water, w/(m-°C) Sc volume fraction of the chip phase Sc\ volume fraction of wood material £c2 volume fraction of entrapped liquor £f volume fraction of the free liquor phase Vf dynamic viscosity of the free liquor phase, Pa-S dynamic viscosity of the chip phase, Pa-S stress tensor for the chip phase stress tensor for the free liquor phase Pc density of the chip phase, kg/m3 Pf density of the free liquor phase, kg/m Viscosity, Pa-S TleflF equivalent viscosity, Pa-S A. thermal conductivity, w/(m-°C) 7 shear rate, s~' T shear stress, Pa volume of computation cell, m A viscous friction coefficient, Pa-S/ m Xll l Acknowledgements I would like to express my gratitude to all those who gave me the possibility of completing this thesis. I would like to express my deep and sincere gratitude to my supervisors, Professor Martha Salcudean and Professor Ian Gartshore, whose gracious support, warm encouragement, and timely guidance helped me all the time in this research work. I am also extremely grateful to both my supervisors for spending significant amounts of time reading and correcting my thesis, and for giving me timely feedback even in some difficult times. I am also deeply indebted to Professor Chad Bennington in the Pulp and Paper Center, Chemical and Biological Department, at the University of British Columbia, for offering the opportunity to simulate the liquid flow in his laboratory digester; this work has served as a good basis for this thesis study. Professor Bennington has also given me useful advice and suggestions for improvement. I would like to acknowledge with much appreciation the help and suggestions given by Professor Steven Rogak and Professor Carl Olliver-Gooch. As my research committee members, both professors have given me much assistance in better representing the transport behavior in the porous medium. I would like to give my special thanks to my master program supervisor, Dr. Jerry Yuan, for his kind support, encouragement, and valuable input to my research. Dr. Yuan has always been of great help, and his knowledge and ideas contributed significantly to my research work. Many thanks to Garry Pageau, Senior Process Engineer in Howe Sound Pulp and Paper, who has generously provided me with the mill records and other useful data; Garry has also offered lots of suggestions and advice, which contributed significantly to xiv this thesis work. I am also grateful to Drs. Peter Gorog and Marguerite Gorog from Weyerhaeuser Company, they have given me very useful information and advice. Also, thanks to Konstantin Pougatch for the very helpful discussions of the equations and non-Newtonian fluid model. Thanks to my colleagues, Xiaosi Feng, Suqin Dong, Zheqiong Wang, Mohammand Shariati for their help and friendship, which made my student life more memorable. Finally, I would thank my parents and my wife, Yamei, whose tremendous love, understanding, and encouragement enabled me to complete this work. xv Chapter 1 Introduction 1.1 Paper making and kraft pulping In the paper manufacturing, wood timber is transformed into wood chips which are then pulped through one of the following three pulping processes: sulphate or kraft, sulphite and Chemi-Thermo-Mechanical-Pulp (CTMP). The kraft pulping process involves the digesting or cooking of wood chips at elevated temperature and pressure in "white liquor", which is an aqueous solution of sodium sulfide and sodium hydroxide. The white liquor chemically dissolves the lignin that binds the cellulose fibers together. The kraft pulping process is used to produce the majority of pulp worldwide. The aim of kraft pulping is to chemically remove enough lignin so that the fibres are free and to give them some of the required characteristics for paper-making at the lowest possible cost. To evaluate the quality of the produced fiber solution or pulp, a large number of properties are measured or tested, among which are kappa number and yield. The kappa number and yield are important control variables; they are also chief variables in our model predictions; the definitions of these two variables can be found in Michelsen's thesis [1]: "Kappa Number is a term used to define the degree of delignification. It is commonly measured from laboratory tests, typically with one hour intervals. In some modern mills the kappa number is measured on-line in the blow line. The measurement is based on reaction with acidic permanganate solution (KMn04), giving the amount of lignin components in a sample from the pulp." "The yield expresses the total residual content of solid material in the pulp. It is defined simply as the sum of the residual content of lignin and carbohydrates. Hence, 1 Chapter 1 Introduction high yield means that a large part of the wood remains in the pulp, and the consumption of wood is low." Rydholm [2] defined the kappa number (K) mathematically by Equation (1.1): K = (%(aligl(alig+acarh\))lcK (1.1) The factor cK is 0.14 according to Christensen [3]. In most pulping wood, the initial mass fraction of lignin and carbohydrates ranges from 0.95 to 0.97; therefore equation (1.1) can be approximated by, /cx0.15 = % « % For instance, when wood chips are cooked to a kappa number of 30, the corresponding lignin content is 4.5%. The yield is expressed as: Y = %(an+acarh) (1.2) 1.2 Continuous digesters 1.2.1 Batch cooking and continuous cooking Kraft pulping occurs in either a batch or continuous digester system. For batch cooking, a mill normally uses several vessels. Wood chips and cooking liquor are added simultaneously to the digester. Then the digester is sealed and raised to a target operating pressure and temperature (approximately 800 kPa, 170°C) for one to two hours; the pulp is then blown into an atmospheric tank, from which it is pumped to the brown stock washing system. For continuous cooking, the most common continuous digester is the vertical down flow type. 2 Chapter 1 Introduction White liquor Wash liquor Figure 1.1 Two-vessel hydraulic continuous cooking system (conventional) Figure 1.1 shows a conventional two-vessel hydraulic continuous cooking system. The wood chips are preheated in a steaming vessel to remove air and some of the volatile wood constituents before the chips enter the digester. The chips are mixed with the cooking liquor and are fed into an impregnation vessel, where the chips are supposed to be fully impregnated. Chips are then sent to the top of the digester vessel so that they move downward by gravity. The hydraulic pressure in the digester vessel is kept at approximately 1,140 kPa, and the temperature is approximately 160 to 170°C. For typical softwood wood chips, one to one and a half hours of residence time is needed to ensure complete pulping reaction. The reaction is stopped by introducing a low temperature liquor flow (which is called quench flow) in the lower region of the digester, where pulp washing is carried out, normally using a countercurrent flow of the filtrate from the first 3 Chapter 1 Introduction brown stock washer. The pulp is blown from the bottom of the digester to a tank at atmospheric pressure. A batch pulping system enables a mill to pulp several different grades at once by using different digesters for different pulping conditions or fiber furnishes. The continuous process that produces pulp at a consistent rate with simple operations is more space efficient than the traditional batch process. Continuous cooking has gradually become the dominating pulp producing process. Today, over 65% of the kraft pulp produced in the world comes from continuous digesters, and most new digester installations are now of a continuous type. 1.2.2 Modifications of continuous digesters Until the early 1980s, most pulping processes used a "conventional" cooking mode, but modified cooking practices became available thereafter. Conventional cooking in continuous digesters means adding all the required white liquor within the chip feeding system, thus impregnation and cooking occur entirely in the cocurrent cooking zones. To decrease the bleaching chemical cost and lower the environmental impact, it is important to dissolve the lignin as much as possible in the cooking process. However, conventional cooking in very low lignin concentrations will cause severe cellulose degradation and thus decreases the yield and pulp quality. "Extended Cooking" originally proposed by the Swedish Forest Products Research Institute (STFI) in the later 1970s gained increasing popularity since the late 1980s [4]. With this technique, lower lignin content can be achieved with retained strength and pulp yield. Figure 1.2 shows a typical installation of a two-vessel hydraulic Extended Modified Continuous Cooking (EMCC) digester developed by Kamyr AB. 4 Chapter 1 Introduction Figure 1.2 Two-vessel hydraulic Kamyr E M C C continuous digester As compared with conventional cooking, the cook starts at a reduced concentration of effective alkali and ends at a higher concentration of effective alkali. The lignin concentration in the chip phase is reduced at the end of the cook. This is accomplished by splitting the alkali charge into different insertion points and by ending the cook in the countercurrent zone where the liquor flows in an opposite direction to the chips. As figure 1.2 shows, the cooking phase in the digester is divided into four different zones. The white liquor charge, comprised of sodium hydroxide and sodium sulfide, is added at four different points in the digester, instead of at just one point as with traditional continuous cooking. Two benefits are gained by adding white liquor to the countercurrent zone. Firstly the initial hydroxide concentration is lowered, which prevents fibers being partially damaged; secondly, some cooking is performed in the countercurrent zone where the dissolved solids are lower and cooking proceeds more effectively. 5 Chapter 1 Introduction 1.3 Objective of this thesis In kraft pulping, digesters play a critical role in shaping pulp and paper qualities. An ideal cooking process leads to reduced consumption of cooking liquor and external heating, thus the recovery load is reduced; it also minimizes the kappa number variation of the blowing pulp and improves pulp washing. However, the physical phenomena in continuous cooking process are quite complex involving chemical reactions and transport phenomena like multiphase flow, mass, momentum and heat transfer. In a commonly used EMCC digester, the counter-current liquor flow combined with a number of liquor circulations retards the motion of the solid chip mass, provides the chips appropriate reaction conditions such as high cooking temperature, enough chemicals and appropriate residence time. The cooking performance is essentially a combined result of many variables. Although industries have achieved great progress in continuous cooking, many underlying physical and chemical mechanisms are still not well understood. Recent studies [5,6] have shown that there is a need to improve further the cooking performance, and therefore efforts are still required to optimize the process and equipment design, operation and control. Some experimental work has been carried out to understand the liquor and chip flow in digesters either in plants or laboratories, but this is still very limited since making measurements in a digester is very difficult. Moreover, these measurements cannot fully represent the process. On the other hand, Computational Fluid Dynamics (CFD) has become a powerful and widely used tool in the investigation of industrial processes. The growing success in CFD technology has sparked an increasing interest in a broad range of industries. In the pulp and paper industry, kraft pulping has also drawn attention of CFD researchers; important CFD work on kraft pulping was carried out by Harkonen [7], Michelsen [1,8], He and Salcudean [5], Pougatch and Salcudean [6,9], Miyanishi and Shimada [10]. The objective of the present work is to investigate computationally the continuous digester in order to improve process understanding and develop a tool that can be used for process optimization and apply the tool to an industrial digester. 6 Chapter 1 Introduction The thesis presents the research undertaken to achieve this objective as follows: In chapter 2, the two-phase two-fluid models to describe the flow process in digesters are presented; issues related to chip compressibility and the wall stress model are also addressed. A general pulping reaction model is derived in chapter 3. Compared with most existing models, the correlation does not appear in a piece-fitting form, thus no transitional parameters are needed; besides, this model is obtained from extensive experimental data and thus gives more accurate predictions over a wide range of cooking conditions. Based on the governing equations, a general 3-D two-phase computational code is developed in chapter 4. The code uses non-staggered mesh and has a second-order of accuracy. The discretization of the governing equations is described in detail, and various validation cases are also given in this chapter. Chapter 5 gives the simulation of the liquid flow in a laboratory digester where stationary wood chips are added with compaction. The simulation results are then checked using Electrical Resistance Tomography images. In chapter 6, simulations of full-scale industrial digesters are presented. The velocity profiles, pressure distributions etc. are discussed; the impacts on cooking performance caused by production rate, chip type and cold blow rate are also investigated. A brief summary of the thesis is given in chapter 7, important conclusions are also summarized in this chapter. Chapter 8 contains the recommendations for future work. Chapter 2 Mathematical Models of the Pulping Process This chapter gives a brief review of previous work in digester modelling. The mathematical model used in the thesis is then described; and the main assumptions on which the model is based are stated. 2.1 Literature review of digester modelling As mentioned in Chapter 1, in the kraft pulping process, strong interactions through mass, momentum and energy exchanges occur between the cooking liquor and the downward moving chip phase; there are also multiple flow circulations and various disturbances from inflows. Awareness of the complexity of the pulping process and its importance to the pulp and paper industry, prompted work on digester modelling over the years lead to two different methods. In the first approach, both liquor and chips are stationary and chip compressibility is usually neglected, whereas diffusion and chemical reactions occurring in the pulping process are of particular interest. One of the earliest kinetic models was developed by Vroom [11], who established an Arrhenius type of relation in the form of %- = -k(T) (2.1) at where k(T) = ae-plT Here, L is lignin content and T is the temperature, a and P are positive constants. Vroom then further derived the concept of " H " factor, 8 Chapter 2 Mathematical Models of the Pulping Process H-factor= \ e { A x m i T ) d t (2.2) By the above integration, the H-factor combines temperature and time into a single variable representing the extent of the cooking; it has been widely used in the control schemes in today's pulping industry. In Kerr's model [12], the pulping process consists of five reaction steps: (1) Transport of chemicals to the site of reactions by flow of liquor into a chip (penetration) and molecular diffusion; (2) Adsorption of chemicals on the reaction site; (3) Chemical reactions; (4) Desorption of reaction products; (5) Transport of reaction products away from the site by diffusion. Kerr expresses the delignification rate as a function of effective alkali, lignin content and temperature, which is more complicated than Vroom's model. By virtue of Vroom's method, the connection between lignin content and " H " factor is established through a single equation in which numerous pulping variables are incorporated. In this model, the relationship between effective alkali concentration and pulp lignin content in initial and bulk phases was approximated by two different expressions. This makes the model very complex. Kerr later found that the relationship could be well approximated by a single expression, and this greatly simplifies the model, while the predictions are still accurate [13]. Gustafson [14] provided a good review of early pulping models, and presented a comprehensive model to study the combined diffusion and kinetics within a wood chip during kraft pulping. Assumptions made in his work include fully impregnated wood chips, well-stirred homogeneous bulk phase, constant sulfide concentration and so on. Kinetic models and data were carefully selected, and the model predictions were shown to be remarkably good when compared to some available laboratory cooking experiments. 9 Chapter 2 Mathematical Models of the Pulping Process The other approach, mainly applied in the continuous cooking process, recognizes the importance of fluid motion on the pulping process performance, and attempts to simulate the overall pulping performance resulting from liquor and chip flows, cooking reactions, heat transfer etc. By assuming plug flow for both chips and alkali solution, Christensen et al. [15] proposed a comprehensive one-dimensional model for continuous digesters. To describe the pulping chemistry, wood is considered as a mixture of four families of compounds: cellulose, hemicellulose, lignin and extractive compounds. In this study, entrance effects, exit effects and cross flows near the extraction ports were assumed to be of minor importance; in addition, it was assumed that temperature and concentration gradients occur only in the vertical direction. A good attempt at simulating the liquor flow two-dimensionally in a digester was made by Harkonen [7]. In his study, a chip is composed of soluble and insoluble materials, and the chip mass is treated as an isotropic porous continuum. Mass and momentum conservation equations are written separately for the solid phase and liquid phase and inertial terms are ignored. Interphase drag is calculated using Ergun's Equation. An energy balance in both phases includes the convective and conductive terms, reaction heat and heat exchange between the chip and liquor phases. A test digester with a diameter of 190 mm and a height of 360 mm was constructed to investigate experimentally the chip compressibility and flow resistance through the chip mass. In 1995, Michelsen [8] proposed a more comprehensive one-dimensional two-phase model, where convection, reaction, diffusion and dispersion were all included. In this model, wood chips are assumed to consist of lignin and carbohydrates, with the free liquor consisting of effective alkali, dissolved solids and water. The first order pulping reaction model given by Christensen et al. [15] was adopted. The inter-phase diffusion rate coefficients for the alkali and dissolved solids were adopted from Gustafson et al. [14] and Christensen et al. [15]. In steady state situations, the model gives reasonable predictions to the compaction and pressure of the chip mass, the velocities of free liquor and chip plug, kappa number etc. However there are limitations due to the one-10 Chapter 2 Mathematical Models of the Pulping Process dimensional approach as there could be significant radial variations in flow and associated phenomena, especially in modern digesters with large dimensions. The dynamic behavior was also studied for particular variations that occur in Michelsen et al's work. Miyanishi and Shimada [10] simulated the Low-solids cooking process and made some comparison with conventional cooking. The computational results are impressive because even the cellulose degradation has been well captured. The previous work cited has made a valuable contribution to the understanding of the kraft pulping process. A general deficiency of these models lies in their one-dimensional basis, so that many important cooking features, like secondary flow, cross-flow diffusion and non-uniformity cannot be predicted. Harkonen's two-dimensional simulation is based on a potential flow treatment and does not address sufficiently the chemical kinetics. On the other hand, the constantly increasing need for a high quality pulp at low cost always requires better understanding and more precise control. Over the past years, continuous digesters have been greatly enlarged and have incorporated novel technologies, such as E M C C and Low-Solids cooking. The increased process complexity has raised many practical concerns, such as safe operation and uniformity of cooking. Hence, it is important to develop a more comprehensive model capable of simulating the pulping process that occurs in digesters. It would be worthwhile to compute the flow three-dimensionally as in some cases the liquor injection or extraction does not take place axi-symmetrically and that can result in circumferential non-uniformity of the pulp. Preliminary work [5,6] has demonstrated the possibilities for numerical modelling of digesters but resulted in a numerical code that was not robust and exceedingly difficult to converge. The assumption that the chip phase was a Newtonian fluid is a weakness of this model that resulted in some unphysical behavior. The work in this thesis starts with a fresh and different approach based on a collocated grid (as opposed to the staggered grid used in the earlier approachs). The present work is also focused on the simulation of a modern industrial digester (Howe Sound Pulp and Paper M i l l at Port Mellon BC) that would allow extensive testing and 11 Chapter 2 Mathematical Models of the Pulping Process comparisons. Currently, after significant development of the preliminary numerical scheme, both codes have the capability of performing computations of a digester. The earlier approach has now been developed successfully and a description of its use and application has been accepted for publication [6]. A comparison between sets of numerical results produced by the two different codes for a single application has been made. This will be described in Chapter 4 of the present thesis. In the following sections, the assumptions and differential equations used in the present thesis are described; sub-models of inter-phase forces and chip compressibility are also discussed. 2.2 Modelling assumptions Wood chips used in the typical industrial cooking process are very small, with a thickness ranging from 1 to 10 mm. Practical industrial digesters, on the other hand, have large dimensions; for a digester with a 1,000 ton/day capacity, typical dimensions are 50 meters in height, and 6 meters in diameter. Even a small digester has dimensions which are much bigger than wood chips, so that detailed CFD computations at a scale of a single wood chip are obviously computationally prohibitive. It is therefore necessary to simulate the process from a "macroscopic" point of view. The cooking conditions in this industrial pulping process must be made as uniform as possible to minimize the non-uniformities in the pulp. In practice, digesters are operated carefully to ensure a long impregnation, enhanced diffusion, low gradients in temperature and alkali concentration. This makes a macroscopic approach more acceptable. In this approach, the mesh size is so big that a typical control volume or representative element volume (REV) contains a large number of chips, as schematically shown in Figure 2.1. The CFD calculation will provide the mean quantities such as velocities and the pressure in such an REV. 12 Chapter 2 Mathematical Models of the Pulping Process Liquor phase refers to the liquor that can move freely in the interstitial space of wood chips Chip phase includes solid wood chips and liquor entrapped in the pores of wood chips Figure 2.1 Chip and free liquor phase in a REV As shown in Figure 2.1, wood chips and the entrapped liquor are defined as the chip phase, while the liquor outside the chips is called free liquor phase or liquor phase, for simplicity. To simulate mathematically the flow in a continuous digester, we make the following assumptions: 1. Air is completely removed from the wood chips; as a result, the voids or cavities in the wood are totally filled with the liquid. We further assume that this part of the liquid (bound liquid) will have the same bulk velocity as the wood chips, so that we will describe the bound liquor and the wood chips as the chip phase. The remaining liquor, often referred to as free liquor or the free liquor phase, will interact physically and/or chemically with the chip phase. In Figure 2.1, sj refers to the volume fraction of the free liquor phase, and s c refers to the volume fraction of the solid/chip phase, which includes two parts, one is solid material s ci and the other is entrapped liquor sC2. Specifically, we have, V, free liquor V chip phase (2.3) REV V solid material V entrapped liquor V v REV (2.4) REV 13 Chapter 2 Mathematical Models of the Pulping Process Immediately, we can get the following relationships: ec + ef = 1 (2.5) (2.6) These relationships are with reference to a control volume. 2. Wood chips have no specific shape. The chip mass is assumed to be isotropic and homogeneous. Individual chips, as we know, are neither isotropic nor homogeneous. However the chip mass, namely a collection of chips, tends to be distributed randomly when compacted in typical pressure range, therefore assuming the chip mass to be isotropic and homogeneous is usually a good approximation. This assumption is supported by the observation that no significant difference was found from a macroscopic point of view, when measuring diffusivities in the three primary directions [14,16]. 3. Solid material in wood chips consists of lignin, carbohydrates, and insoluble material. When lignin and carbohydrates dissolve and enter into the free liquor phase, the corresponding volume will be simultaneously filled by free liquor. 4. The chip mass is compressible, but each individual chip retains its volume throughout the whole cooking process. In other words, when solid pressure (stress) is imposed, the interstitial space in the chip mass will decrease; each single chip may change its shape but not change its volume. This proves to be a valid assumption from pilot plant experiments [8]. 5. The chip mass will be regarded as a moving porous medium.With this assumption, the digestion process can be considered as a continuous two-phase flow, where the free liquor (or unbound liquor) flows through a compressible moving chip bed. 6. The cooking liquor is assumed to be a solution of sodium hydroxide (NaOH) and sodium sulfide (Na2S), with the addition of a certain amount of dissolved wood substance. The main feature of the cooking liquor is that its ingredients (or species concentrations) and density are changing as chemical reactions proceed. The chemical 14 Chapter 2 Mathematical Models of the Pulping Process kinetics are characterized in terms of effective alkali concentration and sulfite concentration. 7. The flow is laminar. A continuous digester usually involves several liquor circulations, which supply the liquor to the digester typically through central downcomers. These inflows are basically liquor jets, where turbulence tends to occur even at low Reynolds numbers. However the flow in the digester has a very high solid fraction and the occurrence of turbulence over large regions is unlikely. Additionally, studies on fluid flow through porous media could give some idea about the flow pattern in a digester. Michelsen [8] estimated the Reynolds number based on the average grain diameter and found it to be in the range of 0.8-15 for typical digester operations; he then concluded that the flow is mainly in the laminar domain. 8. Pulping reactions are irreversible. Alkali diffusivity depends on temperature and the ingredients of the cooking solution. 9. Chips are in thermal equilibrium with the local cooking liquor, i.e., chips and cooking liquor in the same REV always have the same temperature. Because the temperature in a continuous digester does not change sharply in space, and solid chips only move on the order of millimeters per second, this assumption should be acceptable. 2.3 Governing equations Based on the previous work [5], the macroscopic equations for the two-phase flow are the following. 2.3.1 Continuity equation Mass conservations are expressed as: (2.7) —m ex (2.8) dt 15 Chapter 2 Mathematical Models of the Pulping Process Please note, in this thesis, subscript "f' and "c" refer to free liquor phase and chip phase respectively. In Equation (2.7) and (2.8), mcx denotes the mass exchange between the chip phase and the free liquor, and will be represented by the chemical kinetics and mass transfer models. 2.3.2 Momentum equations The single-phase momentum equation can be generally written as, f dp dr N + dx, dx (2.9) ./ J where p is the static pressure, r,y is the stress tensor (described below), pgj and Fj are the gravitational force and other body forces in the i direction, respectively. du, du. v dxj dx,. 2 du, M— s a 3 dx, " (2.10) where p is the molecular viscosity and the second term on the right hand side is the effect of volume dilation. Based on the assumptions made in section 2.2, a two fluids two-phase model will be used to predict the flow in the digester. Figure 2.2 schematically shows the inter-phase forces in the x direction for the free liquor and chip phases, where Ff denotes the friction force per unit volume between the two phases, and is considered as a body force. F, f <: (a) Free liquor phase (b) Chip phase Figure 2.2 Inter-phase forces acting on different phases (within the REV, x direction only) 16 Chapter 2 Mathematical Models of the Pulping Process From Equation (2.9) and the forces shown in Figure 2.5, the momentum equations for each phase can be written as, Jt (Pf£fufj) + -^T (Pf£fufJufJ ) = e. V dx: dx + EfPfgt + Ff,i + mexuc, (2.11) J J -r (Pc£c"cj ) + -r— (pcecuCJuCJ) = dt dx, dx + £cPcgi - Ffj - sc J J dx, (2.12) where Ffj is the inter-phase friction force acting on the fluid per unit volume in the i direction, and is actually the pressure gradient caused by the inter-phase friction force. This will be further discussed in a subsequent section. The stress rl is defined in the same way as in (2.10), and zctj is defined by Equation (2.14), as the chip phase is assumed to be a Bingham plastic fluid. K dXj dx, y 2 dulf rduhC | dUj^ Kdxj dx, ; 2 dulc (2.13) (2.14) 2.3.3 Species transport equation In general, every species equation is determined by chemical kinetics as well as diffusion. Early models in kraft pulping assume that the cooking process is kinetically controlled. This assumption is acceptable for batch cooking. But for continuous cooking, where both chip mass and liquor flow continuously and liquor circulations are often used, ignoring diffusion is no longer reasonable. Apart from this, many of the commercial chips used today are quite large, making the overall reaction rates essentially dominated by both chemical kinetics and diffusion. 17 Chapter 2 Mathematical Models of the Pulping Process As a matter of fact, efficient diffusion is desirable for industrial digesters because it results in fast mixing of chemical reactants giving high production; high diffusivity also helps to distribute products evenly, ensures complete reaction and therefore decreases non-uniformity of cooking. In the history of the development of continuous digesters, the introduction of a central pipe was considered to be a major advance as it enhances the diffusion and greatly decreases the pulp non-uniformity. Therefore, solving general transport equations for some important species or components is necessary. In the free liquor phase, transport equations for effective alkali or sulfite can be written generally as, er 1 '' dx,y ' 1 l j dx j J + R, (2.15) ./ J where Yj is the mass fraction of effective alkali or sulfite ions over the total mass of the free liquor; Rf is the mass rate of creation or depletion of each species by the chemical reaction model, which will be described in detail in chapter 3; Dj is the equivalent diffusion coefficient, which is considered as the sum of a contribution from molecular diffusion and from hydraulic dispersion: D = D0+Ddls (2.16) In Equation (2.16), D0 is the molecular diffusion coefficient, Ddis is the dispersivity from hydraulic dispersion. Scheiddegger [17] found that Ddis depends linearly on the velocity and a characteristic of the porous media referred to as dispersion coefficient; he proposed the following model: Ddls=a\U\ (2.17) where U is superficial velocity in m/s and a is the dispersion coefficient of the porous media in meter. The average grain diameter is considered as a typical value for the dispersion coefficient of homogeneous sands [18]. However, the dispersion coefficient for wood chips has not been investigated yet. Huang et al. [19] computationally studied 18 Chapter 2 Mathematical Models of the Pulping Process the diffusion of salt water in soil, using a dispersion coefficient of 0.01 meter. In Huang's computations, the porosity was in the range of 0.35 to 0.62, which is typical for our digester simulations. Since there is no information available in literature for wood chips, we assume in our study a dispersion coefficient of 0.01 meter. However, with the same porosity as the soil, compacted wood chips have different mean grain diameter, thus different mean free path, and therefore the dispersion coefficient is expected to be different as well. The superficial velocity difference in a digester is in the order of 10"3 m/s. Based on equation (2.17), the diffusivity caused by hydraulic dispersion is expected 5 2 to be of the order of 10" m/s. A sensitivity study in digester modeling has shown that when a is less than 0.1 meter, the variation of a does not influence significantly the distribution of species concentrations. However, large values of a result in noticeable difference in the distribution of species concentrations. The mass fraction Yj in Equation (2.15) can be easily converted to its molar concentration by the following expression: (2.18) Pf And equivalently, Equation (2.15) can be written as dip f[X,]) dt dx; v 7 /L , J ; dx dXj j Pf + ^ - L (2.19) W, The species equations of lignin and carbohydrates are quite similar to Equation (2.15). Because lignin and carbohydrates do not diffuse freely in the solid phase, i.e., D, =0, the corresponding species equations can be simplified to, (^Pc^ )+/-U"/jr/)=^ (2-20) dt dx , 2.3.4 The energy equation Chemical reactions occurring in cooking processes are temperature sensitive. In normal cooking conditions, the delignification rate approximately doubles for every 8-19 Chapter 2 Mathematical Models of the Pulping Process 10°C increase in temperature, and the diffusivity of the reaction agents also strongly depends on temperature [20]. Therefore accurately predicting the temperature distribution in the digester is very important. The heat released in the cooking reactions is not very high, and therefore in normal cooking conditions, the chip phase and free liquor are assumed to have the same temperature within each REV. Thus the following energy equation applies to both phases: —(pfc„ ,T)+—{pfUf ,c„ — dtv y p j ' cbc, v ' f j p J ' dxj . dT f Dp, d{u,T,,) + ^LL + ^±JL> + Q (2.21) Dt dxj (dl dC\ Q = CR\— + — (2.22) \dt dt) here, CR is an empirical constant, with a value of 0.218 MJ/kg/m for dry wood as suggested by Gullichsen [21]. On the right hand side of Equation (2.20), the first term is the heat transport by conduction; the second term accounts for the heat converted from flow work; the third term is the heat generated by viscous dissipation; and the last term Q includes the heat of chemical reaction and external heat supply. The second and third terms are negligible for the digester computations. 2.4 Sub-models of the liquid phase 2.4.1 Flow resistance or inter-phase forces model Flow through a porous medium is very complex, and extensive experimental work has been conducted for different combinations of solids and fluids. To predict the pressure drop resulting from flow through a porous medium, various models have been published in the literature, often in the following two popular forms [22]. A. Ergun and modified Ergun equation [23] form (M-E form) Fk.— = B\ A(l~£) (2.23) where 20 Chapter 2 Mathematical Models of the Pulping Process dPIdL pV02 1 NM = (2.25) M where s is the porosity of the porous medium and D e q is an appropriate characteristic length for the medium. B. Ahmed and Sunada equation [24,25] form (A-S form) dP/dL _ pV0 = a + p—°- (2.26a) pV0 p dP Or: = apV0+fipV02 (2.26b) dL where dP/dL is the pressure gradient in the direction of the flow; Vo is the superficial velocity and is given by Vo=Q/A, where Q is volumetric flow rate and A is the total cross-sectional area perpendicular to the flow direction; a and P are model parameters to be established empirically. The Modified Ergun (M-E) equation [23] is well known as it includes the quantity porosity, which characterizes the porous medium. Also, it is believed that the M-E equation has a good theoretical basis. This equation does give very good results when the porous medium can be characterized fairly accurately. The A-S equation, on the other hand, does not explicitly include a parameter related to the porous medium. Therefore the parameters a and P must be functions of the medium and accordingly must be empirically established for each different medium. However, The A-S equation has the advantage of being simple and appears to be more accurate when fitted by a large number of experimental data, as concluded by MacDonald [22], especially when the porous medium is not well characterized. It is probably for this reason that most flow resistance models for liquor flow through wood chips are given in the A-S form. 21 Chapter 2 Mathematical Models of the Pulping Process Harkonen [7] investigated experimentally the flow resistance in the presence of wood chips and established the following correlation: dP dL f \ v sf £f J (2.27) For pine Ri=4.6xl02, R2=3.9xl06. The pressure drop strongly depends on velocity difference and solid fraction, as plotted in Figure 2.3. 5.0 r E Q. Q. < CL 2 •D £> 3 at to £ CL 2 4 6 8 10 12 14 16 Velocity difference, AV (mm/s) Figure 2.3 Pressure drop over wood chips [7] When both phases are in motion, such as in the digester, the formula should change accordingly. Michelsen extended the use of Equation (2.27) in his dissertation [8] and gave, dP dL sf v 6 n ef f s2 V-V, t - v , ) (2.28) where A is an expression of the viscous friction coefficient, and is a function of the porosity of the porous media and the velocity difference between the two phases. Here, both velocities are physical velocities or local velocities instead of superficial velocities which are also commonly used in literature. Equations (2.11) and (2.12) result in, dP_ ~dL Ji A R ^ + R2sc uc. •(ueJ-ufJ) (2.29) 22 Chapter 2 Mathematical Models of the Pulping Process 2.4.2 Compressibility equation for the chip mass As stated before, compressibility is one of the important features of the chip mass. The chips become weaker as lignin is dissolved from the cell walls. Thus compressibility is a function not only of solid pressure, but also of lignin content in the chip mass: the less lignin there is, the softer the wood chip becomes, and the higher the chip compaction is. Harkonen [7] has established the following correlation based on his experiments, and plotted in Figure 2.4. (2.30a) (2.30b) Solid pressure (kpa) Figure 2.4 Solid fraction versus solid pressure [7] Further work [29, 30, 31] on chip compressibility and flow resistance over wood chips was conducted and correlations of chip compressibility and flow resistance in the form of Equations (2.27) and (2.30a) were obtained. Brief description of the chips tested and the resulting correlations are listed in Table 2.1. 23 Chapter 2 Mathematical Models of the Pulping Process Table 2.1 F l o w resistance through packed w o o d chips (Uni ts : P c — Pa , d P / L — k P a / m , superficial ve loc i ty u - mrn/s ) R e f Raw material Equat ions Pine, size distribution not specified ef = 0.644 + ( P c / 1 0 4 ) ° 5 9 x ( -0 .831 + 0.139ln(x-)) — = 0.0046 u + 0.0039 ±L u2 29 Scandanavian pine, (Pinus silvestrus) Reference Chips: Fraction of accepts chips of sawmill chips (chips passed through a 13 mm round hole but retained on a 7 mm round hole, no other details) " N e w " Chips: 4 mm thick, 40 mm in length f ' new dP L ref dP 0.663 + (Pc 1104)°56 x ( -0 .788 + 0.133 ln(*r)) = 0.840 + (Pc 1104)°39 x ( -1 .092 + 0.161 ln(*r)) (l-sf)2 (\-£f) 2 = 0 . 0 5 2 - l—u + 0 . 0 0 1 5 - ^-u2 (l-ef)2 (l-£f) 2 0.082 1— u - 0.00011 ^ u2 30 Pinus silvestrus (Scandinavian) M i x 1 4~6mm:22.1%; 2~4mm:44.2% 2~3mm:29.1%; <2mm: 4.6% Chip bulk Density=165kg/m 3 M i x 2 4~6mm:23.2%; 2~4mm:46.3% 2-3mm: 30.5%; <2mm: 0% Chip bulk Density =163kg/m 3 M i x 3 4~6mm: 0%; 2~4mm: 60.3% 2~3mm: 39.7%; <2mm: 0% Chips: p =157 kg/m 3 Basic wood: p =401 kg/m 3 s'T = 0.604 + (Pc 1104)°63 x ( -0 .956 + 0.191 ln(/c)) = 0 .615+ ( P C / 1 0 4 ) 0 7 9 x ( -0 .910 +0.181 ln(O) = 0.647 + (Pc / 1 0 4 ) 0 7 4 x ( -1 .021 + 0.205 ln(O) = 0 . 0 2 8 ^ ^ , - 1 . 3 x l O - 4 ^ - ^ W 2 c»iix2 •v dP dP L = 0.051 K 1 - w - 3 . 5 x l 0 " 4 V , J - 2 'mix! dP = 0.0055 w - 5 . 7 x l 0 " ^ u 2 31 Eucalyptus camaldulensis >45mm: 0.6%; 8~45mm: 9.2% 8~13mm:51.4%; <3itim: 1.1% 7~13mm:31.2%; 3~7mm:6.4% Scandinavian b i rch >45mm: 1.1%; 8~45mm: 4.8% 8~13mm:78.8%;<3mm: 0.5% 7~13mm:12.6%;3~7mm:2.2%; Chip bulk Density=175 kg/m 3 Basic wood Density =556kg/m 3 hir f dP L,„„. dP -- 0.591 + (Pc 1104)°56 x ( -0 .645 + 0.148 ln(/c)) 0.630 + (Pc 1104)°64 x ( -0 .697 + 0.151 ln(*-)) 0 . 5 5 2 - l—u + 0 . 0 0 0 7 4 - ^—u2 ( I - * / ) 2 0 - f f f ) 2 0.28 u + 0.0012 ^ - u2 24 Chapter 2 Mathematical Models of the Pulping Process 2.4.3 Concept of "equivalent diameter" of wood chips The flow resistance through the chip phase could be very big and is usually a dominant factor that affects the chip compaction and hence flow in typical digester operations. Experiments have been conducted to measure the flow resistances over different wood chips and empirical correlations have been proposed as listed in Table 2.1. As can be seen, the correlations are all in the same form, where the flow resistance is linked to the velocity and the porosity of the chip bed; the influence of chip size, however, is not explicitly established in these correlations. Although wood chips have irregular shapes and non-uniform sizes, it is possible to derive an average diameter to characterize the whole chip mass; with this diameter, some popular semi-empirical models that describe the resistance of flow through porous media could also be used. In this section, the "equivalent diameter" of wood chips is proposed and determined so that the modified Ergun equation can also be used to calculate the flow resistance through wood chips. Suppose porous medium A] consists of uniform spherical particles and its flow resistance is given by the modified Ergun equation (2.23). By some mathematical manipulations, the following equation can be obtained, dP f.iA dL D eq 1-g V e ) + pB 1-e (2.31) v e J Substitutions of Vf = ef = s and sc = 1 - sf into (2.31) further results in, dP pA dL Dt f \ \£f J 1 D eq r \ \£f J •V (2.32) Now, if chip column A 2 has the same flow resistance as the porous medium A], Equation (2.27) and (2.32) will give identical results. Comparing these two expressions, we have, 25 Chapter 2 Mathematical Models of the Pulping Process juA ' ~ D2 eq R, = pB_ (2.33a) (2.33b) Theoretically, when the constants A and B are chosen, R] and R 2 are experimentally determined, both equation (2.33a) and (2.33b) can then be used to determine D e q . D e q is regarded as the equivalent diameter of chip column B. Macdonald [22] has found a correlation between A and D e q which works very well for Dudgeon's [32] and Gupte's [33] data with porosities ranging from 0.366 to 0.64. When D e q is expressed in meters, the correlation can be written as, A = 75.1 +1.148 xl05£> eq (2.34) Substituting (2.34) into (2.33a) and assuming a liquor viscosity of u=0.001kg/m-s, the equivalent diameter of the chip column is, Deq = 114.8 + 7114.82 + 0.3/?, 2R, (2.35) With the calculated D e q and known constant R 2, the constant B can be determined from equation (2.33b). The equivalent diameters of all the tested wood chips in Table 2.1 are calculated by (2.35) and listed on Table 2.2. Table 2.2 The equivalent diameters of tested wood chips Authors Wood material D e q, mm Reference Harkonen Scandinavian Pine 3.033 7 Wang Scandinavian Pine (Reference chips) 2.735 29 Wang Scandinavian Pine (New chips) 1.885 29 Lindqvist Scandinavian Pinus silvestrus, Mix 3 1.394 30 Lindqvist Scandinavian Pinus silvestrus, Mix 2 1.069 30 Lindqvist Scandinavian Pinus silvestrus, Mix 1 0.838 30 Lammi Scandinavian Birch 0.762 31 Lam mi Eucalytus Camaleulensis 0.488 31 26 Chapter 2 Mathematical Models of the Pulping Process As can be observed here, the equivalent diameters differ significantly, from 0.488 to 3.033 mm as seen in the Table; even for chips from the same wood species, the equivalent diameter could be quite different. For instance, for the two types of pine chips tested by Wang [29], the equivalent diameters are 1.885 and 2.735 mm respectively. This indicates that the size distribution of wood chips could significantly affect the equivalent diameter and therefore the flow resistance of wood chips. In Chapter 6, the effect of the equivalent diameter of wood chips on pulping performance is computationally investigated. 2.5 Special treatment of the chip phase 2.5.1 Non-Newtonian fluid The chip mass has been assumed to be a moving porous medium. The motion of the solid chip phase cannot be well described by a Newtonian fluid model, because some onsite observations indicate that the velocity profile of the chip phase is quite flat near the digester periphery. The chip phase is therefore assumed to be a non-Newtonian fluid, whose behavior can be generally described by the Cross model [27], UlZJL = (Ky)m (2.36) where r|o and r|oo refer to the asymptotic values of viscosity at very low and very high shear rates respectively, K is a constant parameter with the dimensions of time and m is a dimensionless constant. With some approximations, the Cross model can be reduced to some other popular viscosity models. For example, for i"|«r|o and r |»r |oo , the Cross model reduces to 27 Chapter 2 Mathematical Models of the Pulping Process V (Ky)" (2.37) With the definitions of K2 = and n = 1 - m , we further obtain, K V = 7o (Ky)" (2.38) which is the widely used power-law model and n is called the power-law index, K 2 is called the 'consistency'. By intuition, the chip phase internal structure should be rigid enough to prevent movement for values of shear stress less than a specific value, and when the stress exceeds this value, the chip phase will exhibit some shear rate like a Newtonian fluid. This characteristic can be best approximated by a Bingham plastic fluid model; the stress-strain relationships of a Bingham fluid and a Newtonian fluid are shown in Figure 2.5. at P re o Newtonian fluid Bingham fluid shear rate, y Figure 2.5 Flow curves for a Bingham plastic and Newtonian fluids The rheological equation for a Bingham plastic may be written T-T0=7jy (r>ry) (2.39) 28 Chapter 2 Mathematical Models of the Pulping Process However, an equivalent viscosity is desirable for computational purposes, that is, r = Vcffy (2.40) The Bingham plastic behavior given by equation (2.39) can be approximated by the Cross model, with a specific example given in Figure 2.6. 80 I \ — 60 t V in :\ [Pa 50 r \ , 40 - \ /• - V Viscosil 30 _ ' \ - / \ 4 \ Viscosil 20 10 Shear stress (Bingham plastic fluid) Shear Stress (Cross Model) Viscosity (Cross Model) .j i i I i i i u 0.8 H 0.4 H0.3 0.2 Ho.1 re CL in m 2 n CD 0.02 0.04 0.08 shear rate, y (s"1) Figure 2.6 Bingham plastic shear stress approximated by the Cross model [27] Olo=100.Pa-S, 7ra=l,K=150.0, m=1.05Pa) In some situations, to in Equation (2.39) may also change. To better describe the chip phase motion in a digester, Pougatch et al. [6] assumed that To is a function of both shear rate and chip pressure, with the following form, r0=Kpe(l-e-mf) (2.41) The effective viscosity is then expressed as, — = ri + r r -(\-e-'"r) (2.42) Similarly, the Bingham plastic behavior can be approximated by Pougatch et al's expression (2.42); a typical example is given in Figure 2.7. 29 Chapter 2 Mathematical Models of the Pulping Process 60 55 50 45 «T 40 i 35 £ 30 8 25 u * 20 15 10 5 0 0.02 0.04 0.06 0.08 0.10 shear rate, y (s 1) Figure 2.7 Bingham plastic shear stress approximated by Pougatch et al model [6] (110=1.5^-5^=100.0, K=0.001, Pc=500 Pa) In Chapter 4, a power law fluid flowing through round tubes is simulated and computed velocity profiles are compared with available analytical solutions. 2.5.2 Slippage on the wa l l boundary The no-slip condition always applies at wall boundaries for Newtonian fluid flow; in non-Newtonian fluid flow, however, relative slippage between the boundary and the fluid may exist, e.g. when squeezing tooth-paste in its tube. In Figure 2.8, the no-slip and slip conditions are illustrated at the fixed bottom wall. In Figure 2.8, the fluid velocity has a finite value USIJp at the bottom wall, which is different from that in 2.8. The slip condition is typically used for micro scale or "free" surfaces at macroscopic scale [34]. 30 Chapter 2 Mathematical Models of the Pulping Process (a) No-slip condition (b) Slip condition Figure 2.8 No-slip and slip conditions at the bottom wall To better represent the motion of the phase in the digester, slip conditions are used at both the downcomer and digester walls. A generalized slip flow model [34] is used here, sl'P wall v -L ^ vwall ~ -OX (2.43) wall where Aw is the velocity jump at the wall, L s is the slip length, w c is the velocity of the chip phase. Thompson and Troian [35] have experimentally observed that the slip length is a function of shear rate in the liquid flow at the solid surface. By analogy, we assume that the slip length for the chip phase flow at the digester wall is a function of chip pressure and has the following form: 4 = 4 + 1 a = i - 2 ( 1 ^ a r c f a n ^ J (2.44a) (2.44b) 71 where Ls is the characteristic length scale represented by the radius of the downcomer and the digester; P c is the chip pressure on the corresponding wall boundary in Pa; Aw is the surface area between the chip phase and the wall in m 2; a and (3 are dimensionless parameters. 31 Chapter 2 Mathematical Models of the Pulping Process From Equation (2.43) and (2.44), the slip velocity increases with the velocity gradient at the wall and decreases with the chip pressure at the wall. 2.6 Thermal properties of liquor 2.6.1 Density Liquor density depends on the solid concentration and temperature. A good approximation [21] for the density up to 50% dry solids is the following: where p is density of black liquor, kg/m3; T is the temperature, °C; and X is the dry solids concentration, kg dry solids/kg total suspension. 2.6.2 Viscosity The viscosity of black liquor (mixture of white liquor and dissolved solid material) depends on several factors, particularly the temperature and solids content. The source of the black liquor is also important. The following equations provide an estimate of the kinematic viscosity [21]: p25 = 997 + 649X (2.45) p/p25=l .008 - 0.2377 /1000 -1.94(7 /1000)2 (2.46) \nu = A + BIT (2.47) A = -2.4273 + axX + a2X2 + a3X3 (2.48) 73 = 6.1347 xlO 7 +blX + b2X2 +b,X3 where v is the kinematic viscosity, in mm2/s; T is the temperature in degrees, K. (2.49) Constants a,,&, are shown in Table 2.3. 32 Chapter 2 Mathematical Models of the Pulping Process Table 2.3 Coefficients for viscosity equations [21] Softwood Hardwood Tropical a, 9.1578 3.3532 10.482 a2 -56.723 3.7654 -54.046 a3 72.666 -2.4907 61.933 6. -42.178xl07 -5.442x107 -40.165xl0 7 b2 335.12xl0 7 21.915xl0 7 300.55xl0 7 h -349.23xl0 7 17.042xl07 -266.47x107 2.6.3 Specific heat capacity and thermal conductivity The following equations by Masse et al. [28] estimate the black liquor heat capacity, and thermal conductivity. cp = A.l 16(1 - X) + (1.675 + 3.31771000)X + (4.87 - 207 /1000)(1 - X)X3 (2.50) A = AH20(l-X) + aX + bX2 (2.51) <3 = 0.3176 + 2.268xl0" 3r (2.52) b =-1.394x 10"2 -3.069xl0" 37 (2.53) where C p is the black liquor specific heat, kJ/kg • °C; T is the temperature of the black liquor, °C; and X is the thermal conductivity, W/m°C. 2.6.4 Molecular diffusivity Development of continuous digesters has been focused on improvements in production capacity and pulp quality. The largest continuous digester currently in operation has a capacity of over 2000 tons per day. Obviously, the dimension of the digester increases with its capacity. With large dimensions, solid pressure/stress in the chip mass tends to be high, making the chip mass more compacted, and efficient diffusion becomes a key issue in practice. Therefore it is important to have good estimates of the diffusivity. 33 Chapter 2 Mathematical Models of the Pulping Process To estimate the molecular diffusivity o f alkali, McKibbins [16] measured the diffusivity of sodium in kraft cooked chips for a variety of temperatures and presented the following result: D0=3Ax\0-2TU2e-W70/(HT) (2.54) where Do is the molecular diffusivity, cm 2 /min; T is the temperature, K ; R is the gas constant (cal I mol • K). Because McKibb ins ' work was done with cooked chips, the result needs to be corrected with respect to the p H value and lignin content. Taking these into account, Gustafson [14] proposed the following expression for diffusivity: D0 = 5 . 7 x l O " 2 r l / 2 e " 4 8 7 0 / ( / " ' ) ( - 0 . 0 2 Z + 0.13[(9//f 5 5 +0.58) (2.55) where L is the lignin content in the over-dry wood, %; [OH] is the mole concentration of hydroxyl in the liquor solution, mol/L. 34 Chapter 3 The Kraft Pulping Reaction Model As introduced in the previous chapter, the whole pulping process consists of many sub-processes, including reaction between lignin, carbohydrates and cooking liquor, diffusion of active components in the cooking liquor (mainly OH" and HS") through wood cells, associated mass and heat transfer and etc. In chapter 2, the governing equations of the pulping process are described. This chapter will focus on discussing the pulping kinetics. 3.1 Representative reaction models Numerous pulping reaction models have been documented in literature. These models fall into the following three categories. 3.1.1 Single variable model A classical kraft pulping reaction model is based on the so-called H-factor, which combines temperature and time into a single variable representing the extent of the cooking. The delignification is assumed to be one single reaction, and the reaction rate k has the Arrhenius form, From some experimental data and setting a relative reaction rate k as 1 at 100°C, Vroom [11] gives, k = A-e EI(RT) (3.1) ln/t = 43.181 - 16113/T (3.2) And, defines the H-factor = \k{t)dt ={(43.181-16113/7^ (3.3) 35 Chapter 3 The Kraft Pulping Reaction Model Once the H-factor versus yield relationship is determined for a particular set of pulping conditions, the yield for a given H-factor will be known. This simple model finds widespread use in industry, and the H-factor is an extremely important operation parameter in modern process control. However, the H-factor only provides an accurate fit for data with similar cooking chemical charges and wood compositions. It does not include other important factors, like the concentration of alkali or sulfide. In addition, it does not account for any change of chip species. Nevertheless, a single variable method still has its attraction. Due to its simplicity, it serves widely in the control scheme of modern digesters. There are several other similar methods: Edwards and Norberg [35] developed an extension of the H-factor called the T factor, which combines the effects of time, temperature, initial effective alkali and liquor to wood ratio. Kubes et al. [37] developed a G-factor that works for any alkali pulping method. Paulonis and Krishnagopalan [38] discussed a sophisticated process control system for kraft batch digesters that uses liquid density, electrical conductivity, and other techniques in addition to the usual parameters in order to predict kappa number despite variations in wood chip supply and digester operation. 3.1.2 Three phase model: Gustafson's model In the second category of approaches, delignification reactions in kraft cooking are divided into three distinct phases: initial delignification, bulk delignification and residual delignification (Gullichsen [21]). The dissolving rates of lignin and carbohydrates are experimentally obtained for each phase. The following method describing the chemical kinetics has been proposed by Gustafson et al. [14]. Initial delignification: (Lignin>22.5%) dL dt ^ ( 1 7 . 5 - 8 7 6 0 / 7 ^ (3.4) 36 Chapter 3 The Kraft Pulping Reaction Model ^ = kic[OHfU^ dt lcV J dt Bulk delignification: (2.2%<Lignin<22.5%) (3.5) dL , ( 3 5 . 5 - 1 7 2 0 0 / ' / ' dt = k eK b\ rioH-y+kh2e(29A-mo° tT>[0H-]*[HS-]A L dC _ dL dt bc dt (3.6) (3.7) Residual delignification: (Lignin<2.2%) ]pH~f7L dL dt _ i (19.64-10804/77) - Krle dC _ , dL dt dt (3.8) (3.9) In the above equations, L is the lignin content, C is the carbohydrates content, [OH~] is the hydroxyl ion concentration, and [HS~] is the hydrosulfide ion concentration. The parameters ka , kbx , kb2, kr/, kic, kbc and krc are species related constants. For Pinus Silvestris (lignin 27%, carbohydrates 67%, extractives and ash 6%), these constants are listed in Table 3.1. Table 3.1 Constants fitted to pulping reactions [14] Phase Equation Constant Value Delignification Reactions Initial (3.4) ka 1.0 Bulk (3.6) kb\ 0.15 (3.6) kb2 1.65 Residual (3.8) Ki 2.2 Carbohydrate Dissolution Initial (3.5) kic 2.53 Bulk (3.7) kbc 0.47 Residual (3.9) krc 2.19 37 Chapter 3 The Kraft Pulping Reaction Mode l 3.1.3 "Two lignin" model Smith and Williams [39,40] proposed a comprehensive kinetic model for the conventional kraft process. This model, often referred to as the Purdue model, has been a basis for several other digester studies. In this model, the wood components are high- and low- reactive lignin (also called fast and slow lignin), cellulose and the hemi-cellulose components galactoglucomannan and araboxylan. The cooking liquor consists of sulphide and effective alkali. The reaction rates of all these components are described by first order equations. This model considers only the initial and bulk stages of deliginification. One should also mention that the two classes of lignin have different reaction rates, making the calculations more complicated. 3.2 Motivation for deriving a new reaction model Unfortunately, all existing reaction models have some drawbacks. For instance, computations using the three-stage model given by Gustafson etc. [14] exhibit very good agreement with the experimental data. However, the calculation of alkali reaction rate requires the knowledge of the degradation rate of the acetyl group contained in wood chips and the transition points characterizing the reaction stages may change unexpectedly. Therefore, work is continuing on developing an improved kinetic model. Recognizing the effect of sulfide concentrations on the delignification rate, Aurora Santos et al. [41] established a delignification model where hydrosulfide concentration is included. This model is derived for Eucalyptus globulus; when species or operating conditions change, a so-called "stoichiometric coefficient" needs to be determined. Ling Yang et al. [42] presented a mechanistic pulping model where three key kinetic steps are considered: adsorption of hydroxide and hydrosulfide ions on the lignin; chemical reaction on the solid surface to produce degraded lignin products; and desorption of degradation products from the solid surface. This is physically more acceptable. However, the empirical constants in the model are also dependent on the operating conditions and need to be determined case by case, narrowing its application in some situations. 38 Chapter 3 The Kraft Pulping Reaction Model Obviously, further CFD work is still needed to develop an easily applicable, more reliable and accurate pulping reaction model. From the computational modeling point of view, a good reaction model should give the reaction rate for each involved species, with only local chemical and physical conditions specified. In continuous cooking, this is more important because both temperature and chemical concentrations may exhibit significant spatial variation. A common drawback of most pulping reaction models is that they all need some modifications to fit pulping processes under different conditions. It is outside the scope of this research to address issues of fundamental chemical kinetics. Instead, we will take advantage of published work and use the available experimental data to derive an improved and easier to apply reaction model by regression analysis. As will be shown, this model has the capability of predicting pulping kinetics in various situations with good accuracy, although it does not give much physical insight. 3.3 Regression method The following considerations have been taken into account in deriving the general pulping model. 1. Take full advantage of the available experimental data So far, all the existing empirical pulping models were based on one set of specific experimental data. This is probably part of the reason why these models need appropriate modifications when used in other situations. An expanded database of experiments can better reflect the physics and lead to a more extensively applicable kinetic pulping model. In view of this, a relatively large number of data were used in the model regression, with cited data including Aurora Santos [41], Aurell et al. [44], Bugajer [45], Matthews [47], Kleinert [48], Wilder and Daleski [49,50] and Chiang et al. [51]. 2. Select suitable regression variables Prior to discussion, a few definitions are introduced as follows. Residual lignin content is defined as the mass percentage of undissolved lignin in wood chips per unit mass of wood chips (oven-dry basis) in a macroscopic view. The 39 Chapter 3 The Kraft Pulping Reaction Model mass fraction of residual lignin is defined as the percentage of residual lignin over the total lignin, which consists of dissolved and undissolved lignin. The residual carbohydroxide and its mass fraction are defined similarly. Delignification is a complex reaction process. As it proceeds, the dissolved solid material will gradually diffuse into the surrounding cooking liquor. The relatively slow diffusion on the other hand will result in a high concentration of dissolved lignin near the chip surface, which retards further removal of lignin from solid chip surfaces. Therefore, even though residual lignin contents are the same, the one with less dissolved lignin will host a faster lignin reaction. From this point, it is more appropriate to express the delignification rate in terms of the mass fraction of residual lignin. Likewise, the mass fraction of the residual carbohydrates has a closer connection to carbohydrate degradation. Therefore, the mass fractions of residual lignin and carbohydrates, instead of the traditionally used residual lignin and carbohydrate contents are chosen as regression variables. Hydrosulfide concentration is another regression variable as its influence on delignification was documented by Bugajer [45]. 3. Avoid using global variables In some pulping models, global variables such as liquor to wood ratio were used. However, chemical reactions always occur locally; therefore, they are dominated by the local temperature and chemical concentrations. In large industrial digesters, both temperature and chemical concentrations vary from place to place; the fixed quantity liquor to wood ratio may not be a good representation. Further, in continuous cooking, liquor to wood ratio is hardly used as a parameter. 4. Provide sufficient data to ensure model accuracy In normal operating conditions, delignification proceeds quite rapidly/efficiently in the bulk stage but slowly in initial and residual stages. This makes the lignin content versus time a very sharp curve in the middle of the range, while quite flat on both ends. 40 Chapter 3 The Kraft Pulping Reaction Model To characterize this feature in regression, more experimental data are needed, but in practice, only limited measurements are taken. It is therefore very hard to obtain an accurate regression model. Because of this difficulty, all the experimental data to be used in regression are fitted by piece-wise functions, from which a desired number of data are then extracted. This way, a reasonably large number of data are generated for each experiment. Thus the whole database is expanded. This strategy turned out to be very helpful in leading to a higher accuracy regression model. 5. Abandon the "three stage" regression form An extensive review given by Santos [41] showed that the majority of the pulping models divide the delignification process into three distinct stages. Accordingly, three different expressions are used to calculate the delignification rates. In these "three stage" models, the transitional points are considered to be constant. However experiments indicate that these transitional points may change in different cooking conditions. In this work, the "three stage" idea will not be used in order to eliminate the need for transitional points. 3.4 Regression procedures and the pulping model In typical cooking conditions, Santos [41] has experimentally showed that both the effective alkali and sulfide concentrations appear as only functions of the mass fraction of residual lignin content, as shown in Figure 3.1. If we check some previous pulping models, the carbohydrate degradation rate can also be expressed solely by the rate of lignin removal. This is quite encouraging. It clearly indicates that once the delignification rate is known, the calculation of the concentrations of other important species becomes fairly simple. The accuracy of the model prediction therefore depends to a large degree on the computation of the deliginification rate. These experimental results also imply that 41 Chapter 3 The Kraft Pulping Reaction Model the mass fraction of residual lignin is a better variable than the residual lignin for characterizing the reaction process as mentioned earlier. Functions of effective alkali and sulfide concentrations related to mass fraction of residual lignin are obtained using Santos's experimental data. They can be well represented by polynomials: [OH] = 0.4429 + 3.262Zr -17.19L) +42.3SL3r -47.49/4 + 20.04 +{OH\ -1-40 [SH] = 0.04765 + 0.15834 -0.3236Z2. +0.2616Z? + [SH~\ -0.144 where Lr = LI L0, represents the mass fraction of the residual lignin. (3.10) (3.11) The overall F statistics are significant, with F=4.04, p<0.00035 and R-square=0.918 for equation (3.10) and F=8.23, /?<0.00004 and R-square =0.827 for equation (3.11). Calculated concentration profiles are shown in Figure 3.2. 1 p. 1.20 \-0.60 0.20, qp -fc x+ V Exp. Case 1 x Exp. Case 2 • Exp. Case 3 O Exp. Case 4 * Exp. Case 5 • Exp. Case 6 A Exp. Case 7 O Exp. Case 8 + Exp. Case 13 t> Exp. Case 14 v C* x > • *<& 4 0.2 0.4 0.16 _l 0.12 h o E I 0.08 h • A o Exp. Case 1 Exp. Case 2 Exp. Case 3 Exp. Case 4 Exp. Case 5 Exp. Case 6 Exp. Case 7 Exp. Case 8 * A A * 0.0 _l I I L-J I I I I I I I L-0.2 Lr 0.4 0.6 Lr (a) Effective Alkali (b) Sulfide Figure 3.1 Effective alkali and sulfide consumptions versus lignin removal (liquor to wood ratio h=5) [41] 42 Chapter 3 The Kraft Pulping Reaction Mode l E I 1.60 1.40 1.20 0.80 0.40 Exp. Case 1 Exp. Case 2 Exp. Case 3 Exp. Case 4 Exp. Case 5 Exp. Case 6 Exp. Case 7 Exp. Case 8 Exp. Case 13 Exp. Case 14 Polynomial fit R!=0.918 0.12 I Exp. Case 1 Exp. Case 2 Exp. Case 3 Exp. Case 4 Exp. Case 5 Exp. Case 6 Exp. Case 7 Exp. Case 8 Polynomial fit, R!=0.827 0.6 Lr (a) Effective Alkali (b) Sulfide Figure 3.2 Fitted effective alkali and sulfide profiles [41] 0.8 1.0 After developing the correlations between effective alkali, sulfide and the mass fraction of residual lignin, efforts were directed towards obtaining the delignification rate. From theory of chemical kinetics [45], we assume the delignification rate can be given in the following form: dL dt Or equivalently = -ec^ir[OHf[SHfL In dLr dt = C 0 -jCy + \n[OH\C2 + \n[SH]-C3 + ln(4) • C 4 (3.12) (3.13) When sufficient data sets of ln dLr dt ~ , m[OH], \n[SH] and ln(Zr) are provided, the constants Ct, i- 0,1,2,3,4 can be determined by multi-variable linear regression. In typical cooking experiments, the lignin content and temperature at different times are measured, but measurements are usually not taken that often. As a remedy, the continuous form of function Lr(i) in each individual cook can be fitted with the measured data, then ln dLr dt and ln(Zr)can be calculated at any desired time, with ln[0//] and 43 Chapter 3 The Kraft Pulping Reaction Model ln[57/] being determined similarly. As the delignification rate can be affected substantially by heating conditions, cited experimental data are classified as normal heating and fast heating data. Although there is no specific rule, we suggest cases with heating time less than 80 minutes belong to the fast heating group, while the rest belong to the normal heating group. Therefore, data from Bugajer [45] and Aurora Santos [41] fall into the first category; data from Aurell [44], Mathews [47] and Wilder belong to the latter one, and corresponding conditions are listed on Table 3.2. Accordingly, two different regression expressions are obtained for the fast heating and normal heating conditions respectively. dLr dt 17.19-7201 /7'[O/7]-,041[5/7]2 0 7 5 Z° 9 2 5 (Fast heating condition) (3.14) dl\ = _ e 8 6 „ - 5 6 6 3 / r ^ 1 0 373 ^ ^ 0 , 4 8 ^ 0 . 7 6 , c o n d i t i ( m ) (3.15) dt R-squares for these two expressions are 0.927 and 0.965, indicating successful regressions. Table 3.2 Cooking conditions of cited experimental data3 Reference Case L0 Co E A as Sulfidity L / W [NaOH] [Na 2S] th tc 7 C % % [NaOH] mol/L as [SH], mol/L g/L g/L min min °C Aurell 1965 1 27.3 67.7 0.937 0.134 4 32.14 10.44 120 120 170 Aurell, 1965 2 27.3 67.7 1.562 0.223 4 53.57 17.40 120 120 170 Matthews, 27.2 66.5 0.92 0.114 5 32.25 8.86 90 90 170 1974 Wilder, 1965 28.6 71.4 1.50 0.375 200 45.0 24.0 - 205 150 Bugajer, 1984 1 28.77 71.23 1.129 0.161 4 38.71 12.58 65 100 159 Bugajer, 1984 2 28.77 71.23 1.129 0.161 4 38.71 12.58 65 100 168 Santos, 1997 Cl -4 22.1 73.6 1.20 0.18 5 40.80 14.04 45 150 150,160, 170,180 Santos, 1997 C9 22.1 73.6 1.8 0.18 5 64.80 14.04 45 150 160 Santos, 1997 C10 22.1 73.6 0.7 0.18 5 20.80 14.04 45 150 160 Santos, 1997 C13 22.1 73.6 1.2 0.43 5 30.79 33.54 45 150 160 Santos, 1997 C14 22.1 73.6 1.2 0.07 5 45.20 5.46 45 150 160 Chiang etc., 26.8 70.0 1.193 0.211 5 39.31 16.42 100 120 170 1988b Kleinert, C l -3 28.0 72.0 1.226 0.216 10 40.38 16.87 - 60 170,175,180 1966b Note: " 7c, cooking temperature; th heating time; tz cooking time. b Data not used in the model regression, but will be used to check the model prediction. 44 Chapter 3 The Kraft Pulping Reaction Model Unfortunately, experimental data of carbohydrates degradation are not well documented. The degradation rate of carbohydrates given by Aurora [41] is used directly here: dCr — - = < dt 0 . 2 5 ^ while Lr < 0.4 dl ( 3 - 1 6 ) 0.14—^ while L > 0.4 dt From Equation (3.10), the reaction rate of effective alkali (EA) can be further derived. Because the mole concentrations of hydroxyl and sulfide appear in the expression of lignin reaction rate (3.14) and (3.15), the EA reaction rate in its mole concentration form is considered here. To be consistent with Santo's expression (3.10), EA is defined by, EA = (NaOH + 0.5Na2S)as N a 0 H g/L (3.17) Note: In North America, EA is defined as (NaOH + 0.5Na2S)as N a l o g/L. The mole concentration of EA in the free liquor is then, FA FA [EA] = -=^— = — mol/L (3.18) W 40 Also, the mole concentration of EA is essentially the mole fraction of [OH] with the definitions of (3.17) and (3.18), i.e., [EA] = [OH] ( 3 1 9 ) Manipulations of Equations (3.10) and (3.19) leads to, d\EA\ = dyOH] = / _ + 2 _ 1 Q 0 Z 3 + 1 0 0 Z 4 ) ^ ( 3 21) dt dt V r r J dt Similarly, if we define, Sulf = Sulfidity = Na2S NaOH + Na2S ax NaOH (3.22) The derived reaction rate of sulfidity will be, dSu£ = d[HS] f + 8 ^ ( 3 2 3 ) dt dt v ' dt 45 Chapter 3 The Kraft Pulping Reaction Model 3.5 Model validation With Equations 3.12, 3.16, 3.21 and 3.23, all of the four species Lr, C r, OH and SH can be computationally determined. Naturally, the accuracy of the model prediction is a major concern. In view of this, a computational code was developed to simulate the batch cooking experiments. A four-stage Runge-Kutta time advance scheme [43] was used to solve the four differential equations 3.12, 3.16, 3.21 and 3.23. Time steps are reasonably small and vary case by case; measured temperatures are directly used in the computation. As discussed before, the prediction of lignin content is of great importance to the model accuracy. Also, concentration measurements of the other species are not available, and therefore only the lignin content will be checked here. Figure 3.3, 3.4 and 3.5 show the variation of lignin content as time proceeds. Al l lines represent model predictions, while symbols represent experimental measurements, by Aurell et al. [44] and Bugajer [45] in figure 3.3 and Wilder [49,50], Matthews [47] in figure 3.4. As can be seen, fairly good agreement is achieved generally, but the model underestimates Aurell's experiments in the later cooking period. Nevertheless, the biggest predicting error of lignin content on wood is about 11%. 30 • Aurell case 1 • Aurell case 2 • Bugajer case 1 O Bugajer case 2 Prediction A.1 Prediction A.2 Prediction B.1 Prediction B.2 200 250 : 300 time, min Figure 3.3 Lignin profiles by model predictions and experiments Experiments by Aurell et al. [44] and Bugajer [45] 46 Chapter 3 The Kraft Pulping Reaction Model 35 TJ O O c o c CI) Mathews Wilder and Daleski Prediction Mathews — Prediction Wilder _! i i i I i i i i L 50 100 150 200 time, min Figure 3.4 Lignin profiles by model predictions and experiments Experiments by Wilder [49,50] and Matthews [47] •a o o 5 c 'E • o Kleinert170°C Kleinert175°C Kleinert180°C A Chiang etal. Prediction L.170 Prediction L.175 Prediction L.180 N Prediction Chiang \ A \ \ \ A \ \ A \ N A ^ , 100 150 time, min Figure 3.5 Model prediction versus experiments by Kleinert [48] and Chiang et al. [51] In Figure 3.3 and 3.4, however, all the experimental data have been used in the model regression. To further check the quality of the model predictions, some other cases 47 Chapter 3 The Kraft Pulping Reaction Model not previously used in the model regression were modeled and the simulated results are given in Figure 3.5. Again, the model gives very good prediction to Chiang's experiments, with the error mostly limited to 1%. For experiments done by Kleinert [48], the simulations appear promising in the early cooking stage, but under-estimate the measurements thereafter, with the biggest errors in the range of 2-3% occurring in the bulk cooking period. As listed in Table 1, Kleinert's experiments were conducted with high liquid-to-wood ratio and relatively higher cooking temperature, this is probably the reason why model predictions give relatively larger errors. 3.6 Conclusion Based on extensive experimental data, a kinetic pulping model is obtained exclusively by mathematical regression. In this model, two different expressions are used to calculate the delignification rate respectively for normal heating and fast heating conditions; Santos equation of carbohydrates degradation rate is directly employed because of lack of experimental data. By introducing the mass fractions of lignin and carbohydrates as regression variables, better utilizing the experimental results and choosing suitable regression form, the resulting pulping model predicts cooking experiments conducted under various conditions by different investigators, without the need to determine some operating conditions related model constants. In addition, the model is not based on the "three cooking stages", thus no transitional points are needed in the computation. Simulations using this model exhibit good agreement with the experimental data, with prediction errors less than 11% for all the tested cases. It is also expected that the model can be used in the simulation of continuous cooking. 48 Chapter 4 Computational Methods In Chapter 2, the governing equations describing the two-phase flow, energy conservation and chemical reactions in the cooking process have been introduced. In this chapter, the general discretization schemes and computational procedures will be presented. As discussed in previous chapters, the solid chip phase is treated as a fluid, thus Navier-Stokes equations apply for both phases. Instead of giving the detailed discretization techniques for each individual equation, we will first give the approach to solve the single phase Navier-Stokes equations, and then address some additional issues in order to model the two phase flow problem. Examples of imposing typical boundary conditions are also provided in the later part of this chapter. 4.1 Solving Navier-Stokes equations by SIMPLE method In compressible flows, the continuity equation can be used to determine the density and the pressure is calculated from an equation of state. Unfortunately, this approach is generally not appropriate for incompressible flows or low Mach number flows, as density remains constant. By adding a time derivative of pressure into the continuity equation, the compressible computation method can be used in incompressible situations. This is called artificial compressibility method which was proposed by Chorin [52]. Usually, this method is used in conjunction with the collocated variable arrangement, and the resulting pressure oscillations can be lowered numerically. However, for fast convergence, this method requires a good compressibility parameter, which varies significantly from case to case. Selecting such a parameter is not always easy, and that proves to be a major deficiency of the artificial compressibility method. Methods based on the SIMPLE algorithm [53] are fairly efficient for solving both steady and unsteady problems, and this kind of method has been used widely in incompressible flow computations. There are 49 Chapter 4 Computational Methods many variants of this method available; the SIMPLE and SIMPLEC methods have been implemented in the code. In this section, the basic SIMPLE method will be explained. 4.1.1 Navier-Stokes equations We consider a 2-D, steady, laminar, incompressible flows of a Newtonian fluid. In a Cartesian coordinate system, the conservation laws for mass and momentum can be written as: 8(pu) | d(pv) | d(pw) = Q dx dy dz d ( \ d , \ d , x dP — (puu) + — (pvu) + — (pwu) = - — dx dy dz dx d ( \ d ( \ d ( N dP — (puv) + — (pvv) + — (pwv)=- — dx dy dz dy dx dy dz dz dzrr dvyx dz, dx dx dx y \JJ\, KJJ\, KJJ\. J ^dzr„ dz„„ dz ^ xy + - yy + - zy • + dy dy dy + PBX + pBv dz„ | dzyz | dz. dz dz dz + pBz (4.1) (4.2a) (4.2b) (4.2c) For incompressible flow, we have, du XX * o dx ryy=2p Vr- = 2// 5v dy dw ~dz xv vx A* du dv^ dy dx X — Z = u vz zv r* *:x=*x:=P dv dw dz dy 'du dws ydz dx j (4.3) With these relations, the momentum equation can be arranged into its velocity form: ~ (puu) + -|- (pvu) + -|- {pwu) = - ^- + p dx dy dz dx d , \ d i \ d ( \ dP — (puv) + — (pvv)+—(pwv)=- — + p dx dy dz dy d dx du dx dy dv du dx dy dx du dv dy dx f + • dy dv_ dz dw ^du V dx dz j + • dz dv dw dz dy + pBy 50 Chapter 4 Computational Methods d ( \ d , \ d i x dP — (puw) + — (pvw) + —(pww)=- — + p ox dy dz dz dx dw du — + — dx dz + • dy dv dw dz dy + dz V dz J (4.4) + After some mathematical manipulation, the momentum equations can be further arranged into the following uniform form: V -{pvt- pV</>)=-VP + S, (4.5a) In Cartesian coordinates, (4.5a) can be written as, dx dx ^ d + — PV0-JU d0_ dy dz pwcj) - JU d£ dz dP_ dx, where xt• = x, y, z; (/> = u, v, w r dx d du M — v dxj f du^ M — v dy J d ( du^ M — ^ dz j + • dy dv' M — dx + • dx + • dy dv) dy dx + d ( dv^ dy V dz dz d_ dz d_ dz dw dx dw M — v dy j V dz + pBx when 0 = u + pBy when 0 = v + pBz when <j> = w (4.5b) (4.6) In case of constant viscosity p, the viscous terms in S4, are eliminated by continuity. 4.1.2 Collocated grids While solving the Navier-Stokes equations, one can store the velocity variables and pressure at the same grid points. Such a grid system is called a collocated arrangement. As many of the terms in the transport equations are identical, the number of discretization coefficients is minimized and the programming is simplified. For problems with more variables or a more complex geometry, or when the multigrid approach is used, the advantages of a collocated arrangement become more apparent. However, when used in incompressible flows, the collocated arrangement has difficulty handling the pressure-velocity coupling resulting in pressure oscillations. On the other hand, the staggered arrangement, first introduced by Harlow and Welsh [54], offers a strong coupling 51 Chapter 4 Computational Methods between the velocities and the pressure, and thus successfully solves the incompressible Navier-Stokes equations without oscillations. It is for this reason that the collocated arrangement was seldom used in the early days of CFD development. In the 1980s, improved algorithms were developed to solve the pressure and velocity equations, and the collocated arrangement began to win popularity. Today, it is adopted in many commercial codes, such as FLUENT. In the present code, the collocated variable arrangement is used based on the scheme of Ferziger and Peric [55,56]. In collocated grids, only one set of control volumes is required, the computational points (storage locations for variables) are then located in the center of the control volumes, an example of such arrangement in a 2-dimensional case is given in Figure 4.1. o °N l i i o o W n • m i l A w ; p e s ; se o E o i i *S o Figure 4.1 Collocated grids in 2-D situation In the collocated arrangement, careful consideration should be given to the pressure gradient calculation. Ferziger and Peric [55] showed that if pressure derivatives at cell faces are evaluated by interpolation of cell center derivatives, unrealistic pressure distributions could be produced. The main reason is the use of the 2Ax approximation for the first derivatives. To cure this problem, they suggested an effective scheme where pressure derivatives at cell faces are evaluated by central differences and Ax spacing, thus introducing the need for interpolation. 4.1.3 Discretization of the equations by control volume method 52 Chapter 4 Computational Methods For simplicity, the discretization scheme is described below only for a two-dimensional situation, although the code is developed three-dimensionally. Finite Volume Method starts from the conservation equation in integral form, based on equation (4.5a). The integral form of the generic conservation equation can be written as, jpv</> • MS = jfNj • MS + JVPJQ + J S ^ Q (4.7) s s n n Integrating the equations (4.7) over the control volume leads to the following discrete equation: 4 - 4 + 4-4= j f dP ^ AV dx, 4 W = SAV (4.8) where the I's represent the average fluxes through the control volume faces in the direction of increasing coordinates, and SAV is the volume integral of the right hand side of Equation (4.5a). The terms appearing in Equation (4.8) have to be evaluated differently for the mass and momentum conservation equations; this will be outlined below. For convenience, we define Cx = (PUAX) C„ = (PUAXI C„v = (pUAx )w (4.9) Cy=(pVAy) Cyn=(pVAy\ Cys=(pVAy\ (4.10) Performing integration over the control volume for equations (4.5b) yields the discrete equations. 4.1.3.1 Continuity equation In the continuity equation (4.1), the fluxes are the mass flow rates through the control volume faces, and are evaluated as: Ie = Cxe = {pUAx\ In = Cyn = (pVAy\ (4.11) The discrete continuity equation is then pUeAx - PUWAX + pVn Ay - pVs Ay=0 (4.12) 4.1.3.2 Momentum equations 53 Chapter 4 Computational Methods For simplicity, the convection fluxes are approximated with upwind difference scheme (UDS), while the diffusive part is estimated by central difference scheme: I" = r U -I" =c u ju — V Sx f A ^ 8y (UE-UP) = CxeUe-De(UE-UP) (UN-UP) = CynUn-Dn(UN-UP) The pressure gradients are also evaluated with central difference approximations. CxeUe-De(UE-UP)]-[C^UW-DW(UP-UIV)] + CynUn-Dn(UN-UP)]-[CysUs-Ds(UP-Us)] = -Ax^PE-Pw) where cell face velocities are evaluated by upwind scheme, for instance, \UP, if Cxe>0 u„ = u if Cxe<0 (4.13) (4.14) (4.15) (4.16) After some rearrangements and applying Successive Over-Relaxation (SOR) scheme [58], the discrete U-momentum equation can be written as: apUP = Ia„ 6£/„* + (1 - au)U°P - Ax,P(Pe - PJ nb where the index "nb" runs over all neighboring points E, W, N and S; <*E = A + maxC-C^.O) aw =Dw+max{Cxw,0) ®N = D n + max(-C>,0) as = Ds + max(C>,0) ar=-rYuan nb (4.17) (4.18a) (4.18b) (4.18c) (4.18d) (4.18e) where ocu is the over-relaxation parameter, having a typical value of 0.7. Similarly, the V-momentum equation is, apVP = I anbVnh +{l-au )V"P - Ay,P (Pn - Ps) nb The coefficients in equation 4.19 are similar with 4.17, as listed in formula 4.18. (4.19) 54 Chapter 4 Computational Methods 4.1.3.3 Pressure-velocity-coupling The U- and V-momentum equations are solved using guessed values for the pressure and mass fluxes. The velocity components U* and V* calculated with these guessed values will in general not satisfy the continuity equation, thus: P U e A - PUWAX + pVn Ay - pVs Ay = Sm The discretized U-momentum equations for nodes P and E can be written as: dp I nb u;= — X anbU*nb +(l-au )U°P - AXtP (P*e - Pw*) = UP+UP (4.20) (4.21a) UE = 1 X anbU*nb +(l-au )U°P - AX,E (P*e - P*) nb = UE+UE (4.21b) The cell face value is obtained by linear interpolation of the terms in equations 4.2 with the exception of the pressure difference, which is evaluated as in the staggered approach. This yields: ^ l u \ a p j 2ZanbUnb+(l-au)U nb f 1 ^ AxJPE-PP)e=Ue+Ue (4.22) With equation 4.21, we have U] f 1 ^ \ A P JE X^t/;+(i-«„)t/° nb j£ ' 1 ^ \ Q P J (4.23a) Uf r 1 ^ \ A P Jplnb T^bKb+a-ccM f 1 ^ KAP Jp Axp(Pe Pw)p (4.23b) So, U = average V a. u: a. Jp. K A P J dx \oxjE \Sxjp = U + A x 4 v f A~A A f Sp^ Sx K ap J (4.24a) 55 Chapter 4 Computational Methods U„ =U„ (Ps-Pp)e=u:-rAx-Ax^ v P J 'Sp) (Sp = Ue+Ue (4.24b) where the overbar indicates linear interpolation. Neglecting the interpolated pressure gradient in equation 4.22 and 4.24b, the velocity correction Ue is linked to pressure by the so-called velocity correction equation, Ax • Ar V a P J Sp Ax • A, \ S x J e \ Q P J Sp_ ' 1 N \ap j (4.25) Substituting equation 4.24b into equation 4.20, and imposing mass conservation, results in a pressure correction equation: a,,PP = Y u a ^ p n b - S m nb with the coefficients defined by, aE = pA ' 1 ^ x,e \ a P J f 1 > as = pA y,s \ Q P J S ap=Y^anb r , A aw = PA AN = PAy,n KQP J f 1 A (4.26) (4.27) (4.28) (4.29) nb As can be seen, the SIMPLE algorithm with collocated arrangements is very similar to the staggered approach, but requires some interpolation. After discretization, the differential equations are converted to a system of algebraic equations. For structured grids, the resulting algebraic equations always have very sparse coefficient matrices and they can be efficiently solved by some iteration methods, such as strongly implicit procedure (SIP) proposed by Stone [56] and alternating direction implicit (ADI) method introduced by Hageman and Young [58]. 56 Chapter 4 Computational Methods 4.1.4 Solution procedures for steady-state problems 1. Start the calculation of the fields at new iteration step using the latest solution u" and p" as starting estimates for u"+l and p"^ . 2. Assemble and solve the algebraic equation systems for the velocities to obtain uj . 3. Assemble and solve the pressure-correction equation to obtain p . 4. Correct the velocities and pressure to obtain the velocity field w,, which satisfies the continuity equation, and the new pressure p 5. Return to step 1 and repeat, using updated w, and p as improved estimates for u"+x and , until all corrections are negligibly small. 4.1.5 Curvilinear grids In case of curvilinear grids shown in Figure 4.2, the discretization schemes and the approach solving Navier-Stokes equation introduced above still apply, but some modifications are needed when evaluating the fluxes, pressure term, and some source terms. The discretization of corresponding terms in the liquor phase equations is outlined in this section. Figure 4.2 Two-dimensional curvilinear grids 4.1.5.1 Mass flux in continuity equation On the east face of a control volume, the mass flux is evaluated as: 57 Chapter 4 Computational Methods me = }pV • nedS (4.30) The mid-point rule approximation of the above integral leads to, me « (pV •n)eAe= (puAx\ + (pvAy \ 4.1.5.2 Convective fluxes in momentum equation (4.31) Having evaluated the mass flux in continuity equation (4.1), we can further evaluate the convective flux term in the momentum equation (4.5a). On the east face of a control volume, the convective flux integral is, Fec = \p</>V • MS « (p<t>V • n\Ae = mje (4.32) The computation of F° varies in different discretization schemes. For example, in explicit central difference scheme (CDS), m e and <j)e are estimated by interpolation with the values from the last iteration; in implicit upwind difference scheme (UDS), however, Ff is formulated as, F. c,implicit max(me90)up +min(me,0)uE (4.33) 4.1.5.3 Diffusive fluxes in momentum equation As an illustration, we will again evaluate the diffusive flux integral on the east face of a control volume. By using central difference scheme (CDS), the explicit and implicit diffusive flux integral are calculated as, pd^ucu = jMs70.MS*pe Se d±Ax + d±A, dx x dy y V LEP J % -nAe={pA)e(§ -n ?J <t>E-<t>F L (4.34) (4.35) EP J 58 Chapter 4 Computational Methods When the line connecting nodes P and E is orthogonal to the cell face, E, • n = 1; if not, this discrepancy will be corrected by deferred correction, which will be described shortly. 4.1.5.4 Approximation of pressure and source terms The pressure can be treated as conservative forces on CV surfaces, applied with Gauss theory, the pressure integral term can be calculated as, Fp = - J — dn = $P-ndS*YPc4? > (x,=x,y,z) (4.36) n s c On the other hand, the pressure can also be regarded as a body force, so we can compute it another way, Fp=-\^dQ J Ply ndxi r8P) AQ (4.37) Generally speaking, the first approach is fully conservative, but the second does not always remain so. The source term is basically a volume integral. The midpoint rule is commonly used to approximate it, which approximately gives second order of accuracy. F s ' = - { 5 ^ Q « ^ p A Q (4.38) 4.1.5.5 Deferred correction scheme As discussed earlier the adoption of collocated grids produces an oscillatory solution that is not real. A remedy for this is the so-called deferred correction scheme with appropriate discretization schemes; in deferred correction schemes, both convective and diffusive fluxes are computed in the following way: F = Fimp + [Fexp - F'"'p]M (4-39) 59 Chapter 4 Computational Methods where 'imp' and 'exp' denote the implicit and the explicit flux approximation, 'old' means value from previous iteration. In Muzaferija's approach [59], all the fluxes are evaluated by CDS scheme except the implicit convective flux, which is computed by UDS scheme. Then the approximation of convective and diffusive flux integrals at the east face of a control volume is written as: F° = max(me,O)0P + min^,0)^ / ;. + \me</>e -max(mc,0)^p -min(me,0)(/)E]" + ("4 = (4 EP J \ (d^ d£ old </>E-0 V I'EP J (4.40) (4.41) 4.2 Additional treatment of two-phase flow equations 4.2.1 Evaluate the mass flux In continuous digesters, the free liquor and chip phases occupy each individual control volume and flow through the same control volume faces. To evaluate the mass fluxes through the control volume faces in the continuity equations (2.7,2.8), the volume fraction of each phase has to be taken into account; for instance, the liquor phase mass flux through the east face of a control volume is approximated as: me « {pjSjV • n)Ae = (pfefufAx) + (pfefvfAy) 4.2.2 Discretization of the inter-phase interaction terms (4.42) In the presence of wood chips, the liquor phase momentum equation (see equation 2.11) has two extra terms, representing the momentum exchanges caused by the inter-phase friction force and mass transfer respectively. Jt (pf£fufj) + (pf£fuf,,ufj) = £f ' d P f + d r ^ v dx + EfPfgi + FfJ + mnucj (4.43) J J 60 Chapter 4 Computational Methods where mu is determined by chemical kinetics andF ;. is given in Equation (2.15), f s2 , , 4 — + R2£c\UcJ -"/,/ v 4 (4.44) The integral of these two terms over the control volume can be evaluated linearly, which gives, \(F/j + m]2uci)dQ « {FfJ + mnuci)pAQ = (Sp • ufJ + S0)AQ (4.45) where 1 ( s2 c £f V £f 1 ( s2 R\ c £f £f v 2 c c "cy " / , / (4.46) (4.47) The extra terms in the momentum equations of the solid chip phase are treated similarly. 4.3 Discretization of scalar equations The energy and species equations are the scalar equations in our modeling. Al l the scalar equations are solved after the flow fields for both phases are obtained. The discretization methods of these scalar equations are similar; as an example, the discretization of the steady state energy equation is outlined below. Referring to equation 2.17, and neglecting the heat generated by work and viscous dissipation, the steady state energy equation becomes, dx -(PfU/JT) = 8X : Af dT \CPJ dx.iJ + Q (4.48) Using the deferred correction scheme to approximate both left hand side convective and right hand side diffusive term gives, 61 Chapter 4 Computational Methods p c _ pcjmp _|_ j ^ c . e x p _ pc,impl^U p i t _ pdjmp _|_ ^ ( / , e x p _ pdjmp j ° (4.49) (4.50) On the east face of a control volume, fluxes evaluated by explicit and implicit schemes are, (4.51) F e c " e x p = mX xp F d.exp _ CPJ ydx j Aex + A ey F°Jmp = max(me,0)Tp +mm(me,0)TE Fr djrnp C p , f L E P AxEP + rdT\ Ay EP De-(TE-TP) *f Ae D„ = -°P,f L E P (4.52) (4.53) (4.54) (4.55) Similar work is conducted on all other control surfaces, and the following discrete equation is finally obtained. flpT^X^+S, (4.56) nb where the index "nb" is taken over all neighboring points E, W, N , S. aE = Dc + max(-mc ,0) aw = Dw +max(/ne,0) aN = Dn + max(-m„,0) as = Ds + max(mn,0) nb (4.57a) (4.57b) (4.57c) (4.57d) (4.57e) 62 Chapter 4 Computational Methods The approximation to the reaction heat term using midpoint rule yields, AQ (4.58) VCPJ)P 4.4 Imposing boundary conditions Typical boundary conditions in CFD computations include wall boundary, inflow and outflow boundary conditions, gradient and symmetrical boundary conditions. Although it is usually not hard to implement these boundary conditions, careful consideration is necessary in some cases. As an example, the way to impose symmetrical boundary conditions for velocity equations is given in the following. Suppose the control volume surface OAB in figure 4.3 is the symmetrical plane of the velocity, then the net mass flow through plane OAB should be zero, further, the velocity gradient along the normal direction of plane OAB should be also zero, i.e., (4.59) dV = 0 (4.60) dn , b n C j Vh(uh,vh,wh) A B Figure 4.3 Illustration of symmetrical boundary conditions The norm of plane OAB can be computed by n = OC = OA x OB 63 Chapter 4 Computational Methods The second equation can be satisfied by the following two expressions: V„ • OA = K • OA (4.61) Vh • OB = V, • OB (4.62) Where Vb is the velocity at the boundary surface and Vj is the velocity at the center of the boundary cell. The above equations can be written in matrix form as shown below. "OA "OB y 0 A "OA 'OB \xoc yoc ZOC J f t . \ XOA yOA ZOA XOB yOB ZOB 0 0 0 (4.63) Thus the velocity at the boundary can be calculated. 4.5 The code and its structure The code was written in Microsoft C, with the capability of solving either single- or two-phase flow, two or three-dimensional problems. It uses collocated mesh with curvilinear grids and has a second order accuracy. As shown in the flow chart on Figure 4.4, the code consists of three parts, namely "Grid", "Compute" and "Plot". The function of "Grid" is to generate user-defined mesh; as the main part of the code, "Compute" reads the pre-generated mesh, satisfies the user-required boundary conditions and delivers the converged solution; "Plot" functions as a post-processor, which reads the mesh and solution to produce data files recognizable by TecPlot. 64 Chapter 4 Computational Methods ( Grid ) Geometry definition Generate mesh Save mesh file ( End ) Compute^) I Physics definition Single/Two phase flow, Newtonian/non-Newtonian fluid >H Read mesh file T Set boundary conditions Set iteration parameters and convergence rule Estimate initial solution Parameters check Solve flow field Liquid phase N-S equation Solid phase N-S equation Solve energy & species conservation equations Solve chip compressibility equation ( End ) ( Plot ) Plot definition H Read mesh file H Read solution Generate data files in specific formats ( End ) Figure 4.4 Flow chart of the developed CFD code 65 Chapter 4 Computational Methods 4.6 Code validation Validation is a very important stage in the process of code development. The main purpose of validation is to ensure that the implementations in the code are correct and free of errors. The following section describes some of the validation work that has been carried out. 4.6.1 Flow in a cavity As classical study cases, flows in a cavity, either lid-driven or buoyancy driven, two or three-dimensional, are well documented [59,61]. The geometry and boundary conditions in a two-dimensional situation are shown schematically in Figure 4.5 (a) and (b). / Adiabatic 'A As//////////////////, 1 Cold Adiabatic /777777777777777777777 Moving lid Hot (a) Lid-driven (b) Buoyancy-driven Figure 4.5 Geometry and boundary conditions for 2D cavity flows In Figure 4.5(a), the square cavity has unit height and width, the bottom wall moves with unit speed in the positive x direction, the other walls are all at rest and the flow inside the cavity is assumed to be two-dimensional, steady, incompressible, and laminar with constant transport properties. Figure 4.6 shows the streamline and U-velocity along the vertical centerline when Re number equals 1000, the mesh size used in the 66 Chapter 4 Computational Methods computation is 80x80. In Figure 4.6b, data reproduced from Ferziger and Peric's work [55] using center difference scheme are also plotted, as can be seen, the difference between the two velocity profiles is very small. In Figure 4.6c, the computed U-velocity at the vertical centerline using different mesh sizes are plotted, which indicates the solution is grid-independent. The convergence history of a typical computation given in Figure 4.6d shows that the maximal residual of P, V and V-equations drops below 10"6 after 700 iterations and then remains almost constant. This is because single precision data are used exclusively in the code. (a) Streamlines -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 U (b) U-velocity along vertical centerline (c) Computed U-velocity at vertical centerline using different mesh size 600 900 1200 1500 iteration times (d) Convergence history with mesh size of 40x40 Figure 4.6 Flow field for 2D lid driven problem (Re=1000) 67 Chapter 4 Computational Methods In Figure 4.7 and 4.8, velocity profiles of lid-driven flow in a three-dimensional cavity are plotted, where Re numbers are 100 and 400 respectively. The results are very close to those given by Ching et al. [59]. u x (a) U-velocity at plane: x=y=0.5 (b) V-velocity at plane: y=z=0.5 Figure 4.7 Velocity profiles in 3D lid driven problem (Re=100, mesh size: 25x25) u x (a) U-velocity at plane: x=y=0.5 (b) V-velocity at plane: y=z=0.5 Figure 4.8 Velocity profiles in 3D lid driven problem (Re=400, mesh size: 25x25) In buoyancy-driven cases, as shown in Figure 4.5b, the bottom wall, at a temperature T=T.O, moves with unit speed in the positive x direction. The other walls are at rest with 68 Chapter 4 Computational Methods the temperature T=0.0. Chen et al. [60] has simulated this flow using the finite element method, with contour plots of stream function and temperature solutions for Re=100 shown in Figure 4.9. For comparison, the same problem is solved by the present code and corresponding results are provided in Figure 4.10. (c) Temperature for Re=100, Pr=1.0 (d) Temperature for Re=100, Pr=10 Figure 4.9 Stream function and temperature contours from Chen et al. [60] (Re=100) 69 Chapter 4 Computational Methods x x (c) Temperature for Re=100, Pr=l .0 (d) Temperature for Re=100, Pr=10 Figure 4.10 Stream function and temperature contours by this work (Re=100) To make sure the code also works correctly for curvilinear grids, the two-dimensional flow in an inclined cavity is also simulated. This time, the upper wall moves to the right with a unit speed while all the other walls are at rest. The streamline when the cavity is inclined to 45° and the Re=1000 is given in Figure 4.11; the computational results of the same problem given by Ferziger and Peric's [55] are shown in Figure 4.12. Again, the differences between these two simulations are very small. 70 Chapter 4 Computational Methods x Figure 4.11 Streamlines in 2D lid-driven cavity with both side walls inclined at 45° (Re=1000) Figure 4.12 Streamlines in 2D lid-driven cavity by Ferziger and Peric [55] (with both side walls inclined at 45°, Re=1000) 4.6.2 Liquid flow through porous media along a flat plate We have presented the validation of a number of single-phase flow problems with or without heat transfer, occurring in cavities. However, the cooking process is essentially a two-phase flow problem. Although both phases are described by Navier-Stokes equations, there are inter-phase interaction terms contained in the governing equations for both the chip and free liquor phase. It would be very helpful to further validate the code if these terms can be tested. Lauriat and Vafai [63] have obtained the analytical solution to the 71 Chapter 4 Computational Methods fully developed boundary layer flow when liquid passes through a porous media along a plate. The velocity profile is found to be u = l--^-sec/zV / 2[(v* +c,)/2] c, is the integration constant given by c, = 2cT 1 / 2 sec/T'^a//?) where u =ulux; y' = yI4KIS a = (l + 2F-Re p ) ; B = 2F-Rep/3 Re, (4.68) (4.69) (4.70) (4.71) (4.72) As an example, when p=1.0, p=2xl0"3, «M =1.0, the porosity (or void fraction) 8=0.53, the permeability K=4xl0"4, and the inertial coefficient F=0.55, from the above expressions, Rep=l,a=2.1, (3=0.3667, c, =7.617, thus it is easy to compute the dimensionless velocity u versus the dimensionless height y*. The same problem with identical parameters is also computed with our code, with velocity profiles of a typical case shown in Figure 4.13. As can be seen, the fully developed velocity profiles are obtained downward at a certain distance, where the entrance effect is negligible. N 0.03 0.05 0.07 0.09 0.13 Figure 4.13 Boundary layer flow passing through porous medium along a flat plate (p=1.0, p=2xl0"2, s=0.53, k=4xl0 , Rep=l) 72 Chapter 4 Computational Methods The code provides very good predictions, as the results compare well with the analytical solutions in a wide range of Re numbers (based on pore size), as exhibited in Figure 4.14. As pore size based Re number increases, the flow becomes closer to the plug flow, and the code prediction is more accurate. 5- 2 A Rep=1 (Analytical) O Rep=10 (Analytical) • Rep=100 (Analytical) Rep=1 (Computed) Rep=10 (Computed) Rep=100 (Computed) 'fj A/<fr / A/ <i> / A / A /• r* A^.' Oil A ,y o /" A o / II 0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 U* Figure 4.14 Velocity profiles by code prediction and analytical solutions 4.6.3 Power-law non-Newtonian fluid flowing through a tube Skelland [62] studied the velocity profiles for power law fluid flowing through round tubes and found that under laminar flow conditions, the fully developed dimensionless velocity can be written in the form: u V 3» + l n + l 1 (4.73) where u is the velocity, V is the average velocity, n is the power law index, R is the radius of the tube. To ensure the non-Newtonian models are properly implemented in the code, the power law fluid flowing through a round tube is simulated and compared with equation 73 Chapter 4 Computational Methods (4.73). Figure 4.15 shows the comparison of the velocity profiles predicted with the code and the analytical solution (4.73). 1.0 0.8 0.6 0.4 0.2 0.0, Analytical solution — Numerical solution 0.0 0.5 1.0 1.5 2.0 2.5 3.0 u/V (a) Power law index n=l/9 Analytical solution Numerical solution 0.0 0.5 1.0 1.5 2.0 2.5 3.0 u/V (c) Power law index n=1.0 0.2 0.0 Analytical solution Numerical solution I • • - • l 0.0 0.5 1.0 1.5 2.0 2.5 3.0 u/V (b) Power law index n=l/3 1.0 Analytical solution Numerical solution 0.0 0.5 1.0 1.5 2.0 2.5 3.0 u/V (d) Power law index n=3.0 Figure 4.15 Velocity distributions for power law fluids in laminar flow through round tubes (zero-slip between the fluid and the wall, power law index n varies from 1/9 to 3.0) As can be seen in Figure 4.15, over a large range of power law index n (from 1/9 to 3.0), the code can predict the velocity profile reasonably well; the largest error occurs at the tube center when n equals 3.0; also, fluids with low index n approach plug flow over a large portion of the tube diameter. 74 Chapter 4 Computational Methods For comparison purpose, the code is also applied to an industrial digester (in Kamloops, BC) which was simulated and well documented by Pougatch et al [9]; the simulation results given by the current code are then compared with those in the literature. The comparison shows that there are no essential differences in the velocity distributions of both chip and free liquor phases, and the differences in the chip phase pressure and chip compaction are minimal; however, there are relatively large differences in distributions of temperature and alkali concentration. Discrepancies also exist in the distributions of residual carbohydrates and lignin contents. That is expected, as in the two computational codes, different chemistry models are used, and some of the empirical correlations used are also significantly different. Although extensive validation work has been conducted, on-site experimental data will also be used to assist the validation for the overall cooking process in later chapters. 75 Chapter 5 Modelling of Liquid Flow in a Laboratory Digester In this chapter, the liquid flow in a laboratory digester filled with compacted wood chips is computationally modelled under various conditions. The experimental observations through electrical resistance tomography (ERT) images generally agree with the calculated results, but detailed comparison cannot be performed due to difficulties both in CFD modelling and experimental work. 5.1 Introduction Since the early 1980s continuous cooking has been growing rapidly. Many new cooking techniques have also been emerging. Some of them, including Extended Modified Continuous Cooking (EMCC) and Iso-Thermal Cooking (ITC) [21], have quickly found widespread use in the industry. The cooking performance is essentially a combined result of many variables. Although industries have made great progress in continuous cooking, many underlying physical and chemical mechanisms are still not well understood; continued efforts are still required to further improve the cooking performance, optimize the equipment design and get better control of the operation. The liquid flow in a digester is of significant interest as it affects the motion of chip column, and helps to distribute the chemicals and maintain a high enough cooking temperature. It is expected that a better understanding of the liquor flow will lead to an improved digester process. 76 Chapter 5 Modelling of Liquid Flow in a Laboratory Digester In the past several years, a model digester was built by Prof. Bennington's group (Department of Chemical and Biological Engineering, UBC) and the Electrical Resistance Tomography (ERT) technology was implemented and used to visualize the liquid flow in this digester both with and without wood chips [64,65]. However, the ERT technology can only provide flow information in a limited number of digester cross-sections and therefore, it is sometimes very difficult to get a good understanding of the flow structure. Significant progress has been made in the past decade in Computational Fluid Dynamics that allowed it to become a useful tool in achieving a better understanding of industrial processes. In this chapter, the liquid flow in the model digester will be examined numerically with chips in place. 5.2 The model digester and electrical resistance tomography 5.2.1 The laboratory digester Figure 5.1 schematically shows the basic configuration of the laboratory digester [64]. The cylindrical digester has a volume of 40x10"3 m 3 with an inner diameter of 0.28 m and a height of 0.75 m. Two water circulation loops can be established in the digester. In the first loop (shown in blue arrows), water pumped from a reservoir moves up through the digester bottom, and flows out via three nozzles located on the top of the digester periphery. In the second loop (shown in red arrows), water is extracted through the circumferential screen and sent back into the digester through the downcomer. The screen is partitioned into sections by eight equally located tubes through which the extracted water flows. 77 Radial Flow Downcomer Screen Section Distribution Plate Axial Flow 0.75 m Figure 5.1 Configuration of the model digester To model the flow patterns in real digesters, two fluid circulation loops were created. The first one provides axial upward flow Q A , which is injected from the bottom; while the second one gives radial flow Q R which is injected through the central downcomer. The liquid flows out from the extraction screen at the top of the digester. The downcomer can move up and down so as to model different flow situations. Cooked or uncooked wood chips (usually collected from mills) can be added and compacted against the perforated distribution plate near the bottom of the digester. The flow rates Q A and QR can be set at different values. An industrial Electrical Resistance Tomography (ERT) system was installed to generate an image of the flow field. A number of different experiments have already been conducted with this device [64, 65]. 5.2.2 Electrical Resistance Tomography There are numerous methods/ways to make measurements in fluid flows. The usual experimental probes, such as pitot tube, the hot wire and the more advanced laser-doppler 78 Chapter 5 Modelling of Liquid Flow in a Laboratory Digester equipment, can be characterized as point measurements, i.e., they only give information in a single spatial point in the flow field. Although useful, the information is quite limited, especially when flow phenomena are dominated by spatial structures. Newly developed techniques such as particle tracking velocimetry (PTV), particle image velocimetry (PIV) and laser induced fluorescence (LIF) are capable of obtaining spatial information, but they are quite expensive, and usually are used in gaseous and liquid flows. Electrical Resistance Tomography (ERT) belongs to a group of imaging methods collectively termed electrical tomography techniques, which have been extensively reviewed [66,67,68,69,70]. ERT visualizes the flow by imaging the spatial difference in conductivity through the sample volume. In other words, ERT images are generated from measurements of electrical resistance. The ERT technique has greatly progressed since it was developed in the 1980s and offers many benefits, including relative low-cost, easy data collection, high sensitivity, and fast image rates; it can even be used online to monitor and control the operation of existing process equipment with careful design. ERT has been used extensively in various aqueous flows, and recently, it has been employed in studying gas/liquid two-phase flow [70]. A typical ERT system is composed of three main parts: sensors, data acquisition system and image reconstruction system/host computer. As shown in Figure 5.2, in the model digester, there are eight sensor planes equally spaced along the digester height, each sensor plane consists of 16 electrodes evenly configured in the vessel periphery. 79 Chapter 5 Modelling of Liquid Flow in a Laboratory Digester Figure 5.2 The laboratory digester and the ERT sensor planes The ERT system uses a P2000 data acquisition system to inject current between any pair of electrodes and then measures the potential difference. For the configuration of 16 electrodes, 104 independent measurements are made. Afterwards, the image reconstruction by linear backprojection algorithm is conducted on that plane. The ERT generates an image of conductivity distribution in the cross section by measuring the electrical resistance. In the experiments conducted in the model digester, two techniques were used to image fluid flow. One technique uses controlled pulse of the tracer, the other one uses a step change in water conductivity. The higher conductivity solution can be stored in one of the two storage tanks and introduced into the system as desired. 5.3 Computational method In the laboratory digester, wood chips are stationary; cooking reactions do not occur as cold water instead of hot cooking liquor is used. Therefore, inter-phase mass transfer is not involved, as this is basically a liquid flow through a compacted porous media - chip phase. 80 Chapter 5 Modelling of Liquid Flow in a Laboratory Digester The porosity S f is measured when wood chips are loaded. Then the water flow is introduced. As an approximation, the porosity is considered to be a constant, although the water flow does have some impact on the chip column compaction. To model the liquid flow with chips in place, equation 2.7 and 2.11 are used, with the term accounting for mass transfer eliminated. Because £ f is treated as a given constant, the compressibility equation is not solved. The inter-phase friction force determined by Wang [29] (for reference pine chips) is used in the modelling. Figure 5.3 gives the flow resistance over pine chips. Q. 5 Q. < a e £ 3 (A in £> CL 2 4 6 8 10 12 14 16 18 20 Velocity difference AV, mm/s Figure 5.3 Flow resistance over chip column, determined by Wang [29], for reference pine chips In this study, purely liquid flow, i.e., without the presence of wood chips is not simulated, as in this case, the flow is in the laminar/turbulent transition regime. 81 Chapter 5 Modelling of Liquid Flow in a Laboratory Digester 5.4 Computational and experimental results 5.4.1 Approximation of boundary conditions As already described, there are two liquid loops in the model digester. One is injected through the perforated holes located in the downcomer periphery; the other one is introduced from the bottom of the digester. The injection port on the downcomer consists of many densely distributed small open holes (3mm in diameter); for simplicity, the entire perforated region is considered as a single inlet. To check this assumption, two different cases were computed. In the first case the plane water jets exiting from three slots were modelled with the dimensions and flow rate being typical of those occurring in the downcomer. In the second case an equivalent rectangular geometry was used. The port geometry and velocity distributions are shown in Figure 5.4 and 5.5. 3 mm 1mm l 0 0 l s i I i I i j i i i j i i i j i i i I i i i i i I i I i i i i i 0.4 0.8 1.2 1.6 X, cm (a) Port geometry (b) Calculated velocity distribution Figure 5.4 Plane jet: water injected from triple holes, m0=0.33kg/s 82 Chapter 5 Modelling of Liquid Flow in a Laboratory Digester 11m 0.011 1 1 1 1 1 1 1 ' 1 ' 1 s 1 ' 1 5 1 ' 1 ' 1 ' 1 5 1 5 1 = 1 5 1 0.4 0.8 1.2 1.6 X, cm (a) Port geometry (b) Calculated velocity distribution Figure 5.5 Plane jet: water injected from a single port, m0=0.33kg/s As can be seen in Figure 5.4, the velocity between the slots can quickly increase to the value in the main stream; when the distance is larger than twice the slot height, i.e., x>6mm, the velocity profiles in Figure 5.4 and 5.5 are almost identical. So treating the whole perforated area as a single open port will not cause a large error, especially further downstream from the injection port. Because the downcomer has a dead end at the bottom, the inside water pressure will attain its highest value there. Thus the injection pressure is high in the lower portion and low in the upper portion of the injection port. This leads to higher injection velocities in the lower part of the perforated area. This effect has been observed in the experiments. For the bottom inlet, as a distribution plate is installed, uniform inlet velocity can be assumed. Further, we assume that the outflow through the screen section is uniform. The exit flow from the top of the digester is imposed as an outflow boundary condition (zero gradient condition). 83 Chapter 5 Modelling of Liquid Flow in a Laboratory Digester As an approximation, the flow is assumed to be symmetrical, thus only two-dimensional computation is needed to model the liquid flow; the corresponding geometry and mesh are given in Figure 5.6. Downcomer Figure 5.6 Computational domain, mesh, inlet and outlets 5.4.2 Liquid flow in the digester Vlaev and Bennington [65] investigated experimentally the liquid flow in the model digester in the presence of wood chips. In these experiments, chips can be cooked or uncooked, compressed or uncompressed, and the flow can be symmetrical or asymmetrical. In the following, both the experimental and computational results performed on the symmetrical cases with uncooked chips are presented. Figure 5.7, 5.8 and 5.9 show the ERT images and modelled velocity fields for three different flow conditions with the ratios of axial to radial flows of 0.5, 1.0 and 2.0 respectively. In all the cases, the volume fraction of chip phase is 0.53; tracer is injected from the bottom of the digester. In all the ERT images, the red color represents the highest concentration of the tracer, while the blue indicates the lowest. From the color images, one can 84 Chapter 5 Modelling of Liquid Flow in a Laboratory Digester qualitatively infer the flow path and the velocity magnitude. However, it is generally hard to obtain the exact velocity field because of the low resolution and two-dimensional nature of the ERT images and therefore detailed comparison with the computations is difficult; moreover, the ERT images may not represent accurately the real flow if the images are obtained too long after the injection, especially when the tracer diffuses quickly. r, cm r, cm (a) ERT image, after 95 sec (b) Velocity (c) Streamline Figure 5.7 Flow field given by ERT image and CFD modelling QA=125xlO"6 m3/s, QR=250 xlO"6 m3/s From Figure 5.7 (c), all the liquid injected at the digester bottom exits from the screen section. As part of the bottom inflow, the majority of tracer should flow out through the screen section as well, but there might be a small amount of tracer diffused. However, this upflowing tracer will be likely restrained below the liquid injection level from the downcomer, as part of this injection flows down- and outward and exits from the screen section. This is reflected in the ERT image given in Figure 5.7 (a), where the tracer concentration remains high in the lower portion of the digester, while very low in the upper portion above sensor plane 4. 85 Chapter 5 Modelling of Liquid Flow in a Laboratory Digester (a) ERT image, after 229sec 2 5 8 11 14 r, cm (c) Streamline (b) Velocity Figure 5.8 Flow field given by ERT image and CFD modelling QA=125xlfJ6 m3/s, QR=125 xlO"6 m3/s Figure 5.8 (b) and (c) show the velocity and streamline when a flow rate of 125xl0"6 m /s is used in both the axial and radial loops. The flow field is similar with the previous case shown in Figure 5.7, but most of the liquid injected from the downcomer flows upward and leaves the digester from the top exit. In Figure 5.8 (a), relatively high tracer concentrations are observed up to the sensor plane 4. Again, the tracer concentrations in the upper sensor planes are low, indicating most liquid injected from the bottom exits from the extraction screen. 5 0 P " ' ' 40 | s imm I £ ; '3 1 0 k •• 2 5 8 11 14 r, cm (b) Velocity 2 5 8 11 14 r, cm (c) Streamline (a) ERT image, after 95 sec Figure 5.9 Flow field given by ERT image and CFD modeling, with flow rates QA=250xl0"6 m3/s, QR=125 xlO"6 m3/s 86 Chapter 5 Modelling of Liquid Flow in a Laboratory Digester Figure 5.9 shows the flow field for a higher axial flow rate. From the streamline in Figure 5.9 (c), part of the tracer flows up above the extraction screen, moves up continuously near the digester wall, and leaves the digester from the top exit. After a short time (95 seconds), the tracer is detected as high as the top sensor plane as shown in Figure 5.9 (a); in addition, higher tracer concentrations are found near the digester wall. In all the three cases, both experimental ERT images and computational flow fields show that beyond the circulation region, the upward liquid velocity near the digester periphery is higher than that in the center area, indicating that the flow is not uniform, but after some distance away from the screen section, the liquid flow gradually become uniform. Also, as the flow ratio of Q A : Q R increases, the flow pattern in the digester changes as well. In Figure 5.7 and 5.8, all the fluid injected from the bottom exits from the screen section. In Figure 5.9, roughly half of the bottom liquid exits from the extraction screen; the rest of the liquid flows upward above the screen zone near the digester wall and finally exits from the top of the digester. The upward axial velocity in the lower part of the digester increases from case (1) to (3). In case (1), the upward axial velocity is so small that the tracer is experimentally detected mainly in the bottom four plates, as shown in Figure 5.7. In case (2), the tracer can reach the top plate 7, but takes a very long time; the reason is that the bottom flow is constrained to the lower part of the vessel. When compared with (1), the magnitude of the axial velocity in the lower part of the digester increases, and the tracer can be easily transported to plate 3 and 4, as afterwards, it spreads further by diffusion and convection; therefore a moderate tracer concentration is observed, as indicated by the green color shown in Figure 5.8 (a). In case (3), computational results shows that a portion of the bottom flow is pushed to the digester wall and continues to the digester top when it spreads beyond the screen section. On the other hand, the tomographic images in Figure 5.9 (a) indicate that the tracer can rapidly reach the top plate 7. Images in plate 5, 6, and 7 in Figure 5.9 (a) exhibit higher tracer concentrations in the outer layer and this agrees with the computational results in Figure 5.9. 87 Chapter 5 Modelling of Liquid Flow in a Laboratory Digester Similar circulation loops exist in a number of different places in industrial digesters. Knowing the above flow patterns may bring a better understanding that can lead to a better design or operation. Taking wash circulation in an EMCC digester as an example, the flow pattern shown in Figure 5.8 is probably the best in terms of cooking and heating. In the case of Figure 5.7, part of the heated wash liquor and added white liquor are extracted and sent back to the wash liquor circulation loop again; in the flow conditions of Figure 5.9, however, only a part of the wash liquor is extracted and heated; the temperature and white liquor distribution at the bottom of the digester are not uniform, which may undermine the cooking uniformity. The distributions of the liquid dynamic pressure for the three cases are given in Figure 5.10. r, cm r, cm r, cm (a) (b) (c) Figure 5.10 Calculated liquid dynamic pressure (with chips in place) (a) QA=125xlO"6 m3/s, QR=250xlO"6 m3/s; (b) QA= QR=125xlO"6 m3/s; (c) QA=250xlO"6 m3/s, QR=125xlO"6 m3/s; A reference pressure of 22 Pa is set to the digester top. From Figure 5.10, the liquid pressure is high at the injection places and low at the exits. In addition, the circulation ratios have a significant effect on the pressure distributions: in Figure 5.7, the radial flow is large and thus the pressure at the injection port of the downcomer is high; in Figure 88 Chapter 5 Modelling of Liquid Flow in a Laboratory Digester 5.9, the axial flow is twice as large as the radial one, so the upward velocity is large and therefore the liquid pressure at the digester bottom is high. 5.5 Conclusions The liquid flow in a laboratory digester where chips are filled in is modelled computationally. Calculated results agree reasonably well with the experimental observations obtained by ERT images. The computations have led to the following conclusions: In the presence of compacted wood chips, the momentum of the liquid diffuses very efficiently and the liquid flow tends to be uniform away from the injection and extraction places. Depending on the liquor circulation ratios, different flow patterns can occur in a digester. Appropriate flow rates will be beneficial to uniform cooking. As expected, a higher circulation rate requires a higher injection pressure. Because of the low resolution of the ERT images, detailed comparison between computation and experiments is difficult. 89 Chapter 6 Simulation of an Industrial Continuous Digester As pointed out in previous chapters, continuous cooking is a very complex process. It is also a process hard to control mainly for the following reasons [71]. - In practice, there are usually big variations in chip quality, such as moisture content, size distribution etc; sometimes, the wood species also changes; - It is almost impossible to make online measurements inside a digester because of the presence of moving solid chips and the hostile environment; - There is a long process delay; it can take even 5 hours for an individual chip to move from the chip bin to the digester blow line. To maximize the economical benefits, better knowledge and control of the pulping process are required. CFD technology can help meet this demand in two ways. CFD technology can help to identify some key parameters and assess their influence on the pulping process. This is very useful, as for instance, predicting with reasonable accuracy the impact of upstream variations will allow operators to make adequate adjustments in a timely manner. CFD can also help to achieve a better pulping process by predicting the pulping performance before any desired changes are actually made. This helps to give a good retrofit or design that could significantly decrease the capital cost, and improves the pulping quality. This chapter introduces the computational work performed on an EMCC continuous digester operated by Howe Sound Pulp and Paper Ltd. (HSPP). The cold blow is assumed to be a uniform up-flow through the digester bottom; all boundary conditions are derived from the mill operating data. A typical simulation produces the velocity profiles, pressure distributions and volume fractions for both phases; the temperature, effective alkali concentration and kappa number are also obtained. The impacts of variations in production rate, chip type and cold blow rate on cooking performance are also 90 Chapter 6 Simulation of an Industrial Continuous Digester computationally evaluated. In this Chapter, a three-dimensional simulation for the continuous digester is also presented to evaluate the non-uniformity of cooking when the liquor extractions are made in a special manner. 6.1 Introduction 6.1.1 Two-vessel hydraulic system in HSPP Howe Sound Pulp and Paper Ltd. (HSPP) produces both mechanical and kraft pulp. The kraft pulp in HSPP is sold directly in the world market, instead of being transformed into different paper products. Currently, HSPP produces 1100 tonnes of kraft pulp each day. The cooking equipment used in HSPP is a two-vessel hydraulic continuous type, which has a separate impregnation vessel. A technique called Extended Modified Continuous Cooking (EMCC) is being used. Chips, Liquor Chips 1. Steam Chip Bin Steam J High Pressure Feeder Cold Blow Brown Stock Figure 6.1 Two-vessel hydraulic digester system in HSPP 91 Chapter 6 Simulation of an Industrial Continuous Digester Figure 6.1 is a schematic diagram illustrating this cooking system. In this system, the white liquor charge is divided into four parts. The first charge is delivered to the impregnation vessel. The second one is delivered to the top of the digester where the chips are transferred from the impregnation vessel into the digester. A third charge is introduced to the loop of modified cooking liquor circulation, while the last one is introduced to a place above the washing screen through the central downcomer. In the digester, chips move vertically downward continuously, the spent cooking liquor and washing liquor are extracted through a circumferential screen located in the middle of the digester, thus a co-current flow zone is formed in the upper portion of the digester while a counter-current flow occurs in the lower part of the digester. There are also several heating circulations to maintain the required reaction temperatures. The digester can be divided into four zones, namely cooking zone, MCC (Modified Continuous Cooking) zone, EMCC (Extended Modified Continuous Cooking) zone and washing zone. As compared with conventional cooking, EMCC has several advantages: (1) Because part of the cooking is conducted in the countercurrent zone where the liquor has a much lower solid content, the lignin can be further removed without much cellulose degradation. Therefore EMCC has the capacity to cook chips to very low kappa numbers while retaining pulp strength and yield. It was reported that with EMCC, the kappa number could be lowered to 22-24 for softwood and to 14-16 for birch (hard wood). (2) The alkali is charged into several different places to make the alkali content more uniform during the cook. Consequently, the time for delignification becomes longer, but the non-uniformity of cooking decreases. (3) The digester volume is more effectively utilized as cooking reactions take place in the lower part of the digester as well. (4) This process reduces the need for bleaching chemicals and consequently can make a major contribution to the reduction of emissions of chlorinated organic compounds and some other pollutants in the pulp and paper industry. 92 Chapter 6 Simulation of an Industrial Continuous Digester 6.1.2 Mill records of digester operations To facilitate our computational study, HSPP has generously provided us with comprehensive data sets on a number of operations of the EMCC digester. It was a serious effort to keep all operations stable for a minimum of 6 hours; the time-averaged values are obtained from the hourly measured data. Appendix A shows the mill records for two typical cases; Appendix B lists the averaged data for a total of 15 cases. In all computations, only the averaged data are used to determine the boundary conditions. 6.2 Blended chip species Wood pulp consisting of longer fibers is stronger than that of short fibers; so many pulp mills blend different wood species, usually two at a time, to ensure that the produced pulp meets the target strength required by customers. Table 6.1 gives some laboratory measurements of the average fiber length of different wood species in British Columbia [72]. Table 6.1 Average fiber length of different wood species [72] Species Average length, mm BC interior varieties Lodgepole pine 2.3 Spruce 1.9 Douglas fir 2.0 BC coastal varieties Hemlock 3.0 True fir 3.3 Douglas fir 3.4 Western red cedar 2.4 Yellow cedar 2.5 For the simulations, the oven-dry wood density is required. For blended wood chips, the oven-dry wood density can be evaluated by the weighted sum of all the species involved, or mathematically, 93 Chapter 6 Simulation o f an Industrial Continuous Digester PoDw=YXarPoDWj) C6-1) Here, a\ is the mass fraction of the i-th species in the mixture, pODW, is the oven dry wood density of the i-th species. In Howe Sound Pulp and Paper Ltd., four main species are used, namely hemlock, cedar, spruce and Douglas fir. The average specific gravities of these four species are given in Table 6.2. Table 6.2 Specific gravity of some common woods [72] Wood species Specific Gravity Western hemlock 0.44 Western red cedar [74] 0.31 Engelmann spruce 0.36 Douglas-fir, coast-type 0.51 6.3 Computational domain The schematic of the digester used in HSPP is given in Figure 6.2. The digester has a total height of 51.2 meters and a radius of 4.228 meters of the bottom shell. The digester produces around 1300 tons of air-dry pulp per day, with a target kappa number of 28. The chips form a free surface (which is called digester chip level) at a short distance from the digester top. In the digester operation, a fixed chip level is desired so as to maintain a stable operation. In practice, the chip level changes in a certain range. Above this level, wood chips and cooking liquor fall down by gravity. To avoid the difficulty of modelling the complex flow in the top part of the digester, we start the computation below the chip level and assume that chips and cooking liquor enter the digester with uniform velocities. 94 Chapter 6 Simulation of an Industrial Continuous Digester While developing the energy equation in Chapter 2, we assumed that the chips and liquor are in thermal equilibrium, namely the temperatures are identical for both phases. This is a commonly used assumption for liquid flowing through porous media. Then we solve only one energy equation. However, for the digester bottom portion where blowing occurs, the temperatures of the cold blow liquor and chips are significantly different. In addition, the cold blowing liquor is usually injected upward through a number of nozzles and that is difficult to simulate in detail. A s an approximation, the C F D computational domain ends 1 meter below the wash liquor screen, and the cold blow is assumed to be a uniform upflow at the bottom plane; this dramatically simplifies the blowing process. In Figure 6.2, the computational domain is shown in the enclosed box with dotted line. Chips & Liquor •7T>— B C L Cold Blow Trim liquor (not used) M C L W L Brown Stock Figure 6.2 Illustration of the computational domain of the digester (BCL: bulk cooking liquor; B L : black liquor; M C L : modified cooking liquor; WL: wash liquor; Trim liquor: liquor is drawn from chips and circulated to bottom circuit.) 95 Chapter 6 Simulation of an Industrial Continuous Digester 6.4 Boundary conditions 6.4.1 Boundary conditions for velocity In Chapter 2, the chip phase is defined as wood chips and the entrapped liquor, and liquor phase is defined as free liquor. As the Navier-Stokes equations are solved for both the liquor and chip phases, boundary conditions for the velocity should be provided separately for each phase. For the liquor phase, the bottom plane is considered as an outflow boundary, while all the rest of the ports are treated as inflow boundaries. For the chip phase, the bottom plane is also considered as an outflow boundary and the top plane is considered as an inflow boundary; in addition, the slip condition is used for both the downcomer and the digester walls. For all inflow boundaries, the flow rate of the corresponding phase is needed; however, this is not so obvious, for instance, at the digester top, both solid and liquid flow through the same port, but neither of the flow rates is measured there. To compute the velocity at the boundaries for each phase, the flow rate of each phase and the flow area occupied by each phase have to be determined. In this section, the way of estimating the velocities at the digester top and bottom is presented. 6.4.1.1 The top inlet (1) Compute flow areas Both solid and liquid are introduced at the digester top inlet. Based on the assumptions made earlier, the chip phase is a porous matrix; and the free liquor flows through the interstitial space of the porous matrix. When the porosity s is determined, the flow area of each phase can be calculated as, In Equation (6.1) and (6.2), A is the total area through which the chip phase and free liquor phase flow. The chip phase pressure at the chip level is assumed to be zero, i.e., Pc=0. Then from the chip column compaction equation given by Harkonen, (6.2) Ac=ecA = (\-e,)A (6.3) 96 Chapter 6 Simulation of an Industrial Continuous Digester ef = 0.644 + V10 4y (-0.831+ 0.139 In*-) (2.27a) the volume fraction of liquid is st=0.644, thus the volume fraction of solid is sc=T-Sf=0.356. Further, we can calculate the volume fraction of wood material. Suppose the wood material density is /^ lSOOkg/m 3 If the blended chips consist of 43.4% cedar and 56.6% spruce, according to Equation (6.1), the oven dry wood density can be calculated as, p0DW =0.434x310+0.566x360= 338.3 kg/m3 (In case of blended chips, the oven-dry wood density can be calculated by equation 6.1). We have, = Poo, (6-4) £c = P0DW£c = Q 0 8 0 3 ( 6 5 ) Pp. =0.276 (6.6) The volume fractions of the wood material ep w and the entrapped liquor sen will be used to determine the chip phase density in the subsequent section. (2) Calculate flow rate for each phase The digester has a complex inflow/outflow configuration and requires careful consideration when determining the velocity at its boundaries. In the two-vessel system (Refer to Figure 6.1), the total inflow via the digester top includes chips (with moisture), white liquor, consumed steam in the chip bin and steaming vessel, cold blow to the 97 Chapter 6 Simulation of an Industrial Continuous Digester impregnation vessel, and the bulk cooking circulation. (Some liquor is extracted near the digester top, its temperature is increased and then it is directed to the bottom of the impregnation vessel so as to carry wood chips to the main digester vessel; this is called the bulk cooking circulation. This is shown in Figure 6.1). For instance, selected records on 09-MAR-2004 from HSPP digester operations (Refer to Appendix B ) are listed on Table 6.3. Table 6.3 Selection of recorded data from HSPP digester (09-MAR-2004) C o l u m n i n A p p e n d i x B Descr ip t ion Recorded data U n i t A E Weighted average moist 54.26 % A F B l o w Consistency C b | o w 10.00 % A Y Total chips 2428 B D t / d L Digester B low F low 131.9 L/s B E Upper Slice F L 121.6 L/s B D Bulk Cooking Circulation 150.4 L/s 0 H P F W L P U R G E 3.94 L/s P W L Flow to W A S H 4.69 L/s Q W L F low to B C 1.74 L/s R W L Flow to M C 4.68 L/s S W L F low to I V 45.99 L/s T C B F low to IV 0.41 U Cold B l o w to Digester Bottom I 47.78 L/s V Cold B l o w to Digester Bottom 11 10.91 L/s X Cold B l o w to Digester Bottom III 39.54 L/s B F M C Extraction 120.1 L/s W Black Liquor Extraction 139.5 L/s Z B l o w Rate 1267.6 A D t / d 98 Chapter 6 Simulation of an Industrial Continuous Digester The solid/liquor phase flow rate through the top inlet can be obtained by the following procedure: (1) Compute the total inflow M j n (from the digester top and bottom) At steady state flow conditions, the total inflow equals the total outflow, and the total outflow is the sum of the extracted black liquor and the digester blow rate. In Table 6.4, the computation of the total outflow is outlined in rows 1 to 3. (2) Determine the total inflow from the digester top Mi n _to P Wood chips enter the digester from the top; and cooking liquor enters the digester either from the top or the bottom plane. In the bottom, the incoming flow is the cold blow, in the form of three splits. When the cold blow is subtracted from the total inflow, the inflow at the digester top is obtained. This is shown in rows 4 to 5 in Table 6.4. (3) Calculate the flow rate of the solid/liquor phase The inflow of solid chips is known; however, this is not the flow rate of the chip phase. To determine the flow rates of chip and liquor phases, we assume that the chips are fully impregnated. In other words, the interstitial space of wood chips is completely occupied by the cooking liquor. In this situation, the chip phase flow rate will be, m c . p = mpW + men = mpw j 1 Pliq ^ en = 28.10 P £ r pw pw 1 + 1050x0.276 1500x0.0803 =95.71 kg/s And the chip phase density is, r. = P p w S p w + P " " £ c n =1151 kg/m3 Pc Then the flow rate of the incoming liquor phase is, mfp = mm -mcp = 183.35-95.71=87.64 kg/s 99 Chapter 6 Simulation of an Industrial Continuous Digester Table 6.4 Calculation of solid/liquor phase flow rate through the top inlet Row Description Vi Unit Cited data in Table 6.3 1 Black liquor extracted 147.5 kg/s px(W) 2 Digester blow rate 140.5 kg/s px(L) 3 Total mass exiting from the digester M O U T 288.0 kg/s 4 Total cold blow from the bottom M C B 104.6 kg/s px(U+V+X) 5 Total mass inflow via the top M I N TOP 183.4 kg/s V3-V4 6 Total solid inflow (chip flow rate) 28.10 kg/s (AY) 7 Chip phase flow rate via the top M C P 95.7 kg/s 8 Net liquor phase flow rate via the top M F P 87.6 kg/s 9 Blow rate M b i o w 13.2 kg/s (Z) 10 Total flow area A (for both phases) 37.7 m 2 Note: Suppose pi,q=1050, pwater=1000, PcoidBiow=1065 kg/m3 (3) Calculate the velocity at boundary The average velocities of solid and liquor phase are given by, m Wc=—c-*- (6.7) Pc£cA YYl Wf = (6.8) pfefA 6.4.1.2 Bottom exit The liquor flow at the digester bottom is complex. As illustrated in Figure 6.3, the upflowing cold blow is split into three parts, which enter the digester via different channels. A cap and a distributing device are set immediately above the digester exit to 100 Chapter 6 Simulation o f an Industrial Continuous Digester distribute the liquor flow. In practice, the liquor flow in this part could be very complex, for instance, a high injection flow-rate from the side tubes could result in channeling near the digester wall . Consequently, most of the cold blow would move up rapidly through the channeling. In the current stage of our computation, however, this occurrence cannot be captured and assumptions are made to approximate the liquor flow at the digester bottom. A s shown in Figure 6.3, the cold blow is assumed to be a uniform upflow through the bottom plane A i A ] of the computational domain. Brown Stock Cold B l o w Figure 6.3 Illustration of the liquor flow at the digester bottom It is important to note, however, that not all the cold blow moves upward because a portion of the liquor exits along with the chip phase and dilutes the resulting fibers in the subsequent blowing process. To estimate the net upflow of the free liquor phase, we start with the mass balance of the liquor phase, m liq,top 1 7 1 liq,bottom + m\2 m extract (6.9) Equation (6.9) states that the incoming free liquor from the digester top and bottom plane plus the mass transferred from the chip phase equals the extracted black liquor. In 101 Chapter 6 Simulation o f an Industrial Continuous Digester the previous calculation, the free liquor via the digester top and the extracted black liquor are obtained, as given in Table 6.4 row 1 and 8. From equation (6.9), if the mass transfer rate mn is known, the free liquor flow from the digester bottom plane can be calculated. In column Z of Appendix B (also see Table 6.3), the digester blow rate is 1267.55 air-dry tones per day, which is equivalent to 13.2kg/s bone-dry pulp; in Table 6.4, the total solid inflow (chip flow rate) is 28.10 kg/s, so the dissolved solid material is 14.9 kg/s, which corresponds to a dissolved volume AQ. If one assumes that this volume is fully replaced by the surrounding cooking liquor, then the net mass transfer from the chip phase into free liquor phase mn can be approximated by, rhn = (pw - A,, )AQ = (1500-1065) x i ^ = 4.32 kg/s (6.10) Using Equation (6.9), the upflow of the free liquor is then = - mlig,lop ~ ™n = 147.5-87.64-4.32=55.54 kg/s (6.11) 6.4.2 Boundary values for species equation (Lignin, Carbohydrates and EA) In Chapter 3, the proposed pulping reaction model involves four species: lignin, carbohydrates, effective alkali, and sulfidity. However, the sulfidity is usually thought to be a constant through the entire digester and is only measured in the fresh cooking liquor stream. Therefore the measured sulfidity is given as a constant in the digester modelling, instead of solving the conservation equation of sulfidity. Obviously, boundary values are required while solving the rest of the three species conservation equations. 6.4.2.1 Lignin and Carbohydrates The lignin content is directly related to the kappa number, which is an important factor in determining the compressibility and porosity of wood chips. The final contents 102 Chapter 6 Simulation of an Industrial Continuous Digester of lignin and carbohydrates also determine the yield of the pulping process. Therefore, it is important to predict accurately the contents of lignin and carbohydrates. In the species equations given in Chapter 2, the mass fractions, instead of the actual contents of the lignin and carbohydrates, are used. For lignin and carbohydrates, the mass fraction is defined as the residual lignin/carbohydrates content in wood chips over the initial lignin/carbohydrates content. Based on the definition, the mass fractions of lignin and carbohydrates at the digester top inlet should both be 1.0 when a single vessel digester is used. In HSPP, however, a two-vessel hydraulic digester system is used. Thus the initial delignification reactions taking place in the impregnation vessel have to be considered. It was estimated that 5-15% lignin is consumed in the impregnation process. In our computations, it is assumed that 10% of lignin and carbohydrates are dissolved since the exact amounts cannot be determined. Therefore, the mass fractions of lignin and carbohydrates at the digester top inlet are both 0.9. Solving the species equations will give the mass fractions of residual lignin and carbohydrates; in order to obtain the actual contents of residual lignin and carbohydrates, initial contents of lignin and carbohydrates should be provided in a specific pulping process. The contents of lignin and carbohydrates in wood species are affected by many variables, including geographical location, soil and weather conditions, and location of the wood within a given tree. Table 6.5 and Table 6.6 give the cellulose and lignin contents of some soft wood species by different authors. From these two Tables, the cellulose contents among different soft wood species vary only about 3%, while the lignin contents differ significantly; for instance, the lignin content of western hemlock is 11.4%) higher than that of eastern hemlock. Both authors have provided the cellulose and lignin contents of Douglas-fir with some differences, which might be a consequence of different testing methods and experimental errors, or a difference in the chosen samples. 103 Chapter 6 Simulation of an Industrial Continuous Digester Table 6.5 Chemical composition of some soft wood species [75] (% of oven-dry wood) Species Cellulose Lignin Norway spruce 43 29 White spruce 45 27 Black spruce 44 30 Red spruce 44 28 Douglas-fir 44 32 Northern white cedar 43 31 Eastern hemlock 44 33 Table 6.6 Chemical composition of some soft wood in North America [76] (% of oven-dry wood) Species Cellulose Lignin Western hemlock 42.9 21.6 Western red cedar 45.1 20.5 Engelmann spruce 43.6 26.5 Douglas-fir 42.2 29.6 However, the contents of lignin and carbohydrates in wood chips are not measured in HSPP. Because HSPP has a very stable chip source, variations of chemical compositions in the same species are expected to be small. Here, we simply use the data of Table 6.6 to determine lignin contents since these wood species are from North America. The hemi-cellulose content is fixed at 27% of oven-dry wood as an averaged value for softwood [76]. As pointed out earlier, blended wood chips are used in HSPP; the initial contents of lignin and carbohydrates are then computed by, Li8o =TJ(ai-LiSo,,) (6.12) 104 Chapter 6 Simulation of an Industrial Continuous Digester Carb0 = (at • Carb0 t) (6.13) Here, LigQ ,. and Carb0, are the initial contents of lignin and carbohydrates in the i-th species, and the carbohydrates are considered as the sum of cellulose and hemi-cellulose. 6.4.2.2 Effective Alkali concentration In digester operations, effective alkali (EA) concentration is only measured in the white liquor, so EA concentrations in the digester inlet, outlet and extraction ports need to be computed. In the proposed pulping reaction model, mole concentration of effective alkali is used. It is therefore convenient to solve EA concentration equation in its mole concentration form, i.e., Equation (2.18). In the following, the way to calculate EA mole concentrations at the top inlet and at the MC extraction/injection ports is illustrated. 1. The top inlet The white liquor flows to the high-pressure feeder (O), impregnation vessel (S) and bulk cooking circulation and eventually enters the digester via its top plane, and the total amount is, MWL_TOP=3.94+45.99+1.74=51.67 L/s (6.14) The composition of white liquor used in HSPP is given in Table 6.7. Table 6.7 Composition of white liquor in HSPP (g/L) NaOH 69.6 Na 2 S 2 0 3 2.0 Na2S 33.6 EA as NaOH 86.84 N a 2 C 0 3 16.8 A A 103.2 N a 2 S 0 4 4.9 TTA 120 N a 2 S 0 3 1.0 Sulphidity as NaOH 3 2 . 6 0 % As listed in the Table, EA as NaOH is 86.4 g/L, so the mole concentration of EA in the white liquor is, 105 Chapter 6 Simulation o f an Industrial Continuous Digester EA 86.84 2.171 mol/L (6.15) NaOH 40 In Table 6.4 rows 5 and 6, the total mass coming in via the digester top is 183.35 kg/s, among which the solid inflow is 28.10 kg/s, so the liquid inflow via the digester top is 155.24 kg/s. From Figure 6.1, the bulk cooking circulation liquor QBC in the amount of 150.44 L/s (see column B D on Appendix B ) is also introduced at the digester top, with a residual alkali of 12.79 g/L (see column C in Appendix B ) , or equivalently, 0.3198 mol/L. So the mole concentration of EA at the digester top inlet is, 2. M C extraction port/injection ports Figure 6.4 schematically shows the modified cooking circulation loop. A portion of the cooking liquor is extracted from the M C screen. After being mixed with some added white liquor and heated up, it is re-circulated into the digester through the injection port located on the downcomer. [EA] T o p = Total EA (mol) 51.67x2.171 + 150.44x0.3198 = 0.499 mol/L (6.16) Total liquid (liter) 155.24 /1.05 + 150.44 V <$> M C heater M C screen A M C pump „ White liquor - addition Figure 6.4 The extraction/injection ports of the modified cooking circulation ((g) ~ extraction port; (§) - injection port) 106 Chapter 6 Simulation of an Industrial Continuous Digester The selection of the relevant data from Appendix B is given on Table 6.8. Table 6.8 Data relevant to the modified cooking circulation (09-MAR-2004) Column in Appendix B Description Recorded data Unit R White Liquor Flow to MC, Q W L MC 4.68 L/s BF Flow rate of MC Extraction, Q M C 120.1 L/s BP MC Residual Alkali, EA M c 14.2 g/L Note: Suppose pi;q=1050 The mole concentration of EA at the extraction port is then, 14.16 FA [ E A ] = J ^ M L W n ISA 40 = 0.354 mol/L (6.17) The mole concentration of EA in the white liquor [EA]0 is given in Equation (6.15); at the injection port of MC circulation loop, we have, QMC[EA]2+QWL_MC[EA]0 120.08x0.354 + 4.68x2.171 [EA]3 = QMC+Q, WL MC 120.08 + 4.68 = 0.422 mol/L (6.18) 6.4.3 Boundary conditions for the energy equation To compute the energy equation, an outflow boundary condition is set at the bottom exit, while the inflow boundary conditions are imposed at all the remaining. The downcomer and digester walls are assumed to be adiabatic, thus zero heat flux condition is used. 107 Chapter 6 Simulation of an Industrial Continuous Digester 6.5 Simulation results 6.5.1 Boundary values and model constants Because of the symmetrical geometry and operating conditions, the HSPP digester can be well simulated by two-dimensional computations. In Appendix B, 15 sets of data obtained in stable operations of HSPP digester are listed. For each of these cases, the boundary values are determined with the procedures described in Chapter 6, section 4. As an example, Table 6.9 summarizes the computed flow rates of the liquor phase for the digester operation run on March 9, 2004; the temperature and effective alkali concentrations in the corresponding streams are also shown. Table 6.9 Computed boundary values for the case run on March 9, 2004 Liquor St ream Liquor f low (kg/s) Temp (°C) EA concentrat ion (mol/L) Liquor flow via top inlet (in) 245.8 157.3 0.4998 Liquor upflow via bottom port (in) 58.6 151.2 -Wash liquor v ia downcomer (in) 105.0 160.6 0.254 M C circulation via downcomer (in) 131.0 162.6 0.422 Extracted wash liquor (out) 100.1 153.2 0.159 Extracted M C circulation (out) 126.1 159.9 0.354 Black liquor v ia screen (Out) 146.5 159.8 0.260 Bu lk cooking liquor circulation (out) 158.0 157.7 0.320 Other conditions and assumptions Chip phase inflow: 107.43 kg/s; Production rate: 1267.55 ADT/Day [ * ] ; Lignin/carbohydrate contents: 23.93%, 71.07%; Assume 10% of lignin and carbohydrates are dissolved in the impregnation vessel. [*] Note: The term bone-dry (BD) or oven-dry (OD), means pulp from which all water has been removed by heating, but this condition is seldom obtained or sought in the pulp mi l l . Air-dry ( A D ) pulp is understood to contain 10 percent moisture in mass; " A D T / D a y " means tons of air-dry pulp per day. For the chip phase, the computed inflow rate via the top of the digester is 107.43 kg/s. The slip condition is applied on both inner and outer walls with p=0.85. The details of the "slip model" are given in Equation (2.43) and (2.44). To describe the non-Newtonian behavior of the chip phase, the power law model in Equation (2.38) is used with K2 and n equal to 300 Pa • 5 1 2 and 0.8 respectively. 108 Chapter 6 Simulation of an Industrial Continuous Digester 6.5.2 The mesh As discussed before, two-dimensional computations are adequate for the HSPP digester as no circumferential flow is expected under normal operating conditions. To model the axisymmetrical flow, one can either solve the governing equations in cylindrical coordinates with a plane mesh or solve the equations in Cartesian coordinates with a special mesh. Here, the second approach is used. A structured two-dimensional mesh of 40x200 is generated for a small sector of the digester by the developed code, as shown in Figure 6.5. The angle of the sector can be arbitrary, but should be reasonably small so that the curvature of the digester wall can be well fitted by the mesh cells. Using such a mesh is advantageous, as it eliminates the need to solve equations in the cylindrical coordinates. This approach has also been used in FLUENT, a commercial CFD software. In 6.5, the mesh is shown in reduced size for the viewing purpose. | Chips & Liquor Figure 6.5 The structured 2-D mesh shown in reduced mesh size for better viewing (The original mesh size is 40x200) 109 Chapter 6 Simulation of an Industrial Continuous Digester 6.5.3 Computational results With the determined boundary conditions and the generated mesh, the computed results are given in Figure 6.6. uiimiMiuuuu •umii milium Miimiii tuuui 1 5 f c 1 2 3 4 r, m (a) r, m '5 (b) r, m (C) (d) (e) 110 Chapter 6 Simulation of an Industrial Continuous Digester 451 'as 161.9 161.5 161.1 160.7 160.4 160.0 160.0 159.6 159.3 158.9 158.5 158.1 157.8 157.4 157.0 156.6 156.3 155.9 155.5 155.2 154.8 148.0 132.3 122.2 102.4 EA (g/L) 'a 47.3 46.1 44.8 43.5 42.3 41.0 39.8 38.5 37.3 36.0 34.8 33.5 32.2 31.0 29.7 28.5 27.2 26.0 24.7 23.5 22.2 21.0 O) kappa # •— 111.0 — 106.4 — 101.7 97.1 9Z5 87.8 83.2 78.6 73.9 69.3 64.7 60.0 55.4 50.8 46.1 41.5 36.8 32.2 27.6 12 3 4 r, m r, m r, m (f) (g) (h) Figure 6.6 Computational results for the operating case from HSPP on March 9, 2004 (a) Liquor phase velocity; (b) Liquor phase streamline; (c) Liquor phase dynamic pressure; (d) Chip phase pressure; (e) Volume fraction of the chip phase; (f) Temperature; (g) Effective alkali; (h) kappa number distribution. 6.5.3.1 Distributions of pressure and volume fraction Figure 6.6 (a), (b) and (c) show the distributions of the velocity, streamline and pressure of the liquor phase; (d) gives the pressure distribution in the chip phase; and (e) shows the volume fraction of the chip phase. Figure 6.6 (f) shows the temperature distribution for both phases, as thermal equilibrium is assumed between the chip phase and liquor phase. Figure 6.6 (g) shows the mole concentration of effective alkali (over the total liquid in the unit of mol/L), which is high in the central downcomer area and low near the digester periphery. Figure 6.6 (h) shows the distribution of kappa number, thus lignin content, as these two quantities are directly proportional to each other. I l l Chapter 6 Simulation of an Industrial Continuous Digester As can be seen, there is a relatively larger pressure gradient near the liquor injection and extraction ports in both free liquor and chip phases. This indicates a strong interaction between the liquor and chip phase exists. To see this clearly, the liquor velocity and chip phase pressure near the modified cooking circulation level are plotted in Figure 6.7. E a r, m r, m (a) Liquor phase velocity (b) Chip phase pressure Figure 6.7 The liquor phase velocity and the chip phase pressure near the modified cooking circulation level As shown in Figure 6.7 (a), the free liquor has a very high velocity near the injection port due to a small surface area of open nozzles in the downcomer; it spreads quickly in the presence of the compacted wood chips. Figure 6.6 (a) and (b) indicate that a large portion of the injected liquor moves down and outward, being extracted via the modified cooking (MC) screen; the rest of the injected liquor moves upward and is extracted through the black liquor extraction screen. While flowing out through the porous chip phase, the free liquor gradually loses its pressure; the resulting pressure gradient in the liquor phase then raises the chip phase pressure along that direction. Therefore, in the horizontal direction of the injection level, the chip phase pressure is higher on the right end and lower in the injection port, as shown in Figure 6.7 (b). A similar situation happens near the black liquor extraction port, as shown in Figure 6.8. It can also be noted 112 Chapter 6 Simulation of an Industrial Continuous Digester from Figure 6.8 (a) that the compacted wood chips greatly enhances the momentum transfer between phases, resulting in a very flat velocity profile even in the wall vicinities. 36 35 34 33 F DUiimiiiiiiiiuiiiii p i i m u i u i u i u i u i i i in im siiiuuimiiuuiiP niuiiimmiimupi niiimmmumiuw uuuumwuuuup B M i m u m m u u w u 32 F 31b 30 F 29 tumww\ \ \\ • i n / / / ; / / / / / / / ////////ML tmtuin ttttffltmk ?; 'ttlllttttttttlJI m m r, m r, m (a) Liquor phase velocity (b) Chip phase pressure Figure 6.8 The liquor phase velocity and the chip phase pressure near the black liquor extraction port As introduced earlier, wherever a liquid flows across the porous media - the chip phase - there is a pressure drop. It is easy to understand that the pressure of the liquor phase is high in the top and bottom, low in the middle; high near the injection ports, low near the extraction screens. The chip phase moves downward continuously, while the liquor phase flows downward in the upper part, and moves mostly upward in the lower part of the digester, the velocity difference between the liquor and chip phases in the lower part is larger than that in the upper part, and therefore the interphase drag is larger. Also the upflowing liquor flows through a longer path before being extracted. Therefore the liquor phase pressure at the bottom is much higher than at the top of the digester, as shown in Figure 6.6 (c). 113 Chapter 6 Simulation of an Industrial Continuous Digester The chip phase pressure in the lower part of the digester is decreased significantly by the counter-current liquor flow from the bottom. Nevertheless, the chip phase pressure is still high in the bottom while low at the top due to the body force, as shown in Figure 6.6 (d). Consequently, the volume fraction of the chip phase is also high in the bottom, as given in Figure 6.6 (f). 6.5.3.2 Velocity of the chip phase Figure 6.9 presents the chip phase velocity in vector form. As can be seen, the chip phase exhibits close to a plug-flow behavior. In the upper part of the digester, the chip phase volume fraction is small, as a result, the cross section area that the chip phase flows through is small, and therefore the downward velocity of the chip phase is large; in the lower part of the digester, the chip phase volume fraction is large so the downward velocity of the chip phase is small. In Figure 6.9 (b) and (c), the detailed views of the chip phase velocity in two different locations are shown. Obviously, the chip phase has a slippage on the digester wall. 35IIIIUII Ul 'a X 20 1IIIUIIIIIIIUIHI iiimiiiiiuiaj •imiiiiiiiiu •mmiuuu l i imiui i iu iBl iiiuumtmuiB i n i m i i m i w i aiiiiiinuuiKi.i muuiu liiiiiimiiumiw wiiiiiiiiiiiiiiimw iiiiiiiiiiniiiiiiiii'ii iiiiiiiiiiiiiuiiiimii UWIIIIIIIIlllllllllll aumiiiuiiiiiiiiiiii WllUlllllllllllllllll UlliltlilHHIIllllllll mum milium lllllllli liniiiiii iiiiiiinii r, m 24 23 E 22 •J 21 a> T 20 19 18 wwwm in1111 ii i i i i inii i l l l l l i l i m u m III 11111 iinmiiiil 1 2 3 4 r, m (b) High chip pressure regime 4 | 3 E s £ 2 at 1 2 3 4 r, m (a) Overall chip phase velocity (c) Near the digester bottom Figure 6.9 Velocity distribution of the chip phase 114 Chapter 6 Simulation of an Industrial Continuous Digester In Figure 6.10, the predicted velocity profiles of the chip phase at three different digester levels are shown. As can be seen, the averaged downward velocity is decreased by about 20% from the top portion to the middle part of the digester; this is mainly due to an increased chip phase volume fraction as shown in Figure 6.6 (f). In the lower portion of the digester, the volume fraction of the chip phase further increases, and the mass of the chip phase decreases because of the increased delignification rate. As a result, the averaged downward velocity of the chip phase further decreases, being 0.2mm/sec lower at the bottom than at the middle of the digester. 1 — Top (Z = 35.1 m) - 2 . 5 - 2 — Middle (Z = 15.8 m) 3 — Bottom (Z = 5.0 m) L_i i i l i i i i l i i i i l i i i i LJ A 1.0 2 . 0 3 .0 4 .0 ' r, m Downcomer Figure 6.10 Chip phase velocity profiles at different digester levels (Z — digester level) In the upper portion of the digester, the chip phase velocities near the outer wall of the digester shell are slightly smaller than those in the center region, as indicated by curve 1 in Figure 6.10; in the middle part of the digester, this trend is augmented (curve 2); but in the lower portion of the digester, the velocity curve 3 becomes flat. This is due to the friction force between the chip phase and the digester wall. As plotted in Figure 6.5 (e), 115 Chapter 6 Simulation of an Industrial Continuous Digester the chip phase pressure is relatively high in the upper portion and becomes higher in the middle of the digester. Therefore the friction force at the digester shell retards the chips near the wall, causing the chip phase velocities to decrease in the vicinity of the wall. In the lower portion of the digester, the chip phase pressure is significantly decreased by the upflowing liquor and therefore the friction force at the wall is almost negligible. Further, chips near the digester shell gradually accelerate to the velocities of the chips in the middle entrained by the shear stress caused by the neighboring chips. Near the downcomer, the situation is essentially the same, but the surface area of the downcomer outer wall is much smaller than the digester wall. Consequently, the friction force is small and the velocity difference is minimal. 6.5.3.3 Kappa number variation at the bottom exit Figure 6.11 presents the prediction of kappa number distribution at the bottom exit of the digester. The kappa number near the downcomer and the digester wall are lower than that in the middle, indicating a non-uniform cook; the total kappa number variation is about 9. 24 I. i 11 i i 11 * i i 111 i i i 11 i i 111 i i 111 i i 11 11 i 11 i l 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 r, m Figure 6.11 Kappa number variation at the bottom exit of the digester 116 Chapter 6 Simulation of an Industrial Continuous Digester The cooking extent is determined by temperature, liquor concentration and residence time. Since all the circulating liquor was injected from the central downcomer, the concentration of effective alkali and temperature are high in the downcomer vicinity and low near the digester wall, as can be seen in Figure 6.6 (d) and (h). However, Figure 6.11 indicates that the chips near.the digester wall have lower kappa numbers and are thus over-cooked. Therefore the dominating factor of the non-uniform cook in this particular case is probably the residence time. To confirm this assumption, the total residence time of the chip phase is calculated based on its flow field and plotted against the radial locations; as shown in Figure 6.12, the chips near the digester wall have much longer residence time which allows further delignification and results in lower kappa number. 3 o £ 0) E o o c 0) •a '53 £ a. !E o 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 r, m Figure 6.12 Chip phase residence time versus chip locations 117 Chapter 6 Simulation of an Industrial Continuous Digester 6.6 Discussion of simulation results 6.6.1 Mesh convergence check To ensure a mesh converged solution, three meshes are used for the same cooking case and no obvious difference is found either in the resulting flow fields or the scalar distributions. For simplicity, only the downward velocities of the liquor phase along the vertical line r=2m are shown, as in Figure 6.13. As can be seen, the increased mesh size does not significantly influence the liquor phase velocity, and computations using a mesh of 40x200 should give reasonably good predictions. Liquor velocity, mm/s Figure 6.13 Computed downward liquor phase velocities along the vertical line r=2m using different meshes 6.6.2 Predictions of kappa number and yield To check the code's prediction capability, two-dimensional computations have been conducted in 15 available cooking cases, with corresponding cooking conditions listed in Appendix B. In these cases, the production rates range from about 1000 to 1547 tons air-dry-pulp per day; other operating parameters change from case to case. Operations in all these cases are fairly stable over a minimum period of 6 hours. In Table 6.10, the computed kappa number and yield, and the corresponding mill records are listed. As can be seen, the prediction errors of kappa number and yield are less than 10%. 118 Chapter 6 Simulation of an Industrial Continuous Digester Table 6.10 Computed kappa # and yield versus mill records No Production Rate (ADT/d) Chip Meter (RPM) kappa # (Record/Comp.) Error (%) Yield, (« (Record/ Comp) Vo) Error (%) 1 12686 15.06 29.63 /31.40 5.97 46.4/48.2 3.88 2 1313 15.01 27.99 / 29.05 3.79 48.2/50.7 5.19 -> 888 10.40 28.35 /30.86 8.85 46.7/48.7 4.28 4 1317 15.66 27.99/26.74 -4.47 46.4/48.2 3.21 5 1237 15.37 30.99/32.65 5.36 44.4 / 46.4 4.50 6 1223 14.95 29.35 /31.10 5.96 45.1 /47.8 5.99 7 1265 14.53 29.25 /31.31 7.04 48.0/47.2 -1.67 8 1410 14.89 28.26/29.80 5.45 52.2/49.0 -6.13 9 1128 13.79 35.66/33.95 -4.80 45.0/48.9 8.67 10 1100 12.85 40.94/37.36 -8.74 47.1 /49.1 4.25 11 937 11.07 29.42/30.38 3.26 46.4/48.9 5.38 12 1058 11.95 31.97/32.01 0.12 48.6/52.0 7.0 13 1398 14.79 32.27/31.71 -1.74 52.1 /50.3 -3.45 14 1084 13.33 32.29/29.06 -10.0 44.7/43.2 3.36 15 1190 15.38 35.07/33.13 -5.53 42.7/44.0 3.04 6.6.2 Variation of production rate In daily operations, the production rate of a digester may change as desired; in some pulp mills, digesters are often operated well above their original design capacities. To know how the variations of production rate affect the exit kappa number and yield is helpful for operational personnel. Figure 6.14 shows the simulated results with production rate increasing from 1075 to 1455 tons/day. As plotted, both kappa number and yield increase as production rate increases. In the computed production range, the kappa number increases from 25 to 43, while the yield increases from 47.5% to about 51.7%). The chip phase pressure at digester bottom exit drops from 8.8 to 4.7 kPa, as shown in Figure 6.14 (b). 119 Chapter 6 Simulation of an Industrial Continuous Digester 1100 1200 1300 1400 Production rate, ton/day (a) Yield and kappa number 1100 1200 1300 1400 Production rate, ton/day (b) Pc at the bottom exit plane Figure 6.14 Cooking performance affected by variation in production rate (Except production rates, data from HSPP, March 9, 2004) 6.6.3 Cooking performance affected by the chip type In the coastal area of British Columbia, many pulp mills use different types of wood chips or blended chips to meet the diverse needs of the pulp market. Shifting from one type of chips to another often causes significant variation in the exit kappa number. Without adequate knowledge, it is very difficult to make an appropriate and quick adjustment in order to meet the target kappa number. For different wood chips, the lignin content is different, the chip compressibility and flow resistance over different chips are also different, possibly leading to differences in cooking. Table 6.11 lists typical chip combinations used in HSPP; the corresponding contents of lignin and carbohydrates are computed based on data in Table 6.6. 120 Chapter 6 Simulation of an Industrial Continuous Digester Table 6.11 Typical combinations of incoming wood chips in HSPP Chip combination Hemlock (%) Cedar (%) Spruce (%) Fir (%) Ligo Carbo 1 60 40 21.16 70.78 2 45 55 23.80 71.28 3 50 50 24.05 70.25 4 40 60 24.54 70.32 5 50 50 25.60 69.55 6 30 70 27.20 69.41 Simulated cooking effects are shown in Figure 6.15 when the above six chip combinations are used. In these simulations, the same operating conditions from March 9, 2004 in HSPP are used; Harkonen's correlations [7] are used to calculate the chip compressibility and flow resistance. .a E 3 C re o. a. ns '>< in 40 38 36 34 32 30 28 26 24 22 20 • Kappa # A Yield(%) J I I I I L 2 3 4 5 6 Chip combinations Figure 6.15 Cooking performance affected by the chip type (See Table 6.11 for chip combinations) From Table 6.11 and Figure 6.15, the initial lignin content increases from chip combination 1 to 6, and the exit kappa number also increases, from about 24.5 to 33; the yield on the other hand, decreases from about 47.5 to 45. In past years, extensive work [7,29,30,31] has been published on chip compressibility. It is of interest to investigate the effect of chip compressibility on cooking performance. 121 Chapter 6 Simulation of an Industrial Continuous Digester In Chapter 2, the compressibility of different types of wood chips is listed; the concept of "equivalent diameter" characterizing the chip compressibility is also given. To address the effect of chip compressibility on the cooking performances, simulations have been carried out for different cases, where all the cooking conditions are the same as in the case from March 9, 2004. The six softwood chips listed in Table 2.2 are used with the equivalent diameter ranges from 0.85 to 3.0 mm. As seen in Figure 6.16 (a), when the equivalent diameter increases, both kappa number and yield increase at the beginning and decrease when the equivalent diameter exceeds about 1.4mm; the averaged residence time for wood chips exhibits the opposite trend as plotted in Figure 6.16 (b). Because longer residence time results in more complete cooking, lower kappa number of wood chips, the residence time is probably the immediate cause for kappa number changes. From Figure 6.15 and 6.16, variation in chip combination in HSPP and the equivalent diameter affect the kappa number almost to the same extent, with both leading to kappa number variation ranging from about 24.5 to 33. • K a p p a # • Y i e l d (%) • • • • 49 2 o >-1.0 1.5 2.0 2.5 3.0 Equivalent diameter, mm (a) kappa number and yield H 47 46 1.0 1.5 2.0 2.5 3.0 Equivalent diameter, mm (b) Averaged residence time Figure 6.16 Cooking affected by equivalent diameter of wood chips To illustrate how the chip residence time is affected by the equivalent diameter of wood chips, the chip phase pressure distributions corresponding to three different 122 Chapter 6 Simulation of an Industrial Continuous Digester equivalent diameters are shown in Figure 6.17; the chip phase volume fractions are given in Figure 6.18. When the equivalent diameter increases, the flow resistance over wood chips decreases. Consequently, the upflowing liquor decreases the solid pressure less. Therefore with bigger equivalent diameters in cases of Figure 6.17 (b) and (c), the chip phase pressure in the lower part of the digester is higher than that in the case of Figure 6.17(a). The volume fractions of the chip phase shown in Figure 6.18 are determined by the chip phase pressure and chip compressibility. Higher chip phase volume fraction means larger surface area that the chip phase can flow through, thus lower chip phase velocity. Compared with case (b), case (a) has higher chip phase volume fraction in the upper part of the digester, while case (c) has higher chip phase volume fraction in the lower part. As a result, chips in case (c) have the longest residence time, while chips in case (b) have the shortest residence time. 123 Chapter 6 Simulation of an Industrial Continuous Digester 1 2 3 4 1 2 3 4 1 2 3 4 r, m r, m r, m (a) D e q = 0 . 8 3 8 mm (b) D e q =l .394 mm (c) Deq=3.03 mm Figure 6.17 Chip phase pressure affected by the equivalent diameter . 2 5 "3 I 20 0.546 0.540 0.535 0.529 0.524 0.518 0.513 0.507 0.501 0.496 0.490 0.485 0.479 0.474 0.468 0.463 0.457 0.451 0.446 0.440 0.435 0.429 0.424 0.418 ffl X 2 0 | • 0.546 I 0.540 I 0.535 W 0.529 ^ 0.524 •-' 0.518 0.513 0.507 0.501 0.496 0.490 0.485 0.479 0.474 0.468 0.463 0.457 0.451 0.446 0.440 0.435 0.429 0.424 0.418 0.546 0.540 0.535 0.529 0.524 0.518 0.513 0.507 0.501 0.496 0.490 0.485 0.479 0.474 0.468 0.463 0.457 0.451 0.446 0.440 0.435 0.429 0.424 0.418 r, m r, m r, m (a) D e q = 0 . 8 3 8 mm (b) Deq=l.394 mm (c) D e q = 3 . 0 3 mm Figure 6.18 Chip phase volume fraction affected by the equivalent diameter 124 Chapter 6 Simulation of an Industrial Continuous Digester 6.6.4 Effect of the cold blow rate In HSPP, the cold blow rate is adjusted to maintain the free liquor level and control the pulp consistency (the percentage of fiber in the pulp solution) at the blow line. Figure 6.19 shows the variation of the kappa number and yield when the cold blow rate increases; in the simulations, the black liquor extraction increases accordingly so as to maintain steady state operations. re a. o. re 36 34 32 30 28 26 Kappa # Yield (%) J i L J i l_ 52 51 50 2 49 re a. *X 'S a> E & o .Q +J re a." 48 55 65 75 85 95 105 Cold blow rate, kg/s Cold blow rate, kg/s (a) Kappa number and yield (b) Chip phase pressure at bottom exit Figure 6.19 Cooking performance affected by the cold blow rate (Except cold blow rate, data from HSPP, March 9, 2004) As plotted in Figure 6.19 (a), both the kappa number and the yield increase with the cold blow rate; the chip phase pressure at the digester bottom decreases considerably, from about 7 to 2 kPa, as shown in Figure 6.19 (b). If the cold blow rate increases further, the chip phase pressure at the digester bottom may approach zero, which means wood chips lose contact with each other. In this situation, the whole chip column moves quite slowly or even stops, leading to the occurrence called "hang-up" that should be prevented. 125 Chapter 6 Simulation of an Industrial Continuous Digester 6.7 Three-dimensional simulation In large digesters, chips may plug the extraction screen when the extraction flow is suddenly changed or the flow rate is too high. As a remedy, the black liquor and modified cooking liquor extraction screens are divided into two parts and switched periodically from one part to the other. The extractions can be switched from either "up-and-down" or "left-and-right" mode. In the first mode, the upper half and the lower half of the screens are used alternatively; in the second mode, the extraction screens are divided into a number of sections and the liquor extractions are shifted in the circumferential direction. As a hypothetic case given in Figure 6.20, both the black liquor and modified cooking liquor extraction screens are divided into 12 equal sectors; the liquor extractions are shifted in the circumferential direction. In this situation, the uneven extraction of liquor may result in different compaction on the chip phase; it may also affect the temperature distribution near the extraction screens, and thus may contribute to a non-uniform cooking. To investigate this effect, a three-dimensional simulation is performed for the digester in HSPP. (a) Extraction at time T (b) Extraction at time T+At Figure 6.20 Liquor extraction shifted between screen sectors In the steady state computations, extractions using half of the screen are modelled; the dynamic behavior when extraction screens are shifted is not considered. Because of the symmetrical condition, the computation only considers one sector enclosed by OA 126 Chapter 6 Simulation of an Industrial Continuous Digester and OB, but the results in a cross section view are extended to the region enclosed by OA and OC, as shown in Figure 6.20 (b). In Figure 6.21, the structured mesh used in the computation is shown; the liquor flow configuration is also illustrated. (a) Overall mesh and liquor flow (b) Top view at black liquor extraction plane Figure 6.21 The structured mesh shown in reduced size (original size:30x20xl80) (BCL: bulk cooking liquor; BL: black liquor; MCL: modified cooking liquor; WL: wash liquor.) Again, the computation is performed for the operating case on March 9, 2004 from HSPP. Figure 6.22 presents the velocity and pressure distributions for the liquor phase. As can be seen in Figure 6.22 (b) and (c), higher velocity always corresponds to a larger pressure drop; when the circulation liquor is injected from the downcomer at high velocity, the liquor phase pressure drops significantly. Therefore, the total pressure drop in the liquor phase at the modified cooking liquor injection level is much larger than that at the black liquor extraction level. 127 Chapter 6 Simulation of an Industrial Continuous Digester (c) (c) Figure 6.22 Velocity and pressure distributions for the liquor phase (a) Front view in the centered symmetrical plane; (b) cross-section view at the black liquor extraction level (z=31.5m); (c) cross-section view at the modified cooking liquor injection level (z=21.95m). 128 Chapter 6 Simulation of an Industrial Continuous Digester Figure 6.23 shows the streamline in the liquor phase and the chip phase pressure. In Figure 6.23 (b), the chip phase pressure distribution in the central symmetrical plane is very similar to the distribution obtained in a two-dimensional simulation, as given in Figure 6.6 (d). From the streamline of the free liquor phase and the pressure contours in Figure 6.23 (c) and (d), the liquor velocity contributes greatly to the chip phase pressure distribution. Near the black liquor extraction screen, the extracted liquor tends to push the chip column outward, causing a high chip phase pressure near the extraction port. At the injection level of the modified cooking liquor, the free liquor moves down- and outward, resulting in a continuous rise in the chip phase pressure. The magnitude of the injection velocity is high but decreases quickly; therefore the chip phase pressure increases quickly near the injection port and increases slowly afterwards, as given in Figure 6.23 (d). (a) (b) Figure 6.23 Liquor phase streamline and the chip phase pressure (a), (b) velocity and pressure in the center symmetrical plane; (c), (d) pressure in cross-section view at the black liquor extraction level (z=31.5m) and the modified cooking liquor injection level (z=21.95m). 129 Chapter 6 Simulation of an Industrial Continuous Digester Figure 6.24 contains the chip phase velocity distributions. In Figure 6.24 (a), the chip phase exhibits a plug behavior as a whole, but near the digester wall, chip motion is retarded because of the friction force from the digester wall. Chip phase velocity profiles at the black liquor extraction level and wash screen level are given in Figure 6.24 (b). The chip phase velocity at the digester wall has a larger gradient at the black liquor extraction level than at the wash screen level, as the chip phase pressure is higher in the former case. 30 25 JZ CT '<U X 20 10 lilt ll»\ mi. . nil: 35 l^llll!IIU!IH. iiiuimnuA llllllllllimf iiiiiiiiiinimii Liilllllliianiii IIIIIIIIIIIIIIII .IIIIIIIIIIIIIIII iiiiiiiiiuiiii l l l l l l l l l lU l J I I I I M 11 i i r Illllll 111 nun I I I I I I I I I U I I I I I I I I I I I I I I I I I I I I I 15)111111111 iiniii m i n i m i m i n i 11 l imn I I I I I I I I I U I I I I I I I I I I I I I I U I I I I I I I I I I I I 11 i inii i I I I I I I I 11 I I I I I I I I I I I I I I 11 I I I I I I I l l i i i i i i 11 min i I I I I I I I I I I I I I I I I I I I I I I I 11 iiiiiui I I I I I I I i I I I I I I I I I I I I I I J I I I . I I I I , -4.0 1 2 3 4 r, m 1 — black liquor extraction (Z = 31.70 m) 2 — wash screen (Z =2.05 m) (a) (b) Figure 6.24 The velocity of the chip phase (a) velocity in the center symmetrical plane; (b) Axial velocity profiles at two different places. Figure 6.25 presents the distributions of volume fraction of the chip phase, temperature and kappa number. In the center symmetrical plane, the results in Figure 6.25 (1) are all close to those obtained in a two-dimensional computation, as given in Figure 6.6 (f), (g) and (h). The contour plots in the cross section views show that for compaction, temperature and kappa number, the difference in the radial direction is dominant while the difference in the circumferential direction is negligible. This indicates that the non-uniform extraction causes only local variations, and therefore is not a main factor in non-uniform cooking. 130 Chapter 6 Simulation of an Industrial Continuous Digester <°C) Kappa # '5 201 0.513 0.508 0.503 0.499 0.494 0.489 0.485 0.480 0.475 0.470 0.466 0.461 0.456 0.452 0.447 0.442 0.438 0.433 0.428 0.424 0.419 0.414 0.409 0.405 I 20 r, m (a) 161.9 161.5 161.1 160.7 160.4 160.0 160.0 159.6 159.3 158.9 158.5 158.1 157.8 157.4 157.0 156.6 156.3 155.9 155.5 155.2 154.8 148.0 132.3 122.2 102.4 r, m (b) (1) In the center symmetrical plane T(°C) r, m (c) Kappa # 96.32 95.81 95.29 94.76 94.27 93.75 93.24 92.72 92.21 91.70 91.18 90.67 90.15 (a) (b) (c) (2) Cross section view at the black liquor extraction level (z=31.5m) T(°C) 1.0 2.0 3.0 4.0 Kappa # 30.95 30.05 29.15 28.25 27.35 26.45 25.55 24.65 23.74 22.84 21.94 21.04 20.14 (b) (3) Cross section view at the bottom exit (c) Figure 6.25 Computed distributions of scalar variables (a) chip phase volume fraction; (b) temperature; (c) Kappa number. 131 Chapter 6 Simulation of an Industrial Continuous Digester 6.8 Summary This chapter presents the computational study on an industrial digester in Howe Sound Pulp and Paper Ltd. Firstly, the boundary conditions for both phases are described; boundary values, including inlet flow rates and species concentrations are determined. Two-dimensional simulations are then conducted for 15 data sets for stable operations of the digester in Howe Sound Pulp and Paper Ltd. The predicted kappa number and yield agree reasonably well with the mill records. The computed distributions of pressure, velocity, volume fraction etc. from a typical simulation are shown. The interactions between the free liquor and chip phases are discussed. Further efforts were aimed at the investigation of the influence of different parameters on the pulping process. The following was concluded: 1. Both yield and kappa number increase with production rate. However the chip phase pressure at the digester bottom exit decreases very significantly when the production rate increases. 2. Wood chips differ considerably in lignin contents but not in carbohydrates. When wood species changes, cooking conditions should be changed accordingly. Chips with higher lignin content require longer residence time. 3. Chip compressibility greatly affects the cooking performance. Based on simulations in typical cooking conditions, chips with an equivalent diameter of 1.4 mm have the shortest residence time in the digester and are least cooked. Smaller chips have higher compaction in the upper part, while bigger chips have higher compaction in the lower part of the digester, therefore both have a longer residence time and are more cooked, resulting lower yields and kappa numbers. 4. Increasing the cold blow results in higher kappa number and yield but lower chip phase pressure at the digester bottom. Because of this, a too high cold blow could result in "hanging chip column". 132 Chapter 6 Simulation of an Industrial Continuous Digester 6.8 Summary This chapter presents the computational study on an industrial digester in Howe Sound Pulp and Paper Ltd. Firstly, the boundary conditions for both phases are described; boundary values, including inlet flow rates and species concentrations are determined. Two-dimensional simulations are then conducted for 15 data sets for stable operations of the digester in Howe Sound Pulp and Paper Ltd. The predicted kappa number and yield agree reasonably well with the mill records. The computed distributions of pressure, velocity, volume fraction etc. from a typical simulation are shown. The interactions between the free liquor and chip phases are discussed. Further efforts were aimed at the investigation of the influence of different parameters on the pulping process. The following was concluded: 1. Both yield and kappa number increase with production rate. However the chip phase pressure at the digester bottom exit decreases very significantly when the production rate increases. 2. Wood chips differ considerably in lignin contents but not in carbohydrates. When wood species changes, cooking conditions should be changed accordingly. Chips with higher lignin content require longer residence time. 3. Chip compressibility greatly affects the cooking performance. Based on simulations in typical cooking conditions, chips with an equivalent diameter of 1.4 mm have the shortest residence time in the digester and are least cooked. Smaller chips have higher compaction in the upper part, while bigger chips have higher compaction in the lower part of the digester, therefore both have a longer residence time and are more cooked, resulting lower yields and kappa numbers. 4. Increasing the cold blow results in higher kappa number and yield but lower chip phase pressure at the digester bottom. Because of this, a too high cold blow could result in "hanging chip column". 132 Chapter 6 Simulation of an Industrial Continuous Digester In addition, a three-dimensional simulation is carried out. The computation shows that the uneven liquor extraction investigated and used in many mills does not cause significant non-uniformity in cooking. 133 Chapter 7 Summary and Conclusion Chapter 7 Summary and Conclusion In this thesis, a general computational code is developed to simulate the complex cooking process in industrial continuous digesters. In the following, the main contributions and conclusions are listed. A set of comprehensive mathematical equations governing the two-phase flow, heat transfer and chemical reaction in continuous digesters are presented; sub-models accounting for the inter-phase friction, thermal dispersion and chip compressibility are also discussed. In addition, a new, more reliable cooking reaction model is derived based on extensive experimental data. As a main advantage, this model can be used in various cooking conditions without the need to pre-determine some model constants, which in other models are dependent on operating conditions. Unlike most existing pulping models, this model abandons the piece-fitting method, thus parameters characterizing the cooking stages are no longer used, which is another attractive feature of this model. Despite its simple form, the model predictions agree well with a wide range of cooking experimental data. A three-dimensional two-phase CFD code based on a Finite Volume Method is then developed to calculate the flow fields for both phases; the volume fraction of solid phase, the concentrations of chemicals and the residual contents of lignin and carbohydrates are also calculated. The code is based on collocated curvilinear grids and is of second order accuracy; it has the capability of solving many different steady state laminar flow problems with or without the involvement of heat transfer or chemical reactions. The liquid flow in a laboratory digester with chips filled in is simulated using the developed code. Computations indicate that different flow patterns can be formed depending on the circulation rate. Simulations show that the compacted wood chips greatly 134 Chapter 7 Summary and Conclusion enhance the momentum dispersion of the liquor phase; therefore the liquid flow tends to be uniform away from the injection and extraction places. Simulations are performed for a full-scale industrial digester for numerous operating conditions. The predictions of yield and kappa number at the digester exit are fairly accurate for various cooking conditions. Further simulations lead to the following conclusions: 1. With increased production rate, yield and kappa number increase; however the chip phase pressure at the digester bottom exit decreases very significantly; wherever the chip phase pressure approaches zero, chips lose contact with each other and "hanging chip column" may occur. 2. Increased cold blow results in higher kappa number and yield, but lowers the chip phase pressure at the digester bottom. 3. Chip compressibility affects significantly the cooking performance. Chips with an equivalent diameter of around 1.4 mm exhibit a minimum residence time and are least cooked, therefore giving the highest yield and kappa number. A three-dimensional simulation is carried out to investigate the circumferential variation of different parameters. The computation shows that uneven liquor extraction investigated and used in some mills does not cause significant non-uniformity in cooking. The developed three-dimensional tool constitutes a reliable and robust tool to investigate industrial digesters. It is a tool that can be used to assist the design and optimization of the pulping process. 135 Chapter 8 Recommendation for Future Work The mathematical model, as well as the computational code developed in this thesis is a useful tool for the investigation and optimization of the pulping process. The tool can be further developed to improve its predictive capability and widen its use. Future work is suggested along the following lines: • An improved model accounting for the cold blow flow injected into the digester should be developed. The actual geometry is very complex and the simplification used in this thesis allows for good over-all simulation but might miss some local details if the injection is non -uniform or occurs with very high velocity, possibly leading to channeling. • The energy equation should be solved separately for the chip phase and the free liquor phase in the lower part of the digester. The temperature difference between the phases could be significant locally and the evaluation of this effect would be useful. • Chip compressibility and flow resistance are very important for the predictive capability of the model. Therefore further work is necessary for better definition of these two important parameters, especially when blended wood chips are used. • The model should be extended to batch digesters. They are widespread in industry, have been enlarged and incorporate new cooking techniques, and therefore an improved understanding of the process would be a useful contribution. 136 Bibliography 1. Michelsen, F.A., A Dynamic Mechanistic Model and Model-Based Analysis of a Continuous Kamyr Digester, PH.D. Thesis, Dept. of Eng. Cybernetics, The Norwegian Institute of Technology, University of Trondheim, 1995. 2. Rydholm, S., Pulping Processes, Interscience Publishers, New York, 1967. 3. Koch Christensen, P., Treforedlingskjemi, Department of Chemical Engineering, Norwegian Institute of Technology, Trondheim, 1991. 4. Harder, N . "Extended Delignification in Kraft Cooking." Svensk Paperstidning. 81(15):483, 1978. 5. He, P., Salcudean,M., Gartshore I., Bibeau E., Modeling of Kraft Two-Phase Digester Pulping Processes, TAPPI Conference, Anaheim, Sept. 1998. 6. Pougatch, K., Salcudean, M., Gartshore, I., "A numerical model of the reacting multiphase flow in a pulp digester", submitted to Applied Mathematical Modelling, (personal communications) 7. Harkonen, E.J. (1984), A Mathematical Model for Two-Phase Flow, Acta Polytechnica Scandinavia, Mechanical Engineering Series No. 88. 8. Michelsen, F. A., Christensen, R. T., Lunde, G. G., Lundman, G. and Johansson, K., Model Predictive Control of a Continuous Kamyr Digester at SCA-Nordliner Munksund, Sweden. Preprints Control Systems 1992, Dream vs. Reality: Modern Process Control in the Pulp and Paper Industry, Sept.-Oct.1992, Whistler, B.C. Canada, 1992. 9. Pougatch, K., Salcudean, M. , and Gartshore, I., "Mathematical Modelling of Digesters", presented at 52nd Canadian Chemical Engineering Conference, Vancouver, October 20-23, 2002. 137 10. Miyanishi, T. and Shimada, H., Improvement of Pulp Strength and Yield by Computer Simulation of Lo-solids Kraft Cooking, 2001 Tappi Journal Peer Reviewed Paper. 11. Vroom, K.E., The " H " factor: A Means of Expressing Cooking Times and Temperatures as a Single Variable, Pulp and Paper Magazine of Canada, 58(3): 228 (1957). 12. Kerr, A.J., The Kinetics of Kraft Pulping, Appita 24(3): 180-188,1970. 13. Kerr, A.J. and Uprichard, J.M., The Kinetics of Kraft Pulping - Refinement of a Mathematical Model. Appita, 30:48-54, 1976. 14. Gustafson, R., Sleicher,C, Mckean,W., et. al., Theoretical Model of the Kraft Pulping Process, Ind. Eng. Chem. Process Design & Development, 22(l):87-96, 1983. 15. Christensen, T., Albright, L. and Williams,T., TAPPI 1983, Annual Meeting, Atlanta, Georgia, March 2-4, 239(1983). 16. McKibbens, S.W., Application of Diffusion Theory to the Washing of Kraft Cooked Wood Chips. TAPPI Journal, 43(10):801-805, 1960. 17. Scheiddegger, A.E., Statistical hydrodynamics in porous media, J. Appl. Phys. 25(8), 994-1001, 1954. 18. Spitz, K. and Moreno J., A practical guide to groundwater and solute transport modeling, ISBN: 0-471-13687-5, by John Wiley & Sons. Inc., New York, 1996. 19. Huang, H., Li, Z. and Usmani,A.S., Finite element analysis of non-Newtonian flow: theory and software, ISBN: 1852330244. London ; New York : Springer, 1999. 20. Biermann, C.J., Essentials of Pulping and Papermaking, ISBN 0-12-097360-X, Academic Press, 1993. 21. Gullichsen, J., Chemical Pulping, Papermaking Science and Technology, Book 6A, ISBN952-5216-06-03, Fapet Oy, Finland, 1999a. 138 22. Macdonald, I.F., et al., "Flow through Porous Media - the Ergun Equation Revisited", Ind. Eng. Chem. Fundam., Vol.18, No.3, 1979, pi99-208. 23. Ergun, S. "Flow through packed columns", Chem. Eng. Progress., Vol.48, No.2, p89, 1952. 24. N . Ahmed and D. K. Sunada, Non-linear Flow in Porous Media, Journal of Hydr Division, ASCE, Vol.95, No.6, 1969, pi847. 25. N . Ahmed and D. K. Sunada, Non-linear Flow in Porous Media, Journal of Hydr Division, ASCE, Vol.97, 1971, pl233. 26. Barnes, H.A., Hutton J. F. and Walters, K., An Introduction to Rheology, ISBN: 0444871403, Elsevier Science Publishers B.V., 1989. 27. Cross, M.M., Rheology of non-Newtonian fluids: a new flow equation for pseudo-plastic systems, J. Colloid Sci., (20), 417-437. 28. Masse, M.A., Kiran, E., and Fricke, A.L., A Thermodynamical Model of the Heat Capacity of Compositionally Complex Multicomponent Polymer Solutions: Kraft Black Liquor, Chemical Engineering Communications, 50: 81-91 (1987). 29. Wang, Z. and Gullichsen, J., Flow Compressibility and Resistance of Chips Made with a new Chipping Technique, Pro. 1998 TAPPI Pulping Conf, Montreal (Oct25-29), Vol2, p.753-760. 30. Lindqvist, J., The Effects of Flow Conditions in Kraft Cooking of Wood Chips by Various Chip Size Distribution, Masters Thesis, Helsinki University of Technology, 1994. 31. Lammi, L., The Utilization of Eucalyptus Camaldulensis-species with the Superbatch process, Masters Thesis, Helsinki University of Technology, 1996. 32. Dudgeon, C. R., An Experimental Study of Flow of Water through Coarse Granular Media, La Houllie Blanche, France, Vol.21, No.7, p.785, 1966. 139 33. Gupte, A. R., Dissertation, Fluid-Dynamic characteristics of two-phase flow, Karlsruhe, 1970. 34. Nam-Trung Nguyen, Steven T. Wereley, Fundamentals and applications of microfluidics, ISBN: 1580533434, Artech House, Boston, MA, c2002. 35. Thompson, P. A., and Troian, S.M., "A general boundary condition for liquid flow at solid surfaces", Nature, Vol.389, 1997, p. 360-362. 36. Edwards, L. and Norberg, S.E., Alkaline delignification kinetics, A General Model Applied to Oxygen Bleaching and Kraft Pulping, TAPPI J. 56(11): 108-111 1973. 37. Kubes, G.J., Fleming, B.I. and MacLeod, J.M., Viscosities of Unbleached Alkaline Pulps, II. The G factor, J.Wood Chem. Tech. 3(3): 313,1983. 38. Paulonis, M.A. et. al, Adaptive Inferential Control of Kraft Batch Digesters as Based on Pulping Liquor Analysis, TAPPI J. 74(6): 169-175,1991. 39. Smith, C.C. and Williams T.J., A Dynamic Model of a Kamyr Digester and its Application to Pulp Washing and Pulping Reaction Control, Modeling & Control Kraft Production Systems. Proc. ISA Symp. Milwaukee, Oct., 1975. 40. Smith, C.C. and Williams, T.J., "Mathematical Modeling, Simulation, and Control of the Operation of a Kamyr Continuous Digester for the Kraft Process", Report No.64, Purdue Laboratory for Applied Industrial Control, Purdue University, Lafayette, Indiana, 1974. 41. Aurora Santos, Francisco Rodriguez et. al, Kinetic Modeling of Kraft Delignification of Eucalyptus Globulus, Ind. Eng. Chem. Res, 1997, 36, 4114-4125. 42. Ling Yang, Keng Chuang and Shijie Liu, A mechanistic Model for Kraft Pulping, PAPTAC 90th Annual Meeting, 2004. 140 43. Hairer E., Lubich, C. and Roche M., The numerical solution of differential-algebraic systems by Runge-Kutta methods, ISBN: 3540518606, Berlin; New York: Springer-Verlag, cl989. 44. Aurell, R., Harder, N. , Kraft pulping of pine. I. Changes in the compounds of the wood residue during the cooking process. Sven. Papperstidn. 1965, 68 (3), 59-68. 45. Kuo, K. K , Principles of Combustion, Chapter 2, pl09, ISBN 0-471-09852-3, John Wiley & Sons, Inc., 1986. 46. Bugajer, S. Kinetics of delignification reactions in the Eucalyptus pulping process (in Portuguese). Ph.D. Thesis, Polytechnic School, University of Sao Paulo, Brazil, 1984. 47. Matthews, C.H., Carbohydrate losses at high temperature in kraft pulping. Sven. Papperstidn. 1974, 77(17), 629-635. 48. Kleinert, T.N., Mechanisms of alkaline delignification. I. The overall reaction patern. Tappi 1966, 49 (2), 53-57 49. Wilder, H.D., Daleski, E.J., Kraft Pulping Kinetics. I. Literature Review and Research Program. Tappi 1964, 47 (5), 270-275. 50. Wilder, H.D., Daleski, E. J. Delignification rate studies. Part II of a series on kraft pulping kinetics. Tappi 1965, 48 (5), 293-297. 51. Chiang, V.L., Puumala, R.J., Takeuchi, H., Eckert, R. E., Comparison of softwood and hardwood kraft pulping. Tappi J.1988, 71 (9), 173-176. 52. Chorin, A.J., A numerical method for solving incompressible viscous flow problem, J. Comput. Phys., (2) 12, 1967. 53. Patankar Suhas V., Numerical Heat Transfer and fluid flow, Hemisphere Pub. Corp.; New York, McGraw-Hill, P37, 1980. 141 54. Harlow, F.H., Welsh, J.E., Numerical Calculation of Time Dependent Viscous Incompressible Flow with Free Surface. Phys. Fluids, (8) 2181-2189,1965. 55. Ferziger, J.H. and Peric, M. , Computational methods for fluid dynamics, second Edition, Springer, 1999. ISBN3-540-65373-2. QA911 .F434 1999 c l 56. Peric, M. , et al. Comparison of finite-volume numerical methods with staggered and collocated grids, Computers and Fluids, 16(4):389-403, 1998 57. Stone, H.L., Iterative solution of implicit approximations of multi-dimensional partial differential equations, SIAM Journal of Numerical Analysis, (5) 530-558, 1968. 58. Hageman, L.A. and Young, D.M., Applied iterative methods, Wiley, New York, 1981. 59. Muzaferija, S., Adaptive finite volume method for flow predictions using unstructured meshes and multigrid approach. PhD Thesis, Univeristy of London, 1994. 60. Chen, C.J., Bernatz, R., et al, Finite Analytic Method in Flows and Heat Transfer, ISBN 1-56032-898-3, Taylor and Francis, 2000. 61. Hortmann, M. , Peric, M. , Scheurer, G., Finite Volume Multigrid Prediction of laminar natural convection: bench-mark solutions. International Journal of Numerical Methods in Fluids, No.l 1, 189-207, 1990. 62. Skelland, A.H.P, Non-Newtonian Flow and Heat Transfer, John Wiley & Sons, Inc. New York, 1967. 63. Lauriat G. and Vafai K., Forced convective flow and heat transfer through a porous medium exposed to a flat plate or a channel, Convective and Mass Transfer in Porous Media, p289-327, 1991. 64. Vlaev, D.S. and Bennington, C , Using Electrical Resistance Tomography to Image Liquor Flow in a Model Digester, to be published. 142 65. Vlaev, D.S. and Bennington, C , Flow Uniformity in a Model Digester Measured with Electrical Resistance Tomography, 3 r d World Congress on Industrial Process Tomography, Banff, Canada. 66. Tapp, H.S. and Wilson R.H., Developments in low-cost electrical imaging techniques, Process Control Quality 9 (1-3) (1997) 7-16. 67. Dickin, F., Wang, M. , Electrical resistance tomography for process applications, Meas. Sci. Technol. 7 (3) (1996) 247-260. 68. Ma Y., Zheng Z. et al., Application of electrical resistance tomography system to monitor gas/liquid two-phase flow in a horizontal pipe, Flow Measurement and Instrumentation 12, (2001) 259-265 69. Peyton A. J., Yu, Z. Z., Lyon G., AlZeibak S., et al., An overview of electromagnetic inductance tomography: description of three different systems, Meas. Sci. Technol. 7 (3) (1996) 261-271. 70. Reinecke, N . and Mewes, D., Recent developments and industrial/research applications of capacitance tomography, Meas. Sci. Technol. 7 (3) (1996) 233-246. 71. Kenneth E. Smith, Process control for pulp and paper mills, Miller Freeman Publications, ISBN: 087930149X, 1983. 72. Graham Vandegriend, pure species LWAFL and coarseness values, ECOLO Talk, No.4, 1999. 73. Wood Handbook: Wood as an Engineering Material, Agriculture Handbook 72, U.S. Department of Agriculture, Forest Service, Forest Products Laboratory, Madison, Wise. Revised, 1987. 74. Robert H. White and Erik V. Nordheim, Charring Rate of Wood for ASTM E 119 Exposure, Fire Techology,Vol.28, No.l , 1992. 75. R. Pettersen, The Chemical Composition of Wood, in R. Rowell, ed. The Chemistry of Solid Wood. American Chemical Society, Washington, D. C , 1984. 143 76. Chen H., Harmon, M. E., Sexton J., and Fasth, B., Fine-root decomposition and N dynamics in coniferous forests of the Pacific Northwest, U.S.A., Can. J. For. Res. 32: 320-331 (2002). 77. Blum L., The Production of Bleached Kraft Pulp, Environmental Defence Fund, 1996. 144 Appendix A Mill record of digester operation at HSPP* A B C D E F G H I J K L Tagname 411AI9129 411AT9606 411AYI9125 A 411AYI9125 B 411AYI9125 C 411AYI9125 D 411AYI9125 E 411AYI9125 F 411CT9036 411DRPM 411FI9121 Engineering Units KAPPA g/L H-Fctr H-Fctr H-Fctr H-Fctr H-Fctr H-Fctr G/L RPM US Description Dig Blowline Kappa BC RESIDUAL ALKALI H-Factor to Wash H-Factor to Mid EMCC H-Factor to End MCC H-Factor to Mid MCC H-Factor to Extraction H-Factor in IV WL STRENGTH CHIP METER RPM DIG BLOW FLOW(F12A/B) 9-Mar-2004 15:00 30.1 12.8 1494.0 1255.4 852.8 645.7 355.8 0.0 84.2 15.7 133.8 9-Mar-2004 16:00 29.9 12.8 1513.2 1228.1 826.5 656.0 382.8 0.0 84.4 15.3 134.0 9-Mar-2004 17:00 29.1 12.8 1467.4 1210.9 862.6 703.3 409.7 0.0 84.5 15.4 137.7 9-Mar-2004 18:00 29.7 12.8 1500.6 1264.7 914.8 713.2 375.5 0.0 84.6 13.4 130.9 9-Mar-2004 19:00 29.4 12.8 1565.6 1317.1 883.3 713.3 420.8 0.0 84.8 15.4 129.3 9-Mar-2004 20:00 29.6 12.8 1568.7 1268.1 905.2 662.9 336.6 0.0 84.9 15.2 125.6 Min 29.1 12.8 1467.4 1210.9 826.5 645.7 336.6 0.0 84.2 13.4 125.6 Avg 29.6 12.8 1518.2 1257.4 874.2 682.4 380.2 0.0 84.6 15.1 131.9 Max 30.1 12.8 1568.7 1317.1 914.8 713.3 420.8 0.0 84.9 15.7 137.7 (max-min)/avg 23% 3% 1% 7% 8% 10% 10% 22% 1% 15% 9% 11-Mar-2004 14:00 28.4 12.8 1423.9 1173.2 794.1 625.7 349.9 0.0 88.6 14.8 130.8 11-Mar-2004 15:00 27.9 13.1 1424.3 1168.6 791.5 629.8 362.1 0.0 87.1 14.7 136.5 11-Mar-2004 16:00 28.1 12.3 1405.1 1150.2 805.1 647.9 375.6 0.0 87.1 15.4 120.8 11-Mar-2004 17:00 27.1 12.4 1387.3 1154.0 832.5 645.6 353.7 0.0 87.1 15.3 128.5 11-Mar-2004 18:00 27.7 12.5 1410.2 1181.9 800.9 633.7 373.1 0.0 87.0 14.8 141.7 11-Mar-2004 19:00 28.7 12.3 1421.3 1162.0 823.7 659.4 388.5 0.0 87.0 14.9 142.7 Min 27.1 12.3 1387.3 1150.2 791.5 625.7 349.9 0.0 87.0 14.7 120.8 Avg 28.0 12.6 1412.0 1165.0 808.0 640.3 367.2 0.0 87.3 15.0 133.5 Max 28.7 13.1 1424.3 1181.9 832.5 659.4 388.5 0.0 88.6 15.4 142.7 (max-min)/avg 18% 5% 6% 3% 3% 5% 5% 11% 2% 4% 16% *HSPP is an abbreviation of "Howe Sound Pulp and Paper". 145 Appendix A Mill record of digester operation at HSPP (continued) A M N O P Q R S T u V W Tagname 411FT9011 411FT9021 411FT9031 411FT9032 411FT9037 411FT9038 411FT9039 411FT9041 411FT9111 411FT9141 411FT9160 Engineering Units T/H T/H L/S L/S L/S L/S L/S L/S L/S L/S L/S Description HI PRESS STEAM FLOW LP STEAM FLOW F3A HPF WL PURGE F 3EWL FLOW TO WASH F3B WL FLOW TO BC F 3D WL FLOW TO MC F3F WL F4 CB FLOW TO IV FLOW TO IV F11 CB FOR F14CBT O DIG PRES DIG OD F16AEXTR SCRN FLOW 9-Mar-2004 15:00 40.0 27.4 3.9 5.0 2.5 4.9 48.3 0.4 50.0 11.6 151.0 9-Mar-2004 16:00 36.3 28.1 3.9 4.9 1.3 4.9 47.1 0.4 46.4 10.4 136.6 9-Mar-2004 17:00 31.7 27.1 3.9 4.7 1.5 4.7 47.2 0.4 47.5 10.8 133.2 9-Mar-2004 18:00 32.1 24.0 3.9 4.5 1.7 4.5 40.5 0.4 47.6 11.0 125.2 9-Mar-2004 19:00 36.8 29.7 3.9 4.5 2.1 4.5 47.1 0.4 49.0 11.1 148.5 9-Mar-2004 20:00 33.3 27.8 4.0 4.5 1.4 4.5 45.8 0.4 46.3 10.6 142.3 Min 31.7 24.0 3.9 4.5 1.3 4.5 40.5 0.4 46.3 10.4 125.2 Avg 35.0 27.3 3.9 4.7 1.7 4.7 46.0 0.4 47.8 10.9 139.5 Max 40.0 29.7 4.0 5.0 2.5 4.9 48.3 0.4 50.0 11.6 151.0 (max-min)/avg 23% 23% 21% 1% 12% 10% 17% 9% 8% 11% 18% 11-Mar-2004 14:00 38.7 26.1 3.5 4.7 4.4 4.7 43.7 0.3 50.8 5.7 146.1 11-Mar-2004 15:00 35.7 25.9 3.6 4.8 0.0 4.7 44.3 0.3 52.8 5.8 137.4 11-Mar-2004 16:00 32.9 26.9 3.5 4.0 0.0 4.0 46.6 0.3 47.5 5.6 144.2 11-Mar-2004 17:00 39.2 27.5 3.5 4.0 0.0 4.0 46.3 0.3 50.9 6.0 145.2 11-Mar-2004 18:00 38.6 25.7 3.5 4.5 0.0 4.5 44.7 0.3 54.6 5.8 136.9 11-Mar-2004 19:00 38.7 26.4 3.4 4.8 0.0 4.8 45.0 0.3 55.4 6.4 140.3 Min 32.9 25.7 3.4 4.0 0.0 4.0 43.7 0.3 47.5 5.6 136.9 Avg 37.3 26.4 3.5 4.4 0.8 4.5 45.1 0.3 52.0 5.9 141.7 Max 39.2 27.5 3.6 4.8 4.4 4.8 46.6 0.3 55.4 6.4 146.1 (max-min)/avg 18% 17% 7% 4% 18% 18% 6% 0% 15% 13% 6% 146 Appendix A Mill record of digester operation at HSPP (continued) A X Y Z AA AB AC AD A E A F AG AH Tagname 411FT9181 411FYI9025 411FYI9121 411HI9022B 411HT9022B 411LY9102 411LY9752 411MI9022 411NIC9183 411TT9085 411TT9165 Engineering Units L/S T/H ADT/D FURNSH N/A % % % % BD DEG C DEG C CALC STEAM DIGESTER DIGESTER TRIM LQR Description F18 C B T O FLOW TO DIGESTER EXIT PRODUCTIO CHIP IV CHIP WEIGHTED PD18 BLOWTEMPERAT EXTRACTI DIG BOT BIN BLOW RATE FURNISH N GRADE LEVEL LEVEL AVE MOIST CONS MV URE ON TEMP 9-Mar-2004 15:00 41.8 29.2 1286.4 400.0 400.0 29.3 8.8 54.0 10.0 157.9 160.3 9-Mar-2004 16:00 38.3 30.1 1289.6 400.0 400.0 29.8 31.1 54.4 10.0 158.8 160.6 9-Mar-2004 17:00 39.0 28.6 1325.2 400.0 400.0 33.2 27.3 54.2 10.0 158.6 161.6 9-Mar-2004 18:00 39.2 26.1 1259.1 400.0 400.0 24.3 16.5 54.0 10.0 158.3 161.2 9-Mar-2004 19:00 40.7 30.6 1238.5 400.0 400.0 27.8 3.9 54.4 10.0 157.2 161.3 9-Mar-2004 20:00 38.2 28.5 1206.6 400.0 400.0 24.0 19.9 54.5 10.0 156.3 159.7 Min 38.2 26.1 1206.6 400.0 400.0 24.0 3.9 54.0 10.0 156.3 159.7 Avg 39.5 28.8 1267.6 400.0 400.0 28.1 17.9 54.3 10.0 157.9 160.8 Max 41.8 30.6 1325.2 400.0 400.0 33.2 31.1 54.5 10.0 158.8 161.6 (max-min)/avg 23% 9% 16% 9% 0% 0% 1% 1% 2% 1% 11-Mar-2004 14:00 40.0 27.6 1299.8 400.0 400.0 40.0 33.7 54.5 10.3 156.5 159.8 11-Mar-2004 15:00 41.7 27.7 1346.6 400.0 400.0 39.6 35.0 54.5 10.3 157.4 159.6 11-Mar-2004 16:00 37.1 28.3 1189.8 400.0 400.0 24.8 14.0 54.5 10.2 157.5 160.2 11-Mar-2004 17:00 39.4 29.1 1264.8 400.0 400.0 33.0 16.7 54.5 10.2 157.4 159.8 11-Mar-2004 18:00 42.8 27.6 1378.4 400.0 400.0 42.2 33.4 54.5 10.1 157.4 159.8 11-Mar-2004 19:00 43.5 28.3 1397.7 400.0 400.0 38.3 27.4 54.5 10.2 157.7 160.2 Min 37.1 27.6 1189.8 400.0 400.0 24.8 14.0 54.5 10.1 156.5 159.6 Avg 40.8 28.1 1312.9 400.0 400.0 36.3 26.7 54.5 10.2 157.3 159.9 Max 43.5 29.1 1397.7 400.0 400.0 42.2 35.0 54.5 10.3 157.7 160.2 (max-min)/avg 18% 16% 5% 16% 0% 0% 0% 2% 1% 0% 147 Appendix A Mill record of digester operation at HSPP (continued) A A l A J AK AL AM AN AO AP AQ AR AS Tagname 411TT9195 411TT9197 411TT9205 411TT9207 411TT9600 411TT9605 411TT9607 411TT9700 411TT9705 411TT9710 411TT9711 Engineering Units DEG C DEG C DEG C DEG C DEG C DEG C DEG C DEG C DEG C DEG C DEG C Description T19C MC COLD TEMP T19H MC HOT TEMP T20C WASH COLD TEMP T20H WASH HOT TEMP T60C BC COLD TEMP T60H BC HOT TEMP T60L BC H TEMP SP IV UPPER TEMPERATURE IV LOWER TEMPERATURE TRIM ZONE TEMP EXTR ZONE 1 TEMP 9-Mar-2004 15:00 161.9 161.0 150.3 159.1 158.4 177.4 158.3 122.7 129.6 157.2 160.8 9-Mar-2004 16:00 161.8 161.0 153.6 159.1 159.2 177.6 159.0 125.7 129.7 157.8 160.5 9-Mar-2004 17:00 161.7 161.2 157.4 159.2 158.0 175.4 157.4 125.3 130.0 157.9 160.6 9-Mar-2004 18:00 161.4 160.1 152.1 159.4 159.4 174.5 158.8 124.1 130.2 157.7 160.1 9-Mar-2004 19:00 161.8 161.2 149.6 158.9 156.7 174.2 156.4 123.1 130.1 158.1 160.0 9-Mar-2004 20:00 161.4 160.7 156.2 158.9 158.5 176.0 158.3 123.3 129.8 157.3 160.3 Min 161.4 160.1 149.6 158.9 156.7 174.2 156.4 122.7 129.6 157.2 160.0 Avg 161.6 160.9 153.2 159.1 158.4 175.9 158.0 124.0 129.9 157.6 160.4 Max 161.9 161.2 157.4 159.4 159.4 177.6 159.0 125.7 130.2 158.1 160.8 (max-min)/avg 23% 0% 1% 5% 0% 2% 2% 2% 2% 0% 1% 0% 11-Mar-2004 14:00 161.1 160.4 147.0 158.6 157.0 176.2 157.0 123.6 129.2 156.4 160.1 11-Mar-2004 15:00 160.7 159.9 148.6 158.2 158.0 176.1 157.9 124.6 129.8 157.0 159.6 11-Mar-2004 16:00 160.7 160.1 150.2 157.4 157.7 175.3 157.2 123.0 129.8 157.3 159.1 11-Mar-2004 17:00 160.5 160.0 142.8 157.8 157.8 176.3 157.2 123.1 129.1 157.0 159.4 11-Mar-2004 18:00 160.6 160.0 144.0 158.1 158.2 176.7 157.8 123.1 128.7 157.1 159.7 11-Mar-2004 19:00 160.5 159.9 143.8 158.2 158.3 176.7 157.8 124.1 128.8 157.5 159.4 Min 160.5 159.9 142.8 157.4 157.0 175.3 157.0 123.0 128.7 156.4 159.1 Avg 160.7 160.1 146.1 158.1 157.9 176.2 157.5 123.6 129.2 157.0 159.5 Max 161.1 160.4 150.2 158.6 158.3 176.7 157.9 124.6 129.8 157.5 160.1 (max-min)/avg 18% 0% 0% 5% 1% 1% 1% 1% 1% 1% 1% 1% 148 Appendix A Mill record of digester operation at HSPP (continued) A AT AU AV AW AX AY AZ BA BB BC BD BE Tagname 411TT9712 411TT9713 411TT9714 411TT9715 411UI9160 411WI9022411XI9022A411XI9022B 411X19022 C 411X19022 D 411FT9601 411FT9611 Engineering Units DEG C DEG C DEG C DEG C % T/DBD % % % % L/S L/S Description EXTR ZONE 2 TEMP EXTR ZONE 3 TEMP MCC ZONE TEMPERAT URE WASH ZONE TEMP CALC DILU FACTOR TOTAL CHIPS TO CB HEMLOCK CEDAR SPF FIR FURNISH FURNISH FURNISH FURNISH FRACTION FRACTION FRACTION FRACTION F60 BC FLOW F61A UPPER SLUICE FL 9-Mar-2004 15:00 161.0 161.1 161.9 160.8 0.0 2523.9 0.0 42.6 57.3 0.0 156.3 115.0 9-Mar-2004 16:00 160.6 160.9 161.7 162.0 -0.7 2659.5 0.0 47.3 52.7 0.0 152.8 125.3 9-Mar-2004 17:00 160.8 161.0 162.9 161.0 -0.8 2348.0 0.0 46.4 53.7 0.0 150.6 126.1 9-Mar-2004 18:00 160.7 160.7 161.6 160.0 -0.5 2253.2 0.0 30.5 68.7 0.0 148.2 125.8 9-Mar-2004 19:00 160.5 160.6 161.3 160.7 0.1 2481.0 0.0 47.2 52.7 0.0 147.4 117.8 9-Mar-2004 20:00 160.7 160.6 162.5 160.6 -0.2 2302.7 0.0 46.5 55.1 0.0 147.2 119.4 Min 160.5 160.6 161.3 160.0 -0.8 2253.2 0.0 30.5 52.7 0.0 147.2 115.0 Avg 160.7 160.8 162.0 160.9 -0.4 2428.1 0.0 43.4 56.7 0.0 150.4 121.6 Max 161.0 161.1 162.9 162.0 0.1 2659.5 0.0 47.3 68.7 0.0 156.3 126.1 (max-min)/avg 23% 0% 0% 1% 1% 17% 6% 9% 11-Mar-2004 14:00 160.3 160.3 160.3 160.4 0.1 2304.8 0.0 45.9 53.1 0.0 138.4 121.1 11-Mar-2004 15:00 159.9 159.8 160.1 161.5 0.0 2338.5 0.0 46.5 53.3 0.0 136.5 122.3 11-Mar-2004 16:00 159.2 159.5 160.0 160.8 0.1 2608.2 0.0 47.1 52.9 0.0 135.6 126.0 11-Mar-2004 17:00 159.4 159.6 160.0 159.9 0.1 2558.6 0.0 47.0 53.9 0.0 135.3 120.0 11-Mar-2004 18:00 159.7 159.9 160.5 154.2 0.0 2175.3 0.0 47.0 53.4 0.0 136.4 120.2 11-Mar-2004 19:00 159.7 159.4 159.9 156.0 0.1 2213.6 0.0 47.2 52.8 0.0 136.9 122.3 Min 159.2 159.4 159.9 154.2 0.0 2175.3 0.0 45.9 52.8 0.0 135.3 120.0 Avg 159.7 159.7 160.1 158.8 0.1 2366.5 0.0 46.8 53.2 0.0 136.5 122.0 Max 160.3 160.3 160.5 161.5 0.1 2608.2 0.0 47.2 53.9 0.0 138.4 126.0 (max-min)/avg 18% 1% 1% 0% 5% 18% 2% 5% 149 Appendix A Mill record of digester operation at HSPP (continued) A BF BG BH B l B J BK BL BM BN BO BP BQ Tagname 411FT9191 411FT9201 411FT9614 411TT9125 411TT9245 411NT9126 411IT0024 411PDT9183 411FT9871 411AT9158 411AT9190 421AI9129 Engineering Units L/S L/S L/S DEG C DEG C % A KPA L/S g/i g/i KAPPA Description F19AMC UPPER FLOW WASH CIRC FLOW CONT BC Heater Bypass Flow DIF INLET TEMP CB COOLER TEMP DIG BLOW CONSIST OUTLET DEVICE (DIG) PD18 BLOW BLOW LINE CONSIST DILN Extraction Resid Alkali MC Resid Alkali Decker Vat Kappa 9-Mar-2004 15:00 120.1 95.7 13.1 80.7 56.6 10.9 372.2 83.4 0.0 10.3 14.2 32.4 9-Mar-2004 16:00 119.9 95.6 13.7 83.8 55.5 10.9 373.2 81.8 0.1 10.3 14.3 30.9 9-Mar-2004 17:00 120.2 95.5 15.5 81.4 53.0 10.7 343.0 82.6 0.0 10.5 14.2 31.5 9-Mar-2004 18:00 120.1 95.2 16.9 81.9 53.3 10.9 384.0 82.7 0.0 10.5 14.1 32.3 9-Mar-2004 19:00 120.1 95.2 18.4 81.2 57.7 10.9 374.5 83.0 0.0 10.2 14.1 31.4 9-Mar-2004 20:00 120.0 94.8 15.5 81.2 57.7 10.8 349.1 80.9 0.0 10.5 14.0 31.3 Min 119.9 94.8 13.1 80.7 53.0 10.7 343.0 80.9 0.0 10.2 14.0 30.9 Avg 120.1 95.3 15.5 81.7 55.6 10.9 366.0 82.4 0.0 10.4 14.2 31.6 Max 120.2 95.7 18.4 83.8 57.7 10.9 384.0 83.4 0.1 10.5 14.3 32.4 (max-min)'avg 23% 0% 1% 4% 8% 2% 11% 3% 3% 2% 5% 11-Mar-2004 14:00 120.1 89.4 46.6 82.1 58.9 11.4 410.1 93.8 0.0 11.6 13.6 28.1 11-Mar-2004 15:00 120.2 89.2 44.1 82.2 58.1 11.1 409.4 92.3 0.0 11.4 14.6 28.0 11-Mar-2004 16:00 120.1 89.3 42.3 80.7 61.0 10.9 399.4 90.4 0.0 10.6 14.0 27.4 11-Mar-2004 17:00 119.9 89.4 44.0 83.3 58.5 11.0 410.8 90.6 0.0 10.3 13.8 27.4 11-Mar-2004 18:00 119.9 89.1 45.5 82.2 55.6 10.8 394.6 87.1 0.6 10.3 13.7 28.7 11-Mar-2004 19:00 120.1 89.1 43.3 81.5 58.6 10.9 407.5 89.5 0.0 10.3 14.0 31.3 Min 119.9 89.1 42.3 80.7 55.6 10.8 394.6 87.1 0.0 10.3 13.6 27.4 Avg 120.1 89.3 44.3 82.0 58.5 11.0 405.3 90.6 0.1 10.7 13.9 28.5 Max 120.2 89.4 46.6 83.3 61.0 11.4 410.8 93.8 0.6 11.6 14.6 31.3 (max-min)/avg 18% 0% 0% 3% 9% 5% 4% 7% 12% 7% 14% 150 Appendix B Average mill records of digester operation at HSPP * In all the cases below, efforts have been made to keep the operations stable for a minimal of 6 hours. A A l B C D E F G H I J K L Tagname Deviation 411AI9129 411AT9606 411AYI9125 A 411AYI9125 B 411AYI9125 C 411AYI9125 D 411AYI9125 E 411AYI9125 F 411CT9036 411DRPM 411FI9121 Engineering Units % KAPPA g/L H-Fctr H-Fctr H-Fctr H-Fctr H-Fctr H-Fctr G/L RPM L/S Description (max-min)/ avg Dig Blowline Kappa BC RESIDUAL ALKALI H-Factor to Wash H-Factor to Mid EMCC H-Factor to End MCC H-Factor to Mid MCC H-Factor to Extraction H-Factor in IV WL STRENGTH CHIP METER RPM DIG BLOW FLOW(F12A IB) 09-Mar-04 23 29.6 12.8 1518.2 1257.4 874.2 682.4 380.2 0.0 84.6 15.1 131.9 11-Mar-04 18 28.0 12.6 1412.0 1165.0 808.0 640.3 367.2 0.0 87.3 15.0 133.5 26-Mar-04 27 28.4 13.5 1517.6 1338.9 975.4 767.1 413.0 0.0 85.9 10.4 88.2 12-Apr-04 16 28.0 14.6 1438.6 1179.5 809.4 633.5 349.8 0.0 84.0 15.7 135.0 06-May-04 16 31.0 13.5 1344.1 1123.9 799.3 640.1 369.1 0.0 89.7 15.4 125.7 02-Jun-04 225 29.3 14.7 1397.5 1130.8 747.9 573.0 302.8 0.0 84.0 15.0 119.0 05-Jun-04 63 29.2 13.8 1293.1 1106.8 807.0 653.6 377.7 0.0 86.7 14.5 124.3 15-Jun-04 24 28.3 13.5 1437.9 1185.8 824.4 651.6 358.5 0.0 87.0 14.9 139.6 23-Jul-04 19 35.7 12.9 1526.9 1251.5 853.5 665.1 367.5 0.0 90.4 13.8 112.9 07-Aug-04 24 40.9 12.8 1596.3 1291.4 868.7 675.4 369.1 0.0 87.6 12.9 108.9 08-Aug-04 41 29.4 13.6 1632.9 1336.6 912.1 712.0 402.4 0.0 90.0 11.1 92.9 09-Aug-04 36 32.0 13.5 1636.6 1317.4 869.0 672.8 362.8 0.0 90.7 12.0 104.8 11-Aug-04 24 32.3 13.5 1444.4 1171.3 783.6 606.3 329.2 0.0 92.1 14.8 140.0 08-Sep-04 28 32.3 13.3 1217.0 1023.4 734.4 595.1 353.8 0.0 89.6 13.3 108.6 13-Sep-04 19 35.1 13.5 1198.6 1012.3 729.2 592.2 349.1 0.0 85.7 15.4 121.5 151 Appendix B Average mill records of digester operation at HSPP (continued) A M N O P Q R S T U V W Tagname 411FT9011 411FT9021 411FT9031 411FT9032 411FT9037 411FT9038 411FT9039 411FT9041 411FT9111 411FT9141 411FT9160 Engineering Units T/H T/H L/S L/S L/S L/S L/S L/S L/S L/S L/S Description HI PRESS STEAM FLOW LP STEAM FLOW F3A HPF WL PURGE F 3EWL FLOW TO WASH F3BWL FLOW TO BC F 3D WL FLOW TO MC F3FWL FLOW TO IV F4 CB FLOW F11 CB FOR TO IV DIG PRES F14CBTO DIG OD F16AEXTR SCRN FLOW 09-Mar-04 35.0 27.3 3.9 4.7 1.7 4.7 46.0 0.4 47.8 10.9 139.5 11-Mar-04 37.3 26.4 3.5 4.4 0.8 4.5 45.1 0.3 52.0 5.9 141.7 26-Mar-04 24.6 18.0 4.0 3.3 1.9 3.3 30.9 0.3 31.9 3.7 109.2 12-Apr-04 38.1 22.7 4.0 4.4 0.2 4.4 48.3 5.0 52.2 11.4 147.4 06-May-04 32.3 20.7 3.9 4.7 0.8 4.6 43.3 0.4 47.3 15.5 133.6 02-Jun-04 27.3 26.2 4.1 3.9 0.0 3.9 48.7 1.2 30.8 7.3 109.4 05-Jun-04 28.2 33.1 4.0 3.8 2.0 3.7 42.3 4.0 36.8 8.2 110.5 15-Jun-04 29.0 21.4 3.7 4.0 0.5 4.0 44.0 0.5 44.6 10.0 115.9 23-Jul-04 29.4 15.9 3.6 3.9 0.8 3.5 41.2 8.0 45.7 10.7 120.2 07-Aug-04 27.4 16.0 4.1 4.6 1.5 4.4 38.7 6.0 42.2 10.3 115.1 08-Aug-04 26.2 13.1 4.3 3.9 2.1 3.9 32.4 6.0 42.1 10.2 116.8 09-Aug-04 26.1 14.4 4.3 4.3 2.9 4.3 34.9 6.0 40.3 9.7 107.7 11-Aug-04 29.1 17.4 4.3 3.8 0.4 4.2 43.4 8.0 47.8 7.0 120.6 08-Sep-04 25.7 20.4 4.2 3.4 0.4 3.5 38.3 9.0 39.6 11.3 110.4 13-Sep-04 30.6 22.8 3.9 4.5 0.9 4.5 44.6 10.0 40.6 9.7 120.2 152 Appendix B Average mill records of digester operation at HSPP (continued) A X Y Z AA AB AC AD AE AF AG AH Tagname 411FT9181 411FYI9025 411FYI9121 411HI9022B 411HT9022B 411LY9102 411LY9752 411MI9022 411NIC9183 411TT9085 411TT9165 Engineering Units L/S T/H ADT/D FURNSH N/A % % % % BD DEG C DEG C Description F18 CBTO DIG BOT CALC STEAM FLOW TO BIN DIGESTER BLOW RATE DIGESTER EXIT FURNISH PRODUCTIO N GRADE DIGESTER CHIP LEVEL IV CHIP LEVEL WEIGHTED AVE MOIST TRIM LQR PD18 BLOW TEMPERATU EXTRACTIO CONS MV RE N TEMP 09-Mar-04 39.5 28.8 1267.6 400.0 400.0 28.1 17.9 54.3 10.0 157.9 160.8 11-Mar-04 40.8 28.1 1312.9 400.0 400.0 36.3 26.7 54.5 10.2 157.3 159.9 26-Mar-04 29.7 20.9 888.7 400.0 400.0 38.7 29.9 54.5 10.5 153.7 158.3 12-Apr-04 48.3 24.4 1316.6 400.0 400.0 29.7 15.0 52.5 10.2 157.8 160.6 06-May-04 43.2 22.8 1236.6 400.0 400.0 30.4 15.5 50.5 10.2 158.1 160.7 02-Jun-04 29.0 23.9 1223.1 400.0 400.0 28.4 0.7 49.4 10.7 155.6 158.5 05-Jun-04 30.2 22.1 1265.1 400.0 400.0 34.0 10.6 49.5 10.6 157.2 160.7 15-Jun-04 34.6 23.0 1409.7 400.0 400.0 33.0 33.0 49.5 10.5 157.4 160.4 23-Jul-04 37.6 18.9 1127.8 400.0 400.0 29.9 16.9 46.3 10.4 156.6 159.5 07-Aug-04 37.9 18.6 1100.2 400.0 400.0 28.1 16.4 47.2 10.5 155.6 158.7 08-Aug-04 37.7 15.4 937.1 400.0 400.0 32.1 22.8 46.8 10.5 154.8 158.3 09-Aug-04 35.8 16.6 1057.6 400.0 400.0 32.5 27.9 47.1 10.5 154.4 157.5 11-Aug-04 35.6 19.2 1397.7 400.0 400.0 28.3 19.8 47.1 10.4 156.0 159.4 08-Sep-04 32.1 20.1 1083.9 400.0 400.0 30.2 15.7 46.6 10.4 156.4 157.9 13-Sep-04 35.7 24.1 1189.6 400.0 400.0 30.4 12.8 46.5 10.2 158.2 159.7 153 Appendix B Average mill records of digester operation at HSPP (continued) A A l A J AK A L AM AN AO AP AQ AR AS Tagname 411TT9195 411TT9197 411TT9205 411TT9207 411TT9600 411TT9605 411TT9607 411TT9700 411TT9705 411TT9710 411TT9711 Engineering Units DEG C DEG C DEG C DEG C DEG C DEG C DEG C DEG C DEG C DEG C DEG C Description T19C MC COLD TEMP T19H MC HOT TEMP T20C WASH COLD TEMP T20H WASH HOT TEMP T60C BC COLD TEMP T60H BC HOT TEMP T60L BC H TEMP SP IV UPPER TEMP IV LOWER TEMP TRIM ZONE TEMP EXTR ZONE 1 TEMP 09-Mar-04 161.6 160.9 153.2 159.1 158.4 175.9 158.0 124.0 129.9 157.6 160.4 11-Mar-04 160.7 160.1 146.1 158.1 157.9 176.2 157.5 123.6 129.2 157.0 159.5 26-Mar-04 159.9 158.1 146.7 154.0 154.2 166.4 153.6 116.8 127.2 153.7 159.0 12-Apr-04 161.4 160.9 143.8 159.1 157.8 170.4 157.3 122.5 129.1 157.2 160.2 06-May-04 159.8 159.4 145.9 157.1 158.2 170.0 157.5 126.2 131.1 157.9 158.6 02-Jun-04 161.9 160.6 157.1 159.0 155.8 174.5 . 155.4 122.3 128.0 155.0 161.7 05-Jun-04 159.5 157.6 149.9 153.9 157.6 171.7 157.0 122.6 128.5 157.4 159.1 15-Jun-04 161.3 159.9 147.4 157.5 157.4 169.0 156.9 124.8 131.6 157.3 161.0 23-Jul-04 160.5 160.5 145.9 159.1 156.6 168.0 156.1 125.8 133.3 156.4 159.2 07-Aug-04 160.4 160.1 144.9 159.0 155.9 164.5 155.4 125.5 133.0 155.3 159.3 08-Aug-04 159.1 158.6 137.8 157.0 154.4 159.8 153.5 126.0 133.5 154.4 157.1 09-Aug-04 159.6 159.5 144.7 158.2 154.9 161.9 154.4 126.8 133.6 153.8 158.6 11-Aug-04 161.3 160.6 145.9 159.0 156.3 167.7 155.7 126.8 132.5 155.8 160.1 08-Sep-04 155.7 156.0 145.9 153.5 156.5 165.5 155.6 125.3 130.5 155.4 154.2 13-Sep-04 158.2 157.9 148.5 155.3 158.0 170.7 157.1 123.9 128.4 157.1 156.9 154 Appendix B Average mill records of digester operation at HSPP (continued) A AT AU AV AW AX AY AZ BA BB BC BD BE Tagname 411TT9712 411TT9713 411TT9714 411TT9715 411U19160 411WI9022 411XI9022A411XI9022B411XI9022C411XI9022D 411FT9601 411FT9611 Engineering Units DEG C DEG C DEG C DEG C % T/DBD % % % % L/S L/S Description EXTR ZONE 2 TEMP EXTR ZONE 3 TEMP MCC ZONE TEMPERAT URE WASH ZONE TEMP CALC DILU FACTOR TOTAL CHIPS TO CB HEMLOCK FURNISH FRACTION CEDAR FURNISH FRACTION SPF FURNISH FRACTION FIR FURNISH FRACTION F60 BC FLOW F61A UPPER SLUICE FL 09-Mar-04 160.7 160.8 162.0 160.9 -0.4 2428.1 0.0 43.4 56.7 0.0 150.4 121.6 11-Mar-04 159.7 159.7 160.1 158.8 0.1 2366.5 0.0 46.8 53.2 0.0 136.5 122.0 26-Mar-04 159.5 159.5 159.4 155.8 0.0 1909.0 0.0 47.2 52.8 0.0 136.3 89.3 12-Apr-04 160.7 160.0 160.0 158.8 0.5 2517.0 0.0 47.1 52.8 0.0 192.2 140.1 06-May-04 159.0 158.6 157.8 159.1 0.7 2409.0 0.0 46.9 52.9 0.0 164.9 155.0 02-Jun-04 161.8 161.8 161.6 159.7 -1.6 2691.1 0.2 45.7 55.1 0.0 130.0 121.9 05-Jun-04 159.3 159.3 158.1 157.4 -1.3 2527.6 0.0 46.3 53.1 0.2 164.8 127.6 15-Jun-04 160.5 161.3 160.6 160.7 -1.0 2717.8 0.0 47.3 53.1 0.0 167.3 140.9 23-Jul-04 159.4 159.3 159.8 159.7 0.6 2053.4 0.0 43.6 56.2 0.0 149.4 140.1 07-Aug-04 159.6 159.3 159.6 159.4 0.7 1822.6 0.0 49.2 51.1 0.0 161.2 136.4 08-Aug-04 157.5 157.4 157.2 153.5 2.3 1629.0 0.0 47.6 52.9 0.0 165.3 113.1 09-Aug-04 158.7 158.6 158.7 158.2 0.6 1736.0 0.0 48.0 51.5 0.0 166.7 127.3 11-Aug-04 160.5 160.1 160.1 159.9 0.4 2086.2 0.0 47.4 52.6 0.0 144.6 151.3 08-Sep-04 154.4 154.5 154.8 153.0 -0.2 2064.6 0.0 44.9 55.2 0.0 193.2 122.6 13-Sep-04 157.2 157.1 156.7 156.0 -0.8 2361.3 0.0 45.0 54.5 0.0 184.9 137.4 155 Appendix B Average mill records of digester operation at HSPP (continued) A B F BG BH B I B J BK BL BM BN BO BP BQ Tagname 411FT9191 411FT9201 411FT9614 411TT9125 411TT9245 411NT9126 411IT0024 411PDT918 3 411FT9871 411AT9158 411AT9190 421AI9129 Engineering Units L/S L/S L/s DEG C DEG C % A KPA L/S g/i g/i KAPPA Description F19AMC UPPER FLOW WASH CIRC FLOW CONT BC Heater Bypass Flow DIF INLET TEMP CB COOLER TEMP DIG BLOW CONSIST OUTLET DEVICE (DIG) PD18 BLOW CONSIST BLOW Extraction LINE DILN Resid Alkali MC Resid Alkali Decker Vat Kappa 09-Mar-04 120.1 95.3 15.5 81.7 55.6 10.9 366.0 82.4 0.0 10.4 14.2 31.6 11-Mar-04 120.1 89.3 44.3 82.0 58.5 11.0 405.3 90.6 0.1 10.7 13.9 28.5 26-Mar-04 119.8 89.7 47.5 82.3 56.7 11.2 408.7 99.0 0.0 10.4 13.3 30.3 12-Apr-04 120.0 96.7 0.1 82.0 61.0 11.0 378.4 87.6 0.0 11.5 13.8 30.1 06-May-04 120.0 95.3 0.0 81.9 63.0 10.8 402.8 90.7 0.0 10.0 15.2 30.2 02-Jun-04 110.1 92.5 41.5 82.0 40.7 10.9 409.4 105.4 9.0 11.3 20.2 29.1 05-Jun-04 110.1 86.3 -0.1 79.7 44.0 10.7 409.4 102.3 8.9 10.3 20.3 29.7 15-Jun-04 110.1 86.7 0.0 74.4 43.0 10.9 409.1 99.0 8.6 10.1 12.7 28.6 23-Jul-04 119.9 88.2 20.9 84.2 64.3 11.3 403.5 95.2 0.0 9.4 12.6 34.1 07-Aug-04 120.2 84.6 0.2 86.3 64.2 11.2 408.8 98.4 0.0 8.6 12.5 39.9 08-Aug-04 120.2 83.3 0.2 85.1 69.7 11.3 408.5 99.0 0.2 9.1 12.6 28.0 09-Aug-04 120.2 83.1 0.0 86.2 63.9 11.3 407.8 98.3 0.3 9.0 12.8 30.9 11-Aug-04 120.1 90.1 21.0 86.0 65.3 10.1 407.4 95.7 10.3 9.7 13.1 31.5 08-Sep-04 120.0 89.4 -0.6 80.6 57.7 11.2 407.5 95.1 0.2 8.7 12.8 32.9 13-Sep-04 120.1 86.2 3.0 82.6 52.6 10.9 392.8 89.5 0.0 8.5 13.3 35.6 156
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerical investigation of industrial continuous digesters
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerical investigation of industrial continuous digesters Fan, Yaoguo 2005
pdf
Page Metadata
Item Metadata
Title | Numerical investigation of industrial continuous digesters |
Creator |
Fan, Yaoguo |
Date Issued | 2005 |
Description | Digesters, either of batch or continuous type, are large chemical cooking vessels used to convert wood chips into fiber solution or pulp, from which paper products are eventually produced. In the continuous digesters, the complex liquor flow, which usually has a number of circulation loops, interacts with the downward moving wood chips through mass, momentum, and energy exchanges and therefore has a major effect on the chip motion and the cooking reactions. A set of comprehensive differential equations governing the two-phase flow, heat transfer, and chemical reaction is developed, with the purpose of computationally investigating the cooking process in continuous digesters. A new cooking reaction model is also derived, based on extensive experimental data. To model the chip phase motion, the chip phase is assumed to be a non-Newtonian fluid with a slip boundary condition on the digester walls. A number of sub-models are examined and used to account for the compressibility of chip column and the inter-phase friction force. Based on the governing equations and mathematical models, a general two-phase three-dimensional CFD code has been developed to calculate the flow fields for both phases, solid fraction, temperature, and the concentrations of chemicals. The developed CFD code also has the capability of solving many other potential multiphase reactive flow problems. Liquid flow in a laboratory digester filled with compacted chips is simulated by the developed code. The computations compare reasonably well with the experimental observations by ERT (Electrical Resistance Tomography). Currently, the code cannot predict the tracer motion. Therefore detailed comparison between CFD simulation and ERT images cannot be made. Simulations are also performed in a full-scale industrial digester. The predictions of yield and kappa number at the digester exit are accurate for various cooking conditions. The distributions of pressure and chip phase volume fraction are discussed, the chip phase velocity profiles at different digester levels are plotted, and the kappa number variation in the digester exit plane is also addressed. Simulations also show that an increased production rate results in a higher kappa number and a lower solid pressure at the bottom exit, and therefore the digester becomes more difficult or impossible to operate. The impact of variations of chip type and cold blow flow rate on cooking performance is also investigated computationally. The extensive investigation of the parameter variation impact on the process has shown that the model developed constitutes a reliable tool to predict the complex process in industrial digesters. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2009-12-21 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080759 |
URI | http://hdl.handle.net/2429/16946 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2005-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2005-104880.pdf [ 11.7MB ]
- Metadata
- JSON: 831-1.0080759.json
- JSON-LD: 831-1.0080759-ld.json
- RDF/XML (Pretty): 831-1.0080759-rdf.xml
- RDF/JSON: 831-1.0080759-rdf.json
- Turtle: 831-1.0080759-turtle.txt
- N-Triples: 831-1.0080759-rdf-ntriples.txt
- Original Record: 831-1.0080759-source.json
- Full Text
- 831-1.0080759-fulltext.txt
- Citation
- 831-1.0080759.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080759/manifest