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 P A R T I A L F U L F I L L M E N T OF THE REQUIREMENTS OF THE D E G R E E 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 © 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 ; Nomenclature vii xi Acknowledgements xiii Chapter 1 Introduction 1 Chapter 2 Literature Review 4 2.1 Headbox 2.1.1 2.1.2 2.2 4 Headbox flow simulations Fiber orientation 5 6 Numerical methods 2.2.1 2.2.2 2.2.3 7 Direct Numerical Simulation (DNS) R A N S models Large Eddy Simulation Chapter 3 8 9 10 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 3.2.2 Subgrid model 3.2.2.1 Dynamic subgrid model 3.2.2.2 Approximate deconvolution model ( A D M ) 3.3 Numerical methods 21 3.3.1 The flow solver 3.3.2 Numerical methods for solution of the Poisson equation 3.3.2.1 Direct method 3.3.2.2 Conjugate Gradient method and C G S T A B 3.3.3 Parallel computation for solution of the Poisson Equation 3.3.3.1 Parallel computation using the direct method 3.3.3.2 Parallel computation for C G S T A B method 3.3.4 Di screte test filters 3.4 Numerical results 21 26 26 26 27 28 28 30 31 3.4.1 Reynolds number R e = 395 3.4.2 Results for Reynolds number R e = 1655 3.5 13 14 16 18 32 r r Summary 42 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 4.2.2 4.2.3 4.3 Experimental setup Numerical simulation of the converging section Alternative methods for the generation of inflow data Discussion 4.3.1 4.3.2 4.4 69 Dependence on grid resol ution Effects of the different subgrid models Summary Chapter 5 50 52 62 70 74 77 Fiber Concentration and Orientation 80 5.1 Introduction 80 5.2 Fiber concentration 81 5.2.1 5.2.2 5.3 Fiber model and wall model Numerical results 81 83 Fiber orientation in the headbox 85 5.3.1 5.3.2 5.3.3 5.4 Initial conditions Numerical results Discussion 85 86 98 Summary Chapter 6 99 Simulation of a Typical Industrial Headbox 101 6.1 Introduction 101 6.2 Flow simulation 101 6.2.1 6.2.2 6.2.3 6.3 Computational geometry Inflow conditions Numerical results 101 102 108 Summary Chapter 7 122 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 Re TABLE 3 . 3 GRID PARAMETERS FOR 3 3 = 395 r Re = T 1655 4 3 T A B L E 4 . 1 T H E GEOMETRY OF THE CONVERGING SECTION TABLE 4 . 2 GRID PARAMETERS. Nx, Ny 5 2 DENOTE THE NUMBER OF MESHES IN THE STREAMWISE, ,Nz SPANWISE AND WALL-NORMAL DIRECTIONS, RESPECTIVELY. S IS HALF OF THE INLET HEIGHT OF THE CONVERGING SECTION T A B L E 5 . 1 T H E FIRST MOMENT ( ^ or,.)A« .), ( SECOND MOMENT (afp(a j )Aa i ) / i ( o f p{a )A.a ) OF THE PROBABILITY j i T A B L E 5 . 2 T H E FIRST MOMENT ( ^ (af j DISTRIBUTION IN X-Z PLANE 9 4 (X p{jDt )Aor.), SECOND MOMENT (]T a)p{a )Aor,) AND i p{a )Aa ) OF THE PROBABILITY i 7 1 AND THIRD MOMENT i t DISTRIBUTION IN X-Y PLANE VI THIRD MOMENT 9 5 List of Figures FIGURE 1.1 F I G U R E 3.1 T H E S C H E M A FOR T H E PROCESS OF APAPER M A C H I N E ( F R O M SHARIATI 2002) 1 S C H E M A T I C OF E N E R G Y S P E C T R U M ILLUSTRATING T H E RELATIONSHIP A M O N G R E S O L V E D SUBFILTER-SCALE, A N D SUBGRID-SCALE MOTIONS. RESOLVED, CO is T H E W A V E N U M B E R . IS T H E CQ C GRID CUTOFF 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 2003) 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 FOR T H E ^ - M O M E N T U M EQUATION 2 4 F I G U R E 3.4 TRAPEZOIDAL F I G U R E 3.5 S C H E M A OF T H E C H A N N E L F L O W FILTER IN 3 D F I G U R E 3.6 D I M E N S I O N L E S S S M A G O R I N S K Y F I G U R E 3.7 M E A N STREAMWISE VELOCITY yju'u lu U: WITH T H E A D M T M O D E L A N D T H E 3 7 R M S FLUCTUATION OF V : -JV'V WITH T H E A D M Iu r M O D E L A N D T H E M O D E L R M S FLUCTUATION OF y/w'w' I u W: WITH T H E A D M r M O D E L A N D WITH T H E A D M M O D E L A N D T H E M E A N STREAMWISE VELOCITY PROFILES F R O M T H E CASE 2 F I G U R E 3.12 DIMENSIONLESS R M S FLUCTUATION OF u: yju'u lu F I G U R E 3.13 DIMENSIONLESS R M S FLUCTUATION OF v.yjv'v' I u FIGURE 3.14 DIMENSIONLESS R M S FLUCTUATION O F F R O M T H E CASE 2 T yjw'w' I u W. F R O M T H E C A S E 2 T F I G U R E 3.15 S G S F I G U R E 3 . 1 6 RATIO OF S G S E D D Y VISCOSITY TO M O L E C U L A R VISCOSITY: M O D E L A N D T H E A D M M O D E L (IN C A S E 2) F I G U R E 3.17 M O D E L COEFFICIENTS F R O M T H E CASE 2 T DIMENSIONLESS FIGURE3.18 DIMENSIONLESS F I G U R E 3.19 DIMENSIONLESS R M S FLUCTUATION OF 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 FIGURE3.21 DIMENSIONLESS F I G U R E 3 . 2 2 DIMENSIONLESS yju'u' lu R M S FLUCTUATION OF U: V.JV'V lu 1 r AT Re T yju'u lu R M S F L U C T U A T I O N OFu: R M S FLUCTUATION OF A T Iu T v: yjv'v I u T r Re = A T T WITH T H E T Re r r 3 9 3 9 4 0 4 0 S M A G O R I N S K Y 4 2 = 1655 4 4 1655 4 5 = 1655 Re 8 4 1 v lv AT 3 4 1 WITH T H E M E D I U M M E S H CASE (CASE 2) M E A N STREAMWISE VELOCITY PROFILES 8 S M A G O R I N S K Y M O D E L DIMENSIONLESS 3 T H E M O D E L 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\ F I G U R E 3.11 2 T H E 3 7 R M S FLUCTUATION OF F I G U R E 3.9 D I M E N S I O N L E S S S M A G O R I N S K Y M O D E L A N D M O D E L F I G U R E 3.8 D I M E N S I O N L E S S S M A G O R I N S K Y PROFILES WITH T H E A D M M O D E L DIMENSIONLESS S M A G O R I N S K Y 3 1 3 4 5 = 1655 4 6 WITH 96X128X96 M E S H 4 6 WITH 96x128x96 M E S H 4 7 F I G U R E 4.1 F L O W LOOP IN T H E E X P E R I M E N T ( F R O M SHARIATI 2002) 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 SHARIATI 2002) 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 FOR T H E C O N V E R G I N G SECTION. O N L Y ASUBSET OF T H E A C T U A L GRID PLOTTED F I G U R E 4.4 M E A S U R E M E N T LOCATIONS A L O N G T H E CENTERLINE OF T H E C O N V E R G I N G SECTION (LINE A N D T H E VERTICAL LINE AT 0 . 0 3 M AFTER INLET (LINE D E ) F I G U R E 4.5 F I G U R E 4.6 IS 5 5 5 6 C O M P A R I S O N OF T H E M E A N VELOCITY AT T H E CENTERLINE OF T H E C O N V E R G I N G SECTION C O M P A R I S O N OF R M S V A L U E OF T H E VELOCITY vii A B C ) FLUCTUATION u N O R M A L I Z E D B Y M ^ 5 6 5 7 F I G U R E 4.7 C O M P A R I S O N OF R M S V A L U E OF T H E VELOCITY FLUCTUATION V N O R M A L I Z E D BY F I G U R E 4.8 C O M P A R I S O N OF R M S V A L U E OF T H E VELOCITY FLUCTUATION W F I G U R E 4.9 C O M P A R I S O N OF T U R B U L E N C E KINETIC E N E R G Y k NONDIMENSIONALIZED BY 57 u' jn N O R M A L I Z E D B Y 58 u' n k AT T H E INLET OF THE POSITION 0 . 0 0 6 M AFTER T H E C O N V E R G E N C E STARTS FIGURE 4.10 C O M P A R I S O N OF T H E FLUCTUATING VELOCITY 58 U AT INLET FOR O N E C H A N N E L INFLOW CONDITION WITH T H E M E A S U R E D V A L U E AT T H E VERTICAL LINE ( D E LINE IN F I G U R E 4.4) AT T H E LOCATION 0.03M F R O M T H E INLET 59 F I G U R E 4.11 T H E D I A G R A M OF " T W O - C H A N N E L " INFLOW CONDITION 59 FIGURE 4.12 C O M P A R I S O N OF T H E M E A N VELOCITY AT T H E CENTERLINE OF T H E C O N V E R G I N G SECTION WITH DIFFERENT INFLOW CONDITIONS FIGURE 4.13 6 0 C O M P A R I S O N OF R M S V A L U E OF T H E VELOCITY FLUCTUATION U N O R M A L I Z E D BY M ^ W I T H DIFFERENT INFLOW CONDITIONS FIGURE 4.14 60 C O M P A R I S O N OF R M S V A L U E OF T H E VELOCITY FLUCTUATION N O R M A L I Z E D BY V WITH u' n DIFFERENT INFLOW CONDITIONS FIGURE 4.15 61 C O M P A R I S O N OF R M S V A L U E OF T H E VELOCITY FLUCTUATION W N O R M A L I Z E D B Y w' wiTH n DIFFERENT INFLOW CONDITIONS FIGURE 4.16 61 S C H E M A T I C OF N E W INFLOW GENERATION T E C H N I Q U E FIGURE 4.17 C O M P A R I S O N OF R M S V A L U E OF T H E VELOCITY FLUCTUATION 64 AT INLET WITH T H E U 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 LINE (LINE D E IN 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 INLET FIGURE 4.18 C O M P A R I S O N OF STREAMWISE C O M P O N E N T OF T H E M E A N VELOCITY A T T H E INLET WITH 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 LINE (IN F I G U R E 4.4) FIGURE 4.19 65 T H E L O C A T E D 0 . 0 3 M F R O M T H E INLET M E A N STREAMWISE VELOCITIES C A L C U L A T E D WITH T H E N E W INFLOW CONDITION AT 65 T H E V E R T I C A L L I N E ( D E LINE IN F I G U R E 4.4) C O M P A R E D WITH M E A S U R E D V A L U E S A T T H E S A M E LOCATION. FIGURE 4.20 C O M P A R I S O N OF R M S V A L U E OF T H E VELOCITY FLUCTUATION U AT T H E VERTICAL LINE ( D E ..66 LINE IN 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 F I G U R E 4.21 66 M E A N VELOCITY AT T H E CENTERLINE OF T H E C O N V E R G I N G SECTION C A L C U L A T E D WITH N E W INFLOW CONDITION C O M P A R E D TO M E A S U R E D V A L U E S AT T H E S A M E LOCATION FIGURE 4.22 C O M P A R I S O N OF R M S V A L U E OF T H E VELOCITY FLUCTUATION 67 U N O R M A L I Z E D BY u' V N O R M A L I Z E D BY u' jn C A L C U L A T E D WITH T H E N E W INFLOW CONDITION FIGURE 4.23 67 C O M P A R I S O N OF R M S V A L U E OF T H E VELOCITY FLUCTUATION n C A L C U L A T E D WITH T H E N E W INFLOW CONDITION FIGURE 4.24 68 C O M P A R I S O N OF R M S V A L U E OF T H E VELOCITY FLUCTUATION 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 WITH T H E N E W INFLOW CONDITION FIGURE 4.25 68 C O M P A R I S O N OF T U R B U L E N C E KINETIC E N E R G Y K NONDIMENSIONALIZED B Y A: A T T H E I N L E T OF T H E POSITION 0.006M A F T E R T H E C O N V E R G E N C E STARTS WITH T H E N E W INFLOW CONDITION FIGURE 4.26 C A S E 69 C O M P A R I S O N OF T H E M E A N VELOCITY AT T H E CENTERLINE OF T H E C O N V E R G I N G SECTION WITH 1 A N D C A S E 2 FIGURE 4.27 71 C O M P A R I S O N OF R M S V A L U E OF T H E VELOCITY FLUCTUATION U N O R M A L I Z E D B Y u' WITH C A S E W N O R M A L I Z E D B Y u WITH C A S E n 1 A N D C A S E 2 FIGURE 4.28 72 C O M P A R I S O N OF R M S V A L U E OF T H E VELOCITY FLUCTUATION jn 1 A N D C A S E 2 FIGURE 4.29 C A S E I F I G U R E 4.30 72 C O M P A R I S O N OF T H E M E A N VELOCITY AT T H E CENTERLINE OF T H E C O N V E R G I N G SECTION WITH A N D C A S E 3 73 C O M P A R I S O N O F R M S V A L U E OF T H E VELOCITY FLUCTUATION U N O R M A L I Z E D B Y u' W N O R M A L I Z E D B Y WITH n A N D C A S E 3 F I G U R E 4.31 C A S E I 73 C O M P A R I S O N OF R M S V A L U E OF T H E VELOCITY FLUCTUATION C A S E I , A N D C A S E 3 u' jn WITH 74 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 WITH DIFFERENT SUBGRID M O D E L S F I G U R E 4.33 75 C O M P A R I S O N O F R M S V A L U E O F T H EVELOCITY 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 WITH DIFFERENT SUBGRID 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' WITH H DIFFERENT SUBGRID 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 WITH in DIFFERENT SUBGRID 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 : V Iv W I T H T H E t 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 F I G U E R 5.3 C O N C E N T R A T I O N 83 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 FIBER L E N G T H WITH 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 INITIAL FIBER ORIENTATION DISTRIBUTION, (A) IN X - Z PLANE, (B) IN 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 . 0 4 5 M , ( 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 . 1 2 2 M , ( 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 , 92 (A) IN X - Z PLANE, (B) IN X - Y P L A N E 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) IN X - Z PLANE, (B) IN 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 . 2 2 7 M , ( 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 . 2 6 2 M , ( 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 a^a^Acc,) F I G U R E 5.14 T H E FIRST M O M E N T ( ^ FIGURE 5.15 T H E FIRST M O M E N T F I G U R E 5.16 T H E S E C O N D M O M E N T ( ^ F I G U R E 5.17 T H E S E C O N D F I G U R E 5.18 T H E T H I R D M F I G U R E 5.19 T H E T H I R D M (^ IN X - Z P L A N E 95 DIFFERENT STATIONS IN X - Y P L A N E 96 A T DIFFERENT STATIONS )Aflf .) A 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 i M O M E N T (^jT a]p(a )t\a ) 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 O M E N T ( o f p{a )l!s,a ) 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 i O M E N T ( ^ ofp{a)Aa) 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 t i i i i i 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 97 97 98 I N X - Y P L A N E W I T H 8000, 16,000 A N D 24,000 FIBERS F I G U R E 6.1 96 99 T H EG E O M E T R Y O F T H EC O N V E R G I N G SECTION 102 F I G U R E 6.2 O N E - C H A N N E L INFLOW CONDITION F O R T H EC O N V E R G I N G SECTION 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 CONDITION S H O W I N G T H E DETAILS O F O N E OF T H E FIVE IDENTICAL C H A N N E L S F I G U R E 6.4 105 C O M P A R I S O N O F T H EM E A N VELOCITY A T T H E INLET 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 F I G U R E 6.6 104 C O M P A R I S O N O F R M S V A L U E O F T H EVELOCITY 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 EVELOCITY F I G U R E 6.8 C O M P A R I S O N POSITIONS A L O N G T H E H E A D B O X FLUCTUATION F L U C T U A T I O N 106 U A T T H E INLET 106 V A T T H E INLET 107 W ' A T T H E INLET 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 ix 107 110 111 FIGURE 6 . 1 0 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION U AT THE FIRST STATION A B . . . FIGURE 6 . 1 1 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION V FIGURE 6 . 1 2 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION W AT THE FIRST STATION A B . 111 AT THE FIRST STATION A B . . . . 1 1 1 ..112 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 FIGURE 6 . 1 5 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION V AT THE SECOND STATION C D . FIGURE 6 . 1 6 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION W AT THE SECOND STATION C D . AT THE SECOND STATION C D . 113 113 114 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 . . . . FIGURE 6 . 2 0 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION W AT THE THIRD STATION E F . ..116 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 FIGURE 6 . 2 7 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION V AT THE OUTLET WITH THE SMAGORINSKY MODEL FIGURE 6 . 2 8 1 1 9 1 1 9 COMPARISON OF RMS VALUE OF THE VELOCITY FLUCTUATION W' AT THE OUTLET WITH THE SMAGORINSKY MODEL FIGURE 6 . 2 9 COMPARISON OF THE RATIO v lu FIGURE 6 . 3 0 COMPARISON OF THE RATIO vJ lu FIGURE 6 . 3 1 COMPARISON OF THE RATIO v I u CONDITIONS FIGURE 6 . 3 2 115 COMPARISON OF THE RATIO w lu 1 2 0 AT DIFFERENT POSITIONS 1 2 0 AT DIFFERENT POSITIONS 1 2 1 AT THE OUTLET WITH THE TWO DIFFERENT INFLOW 1 2 1 AT THE OUTLET WITH THE TWO DIFFERENT INFLOW CONDITIONS 1 2 2 X Nomenclature C Smagorinsky constant G grid filter function G test filter k kinetic energy of turbulence n unit normal p filtered pressure s R Reynolds number based on bulk velocity b R Reynolds number based on friction velocity S, Stokes number T - 1 du 3w, S =—(—'- + -—-) filtered rate of strain tensor u 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 U bulk velocity V velocity vector x, y, z distance parameter y non-dimensional distance ( h + Greek Letters xi Kronecker delta £ dissipation rate of turbulence M dynamic viscosity P density o Reynolds-stress tensor v eddy viscosity A grid filter width A test filter width T T V Kolmogorov microscale Superscripts 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 pulpfibersin 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-£ 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- £ can not simulate the turbulence well. There are two ways to solve this problem. One way is to modify the k£ equation model (using for example, a non-linear k-£ 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 2.1 Literature Review 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: • Spread stock evenly across the width of the machine. • Level out cross-currents and consistency variations. • Level out machine direction velocity gradients. • Create controlled turbulence to eliminate fiber flocculation. • 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 fiberflocculationin 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- £ 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 onfiberorientation than that of the flow rate. Thefibershave 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: • Fiber orientation distribution is independent of the flow rate through the headbox. • Fibers are more strongly oriented in the plane of the contraction than in the plane of the paper. • The shape of the headbox does not affect the final fiber orientation distribution leaving the headbox. • 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 onedimensional 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 2.2.1 8 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 / £ ) , 4 where V is kinematic viscosity, and £ 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 , which 9 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„ 12,300 30,800 61,600 230,000 N (grid number) 6.7X10 4.0X10 7 1.5x10 2.1X10 9 s 6 Chapter 2 Literature Review 9 U H Here R e = —-—, U H m 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,=u + u' t 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. _ du t —- + u —•- = dt dXj 1 dp - i -+ r-2- vV w, p dx., 1 ^ ij T — p dXj L f^ = 0, (2.1) (2.2) OX, where T = p u Uj is the Reynolds stress tensor. tj t 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 where u v k (23) v is referred to as the effective turbulent viscosity, 8 is the Kronecker delta and T tJ k = u. u 12 . The shear-stress model, however, computes the Reynolds stress directly, without t 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 — £, k — 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 — £ 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, can be 2 2 computed from the assumption of «, =— k — V —— in (2.3). With this assumption used in a 2 T 3 converging section of a headbox (where ox, is very large near the exit), w, can become negative 2 OJC, which of course is unphysical. This leads to the production term not being calculated accurately. The production term in the simpler k — £ 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. Domaindecomposition 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 £=0 (3.1) ox. ( du. dt | diijUj _ dxj where the indices i,j-1,2,3 1 o?J /Jctc, | v d\ dxdxj 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) <f> respectively, as <f) = <j)+<t>' where Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 14 (j) = \G(x-x*, A)<p(x")d x*, 3 D A is filter width, G is a filter function with property: J G(x — x*, A)d x* = 1 3 D A finite volume discretization involves a box-filter with G(\x-£\,A) Ax,Ax Ax = 2 Q 3 1 ' 1 otherwise « / = 1,2,3. After applying the above box-filter to equations (3.1) and (3.2), we can get the filtered Navier-Stokes equations: OXj ^+M _I^ 3/ = dxj + v pax, JS_i?i oXjdXj 3Xj (3 .4) where the subgrid-scale (SGS) Reynolds stress tensor is written as T =u v iUj - ,Uj u (3.5) In order to solve the equations (3.3) and (3.4), the SGS Reynolds stress must be modelled with subgrid scale models. 3.2.2 Subgrid model The most commonly used SGS model is the Smagorinsky model (1963). The model assumes the SGS stress follows a gradient-diffusion process, similar to molecular motion. Consequently, T tj is given by 15 Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow * -\*Aj=-lCA \S\S = -2v S 2 ij lj T (3.6) ij 1 du. du where v is the eddy viscosity related only to the smaller scale motions, S =—(——+ -—^-) is the 2 dXj ox, T u filtered rate of strain tensor, | S \= ^ISySy is the magnitude of the strain tensor and C is the s i Smagorinsky constant (0.1 < C < 0.24, Rogallo and Moin (1984)). Here A = (Ax, * Ax * Ax ) is the 3 s 2 3 filter width, Ax,.(/ = 1,2,3) is /-direction grid spacing. In the region close to the wall, the eddy viscosity has to be reduced, which is usually achieved by using a Driest damping function: v =C,[(l-e- / / 2 5 r )A] |5| (3.7) 2 Lilly (1966) postulated that V = C Aq T (3.8) L where q is the SGS kinetic energy and C is a closure coefficient. An equation for q can be 1 2 L derived from a moment of the Navier-Stokes equation, which involves several terms that must be modelled. It is difficult to conclude that any significant improvement over the Smagorinsky model can be obtained with such a model. The most complicated SGS model was created by Deardorff (1973) for application to the atmospheric boundary layer. The model consists of 10 partial differential equations and bears a strong resemblance to a second-order closure model. The computational effort was increased using this model, but almost no improvement in the results was obtained. Several alternative models for the subgrid scale turbulence have been proposed: Yakhot et al. (1989) proposed a subgrid-scale model based on renormalization group (RNG) theory for the study of channel flow. Lesieur and Metais (1996) showed that Kraichnan's spectral eddy viscosity can be implemented in physical space, yielding the so-called structoe-function model. 16 Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 3.2.2.1 Dynamic subgrid model A more flexible model for the subgrid turbulence has been devised by Germano et al. (1991), the dynamic SGS model. Their formulation begins with the Smagorinsky eddy-viscosity approximation. However, rather than assuming a fixed value of C a priori, they permit it to be s computed as the LES proceeds. This is accomplished by using two filters (test filter G, and grid filter G). The characteristic length scale of G, which is A , is larger than that of the grid filter G, A . The value of a = A / A depends on the filter type. From Najjar and Taffi (1996), Lund (1997), it is known that a = for a box-filter with a trapezoidal rule for a numerical test filter, a = 2 for a box-filter with Simpson's rule for a numerical test filter. The grid filter and test filter are applied to the momentum equations (3.2) to produce the following equations: m. 1 dp 6*u. cT., Lj ' — £_j-i/ ' L dt dxj pdx dXjdXj fiXj (3.9) t where the subtest stress is T (3.10) ij= i j- i ju u U u Now the test filter G is applied to the filtered momentum equations (3.4) du t dt dUjUj _ cbCj 1 dp p dx, &% dfy dXjdXj Sxj de v (3.11) dx. (3.12) Using equation (3.9) and (3.11), the resolved turbulent stress l can be expressed with v (3.13) The subgrid-scale and subtest-scale stress are then modelled by the Smagorinsky model: •2CA |5|5,=-:•2Cfi„ 2 (3.14) Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow T, = -2CA 1f1 f, = - 2 C a , 2 17 (3.15) where^=A |S|^.,a,=A |f|!,. 2 2 Substitution of (3.14) and (3.15) into (3.13) gives ^-\^S =-2C{a -p ). 9 v (3.16) v Equation (3.16) represents a set of five independent equations that unfortunately cannot be solved explicitly for C. Lilly (1992) suggested computing the optimum in the least squares sense. The error is Q = V,-\l*S +2CMtf (3.17) 9 where M,=A |f | f , - A 2 2 | ? ^ (3.18) is minimized by 8QI dC = 0 , which yields: £M '' 2M M u C(x,y,z,t) = - a 9 V (3.19) V The numerator tyMy can have both negative and positive values. This leads to eddy viscosities of both signs. A negative eddy viscosity is considered a way of representing backscatter, which takes place as energy is transferred from the small scales to large scales in turbulent flows. But negative viscosity may cause a numerical instability. One way to avoid this is to employ averaging. It is a common procedure to average the denominator and the numerator in homogeneous directions (Germano et al. 1991, Najjar and Taffl 1996) in turbulent channel flows. This results in: (£ M ) C ( y ' t ) = -nA, MS ( 3 - 2 0 ) Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow where the brackets ( ( ) ) represent an average over the spatial region, () 18 is a plane averaging operation in the x and z directions, which are usually homogeneous in a turbulent channel flow. For fully inhomogeneous flows, Ghosal et al. (1995) proposed a dynamic localization model based on a constrained variational formulation. Both the dynamic subgrid model and the simpler Smagorinsky model for subgrid turbulence will be used to produce the numerical results described in the following sections. 3.2.2.2 Approximate deconvolution model (ADM) For the finite volume method, a box-filter is considered 'implicit' filtering of the NavierStokes equations. This approach has the advantage that it does not need to define a filter function. There are no explicit filter operations. For the dynamic Smagorinsky model, the only test filter is constructed to model the subgrid scales. In contrast, a filter is defined when explicit filtering is used. The flow is divided into resolved and subfilter-scale (SFS) motions. The SFS motions can be further divided into resolved SFS (RSFS) and unresolved SFS (USFS) (Zhou et al. 2001, Carati et al. 2001). The USFS is commonly called SGS. The RSFS motions can be reconstructed from the resolved motions. The SGS (USFS) can not be represented directly and must be modelled. The schematic of energy spectrum is shown in Figure 3.1. Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 19 E Energy (O Figure 3.1 Schematic of energy spectrum illustrating the relationship among resolved, resolved subfilter-scale, and subgrid-scale motions. CO is the wavenumber. COc is the grid cutoff wavenumber (Gullbrand and Chow 2003). Stolz et al. (2001) used u to approximate the unfiltered solution u . The approximate deconvolution u is computed by applying the deconvolution operator Q to the filter quantity u N u=Q *u (3.21) N The operator Q is obtained by expanding the inverse of the filter G as an infinite series and N truncating it at finite N. This leads to Q as an approximation of G , _1 N ^ = £ ( ! - ( ? ) ' ( 3 . 2 2 ) where I is the identity operator (Stolz et al. 1999). It is found that the deconvolution order N = 5 is sufficient for numerical test cases (Stolz et al. 2001). Using (3.22), u can be computed by repeated filtering of u from u =Q u =u+(u-u) N + (u-2u +u) + ... (3.23) The subgrid-scale Reynolds stress tensor T, in (3.5) can be divided into two terms: SGS and RSFS Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 20 T.. = (upj - u Uj) + («,. Uj - u Uj ) t (3.24) t SGS(USFS) RSFS To model the transfer of energy to small scales, represented by U(Uj — w. Uj , a relaxation term ZU(I~QN *G)*U is added for the A D M model. The RSFS term in (3.24) is moved to the left of the filtered Navier-Stokes equation (3.4). The equation (3.4) with the A D M model can be written into at dXj pax, oXjOXj where % >0 is a relaxation coefficient. u The relaxation coefficient % can be determined dynamically (Stolz et al. 2001). The secondu order structure function (Lesieur and Metais 1996) is applied to (I — Q *G)*u N . The structure function F is defined as 2 F (^) = \ \ ^ + r,t)-^j)l^ 2 where h is the grid size in computational (3.26) space £ = (^,£ ,£ ) 2 3 , <p u with 4=(/-fi**G)*fi;. To estimate the relaxation coefficient % , the equation (3.25) is solved two times in one time u step, once using % = X u or a positive F (£, t" + At) | 2 a n t J < scales within At. n c —F 2 e using % = 0. % is the % value from the previous time step u constant = 0 o Ua at time t=0. u The difference of the structure function t") represents an estimate of the energy generation of non-resolved The difference F (£, t" + At) | 2 0 -F (£, t" + At) \ 2 energy would be dissipated by the relaxation term using X ~ X u Uo estimates how much • T h e n ^ can be computed from Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 21 t ' - ^ F t f s + w ^ s + m ^ ' This A D M model will be used in the following chapters to compare with the Smagorinsky model. Another alternative is using the dynamic Smagorinsky model to compute the subgrid scales instead of using the relaxation term in the A D M model. This leads to the dynamic reconstruction model (DRM) which combines the A D M model with the dynamic Smagorinsky model (Gullbrand and Chow 2003). 3.3 Numerical methods 3.3.1 Theflowsolver The equations governing incompressible and isothermal flow with constant density are written into a conservative form: J V * ndS-0 (3.28) |- \udV+ \[uV at * £ n+ P-n -v(n ^- + n ^- + n ^-)}dS = 0 p dx dy dz (3.29) ?-ivdV+[[vV dt • • n+ P-n -v(n ^ + n ^ + n ^)]dS = 0 p dx dy dz (3.30) x x y y x y z z |- \wdV+ f [ w V n+^-n -v(n ^ + n ^ + n ^)]dS = 0 dt J * p dx dy dz z x (3.31) 2 where V is the velocity vector, dV is a volume element, dS is an area element and n = n i + n j x y + n k is a unit normal of the control volume face. z The scheme used in this thesis is based on a fraction-step method (Chorin 1968, Kim and Moin 1987). The pressure gradient and the incompressibility constraint are integrated implicitly in time. The convective and diffusive terms are treated explicitly with the second-order Adams- 22 Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow Bashforth scheme in time. For finite volume discretization, the flux is computed by a second-order accurate averaging which is therefore equivalent to second-order central difference scheme. A general nonorthogonal coordinate system ( , £ £ ) is defined by 2 where r = (x, y,z) r 3 is the Cartesian coordinate system (He and Salcudean 1994, Rosenfeld et al. 1991). The computational domain (£ ,£ x £3) i divided into uniform general cells (or scalar control s 2 cells) with mesh size A£ . = 1, and the center of each general cell corresponds to the indices i,j,k. ( The general cell is illustrated in Figure 3.2. The surface area vectors, S (i = 1,2,3), are 1 defined at the centers of the control volume surfaces, and are locally parallel to the coordinate line i 2 ,j 2 ,* +— 2 ri 1 + — • 2J 1 . 1 , 1 1 2 2 , ' 1 \ 2 \ 1 ,k + — 2 2 • i,j,k i - 1 . 1 —. 1 2 1 . 1 1 + —, j 2 , 2 1 ,k s 1 • 2 1 . 2 J 1 , 1 / + —, ; + —, k Figure 3.2 Illustration of the primary cell 2 2 , 1 2 2 + —, k Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 23 Integrating the equation (3.28) over a primary cell, as shown in Figure 3.2, we get ( S ' . V W (S'.V),,/^ (S .V) 2 - (S .V)j. 2 j+1/2 l/2 + (S .V) i/2-(S .V) . = £ S • V = 0 3 3 (3.32) 1 k+ k 1/2 faces where ^ implies summation (with the proper signs) over all the faces of a primary cell. Each term of (3.32) represents the volume flux over the corresponding face. The volume fluxes are selected as the velocity unknowns. Let u = s'-v u = s -v x 2 2 = s -v 3 u 3 where U , U and U are the volume fluxes defined at the center of £ , £ and £ faces in a l 2 3 t 2 3 primary cell. Then the equation of (3.32) can be written as f/ _f/' +TJ -T/ + f/ _f/ ^ 1+1/2 '-'i-l/2 y+1/2 j-l/2 ^ 4+1/2 ^ k-l/2 1 2 T or L J 3 2 u 3 u =0 " D (U')=0, (3.33) iv where U = (U ,U ,U ). l ] 2 3 The summation operator D is a discrete divergence-like operator (the jv divergence operator itself is ^ D , where V is the cell volume). iv Different computational cells are used for the discretization of the momentum equations. Each cell has the element size of A£, X A£ x A£ in the computational space, but centers are located 2 3 differently at (i + ^,j,k ), (i,j + ^,k ), and (i,j,k + ^ ) for the U\U and U momentum 2 3 equations. This grid uses the staggered-grid convention. The control volume used for the discretization of £/' equation is shown in Figure 3.3. Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow Figure 3.3 Control volume for the 24 -momentum equation. Dotting S V i / 2 with (3.29), (3.30) and (3.31) produces (+1/2 — S i+l/2 i+i/2 F, *i / dt + 1 (3.34) 2 where F / is the total flux through the computational cell of the ^ -momentum equation, and V i+1 2 MI2 is the volume of the control cell. The equation (3.34) also can be written as dUJ +U2 1+1/2 dt _. " A (3.35) where L, is the right side of equation (3.34). £, can be split into two parts, £ , = # , ( £ / ' ) + *, CP), (3.36) where H (U') is an operator representing the discretized convective, diffusive terms , and i?, (P) is t the operator of discretized pressure terms. In the first step, the pressure terms are dropped from the equation (3.35). An AdamsBashforth scheme is used to advance the time step. Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow C/2-^"/2 ^ -(3^ (C/ )-// (C/' - )) + £ / n n 1 1 1 25 (3.37) where U)* is an intermediate velocity which does not generally satisfy the continuity equation. xl2 Substitution of (3.37) into (3.35) yields U^=Ul' -^R (p^) ]/2+ (3.38) t "i+l/2 In a similar way, the £ - and £ -momentum equation, for which the control cell centers are 2 3 located at (/, j + , k), and (/, j, k + -j ) respectively, can be discretized as Kvi = ( ^ 2 (U' ) - H (U' )) + 3 n 2 U j ^ U ^ - ^ R , ^ ) Ui; =Ul: m Ul^ = Ui; 2 where H (U ) l 2 terms , and (3.40) -£t-(3H,(U'")-H (U>"- )) 1 u2 + m + (3.39) 3 -^R (p^) 3 (3.41) (3.42) and H (U ) are the operators representing the discretized convective, diffusive l 3 (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 -D (U' ) = AtD X- ± - -) t iv J E J i "m+1/2 or (3.43) Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 26 A v ( 4 ^ ) =- ^ A v ( 0 where m-i,j, (3-44) 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 {u ,u ,u ,...} 0 x 2 that end when \\r \\ =\\Q- Au \\ is small enough. The effectiveness of this method for sparse k 2 k 2 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 nonorthogonal 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 (BICGSTAB) 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 3.3.3.1 28 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 — 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 subdomain 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): „. . wall clock time for a run on one processor S(p) = £ , 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. 30 Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 1 4 8 Cases 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) jc—Ax The box-filter then becomes G(x-x) =— 2Ax G(x-x') = 0 if \x-x\<Ax (3.46) 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)). 4 (3.48) 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(a Tj )/a 0 + l) Q = 0,\,...,Ny-l). J (3.49) where Ny is the total number of grid points in y direction. The quantity a is a stretching parameter ( 0 < a <1), a =^log((l + a)/(l-a)), 0 Tjj = -l + 2jAn and An = \/(Ny-\). Here a = 0.98846 is used to get finer meshes near the wall. 3.4.1 Reynolds number Re = 395 r The parameters of the first test case have the following values: Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 33 H = 0.075 m, U = 0.204m/s and v - l x l O " m /s, 6 2 b Where U is the bulk velocity. The Reynolds number based on friction velocity is b Re, = ^ - 3 9 5 v where u = ^Jr I p is the friction velocity , T is the wall shear stress and 8 is half of the channel T w w height H. The Reynolds number based on bulk velocity is Re,=-^«7650. v The flow is turbulent since Re^ is greater than Re^ crit ~ 575 (White 1998). The simulations are performed with three meshes. Table 3.2 gives the grid parameters. Case Nz Nx Ny AJC Ay . + + Av + /max 1 32 32 32 65.0 32.6 1.5 56.0 2 64 64 64 37.8 18.9 0.8 32.8 3 96 96 96 25.6 12.8 0.5 22.2 Table 3.2 Grid parameters for R e = 395 r The initial mean velocity profile was generated by interpolating data obtained from a computation with the common k - £ model. Random perturbations were added to the mean velocity distribution. The simulations 1, 2 and 3 correspond to the coarse, medium and fine mesh cases. The CFL number (CFL = v ———;— j 1 dtconv dtdif , where dtconv = min(| —1,1 —1,1 — |) C/ K dtdif = min(—— ,—— ,—— ) ) is 0.15 for all three cases. The time steps are 2(v+v )2(v + v )2(v + v ) T T T Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 0.002T, 0.0018T and 0.0012T, respectively, where 1-HIU . b 34 The fluctuating velocities, u', are computed from: =",-<", >, where u is the instantaneous velocity and < u > is the tempo-spatial average of the velocity field. i j 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 C s is set to 0.1. The dynamic model computes the constant from (3.20). The relaxation coefficient % is u 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: v tol =max(0,v+v ). T 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]u u Iu with two different models. Figure 3.10 shows the Reynolds stress u'v I u . 2 j i T 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 A D M (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 — 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 apartfromwall regions. In Figure 3.15, the dynamic coefficient yJC(y) and the van Driest damping on the Smagorinsky coefficient C ( l — e~ ) are compared from the medium mesh case. The peak value y 125 s 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 samefromFigures 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. 37 Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow Figure 3.7 Dimensionless rms fluctuation of u: yjuu lu with the ADM model and the T Smagorinsky model. 38 Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow Figure 3.8 Dimensionless rms fluctuation of v: V v V lu with the ADM model and the T Smagorinsky model. Figure 3.9 Dimensionless rms fluctuation of w: -Jw'w' I u with the ADM model and T the Smagorinsky model. Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 64*64*64 ADM _-i I • i 0 •i Ii ••• I 50 100 i •i • I•> • .I•• . • I•.•• I.. •• I.•. • I 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 DNS MKM Figure 3.13 Dimensionless rms fluctuation of v: V v V lu from the case 2. T 41 Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow _ _ _ _ 1.4 Smagorinsky Dynamic model ADM model N=5 ADM model N=10 DNS MKM • ••••••••• 1.2 tl I ) I—I 0 ——— — 1 1 1 1 100 I I I I I 200 y+ I I I I I I 300 I I I 400 Figure 3.14 Dimensionless rms fluctuation of w: yjww I u from the case 2. T Figure 3.15 SGS model coefficients with the medium mesh case (case 2). 42 Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow 0.6 r • O Smagorinsky model ADM model 04 o • • • 0.2 o • • o p 8 • o • • O n o o ° n D oo 3© 0.4 0.2 y/H 0.6 0.8 Figure 3.16 Ratio of SGS eddy viscosity to molecular viscosity: V lv with the Smagorinsky T model and the ADM model (in case 2). 3.4.2 Results for Reynolds number Re = 1655 r The parameters of the second test case have the following values: H = 0.075 m, U =0.204m/s and v = 0.2441xl0" m /s, 6 2 h Where U is the bulk velocity. The Reynolds number based on friction velocity is b R e r = ^l655 v where u = ^T I p is the friction velocity , T is the wall shear stress and S is half of the channel T w w height H. The geometry is the same as that with Re = 395 . The initial distribution comes from r the previous Re = 395 case. The grid parameters are shown in Table 3.3. r 43 Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow Case Nz Nx 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 = 1655 T 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 wallregion, is used as the reference profile, with K = 0.40, B = 5.5 : y + = + u +e , K ) 2 2 K ) 6 3 ~ 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 = 1655. The results from the Smagorinsky model are a little better than those from the r 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 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 u varies between the different cases and its vaue is very small, using the shear velocity as a T non dimensional parameter somewhat exaggerates the differences between different cases. Figure 3.17 Dimensionless mean streamwise velocity profiles at Re = 1655. r Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow Figure 3.19 Dimensionless rms fluctuation of w.yfvV lu at Re = 1655. T r 46 Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow G -1 I 0 I I 1 1 1 1—I 500 I I y+ I L_l 1000 I i I i i 1500 Figure 3.20 Dimensionless Reynolds stress: u'v lu at Re = 1655. 2 T 0 1 — i — i — i — i — I — i — i — i — i 0 500 y+ I r i i 1000 i—J I i— 1500 Figure 3.21 Dimensionless rms fluctuation of u: yju'u lu with 96x128x96 mesh. T 47 Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow Figure 3.22 Dimensionless rms fluctuation of v: yjv'v lu with 96x128x96 mesh. x 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 or t n e p i i e l computation of a large eddy simulation. ara Chapter 3 Large Eddy Simulation for Fully Developed Channel Flow The channel flow at R e T 48 = 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 = 1655, the channel flow was computed by using two r 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 nonuniform 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 — £ model and a Reynolds stress model (RSM) to his experimental measurements. He found that the turbulence kinetic energy predicted with the k — £ 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 L E S 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 L E S 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 A D M model. A l 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 53 Chapter 4 L E S for Turbulent Flow through a Converging Section 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. U is set to the mean streamwise velocity at the exit plane. This condition (4.1) c 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 nondimensionalized with the bulk velocity at the inlet of the converging section, U* = UIU . Here b X* = xlLx is the x position normalized with the total length of the converging section L . The x 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 (u ) are presented in Figures 4.6 - 4.8. The jn 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 L E S 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 — £ model in most of regions is close to the measured values, but the computed values with the standard k — £ model increase rapidly towards the exit. The reason is that the standard k — £ 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 — £ 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 "twochannel" 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 (u ) with two inflow conditions. The results of the "two-channel" inflow conditions are better in 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. 56 Chapter 4 LES for Turbulent Flow through a Converging Section 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 o o 9 8 7 F> • LES r~\ U* measurement • 6" 5 4r> O O J I I I I 0.25 I I ! I I 0.5 X* (x/Lx) I I I n n I I 0.75 Figure 4.5 Comparison of the mean velocity at the centerline of the converging section. Chapter 4 L E S for Turbulent Flow through a Converging Section 1.2 1.1 O 0.9 • o u'/u' in measurement o • r • 0.8 • o 3 LES • 1 • o 0.6 = 0.5 •• o •• • • 0.4 o D D oO • 0.3 n 0 0.2 • 0.1 01—I_I——I i 0 0.25 i i i I 0.5 1 1 > 1 X* (x/Lx) 1 1 1 1 0.75 1 p • 1 1 Figure 4.6 Comparison of rms value of the velocity fluctuation u normalized by 1.2 1.1 LES 1 v7u' in measurement O 0.9 0.8 •o 0.7 3 5 o o o 0.6 o o o°o ••••••••• 0.5 •• 0.4 0.3 0.2 £ ] • • • • • • • •• • • •• 0.1 r . ' 0.25 1 1 1 1 • 1 0.5 ' X* (x/Lx) 1 J 0.75 I L_L. Figure 4.7 Comparison of rms value of the velocity fluctuation v normalized by u' n Chapter 4 L E S for Turbulent Flow through a Converging Section 1.2 1.1 o LES • 1 O W/u' in measurement coo • 0.9 0.8 c 0 7 3 0.6 ; o o Q *0.5 • 0.4 •• 0 o° a,oo • 0.2 0.1 J I 0.25 I I I I 0.5 I I I L X* (x/Lx) _L J 0.75 i i i_ Figure 4.8 Comparison of rms value of the velocity fluctuation w normalized by u. 7|A LES • O A k measurement A A 2h A A -i i i O O O O 1 • 0.25 1 i i_ 0.5 X* (x/Lx) 0.75 j I i u 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 L E S for Turbulent Flow through a Converging Section u' for one channel inflow condition O • 1 u ' measurement at 3cm from inlet plane • • • • • • 0 oQ 0.9 0.7 • 0.6 Z* O O O f 0.8 • • • • 0.5 0.4 h 0.3 f a. • o o o o o o o • • • 0.2 0.1 0 t M B - JL 1 0.01 0.02 •• I •,£.• 0.03 u' o 1 1 1 1 1 1 1 0.04 0.05 0.06 o o o o o c II I 0.07 0.08 (m/s) 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. interolation H/2 H/2 t 0.001m interolation Figure 4.11 The diagram of" two-channel " inflow condition. 60 Chapter 4 L E S for Turbulent Flow through a Converging Section 11 D 11111 10 LES with one channel inflow condition • o i> 9 8 U* measurement LES with two channels inflow condition O o s 7 6 5 3 o > o 1 1 1 1 2 1\ m 1 1 n q• o - 4 1 1 J 1 I 0.25 I L O 0.5 O Q p P 0.75 X* (x/Lx) Figure 4.12 Comparison of the mean velocity at the centerline of the converging section with different inflow conditions. AA " J= ^•••••Q 0.9" 0.8 . 0.7 ° • LES with one channel Inflow condition uVu'In measurement A - LES with two channels Inflow condition O L A • — nA j=j 6 nO ^ a A O i 0.6 o ' 0.5 O o AO 0.4 A 0.3 0.2 Q Q r A A 0.1 0.25 0.5 X* (x/Lx) 0.75 Figure 4.13 Comparison of rms value of the velocity fluctuation u normalized by u' with different n inflow conditions. Chapter 4 L E S for Turbulent Flow through a Converging Section 61 1.2 1.1 c •o 0.9 0.8 LES with one channel Inflow condition O v 1 v7u' in measurement LES with two channels Inflow condition O nV?) V V o 0.7 Q c O V "a 0.6 5 0.5 0.4 V, V v ••••••n n D D VV V V o 0.3 0.2 0.1 0 0.2 0.4 0.6 0.8 X* (x/Lx) Figure 4.14 Comparison of rms value of the velocity fluctuation V normalized by u' with n different inflow conditions. 1.2 1.1 1 0.9 O LES with one channel inflow condition O w'/u' in measurement <] LES with two channels inflow condition Q Q 0.8 c 0.7 IPO ° "= 0.6 *0.5 0.4 03 <L O =0 ^ A O 0.2 0.1 0 0.2 0.4 0.6 X* (x/Lx) 0.8 Figure 4.15 Comparison of rms value of the velocity fluctuation w normalized by u' with different in inflow conditions. Chapter 4 L E S for Turbulent Flow through a Converging Section 4.2.3 62 Alternative methods for the generation of inflow data It appears that the turbulence levels near the inlet of the channel are lower than measured values so that some new method for prescribing inlet conditions needs to be considered. The most straightforward way to get the appropriate inflow condition for the converging section is to include the entire length of rectifier tubes in the calculation. This procedure would be a very costly venture in terms of computational effort. Furthermore, the inflow condition for the rectifier tubes is also needed to start this calculation. The simplest way to generate turbulent inflow conditions is to superpose random fluctuations on a desired mean velocity profile. The amplitude of the fluctuations can be estimated from a desired set of one-point second-order turbulence statistics. This approach has been employed to generate the inflow condition in, for example, Lee et al. (1992) and Rai and Moin (1993). This method requires a fairly lengthy straight section to allow for the development of the appropriate organized turbulent motion, which is not available in this converging section application. Bonnet et al. (1997) propose an approach to recover the two-point correlations of the inflow by using linear stochastic estimation. But this requires knowledge of the time history of the velocity at a few reference points and of the two-point spatial correlations tensor. These requirements can not be met from the experimental data. Lund et al. (1998) developed an extracting/rescaling technique to generate turbulent inflow data for spatially developing boundary layer simulations. Their method consists of the following steps: 1. Extract mean values and fluctuations from a down stream plane. 2. Rescale the preceding values according to boundary-layer similarity. 3. Reintroduce the rescaled values at the inlet plane. This method is shown to be highly accurate for spatially developing boundary layer simulations from Lund et al. (1998). For the converging section, there is no boundary layer similarity theory. Nevertheless, a modified Lund's method is constructed to generate the inflow data in this case. The geometry is the same as the channel flow geometry. The inflow obtained using the random Chapter 4 L E S for Turbulent Flow through a Converging Section 63 fluctuation method (see appendix C) is set to have the similar mean profile and Reynolds stress tensor distributions as those produced by Shariati (2002) using the RSM model for the inlet condition of the converging section. As before, the outflow condition is the convective boundary condition (4.1). The simulation is run to a statistically stationary state. The fluctuation values are extracted from a rescaling station (middle plane of the channel) shown in Figure 4.16. The rescaling values are not obtained from boundary4ayer similarity as was done by Lund et al. (1998). In the present case, they all have a common factor, the ratio of the computational rms velocity fluctuation u and the experimental value at the centerline of the converging section. The rescaled values are put back as the inflow data. After many time steps, the instantaneous velocity distribution on the exit plane is written to disk as the inflow data for the converging section. The new inflow rms velocity fluctuation u and the measured value at 0.03m from the inlet plane (line DE in Figure 4.4) are shown in in Figure 4.17. The result is much better than that in Figure 4.10. Figure 4.18 shows the new inflow mean streamwise velocity and the measured value of the same velocity along a vertical line near the channel entrance (line DE in Figure 4.4). Note that the measured values are for a location slightly downstream of those used for the inlet conditions but that there would be essentially no change in the streamwise direction if the channel were long and had fully developed flow. There is a reasonable agreement between the new inflow values and their measured counterparts. Using the new inflow data, the calculation of the converging section in section 4.2.2 is repeated. Figure 4.19 shows the new mean streamwise component of the velocity along the vertical line (line DE in Figure 4.4). The computational velocity profile is similar to the measured profile. The values of the streamwise fluctuating component of velocity at the same location as Figure 4.19 are presented in Figure 4.20. The computational results are in reasonable agreement with the measured data. It is observed that there are still differences between the computational values and the experimental results in the middle region. One reason for this difference could be that the mesh is 64 Chapter 4 L E S for Turbulent Flow through a Converging Section coarse in the middle region because the streching mesh is used in the wall direction. The mean streamwise velocity along the centerline of the converging section is shown in Figure 4.21. 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 (u ) are presented in Figures 4.22 —4.24. in The Figures 4.22 and 4.24 are almost the same as those computed from the earlier one channel inlet condition and plotted in Figures 4.6 and 4.8. From Figure 4.23, the resolved rms fluctuation v with the new inflow condition is better than that using the one channel inflow condition of Figure 4.7. The difference between w'from LES and the measured values remains quite large near the inlet region. The inlet conditions need further improvement. The turbulence kinetic energy k non- dimensionalized by k at the inlet of the position 0.006m after the convergence starts with the new inflow condition is shown in Figure 4.25. There is little difference between the new inlet condition and one-channel inlet condition. random inflow data rescaling station scaling of the rms values according to the experimental data Figure 4.16 Schematic of new inflow generation technique. 65 Chapter 4 L E S for Turbulent Flow through a Converging Section • 1 F 0.9 0.8 • u'at Inlet from new inflow condition O u' measurement at 3cm from inlet plane 0.7 0.6 Z* • • • o o o o o o o • 0.5 • 0.4 a • 0.3 • 0.2 0.1 ' 0 ', 1 1 1 0.01 ' •' ' ' 0.02 0.03 o unJS! 0.04 u' (m/s) • • • 0.05 o o o o o o o * 0.06 0.07 0.08 Figure 4.17 Comparison of rms value of the velocity fluctuation u at inlet with the measured value at the vertical line (line DE in Figure 4.4) of the location 0.03m from the inlet. 1 r • 0.9 0.8 r CD - U measurement at 3cm from inlet plane Q o 0.7 '- o 0.6 j - • Z* 0.5 7 ,o 0.4 r 0.3 7 0.2 7 • • • U at inlet from new inflow condition Q] oon o 0 o o O • • • • a 0.1 "- -e—i—L a 1 0.1 J o U (m/s) 0.2 0.3 Figure 4.18 Comparison of streamwise component of the mean velocity at the inlet with the measured value at the vertical DE line (in Figure 4.4) located 0.03m from the inlet. Chapter 4 L E S for Turbulent Flow through a Converging Section 1 C D 0.9 Q 0.8 ° ° 0.7 O O • U L E S at 3 c m H n * D E O U m e a s u r e m e n t at 3 c m line D E 0.6 • 8 0.5 D 0.4 tr 0.3 T 0.2 r 0.1 •0 oC P p r • _QL 0 D D D o • o • o • ° O n • n ° RD 1 0.1 o n 1 1 1 0.3 0 . 2 U (m/s) Figure 4.19 Mean streamwise velocities calculated with the new inflow condition at the vertical line (DE line in Figure 4.4) compared with measured values at the same location. 1 • § 0.9 0.8 0.7 • u" LES at 3cm line OE u* measurement at 3cm One DE • • • • Z* 0.5 Z • O O O 0.6 oO 0.4 • O c • <& • • o c o o • o CO 0.3 •o 0.2 0.1 J QJ. 0.025 -GL P 0.05 P • •• u' (m/s) 0.075 ' ' 1 l_ 0.1 Figure 4.20 Comparison of rms value of the velocity fluctuation u at the vertical line (DE line in Figure 4.4) with new inflow condition. Chapter 4 LES for Turbulent Flow through a Converging Section 11 67 r- • 10 - • 9 - o o LES Measurement O 8 - 9 7 6 • 5 OQ o • 4 3 O 2 1 o o ILLLLTI O • o • I 0 0.4 0.6 X* (x/Lx) 0.2 -J • • I I • 1 0.8 Figure 4.21 Mean velocity at the centerline of the converging section calculated with new inflow condition compared to measured values at the same location. 1.4 i- 1.2 h o • ° • 'F • o 0.8 o o 0.6 o 0.4 0.2 r- • LES O Measurement 0.2 _l I I 0.4I I L.0.6 X* (x/Lx) o o 0.8 Figure 4.22 Comparison of rms value of the velocity fluctuation u normalized by u calculated in with the new inflow condition. Chapter 4 LES for Turbulent Flow through a Converging Section 68 1.2 1.1 1 • LES O Measurement 0.9^ 0.8 c 3 0 o ! fa 7 o 0.6 E- • ••• ^o.sE0.4 r 0.3 r 0.2 0.1 r o, _i 0.4 0.2 i 0.6 i i X* (x/Lx) i * 0.8 ' 1 i_ Figure 4.23 Comparison of rms value of the velocity fluctuation v normalized by u' calculated in with the new inflow condition. 1.2 1.1 • 1 o o LES Measurement 0.9 0.8 = 0.6 too o o '*0.5 • 0.4 0. o, So o • 3[ 0.2 0.1 0 0.2 0.4 0.6 X* (x/Lx) 0.8 Figure 4.24 Comparison of rms value of the velocity fluctuation w normalized by u' in calculated with the new inflow condition. 69 Chapter 4 L E S for Turbulent Flow through a Converging Section • O A LES Measurement ^ k-e 4|- A 1 3 A 2 : or i i i 0 A i i i i i i i 0.2 0.4 i i i 0.6 X* (x/Lx) i i i i i 0.8 1 Figure 4.25 Comparison of turbulence kinetic energy k nondimensionalized by k at the inlet of the position 0.006m after the convergence starts with the new inflow condition. 4.3 Discussion One possible reason for the apparent difference between measured and predicted results is that an effectively two-dimensional periodic condition is imposed in the spanwise direction while in the experimental apparatus there are walls on both sides of the channel. Another explanation involves the accuracy of the experimental data. In the experiment, there is the possibility of different flow rates in the adjacent rectifier tubes. From Shariati (2002), it is known that the velocity profile is different for the top and the bottom rectifier tube in the experimental measurements. This will lead to inflow conditions that are unique to the apparatus and impossible to model with precision. In this section, the sensitivity of the simulation results with respect to grid size and the subgrid model is tested. This is necessary to assess whether the computational results represent a realistic converged solution. Chapter 4 L E S for Turbulent Flow through a Converging Section 4.3.1 70 Dependence on grid resolution For the present configuration questions of adequate resolution and proper representation of the inflow conditions are closely related. From chapter 3 for fully developed channel flow, LES requires rather fine mesh spacing to resolve the large scales of motion. The Smagorinsky model with 64x64x64 mesh can provide a good inflow condition for the converging section, but it still needs a finer mesh to resolve the core region of the channel flow. Here we use one case: 64x64x96 to test the effect of grid refinement in the normal direction to the wall. It is known that 64 grids in the spanwise direction are enough to resolve the large scales for fully developed channel flow in chapter 3. The spanwise length of the converging section is the same as that of the channel flow. We can therefore expect that 64 grids in the spanwise direction will be sufficient for the converging section. The grid parameters for different cases are shown in Table 4.2. The comparison of mean streamwise velocity along the centerline of the converging section between Case 1 and Case 2 is shown in Figure 4.26. Comparisons of resolved rms fluctuations u and w normalized by u at the inlet (u ) between Case 1 and Case 2 are presented in Figures 4.27 -4.28. The meshes of Case 2 in !n the x direction are double those of Case 1. The results of Case 2 are almost the same as those of Case 1. A comparison of the mean streamwise velocity along the centerline of the converging section for Case 1 and Case 3 is shown in Figure 4.29. The comparisons of resolved rms fluctuations w'and w normalized by w'at the inlet («,„) for Case 1 and Case 3 are presented in Figures 4.304.31. The results from a fine mesh Case 3, 64x64x96 are slightly closer to the experimental data than those of the coarser mesh case but there is not a great change or improvement from a refined mesh. If a finer mesh, such as 64x64x128,64x64x256, would be used, the results would be expected to improve. However the computational cost would become prohibitive. 71 Chapter 4 L E S for Turbulent Flow through a Converging Section Case Nx ~~T~ ~64 Ny Nz Ax/S Ay/S Az • 18 Az IS 64 64 0.03125 0.025 0.002531 0.01277 2 128 64 64 0.015625 0.025 0.002531 0.01277 3 64 64 96 0.03125 0.025 0.00166 0.008517 Table 4.2 Grid parameters. Nx, Ny, Nz denote the number of meshes in the streamwise, spanwise and wall-normal directions, respectively. 5 is half of the inlet height of the converging section. Figure 4.26 Comparison of the mean velocity at the centerline of the converging section with Case 1 and Case 2. Chapter 4 L E S for Turbulent Flow through a Converging Section 1.41- 1.2 A A • A • E .0.8 c A 3 0.6 A 0.4 • LES 64*64*64 A LES 128*64*64 • A A 0.2 _l I I I 0.2 I I L I 0.4 I I I 0.6 X* (x/Lx) J I I I 0.8 I I l_ Figure 4.27 Comparison of rms value of the velocity fluctuation u normalized by u' with in Case 1 and Case 2. Figure 4.28 Comparison of rms value of the velocity fluctuation w normalized by u' with n Case 1 and Case 2 73 Chapter 4 LES for Turbulent Flow through a Converging Section 10 • A 64*64*96 LES 64*64*64 LES • 300 i i i I 0.2 i i i LJ I 0.4 • •• I 0.6 I I X* (x/Lx) I 0.8 I I 1_ Figure 4.29 Comparison of the mean velocity at the centerline of the converging section with Casel and Case 3. 1.4 1.2 On • • A^A • A A • • A • • 0.6 A • A A • • o o o 0.4 h A • o 3 0.2 A o .0.8 a A A • A 64*64*96 LES • • 64*64*64 LES A Measurement o J 0.2 I I I I 0.4 I L ; I j 0.6 X* (x/Lx) I_I l i 0.8 i A • u Figure 4.30 Comparison of rms value of the velocity fluctuation u normalized by u' with Casel n and Case 3. Chapter 4 L E S for Turbulent Flow through a Converging Section 1.2 74 r • A O o 64*64*96 LES 64*64*64 LES Measurement .3 0.8 ICO = 0.6 o 0.2 "0 o o 8 0.4 A • A • A • A • A A A A h 0.2 0.4 0.6 X* (x/Lx) 0.8 Figure 4.31 Comparison of rms value of the velocity fluctuation w normalized by u' with n Casel, and Case 3. 4.3.2 Effects of the different subgrid models To consider the effect of the different subgrid models for the converging section, the ADM model is used with the 64 x 64 x 64 mesh. The comparison of mean streamwise velocity along the centerline of the converging section among the different subgrid models is shown in Figure 4.32. The comparisons of resolved rms fluctuations u, v and w normalized by u at the inlet (u ) between in two different subgrid models are presented in Figures 4.33-4.35. From the results, there is not much difference for two different subgrid models. It is noted that there is a little difference from w in Figure 4.35. The reason is that different filter operations are used in non-uniform meshes at z direction. For the Smagorinsky model, there is an implicit filter which involes one cell. The explicit filter operation, which involes five cells in each direction, is used for the ADM model. The Smagorinsky model is therefore adequate for the converging section. Chapter 4 LES for Turbulent Flow through a Converging Section 75 The ratio of the SGS eddy viscosity to the molecular viscosity, V IV , is shown with the T Smagorinsky model along the centre of the converging section in Figure 4.36. The value is less than 6 in most of region. The maximum value reaches around 13 near the exit. This indicates that the effect of the different subgrid models is not large. 12 Smagorinsky model 10 •r-? ADM modal O Experiment ("j o h 6 • 2h 0.4 0.6 -i X*(x/Lx) 0.2 1 1 1 0.8 L. 1 Figure 4.32 Comparison of the mean velocity at the centerline of the converging section with different subgrid models. 12 . hO G • O o [~] 0.2 o o - Smaoorlnsky modal • ADM modal Exparimant 0.4 X*(x/Lx) 0.6 0.8 Figure 4.33 Comparison of rms value of the velocity fluctuation u normalized by u with different in subgrid models. 76 Chapter 4 LES for Turbulent Flow through a Converging Section 1.4 r- Smagorinsky model • ADM model Experiment v vwvvw ^ V \ h 01 • ' • 0.2 • 0.4 0.6 -i X*(x/Lx) i i 0.8 i i i_ Figure 4.34 Comparison of rms value of the velocity fluctuation v normalized by u with different in subgrid models. 1.4 r- 1.2 Smagorinsky modal 5 ADM modal 3 O • 0.8 * 0.6 o • Experiment (CO ° sf7 v O © VW i v: V n 0.4 1 I 11 I 0.2 0.2 0.4 0.6 0.8 X*(x/Lx) Figure 4.35 Comparison of rms value of the velocity fluctuation w normalized by u with different in subgrid models. Chapter 4 LES for Turbulent Flow through a Converging Section 77 o - o ratio of eddy viscosoty with the molecular viscosity o o v lv T o - cP° OOOOOOCDOOO 1 0 0.2 , , 0.4X* 0.6 (x/Lx) 0.8 1 Figure 4.36 Ratio of SGS eddy viscosity to molecular viscosity: V IV with the Smagorinsky T model. 4.4 Summary The large eddy simulation method has been used to compute the turbulent flow through an asymmetric converging section. An important feature of this flow is that the turbulence is extremely non-isotropic near the exit where the standard k - £ model fails to predict the correct kinetic energy. The proper inflow condition is challenging for present quantitative predictions. The most straightforward approach to get the inflow condition for the converging section is to start the calculation by completely modelling the rectifier tubes. This procedure would be very time consuming. A one channel inflow condition is first employed to compute the turbulent flow in the converging section. This approach is also used by Kaltenbach et al. (1999), where a channel flow is used to generate inflow conditions for a plane diffuser. The computed values of the mean velocity agree well with the measured values. For rms fluctuations, the trend of the calculation results is similar to that of the experimental results, but there are some differences between the numerical Chapter 4 L E S for Turbulent Flow through a Converging Section 78 results and the measured results from the inlet of the converging section to the first third of the way through the converging section. The use of a "two-channel" to represent the inflow condition is also tested. The different inflow conditions can indeed affect the fluctuation results. The "two-channel" inflow condition is better than the one channel inflow condition, but one problem evident in the results using the "two-channel" inflow conditions is that the computational streamwise fluctations are only half of the measured values at 0.03m from the inlet of converging section. An alternative generation of inflow conditions method is therefore proposed. The idea is based on the extracting/rescaling technique developed by Lund et al. (1998). The computed streamwise fluctations and mean velocity using this new inflow data agree reasonably well with the measured values at 0.03m from the inlet of converging section. For the results along the centerline of the converging section, they are almost the same as those with one channel inflow condition. The v' with the new inflow condition is better than that with one channel inflow condition. The grid dependence study also reveals a sensitivity to the inflow condition. The more fine mesh in the wall-normal direction can make the predicted rms values of the velocity fluctuations along the centerline of the converging section more similar to the measured values. The reason is that the value of the velocity fluctuation w at the inlet is closer to the measured values with a fine mesh in the wall-normal direction. The agreement of rms values of velocityfluctuationsbetween the computational results and the measured values is less satisfactory. The remaining disagreement is likely to be partly related to the coarse meshes in the wall-normal direction. As indicated by our grid dependence study, a finer mesh would probably improve the numerical results somewhat. Another possible reason for the apparent difference between measured and predicted results is that we use an effectively twodimensional periodic condition in the spanwise direction while in the experimental apparatus there are walls on both sides of the channel. Another explanation involves the accuracy of the experimental data. In the experiment, there is the possibility of different flow rates in the adjacent Chapter 4 LES for Turbulent Flow through a Converging Section 79 rectifier tubes. From Shariati (2002), we know that the velocity profile is different for the top and the bottom rectifier tube in the experimental measurements. This will lead to inflow conditions that are unique to the apparatus and impossible to model with precision. The numerical code has been used to simulate the turbulence through a converging section. The objective 3, described in chapter 1, is achieved. Chapter 5 5.1 Fiber Concentration and Orientation Introduction The numerical flow simulations have been performed with LES as described in chapters 3 and 4 for fully developed channel flow and turbulent flow through a converging section. In this chapter, a fiber model is coupled with LES to study the fiber concentration in turbulent channel flow and fiber orientation in a headbox. Screening is an important process in both chemical and mechanical pulping (Kumar 1991, Olson 1996, Gooding 1986, 1996, Gooding et al. 2001, Dong 2002). Understanding the fiber motion and fiber concentration in the flow approaching the screen plate is a key part of predicting the screen fractionation efficiency. In section 5.2, numerical results are given to simulate numerically the fiber concentration in a 3D turbulent channel flow that is similar in many ways to the flow upstream of the screen. Fiber orientation in the paper determines the quality of the paper and the study of fiber orientation in the paper machine headbox becomes more and more important. Zhang (2001) measured the fiber orientations along the centerline of an asymmetric headbox. He also predicted the fiber orientation in both symmetric and asymmetric headboxes using a fiber model, but only using the mean flow field obtained from the K - E model. His comparison between the experiment and predicted results showed that the simulated orientations tend to align more with the flow than the experimentally observed orientations because the turbulence was not taken into consideration in his calculation. Dong et al. (2002) partially considered the turbulence effect with a single fixed velocity field found through LES simulation. In section 5.3, the fiber calculation executes one step after which the flow is updated with LES. This means the flow distribution for fibers is changed at each 80 Chapter 5 Fiber Concentration and Orientation 81 time step to consider more realistic turbulence effects in the simulation. The summary is given in section 5.4. 5.2 Fiber concentration Olson (1996) experimentally determined the detailed distribution of fiber concentration for a range of fiber lengths. The channel velocity in the experiment is 7.1 m/s and the channel Reynolds number is about 71,000 based on the channel half-height and bulk velocity. Dyed nylon fibers of length 0.0031m, 0.0020m and 0.001 lm with diameter 40 microns were used. The fiber concentration was less than 5000 fibers/litre so that the fiber-fiber interaction was negligible. Since the fiber suspension is dilute, the fiber-fiber interaction and the effects of the fiber on the turbulent flow are neglected in the present simulation. If a particle with density p and radius a moves in a fluid with dynamic viscosity ft and friction velocity u , the Stokes number (Crowe et al. T 1985) can be defined as In T = p St = ^jh which is essentially the ratio of the particle response time 2 and some fluid time scale T =^(Rouson and Eaton, 1994, 8 is the channel halff height). If the Stokes number is very small, then the particle follows the flow exactly. If it is very large, the particle will not be affected by turbulent fluctuations. If the Stokes number is in the order of one, the particle will be affected by the turbulence and will not follow the mean flow exactly (Fessler, et al, 1994). In the current simulation, if a is taken as the typical long radius of the spheroid used to represent the fiber, the Stokes number St is about 1.9. So the fiber is then expected to be affected by the turbulence and may not exactly follow the mean flow. 5.2.1 Fiber model and wall model The first investigation of the motion of a rigid, neutrally buoyant, ellipsoidal particle in a Stokes flow was conducted by Jeffery (1922). The Jeffery theory could be used to describe the Chapter 5 Fiber Concentration and Orientation 82 motion of rigid fibers (Anczurowski and Mason 1967). Pulp fibers have high aspect ratios and can have considerable flexibility which could not be modelled well by Jeffery's theory (Mason 1954). A flexible fiber model based on Euler-Bernoulli beam bending theory was introduced in Lawryshyn (1996). The immersed boundary method was used in Stockie and Green (1998) to simulate the motion of flexible pulp fibers. Their work was restricted to two-dimensional simulations. A flexible fiber model was proposed by Ross and Klingenberg (1997). It was further developed by Dong (2002) to account for fiber-wall interaction. In this model, each fiber consists of N rigid spheroids connected through N — l joints (see Figure 5.1). The lengths a, b(b < a) are the major and minor axes of each spheroid. The fiber can bend and twist much like a real fiber because of the rotational freedom in each joint. The motion of the fiber is determined by solving each spheroid's translation and rotation equations which are derived from Newton's second Law and the law of moment of momentum. Fibers frequently touch the wall in pulp and paper equipment, so that a wall model that can efficiently deal with the fiber-wall interaction is needed. For a smooth wall, a two-dimensional wall model was developed in Olson (1996). The idea is that a reaction force normal to the wall is exerted on the fiber to stop the fiber passing through the solid wall, and the friction force tangential to the wall is proportional to the normal force on the fiber. Dong (2002) extended the Olson's model to a three-dimensional universal wall model, which can deal with the fiber interaction in any wall geometry. We will use the fiber model and the wall model developed by Dong (2002) in the following simulations. Chapter 5 Fiber Concentration and Orientation 83 Figure 5.1 A fiber consisting of spheroids connected through ball and socket joints (from Ross and Klingenberg(1997)). 5.2.2 Numerical results The schema of the channel flow is shown in Figure 5.2. The height of the channel is 0.02m (z direction), which is the about the same as that of Olson (1996). The numerical methods have been described in chapter 3. The comparison between numerical and experimental results for channel flow can be found in Dong, et al. (2003). In order to record the fiber positions, several stations are set along the channel length. The channel height is divided into equal height intervals. The concentration ratio CICav is calculated by the ratio of the number of fibers at each height interval at each station to the average number of the fibers across the channel height. 8000 rigid fibers of 0.001m, 0.002m, 0.003m are initially chosen across the channel with random height and random orientations such that the number of fibers is about the same at each height interval. Figure 5.2 Schema of the channel flow 84 Chapter 5 Fiber Concentration and Orientation After the flow field from LES is statistically steady, the fibers are put in. Then each LES time step calculation is followed by one step of the fiber motion. Figure 5.3 shows the comparison of concentration for 0.001m, 0.002m, 0.003m fibers. We can see that the concentration curve of 0.002m and 0.003m fibers are very close to the Olson's fitted curve. Olson's fitted concentration curve (Olson 1996) is z/I<l/3.22, (5.1) z/L>l/3.22. where L is the fiber length, z is the fiber center position in normal wall direction. The 0.001m fiber does not have the same curve. The grid size in x direction (0.0006m) is probably too coarse for 0.001m fiber simulation. 3 2.5 — — — - > O O o L E S : 2 m m fiber L E S : 3 m m fiber L E S : 1 m m fiber O l s o n ' s fitting c u r v e for his e x p e r i m e n t data 2 c 1.5 o 2 c s c o o 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 85 Chapter 5 Fiber Concentration and Orientation 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 (x ,y ,z ), polar angle c c c 0 < 6 < 7t, which is the angle between the fiber main axis and Z-axis, and the azimuthal angle 0 < <p < 7t, which is the angle between the Y-axis and the projection of the fiber main axis on X Y plane (see Figure 5.4). It is noted that one end of the fiber is not distinguishable from the other. Z e o X Figure 5.4 A fiber's initial position. Y Chapter 5 Fiber Concentration and Orientation 86 Eight thousand (8,000) rigid fibers of 0.003m length are initially chosen across 0.04 m region (0.0175 m<z< 0.0575 m) in the center of the inlet height with random height and random orientation. Random orientation is implemented by choosing the fiber angles 0 and (/) with random number generators. The angles 0 and 0 are not chosen randomly from uniform distributions 0e [0,?r] and 0e[O,;r], since the area element dQ.-s\n0d(pd0 is a function of 0 on the surface of a unit sphere. So the angles are selected with the following formula (Zhang 2001): f? = cos (2fl,-l) _1 <f) = na (5.2) (5.3) 2 where a and a are random variables between [0, l ] . x 5.3.2 2 Numerical results The projections of the orientation of fibers can be obtained in three different planes. The fiber orientation was only measured in both the plane of the paper (x-y) and in the plane of contraction (x-z) at the central streamline of the converging section along the axis of the headbox in Zhang (2001). The measurements were taken at several points along the headbox as shown in Figure 5.5. So the projections on two of the planes, x-y and x-z planes, are considered here. The fiber K orientation angle K CC(—— <GC< —), either on the x-y or on the x-z projection plane, is defined to be the angle between the projection of the fiber axis on that plane and the machine direction (x axis). Figure 5.6 shows the initial fiber orientation distributions in the computation. The horizontal axis represents the orientation angle OC. The vertical axis stands for the statistical probability density p(a), such that: n_ 2 J p(a)da = 1 (5.4) Chapter 5 Fiber Concentration and Orientation 87 If the fiber orientations are uniform, then p((X) = l//r(=0.318). It can be seen that the initial random fiber orientation distribution is almost uniform from Figure 5.6. Figures 5.7-5.13 show the fiber orientation distributions at different stations. Most of the alignment occurs near the end of converging section, beyond x = 0.221m, where the velocity gradient is the highest. The LES results predict the fiber orientation reasonably well when compared to the experimental data from Zhang (2001). From Figures 5.10-5.12, the k-£ simulated fiber orientations tend to align more with the flow compared to experimental and LES results. At the exit where x = 0.3 Im, both k — £ and LES predict that the fibers are highly aligned in the flow direction in x — z plane as shown in Figure 5.13. But the peak value of LES is lower than that of k — £. Unfortunately there is no experimental data of fiber orientation at the exit in the x — z plane, because the channel is too narrow and the flow speed is too high at that location to get clear fiber images (Zhang 2001). From Figure 5.13, the fiber alignments in the x - z plane are stronger than that in the x — y plane. This phenomenon agrees with the observation in Zhang (2001) and Olson (2002). The statistical probability density pice) gives little information except for the trend of the fiber alignments in the flow direction. The first moment (^flr,./?(ar.)A<2r,.), second moment (^ocf/?(a,)Aa,.) and third moment i < ( ^afp(Cr;,)Acir ) of the probability distribution from numerical results and experimental i measurements are shown in Table 5.1 and Table 5.2. The LES results are somewhat closer to the experimental data than k — £ results in the x — z plane from Table 5.1. The first moment ( ^\ a p{(X^)tS.a ) at different stations in the x-z y i i i plane and the x — y plane are presented in Figures 5.14-5.15. The first moment indicates the mean angle value of the probability density pipe). The LES results of the first moment are closer to the experimental data Chapter 5 Fiber Concentration and Orientation 88 than the k — £ results near the exit in Figure 5.14. The first moment values from k — £ simulations in the x — y plane are better than thosefromLES simulations in Figure 5.15. The second moment (^ccf ) Aof .) at different stations in the x - z plane and the ( x— vplane are shown in Figures 5.16-5.17. The second moment values from LES simulations are closer to the experimental data than k — £ simulations in the x — zplane and thex— vplane from Figures 5.16-5.17. The second moment indicates that the fiber is whether alignment in the flow direction or not. This is important for the design of a headbox. The second moment is decreasing along the converging section from station 0.192m to the exit. The fiber is close to random at the entrance and that gives a large value of the second moment. As the fiber gets closer to the exit it gets aligned with the flow. Turbulence is counteracting this alignment. The fiber angle is affected by the du strain rate in the flow. The normal strain rate —— is large near the exit of the converging section. It is ox not surprising that the fibers tend to align themselves in the flow direction. There are only mean dU quantities in the k — £ simulations. The mean strain rate dominates in the flow along the center dx of the converging section near the exit. For LES simulations, the flow contains two parts: the mean dw part and the resolved fluctuation part. The fluctuation velocity gradient —— is large as seen in dx dw du Figure 4.24 from chapter 4. The shear strain rate —— + —— is not a negligible value. This shear strain dx dz du rate reduces the effect of the normal strain rate ——. That works against aligning the fibers in the dx flow direction. The turbulence tends to randomize the fiber orientation. So the LES simulations give the better comparison with experimental data for the second moment than the k — £ simulations. Chapter 5 Fiber Concentration and Orientation 89 The third moment ( ^T#f/?(«,. )Aar .) at different stations in the x - z plane and the ( x — y plane are presented in Figures 5.18-5.19. The third moment value indicates the skewness of the probability density p((X) . If the third moment is large, the p{(X) is very asymmetric. The third moment has small variation along the length of the converging section. There is some asymmetry but not very significant. The geometry and the flow are not completely symmetric. From Figures 5.18 and 5.19, the third moments from LES are close to experiment data at some stations. But there are some differences between LES results and experimental data at some stations. One obvious difference occurs at 0.262m station. One possible reason is that the fiber number is not large enough in the LES simulations near the exit. It is difficult to compare the results accurately. One reason is that the number of fibers used in the experiments and computations is not large enough. The initial 8,000 fibers are put into more than half channel height. The actual fibers at the center of converging section are less than 500. Much more computational time is needed with more fibers. From Zhang (2001), the experimental fiber number is less than 1500 at some stations. Another very likely reason is the lack of accuracy of flow measurement data as mentioned earlier in chapter 4. The problems in the experimental measurements discussed in chapter 4 have an influence on the fiber calculations as well. The changes in alignment near the exit are very abrupt. Accurate measurements near the exit of the converging section are difficult to make and the precise location is hard to define. Chapter 5 Fiber Concentration and Orientation 0.5 V) 0.6 InitJal fiber orientation 1/0.314 r- Initial fiber orientation 1/0.314 I0.3P s h-B-a FhB^ji R-SsSL f)-e- a 01 . h -1.5 -0.5 0 0.5 angle(radians) 1 1.5 -0.5 0 0.5 angle(radians) (a) (b) Figure 5.6 Initialfiberorientation distribution, (a) in x-z plane, (b) in x-y plane Chapter 5 Fiber Concentration and Orientation (a) (b) Figure 5.8 Fiber orientation distribution at x = 0.122m, (a) i n x-z plane, (b) in x-y plane Chapter 5 Fiber Concentration and Orientation (a) (b) Figure 5.10 Fiber orientation distribution at x = 0.192m, (a) in x-z plane, (b) i n x-y plane Chapter 5 Fiber Concentration and Orientation (a) (b) Figure 5.12 Fiber orientation distribution at x = 0.262m, (a) in x-z plane, (b) in x-y plane 94 Chapter 5 Fiber Concentration and Orientation angle(radians) angle(radians) (a) (b) Figure 5.13 Fiber orientation distribution at the converging section exit, (a) in x-z plane, (b) in x-y plane First Second Third First Second Third First Second Third moment moment moment moment moment moment moment moment (experiment) (experiment) (LES) (LES) (LES) (k-e) (k-e) (k-e) 0.045m -0.093 0.753 -0.125 0.027 0.769 0.029 0.122m -0.094 0.689 -0.118 -0.032 0.775 0.009 0.157m -0.039 0.670 -0.021 -0.050 0.675 -0.052 0.192m -0.025 0.623 -0.008 -0.067 0.652 -0.088 -0.001 0.663 -0.043 0.227m -0.074 0.594 -0.053 -0.028 0.611 -0.002 0.001 0.481 -0.035 0.262m -0.094 0.482 -0.085 -0.070 0.584 -0.205 0.006 0.295 -0.026 0.090 0.262 -0.016 0.033 0.082 -0.013 moment (experiment) exit Table 5.1 The first moment (^or./?(<2f .)Aa ), second moment ( t moment j/^af p{a )Aar ) of the probability distribution in x-z plane. i ( />(«;.)Aa;.) and third 95 Chapter 5 Fiber Concentration and Orientation Third First Second Third First Second Third moment moment moment moment moment moment moment moment (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 0.325 0.001 -0.027 0.538 -0.006 0.002 0.161 -0.004 First Second moment (experiment) exit -0.007 Table 5.2 The first moment (^ CC p{oc ) Acc ), second moment (^afp(a ) A or,) and third moment J I t t t i i ( >J oc]pia^Aot,) of the probability distribution in x-y plane. 0.1 r • LES O Experiment k-e 0.05 h • c o I • -0.05 h A A o • o • c -0.1 A 0.05 0.1 0.15 o 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 • LES A k-e O Experiment • 0.05 • c E 1 0 G o A o 0 • o o A o A A c • • • -0.05 Q A I I I I I I I I I I I I I I 0.1 0.05 0 I I I I I I I I I I I 0.2 0.15 I I I I I I I 0.3 0.25 Figure 5.15 The first moment (]T GC p(0C )Aa ) at different stations in x-y plane. i r • A O 0.8 l t LES k-e Experiment • 0 | 0.6 O £ •o 8 0 A O 0.4 0.2 A 1 1 1 1 1 1 0.05 * 0.1 1 1 1 1 I 1 1 1 1 0.15 * 1 1 0.2 1 1 I 1 1 1 0.25 1 * 0.3 X Figure 5.16 The second moment ( ^ a]pipe,)Aa ) at different stations in x-z plane. i Chapter 5 Fiber Concentration and Orientation r • : • 0.8 LES A k-e Q Experiment o : • c • • | 0.6 O c c A • • E -o 8 0.4 9 A (0 o 0.2 1 1 1 1 1 1 1 1 1 1 1 1 0.05 1 1 1 1 1 0.1 0.15 ' 0.2 ' 1 1 1 1 1 1 1 1 1 11 0.25 0.3 X Figure 5.17 The second moment (^T a]p{a )Aa ) t at different stations in x-y plane. t • 0.1 | - LES k-e A Experiment o • O | f A A n o A i-o,- i - -0.2 -0 31 ' 0 • i i i i I i 0.05 ' I' I i 0.1 i i i I i 0.15 ' i i I 0.2 i i I i 0.25 i i i I i 0.3 i X Figure 5.18 The third moment ( ^ a f p(a )Aa ) i i at different stations in x-z plane. 98 Chapter 5 Fiber Concentration and Orientation 0.04 • 0.02 O G C o O o o o 0 A S-0.02 E o E-0.04 A A • 1 ~-0.06 r- • A -0.08 - o -0.1 -0.12, ' 0.05 k-e Experiment • 0.1 0.15 0.2 0.25 0.3 Figure 5.19 The third moment ( ^ o f p{CC )AflT ) at different stations in x-y plane. i 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 deviationfromthe 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 predicted orientation distributions with LES are better than the k-e 100 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 6.1 Simulation of a Typical Industrial Headbox 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 L E S 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 6.2.1 Flow simulation 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 102 is essentially two-dimensional. 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 L E S 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 m / s (per unit width) which is close to the maximum 3 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 "fivechannel" 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 the "five-channel" inflow condition showing the details of one of the five identical channels 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 velocityfluctuationu at the inlet. Chapter 6 Simulation of a Typical Industrial Headbox w" (m/s) Figure 6.7 Comparison of rms value of the velocity fluctuation w'at the inlet. 107 Chapter 6 Simulation of a Typical Industrial Headbox 6.2.3 108 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 stationfromFigures 6.10-6.12. Thefluctuationsfrom"five-channel" inflow conditions grow more quickly, presumably due to the energy produced by the more nonuniform inlet mean velocity profile that is present in the "five-channel" case. The mean streamwise velocity profilesfromthe two inflow conditions are shown in Figure 6.13 at the second station CD derivedfromuse of the ADM model. Not surprisingly, there is little difference between these two distributions. The resolved rmsfluctuationsu , v and w are presented in Figures 6.14-6.16 at the second station CD with the ADM model. Thefluctuationsfrom "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 lengthfromthe "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 profilesfromthe two inflow conditions are shown in Figure 6.21 at the outlet. Essentially identical mean velocity profiles are obtained at the outletfromthe two different inflow conditions. The resolved rmsfluctuationsu , 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 thatfromone 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 interpolatedfromthe 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 inlet A D C E Figure 6.8 Comparison positions along the headbox. U (m/s) outlet Chapter 6 Simulation of a Typical Industrial Headbox Figure 6.9 Comparison of the mean velocity at the first station AB. n o D © D D o 0.8 • 0.6 One channel inflow ADM N Five channels Inflow ADM 0.4h 0.2 h 0.05 0.1 0.15 u" (m/s) c? F c i f t ^ P . D 0.2 0.25 Figure 6.10 Comparison of rms value of the velocity fluctuation u at the first station AB. 1 r- CP O 0.8 0.6 • • N fj One channel inflow ADM O Five channels inflow ADM 0.4 CK ^ 0.2 o 0.02 • I i?^Q*A' i l l rr 0X40.06 0.08 ^ 0.1 • (m/s) I I 0.12 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 o 0.8 0.6 • N 0.4 r- 0.2 h One channel Inflow ADM Five channels inflow ADM o 0 J • -Qo Q3 0.02 0.04 u" 0.06 0.08 (m/s) 0.1 Figure 6.14 Comparison of rms value of the velocityfluctuationu at the second station CD. 1 r- 0.8 0.6 • O M One channel inflow ADM Five channels Inflow ADM 0.4 0.2 oi fJ i i < i 0.02 I Q i i 0.04 ii i /<mv v" (m/s) 0.06 0.08 i i i i Ol i . i i 012 i i i i 014 Figure 6.15 Comparison of rms value of the velocityfluctuationv at the second station CD. Chapter 6 Simulation of a Typical Industrial Headbox Figure 6.17 Comparison of the mean velocity at the third station EF. 1 Chapter 6 Simulation of a Typical Industrial Headbox 0.8 0.6 One channel Inflow ADM O N Five channels Inflow ADM 0.4 0.2 -i i i i 0.02 i i i i 0.04 i nn "1X06 H06 1 i i i u' (m/s) © Q i 0.08 L 0.1 Figure 6.18 Comparison of rms value of the velocity fluctuation u at the third station EF. o• 11- 0.8 c oo O Q 0° D 0.6 ,8 N ° One channel inflow ADM O Five channels inflow ADM 0.4 0.2 pl fj I i irlni 002 i i I i O04 i ClaCl (j:06 J 0.08 v*(m/s) I I I I I L_I_J I 0.1 0.12 l_ 0.14 Figure 6.19 Comparison of rms value of the velocity fluctuation v at the third station EF. Chapter 6 Simulation of a Typical Industrial Headbox Figure 6.20 Comparison of rms value of the velocity fluctuation w at the third station EF. U (m/s) Figure 6.21 Comparison of the mean velocity at the outlet. Chapter 6 Simulation of a Typical Industrial Headbox u" (m/s) Figure 6.22 Comparison of rms value of the velocity fluctuation u at the outlet. v" (m/s) Figure 6.23 Comparison of rms value of the velocity fluctuation v at the outlet. Chapter 6 Simulation of a Typical Industrial Headbox w' (m/s) 118 2 Figure 6.24 Comparison of rms value of the velocityfluctuationw at the outlet. U (m/s) Figure 6.25 Comparison the mean velocity profiles at the outlet with the Smagorinsky model. Chapter 6 Simulation of a Typical Industrial Headbox t (m/s) Figure 6.26 Comparison of rms value of the velocity fluctuation u at the outlet with the Smagorinsky model. Figure 6.27 Comparison of rms value of the velocity fluctuation v at the outlet with the Smagorinsky model. Chapter 6 Simulation of a Typical Industrial Headbox w" (m/s) Figure 6.28 Comparison of rms value of the velocity fluctuation W at the outlet with the Smagorinsky model. v'/u" Figure 6.29 Comparison of the ratio v I u at different positions. Chapter 6 Simulation of a Typical Industrial Headbox w'/u' Figure 6.30 Comparison of the ratio w' lu' at different positions. v7u Figure 6.31 Comparison of the ratio v I u at the outlet with the two different inflow conditions. Chapter 6 Simulation of a Typical Industrial Headbox w'/u 1 2 2 1 Figure 6.32 Comparison of the ratio w I u at the outlet with the two different inflow conditions. 6.3 Summary Large eddy simulations are carried out for a commercial sized headbox. Two inflow conditions for the converging section, a one-channel inlet and a "five-channel" inlet, are investigated and two alternative subgrid models are compared. The two inflow conditions are significantly different at the inlet of the converging section. At one-sixth of the whole length position, this difference is reduced. The turbulence generated from "five-channel" inflow conditions is greater than that from one-channel inflow conditions at one third of the whole length position. At the outlet of the headbox, more turbulence is indeed generated using the "five-channel" inflow. However, the difference at the outlet between the two inflow conditions is not large. From the simulations, it is clear that the turbulence is close to isotropic at the first and second stations in most regions, but the anisotropy of the turbulence grows quickly near the outlet. At the outlet, the turbulence is very anisotropic in most regions. Chapter 6 Simulation of a Typical Industrial Headbox 123 One limitation of our calculations is that the periodic boundary condition is applied in the spanwise direction for all calculations. In the real commercial headbox, there are walls placed at regular intervals in the y direction in the five channel turbulence-producing tubes that form the entry for the "five-channel" condition in the present simulations. Since the length in the y direction is very large (8.64 m) compared to the length in the streamwise direction (0.665 m ), the effects of wall conditions for mean flow and turbulence would be expected to be small at the outlet. Removal of the periodic boundary conditions in the spanwise direction, while certainly possible, would require a major change in the calculation procedure and so is being left for future work. The turbulence generated in the tubes decays significantly in the converging section. Even the low level of turbulence intensity observed toward the exit of the converging section seems to be very important for the fiber distribution. This conclusion was reached in chapter 5, when computing fiber trajectories considering only the main flow as compared with computing them with turbulence fluctuations accounted for in a laboratory headbox calculation. The fibers move essentially with the same mean speed as the fluid so that the absolute turbulence intensity is of importance to them, not the apparently small relative turbulence intensity. The objective 5, described in chapter 1, is achieved. Chapter 7 Summary and Conclusions In this thesis, a numerical tool has been developed to simulate turbulence in a headbox. Using this tool, the flow through the headbox converging section has been investigated with parallel computing techniques and the large eddy simulation turbulence model. In addition, the statistical orientation of fibers passing through a headbox has been simulated by coupling the flow calculation with a suitable fiber model. The major findings of this thesis are summarized below: 7.1 Parallel computing In chapter 3, two methods, the direct method and the CGSTAB method, are employed to solve the Poisson equation. The MPI library is used to implement the message passing between processors. Parallel computation with the direct method is straightforward using ALL-TO-ALL (MPI command) to exchange the data. For the CGSTAB method, incomplete L U (ILU) is used as a preconditioner. A staircase parallelization algorithm is therefore used for the parallel implementation with the CGSTAB method. The advantage of this algorithm is that the serial and parallel versions have the same behavior with respect to convergence. The CGSTAB method is slower than the direct method, but it can be applied in more general cases. 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 or the parallel computation of a large eddy simulation. The computational time can be reduced much more with more processors. This makes numerical simulation of complex problems possible using parallel computation. 124 Chapter 7 Summary and Conclusions 7.2 125 LES for turbulentflowthrough a converging section 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. CGSTAB method was used to solve the Poisson equation for the converging section flow. The channel flow at Re = 395 was modelled first with three different SGS models T 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 Re = 1655, the channel flow was computed by using two 7 different mesh sizes. The trend in the calculations was similar to that of the low Reynolds number case. The Smagorinsky model with 64 x 64 x 64 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. Chapter 7 Summary and Conclusions 126 The proper inflow condition is challenging for the present quantitative predictions. The most straightforward approach to obtain the inflow conditions for the converging section is to start the calculation by completely modeling the upstream rectifier tubes. This procedure would be very time consuming. One channel inflow condition was employed to compute the turbulent flow in the converging section. This approach was also used by Kaltenbach et al. (1999), where a channel flow was used to generate inflow conditions for a plane diffuser. The computed values of the mean velocity agreed well with the measured values. For rms fluctuations, there were some differences between the numerical results and the measured results over the first third of the converging section length. The use of a "two-channel" inflow condition was also tested. Different inflow conditions can indeed affect the fluctuation results. The "two-channel" inflow condition was better than one channel inflow condition, but one evident problem in the results using the "two-channel" inflow condition was that the computational streamwise fluctations were only half of the measured values at 0.03m from inlet of converging section. A new generation of inflow conditions method was proposed. The idea was based on the extracting/rescaling technique developed by Lund et al. (1998). The computed streamwise fluctations and mean velocity using this new inflow data agreed reasonably well with the measured values at 0.03m from inlet of converging section. For the results along the centerline of the converging section, they were almost the same as those with one channel inflow condition. The values of v with the new inflow condition were better than those with the one channel inflow condition. The grid dependence study also revealed a sensitivity to the inflow condition. A finer mesh in the wall-normal direction can provide better agreement between the predicted and measured rms values of the velocity fluctuations along the centerline of the converging section. This is because the values of the velocity fluctuations V and w at the inlet were closer to the measured values when a fine mesh was used in the wall-normal direction. Chapter 7 Summary and Conclusions 127 The agreement of rms values of velocity fluctuations between the computational results and the measured values was less satisfactory. The disagreement can be partially related to the coarse mesh in the wall-normal direction. The grid dependence study suggested that a finer mesh would probably improve the numerical results. Another possible reason for the apparent difference between measured and predicted results was that a periodic condition in the spanwise direction was used while in the experimental apparatus there was a wall. Finally, it must be noted that in the experiment, there was also the possibility that flow rates in the adjacent rectifier tubes were different. From Shariati (2002), it is known that the velocity profile is different for the top and the bottom rectifier tube in the experimental measurements. This will lead to inflow conditions that are unique to the apparatus and impossible to model with precision. 7.3 Fiber orientation 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 were included in the fiber motion simulations. Random initial fiber orientations were set at the inlet of the channel. The statistical distribution for the orientation of a large number of fibers was evaluated at each station by computing the orientation of each single fiber along the central streamline. The LES results predicted the fiber orientation reasonably well when compared to the experimental data. The predicted orientation distributions near the exit of the converging section were better with LES than with those using the k — £ model for turbulence. 7.4 Simulation of a typical industrial headbox High Reynolds number large eddy simulations were carried out for a commercial sized headbox. Two inflow conditions for the converging section were investigated. Chapter 7 Summary and Conclusions 128 The two inflow conditions were significantly different at the inlet of the converging section. At the sixth of the whole length position, this difference was reduced. The turbulence generated from "five-channel" inflow conditions was greater than that from one channel inflow conditions at one third of the whole length position. At the outlet of the headbox, more turbulence was indeed generated using the "five-channel" inflow. However, the difference at the outlet between the two inflow conditions was not large. From the simulations, it was known that the turbulence was close to isotropic at the first and second stations in most regions, but the anisotropy of the turbulence grew quickly near the outlet. At the outlet, the turbulence was very anisotropic in most regions. The turbulence generated in the tubes decayed significantly in the converging section. Even the low level of turbulence intensity observed toward the exit of the converging section seemed to be very important for the fiber distribution. This conclusion was reached in chapter 5, when computing fiber trajectories considering only the mean flow as compared with computing them with turbulence fluctuations accounted for in a laboratory headbox calculation. From the above summary, the most important conclusions are listed below: • The techniques for parallel computation of a large eddy simulation are necessary and successful. The parallel efficiency can be higher than 80%. This makes numerical simulation of this complex problem possible. • For the channel flow, the dynamic model is very close to the Smagorinsky model. The Smagorinsky model with 64 x 64 x 64 mesh can provide a good representation of the channel flow. 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. Chapter 7 Summary and Conclusions • 129 Computations in an asymmetrical converging section show that the turbulence is very anisotropic, especially toward the exit of the converging section. Turbulence which survives the convergence will be strongly non-isotropic near the exit of any strongly converging section. • The standard k — £ model fails to predict the values of turbulence kinetic energy and therefore it should not be used in strongly converging sections. • For inflow conditions mimicking both one and two upstream channels, the mean flow computations agree with the experimental values. The computed turbulence fluctuations reproduce the experimental trends but the quantitative agreement is less good especially in the first third of the converging section length. • The LES coupled with a suitable fiber model to predict the orientation of nylon "fibers" in the converging section, gives statistical results similar to experimental data. • For the industrial headbox, multiple channels give somewhat higher turbulence at the exit compared to a single inlet channel, but the inlet turbulence has a limited effect on the turbulence at the exit. Chapter 8 Recommendations for Future Work This thesis has developed a parallel numerical tool to simulate the large scales of turbulence in a headbox using the LES model for turbulence. The fiber orientation in a headbox has been studied using this tool when coupled with a suitable fiber model. Further development can be done to provide a more accurate simulation of the turbulence inside the headbox. Some recommendations are as follows: The inflow condition for the headbox is still a challenging topic. In the present thesis a modified method that is similar to the random method is used to generate the inflow field. But the computational inflow condition still cannot represent the experimental inflow condition accurately. This may result in a difference between numerical and experimental results near the exit of a headbox. It would be better to simulate the complete rectifier tubes upstream of the converging section together with the converging section in the future. There may also be significant effects of non-uniform flow entering the rectifier tubes; there is some suggestion that the present experimental data is contaminated by such an effect. Sample calculations which introduce known non-uniform inflow conditions could quantify the downstream effects of upstream flow spatial non-uniformities. In the converging section simulations, the boundary condition assumed in the cross-machine direction is periodic. This assumption is economical and reasonable for focusing on the middle region of the converging section. The wall boundary condition in the spanwise direction in the incoming turbulence tubes and possibly in the converging section itself must be used if we consider the actual headbox converging section. In this situation, the fast FFT cannot be used. The parallel CGSTAB method in chapter 3 could be used to solve the Poisson equation. This would require more CPU time however, which may be practical in the near future. 130 Bibliography Aidun, C. and Kovacs, A.E., Hydrodynamics of the forming section: the origin of non-uniform fiber orientation, TAPPI Journal, 78(11), pp. 97-106, 1995. Anczurowski, E. and Mason, S.G., The kinetics of flowing dispersions. II. Equilibrium orientations of rods and discs (experimental), J. Colloid Interface Sci. 23, pp.533-546, 1967. Banhakavi, V. S. and Aidun, C. Analysis of turbulent flow in the converging zone of a headbox, TAPPI Engineering/Process and Product Quality Conference & Trade Fair, pp.1135-1154, 1999. Bastian, P., and Horton, G., Parallelization of robust multigrid methods: ILU factorization and frequency decomposition method, SIAM. J. Stat. Comput. 12, pp.1457-1470, 1991. Bonnet, J.P., Delville, J., Druault, P., Sagaut, P. and Grohens, R., Linear stochastic estimation of LES inflow conditions. (Advances in DNS/LES, Liu, C. and Liu, Z. eds.) Greyden Press, pp. 341348, 1997. Britt, K. W., Handbook of pulp and paper technology, Van Nostrand Reinhold, 1970. Burden, R.L. and Faires, J.D., numerical analysis, (5 Ed.). PWS Publishing Company, 1993. th Carati, D., Winckelmans, G. and Jeanmart, H . On the modelling of the subgrid-scale and filteredscale stress tensors in large-eddy simulation, J. Fluid Mech. 441, pp.119-138,2001. Chorin, A.J., Numerical solution of the Navier-Stokes equations, Math. Comput., 23, pp. 341-354, 1968. Chowdhury, D., Introduction to the Renormalization Group method and turbulence modelling, Fluent inc., TM-107 (cross-reference), 1993. Clark, J.H., Ferziger, J. H . and Reynolds, W. C , Evaluation of subgrid-scale models using an accurately simulated turbulent flow, J. Fluid. Mech., 91, pp. 1-6, 1978. Crowe, C.T., Gore, R. A . and Troutt, T. R., Particle dispersion by coherent structures in free shear flows, Part. Sci. Tech. 3, pp.149-158, 1985. 131 Bibliography 132 Dahl, H. K. and Weiss, H . G., A new hydraulic principle for headboxes, TAPPI Journal, 58(11), pp. 72-77, 1975. Deardodff, J.W. Three dimensional numerical study of turbulence in an entrained mixing layer. In AMS Workshop in Meteorology, American Meteorological Society, 1973. Deardorff, J.W., A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers. J. Fluid Mech. 41, 1970. Dong, S., Modelling of Fiber motion in pulp and paper equipment, Ph.D. thesis, University of British Columbia, 2002. Dong, S., Feng, X., Salcudean M . and Gartshore, I., Concentration of pulp fibers in 3D turbulent channel flow, Inter. J. Multiphase Flow, vol. 29, 1, pp. 1-21, 2003. Dong, S., Feng, X., Salcudean, M . , Gartshore, I. and Shariati, M . , Turbulence and Fiber Orientation in the Converging Section of a Paper-machine Headbox, The 4th ASME/JSME/KSME Symposium on Computational Techniques for Fluid/Thermal/Chemical Systems with Industrial Applications, Vancouver, Canada, 2002. Ferziger, J. and Peric, M . Computational methods for fluid dynamics, Springer Verlag, Second Edition, 1999. Fessler, J. R., Kulick, J. D. and Eaton, J. K., Preferential concentration of heavy particles in a turbulent channel flow, Phys. Fluids 6, 11, pp.3742-3749, 1994. Gavelin, G., Paper machine design and operation- descriptions and explanations-, Angus Wilde Publications Inc., 1998. Germano, M . , Poimelli, U., Cabot, W.H., and Moin, P. A dynamic subgrid scale eddy viscosity model. Phys. Fluids A, 3(7), pp.1760-1765, 1991. Ghosal, S., An analysis of numerical errors in large-eddy simulation of turbulence, J. comput. Phys., 125, pp. 187-206, 1996. Bibliography 133 Ghosal, S., Lund, T., Moin, P. and Akselvoll, K., A dynamic localization model for large eddy simulation of turbulent flows. J. Fluid Mech., 286, pp. 229-255, 1995. Gooding, R. W., The passage of fibres through slots in pulp screening, M. A. Sc thesis, University of British Columbia, 1986. Gooding, R. W., Flow Resistance of Screen Plate Apertures, Ph. D thesis, University of British Columbia, 1996. Gooding, R. W., Kerekes, R. J. and Salcudean, M., The Flow Resistance of Slotted Apertured in Pulp Screens, The Science ofPapermaking, 1, pp. 287-338, 2001. Golub, G., H. and van Loan, C. F., Matrix Computations (The Johns Hopkins University Press, New York), 1993. Gullbrand, J. and Chow, F. K., The effect of numerical errors and turbulence models in large-eddy simulations of channel flow, with and without explicit filtering, J. Fluid Mech., 495, pp. 323-341, 2003. Gustafsson, I., A class of first order Factorization methods, BIT, 18, pp. 142-156, 1978. Hamalainen, J., Mathematical modelling and simulation of fluid flows in the headbox of paper machines. Ph.D. Thesis, University of Jyvaskyla, 1993. He, P. and Salcudean, M., A numerical method for 3D viscous incompressible flows using nonorthogonal grids, Int. J. Numer. Methods in Fluids, 18, pp. 449-469, 1994. He, P., Bibeau, E. L., Hua, L., Salcudean, M . and Gartshore, I., Fluid dynamics of the flow distribution in a headbox, CPPI conference, Montreal, 1998. Horiuti, K. Comparison of conservative and rotational forms in large eddy simulation of turbulent channel flow. J. Comp. Phys. 71, 1987. Hackbusch, W., Multi-grid methods and applications, Springer, Berlin, 1985. Bibliography 134 Hua, L, Bibeau, E. L., He, P., Salcudean. M. and Gartshore, I., Flow distribution in a hydraulic headbox." TAPPI conference, Anaheim, 1999. Hua, L., He, P., Salcudean, M., Gartshore, I. and Bibeau, E. L., Turbulent flow in hydraulic headbox, TAPPI conference, Vancouver, 2000. Jeffery, G. B., The motion of ellipsoidal particles immered in a viscous fluid. Proc. Royal Soc, A102, pp.161-179, 1922. Joubert, W., Oppe, T., Janardhan, R. and Dearholt, W., Fully Parallel Global M/ILU Preconditioning for 3-D Structured Problems. 2000. Kaltenbach,H., Towards a near-wall model for LES of a separated diffuser flow. Annual Research Briefs, Center for Turbulence Research, NASA Ames/Stanford Univ., pp.255-265, 1998. Kaltenbach,H., Fatica, M., Mittal, R., Lund, T.S. and Moin, P., Study of flow in a planar asymmetric diffuser using large-eddy simulations, J. Fluid Mech., 390, pp. 151-185, 1999. Kim,J. and Moin,P., Application of fractional-step method to incompressible Navier-stokes equations. J. Comp. Phys., 59, pp.308-323, 1985. Kim,J., Moin,P. and Moser, R. D., Turbulence statistics in fully developed channel flow at low Renolds number, J. Fluid Mech., 177, pp. 133-166, 1987. Kumar, A., Passage of fibers through screen apertures. Ph.D. thesis, University of British Columbia, 1991. Lawryshyn, Y. A., Statics and dynamics of pulp fibers, Ph.D. thesis, University of Toronto, 1997. Lee, S., Lele, S. K. and Moin, P., Simulation of spatially evolving turbulence and the applicability of Taylor's hypothesis in compressible flow, Phys. Fluids A 4, pp. 1521-1530, 1992. Lesieur,M. and Metais, O. New trends in large-eddy simulations of turbulence, Annu. Rev. Fluid. Mech. 28,pp.45-82,1996. Lilly,D.K. On the application of the eddy viscosity concept in the inertial sub-range of turbulence, NCAR Manuscript No. 123,1966. Bibliography 135 Lilly,D.K., A proposed modification of the Germano subgrid scale closure method, Phys. Fluids A 4, pp.663-665, 1992. Lund, T. S., On the use of discrete filters for large eddy simulation. Annual Research Briefs, Center for Turbulence Research, NASA Ames/Stanford Univ., pp.83-95, 1997. Lund, T. S., Wu, X., and Squires, K. D., generation of turbulent inflow data for spatially-developing boundary layer simulations, J. Comp. Phys., 140, pp.233-258, 1998. Loewen, S. R., Fiber orientation optimization, Pulp Paper Can., 98(10), pp67-71, 1997. Mason, S.G., Fiber motions and flocculation, TAPPI Journal, 37(11), pp.494-501, 1954. Moin, P. and Kim, J. Numerical investigation of turbulent channel flow, J. Fluid Mech. 118, pp.341377, 1982. Moser, R. D., Kim, J. and Mansour, N. N., Direct numerical simulation of turbulent channel flow up to R e = 590, r Phys. Fluids, 11, pp.943-945, 1999. Najjar, F.M. and Tafti, D. K., Study of discrete test filters and finite difference approximations for the dynamic subgrid-scale stress model, Phys. Fluids, 8, pp. 1076-1088, 1996. Olson, J. A., The effect of fiber length on passage through narrow apertures, Ph.D. thesis, University of British Columbia, 1996. Olson, J. A., Analytic estimate of the fiber orientation distribution in a headbox flow, Nordic Pulp Paper Res. J. ,17 (3), pp. 302-306, 2002. Olson, J. A., Frigaard, I., Chan, C. and Hamalainen, J., Modelling a turbulent fiber suspension flowing in a plannar contraction: The one-dimensional headbox, Inter. J. Multiphase Flow, vol. 30, 1, pp. 51- 66, 2004. Parsheh, M . and Dahlkild, A.A. Modelling the flow around elastic guiding vanes and turbulence in a two-dimensional contraction, TAPPI Engineering/Process and Product Quality Conference & Trade Fair,pp.l433-1452, 1999. Bibliography 136 Peel, J. D., Paper science and paper manufacture, Angus Wilde Publications, 1999. Piomelli, U., Moin, P. and Ferziger, J.H., Model consistency in large eddy simulation of turbulent channel flows. Phys. Fluids. 31, 1988. Pougatch, K., Salcudean, M . , Gartshore, I. and Abdullah, Z., A computational analysis of a headbox flow: a case study, Practical Papermaking conference, Milwaukee, WI, 2005. Rai, M . M . and Moin, P., Direct numerical simulation of transition and turbulence in a spatially evolving boundary layer, J. Comp. Phys. 109, pp. 169-192, 1993. Reynolds, W. C , The potential and limitations of Direct and Large Eddy Simulations. Whither Turbulence? Turbulence at the Crossroads, J. L . Lumley, Ed., vol.357, Lecture Notes in Physics, Ithaca, New York, pp. 313-343, 1989. Rogallo,R.S. and Moin, P. Numerical simulation of turbulent flows. Ann. Rev. Fluid Mech. 16, pp.99-137, 1984. Rosenfeld, M . , Kwak, D. and Vinokur, M . , A fractional step solution method for the unsteady incompressible Navier-Stokes equations in generalized coordinate systems, J. Comp. Phys. 94, pp. 102-137, 1991. Ross, R. F. and Klingenberg, D. J., Dynamic simulation of flexible fibers composed of linked rigid bodies, J. Chem. Phys. 106 (7), pp.2949-2960, 1997. Rouson, D. W. I. and Eaton, J. K., Direct numerical simulation of turbulent channel flow with immersed particles, The 1994 A S M E Fluids Engineering Division Summer Meeting, Lake, 1994. Sagaut, P., Large eddy simulation for incompressible flows, Springer, 2001. Schumann, U., Subgrid scale models for finite difference simulations of turbulent flows in plane channels and annuli, J. Comp. Phys. 18, pp.376-404, 1975. Shariati, M . R., Bibeau, E., Salcudean M . and Gartshore, I., Numerical and Experimental Models of Flow in the Converging Section of a Headbox, TAPPI Papermakers Conference, Vancouver, Canada, 2000. Bibliography 137 Shariati, M . R., Bibeau, E., Salcudean, M . , Gartshore, I., Zhang, X. and Abdullah, Z. Computational Investigation of the Flow through a Headbox, TAPPI Engineering and Finishing Conference, San Antonio, 2001. Shariati, M . R., Experimental and mathematical modelling of flow in headboxes, Ph.D. thesis, the University of British Columbia, (2002). Shimizu, T. and Wada,K. Computer simulation of measurement of flow in a headbox. Proceedings of the Pan Pacific Pulp and Paper technology conference, pp. 157-165, 1992. Smagorinsky , J. General circulation experiments with the primitive equations. Mon. Weather Rev. 91(3), pp.99-164, 1963. Smook, G. A., Handbook for pulp and paper technologists, Angus Wilde Publications, 1992. Spalart, P. R., Theoretical and numerical study of a three-dimensional turbulent boundary layer, J. Fluid Mech., 205, pp.319-340, 1989. Speziale, C. G , On nonlinear k-l and k-e models of turbulence, J. Fluid Mech., 178, pp.459-475, 1987. Stockie, J. M . and Green, S. I., Simulating the motion of flexible pulp fibers using the immersed boundary method, J. Comp. Phys. 147, pp.147-165, 1998. Stolz, S., and Adams, N.A., An approximate deconvolution procedure for large-eddy simulation, Phys. Fluids 11, pp.1699-1701, 1999. Stolz, S., Adams, N.A., and Kleiser, L . , An approximate deconvolution model for large-eddy simulation with application to incompressible wall-bounded flows, Phys. Fluids 13, pp.997-1015, 2001. Ullmar, M . and Norman, B., Observation of fiber orientation in a headbox nozzle at low consistency, TAPPI Proceedings, Engineering and Papermakers Conference, Anaheim, pp.865-869, 1997. Ullmar, M . , On fiber orientation mechanisms in a headbox nozzle, Master's thesis, Royal Institute of Technology, Stockholm, Sweden, 1998. Bibliography 138 Van der Vorst, H . A., A fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems, SIAMJ. Sci. Statist. Comput. 13, pp.631-644, 1992. Vuik, C , van Nooyen, R.R.P., and Wesseling, P., Parallelelism in ILU-preconditioned GMRES, Parallel Computing, 24, pp. 1927-1946, 1998 Wei, T. and Willmarth, W. W., Reynolds-number effects on the structure of a turbulent channel flow. J. Fluid Mech., 204, pp.57-95, 1989. White, F. M., Fluid Mechanics, forth edition, McGraw-Hill, Inc., 1998. Wilcox, D. C , Turbulence modelling for CFD, Second Edition, DCW Industries, La Canada, CA, 1998. Yakhot, A., Orszag, S. A., Yakhot, V. and Israeli, M . , Renormalization group formulation of largeeddy simulations, J. Sci. Comput. 4, p. 139, 1989. Zang, Y , Street, R. L . and Koseff, J. R., A dynamic mixed subgrid-scale model and its application to turbulent recirculating flows, Phys. Fluids A 5, pp. 3186-3196, 1993. Zhang, X., Fiber orientation in Headbox, M . A. Sc thesis, University of British Columbia, 2001. Zhou, Y , Brasseur, J. and Juneja, A., A resolvable subfilter-scale model specific to large-eddy simulation of under-resolvable turbulence, Phys. Fluids 13, pp.2602-2610, 2001. Appendix A CG and CGSTAB algorithms The C G algorithm is presented as follows: • Initialize by setting: k = 0, u = u , r = Q - Au • while || r ||> e do 0 jn 0 ,p =0,s -\O ,e i0 jn 0 0 k k= k+\ Pk= klh-i s P = k+PkP -\ z k k <* = k Kp - Pk) s A k k u =u ^+a p k k k k k= k-\~ k Pk r r a A end while The CGSTAB algorithm is given follows: Initialize by setting: k = 0, u =u ,r =Q-Au , 0 in 0 in «o = Po = h = 1 • while || r \\> e do k k=k+l Pk= o- k-x r r 139 p = 0, 0 = 0, s Appendix <°k = (Aft-iV(a*-iA-i) Pk k-\ + k (Pk-i ~ k-\A-\) = u a a solve the system: Mz = </> =Az k Yk w = =PkW -r ) k 0 k-x-rA r i f || w | | < £ 2 k= k-x-Yk r r z return end i f solve the system: My = w v = Ay a ={v-r )l{vv) k k k = k-i+r u z u k r =w-a v k end while k + ky a 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 L L W = A' , L =A L W L' -A'-L' UN - AN L S U'~ -L'U'- NJ IL, P S U ] E - A I L. E P 142 Appendix 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 Ry(v) =< u'jUj > , , where the operation < z > zJ be c/(v) and Reynolds stress tensor by is an average on the spanwise direction and time. Three random numbers u (z), i = 1,2,3 are generated. They satisfy the following condition: each i 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 ~ By not listed above are set zero. ^ > A>33 - V A 3 3 ~P 31 -Pl2 •
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Large eddy simulation of the turbulence and fiber motion...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Large eddy simulation of the turbulence and fiber motion in a headbox Feng, Xiaosi 2005
pdf
Page Metadata
Item Metadata
Title | Large eddy simulation of the turbulence and fiber motion in a headbox |
Creator |
Feng, Xiaosi |
Date Issued | 2005 |
Description | 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. 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. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2009-12-23 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080757 |
URI | http://hdl.handle.net/2429/17095 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2005-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2005-104899.pdf [ 14.27MB ]
- Metadata
- JSON: 831-1.0080757.json
- JSON-LD: 831-1.0080757-ld.json
- RDF/XML (Pretty): 831-1.0080757-rdf.xml
- RDF/JSON: 831-1.0080757-rdf.json
- Turtle: 831-1.0080757-turtle.txt
- N-Triples: 831-1.0080757-rdf-ntriples.txt
- Original Record: 831-1.0080757-source.json
- Full Text
- 831-1.0080757-fulltext.txt
- Citation
- 831-1.0080757.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080757/manifest