"Applied Science, Faculty of"@en . "Mechanical Engineering, Department of"@en . "DSpace"@en . "UBCV"@en . "Feng, Xiaosi"@en . "2009-12-23T00:36:25Z"@en . "2005"@en . "Doctor of Philosophy - PhD"@en . "University of British Columbia"@en . "This thesis describes numerical simulations of the mean flow and the large-scale turbulence\r\nstructures in a typical paper machine headbox. Parallel computational methods and a large eddy\r\nsimulation (LES) model are used. Turbulent flow through the converging section is modelled and is\r\nthen used, with a suitable fiber model, to make predictions of the statistical fiber orientation in the\r\nheadbox.\r\nThe present large eddy simulations have been made with a new parallel computational code\r\ndeveloped by the author for this purpose. Techniques for parallel computation of LES are discussed.\r\nIn this code, the incompressible Navier-Stokes equations are solved using a staggered finite volume\r\nmethod for a generalized curvilinear coordinate system. The fractional step method is used to solve\r\nthe momentum equations. A Poisson equation for pressure is obtained to satisfy the continuity\r\nequation. The Smagorinsky constant, the dynamic subgrid scale models and the approximate\r\ndeconvolution model (ADM) are used in the large eddy simulations to account for the effects of the\r\nunresolved turbulent motions. Satisfactory validation of the new code is obtained through\r\ncomparisons with experiments and computations of fully developed channel flow. The parallel\r\nefficiency is higher than 80% for the present parallel computation of LES.\r\nFor the converging section flow representing the headbox, the computed values of the mean\r\nvelocity agree well with the measured values. For the rms fluctuations, there are some differences\r\nbetween the numerical and the measured results over the first third of the converging section. These\r\ndifferences appear to be related to the impossiblity to set accurately the flow in the experiment and\r\nto the inflow data used in the simulations. A modified method for the simulation of inflow conditions\r\nis proposed, and used, and better agreement with experimental data is obtained with the new inflow\r\nconditions. The LES computer code is coupled with a suitable fiber model to predict the statistical\r\norientation of nylon \"fibers\" in the converging section. The numerical methods give statistical results\r\nwhich are similar to existing experimental data for the fiber orientation.\r\nFinally, high Reynolds number large eddy simulations are carried out for a commercial sized\r\nheadbox of generic geometry. Two inflow conditions for the converging section are investigated and\r\nsome conclusions are drawn."@en . "https://circle.library.ubc.ca/rest/handle/2429/17095?expand=metadata"@en . "Large Eddy Simulation of the Turbulence and Fiber Motion in a Headbox By Xiaosi Feng B. Sc., Zhengzhou University, China, 1986 M . Sc., The Institute of Applied Physics and Computational Mathematics, China, 1989 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF THE REQUIREMENTS OF THE DEGREE OF DOCTOR OF PHILOSOPHY in THE F A C U L T Y OF G R A D U A T E STUDIES (Mechanical Engineering) THE UNIVERSITY OF BRITISH C O L U M B I A August 2005 \u00C2\u00A9 Xiaosi Feng, 2005 Abstract This thesis describes numerical simulations of the mean flow and the large-scale turbulence structures in a typical paper machine headbox. Parallel computational methods and a large eddy simulation (LES) model are used. Turbulent flow through the converging section is modelled and is then used, with a suitable fiber model, to make predictions of the statistical fiber orientation in the headbox. The present large eddy simulations have been made with a new parallel computational code developed by the author for this purpose. Techniques for parallel computation of LES are discussed. In this code, the incompressible Navier-Stokes equations are solved using a staggered finite volume method for a generalized curvilinear coordinate system. The fractional step method is used to solve the momentum equations. A Poisson equation for pressure is obtained to satisfy the continuity equation. The Smagorinsky constant, the dynamic subgrid scale models and the approximate deconvolution model (ADM) are used in the large eddy simulations to account for the effects of the unresolved turbulent motions. Satisfactory validation of the new code is obtained through comparisons with experiments and computations of fully developed channel flow. The parallel efficiency is higher than 80% for the present parallel computation of LES. For the converging section flow representing the headbox, the computed values of the mean velocity agree well with the measured values. For the rms fluctuations, there are some differences between the numerical and the measured results over the first third of the converging section. These differences appear to be related to the impossiblity to set accurately the flow in the experiment and to the inflow data used in the simulations. A modified method for the simulation of inflow conditions is proposed, and used, and better agreement with experimental data is obtained with the new inflow conditions. if The LES computer code is coupled with a suitable fiber model to predict the statistical orientation of nylon \"fibers\" in the converging section. The numerical methods give statistical results which are similar to existing experimental data for the fiber orientation. Finally, high Reynolds number large eddy simulations are carried out for a commercial sized headbox of generic geometry. Two inflow conditions for the converging section are investigated and some conclusions are drawn. Table of contents Abstract ii Table of contents iv List of Tables vi List of Figures ; vii Nomenclature xi Acknowledgements xiii Chapter 1 Introduction 1 Chapter 2 Literature Review 4 2.1 Headbox 4 2.1.1 Headbox flow simulations 5 2.1.2 Fiber orientation 6 2.2 Numerical methods 7 2.2.1 Direct Numerical Simulation (DNS) 8 2.2.2 RANS models 9 2.2.3 Large Eddy Simulation 10 Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 12 3.1 Introduction 12 3.2 Governing equations 13 3.2.1 Filtered Navier-Stokes equations 13 3.2.2 Subgrid model 14 3.2.2.1 Dynamic subgrid model 16 3.2.2.2 Approximate deconvolution model (ADM) 18 3.3 Numerical methods 21 3.3.1 The flow solver 21 3.3.2 Numerical methods for solution of the Poisson equation 26 3.3.2.1 Direct method 26 3.3.2.2 Conjugate Gradient method and C G S T A B 26 3.3.3 Parallel computation for solution of the Poisson Equation 27 3.3.3.1 Parallel computation using the direct method 28 3.3.3.2 Parallel computation for C G S T A B method 28 3.3.4 Di screte test filters 30 3.4 Numerical results 31 3.4.1 Reynolds number R e r = 395 32 3.4.2 Results for Reynolds number R e r = 1655 42 3.5 Summary 47 iv Chapter 4 LES for Turbulent Flow through a Converging Section 49 4.1 Introduction 49 4.2 Flow simulation 50 4.2.1 Experimental setup 50 4.2.2 Numerical simulation of the converging section 52 4.2.3 Alternative methods for the generation of inflow data 62 4.3 Discussion 69 4.3.1 Dependence on grid resol ution 70 4.3.2 Effects of the different subgrid models 74 4.4 Summary 77 Chapter 5 Fiber Concentration and Orientation 80 5.1 Introduction 80 5.2 Fiber concentration 81 5.2.1 Fiber model and wall model 81 5.2.2 Numerical results 83 5.3 Fiber orientation in the headbox 85 5.3.1 Initial conditions 85 5.3.2 Numerical results 86 5.3.3 Discussion 98 5.4 Summary 99 Chapter 6 Simulation of a Typical Industrial Headbox 101 6.1 Introduction 101 6.2 Flow simulation 101 6.2.1 Computational geometry 101 6.2.2 Inflow conditions 102 6.2.3 Numerical results 108 6.3 Summary 122 Chapter 7 Summary and Conclusions 124 7.1 Parallel computing 124 7.2 LES for turbulent flow through a converging section 125 7.3 Fiber orientation 127 7.4 Simulation of a typical industrial headbox 127 Chapter 8 Recommendations for Future Work 130 Bibliography 131 Appendix A 139 Appendix B 141 Appendix C 142 v List of Tables T A B L E 3 . 1 EXECUTION TIME IN SECONDS, SPEED UP AND EFFICIENCY FOR DIFFERENT NUMBERS OF PROCESSORS ON 64X64X64 GRID 3 0 T A B L E 3 . 2 GRID PARAMETERS FOR Rer = 395 3 3 T A B L E 3 . 3 GRID PARAMETERS FOR ReT= 1 6 5 5 4 3 T A B L E 4 . 1 T H E GEOMETRY OF THE CONVERGING SECTION 5 2 T A B L E 4 . 2 GRID PARAMETERS. Nx, Ny ,Nz DENOTE THE NUMBER OF MESHES IN THE STREAMWISE, SPANWISE AND WALL-NORMAL DIRECTIONS, RESPECTIVELY. S IS HALF OF THE INLET HEIGHT OF THE CONVERGING SECTION 7 1 T A B L E 5 . 1 T H E FIRST MOMENT ( ^ or , . )A\u00C2\u00AB ( .), SECOND MOMENT (afp(a j)Aa i) AND THIRD MOMENT / i ( o f p{aj)A.ai) OF THE PROBABILITY DISTRIBUTION IN X-Z PLANE 9 4 T A B L E 5 . 2 T H E FIRST MOMENT ( ^ (Xip{jDti)Aor.), SECOND MOMENT (]T a)p{at)Aor,) AND THIRD MOMENT (af p{ai)Aaj) OF THE PROBABILITY DISTRIBUTION IN X-Y PLANE 9 5 VI List of Figures F I G U R E 1 . 1 T H E S C H E M A F O R T H E P R O C E S S O F A P A P E R M A C H I N E ( F R O M S H A R I A T I 2 0 0 2 ) 1 F I G U R E 3 . 1 S C H E M A T I C O F E N E R G Y S P E C T R U M I L L U S T R A T I N G T H E R E L A T I O N S H I P A M O N G R E S O L V E D , R E S O L V E D S U B F I L T E R - S C A L E , A N D S U B G R I D - S C A L E M O T I O N S . CO is T H E W A V E N U M B E R . CQC I S T H E G R I D C U T O F F W A V E N U M B E R ( G U L L B R A N D A N D C H O W 2 0 0 3 ) 1 9 F I G U R E 3 . 2 I L L U S T R A T I O N O F T H E P R I M A R Y C E L L 2 2 F I G U R E 3 . 3 C O N T R O L V O L U M E F O R T H E ^ - M O M E N T U M E Q U A T I O N 2 4 F I G U R E 3 . 4 T R A P E Z O I D A L F I L T E R I N 3 D 3 1 F I G U R E 3 . 5 S C H E M A O F T H E C H A N N E L F L O W 3 2 F I G U R E 3 . 6 D I M E N S I O N L E S S M E A N S T R E A M W I S E V E L O C I T Y P R O F I L E S W I T H T H E A D M M O D E L A N D T H E S M A G O R I N S K Y M O D E L 3 7 F I G U R E 3 . 7 D I M E N S I O N L E S S R M S F L U C T U A T I O N O F U : yju'u luT W I T H T H E A D M M O D E L A N D T H E S M A G O R I N S K Y M O D E L 3 7 F I G U R E 3 . 8 D I M E N S I O N L E S S R M S F L U C T U A T I O N O F V : -JV'V I ur W I T H T H E A D M M O D E L A N D T H E S M A G O R I N S K Y M O D E L 3 8 F I G U R E 3 . 9 D I M E N S I O N L E S S R M S F L U C T U A T I O N O F W : y/w'w' I ur W I T H T H E A D M M O D E L A N D T H E S M A G O R I N S K Y M O D E L 3 8 F I G U R E 3 . 1 0 D I M E N S I O N L E S S R E Y N O L D S S T R E S S : u'v' lu\ W I T H T H E A D M M O D E L A N D T H E S M A G O R I N S K Y M O D E L 3 9 F I G U R E 3 . 1 1 D I M E N S I O N L E S S M E A N S T R E A M W I S E V E L O C I T Y P R O F I L E S F R O M T H E C A S E 2 3 9 F I G U R E 3 . 1 2 D I M E N S I O N L E S S R M S F L U C T U A T I O N O F u: yju'u luT F R O M T H E C A S E 2 4 0 F I G U R E 3 . 1 3 D I M E N S I O N L E S S R M S F L U C T U A T I O N O F v.yjv'v' I uT F R O M T H E C A S E 2 4 0 F I G U R E 3 . 1 4 D I M E N S I O N L E S S R M S F L U C T U A T I O N O F W . yjw'w' I uT F R O M T H E C A S E 2 4 1 F I G U R E 3 . 1 5 S G S M O D E L C O E F F I C I E N T S W I T H T H E M E D I U M M E S H C A S E ( C A S E 2 ) 4 1 F I G U R E 3 . 1 6 R A T I O O F S G S E D D Y V I S C O S I T Y T O M O L E C U L A R V I S C O S I T Y : vTlv W I T H T H E S M A G O R I N S K Y M O D E L A N D T H E A D M M O D E L ( I N C A S E 2 ) 4 2 F I G U R E 3 . 1 7 D I M E N S I O N L E S S M E A N S T R E A M W I S E V E L O C I T Y P R O F I L E S A T Rer = 1 6 5 5 4 4 F I G U R E 3 . 1 8 D I M E N S I O N L E S S R M S F L U C T U A T I O N O F U : yju'u' luT A T Rer = 1 6 5 5 4 5 F I G U R E 3 . 1 9 D I M E N S I O N L E S S R M S F L U C T U A T I O N O F V.JV'V I ur A T Rer = 1 6 5 5 4 5 F I G U R E 3 . 2 0 D I M E N S I O N L E S S R E Y N O L D S S T R E S S : u'v lu1 A T ReT = 1 6 5 5 4 6 F I G U R E 3 . 2 1 D I M E N S I O N L E S S R M S F L U C T U A T I O N O F u : yju'u luT W I T H 9 6 X 1 2 8 X 9 6 M E S H 4 6 F I G U R E 3 . 2 2 D I M E N S I O N L E S S R M S F L U C T U A T I O N O F v: yjv'v I uT W I T H 96x128x96 M E S H 4 7 F I G U R E 4 . 1 F L O W L O O P I N T H E E X P E R I M E N T ( F R O M S H A R I A T I 2 0 0 2 ) 5 1 F I G U R E 4 . 2 E X P E R I M E N T A L H E A D B O X D R A W I N G ( F R O M S H A R I A T I 2 0 0 2 ) 5 1 F I G U R E 4 . 3 C O M P U T A T I O N A L D O M A I N F O R T H E C O N V E R G I N G S E C T I O N . O N L Y A S U B S E T O F T H E A C T U A L G R I D I S P L O T T E D 5 5 F I G U R E 4 . 4 M E A S U R E M E N T L O C A T I O N S A L O N G T H E C E N T E R L I N E O F T H E C O N V E R G I N G S E C T I O N ( L I N E A B C ) A N D T H E V E R T I C A L L I N E A T 0 . 0 3 M A F T E R I N L E T ( L I N E D E ) 5 6 F I G U R E 4 . 5 C O M P A R I S O N O F T H E M E A N V E L O C I T Y A T T H E C E N T E R L I N E O F T H E C O N V E R G I N G S E C T I O N 5 6 F I G U R E 4 . 6 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N u N O R M A L I Z E D B Y M ^ 5 7 vii F I G U R E 4 . 7 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N V N O R M A L I Z E D BY u'jn 5 7 F I G U R E 4 . 8 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N W N O R M A L I Z E D B Y u'n 5 8 F I G U R E 4 . 9 C O M P A R I S O N O F T U R B U L E N C E K I N E T I C E N E R G Y k N O N D I M E N S I O N A L I Z E D B Y k A T T H E I N L E T O F T H E P O S I T I O N 0 . 0 0 6 M A F T E R T H E C O N V E R G E N C E S T A R T S 5 8 F I G U R E 4 . 1 0 C O M P A R I S O N O F T H E F L U C T U A T I N G V E L O C I T Y U A T I N L E T F O R O N E C H A N N E L I N F L O W C O N D I T I O N W I T H T H E M E A S U R E D V A L U E A T T H E V E R T I C A L L I N E ( D E L I N E I N F I G U R E 4 . 4 ) A T T H E L O C A T I O N 0 . 0 3 M F R O M T H E I N L E T 5 9 F I G U R E 4 . 1 1 T H E D I A G R A M O F \" T W O - C H A N N E L \" I N F L O W C O N D I T I O N 5 9 F I G U R E 4 . 1 2 C O M P A R I S O N O F T H E M E A N V E L O C I T Y A T T H E C E N T E R L I N E O F T H E C O N V E R G I N G S E C T I O N W I T H D I F F E R E N T I N F L O W C O N D I T I O N S 6 0 F I G U R E 4 . 1 3 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N U N O R M A L I Z E D B Y M ^ W I T H D I F F E R E N T I N F L O W C O N D I T I O N S 6 0 F I G U R E 4 . 1 4 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N V N O R M A L I Z E D B Y u'n W I T H D I F F E R E N T I N F L O W C O N D I T I O N S 6 1 F I G U R E 4 . 1 5 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N W N O R M A L I Z E D B Y w ' n w i T H D I F F E R E N T I N F L O W C O N D I T I O N S 6 1 F I G U R E 4 . 1 6 S C H E M A T I C O F N E W I N F L O W G E N E R A T I O N T E C H N I Q U E 6 4 F I G U R E 4 . 1 7 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N U A T I N L E T W I T H T H E M E A S U R E D V A L U E A T T H E V E R T I C A L L I N E ( L I N E D E I N F I G U R E 4 . 4 ) O F T H E L O C A T I O N 0 . 0 3 M F R O M T H E I N L E T 6 5 F I G U R E 4 . 1 8 C O M P A R I S O N O F S T R E A M W I S E C O M P O N E N T O F T H E M E A N V E L O C I T Y A T T H E I N L E T W I T H T H E M E A S U R E D V A L U E A T T H E V E R T I C A L D E L I N E ( I N F I G U R E 4 . 4 ) L O C A T E D 0 . 0 3 M F R O M T H E I N L E T 6 5 F I G U R E 4 . 1 9 M E A N S T R E A M W I S E V E L O C I T I E S C A L C U L A T E D W I T H T H E N E W I N F L O W C O N D I T I O N A T T H E V E R T I C A L L I N E ( D E L I N E I N F I G U R E 4 . 4 ) C O M P A R E D W I T H M E A S U R E D V A L U E S A T T H E S A M E L O C A T I O N . . . 6 6 F I G U R E 4 . 2 0 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N U A T T H E V E R T I C A L L I N E ( D E L I N E I N F I G U R E 4 . 4 ) W I T H N E W I N F L O W C O N D I T I O N 6 6 F I G U R E 4 . 2 1 M E A N V E L O C I T Y A T T H E C E N T E R L I N E O F T H E C O N V E R G I N G S E C T I O N C A L C U L A T E D W I T H N E W I N F L O W C O N D I T I O N C O M P A R E D T O M E A S U R E D V A L U E S A T T H E S A M E L O C A T I O N 6 7 F I G U R E 4 . 2 2 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N U N O R M A L I Z E D B Y u'jn C A L C U L A T E D W I T H T H E N E W I N F L O W C O N D I T I O N 6 7 F I G U R E 4 . 2 3 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N V N O R M A L I Z E D B Y u'n C A L C U L A T E D W I T H T H E N E W I N F L O W C O N D I T I O N 6 8 F I G U R E 4 . 2 4 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N W N O R M A L I Z E D B Y u'n C A L C U L A T E D W I T H T H E N E W I N F L O W C O N D I T I O N 6 8 F I G U R E 4 . 2 5 C O M P A R I S O N O F T U R B U L E N C E K I N E T I C E N E R G Y K N O N D I M E N S I O N A L I Z E D B Y A : A T T H E I N L E T O F T H E P O S I T I O N 0 . 0 0 6 M A F T E R T H E C O N V E R G E N C E S T A R T S W I T H T H E N E W I N F L O W C O N D I T I O N 6 9 F I G U R E 4 . 2 6 C O M P A R I S O N O F T H E M E A N V E L O C I T Y A T T H E C E N T E R L I N E O F T H E C O N V E R G I N G S E C T I O N W I T H C A S E 1 A N D C A S E 2 7 1 F I G U R E 4 . 2 7 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N U N O R M A L I Z E D B Y u'n W I T H C A S E 1 A N D C A S E 2 7 2 F I G U R E 4 . 2 8 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N W N O R M A L I Z E D B Y ujn W I T H C A S E 1 A N D C A S E 2 7 2 F I G U R E 4 . 2 9 C O M P A R I S O N O F T H E M E A N V E L O C I T Y A T T H E C E N T E R L I N E O F T H E C O N V E R G I N G S E C T I O N W I T H C A S E I A N D C A S E 3 7 3 F I G U R E 4 . 3 0 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N U N O R M A L I Z E D B Y u'n W I T H C A S E I A N D C A S E 3 7 3 F I G U R E 4 . 3 1 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N W N O R M A L I Z E D B Y u'jn W I T H C A S E I , A N D C A S E 3 7 4 viii F I G U R E 4.32 C O M P A R I S O N O F T H E M E A N V E L O C I T Y A T T H E C E N T E R L I N E O F T H E C O N V E R G I N G S E C T I O N W I T H D I F F E R E N T S U B G R I D M O D E L S 75 F I G U R E 4.33 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N U N O R M A L I Z E D B Y u'in W I T H D I F F E R E N T S U B G R I D M O D E L S 75 F I G U R E 4.34 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N V N O R M A L I Z E D B Y W ' H W I T H D I F F E R E N T S U B G R I D M O D E L S 76 F I G U R E 4.35 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N W N O R M A L I Z E D B Y U W I T H in D I F F E R E N T S U B G R I D M O D E L S 76 F I G U R E 4.36 R A T I O O F S G S E D D Y V I S C O S I T Y T O M O L E C U L A R V I S C O S I T Y : Vt I v W I T H T H E S M A G O R I N S K Y M O D E L 77 F I G U R E 5.1 A F I B E R C O N S I S T I N G O F S P H E R O I D S C O N N E C T E D T H R O U G H B A L L A N D S O C K E T J O I N T S ( F R O M R O S S A N D K L I N G E N B E R G (1997)) 83 F I G U R E 5.2 S C H E M A O F T H E C H A N N E L F L O W 83 F I G U E R 5.3 C O N C E N T R A T I O N F O R 0 .001M, 0 .002M, 0 .003M F I B E R V S . C H A N N E L H E I G H T S C A L E D B Y F I B E R L E N G T H W I T H L E S 84 F I G U R E 5.4 A F I B E R ' S I N I T I A L P O S I T I O N 85 F I G U R E 5.5 M E A S U R E M E N T P O I N T S A L O N G T H E H E A D B O X , U N I T M 90 F I G U R E 5.6 I N I T I A L F I B E R O R I E N T A T I O N D I S T R I B U T I O N , ( A ) I N X - Z P L A N E , ( B ) I N X - Y P L A N E 90 F I G U R E 5.7 F I B E R O R I E N T A T I O N D I S T R I B U T I O N A T x = 0 .045M, ( A ) I N X - Z P L A N E , ( B ) I N X - Y P L A N E 91 F I G U R E 5.8 F I B E R O R I E N T A T I O N D I S T R I B U T I O N A T X = 0 .122M, ( A ) I N X - Z P L A N E , ( B ) I N X - Y P L A N E 91 F I G U R E 5.9 F I B E R O R I E N T A T I O N D I S T R I B U T I O N A T X = 0.1 5 7 M , ( A ) I N X - Z P L A N E , ( B ) I N X - Y P L A N E 92 F I G U R E 5.10 F I B E R O R I E N T A T I O N D I S T R I B U T I O N A T X = 0.1 9 2 M , ( A ) I N X - Z P L A N E , ( B ) I N X - Y P L A N E 92 F I G U R E 5.11 F I B E R O R I E N T A T I O N D I S T R I B U T I O N A T X = 0 .227M, ( A ) I N X - Z P L A N E , ( B ) I N X - Y P L A N E 93 F I G U R E 5.12 F I B E R O R I E N T A T I O N D I S T R I B U T I O N A T X = 0 .262M, ( A ) I N X - Z P L A N E , ( B ) I N X - Y P L A N E 93 F I G U R E 5.13 F I B E R O R I E N T A T I O N D I S T R I B U T I O N A T T H E C O N V E R G I N G S E C T I O N E X I T , 94 F I G U R E 5.14 T H E F I R S T M O M E N T ( ^ a^a^Acc,) A T D I F F E R E N T S T A T I O N S I N X - Z P L A N E 95 F I G U R E 5.15 T H E F I R S T M O M E N T ( ^ )Aflf(.) A T D I F F E R E N T S T A T I O N S I N X - Y P L A N E 96 F I G U R E 5.16 T H E S E C O N D M O M E N T ( ^ a)/?(a,.)Aor,.) A T D I F F E R E N T S T A T I O N S I N X - Z P L A N E 96 i F I G U R E 5.17 T H E S E C O N D M O M E N T ( j^T a]p(at )t\ai) A T D I F F E R E N T S T A T I O N S I N X - Y P L A N E 97 F I G U R E 5.18 T H E T H I R D M O M E N T ( o f p{ai)l!s,ai) A T D I F F E R E N T S T A T I O N S I N X - Z P L A N E 97 i F I G U R E 5.19 T H E T H I R D M O M E N T ( ^ ofp{ai)Aai) A T D I F F E R E N T S T A T I O N S I N X - Y P L A N E 98 F I G U R E 5.20 F I B E R O R I E N T A T I O N D I S T R I B U T I O N A T X = 0.1 5 7 M I N X - Y P L A N E W I T H 8000, 16,000 A N D 24,000 F I B E R S 99 F I G U R E 6.1 T H E G E O M E T R Y O F T H E C O N V E R G I N G S E C T I O N 102 F I G U R E 6.2 O N E - C H A N N E L I N F L O W C O N D I T I O N F O R T H E C O N V E R G I N G S E C T I O N 104 F I G U R E 6.3 S C H E M A T I C D I A G R A M O F T H E \" F I V E - C H A N N E L \" I N F L O W C O N D I T I O N S H O W I N G T H E D E T A I L S O F O N E O F T H E F I V E I D E N T I C A L C H A N N E L S 105 F I G U R E 6.4 C O M P A R I S O N O F T H E M E A N V E L O C I T Y A T T H E I N L E T 106 F I G U R E 6.5 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N U A T T H E I N L E T 106 F I G U R E 6.6 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N V A T T H E I N L E T 107 F I G U R E 6.7 C O M P A R I S O N O F R M S V A L U E O F T H E V E L O C I T Y F L U C T U A T I O N W ' A T T H E I N L E T 107 F I G U R E 6.8 C O M P A R I S O N P O S I T I O N S A L O N G T H E H E A D B O X 110 F I G U R E 6.9 C O M P A R I S O N O F T H E M E A N V E L O C I T Y A T T H E F I R S T S T A T I O N A B 111 ix FIGURE 6 . 1 0 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION U AT THE FIRST STATION A B . . . 1 1 1 FIGURE 6 . 1 1 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION V AT THE FIRST STATION A B . . . . 1 1 1 FIGURE 6 . 1 2 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION W AT THE FIRST STATION A B . . . 1 1 2 FIGURE 6 . 1 3 COMPARISON OF THE MEAN VELOCITY PROFILES AT THE SECOND STATION C D 1 1 2 FIGURE 6 . 1 4 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION U AT THE SECOND STATION C D . 1 1 3 FIGURE 6 . 1 5 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION V AT THE SECOND STATION C D . 1 1 3 FIGURE 6 . 1 6 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION W AT THE SECOND STATION C D . 1 1 4 FIGURE 6 . 1 7 COMPARISON OF THE MEAN VELOCITY AT THE THIRD STATION E F 1 1 4 F I G U R E 6 . 1 8 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION U AT THE THIRD STATION E F . . . . 1 1 5 FIGURE 6 . 1 9 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION V AT THE THIRD STATION E F . . . . 1 1 5 FIGURE 6 . 2 0 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION W AT THE THIRD STATION E F . . . 1 1 6 FIGURE 6 . 2 1 COMPARISON OF THE MEAN VELOCITY AT THE OUTLET 1 1 6 FIGURE 6 . 2 2 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION U AT THE OUTLET 1 1 7 FIGURE 6 . 2 3 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION V AT THE OUTLET 1 1 7 FIGURE 6 . 2 4 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION W AT THE OUTLET 1 1 8 FIGURE 6 . 2 5 COMPARISON THE MEAN VELOCITY PROFILES AT THE OUTLET WITH THE SMAGORINSKY MODEL. 1 1 8 FIGURE 6 . 2 6 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION U AT THE OUTLET WITH THE SMAGORINSKY MODEL 1 1 9 FIGURE 6 . 2 7 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION V AT THE OUTLET WITH THE SMAGORINSKY MODEL 1 1 9 FIGURE 6 . 2 8 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION W' AT THE OUTLET WITH THE SMAGORINSKY MODEL 1 2 0 FIGURE 6 . 2 9 COMPARISON OF THE RATIO v lu AT DIFFERENT POSITIONS 1 2 0 FIGURE 6 . 3 0 COMPARISON OF THE RATIO vJ lu AT DIFFERENT POSITIONS 1 2 1 FIGURE 6 . 3 1 COMPARISON OF THE RATIO v I u AT THE OUTLET WITH THE TWO DIFFERENT INFLOW CONDITIONS 1 2 1 FIGURE 6 . 3 2 COMPARISON OF THE RATIO w lu AT THE OUTLET WITH THE TWO DIFFERENT INFLOW CONDITIONS 1 2 2 X Nomenclature Cs Smagorinsky constant G grid filter function G test filter k kinetic energy of turbulence n unit normal p filtered pressure R b Reynolds number based on bulk velocity RT Reynolds number based on friction velocity S, Stokes number - 1 du 3w, Su =\u00E2\u0080\u0094(\u00E2\u0080\u0094'- + -\u00E2\u0080\u0094-) filtered rate of strain tensor 2 dXj dXj | S |= yJ2S~S~ magnitude of the strain tensor t time zT(. filtered velocity U,V,W mean velocity in x, y, z direction, respectively Uh bulk velocity V velocity vector x, y, z distance parameter y + non-dimensional distance Greek Letters x i \u00C2\u00A3 M P To vT A A V Superscripts Kronecker delta dissipation rate of turbulence dynamic viscosity density Reynolds-stress tensor eddy viscosity grid filter width test filter width Kolmogorov microscale nondimensionalized value xii Acknowledgements I would like to express my deepest appreciation to my research supervisors, Dr. Martha Salcudean and Dr. Ian Gartshore for their invaluable guidance and kind support that made this work possible. I would also like to thank my research committee member, Dr. W. Kendal Bushe for his technical guidance and helpful suggestions. I gratefully acknowledge Dr. Mark Martinez for participating in my research committee and for his help. Special thanks are due to Dr. Ned Djilali and Dr. Brian Wetton for their inspiration and useful discussions. I would like to thank to my colleagues at the University of British Columbia (UBC) specially Mohammad R. Shariati and Xun Zhang for providing experimental data for comparison. I would like to thank Dr. Matthew W. Choptuik and Dr. Doug Phillips for their generosity and technical supports in the use of the PIII/Linux and P4-Xeon clusters at UBC and the MACI Alpha cluster at the University of Calgary. I also wish to acknowledge the financial assistance from FRBC and the Natural Science and Engineering Research Council of Canada. Finally, I am deeply indebted to my wife Suqin, my daughter Manning, my parents and my brothers and sisters for their love, infinite patience, unconditional support and help throughout all the years of work. xiii Chapter 1 Introduction The headbox is the first component of a paper machine in the papermaking system (Figure 1.1). The fiber orientation and basis weight profiles of the paper depend on the fluid flow in the headbox, which makes a headbox critical to a successful papermaking system. Drainage Figure 1.1 The schema for the process of a paper machine (from Shariati 2002). The headbox can be several meters wide in the cross machine direction and has a contraction ratio (inlet area to outlet area) of 10 or more in a streamwise length of less than a meter. A dilute concentration of pulp fibers in water flows steadily through the headbox. The detailed introduction of the headbox will be given in chapter 2. One of its functions is to provide a discharge of stock (fibers and other substances in fluid suspension) at its outlet that is uniform in both the machine direction and across the width of the paper web (the cross direction). The stock is fed onto the forming section wire mesh through a slice (exit of the headbox) opening. The forming section will drain water from the stock, leaving behind the fibers that naturally bond together into a sheet. A common way to speed up the papermaking process is to use a twin-wire 1 Chapter 1 Introduction 2 former that squeezes the water off simultaneously from above and below. Following the forming section, the paper web then passes through the presses, in which the rolls squeeze out excess water. After that, the web proceeds to the drying section, which consists of a set of steam-heated drums. Then the paper web goes through a calender. The calender is composed of a set of high-pressure rolls to give the web smoothness and increase its density. Finally the web is carried to the paper machine reel to be made into jumbo rolls. To achieve good paper formation, there must be no fiber flocculation in the jet exiting the headbox. The fibers must also have an even distribution in the cross direction of the machine. Thus the flow in the headbox should have a minimal component of mean velocity in the cross-machine direction. At the same time, the maintenance of high intensity turbulence within the stock is required in headbox designs to prevent fiber flocculation. The turbulence characteristics are critically important, since the turbulence directly influences fiber flocculation and orientation. Turbulence is a complex but common form of fluid motion. In many practical numerical simulations, the RANS (Reynolds averaged Navier-stokes) model is used to simulate the turbulent flow, together with a specific turbulence model such as the k-\u00C2\u00A3 equation model. This provides results that agree well with the experimental data in simple flows. But for complex flows, such as the flow through the converging section of a headbox with its large and rapid contraction, the turbulence may become very anisotropic. In these cases, conventional turbulence models such as k- \u00C2\u00A3 can not simulate the turbulence well. There are two ways to solve this problem. One way is to modify the k-\u00C2\u00A3 equation model (using for example, a non-linear k-\u00C2\u00A3 model such as that of Speziale (1987)) or using other models for the RANS equation. These modifications are not universal and require a test for each modified model. The other way is to use large eddy simulation methods. These are much more universal because they simulate more of the physical features of turbulence and use fewer model constants than RANS models. Chapter 1 Introduction 3 The main objectives of the thesis are to investigate the turbulent flow and fiber orientation in the converging section of the headbox and to develop a computational tool necessary to carry out the investigation. These objectives are listed below in the order they have been approached: 1. To develop a practical numerical tool based on LES which is capable of simulating the large scales of turbulence in a headbox. 2. To study parallel computing algorithms and develop an appropriate LES parallel code, that will allow simulations of the turbulent flow in the converging section of the headbox to be made. 3. To simulate the turbulent flow in the converging section. 4. To study the fiber orientation in a headbox coupling LES with the fiber model. 5. To simulate the turbulent flow in an industrial headbox. Following this introductory chapter, a literature review is given in chapter 2. The headbox and the numerical methods to simulate the turbulence are also described in chapter 2. Chapter 3 describes the numerical results of LES for fully developed channel flow. In chapter 4, the results of LES for the turbulent flow through the converging section are presented. Chapter 5 describes the fiber motion through the converging section and the fiber concentration calculation for fully developed channel flow. The LES simulation for a typical industrial-scale headbox is described in chapter 6. The summary and conclusions are presented in chapter 7. Finally, recommendations for future work are given in chapter 8. Chapter 2 Literature Review 2.1 Headbox The function of the headbox is to supply a continuous flow of fibers and other substances in a fluid suspension at constant velocity both across the width of the paper machine and with respect to time. The formation and uniformity of the final paper product are dependent on the uniformity of flow through the headbox, so a good understanding of the headbox is necessary for a successful papermaking system (Smook 1992, Britt 1970 and Gavelin 1998). The operational objectives of the headbox system are listed by Smook (1992) as follows: \u00E2\u0080\u00A2 Spread stock evenly across the width of the machine. \u00E2\u0080\u00A2 Level out cross-currents and consistency variations. \u00E2\u0080\u00A2 Level out machine direction velocity gradients. \u00E2\u0080\u00A2 Create controlled turbulence to eliminate fiber flocculation. \u00E2\u0080\u00A2 Discharge fluid and fibers evenly from the slice opening so that they impinge on the forming fabric at the correct location and angle. Headboxes can be classified as open, air-cushioned and hydraulic (Smook 1992). In the open headbox, the stock surface is open to the atmosphere and the flow at the slice is varied through adjustments to the level of stock in the headbox. Open headboxes were used in early paper machines. As paper machine speeds increased, it became impractical to increase the height of stock, so the air-cushioned headbox was developed. 4 Chapter 2 Literature Review 5 Air-cushioned headboxes use rotating hollow cylinders with large perforations in their walls to provide a high turbulence and to smooth velocity gradients. These boxes are pressurized by the hydraulic feed. A pressurized air chamber is incorporated at the top of the box (Peel 1999). Hydraulic headboxes usually have no air chamber while the air-cushioned types do. They comprise most of the modern headboxes. Practically all headbox research today deals with hydraulic headboxes without air-cushions (Gavelin 1998). This thesis will therefore focus on hydraulic headboxes. Headboxes can be divided into three sections according to the principal flow processes involved (Britt 1970): flow distribution, flow rectification and jet development. In the first section, the flow of stock is distributed across the width of the machine through a manifold distributor. In the second section, the flow from the distributor is rectified. In a hydraulic headbox, a tube bank is often used in the second section. The wall friction in the tubes dampens flow disturbances originating in the stock approach system, provides a resistance to flow which creates a more even flow distribution across the section and creates turbulence to prevent fiber flocculation in the paper-machine forming zone. The processes in the tube bank may include mixing and blending of separate flows from a distributor, eliminating undesirable cross-flows and eddies, improving the velocity profile and developing turbulence of desired scale and intensity (Dahl and Weiss 1975). The third section is responsible for developing the jet and delivering the stock to the sheet forming section. An ideal headbox should produce a uniform and stable jet over the width of the machine. 2.1.1 Headbox flow simulations Shimizu and Wada (1992) have calculated the flow of water in a generic headbox using the k- \u00C2\u00A3 model. They calculated the relationship between the varying height in the contracting part and the flow velocity distribution at the exit. The flow distribution was investigated in two dimensions and the jets from the diffuser tubes were modelled three-dimensionally. Hamalainen (1993) used a Chapter 2 Literature Review 6 finite element method to model the manifold, the turbulence generating section, and the slice in two dimensions. Aidun and Kovacs (1995) have focussed on the study of secondary flows in the headbox and their effects on non-uniform fiber orientation and mass. A non-linear k - f model developed by Speziale (1987) was employed. They reported surprisingly large secondary flows at the headbox exit. Bandahakavi and Aidun (1999) investigated the turbulent flow in the converging section of a headbox using the RNG k - f model and the Reynolds Stress Model (RSM). The RNG k - f model is an improved version of the standard k - f model (Chowdhury 1993). They found that the results obtained from the RSM model are superior to those from the RNG k - f model. Parsheh and Dahlkild (1999) compared the results of modeling the turbulence with different models. It was observed that the turbulence energy computed by the k - f (linear) model is exaggerated. He et al. (1998), Hua et al. (1999, 2000) and Pougatch et al. (2005) have done a complete modeling of the manifold along with all the individual tubes with the k - f model in three dimensions. Shariati et al. (2000, 2001) and Shariati (2002) studied the headbox flow both experimentally and numerically exploring various models for the turbulence (not including LES). 2.1.2 Fiber orientation Fiber orientation refers to the angular distribution of fibers relative to the paper-machine direction (MD). Fiber orientation, both in the plane of the paper and in the paper thickness direction, significantly affects the paper quality (Loewen 1997). Ullmar and Norman (1997) and Ullmar (1998) investigated experimentally the effect of the headbox contraction ratio on fiber orientation. They found that the effect of the contraction ratio is more significant on fiber orientation than that of the flow rate. The fibers have been found to be more strongly orientated in the machine direction for higher contraction ratios. Zhang (2001) measured the orientation of dyed nylon fibers moving in a pilot plexiglass headbox and conducted the numerical simulation of fiber orientation with a fiber motion model. The fiber orientation distribution was Chapter 2 Literature Review 7 calculated using the mean flow only. The fibers were seen to be more strongly oriented in the machine direction by the numerical predictions than was observed in the experiments. He concluded that the anisotropy of the fiber orientation in the headbox flow was caused not only by the mean flow field characteristics, but also by the turbulence characteristics, and the explicit effects of the turbulence should be included into the fiber model. Olson (2002) provided an analytic expression for the fiber orientation distribution in a headbox flow, neglecting the effect of turbulent dispersion. His analysis showed: \u00E2\u0080\u00A2 Fiber orientation distribution is independent of the flow rate through the headbox. \u00E2\u0080\u00A2 Fibers are more strongly oriented in the plane of the contraction than in the plane of the paper. \u00E2\u0080\u00A2 The shape of the headbox does not affect the final fiber orientation distribution leaving the headbox. \u00E2\u0080\u00A2 The contraction ratio is the only geometrical factor affecting fiber orientation distribution at the slice. Olson et al. (2004) proposed an Eulerian model of a turbulent fiber suspension through a one-dimensional headbox to predict the fiber orientation distribution. The contributions discussed earlier have helped considerably to understand better the flow and fiber motion in headboxes. However, a better representation of the turbulence in the headbox and its effect on fiber motion and distribution is necessary to advance the knowledge of paper making process and equipment design. The thesis objective is to address this issue. 2.2 Numerical methods Generally, there are three numerical methods to simulate turbulence: direct numerical simulation, RANS models and large eddy simulation. Chapter 2 Literature Review 8 2.2.1 Direct Numerical Simulation (DNS) A direct numerical simulation means a complete time dependent solution of the Navier-Stokes and continuity equations. The main idea of this method is to assume that the size of the smallest eddies in any turbulent flow is of the order of the Kolmogorov microscale (7]~ (v / \u00C2\u00A3) 4 , where V is kinematic viscosity, and \u00C2\u00A3 is the energy dissipation rate of the turbulent flow). A huge number of grids are required in space in order to resolve these eddies. Unfortunately, DNS can only be applied to very small Reynolds number problems, because the necessary number of cells is roughly given by (Reynolds 1989) 9 N=0(Re 4 ), where Re is defined by the bulk velocity and the characteristic length of the flow. For a typical headbox, this requirement suggests the need for a very large number of cells, more than 1 0 9 , which is impossible using today's computational systems. In spite of this restriction, DNS has been applied to many relatively simple problems, such as homogenous isotropic turbulent flow and turbulent boundary layer and channel flow (Clark et al. 1978, Kim et al. 1987, Spalart 1989, Moser et al. 1999). For channel flow, the grids required by DNS are Re\u00E2\u0080\u009E N (grid number) 12,300 6.7X106 30,800 4 . 0 X 1 0 7 61,600 1.5x10s 230 ,000 2 . 1 X 1 0 9 Chapter 2 Literature Review 9 U H Here R e H = \u00E2\u0080\u0094-\u00E2\u0080\u0094, Um is the average velocity across the channel and 2 H is the height of the channel. DNS techniques are very accurate (assuming the numerical error to be very small by using a spectral method or other high accuracy methods). The results can be used like experimental data to calibrate other models. But DNS is restricted to low Reynolds numbers. 2.2.2 RANS models From a practical viewpoint, the RANS model provides a simple and powerful tool to obtain solutions to many flows. For incompressible turbulent flow, the velocity w, and kinematic pressure p can be decomposed into ensemble means and fluctuation parts as: u,=ut + u' P = P+p', where an overbar is an ensemble mean ( a time average can be used in statistically steady turbulent flows, and spatial averages can be used in homogeneous turbulent flows ). The mean velocity and pressure are solutions of the Reynolds averaged Navier-Stokes and continuity equations which are given for constant density and viscosity by 3w. _ dut 1 dp r - 2 - 1 ^Tij \u00E2\u0080\u0094- + u \u00E2\u0080\u0094\u00E2\u0080\u00A2- = - i - + vV w, \u00E2\u0080\u0094 L (2.1) dt dXj p dx., p dXj f^ = 0, (2.2) OX, where Ttj = p utUj is the Reynolds stress tensor. The Reynolds stress tensor needs a model or another equation (for higher order closure). The additional model is based on a combination of physical reasoning, dimensional analysis and empiricism. RANS models can be divided into two groups, eddy-viscosity models and shear-stress models, respectively. In eddy-viscosity models, the Reynolds stresses are modelled by: Chapter 2 Literature Review 10 -u-v^;^-^k ( 2 3 ) where vT is referred to as the effective turbulent viscosity, 8tJ is the Kronecker delta and k = u. ut 12 . The shear-stress model, however, computes the Reynolds stress directly, without assuming the existence of an eddy viscosity. There are many models: zero-equation model (such as the mixing length model), one equation model (k, / , etc.), two-equation models (k \u00E2\u0080\u0094 \u00C2\u00A3, k \u00E2\u0080\u0094 CO, etc.) and second-order closure models. A comprehensive overview of these models can be found in Wilcox (1998). The RANS model is fairly efficient for many simple flows, so that results can often be obtained relatively quickly. But there are empirical choices of model constants for different flows. For complex flows, these constants are difficult to determine. The simpler k \u00E2\u0080\u0094 \u00C2\u00A3 model essentially assumes isotropy of the turbulence by assuming a scalar effective viscosity that is independent of direction and k is calculated based on this assumption. The normal stresses like w,2 can be 2 computed from the assumption of \u00C2\u00AB,2 =\u00E2\u0080\u0094 k \u00E2\u0080\u0094 VT\u00E2\u0080\u0094\u00E2\u0080\u0094 in (2.3). With this assumption used in a 3 ox, converging section of a headbox (where is very large near the exit), w,2 can become negative OJC, which of course is unphysical. This leads to the production term not being calculated accurately. The production term in the simpler k \u00E2\u0080\u0094 \u00C2\u00A3 model is highly over estimated in comparison with the similar term in other turbulence models (Shariati 2002). 2.2.3 Large Eddy Simulation From the above discussion, the major disadvantage of DNS is that it can only deal with low Reynolds number problems using currently available computers. The major defect of RANS models is that they are not truly universal for different turbulent flows and are even inapplicable to Chapter 2 Literature Review 11 some flows. LES, however, has advantages over both DNS and RANS. It employs a mesh resolution coarser than that used by DNS but fine enough to resolve the large-scale eddies and uses a turbulence model only for those eddies whose size is smaller than the mesh spacing. This model is called the SGS (Subgrid Scale) model. LES reduces computing time compared to DNS so that it can be applied to relatively high Reynolds number problems. The results of LES research can certainly be used to help improve engineering models of turbulence. Although LES is more economical than DNS (typically requires 5-10% of the CPU time need for DNS), the method still requires large computer resources. Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 3.1 Introduction Large eddy simulation is used to model the flow in a headbox because the contraction in the headbox produces unusual anisotropic turbulence that has an important effect on the fiber distribution leaving the headbox. The overall idea in the LES method is to filter out the smaller scales of turbulence, and to solve explicitly for the detailed unsteady motion in the remaining larger scale motions, using a simplified isotropic assumption for the effects of the smaller scales. Domain-decomposition is used to handle geometrically complicated domains such as those occurring in the headbox contraction. The finite volume method and a staggered grid are applied for discretization in space. A combination of a two-stage Adams-Bashforth scheme and the cyclic reduction method for the Poisson-like pressure equation is used to advance the solution in time. Fully developed turbulent channel flow was tested as a code validation case. It is an ideal model flow for simulations of wall-bounded turbulence. Since the flow is fully developed, it is statistically homogeneous in the streamwise direction. Thus, periodic boundary conditions can be used in that direction, avoiding the problems associated with inflow and outflow boundary conditions. Fully developed channel flow has been studied numerically by other people using LES, for example, Deardorff (1970), Schumann (1975), Moin and Kim (1982), Horiuti (1987) and Piomelli, et al. (1988). 12 Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 13 The filtered Navier-Stokes equations and subgrid models that are used are described in section 3.2. The numerical methods are presented in section 3.3. In section 3.4, the numerical results are provided. A summary is given in section 3.5. 3.2 Governing equations 3.2.1 Filtered Navier-Stokes equations The incompressible Navier-Stokes equations for constant viscosity v can be written in the following form \u00C2\u00A3 = 0 (3.1) ox(. du. | diijUj _ 1 o?J | v d\ dt dxj /Jctc, dxdxj where the indices i,j-1,2,3 refer tox,y and z. Einstein's summation convention is applied, and Uj is the velocity in i -direction. A filtering operation is performed on the Navier-Stokes equations to remove the small spatial scales. The resulting equations that describe the space-time evolution of the \"large eddies\" contain the subgrid scale stress tensor that describes the effect of the unresolved small scales on the resolved larger scales. Any flow variable (f>, in the fluid domain D , can be decomposed into a large-scale part (resolvable-scale filtered variables) (/> and a small-scale part (subgrid scale variables) respectively, as ' where Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 14 (j) = \G(x-x*, A)0 is a relaxation coefficient. The relaxation coefficient %u can be determined dynamically (Stolz et al. 2001). The second-order structure function (Lesieur and Metais 1996) is applied to (I \u00E2\u0080\u0094 QN*G)*u . The structure function F2 is defined as F2(^) = \ \ ^ + r,t)-^j)l^ (3.26) where h is the grid size in computational space \u00C2\u00A3 = (^,\u00C2\u00A32,\u00C2\u00A33) , \"-1)) (3.41) Ul^2 = Ui;m + -^R3(p^) (3.42) where H2(Ul) and H3(Ul) are the operators representing the discretized convective, diffusive terms , and (P) and R^ (P) are the operators of the discretized pressure terms. Substituting (3.38), (3.40) and (3.42) into (3.33) gives the following discrete Poisson equation -Div(U't) = AtDiX-J\u00C2\u00B1E-J-) (3.43) \"m+1/2 or Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 26 A v ( 4 ^ ) = - ^ A v ( 0 (3-44) where m-i,j, or k for / = 1,2 or 3, respectively. In the next step, the pressure field is obtained by solving the Poisson equation (3.44). Finally, the new velocities at the (n +1) time step are computed from (3.38), (3.40) and (3.42). 3.3.2 Numerical methods for solution of the Poisson equation 3.3.2.1 Direct method If there is one periodic direction, we can use fast a Fourier transform (FFT) in that direction, resulting in a two-dimensional Helmholtz equations. The cyclic reduction method (Golub and Van Loan (1993)) is used to solve the Helmholtz equations. Finally an inverse FFT is then applied to obtain the three-dimensional solution. In our code, V F F T P A C K and FISHPACK libraries from netlib were used to implement the fast Fourier transform and cyclic reduction methods. This approach is very fast. The exact solution is found in a finite number of operations. The pressure Poisson equation in the channel flow computation is solved with this approach. Note however that this solution method is restricted to special cases: those having one direction of periodicity and orthogonal computational meshes. 3.3.2.2 Conjugate Gradient method and CGSTAB The Poisson equation (3.44) can be written into a linear system, AU = Q, where the matrix A is sparse. The method of Conjugate gradients (CG) is often used for solving sparse positive definite systems. This method is iterative in that it generates a sequence of vectors {u0,ux,u2,...} that end when \\rk \\2=\\Q- Auk \\2 is small enough. The effectiveness of this method for sparse matrices lies in the fact that only the computation of Ap for some vector p is required. This can be Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 27 done quickly since A contains only a few non-zero entries in each row. Most often, CG is used in combination with some kind of preconditioning to reduce the total number of iterations (Golub and Van Loan (1993)). The CG algorithm is presented in Appendix A. The CG method is applicable only to symmetric positive definite systems. To apply the CG method to systems of equations that are not symmetric, such as the pressure Poisson equation in non-orthogonal meshes, we need to use a different method. One of the most stable and robust conjugate gradient methods for solving non-symmetric systems, Bi-Conjugate Gradient Stabilized (BI-CGSTAB) is used. It is also known as CGSTAB (Van den Vorst 1992). The CGSTAB algorithm can also be found in Appendix A. The CGSTAB algorithm performs approximately twice as many operations as the preconditioned CG algorithm. Although the fast multi-grid method (Hackbusch 1985) can also be used to solve the Poisson equation, it is difficult to use with non-orthogonal meshes. 3.3.3 Parallel computation for solution of the Poisson Equation The rapid growth in the capability of single processor computers has slowed in recent years. It now appears that further increases in speed will require multiple processors, i.e. parallel computers. One of the characteristics of modern distributed memory parallel computers is very high computational potential. This is very attractive for the simulation of complicated turbulent flow problems since these problems, when posed in complex geometries can be very large and time consuming. Furthermore, these problems cannot be solved easily with high grid resolution due to limitations of memory on a single workstation. Parallel computation can reduce the computing time and memory requirements per processor. Message-Passing Interface (MPI), the de facto message-passing standard, is a library specification for message passing. It is designed for high performance on both massively parallel machines and on workstation clusters. It is an efficient match to the hardware available at UBC. Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 28 3.3.3.1 Parallel computation using the direct method For parallel calculation, the processors are decomposed into one direction. First, FFT is used in the periodic direction. Then ALL-TO-ALL (MPI command) is used to exchange the data. The rest is the same as the sequential calculation. 3.3.3.2 Parallel computation for CGSTAB method For the CGSTAB method, the Incomplete LU (ILU) is used as a preconditioner (Gustafsson 1978, Ferziger and Peric 2001). These preconditioners, which have proven their effectiveness on serial machines, have an inherent sequentiality which makes them difficult to parallelize. The ILU preconditioner has a known factorization of M \u00E2\u0080\u0094 LU where M, which is known as the preconditioning matrix, L is a lower triangular matrix with the same non-zero structure as the lower triangular part of A, and U is an upper triangular matrix U of similar form (see Appendix B). The solution of a linear system involving M requires a forward and backward substitution in solving lower and upper triangular systems respectively. Parallelization of the construction of L, U and solution of the triangular systems is comparable. The parallel implementation of Lx = b is only considered here. The dependencies can be seen immediately. In the first step, the equation is solved with no dependencies because the first row of L contains only one element. In the second step, the results can be computed once the first solution is known and used as a substitution. Each of those solutions can be substituted into two (in 2D) or three (in 3D) other dependent equations, and so on. A staircase parallelization algorithm is used for the parallel implementation (Bastian and Horton 1991, Vuik et al., 1998). One advantage of this method is that the serial and parallel versions of this method have the same behavior with respect to convergence. The rectangular computational domain is decomposed into p strips parallel to the y-axis. For a 5-point stencil, unknowns (i, j) depend only on and The algorithm is shown as follows: in the first step, Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 29 processor 1 calculates unknown (7,1) and sends the boundary value to processor 2. The processors 2 and 3 are idle. Then processor 2 calculates (/,1) and sends the boundary value to processor 3. At the same time, processor 1 calculates (/, 2) in parallel and exchanges the boundary data to processor 2. Then processor 3 calculates (/, 1). After some start-up time, all processors are running. For the 3D case, Joubert et al., (2000) extend the staircase parallel strategy. The basic concept of this algorithm is based on imposing 2D processor decomposition on the 3D grid and changing the computation order so that the processors can be used with computations as quickly as possible while the processor idle time is minimized. This algorithm is used to parallel CGSTAB method in our code. For parallel computation with LES, the domain decomposition approach is used to partition the data. Each sub-domain data is assigned to one processor to calculate the quantities and exchange the boundary data with neighbors. Two extra rows of cells are added for each block as ghost cells (sub-domain boundaries). In one time step, each processor solves the same equations with each sub-domain data and passes the boundary information to its neighbors. It then proceeds to the next time step. Table 3.1 shows the execution time of the LES for the converging section in 100 time steps. The computational time reduces significantly with 8 processors. For measuring the performance of the two different systems, the following true speedup and efficiency definitions are used (Vuik et al., 1998): \u00E2\u0080\u009E. . wall clock time for a run on one processor S(p) = \u00C2\u00A3 , wall clock time for a run on p processors E(p) = S(p)/p. From the Table 3.1, the efficiency is higher than 80%. Furthermore, the parallel approach for LES reduces the memory usage with processors. It is therefore possible to compute large sized problems. Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 30 Cases 1 4 8 LES for converging section (100 time steps) 791.79 232.11 122.56 Speedup S(p) 1 3.41 6.46 Efficiency E(p) 1 0.85 0.80 Table 3.1 Execution time in seconds, speed up and efficiency for different numbers of processors on 64X64X64 grid 3.3.4 Discrete test filters The filtering operation in one dimension is defined as f(xj= j G(x-x')f(x')dx'. (3.45) j c \u00E2\u0080\u0094 A x The box-filter then becomes G(x-x) = \u00E2\u0080\u0094 if \x-x\Ax. (3.47) The trapezoidal rule is applied to the numerical integration (Burden and Faires (1993)). The result is the following second-order accurate trapezoidal filter fix) = - (f(x - Ax) + 2/(x) + f(x + Ax)). (3.48) 4 Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 31 Trapezoidal filter Figure 3.4 Trapezoidal filter in 3D. The filters are used serially in each direction for three dimensions. The three dimensional trapezoidal filter is shown in Figure 3.4. 3.4 Numerical results The lengths of the computational domain in the streamwise (x) and spanwise (z) direction are 3.OH and 1.5H respectively where H is the channel height. Periodic boundary conditions are imposed in the x and z directions. The schema of the channel flow is shown in Figure 3.5. Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 32 Figure 3.5 Schema of the channel flow There is no wall model used in the simulations to be presented here. For the Smagorinsky model, a Driest damping function is used in (3.7). A no-slip wall condition is imposed in y direction. To resolve the near wall region, a nonuniform grid is used in the wall direction as in Moin and Kim (1982). The y coordinate is then: y j =^H(pmh(a0TjJ)/a + l) Q = 0,\,...,Ny-l). (3.49) where Ny is the total number of grid points in y direction. The quantity a is a stretching parameter (0, where ui is the instantaneous velocity and < uj > is the tempo-spatial average of the velocity field. Since x and z are two homogeneous directions, the data are collected at each xz -plane normal to the walls and averaged in each time step. The computations are conducted with three SGS models: a constant Smagorinsky, a dynamic model and the A D M model. The Smagorinsky model uses the damping function of (3.7), where Cs is set to 0.1. The dynamic model computes the constant from (3.20). The relaxation coefficient %u is computed from (3.27) in the A D M model. The value of deconvolution order N is 5 if not specified. As discussed in 3.2.2.1, the dynamic model involves two filters: the test filter and the grid filter. The trapezoidal rule is applied for test filter, which is illustrated in Figure 3.4. The total viscosity is restricted to positive values: vtol =max(0,v+vT). Ghosal, et al. (1995) and Zang, et al. (1993) found that this restriction improved the numerical stability. The dimensionless mean streamwise velocity profiles are shown in Figure 3.6 with the A D M model and the Smagorinsky model. The coarse mesh (Case 1) does not predict the logarithmic layer well. It shows some differences when compared with the direct numerical simulation (or DNS) data conducted by Moser, Kim and Mansour (1999), hereafter denoted by M K M . The medium mesh (Case 2) and fine mesh (Case 3) agree well with the profile of M K M . The results of the A D M model are better than those of the Smagorinsky model. Figures 3.7-3.9 present the resolved rms fluctuations: y]ujui IuT with two different models. Figure 3.10 shows the Reynolds stress u'v I u2 . Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 35 Figures 3.7-3.9 show that the fluctuations from Case 1 are significantly different from those of M K M . This can be arributed to an inadequate gird resolution. The near wall region cannot be resolved with such a coarse mesh. The effective shear stress on the wall is reduced and so the streamwise non-dimensional fluctuations are overpredicted. The resolved fluctuations normal to the wall are underpredicted. The fluctuations found in Case 3 from two different subgrid models are close to those of M K M . The results from Case 2 are closer to those of Case 3. Figure 3.11 shows the dimensionless mean streamwise velocity profiles from case 2 using three different SGS models. The improvement from the A D M model is clear in the comparison of the mean velocity profiles in Figure 3.11. The incremental improvement of mean velocity between A D M (TV = 5) and ADM (iV = 10) is small, indicating that good reconstruction of the unfiltered velocity is likely obtained with N = 5. The resolved rms fluctuations from the medium mesh case are presented in Figures 3.12 \u00E2\u0080\u0094 3.14. In the near wall region, the maximum value of the fluctuation urms for the A D M model is close to DNS value in Figure 3.12. The differences between the A D M model and the Smagorinsky model are small apart from wall regions. In Figure 3.15, the dynamic coefficient yJC(y) and the van Driest damping on the Smagorinsky coefficient C s ( l \u00E2\u0080\u0094 e~y 125) are compared from the medium mesh case. The peak value of the dynamic coefficient is 10% higher than the constant Smagorinsky coefficient. In most of the computational region, the difference between the dynamic model coefficient and the Smagorinsky coefficient is less than 2%. So the results of the dynamic model and the Smagorinsky model are almost the same from Figures 3.11-3.14. The ratio of the SGS eddy viscosity to the molecular viscosity in the direction perpendicular to the wall is shown in Figure 3.16 from case 2. The Smagorinsky model and the A D M model are used in this case. As the value of the ratio is less than 0.5 one can conclude that the effect of Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 36 different subgrid models is not very important for fully developed channel flow. It can be observed that the ratio of A D M model eddy viscosity to the molecular viscosity is not symmetric at ylH = 0.5. Numerical errors are most likely reponsible for the lack of simmetry observed. The filter size near the wall region is small. The smallest value of y+ is less than 1. The wall can be resolved well with the small mesh size near the wall as shown in Figures 3.6-3.14 except for the coarse mesh case. On the other hand, the filter size in the channel middle region is large. The value of y + can be 40 times larger than the smallest one. It is not surprising to find that the fluctuations are not well resolved in the core flow. The reason is that a streching mesh is used in the direction perpendicular to the wall, which increases the mesh size in the middle of the channel. More grid points and coarse meshes near the wall can be used with the A D M model if the flow is expected to be well resolved in the middle region of the channel. It is known that the A D M model can alleviate the mesh restriction near the wall region from Stolz et al. (2001). Another alternative is using the wall model to remove the mesh restriction near the wall region (Sagaut 2001). But the wall model contains an additional source of error because a model is used for part of the flow dynamics. Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 37 Figure 3.7 Dimensionless rms fluctuation of u: yjuu luT with the ADM model and the Smagorinsky model. Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 38 Figure 3.8 Dimensionless rms fluctuation of v: V v V luT with the ADM model and the Smagorinsky model. Figure 3.9 Dimensionless rms fluctuation of w: -Jw'w' I uT with the ADM model and the Smagorinsky model. Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 64*64*64 ADM _-i I \u00E2\u0080\u00A2 i \u00E2\u0080\u00A2 i I i \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 I i \u00E2\u0080\u00A2 i \u00E2\u0080\u00A2 I \u00E2\u0080\u00A2 > \u00E2\u0080\u00A2 . I \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 . \u00E2\u0080\u00A2 I \u00E2\u0080\u00A2 . \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 I . . \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 I . \u00E2\u0080\u00A2 . \u00E2\u0080\u00A2 I 0 50 100 150 200 250 300 350 400 Figure 3.10 Dimensionless Reynolds stress: u'vlu] with the ADM model and the Smagorinsky model. 4 Figure 3.11 Dimensionless mean streamwise velocity profiles from the case 2. Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow Smagorinsky Dynamic model ADM model N=5 ADM model N = 10 \u00C2\u00BB\u00C2\u00AB\u00C2\u00AB\u00C2\u00BB \u00E2\u0080\u0094 \u00E2\u0080\u0094 \u00E2\u0080\u0094 \u00E2\u0080\u0094 DNS MKM Figure 3.13 Dimensionless rms fluctuation of v: V v V luT from the case 2. Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 41 1.4 1.2 tl I _ _ _ _ Smagorinsky Dynamic model \u00E2\u0080\u00A2 \u00E2\u0080\u00A2\u00E2\u0080\u00A2\u00E2\u0080\u00A2\u00E2\u0080\u00A2\u00E2\u0080\u00A2\u00E2\u0080\u00A2\u00E2\u0080\u00A2\u00E2\u0080\u00A2\u00E2\u0080\u00A2 ADM model N=5 ADM model N=10 \u00E2\u0080\u0094 \u00E2\u0080\u0094 \u00E2\u0080\u0094 \u00E2\u0080\u0094 DNS MKM ) I \u00E2\u0080\u0094 I 1 1 1 1 I I I I I I I I I I I I I I 0 100 200 300 400 y+ Figure 3.14 Dimensionless rms fluctuation of w: yjww I uT from the case 2. Figure 3.15 SGS model coefficients with the medium mesh case (case 2). Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 42 0.6 r 0 4 0.2 \u00E2\u0080\u00A2 Smagorinsky model O ADM model \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 o \u00E2\u0080\u00A2 8 o o \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 o \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 p \u00C2\u00B0 o o 3\u00C2\u00A9 O n o n o D 0.2 0.4 0.6 y/H 0.8 Figure 3.16 Ratio of SGS eddy viscosity to molecular viscosity: VTlv with the Smagorinsky model and the ADM model (in case 2). 3.4.2 Results for Reynolds number Re r= 1655 The parameters of the second test case have the following values: H = 0.075 m, Uh=0.204m/s and v = 0.2441xl0\"6 m2/s, Where Ub is the bulk velocity. The Reynolds number based on friction velocity is R e r = ^ l 6 5 5 v where uT = ^ Tw I p is the friction velocity , Tw is the wall shear stress and S is half of the channel height H. The geometry is the same as that with Re r = 395 . The initial distribution comes from the previous Re r = 395 case. The grid parameters are shown in Table 3.3. Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 43 Case Nx Nz Ny Ax + Az + Av + /min Ay + /max 1 96 96 96 101.1 51.6 2.1 89.4 2 96 96 128 101.1 51.6 1.2 51.6 Table 3.3 Grid parameters for Re T= 1655 The calculations are conducted with two mesh sizes: 96x96x96 and 96x128x96. The dimensionless mean velocity profiles are shown in Figure 3.17. There is no DNS data available at this high Reynolds problem. Spalding's (1961) formula quoted below, which covers the entire wall-region, is used as the reference profile, with K = 0.40, B = 5.5 : , K ) 2 K ) 3 ~ 2 6 It is shown that the A D M model provides an obvious improvement for the mean velocity profiles. The reason is that the A D M model uses an explicit filter and recovers a part of small scales from the large scales directly. The dimensionless resolvable turbulence intensities and Reynolds stress are presented in Figures 3.18-3.20. WW (1989) is the reference experimental data measured by Wei and Willmarth (1989) at Re r = 1655. The results from the Smagorinsky model are a little better than those from the A D M model. From above results, only the resolved scales quantities are shown. If the small scale is added, the computational results are closer to the experimental data, but the difference is less than 10%. The dimensionless resolvable turbulence intensities plus the small scale part (from the model) with A D M model are shown in Figure 3.21-3.22. The difference between numerical results and the experimental data can be seen clearly above y* >200 region. One reason is that the nonuniform mesh from (3.49) is used, where the mesh y+ = u + +e Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 44 ratio is large away from the wall region. It is expected to get better solutions with more finer mesh or a higher accuracy scheme. It is noted that the shear velocity is used as a non dimensional parameter. Since uT varies between the different cases and its vaue is very small, using the shear velocity as a non dimensional parameter somewhat exaggerates the differences between different cases. Figure 3.17 Dimensionless mean streamwise velocity profiles at Re r = 1655. Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow Figure 3.19 Dimensionless rms fluctuation of w.yfvV luT at Re r = 1655. Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 46 G -1 I I I 1 1 1 1\u00E2\u0080\u0094I I I I L _ l I i I i i 0 500 1000 1500 y+ Figure 3.20 Dimensionless Reynolds stress: u'v lu2T at Re r = 1655. 0 1 \u00E2\u0080\u0094 i \u00E2\u0080\u0094 i \u00E2\u0080\u0094 i \u00E2\u0080\u0094 i \u00E2\u0080\u0094 I \u00E2\u0080\u0094 i \u00E2\u0080\u0094 i \u00E2\u0080\u0094 i \u00E2\u0080\u0094 i I i i i \u00E2\u0080\u0094 J I i\u00E2\u0080\u0094 0 500 1000 1500 y+ Figure 3.21 Dimensionless rms fluctuation of u: yju'u luT with 96x128x96 mesh. Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 47 Figure 3.22 Dimensionless rms fluctuation of v: yjv'v lux with 96x128x96 mesh. 3.5 Summary The Navier-Stokes equations were solved using a staggered finite volume method for a generalized curvilinear coordinate system. The fractional step method was used to solve the momentum equations. A Poisson equation for pressure was obtained to satisfy the continuity equation. For the channel flow, the fast Fourier transform and the cyclic reduction method were presented to solve the Poisson equation. For parallel computation of a large eddy simulation, the domain decomposition approach is used to partition the data. Each sub-domain data is assigned to one processor to calculate the quantities and exchange the boundary data with neighbors. The parallel efficiency is still higher than 80% f o r t n e p a r a i i e l computation of a large eddy simulation. Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 48 The channel flow at Re T = 395 was modelled first with three different SGS models utilizing three alternative mesh sizes. The coarse mesh case did not predict the logarithmic layer well, which can be arributed to an inadequate grid resolution. The fine mesh case was close to the reference DNS result. The medium case was close to the fine mesh case. The Smagorinsky SGS model shows little difference from results obtained with the dynamic SGS model. The peak value of the dynamic coefficient was 10% higher than the Smagorinsky coefficient. In most of the computational region, the difference between the dynamic model coefficient and the Smagorinsky coefficient is less than 2%. The A D M model improves the mean velocity profile compared to the Smagorinsky model. The reason is that the A D M model uses an explicit filter and recovers a part of small scales from the large scales directly. The resolved turbulence intensities with the A D M model are close to the Smagorinsky model results. At the higher Reynolds number R e r = 1655, the channel flow was computed by using two mesh sizes. The trend in the calculations was similar to that of the low Reynolds number case. The Smagorinsky model with 64x64x64 mesh can provide a good representation of channel flow and is used to provide inflow conditions for the modelling of the converging section as described in chapter 4. The A D M model can improve the mean velocity profile using the explicit filter. But the A D M model needs much more computational time than the Smagorinsky model, and so it will be used only for model comparisons in chapter 4 and chapter 6. The objectives 1 and 2, described in chapter 1, are achieved. Chapter 4 LES for Turbulent Flow through a Converging Section 4.1 Introduction As described in chapter 1, the rapidly converging section of a paper-machine headbox carries a dilute concentration of pulp fibers in water to the wire mesh where the fibers are dried to become paper. Ideally, the mean velocity distribution in the stock (fibers and other substances in fluid suspension) leaving the converging section (or slice) should be uniform over the paper thickness direction and across the entire span of the converging section exit. Otherwise, the non-uniform stock can result in defects in the paper being produced by the machine. To achieve good paper formation, there must be no fiber flocculation in the jet exiting the headbox. The turbulence produced in the converging section does not affect the mean flow distribution significantly, but it is critically important in preventing unwanted fiber flocculation and in providing a degree of dispersion and preferential orientation of the fibers. A detailed understanding of this turbulence is therefore essential in order to study the fiber motion and to predict the quality of the paper produced. Shariati (2002) has studied the flow in the converging section of a headbox both experimentally and numerically. The experimental headbox that he used is a laboratory scale model of a typical headbox with the size reduced by a factor of 5. Shariati made L D V (Laser Doppler Velocimeter) measurements of mean and fluctuating components of the velocity in the headbox. He also compared results using the k \u00E2\u0080\u0094 \u00C2\u00A3 model and a Reynolds stress model (RSM) to his experimental measurements. He found that the turbulence kinetic energy predicted with the k \u00E2\u0080\u0094 \u00C2\u00A3 model was much too high near the slice exit. The R S M model provided better agreement with the 49 Chapter 4 LES for Turbulent Flow through a Converging Section 50 experimental measurements. However, the RSM model cannot provide the instantaneous velocity distribution which is required to predict the fiber motion in the converging section. In this chapter, a large eddy simulation computation of the converging section is presented. The geometry of the converging section is the same as that of the experimental headbox used by Shariati (2002). The calculated time-averaged turbulence components predicted by the LES are compared to Shariati's measured values. Different inflow conditions are explored. Qualitatively, the calculated and observed turbulence distributions follow similar trends. More exact inflow conditions for the turbulence must be introduced into the calculated LES results if better agreement is to be obtained. This aspect of the present work needs further investigation. The results of flow simulation are presented in section 4.2. Some discussions are presented in section 4.3. A summary is given in section 4.4. 4.2 Flow simulation 4.2.1 Experimental setup The experiments were performed by Shariati (2002) at the University of British Columbia. The experimental set-up used a closed flow system, diagrammatically shown in Figure 4.1. Chapter 4 LES for Turbulent Flow through a Converging Section 51 Figure 4.2 Experimental headbox drawing (from Shariati 2002). The headbox drawing is shown in Figure 4.2. In the flow loop, water is pumped from the reservoir tank to the headbox through the pipes and rectifier tubes. The rectifier tubes (or diffusers) Chapter 4 LES for Turbulent Flow through a Converging Section 52 are round at the inlet and rectangular at the outlet with slowly increasing cross sectional areas. They are used to provide turbulence energy for the flow and to generate a fairly uniform velocity profile at the converging section inlet. In this headbox model, there are 40 rectifier tubes, two rows in the vertical or height direction and 20 each in the spanwise or horizontal direction. The headbox section starts with a rectangular channel, which remains constant in the cross section area until the length reaches 0.075 m. Table 4.1 shows the geometry of the converging section. By means of LDV, Shariati (2002) measured the mean velocity and the rms velocity along the centerline of the converging section. Parameters Values (m) Width 0.75 (constant) Inlet height 0.075 Slice height 0.0075 Contraction ratio 10 Length 0.3 Table 4.1 The geometry of the converging section. 4.2.2 Numerical simulation of the converging section The numerical method is the same as that of chapter 3 for fully developed channel flow. The Poisson equation is solved with FFT and CGSTAB method. The subgrid model uses the constant Smagorinsky model or the ADM model. Al l the calculations are performed with parallel computation, as described in chapter 3, on the P4-Xeon Cluster with a Myrinet network at the University of British Columbia. The computational domain is shown in Figure 4.3. In the z direction, no-slip wall boundary conditions are imposed along the solid walls. The spanwise direction, y, is treated as periodic. At the inlet plane the unsteady distributions for the velocity components are generated by a separate LES of Chapter 4 LES for Turbulent Flow through a Converging Section 53 fully developed channel flow. Different inflow conditions that simulate the turbulence produced by the rectifier tubes are discussed later. At the exit plane, one possibility is to set all derivatives in the direction normal to the boundary to zero. This condition is often used in steady flow. For unsteady flow, the convective boundary condition (Kaltenbach 1998, Kaltenbach et al. 1999) is employed. Uc is set to the mean streamwise velocity at the exit plane. This condition (4.1) eliminates most of the problems caused by pressure waves reflecting from the outflow condition. The dimensions of the computational domains are the same as the dimensions of the experimental apparatus shown in Figure 4.2 except for the width in the spanwise direction. The computational domain spanwise width is approximately three inlet heights. The 64x64x64 mesh is used. The meshes are uniform in the streamwise direction and the spanwise direction. A non-uniform grid is used in the normal wall direction to resolve the near wall region as in chapter 3 (formula 3.49). The measured and calculated mean streamwise velocities along the centerline of the converging section (line ABC in Figure 4.4) are shown in Figure 4.5. The velocity is non-dimensionalized with the bulk velocity at the inlet of the converging section, U* = UIUb. Here X* = xlLx is the x position normalized with the total length of the converging section Lx. The computed values of the mean velocity agree well with the measured values. The resolved rms fluctuations u,v and w normalized by u at the inlet (ujn) are presented in Figures 4.6 - 4.8. The rms fluctuation of the velocity in the x direction, u diminishes towards the exit, while w in the z direction increases. The trend of the calculated results from u and w is the same as that of the experimental results, but there are some differences between them, particularly in the region between the inlet of the converging section and the first third of the converging section. The computed v increases along the center line, while the measured v decreases in the same region. These Chapter 4 LES for Turbulent Flow through a Converging Section 54 differences can be attributed to the difference between the measured and assumed inflow conditions. Also the expermental flows settings have some errors that produce secondary flows. From the experiment, the two rows of tubes can generate much higher turbulence than the single channel which is used as the model in above calculations. It can also be observed that the computational inlet turbulence intensity is less than one half of the measured value at 0.03m from the inlet plane (DE line in Figure 4.4) in Figure 4.10. Figure 4.9 shows the turbulence kinetic energy k nondimensionalized by k at the inlet of the position 0.006m after the convergence starts. The turbulent kinetic energy predicted by LES and the standard k \u00E2\u0080\u0094 \u00C2\u00A3 model in most of regions is close to the measured values, but the computed values with the standard k \u00E2\u0080\u0094 \u00C2\u00A3 model increase rapidly towards the exit. The reason is that the standard k \u00E2\u0080\u0094 \u00C2\u00A3 model includes a turbulence generation term that depends on Wldx. That term becomes very large close to the exit of the converging section and causes a break down of the model. Also, the standard k \u00E2\u0080\u0094 \u00C2\u00A3 often fails to predict the turbulence kinetic energy where the turbulence is significantly anisotropic such as the case in the region close to the exit of the converging section. To consider the effect of the different inflow conditions, a \"two-channel\" inflow condition is constructed. The inflow condition is linearly interpolated from two half-inlet heights of the converging section channels where one channel is located above the other channel (Figure 4.11). Each channel is then 0.037 m high and the wall thickness is 0.001 m, the total of the two channels and wall thickness making up the 0.075m inlet height of the headbox. The wall thickness is explicitly considered in the generation of the inflow condition. The comparison of the mean streamwise velocity along the centerline of the converging section is shown in Figure 4.12 with the \"two-channel\" inflow conditions. There is no significant difference between the two velocity variations. Figures 4.13-4.15 present the comparison of rms fluctuations w',v'and w normalized by w'at the inlet (uin) with two inflow conditions. The results of the \"two-channel\" inflow conditions are better Chapter 4 LES for Turbulent Flow through a Converging Section 55 than that of the one channel inflow condition, matching the trend of the measured values from u and w but not reproducing them with real accuracy. There is no improvement with \"two-channel\" inflow conditions for v . This appears to be related to the impossiblity to set accurately the flow in the experiment and to the inflow data used in the simulations. The inflow turbulence condition clearly affects the numerical result and better agreement could probably be obtained i f the turbulence at the inflow section could be more precisely modelled in the calculations. 0.25 -0.2 \" Zo.15 \" 0.1 -X Figure 4.3 Computational domain for the converging section. Only a subset of the actual grid is plotted. Chapter 4 LES for Turbulent Flow through a Converging Section 56 E D Figure 4.4 Measurement locations along the centerline of the converging section (line ABC) and the vertical line at 0.03m after inlet (line DE). 11 10 9 8 7 F> \u00E2\u0080\u00A2 L E S r~\ U* measurement 6 \" 5 \u00E2\u0080\u00A2 o o \u00E2\u0080\u00A2 4r> O n O n J I I I I I I ! I I I I I I I 0.25 0.5 0.75 X* (x/Lx) Figure 4.5 Comparison of the mean velocity at the centerline of the converging section. Chapter 4 LES for Turbulent Flow through a Converging Section 1.2 1.1 1 0.9 0.8 3 0.6 = 0.5 0.4 0.3 0.2 0.1 O \u00E2\u0080\u00A2 L E S \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 r o u'/u' in measurement o \u00E2\u0080\u00A2 o \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 o \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 o \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 o 0 p n DD o O \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 01\u00E2\u0080\u0094I_ I \u00E2\u0080\u0094 \u00E2\u0080\u0094 I i i i i I 1 1 1 > 1 1 1 1 1 1 0 0.25 0.5 0.75 1 X* (x/Lx) Figure 4.6 Comparison of rms value of the velocity fluctuation u normalized by 1.2 1.1 1 0.9 0.8 0.7 3 0.6 5 0.5 0.4 0.3 \u00E2\u0080\u00A2o L E S O v7u' in measurement o o o o o o \u00C2\u00B0 o 0.2 \u00C2\u00A3 ] \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 0.1 r \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 . ' 1 1 1 1 \u00E2\u0080\u00A2 1 ' 1 0.25 0.5 0.75 X* (x/Lx) J I L _ L . Figure 4.7 Comparison of rms value of the velocity fluctuation v normalized by u'n Chapter 4 LES for Turbulent Flow through a Converging Section 1.2 1.1 1 0.9 0.8 c 0 7 3 0.6 * 0 . 5 0.4 0 0.2 0.1 \u00E2\u0080\u00A2 LES o O W/u' in measurement coo \u00E2\u0080\u00A2 ; o o o \u00C2\u00B0 Qa ,oo \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 J I I I I I I I I L _L 0.25 0.5 0.75 X* (x/Lx) J i i i_ Figure 4.8 Comparison of rms value of the velocity fluctuation w normalized by u. 7 | -2 h \u00E2\u0080\u00A2 O A LES k measurement A A A A O O O O A - i i i 1 \u00E2\u0080\u00A2 1 i i_ j I i u 0.25 0.5 X* (x/Lx) 0.75 Figure 4.9 Comparison of turbulence kinetic energy k nondimensionalized by k at the inlet of the position 0.006m after the convergence starts. Chapter 4 LES for Turbulent Flow through a Converging Section 1 0.9 0.8 f 0.7 0.6 Z* 0.5 0.4 h 0.3 f 0.2 0.1 t 0 u' for one channel inflow condition O u ' measurement at 3cm from inlet plane \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 a. \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 M B - 1 JL \u00E2\u0080\u00A2 0.01 0.02 \u00E2\u0080\u00A2\u00E2\u0080\u00A2 I \u00E2\u0080\u00A2,\u00C2\u00A3.\u00E2\u0080\u00A2 0 o o o o o o Q O O O o o o o o o o o c 1 1 1 1 1 1 1 I I I 0.03 0.04 0.05 u' (m/s) 0.06 0.07 0.08 Figure 4.10 Comparison of the fluctuating velocity u at inlet for one channel inflow condition with the measured value at the vertical line (DE line in Figure 4.4) at the location 0.03m from the inlet. H/2 H/2 interolation t 0.001m interolation Figure 4.11 The diagram of\" two-channel \" inflow condition. Chapter 4 LES for Turbulent Flow through a Converging Section 60 11 10 11111 9 8 7 6 5 4 -3 2 > 1111 1\ m n 1 \u00E2\u0080\u00A2 o i> LES with one channel inflow condition U* measurement LES with two channels inflow condition D O o s o o o O Q O p P q \u00E2\u0080\u00A2 1 1 1 1 J I I L 0.25 0.5 X* (x/Lx) 0.75 Figure 4.12 Comparison of the mean velocity at the centerline of the converging section with different inflow conditions. A A L \u00E2\u0080\u0094 j=j - LES with two channels Inflow condition \u00C2\u00B0 6 n ^ O a A J = \" \u00E2\u0080\u00A2 LES with one channel Inflow condition ^\u00E2\u0080\u00A2\u00E2\u0080\u00A2\u00E2\u0080\u00A2\u00E2\u0080\u00A2\u00E2\u0080\u00A2Q A O uVu'In measurement \u00E2\u0080\u00A2 n A A 0.9\" 0.8 . 0.7 i 0.6 ' 0.5 0.4 0.3 0.2 0.1 O o A O O o Q Q A r A A 0.25 0.5 0.75 X* (x/Lx) Figure 4.13 Comparison of rms value of the velocity fluctuation u normalized by u'n with different inflow conditions. Chapter 4 LES for Turbulent Flow through a Converging Section 61 1.2 1.1 1 0.9 0.8 0.7 c \"a 0.6 5 0.5 0.4 0.3 0.2 0.1 \u00E2\u0080\u00A2o O c O v o LES with one channel Inflow condition v7u' in measurement LES with two channels Inflow condition Q O n V ? ) V V v V V, V V V V V \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 n D n D o 0 0.2 0.4 0.6 X* (x/Lx) 0.8 Figure 4.14 Comparison of rms value of the velocity fluctuation V normalized by u'n with different inflow conditions. 1.2 1.1 1 0.9 0.8 0.7 c \"= 0.6 * 0 . 5 0.4 0.2 0.1 O LES with one channel inflow condition Q O w'/u' in measurement Q <] LES with two channels inflow condition IPO \u00C2\u00B0 0 3l/3.22. (5.1) 3 2.5 L E S : 2mm fiber \u00E2\u0080\u0094 \u00E2\u0080\u0094 \u00E2\u0080\u0094 - L E S : 3mm fiber L E S : 1mm fiber > O l s o n ' s fitting curve for his exper iment data O O o 2 c o 1.5 o o 2 c s c 0.5 2 3 4 C h a n n e l height z / L 5 6 Figuer 5.3 Concentration for 0.001m, 0.002m, 0.003m fiber vs. channel height scaled by fiber length with LES Chapter 5 Fiber Concentration and Orientation 85 5.3 Fiber orientation in the headbox In addition to LDV measurements, fiber experiments had been conducted by Zhang (2001) who measured the orientation of dyed nylon fibers moving in the Shariati's headbox (Shariati 2002). The fibers are made of nylon and have a nominal length of 0.003m and the diameter of 44 pm. The suspension is well within the dilute regime, with the consistency of no more than 0.001%, which means there is little interaction between fibers. In the experiment, video pictures were taken of fibers in motion at several locations along the center of the converging section. 5.3.1 Initial conditions The fiber's position can be defined by three variables: fiber center (xc,yc,zc), polar angle 0 < 6 < 7t, which is the angle between the fiber main axis and Z-axis, and the azimuthal angle 0 <

(\u00C2\u00AB;.)Aa;.) and third moment j/^af p{ai )Aar() of the probability distribution in x-z plane. Chapter 5 Fiber Concentration and Orientation 95 First Second Third First Second Third First Second Third moment moment moment moment moment moment moment moment moment (experiment) (experiment) (experiment) (LES) (LES) (LES) (k-e) (k-e) (k-e) 0.045m -0.001 0.794 0.020 0.052 0.753 0.045 0.122m -0.005 0.776 0.012 0.027 0.689 0.046 0.157m 0.032 0.737 0.029 -0.025 0.695 -0.027 0.192m 0.020 0.745 0.034 0.045 0.727 0.042 -0.011 0.571 -0.017 0.227m 0.001 0.659 0.012 0.026 0.710 0.038 -0.008 0.468 -0.014 0.262m 0.004 0.611 0.013 -0.020 0.574 -0.108 -0.008 0.350 -0.012 exit -0.007 0.325 0.001 -0.027 0.538 -0.006 0.002 0.161 -0.004 Table 5.2 The first moment (^JCCIp{oct) Acct), second moment (^afp(at) A or,) and third moment i i ( >J oc]pia^Aot,) of the probability distribution in x-y plane. 0.1 r 0.05 h c o I --0.05 h -0.1 \u00E2\u0080\u00A2 LES k-e O Experiment \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 o \u00E2\u0080\u00A2 A A o \u00E2\u0080\u00A2 c o A \u00E2\u0080\u00A2 0.05 0.1 0.15 0.2 0.25 0.3 Figure 5.14 The first moment ( ^ OCjPia,)A*3r.) at different stations in x-z plane. Chapter 5 Fiber Concentration and Orientation 0.1 0.05 c E 1 0 -0.05 \u00E2\u0080\u00A2 LES A k-e O Experiment \u00E2\u0080\u00A2 o \u00E2\u0080\u00A2 o o \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 0 o G A A A o A c \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 Q A I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 0 0.05 0.1 0.15 0.2 0.25 0.3 Figure 5.15 The first moment (]T GCip(0Cl)Aat) at different stations in x-y plane. r 0.8 | 0.6 O \u00C2\u00A3 \u00E2\u0080\u00A2o 8 0.4 0 0.2 1 1 1 1 1 1 \u00E2\u0080\u00A2 LES A k-e O Experiment \u00E2\u0080\u00A2 0 A O A * 1 1 1 1 I 1 1 1 1 * 1 1 1 1 I 1 1 1 1 * 0.05 0.1 0.15 0.2 0.25 0.3 X Figure 5.16 The second moment ( ^ a]pipe,)Aai) at different stations in x-z plane. Chapter 5 Fiber Concentration and Orientation r 0.8 c | 0.6 O E -o 8 0.4 9 (0 0.2 : \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 LES A k-e Q Experiment o \u00E2\u0080\u00A2 : \u00E2\u0080\u00A2 A 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 \u00E2\u0080\u00A2 c c \u00E2\u0080\u00A2 A \u00E2\u0080\u00A2 o 1 1 ' 1 1 1 ' 1 1 1 1 1 1 1 1 0.05 0.1 0.15 0.2 0.25 0.3 X Figure 5.17 The second moment (^T a]p{at)Aat) at different stations in x-y plane. \u00E2\u0080\u00A2 LES 0.1 | - A k-e o \u00E2\u0080\u00A2 -0.2 Experiment O A A | f n A o i - o , -i -\u00E2\u0080\u00A2 -0 31 i i i i I i ' I ' I i i i i I i ' i i I i i I i i i i I i i ' 0 0.05 0.1 0.15 0.2 0.25 0.3 X Figure 5.18 The third moment (^af p(ai)Aai) at different stations in x-z plane. Chapter 5 Fiber Concentration and Orientation 98 0.04 0.02 0 S-0.02 E o E-0.04 1 ~-0.06 r--0.08 --0.1 -0.12, \u00E2\u0080\u00A2 O \u00E2\u0080\u00A2 A o G o o \u00E2\u0080\u00A2 k-e Experiment ' C O A o o A A \u00E2\u0080\u00A2 0.05 0.1 0.15 0.2 0.25 0.3 Figure 5.19 The third moment ( ^ o f p{CCi)AflT() at different stations in x-y plane. 5.3.3 Discussion In the above calculations, 8,000 fibers are used for the fiber orientation distribution calculation. In order to verify that 8000 fibers represent the fibers with sufficient accuracy, the use of 16,000 and 24,000 fibers are also tested. Figure 5.20 shows a comparison of fiber orientation distribution at x = 0.157m in x-y plane with 8000, 16,000 and 24,000 fibers where the error bars stand for the standard deviation from the mean of these three sets of data. The result shows that the three sets of data produce no significant difference in the x-y plane. Consequently 8,000 fibers are used for the calculation. Chapter 5 Fiber Concentration and Orientation 99 Figure 5.20 Fiber orientation distribution at x = 0.157m in x-y plane with 8000, 16,000 and 24,000 fibers 5.4 Summary Using the LES method with different length fibers, the concentration curves of the 0.002m and 0.003m fibers are similar to the experiments of Olson (1996). The LES and a fiber motion model have been coupled for predicting the orientation of rigid fibers in dilute suspensions. The effects of turbulence that would tend to randomize the fibers are included in these simulations of fiber motion. Random initial fiber orientations are set at the inlet of the channel. The statistical expressions of the orientation of a large number of fibers are evaluated at each station by computing the orientation of each single fiber along the central streamline. The LES results predict the fiber orientation reasonably well when compared to the experimental data. The Chapter 5 Fiber Concentration and Orientation 100 predicted orientation distributions with LES are better than the k-e simulated fiber orientations near the exit of the converging section. The numerical results show that the LES calculation scheme can provide a useful method for the simulation of the interaction between fibers and complex fluid flow patterns. The LES code is coupled with a suitable fiber model to study the statistical orientation of nylon \"fibers\" in the converging section. The objective 4, described in chapter 1, is achieved. Chapter 6 Simulation of a Typical Industrial Headbox 6.1 Introduction A numerical tool to simulate turbulence in a headbox has been developed in chapter 3, where its use for simulating the flow through the converging section of a laboratory headbox has been described. However, the Reynolds number for the laboratory headbox is much lower than that for an industrial headbox, and it is therefore still challenging to simulate the turbulence in a typical industrial case. As far as the LES method used for the present simulations is concerned, the higher Reynolds number should be an advantage, because the spectral separation between inertial and dissipative scales, being greater in the high Reynolds number case, allows for very clear and distinct simulation of the larger scales. In this chapter, the flow through the converging section of a typical generic headbox is simulated with LES. Two inflow conditions are investigated in section 6.2. Some observations are presented in section 6.3. 6.2 Flow simulation 6.2.1 Computational geometry The geometry of the converging section is shown in Figure 6.1. The length of the converging section in the streamwise or x direction is 0.665 m . The inlet height of the converging section is 0.105 m . The contraction ratio is 7 . The actual length in the cross-stream or y direction is 8.64 m. In the numerical simulation presented here, the dimension in the y direction is set to 0.1575 m because periodic boundary conditions are used in the y direction, assuming that the flow 101 Chapter 6 Simulation of a Typical Industrial Headbox is essentially two-dimensional. 102 Figure 6.1 The geometry of the converging section. 6.2.2 Inflow conditions In the z direction, which is normal to the walls at entry, no-slip wall boundary conditions are imposed along the solid walls. As noted above, the spanwise direction, y , is treated as periodic. At the exit plane, the convective boundary condition (Kaltenbach 1998, Kaltenbach et al.1999) is imposed. At the inlet plane, two inflow conditions are tested. One is to use the unsteady distribution for the velocity components generated by a separate LES of fully developed channel flow (Figure 6.2). The other is a \"five-channel\" inflow condition, simulating the flow from the outlets of a series of five equal smaller channels, typical of the turbulence generators often used in commercial headboxes. The inflow condition is linearly interpolated from five small height channels (Figure Chapter 6 Simulation of a Typical Industrial Headbox 103 6.3). Each channel is then 0.019 m high and each wall separating the five channels is 0.0025 m thick, the total of the five channels and four intervening walls making up the 0.105 m inlet height of the headbox. The wall thickness is explicitly considered in the generation of the inflow condition. The mean velocity profiles at the inlet for the one channel inflow and \"five-channel\" are shown in Figure 6.4. Z* is the z position from the bottom plane normalized with the total height at the same location. The corresponding resolved rms fluctuations u , v and w at the inlet are presented in Figures 6.5-6.7. From these graphs, some zero values for the \"five-channel\" inflow condition can be seen due to the wall thickness. The volume flow rate is 0.208 m3 / s (per unit width) which is close to the maximum volume flow rate expected for a commercial headbox. The same mass flow rate is imposed for the two inflow conditions. As already noted, the wall thickness between small channels for the \"five-channel\" inflow condition is 0.0025 m and each small channel height is 0.019 m . The numerical method is the same as that of chapter 4. The ADM model and the constant Smagorinsky model are used to consider the effect of subgrid scales. All the calculations are performed with parallel computation, as described in chapter 3, on the P4-Xeon Cluster with a Myrinet network at the University of British Columbia. Chapter 6 Simulation of a Typical Industrial Headbox Figure 6.2 One-channel inflow condition for the converging section. Chapter 6 Simulation of a Typical Industrial Headbox Figure 6.3 Schematic diagram of details of one of the five identical channels the \"five-channel\" inflow condition showing the Chapter 6 Simulation of a Typical Industrial Headbox 1 U (m/s) Figure 6.4 Comparison of the mean velocity at the inlet. Figure 6.5 Comparison of rms value of the velocity fluctuation u at the inlet. Chapter 6 Simulation of a Typical Industrial Headbox 107 w\" (m/s) Figure 6.7 Comparison of rms value of the velocity fluctuation w'at the inlet. Chapter 6 Simulation of a Typical Industrial Headbox 108 6.2.3 Numerical results The flow characteristics are taken at several positions along the converging section. The first position is located at one sixth of the full length where the converging section starts (AB line in Figure 6.8). The second position is at one third of the whole length (CD line in Figure 6.8). The third position is at two thirds of the whole length (EF line in Figure 6.8). The last position is at the outlet. The 96x96x128 mesh is used. The meshes are uniform in the streamwise direction x and the spanwise direction y . A non-uniform grid is employed in the normal wall direction z to resolve the near wall region as described in chapter 3. The mean streamwise velocity profiles computed from two inflow conditions are shown in Figure 6.9 at the first station AB with the ADM model. They are very similar even at this first station. The resolved rms fluctuations u , v and w are presented in Figures 6.10-6.12 at the first station AB with the ADM model. From Figures 6.5-6.7, the fluctuations from the two inflow conditions have a significant difference at the inlet. The differences from two inflow conditions are reduced a great deal at the first station from Figures 6.10-6.12. The fluctuations from \"five-channel\" inflow conditions grow more quickly, presumably due to the energy produced by the more non-uniform inlet mean velocity profile that is present in the \"five-channel\" case. The mean streamwise velocity profiles from the two inflow conditions are shown in Figure 6.13 at the second station CD derived from use of the ADM model. Not surprisingly, there is little difference between these two distributions. The resolved rms fluctuations u , v and w are presented in Figures 6.14-6.16 at the second station CD with the ADM model. The fluctuations from \"five-channel\" inflow conditions are slightly greater than those from the one channel inflow conditions at this second station, showing that more turbulence is generated at one third of the whole length from the \"five-channel\" inflow conditions. Chapter 6 Simulation of a Typical Industrial Headbox 109 The mean streamwise velocity profiles from the two inflow conditions are presented in Figure 6.17 at the third station EF. The resolved rms fluctuations u , v and w are shown in Figures 6.18-6.20 at the third station EF. The trend is the same as that at the second station. The mean streamwise velocity profiles from the two inflow conditions are shown in Figure 6.21 at the outlet. Essentially identical mean velocity profiles are obtained at the outlet from the two different inflow conditions. The resolved rms fluctuations u , v and w are presented in Figure 6.22-6.24 at the outlet. The fluctuations u from the two inflow conditions show little difference here. However, the fluctuations v and w from the \"five-channel\" inflow conditions are greater than that from one channel inflow conditions at the outlet, implying that more turbulence is indeed generated using the \"five-channel\" inflow conditions, and that this increased turbulence lasts to the outlet of the headbox. As an alternative to the results produced by the ADM model, some simulations are repeated with the Smagorinsky model simulating the small scales of turbulence. The mean streamwise velocity profiles from the two inflow conditions are shown in Figure 6.25 at the outlet with the constant Smagorinsky model. The mean velocity profile is not symmetric near the upper wall region. This can be a result of the inflow condition that was interpolated from the asymmetrical channel flow. The resolved rms fluctuations u , v and w at the outlet are presented in Figures 6.26-6.28. The fluctuations u from the two inflow conditions have little difference at the outlet. However, as before at the outlet, the fluctuations v and w from \"five-channel\" inflow conditions are greater than that from the one channel inflow conditions, showing that essentially similar results are obtained when using either model of small-scale turbulence. The ratios of resolved rms fluctuations v I u and \v I u are presented in Figures 6.29-6.30 at different positions with the ADM model. From these graphs, the turbulence is close to isotropic at the first and second stations in most regions. The anisotropy of the turbulence grows at the third Chapter 6 Simulation of a Typical Industrial Headbox 110 position. At the outlet, the turbulence is very anisotropic in most regions. The maximum value of v7 u can be up to 3.5 at the outlet of the converging section. Figures 6.31-6.32 show the ratio of v I u and w'/w'at the outlet with the two inflow conditions. There is little difference in these plots for the two inflow conditions. B D inlet A C E outlet Figure 6.8 Comparison positions along the headbox. U (m/s) Chapter 6 Simulation of a Typical Industrial Headbox Figure 6.9 Comparison of the mean velocity at the first station AB. 0.8 0.6 N 0.4h 0.2 h n o D \u00C2\u00A9 D D o \u00E2\u0080\u00A2 One channel inflow ADM Five channels Inflow ADM c?D Fci f t^P. 0.05 0.1 0.15 0.2 0.25 u\" (m/s) Figure 6.10 Comparison of rms value of the velocity fluctuation u at the first station AB. 1 r-0.8 0.6 N 0.4 0.2 CP O \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 C K ^ fj One channel inflow ADM O Five channels inflow ADM \u00E2\u0080\u00A2 o I i?^Q*A' i l l rr I I 0.02 0 X 4 0 . 0 6 0.08 ^ 0.1 0.12 \u00E2\u0080\u00A2 (m/s) Figure 6.11 Comparison of rms value of the velocity fluctuation v at the first station AB. Chapter 6 Simulation of a Typical Industrial Headbox u.uu w'(m/s) Figure 6.12 Comparison of rms value of the velocity fluctuation w at the first station AJ3. Figure 6.13 Comparison of the mean velocity profiles at the second station CD. Chapter 6 Simulation of a Typical Industrial Headbox 113 0.8 0.6 N 0.4 r-0.2 h o J0 \u00E2\u0080\u00A2 \u00E2\u0080\u00A2 o One channel Inflow ADM Five channels inflow ADM -Qo Q3 0.02 0.04 0.06 0.08 0.1 u\" (m/s) Figure 6.14 Comparison of rms value of the velocity fluctuation u at the second station CD. 1 r-0.8 0.6 M 0.4 0.2 \u00E2\u0080\u00A2 One channel inflow ADM O Five channels Inflow ADM o i i i < i I Q i i i i i / e do k = k + \ Pk=sklh-i Pk =zk+PkPk-\ <*k =sk Kpk-APk) uk=uk^+akpk rk=rk-\~akAPk end while The CGSTAB algorithm is given follows: Initialize by setting: k = 0, u0=uin,r0=Q-Auin, p0 = 0, = 0, s \u00C2\u00ABo = Po = h =1 \u00E2\u0080\u00A2 while || rk \\> e do k = k + l Pk=ro-rk-x 139 Appendix <\u00C2\u00B0k = (Aft- iV(a*-iA-i) Pk = uk-\ + ak (Pk-i ~ ak-\A-\) solve the system: Mz = k=Az Yk =PkWk-r0) w = rk-x-rA i f || w | | 2 < \u00C2\u00A3 rk=rk-x-Ykz return end i f solve the system: My = w v = Ay ak={v-rk)l{vv) uk =uk-i+rkz + aky rk=w-akv end while Appendix 141 Appendix B Incomplete LU Decomposition The preconditioner M = LU should be a good approximation to A. In Incomplete L U (ILU), the L and U matrices have non-zero elements only on diagonals on which A has non-zero elements. Every element on the main diagonal of U is set to unity to make the matrices unique. The L and U have the following schematic presentation: where LLW = A'W , LLS = ALS L' -A'-L' U'~NJ -L'U'-] UN - AN I LP, UE - AE I LP. Appendix 142 The meaning of one-dimensional storage location index / and Nj can be found in Table 3.1 of Ferziger and Peric (1999). Appendix C Random Fluctuation Inflow Generation Method A desired mean flow and Reynolds stress tensor are matched for the random fluctuation method. Let the target mean flow profile be c/(v) and Reynolds stress tensor by Ry(v) =< u'jUj > z , , where the operation < > z J is an average on the spanwise direction and time. Three random numbers ui (z), i = 1,2,3 are generated. They satisfy the following condition: each sequence has zero mean and unit variance. The velocity field is constructed according to u,(y,z) = U,(y) + 0y Sj(z) where Pu = V^ 22 ~ fill ' fill = Al P32 ~ ^ > A>33 - V A 3 3 ~P 31 -Pl2 \u00E2\u0080\u00A2 By not listed above are set zero. "@en . "Thesis/Dissertation"@en . "2005-11"@en . "10.14288/1.0080757"@en . "eng"@en . "Mechanical Engineering"@en . "Vancouver : University of British Columbia Library"@en . "University of British Columbia"@en . "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."@en . "Graduate"@en . "Large eddy simulation of the turbulence and fiber motion in a headbox"@en . "Text"@en . "http://hdl.handle.net/2429/17095"@en .