C F D S I M U L A T I O N O F M I X I N G D Y N A M I C S IN A G I T A T E D PULP S T O C K CHESTS by CLARA FORD Honours in Applied Science., Universidad Pontificia Bolivariana, Medellin, Colombia, A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF CHEMICAL AND BIOLOGICAL ENGINEERING THE UNIVERSITY OF BRITISH COLUMBIA December 2004 © Clara Ford, 2004 2002 11 ABSTRACT Agitated pulp chests function as low-pass filters to reduce high frequency variability in pulp properties (mass concentration, freeness etc.) ahead of many pulping and papermaking operations. Tests made on both industrial and scale-model chests have shown that their dynamic performance is far from ideal, with a significant extent of non-ideal flow (short circuiting, recirculation and stagnation) possible (Ein Mozzafari et al., 2004). In this work, the flow field of a 1/11 scale-model pulp chest was modeled using commercial CFD software th (Fluent) with the pulp suspension treated as a Bingham plastic. A multiple reference frame approach was used with coupling between reference frames made using a velocity transformation. The flow profiles predicted by the simulation agreed qualitatively with the observed profiles in the experiments. The power input predicted by the simulations was slightly higher ( + 12%) than the measured power. The velocity field obtained from the CFD model was used to obtain the system's dynamic response to a frequency-modulated random binary input signal. This data was then used as input to a dynamic model developed by Ein Mozzafari (2002) that treats flow within the chest as following two streams: one that bypasses the mixing zone and one that enters it. For both streams, the fraction of suspension passing through each zone was determined and a time constant and delay time computed. These parameters were then compared with parameters measured experimentally under identical operating conditions and the parameters associated with the mixing zone found to agree to within ±25%. The CFD simulation provides detailed information on the velocity profile within the chest and allows the location(s) of poor mixing regions to be identified. Hence, design methods and trouble shooting of existing stock chests can be significantly improved. iii T A B L E OF CONTENTS ABSTRACT ii List of Tables iii List of Figures v Acknowledgements ix 1 INTRODUCTION 1 2 LITERATURE REVIEW. 7 2.1 Introduction 7 2.2 Pulp suspension rheology 8 2.3 Mixing scales 14 2.4 Performance and Design Considerations of the Mixing Unit 15 2.5 2.6 2.4.1 Chest Shape and Impeller Types 15 2.4.2 Measures of Mixer effectiveness 18 2.4.3 Power Consumption for Newtonian Fluids 18 2.4.4 Power Consumption for Non- Newtonian Fluids 20 2.4.5 Power Consumption for Pulp Suspension 22 2.4.6 Design Methods for Agitated Pulp Stock Chests 24 Dynamics of Agitated Pulp Stock Chests Computational Fluid Dynamics 26 31 2.6.1 C F D Modeling of Stirred Tanks 33 2.6.2 C F D Modeling of Pulp Suspensions 36 iv 2.6.3 3 4 38 EXPERIMENTAL 39 3.1 Experimental setup 39 3.2 Impeller specifications 40 3.3 Materials 42 3.4 Test procedure 43 3.5 Experimental conditions 43 3.6 Power Measurement 44 CFD M O D E L DEVELOPMENT 46 4.1 Introduction 46 4.2 General Transport Equations and Discretization 46 4.2.1 Solution Convergence 49 4.2.2 Discretization of the Momentum Equation 50 4.2.3 Discretization of the Continuity Equation 51 4.3 Flow in a Rotating Frame 52 4.4 Scale Mixing Chest Model 54 4.5 5 Objectives of this work 4.4.1 Geometry 54 4.4.2 Meshing the Geometrical Model 56 4.4.3 Problem Setup for Multiple Reference Frames in F L U E N T 59 Mixing Dynamics Study with Tracer Analysis 66 CFD SIMMULATION RESULTS 69 5.1 Introduction 69 5.2 Newtonian Fluid Simmulation 69 5.3 Non-Newtonian Simmulation for Pulp Suspension 73 V 5.4 Power Input Comparison 5.5 Numerical Dynamic Test Results 6 84 86 CONCLUSIONS AND RECOMMENDATIONS FOR F U T U R E W O R K 101 6.1 Conclusions 101 6.2 Recommendations for future work 104 NOMENCLATURE 106 BIBLIOGRAPHY 108 APPENDIX A: UPWIND S C H E M E S AND T E M P O R A L DISCRETIZATION 116 A P P E N D I X B: T U R B U L E N C E M O D E L L I N G : R N G k - s M O D E L 120 A P P E N D I X C : E X C I T A T I O N S I G N A L UDF'S 124 A P P E N D I X D: SENSITIVITY A N A L Y S I S D A T A 126 vi List of Tables Table 3.1: Maxflo impeller specifications 42 Table 3.2: Summary of operating conditions studied in the experiments 44 Table 4.1: 3-D coordinates of the points along the blade wire frame 55 Table 5.1: Comparison of the power input calculated from C F D and measured experimentally 85 Table 5.3: Comparison of the parameters estimated from experimental and CFD data for F B K pulp suspension at C = 3.3% at Q = 7.9 L/min and Q = 21.1 L/min m for an impeller speed of 1031 rpm 92 Table 5.4: Comparison of the parameters estimated from experimental and CFD data for F B K pulp suspension at C - 3.3% at Q = 37.1L/min and Q = 7.9. L/min m for different impeller speeds 95 Table 5.5: Comparison of the parameters estimated from experimental and C F D data for F B K pulp suspension at C = 3.3% at Q = 7.9 L/min and different impeller m speeds for Config.2 97 Table D . l : Comparison of the parameters estimated from experimental and C F D data for F B K pulp suspension at C = 3.3% ± 95% confidence interval of z m at Q = 7.9 L/min and 7 V = 1031 rpm y 125 vii Table D.2: Comparison of the parameters estimated from experimental and C F D data for F B K pulp suspension C = 3.3%, at Q = 7.9 L/min and N = 1031 rpm m for different yielding viscosity values 126 Table D.3: Comparison of the parameters estimated from experimental and CFD data for F B K pulp suspension C = 3.3% , at Q = 7.9 L/min, N = 1031 rpm. m and ju =100 kg/m.s, for different flow index values 0 128 Vlll. List of Figures Figure 1.1: Step response of an industrial blend chest 2 Figure 1.2: The frequency response measured for the industrial chest compared with its ideal response 3 Figure 2.1: Shear stress vs. shear rate according to the two viscosity Herschel-Bulkley model 11 Figure 2.2: Torque vs. rotational speed response of pulp in dynamic shear test Gullichsen and Harkonen (1981) 12 Figure 2.3: Flow pattern for an axial-flow impeller: (a) top-entry (b) side-entry 18 Figure 2.4: Dynamic Model developed by Ein-Mozaffari et al. (2003) 30 Figure 3.1: Schematic of the experimental set up 40 Figure 3.2: Scale model chest dimensions 41 Figure 3.3: Maxflo impeller used in the scale model chest ...42 Figure 3.4: Input-Output configurations studied 44 Figure 4.1: Maxflo Mark II (6.5") impeller wire-frame in GAMBIT 54 Figure 4.2: Mixing chest geometry in GAMBIT 56 Figure 4.3: Multiblock division of the computational domain 57 Figure 4.4: Enlarged view of the mesh and size function used around the impeller 58 Figure 4.5: Three dimensional mesh of the scale mixing chest 59 Figure 4.6: Overview of the segregated solution method 62 ix Figure 4.7: The boundary conditions specified in the domain 64 Figure 4.8: U D F code of thefrequency-modulatedbinary input signal (for Q =37.1 L/min).. 68 Figure 5.1: Velocity vectors colored by velocity magnitude (m/s) for water wth Q =30 L/min and TV =1250 rpm 70 Figure 5.2: Velocity vectors colored by velocity magnitude (m/s) on the plane x= 0 cm, for water with Q =30 L/min and TV =1250 rpm 71 Figure 5.3: Velocity vectors colored by velocity magnitude (m/s) on the impeller wall z = 0 cm, for water with Q =30 L/min and TV =1250 rpm 72 Figure 5.4: Velocity vectors colored by velocity magnitude (m/s) on the wall opposite to the impeller, z = 53 cm, for water with Q = 30 L/min and TV=1250 rpm .73 Figure 5.5: Velocity vectors colored by velocity magnitude (m/s) on the side wall x= 20 cm, for water with Q = 30 L/min and TV =1250 rpm 74 Figure 5.6: Velocity vectors colored by velocity magnitude (m/s) on the top surface for water with Q = 30 L/min and TV =1250 rpm Figure 5.7: Stress vs. rotor angular displacement for a 1.6% C 74 m semi bleached kraft Pulp tested with the Haake RV12. Rotor MVIP, Housing RVH1 1988) (Bennington, 1988) 77 Figure 5.8: Converge history of the scaled residuals for a typical step by step simulation 78 Figure 5.9: Velocity vectors colored by velocity magnitude (m/s) on a surface aligned with the blade tip at x = 9 cm. Pulp suspension at C =3.3%, m X Q = 37.1 L/min, 7Y= 1203 rpm 79 Figure 5.10: Velocity vectors colored by velocity magnitude (m/s) on the side wall surface. Pulp suspension at C =3.3%, Q = 37.1 L/min, N= 1203 rpm m 80 Figure 5.11: Velocity vectors colored by velocity magnitude (m/s) on the top surface. Pulp suspension at C =3.3%, Q = 37.1 L/min, N = 1203 rpm m 81 Figure 5.12: Velocity vectors colored by velocity magnitude (m/s) on the impeller wall surface. Pulp suspension at C =3.3%, Q = 37.1 L/min, N- 1203 rpm m 82 Figure 5.13: Velocity vectors colored by velocity magnitude (m/s) on the opposite impeller wall surface. Pulp suspension at C =3.3%, m Q = 37.1 L/min, N= 1203 rpm 83 Figure 5.14: Grid independence verified by comparing the velocity in front of the impeller for 3 different levels of adaption 84 Figure 5.15: Calculated and measured power input vs. impeller rotational speed C =3.3%; Q = 37.1 L/min 86 m Figure 5.16: Trajectories of massless tracer particles released from the velocity inlet 87 Figure 5.17: Single step input signal and output response obtained from the C F D simulation. C =3.3%, Q = 7.9 L/min and N = 1456 rpm m 88 Figure 5.18: Input/Output signals from the numerical dynamic test and the scale model chest for different operating conditions 91 Figure 5.19: Velocity vectors colored by velocity magnitude (m/s) on the plane x = 0 cm. Pulp suspension at C =3.3%, Q = 37.1 L/min, N = 1203 rpm m Config.2 Figure 5.20: Trajectories of massless tracer particles released from the velocity inlet 98 xi for Config.2 99 Figure 5.21: Input/Output signals from the numerical dynamic test and the scale model chest for Config.2. operating at Q = 7.9 L/min, N = 1203 rpm 99 Figure D . l : Input/Output signals from the numerical dynamic test and the scale model chest for C =3.3%, Q =7.9 L/min, N = 1031 rpm. Comparison m of the CFD calculated response for different yielding viscosity values 127 Figure D.2: Input/Output signals from the numerical dynamic test and the scale model chest for C =3.3%, Q = 37.1 L/min, N = 1456 rpm. Comparison m of the CFD calculated response for different yielding viscosity values 127 Xll Acknowledgements I would first like to express my sincere gratitude and appreciation to my supervisor Dr. Chad PJ. Bennington for his guidance and encouraging enthusiasm during all times of this work. I also would like to give special thanks to Dr. Fariborz Taghipour for always having time to answer my questions and giving a great input to this work. I gratefully acknowledge Dr. Farhad Ein-Mozaffari for his innumerable assistance and for all the interesting and fruitful discussions on this topic. I acknowledge the assistance of all the staff in the Pulp and Paper center and the Chemical and Biological Engineering Department at UBC. Financial support from the University of British Columbia (UGF Scholarship) is gratefully acknowledged. To my parents and Matthew, For their love and support. Chapter 1: Introduction 1 INTRODUCTION Agitated pulp stock chests are used for a number of purposes in the pulp and paper manufacturing process. They act as storage vessels, allowing continuous operation of processes in the event of a breakdown or upset elsewhere in the system. They can also act as blend chests to mix different pulp stocks and /or pulp suspension with process chemicals. They keep the pulp in motion and prevent the suspension from dewatering. Their chief purpose though, is to reduce high-frequency disturbances in pulp properties (mixture composition, fibre mass concentration, freeness, etc.) ahead of many pulping and papermaking unit operations. In this capacity, stock chests complement the action of control loops that are only able to attenuate low-frequency variability below their cut off frequency (Bialkowski, 1992). In spite of the importance of their applications, the problem of designing and scaling up agitated stock chests has been tackled mainly by means of semi-empirical methods. One common design method has been summarized by Yackel (1990) and is based on matching the momentum flux generated by an impeller with that needed to provide complete motion in the vessel. Yackel defined complete motion as occurring when the suspension was in motion over the entire surface of the vessel. Figure 1.1 shows the measured response of an industrial stock chest to a step input change at t = 200 s (Ein-Mozaffari et al., 2003). The upper curve is the input signal and the lower curve is the output response. The chest response is clearly not ideal. The output shows an initial bump at about t = 300 s, which is due to shortcircuiting within the chest and was not stabilized until t - 450 s. The response also shows a series of first-order exponential responses, each one 2 Chapter 1: Introduction getting progressively smaller. From this it can be concluded that a portion o f the stock is recirculated within the chest (Ein-Mozzafari et al., 2003). Blend Chest Mixing Test - Residence Time=14 minutesmin 4.6 4.4 Steady State in 1400 sec, Appr. Time Constants50 sec=5.8 min i r — " 1 1 1 1 — B l a id Che Bt hiletifTTW) Consistency B end CI lest Outlet Consister icy 4.2 4.0 3.8 1 3.6 ; ; ; ; 3.4 200 400 600 800 1000 1200 1400 1600 1800 200(] Time (sec) Figure 1.1: Step response of an industrial blend chest. The output signal has been shifted by -0.2 so as not to overlap with the input signal. The frequency responses calculated for this industrial chest and for a perfectly mixed chest having the same volume and flow rate are plotted i n Fig. 1.2. The degree o f upset attenuation for the industrial chest is considerably worse than that o f the perfectly mixed chest, especially for frequencies from 0.01 to 0.1 rad/s. These frequencies are higher than the cut-off frequencies o f paper machine control loops, which are usually between 0.002 and 0.005 rad/s (Bialkowski,1992). Therefore, disturbances at these frequencies would not be attenuated and w i l l pass through the process. Consequently paper quality and machine runnability are reduced. 3 Chapter 1: Introduction Frequency (rad/s) Figure 1.2: The frequency response measured for the industrial stock chest compared with its ideal response Ein-Mozzafari et al., (2003) developed a dynamic model that described the observed response of an industrial chest. In this model the pulp feed to the chest splits into two streams: one that enters a well mixed zone with the possibility of suspension recirculation in it and one that bypasses this mixed zone. Two first-order transfer functions with time delays are used to describe the mixing occurring in both the agitated and the bypassing zones. The time constant and time delay for the bypassing zone are always smaller relative to those of the mixed zone because the amount of mixing in the bypass flow is generally very limited when compared to the agitated zone. In order to study dynamic behaviour of stock chests under realistic (less than ideal) mixing conditions, Ein-Mozzafari (2002) constructed a flow loop with a laboratory chest that is a geometric-scale down (1/11 scale) of the industrial chest shown in Fig 1.1. The th system was excited and the dynamic model parameters were calculated from input-output data. Chapter 1: Introduction . 4 From the dynamic tests made on the scale-model chest, Ein-Mozzafari et al. (2003) found that the extent of non-desirable flows for good mixing such as channeling, recirculation and stagnation was significant even when complete surface motion was achieved and therefore the power input calculated from existing design criteria did not completely eliminate poorly mixed zones or bypassing in the chest. A significant improvement in the design of agitated stock chests may be obtained by developing models which take into account the actual flow field inside the vessel. Computational Fluid Dynamics (CFD) makes the development of these models possible by mathematically solving the momentum and mass conservation equations that govern the flow in the agitated chest. A great effort has been made during the last years towards the improvement of different factors related with the development of a reliable C F D model for mixing, including computer speed and memory, grid refinement techniques and new methods to accurately model rotating flows. In particular, the development of sliding mesh and multiple reference frame methods now make it possible to perform numerical simulations that include the impeller region in the computational domain instead of imposing boundary conditions at the impeller location (black box method) which at best can only be extended to geometries close to that for which the experimental data were obtained (Brucato et al., 1998). In this work, commercial C F D software (FLUENT 6.1, Fluent Inc, Lebanon, NH) was used to simulate the flow field inside the lab-scale model agitated chest designed by Ein-Mozzafari (2002). A steady-state approximation method called "Moving Reference Frame" (MRF) was used to solve governing equations. In order to implement this method, the vessel volume was divided in two regions, the inner region rotated with the impeller while the outer region was 5 Chapter 1: Introduction stationary. The conservation of mass and momentum equations were solved in each reference frame, with interface coupling performed by velocity transformation. The pulp suspension rheology was approximated with a Herschel-Bulkley model, and the equations of motion were solved in the laminar regime. The flow field obtained agreed qualitatively with the patterns observed in the experiments. The power input calculated from the CFD simulations was compared to the measured power and found to agree to within ± 12%. The dynamic experiments performed by Ein-Mozzafari model using tracer distribution analysis. (2002) were simulated in the C F D The velocity information previously acquired was used to obtain the dynamic response of the system to a simulated frequency modulated random binary input signal designed by Ein-Mozzafari (2002). Input-output data were collected from the simulations for operating conditions identical to the ones studied experimentally by EinMozzafari (2002). The dynamic data obtained numerically and experimentally were directly compared and we concluded that the C F D predicted response agreed fairly well with the experimental response but specific details of the measured signal were not accurately captured by the simulation. This technique provides a valuable tool for validation of the CFD model. A better understanding of the dynamic behaviour predicted by the C F D simulation was obtained by estimating the parameters of Ein-Mozzafari's dynamic model (Ein-Mozzafari, 2002) from the simulated input-output data. These parameters were compared to the ones experimentally determined under identical conditions and in particular, the percentage of channeling and the time constant and time delay associated with the agitated zone found to 6 Chapter 1: Introduction agree to within ±25%. The development of a more descriptive model for the pulp rheology may improve the accuracy of the CFD solution. The C F D simulations provide extensive information of the flow occurring inside the chest allowing the location of poor mixing regions. Thus, they can be effectively used to aid the design and trouble shooting pulp mixing chests. C O N T R I B U T I O N S O F THIS W O R K The main contributions of this work are: 1. The development of a C F D model using F L U E N T for mixing pulp suspensions in a laboratory scale model chest using a steady state approximation. 2. The acquisition and analysis of dynamic information from the steady-state CFD model using a pseudo random binary input signal. 3. The comparison of model generated and experimentally obtained data from identical mixing chests to validate the CFD model. Chapter 2: Literature Review 7 2. L I T E R A T U R E R E V I E W 2.1. Introduction Mixing is a central feature in many processes in all the chemical industries. In the pulp and paper industry, effective mixing is critical for the success of several operations including blending of different pulp stocks in papermaking, consistency control ahead of a paper machine, steam addition and the contacting of chemicals with pulp in a bleaching stage. Agitated pulp stock chests function as low-pass filters to reduce high frequency variability in pulp properties like mass concentration, freeness etc. ahead of many pulping and paper making operations. This complements the action of control loops which can only control low- frequency variability below the control loop cut-off frequency (Bialkowski, 1992). It is important to ensure that such chests are properly designed in order to achieve the desired degree of upset attenuation. Information on mixing pulp suspensions remains limited due in part to the complex rheology involved which causes non ideal phenomena in the stock chest such as channeling, recirculation and the presence of stagnant regions. These factors are not accounted for in the design criteria typically used for industrial stock chests and thus the performance of these mixing units has been found to be far from ideal (Ein Mozzafari, 2003). Computational fluid dynamics is a powerful design tool that allows us to mathematically solve the flow occurring in a particular system. Compelling advantages of this approach are its low Chapter 2: Literature Review 8 cost, the ability to access detailed and complete information about the relevant variables in the domain of interest, and its versatility for design purposes. 2.2. Pulp Suspension Rheology Rheology is the study of flow and deformation. In principle, rheology includes everything dealing with flow behavior. However, in practice rheology has been restricted to the study of the fundamental relations between force and deformation in materials. The simplest case is when the stress is proportional to the deformation (Hooke's law), which is called the ideal elastic solid. Fluids often obey Newton's law where the stresses are proportional to the rate of deformation instead of the deformation itself. If there is a non-linear dependence between stress and rate of deformation, the fluid is called non-Newtonian. Shear thinning or thickening viscosity and normal stresses in shear are some of the most common phenomenon of non Newtonian fluids. Pulp suspensions, consisting primarily of water and fibre, are non-Newtonian fluids. They are continuous fibre networks that possess structure and strength resulting from the interaction between neighbouring fibres. Suspensions having fibre mass concentration greater than 0.5% exhibit cohesive strength that arises from mechanical forces caused by the bending and hooking of the fibres (Bennington, 1996). As the fibre mass concentration increases, the number of fibre/fibre interactions increases which results in increased network strength. In order to cause motion throughout the suspension, the shear stresses imposed have to be greater than the strength of the pulp network. This network strength is known as the yield stress of the pulp suspension. Chapter 2: Literature Review 9 M a n y authors h a v e r e p o r t e d different m e a s u r e m e n t s o f p u l p n e t w o r k strength, ( K e r e k e s et al, 1983). D e s p i t e the differences a m o n g these results, they a l l s h o w a m a r k e d d e p e n d e n c y o n consistency a n d c a n be represented b y r = aC (2.1) h m where x is the strength o f the p u l p n e t w o r k i n P a s c a l s , C m a n d a a n d b fitted constants. the m a s s c o n c e n t r a t i o n i n percent T h e s e parameters d e p e n d o n m a n y factors, i n c l u d i n g test procedure, the p u l p t y p e a n d degree o f p u l p treatment (beating) p r i o r to testing. V a l u e s o f the constant a b e t w e e n 1.18 a n d 24.5 a n d values o f the e x p o n e n t b b e t w e e n 1.25 a n d 3.02 h a v e b e e n reported i n the literature ( B e n n i n g t o n et al, 1990). T h e initiation o f m o t i o n i n a p u l p s u s p e n s i o n first o c c u r s at the w e a k z o n e s o f the suspension, that is, b e t w e e n floes. A s the a p p l i e d shear is increase, the relative m o t i o n b e t w e e n increases a n d a r e d u c t i o n i n floe size o c c u r s ( K e r e k e s , 1983). floes A further increase i n shear causes m o r e intense m o v e m e n t a n d leads to turbulent f l o w . A c c o r d i n g to the d e s c r i p t i o n a b o v e , p u l p suspensions m i g h t be treated as plastic materials. T h e y s h o w little d e f o r m a t i o n u p to certain levels o f stress a n d a b o v e this y i e l d stress the material f l o w s r e a d i l y . O n e r o u g h a p p r o x i m a t i o n o f the r h e o l o g i c a l b e h a v i o u r o f p u l p s u s p e n s i o n m a y be o b t a i n e d u s i n g the B i n g h a m m o d e l w h i c h is the s i m p l e s t m o d e l for a fluid possessing a y i e l d stress: H o o k e a n b e h a v i o u r at stresses b e l o w y i e l d a n d N e w t o n i a n b e h a v i o u r above. ( M a c o s k o , 1994). F o r one d i m e n s i o n a l d e f o r m a t i o n s r = Gy for T < x (2.2) x = x + rjy for x > x y y Chapter 2: Literature Review 10 Where G is the elastic modulus (solid-like behaviour), n is the non-Newtonian viscosity, r is the stress, y is the strain and y is the strain rate. Measurements in pulp suspensions indicate that they behave according to the Bingham model when the shear rate is low (Gullichsen and Harkonen, 1981). Besides the yield stress, fibre suspensions also exhibit non-linear viscoplastic properties such as shear thinning (a decrease in the plastic viscosity as the shear rate increases) (Wikstrom, 2002). These two features are incorporated in the Herschel and Bulkley model (Herschel and Bulkley, 1926), which is a generalization of the Bingham model that takes into account changes in the non newtonian viscosity n (Eqn. 2.2) with the applied shear rate through a power law behaviour: n = Ky"~ . x In numerical modeling, an inherent difficulty of the Bingham and Herschel-Bulkley model is that they are discontinuous because at vanishing shear rates, the non-Newtonian viscosity becomes unbounded. This introduces significant difficulties in the simulations of complicated problems that are only amenable to numerical analysis. However, studies made on viscoplastic fluids, covering a very wide shear rate range with different types of suspensions indicate that there is a lower Newtonian rather than a Hookean regime. (Macosko, 1994). This suggests that a two viscosity model could be used to overcome the discontinuity issues as follows: For low strain rates, the "rigid" material behaves like a very viscous fluid with viscosity fj. . As the g strain rate increases and the yield stress threshold, r , is passed, the fluid behaviour is described y by a power law: Chapter 2: Literature Review 11 n T +k Y ~ (r n y {/Mo) (2.3) where k is the consistency factor, and n is the power-law index. The shear thinning behaviour is obtained when the power index is set to less than unity. Figure 2.1 shows how shear stress ( r ) varies with shear rate (y) for the two viscosity Herschel-Bulkley model. n>1 y i o y Figure 2.1: Shear stress vs. shear rate according to the two viscosity Herschel-Bulkley model Gullichsen and Harkonen (1981) studied the conditions under which pulp suspensions became turbulent in a concentric cylinder apparatus. Figure 2.2 shows torque vs. rotational speed for the bleached pine kraft pulp tested in their work. This figure is a clear illustration of the behavior of pulp suspensions under shear. The yield stress must be exceeded before motion can be initiated, and this yield stress increases as the fibre mass concentration increases. The points on the curves mark the instant at which the vessel contents were observed to come to complete Chapter 2: Literature Review turbulence. 12 The pulp suspension flow curves exhibited a discontinuity as they neared the flow curve for water. Once motion had been initiated, the shear stress increased or remained approximately constant as the rotational speed increased. The creation of turbulence within a suspension is referred to as fluidization. This motion leads to energy dissipation, thus power and energy expenditure have been used to quantify fluidization. REVOLUTIONS, rpm Figure 2.2: Torque vs. rotational speed response of pulp in dynamic shear test. Gullichsen and Harkonen(1981) Wahren (1980) obtained the following correlation for the estimation of power dissipation per unit volume for the onset of fluidization (s ): f Chapter 2: Literature Review 13 s =-JL = 1 . 2 x 1 0 C * 4 (2.4) 3 M where r is the p u l p s u s p e n s i o n y i e l d stress, C m percent a n d ju is the v i s c o s i t y o f water. is the fibre m a s s c o n c e n t r a t i o n g i v e n as a G u l l i c h s e n and H a r k o n e n (1981) defined fluidization as the onset o f turbulent m o t i o n throughout the vessel. T h i s p o i n t w a s c h a r a c t e r i z e d b y a sharp increase i n the torque o n the rotor (figure 2.2). B e n n i n g t o n ( 1 9 8 8 ) o b t a i n e d the f o l l o w i n g e x p r e s s i o n for e from G u l l i c h s e n a n d H a r k o n e n ' s data: f £ =3.4xl0 C^ 2 (2.5) 4 / T h e noticeable d i s c r e p a n c y b e t w e e n the t w o correlations c a n be d u e t o d i f f e r i n g definitions o f f l u i d i z a t i o n ( B e n n i n g t o n et al, 1991). Another likely reason was given b y Bennington and K e r e k e s (1996). T h e y p o i n t e d out that W a h r e n ' s use o f the v i s c o s i t y o f water to characterize the v i s c o s i t y o f a p u l p s u s p e n s i o n i n E q n . (2.4) w a s inappropriate. " O n c e a p u l p s u s p e n s i o n is ' f l u i d i z e d ' , it c a n be a s s i g n e d a n apparent viscosity, w h i c h is m o s t certainly greater than that o f water". B e n n i n g t o n a n d K e r e k e s ( 1 9 9 6 ) p e r f o r m e d tests u s i n g a c o n c e n t r i c c y l i n d e r d e v i c e a n d d e v e l o p e d a n e x p r e s s i o n to estimate the e n e r g y d i s s i p a t i o n r e q u i r e d for f l u i d i z a t i o n i n it: - r 2 . 5 / r> s =4.5x\0 C\ 4 '(D /Dy ;/ 4 f 2 /5 m 23 l.3<D /D<3A;l.0<C T where D T r»\-2.3. T m (2.6) < 12.6 is the d i a m e t e r o f the h o u s i n g , D is the d i a m e t e r o f the rotor, a n d C m is the fibre m a s s c o n c e n t r a t i o n , i n percent. E q u a t i o n (2.6) c a n b e extrapolated to a z e r o g a p size (D -» D ), to o b t a i n T Chapter 2: Literature Review 14 £ =4.5x\0 C J 4 2 f (2.7) Since the power dissipation depends on geometrical parameters of the mixer, the values given by Eqns. (2.7) and (2.5) are different. As a result of this rheological behaviour, fibre suspensions are extremely difficult to agitate. To provide motion through the whole tank the shear stress has to exceed the yield stress every where in the suspension. Since the stress is not uniform over the whole volume, there will be regions in the fluid close to the impeller where the fibre network is completely disrupted and the flow can attain turbulence, while other regions will be laminar or even stagnant in other locations in the tank. 2.3 Mixing Scales Mixing is achieved by the movement of material between various parts of a bulk mass to reduce the non-uniformities or gradients of the property of interest to an acceptable level. Non uniformities scales may span a wide spectrum of sizes. A useful way to define scales of non uniformities is the distance of relative motion required to even out the variations (Bennington, 1996). Macroscale variations may be defined as those requiring substantial backflow, i.e. bulk movement of fluid over large distances. Respectively, macroscale mixing refers to large scale flow characteristics (e.g. convection, turbulent dispersion) that are responsible for large-scale distributions in the system (e.g. residence time distribution, age distribution, etc.) features are characterized by flow pattern, power input and circulation time. These Chapter 2: Literature Review 15 Microscale mixing is concerned with the features that cause attainment of homogeneity at the small scale and molecular level. Microscale variations can be eliminated by shear and by diffusion. As fibre suspensions are discontinuous over small distances it is useful to define two small scales of mixing: the fibre scale to represent non-uniformities on the fibre and floe level, and the microscale to represent non-uniformities approaching the molecular level (Bennington, 1988). 2.4 Performance and Design considerations of the mixing unit The mixing process essentially consists of a combination of the mechanical movement of an agitator and the resultant induced deformation and flow. The mechanisms of mixing can be summarized as diffusion, convection and bulk movement. A mixing device tries to incorporate all these mechanisms together in the most economical manner. Leaving aside capital expenditure, this means achieving the desired result in the shortest time using the least power. This, in turn, means achieving the right balance between diffusion, convection and bulk movement, keeping in mind the properties of the materials to be mixed and the required final product (Sweeney, 1978). 2.4.1. Chest Shape and impeller types Stock chests are often designed based on proprietary information that has been refined with experience. Simple vortexing, the primary flow generated by a solid rotor in a fluid medium, is not sufficient to properly agitate pulp suspension, as each fluid particle follows a small closed loop and because of the yield stress of the suspension, a large stagnation zone is expected away from the rotor (Silvester, 1985). Secondary motion in the axial direction can be generated by pitching the impeller blades or driving the impeller at very high speed. The greatest percentage Chapter 2: Literature Review 16 of all pulp agitation devices employ some type of axial flow impeller and use the side-insert configuration (Yackel, 1990). Most of these impeller designs are proprietary, with each manufacturer offering unique designs. Most of these impellers have three or four blades that are bolted or welded to a central hub that mounts on the shaft. Axial flow impellers (marine type impeller, pitched blade turbine, Maxflo, etc) produce high flow with low turbulence. They also produce more flow per unit power than radial impellers and are more cost effective in flow controlled operations such as solid suspension and blending (Dickey, 2001). Figure 2.3: Flow pattern for an axial-flow impeller: (a) top-entry (b) side-entry The principal locus of flow occurs along the axis of the impeller. The fluid is drawn from one side and discharged on the other side. Entrainment occurs after discharging and even though the velocity decreases, the flow is increased, resulting in a high volume turnover. The flow pattern of an axial flow impeller is illustrated in figure 2.3. Chapter 2: Literature Review 17 Pulp and paper mills use both cylindrical and rectangular chests. Yackel (1990) suggested that to achieve complete motion of a pulp suspension with minimum power input in a cylindrical chest, the ratio of stock height (Z) to the chest diameter (7) should be 0.8 (Z/T = 0.8). Even though cylindrical chests are wasteful of space, they are less expensive to construct and the hoop stress in cylindrical vessels allow thinner walls. However, if several rectangular chests are nested together that factor becomes negligible, space is saved and the cost is the same as that of equivalent cylindrical chests (Yackel, 1990). Yackel (1990) also found that the optimum shape for a rectangular pulp stock agitated chest in order to achieve complete motion was a cube defined by L/W =1.0 and Z/W = 1.0, where L and W are chest length and chest width, respectively. He pointed out that a cubical chest requires less power input than a cylindrical chest (Z/T = 0.8) with equivalent volume. Since all chests can not be perfect cubes, general rules of thumb exist for their design one of which is that the L/W ratio should not exceed 1.5. Chests with L/W ratio greater than 1.5 to 2.0 will require the use of two agitators. The two units would be sized as though each were in a separate chest, with L/2 used as the effective L for each mixing zone, with each mixing zone designed for an L/W= 1.0 and Z/W = 1.0 (Reed, 1995). Oldshue (1971) compared the performance of side-entering agitators and top entering agitators. He showed that side entering agitators in a chest with Z/T greater than 0.4, are more economical than top-entering agitators for the same total mixed volume. The performance of top entering agitators in chests with Z/T ratios between 0.4 and 0.25 was very similar to that of side entering agitators. However, this Z/T ratio is not appropriate in terms of required power for complete motion. Chapter 2: Literature Review 18 2.4.2. Measures of M i x e r effectiveness The general measure of mixer performance is usually given by the time required to achieve a desired blend, the power load, and the properties of the final product. Because of lack of quantitative measurements of mixer performance, some highly arbitrary criteria exist, especially for suspensions and pastes. With impeller type mixers, usually used in liquid-liquid mixtures, solid suspensions, gas-inliquid dispersions, etc., the power per unit volume has been frequently used as a measure of mixer effectiveness. In general, the best mixer is the one which does the desired job with the smallest amount of power (Sweeney, 1978). Other parameters used to assess the mixer effectiveness are the mixing time (i.e. the time to achieve homogeneity), and the mixing rate (i.e. the rate at which uniformity is approached) and measurements of the degree of uniformity of the product. 2.4.3 Power Consumption for Newtonian Fluids Early investigators used dimensional analysis to derive an equation for agitator power. They considered the impeller power to be a function of the geometry of the impeller and the tank, the properties of the fluid, the impeller speed and the gravitational force. Buckingham's pi theorem gives the following general dimensionless equation for the relationship of variables (Ruhston et al, 1950, Uhl and Gray, 1966) ' pD ^ ^ DN ^ 2 2 V 8 =0 pD N 5 3 (2.8) Chapter 2: Literature Review 19 where D, T, Z , c, w, p, n, I, p, ju, P, N and g are impeller diameter, thank diameter, liquid depth, clearance of impeller off vessel bottom, blade width, pitch of blades, number of blades, blade length, fluid density, fluid viscosity, power, impeller rotational speed, and gravitational acceleration, respectively. For geometrically similar systems Eqn.(2.8) can be reduced to: pD N^ f 2 pD N 5 3 =0 (2.9) The following dimensionless groups will be introduced at this point: Power Number, Po; is a non dimensional measure of the power drawn by the agitator P_ pD N' 5 Froude Number, Fr, represents the ratio of applied forces to gravitational forces. Fr is only considered if there are significant gravitational effects present. = (2,!) g Reynolds Number, Re; indicates the relative importance of the inertial forces to the viscous forces. Re = ^ (2.12) M Equation (2.9) may be written in the form: Po = K(Re) (Fr) a h (2.13) Chapter 2: Literature Review 20 where K, a and b are equation constants. For a fully baffled tank, where surface waves and vortices are absent, the exponent b generally equals 0 and the power number is only a function of the Reynolds number. The power number decreases as inverse of Reynolds number (Re ") through the laminar flow regime, with the steepness of the descent depending upon K, which is a function of the system geometry. At high Reynolds numbers, greater than about 10 , the flow is turbulent and the 4 power number is essentially constant, i.e. Po = B . In between the laminar and the turbulent region there exists a gradual transition zone in which no simple mathematical relationship exists between the power number and the Reynolds number (Harnby et al., 1985). High flow and low turbulence is produced in agitated pulp stock chests equipped with axial flow impellers, thus they generally operate in the laminar and transition regimes. 2.4.4. Power Consumption for non-Newtonian fluids The power requirements for mixing of non-Newtonian fluids have been studied over a wide range of Reynolds number. Blasinski and Rzyski (1972) studied the subject for turbine, propeller and paddle agitators and quoted different forms of the Reynolds number based on a different correlation for each individual geometry. The Reynolds number (Re) for a non-Newtonian fluid is commonly determined in the laminar flow regime by the method introduced by Metzner and Otto (1957). This method assumes that a representative near-impeller average shear rate (y ) can be calculated as a function solely of avg the impeller rotational speed Chapter 2: Literature Review 21 r avg = K,N (2.14) where N (rotations per second, rps, s") is the impeller rotational speed and K was originally s noted (Metzner and Otto, 1957) to be a constant that depends only on impeller geometry. This shear rate is then used to estimate a near-impeller average or effective viscosity, ju . ejr * constant (2.15) For the disk flat-blade turbine used in their study, a value of 13 was proposed for the proportionality constant K s and their approach was found to be a workable one for pseudoplastic fluids. Calderbank and Moo-Young (1959) studied various sizes and types of impellers in Bingham, pseudoplastic and dilatant fluids, and their results agreed with the conclusion of Metzner and Otto (1957). However, they found K to be 10 instead of 13 for pseudoplastic fluids. They also s indicated that this method was least effective for axial flow impellers operating in the transitional flow regime. The power number in this situation was over-predicted by 20-30%, especially for highly shear-thinning fluids. In the transitional flow regime, the impeller lift force becomes significant, and the viscous forces become less important. This flow regime is distinctly different from the laminar-flow situation (Re < 25), for which the Metzner and Otto method was primarily intended. The limitations of the Metzner and Otto Method have been documented by several authors (Posarac and Watkinson, 2000). Some researchers have proposed alternative methods of estimating the power number for a shear thinning fluid in a mixing tank (Rieger and Novak, Chapter 2: Literature Review 1973). 22 Most shear thinning fluids exhibit a yield stress and/or have zero shear viscosity. Bertrand and Tanguy (1996) modified the Metzner and Otto approach, by including the fluid yield stress in the estimation of the impeller Re. From Bertrand and Tanguy (1996), the magnitude of the Bingham number given in Eqn. (2.16) determines whether the fluid yield stress (T ) would contribute significantly to the power consumption of the impeller. Bi = 4 ~ (2-16) L •ff They found that the effect of the fluid yield stress on K and consequently y s avg through Eqn. (2.14), is negligible for Bi < 7500. 2.4.5. Power Consumption for Pulp Suspension Various researchers (Gibbon and Attwood, 1962; Blasinski and Rzyski 1972) have shown that in the mixing of pulp suspension, the power number is a function of Reynolds number and Reynolds number is a function of ^®/ m where N and D are impeller speed and impeller diameter, respectively, m is a parameter introduced by Head and Durst (1957) and it defines the resistance of the fibre suspension to shear. m=P V P (2.17) where p and r are the denstity and yield stress, respectively. Blasinski and Rzyski (1972) assumed that pulp suspension behaved as a Bingham fluid (Eqn. 2.2). The apparent or effective viscosity, introduced by Metzner and Otto (1957) was then given as: Chapter 2: Literature Review 23 ^.L.hlM. r (2,8) r After applying the Metzner and Otto relationship (y = 107V), we get r+10/Vn Head and Durst (1972) showed that the slope of the flow curve for pulp suspension is small, i.e. the plastic viscosity n p was very small. Therefore, T ^>\0Nrj y p and Eqn. (2.15) may be written as Me,r*— (2-20) 107V ejf Substituting Eqns. (2.19) and (2.16) in to the Reynolds number gives: R e = p D W = 1 0 r / v D V (2.21) m j The approach above was the one used by Gibbon et al. (1962) and Blansinski et al. (1972) in their studies of agitation of pulp suspensions in cylindrical chests equipped with top entering impellers. They used a Bingham rheological model for the pulp suspension and their study confirmed the dependence of the power number on the Reynolds number which was defined as a function of ^ ^ / according to Eqn. (2.17) and the Metzner and Otto's correlation given in m Eqn. (2.14). Chapter 2: Literature Review 24 Attwood and Gibbon (1963) obtained the following correlation for the power required to attain motion across the bottom and turn over the entire contents in the agitated stock chest; defined as complete motion. P = CD-m 2 6 29 0 (2.22) where P is the required power for complete motion in the chest, D is the impeller diameter, and m is defined by Eqn. (2.17). Co is a constant that depends only on the geometry of the mixer. 2.4.6. Design methods for agitated pulp stock chests For any impeller, operating freely in a vessel, the volumetric discharge is at its minimum at the impeller while the fluid velocity is at its maximum at the same point. As the fluid stream moves away from the impeller, the boundary layer between the fast moving fluid and the stagnant fluid in the vessel begins to entrain additional flow. At the same time, the velocity of the stream begins to decay (Yackel, 1990). Based on the fact that as the flow increases, velocity decreases in the same ratio, investigators created the concept of the conservation of momentum and used it to describe the capacity of specific impellers. Fox and Gex (1956) found that the mixing quality is a function of the momentum flux generated by the impeller rather than the power input. The same mixing results can be achieved with different power input by changing the impeller speed and diameter, i.e. a large, low speed impeller can produce equivalent mixing results with less power than is required with smaller, high speed impellers. One of the methods used for designing agitated pulp stock chests is Yackel's method (Yackel and Mahony, 1976; Yackel, 1983; Yackel, 1990 and Reed, 1995) and is based on matching the momentum flux generated by an impeller with that needed to provide complete motion in the 25 Chapter 2: Literature Review chest which was defined as that occurring when smooth surface motion over the entire chest surface was observed. The Impeller momentum flux (Mo) is defined by: (2.23) Mo = Qv where Q is the impeller pumping rate and v is the fluid velocity leaving the impeller. Since Q is proportional to ND and v is proportional to ND, Mo can be rewritten as: 3 (2.24) Mo = CN D 2 4 where C is a constant that depends on the impeller geometry and type (the C factor of the impeller). The momentum flux required to fully agitate a pulp stock chest depends on many factors, many of which have been correlated by Yackel (1990). By determining the required momentum flux for a particular installation using the graphical correlations, a suitable impeller may be selected provided the impeller momentum flux for the impeller is known. Yackel (1990) studied the effects of chest configuration, stock type and consistency, temperature and retention time on the momentum number and presented the following relationship between momentum number and the power required for complete motion in the chest (P): CP 0S Mo = Qy = (2.25) 0S 0A Np N In this work (Yackel, 1990) he also introduced a scale up criterion for agitated pulp stock chests. Level momentum (LMo) is defined as: Chapter 2: Literature Review 26 (2.26) where V refers to the fluid volume. Based on this criterion, level momentum must be equal for the model and the prototype. Later, McDonough (1992) showed that the bulk velocity (v buik ) in an agitated tank is directly proportional to the square root of level momentum developed by the impeller: v bulk oc V Z \o M (221) Equal level momentums at various scales yield equal scaled bulk velocities and equal scaled velocities at corresponding points in vessels with different size. Equal velocities in two different scales will produce similar mixing results for macroscale controlled applications such as agitated pulp stock chest (McDonough, 1992). 2.5. Dynamics of agitated pulp stock chests Pulp stock chests are used for a number of purposes in the pulp and paper manufacture. They act as storage vessels, allowing continuous operation in the event of a breakdown elsewhere in the system. They can also be used as blend chests to mix different pulp stocks and/or pulp suspension with process chemicals. They function as low-pass filters to reduce high frequency variability in mixture composition, fibre mass concentration, freenes and other pulp quality factors ahead of many pulping and papermaking unit operations. In this capacity, stock chests compliment the action of control loops that are only able to attenuate low frequency variability below their cut off frequency (Bialkowski, 1992). Chapter 2: Literature Review 27 Perfect mixing for a Continuous Stirred Tank (CST) describes mixing where the input stream is instantaneously dispersed throughout the vessel volume, with the composition of the outlet stream equal to the uniform composition throughout the tank (Oldshue, 1956). For a perfectly mixed chest, assuming that the flows into and out of the chest are equal, the fibre mass concentration at any point in the mixed volume is the same and the mixed volume equals the volume of the suspension in the chest. Under these conditions the fibre mass concentration at the outlet equals the fibre mass concentration inside the chest. A mass balance for the fibres yields: pQC , -pQC , = v^f± m in m 0UI P (2.28) M r . , dC dt r m,in where t, p, Q, V, and T = ^/Q a r mout m,out time, suspension density, pulp flow rate, suspension volume e in the chest, and chest time constant, respectively. At steady state C m in vv =C m out s s and Equation (2.27) can be written as: \-(c \ where m,in m,in,ss C' ,„ = C ,„ - C m m mi n s s / and \ -C m,out C moul \ — r ^^ '" m m,out,ss =C moul J ~ ) ut ^ -C moutss (2.29) are deviation variables. Taking the Laplace transform of both sides of the Eqn. (2.29) G(s) C' , (s)_ m oul CM m 1 TS + \ (2.30) Chapter 2: Literature Review 28 Thus, the dynamic behaviour of such a perfectly mixed chest is represented by a first order transfer function. The assumption of ideal mixing has often been used to study the impact of pulp mixing vessels in pulping processes. Walker and Cholette (1958) calculated the damping factors for an ideally mixed stock chest having a single feed for various input disturbances, including: continuous sinusoidal, continuous square wave and single square wave inputs. They concluded that the degree of upset attenuation is dependent on the disturbance frequency and chest residence time, but they did not make any comparison of their numerical results with experimental data. Reynolds et al. (1964) studied the degree of upset attenuation for a stock chest as a function of the chest residence time and upset frequency. Their study showed that an additional smoothing of variability over that attained only by agitation could be achieved by externally recirculation part of the output stream or by splitting the feed to the top and bottom of the chest. Brown (1968) also assumed ideal mixing in his study of the dynamic behaviour of a paper mill stock chest. He found that for a perfectly mixed chest the frequency at which attenuation begins is fixed by the residence time of the chest. The complex rheology of a pulp suspension can create considerable deviations from ideal mixing. Since ignoring non-ideal flow can lead to errors in system design (Levenspiel, 1998) it is necessary to study the dynamic behaviour of stock mixing under realistic conditions. Ein-Mozafari et al. (2001) studied the dynamics of agitated pulp stock chests under industrial conditions. They attributed the deviations from ideal mixing to a number of factors including bypassing, recirculation and the presence of dead zones in the chest. The dead zones represent Chapter 2: Literature Review unused volume inside the chest. 29 Such zones can arise from the interaction between the circulation patterns generated by the impeller, the suspension flow through the vessel and the chest geometry (Ein-Mozaffari et al, 2004). The flow within a dead zone need not be completely stagnant - the lack of access to the chest outlet is all that is required. (Silvester, 1985). The dynamic tests made on a pilot scale stock chest by Ein Mozzafari et al, (2003), show that the extent of non-ideal flow can be significant. At high input flow rates and low impeller speeds, the percentage of channeling was found to be as high as 90% of the pulp feed. Even when the scale chest was operated under conditions where the whole suspension in the tank was in motion, the percentage of channeling and dead volume was still 26% and 16%, respectively (Ein Mozzafari 2002). Non ideal flow reduces the degree of upset attenuation produced by the chest. Since some disturbances occur at frequencies higher than the cut-off frequencies of paper machine control loops they are not fully attenuated and pass through to the process. Consequently paper quality and machine runnability are reduced (Ein-Mozzafari et al, 2004). Ein-Mozzafari et al. (2003) developed a dynamic model for an agitated pulp stock chest that included a well-mixed zone, with the possibility of suspension recirculation within the mixing zone and short-circuiting past the mixing zone. The continuous dynamic model is shown in figure (2.5) and was constructed from basic principles, where the system parameters have a physical interpretation. / is the portion of stock that bypasses the mixing zone. A limited amount of mixing can occur in the bypass flow and is given by a first order transfer function (Gi). (1-f) is the fraction of pulp that enters the mixing Chapter 2: Literature Review 30 zone and has a first order transfer function G^. A portion of stock, R, exiting the mixing zone can be recirculated within the mixing zone, i i and 12 are the time constants for the bypassing and agitated zones respectively, and typically T I < 12 . Ti and T2 are the time delays for the bypassing and agitated zone respectively. Again, typically Ti < T2. f y + u - t o - Q 2 _ l-R R Figure 2.5: Dynamic Model developed by Ein-Mozzafari et al. (2003) To estimate the parameters for this model, the system has to be excited. The model parameters are then estimated from the input-output data using a numerical method presented by Kammer et al. (2001). To calculate the effective fully mixed volume of the agitated zone the effect of recirculation is eliminated (Ein-Mozaffari et al, 2002). For R = 0 the time constant for the agitated zone (12) is given by fully mixed 0-/)G (2.31) where Q is the pulp flow rate through the chest and T2 and / are estimated from analysis of the input-output data (this allows calculation of the fully mixed volume using Eqn. (2.31)). Chapter 2: Literature Review 31 From the dynamic tests made on the laboratory scale model chest, Ein-Mozaffari et al. (2002) found that the power calculated based on achieving smooth surface motion, or even the power required for onset of complete motion, did not completely eliminate poorly mixed zones or bypassing in the chest. If the purpose of the chest was to attenuate high frequency variability in fibre mass concentration i.e. to act as low pass filter, additional power, typically four times above that required for complete surface motion, was needed to create an ideal dynamic response in the chest. The same characteristic behaviour was observed in several industrial pulp chests studied later by Ein-Mozaffari et al. (2004). In this work, they conducted a bump test to measure the dynamics of different industrial chests and the input-output data was used to identify the extent of non ideal flow in the chests using the dynamic model shown in figure (2.5). In a number of cases they found that the percentage of the chest involved in active mixing was only between 20% to 40%. From comparison between frequency responses of the real industrial chests and a perfectly mixed chest having the same volume and flow rate, they showed that due to stagnant and poor-mixing zones the degree of variability reduction for the industrial chests was considerably less than that of a perfectly mixed chest. For all the industrial chests studied the fully mixed volumes were less than the theoretical value (100%) despite installations where power input was up to three times greater than that calculated based on creation of complete surface motion (Yackel's criterion). This confirmed the dynamic results obtained earlier from the scale model chest (Ein Mozaffari et al, 2003). 2.6. Computational Fluid Dynamics Computational fluid dynamics (CFD) is the science of predicting fluid flow, heat transfer, mass transfer, chemical reactions, and related phenomena by solving the mathematical equations that govern these processes using a numerical algorithm. The results of C F D analyses are relevant Chapter 2: Literature Review 32 engineering data used in conceptual studies of new designs, detailed product development, troubleshooting and optimization of existing devices (Bakker, 2001). To apply CFD, the geometry of interest is first divided, or discretized, into a number of computational cells. Discretization is the method of approximating the differential equations by a system of algebraic equations for the variables of interest at some set of discrete locations in space and time. The discrete locations are referred to as the grid or the mesh. The continuous information from the exact solution of the Navier Stokes partial differential equations is now replaced with discrete values. The number of cells can vary from a few thousand to millions for large and detailed problems. The computational cells can have a variety of shapes. Triangular and quadrilateral cells are generally used for two dimensional (2-D) problems, in which the flow depends only on two spatial coordinates. For three dimensional (3-D) problems, where the flow depends on all three spatial coordinates, hexahedral, tetrahedral, pyramidal or prismatic shaped cells can be used. Geometries are usually created using computer aided design (CAD) software. The geometry, either a wireframe or solid model, is exported to the grid generation software program to create the C F D quality grid. A few packages have combined both functions of C A D geometry creation and mesh generation into a single interface. This phase of the CFD is referred to as preprocessing. Once the grid has been created, boundary conditions need to be applied. Pressures, velocities, mass flows, and scalars such as temperature may be specified at inlets; temperature, wall shear rates, or heat fluxes may be set at walls; and pressure or flow-rate splits may be fixed at outlets. Chapter 2: Literature Review 33 The component material transport properties, such as density, viscosity, and heat capacity, need to be prescribed as constant or selected from a database and can be functions of temperature, pressure or any other variable of state. Fluids can be modeled as being either incompressible or compressible. The viscosity of the fluid can be either Newtonian or non-Newtonian. In mass or heat transfer applications, binary diffusivities and thermal properties need to be defined as well. With the grid created, the boundary conditions and the physical properties specified, the calculations can start. The code will solve the appropriate conservation equations for all grid cells using an iterative procedure. Typical problems involve solving for mass conservation, momentum conservation along with a turbulence model, enthalpy, chemical species concentrations, local reaction rates and local volume fractions for multiphase problems. 2.6.1 C F D Modeling of stirred tanks In spite of the importance and variety of their applications, the problem of designing and scaling-up stirred tanks has been tackled mainly by means of semi-empirical methods. It is believed that a significant improvement on stirred tank design can be obtained by developing models which take into account the actual flow field inside the vessel. A n effort has been made during the last years towards the development of predictive methods based on computational fluid dynamics and capable of providing detailed information on flow and turbulence fields. At the same time, more and more accurate experimental data have been obtained on the mean flow field of mixing tanks, to serve as empirical input to predictive models or as validation benchmarks for the computational methods (Brucato et al, 1998) Chapter 2: Literature Review 34 In the past, the problem of modeling stirred tanks has been dealt with by excluding the impeller region from the computational domain and prescribing, on the basis of simple models and/or experimental information, the flow field near the impeller as a stationary boundary condition for the remainder of the computational domain (Middleton et ah, 1986; Brucato et al. 1989, 1990; Gosman et al, 1992; Ranade and Joshi, 1990; Ranade et al, 1989; Fokema et al, 1994) . These methods are called black box methods. As an alternative, the impeller has been modeled as a distributed source of momentum and turbulence quantities (Hutchings et al., 1989 ) The main limitation of the above approaches is their lack of generality. The near-field flow produced by a given impeller may depend on the overall tank geometry, so that the impeller can not be uniquely characterized. Also, different choices of impeller boundary conditions lead to significantly different predicted flow fields. The conclusion is that numerical simulations using this approach can be extended, at best, to geometries close to that for which the experimental data were obtained. A number of approaches have been proposed which allow the explicit simulation of the whole flow field without any recourse to empirical data. Takeda et al. (1993) used a multi-block technique in which the mesh block associated with the impeller rotated with it, while the outer block associated with the baffles was stationary. interpolation. Inter-block coupling was performed by However, few computational details, and no quantitative comparison with experimental data, were provided. Luo et al. (1993) solved the time-dependent form of the governing transport equations in two domains which were fixed to the respective frames of reference, the outer one being stationary Chapter 2: Literature Review 35 and the inner one rotating with the impeller. At the interface, the finite volume mesh was allowed to shear and/or slide to accommodate the relative motion. Later, Luo et al. (1994) observed that, although the interface conditions between inner and outer domain are generally unsteady, there is a particular radial location at which, for all practical purposes, steady flow conditions can be assumed. By using this fact, they computed the flow in a whole tank without paying the price of fully time dependent calculations; the steady state governing equations were solved in each reference frame, with interface coupling performed by velocity transformation. Results for the same reference geometry compared well with the experimental data. This approach is now available as a standard option in commercial CFD software such as F L U E N T and is called Multiple Reference Frame (MRF) technique. Murthy et al. (1994) described an alternative, sliding-mesh technique that has also been implemented as a standard option in F L U E N T called sliding meshes. Two grids are employed, one moving with the impeller and the other fixed to the tank. The two meshes interact along a surface of slip. The moving grid is allowed to slide relative to the stationary one, with no mesh distortion, and a conservative interpolation is used to obtain flow variables and face fluxes across the slip surface. Simulations are conducted in a time-dependent fashion and the laboratory reference frame was used for both grids. This technique is appropriate when a timeaccurate solution for rotor-stator interaction (rather than a time averaged solution) is desired. The sliding mesh model is the most accurate method for simulating flows where there is a strong interaction between stationary and moving parts (transient rotor stator interaction, impeller-baffles interaction). Since this model requires an unsteady numerical solution, it is computationally more demanding than the M R F model described above. Chapter 2: Literature Review 36 2.6.2. C F D Modeling of Pulp Suspensions The Theological model implemented in a CFD simulation must describe the behaviour of the pulp suspension reasonably well in order to be able to compare the calculated flow field with experimental data. Bakker (1993) was the first to model the fiber suspension flow in a pulp chest with side entering impeller using a black box method (impeller boundary condition method) and a combination of turbulent and laminar flow regimes in the tank. F L U E N T V4 (Fluent Inc.,Lebanon, NH). The C F D solver used was In this work, for every computational cell the computations were first performed as if the flow was turbulent using a k- e turbulence model. After that, it was checked if the total shear stress was indeed larger than the yield stress r . y If this was not the case, then the calculations for that cell were repeated as if the flow was laminar. The local apparent viscosity was then calculated from experimental shear rate - shear stress curves and the local shear rate was calculated by the solver in that particular cell. The shear rate-shear stress curve used in order to model the rheology of pulp suspension was not specific of any pulp type or particular consistency. Further more, numerical values for the yield stress or shear stress threshold were not explicitly given. The results of this CFD model were only compared to flow field visualizations for which the agreement was satisfactory. Later on, Roberg (1997) simulated the flow of paper pulp in a cylindrical stirred tank equipped with a side-entering mixer. They assumed laminar flow in the whole tank and used both a Newtonian and a non-Newtonian model for the rheology of pulp. A power law Chapter 2: Literature Review 37 equation/^ = Ky"~ was used to define the relationship between the shear rate and the viscosity l of the pseudoplastic model. The values of K and n used were 5000 mPas and 0.25, respectively, from experimental data of pulp suspension at 4% consistency. There was no specific information in this paper that describes how these values were measured. The propeller was simulated as a source of momentum, using estimates based on measurement of the stirring power and impeller speed. Commercial CFD software packages PHOENICS and FLOW3D were used. Velocity fields obtained from the Newtonian and non-Newtonian formulation were compared against each other but not against any experimental measurements. They concluded that the Newtonian flow formulation could only be accepted as an approximation or initial guess for the non-Newtonian formulation in order to facilitate convergence. Wikstrom and Rasmuson (1998) studied the agitation of pulp suspension in a 24 m cylindrical 3 tank equipped with a jet nozzle agitator using commercial C F D softwatre (CFX 4, ANSYS Europe Ltd, Oxfordshire) to calculate the flow field in the agitated pulp chest. They didn't explicitly simulate the impeller and its effects were modeled by imposing boundary conditions to the external flow: The outlet of the flow jet agitator was treated as a velocity boundary condition, and the inlet to the agitator as a mass flow boundary condition. This treatment o f the nozzle boundary condition allowed the fluid to circulate in the tank. The pulp suspension was described as Bingham fluid using the yield stress measured, r = 180 Pa, and the flow inside the chest was modeled as being laminar. Experimental velocity measurements in the tank were taken using a sonar Doppler velocity meter at four different positions. They showed that the flow field obtained from the C F D calculations increasingly deviated from the velocities measured (with sonar doppler velocimetry) as the distance from the impeller increased. Since Chapter 2: Literature Review 38 the C F D calculation underestimated the velocity far away from the agitator, they included a shear-thinning parameter through a Herschel-Bulkley model. The result, however, showed no significant deviation from the Bingham model and they suggested that these models were too rigid to describe the real behaviour of the yield stress along a stream line which could explain the deviations. 2.6.3 Objectives of this work The objective of this research is to simulate pulp mixing in a scale model stock chest using CFD commercial software (Fluent 6.1) and a simple model for the pulp suspension rheology using measured values of the suspension yield stress. The calculated velocity field will then be used to obtain the dynamic response of the simulated chest to a random binary input signal for which experimental responses have been previously measured. These dynamic data and measurable variables such as the power input will be then compared. This procedure will allow validation of the CFD model and indicate further avenues for improvement and refinement of CFD simulations. Chapter 3: Experimental 39 3. E X P E R I M E N T A L 3.1. Experimental Setup The experimental setup used is the one designed and built by Ein Mozaffari (2002) to study macroscale mixing and the dynamic behaviour of agitated pulp stock chests. A schematic diagram of this setup is shown in Figure 3.1. Pulp suspension was pumped from the feed tank (1 m ) to the plexiglas stock chest (0.2 m ) and back to the discharge tank (1 m ). Both the feed and discharge pumps are progressive cavity pumps (Moyno Industrial Products, Springfield, OH) equipped with variable frequency drives. The pulp suspension inside the feed tank was agitated using a 63.5 cm (25") top-entering Chemineer (Dayton, OH) A310 impeller driven with a 5-hp variable speed drive. The stock chest is rectangular and its dimensions are: width, W=40cm, height, Z=70cm, and length, L=70cm. The length of the chest could be changed to 40 cm, 53 cm or 70 cm. A 3-hp variable speed motor was used to drive the impeller. The impeller used is a scaled version of Maxflo impellers (Chemineer Inc.) with a diameter of 16.5cm (6.5") identical to those used in industry for stock agitation. The impeller torque and speed were measured using an inductive-rotary torque transducer with an encoder disk (model 0411E50, Staiger Mohilo, Germany). To study the dynamic response of the stock chest, a saline solution was injected through a computer controlled solenoid valve (8300D9u, Pacific Controls Ltd. Vancouver, BC) into the pulp stream using a metering pump (C121-362SI, LMI Milton Roy, Acton, MA). Conductivity variations in the input and output streams were measured using flow-through conductivity sensors (2220254, RoseMount Analytical, Irvine, CA). The pipe diameter was 5.1cm (2") for pulp suspension handling and 1 cm (3/8") for tracer injection. Input/Output conductivity data Chapter 3: Experimental 40 was obtained using a 16 channel data acquisition card (PCI-6024E, National Instruments, Austin, TX) and LabVIEW software. Liquid Tracer • —O' 1 1 r~~l i *1 5 A Discharge Tank (V=1m ) 3 s, 7 1 : Transmitter 2 : Indicator 3: TacbojTWtw 4 : torque Transducer 5: Magnetic Flowmeter 6: Conductivity Sensor 7 : Pump 8 : Solenoid Valve Figure 3.1: Schematic of the experimental setup. 3.2 Impeller specifications The impeller used is a Maxflo Mark II with a diameter of 16.5 cm. Figure 3.3 shows a photo of the impeller. Table 3.1 summarizes the impeller specifications (Pitch, C factor, power number and flow number). The power number given is the one for fully turbulent flow calculated as N„ P pN D 3 5 = 0.276 (3.1) Chapter 3: Experimental 41 where P, p, N and D are power, fluid density, impeller speed and impeller diameter respectively. This value is in the range of typical values of N PJURB (0.1-0.5) for axial flow impellers (AIChE Equipment Testing Procedure, 2001). N«m • ——~—> Figure 3.2: Scale model chest dimensions The dimensionless flow number (NQ) for this Maxflo impeller is N = 0 6 4 3 ND ( 3 V 2 ) J where Q peu is the pumping rate (Uhl and Gray, 1966). This value is in the range of typical im er values (0.4-0.8) for axial flow impellers (AIChE Equipment Testing Procedure, 2001). The C factor is the constant that relates the impeller momentum flux Mo, with the rotational speed N and the impeller diameter D (see equation 2.24) and is 0.41. 42 Chapter 3: Experimental Table 3.1: Maxflo impeller specifications Diameter 16.5 cm Pitch C factor Power number Flow number 0.44 0.41 0.276 0.643 Figure 3.3: M a x f l o impeller used in the scale model chest 3.3. Materials The pulp suspensions were prepared from fully bleached kraft pulp (FBK) with a fiber weighted averaged length of 2.47mm. The pulp was obtained in wet sheet form (about 50% water) from Pacifica Paper,BC These sheets were soaked in water and defibred in a disintegrator and prepared to the desired mass concentration of 3.3%. Table salt from Galloway's Wholesale (Richmond, BC) was used to prepare the saline solution. Chapter 3: Experimental 43 3.4. Test procedure The test procedure for obtain dynamic data was developed by Ein Mozafari (2003). Each test was performed by exciting the system and observing its dynamic response over a time interval typically five times longer than the nominal residence time of the chest at the respective flow rate. A saline solution was injected through a computer controlled solenoid valve into the pulp feed stream and the conductivity variations in input and output streams were measured using flow-through conductivity sensors. These input and output signals were recorded in a computer. The excitation provided was limited to binary sequences as the input signal was controlled by an on-off solenoid valve. The system was excited with a frequency-modulated random binary signal designed by Ein Mozzafari (2003) to excite the system's dynamics at the appropriate frequencies and therefore obtain a more accurate estimation of the parameters of the dynamic model developed for this system. The control of the solenoid valve and data acquisition (conductivity variations in input and output streams) was accomplished using LabVIEW software. Input/Output conductivity measurements were recorded every 1-5 seconds. For instance, the sampling time was one second sXQ- 37.1 L/min and it was 5 seconds at Q = 7.9 L/min while the suspension volume in the chest was kept constant (V= 106 L). Impeller torque and speed were recorded for each test. 3.5. Experimental conditions Various tests were conducted by Ein-Mozzafari (2002) at the different operating conditions listed in table 3.2 for the pulp suspension atC = 3.3% . A subset of these tests was also used m for validation of the C F D model. The effect of pulp feed and exit locations, was also tested for the two chest configurations shown in figure 3.4. Chapter 3: Experimental 44 Table 3.2. Summary of operating conditions studied in the experiments Operating condition Experimental range studied Impeller speed 700 rpm - 1750 rpm Pulp flow rate 7.9 L/min-37.1 L/min Fibre mass concentration 3.3% Figure 3.4: Input-Output configurations studied 3.6. Power Measurement Impeller torque and speed were measured using an inductive-rotary torque transducer with an encoder disk (model 0411IE50, Staiger Mohilo, Germany) mounted on the rotor shaft with two flexible couplings at each end to accommodate radial, angular, and axial misalignment. The Chapter 3: Experimental 45 bearing friction was measured by operating the system with the chest empty. This friction torque was subtracted from all measured torques. The power input to the impeller UP) is obtained from torque (M) and impeller speed (AO measurements using P = 2TTNM and was accurate to ± 0.05 Nrh. (3.3) Chapter 4: CFD Model Development 4. C F D M O D E L 46 DEVELOPMENT 4.1. Introduction C o m p u t a t i o n a l F l u i d D y n a m i c s ( C F D ) is a sophisticated analysis t o o l that a l l o w s us to m a t h e m a t i c a l l y s o l v e the f l o w i n a g i v e n f l u i d process. In this chapter, the d e v e l o p m e n t o f the C F D m o d e l o f the scale m o d e l chest p r e v i o u s l y d e s c r i b e d w i l l b e o u t l i n e d i n detail. T h e c o m m e r c i a l software F L U E N T 6.1 (Fluent Inc., L e b a n o n , N H ) w a s u s e d to n u m e r i c a l l y solve the N a v i e r - S t o k e s equations. GAMBIT, the m e s h i n g t o o l p a c k a g e d w i t h F L U E N T C F D software, w a s u s e d to generate a 3 D g e o m e t r i c a l m o d e l a n d m u l t i b l o c k m e s h o f the m i x i n g chest. T h e p u l p s u s p e n s i o n w a s treated as a B i n g h a m plastic. A m u l t i p l e reference frame a p p r o a c h w a s u s e d (the equations o f m a s s a n d m o m e n t u m c o n s e r v a t i o n are s o l v e d i n different frames o f reference) w i t h c o u p l i n g b e t w e e n frames m a d e u s i n g a v e l o c i t y transformation d e s c r i b e d a h e a d i n this chapter (session 4.3). T r a c e r analysis w a s u s e d to obtain the d y n a m i c response o f the s y s t e m as w a s d o n e i n the experiments. In o r d e r to d o this w i t h the n u m e r i c a l m o d e l , a user d e f i n e d f u n c t i o n that defines the excitation s i g n a l w a s l i n k e d to the solver. T r a n s i e n t species c a l c u l a t i o n w a s p e r f o r m e d i n order to a c q u i r e the specific values o f tracer c o n c e n t r a t i o n at the v e l o c i t y inlet a n d the outflow, a l l o w i n g the c o l l e c t i o n o f input/output data for situations i d e n t i c a l to the e x p e r i m e n t s p e r f o r m e d i n the laboratory. 4.2. General Transport Equations and Discretization T h e general transport e q u a t i o n f o r a v a r i a b l e O i n a C a r t e s i a n c o - o r d i n a t e s y s t e m c a n b e written as: |(^) A(p / *)-A[ .|* ) , + t 1 r + S . (4.,) Chapter 4-.CFD Model Development 47 where T$ is the diffusion coefficient (equal to the viscosity in the momentum equations). The variable O can be any field variable, e.g. velocity components. The first term is the unsteady change of O, the second is the convection of <D, the third is the diffusion of O and the last is the source term of O. Equation (4.1) is then integrated in space and discretized according to the differential scheme chosen. F L U E N T uses a control volume based technique to convert the governing equations to algebraic equations that can be solved numerically. This technique consists of integrating the equations about each control volume, yielding discrete equations that conserve each quantity on a control-volume basis. The steady-state conservation equation for transport of a scalar quantity O can be written in integral form for an arbitrary control volume as follows: (4.2) v where p = density, v = velocity vector, A = surface area vector, F = diffusion coefficient for 0 O , V ® = gradient of O and - source of O per unit volume. Equation (4.2) is applied to each control volume, or cell, in the computational domain. Discretization of equation (4.2) for a given cell yields fpZ f / f O f •A = f / f r 0 (vo) n • A + SV f 0 (4.3) Chapter 4:CFD Model Development Where N faces p Vf-Af f 48 = number of faces enclosing the cell, <& = value of ® convected through face / f = mass flux through the face, Aj = are of face / , ( V ® ) = magnitude of V ® normal n to face/and V= cell volume. The equations solved by F L U E N T take the same general form as the one given above. By default, F L U E N T stores discrete values of the scalar ® at the cell centers. However, face values ® y are required for the convection terms in equation (4.3) and must be interpolated from the cell center values. This is accomplished using an upwind scheme (the face value is derived from quantities in the cell upstream, relative to the direction of the normal velocity v ). n The diffusion terms in equation (4.3) are central differenced and are always second-order accurate. A more detailed definition of the upwind schemes used can be found in Appendix A. Equation (4.3) will, in general, be non-linear with respect to an unknown scalar variable ® at the cell center as well as the unknown values in surrounding neighbor cells. A linearized form of equation (4.3) can be written as: V*>/.=2>,A+6 (-) 4 4 nb where the subscript 'nb' refers to neighbor cells, and a and a are the linearized coefficients p nh for O p and ® „ 6 . Similar equations can be written for each cell in the grid. This result in a set of algebraic equations with a sparse coefficient matrix that is solved by F L U E N T using a point implicit (Gauss-Seidel) linear equation solver in conjunction with an algebraic multigrid method (AMG), and a velocity field is obtained (Ferzieger and Peric, 1996; F L U E N T Inc, 2003). 49 Chapter 4:CFD Model Development Because of the non linearity of the equation set being solved, it is necessary to control the change of 0. This is achieved by under-relaxation, which reduces the change of O in each iteration. The new value of O within a cell depends upon the old value & id, 0 the computed change in O, AG>, and the under relaxation factor, a, as follows: (4.5) When unstable or divergent behaviour is observed, the under-relaxation factors can be reduced in order to control the update of computed variables at each iteration. Lower under-relaxation factors can facilitate stability of the solution but this means that more itereations will be required to reach the final solution. 4.2.1. Solution Convergence One of the most used ways to monitor convergence of the solution is checking the residuals. At the end of each solver iteration, the residual sum for each of the conserved variables is computed and stored, thus recording the convergence history. On a theoretical computer with infinite precision these residuals would go to zero as the solution converges. On an actual computer, the residuals decay to some small value and then level out. The residual 7?* computed by FLUENT's segregated solver is the imbalance in equation (4.4) summed over all the computational cells P. This is referred to as the "unsealed" residual. It may be written as (4.6) Chapter 4:CFD Model Development \ 50 In general, it is difficult to judge convergence by examining the residuals defined by equation (4.6) since no scaling is employed. F L U E N T scales the residual using a scaling factor representative of the flow rate of <D through the domain. This "scaled" residual is defined as cells? nb (4.7) cellsP For the momentum equations the denominator term a <3? is replaced by a u p p p p , where u is the p magnitude of the velocity at cell P. This residual is the default displayed by F L U E N T . For the continuity equation, the unsealed residual for the segregated solver is defined as R = \rate of mass creation in cell P\ c (4.8) cells P and the segregated solver's scaled residual is defined as R° /A Q\ iteration N iteration 5 The denominator is the largest absolute value of the continuity residual in the first five iterations (FLUENT Inc, 2003) 4.2.2. Discretization of the momentum equation The x-momentum equation, for example, can be obtained by setting 0=w: A au = J_a u P nh nb + £ p A -i + S f (4.10) nb Once the pressure field and face mass fluxes are obtained as part of the solution procedure, equation 4.10 can be solved. F L U E N T uses a co-located scheme, whereby pressure and Chapter 4:CFD Model Development 51 velocity are both stored at cell centers. However, equation (4.10) requires the value of the pressure at the face between adjacent cells. Therefore, an interpolation scheme is required to compute the face values of pressure from the cell values. 4.2.3. Discretization of the continuity equation The steady state continuity equation in integral form: Q p v- d A = 0, may be integrated over a control volume to yield the following discrete equation zZ f f=° J (-) A where J is the mass flux through face f, 4 pv . f n In the solution procedure, the continuity n equation is used as an equation for pressure. However, pressure does not appear explicitly in equation (4.11). The SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) family of algorithms (Pantankar, 1980) is used for introducing the pressure into the continuity —> equation. It is also necessary to relate the face values of the velocity v„ to the stored values of the velocity at the cell centers. Momentum-weighted averaging, using weighting factors based on the a p coefficient from equation (4.10) is performed (FLUENT Inc, 2003). Using this procedure the face flux, J/ is given by: J =J f where p c0 and p cX f + d (p -p A f c0 c (4.12) are the pressures within the two cells on either side of the face, and J contains the influence of velocities in these cells. f The term df is a function of a , the average of the momentum equation a coefficients for the cells on either side of face/ p p Chapter 4:CFD Model Development 52 4.3. Flow in a Rotating Frame The flow in the scale mixing chest is unsteady in an inertial frame of reference because the rotor/impeller blades sweep the domain periodically. However, in the absence of stators or baffles, it is possible to perform calculations in a domain that moves with the rotating part. In this case, the flow is steady relative to the rotating (non-inertial) frame, which simplifies the analysis (Luo et al., 1994) F L U E N T has the ability to model flows in an accelerating reference frame. In this situation, the acceleration of the coordinate system is included in the equations of motion describing the flow. In the implementation of the multiple reference frame feature, the calculation domain is divided into subdomains, each of which may be rotating with respect to the inertial frame. The governing equations in each subdomain are written with respect to that subdomain's reference frame. Thus, the flow in stationary subdomains is governed by the continuity and momentum equations as follows: ^ ot + V-(pZ) = S (4.13) m —(pu) + V-(puo) = -Vp + V-(r) + pg+F (4.14) dt = -> -> where p is the static pressure, r is the stress tensor (described below), and p g and F are the gravitational force and external body forces respectively. dependent source terms. The stress tensor is given by F can also contain other model- Chapter 4-.CFD Model Development 53 — T = JU V y + V u V J 3 V-uI (4.15) where ju is the molecular viscosity, I is the unit tensor, and the second term on the right hand side is the effect of volume dilation. The flow in the rotating subdomain is governed by the "modified" (additional terms include the acceleration of the fluid) continuity and conservation of momentum equations as follows: (4.16) dt — (p o) + V • (p v v) + p(Q xu) = - V > + V • ( r ) + p g+ F dt r v r (4.17) refers to the relative velocity and can be related to the absolute velocity by the following equation: u,=u-(Qxr) (4.18) Here, Q is the angular velocity vector (that is, the angular velocity of the rotating frame) and r is the position vector in the rotating frame. At the boundaries between two subdomains, the diffusion and other terms in the governing equations in one subdomain require values for the velocities in the adjacent subdomain. F L U E N T enforces the continuity of the absolute velocity, v, to provide the correct neighbor values of velocity for the subdomain in consideration. Chapter 4:CFD Model Development 54 4.3. Scale Mixing Chest Model 4.3.1. Geometry A mixing chest with identical geometry than the one used in the experiments was modelled in GAMBIT . The impeller geometry was created based on measurements taken from the Maxflo 1 Mark II impeller with diameter 10.2 cm (4"). These measurements were scaled up by a factor of 1.62 in order to properly define the Maxflo Mark II impeller with diameter 16.5 cm (6.5") which is installed in the experimental set up. Table 4.1 shows the coordinates of the points that define the blade curvature and specific parameters of orientation. The edges of the blade were obtained from spline interpolations on every 3 points of the 11 points (3D coordinates) measured along the blade. The wire frame was then imported to GAMBIT and the blade face was created. The same blade was copied and rotated 120° each time. Figure 4.1 shows the impeller wire-frame in GAMBIT. Figure 4.1: Maxflo Mark II (6.5") impeller wire-frame in GAMBIT. GAMBIT, the meshing tool packaged with F L U E N T , is a software package designed to build and mesh models for computational fluid dynamics (CFD) and other scientific applications. Chapter 4:CFD Model Development 55 For a mixing tank, with a single impeller, one can define a rotating frame that encompasses the impeller and the flow surrounding it, and use a stationary frame for the flow outside the impeller region. Thus, in the geometrical model, the impeller is located in a cylinder that includes the fluid that is in the immediate vicinity of the hub and blades and this sub-volume corresponds to the rotating subdomain. The rest of the domain is stationary. Figure 4.2 illustrates the geometry of the scale model mixing chest and the two subdomains used for the numerical solution of the governing equations Table 4.1. 3-D coordinates of the points along the blade wire frame (origin of the coordinate system is at the center of the hub) Impeller diameter =10.2 cm Impeller diameter = 16.5 cm Points X mm Z mm Points 1 16.720 12.165 0.000 1 27.086 19.707 0.000 2 25.570 11.860 -0.600 2 41.423 19.213 -0.972 3 32.200 17.050 -0.350 3 52.164 27.621 -0.567 4 34.630 25.360 0.100 4 56.101 41.083 0.162 5 31.680 33.340 1.200 5 51.322 54.011 1.944 6 25.500 40.000 2.500 6 41.310 64.800 4.050 7 16.600 45.250 4.500 7 26.892 73.305 7.290 8 8.500 46.000 9.500 8 13.770 74.520 15.390 9 0.000 42.000 6.900 9 0.000 68.040 11.178 10 2.000 30.500 8.000 10 3.240 49.410 12.960 11 4.840 20.000 5.600 11 7.841 32.400 9.072 HUB Di.imeter y mm 42mm. I hickncv. - Idmm HI 1$ X mm y mm ^ mm Diameter ^ 26mm,'Thickness = 68.05mm The entrance and exit pipes in this geometrical model are different from the real experimental setup. This discrepancy arises from the boundary conditions selected in F L U E N T . The exit pipe is defined as an outflow, which definition in F L U E N T implies a zero diffusion flux for all Chapter 4:CFD Model Development 56 flow variables at the boundary and this is only physically approached in a situation of fully developed flow. Thus the pipe is modeled with an L/D >10. Figure 4.2: Mixing chest geometry in GAMBIT. The rotating sub-domain is shaded in the views above 4.3.2 Meshing the Geometrical Model A multiblock three dimensional mesh was generated in GAMBIT. The computational domain was divided into seven geometrical blocks. Five of these blocks were meshed using a structured mesh, consisting of hexahedral elements logically connected. The impeller is located Chapter 4:CFD Model Development 57 in a cylinder that encompasses the impeller hub and blades and the flow surrounding it. The computational domain (where the governing equations are solved) is defined by the volume occupied by the fluid in the mixing chest. The impeller hub and blades are modeled as rotating walls. The geometry in this cylinder is quite irregular so it places a high demand on geometrical flexibility of the mesh and thus was meshed using an unstructured mesh, consisting of tetrahedral elements. In order to create a conformal mesh (i.e., the grid node locations are identical at the boundary where two zones meet), the block that is adjacent to the impeller block was also meshed using an unstructured mesh. The multiblock division of the domain is illustrated in Figure 4.3. Block-structured mesh Block-structured mesh Block-structured mesh Block-structured Block-unstructured mesh mesh Y Figure 4.3: Multi-block division of the computational domain Chapter 4:CFD Model Development 58 A finer mesh was generated around the impeller, entrance and exit pipes and tank walls using boundary layers created in GAMBIT (increased mesh density near the walls). The most critical part was the mesh around the impeller. In order to have a very refined mesh in the immediate vicinity of the blades, the sufficient amount of nodes that properly define the curvature of the blades was created on the impeller edges and a size function was used to control the mesh density. This feature allows the mesh elements to grow slowly as a function of the distance from the impeller blades. The mesh around the impeller is shown in figure 4.4. Figure 4.4: Enlarged view of the size function used around the impeller. The mesh density in the block adjacent to the rotating subdomain was also controlled by means of a size function in which the mesh elements edges grow slowly as a function of the distance to the faces of the cylinder that encompasses the impeller. The final three dimensional mesh of the model had originally 112,324 computational cells and is shown in figure 4.5. Chapter 4:CFD Model Development 59 4.3.3 Problem Setup for Multiple Reference Frames in F L U E N T The MRP model is a steady-state approximation in which a subdomain of the full volume is specified to move at a fixed angular speed with respect to the full domain. While this approach is clearly an approximation, it can provide a reasonable model of the time-averaged flow for applications in which rotor-stator interaction is relatively weak like the mixing chest examined. In this section, the problem set up will be outlined step by step. Figure 4.5: Three dimensional mesh of the scale mixing chest. Chapter 4:CFD Model Development 60 1. Grid Setup The geometry and mesh created in GAMBIT is imported into F L U E N T in order to develop the numerical solution. Once the grid information was acquired by F L U E N T it was then scaled to the actual geometrical dimensions of the model and it could be improved by automated procedures called smoothing and swapping. These tools usually increase the quality of the final numerical mesh as follows: Smoothing repositions the nodes and face swapping modifies the cell connectivity, both in order to achieve improvements in cell skewness (an indication of how close the cell is to a perfect cube, i. e. skewness = 0, describes a perfect equilateral cube and skewness =1.0 describes a poorly shaped element). F L U E N T tries to move interior nodes to improve the quality of cells with skewness values greater than the specified "minimum skewness". Four iterations were performed using a minimum skewness of 0.8 and then face swapping was performed until the number of faces swapped decreases to 0. The face swapping process finds configurations of three cells sharing an edge and converts them into two cells sharing a face in order to decrease skewness and cell count. After this, the minimum skewness was decreased to 0.6 and the smoothing/swapping procedure repeated. 2. Material(s) Definition F L U E N T has a complete data base with different materials and their physical properties. Initial simulations were performed for liquid water for which the material properties exist in the F L U E N T data base. The Cm = 3.3% pulp suspension was introduced as a new material with non-Newtonian viscosity defined through the Herschel-Buckley model for Bingham plastics (Eqn.4.19). The inputs for this model are the consistency index k, power-law index n, yield stress threshold^ , and yielding viscosity p. . 0 Chapter 4: C F D Model Development 61 r +k y r - {/Mo) (4.19) These parameters were changed as a function of consistency of the suspension. The work done was focused on pulp suspension at C = 3.3% for which the yield stressr = 475 ± 62 Pa was m y measured by Ein-Mozzafari (2004). The yield stress was varied in the 95% confidence interval of this measurement (T = 337 Pa andr^ = 950 Pa). y An initial yielding viscosity value ju = a 100 ±43 kg/ms, based on experimental data from Bennington (1988) for 1.6% C semim bleached kraft pulp was chosen. The solution sensitivity with respect to ju was studied based o on the output response and the dynamic parameters calculated from the CFD simulation data. These data has been collected in appendix D. 3. Numerical method and Viscous Model The numerical scheme selected in FLUENT was the segregated solver. Using this approach, the governing equations are solved sequentially. Because the governing equations are nonlinear (and coupled), several iterations of the solution loop must be performed before a converged solution is obtained. Each iteration consists of the steps illustrated in Figure 4.6. The pressure-velocity coupling method used was SIMPLE. When the calculations were performed with water, the momentum conservation equations were solved in the turbulent regime using the RNG k-e turbulent viscous model (see appendix C). The calculations performed with pulp suspension were done using the laminar viscous model specified in FLUENT 6.2 (Fluent Inc., 2003). 62 Chapter 4:CFD Model Development Update properties. 1 Solve momentum equations. - f Solve pressure-correction (continuity) equation. Update pressure, face mass flow rate. Solve energy,, species,, turbulence, and other scalar equations. Converged? -*\ Stop Figure 4.6: Overview of the segregated solution method 5. Boundary Conditions The boundary conditions were specified at the locations illustrated in Figure 4.7. For both the rotating and stationary subdomains, an angular velocity (fi) in their respective frame of reference and the axis about which it rotates is specified. The rotation axis origin is at the center of the impeller hub and the orientation is the positive Z direction. The impeller speed is specified in the M R F subdomain and the impeller blades are treated as walls that move with zero relative velocity. Likewise, the chest walls are given a velocity of zero in the absolute reference frame. Chapter 4-.CFD Model Development 63 The velocity inlet is defined in the absolute frame as a flat profile with magnitude Ivl = ^0/ 1 1 / 2 KD and negative Y orientation. The exit is defined as an outflow boundary condition. The boundary conditions used by FLUENT at outflow boundaries are a zero diffusionfluxfor all flow variables, which means that the conditions of the outflow plane are extrapolated from within the domain and have no impact on the upstream flow, and an overall mass balance correction. This condition is approached physically in fully developed flows, thus the exit pipe was modeled with an L/D ratio of 10. 6. Solution Strategies A high degree of rotation introduces a large radial pressure gradient which drivesflowin the axial and radial directions, thereby setting up a distribution of the swirl in the field. This coupling may lead to instabilities in the solution process, and hence require special solution techniques to obtain a converged solution. Some techniques used include the following: > The pressure interpolation scheme selected was PRESTO! (PREssure Staggering Option) which is well-suited for steep pressure gradients involved in rotatingflows.The scheme uses the discrete continuity balance for a "staggered" control volume about the face to compute the "staggered" face pressure (Patankar, 1980). > The calculations were started with a low rotational speed and gradually increase it in order to reach the final desired operating condition. The solution was initialized from the velocity inlet. 64 Chapter 4: C F D Model Development Velocity Inlet Stationary subdomain axis of rotation Stationary wall with zero shear Moving wall at zero Outflow .relative-Velocity Moving wall at zero-absolute velocity MRF angular velocity axis of rotation t Figure 4.7: The boundary conditions specified in the domain > Convergence was monitored dynamically by checking the residuals for each conservation variable: mass, x-velocity, y-velocity and z-velocity for the laminar simulations and all the above plus k and s for the turbulent simulations. The solutions were considered converged when the scaled residuals were below 1 x 10" . 5 > The under-relaxation factors for the velocities were initially reduced to 0.5 and depending upon convergence of the residuals they were further reduced up to 0.2. > A first order upwind scheme was initially used for the momentum equations (the face value 0/is set equal to the cell-center value of 0 in the upstream cell). After convergence was Chapter 4:CFD Model Development 65 reached, the solution was used as the initial conditions for calculations made with a second order upwind scheme (higher-order accuracy is achieved at cell faces through a Taylor series expansion of the cell-centered solution about the cell centroid). While the first-order discretization generally yields better convergence than the second-order scheme, it will usually yield less accurate results, especially on unstructured grids (FLUENT, 2003). > In order to ensure that the mesh was sufficiently refined in regions with large gradients of pressure and swirl velocity, a solution-adaptive refinement of the mesh was incorporated. Assuming the greatest error occurs in high-gradients regions, the available physical features of the evolving flow field were used to drive the grid adaption process in an automated dynamic refinement of the mesh on the cell zones marked with high velocity gradients. 7. Power input calculation The power input to the impeller (P) is obtained from the calculated torque (M) and impeller speed (TV) defined using P = 2nMN (4.20) The moment vector about a specified center for the impeller wall zone is computed by summing the product of the force vectors for each face with the moment vector /. e., summing the forces on each face about the moment center. The moment center is defined in fluent as the center of the impeller hub, a point with coordinates x = 0 cm, y = 12 cm and z =7 cm. The actual components of the pressure, viscous and total moment can be automatically calculated in FLUENT. For example, the net pressure force vector is computed as the vector sum of the individual force vectors for each face as follows: 66 Chapter 4:CFD Model Development (4.21) A where « is the number of faces, A is the area of the face, and n is the unit vector normal to the face. For closed domains, as the mixing chest modeled in this case, the additional term introduced by the reference pressure cancels. 4.4. Mixing Dynamics study with Tracer Analysis The flow field calculated in the CFD simulations was then used to obtain the system's dynamic response in the same manner as it was done in the experiments. When a tracer is added to a fluid in a mixing tank, the transient calculation can be done exclusive of the flow field calculation, given that the properties of the tracer and the background liquid are similar. Once the flow field for the background fluid is satisfactorily converged, the tracer can be introduced and transient species solution for the tracer is performed. Since the properties of the pulp suspension in the chest will not show a significant change with the addition of the tracer, all other equations can be disabled and only species transport of tracer can be solved. The conservation equation solved by F L U E N T takes the general form: 8_ dt (pY ) + V-(puY ) l i = -V-J R +S i+ i i (4.22) where 7, is the local mass fraction of species i, R is the net rate of production of species i by t chemical reaction and S is the rate of creation by addition from the dispersed phase plus any t user-defined sources. Since there is no chemical reaction between the tracer and the pulp -> suspension nor continuous generation of the tracer, R and S were both set to zero. t t J, is the 67 Chapter 4:CFD Model Development diffusion flux of species i, which arises due to concentration gradients. In laminar flows, F L U E N T uses the dilute approximation, under which the diffusion flux can be written as X=-PD,^Y, Here D i>m (4.23) is the diffusion coefficient for species i in the mixture. In the calculations, the diffusion term was neglected and the tracer distribution was obtained based on the convection defined through the flow field previously obtained, thus D im was set to zero. The convection component of the net transport of species at the inlet is fixed by the inlet species mass fraction that is specified in the velocity inlet boundary condition definition. At the walls, a zero flux boundary condition was applied for all species. In order to simulate transient chemical mixing, F L U E N T solves the species conservation equation in time-dependent form. For transient simulations, the governing equations must be discretized in both space and time. Temporal discretization involves the integration of every term in the differential equations over a time step At. The details concerning the integration of transient terms can be found in appendix A . The time step size was chosen so the number of iterations needed to reach convergence (scaled residuals were below lxlO" ) at each time step 5 was about 20-30. The concentration of the tracer at inlet and outlet was monitored with time in order to obtain an input/output data set identical to the one collected in the experimental runs. Initially, the system's dynamic response to a step change in tracer concentration was obtained. After this, in order to specify tracer injection only for some particular time intervals, a User Defined Function (UDF) like the one shown in figure 4.8 was linked to the F L U E N T solver. Chapter 4:CFD Model Development 68 : This function defines the frequency-modulated random binary input signal used in the experimental runs. Jinclude ^udf.h>"^^ DEFINE_PROFILE ( t r a e e r _ m f _ p r o f l i e , t, i) t real . : flow_time face_t = CURRENT_TIME; f; taegm_f_loop if (f, t) (flow_tlme > 932) F_PROFILE(f, else if else if else if else it else if else if else if else if (flow time > t, F_PROFILE(f, t, If low t i m e 843) > F_PROFILE(f, (flow time > F_PROFILE(f, (flow time > F_PROFILE(f, (flow time. > F_PROFILE(f, (flow time > F_PROFILE(f, (flow time > F_PROFILE(f, (flow time i) = 0 .2 7;. i) = q 165; i) = o 2 7;. 0 165; 0 27; 0 165; .0 27; 165; 907) t, 7 68)'.. t, i) 610) t, i) = 474) t, i) 3 87) t, i) 162). t, i) - 0 0. : 2 7 ; > 2 9) F_PROFILE(f, t, i) = F t, i) = 0 else PROFILE(f, 165; } e n d _ f _ l o o p ( f , t) } • Figure 4.8: U D F code of the frequency-modulated random binary input signal (for Q = 37.1 L/min) 69 Chapter 5: CFD Simulation Results 5. C F D S I M U L A T I O N R E S U L T S 5.1. Introduction In the CFD model, the geometry of the experimental mixing chest was captured precisely and the governing equations of the flow occurring in the system were solved as outlined in the previous chapter. Qualitatively, the flow patterns obtained from the numerical solution show a good agreement with the ones observed in the Plexiglas mixing tank. Detailed velocity information and the location of poor mixing zones within the agitated chest became available through the C F D simulation. The power input calculated in the simulations was directly compared with the one measured in the experiment and found to agree to within ±12%. The velocity field obtained from the simulation was used to obtain the system's dynamic response to a frequency-modulated random binary input signal. Input/output conductivity data sets collected in the experimental runs and obtained from the numerical solutions were plotted together and showed good agreement. These input/output data sets were also used as input to the dynamic model developed by Ein-Mozzafari (2003) and the parameters calculated were then compared with the parameters measured experimentally under identical conditions. 5.2. Newtonian Fluid Simulation (Water) Initial simulations were performed with water, a fluid with Newtonian rheology. The experimental conditions simulated are the same ones used in the experiments with the pulp suspension. In this case, the flow within the mixing chest is expected to be turbulent; (the Reynolds number calculated based on Eqn.(2.12) was Re > 10 ) therefore a turbulence model 4 Chapter 5: CFD Simulation Results 70 was used to solve the governing equations. The turbulence model chosen was the RNG-based k-e turbulence model. This is a semi-empirical model based on the transport equations for the turbulence kinetic energy (k) and its dissipation rate (s). Detailed information about the model equations can be found in Appendix B . The R N G model is more responsive to the effects o f rapid strain and stream line curvature, which explains its superior performance for rotating flows such as the mixing chest in consideration. The convergence criteria used was to iterate until the scaled residuals for each transport equation were below 1 x 10~ . The simulation results are shown in the following figures for one 4 of the operating conditions used in the experimental runs: Input flow, Q=30 L/min, and rotational speed, N=1250 rpm. Figure 5.1: Velocity vectors colored by velocity magnitude (m/s), for water with £>=30L/min and N=1250 rpm Chapter 5: CFD Simulation Results /1 Figure 5.1 shows the velocity vectors colored by velocity magnitude in m/s in the whole domain. The flow field obtained agrees with the axial impeller flow patterns described in chapter 2. The vectors in red have the highest velocity magnitude, and they are located near the impeller blades as expected. There is a turbulent region with high shear near the impeller where the mixing occurs as the high velocity streams come into contact with adjacent low velocity streams. Slow moving fluid (vectors colored in light blue/green) is entrained by the flow of the high velocity streams. The vectors colored in blue have the lowest velocity magnitude and they are located mostly near the walls and at the corners of the mixing chest. These poor mixing zones identified in the simulation were also present in the experiments as determined by observation of significant slower suspension motion at the corners of the Plexiglas tank. I 1.20 c-D 1.146'] 1.09c OQ 1.03c 0 0 9.65c 01 9.05ft-Dl 8.4 6 c 01 7.86c D£5^-5: 7.26c 0T"5 6.66e- D1 6.06c 01 5.46c 01 4.86c 01 4.26c 01 3.67e- D1 3.07c 01 2.47c 01 1.87c Dl 1.27c 01 6.71c 02 7.20cD3 Figure 5.2: Velocity vectors colored by velocity magnitude (m/s) on the plane x = 0 cm, for water with <2=30L/min and N =1250 rpm 72 C h a p t e r 5: CFD Simulation Results More detailed information can be accessed by examining the velocity vectors on different planes in the domain that are not accessible by simple observation. For instance, Fig. 5.2 illustrates the velocity vectors colored by velocity magnitude in m/s in a plane located at the center line of the mixing chest (x = 0 cm). This view shows a predominant axial velocity component that results as the fluid is drawn in from behind the impeller and forced down to the end of the tank (wall opposite to the impeller). When the flow encounters this wall it moves upwards to the surface of the chest due to the impeller pumping action, it then moves back towards the impeller wall creating an inflow condition. It can also be seen that a percentage of the flow entering the chest on the top surface is being channeled when it moves on the surface towards the exit location at the impeller wall with out being forced to the mixing zone by the impeller. J0e-0 .4 Oe-D .D De-D ,6le-D ,2le-0 ,8le-D ,41e-0 .Ole-0 ,62e-D .22e-0 ,82e-D A 2 e- D ,D2e-D ,63e-D .23e-0 .83e-D .43e-D .D3e-D .64e-D .24e-0 .3Oe-02 •• 1 * 1 1 1 1 1 1 ' • * . vi J I I / * *.' ,' • J •' 4 J 4 .' I J .41e-l2 .286-03 Figure 5.3: Velocity vectors colored by velocity magnitude (m/s) on the impeller wall, z = 0 cm, for water with (9=30L/min and N =1250 rpm 73 Chapter 5: CFD Simulation Results The flow patterns at the walls can be easily captured in the experiments and then qualitatively compared with the ones obtained from the CFD simulations at the walls and the top surface of the chest. Figures 5.3 to 5.6 illustrate the velocity vectors at the impeller wall, opposite impeller wall, side wall and top surface. l.OSe-'DD 9.97e-Dl 9.45e-Dl 8.93e-Dl 8.4 0e-0l 7.88e-Dl 7.36e-0l 6.846-01 B.32e-Dl 5.80e-0l 5.28e-Dl 4.76e-0l 4.24e-Dl 3.72e-Dl 3.2De-Dl 2.68e-Dl 2.166-01 1.63e-Dl l.lle-Dl 5.93e-D2 7.2D6-D3 Figure 5.4: Velocity vectors colored by velocity magnitude (m/s) on the wall opposite to the impeller, z = 50 cm, for water with f2=30L/min and N =1250 rpm 5.3. Non-Newtonian Simulation for Pulp Suspension (FBK,C =3.3%) m The simulations developed for water were used to verify that the discretization of the geometry was adequate to resolve the physical phenomenon taking place in the mixing chest. The main interest was to simulate the scale mixing chest with pulp suspension as the mixed fluid. However, a great difficulty arises since it is necessary to describe the rheological behaviour of pulp suspension through a mathematical model that can be used in the numerical solution of the governing transport equations outlined in chapter 4. 74 Chapter 5: CFD Simulation Results 8.05e-Dl 7.65e- D1 7.25e-Ql 6.85e- D1 6.45e-Dl 6.115ft- D1 5.66ft- D1 5.26e-Dl 4.86e-Dl 4.46e-Dl 4.D6ft-Dl 3.66ft- 01 3.26ft 01 2.86e II 2.4 7e Dl 2.D7e Dl 1.67ft Dl 1.27e Dl 8.7 De D2 4.71 e 02 7.2 D e 03 * $ JAM •• • - y N . \ V \ O K v Figure 5.5: Velocity vectors colored by velocity magnitude (m/s) on the side wall, x = 20 cm, for water with <2=30L/min and N =1250 rpm .05e-Di .65e- D1 .25e-0l .85ft- 01 .45e-0l .05e-Dl 66e-Dl 26e-D1 86e46ftD6e66e26e86e47e-0l D7e-Dl 67e-D1 ,27ft-01 70e-02 7le-02 20e-03 Figure 5.6: Velocity vectors colored by velocity magnitude (m/s) on the top surface, (y = 50cm in the scale chest, y =38cm in the C F D model) for water with (9=30L/min and N =1250 rpm 75 Chapter 5: CFD Simulation Results The main feature of pulp suspension's rheology is that they posses yield stress. They also exhibit non-linear viscoplastic properties known as shear thinning, which means that the plastic viscosity decreases as the shear rate increases. However, an accurate value for the flow behaviour index could not be obtained from the available data obtained in the yield stress measurement experiments performed by Bennington (1988), because after yielding, the suspension flow behaviour observed was irregular and showed a great dependence on the geometry of the device in which the suspension was yielded, the fibre type and suspension concentration. For a suspension at C =4%, Roberg (1997) reported n = 0.25, and Petterson m (2004) reported n = 0.08. Different n values were tried in the numerical solution, but the output response and dynamic parameters calculated from the CFD simulations showed no significant deviation from the Bingham model {n = 1.0), which had also been found by Wikstrom (2002). The data from these simulations can be found in appendix D. The Bingham model is the simplest model for a fluid possessing a yield stress, with a linear dependence of the effective viscosity with the shear rate. Thus, the pulp suspension rheology was approximated with Bingham plastic model described through the modified HerschelBulkley viscosity model function available in the F L U E N T solver and previously described in chapter 2 (n= 1.0). For very low strain rates, "rigid" materials behave like very viscous fluids with viscosity /u . As the strain rate increases and the yield stress threshold, r , is passed, the o v fluid behaviour is described by equation (5.1). (5.1) Y Where k is the consistency index, a measure of the average viscosity of the fluid and it was defined as k =0.001 kg/ms (equal to the viscosity of water); n = 1.0 is the power-law index for the Bingham Plastic model. The yield stress value was set tor = 475±62 Pa, which was y measured by Ein-Mozaffari (Ein-Mozaffari et al., 2004) for the softwood pulp suspension at 3.3% consistency used in the experiments. The yielding viscosity / / , which is defined as the 0 high viscosity of the fluid before yield, was calculated based on the experimental data collected by Bennington [1988] for semi-bleached kraft pulp tested with the HaakeRV12 viscometer ( H A A K E Inc.,Saddle Brook, NJ). The stresses vs. rotor angular displacement curves obtained by Bennington are shown in Fig. 5.7. The respective shear rate values can be obtained from the rotor displacement using the equations and calculation factors given for the Haake RV12, and the specific rotor and housing used (rotor MVIP, housing RVH1). ju = 100 0 ± 43 kg/ms was calculated for the 1.6% Cm semi-bleached kraft pulp data shown in Fig. 5.7. This value was used as a first approximation for /u in the CFD model. This parameter is 0 expected to increase as the consistency of the suspension increases but it is not possible to develop an accurate prediction for this dependence with out experimental data (suspensions with different consistency will experience different deformations under the same stresses as shown by Bennington (1988) for the behaviour of pulp suspensions after yield). Since there is no data available in the literature that describes the stress vs. deformation rate behaviour of pulp suspension before yield, /u calculated above was initially used and a sensitivity analysis 0 conducted in order to determine the effect of this parameter in the output response and dynamic parameters calculated from the CFD simulation is given in Appendix D. The results showed that increasing /u resulted in higher channeling, and higher time constants for both the bypass 0 zone and the mixing zone defined in the dynamic model developed by Ein-Mozzafari (2002). Chapter 5: CFD Simulation Results ZZ The best agreement between measured and calculated output signals and model parameters was found for /j = 300 kg/ms and ju = 100 kg/ms. A direct comparison between the CFD predicted Q 0 responses for these yielding viscosity values under two different operating conditions can be found in appendix D. Since the pulp suspension flow observed in the experiments is mostly laminar, the mass and momentum conservation equations were solved in the laminar regime. The convergence criteria used was to iterate until the scaled residuals for each transport equation were belowlxlO" . 5 Figure 5.8 shows the convergence history of the residuals for a typical simulation conducted step by step by increasing the angular velocity by 20% each time until reaching the desired value according to the solution strategies outlined in chapter 4. 0 10 20 30 Rotor Displacement, Degrees Figure 5.7: Stress vs rotor angular displacement for a 1.6% C semi4oleached kraft pulp tested with the m Haake RV12. Rotor MVIP, Housing RVH1. (Bennington, 1988) 78 Chapter 5: CFD Simulation Results 1e+00 Residuals —continuity — x-velocity y-velocity z-velocity 1e-0l ie-D2 1&-03 le-04 J 18-05 A 1e-06 5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 Iterations Figure 5.8: Convergence history of the scaled residuals for a typical step by step simulation. The simulation results are shown in the following figures for one of the experimental conditions simulated: inlet flow rate, Q = 37.1 L/min; impeller angular velocity, N = 1203 rpm. Figures 5.10, 5.11, 5.12, and 5.13 show the velocity vectors colored by velocity magnitude in the surfaces that are visible in the experimental set up. Along with the calculated vectors, a picture of each surface is shown in order to highlight the qualitative agreement between the flow field calculated and the one observed in the experiments. Figure 5.9 show the velocity vectors on a surface aligned with the blade tip which is not visible in the experiments. The empty spaces in these figures indicate that the vectors in this region have a magnitude higher than the maximum magnitude specified in the color scale range shown on the left side of the figure. 79 Chapter 5: CFD Simulation Results 56e-D2 25e-02 05.-D2 85.-D2 ,64e-02 ,44.-02 ,24.-02 ,D4e-D2 ,83e-D2 ,63e-02 ,43e-02 ,23e-02 ,02e-02 ,82e-02 .62e-D2 ,42e-02 ,2le-02 .0le-D2 .1 De-03 .07.-03 .05e-03 ,D3e-D3 .02e-03 Figure 5.9: Velocity vectors colored by velocity magnitude (m/s) on a surface aligned with the blade tip at x = 9cm. Pulp suspension at Cm = 3.3%, Q = 37.1 L/min, N=1203 rpm. The box below the inlet pipe has a higher density mesh. The mesh was adapted dynamically based on velocity gradients and normal distance to the walls. By using solution-adaptive grid refinement, cells were added where they were needed in the mesh, thus enabling the features of the flow field to be better resolved. The resulting mesh is optimal for the flow solution because the solution itself is used to determine where more cells are added. Every time the adaptation is performed, the cells with velocity gradients larger than the specified maximum gradient which was set to 0.8 (80% of the maximum velocity gradient in the domain) are splitted into eight tetrahedra. The cells selected for adaptation are always the one located near the impeller blades. To maintain accuracy, neighboring cells are not allowed to differ by more than one level of refinement. This prevents the adaptation from Chapter S: CFD Simulation Results 80 producing excessive cell volume variations and ensures that the positions of the original and refined cell centroids are similar, reducing errors in the flux evaluations. I • ' ' ' •' 0 0e-D3 73e-D3 56e-D3 3&e-03 20e-D3 02e-03 ,&4e-D3 ,67e-D3 ,49e-D3 ,3le-03 ,l3e-03 ,96e-03 7Se-03 ,6De-03 ,42e-03 .25e-D3 ,D7e-03 ,90e-04 .I2e-D4 ,34e-D4 .56e-D4 ,79e-04 ,D2e-06 -- - - • - «. ^ \ • i M i a y „ •--' ~l*~ r "- \ \ \ \ \ \ \ • S \ * '• K K ! v . : i i . I I • •• 1 1 t t i t i • • - . . • . i 1 , I - Impeller Figure 5.10: Velocity vectors colored by velocity magnitude (m/s) on the side wall surface. Pulp suspension at Cm = 3.3%, Q = 37.1 L/min, /V = 1203 rpm. This adaptation process was used to verify grid independence of the solution by demonstrating that additional grid lines near the impeller surface did not change the calculated velocity magnitude measured at a line located right in front of the impeller (z = 10cm) as illustrated in Chapter 5: CFD Simulation Results 8_1 Fig. 5.14. This location was chosen because it is the most critical location for the numerical solution since the highest velocity gradients and mass imbalance residuals are located in this area. The yellow line shows the calculated velocity with the mesh obtained through dynamic adaptation performed after 500 iterations (202,010 computational cells). The 2 nd and 3 rd adaptation lines correspond to the velocity calculated by performing one and two additional adaptations respectively, on the final dynamic adapted mesh (426,520 and 786,422 cells respectively). The velocity magnitude calculated in front of the impeller with these three different meshes did not change by more than 3%, thus the solution obtained with the dynamic adapted mesh was shown to be grid independent. 2.006-02 1.876-B2 i.78e-D2 1.69e-B2 1.6I.-I2 1.5ie-B2 1.42e-B2 1.33e-B2 1.246-D2 l.l6e-B2 1.076-02 9.78*-19 8.89e-03 3.D0e-03 r.ne-o3 6.22e-03 5.33e-03 4.45e-03 3.56e-03 2.67e-03 1.786-B3 8.90e-O4 i.02e-B6 Figure 5.11: Velocity vectors colored by velocity magnitude (m/s) on the top surface. Pulp suspension at Cm = 3.3%, £ = 37.1 L/min, W = 1203 rpm. Because the mesh used is unstructured and adapted in the impeller region, many more cells are located in this area. There is also a very small square with unstructured mesh located below the entrance pipe because the face mapping required by the structured mesh generation can not be satisfied at the point in the domain at which the mesh lines from the entrance and exit pipes meet each other. Since each cell face that intersects the plane being examined will display one velocity vector, there is a large number of velocity vectors near the impeller and in the unstructured mesh areas compared to the rest of the domain, as can be appreciated in Figs. 5.3, 5.9, 5.10 and 5.12. 12 02 D2 02 02 02 02 D2 02 02 02 Z.44Q 02 2.22e 02 2.0 De D2 1.78e 02 1.56e 02 1.33e 02 1.1 le 02 8.89e 03 e.67e OS 5.D De 4.&7e 4.44e 4.22e A.Ue 3.78e 3.56e 3.33e 3.lie 2.89e 2.67e Figure 5.12: Velocity vectors colored by velocity magnitude (m/s) on the impeller wall surface. Pulp suspension at Cm = 3.3%, Q = 37.1 L/min, JV = 1203 rpm C h a p t e r 5: 83 CFD Simulation Results S.QOe- 04 4 . 6 7 V IK 4.45ft-04 4.22e- 04 4.DD&- Dt 3.7&e-04 3.5S6-04 3.34e- 04 3.116-04 2.89ft- 04 2.67e- 04 2.45e- 04 2.23e- 04 2.0le-04 1.7&e-04 1.56ft- 04 1.34 6-04 1.12ft 04 8.97e D5 6.76e 05 4.546 05 2.32* 05 1.02ft OS <*v*SK I *,V*^ i r 1 t i f f i f f f r rr lf J/ ! 1 , ? r ? ? /,' .f r r ? j Figure 5.13: Velocity vectors colored by velocity magnitude (m/s) on the opposite impeller wall surface. Pulp suspension at Cm = 3.3%, Q = 37.1 L/min, N = 1203 rpm The pictures shown were taken under the same operating conditions as the simulations. Under these circumstances the whole suspension inside the chest was in motion but the experimental observations and the simulation results showed that even when the whole suspension is in motion, poor mixing regions exist at the bottom of the chest along the impeller opposite wall and on the surface of the chest along the impeller wall and impeller opposite wall, where pulp Chapter 5: CFD Simulation Results 84 flows significantly slower than in the bulk motion zone. In these regions pulp remains for a longer time which decreases the percentage of fully mixed volume within the chest. x coordinate (cm) Figure 5.14: Grid independence verified by comparing the velocity in front of the impeller for three different levels of mesh adaptation. 5.4. Power Input Comparison The moment vector M , about the center of the impeller wall zone with coordinates i m' Ten,> cm) x 2 = (0,12,7) is computed by F L U E N T in the manner described in chapter 4. The power input can then be calculated as P = IrrNM. In the experiments, impeller torque and speed were measured using an inductive-rotary torque transducer with an encoder disk mounted on the rotor shaft, as described in chapter 3. Thus, the calculated power input from the CFD simulations can be directly compared with the measured one. Table 5.1 shows the calculated 85 Chapter 5: CFD Simulation Results and the measured power input for different impeller speeds and input flow rates, along with a deviation percentage calculated based on the experimental measurement with equation (5.2) %Deviation •• P 1 CFD -P 1 exp xl00% (5.2) exp Table 5.1: Comparison of the power input calculated from CFD and measured experimentally Q=37.1 L/min \\ N (rpm) P (N.m/s) % Deviation 1031 274.3 232.1 18.0 1203 367.3 330.0 11.3 1456 583.0 519.9 12.1 1556 644.3 630.5 2.1 P C F D( N.m/s) Q=21.1 L/min N (rpm) PcFD(N.m) P(N.m) % Deviation 1031 274.1 232.1 18.0 1456 549.0 590.0 6.9 Q=7.9 L/min N (rpm) P C F D( N.m) P (N.m) % Deviation 1031 274.6 233.2 17.7 1456 525.3 519.9 1.0 The power input predicted by the CFD model matched the experimentally determined values reasonably well, especially at higher impeller speeds where the percentage of deviation obtained was as low as 1%. For all conditions studied, the experimentally determined values are over predicted by the C F D model; however, the shape of the power input vs. angular velocity curves is identical, as illustrated Fig. 5.15. Representative error bars are shown (Fig. 5.15) for the measured power data based on the accuracy of the torque meter (± 0.05 Nm). The discrepancy between the experimental and the CFD predicted power input may be a result of the approximated rheology model used to describe pulp suspension flow, the use of a laminarflow model in a transitional-flow regime near the impeller and/or the CFD model not capturing all of the small local recirculation zones and swirls near the blade tip. However, the Chapter 5: CFD Simulation Results 86 improvement observed in the agreement between the computed and the measured power input at higher impeller speeds suggests that at higher shear rates the non-Newtonian viscosity calculated by the rheology model more closely approaches the actual viscosity of the suspension. 720.0 640.0 560.0 480.0 400.0 °- 320.0 240.0 160.0 - 80.0 - 0.0 -, 500 , 600 - --, 700 , , , , . , , 800 900 1000 1100 1200 1300 1400 1 1500 1 , 1600 1700 N (rpm) Figure 5.15: Calculated and measured power input vs. impeller rotational speed. C =3.3%; 0 = 37.1 L/min. 5.5. Numerical Dynamic Test Results The dynamic test conducted experimentally as outlined in chapter 3 was also simulated in F L U E N T by means of tracer addition at the velocity inlet and the transient solution of the species conservation equation shown in chapter 4. Figure 5.16 show the path lines followed by the tracer according to the flow field previously calculated. Non-desirable flows such as channeling of the pulp feed and recirculation can be clearly appreciated from these particle trajectories. Chapter 5: CFD Simulation Results 87 Figure 5.16: Trajectories of massless tracer particles released from the velocity inlet. Each particle is given a different color. Initially, the simulations were performed with a single-step input signal and the input/output tracer concentrations were recorded every time step. The time step size, At, was initially set to be at least one order of magnitude smaller than the time constant estimated for the mixing system and it was further reduced to ensure that the number of iterations per time step was 2030 (FLUENT Inc., 2003). The dynamic data obtained using different step sizes for each flow rate were compared to verify temporal discretization independency. For example, at Q = 37.1 L/min the model generated dynamic response was compared for At = 1 second (obtained using the method above) and for At = 0.5 seconds. The predicted dynamic responses for both tests overlapped, indicating temporal independence. Similar results were found for tests made at Q 88 Chapter 5: CFD Simulation Results = 7.9 L/min for time step sizes of 5 and 2.5 seconds. It turns out that the sampling interval chosen based on the above procedure agreed with the sampling interval used by Ein-Mozzafari (2002) for acquiring the experimental data. Consequently all other C F D tests were made using the experimentally determined At. These single-step tests also confirmed that the simulated system had unitary gain as can be seen in Fig. 5.17 which shows the input and output signals collected from the CFD simulation for the operating conditions Q = 7.9 L/min, N= 1456 rpm. 0.11 flow time (s) Figure 5.17: Single step input signal and output signal obtained from the C F D simulation. C m - 3.3% , Q =7.9 L/min, N= 1456 rpm Since the choice of input signal has a very substantial influence on the calculation of the parameters for the dynamic model, a frequency modulated random binary signal was designed by Ein-Mozafari (2003) to provide adequate excitation for the estimation parameters of the scale mixing chest model shown in Fig. 2.5. of the dynamic This signal was written in C programming language and linked to F L U E N T through a user defined function (UDF). In order to obtain data that could be directly compared to the experimental raw data, the mass fraction concentration of the tracer was set to values analogous to the conductivity of the random binary excitation signal used in the experiments. The excitation signal UDFs written can be found in appendix C (one was needed for each flow rate). The results obtained from the simulated and experimental dynamic test conducted under identical operating conditions, are plotted together in Fig. 5.18. The blue lines correspond to the input signals and the green lines correspond to the output responses. The darker lines are the imposed input signal and the output response obtained from the CFD simulated dynamic test. The excitation signal used in the simulations did not include the noise present in the measured signals, thus the numerically obtained output is a smooth signal. Leaving aside the noise associated with the experimental data, the agreement between the simulated and the measured response is in general quite satisfactory. Figure 5.18 (a) illustrates the input/output raw data comparison for a low flow rate. In this case, the CFD predicted response agrees fairly well with the measured response particularly for the first step of the random binary input signal. In this case, a gradual increase of the conductivity at the output is followed by both the C F D and the experimental response. A dramatic decrease in the conductivity measured at the outlet can be appreciated immediately following the drop in the excitation signal due to the high degree of channeling in the system. The measured responses to the following steps of the input signal show a remarkable initial bump (sudden increase and decrease in conductivity) that is not predicted by the C F D simulated response. This bump is associated with the short-circuiting taking place within the Chapter 5: CFD Simulation Results 90 chest and was quantified by the / parameter included in the mixing chest dynamic model developed by Ein-Mozzafari (2003). These results suggest that the CFD simulation provides a dynamic response that is closer to that of a perfectly mixed chest. At higher flow rates, the CFD predicted output gives better agreement with the experimental response, as shown in Figs. 5.18 (b) and (c). The CFD simulated response is below the measured output for the lowest flow rate and over measured output for the highest flow rate. A good match is noticeable in particular for the medium flow rate case depicted in Fig. 5.18 (b). In Fig. 5.18 (c), the measured response is below the simulated one for the first step of the input signal, but after the second step the predicted and measured output conductivity values are much closer. 4.5 -I 0 1 1 1 1 1 1 1 1 1 500 1000 1500 2000 2500 3000 3500 4000 4500 flow time (•) (a) Q = 7.9 L/min, #=1031 rpm, C„, = 5.5% 1— 5000 91 Chapter 5: CFD Simulation Results 250 500 1000 750 flow 1500 1250 (b) Q = 21.1 L/min, N = 1031 rpm, C m 100 200 300 400 1750 2000 time (sec) 500 600 flow time (sec) = 3.3% 700 800 900 1000 1100 (c) Q = 37.1 L/min, N= 1456 rpm, C _ = 3.3% Figure 5.18: Input (blue)/Output (green) signals from the numerical dynamic test (dark) and the scale model chest (light) for different operating conditions Chapter 5: CFD Simulation Results 92 Despite the mentioned deviations, the simulated signal approximates the experimental response to a good extent in the majority of cases. Thus, the mixing dynamics in the scale model chest are properly predicted by the velocity field calculated in the CFD simulations. This constitutes a valuable technique for validation of the CFD model developed. The parameters of the dynamic model (Ein-Mozzafari, 2003) can be estimated using the numerical method developed by Kammer et al. (2001), using the input-output response data collected both numerically and experimentally for the agitated stock chest. The parameters calculated based on the CFD simulation can be compared to the ones measured experimentally, and a better understanding of the C F D prediction can be obtained. Estimation of the model parameters from input-output data includes two stages: the estimation of discrete-time parameters performed via least squares and estimation of the continuous-time parameters performed via Sequential Quadratic Programming (SQP) due to the non linearity of the relationship among these parameters (Kammer et al., 2001, Ein-Mozzafari, 2003). The computer program written in M A T L A B by Kammer et al. (2001) was used to calculate the dynamic parameters of the model from experimental and CFD calculated data. Table 5.3 shows the dynamic parameters calculated from the experimental and the C F D input/output data and correspond to the plots shown in Figs. 5.18 (a) and (b). It can be seen that at higher feed rate (Q) the system is prone to a higher degree of channeling (/increases as Q increases when all other factors are kept constant). All the parameters calculated from the CFD seem to be in the same range as the experimental ones, except for the time constant r, of Chapter 5: CFD Simulation Results 93 the transfer function G that represents the limited amount of mixing that can occur in the x bypass flow. Even though, r, is still considerably smaller than r , the value calculated from 2 the C F D data is about five times higher than the one obtained from the experimental data. This suggests that the C F D simulated flow provides more mixing in the bypass flow. Table 5.3: Comparison of the parameters estimated from experimental and CFD data for FBK pulp suspension at C = 3.3% at Q = 7.9 L/min and Q = 21.1 L/min for an impeller speed of 1031 rpm. m Experimental Dynamic C F D Test Test Parameters Q = 7.9 L/mih, N= 1031 rpm 0.44 0.51 f 0 0 R 35 55 Ti (sec) r, (sec) 11.4 65 T (sec) r (sec) 170 230 1150 1245 V f u l l y mixed(L) 74 92 2 2 Q = 21.1 L/min, N =1031 rpm 0.59 0.57 / 0 0 R 20 22 Ti (sec) r, (sec) 5.4 22.5 T (sec) r (sec) 74 96 400 493 V f u l l y mixed(L) 60 71 2 2 The fully mixed volume of the agitated zone can be calculated by eliminating the effect of recirculation (i.e. R = 0 for all tests). This is a reasonable assumption, because this effect was not measured in the dynamic responses of the scale model (Ein-Mozzafari, 2002). Therefore, the time constant for the agitated zone T is given by: 2 Chapter 5: CFD Simulation Results 94 V fully mixed / *•>=—— (l-f)Q ^ \ (5-3) c This allows calculation of the fully mixed volume using equation 5.3. The ideal performance for an agitated pulp stock chest is achieved when / = 0 and V fullymixeJ is equal to the total suspension volume inside the chest. Under these circumstances the dynamic behaviour of the mixing chest is defined only by a first order transfer function. The effect of the impeller speed on the dynamic parameters of the model was also examined. Table 5.4 shows the parameters calculated from CFD and experimental data for different impeller speeds under the same flow rate through the chest. The data show that channeling decreases as impeller speed increases. At higher impeller speeds the fibre/fibre network is disrupted further from the impeller leading to improved mixing by eliminating non-ideal flow such as channeling and dead zones. The parameters calculated from CFD data agreed with the observed experimental trend, but the effect of increasing the impeller speed on/(the fraction of the stock that bypasses the mixing zone) is not as appreciable as it is with the experimental test. For instance, while the experimentally calculated parameters show a 72% reduction i n / when the impeller speed is increased from 1031 rpm to 1456rpm, with Q = 7.9 L/min; the C F D calculated parameters show only a 18% reduction in/for the same circumstances. At Q =31A L/min, the experimentally measured parameters show that the onset of significant channeling can be quite sudden. At this flow rate, the percentage of channeling, f, increases from 35% to 83% when the impeller speed drops only 250 rpm. This sudden jump is not predicted in the parameters calculated from the CFD input/output data. Chapter 5: CFD Simulation Results 95 At 1456 rpm and Q = 7.9 L/min, the whole suspension inside the chest was in motion and no stagnant zones were observed. However, from the calculated parameters from both CFD and experimental data, it can be seen that there is still a considerable percentage of channeling which decreases fully mixed volume. Therefore, non ideal flows can persist even upon achieving complete motion. The calculated parameters from both the CFD and experimental input/output data presented in tables 5.2 and 5.3; show that a decrease in pulp flow rate decreases the channeling which translates in an increase in the fully mixed volume. This leads to better mixing by increasing r , the time constant of the mixing zone. 2 Dynamic tests performed by Ein-Mozzafari (2002) showed that the pulp feed and exit locations significantly affect the degree of upset attenuation and the non-ideal behaviour of the mixing chest. Based on the results previously presented, non-ideal flow decreases with increasing impeller speed and decreasing pulp flow rate. But higher impeller speed requires higher installed power and decreasing the pulp flow rate reduces the production capacity of the chest. Ein-Mozzafari (2002) showed that the dynamic performance for the scale chest could be improved by changing the feed and exit locations to that of Config.2 illustrated in Fig. 3.4. The flow patterns Figs. 5.2, 5.9 and 5.16 for Configl. show that a percentage of the pulp feed can easily channel when it moves on the surface towards the exit location at the impeller wall without being forced to the mixing zone by the impeller. Based on this, Config.2 will force flow into the mixing zone before leaving the chest. 96 Chapter 5: CFD Simulation Results Table 5.4: Comparison of the parameters estimated from experimental and C F D data for F B K pulp suspension at C = 3.3% at Q = 37.1 L/min and Q = 7.9 L/min for different impeller speeds. m Dynamic Experimental Parameters Test C F D Test Q = 37.1 L/min, N. = 1203 rpm 0.66 0.83 / R 0 0 Tj (sec) 14 13 r, (sec) 1.8 14.9 T (sec) 108 91 r 133 324.4 14 68 2 (sec) 2 V f u l l v mixed(L) (9 = 37.1 L/niin,.V= 14 56 rpm / 0.35 0.60 R 0 0 Ti (sec) 10 12 r, (sec) 3.0 14.5 T2 (sec) 25 26 167.0 366.5 67 90 r (sec) 2 V f u l l y mixed(L) Q = l » L/min, V = 125S rpm / 0.35 0.50 R 0 0 Ti (sec) 45 60 r, (sec) 13.3 75.5 r 140 295 966.1 1413.1 78 93 r 2 2 (sec) (sec) Vfullv mixed(L) 0 = 1 )< L/min, N = 1456 rpm f R 0.14 0 0 Tj (sec) 45 30 r (sec) 14 32 r 90 50 814 1221 92 102 x r 2 2 (sec) (sec) Vfullv mixed(L) 0.36 Chapter 5: CFD Simulation Results 97 The velocity vectors calculated on a plane at x = 0 cm, and the tracer particle trajectories for configuration 2 are shown in Fig. 5.19 and 5.20 respectively. It can be seen that the effect of channeling has been completely eliminated and the fully mixed volume has been improved. The input/output data raw data comparison for Config.2 atQ- 7.9 L/min and N = 1105 rpm, is illustrated in Fig. 5.21. The measured output and the CFD predicted response agree fairly well. Both signals describe a progressive increase in the conductivity that follows the changes given by the excitation signal. The C F D response is slightly above of the experimentally measured response throughout the whole test. Under these circumstances the dynamic parameters to be calculated are only T , r and R as shown in table 5.4. 2 2 The time delay and time constant calculated from the simulated data are smaller than the ones calculated from the experimental data. The time constant of the mixing zone is under predicted by 8% for the lower impeller speed, and by 11% for the higher impeller speed. Consequently, the fully mixed volume calculated with the parameters obtained from CFD data is smaller than the one calculated with the experimental parameters. Under the same operating conditions, Q = 7.9 L/min and N =1031 rpm, the time constant and time delay of the mixing zone calculated from both numerical and experimental data are much lower for Config.2 than for Config.l. Thus, the system given by Config.2 responds faster to changes in the input signal. The fully mixed volume predicted by the CFD simulation for Config.2 is in all cases smaller than the experimentally calculated one. As can be seen in Fig.5.19, even though the channeling has been eliminated, the velocity vectors show that there is a large area on the opposite side of Chapter 5: CFD Simulation Results 98 the impeller that is moving at very low velocities thus a higher percentage of the full volume is stagnant. Table 5.4: Comparison of the parameters estimated from experimental and C F D data for FBK pulp suspension at C m = 3.3% at g = 7.9 L/min and different impeller speeds for Config.2. Dynamic Experimental C F D Test Parameters Test Q = 7.9 L/min, N = 1031 rpm 0 0 f 0 0 R 10 40 T (sec) 2 r (sec) 2 mixed(L) Vfuiiy 684 623 90 82 £> = 7.9 L/min, Y = 1105 rpm / 0 0 R 0 0 T (sec) 45 10 r 708 673 93 88 2 (sec) 2 mixed(L) Vfully Q = 7.9 L/min, N= 1557 rpm 0 0 / 0 0 R 12 40 T (sec) r (sec) 694 789 2 2 mixed(L) Vfullv 91 103 The fully mixed volume calculated from the experimentally determined parameters under the same operating conditions (Q = 7.9 L/min and N =1031 rpm) is considerably higher for Config.2 than for Configl. It is important to note that at these conditions, the fully mixed volume predicted by the CFD for Config2. (Vf n i d U y m xe = 82 L) is actually smaller than the one predicted for Config.l (Vf n i d = 95 L). This is because for Config.l, the time constant of u y m xe 99 Chapter 5: CFD Simulation Results the mixing zone is over predicted and the percent of channeling is under predicted. These two factors combined together in Eqn. (5.3) yield a higher fully mixed volume. .0 D e - D 2 .85e-02 .70e-D2 .55e-02 .40e-D2 .25e-02 .106-02 .95e-02 .a0e-02 .65e-D2 .5De-02 .35e-D2 .2De-02 .D5e-02 .0 0 e - 0 2 .5De-D3 .0 0 e - 0 3 .a 0 e - D 3 .00e-D3 .5 0 e - D 3 .196-07 us/Rim?' JI &r~ . //(V.^j-tfv''. 'A'' 1 I i i rmrnn i •Vw:->*~-.w.' r Figure 5.19: Velocity vectors colored by velocity magnitude (m/s) on the plane x = 0 cm, for pulp with <2=7.9L/min and N =1031 rpm. Config. 2. The results have shown that the performance of the chest operating under Config.2 is closer to that of an ideal agitated chest. However, Config. 1 is widely used in industry due to the ease of mechanical installation and maintenance. The noise associated with the measured input signal was also propagated in the measured output response as can be appreciated in Figs. 5.17 and 5.21. However, when better mixing was Chapter 5: CFD Simulation Results 100 achieved, as it was the case for Config.2, the output response showed a better attenuation of this noise (Fig. 5.21). This disturbance however was not included in the CFD simulated input signals and thus it was not accounted for in the predicted responses. Figure 5.20: Trajectories of massless tracer particles released from the velocity inlet for Config.2. Each particle is given a different color. Chapter 5: CFD Simulation Results 101 Figure 5.21: Input (blue) /Output (green) signals from the numerical dynamic test (dark) and the scale model chest (light) for Config.2. operating at Q = 7.9 L/min and N = \ 105 rpm Chapter 6: Conclusions and Recommendations for Future Work 102 6. C O N C L U S I O N S A N D R E C O M M E N D A T I O N S F O R F U T U R E W O R K 6.1. C O N C L U S I O N S > The most encouraging conclusion drawn from this work is that C F D simulation can be effectively used to study mixing dynamics in agitated pulp stock chests. The predicted flow fields describe the observed non desirable flows such as channeling, recirculation and stagnation within the chest. The velocity vectors obtained from the CFD simulations agree qualitatively with the flow patterns observed in the experiment. Furthermore, the power input calculated in the C F D model matches the experimentally determined value reasonably well, especially at higher impeller speeds where the percentage of deviation obtained was as low as > 1%. The use of tracer distribution analysis in CFD allows the simulation of dynamic tests and collection of input/output data that can be directly compared with experimental measurements. The chest's dynamic response to a frequency modulated random binary signal that was specifically designed to provide the best excitation to the system, was obtained from the C F D simulation and compared with the measured response. The agreement of these signals was satisfactory. This not only validates the present model but provides a valuable tool for validation of the mixing dynamics predicted by similar CFD simulations. > The calculation of the parameters of the dynamic model developed by Ein Mozzafari (2002) based on both experimental and CFD data were used to quantify and compare the non ideal Chapter 6: Conclusions and Recommendations for Future Work 103 flows within the chest such as channeling and recirculation. Based on the results obtained, the CFD model significantly over predicts the amount of mixing occurring in the bypass flow. > The effect of increasing the impeller speed and flow rate was examined in the simulations and the experiments. In both cases, the degree of channeling and dead zones decreases as the impeller speed increases and the flow rate decreases. The on-set of significant non-ideal behaviour can be quite sudden. Under these circumstances, the dynamics of the mixing chest can be improved by reducing the pulp flow rate. > The effect of pulp feed and exit locations, was also simulated using CFD. The degree of channeling is eliminated once the exit port is located as shown in configuration 2. Therefore, the total mixed volume increases and the dynamic performance of the chest is closer to that of the ideal agitated chest whose behaviour is defined only by a first transfer function. > A validated CFD model of agitated pulp stock chests provides detailed information of flow profiles at different locations which allows the conditions leading to flow pathologies and poor dynamics to be identified and removed. 104 Chapter 6: Conclusions and Recommendations for Future Work 6.2. R E C O M M E N D A T I O N S F O R F U T U R E W O R K > The Bingham model was used as an approximation of the rheological behaviour of pulp suspension. However, measurements and C F D predictions can be attributed to this approximation of the suspension rheology. much of the deviations observed between experimental A more accurate model of the suspension rheology would likely significantly improve the accuracy of the CFD model developed. > The construction of a structured mesh for the whole computational domain including the impeller zone can provide great improvement in the numerical solutions in terms of accuracy, stability and computational expenses. This requires a tremendous effort specially when dealing with complicated geometries. Commercial software packages such as MIXSIM, which contain structured mesh files for most of the industrially relevant mixers, can be purchased in order to overcome this difficulty. > The CFD model can also be improved by the implementation of the unsteady solution of the governing equations using the sliding mesh technique available in F L U E N T . This will be more computational expensive but will allow us to determine the quality of the steady state approximation used in the MRF technique. > The boundary conditions used in the model can be refined by including the modeling of the free surface. This can be done by including multiphase' modeling using the Volume of Chapter 6: Conclusions and Recommendations for Future Work Fluid tool in F L U E N T 105 in order to include the interface (suspension/air) in the computational domain. > More rigorous validation techniques for the CFD simulations are still required. This can be accomplished by comparison of the model predicted velocity vectors with experimental measurements obtained from experimental techniques such as Particle Image Velocimetry (PIV) and Ultrasonic Velocity Profiling (UVP). The implementation of PIV will require the use of transparent model fluids with rheology similar to that of pulp suspension. > Once validated, the C F D model will allow us to conduct reliable "virtual" studies on the effect of all important variables interrelationships. on mixing effectiveness and to model their Hence, design methods and trouble shooting of existing agitated pulp stock chests can be significantly improved. 106 Nomenclature N O M E N C L A T U R E c impeller C factor, dimensionless C fibre mass concentration, % D impeller diameter, m D housing diameter, m f percentage of channeling, fraction k consistency factor (kg/ms) g gravitational acceleration, m/s Gj transfer function for channeling zone G transfer function for mixing zone L chest length, m LMo level momentum, m/s m suspension property (Eqn. 2.12), m/s M torque, N.m Mo impeller momentum number, m /s N impeller rotational speed, s" M T 2 2 2 4 2 1 N Froude number, dimensionless N power number, dimensionless Fr P Np urb power number for fully turbulent flow (dimensionless) NQ impeller flow number, dimensionless N Reynolds number, dimensionless P power, W tt Re 107 Nomenclature ~1 Q pulp flow rate through the chest, m/s Qimpeller impeller pumping rate, m/s R percentage of recirculation, fraction Re Reynolds number, dimensionless t time, s ts sampling time, s T tank diameter, m T, time delay for the channeling zone, s T time delay for the mixing zone, s V fluid velocity leaving the impeller, m/s Vbulk bulk velocity, m/s V suspension volume inside the chest, m Vfully mixed fully mixed volume, m 2 w z 3 chest width, m suspension height, m Greek letters Q angular velocity, rpm f shear rate, s" Yavg average shear rate in the mixing vessel, s power consumption per unit mass, W/kg sj power dissipation at point of fluidization, W/m 1 avg rj p plastic viscosity, N.s/m Nomenclature . fj, fluid viscosity, kg/m.s jUa apparent viscosity, kg/m.s ju yielding viscosity, kg/m.s p fluid density, kg/m T shear stress, N / m n time constant for the channeling zone, s r time constant for the mixing zone, s r yield stress, N/m OJ frequency, rad/s 0 2 y 2 Abbreviations CFD computational fluid dynamics FBK fully bleached kraft MRF moving reference frame PRESTO pressure staggering option RTD residence time distribution SIMPLE semi-implicit method for pressure linked equations VOF volume of fluid model 108 109 Bibliography BIBLIOGRAPHY Attwood, D. and J.D. Gibbon, "The Agitation of Paper stock", Paper Technology, 4(l):54-62 (1963). Bakker, A . and J.B. Fasano, "A Computational Study of the Flow Pattern in an Industrial Paper Pulp Stock Chest with a Side-Entering Impeller", AIChE Symposium Series, 89(293): 118-124 (1993). Bennington, C.P.J., "Mixing Pulp Suspensions", PhD Thesis, University of British Columbia (1988). Bennington, C.P.J., "Mixers and Mixing" in "Pulp Bleaching: Principles and Practice", C W . Dence and R.W. Reeve (eds), TAPPI Press, Atlanta, G A (1996). Bennington, C.P.J., R.J. Kerekes and J.R. Grace, "The Yield Stress of Fibre Suspensions", The Canadian Journal of Chemical Engineering, 68:748-757 (1990). Bennington, C.P.J., R.J. Kerekes and J.R. Grace, "Motion of Pulp Fibre Suspensions in Rotary Devices", The Canadian Journal of Chemical Engineering, 69:251-258 (1991). Bennington, C.P.J, and R.J. Kerekes, "Power Requirements for Pulp Suspension Fluidization", TAPPI Journal, 79(2):253-258 (1996). Bertrand F. And P. Tanguy, "A New Perspective for the Mixing of Yield Stress Fluids with Anchor Impellers", Journal of Chemical Engineering of Japan, 29(1): 51-58 (1996) Bialkowski, W . L . , "Process Control Related Variability and the Link Performance", Control Systems Conference, Helsinki, Finland, 220-248 (1990). to End-Use no Bibliography Bialkowski, W.L., "Newsprint Variability and Its Impact on Competitive Position", Pulp and Paper Canada, 93(11):T299-T306 (1992). Blasinski, H . and E . Rzyski, "Mixing of Non-Newtonian Liquids. Power Consumption for Fibrous Suspensions", Inz. Chem., 2(1):169-182 (1972). Brown, J.R., "The Dynamic Behavior of Paper Mill Stock Chests", Southern Pulp and Paper Manufacturer, 31(6):103-112 (1968). Brucato A . , M . Ciofalo, F. Grisfi and G. Micale, "Numerical Prediction of Flow Fields in Baffled Stirred Vessels: A comparison of Alternative Modeling Approaches". Chemical Egineering Science. Vol.53 (21):3653-3684 (1998) Brucato, A., M.Ciofalo, F. Grisfi and L . Rizzuti, "Application of a Numerical Fluid Dynamics Software to Stirred Tanks Modeling". Super Computing Tools for Science and Engineering, D. Laforenza and R. Perego, eds. Pp 413-419. Franco Angeli, Milano. (1989) Calderbank, P.H. and M . B . Moo-Young, "The Prediction of Power Consumption in The Agitation of Non-Newtonian Fluids", Trans. Instn. Chem. Engrs., 37:26-33 (1959). Dickey, D.S., "AIChE Equipment Testing Procedures: Mixing Equipment (Impeller Type)", AIChE, New York, 3 edition (2001). rd Ein-Mozzafari, "Macroscale Mixing and Dynamic Behaviour of Agitated Pulp Stock Chests", PhD Thesis, University of British Columbia (2002). Ein-Mozaffari, F., G.A. Dumont and C.P.J. Bennington, "Performance and Design of Agitated Pulp Stock Chests", P A C W E S T Conference, Whistler, B.C., Paper No. 3B6, 1-8 (2001). Ill Bibliography Ein-Mozaffari, F., L . C . Kammer, G.A. Dumont and C.P.J. Bennington, "Dynamic Modeling of Agitated Pulp Stock Chests", Control Systems Conference, Stockholm, Sweden, 194-199 (2002). Ein-Mozaffari, F., L . C . Kammer, G.A. Dumont and C.P.J. Bennington, "Design Criteria for Dynamic Behavior of Agitated Pulp Stock Chests", PACWEST Conference, Jasper, A B , Paper No. 2B2, 1-9 (2002). Ein-Mozaffari, F., C.P.J. Bennington and G.A. Dumont. "Performance and Design of Agitated Pulp Stock Chests". Appita Journal 56(2):127-133(2003) Ein-Mozaffari, F., C.P.J. Bennington and G.A. Dumont. "Dynamic Mixing in Agitated Industrial Pulp Chests". Pulp and Paper Canada, 105(5):T118-T122 (2004). Ein-Mozaffari, F., C.P.J. Bennington and G.A. Dumont. "Suspension Yield Stress and the Dynamic Response of Agitated Pulp Chests". Chemical Engineering Science (accepted) 2004. F L U E N T Inc, Fluent 6.1 User's Guide Documentation. F L U E N T Inc. (2003). Fox, E.A. and V . E . Gex, "Single-phase Blending of Liquids", AIChE Journal, 2(4):539-544 (1956). Gibbon, J.D. and D. Attwood, "Prediction of Power Requirements in The Agitation of Fibre Suspensions", Trans. Inst. Chem. Engrs., 40:75-82 (1962). Gosman, A . D . , R.I. Issa, C. Lekakou, M . K . Looney and S. Politis, "Multidimensional Modeling of Turbulent Two Phase Flow in Stirred Vessels". AIChE Journal, 38:1946- 1956(1992). Gullichsen, J., E . Harkonen, "Medium Consistency Technology", TAPPI Journal, 64(6):69-72 (1981). 112 Bibliography Harnby, N . , M . F . Edwards and A . W . Nienow, "Mixing in the Process Industries", Butterworths, London (1985). Head, V.P. and R.E. Durst, "Stock slurry Hydraulics", TAPPI Journal, 40(12):931-936 (1957). Ferzieger J. L . and M . Peric. "Computational Methods for Fluid Dynamics ". Springer-Verlag, Heidelberg, 1996. Fokema, M . D . , S.M. Kresta and P.E. Wood, "Importance of Using the Correct Impeller Boundary Conditions for C F D Simmulations of Stirred Tanks". Canadian Journal of Chemical Engineering, (72):177-183 (1994) Hutchings, B.J., R.J. Weetman and B.R. Patel, " Computation of Flow Fields in Mixing Tanks with Experimental Verification", Paper TN-4810, A S M E Annual Meeting, San Francisco, U S A (1989). Kammer, L . C . , F. Ein-Mozaffari, G.A. Dumont and C.P.J. Bennington, "Identification of Channeling and Recirculation Parameters of Agitated Pulp Stock Chests", AIChE Conference, Reno, Nevada, Paper No. 328g, 1-12 (2001). Kerekes, R.J., "Pulp Flocculation in Decaying Turbulence: A Literature Review" J. Pulp and Paper Sci., 19(3):TR86-TR91 (1983). Levenspiel, O., "Chemical Reaction Engineering", John Wiley and Sons, 3 edition (1998). rd Luo, J.Y., A . D . Gosman, R.I. Issa, J.C. Middleton and M . K . Fitzgerald, "Full Flow Field Computation of Mixing in Baffled Stirred Vessels". Transactions of the Institute of Chemical Engineers. (71.A):342-344 (1993). 113 Bibliography Luo, J.Y., A . D . Gosman and R.I. Issa, "Prediction of Impeller Induced Flows in Mixing Vessels Using Multiple Frames of Reference". Institute of Chemical Engineers Symposium Series. (136):549-556 (1994). Macosko, C.W., "Rheolgy, Principles, measurements and applications". V C H Publishers, Inc. (1994) McDonough, R.J., "Mixing for the Process Industries", Van Nostrand Reinhold, New York (1992). Metzner, A.B. and R.E. Otto, "Agitation of Non-Newtonian Fluids", AIChE Journal, 3(1):3-10 (1957). Murthy, J.Y., S.R. Mathur and D. Choudhury, "CFD Simmulation of Flows in Strirred Tank Reactors Using a Sliding Mesh Technique", Institute of Chemical Engineers Symposium Series, (136):341-348 (1994) Oldshue, J.Y., "Agitation in Groundwood Blending", Pulp and Paper Magazine of Canada, 72(12):T375-T380 (1971). Oldshue, J.Y., "Fluid Mixing Technology", McGraw-Hill, New York (1983). Oldshue, J.Y., A.T. Gretton, "Performance and Design of Paper Stock Mixers", TAPPI Journal, 39(6):378-390 (1956). Posarac, D and P. Watkinson, " Mixing of a Lignin-Based slurry Fuel. Canadian Journal of Chemical Engineering, (78): 265-270 (2000) Ranade, V . , J.B. Joshi and A . G . Marathe, "Flow Generated by Pitched Blade Turbines II: Simulation Using k-s Model". Chemical Engineering Communications. (81): 197-224 (1989) Bibliography 114 Ranade, V . and J.B. Joshi, "Flow generated by a Disc Turbine: Part II. Mathematical Modeling and Comparison with Experimental Data". Trans. IChemE 68(A), 34-50 (1990) Reed, C.S., "Selecting The Right Equipment for Agitation and Blending, Part 1", TAPPI Journal, 78(6):252-254 (1995). Reed, C.S., "Selecting The Right Equipment for Agitation and Blending, Part 2", TAPPI Journal, 78(7):241-244 (1995) Reynolds, E . , J.D. Gibbon and D. Attwood, "Smoothing Quality Variations in Storage Chests Holding Paper Stock", Trans. Instn. Chem. Engrs., 42:T13-T21 (1964). Rushton, J.H., E.W. Costich and H.J. Everett, "Power Characteristics of Mixing Impellers", Chem. Eng. Progress, Part I: 46(8):395-404 (1950), Part II: 46(9):467-476 (1950). Rieger, F and V . Novak, "Power Consumption Scale-Up in Agitating Non-Newtonian Fluids', Transactions of the Institution of Chemical Engineers, (51): 105-111 (1973) Roberg K . E . , "Numerical Simulations of the Flow of Fibre Suspension in a Side-Entering Mixer", Recents Progres en Genie des Precedes. Vol. 11 (51):203-210. Paris (1997) Silvester, R.S., "Mixing of Non-Newtonian Media: A Technical Review", BHRA Fluid Engineering, England (1985). S. V. Patankar. "Numerical Heat Transfer and Fluid Flow". Hemisphere, Washington, D.C., 1980. Sweeney, E.T., "An Introduction and Literature Guide to Mixing". B H R A Fluid Engineering Series, vol.5, England (1978) 115 Bibliography Takeda H . , K . Narasaki, H . Kitajima, S. Sudoh, M . Onofusa and S. Iguchi, " Numerical Simmulation of Mixing Flows in Agitated Vessels with Impellers and Baffles", Computational Fluids (22):223-228 (1993). Uhl, V.W. and J.B. Gray, "Mixing, Theory and Practice", Academic Press, New York (1966). Walker, O.J. and A . Cholette, "Determination of the Optimum Size and Efficiency of Stock Chests. Part I: The Ideal Chest", Pulp and Paper Magazine of Canada, 59:113-117 (1958). Wahren, D., " Fibre Network Structures in Paper Making Operations", Conference Paper Science and Technology, The Cutting Edge, Institute of Paper Chemistry, Appleton, WI, 112132(14980) Wikstrom, T. and A . Rasmuson, "The Agitation of Pulp Suspension With a Jet Nozzle Agitator", Nordic Pulp and Paper Research Journal, 13(2):88-94 (1998). Yackel, D . C . , "Agitation: Theory, Mechanics and Application", in "Pulp and Paper Manufacture", Vol. 6: "Stock Preparation", R.W. Hagemeyer and D.W. Manson, TAPPI Press, Atlanta, G A (1983). Yackel, D.C., "Pulp and Paper Agitation: The History, Mechanics and Process" TAPPI Press, Atlanta, G A (1990). Yackel, D.C. and L . H . Mahony, "Agitation and Mixing in the Pulp and Paper Industry", TAPPI Engineering Conference, Houston, Texas, 207-216 (1976). Yakhot V and L . M Smith, "The renormalization group, the s expansion and derivation of turbulence models" . J. Sci. Comput. (7):35-61 (1992) 116 Appendix A : Upwind Schemes and temporal discretization APPENDIX A UPWIND S C H E M E S A N D T E M P O R A L D I S C R E T I Z A T I O N A . l . FIRST O R D E R UPWIND S C H E M E When first order accuracy is desired, quantities at cell faces are determined by assuming that the cell center values of any field variable represent a cell-average value and hold throughout the entire cell; the face quantities are identical to the cell quantities. Thus when first order upwinding is selected, the face value Of is set equal to the center value of O in the upstream cell. A.2. S E C O N D O R D E R UPWIND S C H E M E When second order accuracy is desired, quantities at cell faces are computed using a multidimensional linear reconstruction approach (Patankar, 1980). In this approach, higher order accuracy is achieved at cell faces through a Taylor series expansion of the cell centered solution about the cell centroid. Thus when second order upwinding is selected, the face value is computed using the following expression: O . = 0 +VO-AJ (A.1) Where ® and V O are the cell centered value and its gradient in the upstream cell, and AJ is the displacement vector from the upstream cell centroid to the face centroid. This formulation requires the determination of the gradient V O in each cell. This gradient is computed using the divergence theorem, which in discrete form is written as N faces v o =- v y4 - <b A 1 f (A.2) 117 Appendix A: Upwind Schemes and temporal discretization Here the face values <b are computed by averaging O from the two cells adjacent to the face. f The gradient V ® is limited so that no new maxima or minima are introduced. A.3. C E N T R A L D I F F E R E N C I N G S C H E M E The central differencing scheme calculates the face for a variable <D f f,CD as follows: 0 ~ (A.3) Where the indices 0 and 1 refer to the cells that share face / , V ® , . and V O 0 reconstructed gradients at cells 0 and 1, respectively, and r r l are the is the vector directed from the cell centroid toward the face centroid (Patankar, 1980). Central differencing schemes can produce unbounded solutions and non physical wiggles, which can lead to stability problems for the numerical procedure. (FLUENT Inc, 2003). These stability problems can be avoided if a deferred approach is used for the central differencing scheme. In this approach the face value is calculated as follows: (A.4) implicit where UP stands for upwind. explicit The upwind part is treated implicitly while the difference between the central difference and the upwind values is treated explicitly. Provided that the numerical solution converges, this approach leads to pure second-order differencing. Appendix A: Upwind Schemes and temporal discretization 118 A.4. TEMPORAL DISCRETIZATION A generic expression for the time evolution of a variable O is given by: where the function F incorporates any spatial discretization. If the time derivative is discretized using backwards differences, the first order accurate temporal discretization is given by O n + l -O" At = F(<D) (A.6) and the second order discretization is given by 30" + 1 -4<D"+® 2At F(0) (A.7) where O = a scalar quantity; n +1 = value at the next time level, t + At; n = value at the current time level, t; n-l = value at the previous time level, t-At. discretized, a choice remains for evaluating F ( O ) . Once the time derivative has been The method chosen for the simulation of transient chemical mixing in the stock chest was the implicit time integration (FLUENT Inc, 2003). When implicit time integration is chosen, F ( O ) is evaluated at the future time level as follows: _ ^ F(o ) n+l = (A.8) This is referred to as "implicit" integration since O" ' in a given cell is related to O 4 neighboring cells through F(o"): +l n + 1 in 119 Appendix A: Upwind Schemes and temporal discretization (D" =0" + ArF(0" ) +1 (A.9) +l This implicit equation can be solved iteratively by initializing O' to O and iterating the n equation <D'=<D"+A/F(<D') (A. 10) for the first order implicit formulation used in the simulation presented in this thesis. The advantage of the fully implicit scheme is that it is unconditionally stable with respect to time step size. Appendix B: k-e RNG Turbulence Model Description 120 APPENDIX B T U R B U L E N C E M O D E L I N G : R N G k-e Model The R N G k-e turbulence model is a semi-empirical model based on the transport equations for the turbulence kinetic energy (k) and its dissipation rate (e). It is derived from the instantaneous Navier-Stokes equations, using a mathematical technique called "renormalization group" (RNG) methods (Yakhot et al., 1992) The transport equations for the R N G k-e model are: dk + G +G -p -Y S kMeff a k b £ M+ (BA) k and dt ' K ' + -£ (peu ) r l dx/' " d£ ) = -l-{a |f| eMeff dx v £ , ^ „ 4 ^ s „ + C G) - C 3e b £ u 2 p ~ R +S £ E (B.2) J J J ^ + C„ | ( G 8X In these equations, Gk represents the generation of turbulence kinetic energy due to the mean velocity gradients, calculated as — du G = -p u,Uj -^du.. j i L k (B.3) Gb is the generation of turbulence kinetic energy due to buoyancy, in this case Gt, = 0. YM represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate, this effect is only important for high Mach number, compressible flows; thus YM =0 in our case. The quantities a and a are the inverse effective Prandtl numbers for k and k e e, respectively. Sk and S are the source terms for k and e respectively. E Appendix B: k-s RNG Turbulence Model Description 121 The inverse Prandtl numbers, are computed using the following formula derived analytically by the R N G theory (Yakhot et al., 1992) a-1.3929 0.6321 a + 2.3929 0.3679 _ Mmol a + 2.3929 a -1.3929 where a = 1.0 In the high Reynolds number limit y p jp 0 mol (B.4) Meff 0 0 eff « : 1 j , a =a k e «1.393 The scale elimination procedure in R N G theory results in a differential equation for turbulent viscosity: pk 2 d \^i P£ where v = p jp eff ^dv = 1.72 j (B.5) v -1 + C , ; C « 1 0 0 . Equation (B.5) is integrated to obtain an accurate description on v how the effective turbulent transport varies with the effective Reynolds number (or eddy scale), allowing the model to better handle low Reynolds number and near wall flows. In the high Reynolds number limit, Equation (B.5) gives k 2 M, = pC (B.6) M with C = 0.0845, derived using R N G theory. M Turbulence in general, is affected by rotation or swirl in the mean flow. The R N G model in F L U E N T accounts for the effects of swirl or rotation by properly modifying the turbulent viscosity (FLUENT Inc, 2003). The modification takes the following functional form: f Mi - P,of a . v > Q > - (B.7) Appendix B: k-e RNG Turbulence Model Description 122 where ju is the value of turbulent viscosity calculated without the swirl modification using l0 Q is a characteristic swirl number evaluated within the solver, and a equation (B.5). is a s swirl constant that assumes different values depending on whether the flow is swirl-dominated or only mildly swirling. The main difference between the R N G and standard k-e model lies in the additional term in the e equation given by R . ^ ' O - ^ f l i+/V k where rj = S /s,rf k 0 ( B . 8 ) = 4.38,/? = 0.012. The effects of this term in the R N G e equation can be seen more clearly by rearranging equation (B.2) using equation (B.7) and merging the third and fourth terms on the right-hand side. The resulting e equation can be written as d 6r / ( ) p£ v d , x ' + ^-( £ ) P dx^ Ui x " d r =^ - dx a, - \ ds M 9 T V ^D X J 2 +C, y(G C G )-C; .p^E k+ ie h (B.9) t J where C* is given by 2e . c £^Uhpl .,o) (B In regions where rj < rj , the R term makes a positive contribution, and C\ becomes larger than 0 C. le c As a result, for weakly to moderately strained flows, the R N G model tends to give results largely comparable to the standard k-e model. In regions of large strain rate (r/>r/ ), however, the R terms makes a negative contribution, making the value of C \ less than C . The smaller destruction of e augments e, reducing k, 0 2 2e Appendix B: k-e RNG Turbulence Model Description 123 and, eventually, the effective viscosity. Therefore, in rapidly strained flows, the RNG models yields a lower turbulent viscosity and, thus it is more responsive to rapid strain effects. The model constants C^and C .in equation B.2. have values derived analytically by the RNG 2t theory (Yakhot et al., 1992). These values, used by default in FLUENT, are C =1.42, u C =1.68 2B (B.ll) 124 Appendix C : Excitation signal UDF's APPENDIX C E X C I T A T I O N S I G N A L UDF'S The excitation signals written in C programming language and linked to F L U E N T as User Defined Functions (UDF) are presented here for the different flow rates simulated. The conductivity maximum and minimum values were varied according to the experimental data. C . l . Q = 7.9 L/min #include <udf.h" DEFINE_PROFILE (tracer_mf_prof i l e , t, 1) { r e a l flow_time face_t = CURRENTJTIME; f; begin_f_loop(f,t) { i f (flow time > 4 665 F_PROFILE(f, e Ise i f else i f e l s e i f else i f else i f else i f else i f else i f (flow time t, time time t, time time F time time time 0 ,. 0 5 2 ; = 0 .. 0 8 8 ; i) = 0 .. 0 5 2 ; t, i) = 0 ,. 0 8 8 ; = 0 ,. 0 5 2 ; = 0 ., 0 8 8 ; = 0 ., 0 5 2 ; > 2375) t, > i) 1940) t, i) > 8.15) F_PROFILE(f, (flow i) t, PROFILE(f, (flow = > 3055) F_PROFILE(f, (flow i) t, F_PROFILE(f, (flow 0 ,. 0 8 8 ; > 3345) F_PROFILE[f, (flow = > 4220) F_PROFILE(f, (flow i) > 4540) F_PROFILE(f, (flow ) t, > i) 150) F_PROFILE(f, t, i) = 0 ,, 0 8 3 ; F_PROFILE(f, t, i) = 0 ,, 0 5 2 ; else } end_f_loop(f, } t) Appendix C : Excitation signal UDF's 125 C.2. (2 = 21.1 L/min #include <udf.h>. _ DEFINE_PROFILE(tracer_rof_profile, t, _ , i) { r e a l flow_tiitie face_t = CURRENT_TIME; f; becfin_f_loop (f, t) { i f (flow_tinie > 1864 ) F_PROFILE(f, else i f else i f else i f else i f else i f else i f else i f else i f (flow_time t , i) =0 . 6 7 ; >1814) F_PROFILE(f, (flow_time t , i) =0 . 4 6 ; >1686) F_PR0FILE(f, t , (flow_time i) =0.67; > 1 5 3 6) F_PROFILE(f, (flow_time t , i) =0.46; > 1220) F_PROFILE(f, (flow_time t , i) =0.67; > 948) F_PROFILE(f, (flow_time t , i) =0 . 4 6 ; > 774) F_PR0FILE(f, t , (flow_time i) =0 . 6 7 ; > 324) F_PR0FILE(f, t , (flow_time i) =0 . 4 6 ; > 58) F_PR0FILE(f, t , i) =0 . 6 7 ; F_PROFILE(f, t , i) =0 . 4 6 ; else } end f loop(f,t) C.3. Q = 37.1 L/min j& i n c l u d e <udf.h> _ DEFINE_PR0FILE(t r acer_ntf _ p r o f i l e , t , i) { r e a l flow_time face_t = CURRENT_TIHE; f; becfin_f_loop (f,t) I i f (f low_tirne > 932 ) F_PROFILE(f, else i f (flou_time F t , i): = 0.275; > 907) PROFILE(f, t , i) = 0.165; Appendix C : Excitation signal UDF's else if e lse if else if e lse if else if else if e lse if 126 843) t, 1) 768) t, i) ( f l o w t i m e > 610) i) F_PROFILE(f, t, ( f l o w t i m e > 474) F_PROFILE(f, t, i) ( f l o w t i m e > 387) i) F_PROFILE(f, t, ( f l o w t i m e > 162) F_PROFILE(f, t, i) ( f l o w t i m e > 29) (flow time > F_PROFILE(f, (flow time > F_PROFILE(f, = 0 ,,2 7 5 ; = 0 ,, 1 6 5 ; = 0 ., 2 7 5 ; = 0 ., 1 6 5 ; = 0 ., 2 7 5 ; = 0 ., 1 6 5 ; F_PROFILE(f, t, i) = 0 ,, 2 7 5 ; F_PROFILE(f, t, i) = 0 ,. 1 6 5 ; else 124 Appendix C: Excitation signal UDF's APPENDIX C EXCITATION SIGNAL UDF'S T h e excitation signals written i n C p r o g r a m m i n g language a n d l i n k e d to F L U E N T as U s e r D e f i n e d F u n c t i o n s ( U D F ) are presented here f o r the different f l o w rates simulated. T h e c o n d u c t i v i t y m a x i m u m a n d m i n i m u m values were v a r i e d a c c o r d i n g to the e x p e r i m e n t a l data. C . l . Q = 7.9 L/min #include .urit.h DEFINE_PROFILE(tracer_mf_prof i l e , t , i) { r e a l flow_time face_t = CURRENT_TIHE; f; begin_f_loop(f,t) < i f [flow time > 4665 else i f else i f else i f else i f else i f else i f else i f else i f (flow time time time t, time t, time t, time t, t, > F_PROFILE(f, [flow time F_PROFILE (flow 0 052 i) = 0 088 i) = 0 052 i) = 0 088 i) = 0 052 = 0 088 = 0 052 1940) t, i) > 815) (f:,. t , time = > 2 3 75) F_PROFILE(f, (flow i) > 3055) F_PROFILE(f, (flow 0 088 > 3345) F_PROFILE(f, (flow = > 4220) F_PROFILE(f, (flow i) > 4540) F_PROFILE(f, (flow ) t, F_PROFILE(f, > i) 150) F_PROFILE(f, t, i) = 0 088 F_PROFILE(f, t, i.) = 0 052 else } end_f_loop(f,t) } 125 A p p e n d i x C: Excitation sisnal UDF's C.2. 0 = 21.1 L/min #inciude <udf.h> DEFINE_PROFILE(tracer_mf_prof i l e , t, i) { r e a l flow_time face_t = CURRENT_TIHE; f; beo;in_f_loop (f, t) { i f (flow_time > 13 64 ) F_PROFILE(f, else i f else i f else i f else i f (flow t, time t, F_PROFILE(f, (flow time t, time t, time (flow i f t, time i f else i f else i f (flow t, t, i ) == 0 ., 4 6 ; i) : = 0 ., 6 7 ; i) : = 0 ., 4 6 ; i) ••= 0 ., 6 7 ; t, i ) == 0 ., 4 6 ; > 58) time F_^PROFILE 0 ., 6 7 ; > 324) time F_PROFILE(f, (flow i ) == > 774) t ime F_PROFILE(f, (flow 0 ., 4 6 ; > 948) F_PROFILE(f, else i ) == 1220) > F_PROFILE(f, e lse 0 ., 6 7 ; 1536) > F_PROFILE(f, (flow == 1686) > F_PROFILE(f, (flow i) 1814) > (f, = 0 ,, 6 7 ; t, i) t, i ) == 0 ., 4 6 ; i) = 0.275; = 0.165; ; else F_PROFILE(f, } end_f_loop(f,t) } C.3. Q = 37.1 L/min ^include' <udf.h>___ , \ DEFINE_PROFILE(tracer_mf_profile, r e a l flow_time face_t = t, i) CURRENT_TIHE; f; begin_f_loop(f, t) { i f (flow_time > 93 2 ) F_PROFILE(f, else i f (flow_time F t, > 907) PROFILE(f, t, i) 126 Appendix C : Excitation signal UDF's else if else if else if e lse if 843) t, i) 768) t, i) ( f l o w t i m e > 610) F_PROFILE(f, t, i) ( f l o w t i m e > 474) (flow time > F_PROFILE(f, (flow time > F_PROFILE(f, = 0 ., 2 7 5 ; = 0 ., 1 6 5 ; = 0 ., 2 7 5 ; i) PROFILE(f, t, ( f l o w t i m e > 387) i) F_PROFILE(f, t, ( f l o w t i m e > 162) = 0 ,. 1 6 5 ; = 0 .. 2 7 5 ; t, ( f l o w t i m e > 29) F_PROFILE(f, t, i) = 0 ,. 1 6 5 ; i) = 0 ,. 2 7 5 ; i) = 0 ,. 1 6 5 ; F else if else if else if F_PROFILE(f, else F_PROFILE(f, t, 127 Appendix D: Sensitivity analysis data APPENDIX D SENSITIVITY ANALYSIS TABLES D.l. Dynamic data for the 95% confidence interval of the yield stress measurement at Q =7.9 L/min, TV = 1031 rpm and ju = 100 kg/m.s 0 Table D . l : Comparison of the parameters estimated from experimental and C F D data for F B K pulp suspension C = 3.3% ± 95% confidence interval, at Q = 7.9 L/min and 7V = 1031 rpm. m Experimental test r,=337 r =450 r =950 (Pa) (Pa) (Pa) / 0.51 0.46 0.44 0.49 R 0 0 0 0 Ti (sec) 55 30 35 40 r, (sec) 11 84 65 112 Ti (sec) 170 300 230 350 r (sec) 1150 1297 1245 1292 Vfuify Mixed (L) 74 92 92 87 2 y y D.2. Effect of the yielding viscosity on the CFD predicted output response and parameters calculation Table D.2 shows the parameters of the dynamic model calculated from the input/output data collected from the C F D simulations at different values of the yielding viscosity at Q =7.9 L/min, N = 1031 rpm and r = 475 Pa. From this results, the values of // that provide better y 0 Appendix D: Sensitivity analysis data i^o" agreement between the parameters calculated from experimental and numerical data are ^ =100 kg/ms and ju =300 kg/ms. The differences in the output responses predicted by the 0 Q CFD simulations for these two yielding viscosity values, can be appreciated in Figs D . l and D.2 for two different operating conditions. Table D.2: Comparison of the parameters estimated from experimental and C F D data for F B K pulp suspension C m = 3.3% , at Q = 7.9 L/min and TV = 1031 rpm, for different yielding viscosity values Experimental test =10 Mo = M (kg/m.s) f 0.51 R ju = 1000 ju = 100 Mo = (kg/m.s) (kg/m.s) (kg/m.s) (kg/m.s) 0.23 0.29 0.44 0.55 0.67 0 0 0 0 0 0 T] (sec) 55 20 30 35 60 80 r, (sec) 11 42 29 65 111 253 T (Sec) 170 120 150 230 345 640 r (sec) 1150 1012 1192 1245 1296 1482 Vfully Mixed (L) 74 103 111 92 84 64 2 2 Mo 0 300 0 Appendix D: Sensitivity analysis data 4.5 -I 1 0 500 1 1000 i 1500 129 i i 1 1 i 1 1 1 2000 2500 3000 3500 4000 4500 5000 5500 time (s) Figure D. 1: Input (blue)/Output (green) signals from the numerical dynamic test (dark) and the scale model chest (light) for (9=7.9 L/min, # = 1031 rpm, C = 3.3% . Comparison of the C F D calculated m response for different yielding viscosity values: ju = 100 kg/ms (green), / / = 300 kg/ms (orange) 0 400 0 500 600 flow time (sec) 700 800 900 1000 1100 Figure D.2: Input (blue)/Output (green) signals from the numerical dynamic test (dark) and the scale model chest (light) for (2=37.1 L/min, /V = 1456 rpm, C m - 3.3% . Comparison of the C F D calculated response for different yielding viscosity values: ju0= 100 kg/ms (green), /u0 = 300 kg/ms (orange) 130 Appendix D: Sensitivity analysis data D.3. Effect of theflowindex on the parameters calculation from CFD input/output data Table D.3: Comparison of the parameters estimated from experimental and CFD data for FBK pulp suspensionC = 3.3% , at Q = 7.9 L/min and N= 1031 rpm and ju =100 kg/ms, for different flow m 0 index values Experimental « = 0.1 test « = 0.25 n= 1 f 0.51 0.40 0.43 0.44 R 0 0 0 0 Ti (sec) 55 20 30 35 r, (sec) 11 62 60 65 T (Sec) 170 270 240 230 (sec) 1150 1282 1284 1245 74 95 96 92 2 T 2 Vfulfy Mixed (L)
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- CFD simulation of mixing dynamics in agitated pulp...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
CFD simulation of mixing dynamics in agitated pulp stock chests Ford, Clara 2004
pdf
Page Metadata
Item Metadata
Title | CFD simulation of mixing dynamics in agitated pulp stock chests |
Creator |
Ford, Clara |
Date Issued | 2004 |
Description | Agitated pulp chests function as low-pass filters to reduce high frequency variability in pulp properties (mass concentration, freeness etc.) ahead of many pulping and papermaking operations. Tests made on both industrial and scale-model chests have shown that their dynamic performance is far from ideal, with a significant extent of non-ideal flow (short circuiting, recirculation and stagnation) possible (Ein Mozzafari et al., 2004). In this work, the flow field of a 1/11th scale-model pulp chest was modeled using commercial CFD software (Fluent) with the pulp suspension treated as a Bingham plastic. A multiple reference frame approach was used with coupling between reference frames made using a velocity transformation. The flow profiles predicted by the simulation agreed qualitatively with the observed profiles in the experiments. The power input predicted by the simulations was slightly higher ( + 12%) than the measured power. The velocity field obtained from the CFD model was used to obtain the system's dynamic response to a frequency-modulated random binary input signal. This data was then used as input to a dynamic model developed by Ein Mozzafari (2002) that treats flow within the chest as following two streams: one that bypasses the mixing zone and one that enters it. For both streams, the fraction of suspension passing through each zone was determined and a time constant and delay time computed. These parameters were then compared with parameters measured experimentally under identical operating conditions and the parameters associated with the mixing zone found to agree to within ±25%. The CFD simulation provides detailed information on the velocity profile within the chest and allows the location(s) of poor mixing regions to be identified. Hence, design methods and trouble shooting of existing stock chests can be significantly improved. |
Extent | 12815782 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-12-03 |
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.0059039 |
URI | http://hdl.handle.net/2429/16307 |
Degree |
Master of Applied Science - MASc |
Program |
Chemical and Biological Engineering |
Affiliation |
Applied Science, Faculty of Chemical and Biological Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2005-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2005-0038.pdf [ 12.22MB ]
- Metadata
- JSON: 831-1.0059039.json
- JSON-LD: 831-1.0059039-ld.json
- RDF/XML (Pretty): 831-1.0059039-rdf.xml
- RDF/JSON: 831-1.0059039-rdf.json
- Turtle: 831-1.0059039-turtle.txt
- N-Triples: 831-1.0059039-rdf-ntriples.txt
- Original Record: 831-1.0059039-source.json
- Full Text
- 831-1.0059039-fulltext.txt
- Citation
- 831-1.0059039.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-0059039/manifest